C++ compilation error when using slicing in recursive function call

[spinkney]

The following recursive function fails to compile. I believe this happens because if D - 1 == 1 then the index is 1:1 and the container becomes a real. Wrapping gamma[1:D - 1] in to_vector() allows the model to compile.

functions {
vector test2(vector gamma);
 vector test2(vector gamma) {
   int D = num_elements(gamma);

   if (D == 1)
      return rep_vector(D, 0);
    else
      return test2(gamma[1:D - 1]);
  }

}
data {
  int N;
  int times;
}
parameters {
  vector[times] gamma;
}
transformed parameters {
  vector[N] z_hat = test2(gamma);
}

@WardBrian @rok-cesnovar

[WardBrian]

The error specifically mentions

stan/lib/stan_math/lib/eigen_3.3.9/Eigen/src/Core/MatrixBase.h:48:34: fatal error: template instantiation depth exceeds maximum of 900 (use ‘-ftemplate-depth=’ to increase the maximum)
   48 | template<typename Derived> class MatrixBase

and a ton of nested Eigen::Block functions. Wrapping this in to_vector seems to stop the recursive template instantiation.

I'm not sure what we can do to resolve this to be honest

[andrjohns]

A quick fix could be to add -ftemplate-depth=1500 (or whichever limit) to the Math library default compiler flags

[WardBrian]

I think the problem is the code generated by this will genuinely recurse infinitely in the C++ compiler. Each new template expansion adds another wrapping Eigen::Block<…> which then leads to a new template expansion

[andrjohns]

Ahhh of course. Would it be tricky on the Stanc3 side of things to always eval recursive functions? Or having a plain_type_t<> return type maybe?

[WardBrian]

I think it would be tricky, or at least quiet messy. By the time we're doing code generation we don't want to be trying to do function resolution to figure out if a call is recursive (as opposed to an overloaded call, or just a normal function call).

If there was something we could do in terms of the return types or declarations to resolve this it would be great.

For some reason this doesn't have the same problem when the code is test2(gamma + gamma), even though I would expect that to create an infinite chain of CWiseUnaryOps. Any idea what the difference is @andrjohns? Could we try to update the templating of rvalue to avoid this?