StanFitter.cxx
Go to the documentation of this file.
1 #include <string>
2 #include <vector>
3 
4 // Stan's dependencies trigger some warnings...
5 #pragma GCC diagnostic push
6 #if __GNUC__ >= 6
7 #pragma GCC diagnostic ignored "-Wmisleading-indentation"
8 #endif
9 #pragma GCC diagnostic ignored "-Wunused-function"
12 #include "stan/io/reader.hpp"
13 #include "stan/io/writer.hpp"
18 #pragma GCC diagnostic pop
19 
20 #include "CAFAna/Fit/StanFitter.h"
21 
22 #include "CAFAna/Core/EigenUtils.h"
23 #include "CAFAna/Core/Utilities.h"
24 
26 
27 #include "OscLib/IOscCalc.h"
28 #include "OscLib/OscCalc.h"
29 #include "OscLib/OscCalcDMP.h"
30 #include "OscLib/OscCalcPMNS.h"
31 #include "OscLib/OscCalcPMNSOpt.h"
32 #include "OscLib/OscCalcAnalytic.h"
33 
34 #include "Utilities/func/MathUtil.h"
36 
37 // these will come in handy below
40 namespace ana
41 {
42 
43  //----------------------------------------------------------------------
45  std::vector<const IFitVar *> vars,
46  std::vector<const ISyst *> systs)
47  : IFitter(vars, systs),
48  stan::model::prob_grad(fVars.size() + fSysts.size()),
49  fCalc(nullptr),
50  fExpt(expt),
51  fMCMCSamples(vars, systs),
52  fMCMCWarmup(vars, systs)
53  {}
54 
55  //----------------------------------------------------------------------
58  const SystShifts &systSeed) const
59  {
60  std::vector<double> vals;
61  std::vector<std::string> names;
63  std::cout << "Building init context for StanFitter object. Var vals:" << std::endl;
64  for (const IFitVar *v: fVars)
65  {
67  std::cout << " " << v->ShortName() << " = " << v->GetValue(seed) << std::endl;
68 
69  vals.push_back(util::GetValAs<double>(v->GetValue(seed)));
70  names.push_back(v->ShortName());
71  }
72  // One way this can go wrong is if two variables have the same ShortName
73  assert(vals.size() == fVars.size());
74  for (const ISyst *s: fSysts)
75  {
76  vals.push_back(util::GetValAs<double>(systSeed.GetShift<stan::math::var>(s)));
77  names.push_back(s->ShortName());
78  }
79  // One way this can go wrong is if two variables have the same ShortName
80  assert(vals.size() == (fVars.size() + fSysts.size()) && vals.size() == names.size());
81 
82  // in principle you can pass 'array' vars,
83  // and they're read from the vector of doubles
84  // by counting off var_sizes[var idx] doubles from vals.
85  // each of our params stands on its own though so we don't need that complication
86  std::vector<std::vector<std::size_t>> var_sizes(names.size(), std::vector<size_t>{1});
87  stan::io::array_var_context ret(names, vals, var_sizes);
88 
89  return std::move(ret);
90  }
91 
92  //----------------------------------------------------------------------
93  void StanFitter::constrained_param_names(std::vector<std::string> &param_names, bool include_tparams,
94  bool include_gqs) const
95  {
96  // other kinds of Stan models need more here, but we don't.
97  std::vector<std::string> names;
98  get_param_names(names);
99  param_names.insert(param_names.end(), names.begin(), names.end());
100  }
101 
102  //----------------------------------------------------------------------
104  {
105  // default to PMNSOpt
106  if (!seed || dynamic_cast<osc::OscCalcPMNSOpt*>(seed))
107  fCalc = std::make_unique<osc::OscCalcPMNSOptStan>();
108  else if (dynamic_cast<osc::OscCalcPMNS*>(seed))
109  fCalc = std::make_unique<osc::OscCalcPMNSStan>();
110  else if (dynamic_cast<osc::OscCalcDMP*>(seed))
111  fCalc = std::make_unique<osc::OscCalcDMPStan>();
112  else if (dynamic_cast<osc::OscCalcAnalytic*>(seed))
113  fCalc = std::make_unique<osc::OscCalcAnalyticStan>();
114  else
115  {
116  std::cerr << "Unexpected oscillation calculator type: " << DemangledTypeName(seed) << std::endl;
117  std::cerr << "No conversion to Stan-aware calculator available!" << std::endl;
118  std::cerr << "Abort." << std::endl;
119  abort();
120  }
121 
122  if (seed)
123  CopyParams(seed, fCalc.get());
124 
125  }
126 
127  //----------------------------------------------------------------------
128  std::unique_ptr<IFitter::IFitSummary>
130  SystShifts &bestSysts,
131  const SeedList& seedPts,
132  const std::vector<SystShifts> &systSeedPts,
133  Verbosity verb) const
134  {
135  if (seed)
136  {
137  if (!fOscCalcCache)
138  fOscCalcCache = std::make_unique<osc::OscCalc>();
139  *fOscCalcCache = *seed;
140  }
141 
142  // check that state of warmup is what we expected
143  if (fStanConfig.num_warmup > 0 && fMCMCWarmup.NumSamples() > 0)
144  {
145  std::cerr << "You supplied a previous collection of MCMC samples for warmup and also requested num_warmup > 0 in your StanConfig." << std::endl;
146  std::cerr << "Which do you want?" << std::endl;
147  abort();
148  }
149 
150  // const-casts here because we need to initialize the writer interface, which requires a non-const pointer,
151  // but this method is const. prefer not to make the members mutable for this one instance
152  fValueWriter = std::make_unique<MemoryTupleWriter>(fStanConfig.num_samples > 0 ? const_cast<MCMCSamples*>(&fMCMCSamples) : nullptr,
153  fStanConfig.num_warmup > 0 ? const_cast<MCMCSamples*>(&fMCMCWarmup) : nullptr);
154  return IFitter::Fit(seed, bestSysts, seedPts, systSeedPts, verb);
155  }
156 
157 
158  //----------------------------------------------------------------------
159  std::unique_ptr<IFitter::IFitSummary>
161  SystShifts &systSeed,
162  Verbosity verb) const
163  {
164  CreateCalculator(seed);
165 
166  fShifts = systSeed.Copy();
167 
168  // status and other stuff that get passed back & forth between us and Stan
169  stan::callbacks::writer init_writer;
170  samplecounter_callback interrupt(std::size_t(fStanConfig.num_warmup),
171  std::size_t(fStanConfig.num_samples)); // creates a nice CAFAna-style Progress bar
172 
173  std::ostream nullStream(nullptr);
174  std::ostream & diagStream = (fStanConfig.verbosity < StanConfig::Verbosity::kQuiet) ? std::cout : nullStream;
175  std::unique_ptr<stan::callbacks::stream_logger> logger;
176  switch(fStanConfig.verbosity)
177  {
178  // if ultraverbose mode, run all the printouts to the screen...
180  logger = std::make_unique<stan::callbacks::stream_logger>(std::cout, std::cout, std::cout,
182  break;
183 
184  // and work our way down from there
186  logger = std::make_unique<stan::callbacks::stream_logger>(nullStream, std::cout, std::cout,
188  break;
189 
191  logger = std::make_unique<stan::callbacks::stream_logger>(nullStream, nullStream, std::cout,
193  break;
194 
195  // always keep errors & fatal?
197  logger = std::make_unique<stan::callbacks::stream_logger>(nullStream, nullStream, nullStream,
199  break;
200  }
201 
202 
203  // Actually run Stan!
204  // this is the init vals of all the parameters.
205  // (the pointer is because array_var_context does not have a copy constructor,
206  // so it needs to only be initialized once.)
207  std::unique_ptr<stan::io::array_var_context> init_context;
208  if (fMCMCWarmup.NumSamples() < 1)
209  init_context = std::make_unique<stan::io::array_var_context>(BuildInitContext(seed, systSeed));
210  // however if we're reusing MCMC samples from a previous run we take the last point from them
211  if (fMCMCWarmup.NumSamples() > 0)
212  {
213  std::unique_ptr<osc::IOscCalcAdjustable> calc(seed->Copy());
214  for (const auto & v : fMCMCWarmup.Vars())
215  v->SetValue(calc.get(), fMCMCWarmup.SampleValue(v, fMCMCWarmup.NumSamples()-1));
216  auto shifts = systSeed.Copy();
217  for (const auto & s : fMCMCWarmup.Systs())
218  shifts->SetShift(s, fMCMCWarmup.SampleValue(s, fMCMCWarmup.NumSamples()-1));
219  init_context = std::make_unique<stan::io::array_var_context>(BuildInitContext(calc.get(), *shifts));
220  }
221 
222  // id: only needed when running multiple chains to combine.
223  // if the grid var $PROCESS is defined, use that
224  unsigned int procId = 0;
225  const char* process = getenv("PROCESS");
226  if(process)
227  procId = std::stoul(process);
228 
229  // n.b. there are _lots_ more options for ways to call Stan but let's start here.
230  // this is an exploration using the "no-u-turn sampler" (NUTS)
231  // within the Hamiltonian MC algorithm with a simple Euclidian metric
232  // in unbounded parameter space.
233  // this call can be replaced by a call to the stock Stan function if needed (see following),
234  // but we've customized it in order to be able to save the state after warmup and restore it.
235  int return_code;
237  return_code = RunHMC<stan_dense_t>(init_writer, interrupt, diagStream, logger, *init_context, procId);
238  else
239  return_code = RunHMC<stan_diag_t>(init_writer, interrupt, diagStream, logger, *init_context, procId);
240 
241  // Use Stan's wrapper function instead of the above. Can be used for diagnostics as needed.
242 // auto return_code = stan::services::sample::hmc_nuts_diag_e_adapt(*this,
243 // init_context,
244 // fStanConfig.random_seed,
245 // procId,
246 // fStanConfig.init_radius,
247 // fStanConfig.num_warmup,
248 // fStanConfig.num_samples,
249 // fStanConfig.num_thin,
250 // fStanConfig.save_warmup,
251 // fStanConfig.refresh,
252 // fStanConfig.stepsize,
253 // fStanConfig.stepsize_jitter,
254 // fStanConfig.max_depth,
255 // fStanConfig.delta,
256 // fStanConfig.gamma,
257 // fStanConfig.kappa,
258 // fStanConfig.t0,
259 // fStanConfig.init_buffer,
260 // fStanConfig.term_buffer,
261 // fStanConfig.window,
262 // interrupt,
263 // *logger,
264 // init_writer,
265 // *fValueWriter,
266 // diagnostic_writer);
267 
268  // todo: something smarter here? or just more output?
269  // also todo: need to check the Stan diagnostics for divergences, autocorrelation, etc.
270  if (return_code != stan::services::error_codes::OK)
271  std::cerr << "warning: Stan fit did not converge..." << std::endl;
272 
273  auto bestSampleIdx = fMCMCSamples.BestFitSampleIdx();
274 
275  // Store results back to the "seed" variable
276  for (auto & var : fVars)
277  var->SetValue(seed, fMCMCSamples.SampleValue(var, bestSampleIdx));
278 
279  // Store systematic results back into "systSeed".
280  // cast to stan-var so that we don't lose the value here
281  for (const auto & syst : fSysts)
282  systSeed.SetShift(syst, stan::math::var(fMCMCSamples.SampleValue(syst, bestSampleIdx)));
283 
285 
286  return std::make_unique<StanFitSummary>(fMCMCSamples.SampleLL(bestSampleIdx));
287  } // StanFitter::FitHelperSeeded()
288 
289  //----------------------------------------------------------------------
290  void StanFitter::get_param_names(std::vector<std::string>& names) const
291  {
292  names.resize(0);
293  for (const auto & var : fVars)
294  names.push_back(var->ShortName());
295  for (const auto & syst : fSysts)
296  names.push_back(syst->ShortName());
297  }
298 
299  //----------------------------------------------------------------------
300  template <bool propto__, bool jacobian__, typename T>
301  T StanFitter::log_prob(std::vector<T>& params_real,
302  std::vector<int>& params_int, // we don't have any integer parameters, but req'd by interface
303  std::ostream*) const
304  {
305  T logprob = 0;
306 
307  try
308  {
309  // read the model parameters.
310  // they'll come out in the order they're fed in here.
311  // we were careful to feed them in in the same order
312  // as in the loops below, so hopefully no wires
313  // get crossed anywhere...
314  // (boy do I wish there were another way to do this...)
315  stan::io::reader<T> reader(params_real, params_int);
316 
318  std::cout << "-----------" << std::endl;
319 
320  std::size_t numVals = fVars.size() + fSysts.size();
321  std::vector<T> params;
322  params.reserve(numVals);
323  T val;
324  for (const auto & var : fVars)
325  {
326  // note: fVars always contains stan::math::vars...
327  auto constrVar = dynamic_cast<const IConstrainedFitVar*>(var);
328  if (constrVar)
329  {
330  if (jacobian__)
331  // GetValAs() because the version of this function where T is double gets instantiated and is req'd by Stan
332  val = reader.scalar_lub_constrain(util::GetValAs<T>(constrVar->LowLimit()),
333  util::GetValAs<T>(constrVar->HighLimit()), logprob);
334  else
335  val = reader.scalar_lub_constrain(util::GetValAs<T>(constrVar->LowLimit()),
336  util::GetValAs<T>(constrVar->HighLimit()));
337  }
338  else
339  {
340  if (jacobian__)
341  val = reader.scalar_constrain(logprob);
342  else
343  val = reader.scalar_constrain();
344  }
346  std::cout << "var " << var->ShortName() << " has val = " << util::GetValAs<double>(val) << std::endl;
347  params.push_back(val);
348  }
349  // systs are always unbounded
350  for (std::size_t idx = 0; idx < fSysts.size(); idx++)
351  {
352  if (jacobian__)
353  val = reader.scalar_constrain(logprob);
354  else
355  val = reader.scalar_constrain();
356 
358  std::cout << "syst " << fSysts[idx]->ShortName() << " has val = " << util::GetValAs<double>(val) << std::endl;
359  params.push_back(val);
360  }
361 
363  std::cout << "log-prob from re-entering constrained space = " << logprob << std::endl;
364 
365  // reset all the parameters.
366  // the fitted ones will be set according to the fitted values by DecodePars().
367  if (fOscCalcCache)
368  CopyParams(fOscCalcCache.get(), fCalc.get());
369 
370  // applies the parameters to the vars & shifts.
371  // again assumes the ordering is right... (shudder)
372  DecodePars(params, fCalc.get());
373 
374  for (unsigned int i = 0; i < fVars.size(); ++i)
375  {
376  // sadly StanFitSupport<IConstrainedFitVar> is not a derived class StanFitSupport<IFitVar>
377  auto var = dynamic_cast<const StanFitSupport<IFitVar>*>(fVars[i]);
378  auto constrVar = dynamic_cast<const StanFitSupport<IConstrainedFitVar>*>(fVars[i]);
379  assert(var || constrVar);
380  auto prior = constrVar ? constrVar->LogPrior(params[i], fCalc.get())
381  : var->LogPrior(params[i], fCalc.get());
383  {
384  std::cout << " var '" << fVars[i]->ShortName() << "' has value " << params[i]
385  << " (in calc = " << (constrVar ? constrVar->GetValue(fCalc.get()) : var->GetValue(fCalc.get())) << ")"
386  << " and ll from prior = " << prior << std::endl;
387  }
388 
389  logprob += util::GetValAs<T>(prior);
390  }
391  logprob += util::GetValAs<T>(fShifts->LogPrior());
393  std::cout << "log-prob from syst prior " << util::GetValAs<double>(fShifts->LogPrior()) << std::endl;
394 
395  auto ll = fExpt->LogLikelihood(fCalc.get(), *fShifts);
397  std::cout << "log-prob from spectrum comparison = " << ll << std::endl;
398  logprob += util::GetValAs<T>(ll);
399 
400  }
401  catch (const std::exception& e)
402  {
403  // just catch to add a bit more info to make it easier to find
404  std::cerr << "Exception during evaluation of likelihood in Stan fitting (see StanFitter::log_prob()): " << std::endl;
405  std::cerr << e.what();
406  throw;
407  }
408 
410  std::cout << "accumulated prob = " << logprob << std::endl;
411  // Stan will "recover memory" (i.e., clear the var cache)
412  // after we're done here.
413  // Invalidate the cache so that it can't be re-used.
414  fCalc->InvalidateCache();
415  return logprob;
416 
417  } // StanFitter::log_prob()
418 
419  //----------------------------------------------------------------------
420  template <bool propto, bool jacobian, typename T_>
421  T_ StanFitter::log_prob(Eigen::Matrix<T_,Eigen::Dynamic,1>& params_r,
422  std::ostream* pstream) const
423  {
424  std::vector<T_> vec_params_r;
425  vec_params_r.reserve(params_r.size());
426  for (int i = 0; i < params_r.size(); ++i)
427  vec_params_r.push_back(params_r(i));
428  std::vector<int> vec_params_i;
429  return log_prob<propto,jacobian,T_>(vec_params_r, vec_params_i, pstream);
430  }
431 
432  //----------------------------------------------------------------------
434  {
435  // before going ahead, check that the supplied samples are compatible
436  // with how this Fitter was initialized.
437  if (fMCMCWarmup.Vars() != warmup.Vars())
438  {
439  std::cerr << "Warmup to be reused is incompatible with this StanFitter." << std::endl;
440  std::cerr << " StanFitter fit vars: ";
441  for (const auto & v : fMCMCWarmup.Vars())
442  std::cerr << v << " ";
443  std::cerr << std::endl;
444  std::cerr << " Supplied warmup fit vars: ";
445  for (const auto & v : warmup.Vars())
446  std::cerr << v << " ";
447  std::cerr << std::endl;
448 
449  abort();
450  }
451  if (fMCMCWarmup.Systs() != warmup.Systs())
452  {
453  std::cerr << "Warmup to be reused is incompatible with this StanFitter." << std::endl;
454  std::cerr << " StanFitter systs: ";
455  for (const auto & s : fMCMCWarmup.Systs())
456  std::cerr << s << " ";
457  std::cerr << std::endl;
458  std::cerr << " Supplied systs: ";
459  for (const auto & s : warmup.Systs())
460  std::cerr << s << " ";
461  std::cerr << std::endl;
462 
463  abort();
464  }
465 
466  fMCMCWarmup = std::move(warmup);
467  }
468 
469  //----------------------------------------------------------------------
470  template <typename Sampler>
473  std::ostream &diagStream,
474  const std::unique_ptr<stan::callbacks::stream_logger> &logger,
475  stan::io::array_var_context &init_context,
476  unsigned int procId) const
477  {
478  // for now, just write the diagnostics to stdout.
479  // we can do something fancier later...
480  stan::callbacks::stream_writer diagnostic_writer(diagStream, "# ");
481 
482  // *********************************************
483  // originally cribbed from stan::services::sample::hmc_nuts_diag_e_adapt()
484  // (in stan/services/sample/hmc_nuts_diag_e_adapt.hpp)
485  // but 'metric' and 'step size' code rewritten to support reloading from previous warmup
487 
488 
489 
490 
491  boost::ecuyer1988 rng = stan::services::util::create_rng(fStanConfig.random_seed, procId);
492 
493  // continuous parameters
494  std::vector<double> cont_vector = stan::services::util::initialize(*this,
495  init_context,
496  rng,
498  true,
499  *logger,
500  init_writer);
501 
502  Eigen::VectorXd inv_metric;
504  {
505  try
506  {
508  stan::io::var_context &unit_e_metric = dmp;
509  inv_metric = stan::services::util::read_diag_inv_metric(unit_e_metric, num_params_r(), *logger);
511  }
512  catch (const std::domain_error &e)
513  {
515  }
516  }
517  else
518  {
520  }
521 
523  {
524  Sampler sampler(*this, rng);
525 
526  sampler.set_metric(inv_metric);
528  sampler.set_nominal_stepsize(fStanConfig.stepsize);
529  else
530  sampler.set_nominal_stepsize(fMCMCWarmup.Hyperparams().stepSize);
531  sampler.set_stepsize_jitter(fStanConfig.stepsize_jitter);
532  sampler.set_max_depth(fStanConfig.max_depth);
533 
534  sampler.get_stepsize_adaptation().set_mu(log(10 * fStanConfig.stepsize));
535  sampler.get_stepsize_adaptation().set_delta(fStanConfig.delta);
536  sampler.get_stepsize_adaptation().set_gamma(fStanConfig.gamma);
537  sampler.get_stepsize_adaptation().set_kappa(fStanConfig.kappa);
538  sampler.get_stepsize_adaptation().set_t0(fStanConfig.t0);
539 
540  sampler.set_window_params(fStanConfig.num_warmup,
544  *logger);
545 
546  RunSampler(sampler, cont_vector, rng, interrupt, *logger, diagnostic_writer);
547 
548 // stan::services::util::run_adaptive_sampler(sampler,
549 // *this,
550 // cont_vector,
551 // fStanConfig.num_warmup,
552 // fStanConfig.num_samples,
553 // fStanConfig.num_thin,
554 // fStanConfig.refresh,
555 // fStanConfig.save_warmup,
556 // rng,
557 // interrupt,
558 // *logger,
559 // *fValueWriter,
560 // diagnostic_writer);
561 
562  }
563 
564  // *********************************************
565 
566  return return_code;
567  }
568 
569  // see at the top of this file for the type aliases
570  template int StanFitter::RunHMC<stan_diag_t>(stan::callbacks::writer &init_writer,
572  std::ostream &diagStream,
573  const std::unique_ptr<stan::callbacks::stream_logger> &logger,
574  stan::io::array_var_context &init_context,
575  unsigned int procId) const;
576  template int StanFitter::RunHMC<stan_dense_t>(stan::callbacks::writer &init_writer,
578  std::ostream &diagStream,
579  const std::unique_ptr<stan::callbacks::stream_logger> &logger,
580  stan::io::array_var_context &init_context,
581  unsigned int procId) const;
582 
583 
584  // StanFitter::RunHMC()
585 
586  //----------------------------------------------------------------------
587  template <typename Sampler>
588  void StanFitter::RunSampler(Sampler& sampler,
589  std::vector<double>& cont_vector,
590  boost::ecuyer1988& rng,
591  stan::callbacks::interrupt& interrupt,
592  stan::callbacks::logger& logger,
593  stan::callbacks::writer& diagnostic_writer) const
594  {
595  // Stan passes around a templated Model type, but this class is the Model.
596  auto & model = *this;
597 
598  // *********************************************
599  // Cribbed from stan::services::run_adaptive_sampler() (in stan/services/util/run_adaptive_sampler.hpp)
600  Eigen::Map<Eigen::VectorXd> cont_params(cont_vector.data(),
601  cont_vector.size());
602 
603  stan::services::util::mcmc_writer writer(*fValueWriter, diagnostic_writer, logger);
604  stan::mcmc::sample s(cont_params, 0, 0);
605 
606  clock_t start, end;
607  double warm_delta_t = 0;
608 
609  if (fMCMCWarmup.NumSamples() < 1)
610  {
611  sampler.engage_adaptation();
612  try
613  {
614  sampler.z().q = cont_params;
615  sampler.init_stepsize(logger);
616  } catch (const std::exception &e)
617  {
618  logger.info("Exception initializing step size.");
619  logger.info(e.what());
620  return;
621  }
622 
623  // -------------------------------
624  // this is a new bit we added in between the bits copied from Stan
626  // -------------------------------
627 
628 
629  // Headers
630  writer.write_sample_names(s, sampler, model);
631  writer.write_diagnostic_names(s, sampler, model);
632 
633  start = clock();
636  0,
641  true,
642  writer,
643  s,
644  model,
645  rng,
646  interrupt,
647  logger);
648  end = clock();
649  warm_delta_t = static_cast<double>(end - start) / CLOCKS_PER_SEC;
650 
651  sampler.disengage_adaptation();
652  writer.write_adapt_finish(sampler);
653 
654 
655  // -------------------------------
656  // this is a new bit we added in between the bits copied from Stan
657  fValueWriter->SaveSamplerState(sampler, warm_delta_t);
658  }
659  if (fStanConfig.num_samples < 1)
660  return;
661 
663  writer.write_sample_names(s, sampler, model);
664  writer.write_diagnostic_names(s, sampler, model);
665  // -------------------------------
666 
667  start = clock();
674  true,
675  false,
676  writer,
677  s,
678  model,
679  rng,
680  interrupt,
681  logger);
682  end = clock();
683  double sample_delta_t = static_cast<double>(end - start) / CLOCKS_PER_SEC;
684 
685  writer.write_timing(warm_delta_t, sample_delta_t);
686 
687  // *********************************************
688 
689  fValueWriter->SaveSamplerState(sampler, sample_delta_t);
690 
691  }
692 
693  // see above near the top of the file for the type aliases
694  template void StanFitter::RunSampler<stan_diag_t>(stan_diag_t& sampler,
695  std::vector<double>& cont_vector,
696  boost::ecuyer1988& rng,
697  stan::callbacks::interrupt& interrupt,
698  stan::callbacks::logger& logger,
699  stan::callbacks::writer& diagnostic_writer) const;
700  template void StanFitter::RunSampler<stan_dense_t>(stan_dense_t& sampler,
701  std::vector<double>& cont_vector,
702  boost::ecuyer1988& rng,
703  stan::callbacks::interrupt& interrupt,
704  stan::callbacks::logger& logger,
705  stan::callbacks::writer& diagnostic_writer) const;
706 
707  //----------------------------------------------------------------------
709  SystShifts &systSeed) const
710  {
711  if (seed)
712  {
713  if (!fOscCalcCache)
714  fOscCalcCache = std::make_unique<osc::OscCalc>();
715  *fOscCalcCache = *seed;
716  }
717 
718  CreateCalculator(seed);
719  fShifts = systSeed.Copy();
720 
721  // status and other stuff that get passed back & forth between us and Stan
722  stan::callbacks::writer init_writer;
723 
724  stan::callbacks::interrupt interrupt;
725 
726  // always go to screen since this is always a test
728  stan::callbacks::stream_writer diagnostic_writer(std::cout, "# ");
729 
730  // this is the init vals of all the parameters
731  auto init_context = BuildInitContext(seed, systSeed);
732 
733  // this would normally get initialized in Fit(), but we're not fitting
734  fValueWriter = std::make_unique<MemoryTupleWriter>(fStanConfig.num_samples > 0 ? const_cast<MCMCSamples*>(&fMCMCSamples) : nullptr,
735  fStanConfig.num_warmup > 0 ? const_cast<MCMCSamples*>(&fMCMCWarmup) : nullptr);
736  // diagnostic mode, where the model's gradients calculated via Stan's autodiff
737  // are compared to those from finite difference calculations
738  const double error = 1e-6;
739  std::cout << "fValueWriter is at: " << fValueWriter.get() << std::endl;
741  init_context,
743  0,
745  1e-6, // epsilon (finite diff step size) -- see CmdStan manual
746  error,
747  interrupt,
748  logger,
749  init_writer,
750  *fValueWriter);
751  if (return_code != 0)
752  std::cerr << "warning: Stan diagnostic gradient test reports " << return_code
753  << " parameters were not within " << error << " of finite diff expectation..." << std::endl;
754  }
755 
756  //----------------------------------------------------------------------
757  template <typename T>
759  stan::io::writer<double>& writer,
760  const T& var) const
761  {
763  "transform_helper() can only be instantiated with IFitVar* or ISyst* as the type");
764 
765  if (!(context.contains_r(var->ShortName())))
766  throw std::runtime_error("variable '" + var->ShortName() + "' missing from context given to Stan!");
767 
768  std::vector<double> vals_real = context.vals_r(var->ShortName());
769  double val = vals_real[0]; // we always have single-valued variables
770  try
771  {
772  // note: T is either IFitVar* or ISyst*, and systs are always unbounded
773  const IConstrainedFitVar * constrFitVar;
774  if ( (constrFitVar = dynamic_cast<const IConstrainedFitVar*>(var)) )
775  {
777  {
778  std::cout << " transforming variable " << constrFitVar->ShortName()
779  << " from constrained between " << constrFitVar->LowLimit() << " and " << constrFitVar->HighLimit()
780  << " to unbounded space" << std::endl;
781  }
782  writer.scalar_lub_unconstrain(util::GetValAs<double>(constrFitVar->LowLimit()),
783  util::GetValAs<double>(constrFitVar->HighLimit()), val);
784  }
785  else
786  writer.scalar_unconstrain(val);
787 
788  }
789  catch (const std::exception& e)
790  {
791  throw std::runtime_error(std::string("Error transforming variable '" + var->ShortName() + "': ") + e.what());
792  }
793 
794  }
795 
796  // explicitly instantiate for the two relevant cases
799  const IFitVar * const &) const;
802  const ISyst * const &) const;
803 
804  //----------------------------------------------------------------------
806  std::vector<int>& params_int,
807  std::vector<double>& params_real,
808  std::ostream*) const
809  {
810 
811  stan::io::writer<double> stanWriter(params_real, params_int);
812 
813  for (const auto & var : fVars)
814  transform_helper(context, stanWriter, var);
815  for (const auto & syst : fSysts)
816  transform_helper(context, stanWriter, syst);
817 
818  params_real = stanWriter.data_r();
819  params_int = stanWriter.data_i();
820  }
821 
822 
823  //----------------------------------------------------------------------
824  template <typename RNG>
826  std::vector<double>& params_real,
827  std::vector<int>& params_int,
828  std::vector<double>& vars,
829  bool, // these are used in autogenerated models,
830  bool, // but aren't needed in our version
831  std::ostream*) const
832  {
833  vars.resize(0);
834  stan::io::reader<double> paramReader(params_real, params_int);
835  for (std::size_t i = 0; i < params_real.size(); i++)
836  {
837  const IConstrainedFitVar * constrFitVar = nullptr;
838  if (i < fVars.size())
839  constrFitVar = dynamic_cast<const IConstrainedFitVar*>(fVars[i]);
840  double parVal;
841  if (constrFitVar)
842  {
844  {
845  std::cout << " transforming variable " << constrFitVar->ShortName()
846  << " from unbounded space"
847  << " to constrained between " << constrFitVar->LowLimit() << " and " << constrFitVar->HighLimit()
848  << std::endl;
849  }
850  parVal = paramReader.scalar_lub_constrain(util::GetValAs<double>(constrFitVar->LowLimit()),
851  util::GetValAs<double>(constrFitVar->HighLimit()));
853  std::cout << " got value = " << parVal << std::endl;
854  }
855  else
856  parVal = paramReader.scalar_constrain();
857  vars.push_back(parVal);
858  }
859  }
860 
861 
862 } // namespace ana
double delta
Adaptation target acceptance statistic; 0 < delta < 1.
Definition: StanConfig.h:80
int diagnose(Model &model, stan::io::var_context &init, unsigned int random_seed, unsigned int chain, double init_radius, double epsilon, double error, callbacks::interrupt &interrupt, callbacks::logger &logger, callbacks::writer &init_writer, callbacks::writer &parameter_writer)
Definition: diagnose.hpp:43
std::unique_ptr< TMatrixD > invMetric
Definition: MCMCSamples.h:41
virtual double HighLimit() const =0
T scalar_lub_constrain(const TL lb, const TU ub)
Definition: reader.hpp:657
std::size_t NumSamples() const
How many samples do we have?
Definition: MCMCSamples.h:135
virtual stan::math::var GetValue(const osc::IOscCalcAdjustableStan *osc) const =0
virtual _IOscCalcAdjustable< T > * Copy() const =0
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
virtual std::unique_ptr< SystShifts > Copy() const
Definition: SystShifts.cxx:67
void get_param_names(std::vector< std::string > &names) const
Return names of parameters. (Required by Stan interface.)
Definition: StanFitter.cxx:290
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
T scalar_constrain()
Definition: reader.hpp:165
int RunHMC(stan::callbacks::writer &init_writer, StanFitter::samplecounter_callback &interrupt, std::ostream &diagStream, const std::unique_ptr< stan::callbacks::stream_logger > &logger, stan::io::array_var_context &init_context, unsigned int procId) const
Definition: StanFitter.cxx:471
std::string DemangledTypeName(T *obj)
utility method to figure out exactly what kind of object I am
Definition: UtilsExt.h:114
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
double stepsize
Integrator step size used in each simulation. Typically this is adaptive and the default value is a f...
Definition: StanConfig.h:68
std::vector< T > & data_r()
Definition: writer.hpp:70
virtual double LowLimit() const =0
OStream cerr
Definition: OStream.cxx:7
double SampleValue(const IFitVar *var, std::size_t idx) const
Get the value of Var var for sample number .
Definition: MCMCSamples.h:169
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
MCMCSamples fMCMCWarmup
Definition: StanFitter.h:415
const std::vector< const IFitVar * > & Vars() const
Which Vars are sampled in these samples?
Definition: MCMCSamples.h:212
void write_array(RNG &base_rng__, std::vector< double > &params_real, std::vector< int > &params_int, std::vector< double > &vars, bool include_tparams__=true, bool include_gqs__=true, std::ostream *pstream__=0) const
Method for writing out data inherited from Stan interface.
Definition: StanFitter.cxx:825
Verbosity verbosity
How verbose do you want me to be?
Definition: StanConfig.h:89
void DecodePars(const std::vector< T > &pars, osc::IOscCalcAdjustableStan *calc) const
Definition: StanFitter.h:303
void scalar_lub_unconstrain(double lb, double ub, T &y)
Definition: writer.hpp:163
bool save_warmup
Save the warmup steps?
Definition: StanConfig.h:66
stan::io::array_var_context BuildInitContext(osc::IOscCalcAdjustable *seed, const SystShifts &systSeed) const
Helper function to build the initial seed point context for Stan.
Definition: StanFitter.cxx:57
osc::OscCalcDumb calc
virtual std::vector< double > vals_r(const std::string &name) const =0
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
virtual bool contains_r(const std::string &name) const =0
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
MCMCSamples fMCMCSamples
Definition: StanFitter.h:414
void scalar_unconstrain(T &y)
Definition: writer.hpp:101
T__ log_prob(std::vector< T__ > &params_real, std::vector< int > &params_i__, std::ostream *pstream__=0) const
std::vector< const ISyst * > fSysts
Definition: IFitter.h:115
unsigned int term_buffer
Width of final fast adaptation interval.
Definition: StanConfig.h:85
const XML_Char * s
Definition: expat.h:262
double init_radius
Size of the range in unconstrained parameter space where the initial point for un-specified parameter...
Definition: StanConfig.h:62
StanFitter(const IExperiment *expt, std::vector< const IFitVar * > vars, std::vector< const ISyst * > systs={})
Definition: StanFitter.cxx:44
expt
Definition: demo5.py:34
stan::io::dump create_unit_e_diag_inv_metric(size_t num_params)
void RunDiagnostics(const StanConfig &cfg) const
Do some checks on the post-fit samples.
void transform_inits(const stan::io::var_context &context, std::vector< int > &params_int, std::vector< double > &params_real, std::ostream *pstream__) const
Take the seed points and convert them into Stan&#39;s unconstrained internal space.
Definition: StanFitter.cxx:805
const XML_Char int const XML_Char * value
Definition: expat.h:331
virtual std::unique_ptr< IFitSummary > Fit(osc::IOscCalcAdjustable *seed, SystShifts &bestSysts=junkShifts, const SeedList &seedPts=SeedList(), const std::vector< SystShifts > &systSeedPts={}, Verbosity verb=kVerbose) const override
Definition: StanFitter.cxx:129
unsigned int seed
Definition: runWimpSim.h:102
int refresh
Number of iterations between printout updates. Only relevant when fitting in &#39;verbose&#39; mode...
Definition: StanConfig.h:67
std::string getenv(std::string const &name)
void CopyParams(const osc::_IOscCalcAdjustable< T > *inCalc, osc::_IOscCalcAdjustable< U > *outCalc)
Copy parameters from one calculator to another, irrespective of their type.
Definition: IOscCalc.cxx:62
virtual std::unique_ptr< IFitSummary > Fit(osc::IOscCalcAdjustable *seed, SystShifts &bestSysts=junkShifts, const SeedList &seedPts=SeedList(), const std::vector< SystShifts > &systSeedPts={}, Verbosity verb=kVerbose) const
Master fitting method. Depends on FitHelper and FitHelperSeeded.
Definition: IFitter.cxx:69
T GetShift(const ISyst *syst) const
const std::map< std::pair< std::string, std::string >, Variable > vars
Base class for fitters.
Definition: IFitter.h:16
int num_warmup
Number of initial steps in the Markov chain (used to enter the typical set). These are usually discar...
Definition: StanConfig.h:63
void constrained_param_names(std::vector< std::string > &param_names, bool include_tparams=true, bool include_gqs=true) const
Return names of parameters . (Required by Stan interface.)
Definition: StanFitter.cxx:93
const IExperiment * fExpt
Definition: StanFitter.h:412
void validate_diag_inv_metric(const Eigen::VectorXd &inv_metric, callbacks::logger &logger)
void ReuseWarmup(MCMCSamples &&warmup)
Definition: StanFitter.cxx:433
const Hyperparameters & Hyperparams() const
Definition: MCMCSamples.h:90
const std::vector< const ana::ISyst * > & Systs() const
Which Systs are sampled in these samples?
Definition: MCMCSamples.h:209
void RunSampler(Sampler &sampler, std::vector< double > &cont_vector, boost::ecuyer1988 &rng, stan::callbacks::interrupt &interrupt, stan::callbacks::logger &logger, stan::callbacks::writer &diagnostic_writer) const
Definition: StanFitter.cxx:588
std::vector< double > initialize(Model &model, stan::io::var_context &init, RNG &rng, double init_radius, bool print_timing, stan::callbacks::logger &logger, stan::callbacks::writer &init_writer)
Definition: initialize.hpp:68
bool denseMassMx
Should the mass matrix used in HMC be diagonal (default) or dense?
Definition: StanConfig.h:90
void transform_helper(const stan::io::var_context &context, stan::io::writer< double > &writer, const T &var) const
Helper function that actually does the unconstraining Stan needs.
Definition: StanFitter.cxx:758
OStream cout
Definition: OStream.cxx:6
Eigen::MatrixXd EigenMatrixXdFromTMatrixD(const TMatrixD *mat)
Definition: EigenUtils.cxx:10
void domain_error(const char *function, const char *name, const T &y, const char *msg1, const char *msg2)
unsigned int window
Definition: StanConfig.h:86
virtual stan::math::var LogPrior(const stan::math::var &var, const osc::IOscCalcAdjustableStan *calc) const
Definition: IFitVar.h:102
double SampleLL(std::size_t idx) const
Get the LL for sample number .
Definition: MCMCSamples.h:162
virtual stan::math::var LogLikelihood(osc::IOscCalcAdjustableStan *osc, const SystShifts &syst=kNoShift) const
Definition: IExperiment.h:25
const XML_Char * context
Definition: expat.h:434
virtual void info(const std::string &message)
Definition: logger.hpp:47
StanConfig fStanConfig
Configuration passed to Stan for fitting. See the StanConfig struct documentation for ideas...
Definition: StanFitter.h:413
std::unique_ptr< osc::IOscCalcAdjustableStan > fCalc
Definition: StanFitter.h:411
double stepsize_jitter
Randomly "jitter" the step size by this fraction of the step size. (Sometimes can help if the integra...
Definition: StanConfig.h:69
std::vector< const IFitVar * > fVars
Definition: IFitter.h:114
int max_depth
Depth of the binary tree used by the HMC when taking its leapfrog steps. Note: the number of leapfrog...
Definition: StanConfig.h:70
const std::string & ShortName() const
Definition: IFitVar.h:36
double gamma
Adaptation regularization scale; gamma > 0.
Definition: StanConfig.h:82
std::unique_ptr< IFitSummary > FitHelperSeeded(osc::IOscCalcAdjustable *seed, SystShifts &systSeed, Verbosity verb) const override
Implement workhorse method from IFitter interface.
Definition: StanFitter.cxx:160
void CreateCalculator(osc::IOscCalcAdjustable *seed) const
Convert a &#39;normal&#39; calculator into the Stan-aware variant used internally.
Definition: StanFitter.cxx:103
unsigned int init_buffer
Width of initial fast adaptation interval.
Definition: StanConfig.h:84
int num_thin
Number of Markov chain steps between saved samples when sampling.
Definition: StanConfig.h:65
Eigen::VectorXd read_diag_inv_metric(stan::io::var_context &init_context, size_t num_params, callbacks::logger &logger)
Base class defining interface for experiments.
Definition: IExperiment.h:14
std::unique_ptr< MemoryTupleWriter > fValueWriter
Definition: StanFitter.h:419
void TestGradients(osc::IOscCalcAdjustable *seed, SystShifts &systSeed) const
Run Stan&#39;s test of its auto-differentiation (comparing to a finite-differences calculation) ...
Definition: StanFitter.cxx:708
std::unique_ptr< osc::IOscCalcAdjustable > fOscCalcCache
Definition: StanFitter.h:432
assert(nhit_max >=nhit_nbins)
Interface definition for fittable variables.
Definition: IFitVar.h:16
double kappa
Adaptation relaxation exponent; kappa > 0.
Definition: StanConfig.h:81
std::unique_ptr< SystShifts > fShifts
Definition: IFitter.h:116
size_t num_params_r() const
Definition: prob_grad.hpp:36
std::vector< int > & data_i()
Definition: writer.hpp:81
double T
Definition: Xdiff_gwt.C:5
Float_t e
Definition: plot.C:35
double t0
Adaptation iteration offset; t0 > 0.
Definition: StanConfig.h:83
const XML_Char XML_Content * model
Definition: expat.h:151
void generate_transitions(stan::mcmc::base_mcmc &sampler, int num_iterations, int start, int finish, int num_thin, int refresh, bool save, bool warmup, util::mcmc_writer &mcmc_writer, stan::mcmc::sample &init_s, Model &model, RNG &base_rng, callbacks::interrupt &callback, callbacks::logger &logger)
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:80
std::size_t BestFitSampleIdx() const
Retrieve the index of the best-fit sample (lowest LL), which can then be used with the SampleLL() or ...
unsigned int random_seed
Random seed used by Stan internally.
Definition: StanConfig.h:60
boost::ecuyer1988 create_rng(unsigned int seed, unsigned int chain)
Definition: create_rng.hpp:25
int num_samples
Number of steps in the Markov chain retained for analysis (after warmup).
Definition: StanConfig.h:64
enum BeamMode string