stan code parsers successfully but fails to compile due to pstream__ not declared error

[BlackEdder]

Summary:

I tried to compile a stan model, which was successfully parsed, but fails to compile due to:
./sirStrain.hpp:2822:40: error: ‘pstream__’ was not declared in this scope. Compilation was attempted both with rstan and cmdstan

Description:

See Summary

Reproducible Steps:

Try to compile with cmdStan (2.18.0):

make sirStrain

Note that the actual stan code is included below.

Current Output:

$ make sirStrain

--- Translating Stan model to C++ code ---
bin/stanc  sirStrain.stan --o=sirStrain.hpp
Model name=sirStrain_model
Input file=sirStrain.stan
Output file=sirStrain.hpp
DIAGNOSTIC(S) FROM PARSER:
Warning: left-hand side variable (name=healthy) occurs on right-hand side of assignment, causing inefficient deep copy to avoid aliasing.
Warning: left-hand side variable (name=healthy) occurs on right-hand side of assignment, causing inefficient deep copy to avoid aliasing.
Warning: left-hand side variable (name=y0) occurs on right-hand side of assignment, causing inefficient deep copy to avoid aliasing.
Warning: left-hand side variable (name=ili) occurs on right-hand side of assignment, causing inefficient deep copy to avoid aliasing.

Compiling pre-compiled header
g++ -Wall -I . -isystem stan/lib/stan_math/lib/eigen_3.3.3 -isystem stan/lib/stan_math/lib/boost_1.66.0 -isystem stan/lib/stan_math/lib/sundials_3.1.0/include -std=c++1y -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE -DBOOST_DISABLE_ASSERTS
 -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -Wno-unused-function -Wno-uninitialized -I src -isystem stan/src -isystem stan/lib/stan_math/ -DFUSION_MAX_VECTOR_SIZE=12 -Wno-unused-local-typedefs -DEIGEN_NO_DEBUG -DNO_FPRINTF_OUTPUT -pipe  -c -O3
 stan/src/stan/model/model_header.hpp -o stan/src/stan/model/model_header.hpp.gch

--- Linking C++ model ---
g++ -Wall -I . -isystem stan/lib/stan_math/lib/eigen_3.3.3 -isystem stan/lib/stan_math/lib/boost_1.66.0 -isystem stan/lib/stan_math/lib/sundials_3.1.0/include -std=c++1y -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE -DBOOST_DISABLE_ASSERTS
 -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -Wno-unused-function -Wno-uninitialized -I src -isystem stan/src -isystem stan/lib/stan_math/ -DFUSION_MAX_VECTOR_SIZE=12 -Wno-unused-local-typedefs -DEIGEN_NO_DEBUG -DNO_FPRINTF_OUTPUT -pipe   src/c
mdstan/main.cpp  -O3 -o sirStrain -include sirStrain.hpp stan/lib/stan_math/lib/sundials_3.1.0/lib/libsundials_nvecserial.a stan/lib/stan_math/lib/sundials_3.1.0/lib/libsundials_cvodes.a stan/lib/stan_math/lib/sundials_3.1.0/lib/libsundial
s_idas.a
In file included from <command-line>:0:0:
./sirStrain.hpp: In member function ‘void sirStrain_model_namespace::sirStrain_model::get_dims(std::vector<std::vector<long unsigned int> >&) const’:
./sirStrain.hpp:2822:40: error: ‘pstream__’ was not declared in this scope
         dims__.push_back(get_state_dim(pstream__));
                                        ^~~~~~~~~
./sirStrain.hpp:2822:40: note: suggested alternative: ‘pthread_t’
         dims__.push_back(get_state_dim(pstream__));
                                        ^~~~~~~~~
                                        pthread_t
In file included from <command-line>:0:0:
./sirStrain.hpp: In member function ‘void sirStrain_model_namespace::sirStrain_model::constrained_param_names(std::vector<std::__cxx11::basic_string<char> >&, bool, bool) const’:
./sirStrain.hpp:3159:56: error: ‘pstream__’ was not declared in this scope
             for (int k_0__ = 1; k_0__ <= get_state_dim(pstream__); ++k_0__) {
                                                        ^~~~~~~~~
./sirStrain.hpp:3159:56: note: suggested alternative: ‘pthread_t’
             for (int k_0__ = 1; k_0__ <= get_state_dim(pstream__); ++k_0__) {
                                                        ^~~~~~~~~
                                                        pthread_t
./sirStrain.hpp:3164:56: error: ‘pstream__’ was not declared in this scope
             for (int k_1__ = 1; k_1__ <= get_state_dim(pstream__); ++k_1__) {
                                                        ^~~~~~~~~
./sirStrain.hpp:3164:56: note: suggested alternative: ‘pthread_t’
             for (int k_1__ = 1; k_1__ <= get_state_dim(pstream__); ++k_1__) {
                                                        ^~~~~~~~~
                                                        pthread_t
./sirStrain.hpp:3171:56: error: ‘pstream__’ was not declared in this scope
             for (int k_0__ = 1; k_0__ <= get_theta_dim(pstream__); ++k_0__) {
                                                        ^~~~~~~~~
./sirStrain.hpp:3171:56: note: suggested alternative: ‘pthread_t’
             for (int k_0__ = 1; k_0__ <= get_theta_dim(pstream__); ++k_0__) {
                                                        ^~~~~~~~~
                                                        pthread_t
./sirStrain.hpp: In member function ‘void sirStrain_model_namespace::sirStrain_model::unconstrained_param_names(std::vector<std::__cxx11::basic_string<char> >&, bool, bool) const’:
./sirStrain.hpp:3253:56: error: ‘pstream__’ was not declared in this scope
             for (int k_0__ = 1; k_0__ <= get_state_dim(pstream__); ++k_0__) {
                                                        ^~~~~~~~~
./sirStrain.hpp:3253:56: note: suggested alternative: ‘pthread_t’
             for (int k_0__ = 1; k_0__ <= get_state_dim(pstream__); ++k_0__) {
                                                        ^~~~~~~~~
                                                        pthread_t
./sirStrain.hpp:3258:56: error: ‘pstream__’ was not declared in this scope
             for (int k_1__ = 1; k_1__ <= get_state_dim(pstream__); ++k_1__) {
                                                        ^~~~~~~~~
./sirStrain.hpp:3258:56: note: suggested alternative: ‘pthread_t’
             for (int k_1__ = 1; k_1__ <= get_state_dim(pstream__); ++k_1__) {
                                                        ^~~~~~~~~
                                                        pthread_t
./sirStrain.hpp:3265:56: error: ‘pstream__’ was not declared in this scope
             for (int k_0__ = 1; k_0__ <= get_theta_dim(pstream__); ++k_0__) {
                                                        ^~~~~~~~~
./sirStrain.hpp:3265:56: note: suggested alternative: ‘pthread_t’
             for (int k_0__ = 1; k_0__ <= get_theta_dim(pstream__); ++k_0__) {
                                                        ^~~~~~~~~
                                                        pthread_t
make/models:14: recipe for target 'sirStrain' failed
make: *** [sirStrain] Error 1

Expected Output:

A working model

Additional Information:

I tested this both from R (rstan) and cmdStan. In R it failed, but did not provide helpful information. Included here is the cmdStan output.

Current Version:

v2.18.0

The model consists of two files. The main file: sirStrain.stan

functions {
#include sirFunctions.stan
  real approxbin_lpdf(real k, real n, real p) {
    return beta_lpdf(p | k+1, n-k+1);
  }
  real niliForcing(real t, real amplitude, real period, real peakTime) {
    real sigma;
    sigma = 0.25*period;
    return amplitude*(exp(-pow(t-peakTime,2)/(2*pow(sigma,2)))-1);
  }
  real[] ode(real time, real[] state, real[] theta, real[] x_r, int[] x_i) {
    real dydt[get_state_dim()];
    real N;
    real healthy;
    healthy = N;
    N = get_N(x_r);
    for (k in 1:get_K(x_i)) {
      real newr;
      newr = get_gamma(theta)*get_I(state, {k});
      dydt[get_S_id({k})] = -get_beta(theta, {k})*get_S(state, {k})*get_I(state, {k})/N;
      dydt[get_I_id({k})] = get_beta(theta, {k})*get_S(state, {k})*get_I(state, {k})/N -
        newr;
      dydt[get_R_id({k})] = newr;
      dydt[get_PILI_id({k})] = newr*get_thetap(theta, {k});
      healthy = healthy - get_I(state, {k});
    }
    healthy = fmin(healthy, 1.0);
    dydt[get_NILI_id()] = exp(get_thetan(theta) +
        niliForcing(time, get_niliAmplitude(theta), get_niliPeriod(theta), get_niliPeak(theta))
      )*healthy;
    return dydt;
  }
}
data {
  int T;
  int K;
  real N;
  real ts[T];

  real sero_mu[K];
  real sero_sd[K];
  real niliPeak_mu;

  real thetap_alpha[K];
  real thetap_beta[K];

  int Y[T-1,K+1];
  real YILI[T,2];
}
transformed data {
  real x_r[get_xr_dim()];
  int x_i[get_xi_dim()];
  x_r = xr_to_1d(N);
  x_i = xi_to_1d(K);
}
parameters {
  real<lower=-10,upper=-6> I0[K];
  real<lower=0,upper=1> susc[K];
  real<lower=0,upper=1> asc;

  real<lower=0,upper=4> beta[K];
  real<lower=0,upper=1> gamma;
  real<lower=0,upper=1>  thetap[K];
  real<lower=-15,upper=0> thetan;

  real<lower=40,upper=200> niliPeriod;
  real<lower=0,upper=3> niliAmplitude;
  real<lower=14,upper=120> niliPeak;
}
transformed parameters {
  real y0[get_state_dim()];
  real y[T, get_state_dim()];
  real theta[get_theta_dim()];
  real pili[T-1,K];
  real nili[T-1];
  real ili[T-1];
  vector[K+1] ppili[T-1];

  theta = theta_to_1d(beta, gamma, thetap, thetan, niliPeriod, niliAmplitude, niliPeak);

  // initial
  y0[get_NILI_id()] = 0;
  for (k in 1:K) {
    y0[get_I_id({k})] = exp(I0[k])*susc[k]*N;
    y0[get_S_id({k})] = susc[k]*N - y0[get_I_id({k})];
    y0[get_R_id({k})] = (1-susc[k])*N;
    y0[get_PILI_id({k})] = 0;
  }

  y = integrate_ode_rk45(ode, y0, 0, ts, theta, x_r, x_i);

  // new nili and pili cases for last week
  for(t in 2:T) {
    nili[t-1] = get_NILI(y[t,]) - get_NILI(y[t-1,]);
    ili[t-1] = nili[t-1];

    for (k in 1:K) {
      pili[t-1,k] = get_PILI(y[t,], {k}) - get_PILI(y[t-1,], {k});
      ili[t-1] = ili[t-1] + pili[t-1,k];
    }
  }

  for (t in 1:(T-1)) {
    for (k in 1:K) {
      ppili[k, t] = pili[t,k]/ili[t];
    }
    ppili[K + 1, t] = nili[t]/ili[t];
  }
}
model {
  real sheddingPeriod;
  sheddingPeriod = 1/gamma;
  sheddingPeriod ~ normal(4.8, 0.245); // carrat_time_2008

  niliPeak ~ normal(niliPeak_mu + 7, 21); // +7 because we start measuring after one running weekp
  for (k in 1:K) beta[k] ~ normal(1.28*gamma*susc[k], 0.133); // Assume biggerstaff is on Reff
  for (k in 1:K) susc[k] ~ normal(sero_mu[k], sero_sd[k]);
  for (k in 1:K) thetap[k] ~ beta(thetap_alpha[k], thetap_beta[k]);

  asc ~ beta(35.644, 69.314);

  for (t in 1:(T-1)) {
    YILI[t,1] ~ approxbin(ili[t]*YILI[t,2]/N, asc);
    Y[t] ~ multinomial(ppili[t]);
  }
}

and the included file sirFunctions.stan:

int index_to_1d(int[] index, int[] dims, int offset);
int index_to_1d(int[] index, int[] dims, int offset) {
  int n;
  n = size(dims);
  if (n == 1) {
    return offset + index[1];
  } else {
    //return (index[n] - 1) * prod(head(dims, n - 1)) +
    //  index_to_1d(head(index, n - 1), head(dims, n - 1), offset);
    return (index[1] - 1) * prod(tail(dims, n - 1)) +
      index_to_1d(tail(index, n - 1), tail(dims, n - 1), offset);
  }
}

int get_beta_id(int[] index) { return index_to_1d(index, {3} , 0); }
real get_beta(real[] v, int[] index) { return v[get_beta_id(index)]; }
real[] get_beta_array(real[] v) { real arr[3];
              for (id1 in 1:3) arr[id1] = get_beta(v, {id1});
              return arr; }
int get_gamma_id() { return 4; }
real get_gamma(real[] v) { return v[get_gamma_id()]; }
real get_gamma_array(real[] v) { return get_gamma(v); }
int get_thetap_id(int[] index) { return index_to_1d(index, {3} , 4); }
real get_thetap(real[] v, int[] index) { return v[get_thetap_id(index)]; }
real[] get_thetap_array(real[] v) { real arr[3];
              for (id1 in 1:3) arr[id1] = get_thetap(v, {id1});
              return arr; }
int get_thetan_id() { return 8; }
real get_thetan(real[] v) { return v[get_thetan_id()]; }
real get_thetan_array(real[] v) { return get_thetan(v); }
int get_niliPeriod_id() { return 9; }
real get_niliPeriod(real[] v) { return v[get_niliPeriod_id()]; }
real get_niliPeriod_array(real[] v) { return get_niliPeriod(v); }
int get_niliAmplitude_id() { return 10; }
real get_niliAmplitude(real[] v) { return v[get_niliAmplitude_id()]; }
real get_niliAmplitude_array(real[] v) { return get_niliAmplitude(v); }
int get_niliPeak_id() { return 11; }
real get_niliPeak(real[] v) { return v[get_niliPeak_id()]; }
real get_niliPeak_array(real[] v) { return get_niliPeak(v); }
real[] theta_to_1d(real[] beta, real gamma, real[] thetap, real thetan, real niliPeriod, real niliAmplitude, real niliPeak) { real v[11]; v[1:3] = to_array_1d(beta); v[4] = gamma; v[5:7] = to_array_1d(thetap); v[8] = thetan; v[9] = niliPeriod; v[10] = niliAmplitude; v[11] = niliPeak; return v; }
int get_theta_dim() { return 11; }
int get_S_id(int[] index) { return index_to_1d(index, {3} , 0); }
real get_S(real[] v, int[] index) { return v[get_S_id(index)]; }
real[] get_S_array(real[] v) { real arr[3];
              for (id1 in 1:3) arr[id1] = get_S(v, {id1});
              return arr; }
int get_I_id(int[] index) { return index_to_1d(index, {3} , 3); }
real get_I(real[] v, int[] index) { return v[get_I_id(index)]; }
real[] get_I_array(real[] v) { real arr[3];
              for (id1 in 1:3) arr[id1] = get_I(v, {id1});
              return arr; }
int get_R_id(int[] index) { return index_to_1d(index, {3} , 6); }
real get_R(real[] v, int[] index) { return v[get_R_id(index)]; }
real[] get_R_array(real[] v) { real arr[3];
              for (id1 in 1:3) arr[id1] = get_R(v, {id1});
              return arr; }
int get_NILI_id() { return 10; }
real get_NILI(real[] v) { return v[get_NILI_id()]; }
real get_NILI_array(real[] v) { return get_NILI(v); }
int get_PILI_id(int[] index) { return index_to_1d(index, {3} , 10); }
real get_PILI(real[] v, int[] index) { return v[get_PILI_id(index)]; }
real[] get_PILI_array(real[] v) { real arr[3];
              for (id1 in 1:3) arr[id1] = get_PILI(v, {id1});
              return arr; }
real[] state_to_1d(real[] S, real[] I, real[] R, real NILI, real[] PILI) { real v[13]; v[1:3] = to_array_1d(S); v[4:6] = to_array_1d(I); v[7:9] = to_array_1d(R); v[10] = NILI; v[11:13] = to_array_1d(PILI); return v; }
int get_state_dim() { return 13; }
int get_N_id() { return 1; }
real get_N(real[] v) { return v[get_N_id()]; }
real get_N_array(real[] v) { return get_N(v); }
real[] xr_to_1d(real N) { real v[1]; v[1] = N; return v; }
int get_xr_dim() { return 1; }
int get_K_id() { return 1; }
int get_K(int[] v) { return v[get_K_id()]; }
int get_K_array(int[] v) { return get_K(v); }
int[] xi_to_1d(int K) { int v[1]; v[1] = K; return v; }
int get_xi_dim() { return 1; }

EDIT: clarify text

[BlackEdder]

This seems to be caused by the fact that you can't use a function to set vector size. Is this a "deliberate" design choice or a bug?

I was trying to use constants to set the size of my arrays instead of hard coded values (sirFunctions.stan is auto generated). The manual says the way to define constants in stan is:

Constants in Stan are nothing more than nullary (no-argument) functions. For in-
stance, the mathematical constants π and e are represented as nullary functions
named pi() and e(). See Section 41.2 for a list of built-in constants.

So I assumed this would mean you could use such constants when defining vector/array size.

[bob-carpenter]

Thanks for reporting. Models shouldn't translate to C++ and then not compile, so that's a bug.

As far as the language design, the intention is that a function may be used to determine a size for a block variable or local variable as long as the function is constant---it can't change values based on parameters.

The situation that's failing here is using it to define a local variable in a function:

int get_state_dim() { return 13; }
...
real[] ode(real time, real[] state, real[] theta, real[] x_r, int[] x_i) {
    real dydt[get_state_dim()];
...

The problem is that it's not getting a print stream passed in.

In file included from <command-line>:0:0:
./sirStrain.hpp: In member function ‘void sirStrain_model_namespace::sirStrain_model::get_dims(std::vector<std::vector<long unsigned int> >&) const’:
./sirStrain.hpp:2822:40: error: ‘pstream__’ was not declared in this scope
         dims__.push_back(get_state_dim(pstream__));
                                        ^~~~~~~~~

We'll fix this.

[BlackEdder]

Thanks!

Just as a follow-up. This does work as expected in the transformed data section, but not in the transformed parameters section or within a function.

[bob-carpenter]

Here's a smaller example of both where I get the same failure to compile:

functions {
  int sizer() {
    return 5;
  }
  real[] foo() {
    real x[sizer()];
    return x;
  }
}
transformed data {
  real x[sizer()] = foo();
}
transformed parameters {
  real z[sizer()];
  z = x;
}

[rok-cesnovar]

Still an issue in stanc3. Transferring.