Using external software inside Stan

[Ben B]

I think this question has been asked a lot, but I wanted to clarify.

I've got a small custom HMC code that does Bayesian inference for a structural mechanics problem.

I'd like to do inferences on it in Stan so I can take advantage of the Stan samplers and use the Stan modeling language for the priors and stuff. However:

1. The problem, at the end of the day, reduces to a big eigenvalue solve. I need control over the eigenvalue solve beyond what Stan currently offers (generalized eigenvalue solve, only need the few smallest eigenvalues)

2. The actual mechanics are kinda black box too. I could conceivably code up some of the matrix assembly stuff in Stan, but it isn't terribly practical, and in the near future I might be needing some other software to do the assembly for me anyway.

Either way, I can provide my stuff as a black box. So it takes in parameters, and spits out some stuff that corresponds roughly to data I'm measuring. I can also provide the gradients for this as well.

So if I want to implement the 'mechanics_solve' thing in something like:

data ~ normal(mechanics_solve(parameters), sigma);

And I have my forward function evaluation as well as my gradients, where's the best place for me to start? I use PyStan, if that changes anything.

From what I've dug up on the mailing list, I go here: https://github.com/stan-dev/stan/wiki/Contributing-New-Functions-to-Stan, find a function that looks reasonably like what I want in the Stan Math library, and plug away?

Ben

[Bob Carpenter]

If you can write up a C++ function that does what you want,
then we can help you plug it together. If you have analytic
gradients, you can just plug them in using:

stan_math/stan/math/rev/core/stored_gradient_vari.hpp

That is, you want your function to be something like this:

std::vector<var> mechanics_solve(const std::vector<var>& v) {
copy args into double structure
compute values and gradients
use stored_gradient_vari to create return vars
return them
}

You need to use our built in allocation memalloc_.allocate_array<T>(size)
to generate the pointer arrays to pass into the stored gradients so
there are no memory leaks.

I should write this up in the how to write functions.

Then just add it to symbol table and you're good to go (assuming you
can build it into PyStan or use it from CmdStan).

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Andrew Gelman]

Just to follow up on Bob's comment: If you (Ben B) do this for a real application, it would be great to see this written up. I think that if people can see a successful application of this sort by a Stan "outsider," then it could open the floodgates--in a good way!
Andrew

[Ben B]

@Bob/@Andrew, Thanks for the enthusiasm!

I'll get a simple example coded up. It shouldn't take me too long. I'll follow up in this thread when I've got something concrete. What Bob said sounds exactly like what I'm wanting to do.

Ben

>> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

>> To post to this group, send email to stan-...@googlegroups.com.
>> For more options, visit https://groups.google.com/d/optout.
>
> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.


--

You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/PQICLhcmOlg/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Ben Goodrich]

On Monday, March 13, 2017 at 6:32:36 PM UTC-4, Ben B wrote:

1. The problem, at the end of the day, reduces to a big eigenvalue solve. I need control over the eigenvalue solve beyond what Stan currently offers (generalized eigenvalue solve, only need the few smallest eigenvalues)

Stan cannot autodiff through complex generalized eigenvalues, which is why we currently only support eigenvalues of symmetric matrices. Even if your small eigenvalues happen to be real, if complex arithmetic is used, it will probably segfault. I have some stuff on a branch that mitigates some of these problems using what Eigen calls a pseudo-eigendecomposition that only uses real arithmetic but it doesn't do anything partial like finding the few smallest generalized pseudo-eigenvalues.

Ben

[Michael Betancourt]

+1

I’ll also note that statistical models built out of inversions/decompositions

are almost always a mess to fit.  This can be explained by their compromising

the generative structure of the model, or just looking at how ill-posed the

gradients are.

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

[Ben B]

@BGoodrich

Yeah, the eigenvalues are actually real in this problem. I'm pretty sure something in ARPACK will solve my problem (I think that's what Scipy uses under the hood anyway?) -- it's just a matter of knocking the dust off my documentation reading skills.

@Michael

>> I’ll also note that statistical models built out of inversions/decompositions are almost always a mess to fit.

Haha, agreed in my limited experience. I think we got pretty lucky with this model/data. I coded this up in Python with the 2012 Radford Neal HMC paper and things just kinda work, which is why I think there's a chance I could take advantage of all the nice Stan stuff if I can weasel it in. You're right though; it still feels like a coinflip.

Ben

[Ben Goodrich]

On Tuesday, March 14, 2017 at 12:17:39 PM UTC-4, Ben B wrote:

@BGoodrich

Yeah, the eigenvalues are actually real in this problem. I'm pretty sure something in ARPACK will solve my problem (I think that's what Scipy uses under the hood anyway?) -- it's just a matter of knocking the dust off my documentation reading skills.

The eigenvalues not only need to be real; the C++ needs to use real arithmetic exclusively to calculate them. If it is using complex<T> then you are going to have pull the double values out of the stan::math::var types, calculate the small eigenvalues, do the derivatives analytically, and then attach the derivatives to new stan::math::var values. And that is a lot of work.

[Bob Carpenter]

If Ben B calculates the derivatives himself, then it's
no holds barred on both the values and derivatives.

- Bob

[Ben B]

Ah! Yes, it's what Bob is saying. I do the derivatives myself, thanks.

Ben

On Tue, Mar 14, 2017 at 10:35 AM, Bob Carpenter <ca...@alias-i.com> wrote:

If Ben B calculates the derivatives himself, then it's
no holds barred on both the values and derivatives.

- Bob

> On Mar 14, 2017, at 1:18 PM, Ben Goodrich <goodri...@gmail.com> wrote:
>
> On Tuesday, March 14, 2017 at 12:17:39 PM UTC-4, Ben B wrote:
> @BGoodrich
>
> Yeah, the eigenvalues are actually real in this problem. I'm pretty sure something in ARPACK will solve my problem (I think that's what Scipy uses under the hood anyway?) -- it's just a matter of knocking the dust off my documentation reading skills.
>
> The eigenvalues not only need to be real; the C++ needs to use real arithmetic exclusively to calculate them. If it is using complex<T> then you are going to have pull the double values out of the stan::math::var types, calculate the small eigenvalues, do the derivatives analytically, and then attach the derivatives to new stan::math::var values. And that is a lot of work.
>
> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

--

You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/PQICLhcmOlg/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Michael Betancourt]

>> I’ll also note that statistical models built out of inversions/decompositions are almost always a mess to fit.

Haha, agreed in my limited experience. I think we got pretty lucky with this model/data. I coded this up in Python with the 2012 Radford Neal HMC paper and things just kinda work, which is why I think there's a chance I could take advantage of all the nice Stan stuff if I can weasel it in. You're right though; it still feels like a coinflip.

You will likely find that “just kinda works” is really “works for

a small part of parameter space, preventing the sampler from

exploring all of the relevant parts of parameters space and

hence biasing inferences.”  Stan will certainly help you 

diagnose such problems, but it won’t be able to fix them.

[Ben B]

@Bob

I got around to implementing a little C++ version of the mechanics code.

The function signature I settled on was:

void mechanics(double c11, double anisotropic, double c44, // These are the parameters I want to estimate for now

               Eigen::Tensor<double, 4> *dp, Eigen::Tensor<double, 2> *pv, unsigned int nevs, // These are constant bits of data I can initialize and pass into Stan in as vectors or whatever

               MatrixXd **freqs,  // This is the output -- an (nevs x 1) matrix of eigenvalues

               MatrixXd **dfreqs_dc11 = NULL, MatrixXd **dfreqs_da = NULL, MatrixXd **dfreqs_dc44 = NULL) // This is output too -- these are the derivatives of the eigenvalues wrt the parameters

What's the closest thing in the Stan Math library to this (for me to start working on)? I can rearrange things so the input is just a vector and output is just a vector (like you showed in your original e-mail), but for testing I thought this was cleaner. What do I need to do about the constants (the things I won't need to differentiate with respect to)?

I attached a little PDF with some details about the problem as well as the rough Stan model I'd like to write. If the attachment doesn't work, it's available here: https://www.overleaf.com/read/qpcmfzmxfxsf

Basically all the code that I'd need hook up to the Stan Math library is here: https://github.com/bbbales2/modal_cpp/blob/master/mechanics.hpp and https://github.com/bbbales2/modal_cpp/blob/master/polybasis.hpp . The only super necessary external dependency for now is LAPACK.

There's a really simple sampler here (from the 2012 Radford Neal paper) that I tested things with: https://github.com/bbbales2/modal_cpp/blob/master/main.cpp . I've got traceplots from that, if anyone is curious.

@Michael, yeah, we're all pretty curious over here to see what NUTS thinks of the problem. It would be nice to not have to guess and check with timesteps a bunch to get the the thing sampling efficiently too.

Thanks a lot!

Ben

--

You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/PQICLhcmOlg/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Bob Carpenter]

This should all be OK if you're writing a function for
your own use and can deal with the LAPACK linkage. Our
existing functions tend to be much more C++ oriented, with
constant reference arguments and return values; under the hood
we sometimes work with raw pointers when we need to control our
own memory.

Basically, you need to be able to construct the vari objects
that correspond to your results with chain() methods that do
the reverse mode autodiff (multiply the function's result's
adjoints by the partials and add to the corresponding operand's
adjoints). This is all explained in our arXiv paper.

If you just want an example of multivariate functions, all
of the matrix functions are like this, but their code's pretty
involved.

For constants (double arguments, presumably), you don't need
to do anything. You won't be able to run Stan programs
that deal with them as parameters, but you shouldn't need to.
We don't have built-in functions like that.

- Bob

> To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
>
>

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

> <resonance-modes.pdf>

[Ben B]

Damn Bob, I got this e-mail at 9:19PM PST. Aren't you guys EST? Get some sleep my dude

>> This should all be OK if you're writing a function for your own use

Yeah, hopefully I can manage a math library fork. It's either I do that or I maintain all my own samplers and stats stuff. I honestly dunno which'll be harder at the end of the day.

>> This is all explained in our arXiv paper.

Oh dang I better read this then. Looks like it's got what I need.

I'll letchall know if I can this working.

Thanks!

Ben

> To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
>
>
> --

> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
> <resonance-modes.pdf>

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/PQICLhcmOlg/unsubscribe.

To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Ben Goodrich]

On Wednesday, March 29, 2017 at 12:59:21 PM UTC-4, Ben B wrote:

Yeah, hopefully I can manage a math library fork. It's either I do that or I maintain all my own samplers and stats stuff. I honestly dunno which'll be harder at the end of the day.

Hopefully, you do not need a math library fork. If you are using R, see

https://cran.r-project.org/web/packages/rstan/vignettes/external.html

or if you are using CmdStan, see its manual for (slightly different) instructions on how to inject C++ into a Stan program.

Ben

[Ben B]

@BGoodrich

That's pretty neat, thanks! I think I still need to inject into the Math library cause I'm messing around with custom gradients, but there's probably some other things I could use that for.

I read the Arxiv Stan Math paper and I think what I'm doing looks a lot like the log determinant example there. I'm gonna give that a go and see what happens.

Ben

[Ben B]

@Bob/BGoodrich

I made a Stan function for my stuff and tested the autodiff on it. It was a lot easier than I expected haha! Came out looking pretty cool in the end (https://github.com/bbbales2/modal_cpp/blob/master/stan_test.cpp (lines 30-62)). The stored_gradient_varis was exactly what I needed.

Now to use this from a Stan Model I'm guessing I just follow the directions here: https://github.com/stan-dev/stan/wiki/Contributing-New-Functions-to-Stan ?

One more question though, one of my input arguments (dp, of type dpT) is of type Matrix< Matrix<double, 3, 3>, Dynamic, Dynamic >. I was hoping that'd map to a type in Stan like "matrix[N, N] y[3, 3]." Any chance of that? I didn't see an obvious thing in https://github.com/stan-dev/stan/blob/develop/src/stan/lang/function_signatures.h. Or do I need to go ahead and make it a regular Matrix? It isn't really a problem. Just cleaner to have it as a 4d sorta thing.

Also @BGoodrich, I had a look at the cmdStan manual. Is what you were saying that I might be able to just take the function I wrote (that returns a vector of vars with stored_gradient_varis preset) and include it in a Stan model like on page 25/26 with the Bernoulli example? I think I may have misunderstood the first time.

Ben

[Ben Goodrich]

On Thursday, March 30, 2017 at 10:25:30 PM UTC-4, Ben B wrote:

Also @BGoodrich, I had a look at the cmdStan manual. Is what you were saying that I might be able to just take the function I wrote (that returns a vector of vars with stored_gradient_varis preset) and include it in a Stan model like on page 25/26 with the Bernoulli example? I think I may have misunderstood the first time.

Yes

[Bob Carpenter]

> On Mar 30, 2017, at 10:25 PM, Ben B <bbba...@gmail.com> wrote:
>
> @Bob/BGoodrich
>
> I made a Stan function for my stuff and tested the autodiff on it. It was a lot easier than I expected haha! Came out looking pretty cool in the end (https://github.com/bbbales2/modal_cpp/blob/master/stan_test.cpp (lines 30-62)). The stored_gradient_varis was exactly what I needed.

Glad to hear it. The new operands_and_partials Sean is coding
up will be the way to do this in the future with more natural types.
And we hope to add better doc.

> Now to use this from a Stan Model I'm guessing I just follow the directions here: https://github.com/stan-dev/stan/wiki/Contributing-New-Functions-to-Stan ?

Yes.

> One more question though, one of my input arguments (dp, of type dpT) is of type Matrix< Matrix<double, 3, 3>, Dynamic, Dynamic >. I was hoping that'd map to a type in Stan like "matrix[N, N] y[3, 3]." Any chance of that?

That Stan type has the same shape, but it'll produce this

std::vector<std::vector<Matrix<double, -1, -1> > >

on the inside, so you'd have to write a wrapper for your function with
this type and convert it.

> I didn't see an obvious thing in https://github.com/stan-dev/stan/blob/develop/src/stan/lang/function_signatures.h. Or do I need to go ahead and make it a regular Matrix? It isn't really a problem. Just cleaner to have it as a 4d sorta thing.
>
> Also @BGoodrich, I had a look at the cmdStan manual. Is what you were saying that I might be able to just take the function I wrote (that returns a vector of vars with stored_gradient_varis preset) and include it in a Stan model like on page 25/26 with the Bernoulli example? I think I may have misunderstood the first time.

You can also include it directly in a Stan model.

Ben was also saying that you don't need to fork per se so much as
get your own function into the compiel path. It will have to
be in the stan::math namespace for everything to work.

>
> Ben
>
> On Wed, Mar 29, 2017 at 8:41 PM, Ben B <bbba...@gmail.com> wrote:
> @BGoodrich
>
> That's pretty neat, thanks! I think I still need to inject into the Math library cause I'm messing around with custom gradients, but there's probably some other things I could use that for.
>
> I read the Arxiv Stan Math paper and I think what I'm doing looks a lot like the log determinant example there. I'm gonna give that a go and see what happens.
>
> Ben
>
> On Wed, Mar 29, 2017 at 10:14 AM, Ben Goodrich <goodri...@gmail.com> wrote:
> On Wednesday, March 29, 2017 at 12:59:21 PM UTC-4, Ben B wrote:
> Yeah, hopefully I can manage a math library fork. It's either I do that or I maintain all my own samplers and stats stuff. I honestly dunno which'll be harder at the end of the day.
>
> Hopefully, you do not need a math library fork. If you are using R, see
>
> https://cran.r-project.org/web/packages/rstan/vignettes/external.html
>
> or if you are using CmdStan, see its manual for (slightly different) instructions on how to inject C++ into a Stan program.
>
> Ben
>
>
>
>
>
> --
> You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/PQICLhcmOlg/unsubscribe.

> To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
>
>
>

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

[Ben Goodrich]

On Friday, March 31, 2017 at 2:39:49 PM UTC-4, Bob Carpenter wrote:

Ben was also saying that you don't need to fork per se so much as
get your own function into the compiel path.  It will have to
be in the stan::math namespace for everything to work.

For CmdStan, yes. For RStan, it is included in the model's namespace.

[Bob Carpenter]

Neat---that way you don't need the general argument-dependent
lookup we usually rely on.

- Bob

[Andrew Gelman]

It would be great for this all to be a case study, blog, etc., at some point.
A

[Ben Goodrich]

On Friday, March 31, 2017 at 5:32:50 PM UTC-4, Andrew Gelman wrote:

It would be great for this all to be a case study, blog, etc., at some point.

Vignette

[Andrew Gelman]

http://andrewgelman.com/2017/03/31/running-stan-external-c-code/

[Ben B]

@Bob/Ben

Groovy. I'll keep you posted. Thanks for the help.

@Andrew

Yeah, sure thing! I'll make a little write up once all the pieces are together and I get a sense for how the sampler is behaving. I'm with Michael in that I'm not ready to pull out the party hats till I have a big pile of healthy samples in my hands :D.

Ben

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/PQICLhcmOlg/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Ben B]

It's party hat time!

Got a really simple exampled coded up (using the cmdStan stuff BGoodrich recommended). Here's the direct comparison from what I had to Stan (model here: https://github.com/bbbales2/modal_cpp/blob/master/test.stan):

Parameter     HMC (100 samples)  Stan (100 samples), posteriors shown as mean +- 1 standard deviation

c11           1.7 +- 0.019       1.7 +- 0.017

a             1.0 +- 0.0024      1.0 +- 0.0024

c44           0.45 +- 0.00078    0.45 +- 0.0007

sigma         0.43 +- 0.062      0.44 +- 0.069

I told Gmail to use a fixed space font so hopefully that worked.

Not exactly a whole lot of samples (my forward model is pretty slow [200ms/eval] and I wanted to see results), and I made the problem really easy by giving Stan good guesses (the HMC had pretty good guesses too), but woohoo anyway! I'll need to add a few things to put the full mechanics model in and test some harder problems (and see if I can get away without the guesses) but it looks like we have liftoff!

To get everything working, I needed to define a template specialization of my mechanics function that took in vars for parameters as well as one that took doubles. The linker wanted a version of each (code here: https://github.com/bbbales2/modal_cpp/blob/master/stan_mech.hpp).

Does NUTS use the 2nd version (the non-var specialization)? I added some debug print statements and it looks like the variational inference code uses it.

I'll follow up again after I've done more tests and attach a write up.

Thanks folks! This is gonna be super useful.

Ben

[Andrew Gelman]

I'm looking forward to seeing this on the blog!

You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

[Bob Carpenter]

Glad to see the results match.

I'm pretty sure we only need to instantiate a double-only
version of the transformed parameters block and the generated
quantities, not the model block.

- Bob

> To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
>
>
>

> --

> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.