Speed Results

[Ben Goodrich]

I created a wiki page to compare speed the speed models across configurations:

https://code.google.com/p/stan/wiki/SpeedResults

Unlike the test-unit results where we are fishing for bugs, here I think we only need to compare the best versions of each compiler for whatever the optimal level of optimization is (I've done O3 for now, but for me O2 is about the same for execution and somewhat faster for compilation). It would be good if someone could put some Mac results up there. Windows, too, I guess although any compiler other than the one included in Rtools is probably a headache to get working.

Note that the wiki is a separate git repo

https://code.google.com/p/stan/source/checkout?repo=wiki

And rather than using the web interface, it is probably easier to clone, edit that page with a text editor, add, commit, and push.

Ben

P.S. We have got to do something about the schools model in vol2.

[Daniel Lee]

Whoa. That's a lot of work. 

Just wondering, were the seeds fixed?

--
 
 

[Ben Goodrich]

A lot of work for the CPU, not for the human. The file sourced fixes the seed but otherwise pretty much uses the defaults.

Ben

[Jiqiang Guo]

I would try to get some results from Windows.  Actually, I can get results for both Windows and Linux from the same computer. 

For comparing O3 and O2, is package RStan also installed differently? 

--

Jiqiang 

On Wed, Dec 19, 2012 at 5:00 PM, Ben Goodrich <goodri...@gmail.com> wrote:

A lot of work for the CPU, not for the human. The file sourced fixes the seed but otherwise pretty much uses the defaults.

Ben

--

[Ben Goodrich]

On Wednesday, December 19, 2012 8:40:32 PM UTC-5, Jiqiang Guo wrote:

I would try to get some results from Windows.  Actually, I can get results for both Windows and Linux from the same computer. 

For comparing O3 and O2, is package RStan also installed differently? 

I've done O3 vs. O2 on my laptop just through .R/Makevars and not found much difference in execution time with previous versions of RStan. I suppose to be really rigorous about it, I should change the O in rstan/src/Makefile locally and reinstall.

The thing about Windows is that we don't have any concrete experience with any compiler except the one in Rtools, so there is not much to say except for the difference between O2 and O3. If you do manage to get clang, Intel's, MSVC, etc. whatever working with RStan on Windows be sure to write down your steps.

Ben

[Andrew Gelman]

Cool!

> --
>
>

[Bob Carpenter]

> On Dec 19, 2012, at 4:33 PM, Ben Goodrich wrote:

...

>> P.S. We have got to do something about the schools model in vol2.

I'll second that. Very cool.

I'm sort of amazed at the variation in compilers
across models. But then maybe that's random init
differences? Are they all running with the same seed?

Schools is a pretty nasty model! But there's definitely
room for both vectorization in the big data loop at the
top of the model.

And it's time to generalize the Matt trick to multivariate
models for alpha in vol2/schools/schools.stan

transformed data {
vector[3] gamma_mu;
cov_matrix[3] gamma_Sigma;
cov_matrix[3] invR;
invR <- inverse(R);
...
parameters {
...
vector[3] alpha[M];
...
model {

# Priors for random coefficients:
for (m in 1:M) alpha[m] ~ multi_normal(gamma, Sigma);

# Hyper-priors:
gamma ~ multi_normal(gamma_mu, gamma_Sigma); //
Sigma ~ inv_wishart(3, invR);

- Bob

[Ben Goodrich]

On Thursday, December 20, 2012 2:18:35 AM UTC-5, Bob Carpenter wrote:

> On Dec 19, 2012, at 4:33 PM, Ben Goodrich wrote:
I'm sort of amazed at the variation in compilers
across models.  But then maybe that's random init
differences?  Are they all running with the same seed?

Same seed. TFC includes scarcely more than this R function

bench <- function(i) {
  mark <- "elapsed"
  compile <- system.time(fit <- stan_demo(i, chains = 0, iter = 0))[mark]
  execute <- system.time(fit <- stan_demo(i, fit = fit, seed = 12345, refresh = -1))[mark]
  if (length(dim(fit)) == 0) execute <- NA_real_ # failure
  return(c(compile, execute))
}

which is executed on each .stan file. Something could be said for increasing the number of observations so that we have fewer examples that take a tenth of a second to execute, but that would be killer for the examples that take minutes to execute.

And it's time to generalize the Matt trick to multivariate
models for alpha in vol2/schools/schools.stan

We do have stuff in the manual about that. More generally, are the bugs_examples intended to be literal translations of BUGS examples or should we be tricking them out? In a couple of cases, we have both an unvectorized and a vectorized version of the model, but that is not so helpful when the unvectorized version is unreasonably slow.

Ben

[Bob Carpenter]

On 12/20/12 2:59 AM, Ben Goodrich wrote:
> On Thursday, December 20, 2012 2:18:35 AM UTC-5, Bob Carpenter wrote:

> And it's time to generalize the Matt trick to multivariate
> models for alpha in vol2/schools/schools.stan
>
>
> We do have stuff in the manual about that. More generally, are the bugs_examples intended to be literal translations of
> BUGS examples or should we be tricking them out?

The BUGS examples are being used for both evaluation
and pedagogy.

For both purposes, we want to leave the models themselves
unchanged (I'd be OK with reparameterizations where convenient)

The JAGS versions modify the BUGS versions in
some cases. We decided ages ago to take the JAGS versions
as our basis.

For evaluation, it's clear we should do efficient
implementations.

On the pedagogy side, we have the conflicting goals
of clarity and speed.

> In a couple of cases, we have both an unvectorized and a vectorized
> version of the model, but that is not so helpful when the unvectorized version is unreasonably slow.

I see room for two versions -- one clear and one fast.

Pedagogically, that would be the best of both worlds.

- Bob

[Ben Goodrich]

On Friday, December 21, 2012 2:45:03 PM UTC-5, Bob Carpenter wrote:

For evaluation, it's clear we should do efficient
implementations.

On the pedagogy side,  we have the conflicting goals
of clarity and speed.

> In a couple of cases, we have both an unvectorized and a vectorized
> version of the model, but that is not so helpful when the unvectorized version is unreasonably slow.

I see room for two versions -- one clear and one fast.

Pedagogically, that would be the best of both worlds.

It's not clear what clear entails. Or to whom, the average BUGS user or the average R user or ...? To take that model as an example

model {
  for(p in 1:N) {
     Y[p] ~ normal(alpha[school[p], 1] + alpha[school[p], 2] * LRT[p] +
                   alpha[school[p], 3] * VR[p, 1] + beta[1] * LRT[p] * LRT[p] +
                   beta[2] * VR[p, 2] + beta[3] * Gender[p] +
                   beta[4] * School_gender[p, 1] + beta[5] * School_gender[p, 2] +
                   beta[6] * School_denom[p, 1] + beta[7] * School_denom[p, 2] +
                   beta[8] * School_denom[p, 3],  exp(-.5 * (theta + phi * LRT[p])));
  }
  // min.var <- exp(-(theta + phi * (-34.6193))) # lowest LRT score = -34.6193
  // max.var <- exp(-(theta + phi * (37.3807)))  # highest LRT score = 37.3807

  # Priors for fixed effects:
  beta ~ normal(0, 100);
  // for (k in 1:8)  beta[k] ~ normal(0.0, 100);
  theta ~ normal(0.0, 100);
  phi ~ normal(0.0, 100);

  # Priors for random coefficients:

  for (m in 1:M) alpha[m] ~ multi_normal(gamma, Sigma);
  # Hyper-priors:
  gamma ~ multi_normal(gamma_mu, gamma_Sigma);

  Sigma ~ inv_wishart(3, invR);
}

Would it be less clear to instead do the faster

model {
  vector[N] mu;
  vector[N] sigma;
  for(p in 1:N) {
     mu[p] <- alpha[school[p], 1] + alpha[school[p], 2] * LRT[p] +
              alpha[school[p], 3] * VR[p, 1] + beta[1] * LRT[p] * LRT[p] +
              beta[2] * VR[p, 2] + beta[3] * Gender[p] +
              beta[4] * School_gender[p, 1] + beta[5] * School_gender[p, 2] +
              beta[6] * School_denom[p, 1]  + beta[7] * School_denom[p, 2] +
              beta[8] * School_denom[p, 3]; 
     sigma[p] <- exp(-.5 * (theta + phi * LRT[p]));
  }
  Y ~ normal(mu, sigma);


? To me that is no less clear, but it is certainly less BUGS-ish and more R-ish, and similarly for the vectorize vs. for loop choice

  # Priors for fixed effects:
  beta ~ normal(0, 100);
  // for (k in 1:8)  beta[k] ~ normal(0.0, 100);


I guess Matt tricking around a multi_normal or inv_wishart is both less clear, less BUGS-ish, and less R-ish.

Ben

[Andrew Gelman]

In general I think we should translate the bugs example so as to run as fast as we can in Stan. I think it's very appropriate for us to use the Matt trick wherever it makes sense. People will use our examples as templates (as they do with the bugs examples).
A

P.S. Also, I think we should consider putting all our examples in a big pdf file. Yes, I know they are all there on the webpage, but I think there's something about a recipe book which is helpful to people.

> --
>

[Bob Carpenter]

On 12/21/12 4:10 PM, Andrew Gelman wrote:
> In general I think we should translate the bugs example so as to run as fast as we can in Stan. I think it's very appropriate for us to use the Matt trick wherever it makes sense. People will use our examples as templates (as they do with the bugs examples).

Exactly.

That's why we kept bringing up the issue of whether
we should change all the inv_gamma(0.001,0.001) priors :-)

I think it's helpful for people to see the simpler
models (without vectorization) AND the vectorized
models, so they can see how to get from one to the other.

> P.S. Also, I think we should consider putting all our examples in a big pdf file. Yes, I know they are all there on the webpage, but I think there's something about a recipe book which is helpful to people.

Formatting is problematic for this kind of thing.
The attached uses 11pt font on regular notebook
paper (8.5" x 11"), with text width 6.5" and height
9.5". That'll make it challenging to read on an iPad,
but it should print OK.

I also include a small sized model (room for 88 chars)
and a footnote sized model (98 chars, but obviously
too small).

Any suggestions on formatting?

I can just write a script to include all of our models,
each in its own chapter.

Of course, it'd be nice to go through and add references
and maybe some discussions of the models and implementation
choices. (And add all those warnings about copying the
priors!)


- Bob

[Bob Carpenter]

I share your notion of clarity.

The Matt trick's clarity is probably an issue with
habit. In some ways, it's clearer to think of
a bunch of normal variates getting rescaled than to
think about drawing from all sorts of differently scaled
and located distributions.

I think having to create a vector and write into it
for vectorization is less clear for everyone.

The main reason plain old vectorization may be less
clear is that it's not immediately evident to either
the BUGS or R user how we've vectorized. That is,
we allow any of the arguments to be a scalar that's
just repeated. Most R functions won't let you do that,
will they?

- Bob

> --
>
>

[Andrew Gelman]

Yes, that's what I was thinking too. Programming using the Matt trick is as easy as any other version. You just get used to having a lot of N(0,1)'s around.

[Andrew Gelman]

I don't think we should redo all the Bugs models, but gamma(,001,.001) just jumps out at me. I really feel uncomfortable having it in our templates.

[Ben Goodrich]

On Friday, December 21, 2012 6:09:26 PM UTC-5, Bob Carpenter wrote:

I share your notion of clarity.

The Matt trick's clarity is probably an issue with
habit.  In some ways, it's clearer to think of
a bunch of normal variates getting rescaled than to
think about drawing from all sorts of differently scaled
and located distributions.

I'd say the univariate Matt trick is only slightly less clear than just having normal(mu[i], sigma). But the multivariate Matt tricks are not that clear to me, even after writing several of them into the manual.

I think having to create a vector and write into it
for vectorization is less clear for everyone.

The main reason plain old vectorization may be less
clear is that it's not immediately evident to either
the BUGS or R user how we've vectorized.  That is,
we allow any of the arguments to be a scalar that's
just repeated.  Most R functions won't let you do that,
will they?

Most R functions *do* let you do that, e.g.

x <- rnorm(1000)
sum(dnorm(x, mean = rep(0,1000), sd = 1, log = TRUE))

is valid. That is why I was saying that being less BUGS-ish is less clear to someone who only knows WinBUGS and Stata / SPSS / etc. but more clear to someone who only knows R.

Ben

[Ben Goodrich]

On Friday, December 21, 2012 5:55:17 PM UTC-5, Bob Carpenter wrote:

> P.S.  Also, I think we should consider putting all our examples in a big pdf file.  Yes, I know they are all there on the webpage, but I think there's something about a recipe book which is helpful to people.

Formatting is problematic for this kind of thing.
The attached uses 11pt font on regular notebook
paper (8.5" x 11"), with text width 6.5" and height
9.5".  That'll make it challenging to read on an iPad,
but it should print OK.

I also include a small sized model (room for 88 chars)
and a footnote sized model (98 chars, but obviously
too small).

Any suggestions on formatting?

Why PDF as opposed to just text in HTML? The web browser can help you out on sizing and is easier to copy / paste with.

Ben

[Andrew Gelman]

Maybe html is better. I think the key is that it's in one place. The Bugs examples are in html or something similar, I think.

On Dec 21, 2012, at 6:43 PM, Ben Goodrich wrote:

> Why PDF as opposed to just text in HTML? The web browser can help you out on sizing and is easier to copy / paste with.
>
> Ben
>
>

> --
>
>

[Ben Goodrich]

I think HTML could be better, but if we want PDF the attached script (to be put in src/docs/stan-reference/) gets us most of the way there if pdflatex the resulting examples.tex. The problem is that many of the .stan files have line lengths that are too long.

Ben

[Andrew Gelman]

I defer to the rest of you.  I just think that a single document of some sort will be helpful.  The current system is well organized but I think is a bit difficult to browse, and browsing is one way for people to find a template.

--
 
 
<do_examples.sh>

[Bob Carpenter]

On 12/21/12 8:58 PM, Andrew Gelman wrote:
> I defer to the rest of you. I just think that a single document of some sort will be helpful. The current system is
> well organized but I think is a bit difficult to browse, and browsing is one way for people to find a template.
>
> On Dec 21, 2012, at 8:55 PM, Ben Goodrich wrote:
>
>> I think HTML could be better,

We already have basic HTML through Google Code's
source browser, the following link to which is
on the Stan home page as per Andrew's original request
for a pointer to browsable models:

http://code.google.com/p/stan/source/browse/#git%2Fsrc%2Fmodels%2Fbugs_examples

>> but if we want PDF the attached script (to be put in src/docs/stan-reference/) gets us
>> most of the way there if pdflatex the resulting examples.tex. The problem is that many of the .stan files have line
>> lengths that are too long.

:-) Yeah, that's what I meant.

I think the solution here is to modify the .stan files
so they don't have lines longer than 80 chars if at
all possible.

On 12/21/12 6:50 PM, Andrew Gelman wrote:
> Maybe html is better. I think the key is that it's in one place.

You mean not having to click through directory structures?

> The Bugs examples are in html or something

What's funny is that they have PDFs that look
like HTML :-) They're organized by volume, as in:

http://www.mrc-bsu.cam.ac.uk/bugs/winbugs/Vol1.pdf

- Bob

[Andrew Gelman]

I don't know. It might be something as simple as that, with Bugs, the examples are in one place that's clear bounded, while with Stan, it's part of a big ugly directory including all sorts of other things. The googlecode interface is much much better than the interface in win bugs, but it looks a bit more intimidating, it's got buttons all over the place, it looks like a developer tool. Maybe people can handle it, I don't know. There's something about it that looks scary.

It's funny, here I'm taking the opposite position I take in the infovis-vs-stat-graphics debate. I'm saying that this is a superior interface for people who know what they're doing, but it's scarier for outsiders.

Given that we have the html already, I think it would be fine if it were supplemented with a pdf document that just had the R and stan code (but not the data). That way, newcomers could start out by downloading Stan, the Stan manual, and the Stan examples.

A

> --
>

[Bob Carpenter]

I like the pdf -- I find it easier to browse
than having to navigate among a bunch of web
pages.

I'd think the place to go from here would be to
have a format for dumping LaTeX into the header
of the models and then putting that in the doc.
I'd like to keep the model doc as close as possible
to the models.

Somehow my bash on the Mac requires more escapes in
the shell script, so I started trying to debug Ben's
script, then decided it'd be easier for me to just write
a Java program. (One of these days I'll learn shell
scripting, really!)

It's a lot of pages, but not many bytes, so I
attached the output, which is still pretty crude.

I'm just going through the directories recursively
looking for .stan files, the same way Ben's script did.

- Bob

[Ben Goodrich]

On Saturday, December 22, 2012 2:35:24 AM UTC-5, Bob Carpenter wrote:

I like the pdf -- I find it easier to browse
than having to navigate among a bunch of web
pages.

OK, but at least make the table of contents clickable. Also, it would be cool if we could append a function index to it.

Ben

[Bob Carpenter]

Clickable contents obviously and maybe a title.

The function index is a great idea and should be pretty
easy because the model syntax is so easy.

- Bob

> --
>
>

[Bob Carpenter]

Here's a verison with function index and clickable
contents. And a title.

There are false positives in the index because I
just used a regex to pull out function names matching
our variable syntax followed by an open paren. I could
extract a list of all of our functions and match it,
but I think we probably just want to do something
different with comments.

The next thing I'll do is read any initial comments in
the model doc straight into LaTeX.

Then we just need to go through and clean the models up
in terms of line overruns, variable naming, using more
vectorization, using some better parameterizations, etc.

I'll take a first pass at that. I've been meaning to do it
for ages anyway.

I also attached Java code I used to do the generation.
Don't laugh -- being Java, it's pretty clunky compared to
terse shell scripts. On the other hand, it's very portable.

- Bob

On 12/22/12 3:13 AM, Ben Goodrich wrote:

> --
>
>

[Michael Betancourt]

Any tex magic to avoid page breaks within braces? I'm not sure how to use the usual orphan/widow with multiline blocks.

> --
>
>
> <MakeModelTex.java><models-out.pdf>

[Bob Carpenter]

On 12/22/12 6:48 PM, Michael Betancourt wrote:
> Any tex magic to avoid page breaks within braces? I'm not sure how to use the usual orphan/widow with multiline blocks.

I don't think so -- the models are all in verbatim environs.
For one, the basic paragraph marker (a blank line) isn't there.

We cold replace the verbatims with line-based things and try
to do it ourselves, but that's a bit more fine detail than I'd
like to go into. I did that for our company's manual, but it was
a pain to do all the escaping. The upside is that it let me replace
vertical lines in code with half lines and play with fonts
and things inside the code to bold pieces, etc. And also chop
down to rectangles inside of code based on comment breadcrumbs.

- Bob

[Jiqiang Guo]

I hate to say I am giving up on this, trying to get a multiple compiler on windows that works seems too difficult to accomplish for me.  I built a 4.7 gcc that seems to be able to build rstan, but just crashes R when running a model. 

A byproduct is that I am able to set up the experimental clang 3.2 for windows.  But it is 32 bits.  I will post another message so it would be easily searched later.  But this clang 3.2 could not pass all our unit tests.  (And clang 3.1 has a bug when running on windows.)

--

Jiqiang 

[Bob Carpenter]

On 12/26/12 4:44 PM, Jiqiang Guo wrote:
> I hate to say I am giving up on this,

By all means, don't drive yourself crazy with builds.
It's hard to set up multiple dev environments because
you need to either change all the environment variables
and paths to switch between them or specify the path
in the compilation and linking steps. Daniel's makefile
makes the latter pretty easy, but I have no idea what
obstacles R throws down after that.

When Daniel and I are back after the break, we're going
to sit down and try to work through the parsing issues
for 4.7 GCC and all of the C++11 modes.

- Bob