more metaprogram unfolding

[Bob Carpenter]

Here's Matt's version of partials3. We can do some
more metaprogramming for efficiency.

template <typename T1, typename T2, typename T3>
class partials3_vari : public vari {
private:
const static size_t nvars_ = !stan::is_constant<T1>::value +
!stan::is_constant<T2>::value +
!stan::is_constant<T3>::value;
vari* operands_[nvars_];
double partials_[nvars_];

public:
partials3_vari(double value,
T1 operand1, double partial1,
T2 operand2, double partial2,
T3 operand3, double partial3)
: vari(value) {
size_t i = 0;
if (!stan::is_constant<T1>::value) {
operands_[i] = *(vari**)(&operand1);
partials_[i] = partial1;
i++;
}
if (!stan::is_constant<T2>::value) {
operands_[i] = *(vari**)(&operand2);
partials_[i] = partial2;
i++;
}
if (!stan::is_constant<T3>::value) {
operands_[i] = *(vari**)(&operand3);
partials_[i] = partial3;
}
}

void chain() {
for (size_t i = 0; i < nvars_; i++)
operands_[i]->adj_ += adj_ * partials_[i];
}
};

I've been reading more about metaprogramming in Vandevoorde
and Josuttis.

For chain(), we can unroll the loop statically
with a metaprogram:

void chain() {
chain_meta<nvars_>::chain(operands_,partials_,adj_);
}

First the recursive case:

template <int DIM>
class chain_meta {
static void chain(vari** operands, double* partials, double adj) {
(*operands)->adj_ += adj * (*partials);
chain_meta<nvars_-1>::chain(operands+1, partials+1, adj);
}
}

Then the base case specialization:

template <>
class chain_meta<0> {
static void chain(vari**,double*,double) { }
}

Assuming everything gets inlined and you have 3 args, that'll be
like:

operands_[0]->adj_ += adj_ * partials_[0];
operands_[1]->adj_ += adj_ * partials_[1];
operanads_[2]->adj_ += adj_ * partials_[2];

Alternatively, you could implement this statically w/o recursion for
three entries with:

void chain() {
if (nvars_ == 0) return;
operands_[0]->adj_ += adj_ * partials_[0];
if (nvars_ == 1) return;
operands_[1]->adj_ += adj_ * partials_[1];
if (nvars_ == 2) return;
operands_[2]->adj_ += adj_ * partials_[2];
}

We could do the same sort of thing with the constructor to
figure out the index in the array of each variable and thus bypass
the need for incrementing i.

All in all, though, I'm wondering how much clearer this all is
and if it'll be easier to understand and maintain than the
whole-bunch-of-functions approach.

- Bob

[Matt Hoffman]

I'm not sure that there's an efficiency advantage to this
approach—nvars_ is known at compile time, so shouldn't the compiler be
able to make its own decisions about whether or not to unroll that
loop?

Matt

[Bob Carpenter]

On 6/23/12 5:56 PM, Matt Hoffman wrote:
> I'm not sure that there's an efficiency advantage to this

> approach�nvars_ is known at compile time, so shouldn't the compiler be

> able to make its own decisions about whether or not to unroll that
> loop?

If you didn't know nvars_ statically, the metaprogram
unrolling wouldn't work.

I don't know how much loop unrolling the compiler will
do if the bounds are small and known. The metaprogram
involves relying on inlining.

I somehow doubt the compiler would unroll the i in the
constructor, though it's not impossible -- the info is
all there.

- Bob

[Matt Hoffman]

We're just talking about chain(), right? Compilers are pretty hip to the
for(size_t i = 0; i < constant; i++) {...}
construction, and my understanding is that gcc unrolls pretty
aggressively. Glancing at the object code in Instruments, it seems as
though there aren't any jump instructions, at least for the version
that I saw.

Matt

On Sat, Jun 23, 2012 at 6:26 PM, Bob Carpenter <ca...@alias-i.com> wrote:
>
>
> On 6/23/12 5:56 PM, Matt Hoffman wrote:
>>
>> I'm not sure that there's an efficiency advantage to this

>> approach—nvars_ is known at compile time, so shouldn't the compiler be

[Ben Goodrich]

On Saturday, June 23, 2012 5:23:57 PM UTC-4, Bob Carpenter wrote:

All in all, though, I'm wondering how much clearer this all is
and if it'll be easier to understand and maintain than the
whole-bunch-of-functions approach.

My $0.02 is that this partials stuff with metaprogramming is likely to be understood only by computer scientists with an above average understanding of C++. But, it seems that it would substantially lower the barriers to entry for a typical statistician who might be able to add a new density to stan. So, even if there is no net improvement in speed or no net improvement in maintainability, I think if there is a gross improvement on a dimension where a larger pool of people might be interested in contributing to stan, that is a good thing in the long term.

Ben

[Bob Carpenter]

On 6/23/12 7:57 PM, Ben Goodrich wrote:
> On Saturday, June 23, 2012 5:23:57 PM UTC-4, Bob Carpenter wrote:
>
> All in all, though, I'm wondering how much clearer this all is
> and if it'll be easier to understand and maintain than the
> whole-bunch-of-functions approach.
>
>
> My $0.02 is that this partials stuff with metaprogramming is likely to be understood only by computer scientists with an
> above average understanding of C++.

I'm afraid that's the ante at this point.

It still helps to make things clear. I'm still confused
every time I look at Eigen, but I think that's largely
because of their use of macros more than the C++ itself.

> But, it seems that it would substantially lower the barriers to entry for a typical
> statistician who might be able to add a new density to stan. So, even if there is no net improvement in speed or no net
> improvement in maintainability, I think if there is a gross improvement on a dimension where a larger pool of people
> might be interested in contributing to stan, that is a good thing in the long term.

This is mainly going to matter in the case where someone
wants to unfold the derivative calcs rather than using
auto-diff. It makes a big difference in speed.

At that point, we have a coding style in play that'll
make it pretty clear in a way rather like the Rcpp interface
to the R data structures:

1. pull values out in double-based scalars, vectors, etc.
using value_of() function.

2. compute log probability and derivative, dropping terms
based on include_summand<> metaprograms and dropping
gradient terms based on is_constant<>.

3. return by calling something like simple_var(value,y,dy,mu,dmu,sigma,dsigma);
[that function needs a better name, like partials_var].

What we're talking about now is mostly about how to
save ourselves a combinatorial explosion in writing
simple_var() and the underlying support class in the
face of variable/double/vector of variable/vector of double,
etc. arguments.

- Bob

[Bob Carpenter]

On 6/23/12 6:31 PM, Matt Hoffman wrote:
> We're just talking about chain(), right?

Both chain() and the constructor.

That's great the compiler does that much unfolding
for chains().

I'm writing out the first vector-based normal now.
I still don't see how to neatly combine the
vector and non-vector code without either lots of
if-else switches based on vector/non-vector or
with a simple expression template wrapper, which
defeats some of the optimizations.

I'd love to talk about some of this stuff if you'll
be around next week. (Otherwise, I'm bringing a lot
of napkins to Dinosaur BBQ.)

- Bob

Compilers are pretty hip to the
> for(size_t i = 0; i< constant; i++) {...}
> construction, and my understanding is that gcc unrolls pretty
> aggressively. Glancing at the object code in Instruments, it seems as
> though there aren't any jump instructions, at least for the version
> that I saw.
>
> Matt
>
> On Sat, Jun 23, 2012 at 6:26 PM, Bob Carpenter<ca...@alias-i.com> wrote:
>>
>>
>> On 6/23/12 5:56 PM, Matt Hoffman wrote:
>>>
>>> I'm not sure that there's an efficiency advantage to this

>>> approach�nvars_ is known at compile time, so shouldn't the compiler be

[Matt Hoffman]

I think they already have lots of napkins at DBBQ, although I think
it's for more practical reasons.

I've got what seems to be a working, fully vectorized, fully
double/var specializable normal_log() now. There are no messy
if/thens, but it's a little ugly, and the performance is a little
disappointing, but I think both of those things can be improved. I'll
be around next week and we can talk more about the details.

Matt

>>>> approach—nvars_ is known at compile time, so shouldn't the compiler be