stan-dev / stan

Stan development repository. The master branch contains the current release. The develop branch contains the latest stable development. See the Developer Process Wiki for details.
https://mc-stan.org
BSD 3-Clause "New" or "Revised" License
2.56k stars 366 forks source link

min:max indexing does not support SoA but does not prevent this at compile time? #3244

Closed WardBrian closed 8 months ago

WardBrian commented 8 months ago

Summary:

Originally reported by @avehtari using a larger model. This is breaks sampling behavior when using STANCFLAGS=--O1

Reproducible Steps:

transformed data {
  int N = 10;
  int<upper=N> M = 10;
}
parameters {
  vector[N] x;
}
transformed parameters {
  vector[M] y;
  y = x[1 : M];
}

model {
  y ~ normal(0, 1);
}

When compiled with cmdstan and sampled, this does perfectly well. When you re-compile with STANCFLAGS=--O1, which flips on SoA for x and y, the sampling behavior becomes incredibly slow and the resulting ESS is horrible.

Current Output:

See above

Expected Output:

Either this should work, or if it is not intended to be used with SoA I would expect C++ compilation to fail. There are no uses of index_min_max in https://github.com/stan-dev/stan/blob/develop/src/stan/model/indexing/rvalue_varmat.hpp and no guards against it in https://github.com/stan-dev/stan/blob/develop/src/stan/model/indexing/rvalue.hpp, so those are the signatures being used

Additional Information:

Current Version:

v2.33.0

bob-carpenter commented 8 months ago

Indexing into an SOA matrix should be efficient, so it shouldn't be too expensive to create the non-SOA matrix x[1:M]. Then we can assign it to the SOA matrix y, we just need the chain() on the SOA data type to add all the adjoints of y to the adjoints of x. Not super efficient, but it shouldn't be dreadful.

SteveBronder commented 8 months ago

That's odd, it should work. Let me look at this today it sounds like a bug in the rvalue function for min max indexing

WardBrian commented 8 months ago

I just opened another issue on Math, but worth noting here that replacing the indexing with a call to head or segment also leads to failure (though in that case, it is a segfault rather than junk inference)

SteveBronder commented 8 months ago

Okay yeah digging into this it's an issue in the math library. I'm going to put up a PR in a minute

SteveBronder commented 8 months ago

Closing as this was a bug in math