laplace_marginal_tol with integrate_1d inside likelihood doesn't compile

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));
  }
}