programming user-defined functions (C++, or another language)

[Tamas Papp]

Hi,

1. I need to implement a function that is very complex, and uses a
numerical, iterative algorithm. My understanding is that I need to
program a C++ function for this, that implements the function and its
Jacobian, compile and link it to Stan. However, this is just a broad
outline, and I am wondering if there is a tutorial or an example.

2. Technically, if the function is f:R^N -> R^M, and I can calculate
f(x) and the Jacobian at x, for x which is a vector of doubles, is that
enough?

3. Ideally, I would like to implement the function not in C++, but
Julia, which has a C interface, which I believe that I can somehow patch
into a C++ interface that Stan can call. Is there an a priori problem
with this? I know this is more hassle and I will need to get my hands
dirty, but I am more comfortable with Julia than C++. I will make the
code available but it is for a paper, so I don't aim to integrate it
into Stan (I am just mentioning this to make it clear that the toolchain
nightmare will be mine, not anyone else's).

Best,

Tamas

[Allen B. Riddell]

On 09/23, Tamas Papp wrote:
> 3. Ideally, I would like to implement the function not in C++, but
> Julia, which has a C interface, which I believe that I can somehow patch
> into a C++ interface that Stan can call. Is there an a priori problem
> with this? I know this is more hassle and I will need to get my hands
> dirty, but I am more comfortable with Julia than C++. I will make the
> code available but it is for a paper, so I don't aim to integrate it
> into Stan (I am just mentioning this to make it clear that the toolchain
> nightmare will be mine, not anyone else's).

Embedding interpreters in C++ is something that people do. Here's a blog
post about doing it with Python
http://www.codeproject.com/Articles/11805/Embedding-Python-in-C-C-Part-I,
and I noticed this thread on the Julia list "Embed julia in C/C++"
http://julia-programming-language.2336112.n4.nabble.com/Embed-julia-in-C-C-td3339.html

Others can weigh in but I suspect your best bet would be to fork the
Stan codebase and add your function (writing it in C++ strikes me as
easier than dealing with the embedding problem).

Best,

Allen

[Michael Betancourt]

Most importantly such embeddings will break the
autodiff, so you’d have to be _really_ careful to
isolate them within a var implementation. This
is possible, but not easy.

> --
> 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.

[Daniel Lee]

Just wanted to respond to #1:

On Sep 23, 2015 11:54 AM, "Tamas Papp" <tkp...@gmail.com> wrote:
> 1. I need to implement a function that is very complex, and uses a
> numerical, iterative algorithm. My understanding is that I need to
> program a C++ function for this, that implements the function and its
> Jacobian, compile and link it to Stan. However, this is just a broad
> outline, and I am wondering if there is a tutorial or an example.

You can write the function in Stan. If there's an iterative algorithm, that doesn't necessarily pose a problem. The autodiff will work through the iterative algorithm. (One of the benefits of autodiff compared to symbolic differentiation.)

There might be some speed gains if you implemented it in c++ directly, especially if you have the gradients with respect to all the parameters of the model, but if it were me, I'd start with a full Stan implementation.

[Tamas Papp]

Is there an example somewhere?

Best,

Tamas

[Tamas Papp]

The algorithm is very complicated, and not ideally suited for autodiff:
I find y=f(x) iteratively, but given y there is an (almost) closed-form
for the Jacobian.

Best,

Tamas

[Daniel Lee]

I'm not quite following what you mean by the Jacobian.

For Stan, we need the gradients with respect to all parameters of the Stan program (not necessarily x). So, we really need d f(x) / d theta for all thetas of the model.

Do you have those in closed form?

Daniel

[Ben Goodrich]

On Wednesday, September 23, 2015 at 11:54:33 AM UTC-4, Tamas Papp wrote:

1. I need to implement a function that is very complex, and uses a
numerical, iterative algorithm. My understanding is that I need to
program a C++ function for this, that implements the function and its
Jacobian, compile and link it to Stan. However, this is just a broad
outline, and I am wondering if there is a tutorial or an example.

The Wiki is somewhat out of date. A similar question came up on

http://stackoverflow.com/questions/31862831/calling-stan-routines-from-a-c-program

and that guy managed to get a hack I suggested to work.

Ben

[Tamas Papp]

I thought (correct me if I am wrong) that for each function you would
need the derivatives wrt to its own arguments, the rest come from the
chain rule.

Best,

Tamas

[Daniel Lee]

Yes, absolutely. Ignore my last email.

[Bob Carpenter]

This is exactly what Michael was warning against --- naively
autodiffing an iterative algorithm. And both Michael and Daniel
are right. It can work to just autodiff an iterative algorithm,
but it can have bad numerical properties. So it's better in terms of
both efficiency and numerical stability to write a separate iterative
algorithm for the derivative that takes its numerical properties
into account.

I also want to qualify that for a full Stan implementation of
a function, you need not only first-order derivatives, but also
second- and third-order. Those usually come about through a templated
implementation of fvar<T> where T is either a var or fvar<var>.
But the same logic applies there --- autodiffing the first-order
Jacobian code to get the Hessian may not be stable. But if all you
need to do is HMC, you can get by with just the var implementation,
which only needs the Jacobian.

Yes, the Stan Math Library is full of examples in C++.

I'm guessing the Julia to Stan combo shot will require you to know
more about C++, not less, than implementing it in C++ directly.

- Bob

[Bob Carpenter]

Right --- the Jacobian's just because it may be a function
mapping N inputs to M outputs, so we need M gradients. Often
it's more efficient to do them all at once, as in many of the
matrix implemntations now (many of which used to just autodiff
the algorithms).

- Bob