Type checking for ODE-solvers is broken

[VMatthijs]

Summary:

Some inputs to the ODE-solvers pass the Stan type check and fail on C++ compilation.

Description:

Data only argument status is not propagated correctly through function calls.

Reproducible Steps:

Try to compile

functions {

  real[] two_pool_feedback(real t, real[] C, real[] theta,
                           real[] x_r, int[] x_i) {
    real k1;
    real k2;
    real alpha21;
    real alpha12;

    real dC_dt[2];

    k1 <- theta[1];
    k2 <- theta[2];
    alpha21 <- theta[3];
    alpha12 <- theta[4];
    
    dC_dt[1] <- -k1 * C[1] + alpha12 * k2 * C[2];
    dC_dt[2] <- alpha21 * k1 * C[1] - k2 * C[2];

    return dC_dt;
  }

  real[] evolved_CO2(int N_t, real t0, real[] ts,
                     real gamma, real totalC_t0,
                     real k1, real k2, 
                     real alpha21, real alpha12,
                     real[] x_r, int[] x_i) {

    real C_t0[2];               // initial state
    real theta[4];              // ODE parameters
    real C_hat[N_t,2];          // predicted pool content

    real eCO2_hat[N_t];

    C_t0[1] <- gamma * totalC_t0;
    C_t0[2] <- (1 - gamma) * totalC_t0;

    theta[1] <- k1;
    theta[2] <- k2;
    theta[3] <- alpha21;
    theta[4] <- alpha12;

    C_hat <- integrate_ode_bdf(two_pool_feedback, 
                           C_t0, t0, ts, theta, x_r, x_i);

    for (t in 1:N_t)
      eCO2_hat[t] <- totalC_t0 - sum(C_hat[t]);
    return eCO2_hat;
  }

}
data {
  real<lower=0> totalC_t0;     // initial total carbon

  int<lower=0> N_t;            // number of measurement times

  vector<lower=0>[N_t] eCO2mean; // measured cumulative evolved carbon
  real<lower=0> eCO2sd[N_t];   // measured cumulative evolved carbon
}
transformed data {
  real x_r[0];                 // no real data for ODE system
  int x_i[0];                  // no integer data for ODE system
}
parameters {
  real<lower=0> k1;            // pool 1 decomposition rate
  real<lower=0> k2;            // pool 2 decomposition rate

  real<lower=0> alpha21;       // transfer coeff from pool 2 to 1
  real<lower=0> alpha12;       // transfer coeff from pool 1 to 2

  real<lower=0,upper=1> gamma; // partitioning coefficient

  real<lower=0> sigma;         // observation std dev

  real t0;                     // initial time

  real<lower=t0> ts[N_t];      // measurement times
  vector<lower=0>[N_t] eCO2;   // evolved CO2
}
transformed parameters {
  real eCO2_hat[N_t];
  eCO2_hat <- evolved_CO2(N_t, t0, ts, gamma, totalC_t0,
                          k1, k2, alpha21, alpha12, x_r, x_i);
}
model {
  gamma ~ beta(10,1);         // identifies pools
  k1 ~ normal(0,1);           // weakly informative
  k2 ~ normal(0,1);
  alpha21 ~ normal(0,1);
  alpha12 ~ normal(0,1);
  sigma ~ normal(0,1);

  // likelihood marginalizes out eCO2

  // true eCO2 (unknown) generated from system plus noise
  eCO2 ~ normal(eCO2_hat, sigma); 

  // measurement eCO2mean (observed) of true eCO2 value with
  // measurement noise eCO2sd (observed)
  eCO2mean ~ normal(eCO2, eCO2sd); 

}

It passes stanc, but fails on C++ compilation. Same with any of the other ODE-solvers.

Current Output:

No stanc type checking error.

Expected Output:

A stanc type checking error.

Additional Information:

This should also be fixed in the manual.

Current Version:

v2.18.1

[syclik]

Wow. Thanks for catching this. I've verified that it's a bug and it's been around since at least v2.16.0.

At some point I thought we had done this properly, but clearly not.

[VMatthijs]

Thanks for verifying! It's been fixed in stanc3.

[syclik]

We should fix this on develop if we're able to. I'll cc people later; I'm on my phone right now and can't remember handles.

…

On Fri, Jan 25, 2019 at 8:30 AM Matthijs Vákár ***@***.***> wrote: Thanks for verifying! It's been fixed in stanc3. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#2731 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAZ_F0Uy8eASuh16TLstVy0OAxJGPGi_ks5vGwcAgaJpZM4aRtpZ> .

[syclik]

cc: @bob-carpenter, @mitzimorris, @wds15, @seantalts

[seantalts]

I was in the room when @VMatthijs found this 😂

[mitzimorris]

assigning this myself with the best of intentions.

[bob-carpenter]

I believe there's a PR in progress for letting the initial time be variable.

…

The initial and non-initial time arguments to the ODE-solvers should be data only. This is the only thing that is currently implemented in the Math library.

[bob-carpenter]

This is a language issue. Are we going to keep patching the C++ transpiler or just wait for the OCaml version for all these fixes? There's a pile of these things that Matthijs has found already.

[syclik]

@bob-carpenter, let's have the general discussion on discourse: https://discourse.mc-stan.org/t/stan-language-do-we-patch-existing-c-transpiler-or-do-we-wait-for-the-ocaml-version/7436

Once we have a consistent position, we can report back on whether we'll fix this or not using the existing C++ transpiler.

[wds15]

there is a PR underway for some while now which will allow the initial time and the time vector to be var.

[VMatthijs]

This issue seems to have been solved by the PR stan-dev/math#1072 .