Ode integrator calls from user defined functions

Modeling
1

I’m trying to wrap up a call to integrate_ode_rk45 in a user function to avoid copy pasta for packing of domain parameters into the generic signature, e.g.

real[] solve(real[] t, real[] y, real tau, real k, real[] a, matrix fc) {
  int nn = size(a);
  real ys[size(t), nn];
  real ode_xr[0];
  int ode_xi[1] = { nn };
  real ode_par[2 + nn + nn*nn];
  ode_par[1] = tau;
  ode_par[2] = k;
  ode_par[3:nn+2] = a;
  ode_par[nn+2:nn+2+nn*nn] = to_array_1d(fc);
  return integrate_ode_rk45(de2, y, 0.0, ct, ode_par, ode_xr, ode_xi);
}

but I’m getting a message I find hard to believe,

sixth argument to integrate_ode_rk45 (real data) must be data only and not reference parameters
ERROR at line 38

 36:      ode_par[3:nn+2] = a;
 37:      ode_par[nn+2:nn+2+nn*nn] = to_array_1d(fc);
 38:      return integrate_ode_rk45(de2, y, 0.0, t, ode_par, ode_xr, ode_xi);
                                                                           ^
 39:    }

PARSER EXPECTED: ")"

whereas, if I copy and paste the function body into transformed parameters or generated quantities there’s no complaint.

Any tips beside #include?

2

Ouch, foiled by a length zero argument. You might try making ode_xr an argument to solve with a data real[] signature.

3

Thanks, but just to be sure I understand, I did try

real ode_xr[1] = { 0.0 } ;

and had the same error. I’ll try moving it to the func sig.

4

Don’t define ‘x_r’ as a local variable. Instead pass ‘rep_array(0.0, 0)’ into the sixth argument ‘integrate_ode_*’ directly.

1 Like
5

Thanks, that was unreasonably effective… is this a bug or a feature?

edit I’ve got something I like, so I promise not to complain about this anymore, e.g.

  real[] ode_par(real tau, real k, real[] a, matrix fc) {
    int nn = size(a);
    real par[2 + nn + nn*nn];
    par[1] = tau;
    par[2] = k;
    par[3:nn+2] = a;
    par[nn+2:nn+2+nn*nn] = to_array_1d(fc);
    return par;
  }

  real[,] solve(real[] t, real[] y, real tau, real k, real[] a, matrix fc) {
    return integrate_ode_rk45(
      de2, y, 0.0, t,
      ode_par(tau, k, a, fc),
      rep_array(0.0, 0),
      rep_array(size(a), 1)
    );
  }
}
6

It’s a documented “feature”, but it’s subtle to explain, especially in function contexts.

The core problem is that local variables will be autodiff variables rather than “data” variables if any of the types involved in the context (function arguments, or locals in transformed parameters or model) may be autodiff variables. The solution is tighter static analysis and code generation, where what you have would be inferred to be compatible with a double scalar type.

7

that’s totes reasonable; I would’ve expected the compiler to mention it, that’s all. I noticed there’s an possibly upcoming data qualifier (at least discourse stan syntax highlights it already?) so this will be easier I guess, going forward.

8

The data qualifier has been merged, so that’ll go into the next release (2.18). But that will just qualify arguments. What we need to do here is tighter type inference on the expressions on the right-hand side to infer the minimal type we can use for the declaration of the variable being assigned to. It’s not so easy, as a variable may be assigned to later, so it’s not something that can be done in one pass as our current compiler works.

9

If I might make a suggestion there: it might be useful to have a warning on use of uninitialized variables but not for parameters. For example I make egregious errors like

generated quantities {
  real oob;
  for (t in 1:T)
    if (abs(z[t] - zh[t]) > (2 * eps))
      oob = oob + 1;
}

since oob is not initialized to 0, this has no hope of working correctly, and a warning would be nice. This is a pretty obvious one, but use of uninitialized can produce more subtle problems.

1 Like
10

That’s on our to-do list. We have a bunch of things we’d like to put into a pedantic mode.

By the way, you can now do oob += 1;.

They’re all tricky in that they involve global analysis of the code.

1 Like