Broken build: namespace issues.

[Daniel Lee]

One instance of our issue is in src/prob/distributions/multivariate/continuous/inv_wishart.hpp. There is a using statement for stan::math::determinant, defined in src/stan/math/matrix.hpp. The problem is that the stan::agrad::var version of determinant() exists in stan::agrad, defined in src/stan/agrad/matrix.hpp.

The simple solution is not what we want: adding a using statement for stan::agrad::determinant. That will create a dependency from our probability library to the stan::agrad::var objects.

Should we move the src/stan/agrad/matrix.hpp functions into namespace stan::math?



Daniel

PS. I am not sure how this worked before.

On Tue, Apr 17, 2012 at 11:40 AM, Matt Hoffman <mdho...@cs.princeton.edu> wrote:

I am not on top of it yet, so feel free to look into it.

On Tue, Apr 17, 2012 at 11:34 AM, Daniel Lee <bea...@alum.mit.edu> wrote:
> Thanks for modifying make clean.
>
> The build is still down, though. I'm taking a look at it now, unless you're
> already on top of it, Matt.
>
>
>
>
>
> On Mon, Apr 16, 2012 at 7:20 PM, Matt Hoffman <mdho...@cs.princeton.edu>
> wrote:
>>
>> -L specifies directories to check for libraries, -l specifies the library
>> names to include. (Analogous to -I and -i.)
>>
>> I can modify make clean, unless Daniel wants to do it.
>>
>> Matt
>>
>>
>> On Monday, April 16, 2012, Bob Carpenter wrote:
>>>
>>> Could someone (by that I mean Matt or Daniel)
>>> fix the build file so the clean task gets rid
>>> of these .o files?  Or is there a super clean target
>>> other than "make clean" that will do it?
>>>
>>> More comments inline.
>>>
>>> On 4/16/12 6:17 PM, Matt Hoffman wrote:
>>>>
>>>> You need to add the flags
>>>> -Lbin -lstan
>>>
>>>
>>> Thanks -- that's what I needed.  I can build the
>>> models with
>>>
>>>  > clang++ -O0 -I src -I lib -Lbin -lstan anon_model.cpp
>>>
>>> Is the difference between "-L" and "-l" the diff
>>> between archives (.a) and object files (.o)?  I
>>> didn't see either one with "clang++ --help" or "g++ --help".
>>>
>>>> to the linker call (assuming that libstan.a has been built by saying
>>>> "make bin/libstan.a"). This tells the linker to link against the stan
>>>> library. In the makefile, those flags are in the LDFLAGS variable.
>>>>
>>>> If you don't want to use a makefile at all, you can build libstan.a by
>>>> hand like so:
>>>> g++ -I src -I lib -O3 -g -c -o src/stan/agrad/agrad.o
>>>> src/stan/agrad/agrad.cpp
>>>> g++ -I src -I lib -O3 -g -c -o src/stan/math/matrix.o
>>>> src/stan/math/matrix.cpp
>>>> g++ -I src -I lib -O3 -g -c -o src/stan/agrad/matrix.o
>>>> src/stan/agrad/matrix.cpp
>>>> ar -rs bin/libstan.a src/stan/agrad/agrad.o src/stan/math/matrix.o
>>>> src/stan/agrad/matrix.o
>>>
>>>
>>> I just wanted to understand the steps without
>>> having to reverse engineer them from the makefiles,
>>> which I find very confusing.
>>>
>>> Our users are going to have to run make to get lib/stanc
>>> built anyway.  I'll put these installation instructions
>>> in an appendix of the manual.
>>>
>>> - Bob
>>>
>

[Bob Carpenter]

Argument-dependent name lookup:

http://en.wikipedia.org/wiki/Argument-dependent_name_lookup

The question is why won't this still work?

You still need to instantiate all the templates,
but you should be able to do one batch in one
.cpp and another batch in a different .cpp.

- Bob

On 4/17/12 12:47 PM, Daniel Lee wrote:
> One instance of our issue is in
> src/prob/distributions/multivariate/continuous/inv_wishart.hpp. There is
> a using statement for stan::math::determinant, defined in
> src/stan/math/matrix.hpp. The problem is that the stan::agrad::var
> version of determinant() exists in stan::agrad, defined in
> src/stan/agrad/matrix.hpp.
>
> The simple solution is not what we want: adding a using statement for
> stan::agrad::determinant. That will create a dependency from our
> probability library to the stan::agrad::var objects.
>
> Should we move the src/stan/agrad/matrix.hpp functions into namespace
> stan::math?
>
>
>
> Daniel
>
> PS. I am not sure how this worked before.
>
>
> On Tue, Apr 17, 2012 at 11:40 AM, Matt Hoffman
> <mdho...@cs.princeton.edu <mailto:mdho...@cs.princeton.edu>> wrote:
>
> I am not on top of it yet, so feel free to look into it.
>
> On Tue, Apr 17, 2012 at 11:34 AM, Daniel Lee <bea...@alum.mit.edu
> <mailto:bea...@alum.mit.edu>> wrote:
> > Thanks for modifying make clean.
> >
> > The build is still down, though. I'm taking a look at it now,
> unless you're
> > already on top of it, Matt.
> >
> >
> >
> >
> >
> > On Mon, Apr 16, 2012 at 7:20 PM, Matt Hoffman

> <mdho...@cs.princeton.edu <mailto:mdho...@cs.princeton.edu>>

[Ben Goodrich]

We shouldn't be calling determinant(S) or determinant(W) there anyway, as the underlying Eigen function does not take advantage of the presumed positive definiteness of these matrices. Better to call .llt() and get the determinant and inverse from L.

Ben

[Bob Carpenter]

On 4/17/12 6:37 PM, Ben Goodrich wrote:
> We shouldn't be calling determinant(S) or determinant(W) there anyway,
> as the underlying Eigen function does not take advantage of the presumed
> positive definiteness of these matrices. Better to call .llt() and get
> the determinant and inverse from L.

Sounds sensible.

Any volunteers to follow Ben's advice? I'd do it, but
my linear algebra's still too shaky to be trusted
in this setting (unless someone writes all the formulas
out more directly for me -- I'm happy to be the code monkey).

> On Tuesday, April 17, 2012 12:47:05 PM UTC-4, Daniel Lee wrote:
>
> One instance of our issue is in
> src/prob/distributions/multivariate/continuous/inv_wishart.hpp.
> There is a using statement for stan::math::determinant, defined in
> src/stan/math/matrix.hpp. The problem is that the stan::agrad::var
> version of determinant() exists in stan::agrad, defined in
> src/stan/agrad/matrix.hpp.
>
> The simple solution is not what we want: adding a using statement
> for stan::agrad::determinant. That will create a dependency from our
> probability library to the stan::agrad::var objects.

Right -- we don't want dependencies from stan::prob to stan::agrad.

The using statement could be added at a higher level if
we stuck to header-only usage, but we're not doing that, either.

So what we need is a template specialization declared in
the instantation file in src/stan/agrad/.

Then, if we import all the headers, argument-dependent lookup
should find the version from the right namespace.

> Should we move the src/stan/agrad/matrix.hpp functions into
> namespace stan::math?

That's also OK, but we can't in general have this because
of the std:: functions we implement auto-dif for.

> Daniel
>
> PS. I am not sure how this worked before.

All of the header definitions (not just declarations)
were in place for the single translation unit. At that
point, it's just argument dependent lookup again.

- Bob

[Marcus Brubaker]

For the record, I'm pretty certain this is a linker issue and not a namespace issue.  Specifically, it's not finding stan::math::determinant in libstan.a for some reason.  My library foo is weak, however the symbol list

[mbrubake@rincewind bin]$ nm -gC libstan.a | grep determinant
0000000000005590 T stan::math::determinant(Eigen::Matrix<double, -1, -1, 0, -1, -1>)
0000000000000000 W stan::agrad::var stan::agrad::determinant<stan::agrad::var>(Eigen::Matrix<stan::agrad::var, -1, -1, 0, -1, -1> const&)
0000000000000000 W stan::agrad::var stan::agrad::determinant<double>(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&)
0000000000000000 W Eigen::PartialPivLU<Eigen::Matrix<stan::agrad::var, -1, -1, 0, -1, -1> >::determinant() const

suggests that when the .o files are combined together, something different is happening with the stan::math version than the stan::agrad version.  (See the T vs W prefix.)  This is causing the linker to not find the math versions.  I had a theory that this had something to do with the fact that both object files had the same name but playing in that direction didn't yield any results.

About using LLT here, Bob check out multi_normal_log_chol for an idea.  The determinant of a matrix A = L*L^T is the product of the square of the diagonals.  I.E., det(A) = \prod_i L(i,i)^2 or, more stably, log(det(A)) = 2 * \sum_i log(L(i,i)).

Cheers,
Marcus

[Ben Goodrich]

On Tuesday, April 17, 2012 7:22:21 PM UTC-4, Bob Carpenter wrote:

On 4/17/12 6:37 PM, Ben Goodrich wrote:
> We shouldn't be calling determinant(S) or determinant(W) there anyway,
> as the underlying Eigen function does not take advantage of the presumed
> positive definiteness of these matrices. Better to call .llt() and get
> the determinant and inverse from L.

Sounds sensible.

Any volunteers to follow Ben's advice?  

I'm on it. Wishart needs to get fixed too.

Ben

[Marcus Brubaker]

Ok, turns out this was a much simpler bug than a linker issue...there was a & dropped when stan/math/matrix.cpp was created so that there was no longer a definition of stan::math::determinant.  I've now learned far more about linking than I ever wanted to know... :)

Anyway, I've fixed and pushed this and test-unit now builds and passes for me.

Cheers,
Marcus

[Daniel Lee]

Thanks!

[Ben Goodrich]

What do I need to do to get multiply() to work in this chunk?

 if (include_summand<propto,T_y,T_scale>::value) {
    Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic> I(k,k);
    I.setIdentity();
    L = L.template triangularView<Eigen::Lower>().solve(I).array().square().matrix();
    Eigen::Matrix<T_scale,Eigen::Dynamic,Eigen::Dynamic> d = S.diagonal();
    lp -= 0.5 * stan::math::multiply(L, d).array().sum(); // tr(SW^-1)
}

I get all these errors like

In file included from src/test/agrad/distributions/multivariate/continuous/inv_wishart_test.cpp:3:                                                                                            
src/stan/prob/distributions/multivariate/continuous/inv_wishart.hpp:90:21: error: no matching function for call to 'multiply'

Ben

[Marcus Brubaker]

On Tue, Apr 17, 2012 at 9:53 PM, Ben Goodrich <goodri...@gmail.com> wrote:

On Tuesday, April 17, 2012 8:13:34 PM UTC-4, Ben Goodrich wrote:

On Tuesday, April 17, 2012 7:22:21 PM UTC-4, Bob Carpenter wrote:

On 4/17/12 6:37 PM, Ben Goodrich wrote:
> We shouldn't be calling determinant(S) or determinant(W) there anyway,
> as the underlying Eigen function does not take advantage of the presumed
> positive definiteness of these matrices. Better to call .llt() and get
> the determinant and inverse from L.

Sounds sensible.

Any volunteers to follow Ben's advice?  

I'm on it. Wishart needs to get fixed too.

What do I need to do to get multiply() to work in this chunk?

 if (include_summand<propto,T_y,T_scale>::value) {
    Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic> I(k,k);
    I.setIdentity();
    L = L.template triangularView<Eigen::Lower>().solve(I).array().square().matrix();

    Eigen::Matrix<T_scale,Eigen::Dynamic,Eigen::Dynamic> d = S.diagonal();
    lp -= 0.5 * stan::math::multiply(L, d).array().sum(); // tr(SW^-1)

Just call multiply(L,d) not stan::math::multiply(L,d).  You'll probably need a "using stan::math::multiply;" somewhere too.

Cheers,
Marcus

[Bob Carpenter]

On 4/17/12 10:38 PM, Marcus Brubaker wrote:

Right -- this is related to Daniel's earlier
question. Then the test will need to include
"stan/agrad/matrix.hpp" and everything should be
OK after that.

- Bob

[Ben Goodrich]

OK, almost working now. Is there a customized function to multiply L by d when L is lower-triangular and d is a vector?

Ben
 

[Bob Carpenter]

On 4/17/12 9:03 PM, Marcus Brubaker wrote:
> Ok, turns out this was a much simpler bug than a linker issue...there was a & dropped when stan/math/matrix.cpp was
> created so that there was no longer a definition of stan::math::determinant. I've now learned far more about linking
> than I ever wanted to know... :)

I feel the same way, except that I wanted to know
how it worked and never did until I refactored
the parser into five pieces.

The Vandevoorde and Josuttis book is pretty good on
this front, but I really had to get my hands dirty
before I really understood it.

> Anyway, I've fixed and pushed this and test-unit now builds and passes for me.

Great.

- Bob

[Bob Carpenter]

On 4/17/12 11:07 PM, Ben Goodrich wrote:

> OK, almost working now. Is there a customized function to multiply L by d when L is lower-triangular and d is a vector?

Not in Stan. You could always add one.
It would get a different sized dot product
in the value of the result vector.

I don't know how Eigen deals with
(triangular-matrix * vector) multiplication.
To state the obvious, I'd look for operator*
with a triangular matrix argument.

- Bob

[Marcus Brubaker]

It should work (in theory) much like the triangular solvers work (see mdivide_left_tri).  You take a triangular view of the matrix and then perform your operation based on that.  Eigen has specializations for triangular matrix products so it should call those when needed.  However, that's unlikely to work properly right now if you have mixed types.  For the time being you may be stuck with evaluating out the full thing untoi lwe have better handling of mixed types and are more tightly integrated with Eigen's operations.

Cheers,
Marcus

[Bob Carpenter]

The easiest thing to do is call Stan's multiply().
It will probably be faster than promoting double to var
and calling Eigen's specialized triangular multiply
because of the auto-dif savings.

Better yet would be to write a version of multiply()
for triangular matrices into in agrad/matrix.hpp.
It would have tighter loop bounds and use different
dot_product lengths based on which row you're using
from the matrix. To avoid promotion, this requires
four implementations ({double,var}*{double,var}).

Longer term, of course, we want to have our
expression template cake and eat it, too.

- Bob

On 4/18/12 12:48 AM, Marcus Brubaker wrote:
> It should work (in theory) much like the triangular solvers work (see mdivide_left_tri). You take a triangular view of
> the matrix and then perform your operation based on that. Eigen has specializations for triangular matrix products so
> it should call those when needed. However, that's unlikely to work properly right now if you have mixed types. For the
> time being you may be stuck with evaluating out the full thing untoi lwe have better handling of mixed types and are
> more tightly integrated with Eigen's operations.
>
> Cheers,
> Marcus
>

[Ben Goodrich]

On Tuesday, April 17, 2012 11:07:08 PM UTC-4, Ben Goodrich wrote:

OK, almost working now. Is there a customized function to multiply L by d when L is lower-triangular and d is a vector?

Pushed. Although multiplication involving a triangular matrix would be nice to have, it turns out that I was doing it wrong and didn't need it right now. At the moment, there is an unrelated (I think) build error for src/models/bugs_examples/vol1/dogs/ complaining about malloc and whatnot.

[Bob Carpenter]

On 4/18/12 3:04 PM, Ben Goodrich wrote:
> On Tuesday, April 17, 2012 11:07:08 PM UTC-4, Ben Goodrich wrote:
>
> OK, almost working now. Is there a customized function to multiply L
> by d when L is lower-triangular and d is a vector?
>
>

> Pushed.[Cholesky-factored distributions]

That's great. Now I can look at the code to figure
out how you did it!

> Although multiplication involving a triangular matrix would be
> nice to have, it turns out that I was doing it wrong and didn't need it
> right now. At the moment, there is an unrelated (I think) build error
> for src/models/bugs_examples/vol1/dogs/ complaining about malloc and
> whatnot.

I'll look into that later if I get some time.

There's also a problem in the unit tests, which I'm
about to fix (appended below).

- Bob

clang++ lib/gtest/src/gtest_main.cc
test/prob/distributions/multivariate/continuous/inv_wishart.o -I
lib/gtest/include -I lib/gtest -Lbin -lstan test/libgtest.a -o
test/prob/distributions/multivariate/continuous/inv_wishart
test/prob/distributions/multivariate/continuous/inv_wishart
--gtest_output="xml:test/prob/distributions/multivariate/continuous/inv_wishart.xml"
Running main() from gtest_main.cc
[==========] Running 4 tests from 1 test case.
[----------] Global test environment set-up.
[----------] 4 tests from ProbDistributionsInvWishart
[ RUN ] ProbDistributionsInvWishart.InvWishart
[ OK ] ProbDistributionsInvWishart.InvWishart (0 ms)
[ RUN ] ProbDistributionsInvWishart.Propto
[ OK ] ProbDistributionsInvWishart.Propto (0 ms)
[ RUN ] ProbDistributionsInvWishart.DefaultPolicy
src/test/prob/distributions/multivariate/continuous/inv_wishart_test.cpp:65:
Failure
Expected: inv_wishart_log(Y, nu, Sigma) doesn't throw an exception.
Actual: it throws.
[ FAILED ] ProbDistributionsInvWishart.DefaultPolicy (0 ms)
[ RUN ] ProbDistributionsInvWishart.ErrnoPolicy
[ OK ] ProbDistributionsInvWishart.ErrnoPolicy (0 ms)
[----------] 4 tests from ProbDistributionsInvWishart (1 ms total)

[----------] Global test environment tear-down
[==========] 4 tests from 1 test case ran. (1 ms total)
[ PASSED ] 3 tests.
[ FAILED ] 1 test, listed below:
[ FAILED ] ProbDistributionsInvWishart.DefaultPolicy

1 FAILED TEST

[Ben Goodrich]

On Wednesday, April 18, 2012 3:12:37 PM UTC-4, Bob Carpenter wrote:

There's also a problem in the unit tests, which I'm
about to fix (appended below).

I fixed a couple of those in 8ff07a88ff275c6db1833d752c5dec27d148838d . Do you have that commit? They arise because the test was expecting a no throw due to the degrees of freedom but passing in an unallocated scale matrix. So, now that *wishart_log checks the scale matrix for positive definiteness, it was throwing an error. The fix is to change the test to pass the identity matrix or whatever.

Ben

 

[Bob Carpenter]

Ben -- are you not on the build list to get the
build errors e-mailed from the Jenkins?

The problem is here (with an identical problem for Wishart):

TEST(ProbDistributionsInvWishart,DefaultPolicy) {
Matrix<double,Dynamic,Dynamic> Sigma;
Matrix<double,Dynamic,Dynamic> Y;
double nu;

Sigma.resize(1,1);
Y.resize(1,1);
nu = 1;
EXPECT_NO_THROW(inv_wishart_log(Y, nu, Sigma));

The new implementation is balking on 1x1 covariance matrices.
(Hat-tip to Daniel for writing such extensive error tests
that they caught this boundary condition.)

For now, I added a test to throw a domain error (or
whatever the policy is) if either Wishart or inverse Wishart
get a 1x1 covariance matrix and changed the unit tests
to match and everything's passing agin.

But I see no reason why it can't work for 1 x 1 matrices,
and it's not being checked, so it must be one of the
linear algebra ops throwing an exception.

Ben -- could you fix it so it works for 1 x 1 matrices?
(There's really not a huge rush on this -- it's just nice
to have everything work in full generality -- it's easier
to doc, for one thing.)

If the (inv_)Wishart is really a generalization of
(inv_)chi_square, we could just call those.

Otherwise, we can add a K=1 case to both distributions' code.

if (k == 1) {
// handle 1 x 1 case
}
... rest of code

If the Wishart and its inverse twin are changed
to handle 1x1 matrices we should also have a test
for the values produced on 1x1 matrices if they
are indeed implemented separately in the 1 x 1 case.

- Bob

On 4/18/12 3:04 PM, Ben Goodrich wrote:

[Bob Carpenter]

Luckily, I just saw Ben's message before I
committed an egegrious error (not including that
pun).

I'll wait for Ben's fix.

- Bob

[Ben Goodrich]

On Wednesday, April 18, 2012 3:38:26 PM UTC-4, Bob Carpenter wrote:

Ben -- are you not on the build list to get the
build errors e-mailed from the Jenkins?

I guess not but I can read them from Google Groups.
 

The problem is here (with an identical problem for Wishart):

   TEST(ProbDistributionsInvWishart,DefaultPolicy) {
     Matrix<double,Dynamic,Dynamic> Sigma;
     Matrix<double,Dynamic,Dynamic> Y;
     double nu;

     Sigma.resize(1,1);
     Y.resize(1,1);
     nu = 1;
     EXPECT_NO_THROW(inv_wishart_log(Y, nu, Sigma));

The new implementation is balking on 1x1 covariance matrices.

I think the code should work for 1x1. I think it is balking on the fact that Sigma and Y are not allocated or maybe are allocated with zeros (or maybe are allocated with ones on Linux because it worked for me). So, the fix that I just pushed (which I did in the other NOTHROW case) was to

Sigma.setIdentity()
Y.setIdentity()

which got Jenkins back to normal.

For now, I added a test to throw a domain error (or
whatever the policy is) if either Wishart or inverse Wishart
get a 1x1 covariance matrix and changed the unit tests
to match and everything's passing agin.

But I see no reason why it can't work for 1 x 1 matrices,
and it's not being checked, so it must be one of the
linear algebra ops throwing an exception.

We can certainly fix it more if need be, but I am inclined to say that stan should throw a domain error in the 1x1 case. I could see why the underlying Eigen functions should support 1x1, but I would guess that the majority of the time when this occurs in stan, it is a mistake on the user's part.

Ben