Write dynamically determined number of warmup to csv output

I’m seeking help from cmdstan and I/O people on how to write num_warmup in csv output in a later stage, instead of during initialization. This came up when I experiment on terminating warmup dynamically so the actually num_warmup would be only available at the end of warmup.

Intuition says a rewind & replace on the fstream, but currently I’m butchering it with a system call to awk. What would be a proper way? @mitzimorris @bbbales2 @Stevo15025 @rok_cesnovar

Well the definition of “num_warmup” changed. It’s now the maximum amount of warmup.

I think this’d be the responsibility of whatever is parsing the output to figure out how much warmup is there. Presumably it’d just count the lines.

we can either burden sampler stream writing or interface parsing, either works for me.

I would first remove the num_warmup line from the existing metadata written to the file and then add a similar line

as we have for the mass matrix

# "new warmup algorithm name" finished
# num_warmup = 123

This should be easy to add. Let me know if you need help.
So that would be

# stan_version_major = 2
# stan_version_minor = 18
# stan_version_patch = 1
# model = bernoulli_model
# method = sample (Default)
#   sample
#     num_samples = 1000 (Default)
#     num_warmup = 1000 (Default)
#     save_warmup = 0 (Default)
#     thin = 1 (Default)
#     adapt
#       engaged = 1 (Default)
#       gamma = 0.050000000000000003 (Default)
#       delta = 0.80000000000000004 (Default)
#       kappa = 0.75 (Default)
#       t0 = 10 (Default)
#       init_buffer = 75 (Default)
#       term_buffer = 50 (Default)
#       window = 25 (Default)
#     algorithm = hmc (Default)
#       hmc
#         engine = nuts (Default)
#           nuts
#             max_depth = 10 (Default)
#         metric = diag_e (Default)
#         metric_file =  (Default)
#         stepsize = 1 (Default)
#         stepsize_jitter = 0 (Default)
# id = 0 (Default)
# data
#   file = examples/bernoulli/bernoulli.data.json
# init = 2 (Default)
# random
#   seed = 2057108222
# output
#   file = output.csv (Default)
#   diagnostic_file =  (Default)
#   refresh = 100 (Default)
lp__,accept_stat__,stepsize__,treedepth__,n_leapfrog__,divergent__,energy__,theta
# Adaptation terminated
# Step size = 0.822884
# Diagonal elements of inverse mass matrix:
# 0.417943
-7.1053,0.953597,0.822884,2,3,0,7.96837,0.364082
-7.27881,0.989561,0.822884,3,7,0,7.30662,0.390752
-7.11597,1,0.822884,1,1,0,7.27771,0.365875
-6.98136,1,0.822884,3,7,0,7.08989,0.341094
-6.84576,1,0.822884,1,1,0,6.95638,0.307776
-6.79428,0.903813,0.822884,3,7,0,7.65851,0.289241
-6.88147,0.944843,0.822884,1,3,0,7.24226,0.189408
-6.80644,0.967211,0.822884,3,7,0,7.28921,0.294255
-6.78961,1,0.822884,1,1,0,6.80642,0.287149
-6.75497,1,0.822884,3,7,0,6.78592,0.26492
-6.74848,0.999216,0.822884,2,7,0,6.76339,0.246229
-6.89863,0.946399,0.822884,2,3,0,7.28267,0.322389
-6.91067,0.997765,0.822884,1,1,0,6.93378,0.325363
-7.00688,0.982414,0.822884,2,7,0,7.13223,0.167972
-6.94034,1,0.822884,2,3,0,7.00674,0.178291
-6.88987,1,0.822884,1,1,0,6.94007,0.187666
-7.09596,0.979745,0.822884,2,3,0,7.1828,0.156488
-6.9546,1,0.822884,1,1,0,7.07972,0.175914
-6.82338,0.97417,0.822884,2,3,0,7.29343,0.203688
-6.85475,0.99835,0.822884,3,7,0,6.85771,0.195414
-6.79183,0.989882,0.822884,3,7,0,6.98517,0.288156
-6.7733,1,0.822884,1,1,0,6.78926,0.278784
-7.42424,0.795049,0.822884,1,3,0,8.48405,0.126167
-7.7722,0.969716,0.822884,1,1,0,7.77409,0.104332
-8.11457,0.968137,0.822884,3,7,0,8.604,0.0882366
-7.98461,1,0.822884,1,1,0,8.1808,0.0938735
-8.41624,0.971264,0.822884,1,1,0,8.42639,0.0768887
-8.28891,0.988322,0.822884,2,3,0,9.01846,0.0814136
-8.15577,1,0.822884,1,1,0,8.36528,0.0865526
-6.81704,0.998172,0.822884,3,7,0,8.06708,0.298239
-6.78833,1,0.822884,1,1,0,6.81315,0.286557
-6.98102,0.914352,0.822884,1,3,0,7.47389,0.171763
-6.91594,1,0.822884,1,1,0,6.97896,0.182614
-6.75555,0.989468,0.822884,2,3,0,7.06065,0.234878
-7.00202,0.979195,0.822884,2,7,0,7.00226,0.345258
-7.02305,0.998654,0.822884,2,3,0,7.06316,0.349341
-6.74824,1,0.822884,2,3,0,6.98094,0.247405
-6.80348,0.989164,0.822884,1,3,0,6.84041,0.209986
-6.84824,0.969801,0.822884,1,3,0,7.15788,0.308535
-7.02817,0.878704,0.822884,1,3,0,7.82137,0.165024
-7.02031,0.91386,0.822884,1,3,0,8.23531,0.348818
# "new warmup algorithm name" finished
# num_warmup = 123
-6.97083,1,0.822884,2,3,

Yeah that’s the easiest on Stan side, but requires changes on read_stan_csv, does it?

Well, as long as its written before

# Adaptation terminated
# Step size = 0.822884

you are fine. And I guess this also finishes before adaptation (yes?). The following example returns the correct value if I do:
a<- read_stan_csv("~/Desktop/cmdstan/src/test/test-models/gq_model_output.csv")

a@sim$iter
[1] 1112

# stan_version_major = 2
# stan_version_minor = 18
# stan_version_patch = 1
# model = bernoulli_model
# method = sample (Default)
#   sample
#     num_samples = 1000 (Default)
#     save_warmup = 0 (Default)
#     thin = 1 (Default)
#     adapt
#       engaged = 1 (Default)
#       gamma = 0.050000000000000003 (Default)
#       delta = 0.80000000000000004 (Default)
#       kappa = 0.75 (Default)
#       t0 = 10 (Default)
#       init_buffer = 75 (Default)
#       term_buffer = 50 (Default)
#       window = 25 (Default)
#     algorithm = hmc (Default)
#       hmc
#         engine = nuts (Default)
#           nuts
#             max_depth = 10 (Default)
#         metric = diag_e (Default)
#         metric_file =  (Default)
#         stepsize = 1 (Default)
#         stepsize_jitter = 0 (Default)
# id = 0 (Default)
# data
#   file = examples/bernoulli/bernoulli.data.json
# init = 2 (Default)
# random
#   seed = 2057108222
# output
#   file = output.csv (Default)
#   diagnostic_file =  (Default)
#   refresh = 100 (Default)
lp__,accept_stat__,stepsize__,treedepth__,n_leapfrog__,divergent__,energy__,theta
# num_warmup = 112
# Adaptation terminated
# Step size = 0.822884
# Diagonal elements of inverse mass matrix:
# 0.417943
-7.1053,0.953597,0.822884,2,3,0,7.96837,0.364082
-7.27881,0.989561,0.822884,3,7,0,7.30662,0.390752
-7.11597,1,0.822884,1,1,0,7.27771,0.365875
-6.98136,1,0.822884,3,7,0,7.08989,0.341094
-6.84576,1,0.822884,1,1,0,6.95638,0.307776
-6.79428,0.903813,0.822884,3,7,0,7.65851,0.289241
-6.88147,0.944843,0.822884,1,3,0,7.24226,0.189408
-6.80644,0.967211,0.822884,3,7,0,7.28921,0.294255
-6.78961,1,0.822884,1,1,0,6.80642,0.287149
-6.75497,1,0.822884,3,7,0,6.78592,0.26492
-6.74848,0.999216,0.822884,2,7,0,6.76339,0.246229
-6.89863,0.946399,0.822884,2,3,0,7.28267,0.322389
-6.91067,0.997765,0.822884,1,1,0,6.93378,0.325363
-7.00688,0.982414,0.822884,2,7,0,7.13223,0.167972
-6.94034,1,0.822884,2,3,0,7.00674,0.178291
-6.88987,1,0.822884,1,1,0,6.94007,0.187666
-7.09596,0.979745,0.822884,2,3,0,7.1828,0.156488
-6.9546,1,0.822884,1,1,0,7.07972,0.175914
-6.82338,0.97417,0.822884,2,3,0,7.29343,0.203688
-6.85475,0.99835,0.822884,3,7,0,6.85771,0.195414

fantastic, thanks a lot.

1 Like

What happens if different chains have different number of warmup iterations? It may wreak havoc downstream with interfaces. Nothing can get released until we know this won’t crash RStan, PyStan and CmdStan.

The way that the algorithm I’m implementing stops warmup guarantees same number of warmup iterations. Only that it’s no longer known a priori.

1 Like

This may be a problem for RStan, which allocates ahead of time. It might be a problem with all the comments and the CSV reader. Our CSV format is very brittle.

In that case there’s going to be even more rigging before cross-chain landing(if it can land).

I don’t think there’s any fundamental reason why it couldn’t. The changes shouldn’t be that hard. I think it’ll be worth it to get faster and/or more robust adaptation.

Eventually we’ll have to deal with ragged chains. The only obstacle is that we have a lot of code assuming rectangular structures.

1 Like