Stan  1.0
probability, sampling & optimization
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Pages
command.hpp
Go to the documentation of this file.
1 #ifndef __STAN__GM__COMMAND_HPP__
2 #define __STAN__GM__COMMAND_HPP__
3 
4 #include <cmath>
5 #include <cstddef>
6 #include <boost/math/special_functions/fpclassify.hpp>
7 #include <boost/date_time/posix_time/posix_time_types.hpp>
8 #include <iomanip>
9 #include <iostream>
10 #include <fstream>
11 #include <sstream>
12 #include <vector>
13 #include <boost/random/additive_combine.hpp> // L'Ecuyer RNG
14 #include <boost/random/mersenne_twister.hpp>
15 #include <boost/random/uniform_01.hpp>
16 #include <boost/random/uniform_real_distribution.hpp>
17 #include <stan/version.hpp>
18 #include <stan/io/cmd_line.hpp>
19 #include <stan/io/dump.hpp>
20 #include <stan/mcmc/adaptive_sampler.hpp>
21 #include <stan/mcmc/adaptive_hmc.hpp>
22 #include <stan/mcmc/hmc.hpp>
23 #include <stan/mcmc/nuts.hpp>
24 #include <stan/mcmc/nuts_diag.hpp>
25 #include <stan/model/prob_grad_ad.hpp>
26 #include <stan/model/prob_grad.hpp>
27 #include <stan/mcmc/sampler.hpp>
28 
29 namespace stan {
30 
31  namespace gm {
32 
33 
34  void print_nuts_help(std::string cmd) {
35  using stan::io::print_help_option;
36 
37  std::cout << std::endl;
38  std::cout << "Compiled Stan Graphical Model Command" << std::endl;
39  std::cout << std::endl;
40 
41  std::cout << "USAGE: " << cmd << " [options]" << std::endl;
42  std::cout << std::endl;
43 
44  std::cout << "OPTIONS:" << std::endl;
45  std::cout << std::endl;
46 
47  print_help_option(&std::cout,
48  "help","",
49  "Display this information");
50 
51  print_help_option(&std::cout,
52  "data","file",
53  "Read data from specified dump-format file",
54  "required if model declares data");
55 
56  print_help_option(&std::cout,
57  "init","file",
58  "Use initial values from specified file or zero values if <file>=0",
59  "default is random initialization");
60 
61  print_help_option(&std::cout,
62  "samples","file",
63  "File into which samples are written",
64  "default = samples.csv");
65 
66  print_help_option(&std::cout,
67  "append_samples","",
68  "Append samples to existing file if it exists",
69  "does not write header in append mode");
70 
71  print_help_option(&std::cout,
72  "seed","int",
73  "Random number generation seed",
74  "default = randomly generated from time");
75 
76  print_help_option(&std::cout,
77  "chain_id","int",
78  "Markov chain identifier",
79  "default = 1");
80 
81  print_help_option(&std::cout,
82  "iter","+int",
83  "Total number of iterations, including warmup",
84  "default = 2000");
85 
86  print_help_option(&std::cout,
87  "warmup","+int",
88  "Discard the specified number of initial samples",
89  "default = iter / 2");
90 
91  print_help_option(&std::cout,
92  "thin","+int",
93  "Period between saved samples after warm up",
94  "default = max(1, floor(iter - warmup) / 1000)");
95 
96  print_help_option(&std::cout,
97  "refresh","int",
98  "Period between samples updating progress report print (0 for no printing)",
99  "default = max(1,iter/200))");
100 
101  print_help_option(&std::cout,
102  "leapfrog_steps","int",
103  "Number of leapfrog steps; -1 for no-U-turn adaptation",
104  "default = -1");
105 
106  print_help_option(&std::cout,
107  "max_treedepth","int",
108  "Limit NUTS leapfrog steps to 2^max_tree_depth; -1 for no limit",
109  "default = 10");
110 
111  print_help_option(&std::cout,
112  "epsilon","float",
113  "Initial value for step size, or -1 to set automatically",
114  "default = -1");
115 
116  print_help_option(&std::cout,
117  "epsilon_pm","[0,1]",
118  "Sample epsilon +/- epsilon * epsilon_pm",
119  "default = 0.0");
120 
121  print_help_option(&std::cout,
122  "equal_step_sizes","",
123  "Use same step size for every parameter with NUTS",
124  "default is to estimate varying step sizes during warmup");
125 
126  print_help_option(&std::cout,
127  "delta","+float",
128  "Initial step size for step-size adaptation",
129  "default = 0.5");
130 
131  print_help_option(&std::cout,
132  "gamma","+float",
133  "Gamma parameter for dual averaging step-size adaptation",
134  "default = 0.05");
135 
136  print_help_option(&std::cout,
137  "save_warmup","",
138  "Save the warmup samples");
139 
140  print_help_option(&std::cout,
141  "test_grad","",
142  "Test gradient calculations using finite differences");
143 
144  std::cout << std::endl;
145  }
146 
147  bool do_print(int n, int refresh) {
148  return (refresh > 0)
149  && (n == 0
150  || ((n + 1) % refresh == 0) );
151  }
152 
153  template <class Sampler, class Model>
154  void sample_from(Sampler& sampler,
155  bool epsilon_adapt,
156  int refresh,
157  int num_iterations,
158  int num_warmup,
159  int num_thin,
160  bool save_warmup,
161  std::ostream& sample_file_stream,
162  std::vector<double>& params_r,
163  std::vector<int>& params_i,
164  Model& model) {
165 
166  sampler.set_params(params_r,params_i);
167 
168  int it_print_width = std::ceil(std::log10(num_iterations));
169  std::cout << std::endl;
170 
171  if (epsilon_adapt)
172  sampler.adapt_on();
173  for (int m = 0; m < num_iterations; ++m) {
174  if (do_print(m,refresh)) {
175  std::cout << "Iteration: ";
176  std::cout << std::setw(it_print_width) << (m + 1)
177  << " / " << num_iterations;
178  std::cout << " [" << std::setw(3)
179  << static_cast<int>((100.0 * (m + 1))/num_iterations)
180  << "%] ";
181  std::cout << ((m < num_warmup) ? " (Adapting)" : " (Sampling)");
182  std::cout << std::endl;
183  std::cout.flush();
184  }
185  if (m < num_warmup) {
186  if (save_warmup && (m % num_thin) == 0) {
187  stan::mcmc::sample sample = sampler.next();
188 
189  // FIXME: use csv_writer arg to make comma optional?
190  sample_file_stream << sample.log_prob() << ',';
191  sampler.write_sampler_params(sample_file_stream);
192  sample.params_r(params_r);
193  sample.params_i(params_i);
194  model.write_csv(params_r,params_i,sample_file_stream);
195  } else {
196  sampler.next(); // discard
197  }
198  } else {
199  if (epsilon_adapt && sampler.adapting()) {
200  sampler.adapt_off();
201  sampler.write_adaptation_params(sample_file_stream);
202  }
203  if (((m - num_warmup) % num_thin) != 0) {
204  sampler.next();
205  continue;
206  } else {
207  stan::mcmc::sample sample = sampler.next();
208 
209  // FIXME: use csv_writer arg to make comma optional?
210  sample_file_stream << sample.log_prob() << ',';
211  sampler.write_sampler_params(sample_file_stream);
212  sample.params_r(params_r);
213  sample.params_i(params_i);
214  model.write_csv(params_r,params_i,sample_file_stream);
215  }
216  }
217  }
218  }
219 
220  void write_comment(std::ostream& o) {
221  o << "#" << std::endl;
222  }
223  template <typename M>
224  void write_comment(std::ostream& o,
225  const M& msg) {
226  o << "# " << msg << std::endl;
227  }
228  template <typename K, typename V>
229  void write_comment_property(std::ostream& o,
230  const K& key,
231  const V& val) {
232  o << "# " << key << "=" << val << std::endl;
233  }
234 
235  template <class Model>
236  int nuts_command(int argc, const char* argv[]) {
237 
238  stan::io::cmd_line command(argc,argv);
239 
240  if (command.has_flag("help")) {
241  print_nuts_help(argv[0]);
242  return 0;
243  }
244 
245  std::string data_file;
246  command.val("data",data_file);
247  std::fstream data_stream(data_file.c_str(),
248  std::fstream::in);
249  stan::io::dump data_var_context(data_stream);
250  data_stream.close();
251 
252  // FIXME: configure to use alt to std::cout for model prints
253  Model model(data_var_context, &std::cout);
254 
255  std::string sample_file = "samples.csv";
256  command.val("samples",sample_file);
257 
258  unsigned int num_iterations = 2000U;
259  command.val("iter",num_iterations);
260 
261  unsigned int num_warmup = num_iterations / 2;
262  command.val("warmup",num_warmup);
263 
264  unsigned int calculated_thin = (num_iterations - num_warmup) / 1000U;
265  unsigned int num_thin = (calculated_thin > 1) ? calculated_thin : 1U;
266  command.val("thin",num_thin);
267 
268  bool user_supplied_thin = command.has_key("thin");
269 
270  int leapfrog_steps = -1;
271  command.val("leapfrog_steps",leapfrog_steps);
272 
273  double epsilon = -1.0;
274  command.val("epsilon",epsilon);
275 
276  int max_treedepth = 10;
277  command.val("max_treedepth",max_treedepth);
278 
279  double epsilon_pm = 0.0;
280  command.val("epsilon_pm",epsilon_pm);
281 
282  bool epsilon_adapt = epsilon <= 0.0;
283 
284  bool equal_step_sizes = command.has_flag("equal_step_sizes");
285 
286  double delta = 0.5;
287  command.val("delta", delta);
288 
289  double gamma = 0.05;
290  command.val("gamma", gamma);
291 
292  int refresh = num_iterations / 200;
293  refresh = refresh <= 0 ? 1 : refresh; // just for default
294  command.val("refresh",refresh);
295 
296  unsigned int random_seed = 0;
297  if (command.has_key("seed")) {
298  bool well_formed = command.val("seed",random_seed);
299  if (!well_formed) {
300  std::string seed_val;
301  command.val("seed",seed_val);
302  std::cerr << "value for seed must be integer"
303  << "; found value=" << seed_val << std::endl;
304  return -1;
305  }
306  } else {
307  random_seed
308  = (boost::posix_time::microsec_clock::universal_time() -
309  boost::posix_time::ptime(boost::posix_time::min_date_time))
310  .total_milliseconds();
311  }
312 
313  int chain_id = 1;
314  if (command.has_key("chain_id")) {
315  bool well_formed = command.val("chain_id",chain_id);
316  if (!well_formed || chain_id < 0) {
317  std::string chain_id_val;
318  command.val("chain_id",chain_id_val);
319  std::cerr << "value for chain_id must be positive integer"
320  << "; found chain_id=" << chain_id_val
321  << std::endl;
322  return -1;
323  }
324  }
325 
326  // FASTER, but no parallel guarantees:
327  // typedef boost::mt19937 rng_t;
328  // rng_t base_rng(static_cast<unsigned int>(random_seed + chain_id - 1);
329 
330  typedef boost::ecuyer1988 rng_t;
331  rng_t base_rng(random_seed);
332  // (2**50 = 1T samples, 1000 chains)
333  static boost::uintmax_t DISCARD_STRIDE = static_cast<boost::uintmax_t>(1) << 50;
334  // DISCARD_STRIDE <<= 50;
335  base_rng.discard(DISCARD_STRIDE * (chain_id - 1));
336 
337  std::vector<int> params_i;
338  std::vector<double> params_r;
339 
340  std::string init_val;
341  // parameter initialization
342  int num_init_tries = 1; // up here for printing below
343  if (command.has_key("init")) {
344  num_init_tries = -1;
345  command.val("init",init_val);
346  if (init_val == "0") {
347  params_i = std::vector<int>(model.num_params_i(),0);
348  params_r = std::vector<double>(model.num_params_r(),0.0);
349  } else {
350  std::fstream init_stream(init_val.c_str(),std::fstream::in);
351  stan::io::dump init_var_context(init_stream);
352  init_stream.close();
353  model.transform_inits(init_var_context,params_i,params_r);
354  }
355  } else {
356  init_val = "random initialization"; // for I/O
357  // init_rng generates uniformly from -2 to 2
358  boost::random::uniform_real_distribution<double>
359  init_range_distribution(-2.0,2.0);
360  boost::variate_generator<rng_t&,
361  boost::random::uniform_real_distribution<double> >
362  init_rng(base_rng,init_range_distribution);
363 
364  params_i = std::vector<int>(model.num_params_i(),0);
365  params_r = std::vector<double>(model.num_params_r());
366 
367  // retry inits until get a finite log prob value
368  std::vector<double> init_grad;
369  static int MAX_INIT_TRIES = 100;
370  for (num_init_tries = 1; num_init_tries <= MAX_INIT_TRIES; ++num_init_tries) {
371  for (size_t i = 0; i < params_r.size(); ++i)
372  params_r[i] = init_rng();
373  // FIXME: allow config vs. std::cout
374  double init_log_prob = model.grad_log_prob(params_r,params_i,init_grad,&std::cout);
375  if (!boost::math::isfinite(init_log_prob))
376  continue;
377  for (size_t i = 0; i < init_grad.size(); ++i)
378  if (!boost::math::isfinite(init_grad[i]))
379  continue;
380  break;
381  }
382  if (num_init_tries > MAX_INIT_TRIES) {
383  std::cout << "Initialization failed after " << MAX_INIT_TRIES
384  << " attempts. "
385  << " Try specifying initial values,"
386  << " reducing ranges of constrained values,"
387  << " or reparameterizing the model."
388  << std::endl;
389  return -1;
390  }
391  }
392 
393  bool save_warmup = command.has_flag("save_warmup");
394 
395  bool append_samples = command.has_flag("append_samples");
396  std::ios_base::openmode samples_append_mode
397  = append_samples
398  ? (std::fstream::out | std::fstream::app)
399  : std::fstream::out;
400 
401 
402  std::cout << "STAN SAMPLING COMMAND" << std::endl;
403  if (data_file == "")
404  std::cout << "data = (specified model requires no data)" << std::endl;
405  else
406  std::cout << "data = " << data_file << std::endl;
407 
408  std::cout << "init = " << init_val << std::endl;
409  if (num_init_tries > 0)
410  std::cout << "init tries = " << num_init_tries << std::endl;
411 
412  std::cout << "samples = " << sample_file << std::endl;
413  std::cout << "append_samples = " << append_samples << std::endl;
414  std::cout << "save_warmup = " << save_warmup<< std::endl;
415 
416  std::cout << "seed = " << random_seed
417  << " (" << (command.has_key("seed")
418  ? "user specified"
419  : "randomly generated") << ")"
420  << std::endl;
421  std::cout << "chain_id = " << chain_id
422  << " (" << (command.has_key("chain_id")
423  ? "user specified"
424  : "default") << ")"
425  << std::endl;
426 
427  std::cout << "iter = " << num_iterations << std::endl;
428  std::cout << "warmup = " << num_warmup << std::endl;
429  std::cout << "thin = " << num_thin
430  << (user_supplied_thin ? " (user supplied)" : " (default)")
431  << std::endl;
432 
433  std::cout << "equal_step_sizes = " << equal_step_sizes << std::endl;
434  std::cout << "leapfrog_steps = " << leapfrog_steps << std::endl;
435  std::cout << "max_treedepth = " << max_treedepth << std::endl;;
436  std::cout << "epsilon = " << epsilon << std::endl;;
437  std::cout << "epsilon_pm = " << epsilon_pm << std::endl;;
438  std::cout << "delta = " << delta << std::endl;
439  std::cout << "gamma = " << gamma << std::endl;
440 
441  if (command.has_flag("test_grad")) {
442  std::cout << std::endl << "TEST GRADIENT MODE" << std::endl;
443  return model.test_gradients(params_r,params_i);
444  }
445 
446  std::fstream sample_stream(sample_file.c_str(),
447  samples_append_mode);
448 
449  write_comment(sample_stream,"Samples Generated by Stan");
450  write_comment(sample_stream);
451  write_comment_property(sample_stream,"stan_version_major",stan::MAJOR_VERSION);
452  write_comment_property(sample_stream,"stan_version_minor",stan::MINOR_VERSION);
453  write_comment_property(sample_stream,"stan_version_patch",stan::PATCH_VERSION);
454  write_comment_property(sample_stream,"data",data_file);
455  write_comment_property(sample_stream,"init",init_val);
456  write_comment_property(sample_stream,"append_samples",append_samples);
457  write_comment_property(sample_stream,"save_warmup",save_warmup);
458  write_comment_property(sample_stream,"seed",random_seed);
459  write_comment_property(sample_stream,"chain_id",chain_id);
460  write_comment_property(sample_stream,"iter",num_iterations);
461  write_comment_property(sample_stream,"warmup",num_warmup);
462  write_comment_property(sample_stream,"thin",num_thin);
463  write_comment_property(sample_stream,"equal_step_sizes",equal_step_sizes);
464  write_comment_property(sample_stream,"leapfrog_steps",leapfrog_steps);
465  write_comment_property(sample_stream,"max_treedepth",max_treedepth);
466  write_comment_property(sample_stream,"epsilon",epsilon);
467  write_comment_property(sample_stream,"epsilon_pm",epsilon_pm);
468  write_comment_property(sample_stream,"delta",delta);
469  write_comment_property(sample_stream,"gamma",gamma);
470  write_comment(sample_stream);
471 
472  if (leapfrog_steps < 0 && !equal_step_sizes) {
473  // NUTS II (with varying step size estimation during warmup)
474  stan::mcmc::nuts_diag<rng_t> nuts2_sampler(model,
475  max_treedepth, epsilon,
476  epsilon_pm, epsilon_adapt,
477  delta, gamma,
478  base_rng, &params_r,
479  &params_i);
480 
481  // cut & paste (see below) to enable sample-specific params
482  if (!append_samples) {
483  sample_stream << "lp__,"; // log probability first
484  nuts2_sampler.write_sampler_param_names(sample_stream);
485  model.write_csv_header(sample_stream);
486  }
487  nuts2_sampler.set_error_stream(std::cout); // cout intended
488  nuts2_sampler.set_output_stream(std::cout);
489 
490  sample_from(nuts2_sampler,epsilon_adapt,refresh,
491  num_iterations,num_warmup,num_thin,save_warmup,
492  sample_stream,params_r,params_i,
493  model);
494 
495  } else if (leapfrog_steps < 0 && equal_step_sizes) {
496 
497  // NUTS I (equal step sizes)
498  stan::mcmc::nuts<rng_t> nuts_sampler(model,
499  max_treedepth, epsilon,
500  epsilon_pm, epsilon_adapt,
501  delta, gamma,
502  base_rng, &params_r,
503  &params_i);
504 
505  nuts_sampler.set_error_stream(std::cout);
506  nuts_sampler.set_output_stream(std::cout); // cout intended
507  // cut & paste (see below) to enable sample-specific params
508  if (!append_samples) {
509  sample_stream << "lp__,"; // log probability first
510  nuts_sampler.write_sampler_param_names(sample_stream);
511  model.write_csv_header(sample_stream);
512  }
513 
514  sample_from(nuts_sampler,epsilon_adapt,refresh,
515  num_iterations,num_warmup,num_thin,save_warmup,
516  sample_stream,params_r,params_i,
517  model);
518 
519  } else {
520 
521  // STANDARD HMC
522  stan::mcmc::adaptive_hmc<rng_t> hmc_sampler(model,
523  leapfrog_steps,
524  epsilon, epsilon_pm, epsilon_adapt,
525  delta, gamma,
526  base_rng, &params_r,
527  &params_i);
528 
529  hmc_sampler.set_error_stream(std::cout); // intended
530  hmc_sampler.set_output_stream(std::cout);
531  // cut & paste (see above) to enable sample-specific params
532  if (!append_samples) {
533  sample_stream << "lp__,"; // log probability first
534  hmc_sampler.write_sampler_param_names(sample_stream);
535  model.write_csv_header(sample_stream);
536  }
537 
538  sample_from(hmc_sampler,epsilon_adapt,refresh,
539  num_iterations,num_warmup,num_thin,save_warmup,
540  sample_stream,params_r,params_i,
541  model);
542  }
543 
544  sample_stream.close();
545  std::cout << std::endl << std::endl;
546  return 0;
547  }
548 
549  } // namespace prob
550 
551 
552 } // namespace stan
553 
554 #endif
     [ Stan Home Page ] © 2011–2012, Stan Development Team.