Stan Math functions should not return Eigen expressions

[avehtari]

Summary:

Using cmdstan-2.22.0.tar.gz, fixed stan-dev/cmdstan#804 so that make build works.
EDIT: the same problem with the latest git commit e9075d0
EDIT: by @wds15 suggestion transferred tgo stanc3 repo

Using cmdstanr (but this is clearly an issue in cmdstan or stanc3)

One model compiles and runs successfully, but other one one fails the compilation.

Description:

This Stan model compiles and runs with previous cmdstan and 2.22.0

data {
  int<lower=0> D;           // number of dimensions
}
parameters {
  vector[D] mu; 
}
model {
  mu ~ normal(0, 1);
}

This model compiles and runs with the previous cmdstan, but doesn't compile with 2.22.0

data {
  int<lower=0> n;				      // number of observations
  int<lower=0> d;             // number of predictors
  int<lower=0,upper=1> y[n];	// outputs
  matrix[n,d] x;				      // inputs
  real<lower=0> scale_icept;	// prior std for the intercept
  real<lower=0> scale_global;	// scale for the half-t prior for tau
  real<lower=1> nu_global;	  // degrees of freedom for the half-t priors for tau
  real<lower=1> nu_local;		  // degrees of freedom for the half-t priors for lambdas
                              // (nu_local = 1 corresponds to the horseshoe)
  real<lower=0> slab_scale;   // for the regularized horseshoe
  real<lower=0> slab_df;
}

parameters {
  real beta0;
  vector[d] z;                // for non-centered parameterization
  real <lower=0> tau;         // global shrinkage parameter
  vector <lower=0>[d] lambda; // local shrinkage parameter
  real<lower=0> caux;
}

transformed parameters {
  vector[d] beta;                     // regression coefficients
  {
    vector[d] lambda_tilde;   // 'truncated' local shrinkage parameter
    real c = slab_scale * sqrt(caux); // slab scale
    lambda_tilde = sqrt( c^2 * square(lambda) ./ (c^2 + tau^2*square(lambda)));
    beta = z .* lambda_tilde*tau;
  }
}

model {
  // half-t priors for lambdas and tau, and inverse-gamma for c^2
  z ~ std_normal();
  lambda ~ student_t(nu_local, 0, 1);
  tau ~ student_t(nu_global, 0, scale_global*2);
  caux ~ inv_gamma(0.5*slab_df, 0.5*slab_df);
  beta0 ~ normal(0, scale_icept);
  
  y ~ bernoulli_logit_glm(x, beta0, beta);
}

The start of the (very long) error message:

--- Compiling, linking C++ code ---
g++ -std=c++1y -pthread -D_REENTRANT -Wno-sign-compare      -I stan/lib/stan_math/lib/tbb_2019_U8/include -O3 -I src -I stan/src -I lib/rapidjson_1.1.0/ -I stan/lib/stan_math/ -I stan/lib/stan_math/lib/eigen_3.3.3 -I stan/lib/stan_math/lib/boost_1.72.0 -I stan/lib/stan_math/lib/sundials_4.1.0/include    -DBOOST_DISABLE_ASSERTS      -c  -x c++ -o /u/77/ave/unix/proj/ovarian/bernoulli_logit_glm_rhs.o /u/77/ave/unix/proj/ovarian/bernoulli_logit_glm_rhs.hpp
/u/77/ave/unix/proj/ovarian/bernoulli_logit_glm_rhs.hpp: In instantiation of ‘T__ bernoulli_logit_glm_rhs_model_namespace::bernoulli_logit_glm_rhs_model::log_prob(std::vector<T_l>&, std::vector<int>&, std::ostream*) const [with bool propto__ = false; bool jacobian__ = false; T__ = double; std::ostream = std::basic_ostream<char>]’:
stan/src/stan/model/model_base_crtp.hpp:147:29:   required from ‘double stan::model::model_base_crtp<M>::log_prob(std::vector<double>&, std::vector<int>&, std::ostream*) const [with M = bernoulli_logit_glm_rhs_model_namespace::bernoulli_logit_glm_rhs_model; std::ostream = std::basic_ostream<char>]’
/u/77/ave/unix/proj/ovarian/bernoulli_logit_glm_rhs.hpp:847:1:   required from here
/u/77/ave/unix/proj/ovarian/bernoulli_logit_glm_rhs.hpp:392:32: error: no matching function for call to ‘multiply(double, Eigen::CwiseUnaryOp<stan::math::apply_scalar_unary<F, T, typename std::enable_if<stan::is_eigen<S, void>::value, void>::type>::apply(const T&) [with F = stan::math::square_fun; T = Eigen::Matrix<double, -1, 1>]::<lambda(stan::math::apply_scalar_unary<stan::math::square_fun, Eigen::Matrix<double, -1, 1>, void>::scalar_t)>, const Eigen::Matrix<double, -1, 1> >)’
             elt_divide(multiply(pow(c, 2), square(lambda)),
                                ^
In file included from stan/lib/stan_math/stan/math/prim/fun.hpp:195:0,
                 from stan/lib/stan_math/stan/math/prim.hpp:10,
                 from stan/lib/stan_math/stan/math/rev.hpp:11,
                 from stan/lib/stan_math/stan/math.hpp:148,
                 from stan/src/stan/model/model_header.hpp:4,
                 from /u/77/ave/unix/proj/ovarian/bernoulli_logit_glm_rhs.hpp:3:
stan/lib/stan_math/stan/math/prim/fun/multiply.hpp:29:51: note: candidate: template<int R, int C, class T1, class T2, class> Eigen::Matrix<typename stan::return_type<T1, T2>::type, R, C> stan::math::multiply(const Eigen::Matrix<T1, R, C>&, T2)
 inline Eigen::Matrix<return_type_t<T1, T2>, R, C> multiply(
                                                   ^
stan/lib/stan_math/stan/math/prim/fun/multiply.hpp:29:51: note:   template argument deduction/substitution failed:
/u/77/ave/unix/proj/ovarian/bernoulli_logit_glm_rhs.hpp:392:32: note:   mismatched types ‘const Eigen::Matrix<T1, R, C>’ and ‘double’
             elt_divide(multiply(pow(c, 2), square(lambda)),
                                ^
In file included from stan/lib/stan_math/stan/math/prim/fun.hpp:195:0,
                 from stan/lib/stan_math/stan/math/prim.hpp:10,
                 from stan/lib/stan_math/stan/math/rev.hpp:11,
                 from stan/lib/stan_math/stan/math.hpp:148,
                 from stan/src/stan/model/model_header.hpp:4,
                 from /u/77/ave/unix/proj/ovarian/bernoulli_logit_glm_rhs.hpp:3:
stan/lib/stan_math/stan/math/prim/fun/multiply.hpp:48:51: note: candidate: template<int R, int C, class T1, class T2, class> Eigen::Matrix<typename stan::return_type<T1, T2>::type, R, C> stan::math::multiply(T1, const Eigen::Matrix<T2, R, C>&)
 inline Eigen::Matrix<return_type_t<T1, T2>, R, C> multiply(
                                                   ^
stan/lib/stan_math/stan/math/prim/fun/multiply.hpp:48:51: note:   template argument deduction/substitution failed:
/u/77/ave/unix/proj/ovarian/bernoulli_logit_glm_rhs.hpp:392:32: note:   ‘Eigen::CwiseUnaryOp<stan::math::apply_scalar_unary<F, T, typename std::enable_if<stan::is_eigen<S, void>::value, void>::type>::apply(const T&) [with F = stan::math::square_fun; T = Eigen::Matrix<double, -1, 1>]::<lambda(stan::math::apply_scalar_unary<stan::math::square_fun, Eigen::Matrix<double, -1, 1>, void>::scalar_t)>, const Eigen::Matrix<double, -1, 1> >’ is not derived from ‘const Eigen::Matrix<T2, R, C>’
             elt_divide(multiply(pow(c, 2), square(lambda)),
                                ^
In file included from stan/lib/stan_math/stan/math/prim/fun.hpp:195:0,
                 from stan/lib/stan_math/stan/math/prim.hpp:10,
                 from stan/lib/stan_math/stan/math/rev.hpp:11,
                 from stan/lib/stan_math/stan/math.hpp:148,
                 from stan/src/stan/model/model_header.hpp:4,
                 from /u/77/ave/unix/proj/ovarian/bernoulli_logit_glm_rhs.hpp:3:
stan/lib/stan_math/stan/math/prim/fun/multiply.hpp:73:53: note: candidate: template<int R1, int C1, int R2, int C2, class T1, class T2, class> Eigen::Matrix<typename stan::return_type<T1, T2>::type, R1, C2> stan::math::multiply(const Eigen::Matrix<T1, R1, C1>&, const Eigen::Matrix<T2, R2, C2>&)
 inline Eigen::Matrix<return_type_t<T1, T2>, R1, C2> multiply(
                                                     ^
stan/lib/stan_math/stan/math/prim/fun/multiply.hpp:73:53: note:   template argument deduction/substitution failed:
/u/77/ave/unix/proj/ovarian/bernoulli_logit_glm_rhs.hpp:392:32: note:   mismatched types ‘const Eigen::Matrix<T1, R1, C1>’ and ‘double’
             elt_divide(multiply(pow(c, 2), square(lambda)),

A smaller example with the same (similar) error

data {
  int<lower=0> d;             // number of predictors
}
parameters {
  vector <lower=0>[d] lambda; // local shrinkage parameter
}
transformed parameters {
  vector[d] lambda_tilde;   // 'truncated' local shrinkage parameter
  real c = 1.0; // slab scale
  lambda_tilde = sqrt( c^2 * square(lambda) ./ (c^2 + square(lambda)));
}
model {
}

With an error (just the beginning of a very long message)

--- Compiling, linking C++ code ---
g++ -std=c++1y -pthread -D_REENTRANT -Wno-sign-compare      -I stan/lib/stan_math/lib/tbb_2019_U8/include -O3 -I src -I stan/src -I lib/rapidjson_1.1.0/ -I stan/lib/stan_math/ -I stan/lib/stan_math/lib/eigen_3.3.3 -I stan/lib/stan_math/lib/boost_1.72.0 -I stan/lib/stan_math/lib/sundials_4.1.0/include    -DBOOST_DISABLE_ASSERTS      -c  -x c++ -o /u/77/ave/unix/proj/ovarian/test.o /u/77/ave/unix/proj/ovarian/test.hpp
/u/77/ave/unix/proj/ovarian/test.hpp: In instantiation of ‘T__ test_model_namespace::test_model::log_prob(std::vector<T_l>&, std::vector<int>&, std::ostream*) const [with bool propto__ = false; bool jacobian__ = false; T__ = double; std::ostream = std::basic_ostream<char>]’:
stan/src/stan/model/model_base_crtp.hpp:147:29:   required from ‘double stan::model::model_base_crtp<M>::log_prob(std::vector<double>&, std::vector<int>&, std::ostream*) const [with M = test_model_namespace::test_model; std::ostream = std::basic_ostream<char>]’
/u/77/ave/unix/proj/ovarian/test.hpp:467:1:   required from here
/u/77/ave/unix/proj/ovarian/test.hpp:201:30: error: no matching function for call to ‘multiply(double, Eigen::CwiseUnaryOp<stan::math::apply_scalar_unary<F, T, typename std::enable_if<stan::is_eigen<S, void>::value, void>::type>::apply(const T&) [with F = stan::math::square_fun; T = Eigen::Matrix<double, -1, 1>]::<lambda(stan::math::apply_scalar_unary<stan::math::square_fun, Eigen::Matrix<double, -1, 1>, void>::scalar_t)>, const Eigen::Matrix<double, -1, 1> >)’
           elt_divide(multiply(pow(c, 2), square(lambda)),
                              ^
In file included from stan/lib/stan_math/stan/math/prim/fun.hpp:195:0,
                 from stan/lib/stan_math/stan/math/prim.hpp:10,
                 from stan/lib/stan_math/stan/math/rev.hpp:11,
                 from stan/lib/stan_math/stan/math.hpp:148,
                 from stan/src/stan/model/model_header.hpp:4,
                 from /u/77/ave/unix/proj/ovarian/test.hpp:3:
stan/lib/stan_math/stan/math/prim/fun/multiply.hpp:29:51: note: candidate: template<int R, int C, class T1, class T2, class> Eigen::Matrix<typename stan::return_type<T1, T2>::type, R, C> stan::math::multiply(const Eigen::Matrix<T1, R, C>&, T2)
 inline Eigen::Matrix<return_type_t<T1, T2>, R, C> multiply(
                                                   ^
stan/lib/stan_math/stan/math/prim/fun/multiply.hpp:29:51: note:   template argument deduction/substitution failed:
/u/77/ave/unix/proj/ovarian/test.hpp:201:30: note:   mismatched types ‘const Eigen::Matrix<T1, R, C>’ and ‘double’
           elt_divide(multiply(pow(c, 2), square(lambda)),
                              ^
In file included from stan/lib/stan_math/stan/math/prim/fun.hpp:195:0,
                 from stan/lib/stan_math/stan/math/prim.hpp:10,
                 from stan/lib/stan_math/stan/math/rev.hpp:11,
                 from stan/lib/stan_math/stan/math.hpp:148,
                 from stan/src/stan/model/model_header.hpp:4,
                 from /u/77/ave/unix/proj/ovarian/test.hpp:3:
stan/lib/stan_math/stan/math/prim/fun/multiply.hpp:48:51: note: candidate: template<int R, int C, class T1, class T2, class> Eigen::Matrix<typename stan::return_type<T1, T2>::type, R, C> stan::math::multiply(T1, const Eigen::Matrix<T2, R, C>&)
 inline Eigen::Matrix<return_type_t<T1, T2>, R, C> multiply(
                                                   ^

Another even smaller example with slightly different error

data {
  int<lower=0> d;             // number of predictors
}
parameters {
  vector <lower=0>[d] lambda; // local shrinkage parameter
}
transformed parameters {
  vector[d] lambda_tilde;   // 'truncated' local shrinkage parameter
  lambda_tilde = sqrt( square(lambda) ./ (square(lambda)));
}
model {
}

With an error message (again just the beginning):

--- Compiling, linking C++ code ---
g++ -std=c++1y -pthread -D_REENTRANT -Wno-sign-compare      -I stan/lib/stan_math/lib/tbb_2019_U8/include -O3 -I src -I stan/src -I lib/rapidjson_1.1.0/ -I stan/lib/stan_math/ -I stan/lib/stan_math/lib/eigen_3.3.3 -I stan/lib/stan_math/lib/boost_1.72.0 -I stan/lib/stan_math/lib/sundials_4.1.0/include    -DBOOST_DISABLE_ASSERTS      -c  -x c++ -o /u/77/ave/unix/proj/ovarian/test.o /u/77/ave/unix/proj/ovarian/test.hpp
/u/77/ave/unix/proj/ovarian/test.hpp: In instantiation of ‘T__ test_model_namespace::test_model::log_prob(std::vector<T_l>&, std::vector<int>&, std::ostream*) const [with bool propto__ = false; bool jacobian__ = false; T__ = double; std::ostream = std::basic_ostream<char>]’:
stan/src/stan/model/model_base_crtp.hpp:147:29:   required from ‘double stan::model::model_base_crtp<M>::log_prob(std::vector<double>&, std::vector<int>&, std::ostream*) const [with M = test_model_namespace::test_model; std::ostream = std::basic_ostream<char>]’
/u/77/ave/unix/proj/ovarian/test.hpp:444:1:   required from here
/u/77/ave/unix/proj/ovarian/test.hpp:193:36: error: no matching function for call to ‘elt_divide(Eigen::CwiseUnaryOp<stan::math::apply_scalar_unary<F, T, typename std::enable_if<stan::is_eigen<S, void>::value, void>::type>::apply(const T&) [with F = stan::math::square_fun; T = Eigen::Matrix<double, -1, 1>]::<lambda(stan::math::apply_scalar_unary<stan::math::square_fun, Eigen::Matrix<double, -1, 1>, void>::scalar_t)>, const Eigen::Matrix<double, -1, 1> >, Eigen::CwiseUnaryOp<stan::math::apply_scalar_unary<F, T, typename std::enable_if<stan::is_eigen<S, void>::value, void>::type>::apply(const T&) [with F = stan::math::square_fun; T = Eigen::Matrix<double, -1, 1>]::<lambda(stan::math::apply_scalar_unary<stan::math::square_fun, Eigen::Matrix<double, -1, 1>, void>::scalar_t)>, const Eigen::Matrix<double, -1, 1> >)’
         stan::math::sqrt(elt_divide(square(lambda), square(lambda))),
                                    ^
In file included from stan/lib/stan_math/stan/math/prim/fun.hpp:78:0,
                 from stan/lib/stan_math/stan/math/prim.hpp:10,
                 from stan/lib/stan_math/stan/math/rev.hpp:11,
                 from stan/lib/stan_math/stan/math.hpp:148,
                 from stan/src/stan/model/model_header.hpp:4,
                 from /u/77/ave/unix/proj/ovarian/test.hpp:3:
stan/lib/stan_math/stan/math/prim/fun/elt_divide.hpp:23:44: note: candidate: template<class T1, class T2, int R, int C> Eigen::Matrix<typename stan::return_type<T1, T2>::type, R, C> stan::math::elt_divide(const Eigen::Matrix<T1, R, C>&, const Eigen::Matrix<T2, R, C>&)
 Eigen::Matrix<return_type_t<T1, T2>, R, C> elt_divide(
                                            ^
stan/lib/stan_math/stan/math/prim/fun/elt_divide.hpp:23:44: note:   template argument deduction/substitution failed:

Current Version:

• v2.22.0
• Linux t31300-lr010 4.15.0-74-generic stan-dev/cmdstan#83~16.04.1-Ubuntu SMP Wed Dec 18 04:56:23 UTC 2019 x86_64 x86_64 x86_64 GNU/Linux
• gcc (Ubuntu 5.4.0-6ubuntu1~16.04.12) 5.4.0 20160609

[wds15]

This is a parser related issue as it looks to me. I suggest to use the "Transfer issue" feature from GitHub to move this to stanc3 repo.

You should be able to use the old parser again if you want with some flags.

[avehtari]

Adding STANC2=true to make/local as instructed in CmdStan README.md didn't help.

[rok-cesnovar]

Ok, so that would mean its not a stanc3 issue, but an issue with the Math signatures. I will try and pinpoint this today, as we will have to make a patch release anyways.

[avehtari]

How do I can confirm whether stanc2 or stanc3 is actually used? bin/stanc --version and hpp files just say 2.22.0

[t4c1]

This one is on me. I overlooked this issue, when preparing #1471. I will prepare a fix soon. Meanwhile I think you can maybe work around this by assigning problematic expression to a variable, before using it (not tested).

[avehtari]

I was just testing stanc3 to catch problems before there's public announcement of new release. @t4c1 great that you know how to fix this!

[wds15]

Sorry for blaming stanc3... is this a bigger problem?

[rok-cesnovar]

How do I can confirm whether stanc2 or stanc3 is actually used? bin/stanc --version and hpp files just say 2.22.0

STANC2 should return 2.21.0 on bin/stanc --version.

After adding STANC2=true to make/local you have to run make clean-all and make build again.

Sorry for blaming stanc3... is this a bigger problem?

Not really, its a problem with the PR @t4c1 mentions, the move to general function signatures that accept Eigen expressions. Tadej explained to me that basically some functions still return Eigen expressions and because we still have some functions that dont accept them, that should not be used (and that is the fix he is working on). Once all functions are able to accept Eigen expression that will work and will be the preffered behavior.

We didnt catch it before because we dont have test models that would test all Stan functions in this context (without intermediate variables). Once all functions will accept Eigen expressions this wont be a problem. And using eigen expression is what we are aiming for ( I think we discussed this somewhere, I think on Discourse).

[avehtari]

STANC2 should return 2.21.0 on bin/stanc --version.
After adding STANC2=true to make/local you have to run make clean-all and make build again.

Yes, I had run make clean-all (twice), but I still got 2.22.0. Just in case also edited both cmdstan/make/local and stan/make/local

[rok-cesnovar]

Ok, will try to replicate that. You used the release cmdstan-2.22.0.tar.gz file + the 804 patch right? That is certainly something that should work and needs to be addressed immediately if it doesnt.

[avehtari]

You used the release cmdstan-2.22.0.tar.gz file + the 804 patch right?

I used cmdstan-2.22.0.tar.gz and manually removed that one extra paren. Two other things in the patch shouldn't have effect on this.