Classes | Public Types | Public Member Functions | Protected Member Functions | Protected Attributes | Static Protected Attributes | Private Member Functions | Private Attributes | List of all members
ana::StanFitter Class Reference

Fitter type that bolts the Stan fitting tools onto CAFAna. More...

#include "/cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N20-11-25/CAFAna/Fit/StanFitter.h"

Inheritance diagram for ana::StanFitter:
ana::IFitter stan::model::prob_grad

Classes

class  samplecounter_callback
 
class  StanFitSummary
 

Public Types

enum  Verbosity { kQuiet, kVerbose }
 

Public Member Functions

 StanFitter (const IExperiment *expt, std::vector< const IFitVar * > vars, std::vector< const ISyst * > systs={})
 
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.) More...
 
virtual std::unique_ptr< IFitSummaryFit (osc::IOscCalcAdjustable *seed, SystShifts &bestSysts=junkShifts, const SeedList &seedPts=SeedList(), const std::vector< SystShifts > &systSeedPts={}, Verbosity verb=kVerbose) const override
 
void get_dims (std::vector< std::vector< size_t > > &dims) const
 
void get_param_names (std::vector< std::string > &names) const
 Return names of parameters. (Required by Stan interface.) More...
 
template<bool propto__, bool jacobian__, typename T__ >
T__ log_prob (std::vector< T__ > &params_real, std::vector< int > &params_i__, std::ostream *pstream__=0) const
 
template<bool propto, bool jacobian, typename T_ >
T_ log_prob (Eigen::Matrix< T_, Eigen::Dynamic, 1 > &params_r, std::ostream *pstream=0) const
 Version with Matrix is req'd by Stan interface, but we won't use it. More...
 
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's unconstrained internal space. More...
 
template<typename RNG >
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. More...
 
const MCMCSamplesGetSamples (MemoryTupleWriter::WhichSamples ws=MemoryTupleWriter::WhichSamples::kPostWarmup) const
 Peruse the samples the MCMC generated. If you want to take ownership of them instead, see ExtractSamples() More...
 
void ExtractSamples (MCMCSamples &samples, MemoryTupleWriter::WhichSamples ws=MemoryTupleWriter::WhichSamples::kPostWarmup)
 Extract the MCMC samples from the fitter (you own them now). More...
 
std::unique_ptr< MCMCSamplesExtractSamples (MemoryTupleWriter::WhichSamples ws=MemoryTupleWriter::WhichSamples::kPostWarmup)
 Another way of extracting the MCMC samples (so you own them now) if you prefer unique_ptrs. More...
 
void ReuseWarmup (MCMCSamples &&warmup)
 
void SetStanConfig (const StanConfig &cfg)
 Change the config used for Stan. See the StanConfig struct documentation for ideas. More...
 
void TestGradients (osc::IOscCalcAdjustable *seed, SystShifts &systSeed) const
 Run Stan's test of its auto-differentiation (comparing to a finite-differences calculation) More...
 
void unconstrained_param_names (std::vector< std::string > &param_names, bool include_tparams=true, bool include_gqs=true) const
 Return names of parameters in Stan's unconstrained space. (Required by Stan interface.) More...
 
template<bool propto__, bool jacobian__, typename T >
T log_prob (std::vector< T > &params_real, std::vector< int > &params_int, std::ostream *) const
 
virtual std::unique_ptr< IFitSummaryFit (osc::IOscCalcAdjustable *seed, SystShifts &systSeed, Verbosity verb) const
 Variant with no seedPts. More...
 
virtual std::unique_ptr< IFitSummaryFit (osc::IOscCalcAdjustable *seed, Verbosity verb) const
 Variant with no seedPts and no systematics result returned. More...
 
virtual std::unique_ptr< IFitSummaryFit (SystShifts &systSeed, Verbosity verb=kVerbose) const
 Variant with no oscillations - useful for ND fits. More...
 
std::unique_ptr< SystShiftsGetSystShifts () const
 
size_t num_params_r () const
 
size_t num_params_i () const
 
std::pair< int, intparam_range_i (size_t idx) const
 

Protected Member Functions

std::vector< SeedPtExpandSeeds (const SeedList &seedPts, const std::vector< SystShifts > &systSeedPts) const
 
virtual std::unique_ptr< IFitSummaryFitHelper (osc::IOscCalcAdjustable *seed, SystShifts &bestSysts, const SeedList &seedPts, const std::vector< SystShifts > &systSeedPts, Verbosity verb) const
 Helper for Fit. More...
 
void ValidateSeeds (osc::IOscCalcAdjustable *seed, const SeedList &seedPts, const std::vector< SystShifts > &systSeedPts) const
 Check that the seeds that were specified are compatible with the vars being fitted. More...
 

Protected Attributes

std::vector< const IFitVar * > fVars
 
std::vector< const ISyst * > fSysts
 
std::unique_ptr< SystShiftsfShifts
 
Verbosity fVerb
 
size_t num_params_r__
 
std::vector< std::pair< int, int > > param_ranges_i__
 

Static Protected Attributes

static SystShifts junkShifts = SystShifts()
 

Private Member Functions

stan::io::array_var_context BuildInitContext (osc::IOscCalcAdjustable *seed, const SystShifts &systSeed) const
 Helper function to build the initial seed point context for Stan. More...
 
void CreateCalculator (osc::IOscCalcAdjustable *seed) const
 Convert a 'normal' calculator into the Stan-aware variant used internally. More...
 
template<typename T >
void DecodePars (const std::vector< T > &pars, osc::IOscCalcAdjustableStan *calc) const
 
std::unique_ptr< IFitSummaryFitHelperSeeded (osc::IOscCalcAdjustable *seed, SystShifts &systSeed, Verbosity verb) const override
 Implement workhorse method from IFitter interface. More...
 
template<typename Sampler >
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
 
template<typename Sampler >
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
 
template<typename T >
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. More...
 

Private Attributes

std::unique_ptr< osc::IOscCalcAdjustableStanfCalc
 
const IExperimentfExpt
 
StanConfig fStanConfig
 Configuration passed to Stan for fitting. See the StanConfig struct documentation for ideas. More...
 
MCMCSamples fMCMCSamples
 
MCMCSamples fMCMCWarmup
 
std::unique_ptr< MemoryTupleWriterfValueWriter
 
std::unique_ptr< osc::IOscCalcAdjustablefOscCalcCache
 

Detailed Description

Fitter type that bolts the Stan fitting tools onto CAFAna.

Unfortunately CAFAna and Stan have different case conventions so the capitalization of the methods here is a haphazard mix :(

Definition at line 126 of file StanFitter.h.

Member Enumeration Documentation

enum ana::IFitter::Verbosity
inherited
Enumerator
kQuiet 
kVerbose 

Definition at line 22 of file IFitter.h.

Constructor & Destructor Documentation

ana::StanFitter::StanFitter ( const IExperiment expt,
std::vector< const IFitVar * >  vars,
std::vector< const ISyst * >  systs = {} 
)

Definition at line 44 of file StanFitter.cxx.

47  : IFitter(vars, systs),
48  stan::model::prob_grad(fVars.size() + fSysts.size()),
49  fCalc(nullptr),
50  fExpt(expt),
53  {}
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
MCMCSamples fMCMCWarmup
Definition: StanFitter.h:415
MCMCSamples fMCMCSamples
Definition: StanFitter.h:414
std::vector< const ISyst * > fSysts
Definition: IFitter.h:115
expt
Definition: demo5.py:34
const std::map< std::pair< std::string, std::string >, Variable > vars
const IExperiment * fExpt
Definition: StanFitter.h:412
std::unique_ptr< osc::IOscCalcAdjustableStan > fCalc
Definition: StanFitter.h:411
std::vector< const IFitVar * > fVars
Definition: IFitter.h:114
IFitter(const std::vector< const IFitVar * > &vars, const std::vector< const ISyst * > &systs={})
Definition: IFitter.cxx:21

Member Function Documentation

stan::io::array_var_context ana::StanFitter::BuildInitContext ( osc::IOscCalcAdjustable seed,
const SystShifts systSeed 
) const
private

Helper function to build the initial seed point context for Stan.

Definition at line 57 of file StanFitter.cxx.

References ana::assert(), om::cout, allTimeWatchdog::endl, fStanConfig, ana::IFitter::fSysts, ana::IFitter::fVars, ana::SystShifts::GetShift(), ana::StanConfig::kQuiet, gen_hdf5record::names, runNovaSAM::ret, registry_explorer::v, and ana::StanConfig::verbosity.

Referenced by FitHelperSeeded(), and TestGradients().

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  }
Verbosity verbosity
How verbose do you want me to be?
Definition: StanConfig.h:89
std::vector< const ISyst * > fSysts
Definition: IFitter.h:115
const XML_Char * s
Definition: expat.h:262
OStream cout
Definition: OStream.cxx:6
StanConfig fStanConfig
Configuration passed to Stan for fitting. See the StanConfig struct documentation for ideas...
Definition: StanFitter.h:413
std::vector< const IFitVar * > fVars
Definition: IFitter.h:114
assert(nhit_max >=nhit_nbins)
void ana::StanFitter::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.)

Parameters
param_namesVector of parameter names to be filled.
include_tparamsUnused in this implementation. (Other types of Stan models need it.)
include_gqsUnused in this implementation. (Other types of Stan models need it.)

Definition at line 93 of file StanFitter.cxx.

References get_param_names(), and gen_hdf5record::names.

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  }
void get_param_names(std::vector< std::string > &names) const
Return names of parameters. (Required by Stan interface.)
Definition: StanFitter.cxx:290
void ana::StanFitter::CreateCalculator ( osc::IOscCalcAdjustable seed) const
private

Convert a 'normal' calculator into the Stan-aware variant used internally.

Definition at line 103 of file StanFitter.cxx.

References om::cerr, osc::CopyParams(), ana::DemangledTypeName(), allTimeWatchdog::endl, and fCalc.

Referenced by FitHelperSeeded(), and TestGradients().

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  }
std::string DemangledTypeName(T *obj)
utility method to figure out exactly what kind of object I am
Definition: UtilsExt.h:114
OStream cerr
Definition: OStream.cxx:7
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
std::unique_ptr< osc::IOscCalcAdjustableStan > fCalc
Definition: StanFitter.h:411
template<typename T >
void ana::StanFitter::DecodePars ( const std::vector< T > &  pars,
osc::IOscCalcAdjustableStan calc 
) const
inlineprivate

Unroll the parameters and stuff them into the calculator and/or syst shifts obj. Needs to be templated because Stan wants to instantiate it with <double>, even though that'll never actually be used (sigh)...

Definition at line 303 of file StanFitter.h.

References ana::assert(), MECModelEnuComparisons::i, calib::j, ana::StanFitSupport< VarClass >::SetValue(), febshutoff_auto::val, and PandAna.Demos.tute_pid_validation::var.

Referenced by log_prob().

304  {
306  "DecodePars() can only be used with double or stan::math::var");
307  assert(pars.size() == fVars.size()+fSysts.size());
308 
309  if (fVars.size() > 0)
310  {
311  assert(calc);
312  for(unsigned int i = 0; i < fVars.size(); ++i){
313  auto val = pars[i];
314  // sadly StanFitSupport<IConstrainedFitVar> is not a derived class StanFitSupport<IFitVar>
315  auto var = dynamic_cast<const StanFitSupport<IFitVar>*>(fVars[i]);
316  auto constrVar = dynamic_cast<const StanFitSupport<IConstrainedFitVar>*>(fVars[i]);
317  assert(var || constrVar);
318  if (constrVar)
319  constrVar->SetValue(calc, val);
320  else
321  var->SetValue(calc, val);
322  }
323  }
324 
325  fShifts->ResetToNominal();
326  for(unsigned int j = 0; j < fSysts.size(); ++j){
327  auto val = pars[fVars.size()+j];
328  fShifts->SetShift(fSysts[j], val);
329  }
330  }
std::vector< const ISyst * > fSysts
Definition: IFitter.h:115
const XML_Char int const XML_Char * value
Definition: expat.h:331
const double j
Definition: BetheBloch.cxx:29
std::vector< const IFitVar * > fVars
Definition: IFitter.h:114
std::string pars("Th23Dmsq32")
assert(nhit_max >=nhit_nbins)
std::unique_ptr< SystShifts > fShifts
Definition: IFitter.h:116
std::vector< IFitter::SeedPt > ana::IFitter::ExpandSeeds ( const SeedList seedPts,
const std::vector< SystShifts > &  systSeedPts 
) const
protectedinherited

Definition at line 46 of file IFitter.cxx.

References ana::SeedList::GetSeeds(), ana::SystShifts::Nominal(), gen_flatrecord::pt, runNovaSAM::ret, and seed.

Referenced by ana::IFitter::FitHelper().

48  {
49  std::vector<SeedPt> ret;
50  for(Seed seed: seedPts.GetSeeds()) ret.emplace_back(seed, SystShifts::Nominal());
51 
52  // Now duplicate as many times as required for the syst seeds
53  if(!systSeedPts.empty()){
54  std::vector<SeedPt> newret;
55  for(const SystShifts& s: systSeedPts){
56  for(SeedPt pt: ret){
57  pt.shift = s;
58  newret.push_back(pt);
59  }
60  }
61  ret = newret;
62  }
63 
64  return ret;
65  }
static SystShifts Nominal()
Definition: SystShifts.h:34
const XML_Char * s
Definition: expat.h:262
unsigned int seed
Definition: runWimpSim.h:102
void ana::StanFitter::ExtractSamples ( MCMCSamples samples,
MemoryTupleWriter::WhichSamples  ws = MemoryTupleWriter::WhichSamples::kPostWarmup 
)
inline

Extract the MCMC samples from the fitter (you own them now).

Definition at line 206 of file StanFitter.h.

Referenced by test_stanfit_withsysts().

208  {
210  };
MCMCSamples fMCMCWarmup
Definition: StanFitter.h:415
MCMCSamples fMCMCSamples
Definition: StanFitter.h:414
std::unique_ptr<MCMCSamples> ana::StanFitter::ExtractSamples ( MemoryTupleWriter::WhichSamples  ws = MemoryTupleWriter::WhichSamples::kPostWarmup)
inline

Another way of extracting the MCMC samples (so you own them now) if you prefer unique_ptrs.

Definition at line 213 of file StanFitter.h.

214  {
215  return std::make_unique<MCMCSamples>(std::move(ws == MemoryTupleWriter::WhichSamples::kWarmup ? fMCMCWarmup : fMCMCSamples));
216  }
MCMCSamples fMCMCWarmup
Definition: StanFitter.h:415
MCMCSamples fMCMCSamples
Definition: StanFitter.h:414
virtual std::unique_ptr<IFitSummary> ana::IFitter::Fit ( osc::IOscCalcAdjustable seed,
SystShifts systSeed,
Verbosity  verb 
) const
inlinevirtualinherited

Variant with no seedPts.

Definition at line 61 of file IFitter.h.

References ana::IFitter::Fit().

64  {
65  return Fit(seed, systSeed, {}, std::vector<SystShifts>(1, systSeed), verb);
66  }
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
virtual std::unique_ptr<IFitSummary> ana::IFitter::Fit ( osc::IOscCalcAdjustable seed,
Verbosity  verb 
) const
inlinevirtualinherited

Variant with no seedPts and no systematics result returned.

Definition at line 69 of file IFitter.h.

References ana::IFitter::Fit().

71  {
72  return Fit(seed, junkShifts, {}, {}, verb);
73  }
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
static SystShifts junkShifts
Definition: IFitter.h:19
virtual std::unique_ptr<IFitSummary> ana::IFitter::Fit ( SystShifts systSeed,
Verbosity  verb = kVerbose 
) const
inlinevirtualinherited

Variant with no oscillations - useful for ND fits.

Definition at line 76 of file IFitter.h.

References ana::IFitter::Fit().

77  {
78  return Fit(nullptr, systSeed, {}, std::vector<SystShifts>(1, systSeed), verb);
79  }
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
std::unique_ptr< IFitter::IFitSummary > ana::StanFitter::Fit ( osc::IOscCalcAdjustable seed,
SystShifts bestSysts = junkShifts,
const SeedList seedPts = SeedList(),
const std::vector< SystShifts > &  systSeedPts = {},
Verbosity  verb = kVerbose 
) const
overridevirtual

Override (and call) the base class version, doing some caching of the oscillation parameters to avoid the problem discussed in the doc for fParamCache

Reimplemented from ana::IFitter.

Definition at line 129 of file StanFitter.cxx.

References om::cerr, allTimeWatchdog::endl, ana::IFitter::Fit(), fMCMCSamples, fMCMCWarmup, fOscCalcCache, fStanConfig, fValueWriter, ana::StanConfig::num_samples, ana::StanConfig::num_warmup, ana::MCMCSamples::NumSamples(), and seed.

Referenced by test_stanfit_withsysts().

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  }
std::size_t NumSamples() const
How many samples do we have?
Definition: MCMCSamples.h:135
OStream cerr
Definition: OStream.cxx:7
MCMCSamples fMCMCWarmup
Definition: StanFitter.h:415
MCMCSamples fMCMCSamples
Definition: StanFitter.h:414
unsigned int seed
Definition: runWimpSim.h:102
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
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
StanConfig fStanConfig
Configuration passed to Stan for fitting. See the StanConfig struct documentation for ideas...
Definition: StanFitter.h:413
std::unique_ptr< MemoryTupleWriter > fValueWriter
Definition: StanFitter.h:419
std::unique_ptr< osc::IOscCalcAdjustable > fOscCalcCache
Definition: StanFitter.h:432
int num_samples
Number of steps in the Markov chain retained for analysis (after warmup).
Definition: StanConfig.h:64
std::unique_ptr< IFitter::IFitSummary > ana::IFitter::FitHelper ( osc::IOscCalcAdjustable seed,
SystShifts bestSysts,
const SeedList seedPts,
const std::vector< SystShifts > &  systSeedPts,
Verbosity  verb 
) const
protectedvirtualinherited

Helper for Fit.

Definition at line 117 of file IFitter.cxx.

References ana::assert(), ana::SystShifts::Copy(), osc::_IOscCalcAdjustable< T >::Copy(), ana::IFitter::ExpandSeeds(), ana::IFitter::FitHelperSeeded(), ana::IFitter::fShifts, ana::IFitter::fSysts, ana::IFitter::fVars, MECModelEnuComparisons::i, gen_flatrecord::pt, seed, ana::SystShifts::SetShift(), and registry_explorer::v.

Referenced by ana::IFitter::Fit().

122  {
123  // if user passed a derived kind of SystShifts, this preserves it
124  fShifts = bestSysts.Copy();
125  fShifts->ResetToNominal();
126 
127  const std::vector<SeedPt> pts = ExpandSeeds(seedPts, systSeedPts);
128 
129  std::unique_ptr<IFitSummary> bestFitSummary;
130  std::vector<double> bestFitPars, bestSystPars;
131 
132  for (const SeedPt &pt: pts)
133  {
134  osc::IOscCalcAdjustable* seed = nullptr;
135  if(initseed) seed = initseed->Copy();
136 
137  pt.fitvars.ResetCalc(seed);
138 
139  // be sure to keep any derived class stuff around
140  auto shift = fShifts->Copy();
141  *shift = pt.shift;
142  auto fitSummary = FitHelperSeeded(seed, *shift, verb);
143  if (fitSummary->IsBetterThan(bestFitSummary.get()))
144  {
145  bestFitSummary = std::move(fitSummary);
146  // Store the best fit values of all the parameters we know are being
147  // varied.
148  bestFitPars.clear();
149  bestSystPars.clear();
150  for (const auto *v: fVars)
151  bestFitPars.push_back(v->GetValue(seed));
152  for (const auto *s: fSysts)
153  bestSystPars.push_back(shift->GetShift<double>(s));
154  }
155 
156  delete seed;
157  } // end for pt
158 
159  assert(bestFitSummary);
160  // Stuff the results of the actual best fit back into the seeds
161  for (unsigned int i = 0; i < fVars.size(); ++i)
162  fVars[i]->SetValue(initseed, bestFitPars[i]);
163  for (unsigned int i = 0; i < fSysts.size(); ++i)
164  bestSysts.SetShift(fSysts[i], bestSystPars[i]);
165 
166  return std::move(bestFitSummary);
167  }
virtual _IOscCalcAdjustable< T > * Copy() const =0
std::vector< const ISyst * > fSysts
Definition: IFitter.h:115
const XML_Char * s
Definition: expat.h:262
unsigned int seed
Definition: runWimpSim.h:102
virtual std::unique_ptr< IFitSummary > FitHelperSeeded(osc::IOscCalcAdjustable *seed, SystShifts &systSeed, Verbosity verb) const =0
Helper for FitHelper – actually does fitting. Reimplement in derived classes.
std::vector< const IFitVar * > fVars
Definition: IFitter.h:114
assert(nhit_max >=nhit_nbins)
std::vector< SeedPt > ExpandSeeds(const SeedList &seedPts, const std::vector< SystShifts > &systSeedPts) const
Definition: IFitter.cxx:46
std::unique_ptr< SystShifts > fShifts
Definition: IFitter.h:116
std::unique_ptr< IFitter::IFitSummary > ana::StanFitter::FitHelperSeeded ( osc::IOscCalcAdjustable seed,
SystShifts systSeed,
Verbosity  verb 
) const
overrideprivatevirtual

Implement workhorse method from IFitter interface.

Implements ana::IFitter.

Definition at line 160 of file StanFitter.cxx.

References ana::MCMCSamples::BestFitSampleIdx(), BuildInitContext(), calc, om::cerr, ana::SystShifts::Copy(), osc::_IOscCalcAdjustable< T >::Copy(), om::cout, CreateCalculator(), ana::StanConfig::denseMassMx, allTimeWatchdog::endl, fMCMCSamples, fMCMCWarmup, ana::IFitter::fShifts, fStanConfig, ana::IFitter::fSysts, ana::IFitter::fVars, cet::getenv(), ana::StanConfig::kEverything, ana::StanConfig::kQuiet, ana::StanConfig::kSilent, ana::StanConfig::kVerbose, ana::StanConfig::num_samples, ana::StanConfig::num_warmup, ana::MCMCSamples::NumSamples(), stan::services::error_codes::OK, launch_batch_jobs::process, makeTestPickles::return_code, ana::MCMCSamples::RunDiagnostics(), ana::MCMCSamples::SampleLL(), ana::MCMCSamples::SampleValue(), ana::SystShifts::SetShift(), ana::MCMCSamples::Systs(), registry_explorer::v, PandAna.Demos.tute_pid_validation::var, ana::MCMCSamples::Vars(), and ana::StanConfig::verbosity.

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()
std::size_t NumSamples() const
How many samples do we have?
Definition: MCMCSamples.h:135
virtual _IOscCalcAdjustable< T > * Copy() 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
MCMCSamples fMCMCWarmup
Definition: StanFitter.h:415
const std::vector< const IFitVar * > & Vars() const
Which Vars are sampled in these samples?
Definition: MCMCSamples.h:212
Verbosity verbosity
How verbose do you want me to be?
Definition: StanConfig.h:89
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
MCMCSamples fMCMCSamples
Definition: StanFitter.h:414
std::vector< const ISyst * > fSysts
Definition: IFitter.h:115
const XML_Char * s
Definition: expat.h:262
void RunDiagnostics(const StanConfig &cfg) const
Do some checks on the post-fit samples.
std::string getenv(std::string const &name)
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
const std::vector< const ana::ISyst * > & Systs() const
Which Systs are sampled in these samples?
Definition: MCMCSamples.h:209
bool denseMassMx
Should the mass matrix used in HMC be diagonal (default) or dense?
Definition: StanConfig.h:90
OStream cout
Definition: OStream.cxx:6
double SampleLL(std::size_t idx) const
Get the LL for sample number .
Definition: MCMCSamples.h:162
StanConfig fStanConfig
Configuration passed to Stan for fitting. See the StanConfig struct documentation for ideas...
Definition: StanFitter.h:413
std::vector< const IFitVar * > fVars
Definition: IFitter.h:114
void CreateCalculator(osc::IOscCalcAdjustable *seed) const
Convert a &#39;normal&#39; calculator into the Stan-aware variant used internally.
Definition: StanFitter.cxx:103
std::unique_ptr< SystShifts > fShifts
Definition: IFitter.h:116
std::size_t BestFitSampleIdx() const
Retrieve the index of the best-fit sample (lowest LL), which can then be used with the SampleLL() or ...
int num_samples
Number of steps in the Markov chain retained for analysis (after warmup).
Definition: StanConfig.h:64
void ana::StanFitter::get_dims ( std::vector< std::vector< size_t > > &  dims) const
inline

Return vector of dimensions of each parameter. (Required by Stan interface.)

Parameters
dimsReturn vector of dimensions. Always 1 for each parameter in our application.

Definition at line 156 of file StanFitter.h.

References fetch_tb_beamline_files::dims, and gen_hdf5record::names.

157  {
158  dims = std::vector<std::vector<std::size_t>>(fVars.size() + fSysts.size(), std::vector<size_t>{1});
159  }
std::vector< const ISyst * > fSysts
Definition: IFitter.h:115
std::vector< const IFitVar * > fVars
Definition: IFitter.h:114
void ana::StanFitter::get_param_names ( std::vector< std::string > &  names) const

Return names of parameters. (Required by Stan interface.)

Definition at line 290 of file StanFitter.cxx.

References ana::IFitter::fSysts, ana::IFitter::fVars, and PandAna.Demos.tute_pid_validation::var.

Referenced by constrained_param_names().

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  }
std::vector< const ISyst * > fSysts
Definition: IFitter.h:115
std::vector< const IFitVar * > fVars
Definition: IFitter.h:114
const MCMCSamples& ana::StanFitter::GetSamples ( MemoryTupleWriter::WhichSamples  ws = MemoryTupleWriter::WhichSamples::kPostWarmup) const
inline

Peruse the samples the MCMC generated. If you want to take ownership of them instead, see ExtractSamples()

Definition at line 200 of file StanFitter.h.

std::unique_ptr<SystShifts> ana::IFitter::GetSystShifts ( ) const
inlineinherited

Definition at line 82 of file IFitter.h.

References ana::IFitter::fShifts.

Referenced by plot_shifts().

82 {return fShifts->Copy();}
std::unique_ptr< SystShifts > fShifts
Definition: IFitter.h:116
template<bool propto__, bool jacobian__, typename T__ >
T__ ana::StanFitter::log_prob ( std::vector< T__ > &  params_real,
std::vector< int > &  params_i__,
std::ostream *  pstream__ = 0 
) const

Implement workhorse method that actually calculates the posterior probability from Stan interface

Template Parameters
propto__
jacobian__
T__Template type for "real" parameters. (Stan uses multiple internal types to allow for matrices and other things.)
Parameters
params_realReal-valued parameters in Stan's unconstrained parameter space.
params_i__Integer-valued parameters in Stan's unconstrained parameter space.
pstream__
Returns
template<bool propto, bool jacobian, typename T_ >
T_ ana::StanFitter::log_prob ( Eigen::Matrix< T_, Eigen::Dynamic, 1 > &  params_r,
std::ostream *  pstream = 0 
) const

Version with Matrix is req'd by Stan interface, but we won't use it.

Definition at line 421 of file StanFitter.cxx.

References MECModelEnuComparisons::i.

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  }
template<bool propto__, bool jacobian__, typename T >
T ana::StanFitter::log_prob ( std::vector< T > &  params_real,
std::vector< int > &  params_int,
std::ostream *   
) const

Definition at line 301 of file StanFitter.cxx.

References ana::assert(), om::cerr, osc::CopyParams(), om::cout, DecodePars(), e, allTimeWatchdog::endl, fCalc, fExpt, fOscCalcCache, ana::IFitter::fShifts, fStanConfig, ana::IFitter::fSysts, ana::IFitter::fVars, ana::StanFitSupport< VarClass >::GetValue(), MECModelEnuComparisons::i, compare_h5_caf::idx, ana::StanConfig::kEverything, ana::IExperiment::LogLikelihood(), ana::StanFitSupport< VarClass >::LogPrior(), stan::io::reader< T >::scalar_constrain(), stan::io::reader< T >::scalar_lub_constrain(), T, febshutoff_auto::val, PandAna.Demos.tute_pid_validation::var, and ana::StanConfig::verbosity.

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()
OStream cerr
Definition: OStream.cxx:7
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
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
std::vector< const ISyst * > fSysts
Definition: IFitter.h:115
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
const IExperiment * fExpt
Definition: StanFitter.h:412
OStream cout
Definition: OStream.cxx:6
virtual stan::math::var LogLikelihood(osc::IOscCalcAdjustableStan *osc, const SystShifts &syst=kNoShift) const
Definition: IExperiment.h:25
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
std::vector< const IFitVar * > fVars
Definition: IFitter.h:114
std::unique_ptr< osc::IOscCalcAdjustable > fOscCalcCache
Definition: StanFitter.h:432
assert(nhit_max >=nhit_nbins)
std::unique_ptr< SystShifts > fShifts
Definition: IFitter.h:116
double T
Definition: Xdiff_gwt.C:5
Float_t e
Definition: plot.C:35
size_t stan::model::prob_grad::num_params_i ( ) const
inlineinherited

Definition at line 40 of file prob_grad.hpp.

40  {
41  return param_ranges_i__.size();
42  }
std::vector< std::pair< int, int > > param_ranges_i__
Definition: prob_grad.hpp:20
size_t stan::model::prob_grad::num_params_r ( ) const
inlineinherited

Definition at line 36 of file prob_grad.hpp.

References stan::model::prob_grad::num_params_r__.

Referenced by RunHMC().

36  {
37  return num_params_r__;
38  }
std::pair<int, int> stan::model::prob_grad::param_range_i ( size_t  idx) const
inlineinherited

Definition at line 44 of file prob_grad.hpp.

References compare_h5_caf::idx.

44  {
45  return param_ranges_i__[idx];
46  }
std::vector< std::pair< int, int > > param_ranges_i__
Definition: prob_grad.hpp:20
void ana::StanFitter::ReuseWarmup ( MCMCSamples &&  warmup)

Supply samples and hyperparameters from a previously run warmup of MCMC. Burden is on the user to ensure they are compatible with the sampling required. (Call this prior to calling Fit() in order to use it.)

Definition at line 433 of file StanFitter.cxx.

References om::cerr, allTimeWatchdog::endl, fMCMCWarmup, ana::MCMCSamples::Systs(), registry_explorer::v, and ana::MCMCSamples::Vars().

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  }
OStream cerr
Definition: OStream.cxx:7
MCMCSamples fMCMCWarmup
Definition: StanFitter.h:415
const std::vector< const IFitVar * > & Vars() const
Which Vars are sampled in these samples?
Definition: MCMCSamples.h:212
const XML_Char * s
Definition: expat.h:262
const std::vector< const ana::ISyst * > & Systs() const
Which Systs are sampled in these samples?
Definition: MCMCSamples.h:209
template<typename Sampler >
template int ana::StanFitter::RunHMC< stan_dense_t > ( 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
private

Run Stan with HMC. Lots of copy-paste from stan::services::sample::hmc_nuts_diag_e_adapt() (in stan/services/sample/hmc_nuts_diag_e_adapt.hpp) so that we can customize in order to save the sampler state after warmup is finished.

Parameters
init_writerStream to write initialization info to
interruptObject whose operator()() will be called after every sample
diagStreamStream to write diagnostic info to
loggerGeneral object to write log to
init_contextMCMC initialization info
procIdProcess ID
Returns
Stan return code

Definition at line 471 of file StanFitter.cxx.

References stan::services::error_codes::CONFIG, stan::services::util::create_rng(), stan::services::util::create_unit_e_diag_inv_metric(), ana::StanConfig::delta, stan::math::domain_error(), e, ana::EigenMatrixXdFromTMatrixD(), fMCMCWarmup, fStanConfig, ana::StanConfig::gamma, ana::MCMCSamples::Hyperparams(), ana::StanConfig::init_buffer, ana::StanConfig::init_radius, stan::services::util::initialize(), ana::MCMCSamples::Hyperparameters::invMetric, std::isnan(), ana::StanConfig::kappa, test_ParserArtEvents::log, ana::StanConfig::max_depth, stan::model::prob_grad::num_params_r(), ana::StanConfig::num_warmup, stan::services::error_codes::OK, ana::StanConfig::random_seed, stan::services::util::read_diag_inv_metric(), makeTestPickles::return_code, RunSampler(), ana::MCMCSamples::Hyperparameters::stepSize, ana::StanConfig::stepsize, ana::StanConfig::stepsize_jitter, ana::StanConfig::t0, ana::StanConfig::term_buffer, stan::services::util::validate_diag_inv_metric(), and ana::StanConfig::window.

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  }
double delta
Adaptation target acceptance statistic; 0 < delta < 1.
Definition: StanConfig.h:80
std::unique_ptr< TMatrixD > invMetric
Definition: MCMCSamples.h:41
double stepsize
Integrator step size used in each simulation. Typically this is adaptive and the default value is a f...
Definition: StanConfig.h:68
MCMCSamples fMCMCWarmup
Definition: StanFitter.h:415
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
unsigned int term_buffer
Width of final fast adaptation interval.
Definition: StanConfig.h:85
double init_radius
Size of the range in unconstrained parameter space where the initial point for un-specified parameter...
Definition: StanConfig.h:62
stan::io::dump create_unit_e_diag_inv_metric(size_t num_params)
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 validate_diag_inv_metric(const Eigen::VectorXd &inv_metric, callbacks::logger &logger)
const Hyperparameters & Hyperparams() const
Definition: MCMCSamples.h:90
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
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
StanConfig fStanConfig
Configuration passed to Stan for fitting. See the StanConfig struct documentation for ideas...
Definition: StanFitter.h:413
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
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
double gamma
Adaptation regularization scale; gamma > 0.
Definition: StanConfig.h:82
unsigned int init_buffer
Width of initial fast adaptation interval.
Definition: StanConfig.h:84
Eigen::VectorXd read_diag_inv_metric(stan::io::var_context &init_context, size_t num_params, callbacks::logger &logger)
double kappa
Adaptation relaxation exponent; kappa > 0.
Definition: StanConfig.h:81
size_t num_params_r() const
Definition: prob_grad.hpp:36
Float_t e
Definition: plot.C:35
double t0
Adaptation iteration offset; t0 > 0.
Definition: StanConfig.h:83
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
template<typename Sampler >
template void ana::StanFitter::RunSampler< stan_dense_t > ( Sampler &  sampler,
std::vector< double > &  cont_vector,
boost::ecuyer1988 &  rng,
stan::callbacks::interrupt interrupt,
stan::callbacks::logger logger,
stan::callbacks::writer diagnostic_writer 
) const
private

Run the HMC sampler explicitly. Cribbed from stan::services::run_adaptive_sampler() (in stan/services/util/run_adaptive_sampler.hpp) in order to insert save/restore state code after warmup is done.

Parameters
samplerThe MCMC sampler object to use. (In Stan this is a templated type, but we'll only ever use this version)
cont_vectorStarting parameter values (transformed into Stan's internal working space)
rngThe random number generator being used
interruptObject whose operator()() will be called after every sample
loggerGeneral object to write log to
diagnostic_writerObject to write diagnostics to

Definition at line 588 of file StanFitter.cxx.

References e, febshutoff_auto::end, fMCMCWarmup, fStanConfig, fValueWriter, stan::services::util::generate_transitions(), stan::callbacks::logger::info(), ana::MemoryTupleWriter::kPostWarmup, ana::MemoryTupleWriter::kWarmup, ana::StanConfig::num_samples, ana::StanConfig::num_thin, ana::StanConfig::num_warmup, ana::MCMCSamples::NumSamples(), ana::StanConfig::refresh, ana::StanConfig::save_warmup, and febshutoff_auto::start.

Referenced by RunHMC().

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  }
std::size_t NumSamples() const
How many samples do we have?
Definition: MCMCSamples.h:135
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
MCMCSamples fMCMCWarmup
Definition: StanFitter.h:415
bool save_warmup
Save the warmup steps?
Definition: StanConfig.h:66
const XML_Char * s
Definition: expat.h:262
int refresh
Number of iterations between printout updates. Only relevant when fitting in &#39;verbose&#39; mode...
Definition: StanConfig.h:67
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
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
int num_thin
Number of Markov chain steps between saved samples when sampling.
Definition: StanConfig.h:65
std::unique_ptr< MemoryTupleWriter > fValueWriter
Definition: StanFitter.h:419
Float_t e
Definition: plot.C:35
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)
int num_samples
Number of steps in the Markov chain retained for analysis (after warmup).
Definition: StanConfig.h:64
void ana::StanFitter::SetStanConfig ( const StanConfig cfg)
inline

Change the config used for Stan. See the StanConfig struct documentation for ideas.

Definition at line 224 of file StanFitter.h.

References plotROC::cfg.

Referenced by fit_3flavor_withsysts(), and test_stanfit_withsysts().

224 { fStanConfig = cfg; }
StanConfig fStanConfig
Configuration passed to Stan for fitting. See the StanConfig struct documentation for ideas...
Definition: StanFitter.h:413
void ana::StanFitter::TestGradients ( osc::IOscCalcAdjustable seed,
SystShifts systSeed 
) const

Run Stan's test of its auto-differentiation (comparing to a finite-differences calculation)

Definition at line 708 of file StanFitter.cxx.

References BuildInitContext(), om::cerr, ana::SystShifts::Copy(), om::cout, CreateCalculator(), stan::services::diagnose::diagnose(), e, allTimeWatchdog::endl, fMCMCSamples, fMCMCWarmup, fOscCalcCache, ana::IFitter::fShifts, fStanConfig, fValueWriter, ana::StanConfig::init_radius, ana::StanConfig::num_samples, ana::StanConfig::num_warmup, ana::StanConfig::random_seed, makeTestPickles::return_code, and seed.

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  }
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
OStream cerr
Definition: OStream.cxx:7
MCMCSamples fMCMCWarmup
Definition: StanFitter.h:415
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
MCMCSamples fMCMCSamples
Definition: StanFitter.h:414
double init_radius
Size of the range in unconstrained parameter space where the initial point for un-specified parameter...
Definition: StanConfig.h:62
unsigned int seed
Definition: runWimpSim.h:102
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
OStream cout
Definition: OStream.cxx:6
StanConfig fStanConfig
Configuration passed to Stan for fitting. See the StanConfig struct documentation for ideas...
Definition: StanFitter.h:413
void CreateCalculator(osc::IOscCalcAdjustable *seed) const
Convert a &#39;normal&#39; calculator into the Stan-aware variant used internally.
Definition: StanFitter.cxx:103
std::unique_ptr< MemoryTupleWriter > fValueWriter
Definition: StanFitter.h:419
std::unique_ptr< osc::IOscCalcAdjustable > fOscCalcCache
Definition: StanFitter.h:432
std::unique_ptr< SystShifts > fShifts
Definition: IFitter.h:116
Float_t e
Definition: plot.C:35
unsigned int random_seed
Random seed used by Stan internally.
Definition: StanConfig.h:60
int num_samples
Number of steps in the Markov chain retained for analysis (after warmup).
Definition: StanConfig.h:64
template<typename T >
template void ana::StanFitter::transform_helper ( const stan::io::var_context context,
stan::io::writer< double > &  writer,
const T var 
) const
private

Helper function that actually does the unconstraining Stan needs.

Definition at line 758 of file StanFitter.cxx.

References stan::io::var_context::contains_r(), om::cout, e, allTimeWatchdog::endl, fStanConfig, ana::IConstrainedFitVar::HighLimit(), ana::StanConfig::kEverything, ana::IConstrainedFitVar::LowLimit(), stan::io::writer< T >::scalar_lub_unconstrain(), stan::io::writer< T >::scalar_unconstrain(), ana::IFitVar::ShortName(), string, febshutoff_auto::val, stan::io::var_context::vals_r(), and ana::StanConfig::verbosity.

Referenced by transform_inits().

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  }
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
Verbosity verbosity
How verbose do you want me to be?
Definition: StanConfig.h:89
void scalar_lub_unconstrain(double lb, double ub, T &y)
Definition: writer.hpp:163
virtual std::vector< double > vals_r(const std::string &name) const =0
virtual bool contains_r(const std::string &name) const =0
void scalar_unconstrain(T &y)
Definition: writer.hpp:101
const XML_Char int const XML_Char * value
Definition: expat.h:331
OStream cout
Definition: OStream.cxx:6
StanConfig fStanConfig
Configuration passed to Stan for fitting. See the StanConfig struct documentation for ideas...
Definition: StanFitter.h:413
Float_t e
Definition: plot.C:35
enum BeamMode string
void ana::StanFitter::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's unconstrained internal space.

Definition at line 805 of file StanFitter.cxx.

References stan::io::writer< T >::data_i(), stan::io::writer< T >::data_r(), ana::IFitter::fSysts, ana::IFitter::fVars, transform_helper(), and PandAna.Demos.tute_pid_validation::var.

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  }
std::vector< const ISyst * > fSysts
Definition: IFitter.h:115
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
std::vector< const IFitVar * > fVars
Definition: IFitter.h:114
void ana::StanFitter::unconstrained_param_names ( std::vector< std::string > &  param_names,
bool  include_tparams = true,
bool  include_gqs = true 
) const
inline

Return names of parameters in Stan's unconstrained space. (Required by Stan interface.)

Parameters
param_namesVector of parameter names to be filled.
include_tparamsUnused in this implementation. (Other types of Stan models need it.)
include_gqsUnused in this implementation. (Other types of Stan models need it.)

Definition at line 235 of file StanFitter.h.

References gen_hdf5record::names.

238  {
239  // other kinds of Stan models need more here, but we don't.
240  std::vector<std::string> names;
241  get_param_names(names);
242  param_names.insert(param_names.end(), names.begin(), names.end());
243  }
void get_param_names(std::vector< std::string > &names) const
Return names of parameters. (Required by Stan interface.)
Definition: StanFitter.cxx:290
void ana::IFitter::ValidateSeeds ( osc::IOscCalcAdjustable seed,
const SeedList seedPts,
const std::vector< SystShifts > &  systSeedPts 
) const
protectedinherited

Check that the seeds that were specified are compatible with the vars being fitted.

Definition at line 170 of file IFitter.cxx.

References ana::SeedList::ActiveFitVars(), om::cout, allTimeWatchdog::endl, ana::IFitter::fSysts, ana::IFitter::fVars, gen_flatrecord::pt, and registry_explorer::v.

Referenced by ana::IFitter::Fit().

173  {
174  if(!seed && !fVars.empty()){
175  std::cout << "ERROR: MinuitFitter::Fit() trying to fit oscillation parameters without an oscillation calculator" << std::endl;
176  abort();
177  }
178 
179  for(const IFitVar* v: seedPts.ActiveFitVars()){
180  if (std::find(fVars.begin(), fVars.end(), v) == fVars.end()){
181  std::cout << "ERROR MinuitFitter::Fit() trying to seed '"
182  << v->ShortName()
183  << "' which is not part of the fit." << std::endl;
184  abort();
185  }
186  }
187 
188  for(const SystShifts& pt: systSeedPts){
189  for(const ISyst* s: pt.ActiveSysts()){
190  if(std::find(fSysts.begin(), fSysts.end(), s) == fSysts.end()){
191  std::cout << "ERROR MinuitFitter::Fit() trying to seed '"
192  << s->ShortName()
193  << "' which is not part of the fit." << std::endl;
194  abort();
195  }
196  }
197  }
198  }
std::vector< const ISyst * > fSysts
Definition: IFitter.h:115
const XML_Char * s
Definition: expat.h:262
OStream cout
Definition: OStream.cxx:6
std::vector< const IFitVar * > fVars
Definition: IFitter.h:114
template<typename RNG >
void ana::StanFitter::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 at line 825 of file StanFitter.cxx.

References om::cout, allTimeWatchdog::endl, fStanConfig, ana::IFitter::fVars, ana::IConstrainedFitVar::HighLimit(), MECModelEnuComparisons::i, ana::StanConfig::kEverything, ana::IConstrainedFitVar::LowLimit(), stan::io::reader< T >::scalar_constrain(), stan::io::reader< T >::scalar_lub_constrain(), ana::IFitVar::ShortName(), and ana::StanConfig::verbosity.

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  }
Verbosity verbosity
How verbose do you want me to be?
Definition: StanConfig.h:89
const std::map< std::pair< std::string, std::string >, Variable > vars
OStream cout
Definition: OStream.cxx:6
StanConfig fStanConfig
Configuration passed to Stan for fitting. See the StanConfig struct documentation for ideas...
Definition: StanFitter.h:413
std::vector< const IFitVar * > fVars
Definition: IFitter.h:114

Member Data Documentation

std::unique_ptr<osc::IOscCalcAdjustableStan> ana::StanFitter::fCalc
mutableprivate

Definition at line 411 of file StanFitter.h.

Referenced by CreateCalculator(), and log_prob().

const IExperiment* ana::StanFitter::fExpt
private

Definition at line 412 of file StanFitter.h.

Referenced by log_prob().

MCMCSamples ana::StanFitter::fMCMCSamples
private

Definition at line 414 of file StanFitter.h.

Referenced by Fit(), FitHelperSeeded(), and TestGradients().

MCMCSamples ana::StanFitter::fMCMCWarmup
private

Definition at line 415 of file StanFitter.h.

Referenced by Fit(), FitHelperSeeded(), ReuseWarmup(), RunHMC(), RunSampler(), and TestGradients().

std::unique_ptr<osc::IOscCalcAdjustable> ana::StanFitter::fOscCalcCache
mutableprivate

stan::math::var objects have one 'gotcha' associated with them: after the gradient of the log-prob is calculated, Stan internally 'recovers' the memory associated with every stan::math::var's value by calling recover_memory() (inside log_prob_grad()). This is no problem for any variables that are computed every fit iteration, like all the ones we're fitting. But any unfitted parameters that happen to have stan::math::var type (i.e: unfitted oscillation params) get the rug yanked from under them. This cache is made at the beginning from all those values and any that aren't part of the values passed by the fitter are rewritten at the beginning of each iteration.

Definition at line 432 of file StanFitter.h.

Referenced by Fit(), log_prob(), and TestGradients().

std::unique_ptr<SystShifts> ana::IFitter::fShifts
mutableprotectedinherited
StanConfig ana::StanFitter::fStanConfig
private

Configuration passed to Stan for fitting. See the StanConfig struct documentation for ideas.

Definition at line 413 of file StanFitter.h.

Referenced by BuildInitContext(), Fit(), FitHelperSeeded(), log_prob(), RunHMC(), RunSampler(), TestGradients(), transform_helper(), and write_array().

std::vector<const ISyst*> ana::IFitter::fSysts
protectedinherited
std::unique_ptr<MemoryTupleWriter> ana::StanFitter::fValueWriter
mutableprivate

See MemoryTupleWriter class documentation for more info. Pointer so it can be initialized lazily

Definition at line 419 of file StanFitter.h.

Referenced by Fit(), RunSampler(), and TestGradients().

std::vector<const IFitVar*> ana::IFitter::fVars
protectedinherited
Verbosity ana::IFitter::fVerb
mutableprotectedinherited

Definition at line 117 of file IFitter.h.

Referenced by ana::IFitter::Fit(), and ana::MinuitFitter::FitHelperSeeded().

SystShifts ana::IFitter::junkShifts = SystShifts()
staticprotectedinherited

Definition at line 19 of file IFitter.h.

size_t stan::model::prob_grad::num_params_r__
protectedinherited
std::vector<std::pair<int, int> > stan::model::prob_grad::param_ranges_i__
protectedinherited

Definition at line 20 of file prob_grad.hpp.


The documentation for this class was generated from the following files: