Proposed Var Context Refactor

[Michael Betancourt]

Will everyone take some time to look through the proposed refractor
on feature/refactor_io_reader_writer? There’s a lot more work to be
done but the basic pattern is in place and I don’t want to implement
everything only to have the pattern change later.

The biggest change is that src/stan/io/writer.hpp and src/stan/io/reader.hpp
have been refactored into src/stan/io/transforms, with one file per
transform. This not only cleans up the model code, simplifies testing, etc,
but admits new methods in var_context,

void var_context::get_unconstrained_vector(std::vector<double>& parameters,
const std::string& name,
const std::string& scope,
size_t N,
stan::io::vector_transform& transform) {
vector_t constr_val = vector(name, N);
transform.unconstrain(constr_val, parameters);
}

The new var_context (and its derivations) can be found in
src/stan/interface/var_context. Note also the structured accessors,

double var_context::scalar(const std::string& name)
vector_t var_context::vector(const std::string& name, size_t N)
row_vector_t var_context::row_vector(const std::string& name, size_t N)
matrix_t var_context::matrix(const std::string& name, size_t n_rows, size_t n_cols)

With these changes the model code becomes substantially cleaner.
For example, in the constructor something like

context__.validate_dims("data initialization", "a", "matrix_d", context__.to_vec(5,5));
validate_non_negative_index("a", "5", 5);
validate_non_negative_index("a", "5", 5);
a = matrix_d(5,5);
vals_r__ = context__.vals_r("a");
pos__ = 0;
size_t a_k_mat_lim__ = 5;
for (size_t n_mat__ = 0; n_mat__ < a_k_mat_lim__; ++n_mat__) {
for (size_t m_mat__ = 0; m_mat__ < a_k_mat_lim__; ++m_mat__) {
a(m_mat__,n_mat__) = vals_r__[pos__++];
}
}

stan::math::check_cov_matrix(function__,"a",a);

becomes

validate_non_negative_index("a", "5", 5);
validate_non_negative_index("a", "5", 5);
a = context__.matrix(“a”, 5, 5);
stan::math::check_cov_matrix(function__,"a",a);

Similarly,

if (!(context__.contains_r("beta")))
throw std::runtime_error("variable beta missing");
vals_r__ = context__.vals_r("beta");
pos__ = 0U;
context__.validate_dims("initialization", "beta", "vector_d", context__.to_vec(5));
vector_d beta(5);
for (int j1__ = 0U; j1__ < 5; ++j1__)
beta(j1__) = vals_r__[pos__++];
try {
writer__.ordered_unconstrain(beta);
} catch (const std::exception& e) {
throw std::runtime_error(std::string("Error transforming variable beta: ") + e.what());
}

becomes

try {
context__.get_unconstrained_vector(params_r, “beta”, “initialization”, 5, stan::io::vector_ordered(5));
} catch(const std::exception& e) {
throw std::runtime_error(e.what());
}

We can also change the log_prob code and move from using std::vectors in all of the model
methods to Eigen::Vectors is desired.

It’d be nice to settle on this pattern as soon as possible so I can start the onerous task
of implementing everything. If you don’t like it then suggest an alternative — we’re already
way behind on the refactor. Also, I’ll be at a conference tomorrow so I won’t be online but
I should have plenty of email access.

[Bob Carpenter]

Some comments in line.

We could also do a hangout to discuss more.

> On Jun 29, 2015, at 11:13 AM, Michael Betancourt <betan...@gmail.com> wrote:
>
> Will everyone take some time to look through the proposed refractor
> on feature/refactor_io_reader_writer? There’s a lot more work to be
> done but the basic pattern is in place and I don’t want to implement
> everything only to have the pattern change later.
>
> The biggest change is that src/stan/io/writer.hpp and src/stan/io/reader.hpp
> have been refactored into src/stan/io/transforms, with one file per
> transform.

I think that's a very good idea.

> This not only cleans up the model code, simplifies testing, etc,
> but admits new methods in var_context,
>
> void var_context::get_unconstrained_vector(std::vector<double>& parameters,
> const std::string& name,
> const std::string& scope,
> size_t N,
> stan::io::vector_transform& transform) {
> vector_t constr_val = vector(name, N);
> transform.unconstrain(constr_val, parameters);
> }

I'm afraid I didn't understand the intent of this method.

* The standard convention is to make the return type non-constant args
last in argument lists (we haven't been so good about that, but it's
another thing I try to do and Google code urges us to do).

* What's the scope argument? It's not used in the implementation.

* (See below --- I rethought this after writing, but want to throw it
out there anyway.)

Why is an unconstrain transform an argument rather than being
built into the function named get_unconstrained_vector. Is the idea
that you might pass in either logit + bounds or log + location + sign?
I'd prefer just having mulitple functions exposed to the API client
so that the creating of the transform and getting of the unconstrained vec
was one call:

void get_lb_unconstrained(
const string& name,
double lb,
size_t n,
vector<double>& unconstrained_vec);

And then a whole bunch more. They could call the more generic
thing above if there's a reason to have the transforms not be
static.

> The new var_context (and its derivations) can be found in
> src/stan/interface/var_context. Note also the structured accessors,
>
> double var_context::scalar(const std::string& name)
> vector_t var_context::vector(const std::string& name, size_t N)
> row_vector_t var_context::row_vector(const std::string& name, size_t N)
> matrix_t var_context::matrix(const std::string& name, size_t n_rows, size_t n_cols)

Right --- this is more what I was thinking.
But now I realize we're looking at

no bound
lower bound
upper bound
lower & upper bound

crossed with:

real
vector
row vector
matrix

for each function. I'd just built those up over time in a really
awful way in the generated code itself.

OK, so now I'm not so sure about my suggestion above. But then
I'm also not so sure on how far to go with traits to reduce the
total number of functions (which all else being equal is a good thing).
Or how to encapsulate --- do I pass the matrix readers an unconstrained
reader?

> With these changes the model code becomes substantially cleaner.
> For example, in the constructor something like
>
> context__.validate_dims("data initialization", "a", "matrix_d", context__.to_vec(5,5));
> validate_non_negative_index("a", "5", 5);
> validate_non_negative_index("a", "5", 5);
> a = matrix_d(5,5);
> vals_r__ = context__.vals_r("a");
> pos__ = 0;
> size_t a_k_mat_lim__ = 5;
> for (size_t n_mat__ = 0; n_mat__ < a_k_mat_lim__; ++n_mat__) {
> for (size_t m_mat__ = 0; m_mat__ < a_k_mat_lim__; ++m_mat__) {
> a(m_mat__,n_mat__) = vals_r__[pos__++];
> }
> }
>
> stan::math::check_cov_matrix(function__,"a",a);
>
> becomes
>
> validate_non_negative_index("a", "5", 5);
> validate_non_negative_index("a", "5", 5);
> a = context__.matrix(“a”, 5, 5);
> stan::math::check_cov_matrix(function__,"a",a);

This will be awesome and very easy to change in the generator.
And thanks for working out this detail so I don't have to reason
through it myself.

> Similarly,
>
> if (!(context__.contains_r("beta")))
> throw std::runtime_error("variable beta missing");
> vals_r__ = context__.vals_r("beta");
> pos__ = 0U;
> context__.validate_dims("initialization", "beta", "vector_d", context__.to_vec(5));
> vector_d beta(5);
> for (int j1__ = 0U; j1__ < 5; ++j1__)
> beta(j1__) = vals_r__[pos__++];
> try {
> writer__.ordered_unconstrain(beta);
> } catch (const std::exception& e) {
> throw std::runtime_error(std::string("Error transforming variable beta: ") + e.what());
> }
>
> becomes
>
> try {
> context__.get_unconstrained_vector(params_r, “beta”, “initialization”, 5, stan::io::vector_ordered(5));
> } catch(const std::exception& e) {
> throw std::runtime_error(e.what());
> }

* I'm missing where the assignment to the local variable is, so I think
that should be:

beta = context__.get_unconstrained_vector(...);

or

context__.get_unconstrained_vector(..., beta);

* Is the "initialization" just for error messages? If so, we can perhaps
try/catch around the whole block to add a location if it's helpful so
we don't have to clutter all these up --- though I guess I'm not too bothered
by this

* Whatever we do, we shouldn't have those char* converted to std::string
on the fly if we can help it. We found that to have a pretty high overhead
in tight interfaces (I'm suring it's doing it now, but it'd be nice to
get rid of).

* I think we can go one step further and roll the exception
throwing into the get_unconstrained_vector operation itself

* Having to redundantly specify the 5 is a place programmers
(like me during code gen) can make errors.

Is there a way to make vector_ordered a simple transform that
takes size arguments? The issue is then what to do with transforms
that don't require sizes?

I can see this going by way of templates, with:

context__.get_unconstrained_vector<vector_ordered>(params_r, "beta", 5);

and that's implemented by calling a static function

vector_ordered::transform(5);

But then I'm not quite sure on how to make that more general to deal with
all the permutations:

context__.get_unconstrained_matrix<lower_bound>(params_r, "Sigma", 0, 5, 7);

The lower bound of 0 needs to get to the lower_bound constraint and the 5 x 7
dims needs to be handled by get_unconstrained_matrix.

I'm back to thinking just blowing it all out wouldn't be such a bad idea.

> We can also change the log_prob code and move from using std::vectors in all of the model
> methods to Eigen::Vectors is desired.

Didn't that already happen? Or are we still wrapping? We could perhaps provide
both if we could deal with the abstraction properly --- they both give us size() and
operator[].

> It’d be nice to settle on this pattern as soon as possible so I can start the onerous task
> of implementing everything.

Agreed. If you want help, I can implement some of it.

> If you don’t like it then suggest an alternative — we’re already
> way behind on the refactor. Also, I’ll be at a conference tomorrow so I won’t be online but
> I should have plenty of email access.

This should keep you busy during boring talks.

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

[Michael Betancourt]

I'm afraid I didn't understand the intent of this method.  

It allows this,

try {
 context__.get_unconstrained_vector(params_r, “beta”, “initialization”, 5, stan::io::vector_ordered(5));
} catch(const std::exception& e) {
 throw std::runtime_error(e.what());
}

In particular, it makes partial random inits trivial.  See

src/stan/interface/var_context/random_var_context.hpp

* The standard convention is to make the return type non-constant args
 last in argument lists (we haven't been so good about that, but it's
 another thing I try to do and Google code urges us to do).

We can fix order of headers and args later after the signatures

have been finalized.

* What's the scope argument?  It's not used in the implementation.

It’s used to construct exceptions — see below about cleaning

all of that up.

* (See below --- I rethought this after writing, but want to throw it
  out there anyway.)

 Why is an unconstrain transform an argument rather than being
 built into the function named get_unconstrained_vector.  Is the idea
 that you might pass in either logit + bounds or log + location + sign?  
 I'd prefer just having mulitple functions exposed to the API client
 so that the creating of the transform and getting of the unconstrained vec
 was one call:

   void get_lb_unconstrained(
     const string& name,
     double lb, 
     size_t n,
     vector<double>& unconstrained_vec);

 And then a whole bunch more.  They could call the more generic
 thing above if there's a reason to have the transforms not be
 static.

This couples the transforms and the var_context — every time you create

a new transform it has to be reimplemented (or at least wrapped) in the

var_context.  We have way too much duplication like this is in the code

already and it makes adding new features a pain in the ass.

OK, so now I'm not so sure about my suggestion above.  But then
I'm also not so sure on how far to go with traits to reduce the
total number of functions (which all else being equal is a good thing).
Or how to encapsulate --- do I pass the matrix readers an unconstrained
reader?  

What’s a matrix reader? 

try {
context__.get_unconstrained_vector(params_r, “beta”, “initialization”, 5, stan::io::vector_ordered(5));
} catch(const std::exception& e) {
throw std::runtime_error(e.what());
}

* I'm missing where the assignment to the local variable is, so I think
 that should be:

    beta = context__.get_unconstrained_vector(...);

  or

    context__.get_unconstrained_vector(..., beta);

There is no more local variable.  get_unconstrained grabs the constrained

values, unconstrained them, and then appends the unconstrained values

to the params_r vector.

* Is the "initialization" just for error messages?  If so, we can perhaps
 try/catch around the whole block to add a location if it's helpful so
 we don't have to clutter all these up --- though I guess I'm not too bothered
 by this

We have to be careful with what throws where, but, yes, this was my thinking

and would be implemented as the refactor is implemented in earnest.

* Whatever we do, we shouldn't have those char* converted to std::string
 on the fly if we can help it.  We found that to have a pretty high overhead
 in tight interfaces (I'm suring it's doing it now, but it'd be nice to
 get rid of).

Easy fix once we finalize the signatures.

* I think we can go one step further and roll the exception
 throwing into the get_unconstrained_vector operation itself

As I warned above this will be tricky as the transforms call free/constrain

functions that themselves can throw exceptions.  And we have to be careful

to validate constraints but without validating them too often.

* Having to redundantly specify the 5 is a place programmers
 (like me during code gen) can make errors.

This is all over the current code generator.  Not sure how to simplify it.

 Is there a way to make vector_ordered a simple transform that
 takes size arguments?  The issue is then what to do with transforms
 that don't require sizes?

 I can see this going by way of templates, with:

   context__.get_unconstrained_vector<vector_ordered>(params_r, "beta", 5);

 and that's implemented by calling a static function

   vector_ordered::transform(5);

 But then I'm not quite sure on how to make that more general to deal with
 all the permutations:

   context__.get_unconstrained_matrix<lower_bound>(params_r, "Sigma", 0, 5, 7);

 The lower bound of 0 needs to get to the lower_bound constraint and the 5 x 7
 dims needs to be handled by get_unconstrained_matrix.

 I'm back to thinking just blowing it all out wouldn't be such a bad idea.

See the current transform implementations — they’re stateful classes to store things

like bounds and sizes.  That way everything can be wrapped in a single 

get_unconstrained function.  Because bounds can be dynamic I don’t think there’s

any other way.

We can also change the log_prob code and move from using std::vectors in all of the model
methods to Eigen::Vectors is desired.

Didn't that already happen?  Or are we still wrapping?  We could perhaps provide
both if we could deal with the abstraction properly --- they both give us size() and
operator[].  

Nope, they wrap.  

[Bob Carpenter]

I had completely confused which direction you were
going here.

Let me take a shot at documenting the intent now. It
is a var_context method, right?

/**
* Reads an N-vector of the specified name out of this variable
* context and appends it to the end of the parameters vector,
* after applying the specified transform.
*
* The scope variable just indicates context for constructing an
* error message. There will be an error if the variable in
* this var_context is not of size N (or if the transformed value is
* not of size N?).
*/
void var_context::get_unconstrained_vector(

std::vector<double>& parameters,
const std::string& name,
const std::string& scope,
size_t N,

stan::io::vector_transform& transform);

As you can see, I'm still confused about specifying the sizes twice.
Is the N for the number of constrained or unconstrained values? Same
question for the argument to transform.

This is why I find it easier if the N is from the consumer's point of view.
That is, the model wants an N-simplex, so N is the number of dimensions, which
corresponds to (N - 1) unconstrained values. I don't then know which you
want to specify to the transform.

> On Jun 29, 2015, at 12:45 PM, Michael Betancourt <betan...@gmail.com> wrote:
>
> This couples the transforms and the var_context — every time you create
> a new transform it has to be reimplemented (or at least wrapped) in the
> var_context. We have way too much duplication like this is in the code
> already and it makes adding new features a pain in the ass.

I wouldn't say it's so much duplication as it is composing existing functions
to hide the internals to simplify the external interfaces. But I agree
it's a tradeoff, and as I said, I can live with the redundancy and cut-and-paste
(or cut-and-paste and resize) of sizes --- we'll just have to be very careful in
testing it all.

By "matrix reader", I just meant whatever did the same thing for matrices,
probably named "get_unconstrained_matrix". And there, you have more parameters
to play with that have to be synchronized across the equivalent of size N and
the equivalent of the arguments to the transform that might get messed up.

> As I warned above this will be tricky as the transforms call free/constrain
> functions that themselves can throw exceptions. And we have to be careful
> to validate constraints but without validating them too often.

I had the opposite concern --- I don't want to have to generate try/catch
any more in our model code --- it's just more redundancy that should be
rolled into the transforms (or even below if the transforms throw). So I
think we're on the same page here.

>> Didn't that already happen? Or are we still wrapping? We could perhaps provide
>> both if we could deal with the abstraction properly --- they both give us size() and
>> operator[].
>
> Nope, they wrap.

Can we resize an Eigen vector? Or do we just reallocate and copy? If so,
we won't be able to work efficiently directly with an Eigen vector unless we
knew how big it was and kept an offset along with all the rest of this.

- Bob

[Michael Betancourt]

>
> As you can see, I'm still confused about specifying the sizes twice.
> Is the N for the number of constrained or unconstrained values? Same
> question for the argument to transform.
>
> This is why I find it easier if the N is from the consumer's point of view.
> That is, the model wants an N-simplex, so N is the number of dimensions, which
> corresponds to (N - 1) unconstrained values. I don't then know which you
> want to specify to the transform.

Both refer to the constrained (hence user-facing) size. The double
specification should be possible to eliminate given how the transformations
are currently implemented — in a previous iteration all transformations
derived from a single base class which only had one constrained and
unconstrained size, hence the function needed additional arguments
to decompose the size into rows and columns.

> I wouldn't say it's so much duplication as it is composing existing functions
> to hide the internals to simplify the external interfaces. But I agree
> it's a tradeoff, and as I said, I can live with the redundancy and cut-and-paste
> (or cut-and-paste and resize) of sizes --- we'll just have to be very careful in
> testing it all.
>
> By "matrix reader", I just meant whatever did the same thing for matrices,
> probably named "get_unconstrained_matrix". And there, you have more parameters
> to play with that have to be synchronized across the equivalent of size N and
> the equivalent of the arguments to the transform that might get messed up.

There is only one “get_unconstrained” function — it works or scalars, vectors,
matrices, and arrays.

If you want redundancy then you have to support it — every time we want a
new transformation you have to teach the dev how to plumb it all through.
This lack of modularity has been a bane to introducing new devs for a long
time and I am very much on the side of avoiding it in the future.

> I had the opposite concern --- I don't want to have to generate try/catch
> any more in our model code --- it's just more redundancy that should be
> rolled into the transforms (or even below if the transforms throw). So I
> think we're on the same page here.

The reason for try/catches in the model code is to identify what was calling
the functions. Personally I find it slightly cleaner to not pass “scope” strings
around and have informative catch/throws in the model code, but I am not
wedded to that.

> Can we resize an Eigen vector? Or do we just reallocate and copy? If so,
> we won't be able to work efficiently directly with an Eigen vector unless we
> knew how big it was and kept an offset along with all the rest of this.

I don’t see why. At not point do really really need to resize and refill an input
vector — it’s just as easy to check for consistent sizes and assign individual
elements.

[Bob Carpenter]

> On Jun 29, 2015, at 5:57 PM, Michael Betancourt <betan...@gmail.com> wrote:
>
>>
>> As you can see, I'm still confused about specifying the sizes twice.
>> Is the N for the number of constrained or unconstrained values? Same
>> question for the argument to transform.
>>
>> This is why I find it easier if the N is from the consumer's point of view.
>> That is, the model wants an N-simplex, so N is the number of dimensions, which
>> corresponds to (N - 1) unconstrained values. I don't then know which you
>> want to specify to the transform.
>
> Both refer to the constrained (hence user-facing) size.

Good --- that'll make it much much easier to keep in synch.

> The double
> specification should be possible to eliminate given how the transformations
> are currently implemented — in a previous iteration all transformations
> derived from a single base class which only had one constrained and
> unconstrained size, hence the function needed additional arguments
> to decompose the size into rows and columns.

That'd be great. I can live with some redundancy in spec, given
the huge reduction in wrapper functions.

>> I wouldn't say it's so much duplication as it is composing existing functions
>> to hide the internals to simplify the external interfaces. But I agree
>> it's a tradeoff, and as I said, I can live with the redundancy and cut-and-paste
>> (or cut-and-paste and resize) of sizes --- we'll just have to be very careful in
>> testing it all.
>>
>> By "matrix reader", I just meant whatever did the same thing for matrices,
>> probably named "get_unconstrained_matrix". And there, you have more parameters
>> to play with that have to be synchronized across the equivalent of size N and
>> the equivalent of the arguments to the transform that might get messed up.
>
> There is only one “get_unconstrained” function — it works or scalars, vectors,
> matrices, and arrays.

Nice.

> If you want redundancy then you have to support it — every time we want a
> new transformation you have to teach the dev how to plumb it all through.
> This lack of modularity has been a bane to introducing new devs for a long
> time and I am very much on the side of avoiding it in the future.

I really really don't want redundancy.

Are there more transformations we want?

I think the big thing we're going to need to add in the immediate future
are:

* ragged arrays

* sparse arrays

>> I had the opposite concern --- I don't want to have to generate try/catch
>> any more in our model code --- it's just more redundancy that should be
>> rolled into the transforms (or even below if the transforms throw). So I
>> think we're on the same page here.
>
> The reason for try/catches in the model code is to identify what was calling
> the functions. Personally I find it slightly cleaner to not pass “scope” strings
> around and have informative catch/throws in the model code, but I am not
> wedded to that.

Maybe the best thing to do is catch at a block level,
so it doesn't need to be in the generated code, but can
be in how it's used. Does that make sense? That's a bigger
change, but it might simplify things. I haven't thought it through.

>> Can we resize an Eigen vector? Or do we just reallocate and copy? If so,
>> we won't be able to work efficiently directly with an Eigen vector unless we
>> knew how big it was and kept an offset along with all the rest of this.
>
> I don’t see why. At not point do really really need to resize and refill an input
> vector — it’s just as easy to check for consistent sizes and assign individual
> elements.

There wasn't an argument for position --- where do the functions write into
the parameters vector or Eigen matrix? I'm imagining filling a vector with multiple values,
and with a std::vector, you can just keep push_back()-ing.

- Bob

[Michael Betancourt]

> Are there more transformations we want?
>
> I think the big thing we're going to need to add in the immediate future
> are:
>
> * ragged arrays
>
> * sparse arrays

I don’t think we’d ever want to map sparse arrays to dense arrays,
that would immediately compromise their benefit. But ultimately
I’m not sure what transformations we may come across in the future
and we should try to design the code to be as welcoming to new
research ideas as possible.

> Maybe the best thing to do is catch at a block level,
> so it doesn't need to be in the generated code, but can
> be in how it's used. Does that make sense? That's a bigger
> change, but it might simplify things. I haven't thought it through.

That would be easy. I just think it’s helpful to somehow append
the scope (data, init, log_prob) to the exception to help guide the
user to the problem.

> There wasn't an argument for position --- where do the functions write into
> the parameters vector or Eigen matrix? I'm imagining filling a vector with multiple values,
> and with a std::vector, you can just keep push_back()-ing.

Yeah, the proposal assumed vectors with push_backing, but
it’s easy to add a position argument provided that we decide
on it now rather than later.

[Bob Carpenter]

> On Jun 30, 2015, at 4:57 AM, Michael Betancourt <betan...@gmail.com> wrote:
>
>> Are there more transformations we want?
>>
>> I think the big thing we're going to need to add in the immediate future
>> are:
>>
>> * ragged arrays
>>
>> * sparse arrays
>
> I don’t think we’d ever want to map sparse arrays to dense arrays,
> that would immediately compromise their benefit.

I hadn't even thought about cross assignability, but I don't see
any reason to outlaw it other than that it's almost always a bad idea
in practice.

> But ultimately
> I’m not sure what transformations we may come across in the future
> and we should try to design the code to be as welcoming to new
> research ideas as possible.

It's more just that we might have a sparse array of parameters
with lower bounds of 0, for instance. This won't change the transforms,
just what they're applied to.

>> Maybe the best thing to do is catch at a block level,
>> so it doesn't need to be in the generated code, but can
>> be in how it's used. Does that make sense? That's a bigger
>> change, but it might simplify things. I haven't thought it through.
>
> That would be easy. I just think it’s helpful to somehow append
> the scope (data, init, log_prob) to the exception to help guide the
> user to the problem.

Agreed. I'm trying to update all of our error reporting to make
it more useful. The only biggy left is the way arrays are printed.

>> There wasn't an argument for position --- where do the functions write into
>> the parameters vector or Eigen matrix? I'm imagining filling a vector with multiple values,
>> and with a std::vector, you can just keep push_back()-ing.
>
> Yeah, the proposal assumed vectors with push_backing, but
> it’s easy to add a position argument provided that we decide
> on it now rather than later.

Yes it would.

- Bob