Embedded Laplace - model.hpp gets C++ compiler novel

Developers
1

I don’t know much about GPs, but I’ve always wanted to learn. I have a WHO dataset of PM2.5 air pollution data - the one that Gabry et al analyzed in the “Visualization in Bayesian Workflow” paper.

I wrote a simple regression model which gives reasonable results and tried to code the corresponding regression + GP for spatial smoothing. I gave Claude all the docs, the papers, and the simple model and amazingly enough, she wrote the expanded model which transpiled to C++ on the first go.

Here’s the model:

// Spatial Regression for WHO PM2.5 Data - North America
// Using Embedded Laplace Approximation with Matern 3/2 kernel

functions {
  // Likelihood function for embedded Laplace approximation
  // theta: spatial random effects (length N)
  real spatial_likelihood(vector theta, 
                         vector y_data, 
                         vector mu_fixed_data,
                         real sigma_obs_data) {
    vector[rows(y_data)] mu_total = mu_fixed_data + theta;
    return normal_lpdf(y_data | mu_total, sigma_obs_data);
  }
  
  // Covariance function for spatial process
  matrix spatial_covariance(array[] vector coords_data, 
                           real sigma_spatial, 
                           real length_scale) {
    return gp_matern32_cov(coords_data, sigma_spatial, length_scale);
  }
}

data {
  int<lower=1> N;                    // Number of observations
  vector[N] y;                       // log(PM2.5) observations
  matrix[N, 2] coords;               // Spatial coordinates (km)
  
  array[N] int<lower=0,upper=1> x_urban;      // Urban monitor indicator
  array[N] int<lower=0,upper=1> x_geocoded;   // High-quality geocoding
  array[N] int<lower=1,upper=3> x_country;    // Country, 1 = CAN, 2 = MEX, 3 = USA
}

transformed data {
  // Convert coordinates matrix to array of vectors for GP functions
  array[N] vector[2] coords_vec;
  for (i in 1:N) {
    coords_vec[i] = coords[i]';
  }
}

parameters {
  // Fixed effects (mean structure) - same as original model
  real beta_0;                       // Intercept
  real beta_urban;                   // Urban effect
  real beta_geocoded;                // Geocoding quality effect
  vector[3] beta_country;            // Country effect
  
  // Spatial process parameters
  real<lower=0> sigma_spatial;       // Spatial variance
  real<lower=0> length_scale;        // Spatial correlation length scale (km)
  
  // Observation noise
  real<lower=0> sigma_obs;           // Measurement error standard deviation
}

transformed parameters {
  // Mean function (linear predictor) - same structure as original
  vector[N] mu_fixed = beta_0 + 
    beta_urban * to_vector(x_urban) +
    beta_geocoded * to_vector(x_geocoded) +
    beta_country[x_country];
}

model {
  // Priors - same as original model
  beta_0 ~ normal(2.5, 1);           // log(12) ≈ 2.5, reasonable PM2.5 baseline
  beta_urban ~ normal(0.3, 0.2);     // Urban sites typically higher
  beta_geocoded ~ normal(0, 0.1);    // Geocoding quality (small effect)
  beta_country ~ normal(0, 0.5);
  
  // Priors for spatial parameters
  sigma_spatial ~ normal(0, 1);      // Spatial standard deviation
  length_scale ~ inv_gamma(2, 50);   // Length scale ~50km with some uncertainty
  
  sigma_obs ~ normal(0, 0.5);        // Measurement noise
  
  // Initial guess for spatial random effects (zeros)
  vector[N] theta_init = rep_vector(0, N);
  
  // Embedded Laplace approximation for spatial likelihood
  target += laplace_marginal(spatial_likelihood, 
                           (y, mu_fixed, sigma_obs),
                           theta_init,
                           spatial_covariance,
                           (coords_vec, sigma_spatial, length_scale));
}

generated quantities {
  // Sample spatial random effects from approximate posterior
  vector[N] theta_init = rep_vector(0, N);
  vector[N] spatial_effects = laplace_latent_rng(spatial_likelihood,
                                               (y, mu_fixed, sigma_obs),
                                               theta_init,
                                               spatial_covariance,
                                               (coords_vec, sigma_spatial, length_scale));
  
  // Total mean including spatial effects
  vector[N] mu = mu_fixed + spatial_effects;
  
  // Posterior predictive samples
  array[N] real y_rep = normal_rng(mu, sigma_obs);
}
2

here’s the abridged error novel:

22:52:30 - cmdstanpy - INFO - compiling stan file /Users/mitzi/github/zmorris/embedded-laplace/stan/gp_pm25.stan to exe file /Users/mitzi/github/zmorris/embedded-laplace/stan/gp_pm25
--- Translating Stan model to C++ code ---
bin/stanc --filename-in-msg=gp_pm25.stan --o=/Users/mitzi/github/zmorris/embedded-laplace/stan/gp_pm25.hpp /Users/mitzi/github/zmorris/embedded-laplace/stan/gp_pm25.stan

--- Compiling C++ code ---
clang++ -std=c++17 -Wno-unknown-warning-option -Wno-tautological-compare -Wno-sign-compare -D_REENTRANT -ffp-contract=off -Wno-ignored-attributes      -I stan/lib/stan_math/lib/tbb_2020.3/include    -O3 -I src -I stan/src -I stan/lib/rapidjson_1.1.0/ -I lib/CLI11-1.9.1/ -I stan/lib/stan_math/ -I stan/lib/stan_math/lib/eigen_3.4.0 -I stan/lib/stan_math/lib/boost_1.87.0 -I stan/lib/stan_math/lib/sundials_6.1.1/include -I stan/lib/stan_math/lib/sundials_6.1.1/src/sundials    -DBOOST_DISABLE_ASSERTS          -c -include-pch stan/src/stan/model/model_header.hpp.gch/model_header_16_0.hpp.gch -x c++ -o /Users/mitzi/github/zmorris/embedded-laplace/stan/gp_pm25.o /Users/mitzi/github/zmorris/embedded-laplace/stan/gp_pm25.hpp
In file included from /Users/mitzi/github/zmorris/embedded-laplace/stan/gp_pm25.hpp:1:
In file included from /Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/src/stan/model/model_header.hpp:4:
In file included from /Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/src/stan/model/model_base.hpp:8:
In file included from /Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/stan/math/rev/core.hpp:4:
In file included from /Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/stan/math/rev/core/accumulate_adjoints.hpp:4:
In file included from /Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/stan/math/prim/meta.hpp:72:
In file included from /Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/stan/math/prim/meta/append_return_type.hpp:4:
In file included from /Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/stan/math/prim/fun/Eigen.hpp:31:
In file included from /Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/lib/eigen_3.4.0/Eigen/Dense:1:
In file included from /Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/lib/eigen_3.4.0/Eigen/Core:279:
/Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/lib/eigen_3.4.0/Eigen/src/Core/AssignEvaluator.h:886:3: error: static assertion failed due to requirement 'Eigen::internal::is_lvalue<Eigen::CwiseUnaryOp<(lambda at /Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/stan/math/prim/fun/value_of.hpp:85:28), const Eigen::Matrix<stan::math::var_value<double, void>, -1, 1, 0, -1, 1>>>::value': THIS_EXPRESSION_IS_NOT_A_LVALUE__IT_IS_READ_ONLY
  886 |   EIGEN_STATIC_ASSERT_LVALUE(Dst)
      |   ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/lib/eigen_3.4.0/Eigen/src/Core/util/StaticAssert.h:203:27: note: expanded from macro 'EIGEN_STATIC_ASSERT_LVALUE'
  203 |       EIGEN_STATIC_ASSERT(Eigen::internal::is_lvalue<Derived>::value, \
      |       ~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  204 |                           THIS_EXPRESSION_IS_NOT_A_LVALUE__IT_IS_READ_ONLY)
      |                           ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/lib/eigen_3.4.0/Eigen/src/Core/util/StaticAssert.h:33:54: note: expanded from macro 'EIGEN_STATIC_ASSERT'
   33 |     #define EIGEN_STATIC_ASSERT(X,MSG) static_assert(X,#MSG);
      |                                                      ^
/Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/lib/eigen_3.4.0/Eigen/src/Core/NoAlias.h:43:7: note: in instantiation of function template specialization 'Eigen::internal::call_assignment_no_alias<Eigen::CwiseUnaryOp<(lambda at /Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/stan/math/prim/fun/value_of.hpp:85:28), const Eigen::Matrix<stan::math::var_value<double>, -1, 1>>, Eigen::Product<Eigen::Matrix<double, -1, -1>, Eigen::Matrix<double, -1, 1>>, Eigen::internal::assign_op<double, double>>' requested here
   43 |       call_assignment_no_alias(m_expression, other.derived(), internal::assign_op<Scalar,typename OtherDerived::Scalar>());
      |       ^
stan/lib/stan_math/stan/math/mix/functor/laplace_marginal_density.hpp:535:25: note: in instantiation of function template specialization 'Eigen::NoAlias<Eigen::CwiseUnaryOp<(lambda at /Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/stan/math/prim/fun/value_of.hpp:85:28), const Eigen::Matrix<stan::math::var_value<double>, -1, 1>>, Eigen::MatrixBase>::operator=<Eigen::Product<Eigen::Matrix<double, -1, -1>, Eigen::Matrix<double, -1, 1>>>' requested here
  535 |         theta.noalias() = covariance * a;
      |                         ^
/Users/mitzi/github/zmorris/embedded-laplace/stan/gp_pm25.hpp:986:12: note: in instantiation of function template specialization 'gp_pm25_model_namespace::gp_pm25_model::log_prob_impl<false, false, Eigen::Matrix<stan::math::var_value<double>, -1, 1>, Eigen::Matrix<int, -1, 1>, nullptr, nullptr, nullptr>' requested here
  986 |     return log_prob_impl<propto__, jacobian__>(params_r, params_i, pstream);
      |            ^
/Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/src/stan/model/model_base_crtp.hpp:98:50: note: in instantiation of function template specialization 'gp_pm25_model_namespace::gp_pm25_model::log_prob<false, false, stan::math::var_value<double>>' requested here
   98 |     return static_cast<const M*>(this)->template log_prob<false, false>(theta,
      |                                                  ^
/Users/mitzi/github/zmorris/embedded-laplace/stan/gp_pm25.hpp:219:3: note: in instantiation of member function 'stan::model::model_base_crtp<gp_pm25_model_namespace::gp_pm25_model>::log_prob' requested here
  219 |   ~gp_pm25_model() {}
      |   ^
In file included from /Users/mitzi/github/zmorris/embedded-laplace/stan/gp_pm25.hpp:1:
In file included from /Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/src/stan/model/model_header.hpp:4:
In file included from /Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/src/stan/model/model_base.hpp:8:
In file included from /Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/stan/math/rev/core.hpp:4:
In file included from /Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/stan/math/rev/core/accumulate_adjoints.hpp:4:
In file included from /Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/stan/math/prim/meta.hpp:72:
In file included from /Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/stan/math/prim/meta/append_return_type.hpp:4:
In file included from /Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/stan/math/prim/fun/Eigen.hpp:31:
In file included from /Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/lib/eigen_3.4.0/Eigen/Dense:1:
In file included from /Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/lib/eigen_3.4.0/Eigen/Core:301:
/Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/lib/eigen_3.4.0/Eigen/src/Core/CwiseNullaryOp.h:347:20: error: no viable overloaded '='
  347 |   return derived() = Constant(rows(), cols(), val);
      |          ~~~~~~~~~ ^ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/lib/eigen_3.4.0/Eigen/src/Core/CwiseNullaryOp.h:548:10: note: in instantiation of member function 'Eigen::DenseBase<Eigen::CwiseUnaryOp<(lambda at /Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/stan/math/prim/fun/value_of.hpp:85:28), const Eigen::Matrix<stan::math::var_value<double>, -1, 1>>>::setConstant' requested here
  548 |   return setConstant(Scalar(0));
      |          ^
/Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/lib/eigen_3.4.0/Eigen/src/Core/ProductEvaluators.h:349:9: note: in instantiation of member function 'Eigen::DenseBase<Eigen::CwiseUnaryOp<(lambda at /Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/stan/math/prim/fun/value_of.hpp:85:28), const Eigen::Matrix<stan::math::var_value<double>, -1, 1>>>::setZero' requested here
  349 |   { dst.setZero(); scaleAndAddTo(dst, lhs, rhs, Scalar(1)); }
      |         ^
/Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/lib/eigen_3.4.0/Eigen/src/Core/ProductEvaluators.h:148:37: note: in instantiation of function template specialization 'Eigen::internal::generic_product_impl_base<Eigen::Matrix<double, -1, -1>, Eigen::Matrix<double, -1, 1>, Eigen::internal::generic_product_impl<Eigen::Matrix<double, -1, -1>, Eigen::Matrix<double, -1, 1>, Eigen::DenseShape, Eigen::DenseShape, 7>>::evalTo<Eigen::CwiseUnaryOp<(lambda at /Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/stan/math/prim/fun/value_of.hpp:85:28), const Eigen::Matrix<stan::math::var_value<double>, -1, 1>>>' requested here
  148 |     generic_product_impl<Lhs, Rhs>::evalTo(dst, src.lhs(), src.rhs());
      |                                     ^
/Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/lib/eigen_3.4.0/Eigen/src/Core/AssignEvaluator.h:890:46: note: in instantiation of member function 'Eigen::internal::Assignment<Eigen::CwiseUnaryOp<(lambda at /Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/stan/math/prim/fun/value_of.hpp:85:28), const Eigen::Matrix<stan::math::var_value<double>, -1, 1>>, Eigen::Product<Eigen::Matrix<double, -1, -1>, Eigen::Matrix<double, -1, 1>>, Eigen::internal::assign_op<double, double>>::run' requested here
  890 |   Assignment<ActualDstTypeCleaned,Src,Func>::run(actualDst, src, func);
      |                                              ^
/Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/lib/eigen_3.4.0/Eigen/src/Core/NoAlias.h:43:7: note: in instantiation of function template specialization 'Eigen::internal::call_assignment_no_alias<Eigen::CwiseUnaryOp<(lambda at /Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/stan/math/prim/fun/value_of.hpp:85:28), const Eigen::Matrix<stan::math::var_value<double>, -1, 1>>, Eigen::Product<Eigen::Matrix<double, -1, -1>, Eigen::Matrix<double, -1, 1>>, Eigen::internal::assign_op<double, double>>' requested here
   43 |       call_assignment_no_alias(m_expression, other.derived(), internal::assign_op<Scalar,typename OtherDerived::Scalar>());
      |       ^
stan/lib/stan_math/stan/math/mix/functor/laplace_marginal_density.hpp:535:25: note: in instantiation of function template specialization 'Eigen::NoAlias<Eigen::CwiseUnaryOp<(lambda at /Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/stan/math/prim/fun/value_of.hpp:85:28), const Eigen::Matrix<stan::math::var_value<double>, -1, 1>>, Eigen::MatrixBase>::operator=<Eigen::Product<Eigen::Matrix<double, -1, -1>, Eigen::Matrix<double, -1, 1>>>' requested here
  535 |         theta.noalias() = covariance * a;
      |                         ^
/Users/mitzi/github/zmorris/embedded-laplace/stan/gp_pm25.hpp:986:12: note: in instantiation of function template specialization 'gp_pm25_model_namespace::gp_pm25_model::log_prob_impl<false, false, Eigen::Matrix<stan::math::var_value<double>, -1, 1>, Eigen::Matrix<int, -1, 1>, nullptr, nullptr, nullptr>' requested here
  986 |     return log_prob_impl<propto__, jacobian__>(params_r, params_i, pstream);
      |            ^
/Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/src/stan/model/model_base_crtp.hpp:98:50: note: in instantiation of function template specialization 'gp_pm25_model_namespace::gp_pm25_model::log_prob<false, false, stan::math::var_value<double>>' requested here
   98 |     return static_cast<const M*>(this)->template log_prob<false, false>(theta,
      |                                                  ^
/Users/mitzi/github/zmorris/embedded-laplace/stan/gp_pm25.hpp:219:3: note: in instantiation of member function 'stan::model::model_base_crtp<gp_pm25_model_namespace::gp_pm25_model>::log_prob' requested here
  219 |   ~gp_pm25_model() {}
      |   ^
/Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/lib/eigen_3.4.0/Eigen/src/Core/util/ForwardDeclarations.h:89:65: note: candidate function (the implicit copy assignment operator) not viable: no known conversion from 'const ConstantReturnType' (aka 'const CwiseNullaryOp<scalar_constant_op<double>, Eigen::Matrix<double, -1, 1, 0, -1, 1>>') to 'const CwiseUnaryOp<(lambda at /Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/stan/math/prim/fun/value_of.hpp:85:28), const Matrix<var_value<double>, -1, 1>>' for 1st argument
   89 | template<typename UnaryOp,   typename MatrixType>         class CwiseUnaryOp;
      |                                                                 ^~~~~~~~~~~~

.............   on and on ................


stan/lib/stan_math/stan/math/mix/functor/laplace_marginal_density.hpp:283:15: note: in instantiation of function template specialization 'Eigen::DenseBase<Eigen::CwiseUnaryOp<(lambda at /Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/stan/math/prim/fun/value_of.hpp:85:28), const Eigen::Matrix<stan::math::var_value<double>, -1, 1>>>::swap<Eigen::Matrix<double, -1, 1>>' requested here
  283 |         theta.swap(theta_tmp);
      |               ^
stan/lib/stan_math/stan/math/mix/functor/laplace_marginal_density.hpp:546:11: note: in instantiation of function template specialization 'stan::math::internal::line_search<Eigen::Matrix<double, -1, 1>, Eigen::Matrix<double, -1, 1>, Eigen::CwiseUnaryOp<(lambda at /Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/stan/math/prim/fun/value_of.hpp:85:28), const Eigen::Matrix<stan::math::var_value<double>, -1, 1>>, const gp_pm25_model_namespace::spatial_likelihood_functor__ &, std::tuple<Eigen::Matrix<double, -1, 1> &, Eigen::CwiseUnaryOp<(lambda at /Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/lib/stan_math/stan/math/prim/fun/value_of.hpp:85:28), const Eigen::Matrix<stan::math::var_value<double>, -1, 1>>, const double &> &, Eigen::Matrix<double, -1, -1> &, std::ostream>' requested here
  546 |           line_search(objective_new, a, theta, a_prev, ll_fun, ll_args_vals,
      |           ^
/Users/mitzi/github/zmorris/embedded-laplace/stan/gp_pm25.hpp:986:12: note: in instantiation of function template specialization 'gp_pm25_model_namespace::gp_pm25_model::log_prob_impl<false, false, Eigen::Matrix<stan::math::var_value<double>, -1, 1>, Eigen::Matrix<int, -1, 1>, nullptr, nullptr, nullptr>' requested here
  986 |     return log_prob_impl<propto__, jacobian__>(params_r, params_i, pstream);
      |            ^
/Users/mitzi/.cmdstan/cmdstan-2.37.0-rc1/stan/src/stan/model/model_base_crtp.hpp:98:50: note: in instantiation of function template specialization 'gp_pm25_model_namespace::gp_pm25_model::log_prob<false, false, stan::math::var_value<double>>' requested here
   98 |     return static_cast<const M*>(this)->template log_prob<false, false>(theta,
      |                                                  ^
/Users/mitzi/github/zmorris/embedded-laplace/stan/gp_pm25.hpp:219:3: note: in instantiation of member function 'stan::model::model_base_crtp<gp_pm25_model_namespace::gp_pm25_model>::log_prob' requested here
  219 |   ~gp_pm25_model() {}
      |   ^
8 errors generated.
make: *** [/Users/mitzi/github/zmorris/embedded-laplace/stan/gp_pm25.o] Error 1
rm /Users/mitzi/github/zmorris/embedded-laplace/stan/gp_pm25.hpp

Command ['make', 'STANCFLAGS+=--filename-in-msg=gp_pm25.stan', '/Users/mitzi/github/zmorris/embedded-laplace/stan/gp_pm25']
	error during processing No such file or directory
3

Not related to the error, but you don’t need Laplace if both prior and observation model are normal

4

This specific error was one of the things we tried to address in rc2, does it still happen if you update?

5

thanks.
plenty more datasets to work with.

6

rc2 has indeed fixed the error - model compiles.

unfortunately, it hangs - on my machine, it writes the header - almost all of the header, but stops mid-parameter name. I sent model and data to @stevebronder who will investigate.

in the meantime, considering other candidate datasets and models.