creating user defined (normal mean) function in C++

[Jarrett Barber]

Hi All,

I have a simulator (of tree growth) written in C++ whose returned value I'd like to use as the mean of a Gaussian specification in Stan. The simulator has about 30 parameters plus various other inputs. All parameters spaces are subsets of the reals. I'm hoping that I can simply modify mymodel.cpp, generated by Stan, in the (hopefully) obvious place(s), then link/compile as usual, letting Stan do the rest, e.g., gradients, etc.

Am I off track or can you provide some suggestions on how to proceed?

Best,

Jarrett


Jarrett Barber
School of Mathematical and Statistical Sciences
Arizona State University
P.O. Box 871804
Tempe, AZ 85287-1804
Jarrett dot Barber at asu dot edu

Save Net Neutrality!

[Bob Carpenter]

On 1/15/14, 10:55 PM, Jarrett Barber wrote:
> I have a simulator (of tree growth) written in C++ whose returned value I'd like
> to use as the mean of a Gaussian specification in Stan. The simulator has
> about 30 parameters plus various other inputs. All parameters spaces are subsets
> of the reals. I'm hoping that I can simply modify mymodel.cpp, generated by Stan,
> in the (hopefully) obvious place(s), then link/compile as usual, letting Stan do
> the rest, e.g., gradients, etc.

I think Daniel's writing some general instructions for people who want to do this.

You have to:

1. Make sure you declare the parameters with appropriate constraints. Stan assumes
the support of the log probability function matches the declared constraints on
parameters. Then you leave the parameter transforms in place.

2. In general, you have to code a templated version of your function, with different
template parameters for each argument. You can get by leaving the data or transformed
data variables as doubles. You can test independently of the model by making sure
you can instantiate the arguments that will be parameters, transformed parameters, or
local variables with stan::agrad::var. This will also let you test the gradients, and
it's a good idea to test vs. a finite-difference implementation.

Then you have to call our built-in functions rather than the built-in C++ library functions
to make sure auto-diff will work. But the compiler will let you know where you need to
do that.

If you send us the signature of your function, I can show you what the signature needs
to look like.

An alternative is to add your new function to Stan directly. This is cleaner,
but we don't yet have a way to make it easy to extend, so you have to modify our
source and compile with that. There are some sketchy instructions on how to do this
here:

https://github.com/stan-dev/stan/wiki/Contributing-New-Functions-to-Stan

- Bob

[Jarrett Barber]

Hi Bob,

Here's the function prototype.  Only parameters in the struct sparms are to be sampled, all other arguments are not to be sampled.  The struct sparms prototype follows, in case that's needed.  Thanks for the help.  -- Jarrett

void growthloop(sparms *p, gparms *gp, double *r0, double r[], double h[], double rBH[], int *t);

typedef struct{

  double hmax; ///< Max. tree height

  double phih; ///< Slope at H vs. r curve at r=0

  double eta; ///< Relative crown height

  double etaB; ///< ???

  double swmax; ///< Maximum sapwood width

  double lamdas; ///< ???

  double lamdah; ///< ???

  double rhomax; ///< ???

  double rhomin; ///< ???

  double f2; ///< Leaf area to xylem conducting area ratio

  double f1; ///< Fine root area to leaf area ratio

  double gammac; ///< Max. storage capacity of sapwood cellsx

  double gammaw; ///< (Inverse) density of sapwood structural tissue

  double gammax; ///< Xylem conducting area to sapwood area ratio

  double cgl; ///< Construction costs of producing leaves

  double cgr; ///< Construction costs of producing fine roots

  double cgw; ///< Construction costs of producing sapwood

  double deltal; ///< Labile carbon storage capacity of leaves

  double deltar; ///< Labile carbon storage capacity of fine roots

  double sl; ///< Senescence rate of leaves

  double sla; ///< Specific leaf area

  double sr; ///< Senescence rate of fine roots

  double so; ///< Senescence rate of course roots and branches

  double rr; ///< ???

  double rhor; ///< Tissue density of fine roots

  double rml; ///< Maintenance respiration rate of leaves

  double rms; ///< Maintenance respiration rate of sapwood

  double rmr; ///< Maintenance respiration rate of fine roots

  double drcrit; ///< ??? may vary by tree

  double drinit; ///< radial increment (may vary by tree??)

  double K; ///< ???

  double epsg; ///< ???

  double M; ///< ???

  double alpha; ///< ???

  double R0; ///< ???

  double R40; ///< ???

} sparms;

[Bob Carpenter]

On 1/20/14, 9:46 PM, Jarrett Barber wrote:
> Hi Bob,
>
> Here's the function prototype. Only parameters in the struct sparms are to be sampled, all other arguments are not to
> be sampled. The struct sparms prototype follows, in case that's needed. Thanks for the help. -- Jarrett

> void growthloop(sparms *p, gparms *gp, double *r0, double r[], double h[], double rBH[], int *t);

Which of these are inputs and which are outputs? And are the double* and int* arrays
or just single values? And is the gparms struct just another bunch of double values?

This old-school C is very underspecified at the function level!

There is no notion of struct in Stan, so you'd have to pass in std::vector containers
with the values for sparms elements and for gparms elements and then create
the structs internally.

Then, sparm is going to have to be generalized so that all of the double
values are stan::agrad::var values. And then any function that operates
on them is going to have to be generalized so that it deals with var
values appropriately. We probably have most of the functions you need,
but if your C++ function is complicated, it's going to be a lot of work.

- Bob

> <https://github.com/stan-dev/stan/wiki/Contributing-New-Functions-to-Stan>
>
> - 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.
> For more options, visit https://groups.google.com/groups/opt_out.

[Jarrett Barber]

Hi Bob,

sparms, gparms, double *r0 are inputs.  sparms contains parameters, gparms and r0 are data.  r0 is a single value.  Everything else is output, including int *t which is just a single value.

It may be better for me to change the structures to std::vec to be consistent with Stan, rather than create structures every time the function is called. ??

By "generalize to stan::agrad::var values" you are referring to templates? (I _am_ old school C, so any examples in stan generated C++ that you could point to would be helpful.)

I suspect the time I spend updating my function (which _is_ complicated) will pay off in the future with much shorter chains to convergence using Stan.  (I plan to use/extend the function for years to come.)

Thanks again.

Jarrett

[Bob Carpenter]

You can look at any of the built-in matrix functions in
stan/agrad/rev/matrix, though it'd probably be easiest to
start with stan/agrad/rev/functions, which just take
single inputs.

You can either template out the arguments or write them
for stan::agrad::var. Either way, the inputs should be
constant references, either "const T&" or "const stan::agrad::var&".

The Stan language itself doesn't have any way of
taking mutable arguments that are set, so everything
has to be packed into a vector return.

Eventually, we will probably add R-style-lists/C-style-structs,
but it's not on our short-term to-do list (meaning it's
unlikely to be done in 2014).

- Bob

> stan-users+...@googlegroups.com <javascript:>.
> > For more options, visit https://groups.google.com/groups/opt_out <https://groups.google.com/groups/opt_out>.

[Jarrett Barber]

To check my understanding so far:

In short: template my functions and template my parameter arguments or, alternatively,  “classify” my parameters according to the stan::agrad::var hierarchy.  Other, non-parameter arguments (function data inputs or function outputs, even by reference) do not need to be generalized, correct? Parameters, other inputs and outputs should be done via std::vector that Stan creates when translating from Stan to C++, correct?

Without a firmer grasp of these details, I would prefer to be as Stan-independent as possible to allow for running of my growth model outside of Stan, so this makes me want to choose _not_ to use the stan::agrad::var hierarchy but to use instead the template mechanism.

Does this make sense?

Best,

Jarrett

[Michael Betancourt]

Here’s an example.  Consider the function

double f(double x1, double x2, double y1, double y2) {

  double x_sum = x1 * x1 + x2 * x2;

  double y_sum = y1 * y2;

  return x_sum / y_sum;

}

where x1, x2 are parameters and y1, y2 are data.  You want to template the function as

template <typename T>

T f(T x1, T x2, double y1, double y2) {

  T x_sum = x1 * x1 + x2 * x2;

  double y_sum = y1 * y2;

  return x_sum / y_sum;

}

Because *, +, and / are all defined for our autodiff variables you can then call

double out = f<double>(x1, x2, y1, y2) [ calls regular version of code ]

stan::a grad::var out = f<stan::agrad::var>(x1, x2, y1, y2) [ calls autodiff version ]

You just have to be careful that you don’t use any functions or operations in the definition f

that aren’t defined for out autodiff variables.

[Jarrett Barber]

This helps!  Thanks.

Alas, my tree growth function calls other user-written C functions, which call other user-written C functions…most of which make heavy use of gsl C library functions.  Instead of templafying “everything”, I am hoping to create a sort of wrapper for my “main” tree growth function, so that the wrapper assigns the values of stan::agrad::var typed parameters to standard types that subsequent functions can operate on without (much) further generalization.  Does this sound reasonably doable to you?

Best,

Jarrett

[Michael Betancourt]

Alas, my tree growth function calls other user-written C functions, which call other user-written C functions…most of which make heavy use of gsl C library functions.  Instead of templafying “everything”, I am hoping to create a sort of wrapper for my “main” tree growth function, so that the wrapper assigns the values of stan::agrad::var typed parameters to standard types that subsequent functions can operate on without (much) further generalization.  Does this sound reasonably doable to you?

Not particularly.  It is possible but would require adding a new var implementation which requires a much deeper understanding of the inner workings of the autodiff code.  If you want to proceed then take a look at how other functions have been implemented (you have to define a var and then implement the chainable method in the corresponding vari).

[Bob Carpenter]

You either have to be able to template everything
else out or you have to be able to calculate
derivatives with respect to the parameters. If you
can calculate derivatives, it's actually not that hard
to interface with the derivative code. But presumably
that won't be easy if it's a nest of C calls.

Stan's not really defined to let you plug in arbitrary
functions. Some other packages, such as emcee in Python
(which I haven't used) are better positioned to let you plug
in arbitrary probability functions.

- Bob

On 1/21/14, 11:04 PM, Michael Betancourt wrote:
>> Alas, my tree growth function calls other user-written C functions, which call other user-written C functions�most of
>> which make heavy use of gsl C library functions. Instead of templafying �everything�, I am hoping to create a sort of
>> wrapper for my �main� tree growth function, so that the wrapper assigns the values of stan::agrad::var typed

>> parameters to standard types that subsequent functions can operate on without (much) further generalization. Does
>> this sound reasonably doable to you?
>>
> Not particularly. It is possible but would require adding a new var implementation which requires a much deeper
> understanding of the inner workings of the autodiff code. If you want to proceed then take a look at how other
> functions have been implemented (you have to define a var and then implement the chainable method in the corresponding
> vari).
>

[Jarrett Barber]

After looking at some of the function definitions in the C++ (Stan) source code, and recalling the introductory discussion in the Stan Manual about analytic derivatives, I now understand the main reason for templating (analytic derivatives).  It seems silly now that I could have somehow (implicitly) considered that numerical differentiation could be used and still result in the speediness of Stan.

It does appear that the Stan C++ library does provide all of the functions that I will need to revise my (my colleague's really) tree growth model code.  Encouraging.  But, after my realization of analytic derivatives, I began to think of the chain rule, which, in my case, quickly reminded me that one of the output quantities in my growth model is defined implicitly as the result of a root-finding routine: find the tree radius such that an expression of tree size (allometry) equals an expression of carbon allocation, else risk getting a grossly misproportioned and dysfunctional tree.  Of course, these size and allocation expressions are functions of parameters to be sampled.  (I intend that the tree radius provided by my model to be essentially the mean of a normal distribution specified in Stan language, hence the connection to gradients of the log-posterior.)  More precisely, I have functions f(radius, parameters) and g(radius, parameters), and radius is defined implicitly as the solution to f(radius, parameters) = g(radius, parameters), and f and g are otherwise fairly simple.

Differentiation (wrt parameters) for my implicitly defined tree radius is no problem on paper (Implicit Function Theorem) (I haven't actually done this yet, but I expect this to be no problem, on paper, because my equations are relatively simple.)  But, my lack of experience in C++/Stan innards leaves me wondering: will I be able to exploit a relatively high-level solution using templates, as we have been discussing, or do I have to dig down and "calculate derivatives" and "interface with the derivative code" as Bob suggests?

Anyway, responses are appreciated.  I hope someone finds this discussion at least mildly interesting and useful.

BTW, have been using my own MH MCMC routine for sampling a small subset of parameters, but I anticipate problems in the future as I try to sample a larger number of parameters, use larger data sets, and make model extensions.  I think HMC might be one way to help with anticipated problems, and I think Stan might save me from reinventing the wheel. 

Best,

Jarrett

On Tuesday, January 21, 2014 5:38:23 PM UTC-7, Bob Carpenter wrote:

You either have to be able to template everything
else out or you have to be able to calculate
derivatives with respect to the parameters.  If you
can calculate derivatives, it's actually not that hard
to interface with the derivative code.  But p

[Michael Betancourt]

Differentiation (wrt parameters) for my implicitly defined tree radius is no problem on paper (Implicit Function Theorem) (I haven't actually done this yet, but I expect this to be no problem, on paper, because my equations are relatively simple.)  But, my lack of experience in C++/Stan innards leaves me wondering: will I be able to exploit a relatively high-level solution using templates, as we have been discussing, or do I have to dig down and "calculate derivatives" and "interface with the derivative code" as Bob suggests?

It is possible to automatically differentiate these kinds of solution, by either explicitly differentiating the root finder or appealing to the things like the Implicit Function Theorem.  We have been skirting with these kind of approaches as we work on introducing differential equations into Stan, but nothing has been implemented yet.  So for the moment you would either have to hack the code or incorporate the particular system into a proper Stan function.

 

BTW, have been using my own MH MCMC routine for sampling a small subset of parameters, but I anticipate problems in the future as I try to sample a larger number of parameters, use larger data sets, and make model extensions.  I think HMC might be one way to help with anticipated problems, and I think Stan might save me from reinventing the wheel. 

HMC will make a huge difference going forward, so any efforts to utilize it now will pay off in spades down the line. 

[Bob Carpenter]

You could also use finite differences, though it's likely
to be slow, and also can have problems with accuracy. But
it might work.

Maybe we should think about doing that in general for complex
functions of not-too-many parameters?

- Bob

[Jarrett Barber]

Continued thanks to Bob and Michael for kindly guiding me toward to a solution.

I’m trying to reconcile Michael’s comments, “It is possible to automatically differentiate these kinds of solutions, by either explicitly differentiating the root finder or appealing to the things like the Implicit Function Theorem” with, in an earlier post, “It is possible but would require adding a new var implementation which requires a much deeper understanding of the inner workings of the autodiff code.”  (Perhaps I mistakenly conflate these 2 posts.)

Selfishly assuming that I’m only looking for a solution to my particular problem at the moment, I may use implicit differentiation to get the derivative of my implicit function.  A more careful (considering existence, multiple implicit functions, etc.) approach for more a more general solution of cases with implicit functions may have to appeal to the Implicit Function Theorem (or related), along with a much more ambitious code effort.  But I think I may get by for now with the derivative obtained by implicit differentiation for my particular case.

So, I have (“…explicitly…”, as I take Michael to mean) the analytic derivative in hand.  Then, I must implement “…via a new var implementation…then implement the chainable method in the corresponding vari”.  Once I do this, then, as Michael says, I can “…automatically differentiate…”.  Yes?

Am I missing a possibly more “automatic”, high-level approach that doesn’t require entering the innards of var and vari and the chainable method?

And, yes, I do anticipate huge differences going forward (and Stan's cool), hence my persistence in spite of revealing my neophyte level of C++ know-how. (I can see the power of generalization, but I’m still devising my coup.)

I’ll keep in mind Bob’s suggestion about finite differences.

Somewhat gropingly,

Jarrett

[Bob Carpenter]

You should definitely listen to Michael on this.
He has a much deeper understanding of the calc and
of HMC than I do!

- Bob

[Michael Betancourt]

So, I have (“…explicitly…”, as I take Michael to mean) the analytic derivative in hand.  Then, I must implement “…via a new var implementation…then implement the chainable method in the corresponding vari”.  Once I do this, then, as Michael says, I can “…automatically differentiate…”.  Yes?

Correct, not that it’s an easy task — it took me a while to get my head around our implementation in order to do things like this myself.  We can help you through it if you want to push forward.

Autodiff really just computes the chain rule automatically — it has no idea how to compute partial derivatives, so what our autodiff library does is define partial derivatives for most of the useful functions out there.  If you want to add a new function then you have to tell Stan its partials.