atomic vector 1rst and 2nd derivatives

[Hervé DAKPO]

Hello all, I am estimating a likelihood function that implies the use of the parabolic cylinder function (pcf). For the pcf I translated the fotran code by Shanjie Zhang, Jianming Jin into c++. Unfortunately it contains so many conditional statements that it is incompatible with TMB-AD. I manage to write the gradient and the hessian associated with the pcf. I tried to include that in the TMB though atomic vector but it seems that TMB discards the hessian. Is there a way to make TMB also consider the second order derivative. Here is a sketch of what I have done:

TMB_ATOMIC_VECTOR_FUNCTION(
  pbdv_atomic,
  1,
  {
    // Recover the order from tx.
    // With 2 inputs, tx.size() == 2*(order+1)
    size_t order = tx.size()/2 - 1;

    // --- Forward Mode: Order 0 (function value) ---
    ty[0] = pbdvGPT(asDouble(tx[0]), asDouble(tx[1]));

    // --- Forward Mode: Order 1 (first derivative) ---
    if (order >= 1) {
      double v_val = asDouble(tx[0]);
      double x_val = asDouble(tx[1]);
      // The first-order directional derivatives for v and x are stored in tx[2] and tx[3]
      double v1 = asDouble(tx[2]);
      double x1 = asDouble(tx[3]);

      double dfdv = dpbdv_dv(v_val, x_val);
      double dfdx = dpbdv_dx(v_val, x_val);
      ty[1] = dfdv * v1 + dfdx * x1;
    }

    // --- Forward Mode: Order 2 (second derivative) ---
    if (order >= 2) {
      double v_val = asDouble(tx[0]);
      double x_val = asDouble(tx[1]);
      double v1 = asDouble(tx[2]);
      double x1 = asDouble(tx[3]);

      double d2fdv2 = d2pbdv_dv2(v_val, x_val);
      double d2fdvdx = d2pbdv_dvdx(v_val, x_val);
      double d2fdx2 = d2pbdv_dx2(v_val, x_val);

      ty[2] = d2fdv2 * v1 * v1 + 2.0 * d2fdvdx * v1 * x1 + d2fdx2 * x1 * x1;
    }
  },
  {
    // --- Reverse Mode (Adjoints) ---
    double v_val = asDouble(tx[0]);
    double x_val = asDouble(tx[1]);
    double dfdv = dpbdv_dv(v_val, x_val);
    double dfdx = dpbdv_dx(v_val, x_val);

    px[0] = py[0] * dfdv;  // derivative with respect to v
    px[1] = py[0] * dfdx;  // derivative with respect to x
  }
)

template<class Type>
Type pbdv(Type v, Type x) {
  // Create an input vector with the function inputs.
  CppAD::vector<Type> tx(2);
  tx[0] = v;
  tx[1] = x;

  // When TMB/CppAD calls pbdv_atomic in forward mode, it will automatically
  // extend tx with the necessary first‐order (and higher) directions.
  CppAD::vector<Type> ty = pbdv_atomic(tx);
  return ty[0];
}