# Adding custom C++ function to R package using RStan

**URL:** <https://discourse.mc-stan.org/t/adding-custom-c-function-to-r-package-using-rstan/14289>\
**Category:** RStan\
**Created:** [April 13, 2020, 3:10pm UTC](https://discourse.mc-stan.org/t/adding-custom-c-function-to-r-package-using-rstan/14289 "2020-04-13T15:10:47Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![karimn](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/karimn/32/4974_2.png) [@karimn](https://discourse.mc-stan.org/u/karimn)\
**Post date:** [April 13, 2020, 3:10pm UTC](https://discourse.mc-stan.org/t/adding-custom-c-function-to-r-package-using-rstan/14289/1 "2020-04-13T15:10:48Z")

</div>

Hi,

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](https://mc-stan.org/rstantools/articles/minimal-rstan-package.html) 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](https://mc-stan.org/rstan/articles/external.html) 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.

Thanks,  
Karim

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [April 13, 2020, 4:33pm UTC](https://discourse.mc-stan.org/t/adding-custom-c-function-to-r-package-using-rstan/14289/2 "2020-04-13T16:33:12Z")

</div>

Can’t go into details right now, but we do something very similar in rstanarm

> **[stan-dev/rstanarm](https://github.com/stan-dev/rstanarm/tree/master/inst/include)**
>
> rstanarm R package for Bayesian applied regression modeling - stan-dev/rstanarm

that you can copy the structure of. If you used rstantools to create the package, then it should work.

---

<div class="post-metadata">

**Author:** ![karimn](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/karimn/32/4974_2.png) [@karimn](https://discourse.mc-stan.org/u/karimn)\
**Post date:** [April 13, 2020, 5:08pm UTC](https://discourse.mc-stan.org/t/adding-custom-c-function-to-r-package-using-rstan/14289/3 "2020-04-13T17:08:21Z")

</div>

Thanks, Ben, this is helpful. I did use rstantools.

---

<div class="post-metadata">

**Author:** ![karimn](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/karimn/32/4974_2.png) [@karimn](https://discourse.mc-stan.org/u/karimn)\
**Post date:** [April 14, 2020, 3:32am UTC](https://discourse.mc-stan.org/t/adding-custom-c-function-to-r-package-using-rstan/14289/4 "2020-04-14T03:32:28Z")

</div>

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;
}

```

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [April 14, 2020, 5:40am UTC](https://discourse.mc-stan.org/t/adding-custom-c-function-to-r-package-using-rstan/14289/5 "2020-04-14T05:40:57Z")

</div>

> [@karimn](#):
>
> is there any advantage to writing this in C++ compared to just doing the same loops in Stan function?

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.

---

<div class="post-metadata">

**Author:** ![karimn](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/karimn/32/4974_2.png) [@karimn](https://discourse.mc-stan.org/u/karimn)\
**Post date:** [April 14, 2020, 1:34pm UTC](https://discourse.mc-stan.org/t/adding-custom-c-function-to-r-package-using-rstan/14289/6 "2020-04-14T13:34:37Z")

</div>

Oh, no. I feared as much. Ok, I’ll need to learn how to handle this analytically. Thanks for your help.

BTW, I was going by what’s implemented in the below files. I noticed that it’s different from what is in [https://mc-stan.org/math/](https://mc-stan.org/math/).

> <https://github.com/stan-dev/math/blob/develop/stan/math/prim/fun/csr_matrix_times_vector.hpp>

  

> <https://github.com/stan-dev/math/blob/develop/stan/math/prim/fun/log_sum_exp.hpp>

---

<div class="post-metadata">

**Author:** ![karimn](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/karimn/32/4974_2.png) [@karimn](https://discourse.mc-stan.org/u/karimn)\
**Post date:** [April 14, 2020, 2:02pm UTC](https://discourse.mc-stan.org/t/adding-custom-c-function-to-r-package-using-rstan/14289/7 "2020-04-14T14:02:23Z")

</div>

Also, now finding the wiki here helpful, for anyone who has some of the same questions I have.

> **[stan-dev/math](https://github.com/stan-dev/math/wiki)**
>
> The Stan Math Library is a C++ template library for automatic differentiation of any order using forward, reverse, and mixed modes. It includes a range of built-in functions for probabilistic mode...