Data Access Layer Spec

[Michael Betancourt]

I’ve attached the basic spec for the data access layer (which will replace var_context),

refactored transforms, and how they will be used in the generated model code.  It should

be self consistent.  Any thoughts?

// DATA ACCESS LAYERS

// base_dal

vector_t get_vector(const std::string& name, size_t N) {

  if (!contains_r(name)) {

    std::stringstream msg;

    msg << "Variable " << name

        << " not found in variable context" << std::endl;

    throw std::runtime_error(msg.str());

  }

  std::vector<size_t> dims = dims_r(name);

  if (dims.size() != 1) {

    std::stringstream msg;

    msg << "Requested the vector " << name

        << ", but " << name

        << " is not defined as a vector in variable context"

        << std::endl;

    throw std::runtime_error(msg.str());

  }

  if (dims[0] != N) {

    std::stringstream msg;

    msg << "Requested the vector " << name

        << " with size " << N << ", but " << name

        << " is defined with size " << dims[0] << " in variable context"

        << std::endl;

    throw std::runtime_error(msg.str());

  }

  std::vector<double> vals = vals_r(name);

  vector_t output(N);

  for (vector_t_idx n = 0; n < N; ++n)

    output(n) = vals[n];

  return output;

}

// constr_dal

template <typename T>

void get_constrained_vector(const string& name,

                            const stan::io::transform& transform,

                            vector_t param) {

  param.resize(transform.size());

  param = get_vector(name, transform.size());

  transform.validate(param);

}

// unconstr_dal

template <typename T>

void get_unconstrained_vector(const string& name,

                              const stan::io::transform& transform,

                              vector_t::InnerIterator& it) {

  vector_t constr_val = get_vector(name, transform.size());

  transform.unconstrain(constr_val, it);

}

// EXAMPLE TRANSFORM

template <typename T, typename T_LB>

class vector_lower_bound: public vector_transform<T> {

private:

  T_LB lower_bound_;

public:

  vector_lower_bound(T_LB lower_bound, size_t N):

    lower_bound_(lower_bound), vector_transform<T>(N) {}

  size_t unconstrained_dim() { return N_; }

  void validate(vector_t constr_val) {

    stan::math::check_greater_or_equal("stan::io::vector_lower_bound",

                                       "Lower bound constraint", constr_val,

                                       lower_bound_);

  }

  template <typename Jacobian>

  vector_t constrain(const vector_t& unconstr_val, lp_accumulator<T_LB>& acc) {

    vector_t constr_val;

    for (idx_t n = 0; n < N_; ++n) {

      if (Jacobian) {

        T_LB lp;

        output(n) = stan::prob::lb_constrain(unconstr_val[n], lower_bound_, lp);

        acc.push_back(lp);

      } else {

        output(n) = stan::prob::lb_constrain(unconstr_val[n], lower_bound);

      }

    }

  }

  vector_t constrain(const vector_t& unconstr_val,

                     std::vector<double> constr_val) {

    for (idx_t n = 0; n < N_; ++n)

      constr_val.push_back(stan::prob::lb_constrain(unconstr_val[n],

                                                    lower_bound_));

  }

  void unconstrain(const vector_t& constr_val, vector_t::InnerIterator& it) {

    for (idx_t n = 0; n < N_; ++n)

      *(it++) = stan::prob::lb_free(constr_val[n], lower_bound_);

  }

};

// MODEL USE CASES

// constructor

model_name {

  ...

  vector_lb_transform<input_type, lb_type>

    param_name1_transform__(lb_value);

  data_var_context.get_constrained_vector("param_name1",

                                          param_name1_transform,

                                          param_name1__);

  ...

}

// init

void init(vector_t& unconstr_state) {

  vector_t::InnerIterator it(unconstr_state);

  ...

  vector_lb_transform<input_type, lb_type>

    param_name1_transform__(lb_value);

  get_unconstrained_vector(“param_name1", param_name1_transform, it);

  ...

}

// what write_array was

void constr_state(vector_t unconstr_state,

                  std::vector<double>& constr_state) {

  ...

  vector_lb_transform<input_type, lb_type>

    param_name1_transform__(lb_value);

  param_name1_transform__.constrain(unconstr_state, constr_state);

  ...

}

// log_prob

template <typename T>

T log_prob(vector_t& unconstr_state) {

  lp_accumulator<T> lp_acc__;

  ...

  vector_lb_transform<input_type, lb_type>

    param_name1_transform__(lb_value);

  vector_t param_name1__ =

    param_name1_transform__.constrain<Jacobian>(unconstr_state, lp_acc__);

  ...

}

[Bob Carpenter]

Thanks for writing this up. Overall it looks great.

I have a few comments inline.

> On Dec 16, 2015, at 5:49 PM, Michael Betancourt <betan...@gmail.com> wrote:
>
> I’ve attached the basic spec for the data access layer (which will replace var_context),
> refactored transforms, and how they will be used in the generated model code. It should
> be self consistent. Any thoughts?
>
> // DATA ACCESS LAYERS
>
> // base_dal
> vector_t get_vector(const std::string& name, size_t N) {
> if (!contains_r(name)) {
> std::stringstream msg;
> msg << "Variable " << name
> << " not found in variable context" << std::endl;
> throw std::runtime_error(msg.str());
> }

Should we include std::endl in the error messages or insert
from the catching context? I think we need Ben and Allen's
input w.r.t. RStan and PyStan.

> std::vector<size_t> dims = dims_r(name);

I like having the new data access layer knowing what the
dimensions are.

> if (dims.size() != 1) {
> std::stringstream msg;
> msg << "Requested the vector " << name
> << ", but " << name
> << " is not defined as a vector in variable context"
> << std::endl;
> throw std::runtime_error(msg.str());
> }
>
> if (dims[0] != N) {
> std::stringstream msg;
> msg << "Requested the vector " << name
> << " with size " << N << ", but " << name
> << " is defined with size " << dims[0] << " in variable context"
> << std::endl;
> throw std::runtime_error(msg.str());
> }
>
> std::vector<double> vals = vals_r(name);
> vector_t output(N);
> for (vector_t_idx n = 0; n < N; ++n)
> output(n) = vals[n];

I think we can probably get away without the extra
copy into vals here by just iterating from vals_r(name),
at least in some implementations --- don't know if it's
general enough, though.

> return output;
> }
>
> // constr_dal
> template <typename T>
> void get_constrained_vector(const string& name,
> const stan::io::transform& transform,
> vector_t param) {
> param.resize(transform.size());
> param = get_vector(name, transform.size());
> transform.validate(param);
> }

What is the template parameter T for here? It's not being used.

> // unconstr_dal
> template <typename T>
> void get_unconstrained_vector(const string& name,
> const stan::io::transform& transform,
> vector_t::InnerIterator& it) {
> vector_t constr_val = get_vector(name, transform.size());
> transform.unconstrain(constr_val, it);
> }

Shouldn't that iterator just be templated out? And again, what is T?

>
> // EXAMPLE TRANSFORM
> template <typename T, typename T_LB>
> class vector_lower_bound: public vector_transform<T> {
> private:
> T_LB lower_bound_;
>
> public:
> vector_lower_bound(T_LB lower_bound, size_t N):
> lower_bound_(lower_bound), vector_transform<T>(N) {}
>
> size_t unconstrained_dim() { return N_; }
>
> void validate(vector_t constr_val) {
> stan::math::check_greater_or_equal("stan::io::vector_lower_bound",
> "Lower bound constraint", constr_val,
> lower_bound_);
> }
>
> template <typename Jacobian>
> vector_t constrain(const vector_t& unconstr_val, lp_accumulator<T_LB>& acc) {
> vector_t constr_val;
> for (idx_t n = 0; n < N_; ++n) {
> if (Jacobian) {
> T_LB lp;
> output(n) = stan::prob::lb_constrain(unconstr_val[n], lower_bound_, lp);
> acc.push_back(lp);
> } else {
> output(n) = stan::prob::lb_constrain(unconstr_val[n], lower_bound);
> }
> }
> }

I think it'd be clearer with "if (Jacobian)" on the outside of the loop.

This is problematic with lp declared on the inside --- it won't work
for a double value because the init isn't defined. I think we should do
what you said in person, which is pass the accumulator in. The way they're
written now, it could be done as:

if (Jacobian) {
T_LB lp(0);

for (idx_t n = 0; n < N_; ++n)

output(n) = stan::prob::lb_constrain(unconstr_val[n], lower_bound_, lp);
acc.push_back(lp);
} else {

for (idx_t n = 0; n < N_; ++n)

output(n) = stan::prob::lb_constrain(unconstr_val[n], lower_bound_);
}

and you only get one init for lp, but not the benefit of turning it
into a sum. Alternatively, we could have the ones with lp arguments
take in a reference and set the value, then push back on the outside.

The type of vector_t and T_LB can be different? It might be a double
value for the accumulator. Can we then drop the Jacobian if it's
constant? Or do we need to worry about getting the full density either
way?

I don't think there's a stan::prob namespace any more. It's
all just stan::math (and if it's not, it should be stan::math).

> vector_t constrain(const vector_t& unconstr_val,
> std::vector<double> constr_val) {
> for (idx_t n = 0; n < N_; ++n)
> constr_val.push_back(stan::prob::lb_constrain(unconstr_val[n],
> lower_bound_));
> }

Whenever you have a string of push_backs and know the size, it
helps to use

constr_val.reserve(N_);

up front.

>
> void unconstrain(const vector_t& constr_val, vector_t::InnerIterator& it) {
> for (idx_t n = 0; n < N_; ++n)
> *(it++) = stan::prob::lb_free(constr_val[n], lower_bound_);
> }
> };

For vector_transform<T>, are there any virtual functions? I
didn't see the base class spec.

>
> // MODEL USE CASES

Awesome!

I was thinking we might be able to make the transforms class
member variables, but that won't work when the bounds are
parameters.

> // log_prob
> template <typename T>
> T log_prob(vector_t& unconstr_state) {
> lp_accumulator<T> lp_acc__;
>
> ...
> vector_lb_transform<input_type, lb_type>
> param_name1_transform__(lb_value);
> vector_t param_name1__ =
> param_name1_transform__.constrain<Jacobian>(unconstr_state, lp_acc__);
> ...
> }

Overall, this looks great. It's going to hugely simplify the
generated model code and reduce a ton of code duplication.

- Bob

[Michael Betancourt]

What is the template parameter T for here?  It's not being used.

Sorry, it’s not really a template on the function but the class.

In particular, vector_t is a typedef that uses T.

typedef Eigen::Matrix<T, Eigen::Dynamic, 1> vector_t;

This is problematic with lp declared on the inside --- it won't work
for a double value because the init isn't defined.  I think we should do
what you said in person, which is pass the accumulator in.  The way they're
written now, it could be done as:

 if (Jacobian) {
   T_LB lp(0);
   for (idx_t n = 0; n < N_; ++n)
     output(n) = stan::prob::lb_constrain(unconstr_val[n], lower_bound_, lp);
   acc.push_back(lp);
 } else {
   for (idx_t n = 0; n < N_; ++n)
     output(n) = stan::prob::lb_constrain(unconstr_val[n], lower_bound_);
 }

and you only get one init for lp, but not the benefit of turning it
into a sum.  Alternatively, we could have the ones with lp arguments
take in a reference and set the value, then push back on the outside.

Why can’t lp(0) be used on the inside?

The type of vector_t and T_LB can be different?  It might be a double
value for the accumulator.  Can we then drop the Jacobian if it's
constant?  Or do we need to worry about getting the full density either
way?

vector_t and T_LB can indeed be different as necessary.  I haven’t

thought about dropping any constants, but in any case I think that

logic should be pushed back into the math transforms.

For vector_transform<T>, are there any virtual functions?  I
didn't see the base class spec.

Yeah, all of the functions mentioned above are defined as pure virtual

in a base class (or they will be).

[Bob Carpenter]

> On Dec 20, 2015, at 6:32 PM, Michael Betancourt <betan...@gmail.com> wrote:
>
>
>> What is the template parameter T for here? It's not being used.
>
> Sorry, it’s not really a template on the function but the class.
> In particular, vector_t is a typedef that uses T.
>
> typedef Eigen::Matrix<T, Eigen::Dynamic, 1> vector_t;

Got it.

>> This is problematic with lp declared on the inside --- it won't work
>> for a double value because the init isn't defined. I think we should do
>> what you said in person, which is pass the accumulator in. The way they're
>> written now, it could be done as:
>>
>> if (Jacobian) {
>> T_LB lp(0);
>> for (idx_t n = 0; n < N_; ++n)
>> output(n) = stan::prob::lb_constrain(unconstr_val[n], lower_bound_, lp);
>> acc.push_back(lp);
>> } else {
>> for (idx_t n = 0; n < N_; ++n)
>> output(n) = stan::prob::lb_constrain(unconstr_val[n], lower_bound_);
>> }
>>
>> and you only get one init for lp, but not the benefit of turning it
>> into a sum. Alternatively, we could have the ones with lp arguments
>> take in a reference and set the value, then push back on the outside.
>
> Why can’t lp(0) be used on the inside?

It can be, but it just allocates a new variable on the stack each time.
Passing in the accumulator would be most efficient. So on the inside,
you get

sum(0 + j[1], 0 + j[2], ..., 0 + j[n])

whereas doing it on the outside you get

0 + j[1] + j[2] + ... + j[n]

and with an accumulator you get the ideal

sum(j[1], j[2], ..., j[n])

and if you pass in the global accumulator, it's a bit better still
because those j[1]...j[n] just go in the global accumulator's sum.

>> The type of vector_t and T_LB can be different? It might be a double
>> value for the accumulator. Can we then drop the Jacobian if it's
>> constant? Or do we need to worry about getting the full density either
>> way?
>
> vector_t and T_LB can indeed be different as necessary. I haven’t
> thought about dropping any constants, but in any case I think that
> logic should be pushed back into the math transforms.

Agreed.

>> For vector_transform<T>, are there any virtual functions? I
>> didn't see the base class spec.
>
> Yeah, all of the functions mentioned above are defined as pure virtual
> in a base class (or they will be).

Right --- if T is the class template parameter type, I think that's
the only option. The alternative, that wouldn't use virtual functions, is
to add another template parameter for the transform. I'd go with that,
because it'll be most efficient and given that it's not some complex
template metaprogram, won't adversely affect compile time (no recursion,
just simple matching).

- Bob

[Michael Betancourt]

It can be, but it just allocates a new variable on the stack each time.
Passing in the accumulator would be most efficient.  So on the inside,
you get

 sum(0 + j[1], 0 + j[2], ..., 0 + j[n])

whereas doing it on the outside you get

 0 + j[1] + j[2] + ... + j[n]

and with an accumulator you get the ideal

 sum(j[1], j[2], ..., j[n])

and if you pass in the global accumulator, it's a bit better still
because those j[1]...j[n] just go in the global accumulator's sum.

I’m still confused.  Firstly, you can’t pass in the accumulator because

the math transforms require that you pass in a var to accept the

Jacobian increment.  Unless the math transforms get modified

we’re limited to passing in something local.

My understanding is that on the inside the temporary doesn’t get

connected to lp (i.e. it’s an orphan node in the expression graph)

until it gets pushed back onto the accumulator, so it should achieve

the ideal reduction.  

Right --- if T is the class template parameter type, I think that's
the only option.  The alternative, that wouldn't use virtual functions, is
to add another template parameter for the transform.  I'd go with that,
because it'll be most efficient and given that it's not some complex
template metaprogram, won't adversely affect compile time (no recursion,
just simple matching).

But the transforms are just implementations of this base class.

They are not defined in their own classes so there’s nothing

to template on.

[Bob Carpenter]

> On Dec 20, 2015, at 11:07 PM, Michael Betancourt <betan...@gmail.com> wrote:
>
>> It can be, but it just allocates a new variable on the stack each time.
>> Passing in the accumulator would be most efficient. So on the inside,
>> you get
>>
>> sum(0 + j[1], 0 + j[2], ..., 0 + j[n])
>>
>> whereas doing it on the outside you get
>>
>> 0 + j[1] + j[2] + ... + j[n]
>>
>> and with an accumulator you get the ideal
>>
>> sum(j[1], j[2], ..., j[n])
>>
>> and if you pass in the global accumulator, it's a bit better still
>> because those j[1]...j[n] just go in the global accumulator's sum.
>
> I’m still confused. Firstly, you can’t pass in the accumulator because
> the math transforms require that you pass in a var to accept the
> Jacobian increment. Unless the math transforms get modified
> we’re limited to passing in something local.

Right.

> My understanding is that on the inside the temporary doesn’t get
> connected to lp (i.e. it’s an orphan node in the expression graph)
> until it gets pushed back onto the accumulator, so it should achieve
> the ideal reduction.

The temp is incremented by the transforms, not set. So it
needs to be initialized to zero.

>> Right --- if T is the class template parameter type, I think that's
>> the only option. The alternative, that wouldn't use virtual functions, is
>> to add another template parameter for the transform. I'd go with that,
>> because it'll be most efficient and given that it's not some complex
>> template metaprogram, won't adversely affect compile time (no recursion,
>> just simple matching).
>
> But the transforms are just implementations of this base class.
> They are not defined in their own classes so there’s nothing
> to template on.

If you get rid of the base class, implement what you're thinking
of as specializations as independent classes and template, you lose
the virtual function call. Probably not that big of a deal.

The main reason not to do that is if there's implementation details
in the base class. That'd reduce code duplication at the cost of
virtual function call overheads.

- Bob

[Michael Betancourt]

If you get rid of the base class, implement what you're thinking
of as specializations as independent classes and template, you lose
the virtual function call.  Probably not that big of a deal.

The main reason not to do that is if there's implementation details
in the base class.  That'd reduce code duplication at the cost of
virtual function call overheads.

There is a small amount of implementation in the base class, although

I prefer it more for defining the interface.  We could template out the

DAL methods on the transform types but because they’re used only

once I don’t think it’ll have any affect on performance.

[Bob Carpenter]

> On Dec 21, 2015, at 10:39 AM, Michael Betancourt <betan...@gmail.com> wrote:
>
>
>> If you get rid of the base class, implement what you're thinking
>> of as specializations as independent classes and template, you lose
>> the virtual function call. Probably not that big of a deal.
>>
>> The main reason not to do that is if there's implementation details
>> in the base class. That'd reduce code duplication at the cost of
>> virtual function call overheads.
>
> There is a small amount of implementation in the base class, although
> I prefer it more for defining the interface.

There's a tension in C++ between using templates and
meta-defined "concepts" (meta in the sense that they aren't
defined in the object language) vs. proper object-oriented
programming. I agree that it's easier to understand as
a developer if there's a base class that gets passed rather
than a templated concept. The tension is because the templated
approach is faster. The in-between "solution" is the CRTP, like
Eigen uses for its matrix classes, though we've also seen what
a pain that is in practice.

> We could template out the
> DAL methods on the transform types but because they’re used only
> once I don’t think it’ll have any affect on performance.

Right. The run-time issue is with the unconstraining
transforms. Anyway, this is something that'll be easy to modify
later without affecting the interface.

- Bob

[Michael Betancourt]

> There's a tension in C++ between using templates and
> meta-defined "concepts" (meta in the sense that they aren't
> defined in the object language) vs. proper object-oriented
> programming. I agree that it's easier to understand as
> a developer if there's a base class that gets passed rather
> than a templated concept. The tension is because the templated
> approach is faster. The in-between "solution" is the CRTP, like
> Eigen uses for its matrix classes, though we've also seen what
> a pain that is in practice.

Right, but you don’t _have_ to use the base class with virtual
functions. I just like having it there to define the concepts
concretely and in the code we can always use templates to
avoid the virtual overhead.

> Right. The run-time issue is with the unconstraining
> transforms. Anyway, this is something that'll be easy to modify
> later without affecting the interface.

Yeah, but those get code generated and we can use the top-level
implementation and never worry about virtuals. Virtuals get used
only in the DAL which is called once on construction and once on
initialization.

[Bob Carpenter]

> On Dec 21, 2015, at 12:35 PM, Michael Betancourt <betan...@gmail.com> wrote:
>
>> There's a tension in C++ between using templates and
>> meta-defined "concepts" (meta in the sense that they aren't
>> defined in the object language) vs. proper object-oriented
>> programming. I agree that it's easier to understand as
>> a developer if there's a base class that gets passed rather
>> than a templated concept. The tension is because the templated
>> approach is faster. The in-between "solution" is the CRTP, like
>> Eigen uses for its matrix classes, though we've also seen what
>> a pain that is in practice.
>
> Right, but you don’t _have_ to use the base class with virtual
> functions. I just like having it there to define the concepts
> concretely and in the code we can always use templates to
> avoid the virtual overhead.

Here's an example of how inheritance works in C++ with
virtual functions. Note that you only need to define the
base class method as virtual and then all the inherited class methods
will also be virtual. The behavior doesn't change if B::foo is
declared as virtual.

#include <iostream>

struct A {
virtual int foo() const { return 0; }
};

struct B : public A {
int foo() const { return 1; }
};

void printer(const A& a) {
std::cout << "printer " << a.foo() << std::endl;
}

template <typename T>
void printer2(const T& x) {
std::cout << "templated printer " << x.foo() << std::endl;
}

int main() {
A* a1 = new B();
std::cout << "a1->foo() " << a1->foo() << std::endl;

A a2 = B();
std::cout << "a2.foo() " << a2.foo() << std::endl;

std::cout << "B().foo() " << B().foo() << std::endl;

B b = B();
printer(b);

printer2(b);
}


The output is:

~/temp2$ clang++ foo.cpp

~/temp2$ ./a.out
a1->foo() 1
a2.foo() 0
B().foo() 1
printer 1
templated printer 1


So I don't think we can avoid the virtual call overhead by templating.

Now Java has a clever way of declaring a class final, but
C++ allows linking in new instances, so it can't decide statically
if a class is final (i.e., doesn't allow extensions). So it can't
know how that printer2() function dispatches at compile time.

- Bob

[Michael Betancourt]

When oh when can we move to 11? http://en.cppreference.com/w/cpp/language/final

I see what you mean, though.  We could always avoid making the methods

virtual in the first place and treat the base class as just a specification of the

concept.

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