如何基于Boost Autodiff实现双参数函数的导数计算？

我尝试用Boost Autodiff API为双参数函数double response(double x, double y)计算导数。该函数返回由参数索引的矩阵单元格值，因此无法通过常规方式实现自动微分。

根据Boost Autodiff手册第5.3节的说明，这类场景的解决方案是手动实现导数。我可以预先计算存储导数值的第二个矩阵，这完全可行；但不清楚如何基于Autodiff API编写适配后的微分函数。

手册中提供了一个单变量自定义导数函数的示例（复制如下）：

template <typename RealType, size_t Order>
fvar<RealType, Order> cos(fvar<RealType, Order> const& cr) {
  using std::cos;
  using std::sin;
  using root_type = typename fvar<RealType, Order>::root_type;
  constexpr size_t order = fvar<RealType, Order>::order_sum;
  root_type const d0 = cos(static_cast<root_type>(cr));
  if constexpr (order == 0)
    return fvar<RealType, Order>(d0);
  else {
    root_type const d1 = -sin(static_cast<root_type>(cr));
    root_type const derivatives[4]{d0, d1, -d0, -d1};
    return cr.apply_derivatives(order, [&derivatives](size_t i) { return derivatives[i & 3]; });
  }
}

我需要将这个示例推广到double response(double x, double y)函数，支持对x、y其中一个或两者同时求导，但不知道具体实现方式。

内容的提问来源于stack exchange，提问作者xperroni