## Problem statement

I’m encountering a problem using external C++ code with MPI and `map_rect()`

. Briefly, I can compile and link fine when I avoid `map_rect()`

, but when I switch to using MPI + `map_rect()`

I get linker errors.

An example of the several errors I get is:

```
ok_wheat_pheno_analytic.hpp:(.text._ZN39ok_wheat_pheno_analytic_model_namespace15continuous_dydtIN4stan4math3varEdS3_EEN5Eigen6MatrixIN5boost4math5tools12promote_argsIT_T0_T1_fffE4typeELin1ELi1ELi0ELin1ELi1EEERKNS5_ISA_Lin1ELi1ELi0ELin1ELi1EEERKNS5_ISB_Lin1ELi1ELi0ELin1ELi1EEERKSt6vectorISC_SaISC_EEPSo[_ZN39ok_wheat_pheno_analytic_model_namespace15continuous_dydtIN4stan4math3varEdS3_EEN5Eigen6MatrixIN5boost4math5tools12promote_argsIT_T0_T1_fffE4typeELin1ELi1ELi0ELin1ELi1EEERKNS5_ISA_Lin1ELi1ELi0ELin1ELi1EEERKNS5_ISB_Lin1ELi1ELi0ELin1ELi1EEERKSt6vectorISC_SaISC_EEPSo]+0x325): undefined reference to `boost::math::tools::promote_args<stan::math::var, double, double, double, float, float>::type ok_wheat_pheno_analytic_model_namespace::sigmoid<stan::math::var, double, double, double>(stan::math::var const&, double const&, double const&, double const&, std::ostream*)'
```

I have encountered and solved (thanks to @bbbales2!) a similar problem before. I believe that the linker error is arising because the linker is looking for (and not finding) a function with signature:

```
inline stan::math::var sigmoid(const stan::math::var& x,
double sen,
double mid_pt,
double rate,
std::ostream* pstream__){
```

However, I don’t understand why it would be looking for a function with that signature. The `x`

argument is from data and so should be `double`

, but I believe that because these data are unpacked from a large 1-D array and assigned to a local variable they get “promoted” to `stan::math::var`

(that doesn’t trouble me much, although I’m open to suggestions to avoid it). All the other arguments (`sen`

, `mid_pt`

, and `rate`

) are parameters being estimated and should, I believe, only ever be `stan::math::var`

. Interestingly, the problem only manifests if I use:

```
target += sum( map_rect(parallel_log_likelihood, theta_vec, dummy, x_r_training, x_i_training) );
```

The problem goes away if I replace the above code with the following:

```
for(i in 1:Nshards){
target += parallel_log_likelihood(theta_vec, dummy[i], x_r_training[i,], x_i_training[i,]);
}
```

I could simply implement a C++ version that matches the signature it’s looking for, but this “solution” would gloss over what to me seems to be a more fundamental problem. Namely that `map_rect()`

is somehow “demoting” parameters to `double`

and thereby dropping gradient information.

I tried to create a minimal reproducible example of the problem, but have been unsuccessful so far.

Am I missing something? Am I somehow off-base in this interpretation? Any guidance would be much appreciated.

More details are provided below.

## Implemented external C++ functions

I have defined 3 versions of the `sigmoid()`

function in C++ the signatures of which are given here:

```
inline double sigmoid(double x,
double sen,
double mid_pt,
double rate,
std::ostream* pstream__){
```

```
inline stan::math::var sigmoid(double x,
const stan::math::var& sen,
const stan::math::var& mid_pt,
const stan::math::var& rate,
std::ostream* pstream__){
```

```
inline stan::math::var sigmoid(const stan::math::var& x,
const stan::math::var& sen,
const stan::math::var& mid_pt,
const stan::math::var& rate,
std::ostream* pstream__){
```

All three versions are in a header file, which I include via `USER_HEADER`

when compiling with `cmdstan`

. I don’t think I should even need the first version of the function (all `double`

arguments) given my description above regarding the arguments. Thus, I think these function signatures should be sufficient.

## Stan code

As mentioned briefly above, the only difference in Stan code between the two versions is in the model block where the code without MPI + `map_rect()`

is:

```
for(i in 1:Nshards){
target += parallel_log_likelihood(theta_vec, dummy[i], x_r_training[i,], x_i_training[i,]);
}
```

The code with MPI + `map_rect()`

is:

```
target += sum( map_rect(parallel_log_likelihood, theta_vec, dummy, x_r_training, x_i_training) );
```

I am literally commenting/uncommenting lines between running the compiler, so I’m very confident that this is the only difference.

In both cases, the arguments are the same: `theta_vec`

is a vector that contains transformed parameters (derived using a non-centered parameterization), `dummy`

is a length-`Nshards`

1-D array of length-0 vectors, `x_r_training`

is a `real`

2-D array with `Nshards`

rows, and `x_i_training`

is an `int`

2-D array with `Nshards`

rows.

`parallel_log_likelihood`

is a function that calls `run_model()`

, which calls `euler_integrate_ode()`

, which calls `continuous_dydt()`

, which calls the function in question `sigmoid()`

.