fixed matrix elt-divide + unit tests

[Bob Carpenter]

I fixed the element-wise division so it actually worked
(the matrix case was broken), renamed the misnamed file
Marcus pointed out, then added unit tests. I also switched
the order of all of our loops over rows/cols to be memory
local (Eigen default's to column major, which is what we're
using).

I like git, but I'm still getting the hang of it and
tend to do "commit -a" when I just mean "commit". I read
the chapter of the "Git Pro" book on merges, and it looks
awesome. Now if I can just remember not to do "commit -a"
reflexively, the commits should start to be more coherent.

- Bob

[Andrew Gelman]

Just curious: stan can handle multi-way arrays, I assume. Can it do vector and matrix operations using R-like conventions such as a[i,,] or a[,j,] or a[,j,k]?
A

[Bob Carpenter]

No,

On 3/7/12 7:35 PM, Andrew Gelman wrote:
> Just curious: stan can handle multi-way arrays, I assume.

Yes.

> Can it do vector and matrix operations using R-like conventions such as a[i,,] or a[,j,] or a[,j,k]?

No.

It's on the nice-to-have to-do list.

One thing that makes it less necessary is that loops and
assignments are much faster in Stan than in R.

- Bob

[Andrew Gelman]

Thx, My general opinion on such cool features is that we can wait until we're working on a model where we have a great desire for such a feature, then that will motivate us to implement it.
A

[Marcus Brubaker]

If those things aren't high priority for you guys, I might take a crack at implementing them in the not-too-distant future.  Particularly for the larger models which I'm playing with (on the order of 1e5 parameters) things like that will make a big difference in performance and readability of the model specifications.

[Bob Carpenter]

The nice thing about open source is that if you do the
work, you can implement the cool features you want.

Extending the modeling language with better matrix
notation is great. Especially given that adding
new operator types is pretty painless given the
parsing and type-checking framework in place.

New functions like id_matrix(K) or constant_matrix(M,N,x)
are also pretty easy to add.

- Bob

[Ben Goodrich]

 On Wednesday, March 7, 2012 5:55:49 PM UTC-5, Bob Carpenter wrote:

I fixed the element-wise division so it actually worked
(the matrix case was broken), renamed the misnamed file
Marcus pointed out, then added unit tests.  I also switched
the order of all of our loops over rows/cols to be memory
local (Eigen default's to column major, which is what we're
using).

I'm just now looking over the new elt_* functions in src/stan/agrad/matrix.hpp . These are not checking whether m1 and m2 are conformable for elementwise operations. Is this checked somewhere else?

Also, is there no type promotion hackery that would allow us to take advantage of Eigen's templated stuff like

( m1.array() OP m2.array() ).matrix();

I am guessing not (yet) but life will be so much easier when Eigen implements mixed type ops.

Ben

[Marcus Brubaker]

That should work for the elt_* functions.  I'm using a similar pattern for element-wise exp/log functions.  I'm just finishing adding unit tests and will commit that change.

Cheers,
Marcus

[Ben Goodrich]

On Saturday, March 10, 2012 7:17:32 PM UTC-5, Marcus Brubaker wrote:

( m1.array() OP m2.array() ).matrix();

I am guessing not (yet) but life will be so much easier when Eigen implements mixed type ops.

That should work for the elt_* functions.  I'm using a similar pattern for element-wise exp/log functions.  I'm just finishing adding unit tests and will commit that change.

That is idiomatic Eigen syntax, but Stan typically has compilation problems when m1 consists of parameters and m2 consists of doubles. In fact, I just tried this for elt_multiply() and the tests failed. It might be possible to promote m2 to whatever type m1 is and then do

( m1.array() OP m2.array() ).matrix();

which is basically what stan does with genuine matrix multiplication. All of this type promoting we would rather not have, but Eigen has not yet implemented Bob's feature request to allow binary operations among mixed types.

It has been a couple months since I last tried to trick Eigen into accepting things like the above. Now that we have git (branches), it is easier to take a stab at it.

Ben

[Bob Carpenter]

On 3/10/12 7:17 PM, Marcus Brubaker wrote:
> On Sat, Mar 10, 2012 at 7:00 PM, Ben Goodrich <goodri...@gmail.com
> <mailto:goodri...@gmail.com>> wrote:
>
> On Wednesday, March 7, 2012 5:55:49 PM UTC-5, Bob Carpenter wrote:
>
> I fixed the element-wise division so it actually worked
> (the matrix case was broken), renamed the misnamed file
> Marcus pointed out, then added unit tests. I also switched
> the order of all of our loops over rows/cols to be memory
> local (Eigen default's to column major, which is what we're
> using).
>
>
> I'm just now looking over the new elt_* functions in
> src/stan/agrad/matrix.hpp . These are not checking whether m1 and m2
> are conformable for elementwise operations. Is this checked
> somewhere else?

It's certainly not being checked elsewhere.
We want to catch these indexing errors and throw
exceptions.

> Also, is there no type promotion hackery that would allow us to take
> advantage of Eigen's templated stuff like
>
> ( m1.array() OP m2.array() ).matrix();
>
> I am guessing not (yet) but life will be so much easier when Eigen
> implements mixed type ops.
>
>
> That should work for the elt_* functions. I'm using a similar pattern
> for element-wise exp/log functions. I'm just finishing adding unit
> tests and will commit that change.

Right -- this stuff's no problem as long as you can
make the types consistent. But if m1 is double and
m2 is var, it probably won't work.

- Bob