I’m trying write my own C++ function to add to my model that I’m putting in an R package. I followed the step-by-step guide and everything works great, but now that I’m trying to figure out where to put a new C++ function. I read the instructions on how to interface with external C++ code but I’m not sure how to apply that to package. (BTW, that RStan guide on how to interface with C++ seems to have an error: there’s a line that reads “The sinc.hpp file reads as” and then there is no code after it).

Some background on why I’m trying to do this: I’m want to create my own custom version of csr_matrix_times_vector() where the vector is actually a vector of log probabilities and I need to something like log_sum_exp() when I do the dot product step. I’m looking at the code for these functions on github and I think I understand what is going on but if anyone knows of some pitfalls I should avoid, please let me know.

I managed to successfully build my package having added a new C++ function. I kind of mixed what csr_matrix_times_vector() and log_sum_exp() do. I was wondering if you see anything obviously wrong with this code? This is my first time dabbling with this. This is what I want to do: given a dense n-vector of log probabilities, I want to calculate the log_sum_exp() of m, of their subsets using a sparse m\times n boolean matrix, where each row identifies which log probabilities to calculate over. Also, is there any advantage to writing this in C++ compared to just doing the same loops in Stan function?

Thanks,
Karim

template<typename T>
inline
Eigen::Matrix<T, Eigen::Dynamic, 1>
csr_log_sum_exp(int m, int n,
const std::vector<int>& v,
const std::vector<int>& u,
const Eigen::Matrix<T, Eigen::Dynamic, 1>& b,
std::ostream* pstream__) {
Eigen::Matrix<T, Eigen::Dynamic, 1> result(m);
result.setZero();
for (int row = 0; row < m; ++row) {
int idx = csr_u_to_z(u, row); // u[row + 1] - u[row]
int row_end_in_w = (u[row] - stan::error_index::value) + idx;
Eigen::Matrix<T, Eigen::Dynamic, 1> b_sub(idx);
b_sub.setZero();
int i = 0;
for (int nze = u[row] - stan::error_index::value; nze < row_end_in_w; ++nze, ++i) {
b_sub.coeffRef(i) = b.coeffRef(v[nze] - stan::error_index::value);
}
const double max = b_sub.maxCoeff();
result.coeffRef(row) = max + std::log((b_sub.array() - max).exp().sum());
}
return result;
}

I think it is actually worse. If you are going to the trouble of writing custom C++, then you need to be handling the derivatives analytically rather than letting it autodiff.