Stan issue with adj updates

clang.txt

--- Compiling, linking C++ code ---
clang++ -std=c++1y -Wno-unknown-warning-option -Wno-tautological-compare -Wno-sign-compare -D_REENTRANT -Wno-ignored-attributes      -I stan/lib/stan_math/lib/tbb_2020.3/include    -O3 -I src -I stan/src -I lib/rapidjson_1.1.0/ -I lib/CLI11-1.9.1/ -I stan/lib/stan_math/ -I stan/lib/stan_math/lib/eigen_3.3.9 -I stan/lib/stan_math/lib/boost_1.75.0 -I stan/lib/stan_math/lib/sundials_6.0.0/include -I stan/lib/stan_math/lib/sundials_6.0.0/src/sundials    -DBOOST_DISABLE_ASSERTS          -c -include-pch stan/src/stan/model/model_header.hpp.gch -include /Users/duke/Nextcloud/Work/HiresHMC/bigpbvepinc.hpp -x c++ -o /Users/duke/Nextcloud/Work/HiresHMC/bigpbvep.o /Users/duke/Nextcloud/Work/HiresHMC/bigpbvep.hpp
In file included from <built-in>:1:
user_header.hpp:82:14: error: no viable overloaded '+='
    xz.adj() += g_x;
    ~~~~~~~~ ^  ~~~
user_header.hpp:105:7: note: in instantiation of function template specialization 'bigpbvep_model_namespace::ode_rhs_rev<stan::math::arena_matrix<Eigen::Matrix<stan::math::var_value<double>, -1, 1, 0>>, stan::math::arena_matrix<Eigen::Matrix<stan::math::var_value<double>, -1, 1, 0>>, stan::math::arena_matrix<Eigen::Matrix<stan::math::var_value<double>, -1, 1, 0>>, stan::math::arena_matrix<Eigen::Matrix<double, -1, 1, 0>>, stan::math::var_value<double>, std::vector<int, stan::math::arena_allocator<int>>, std::vector<int, stan::math::arena_allocator<int>>>' requested here
      ode_rhs_rev(time, dxz_a, xz_a, C_nnz, C_w_a, C_v_a, C_u_a, I1, tau0, K_a, eta_a);
      ^
user_header.hpp:415:9: note: in instantiation of function template specialization 'bigpbvep_model_namespace::ode_rhs_c<Eigen::Matrix<stan::math::var_value<double>, -1, 1, 0>, Eigen::Matrix<stan::math::var_value<double>, -1, 1, 0>, Eigen::Map<Eigen::Matrix<double, -1, 1, 0>, 0>, stan::math::var_value<double>, nullptr>' requested here
        ode_rhs_c(time, xz, C_nnz, C_w, C_v, C_u, I1, tau0, K,
        ^
user_header.hpp:506:11: note: in instantiation of function template specialization 'bigpbvep_model_namespace::ode_step<double, double, Eigen::Matrix<stan::math::var_value<double>, -1, 1, 0>, Eigen::Map<Eigen::Matrix<double, -1, 1, 0>, 0>, double, double, stan::math::var_value<double>, Eigen::Matrix<stan::math::var_value<double>, -1, 1, 0>, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr>' requested here
          ode_step(dt, (t * dt),
          ^
user_header.hpp:588:11: note: in instantiation of function template specialization 'bigpbvep_model_namespace::ode_sol<double, Eigen::Matrix<stan::math::var_value<double>, -1, 1, 0>, Eigen::Map<Eigen::Matrix<double, -1, 1, 0>, 0>, double, double, stan::math::var_value<double>, Eigen::Matrix<stan::math::var_value<double>, -1, 1, 0>, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr>' requested here
          ode_sol(dt, nt,
          ^
user_header.hpp:772:10: note: in instantiation of function template specialization 'bigpbvep_model_namespace::partial<stan::math::var_value<double>, double, Eigen::Map<Eigen::Matrix<double, -1, 1, 0>, 0>, double, double, stan::math::var_value<double>, Eigen::Matrix<stan::math::var_value<double>, -1, 1, 0>, stan::math::var_value<double>, Eigen::Map<Eigen::Matrix<double, -1, -1, 0>, 0>, Eigen::Matrix<stan::math::var_value<double>, -1, 1, 0>, Eigen::Map<Eigen::Matrix<double, -1, -1, 0>, 0>, stan::math::var_value<double>, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr>' requested here
  return partial(idx, i0 + 1, i1 + 1, xzs, dt, nt, C_nnz, C_w, C_v, C_u, I1,
         ^
/Users/duke/.cmdstan/cmdstan-2.29.2/stan/lib/stan_math/stan/math/prim/functor/reduce_sum.hpp:215:10: note: in instantiation of function template specialization 'bigpbvep_model_namespace::partial_rsfunctor__::operator()<stan::math::var_value<double>, double, Eigen::Map<Eigen::Matrix<double, -1, 1, 0>, 0>, double, double, stan::math::var_value<double>, Eigen::Matrix<stan::math::var_value<double>, -1, 1, 0>, stan::math::var_value<double>, Eigen::Map<Eigen::Matrix<double, -1, -1, 0>, 0>, Eigen::Matrix<stan::math::var_value<double>, -1, 1, 0>, Eigen::Map<Eigen::Matrix<double, -1, -1, 0>, 0>, stan::math::var_value<double>, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr>' requested here
  return ReduceFunction()(std::forward<Vec>(vmapped), 0, vmapped.size() - 1,
         ^
/Users/duke/Nextcloud/Work/HiresHMC/bigpbvep.hpp:1725:14: note: in instantiation of function template specialization 'bigpbvep_model_namespace::bigpbvep_model::log_prob_impl<false, false, Eigen::Matrix<stan::math::var_value<double>, -1, 1, 0>, Eigen::Matrix<int, -1, 1, 0>, nullptr, nullptr>' requested here
      return log_prob_impl<propto__, jacobian__>(params_r, params_i, pstream);
             ^
/Users/duke/.cmdstan/cmdstan-2.29.2/stan/src/stan/model/model_base_crtp.hpp:95:50: note: in instantiation of function template specialization 'bigpbvep_model_namespace::bigpbvep_model::log_prob<false, false, stan::math::var_value<double>>' requested here
    return static_cast<const M*>(this)->template log_prob<false, false>(theta,
                                                 ^
/Users/duke/Nextcloud/Work/HiresHMC/bigpbvep.hpp:807:3: note: in instantiation of member function 'stan::model::model_base_crtp<bigpbvep_model_namespace::bigpbvep_model>::log_prob' requested here
  ~bigpbvep_model() { }
  ^
/Users/duke/.cmdstan/cmdstan-2.29.2/stan/lib/stan_math/lib/eigen_3.3.9/Eigen/src/Core/DenseBase.h:289:14: note: candidate function not viable: no known conversion from 'const CwiseUnaryOp<Eigen::MatrixBase<Eigen::Map<Eigen::Matrix<stan::math::var_value<double>, -1, 1, 0>, 0>>::adj_Op, const Eigen::Map<Eigen::Matrix<stan::math::var_value<double>, -1, 1, 0>, 0>>' to 'Eigen::MatrixBase<Eigen::CwiseUnaryOp<Eigen::MatrixBase<Eigen::Map<Eigen::Matrix<stan::math::var_value<double>, -1, 1, 0>, 0>>::adj_Op, const Eigen::Map<Eigen::Matrix<stan::math::var_value<double>, -1, 1, 0>, 0>>>' for object argument
    Derived& operator+=(const EigenBase<OtherDerived> &other);
             ^
/Users/duke/.cmdstan/cmdstan-2.29.2/stan/lib/stan_math/lib/eigen_3.3.9/Eigen/src/Core/MatrixBase.h:158:14: note: candidate function template not viable: no known conversion from 'const CwiseUnaryOp<Eigen::MatrixBase<Eigen::Map<Eigen::Matrix<stan::math::var_value<double>, -1, 1, 0>, 0>>::adj_Op, const Eigen::Map<Eigen::Matrix<stan::math::var_value<double>, -1, 1, 0>, 0>>' to 'Eigen::MatrixBase<Eigen::CwiseUnaryOp<Eigen::MatrixBase<Eigen::Map<Eigen::Matrix<stan::math::var_value<double>, -1, 1, 0>, 0>>::adj_Op, const Eigen::Map<Eigen::Matrix<stan::math::var_value<double>, -1, 1, 0>, 0>>>' for object argument
    Derived& operator+=(const MatrixBase<OtherDerived>& other);
             ^
/Users/duke/.cmdstan/cmdstan-2.29.2/stan/lib/stan_math/lib/eigen_3.3.9/Eigen/src/Core/MatrixBase.h:476:46: note: candidate template ignored: could not match 'ArrayBase' against 'Matrix'
    template<typename OtherDerived> Derived& operator+=(const ArrayBase<OtherDerived>& )

functions.stan

functions {

  vector ode_rhs_c(real time, vector xz,
                 int C_nnz, vector C_w, array[] int C_v, array[] int C_u,
                 real I1, real tau0, real K, vector eta);
    
  vector ode_rhs(real time, vector xz,
                 int C_nnz, vector C_w, array[] int C_v, array[] int C_u,
                 real I1, real tau0, real K, vector eta) {
    int nn = rows(eta);
    vector[nn] x = xz[1:nn];
    vector[nn] z = xz[nn+1:2*nn];

    vector[nn] gx = csr_matrix_times_vector(nn, nn, C_w, C_v, C_u, x);
    vector[nn] dx = 1.0 - x.*x.*x - 2.0*x.*x - z + I1;
    vector[nn] dz = (1/tau0)*(4*(x - eta) - z - K*gx);
    return append_row(dx, dz);
  }
}

user_header.hpp

#include "stan/math/prim/meta/is_rev_matrix.hpp"
#include "stan/math/prim/meta/is_var_matrix.hpp"
#include <stan/math.hpp>
#include <stan/math/rev/core/reverse_pass_callback.hpp>

namespace model_namespace {

  // namespaces
  using namespace stan;
  using namespace stan::math;
  using namespace Eigen;
  using namespace std;
  
  // fwd pass
  template <typename V2, typename V3, typename V4, typename S1,
            require_not_var_t<S1>* = nullptr
            >
  VectorXd ode_rhs_c(const double time, const V2& xz,
                const int& C_nnz, const V4& C_w, const std::vector<int>& C_v, const std::vector<int>& C_u,
                const double I1, const double tau0, const S1& K, const V3& eta, std::ostream* pstream) {
    int nn = xz.size() / 2;
    auto x = xz.head(nn).array(), z = xz.tail(nn).array();
    // auto gx = csr_matrix_times_vector(nn, nn, C_w, C_v, C_u, x);
    VectorXd gx(nn);
    for (int i=0; i<C_u.size()-1; i++) {
      gx[i] = 0.0;
      for (int j=C_u[i]; j<C_u[i+1]; j++)
        gx[i] += C_w[j-1] * x[C_v[j - 1] - 1];
    }
    VectorXd dxz(xz.size());
    dxz.head(nn) = 1.0 - x.cube() - 2*x.square() - z + I1;
    dxz.tail(nn) = (1/tau0)*(4*(x - eta.array()) - z - (K*gx).array());
    return dxz;
  }

  // reusable rev
  template <typename V1, typename V2, typename V3, typename V4, typename S1, typename I3, typename I2>
  void ode_rhs_rev(const double time, V1& dxz, const V2& xz,
                const int& C_nnz, const V4& C_w, const I3 & C_v, const I2& C_u,
                const double I1, const double tau0, const S1& K, V3& eta)
  {
    double rtau0 = 1 / tau0;
    int nn = xz.size() / 2;
    VectorXd g = dxz.adj();
    VectorXd g_dx = g.head(nn), g_dz = g.tail(nn);
    VectorXd x = value_of(xz.head(nn));

    eta.adj() += -g_dz.tail(nn)*rtau0*4;

    VectorXd gx(nn), C_g_dz(nn);
    for (int i=0; i<C_u.size()-1; i++) {
      gx[i] = 0.0;
      C_g_dz[i] = 0.0;
      for (int j=C_u[i]; j<C_u[i+1]; j++) {
        gx[i] += C_w[j-1] * x[C_v[j - 1] - 1];
        C_g_dz[i] += C_w[j-1] * g_dz[C_v[j - 1] - 1];
      }
    }
    double ksum = 0.0;
    for (int i=0;i<nn;i++) ksum += g_dz(i)*gx(i);
    K.adj() -= ksum * rtau0;
    VectorXd g_x(xz.size());
    g_x.fill(0.0);
    g_x.head(nn).array() += g_dx.array()*(-3*x.array().square() - 4*x.array());
    g_x.head(nn) += g_dz*4*rtau0 - value_of(K)*rtau0*(C_g_dz);
    g_x.tail(nn) += -g_dx;
    g_x.tail(nn) += -g_dz*rtau0;
    // typename decltype(xz.adj())::foo bar;
    xz.adj() += g_x;
  }

  typedef Matrix<stan::math::var_value<double>, Dynamic, 1> VarVec;

  // rev
  template <typename V2, typename V3, typename V4, typename S1,
            require_var_t<S1>* = nullptr
            >
  VarVec ode_rhs_c(const double time, const V2& xz,
                const int& C_nnz, const V4& C_w, const std::vector<int>& C_v, const std::vector<int>& C_u,
                const double I1, const double tau0, const S1& K, const V3& eta, std::ostream* pstream) {
    VectorXd dxz_val = ode_rhs_c(
      time, value_of(xz), C_nnz, C_w, C_v, C_u, I1, tau0, value_of(K), value_of(eta), pstream);
    VarVec dxz(dxz_val);
    arena_t<VarVec> dxz_a(dxz);
    arena_t<V2> xz_a(xz);
    arena_t<V4> C_w_a(C_w);
    arena_t<const std::vector<int>> C_v_a(C_v.begin(), C_v.end());
    arena_t<const std::vector<int>> C_u_a(C_u.begin(), C_u.end());
    arena_t<S1> K_a(K);
    arena_t<V3> eta_a(eta);
    reverse_pass_callback([time, dxz_a, xz_a, C_nnz, C_w_a, C_v_a, C_u_a, I1, tau0, K_a, eta_a]() mutable {
      ode_rhs_rev(time, dxz_a, xz_a, C_nnz, C_w_a, C_v_a, C_u_a, I1, tau0, K_a, eta_a);
    });
    return dxz;
  }
}