read_cmdstan_csv efficiency issue for large number of parameters

[yizhang-yiz]

Describe the bug
CSV read speed drastically degrades when the number of params is large(>10^5) .

To Reproduce
@billgillespie and I ran into a model with 10^6 params and read_cmdstan_csv becomes stalled . I was able to reproduce it with the following dummy model

parameters {
  real k;
}

transformed parameters {
  real x[100000] = rep_array(k, 100000);
}

model {
  k ~ normal(0, 1);
}

The output file contains only one sample:

./dummy sample algorithm=fixed_param num_warmup=0 num_samples=1 output file=dummy.csv

rstan::read_stan_csv vs cmdstanr::read_cmdstan_csv:

> system.time(fit <- cmdstanr::read_cmdstan_csv("dummy.csv"))
   user  system elapsed 
 55.873  26.429  67.056 
> system.time(fit <- rstan::read_stan_csv("dummy.csv"))
   user  system elapsed 
  5.449   0.341   5.847 

Expected behavior
More realistic reading speed in this scenario.

Operating system
macos mojave

> sessionInfo()
R version 3.6.2 (2019-12-12)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Mojave 10.14.6

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
[1] compiler_3.6.2 tools_3.6.2    tcltk_3.6.2   

CmdStanR version number

> library("cmdstanr")
This is cmdstanr version 0.1.3
- Online documentation and vignettes at mc-stan.org/cmdstanr
- Use set_cmdstan_path() to set the path to CmdStan
- Use install_cmdstan() to install CmdStan

[mike-lawrence]

You're not seeing a message that says "Fast CSV reading with vroom::vroom() failed." are you?

If not, then reading the CSVs already employs the current fastest csv-reading method, vroom. Someday cmdstan may move to a proper binary output format and be ridiculously faster (see here for discussions on this), but for now we're stuck with the current speeds.

[mike-lawrence]

Ah, slight correction: I guess data.table::fread() is slightly faster than vroom::vroom() (see here), and while it's possible that cmdstanr might move to the former for it's lighter dependency requirements, that probably won't yield noticeable speed improvements.

[rok-cesnovar]

Thanks for the report @yizhang-yiz and sorry for the issues. Will take a look this weekend.

[yizhang-yiz]

You're not seeing a message that says "Fast CSV reading with vroom::vroom() failed." are you?

No.

I guess data.table::fread() is slightly faster than vroom::vroom()

@billgillespie just tested it and fread is much faster.

[mike-lawrence]

@billgillespie just tested it and fread is much faster.

Cool! Did you guys do this test in a fork of cmdstanr? If so, feel free to do a pull request!

[dpastoor]

FYI, from vroom itself:

https://vroom.r-lib.org/articles/benchmarks.html#all-numeric-data

All numeric data is really a worst case scenario for vroom. The index takes about as much memory as the parsed data. Also because parsing doubles can be done quickly in parallel and text representations of doubles are only ~25 characters at most there isn’t a great deal of savings for delayed parsing.

For these reasons (and because the data.table implementation is very fast) vroom is a bit slower than fread for pure numeric data.

[rok-cesnovar]

Sorry for the radio silence on this issue. After debugging this issue for some time I took another look at fread and some preliminary tests show that we can just replace vroom with data.table and its fread function. Either data.table was updated recently to play nicer with Stan CSV or I was sloppy in my initial testing of fread.

this will get sorted by replacing vroom with fread which will also close #262 and #198

[rok-cesnovar]

Can someone on a mac try the following:

library(cmdstanr)
library(data.table)

model_code <- "
parameters {
  real k;
}

transformed parameters {
  real x[100000] = rep_array(k, 100000);
}

model {
  k ~ normal(0, 1);
}
"
model_file <- write_stan_file(model_code)
mod <- cmdstan_model(model_file)
fit <- mod$sample(iter_sampling = 1, iter_warmup = 200, validate_csv = FALSE, chains = 4)

a <- Sys.time()
t1 <- rstan::read_stan_csv(fit$output_files())
print(Sys.time()-a)

a <- Sys.time()
t2 <- list()
for (f in fit$output_files()) {
  t2[[f]] <- fread(cmd= paste0("grep -v '^#' ", f))
}
print(Sys.time()-a)

The grep thing is ugly but necessary as fread does not support comment.char. There is an equivalent Windows approach to this. On Linux this is really fast (the bottom approach).

[yizhang-yiz]

Running MCMC with 4 sequential chains...

Chain 1 Iteration:   1 / 201 [  0%]  (Warmup) 
Chain 1 Iteration: 100 / 201 [ 49%]  (Warmup) 
Chain 1 Iteration: 200 / 201 [ 99%]  (Warmup) 
Chain 1 Iteration: 201 / 201 [100%]  (Sampling) 
Chain 1 finished in 0.2 seconds.
Chain 2 Iteration:   1 / 201 [  0%]  (Warmup) 
Chain 2 Iteration: 100 / 201 [ 49%]  (Warmup) 
Chain 2 Iteration: 200 / 201 [ 99%]  (Warmup) 
Chain 2 Iteration: 201 / 201 [100%]  (Sampling) 
Chain 2 finished in 0.2 seconds.
Chain 3 Iteration:   1 / 201 [  0%]  (Warmup) 
Chain 3 Iteration: 100 / 201 [ 49%]  (Warmup) 
Chain 3 Iteration: 200 / 201 [ 99%]  (Warmup) 
Chain 3 Iteration: 201 / 201 [100%]  (Sampling) 
Chain 3 finished in 0.2 seconds.
Chain 4 Iteration:   1 / 201 [  0%]  (Warmup) 
Chain 4 Iteration: 100 / 201 [ 49%]  (Warmup) 
Chain 4 Iteration: 200 / 201 [ 99%]  (Warmup) 
Chain 4 Iteration: 201 / 201 [100%]  (Sampling) 
Chain 4 finished in 0.2 seconds.

All 4 chains finished successfully.
Mean chain execution time: 0.2 seconds.
Total execution time: 0.9 seconds.
> Time difference of 8.561616 secs
> + + > Time difference of 1.311881 secs
> 

[rok-cesnovar]

Cool, thanks.

At least for unix systems then this looks promising. For windows we can use findstr.
Hopefully that will work.

No comment.char in fread is a longstanding issue (Rdatatable/data.table#856 (comment)) and there does not seem to be much traction there, so this grep/findstr thing will have to do.

[spinkney]

fwiw, it was taking a super long time on a recent model for fit$summary() that I just stopped it. Running the above from Rok on my model yields

> a <- Sys.time()
> t1 <- rstan::read_stan_csv(fit_vi$output_files())
Warning messages:
1: In parse_stancsv_comments(comments) :
  line with "Elapsed Time" not found
2: In validityMethod(object) :
  The following variables have undefined values:  log_p__. Many subsequent functions will not work correctly.
> print(Sys.time()-a)
Time difference of 3.083478 mins
> 
> a <- Sys.time()
> t2 <- list()
> for (f in fit_vi$output_files()) {
+   t2[[f]] <- fread(cmd= paste0("grep -v '^#' ", f))
+ }
> print(Sys.time()-a)
Time difference of 20.21756 secs

9x faster! On the base ec2 (2 cores, 8 gb mem) linux ubuntu 20.04