第9章：稀疏数据结构与可微编程

struct OctreeNode {
    union {
        OctreeNode* children[8];  // 内部节点
        float value;              // 叶节点
    };
    uint8_t childMask;           // 子节点存在标记
};

blocks = {block_0, block_1, ..., block_n}
block_map: (i,j,k) -> block_index or NULL

// Array of Structures (AoS)
struct Voxel {
    float density;
    vec3 velocity;
    float temperature;
};
Voxel voxels[N];

// Structure of Arrays (SoA)
struct VoxelData {
    float density[N];
    vec3 velocity[N];
    float temperature[N];
};

struct Block {
    float density[BLOCK_SIZE];
    vec3 velocity[BLOCK_SIZE];
    // 按块组织，块内SoA
};

struct alignas(64) AlignedBlock {
    float data[BLOCK_SIZE];
    char padding[PAD_SIZE];  // 确保大小是64的倍数
};

// 伪代码
for block in active_blocks:
    prefetch(next_block)  // 预取下一个块
    process(block)

根节点 (Root): 32³ = 32768 分支
 └─ 内部节点 (Internal): 16³ = 4096 分支
     └─ 叶节点 (Leaf): 8³ = 512 体素

template<typename ChildT>
class RootNode {
    struct RootData {
        ChildT* child;     // 子节点指针
        Coord origin;      // 子节点原点坐标
    };
    std::map<Coord, RootData> mTable;  // 稀疏哈希表
};

template<typename ChildT>
class InternalNode {
    static constexpr int DIM = 16;
    static constexpr int SIZE = DIM * DIM * DIM;

    union Child {
        ChildT* child;
        ValueType value;
    };

    Child mNodes[SIZE];
    Mask mChildMask;    // 4096位，标记是否为子节点
    Mask mValueMask;    // 4096位，标记是否为活跃值
};

template<typename ValueT>
class LeafNode {
    static constexpr int DIM = 8;
    static constexpr int SIZE = DIM * DIM * DIM;

    ValueT mValues[SIZE];
    Mask mValueMask;    // 512位，标记活跃体素
};

template<size_t SIZE>
class Mask {
    static constexpr size_t WORD_COUNT = (SIZE + 63) / 64;
    uint64_t mWords[WORD_COUNT];

public:
    bool isOn(size_t n) const {
        return mWords[n >> 6] & (uint64_t(1) << (n & 63));
    }

    void setOn(size_t n) {
        mWords[n >> 6] |= uint64_t(1) << (n & 63);
    }

    size_t countOn() const {
        size_t count = 0;
        for (auto word : mWords)
            count += __builtin_popcountll(word);
        return count;
    }
};

class InternalNode {
    bool isTileValue(size_t n) const {
        return !mChildMask.isOn(n) && mValueMask.isOn(n);
    }

    ValueType getTileValue(size_t n) const {
        return mNodes[n].value;  // 直接存储值，无需子节点
    }
};

template<typename NodeT>
class NodePool {
    struct Block {
        NodeT nodes[NODES_PER_BLOCK];
        Block* next;
    };

    Block* mHead;
    std::vector<NodeT*> mFreeList;

public:
    NodeT* allocate() {
        if (mFreeList.empty()) {
            allocateNewBlock();
        }
        NodeT* node = mFreeList.back();
        mFreeList.pop_back();
        return node;
    }
};

void LeafNode::fill(ValueType value) {
    __m256 val = _mm256_set1_ps(value);
    float* data = mValues;

    for (int i = 0; i < SIZE; i += 8) {
        _mm256_store_ps(data + i, val);
    }
}

void LeafNode::add(const LeafNode& other) {
    for (int i = 0; i < SIZE; i += 8) {
        __m256 a = _mm256_load_ps(mValues + i);
        __m256 b = _mm256_load_ps(other.mValues + i);
        __m256 sum = _mm256_add_ps(a, b);
        _mm256_store_ps(mValues + i, sum);
    }
}

template<typename TreeT>
class ValueIterator {
    using LeafIterT = typename TreeT::LeafIterator;
    using VoxelIterT = typename LeafNode::ValueOnIterator;

    LeafIterT mLeafIter;
    VoxelIterT mVoxelIter;

public:
    void operator++() {
        ++mVoxelIter;
        if (!mVoxelIter) {  // 当前叶节点遍历完
            ++mLeafIter;
            if (mLeafIter) {
                mVoxelIter = mLeafIter->beginValueOn();
            }
        }
    }
};

class NodeIterator {
    std::stack<NodePtr> mStack;

public:
    void advance() {
        NodePtr node = mStack.top();
        mStack.pop();

        // 深度优先遍历
        for (auto child : node->children()) {
            if (child) mStack.push(child);
        }
    }
};

template<int RADIUS>
class StencilIterator {
    static constexpr int SIZE = 2 * RADIUS + 1;
    ValueType mCache[SIZE][SIZE][SIZE];
    Coord mCenter;

public:
    void moveTo(const Coord& xyz) {
        // 智能缓存更新：只加载新进入的值
        updateCache(xyz - mCenter);
        mCenter = xyz;
    }
};

// Morton序遍历提高缓存命中率
void mortonTraverse(LeafNode& leaf) {
    static const int morton[8] = {0,1,4,5,2,3,6,7};

    for (int z = 0; z < 8; z += 2) {
        for (int y = 0; y < 8; y += 2) {
            for (int x = 0; x < 8; x += 2) {
                // 2x2x2块内按Morton序访问
                for (int i = 0; i < 8; ++i) {
                    int dx = morton[i] & 1;
                    int dy = (morton[i] >> 1) & 1;
                    int dz = (morton[i] >> 2) & 1;
                    process(leaf(x+dx, y+dy, z+dz));
                }
            }
        }
    }
}

void dilate(Tree& tree, int radius) {
    // 收集所有活跃体素的邻域
    std::vector<Coord> newActive;

    for (auto iter = tree.beginValueOn(); iter; ++iter) {
        Coord xyz = iter.getCoord();

        // 26邻域膨胀
        for (int dz = -radius; dz <= radius; ++dz) {
            for (int dy = -radius; dy <= radius; ++dy) {
                for (int dx = -radius; dx <= radius; ++dx) {
                    if (dx*dx + dy*dy + dz*dz <= radius*radius) {
                        newActive.push_back(xyz + Coord(dx,dy,dz));
                    }
                }
            }
        }
    }

    // 激活新体素
    for (const auto& coord : newActive) {
        tree.setValueOn(coord);
    }
}

void erode(Tree& tree, int radius) {
    std::vector<Coord> toDeactivate;

    for (auto iter = tree.beginValueOn(); iter; ++iter) {
        if (!hasFullNeighborhood(tree, iter.getCoord(), radius)) {
            toDeactivate.push_back(iter.getCoord());
        }
    }

    for (const auto& coord : toDeactivate) {
        tree.setValueOff(coord);
    }
}

void parallelDilate(Tree& tree, int radius) {
    // 1. 按叶节点分组
    std::vector<LeafNode*> leaves;
    tree.getNodes(leaves);

    // 2. 并行处理每个叶节点
    tbb::parallel_for(size_t(0), leaves.size(), [&](size_t i) {
        LeafNode* leaf = leaves[i];
        LeafNode dilated = *leaf;  // 复制

        // 叶节点内部膨胀
        dilateLeaf(dilated, radius);

        // 原子更新
        leaf->merge(dilated);
    });

    // 3. 处理跨叶节点边界
    processLeafBoundaries(tree, radius);
}

size_t countActiveVoxels(const Tree& tree) {
    size_t count = 0;

    // 利用层次结构避免遍历空区域
    tree.visitActiveBBox([&](const CoordBBox& bbox) {
        if (bbox.volume() == 1) {
            count++;  // 单个体素
        } else {
            // 瓦片值区域
            count += bbox.volume();
        }
    });

    return count;
}

CoordBBox computeBBox(const Tree& tree) {
    CoordBBox bbox;

    // 只需检查活跃叶节点的边界
    for (auto iter = tree.beginLeaf(); iter; ++iter) {
        bbox.expand(iter->getNodeBoundingBox());
    }

    return bbox;
}

## 9.3 稀疏MPM实现

将MPM与稀疏数据结构结合，可以显著提高大规模模拟的效率。稀疏MPM的核心挑战在于高效的粒子-网格传输和稀疏线性系统求解。

### 9.3.1 稀疏网格表示

#### MPM中的稀疏性特征

MPM模拟中的稀疏性来源于：

1. **材料分布**：物体只占据部分空间
2. **影响域局部性**：每个粒子只影响周围3×3×3的网格节点
3. **时间相关性**：活跃区域在时间步之间变化缓慢

**活跃节点定义：**
节点 $\mathbf{x}\_i$ 在时间步 $t$ 是活跃的，当且仅当：
$$\exists p \in \mathcal{P}: \|\mathbf{x}\_p^t - \mathbf{x}\_i\| < 2h$$
其中 $\mathcal{P}$ 是粒子集合，$h$ 是网格间距。

#### 稀疏网格数据结构

**分块稀疏网格：**
```cpp
template<int BLOCK_SIZE = 8>
class SparseMPMGrid {
    struct Block {
        // 节点数据（速度、质量等）
        Vec3f velocity[BLOCK_SIZE³];
        float mass[BLOCK_SIZE³];

        // 块元数据
        Vec3i baseCoord;  // 块的全局坐标
        uint64_t activeMask;  // 活跃节点掩码
    };

    std::unordered_map<Vec3i, Block*, GridHasher> blocks;

public:
    // 激活节点
    void activate(const Vec3i& coord) {
        Vec3i blockCoord = coord / BLOCK_SIZE;
        Block*& block = blocks[blockCoord];

        if (!block) {
            block = allocateBlock();
            block->baseCoord = blockCoord * BLOCK_SIZE;
        }

        Vec3i localCoord = coord % BLOCK_SIZE;
        int idx = localCoord.x + BLOCK_SIZE * 
                 (localCoord.y + BLOCK_SIZE * localCoord.z);
        block->activeMask |= (1ULL << idx);
    }
};

void updateActiveNodes(SparseMPMGrid& grid, 
                      const std::vector<Particle>& particles) {
    // 清除旧的活跃标记
    grid.clearActive();

    // 并行标记新的活跃节点
    tbb::parallel_for(size_t(0), particles.size(), [&](size_t p) {
        Vec3f pos = particles[p].position;
        Vec3i base = Vec3i(pos / grid.dx);

        // 标记3x3x3邻域
        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) {
                for (int k = 0; k < 3; ++k) {
                    grid.activate(base + Vec3i(i,j,k));
                }
            }
        }
    });
}

void particleToGrid(const std::vector<Particle>& particles,
                   SparseMPMGrid& grid) {
    // 1. 清零网格数据
    grid.zeroActiveNodes();

    // 2. 并行累积粒子贡献
    tbb::parallel_for(size_t(0), particles.size(), [&](size_t p) {
        const Particle& particle = particles[p];
        Vec3f xp = particle.position;
        Vec3i base = Vec3i(xp / grid.dx);
        Vec3f fx = xp / grid.dx - Vec3f(base);

        // 二次B样条权重
        Vec3f w[3];
        w[0] = 0.5f * (Vec3f(1.5) - fx).square();
        w[1] = 0.75f - (fx - Vec3f(1.0)).square();
        w[2] = 0.5f * (fx - Vec3f(0.5)).square();

        // 计算应力
        Mat3f stress = computeStress(particle);

        // 传输到3x3x3邻域
        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) {
                for (int k = 0; k < 3; ++k) {
                    Vec3i node = base + Vec3i(i,j,k);
                    float weight = w[i].x * w[j].y * w[k].z;

                    Vec3f xg = Vec3f(node) * grid.dx;
                    Vec3f momentum = particle.mass * 
                        (particle.velocity + particle.C * (xg - xp));

                    // 原子操作累积
                    grid.addMomentum(node, weight * momentum);
                    grid.addMass(node, weight * particle.mass);
                    grid.addForce(node, -particle.volume * 
                                       stress * gradWeight(fx, i,j,k));
                }
            }
        }
    });
}

struct LocalBuffer {
    static constexpr int BUFFER_SIZE = 27;  // 3x3x3
    Vec3f momentum[BUFFER_SIZE];
    float mass[BUFFER_SIZE];
    Vec3i indices[BUFFER_SIZE];
    int count = 0;

    void flush(SparseMPMGrid& grid) {
        for (int i = 0; i < count; ++i) {
            grid.addMomentum(indices[i], momentum[i]);
            grid.addMass(indices[i], mass[i]);
        }
        count = 0;
    }
};

void gridToParticle(SparseMPMGrid& grid,
                   std::vector<Particle>& particles) {
    tbb::parallel_for(size_t(0), particles.size(), [&](size_t p) {
        Particle& particle = particles[p];
        Vec3f xp = particle.position;
        Vec3i base = Vec3i(xp / grid.dx);
        Vec3f fx = xp / grid.dx - Vec3f(base);

        // 计算权重（同P2G）
        Vec3f w[3];
        computeWeights(fx, w);

        // 收集网格信息
        Vec3f newVelocity(0.0f);
        Mat3f velocityGradient(0.0f);

        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) {
                for (int k = 0; k < 3; ++k) {
                    Vec3i node = base + Vec3i(i,j,k);
                    float weight = w[i].x * w[j].y * w[k].z;

                    Vec3f vi = grid.getVelocity(node);
                    newVelocity += weight * vi;

                    Vec3f gradW = gradWeight(fx, i,j,k);
                    velocityGradient += outer(vi, gradW);
                }
            }
        }

        // 更新粒子状态
        particle.velocity = newVelocity;
        particle.C = velocityGradient * (4.0f / (grid.dx * grid.dx));
        particle.position += dt * newVelocity;
    });
}

struct SparsePattern {
    std::vector<std::set<int>> adjacency;

    void analyze(const SparseMPMGrid& grid) {
        // 遍历所有活跃节点对
        for (auto& [coord1, node1] : grid.activeNodes()) {
            int idx1 = grid.getIndex(coord1);

            // 检查27邻域
            for (int di = -1; di <= 1; ++di) {
                for (int dj = -1; dj <= 1; ++dj) {
                    for (int dk = -1; dk <= 1; ++dk) {
                        Vec3i coord2 = coord1 + Vec3i(di,dj,dk);
                        if (grid.isActive(coord2)) {
                            int idx2 = grid.getIndex(coord2);
                            adjacency[idx1].insert(idx2);
                        }
                    }
                }
            }
        }
    }
};

class SparseMPMMatrix {
    std::vector<float> values;
    std::vector<int> colIndices;
    std::vector<int> rowOffsets;

public:
    void construct(const SparsePattern& pattern) {
        int nnz = 0;
        rowOffsets.push_back(0);

        for (int row = 0; row < pattern.size(); ++row) {
            for (int col : pattern.adjacency[row]) {
                colIndices.push_back(col);
                values.push_back(0.0f);  // 待填充
                nnz++;
            }
            rowOffsets.push_back(nnz);
        }
    }

    // 稀疏矩阵-向量乘法
    void multiply(const float* x, float* y) const {
        #pragma omp parallel for
        for (int row = 0; row < numRows(); ++row) {
            float sum = 0.0f;
            for (int idx = rowOffsets[row]; idx < rowOffsets[row+1]; ++idx) {
                sum += values[idx] * x[colIndices[idx]];
            }
            y[row] = sum;
        }
    }
};

class SparsePCGSolver {
    // Jacobi预条件
    std::vector<float> diagonal;

public:
    void solve(const SparseMPMMatrix& A, const float* b, float* x,
               float tolerance = 1e-6f, int maxIter = 1000) {
        int n = A.numRows();
        std::vector<float> r(n), z(n), p(n), Ap(n);

        // 初始化：r = b - Ax
        A.multiply(x, Ap.data());
        for (int i = 0; i < n; ++i) {
            r[i] = b[i] - Ap[i];
        }

        // 应用预条件：z = M^{-1}r
        applyPreconditioner(r.data(), z.data());

        // p = z
        std::copy(z.begin(), z.end(), p.begin());

        float rz_old = dot(r.data(), z.data(), n);

        for (int iter = 0; iter < maxIter; ++iter) {
            // Ap = A * p
            A.multiply(p.data(), Ap.data());

            // α = (r·z) / (p·Ap)
            float pAp = dot(p.data(), Ap.data(), n);
            float alpha = rz_old / pAp;

            // x = x + α*p, r = r - α*Ap
            #pragma omp parallel for
            for (int i = 0; i < n; ++i) {
                x[i] += alpha * p[i];
                r[i] -= alpha * Ap[i];
            }

            // 检查收敛
            float residual = norm(r.data(), n);
            if (residual < tolerance) break;

            // z = M^{-1}r
            applyPreconditioner(r.data(), z.data());

            // β = (r·z)_new / (r·z)_old
            float rz_new = dot(r.data(), z.data(), n);
            float beta = rz_new / rz_old;

            // p = z + β*p
            #pragma omp parallel for
            for (int i = 0; i < n; ++i) {
                p[i] = z[i] + beta * p[i];
            }

            rz_old = rz_new;
        }
    }

private:
    void applyPreconditioner(const float* r, float* z) {
        #pragma omp parallel for
        for (int i = 0; i < diagonal.size(); ++i) {
            z[i] = r[i] / diagonal[i];
        }
    }
};

void mpmStep(SparseMPMGrid& grid, std::vector<Particle>& particles) {
    // Level 1: 块级并行
    tbb::task_group tg;

    // P2G阶段
    tg.run([&] {
        // Level 2: 粒子并行
        tbb::parallel_for_each(particles, [&](Particle& p) {
            // Level 3: SIMD
            processParticleP2G_SIMD(p, grid);
        });
    });

    // 网格更新
    tg.run([&] {
        grid.updateVelocities();
    });

    tg.wait();
}

void sortParticlesByMorton(std::vector<Particle>& particles) {
    // 计算Morton码
    std::vector<std::pair<uint64_t, int>> mortonCodes;
    mortonCodes.reserve(particles.size());

    for (int i = 0; i < particles.size(); ++i) {
        uint64_t code = computeMortonCode(particles[i].position);
        mortonCodes.push_back({code, i});
    }

    // 按Morton码排序
    std::sort(mortonCodes.begin(), mortonCodes.end());

    // 重排粒子
    std::vector<Particle> sorted(particles.size());
    for (int i = 0; i < particles.size(); ++i) {
        sorted[i] = particles[mortonCodes[i].second];
    }
    particles = std::move(sorted);
}

float computeAdaptiveTimestep(const SparseMPMGrid& grid,
                             const std::vector<Particle>& particles) {
    float maxVel = 0.0f;
    float maxAccel = 0.0f;

    // 并行计算最大值
    tbb::parallel_reduce(
        tbb::blocked_range<size_t>(0, particles.size()),
        std::make_pair(0.0f, 0.0f),
        [&](const auto& range, auto maxVals) {
            for (size_t i = range.begin(); i < range.end(); ++i) {
                float vel = particles[i].velocity.norm();
                float accel = particles[i].force.norm() / particles[i].mass;
                maxVals.first = std::max(maxVals.first, vel);
                maxVals.second = std::max(maxVals.second, accel);
            }
            return maxVals;
        },
        [](auto a, auto b) {
            return std::make_pair(
                std::max(a.first, b.first),
                std::max(a.second, b.second)
            );
        }
    );

    float dt_vel = grid.dx / maxVel;
    float dt_accel = std::sqrt(grid.dx / maxAccel);

    return 0.4f * std::min(dt_vel, dt_accel);
}

class GPUSparseMPMGrid {
    // 使用统一内存避免显式传输
    Block* blocks;  // cudaMallocManaged
    int* blockMap;  // 块索引映射

    __device__ Vec3f& velocity(int x, int y, int z) {
        int blockIdx = getBlockIndex(x, y, z);
        int localIdx = getLocalIndex(x, y, z);
        return blocks[blockIdx].velocity[localIdx];
    }
};

__device__ void atomicAddVec3(Vec3f* addr, Vec3f val) {
    // 利用warp内线程合作减少原子操作冲突
    int lane = threadIdx.x & 31;
    int leader = __ffs(__ballot_sync(0xffffffff, true)) - 1;

    if (lane == leader) {
        // 只有warp leader执行原子操作
        atomicAdd(&addr->x, val.x);
        atomicAdd(&addr->y, val.y);
        atomicAdd(&addr->z, val.z);
    }
}

struct DiffMPMState {
    // 前向值
    std::vector<Vec3f> positions;
    std::vector<Vec3f> velocities;
    std::vector<Mat3f> deformationGradients;

    // 梯度
    std::vector<Vec3f> dPositions;
    std::vector<Vec3f> dVelocities;
    std::vector<Mat3f> dDeformationGradients;
};

void diffMPMStep(DiffMPMState& state, const MPMParams& params,
                bool computeGradients = false) {
    // 前向传播
    particleToGrid(state, params);
    gridUpdate(state, params);
    gridToParticle(state, params);

    if (computeGradients) {
        // 反向传播
        gridToParticleBackward(state, params);
        gridUpdateBackward(state, params);
        particleToGridBackward(state, params);
    }
}

void diffImplicitSolve(DiffMPMState& state, const MPMParams& params) {
    // 前向：牛顿迭代
    for (int iter = 0; iter < maxIter; ++iter) {
        Mat Jacobian = computeJacobian(state, params);
        Vec residual = computeResidual(state, params);
        Vec delta = solve(Jacobian, -residual);
        state.positions += delta;

        if (residual.norm() < tolerance) break;
    }

    // 反向：通过伴随法计算梯度
    if (state.hasGradients) {
        // λ = (J^T)^{-1} * dL/dx_{t+1}
        Vec lambda = solveTranspose(Jacobian, state.dPositions);

        // dL/dθ = -λ^T * ∂g/∂θ
        state.dParams = -lambda.dot(computeParamDerivative(state, params));
    }
}

x_0 = x
x_1 = f_1(x_0)
x_2 = f_2(x_1)
...
y = x_n = f_n(x_{n-1})

v_n = ∂L/∂y
v_{n-1} = v_n · ∂f_n/∂x_{n-1}
v_{n-2} = v_{n-1} · ∂f_{n-1}/∂x_{n-2}
...
∂L/∂x = v_0 = v_1 · ∂f_1/∂x_0

void p2gBackward(DiffMPMState& state, const MPMParams& params) {
    // 清零粒子梯度
    state.clearParticleGradients();

    for (int p = 0; p < numParticles; ++p) {
        Vec3f xp = state.positions[p];
        Vec3i base = Vec3i(xp / dx);
        Vec3f fx = xp / dx - Vec3f(base);

        // 计算权重及其导数
        float w[3][3];
        float dw[3][3];
        computeWeightsAndDerivatives(fx, w, dw);

        // 从网格收集梯度
        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) {
                for (int k = 0; k < 3; ++k) {
                    Vec3i node = base + Vec3i(i,j,k);

                    // 位置梯度
                    Vec3f dMomentum = state.grid.dMomentum[node];
                    state.dPositions[p] += dMomentum * 
                        (dw[i] * w[j] * w[k] + w[i] * dw[j] * w[k] + 
                         w[i] * w[j] * dw[k]) / dx;

                    // 速度梯度
                    state.dVelocities[p] += dMomentum * 
                        (w[i] * w[j] * w[k]) * state.mass[p];

                    // 变形梯度的梯度（通过应力）
                    Mat3f dStress = computeStressGradient(state, p);
                    state.dDeformationGradients[p] += 
                        dStress : state.grid.dForce[node];
                }
            }
        }
    }
}

bool gradientCheck(const Function& f, const Vec& x, const Vec& grad,
                  float epsilon = 1e-5f, float tolerance = 1e-4f) {
    Vec numGrad(x.size());

    for (int i = 0; i < x.size(); ++i) {
        Vec xPlus = x;
        Vec xMinus = x;
        xPlus[i] += epsilon;
        xMinus[i] -= epsilon;

        float fPlus = f(xPlus);
        float fMinus = f(xMinus);

        numGrad[i] = (fPlus - fMinus) / (2 * epsilon);
    }

    float relError = (grad - numGrad).norm() / 
                     (grad.norm() + numGrad.norm() + 1e-8f);

    return relError < tolerance;
}

struct MaterialOptimizer {
    DiffMPMSimulator simulator;
    std::vector<ObservationData> observations;

    MaterialParams optimize(const MaterialParams& initial,
                          int maxIterations = 100) {
        MaterialParams params = initial;
        AdamOptimizer optimizer(learningRate = 0.01);

        for (int iter = 0; iter < maxIterations; ++iter) {
            float loss = 0.0f;
            MaterialParams gradParams;

            // 计算所有观测点的损失和梯度
            for (const auto& obs : observations) {
                auto simResult = simulator.simulate(params, obs.initialState);

                // 前向：计算损失
                loss += computeLoss(simResult.finalState, obs.targetState);

                // 反向：计算参数梯度
                auto dLoss = computeLossGradient(simResult.finalState, 
                                                obs.targetState);
                simulator.backward(dLoss);
                gradParams += simulator.getParamGradient();
            }

            // 更新参数
            params = optimizer.update(params, gradParams);

            if (loss < tolerance) break;
        }

        return params;
    }
};

void trajectoryOptimization(DiffMPMSimulator& sim,
                          const State& initialState,
                          const State& targetState) {
    std::vector<ControlInput> controls(T);

    for (int optIter = 0; optIter < maxOptIter; ++optIter) {
        // 前向模拟
        std::vector<State> trajectory = sim.rollout(initialState, controls);

        // 计算损失
        float loss = computeTrajectoryLoss(trajectory, targetState);

        // 反向传播（伴随法）
        std::vector<State> adjoints(T+1);
        adjoints[T] = computeTerminalGradient(trajectory[T], targetState);

        for (int t = T-1; t >= 0; --t) {
            // 状态伴随变量
            adjoints[t] = sim.dynamicsTranspose(trajectory[t], controls[t]) 

                         * adjoints[t+1]
                         + computeStateGradient(trajectory[t]);

            // 控制梯度
            Vec dControl = sim.controlJacobian(trajectory[t], controls[t])
                          .transpose() * adjoints[t+1];

            // 更新控制
            controls[t] -= learningRate * dControl;
        }
    }
}

class ShapeOptimizer {
    void optimizeShape(DiffMPMSimulator& sim,
                      const TargetShape& target) {
        LevelSet shape = initializeShape();

        for (int iter = 0; iter < maxIter; ++iter) {
            // 生成粒子
            auto particles = sampleParticles(shape);

            // 模拟变形
            auto result = sim.simulate(particles);

            // 计算形状误差
            float error = computeShapeError(result.finalShape, target);

            // 计算形状梯度
            auto shapeGrad = computeShapeGradient(sim, result, target);

            // 更新水平集
            shape.update(-learningRate * shapeGrad);
            shape.reinitialize();  // 保持符号距离性质
        }
    }
};

// 原始代码
float f(float x, float y) {
    float a = x * y;
    float b = sin(x);
    float c = a + b;
    return c;
}

// 计算图表示
struct Node {
    enum Op { INPUT, MUL, SIN, ADD };
    Op operation;
    std::vector<Node*> inputs;
    float value;
    float gradient;
};

template<typename T>
struct Dual {
    T value;
    T derivative;

    Dual(T v, T d = 0) : value(v), derivative(d) {}

    Dual operator+(const Dual& other) const {
        return Dual(value + other.value, 
                   derivative + other.derivative);
    }

    Dual operator*(const Dual& other) const {
        return Dual(value * other.value,
                   derivative * other.value + value * other.derivative);
    }

    friend Dual sin(const Dual& x) {
        return Dual(std::sin(x.value), 
                   std::cos(x.value) * x.derivative);
    }
};

// 使用示例
Dual<float> f(Dual<float> x, Dual<float> y) {
    return sin(x) + x * y;
}

// 计算 df/dx at (2, 3)
auto result = f(Dual<float>(2, 1), Dual<float>(3, 0));
float dfdx = result.derivative;  // cos(2) + 3