laplace_marginal_tol with integrate_1d inside likelihood doesn't compile

[avehtari]

Not a minimal example, but maybe the reason is obvious anyway? If laplace_marginal_tol likelihood function has integrate_1d inside, the C++ compilation fails. I'm fine if that is not supported, but then it would be good to check that before C++ compilation and report that it's not supported

Errror is quite long, so here only the end

/tmp/RtmpGPlDrK/model-12448b780d70f3.hpp:1504:43:   required from here
stan/lib/stan_math/stan/
math/rev/functor/integrate_1d.hpp:128:29: error: no match for call to ‘(const integrate_1d_adapter<disc_golf_base_nonconst_angle2_dist2_hier_error_la_model_namespace::integrand_functor__>) (const double&, const double&, std::ostream* const&, std::vector<stan::math::var_value<double>, std::allocator<stan::math::var_value<double> > >&, std::vector<stan::math::var_value<double>, std::allocator<stan::math::var_value<double> > >&, std::vector<stan::math::var_value<double>, std::allocator<stan::math::var_value<double> > >&)’
  128 |                     return f(x, xc, msgs, local_args...);
      |                            ~^~~~~~~~~~~~~~~~~~~~~~~~~~~~
In file included from stan/lib/stan_math/stan/math/prim/functor/integrate_1d.hpp:7,
                 from stan/lib/stan_math/stan/math/prim/functor/hcubature.hpp:27,
                 from stan/lib/stan_math/stan/math/prim/functor.hpp:14,
                 from stan/lib/stan_math/stan/math/rev/fun/to_arena.hpp:7,
                 from stan/lib/stan_math/stan/math/rev/core/operator_division.hpp:15,
                 from stan/lib/stan_math/stan/math/rev/core/operator_divide_equal.hpp:5,
                 from stan/lib/stan_math/stan/math/rev/core.hpp:32,
                 from stan/lib/stan_math/stan/math/mix/meta.hpp:6,
                 from stan/lib/stan_math/stan/math/mix.hpp:8,
                 from /tmp/RtmpGPlDrK/model-12448b780d70f3.hpp:3:
stan/lib/stan_math/stan/math/prim/functor/integrate_1d_adapter.hpp:20:8: note: candidate: ‘template<class T_a, class T_b, class T_theta> auto integrate_1d_adapter<F>::operator()(const T_a&, const T_b&, std::ostream*, const std::vector<T_theta>&, const std::vector<double, std::allocator<double> >&, const std::vector<int>&) const [with T_a = T_a; T_b = T_b; T_theta = T_theta; F = disc_golf_base_nonconst_angle2_dist2_hier_error_la_model_namespace::integrand_functor__]’
   20 |   auto operator()(const T_a& x, const T_b& xc, std::ostream* msgs,
      |        ^~~~~~~~
stan/lib/stan_math/stan
/math/prim/functor/integrate_1d_adapter.hpp:20:8: note:   template argument deduction/substitution failed:
In file included from stan/lib/stan_math/stan/math/rev/functor.hpp:15,
                 from stan/lib/stan_math/stan/math/mix/functor/wolfe_line_search.hpp:10,
                 from stan/lib/stan_math/stan/math/mix/functor/laplace_marginal_density_estimator.hpp:5,
                 from stan/lib/stan_math/stan/math/mix/functor/laplace_marginal_density.hpp:5,
                 from stan/lib/stan_math/stan/math/mix/functor/laplace_base_rng.hpp:5,
                 from stan/lib/stan_math/stan/math/mix/functor.hpp:11,
                 from stan/lib/stan_math/stan/math/mix.hpp:10,
                 from /tmp/RtmpGPlDrK/model-12448b780d70f3.hpp:3:
stan/lib/stan_math/stan/math/rev/functor/integrate_1d.hpp:128:29: note:   cannot convert ‘local_args#1’ (type ‘std::vector<stan::math::var_value<double>, std::allocator<stan::math::var_value<double> > >’) to type ‘const std::vector<double, std::allocator<double> >&’
  128 |                     return f(x, xc, msgs, local_args...);
      |                            ~^~~~~~~~~~~~~~~~~~~~~~~~~~~~

make: *** [make/program:77: /tmp/RtmpGPlDrK/model-12448b780d70f3.o] Error 1

And the Stan model code causing that error

// Disc golf putting model with 2D angular uncertainty and distance uncertainty
// Hierarchical version with player-level effects for sigma_angle_base, sigma_angle_slope, sigma_distance
functions {
  real integrand(real z, real notused, array[] real theta,
               array[] real X_j, array[] int data_j) {
    real sigma_epsilon = theta[1];
    real p_j = theta[2];
    int n_j = data_j[1];
    int y_j = data_j[2];
    real p = (z<1) ? exp(gamma_lpdf(z | 1, 1/sigma_epsilon) + binomial_lpmf(y_j | n_j, p_j*(1-z))) : 0;
    return (is_inf(p) || is_nan(p)) ? 0 : p;
  }

  // Log-likelihood function for a single player (group)
  // theta = [eta_1, eta_2, eta_3] (player-level effects on log scale)
  // Data tuple contains player-specific observations
  real player_group_ll(vector theta,
                       int n_g, array[] int y_g, array[] int n_attempts_g,
                       vector x_g, vector threshold_angle_g,
                       real sigma_base, real threshold_distance,
                       real sigma_epsilon, real R_val) {
    real sigma_angle_base_p = exp(theta[1]);
    real sigma_angle_slope_p = exp(theta[2]);
    real sigma_distance_p = exp(theta[3]);
    real p_base = Phi(1.0 / sigma_base);
    real ll = 0;
    for (i in 1:n_g) {
      real sigma_angle = sqrt(square(sigma_angle_base_p) + square(sigma_angle_slope_p * x_g[i]));
      real p_angle = rayleigh_cdf(threshold_angle_g[i] | sigma_angle);
      real p_distance;
      if (x_g[i] <= threshold_distance) {
        p_distance = 1.0;
      } else {
        p_distance = normal_cdf(1.0 | 0.0, sigma_distance_p * (x_g[i] - threshold_distance));
      }
      real pr_i = p_base * p_angle * p_distance;
      // Integrate out epsilon for this observation
      ll += log(integrate_1d(integrand,
                             0,
                             positive_infinity(),
                             append_array({sigma_epsilon}, {pr_i}),
                             {0},
                             append_array({n_attempts_g[i]}, {y_g[i]}),
                             1e-9));
    }
    return ll;
  }

  // Covariance function for player-level effects
  matrix cov_group(vector tau_v, matrix L_Omega_v) {
    return multiply_lower_tri_self_transpose(diag_pre_multiply(tau_v, L_Omega_v));
  }
}
data {
  int N_obs;
  int N_players;
  vector[N_obs] x;              // distance in feet
  array[N_obs] int n;           // number of attempts
  array[N_obs] int y;           // number of successes
  array[N_obs] int player;      // player index
  real R;                       // target radius (half of basket diameter)
}
transformed data {
  // Maximum allowable angle at each distance
  vector[N_obs] threshold_angle = asin(R ./ x);
  int K = 3;  // number of player-level parameters (sigma_angle_base, sigma_angle_slope, sigma_distance)

  // Prepare player-specific data arrays for Laplace approximation
  array[N_players] int n_pg;  // number of observations per player
  array[N_players, N_obs] int y_pg;  // y values per player (padded)
  array[N_players, N_obs] int n_attempts_pg;  // n values per player (padded)
  array[N_players, N_obs] real x_pg_arr;  // x values per player (padded)
  array[N_players, N_obs] real threshold_angle_pg_arr;  // threshold_angle per player (padded)

  // initialize counts
  for (p in 1:N_players) {
    n_pg[p] = 0;
  }
  // populate player-specific arrays
  for (j in 1:N_obs) {
    int p = player[j];
    n_pg[p] += 1;
    y_pg[p, n_pg[p]] = y[j];
    n_attempts_pg[p, n_pg[p]] = n[j];
    x_pg_arr[p, n_pg[p]] = x[j];
    threshold_angle_pg_arr[p, n_pg[p]] = threshold_angle[j];
  }

  // Control parameters for Laplace approximation
  real la_tolerance = 1e-4;
  int la_max_num_steps = 100, la_hessian_block_size = 3, la_solver = 2, la_max_steps_line_search = 100;
}
parameters {
  // Population-level parameters (on log scale for positive constraint)
  vector[K] mu;                     // population means
  vector<lower=0>[K] tau;           // population SDs
  cholesky_factor_corr[K] L_Omega;  // Cholesky factor of correlation matrix
  
  // Non-centered player effects (standard normal)
  matrix[K, N_players] z;
  
  // Pooled parameters
  real<lower=0> sigma_base;         // base angular precision (radians)
  real<lower=0> threshold_distance; // distance threshold
  real<lower=0> sigma_epsilon;
  vector<lower=0, upper=1>[N_obs] epsilon;
}
transformed parameters {
  // Player-level parameters (on log scale)
  matrix[N_players, K] eta;
  
  // Transform from non-centered to centered parameterization
  // eta = mu + diag(tau) * L_Omega * z
  eta = (diag_pre_multiply(tau, L_Omega) * z)';
  for (p in 1:N_players) {
    eta[p] = eta[p] + mu';
  }
  
  // Player-level parameters on natural scale
  vector<lower=0>[N_players] sigma_angle_base;
  vector<lower=0>[N_players] sigma_angle_slope;
  vector<lower=0>[N_players] sigma_distance;
  for (p in 1:N_players) {
    sigma_angle_base[p] = exp(eta[p, 1]);
    sigma_angle_slope[p] = exp(eta[p, 2]);
    sigma_distance[p] = exp(eta[p, 3]);
  }
  
  // Observation-level probabilities
  vector[N_obs] pr;
  for (j in 1:N_obs) {
    int pl = player[j];
    // base probability
    real p_base = Phi(1.0 / sigma_base);
    // probability of correct angle (distance-dependent SD)
    real sigma_angle = sqrt(square(sigma_angle_base[pl]) + square(sigma_angle_slope[pl] * x[j]));
    real p_angle = rayleigh_cdf(threshold_angle[j] | sigma_angle);
    // probability of correct distance
    real p_distance;
    if (x[j] <= threshold_distance) {
      p_distance = 1.0;
    } else {
      p_distance = normal_cdf(1.0 | 0.0, sigma_distance[pl] * (x[j] - threshold_distance));
    }
    // combined probability assuming independence
    pr[j] = p_base * p_angle * p_distance;
  }
}
model {
  // priors on population parameters
  mu[1] ~ normal(-6, 2);            // prior for mean log_sigma_angle_base
  mu[2] ~ normal(-6, 2);            // prior for mean log_sigma_angle_slope
  mu[3] ~ normal(-6, 2);            // prior for mean log_sigma_distance
  tau[1] ~ normal(0, 0.5);            // population SDs
  tau[2] ~ normal(0, 0.5);            // population SDs
  tau[3] ~ normal(0, 0.5);            // population SDs
  L_Omega ~ lkj_corr_cholesky(3);  // LKJ prior on correlation
  
  // non-centered parameterization: z ~ normal(0, 1)
  to_vector(z) ~ std_normal();
  
  // priors on pooled parameters
  sigma_base ~ normal(0, 1);
  threshold_distance ~ normal(40, 10);

  epsilon ~ gamma(1, 1/sigma_epsilon);

  // data model
  y ~ binomial(n, pr .* (1 - epsilon));
}
generated quantities {
  vector[N_obs] p = pr .* (1 - gamma_rng(1, 1/sigma_epsilon));
  vector[N_obs] log_lik;
  vector[N_obs] log_liki;
  for (j in 1:N_obs) {
    log_lik[j] = binomial_lpmf(y[j] | n[j], pr[j] * (1 - epsilon[j]));
    log_liki[j] = log(integrate_1d(integrand,
                              0,
                              positive_infinity(), // this works better than 1!
                              append_array({sigma_epsilon}, {pr[j]}),
			      {0}, // not used, but an empty array not allowed
			      append_array({n[j]}, {y[j]}),
			      1e-9));
  }
  // Player-level log-likelihoods using Laplace approximation
  vector[N_players] log_lik_g;
  for (g in 1:N_players) {
    vector[N_obs] x_pg_vec = to_vector(x_pg_arr[g]);
    vector[N_obs] threshold_angle_pg_vec = to_vector(threshold_angle_pg_arr[g]);

    // Use sampled player effects as starting point for Newton solver
    vector[K] theta_init = eta[g]';

    log_lik_g[g] = laplace_marginal_tol(player_group_ll,
                                     (n_pg[g], y_pg[g], n_attempts_pg[g],
                                      x_pg_vec, threshold_angle_pg_vec,
                                      sigma_base, threshold_distance,
                                      sigma_epsilon, R),
                                     cov_group, (tau, L_Omega),
                                     (theta_init, la_tolerance, la_max_num_steps,
                                      la_hessian_block_size, la_solver, la_max_steps_line_search, 1));
  }

  vector[N_players] p33;
  for (j in 1:N_players) {
    real x33 = 33;
    real threshold_angle33 = asin(R ./ x33);
    int pl = j;
    // base probability
    real p_base = Phi(1.0 / sigma_base);
    // probability of correct angle (distance-dependent SD)
    real sigma_angle = sqrt(square(sigma_angle_base[pl]) + square(sigma_angle_slope[pl] * x33));
    real p_angle = rayleigh_cdf(threshold_angle33 | sigma_angle);
    // probability of correct distance
    real p_distance;
    if (x33 <= threshold_distance) {
      p_distance = 1.0;
    } else {
      p_distance = normal_cdf(1.0 | 0.0, sigma_distance[pl] * (x33 - threshold_distance));
    }
    // combined probability assuming independence
    p33[pl] = p_base * p_angle * p_distance .* (1 - gamma_rng(1, 1/sigma_epsilon));
  }
}

[WardBrian]

Here's the complete error: https://gist.github.com/WardBrian/6f8ae460ca6f00d5af51a15ab3c11e3a

It also mentions

stan/lib/stan_math/lib/eigen_3.4.0/Eigen/src/Core/util/XprHelper.h:851:96: error: static assertion failed: YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY
  851 |   EIGEN_STATIC_ASSERT((Eigen::internal::has_ReturnType<ScalarBinaryOpTraits<LHS, RHS,BINOP> >::value), \
      |                       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~

which seems like it could be the issue

[avehtari]

At least years ago nested integrate_1d did work, so maybe integrate_1d nested inside laplace_marginal could work, too?

[WardBrian]

I don't think there is any inherent reason it shouldn't, this just looks like a bug

[SteveBronder]

This looks like a bug. Ty will investigate. My guess is that we are missing a cast somewhere

[SteveBronder]

Here is the full error without any skips in template instantiations

https://gist.github.com/SteveBronder/11cb838a72ddc5e1c2e919e966ed625e

So this is sort of an Eigen bug and sort of an us bug (but not a laplace bug!). For higher order autodiff and integrate_1d we have to use a finite diff method. The failure comes from the function to aggregate the tangents where Eigen thinks it is doing an operation on two types with fvar<var> and var scalar types. But we are actually doing a multiply of the fvar's adjoint with another adjoint which are both var types

template <typename FuncTangent, typename InputArg,
          require_st_fvar<InputArg>* = nullptr>
inline auto aggregate_tangent(const FuncTangent& tangent, const InputArg& arg) {
  return sum(apply_scalar_binary([](auto&& x, auto&& y) { return x * y.d_; },
                                 tangent, arg));
}

in apply_scalar_binary we call x.binary_expr(y) and so all Eigen sees is that we are trying to operate on two types that have mixed scalar types which they do not support.

I can do a fix for this that just gets the adjoint/tangent generally and that should fix this. since then we can make this over var, var`

template <typename FuncTangent, typename InputArg,
          require_st_fvar<InputArg>* = nullptr>
inline auto aggregate_tangent(const FuncTangent& tangent, const InputArg& arg) {
  return sum(apply_scalar_binary([](auto&& x, auto&& y) { return x * y; },
                                 tangent, tangent_of(arg)));
}

[SteveBronder]

I have a full fix for this here

https://github.com/stan-dev/math/tree/fix/finite-diff-tangents

I ended up writing generic adjoint_of and tangent_of functions as I've been meaning to write those anyway. I need one more test to make sure that Aki's example now compiles then I'll open up a branch.

[avehtari]

I noticed the Laplace doc PR says

3. The operations in the function must support higher-order automatic
differentiation (AD). Most functions in Stan support higher-order AD.
The exceptions are functions with specialized calls for reverse-mode AD, and
these are higher-order functions (algebraic solvers, differential equation
solvers, and integrators), the marginalization function for hidden Markov
models (HMM) function, and the embedded Laplace approximation itself.

@SteveBronder, are you saying your fix will make all these work? If not, then it would be good to give informative error message if some of these are not supported and user still tries to use them (e.g. I didn't remember this)

[WardBrian]

@avehtari the compiler will give you a specific error for the ones that don’t work, e.g. if you try to put reduce_sum in there. Integrate 1d is specifically supposed to thanks to some previous work by Andrew Johnson