stan-dev / stanc3

The Stan transpiler (from Stan to C++ and beyond).
BSD 3-Clause "New" or "Revised" License
141 stars 46 forks source link

Use Eigen Initialize to NaN macro instead of looping over matrix #514

Closed SteveBronder closed 3 years ago

SteveBronder commented 4 years ago

Summary

When we initialize containers of parameters the compiler usually generates code like

      for (size_t sym1__ = 1; sym1__ <= N; ++sym1__) {
        current_statement__ = 10;
        for (size_t sym2__ = 1; sym2__ <= M; ++sym2__) {
          current_statement__ = 10;
          assign(xx,
            cons_list(index_uni(sym1__),
              cons_list(index_uni(sym2__), nil_index_list())),
            std::numeric_limits<double>::quiet_NaN(), "assigning variable xx");
        }
  }

What if instead of doing this ourselves we just used the EIGEN_INITIALIZE_MATRICES_BY_NAN? I'm guessing Eigen would have some tricks to fill these out. It might even be a nice optimization to turn off NaN initialization

bob-carpenter commented 4 years ago

For reverse-mode autodiff types, we need to do it ourselves so that we don't make extra copies. What we want is more like this:

var nan = std::numeric_limits<double>::quiet_NaN();
Matrix<var, -1, -1> a(M, N);
for (int i = 0; i < a.size(); ++i) a(i) = nan;

We could even go further by making sure the vari created for nan doesn't go on the autodiff stack.

We could probably get away with

var nan = std::numeric_limits<double>::quiet_NaN();
Matrix<T, -1, -1> a = Eigen::Matrix<T, -1, -1>::Constant(nan);

The way it's done now is extra bad because those assigns have to do integer arithmetic in there to convert the indexing-by-1 back to indexing-by-0.

Same for intializing vector<var>.

SteveBronder commented 4 years ago

Hmm, I think the first one is v possible. We'd have to do something like in #433 where the value is declared up front

SteveBronder commented 4 years ago

I tried taking a stab at this but it was a lil' over my head. though yeah it looks like you can do almost the same thing as with the flat stuff

nhuurre commented 4 years ago

@SteveBronder can you take a look at #493 ? The latest commit removes the loop in favor of stan::math::fill(). I guess Eigen::Matrix::Constant() would be even better.

SteveBronder commented 4 years ago

Nice! Yeah I think stan::math::fill makes a lot of sense. On the stan math backend imo I think it's a lot easier to write the optimized version