Re: [stan-users] expose Stan model user functions to interfaces (rstan)
102 views
Skip to first unread message
Bob Carpenter
Mar 2, 2015, 11:20:56 PM3/2/15
to stan...@googlegroups.com
[moved to stan-dev]
> On Mar 3, 2015, at 3:06 PM, Ben Goodrich <goodri...@gmail.com> wrote:
>
> On Monday, March 2, 2015 at 6:16:27 PM UTC-5, Bob Carpenter wrote:
>
> > On Monday, March 2, 2015 at 10:39:25 PM UTC+1, Ben Goodrich wrote:
> >
> > It would still entail two compilations because the C++ source used to expose the user-defined functions is a shortened and heavily modified version of the C++ source used to draw from the posterior.
>
> Why does it need to be modified? And how are you modifying it?
> • It ignores everything starting with the class definition for the model (because that isn't relevant to the user-defined functions)
> • It takes away the local namespace (because Rcpp can only expose functions from the global namespace), which entails cutting all the using:: statements from what used to be the local namespace and pasting them inside the user-defined functions (so as not to create conflicts in the global namespace)
That's really bad programming practice on Rcpp's part.
They probably did it for "simplicity".
One alternative would be to write global functions that call the namespace-based
ones. Then you shouldn't need to modify the Stan-generated code at all, just
write a wrapper for it.
> • It gets rid of all the templating (because R can only use double or int)
Do you need to get rid of it? Or is it just that it needs to only be instantiated
that way? Again, calling from the outside should take care of this, too.
> • All the user-defined functions that end in _rng have their arguments rewired to not pass a RNG& (because Rcpp can't deal with it). Instead they get an int& seed = 0 argument and at the first line of the body of the function is inserted as static boost::ecuyer1988 base_rng__(seed);. Thus, that gets initialized once with the initial seed but has global scope.
Why can't Rcpp deal with a PRNG reference? Is it just a problem
of passing in from R? That is, there's no way to instantiate an RNG
in R? Can't we do that by wrapping with Rcpp?
> • All the user-defined functions that end in _rng have their arguments rewired to not pass a T_lp_accum__& (because we can't template). Instead that is deleted and the first line of the body of the function is inserted as stan::math::accumulator<double> lp_accum__;.
I think this is problematic because it won't let you test the lp__ increments
inside the function. Would it again be possible to call this from the outside?
> • All the user-defined functions have their arguments rewired to not pass a std::ostream* pstream__ (because Rcpp can't deal with it). Instead that is deleted and the first or second line of the body of the function is inserted as std::ostream* pstream__ = &Rcpp::Rcout;. Thus, if the user-defined function has print() statements, the Rcpp::Rcout thing makes sure they end up in R somehow.
Why not? I have a feeling this is related to two questions back.
> • All the inline declarations get taken away for some reason I can't remember (probably a Rcpp limitation)
Is that just for the wrapper? Again, calling from the outside would help
rather than modifying the existing function. It's the modifications that
I think will be hard to maintain.
> • The relevant functions (but not their declarations if there are declarations) get marked with the attribute // [[Rcpp::export]] before their definitions
I don't even know what that means. Is that a comment in
the C++ code that Rcpp can look at?
> • A variety of things get deleted from the corresponding functors. Perhaps the functors themselves should be deleted because I don't think they can be called from R anyway.
The functors are there because I need them to pass arguments
to the ODE solver. Is it that they can't be called or that
Rcpp places restrictions on what can even be compiled? I can't
imagine it's the latter or we'd never be able to call a Stan function.
> • Everything that looks like a templated type gets changed to a double, inclusive of all the internal const static bool propto__ = true; being changed to false so that (correct) log-density calculations will match those in R (if Stan is using the same parameterization)
That should be OK, but again I think it'd be better handled from
the outside.
> • There is some special handling of get_lp()
> • Dependencies are added to StanHeaders, BH, RcppEigen, and Rcpp so that it can find its #includes (because the compiler is not provided our ordinary CFLAGS)
> That's all. So, yes it probably will be a nightmere to maintain, but all of that was a result of it first failing to compile without it. So, if that is the price that we have to pay for users to be able to test their user-defined functions (and Sebastian can do his posterior predictive stuff outside of generated quantities) then so be it. It would be good to have a test .stan file with every conceivable kind of valid user-defined function (Sebastian got 90% of the way there by himself), so that Travis or something could alert us when it inevitably gets broken.
I'd rather take the approach of making things more modular.
> Ideally, as far as making something like this a feature that the other Stan interfaces could support, it would probably be better to have a flag that could be passed to stanc so that it would emit C++ code that would obviate the need for as many of the above tweaks as possible. But that would require a significant undertaking on the part of someone who understands stanc whereas my approach just requires a significant undertaking on the part of someone who understands regexes.
That's what I was hoping we could do. It'd be pretty straightforward.
But I think it's going to be hard to coordinate this via Google+ or
e-mail, so maybe we can take it up more seriously when I get back and
you and Daniel and I can sit down in a room and go over it.
- Bob
> On Mar 3, 2015, at 3:06 PM, Ben Goodrich <goodri...@gmail.com> wrote:
>
> On Monday, March 2, 2015 at 6:16:27 PM UTC-5, Bob Carpenter wrote:
>
> > On Monday, March 2, 2015 at 10:39:25 PM UTC+1, Ben Goodrich wrote:
> >
> > It would still entail two compilations because the C++ source used to expose the user-defined functions is a shortened and heavily modified version of the C++ source used to draw from the posterior.
>
> Why does it need to be modified? And how are you modifying it?
> • It ignores everything starting with the class definition for the model (because that isn't relevant to the user-defined functions)
> • It takes away the local namespace (because Rcpp can only expose functions from the global namespace), which entails cutting all the using:: statements from what used to be the local namespace and pasting them inside the user-defined functions (so as not to create conflicts in the global namespace)
That's really bad programming practice on Rcpp's part.
They probably did it for "simplicity".
One alternative would be to write global functions that call the namespace-based
ones. Then you shouldn't need to modify the Stan-generated code at all, just
write a wrapper for it.
> • It gets rid of all the templating (because R can only use double or int)
Do you need to get rid of it? Or is it just that it needs to only be instantiated
that way? Again, calling from the outside should take care of this, too.
> • All the user-defined functions that end in _rng have their arguments rewired to not pass a RNG& (because Rcpp can't deal with it). Instead they get an int& seed = 0 argument and at the first line of the body of the function is inserted as static boost::ecuyer1988 base_rng__(seed);. Thus, that gets initialized once with the initial seed but has global scope.
Why can't Rcpp deal with a PRNG reference? Is it just a problem
of passing in from R? That is, there's no way to instantiate an RNG
in R? Can't we do that by wrapping with Rcpp?
> • All the user-defined functions that end in _rng have their arguments rewired to not pass a T_lp_accum__& (because we can't template). Instead that is deleted and the first line of the body of the function is inserted as stan::math::accumulator<double> lp_accum__;.
I think this is problematic because it won't let you test the lp__ increments
inside the function. Would it again be possible to call this from the outside?
> • All the user-defined functions have their arguments rewired to not pass a std::ostream* pstream__ (because Rcpp can't deal with it). Instead that is deleted and the first or second line of the body of the function is inserted as std::ostream* pstream__ = &Rcpp::Rcout;. Thus, if the user-defined function has print() statements, the Rcpp::Rcout thing makes sure they end up in R somehow.
Why not? I have a feeling this is related to two questions back.
> • All the inline declarations get taken away for some reason I can't remember (probably a Rcpp limitation)
Is that just for the wrapper? Again, calling from the outside would help
rather than modifying the existing function. It's the modifications that
I think will be hard to maintain.
> • The relevant functions (but not their declarations if there are declarations) get marked with the attribute // [[Rcpp::export]] before their definitions
I don't even know what that means. Is that a comment in
the C++ code that Rcpp can look at?
> • A variety of things get deleted from the corresponding functors. Perhaps the functors themselves should be deleted because I don't think they can be called from R anyway.
The functors are there because I need them to pass arguments
to the ODE solver. Is it that they can't be called or that
Rcpp places restrictions on what can even be compiled? I can't
imagine it's the latter or we'd never be able to call a Stan function.
> • Everything that looks like a templated type gets changed to a double, inclusive of all the internal const static bool propto__ = true; being changed to false so that (correct) log-density calculations will match those in R (if Stan is using the same parameterization)
That should be OK, but again I think it'd be better handled from
the outside.
> • There is some special handling of get_lp()
> • Dependencies are added to StanHeaders, BH, RcppEigen, and Rcpp so that it can find its #includes (because the compiler is not provided our ordinary CFLAGS)
> That's all. So, yes it probably will be a nightmere to maintain, but all of that was a result of it first failing to compile without it. So, if that is the price that we have to pay for users to be able to test their user-defined functions (and Sebastian can do his posterior predictive stuff outside of generated quantities) then so be it. It would be good to have a test .stan file with every conceivable kind of valid user-defined function (Sebastian got 90% of the way there by himself), so that Travis or something could alert us when it inevitably gets broken.
I'd rather take the approach of making things more modular.
> Ideally, as far as making something like this a feature that the other Stan interfaces could support, it would probably be better to have a flag that could be passed to stanc so that it would emit C++ code that would obviate the need for as many of the above tweaks as possible. But that would require a significant undertaking on the part of someone who understands stanc whereas my approach just requires a significant undertaking on the part of someone who understands regexes.
That's what I was hoping we could do. It'd be pretty straightforward.
But I think it's going to be hard to coordinate this via Google+ or
e-mail, so maybe we can take it up more seriously when I get back and
you and Daniel and I can sit down in a room and go over it.
- Bob
Ben Goodrich
Mar 2, 2015, 11:59:48 PM3/2/15
to stan...@googlegroups.com
On Monday, March 2, 2015 at 11:20:56 PM UTC-5, Bob Carpenter wrote:
One alternative would be to write global functions that call the namespace-based
ones. Then you shouldn't need to modify the Stan-generated code at all, just
write a wrapper for it.
Possibly. Hopefully.
> • It gets rid of all the templating (because R can only use double or int)
Do you need to get rid of it? Or is it just that it needs to only be instantiated
that way? Again, calling from the outside should take care of this, too.
Yeah, whatever function is being called from R can basically only pass (and return)
- double or int
- std::vector<> with double or int inside the <>
- matrix_d, vector_d, or row_vector_d
> • All the user-defined functions that end in _rng have their arguments rewired to not pass a RNG& (because Rcpp can't deal with it). Instead they get an int& seed = 0 argument and at the first line of the body of the function is inserted as static boost::ecuyer1988 base_rng__(seed);. Thus, that gets initialized once with the initial seed but has global scope.
Why can't Rcpp deal with a PRNG reference? Is it just a problem
of passing in from R? That is, there's no way to instantiate an RNG
in R?
I didn't totally figure out _why_ that was the case but yes, Rcpp wouldn't instantiate a boost::ecuyer1988 based on what was passed from R. I tried for a few days but got stuck on this error message
https://groups.google.com/d/msg/stan-users/wn3gZZvWi8g/uipFTfcQbzkJ
> • All the user-defined functions that end in _rng have their arguments rewired to not pass a T_lp_accum__& (because we can't template). Instead that is deleted and the first line of the body of the function is inserted as stan::math::accumulator<double> lp_accum__;.
I think this is problematic because it won't let you test the lp__ increments
inside the function. Would it again be possible to call this from the outside?
It is kind of limited. Currently, it passes in 0 from R and you can only test that get_lp() is right for that function. If you wanted to test a sequence of function calls, you would have to pass get_lp() back to R and pass that back to the lp__ argument for the next function ending in __lp.
> • All the user-defined functions have their arguments rewired to not pass a std::ostream* pstream__ (because Rcpp can't deal with it). Instead that is deleted and the first or second line of the body of the function is inserted as std::ostream* pstream__ = &Rcpp::Rcout;. Thus, if the user-defined function has print() statements, the Rcpp::Rcout thing makes sure they end up in R somehow.
Why not? I have a feeling this is related to two questions back.
Yes. Mav came up with the current solution
https://groups.google.com/d/msg/stan-users/wn3gZZvWi8g/LDxJfxdTij8J
> • All the inline declarations get taken away for some reason I can't remember (probably a Rcpp limitation)
Is that just for the wrapper? Again, calling from the outside would help
rather than modifying the existing function. It's the modifications that
I think will be hard to maintain.
Still don't remember. Maybe the compiler inlined a function and then R couldn't call it or something.
> • The relevant functions (but not their declarations if there are declarations) get marked with the attribute // [[Rcpp::export]] before their definitions
I don't even know what that means. Is that a comment in
the C++ code that Rcpp can look at?
Yes, that is how Rcpp knows which functions to expose to R. That part is the least of the problems.
> • A variety of things get deleted from the corresponding functors. Perhaps the functors themselves should be deleted because I don't think they can be called from R anyway.
The functors are there because I need them to pass arguments
to the ODE solver. Is it that they can't be called or that
Rcpp places restrictions on what can even be compiled? I can't
imagine it's the latter or we'd never be able to call a Stan function.
They compile if you take away some of the extraneous arguments but I didn't know what they are for. Sebastian's file had some integrate_ode() stuff I think, so I guess it works at the moment.
I'd rather take the approach of making things more modular.
Yes, I think having global wrapper functions is a good idea when you get back.
Ben
Bob Carpenter
Mar 3, 2015, 1:19:38 AM3/3/15
to stan...@googlegroups.com
Great. That all sounds like we can get it under control.
I'm working right now on plumbing a runtime Stan-source line counter
through the code, which will throw another wrench into the works.
Boost Spirit is a real piece of work --- it's a crazy hack in
the grammar to get source line numbers in, something I'd think
almost everyone would want (it's easy for errors, but not for
non-errors). As the answer on Stackoverflow said, "and then we do
a little dance with "raw[ ... ]". I'm going to use "little dance"
from now on instead of "hack".
- Bob
> --
> You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
I'm working right now on plumbing a runtime Stan-source line counter
through the code, which will throw another wrench into the works.
Boost Spirit is a real piece of work --- it's a crazy hack in
the grammar to get source line numbers in, something I'd think
almost everyone would want (it's easy for errors, but not for
non-errors). As the answer on Stackoverflow said, "and then we do
a little dance with "raw[ ... ]". I'm going to use "little dance"
from now on instead of "hack".
- Bob
> You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
Allen B. Riddell
Mar 3, 2015, 8:53:44 AM3/3/15
to stan...@googlegroups.com
An alternative route might be to figure out something with Protocol
Buffers since that has C++, Python, and R support. I'll note that Dirk
Eddelbuettel (Rcpp) is also the author of RProtoBuf. Might be a bit
cleaner and could solve the problems around the sampler arguments (which
need to accept vector/matrix initialization values).
Best,
Allen
On 03/02, Ben Goodrich wrote:
> On Monday, March 2, 2015 at 11:20:56 PM UTC-5, Bob Carpenter wrote:
> Yeah, whatever function is being called from R can basically only pass (and
> return)
>
> - double or int
> - std::vector<> with double or int inside the <>
> - matrix_d, vector_d, or row_vector_d
Buffers since that has C++, Python, and R support. I'll note that Dirk
Eddelbuettel (Rcpp) is also the author of RProtoBuf. Might be a bit
cleaner and could solve the problems around the sampler arguments (which
need to accept vector/matrix initialization values).
Best,
Allen
On 03/02, Ben Goodrich wrote:
> On Monday, March 2, 2015 at 11:20:56 PM UTC-5, Bob Carpenter wrote:
> Yeah, whatever function is being called from R can basically only pass (and
> return)
>
> - std::vector<> with double or int inside the <>
> - matrix_d, vector_d, or row_vector_d
Bob Carpenter
Mar 3, 2015, 9:01:14 AM3/3/15
to stan...@googlegroups.com
Just to be clear, the context here is calling Stan functions, not
doing something like inits. I thought protocol buffers were dead
and JSON won. You can also do that in any of these languages.
I don't really care one way or another other than that it be easily
understood by users.
What's the "C++ support"? They're not built into C++, so any of
these things we use will need a 3rd party library.
Looking at the Google page (I believe they came out of Google), I
have to say I like the typing compared to JSON after a 30s quick
scan:
https://developers.google.com/protocol-buffers/docs/proto
- Bob
doing something like inits. I thought protocol buffers were dead
and JSON won. You can also do that in any of these languages.
I don't really care one way or another other than that it be easily
understood by users.
What's the "C++ support"? They're not built into C++, so any of
these things we use will need a 3rd party library.
Looking at the Google page (I believe they came out of Google), I
have to say I like the typing compared to JSON after a 30s quick
scan:
https://developers.google.com/protocol-buffers/docs/proto
- Bob
Allen B. Riddell
Mar 3, 2015, 10:17:57 AM3/3/15
to stan...@googlegroups.com
I guess I see these things as connected: figuring out structured data
exchange to and from C++. For example, I'd love to figure out a way to
to call the "Stan function" split_rhat from R and Python --- which I
think would be an excellent idea if it's going to change from release to
release. Casting R and Python data into Eigen data structures is too
painful, IMHO.
Doing JSON in C++ does not look like fun at all. That's why I've been
looking at Protocol Buffers. (The new protocol buffers 3.0 also defines
a canonical serialization to JSON, which might be nice if one wants to
move away from rdump.)
Of course, I should add that I have this vision of reducing the RStan
and PyStan codebase by 30-40% by moving just about everything into Stan.
I'm not sure how widely held this ambition is...
Thanks,
Allen
p.s. I'm happy to work on any side of this (R, C++, Python) if there's
consensus about the best way to do this.
exchange to and from C++. For example, I'd love to figure out a way to
to call the "Stan function" split_rhat from R and Python --- which I
think would be an excellent idea if it's going to change from release to
release. Casting R and Python data into Eigen data structures is too
painful, IMHO.
Doing JSON in C++ does not look like fun at all. That's why I've been
looking at Protocol Buffers. (The new protocol buffers 3.0 also defines
a canonical serialization to JSON, which might be nice if one wants to
move away from rdump.)
Of course, I should add that I have this vision of reducing the RStan
and PyStan codebase by 30-40% by moving just about everything into Stan.
I'm not sure how widely held this ambition is...
Thanks,
Allen
p.s. I'm happy to work on any side of this (R, C++, Python) if there's
consensus about the best way to do this.
Bob Carpenter
Mar 4, 2015, 7:36:13 AM3/4/15
to stan...@googlegroups.com
> On Mar 4, 2015, at 2:17 AM, Allen B. Riddell <a...@ariddell.org> wrote:
>
> I guess I see these things as connected: figuring out structured data
> exchange to and from C++.
serialize data to a byte stream. But it's hugely inefficient compared
to just using data structure directly.
I just want to use JSON or protocol buffers or whatever to do
deserialization and serialization to output streams (particularly
files, but it would also work over a network).
Maybe efficiency isn't such a big deal, though.
> For example, I'd love to figure out a way to
> to call the "Stan function" split_rhat from R and Python --- which I
> think would be an excellent idea if it's going to change from release to
> release. Casting R and Python data into Eigen data structures is too
> painful, IMHO.
structures, but it's totally straightforward in R. Would it help
if we provided other data structures for interfaces and did the
conversion on the C++ side? If so, what do you want in the way of
interfaces. C-style arrays, nested C++ std::vector, Boost arrays,
or what?
> Doing JSON in C++ does not look like fun at all.
*** Stan already has a relatively efficient JSON parser written in C++ ***
> That's why I've been
> looking at Protocol Buffers. (The new protocol buffers 3.0 also defines
> a canonical serialization to JSON, which might be nice if one wants to
> move away from rdump.)
to doc.
> Of course, I should add that I have this vision of reducing the RStan
> and PyStan codebase by 30-40% by moving just about everything into Stan.
> I'm not sure how widely held this ambition is...
just that topic in the Stan meeting yesterday. Michael put up
a design doc on the Wiki, so please have a look and comment.
> Thanks,
>
> Allen
>
> p.s. I'm happy to work on any side of this (R, C++, Python) if there's
> consensus about the best way to do this.
- Bob
Allen B. Riddell
Mar 4, 2015, 8:50:43 AM3/4/15
to stan...@googlegroups.com
On 03/04, Bob Carpenter wrote:
>
>
> > On Mar 4, 2015, at 2:17 AM, Allen B. Riddell <a...@ariddell.org> wrote:
> >
> > I guess I see these things as connected: figuring out structured data
> > exchange to and from C++.
>
> JSON and protocol buffers provide a nice cross-language way to
> serialize data to a byte stream. But it's hugely inefficient compared
> to just using data structure directly.
>
> I just want to use JSON or protocol buffers or whatever to do
> deserialization and serialization to output streams (particularly
> files, but it would also work over a network).
>
> Maybe efficiency isn't such a big deal, though.
>
If it's just the data for var_context (initialization and data) then
> >
> > I guess I see these things as connected: figuring out structured data
> > exchange to and from C++.
>
> JSON and protocol buffers provide a nice cross-language way to
> serialize data to a byte stream. But it's hugely inefficient compared
> to just using data structure directly.
>
> I just want to use JSON or protocol buffers or whatever to do
> deserialization and serialization to output streams (particularly
> files, but it would also work over a network).
>
> Maybe efficiency isn't such a big deal, though.
>
this typically will not be that much data, especially if there was a
binary serialization.
> > For example, I'd love to figure out a way to
> > to call the "Stan function" split_rhat from R and Python --- which I
> > think would be an excellent idea if it's going to change from release to
> > release. Casting R and Python data into Eigen data structures is too
> > painful, IMHO.
>
> Where's the pain? I don't know anything about Python data
> structures, but it's totally straightforward in R. Would it help
> if we provided other data structures for interfaces and did the
> conversion on the C++ side? If so, what do you want in the way of
> interfaces. C-style arrays, nested C++ std::vector, Boost arrays,
> or what?
>
good idea?
The pain comes about because the data chains are vectors of doubles and
are not in contiguous memory in R or in Python -- and casting these
vectors to Eigen requires wrapping Eigen.
I actually don't have a good solution to this problem (that doesn't
involve wrapping Eigen). I could use help here. This is probably lower
hanging fruit than the command API refactor (but related).
> > Doing JSON in C++ does not look like fun at all.
>
> Let me stress this again:
>
> *** Stan already has a relatively efficient JSON parser written in C++ ***
>
as somewhat cumbersome vs. Protocol Buffers. I'm happy to use JSON. I
should also clarify that I mean going from a JSON serialized C++ class
(or similar) to the C++ class itself --- that's what's (moderately)
painful.
Bob Carpenter
Mar 5, 2015, 2:13:13 AM3/5/15
to stan...@googlegroups.com
> On Mar 5, 2015, at 12:50 AM, Allen B. Riddell <a...@ariddell.org> wrote:
>
> On 03/04, Bob Carpenter wrote:
>>
>> Maybe efficiency isn't such a big deal, though.
>>
>
> If it's just the data for var_context (initialization and data) then
> this typically will not be that much data, especially if there was a
> binary serialization.
but they're a huge pain to maintain, doc, etc., and have the downside
that users can't inspect them by hand or with other programs without
a custom reader. (Now if protocol buffers have a binary format, that
might help.) There's also portability issues across 32-bit, 64-bit,
little endian and big endian, though I don't know how much of an issue
that is these days and if protocol buffers have binary formats, it's
probably already solved for us.
...
> Fair enough, R does have REigen. But is having interfaces wrap Eigen a
> good idea?
>
> The pain comes about because the data chains are vectors of doubles and
> are not in contiguous memory in R or in Python -- and casting these
> vectors to Eigen requires wrapping Eigen.
>
> I actually don't have a good solution to this problem (that doesn't
> involve wrapping Eigen). I could use help here. This is probably lower
> hanging fruit than the command API refactor (but related).
this discussion.
Why would you need to wrap Eigen? Can't you just
write C++ that just includes the Eigen lib? Do we somehow need
to expose Eigen data structures to users? For the same reason, I'm also
confused why RStan would need to use RcppEigen.
If you don't want to use Eigen, what would you like to use?
std::vector<std::vector<T>> is too cumbersome and many of our interfaces
are explicitly matrix based, which is why we chose to go with Eigen.
I agree that it's not the most portable thing in the world, which is
why I needed a lot of convincing to agree to going with Eigen.
...
> should also clarify that I mean going from a JSON serialized C++ class
> (or similar) to the C++ class itself --- that's what's (moderately)
> painful.
can be anything. We could write specific cases if the data's all
strings or numbers.
- Bob
Allen B. Riddell
Mar 5, 2015, 9:51:54 AM3/5/15
to stan...@googlegroups.com
On 03/05, Bob Carpenter wrote:
>
> > On Mar 5, 2015, at 12:50 AM, Allen B. Riddell <a...@ariddell.org> wrote:
> >
> > On 03/04, Bob Carpenter wrote:
> >>
>
> ...
>
> >> Maybe efficiency isn't such a big deal, though.
> >>
> >
> > If it's just the data for var_context (initialization and data) then
> > this typically will not be that much data, especially if there was a
> > binary serialization.
>
> We've discussed binary formats before. They'd be more efficient,
> but they're a huge pain to maintain, doc, etc., and have the downside
> that users can't inspect them by hand or with other programs without
> a custom reader. (Now if protocol buffers have a binary format, that
> might help.) There's also portability issues across 32-bit, 64-bit,
> little endian and big endian, though I don't know how much of an issue
> that is these days and if protocol buffers have binary formats, it's
> probably already solved for us.
>
Yes, binary and text (JSON in Protocol Buffers 3.0) serialization is
>
> > On Mar 5, 2015, at 12:50 AM, Allen B. Riddell <a...@ariddell.org> wrote:
> >
> > On 03/04, Bob Carpenter wrote:
> >>
>
> ...
>
> >> Maybe efficiency isn't such a big deal, though.
> >>
> >
> > If it's just the data for var_context (initialization and data) then
> > this typically will not be that much data, especially if there was a
> > binary serialization.
>
> We've discussed binary formats before. They'd be more efficient,
> but they're a huge pain to maintain, doc, etc., and have the downside
> that users can't inspect them by hand or with other programs without
> a custom reader. (Now if protocol buffers have a binary format, that
> might help.) There's also portability issues across 32-bit, 64-bit,
> little endian and big endian, though I don't know how much of an issue
> that is these days and if protocol buffers have binary formats, it's
> probably already solved for us.
>
available. (As an aside, I'm not a huge fan of PB in general but for
dealing with C++ interop specifically it seems like the best tool.)
> ...
>
> > I actually don't have a good solution to this problem (that doesn't
> > involve wrapping Eigen). I could use help here. This is probably lower
> > hanging fruit than the command API refactor (but related).
>
> I'm still confused, so I'm probably missing the whole point of
> this discussion.
>
> Why would you need to wrap Eigen? Can't you just
> write C++ that just includes the Eigen lib? Do we somehow need
> to expose Eigen data structures to users? For the same reason, I'm also
> confused why RStan would need to use RcppEigen.
>
> If you don't want to use Eigen, what would you like to use?
> std::vector<std::vector<T>> is too cumbersome and many of our interfaces
> are explicitly matrix based, which is why we chose to go with Eigen.
> I agree that it's not the most portable thing in the world, which is
> why I needed a lot of convincing to agree to going with Eigen.
>
Ok, I think we're on the same page here. Yes, I think the signature
> > I actually don't have a good solution to this problem (that doesn't
> > involve wrapping Eigen). I could use help here. This is probably lower
> > hanging fruit than the command API refactor (but related).
>
> I'm still confused, so I'm probably missing the whole point of
> this discussion.
>
> Why would you need to wrap Eigen? Can't you just
> write C++ that just includes the Eigen lib? Do we somehow need
> to expose Eigen data structures to users? For the same reason, I'm also
> confused why RStan would need to use RcppEigen.
>
> If you don't want to use Eigen, what would you like to use?
> std::vector<std::vector<T>> is too cumbersome and many of our interfaces
> are explicitly matrix based, which is why we chose to go with Eigen.
> I agree that it's not the most portable thing in the world, which is
> why I needed a lot of convincing to agree to going with Eigen.
>
should be std::vector<std::vector<T>> and that one should guarantee (not
hard, IIRC) that within each std::vector<T> memory is contiguous. I
*think* you can cast that to an Eigen matrix (in Stan). This would avoid
having Python or R touch Eigen.
> > should also clarify that I mean going from a JSON serialized C++ class
> > (or similar) to the C++ class itself --- that's what's (moderately)
> > painful.
>
> I think that's pretty much impossible in general, since C++ classes
> can be anything. We could write specific cases if the data's all
> strings or numbers.
easily. Remember what var_context is:
std::map<std::string,
std::pair<std::vector<double>,
std::vector<size_t> > > vars_r_;
std::map<std::string,
std::pair<std::vector<int>,
std::vector<size_t> > > vars_i_;
std::vector<double> const empty_vec_r_;
std::vector<int> const empty_vec_i_;
std::vector<size_t> const empty_vec_ui_;
Serializing and deserializing (binary or text) is the sort of thing
Protocol Buffers was made for as far as I can tell.
Best wishes,
Allen
Bob Carpenter
Mar 5, 2015, 10:03:14 AM3/5/15
to stan...@googlegroups.com
> On Mar 6, 2015, at 1:51 AM, Allen B. Riddell <a...@ariddell.org> wrote:
>
> On 03/05, Bob Carpenter wrote:
>
> Yes, binary and text (JSON in Protocol Buffers 3.0) serialization is
> available. (As an aside, I'm not a huge fan of PB in general but for
> dealing with C++ interop specifically it seems like the best tool.)
Sounds good.
> On 03/05, Bob Carpenter wrote:
>
> Yes, binary and text (JSON in Protocol Buffers 3.0) serialization is
> available. (As an aside, I'm not a huge fan of PB in general but for
> dealing with C++ interop specifically it seems like the best tool.)
...
> Ok, I think we're on the same page here. Yes, I think the signature
> should be std::vector<std::vector<T>> and that one should guarantee (not
> hard, IIRC) that within each std::vector<T> memory is contiguous. I
> *think* you can cast that to an Eigen matrix (in Stan). This would avoid
> having Python or R touch Eigen.
But Eigen::Matrix does. You can't quite cast---you need
to use their Eigen::Map expression template.
We don't want to use vector<vector<T>> for anything compute
intensive because of the lack of contiguity guarantee.
>>> should also clarify that I mean going from a JSON serialized C++ class
>>> (or similar) to the C++ class itself --- that's what's (moderately)
>>> painful.
>>
>> I think that's pretty much impossible in general, since C++ classes
>> can be anything. We could write specific cases if the data's all
>> strings or numbers.
>
> But in the specific case of var_context, Protocol Buffers could do it
> easily.
How was Eigen related to this?
> Remember what var_context is:
>
> std::map<std::string,
> std::pair<std::vector<double>,
> std::vector<size_t> > > vars_r_;
> std::map<std::string,
> std::pair<std::vector<int>,
> std::vector<size_t> > > vars_i_;
> std::vector<double> const empty_vec_r_;
> std::vector<int> const empty_vec_i_;
> std::vector<size_t> const empty_vec_ui_;
a specific implementation of var_context I used for the dump
data format.
I think if we went to something like JSON or protocol buffers,
then following the shape of the data would be more natural
than this mess, which I only introduced for convenience for the
case at hand.
> Serializing and deserializing (binary or text) is the sort of thing
> Protocol Buffers was made for as far as I can tell.
C++ classes, which can have anything in them, like pointers and
references and file handles and ... That's why I was confused
earlier.
Unlike Java, there's no general built-in way to serialize a C++
class. You can just grab the memory, but it might implicitly
be pointing off to the heap, as in both std::vector and Eigen::Matrix.
- Bob
Allen B. Riddell
Mar 5, 2015, 11:55:28 AM3/5/15
to stan...@googlegroups.com
On 03/06, Bob Carpenter wrote:
>
> > On Mar 6, 2015, at 1:51 AM, Allen B. Riddell <a...@ariddell.org> wrote:
> >
> > On 03/05, Bob Carpenter wrote:
> >
> > Yes, binary and text (JSON in Protocol Buffers 3.0) serialization is
> > available. (As an aside, I'm not a huge fan of PB in general but for
> > dealing with C++ interop specifically it seems like the best tool.)
>
> Sounds good.
>
>
> ...
>
> > Ok, I think we're on the same page here. Yes, I think the signature
> > should be std::vector<std::vector<T>> and that one should guarantee (not
> > hard, IIRC) that within each std::vector<T> memory is contiguous. I
> > *think* you can cast that to an Eigen matrix (in Stan). This would avoid
> > having Python or R touch Eigen.
>
> Not quite. vector<vector<T>> has no contiguity guarantee.
> But Eigen::Matrix does. You can't quite cast---you need
> to use their Eigen::Map expression template.
>
> We don't want to use vector<vector<T>> for anything compute
> intensive because of the lack of contiguity guarantee.
What if there's contiguous memory for the inner vector, the vector<T>?
>
> > On Mar 6, 2015, at 1:51 AM, Allen B. Riddell <a...@ariddell.org> wrote:
> >
> > On 03/05, Bob Carpenter wrote:
> >
> > Yes, binary and text (JSON in Protocol Buffers 3.0) serialization is
> > available. (As an aside, I'm not a huge fan of PB in general but for
> > dealing with C++ interop specifically it seems like the best tool.)
>
> Sounds good.
>
>
> ...
>
> > Ok, I think we're on the same page here. Yes, I think the signature
> > should be std::vector<std::vector<T>> and that one should guarantee (not
> > hard, IIRC) that within each std::vector<T> memory is contiguous. I
> > *think* you can cast that to an Eigen matrix (in Stan). This would avoid
> > having Python or R touch Eigen.
>
> Not quite. vector<vector<T>> has no contiguity guarantee.
> But Eigen::Matrix does. You can't quite cast---you need
> to use their Eigen::Map expression template.
>
> We don't want to use vector<vector<T>> for anything compute
> intensive because of the lack of contiguity guarantee.
Can Eigen work with that? Or does the entire thing need to be
contiguous?
>
> >>> should also clarify that I mean going from a JSON serialized C++ class
> >>> (or similar) to the C++ class itself --- that's what's (moderately)
> >>> painful.
> >>
> >> I think that's pretty much impossible in general, since C++ classes
> >> can be anything. We could write specific cases if the data's all
> >> strings or numbers.
> >
> > But in the specific case of var_context, Protocol Buffers could do it
> > easily.
>
> Agreed. I didn't realize that's what we were talking about.
> How was Eigen related to this?
>
serialize, no matter how efficiently one goes about it.
>
> > Remember what var_context is:
> >
> > std::map<std::string,
> > std::pair<std::vector<double>,
> > std::vector<size_t> > > vars_r_;
> > std::map<std::string,
> > std::pair<std::vector<int>,
> > std::vector<size_t> > > vars_i_;
> > std::vector<double> const empty_vec_r_;
> > std::vector<int> const empty_vec_i_;
> > std::vector<size_t> const empty_vec_ui_;
>
> I like to think of the var_context as an interface. That's
> a specific implementation of var_context I used for the dump
> data format.
>
> I think if we went to something like JSON or protocol buffers,
> then following the shape of the data would be more natural
> than this mess, which I only introduced for convenience for the
> case at hand.
>
>
> > Serializing and deserializing (binary or text) is the sort of thing
> > Protocol Buffers was made for as far as I can tell.
>
> Yes, but only for structured data like this --- not for general
> C++ classes, which can have anything in them, like pointers and
> references and file handles and ... That's why I was confused
> earlier.
separately. But I don't think we pass file handles except for the stderr
stream, right?
Just to check, I'm looking at the stan_args configuration class in R --
I think the function and file pointers are already handled separately. I
think the entire thing could be Protocol Buffer/JSON-ized. (In some
respects this seems easier than var_context.)
Data in the C++ class rstan::stan_args (which is, btw, 600 LOC):
unsigned int random_seed;
unsigned int chain_id;
std::string init;
SEXP init_list;
double init_radius;
std::string sample_file; // the file for outputting the samples
bool append_samples;
bool sample_file_flag; // true: write out to a file; false, do not
stan_args_method_t method;
std::string diagnostic_file;
bool diagnostic_file_flag;
union {
struct {
int iter; // number of iterations
int refresh; //
sampling_algo_t algorithm;
int warmup; // number of warmup
int thin;
bool save_warmup; // weather to save warmup samples (always true now)
int iter_save; // number of iterations saved
int iter_save_wo_warmup; // number of iterations saved wo warmup
bool adapt_engaged;
double adapt_gamma;
double adapt_delta;
double adapt_kappa;
unsigned int adapt_init_buffer;
unsigned int adapt_term_buffer;
unsigned int adapt_window;
double adapt_t0;
sampling_metric_t metric; // UNIT_E, DIAG_E, DENSE_E;
double stepsize; // defaut to 1;
double stepsize_jitter;
int max_treedepth; // for NUTS, default to 10.
double int_time; // for HMC, default to 2 * pi
} sampling;
struct {
int iter; // default to 2000
int refresh; // default to 100
optim_algo_t algorithm; // Newton, (L)BFGS
bool save_iterations; // default to false
double init_alpha; // default to 0.001, for (L)BFGS
double tol_obj; // default to 1e-12, for (L)BFGS
double tol_grad; // default to 1e-8, for (L)BFGS
double tol_param; // default to 1e-8, for (L)BFGS
double tol_rel_obj; // default to 1e4, for (L)BFGS
double tol_rel_grad; // default to 1e7, for (L)BFGS
int history_size; // default to 5, for LBFGS only
} optim;
struct {
double epsilon; // default to 1e-6, for test_grad
double error; // default to 1e-6, for test_grad
} test_grad;
} ctrl;
Allen B. Riddell
Mar 5, 2015, 12:58:03 PM3/5/15
to stan...@googlegroups.com
On 03/05, Allen B. Riddell wrote:
> On 03/06, Bob Carpenter wrote:
> > > Ok, I think we're on the same page here. Yes, I think the signature
> > > should be std::vector<std::vector<T>> and that one should guarantee (not
> > > hard, IIRC) that within each std::vector<T> memory is contiguous. I
> > > *think* you can cast that to an Eigen matrix (in Stan). This would avoid
> > > having Python or R touch Eigen.
> >
> > Not quite. vector<vector<T>> has no contiguity guarantee.
> > But Eigen::Matrix does. You can't quite cast---you need
> > to use their Eigen::Map expression template.
> >
> > We don't want to use vector<vector<T>> for anything compute
> > intensive because of the lack of contiguity guarantee.
>
Another idea:
> On 03/06, Bob Carpenter wrote:
> > > Ok, I think we're on the same page here. Yes, I think the signature
> > > should be std::vector<std::vector<T>> and that one should guarantee (not
> > > hard, IIRC) that within each std::vector<T> memory is contiguous. I
> > > *think* you can cast that to an Eigen matrix (in Stan). This would avoid
> > > having Python or R touch Eigen.
> >
> > Not quite. vector<vector<T>> has no contiguity guarantee.
> > But Eigen::Matrix does. You can't quite cast---you need
> > to use their Eigen::Map expression template.
> >
> > We don't want to use vector<vector<T>> for anything compute
> > intensive because of the lack of contiguity guarantee.
>
What if instead of std::vector<std::vector<T>> we just had
std::vector<T> (contiguous) and a std::vector<size_t> for the chain
lengths? Would that allow R and Python to use the rhat and ess services
of Stan without wrapping Eigen?
Bob Carpenter
Mar 5, 2015, 7:03:14 PM3/5/15
to stan...@googlegroups.com
Yes, we could start index fiddling with single arrays.
But it'd be easier to have a vector<vector<double>> interface.
We may wind up having to go to vector<vector<double>> anyway
when we start using ragged arrays for comparing chains of
different lengths.
I don't know what you mean by chains being too large
to serialize. CmdStan writes them out to a file just
using ASCII, which is relatively slow, but not really
a bottleneck for Stan.
I wouldn't think that protocol buffers or JSON would replace
var_context --- they'd just be a way to get things into
memory. What we need is not the memory underlying var_context
but the interface that lets us pull out values based on Stan
data types. I'm sure we could implement that on top of JSON
or protocol buffers the same way it's now implemented on top of
our half-assed Rdump implementation (half assed in that it
doesn't support the same range of features as R itself, which
has a no-holds-barred approach to dump).
Michael and Daniel are working on refactoring the entire
interface. So hopefully those 600 lines of code will be
reduced on both the RStan and PyStan sides. It's just that so
many other issues keep coming up that neither of them have time
to work on it properly.
Thanks for doing the rundown of arguments. You can see
something like that in the CmdStan manual, as well, because
I tried to doc it that way.
In the rundown of rstan::stan_args, the SEXP init_list is
the killer --- that's what we'd need something generic like
an in-memory structure for or would need to somehow marshal inside
of Python and R and pass in via memory. That's where the
var_context comes in. Serializing to and from JSON or protocol
buffers should be much much slower than staying in memory --- the
conversion back and forth from a character stream is very expensive.
As Marcus pointed out, this can dominate computation for some of
the efficient optimization models, like linear regression.
I have no idea what optim_alog_t is --- some kind of enum?
If we only ever deal with files in terms of their file names,
it should be easy --- I don't know if Daniel or Michael are
thinking along those lines or in terms of actual streams.
- Bob
But it'd be easier to have a vector<vector<double>> interface.
We may wind up having to go to vector<vector<double>> anyway
when we start using ragged arrays for comparing chains of
different lengths.
I don't know what you mean by chains being too large
to serialize. CmdStan writes them out to a file just
using ASCII, which is relatively slow, but not really
a bottleneck for Stan.
I wouldn't think that protocol buffers or JSON would replace
var_context --- they'd just be a way to get things into
memory. What we need is not the memory underlying var_context
but the interface that lets us pull out values based on Stan
data types. I'm sure we could implement that on top of JSON
or protocol buffers the same way it's now implemented on top of
our half-assed Rdump implementation (half assed in that it
doesn't support the same range of features as R itself, which
has a no-holds-barred approach to dump).
Michael and Daniel are working on refactoring the entire
interface. So hopefully those 600 lines of code will be
reduced on both the RStan and PyStan sides. It's just that so
many other issues keep coming up that neither of them have time
to work on it properly.
Thanks for doing the rundown of arguments. You can see
something like that in the CmdStan manual, as well, because
I tried to doc it that way.
In the rundown of rstan::stan_args, the SEXP init_list is
the killer --- that's what we'd need something generic like
an in-memory structure for or would need to somehow marshal inside
of Python and R and pass in via memory. That's where the
var_context comes in. Serializing to and from JSON or protocol
buffers should be much much slower than staying in memory --- the
conversion back and forth from a character stream is very expensive.
As Marcus pointed out, this can dominate computation for some of
the efficient optimization models, like linear regression.
I have no idea what optim_alog_t is --- some kind of enum?
If we only ever deal with files in terms of their file names,
it should be easy --- I don't know if Daniel or Michael are
thinking along those lines or in terms of actual streams.
- Bob
Reply all
Reply to author
Forward
0 new messages