# Link Symbolicc++ function to Stan interface

**URL:** <https://discourse.mc-stan.org/t/link-symbolicc-function-to-stan-interface/24216>\
**Category:** Modeling\
**Created:** [September 2, 2021, 7:31am UTC](https://discourse.mc-stan.org/t/link-symbolicc-function-to-stan-interface/24216 "2021-09-02T07:31:56Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![chI-No](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/chi-no/32/12632_2.png) [@chI-No](https://discourse.mc-stan.org/u/chI-No)\
**Post date:** [September 2, 2021, 7:31am UTC](https://discourse.mc-stan.org/t/link-symbolicc-function-to-stan-interface/24216/1 "2021-09-02T07:31:56Z")

</div>

My question is simple. How we can add a function created via symbolic derivative from SymbolicC++ to the stan program in a mathematica interface? Can i load this library indirectly in Stan?

Thanks!

---

<div class="post-metadata">

**Author:** ![andrjohns](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/andrjohns/32/15297_2.png) [@andrjohns](https://discourse.mc-stan.org/u/andrjohns)\
**Post date:** [September 2, 2021, 8:37am UTC](https://discourse.mc-stan.org/t/link-symbolicc-function-to-stan-interface/24216/2 "2021-09-02T08:37:58Z")

</div>

Stan isn’t able to use external packages, unfortunately. You would need to specify your function using Stan syntax in the `functions{}` block of your Stan model

---

<div class="post-metadata">

**Author:** ![chI-No](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/chi-no/32/12632_2.png) [@chI-No](https://discourse.mc-stan.org/u/chI-No)\
**Post date:** [September 2, 2021, 4:56pm UTC](https://discourse.mc-stan.org/t/link-symbolicc-function-to-stan-interface/24216/3 "2021-09-02T16:56:38Z")

</div>

Thats very bad. Then how can i perform a derivative of a function in stan? Cause i need to use a derivative of a complex function to fit my model

---

<div class="post-metadata">

**Author:** ![chI-No](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/chi-no/32/12632_2.png) [@chI-No](https://discourse.mc-stan.org/u/chI-No)\
**Post date:** [September 2, 2021, 5:17pm UTC](https://discourse.mc-stan.org/t/link-symbolicc-function-to-stan-interface/24216/4 "2021-09-02T17:17:06Z")

</div>

Another solution is to perform the symbolic derivative and just copy explicit the result in another archive. Can i at least load an external c++ archive to the stan interface in mathematica?

---

<div class="post-metadata">

**Author:** ![andrjohns](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/andrjohns/32/15297_2.png) [@andrjohns](https://discourse.mc-stan.org/u/andrjohns)\
**Post date:** [September 3, 2021, 12:44am UTC](https://discourse.mc-stan.org/t/link-symbolicc-function-to-stan-interface/24216/5 "2021-09-03T00:44:15Z")

</div>

No that’s not possible. Stan uses custom data types and functions within its c++, so you can’t interface with external c++ libraries very easily.

You need to know the functional form of the derivative and specify that in Stan, unfortunately

---

<div class="post-metadata">

**Author:** ![chI-No](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/chi-no/32/12632_2.png) [@chI-No](https://discourse.mc-stan.org/u/chI-No)\
**Post date:** [September 3, 2021, 2:50am UTC](https://discourse.mc-stan.org/t/link-symbolicc-function-to-stan-interface/24216/6 "2021-09-03T02:50:39Z")

</div>

Yes that was what im trying to say. I mean compute the derivative and obtain the explicit form of the function. Then just copy that and put in stan.

My question is just if i can copy the explicit function in an archive.txt and call this archive with stan, not using any external library just call a explicit function.

---

<div class="post-metadata">

**Author:** ![andrjohns](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/andrjohns/32/15297_2.png) [@andrjohns](https://discourse.mc-stan.org/u/andrjohns)\
**Post date:** [September 3, 2021, 2:57am UTC](https://discourse.mc-stan.org/t/link-symbolicc-function-to-stan-interface/24216/7 "2021-09-03T02:57:44Z")

</div>

The function still needs to be specified in a language that is compatible with Stan. You can either specify the function using c++ and then include that function within Stan, or just write the function in pure Stan syntax;

But that is not always a straightforward process, if the derivative requires any functions outside of the c++ standard library, then you also need to specify the derivative of the derivative function for Stan to use it for sampling. If your derivative function doesn’t need any specialised functions, then it is generally much easier to respecify it in pure Stan.

Out of interest, what is the derivative function that you’re trying to use?

---

<div class="post-metadata">

**Author:** ![chI-No](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/chi-no/32/12632_2.png) [@chI-No](https://discourse.mc-stan.org/u/chI-No)\
**Post date:** [September 3, 2021, 8:48am UTC](https://discourse.mc-stan.org/t/link-symbolicc-function-to-stan-interface/24216/8 "2021-09-03T08:48:47Z")

</div>

Im gonna copy here my program. Im trying to estimate parameters of a cosmological void model with supernovae data. This is the stan code with de definition of all functions that i need to use:

stanCode = "  
functions {  
real OmegaM(real r, // Comovil Position  
real Omegaout, // Omega out  
real Omegain, // Omega in  
real r0, // Void size parameter  
real Deltar // Change parameter  
) {  
return Omegaout+(Omegain-Omegaout)_((1-tanh((r-r0)/(2.0_   
Deltar)))/(1+tanh(r0/(2.0\* Deltar))));  
}

real OmegaK(real r, // Comovil Position  
real Omegaout, // Omega out  
real Omegain, // Omega in  
real r0, // Void size parameter  
real Deltar // Change parameter  
) {  
return 1-OmegaM(r,Omegaout,Omegain,r0,Deltar);  
}

real tBB(real r, // Comovil Position  
real H0 // Local Hubble Rate  
) {  
real c = 1;  
return c/H0;  
}

real Hubble0(real r, // Comovil Position  
real H0 , // Local Hubble Rate  
real Omegaout, // Omega out  
real Omegain, // Omega in  
real r0, // Void size parameter  
real Deltar // Change parameter  
) {  
return (1/tBB(r,H0))_((1/OmegaK(r,Omegaout,Omegain,r0,Deltar))-(  
OmegaM(r,Omegaout,Omegain,r0,Deltar)/sqrt(OmegaK(r,Omegaout,Omegain,  
r0,Deltar))^3))_  
asinh(sqrt(OmegaK(r,Omegaout,Omegain,r0,Deltar)/OmegaM(r,Omegaout,  
Omegain,r0,Deltar)));  
}

real R0(real r // Comovil Position  
) {  
return r;  
}

real M(real r, // Comovil Position  
real H0 , // Local Hubble Rate  
real Omegaout, // Omega out  
real Omegain, // Omega in  
real r0, // Void size parameter  
real Deltar // Change parameter  
) {  
real G=1;  
return (1/(2.0_G))_(Hubble0(r,H0,Omegaout,Omegain,r0,Deltar))^2 \*   
OmegaM(r,Omegaout,Omegain,r0,Deltar)\*(R0(r))^3;  
}

real k(real r, // Comovil Position  
real H0 , // Local Hubble Rate  
real Omegaout, // Omega out  
real Omegain, // Omega in  
real r0, // Void size parameter  
real Deltar // Change parameter  
) {  
return (Hubble0(r,H0,Omegaout,Omegain,r0,Deltar))^2\*(OmegaM(r,  
Omegaout,Omegain,r0,Deltar)-1)\*(R0(r))^3;  
}

real InverseSerie(real y // variable  
) {  
return (3.0/2)^(2.0/3)\* y^(2.0/3)+(3.0/10)\* (3.0/2)^(1.0/3)\*   
y^(4.0/3)-(27.0/700)_y^2+(23.0/3500)_(3.0/2)^(2.0/3)_y^(8.0/3)-(2841.  
0/1247500)_(3.0/2)^(1.0/3)\*y^(10.0/3);  
}

real contrast(real r, // Comovil Position  
real Omegaout, // Omega out  
real Omegain, // Omega in  
real r0, // Void size parameter  
real Deltar // Change parameter  
) {  
return OmegaK(r,Omegaout,Omegain,r0,Deltar)/OmegaM(r,Omegaout,  
Omegain,r0,Deltar);  
}

real a(real r, // Comovil Position  
real t, // Comovil Time  
real H0 , // Local Hubble Rate  
real Omegaout, // Omega out  
real Omegain, // Omega in  
real r0, // Void size parameter  
real Deltar // Change parameter  
) {  
return ((3.0/2)\* Hubble0(r,H0,Omegaout,Omegain,r0,Deltar)\*   
sqrt(OmegaM(r,Omegaout,Omegain,r0,Deltar))\* t)^(2.0/3);  
}

real R(real r, // Comovil Position  
real t, // Comovil Time  
real H0 , // Local Hubble Rate  
real Omegaout, // Omega out  
real Omegain, // Omega in  
real r0, // Void size parameter  
real Deltar // Change parameter  
) {  
return (R0(r)/contrast(r,Omegaout,Omegain,r0,Deltar))_InverseSerie((  
2.0/3)_(contrast(r,Omegaout,Omegain,r0,Deltar)\*a(r,t,H0,Omegaout,  
Omegain,r0,Deltar))^(3.0/2));  
}

real dR\_dr(real r, // Comovil Position  
real t, // Comovil Time  
real H0 , // Local Hubble Rate  
real Omegaout, // Omega out  
real Omegain, // Omega in  
real r0, // Void size parameter  
real Deltar // Change parameter  
) {  
return dR\_dr;  
}

real dR\_drt(real r, // Comovil Position  
real t, // Comovil Time  
real H0 , // Local Hubble Rate  
real Omegaout, // Omega out  
real Omegain, // Omega in  
real r0, // Void size parameter  
real Deltar // Change parameter  
) {  
return dR\_drt;  
}

real dL(real z, // Redshift  
real r, // Comovil Position  
real t, // Comovil Time  
real H0 , // Local Hubble Rate  
real Omegaout, // Omega out  
real Omegain, // Omega in  
real r0, // Void size parameter  
real Deltar // Change parameter  
) {  
return (1.0+z)^2 \* R(r,t,H0,Omegaout,Omegain,r0,Deltar);  
}

real Ev(real z // Redshift  
) {  
return z;  
}

real mb(real z, // Redshift  
real r, // Comovil Position  
real t, // Comovil Time  
real H0 , // Local Hubble Rate  
real Omegaout, // Omega out  
real Omegain, // Omega in  
real r0, // Void size parameter  
real Deltar // Change parameter  
) {  
return 5\*log10(dL(z,r,t,H0,Omegaout,Omegain,r0,Deltar))+Ev(z);  
}

real[] GeodesicEquation(real z, // Redshift  
real[] c, // Comovil Coordinates  
real[] theta, // Parameters  
real[] x\_r, // unused data  
int[] x\_i  
) {  
real r=c[1];  
real t=c[2];

```
real H0 = theta[1];   
real Omegaout = theta[2]; 
real Omegain = theta[3];   
real r0 = theta[4];        
real Deltar = theta[5];    

real dt_dz=(dR_dr)/((1+z)*(dR_drt) // Not defined yet;
real dr_dz=sqrt(1-k(r,H0,Omegaout,Omegain,r0,Deltar))/((1+z)*(dR_drt)\

```

// Not defined yet;  
return {dr\_dz,dt\_dz};  
}  
}  
data {  
int\<lower = 0\> N; // number of measurement  
real zobs[N]; // Redshifts  
real\<lower = 0\> msobs[N]; // Supernovaes magnitude  
real\<lower = 0\> errorm[N]; // Supernovaes error  
}

parameters {  
real\<lower = 0\> H0; // Hubble local rate  
real\<lower = 0\> Omegain; // Local density matter  
real\<lower = 0\> r0; // Void size parameter  
real\<lower = 0\> Deltar; // Change parameter  
real\<lower = 0\> sigma; // Error  
real\<lower = 0\> alpha; // Evolution parameter  
real\<lower = 0\> beta; // Evolution parameter  
}

transformed parameters {  
real c=1;  
real tG0=c/H0;  
real rG0=0;  
real ComovilCoordinates[N,2]  
= integrate\_ode\_rk45(GeodesicEquation, {tG0,rG0}, 0, zobs,   
{H0,1.0,Omegain,r0,Deltar},  
rep\_array(0.0, 0), rep\_array(0, 0),  
1e-5, 1e-3, 5e2);  
}

model {  
H0 ~ normal(1, 0.5);  
Omegain ~ normal(0.05, 0.05);  
r0 ~ lognormal(-1, 1);  
Deltar ~ lognormal(log(10), 1);  
sigma ~ lognormal(log(10), 1);  
alpha ~ lognormal(log(10), 1);  
beta ~ lognormal(log(10), 1); //not true priors

for(i in 1:N)  
msobs[i] ~   
normal(mb(zobs[i],ComovilCoordinates[i,2],ComovilCoordinates[i,1],H0,  
1.0,Omegain,r0,Deltar), sigma);  
}";

You can see that i need te derivatives of the function R (dR\_dr, dR\_drt) to define my differential equation in the z var and obtain the mb function that i need

---

<div class="post-metadata">

**Author:** ![chI-No](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/chi-no/32/12632_2.png) [@chI-No](https://discourse.mc-stan.org/u/chI-No)\
**Post date:** [September 3, 2021, 8:50am UTC](https://discourse.mc-stan.org/t/link-symbolicc-function-to-stan-interface/24216/9 "2021-09-03T08:50:53Z")

</div>

Overall, all the quatities just depend on the functions:

OmegaM(r)=Omegaout+(Omegain-Omegaout) _((1-tanh((r-r0)/(2.0_  
Deltar)))/(1+tanh(r0/(2.0\* Deltar))))

Hubble0(r)= (1/OmegaK(r,Omegaout,Omegain,r0,Deltar))-(  
OmegaM(r,Omegaout,Omegain,r0,Deltar)/sqrt(OmegaK(r,Omegaout,Omegain,  
r0,Deltar))^3))\*  
asinh(sqrt(OmegaK(r,Omegaout,Omegain,r0,Deltar)/OmegaM(r,Omegaout,  
Omegain,r0,Deltar))

All other expressions are just algebraic manipulation of this

---

<div class="post-metadata">

**Author:** ![chI-No](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/chi-no/32/12632_2.png) [@chI-No](https://discourse.mc-stan.org/u/chI-No)\
**Post date:** [September 12, 2021, 2:55am UTC](https://discourse.mc-stan.org/t/link-symbolicc-function-to-stan-interface/24216/10 "2021-09-12T02:55:23Z")

</div>

Just to close this issue, i found that the best solution is to perform the deruvative in mathematica and then write this in c languaje via Cform[] command. It will gives you the function and you can add it to the stan code exporting to a file “func.txt” and then just importing and pasting this text in the stan code like:

“… return” \<\> file \<\> ;  
}

Otherfunction(…"

Hope it helps to someone