R crashes when using rstan

[Michael Braun]

 I wanted to let you know that my first attempt to use Stan (via rstan) has not gone very well.  The result is an Abort trap, and all I wanted to do was try a simple hierarchical normal regression.  

Below, I've pasted the code, and then the verbose output.

Any idea what's wrong?

Thanks,

Michael




Here's the stan file, hmvn.stan

data{

int<lower=1> N;
int<lower=1> T;
int<lower=1> k;

matrix[T,N] Y;
matrix[k,N] X;

vector[k] mu_prior_mean;
matrix[k,k] mu_prior_cov;
real<lower=k+1> nu;
matrix[k,k] S;
}

transformed data {
int<lower=0> idx;
vector[N] xi[k];
for (i in 1:N) {
    xi[i] <- col(X,i);
}      
}

parameters {
vector[N] B[k];
vector[k] mu;
cov_matrix[k] Sig;
}

model {
Sig ~ inv_wishart(nu, S);
mu ~ multi_normal(mu_prior_mean, mu_prior_cov);
B[i] ~ multi_normal(mu, Sig);
for (i in 1:N) {
  for (t in 1:T) {
    Y[t,i] ~ normal(dot_product(xi[i],B[i]),1);
  }
} 
}

Now, the R file:


library(inline)
library(Rcpp)
library(RcppEigen)
library(mvtnorm)
library(plyr)
library(rstan)


N <- 25
T <- 5
sd1 <- .25
sd2 <- .25
k <- 4

## simulate data

X <- rbind(rep(1,N),matrix(rnorm(N*(k-1),0,1),k-1,N))
B <- laply(c(5,0,2,0),function(mu) return(rnorm(N,mu,sd1)))
y.mean <- colSums(X*B)
Y <- t(laply(y.mean, function(mu) return(rnorm(T,mu, sd2))))

## set hyperpriors

mu.prior.mean <- rep(0,k)
mu.prior.cov <- .04*diag(k)
mu.prior.chol.prec <- t(chol(mu.prior.cov))
nu <- as.numeric(k+6)
G.mean <- .10*diag(k)
S <- solve(G.mean/nu)
chol.S <- t(chol(S))


data.list <- list(N=N, T=T, k=k,
                Y=Y, X=X,
                mu_prior_mean=mu.prior.mean,
                mu_prior_cov=mu.prior.cov,
                nu=nu,
                S=S
                )

fit1 <- stan(file="hmvn.stan",
           data=data.list,
           iter=10,
           n_chains=1,
           verbose=TRUE,
           warmup=0
           )

Here's the verbose output

source("hmvn.R")

Loading required package: Matrix
Loading required package: lattice
rstan (Version 1.0.0, packaged: 2012-08-31 00:25:52 UTC, GitRev: 7e72f4e8befe)

TRANSLATING MODEL 'anon_model' FROM Stan CODE TO C++ CODE NOW.
successful of parsing the Stan model 'anon_model'.
COMPILING THE C++ CODE FOR MODEL 'anon_model' NOW.

setting environment variables: 

PKG_LIBS =  /Library/Frameworks/R.framework/Versions/2.15/Resources/library/Rcpp/lib/x86_64/libRcpp.a "/Library/Frameworks/R.framework/Versions/2.15/Resources/library/rstan/libstan/x86_64/libstan.a"
PKG_CPPFLAGS =   -I"/Library/Frameworks/R.framework/Versions/2.15/Resources/library/rstan/include//stansrc"  -I"/Library/Frameworks/R.framework/Versions/2.15/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/2.15/Resources/library/RcppEigen/include//unsupported"  -I"/Library/Frameworks/R.framework/Versions/2.15/Resources/library/rstan/include//stanlib/boost_1.51.0" -I"/Library/Frameworks/R.framework/Versions/2.15/Resources/library/rstan/include"

LinkingTo : Rcpp

CLINK_CPPFLAGS =  -I"/Library/Frameworks/R.framework/Versions/2.15/Resources/library/Rcpp/include" 

Program source :

<<<  Not including 423 lines of source code >>>


Compilation argument:
/Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB file135956348088.cpp 2> file135956348088.cpp.err.txt 
g++ -m64 -std=c++03 -I/Library/Frameworks/R.framework/Resources/include -I/Library/Frameworks/R.framework/Resources/include/x86_64 -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/2.15/Resources/library/rstan/include//stansrc"  -I"/Library/Frameworks/R.framework/Versions/2.15/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/2.15/Resources/library/RcppEigen/include//unsupported"  -I"/Library/Frameworks/R.framework/Versions/2.15/Resources/library/rstan/include//stanlib/boost_1.51.0" -I"/Library/Frameworks/R.framework/Versions/2.15/Resources/library/rstan/include" -I/usr/local/include -I"/Library/Frameworks/R.framework/Versions/2.15/Resources/library/Rcpp/include"   -fPIC  -O2 -mtune=corei7 -msse4.2  -c file135956348088.cpp -o file135956348088.o
g++ -m64 -std=c++03 -dynamiclib -Wl,-headerpad_max_install_names -undefined dynamic_lookup -single_module -multiply_defined suppress -L/opt/intel/composer_xe_2011_sp1.10.328/mkl/lib -lmkl_rt -lpthread -liomp5 -o file135956348088.so file135956348088.o /Library/Frameworks/R.framework/Versions/2.15/Resources/library/Rcpp/lib/x86_64/libRcpp.a /Library/Frameworks/R.framework/Versions/2.15/Resources/library/rstan/libstan/x86_64/libstan.a -F/Library/Frameworks/R.framework/.. -framework R -Wl,-framework -Wl,CoreFoundation
R(4953,0x7fff78edc180) malloc: *** error for object 0x1000000000000000: pointer being freed was not allocated
*** set a breakpoint in malloc_error_break to debug
Abort trap: 6

[Ben Goodrich]

Hi Michael,

On Saturday, September 1, 2012 9:46:28 AM UTC-4, Michael Braun wrote:

Any idea what's wrong?

transformed data {
int<lower=0> idx;
vector[N] xi[k];
for (i in 1:N) {
    xi[i] <- col(X,i);
}      
}

The idx variable appears to be unused, which may cause the abort you are seeing.
 

model {
Sig ~ inv_wishart(nu, S);
mu ~ multi_normal(mu_prior_mean, mu_prior_cov);
B[i] ~ multi_normal(mu, Sig);

This won't compile; did you mean to put B[i] ~ multi_normal(mu, Sig) inside the for(i in 1:N) loop?
 
Ben

P.S. You can attach the .stan and .R files via email or the Google Groups interface.

[Michael Braun]

Ben:

I removed the idx variable (which was in the file I actually ran), and moved the B[i] line (which was not, but was out of order in the email because I got a little fancy with cutting and pasting).  I have attached the two files that are exactly what generate the following:

source("hmvn.R")

rstan (Version 1.0.0, packaged: 2012-08-31 00:25:52 UTC, GitRev: 7e72f4e8befe)

TRANSLATING MODEL 'anon_model' FROM Stan CODE TO C++ CODE NOW.

COMPILING THE C++ CODE FOR MODEL 'anon_model' NOW.

R(54576,0x7fff78edc180) malloc: *** error for object 0xf100000d: pointer being freed was not allocated

*** set a breakpoint in malloc_error_break to debug

Abort trap: 6

The hmvn.stan file does compile using make from a shell prompt.  I cannot try running it from the shell because I cannot find a function that converts an R list of variables into the dump data format that stan seems to want (any suggestions on that?).

Thanks,

Michael

On Saturday, September 1, 2012 9:46:28 AM UTC-4, Michael Braun wrote:

[Jiqiang Guo]

For this purpose, R function dump or function stan_rdump in rstan package can be used.  

--

Jiqiang

[Ben Goodrich]

Hi Michael,

On Saturday, September 1, 2012 3:46:17 PM UTC-4, Michael Braun wrote:

The hmvn.stan file does compile using make from a shell prompt.  I cannot try running it from the shell because I cannot find a function that converts an R list of variables into the dump data format that stan seems to want (any suggestions on that?).

In R, see help(dump). It takes a character vector of object names. On the command-line, be sure to compile with -O0 (no optimization) and -g (debug symbols). Then, you should be able to execute the binary under gdb.

Ben

[Michael Braun]

Fixed the problem.  Once I exported the data, and ran stan from a shell, I got the following helpful message:

Exception: INDEX OPERATOR [] OUT OF BOUNDS; index=5; lower bound=1; upper bound=4; index position=1; xi

Diagnostic information: 

Dynamic exception type: std::out_of_range

std::exception::what: INDEX OPERATOR [] OUT OF BOUNDS; index=5; lower bound=1; upper bound=4; index position=1; xi

This told me that I was interpreting the statement

 vector[N] xi[k]

incorrectly.  I thought it meant a vector of N k-dimensional arrays, when in fact it is a vector of k N-dimensional arrays. Once I fixed that, it all worked.

It might be useful for those messages to pass through the rstan package somehow.

On Saturday, September 1, 2012 9:46:28 AM UTC-4, Michael Braun wrote:

[Ben Goodrich]

On Saturday, September 1, 2012 5:13:45 PM UTC-4, Michael Braun wrote:

It might be useful for those messages to pass through the rstan package somehow.

 

Glad it is working now.  I think those messages "are" being passed to the R console, but in your case, R crashed before it flushed the message through.

Ben

[Jiqiang Guo]

Thanks, Ben.

Jiqiang

On Sat, Sep 1, 2012 at 5:39 PM, Ben Goodrich <goodri...@gmail.com> wrote:

I don't think so because in this case, the guy had his rows and columns reversed. Ideally, there might be some exceptions policy that would cause the sampler to just break out of its loop over the iterations. But we would have to ask Daniel about that.

Ben

On Sat, Sep 1, 2012 at 5:36 PM, Jiqiang Guo <guo...@gmail.com> wrote:

Ben,

I think you are right since sometimes we can see the error message being popped up in R console.  So can we do anything preventing R from crashing in this case? 

Jiqiang

[Michael Braun]

Yes, it was my mistake that caused the crash.  But I think what made it frustrating is that there was no indication as to why it crashed (no error message), no obvious way to debug it that an end-user would be able to understand (not all R users know C++ and gdb), and that it actually *crashed R*.

On Sep 1, 2012, at 5:42 PM, Jiqiang Guo <guo...@gmail.com>

[Ben Goodrich]

We know it is frustrating (I've crashed R from stan() many times myself). But we don't really know what to do about it. In this case, trying to assign out of range to a std::vector causes abort() to get called by code within the std library. Possibly we can catch the SIGABRT and flush the messages, but it is still going to crash R.

Ben

[Bob Carpenter]

On 9/1/12 5:13 PM, Michael Braun wrote:
> Fixed the problem. Once I exported the data, and ran stan from a shell, I got the following helpful message:
>
> Exception: INDEX OPERATOR [] OUT OF BOUNDS; index=5; lower bound=1; upper bound=4; index position=1; xi
>
> Diagnostic information:
> Dynamic exception type: std::out_of_range
> std::exception::what: INDEX OPERATOR [] OUT OF BOUNDS; index=5; lower bound=1; upper bound=4; index position=1; xi
>
> This told me that I was interpreting the statement
>
> vector[N] xi[k]
>
> incorrectly. I thought it meant a vector of N k-dimensional arrays, when in fact it is a vector of k N-dimensional
> arrays. Once I fixed that, it all worked.

Just the kind of thing we were worried people would confuse.
We just didn't want to go all the way to the more verbose
notation:

matrix<rows=M,columns=N> a;
vector<columns=N> xi[k];
row_vector<rows=M> psi[j];

or

matrix<rows=M,columns=N> a;
vector<size=N> xi[k];
row_vector<size=N> psi[j];

> It might be useful for those messages to pass through the rstan package somehow.

Absolutely. We'll work on making sure the errors
propagate appropriately.

- Bob