Comments on Command Refactor

[Michael Betancourt]

If people have a minute can you check out the current progress at
feature/issue-1361-writer_callback? This branch introduces a
writer callback and then refactors a ton of services to use the
callback. Some of that refactoring addresses a bunch of
outstanding issues and substantially pushes command towards
the spec outlined at https://github.com/stan-dev/stan/wiki/Stan-API-Refactor.

If people are cool with it I’ll proceed to the next step, which involves
pushing the writer callback all the way through the sampling and
optimization codes and then updating an incredibly large number of
tests (not to mention adding a bunch of new ones).

-Betancourt

[Bob Carpenter]

I did some major edits on the services and interface proposal:

https://github.com/stan-dev/stan/wiki/Stan-API-Refactor

I tried to flesh out the missing arguments to the command
service function, added a discussion of templating vs.
abstract base classes for the callback option, and added
a discussion of just replacing everything with JSON or
Protocol Buffers. I tried to lay out more of what each
component needs as well as discussing some pros and cons
as I see them without trying to editorialize my own opinions.

My own opinion is that we should go the callback route.
I think the pros outweigh the cons, unless the pain to
Allen and Jiqiang of writing C++ is so great as to outweigh
the other benefits.

Feel free to continue editing --- I'm sure I missed some of
the pros and cons.

We need to dive one level deeper and start fleshing out what
each of the elements is going to be. I think I'll go and drop
in the CmdStan config skeletons, which I think show what's going
on.

This also brings up the issue of whether the interfaces provide
a monolithic god-object for config and single function, or whether they
break things down more modularly. I'm in favor of smaller functions,
which is why I documented CmdStan as I did, whereas I believe Michael's
in favor of the god-object approach, which is why the config is
implemented as it is. So the implementation/interface decisions
can go their own way.

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

[Bob Carpenter]

I just went and added all the commands and a discussion
of some of the extra services provided by RStan and PyStan
that aren't part of the proposed Stan C++ services.

- Bob

[Michael Betancourt]

Has anybody looked at the actual refactoring at feature/issue-1361-writer_callback?

[Bob Carpenter]

I'm looking at it now, but am not sure what I'm supposed
to be seeing. It's hard to read these templated functions
because I can't tell what the types are.

It might be prudent in terms of your time to get some
consensus on what things will look like from the
top down in a way that Allen, Jiqiang and Ben can live with
before doing too much refactoring of what we currently have.

I somehow don't think I'm going to get my way in how this
is all designed because Michael wants it one way and Allen and
Jiqiang want it another, and I want something different again.

I had written a whole discourse on how I'd like to see it
but I think that's just going to make you unhappy if I post
it.

- Bob

[Michael Betancourt]

I think it’s important to remember that because none of the
interfaces use command nothing in this branch has a huge
effect on anything but CmdStan so any related comments
are irrelevant.

Seriously, though, this is becoming a problem. Allen essentially
wants a C-interface which is a gigantic change to our current
setup, including the perspective of Stan as a C++ API on the
whole. I’m not even sure what Jiqiang is advocating, but
the whole point of a unified UX is that the interfaces _don’t_
get to do whatever they want. We’re getting dangerously
close to every interface being completely independent, and
once that happens we lose a tremendous potential for Stan.

[Allen B. Riddell]

Hi everyone,

I find Bob's edits really helpful.

The goal of increasing code reuse trumps my concerns about C++ and
callbacks. While I'm still strongly in favor of reducing the amount of C++
code that needs to be maintained from the R and Python side, if we can
get some shared way of validating and populating-with-defaults the
proposed "stan::command::config" instance, then I'll be happy. Sharing rhat
and ess calculations would be next on the list.

I also support breaking up the interface into different functions
(lbfgs, hmc, ...) -- to what extent does this overlap with the RStan
Wiki discussion
https://github.com/stan-dev/rstan/wiki/Interfaces-3.0-Spec?

Thanks,

-a

[Bob Carpenter]

> On Mar 20, 2015, at 10:16 PM, Michael Betancourt <betan...@gmail.com> wrote:
>
> I think it’s important to remember that because none of the
> interfaces use command nothing in this branch has a huge
> effect on anything but CmdStan so any related comments
> are irrelevant.

So what exactly did you want me to look at?

> Seriously, though, this is becoming a problem.

Hopefully we can find a compromise that everyone can live
with.

> Allen essentially
> wants a C-interface which is a gigantic change to our current
> setup, including the perspective of Stan as a C++ API on the
> whole.

I thought he was advocating something that would pass
data as general interface neutral streams, like protocol
buffers. That'd be independent of C or C++ or Python or R.
Of course we'd code it all in C++ on the Stan side.

> I’m not even sure what Jiqiang is advocating, but
> the whole point of a unified UX is that the interfaces _don’t_
> get to do whatever they want.

I agree --- everything should be unified. But I think
we should design with the interfaces in mind --- they're how
almost everyone uses Stan, so getting the interfaces right
is important.

> We’re getting dangerously
> close to every interface being completely independent, and
> once that happens we lose a tremendous potential for Stan.

I think the key thing is that the language stays the same
and that the interfaces look relatively natural from their
host languages and that we use common computations in Stan C++
for all the heavy lifting.

Everything else is a detail.

- Bob

[Michael Betancourt]

So what exactly did you want me to look at?

The interface_callback/services pattern, in particular

the use of base_writer derivations everywhere and

how command.hpp currently looks.

I think the key thing is that the language stays the same
and that the interfaces look relatively natural from their
host languages and that we use common computations in Stan C++
for all the heavy lifting.

Everything else is a detail.

Except we are now being help up by details.

[Rob Goedman]

Hi,

Yes, Bob's summary is very useful.

As far as Julia and today’s Stan.jl are concerned, for Julia 0.3.x Allen’s earlier comment about following R and Python without C++ is correct. As alternatives folks can also fall back to RCall.jl and/or PyCall.jl.

Once Julia 0.4 + Cxx.jl are in general use (2nd half of this year?) I’ll reconsider if 1) I am convinced I don’t end up in the C++/Rtools version struggles across the 3 platforms and 2) there are clear advantages in using the C++ APIs.

As far as item 1) is concerned, my goal remains to install Stan through build.jl when installing Stan.jl. As far as item 2) is concerned, I do not consider ’shelling out’ CmdStan in Julia’s case a mistake because it keeps using Make simple (and thus averts extra compilations), it gives parallel chains almost for free and as Bob also mentioned I’ve not seen a case yet where i/o was a real bottleneck.

Either way, I’m on board with a callback (API) approach. From the Julia user (interface) point of view, Stan.jl attempts to bring Stan to Julia, not vice versa. Like input data pre-processing, I probably prefer ‘common’ configuration (e.g. dictionary based) and posterior tools for Julia/JuliaStats (today’s Stan.jl is leaning heavily on Mamba.jl).

As I am going to work on on-line help stuff for Stan.jl next week I’ll certainly keep this discussion in mind.

Regards,
Rob

[Michael Betancourt]

Interface devs,

Is a writer method using Eigen like

Writer::operator()(std::string key, Eigen::MatrixXd& m)

okay or is it substantially easier to use vectors, a la

Writer::operator()(std::string key, std::vector<double> m, int n_rows, int n_cols)?

[Bob Carpenter]

Is there a a reason these aren't const references?
Are some of the values getting written to? Or are
they callbacks where the interface implements the
Writer interface to do writing elsewhere?

- Bob

[Michael Betancourt]

I’m being lazy — they’re all const references in the code.

[Allen B. Riddell]

vector (or double *) please! It's easy to work with vector in Rcpp and Cython.

[Michael Betancourt]

If I pass pure double* arrays would people be comfortable
working without all of the STL methods? Right now using
vectors costs a big copy which we won’t be able to remedy
until we upgrade to C++11 and get to use all the sweet rvalue
copies.

[Allen B. Riddell]

Double* is preferable. Double* is the internal representation used by
Python (NumPy) anyway and I suspect it's the same for R. So if you pass
a vector it will have to be converted to double* anyway.

In fact, double* is better because there's no uncertainty about memory
contiguity from the interface's side (IIRC).

Thanks,

Allen

[Michael Betancourt]

What about arrays of strings (such as parameter names)?
Right now I’m using a vector<string> but would a string* array be easier?

[Ido Rosen]

Two other contenders to throw in the mix: heap-allocated vectors (vector<T>*) or heap-allocated strings (string<T>*) to save the copy until C++11?    String ensures contiguity and can be accessed directly via string::data() when needed.  (std::vector is contiguous but not by spec, tbb::concurrent_vector for example is not contiguous, std::string is contiguous by spec & impl.)

Related, an issue that's bitten me with C++14: VLAs go away in C++14 (and maybe come back in C++17), so using C-style arrays may become problematic in the future in a C++14 move.  (double* d will be fine, but double d[] might not.)

[Allen B. Riddell]

vector<string> is fine and works best from the Cython/Python side. A
copy will be made to get the strings into Python but this doesn't seem
like much of problem. We never have truly giant sequences of strings to
worry about, do we?

[Jiqiang Guo]

On Mon, Mar 30, 2015 at 10:05:41AM -0400, Allen B. Riddell wrote:
> Double* is preferable. Double* is the internal representation used by
> Python (NumPy) anyway and I suspect it's the same for R. So if you pass
> a vector it will have to be converted to double* anyway.

I don't know the context here. R does not use double*. R has its own
vector.

Jiqiang

[Allen B. Riddell]

Hi Jiqiang,

I think the question is this: which is easier to convert into an Rcpp::RList?

- Eigen::MatrixXd
- vector<double>
- double*

My sense is that the answer is vector<double>, right? e.g.,

Rcpp::List(sample_recorder.values_.x().begin(), sample_recorder.values_.x().end())

Best,

Allen

[Bob Carpenter]

I believe std::vector is required to be contiguous.
There's a C++11 method to get the data directory (rather
than as &x[0]), documented here:

http://www.cplusplus.com/reference/vector/vector/data/

Agreed that we should stay away from variable-length arrays.

I'd rather stay away with heap allocation that's not encapsulated
in initialization (the RIIA pattern), because I don't want to
have to track down memory leaks.

- Bob

[Allen B. Riddell]

If it's guaranteed to be contiguous then vector<double> is certainly
better.

Maybe I'm remembering something different -- perhaps it was
vector<vector<double>> that I couldn't seamlessly deal with in NumPy.

Best,

Allen

[Krzysztof Sakrejda]

On Wednesday, April 1, 2015 at 8:23:11 AM UTC-4, Allen B. Riddell wrote:
> If it's guaranteed to be contiguous then vector<double> is certainly
> better.
>
> Maybe I'm remembering something different -- perhaps it was
> vector<vector<double>> that I couldn't seamlessly deal with in NumPy.

vector became contiguous by standard with the 03 standard:

http://stackoverflow.com/questions/849168/are-stdvector-elements-guaranteed-to-be-contiguous

I'm pretty sure vector<vector<double>> would not be contiguous, but at that point using Eigen types seems easier. Rcpp (well, RcppEigen) will let rstan handle Eigen matrices easily, they are contiguous, and you can get a pointer to the Eigen matrix data array (and go the double* route) for Python.

Krzysztof

>
> Best,
>
> Allen

[Marcus Brubaker]

One advantage to using double* is that it can generally be passed without requiring a memory copy.  If a call requires std::vector<double> and you have anything *except* that type, you will need to allocate more memory and make a copy in order to make the call.  Using std::vector<double> for the API pretty much ensures that this happens for almost any call that interfaces between R/Matlab/Python/etc and Stan.  While not a huge performance hit for typical models now, if we are serious about scaling to big datasets, it's something to consider.

Cheers,

Marcus

--
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+unsubscribe@googlegroups.com.

[Krzysztof Sakrejda]

On Wednesday, April 1, 2015 at 1:23:54 PM UTC-4, Marcus Brubaker wrote:
> One advantage to using double* is that it can generally be passed without requiring a memory copy.  If a call requires std::vector<double> and you have anything *except* that type, you will need to allocate more memory and make a copy in order to make the call.  Using std::vector<double> for the API pretty much ensures that this happens for almost any call that interfaces between R/Matlab/Python/etc and Stan.  While not a huge performance hit for typical models now, if we are serious about scaling to big datasets, it's something to consider.

Turns out, like all good things, the move constructor for std::vector arrives with c++11. That does mean it will be here in a few months for Stan.

>
>
> Cheers,
> Marcus
>
>
>
> On Wed, Apr 1, 2015 at 10:23 AM Krzysztof Sakrejda <krzysztof...@gmail.com> wrote:
> On Wednesday, April 1, 2015 at 8:23:11 AM UTC-4, Allen B. Riddell wrote:
>
> > If it's guaranteed to be contiguous then vector<double> is certainly
>
> > better.
>
> >
>
> > Maybe I'm remembering something different -- perhaps it was
>
> > vector<vector<double>> that I couldn't seamlessly deal with in NumPy.
>
>
>
> vector became contiguous by standard with the 03 standard:
>
>
>
> http://stackoverflow.com/questions/849168/are-stdvector-elements-guaranteed-to-be-contiguous
>
>
>
> I'm pretty sure vector<vector<double>> would not be contiguous, but at that point using Eigen types seems easier. Rcpp (well, RcppEigen) will let rstan handle Eigen matrices easily, they are contiguous, and you can get a pointer to the Eigen matrix data array (and go the double* route) for Python.
>
>
>
> Krzysztof
>
>
>
> >
>
> > Best,
>
> >
>
> > Allen
>
>
>
> --
>
> 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.

[Marcus Brubaker]

From what I understand, I don't think the move constructor actually makes a difference here, as it only works if you're creating from an existing std::vector.  For instance, if I have an Eigen::Matrix<double,N,1> and I need to call something that takes std::vector<double>, there is no way to do that (as far as I know) without copying the entire contents into a newly allocated std::vector<double> of length N.  Please correct me if I'm wrong though.

Cheers,
Marcus

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

[Michael Betancourt]

Marcus, this is what the move constructors will accomplish.  Once we move to C++11

then we can use whichever type we want (double*, std::vector, Eigen::Matrix) and

translate between them with no copy or allocation cost.  Right now we can hack

around this by using &v[0] and assuming that the implementation is contiguous memory

of the right size without any sneaky alignment tricks.

[Krzysztof Sakrejda]

On Wednesday, April 1, 2015 at 2:45:01 PM UTC-4, Marcus Brubaker wrote:
> From what I understand, I don't think the move constructor actually makes a difference here, as it only works if you're creating from an existing std::vector.  For instance, if I have an Eigen::Matrix<double,N,1> and I need to call something that takes std::vector<double>, there is no way to do that (as far as I know) without copying the entire contents into a newly allocated std::vector<double> of length N.  Please correct me if I'm wrong though.

Nope, you're right, I didn't read carefully enough. It almost makes more sense to pass around Eigen types.

>
>
>
> Cheers,
> Marcus

[Rob Trangucci]

Apparently std::swap is as performant as move: http://stackoverflow.com/questions/12613428/stl-vector-moving-all-elements-of-a-vector

From what I gather, as of C++03 elements of the vector are stored contiguously, and are guaranteed to adhere to &v[n] = &v[0] + n, though this can change if elements are added to the vector as it may be reallocated (paraphrased from here: http://stackoverflow.com/questions/247738/is-it-safe-to-assume-that-stl-vector-storage-is-always-contiguous) 

Rob

[Allen B. Riddell]

PyStan currently doesn't wrap Eigen and there's no equivalent for
RcppEigen. If vector<double> is contiguous it can be converted into a
NumPy array without a copy.

Here's the C++ FAQ entry on the topic:

https://isocpp.org/wiki/faq/containers#vector-is-contiguous

It says the proper way to access the memory is via `vec.data()`.

Best,

Allen

[Michael Betancourt]

Allen,

But Eigen exposes the underlying data (which is contiguous) so you should be able to
do anything you can do with an array, no?

[Allen B. Riddell]

In theory yes (I think) -- but you have to tell Python or R that you're
dealing with an Eigen matrix that has a .data() method. Whereas
everything that supports wrapping C++ already knows about
vector<double>.

If the data in question is already flat, please use vector<double> and
have the Stan library do the mapping between the vector<double>'s memory
and Eigen MatrixXd's memory. That's in keeping with the way things
currently work in any event (e.g., var_context uses vector<double>).

Also, to the extent that we want to make the interface welcoming for
other software (Julia, etc.) it seems like sticking to just C++11 is a
good idea.

-a

[Marcus Brubaker]

Even if it knows about vector<double> it will still require a copy.  E.G., if you get a large array of data from Python as a NumPy array and need to pass it to Stan as a std::vector<double> there's no way to do that without a) an ugly hack potentially involving some weird class derived from std::vector<double> or b) copying the data.  It seems that the C++11 move constructor doesn't solve this issue.

Cheers,
Marcus

> > To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+unsubscribe@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 development mailing list" group.

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+unsubscribe@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 development mailing list" group.

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

[Michael Betancourt]

Firstly, we’re talking about moving data from Stan into the interfaces

so the important part is that we can expose the underlying contiguous

memory so that the interfaces can do what they need to do.

Going the other way we can play a few tricks — C++11 adds a bunch

of direct-data manipulation methods to the STL library so I don’t think

it will require trying to create an std::vector wrapper.

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

[Bob Carpenter]

This is, by the way, the reason I originally used vector<double> instead
of Eigen --- it's just more portable to other interfaces.

Yes, memory is contiguous for vector<double>, but not vector<vector<double> >.

Marcus is absolutely right that there's no way to create a vector<double> from
a double* without copying, even in C++11, because vector<double> wants to manage
the memory. That's the main advantage of using it over a double* (other than
that you also don't have to keep passing lengths around). See:

http://www.cplusplus.com/reference/vector/vector/vector/

Yes, we can copy a vector<double> to another without a deep copy using move
semantics in C++11. Or swap as Rob points out.

To avoid copying a vector<double>, can't we pass by reference, or by
pointer as a last resort? Or we can manage memory using a vector<double>
and pass a reference to the memory, but this is dangerous, as resizing
can cause the contents to move. That is, a vector<double> is guaranteed
to have contiguous memory at any given time, but that location is just
a handle, not a hard reference, so it may move.

We can't create a MatrixXd in Eigen without a copy, but we can use
the Map<MatrixXd> expression template to do the copy lazily or not
at all depending on how it's used.

- Bob

[Allen B. Riddell]

If it's only a vector<double> you can use Map<MatrixXd> to map it onto
an Eigen matrix without a copy. It's the vector<vector<double>> that's
the problem (and a big one, since that's how RStan and PyStan hold their
chains).

Is there any C++11 way to create a vector<vector<double>> that is
contiguous? Can we use an std::array or something? If there was a way to
do this cleanly, RStan and PyStan could use Stan's ess and rhat
routines.

Best,

Allen

> > an email to stan-dev+u...@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 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.
> >
> > --
> > 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.
> >
>
> --
> 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.

[Bob Carpenter]

> On Apr 2, 2015, at 12:04 PM, Allen B. Riddell <a...@ariddell.org> wrote:
>
> If it's only a vector<double> you can use Map<MatrixXd> to map it onto
> an Eigen matrix without a copy.

Map<MatrixXd> is an expression template that wraps raw data.
As soon as it's assigned, there's a copy implicated. So if you do this:

Map<MatrixXd> y(...);
MatrixXd m <- y;

you get a copy. If you just use it in something like an addition,
you can sometimes avoid the copy.

> It's the vector<vector<double>> that's
> the problem (and a big one, since that's how RStan and PyStan hold their
> chains).

OK --- that decisions on the RStan and PyStan side.

> Is there any C++11 way to create a vector<vector<double>> that is
> contiguous?

No. A vector<double> manages its memory independently on the heap.

You need to use an array. Or use the Boost Array package to do it
more generally.

> Can we use an std::array or something? If there was a way to
> do this cleanly, RStan and PyStan could use Stan's ess and rhat
> routines.

std::array is C++11, I'm afraid. It also doesn't do anything for you
other than template out the size so the size is static, which you
don't want.

I don't see anything for you other than to malloc your own 2D array
structure. I'd suggest doing it following std::vector so it uses
the RAII pattern in C++ and thus will ensure memory safety.

- Bob

[Allen B. Riddell]

On 04/02, Bob Carpenter wrote:
> > Can we use an std::array or something? If there was a way to
> > do this cleanly, RStan and PyStan could use Stan's ess and rhat
> > routines.
>
> std::array is C++11, I'm afraid. It also doesn't do anything for you
> other than template out the size so the size is static, which you
> don't want.
>
> I don't see anything for you other than to malloc your own 2D array
> structure. I'd suggest doing it following std::vector so it uses
> the RAII pattern in C++ and thus will ensure memory safety.
>

That sounds right. What if there are different chain lengths? Or can we
safely worry about that some other day?

[Bob Carpenter]

I don't think any of the code handles different chain lengths.

But I'm about to publish my preliminary spec for ragged arrays
this afternoon, and we can probably use whatever we do for that.

- Bob