5 #pragma GCC diagnostic push 7 #pragma GCC diagnostic ignored "-Wmisleading-indentation" 9 #pragma GCC diagnostic ignored "-Wunused-function" 18 #pragma GCC diagnostic pop 34 #include "Utilities/func/MathUtil.h" 45 std::vector<const IFitVar *>
vars,
46 std::vector<const ISyst *>
systs)
51 fMCMCSamples(vars, systs),
52 fMCMCWarmup(vars, systs)
60 std::vector<double> vals;
61 std::vector<std::string>
names;
69 vals.push_back(util::GetValAs<double>(
v->GetValue(seed)));
70 names.push_back(
v->ShortName());
73 assert(vals.size() == fVars.size());
77 names.push_back(
s->ShortName());
80 assert(vals.size() == (fVars.size() + fSysts.size()) && vals.size() == names.size());
86 std::vector<std::vector<std::size_t>> var_sizes(names.size(), std::vector<size_t>{1});
89 return std::move(ret);
94 bool include_gqs)
const 97 std::vector<std::string>
names;
99 param_names.insert(param_names.end(), names.begin(), names.end());
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>();
128 std::unique_ptr<IFitter::IFitSummary>
132 const std::vector<SystShifts> &systSeedPts,
145 std::cerr <<
"You supplied a previous collection of MCMC samples for warmup and also requested num_warmup > 0 in your StanConfig." <<
std::endl;
154 return IFitter::Fit(seed, bestSysts, seedPts, systSeedPts, verb);
159 std::unique_ptr<IFitter::IFitSummary>
173 std::ostream nullStream(
nullptr);
175 std::unique_ptr<stan::callbacks::stream_logger> logger;
186 logger = std::make_unique<stan::callbacks::stream_logger>(nullStream,
std::cout,
std::cout,
191 logger = std::make_unique<stan::callbacks::stream_logger>(nullStream, nullStream,
std::cout,
197 logger = std::make_unique<stan::callbacks::stream_logger>(nullStream, nullStream, nullStream,
207 std::unique_ptr<stan::io::array_var_context> init_context;
209 init_context = std::make_unique<stan::io::array_var_context>(
BuildInitContext(seed, systSeed));
213 std::unique_ptr<osc::IOscCalcAdjustable>
calc(seed->
Copy());
216 auto shifts = systSeed.
Copy();
219 init_context = std::make_unique<stan::io::array_var_context>(
BuildInitContext(
calc.get(), *shifts));
224 unsigned int procId = 0;
227 procId = std::stoul(process);
237 return_code = RunHMC<stan_dense_t>(init_writer, interrupt, diagStream, logger, *init_context, procId);
239 return_code = RunHMC<stan_diag_t>(init_writer, interrupt, diagStream, logger, *init_context, procId);
281 for (
const auto & syst :
fSysts)
294 names.push_back(
var->ShortName());
295 for (
const auto & syst :
fSysts)
296 names.push_back(syst->ShortName());
300 template <
bool propto__,
bool jacobian__,
typename T>
302 std::vector<int>& params_int,
320 std::size_t numVals =
fVars.size() +
fSysts.size();
321 std::vector<T> params;
322 params.reserve(numVals);
333 util::GetValAs<T>(constrVar->HighLimit()), logprob);
336 util::GetValAs<T>(constrVar->HighLimit()));
346 std::cout <<
"var " <<
var->ShortName() <<
" has val = " << util::GetValAs<double>(
val) << std::endl;
347 params.push_back(val);
358 std::cout <<
"syst " <<
fSysts[
idx]->ShortName() <<
" has val = " << util::GetValAs<double>(
val) << std::endl;
359 params.push_back(val);
363 std::cout <<
"log-prob from re-entering constrained space = " << logprob <<
std::endl;
374 for (
unsigned int i = 0;
i < fVars.size(); ++
i)
380 auto prior = constrVar ? constrVar->
LogPrior(params[
i],
fCalc.get())
381 :
var->LogPrior(params[i],
fCalc.get());
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;
389 logprob += util::GetValAs<T>(prior);
391 logprob += util::GetValAs<T>(
fShifts->LogPrior());
393 std::cout <<
"log-prob from syst prior " << util::GetValAs<double>(
fShifts->LogPrior()) << std::endl;
398 logprob += util::GetValAs<T>(ll);
404 std::cerr <<
"Exception during evaluation of likelihood in Stan fitting (see StanFitter::log_prob()): " <<
std::endl;
414 fCalc->InvalidateCache();
420 template <
bool propto,
bool jacobian,
typename T_>
422 std::ostream* pstream)
const 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);
439 std::cerr <<
"Warmup to be reused is incompatible with this StanFitter." <<
std::endl;
444 std::cerr <<
" Supplied warmup fit vars: ";
445 for (
const auto &
v : warmup.Vars())
453 std::cerr <<
"Warmup to be reused is incompatible with this StanFitter." <<
std::endl;
459 for (
const auto &
s : warmup.Systs())
470 template <
typename Sampler>
473 std::ostream &diagStream,
474 const std::unique_ptr<stan::callbacks::stream_logger> &logger,
476 unsigned int procId)
const 502 Eigen::VectorXd inv_metric;
524 Sampler sampler(*
this, rng);
526 sampler.set_metric(inv_metric);
546 RunSampler(sampler, cont_vector, rng, interrupt, *logger, diagnostic_writer);
572 std::ostream &diagStream,
573 const std::unique_ptr<stan::callbacks::stream_logger> &logger,
575 unsigned int procId)
const;
578 std::ostream &diagStream,
579 const std::unique_ptr<stan::callbacks::stream_logger> &logger,
581 unsigned int procId)
const;
587 template <
typename Sampler>
589 std::vector<double>& cont_vector,
590 boost::ecuyer1988& rng,
596 auto &
model = *
this;
600 Eigen::Map<Eigen::VectorXd> cont_params(cont_vector.data(),
607 double warm_delta_t = 0;
611 sampler.engage_adaptation();
614 sampler.z().q = cont_params;
615 sampler.init_stepsize(logger);
618 logger.
info(
"Exception initializing step size.");
619 logger.
info(e.what());
630 writer.write_sample_names(s, sampler,
model);
631 writer.write_diagnostic_names(s, sampler,
model);
649 warm_delta_t =
static_cast<double>(end -
start) / CLOCKS_PER_SEC;
651 sampler.disengage_adaptation();
652 writer.write_adapt_finish(sampler);
663 writer.write_sample_names(s, sampler,
model);
664 writer.write_diagnostic_names(s, sampler,
model);
683 double sample_delta_t =
static_cast<double>(end -
start) / CLOCKS_PER_SEC;
685 writer.write_timing(warm_delta_t, sample_delta_t);
689 fValueWriter->SaveSamplerState(sampler, sample_delta_t);
694 template void StanFitter::RunSampler<stan_diag_t>(
stan_diag_t& sampler,
695 std::vector<double>& cont_vector,
696 boost::ecuyer1988& rng,
700 template void StanFitter::RunSampler<stan_dense_t>(
stan_dense_t& sampler,
701 std::vector<double>& cont_vector,
702 boost::ecuyer1988& rng,
738 const double error = 1
e-6;
753 <<
" parameters were not within " << error <<
" of finite diff expectation..." <<
std::endl;
757 template <
typename T>
763 "transform_helper() can only be instantiated with IFitVar* or ISyst* as the type");
766 throw std::runtime_error(
"variable '" + var->ShortName() +
"' missing from context given to Stan!");
768 std::vector<double> vals_real = context.
vals_r(var->ShortName());
769 double val = vals_real[0];
774 if ( (constrFitVar = dynamic_cast<const IConstrainedFitVar*>(var)) )
779 <<
" from constrained between " << constrFitVar->
LowLimit() <<
" and " << constrFitVar->
HighLimit()
783 util::GetValAs<double>(constrFitVar->
HighLimit()), val);
791 throw std::runtime_error(
std::string(
"Error transforming variable '" + var->ShortName() +
"': ") + e.what());
799 const IFitVar *
const &)
const;
802 const ISyst *
const &)
const;
806 std::vector<int>& params_int,
807 std::vector<double>& params_real,
815 for (
const auto & syst :
fSysts)
818 params_real = stanWriter.
data_r();
819 params_int = stanWriter.
data_i();
824 template <
typename RNG>
826 std::vector<double>& params_real,
827 std::vector<int>& params_int,
828 std::vector<double>&
vars,
835 for (std::size_t
i = 0;
i < params_real.size();
i++)
839 constrFitVar = dynamic_cast<const IConstrainedFitVar*>(
fVars[
i]);
846 <<
" from unbounded space" 847 <<
" to constrained between " << constrFitVar->
LowLimit() <<
" and " << constrFitVar->
HighLimit()
851 util::GetValAs<double>(constrFitVar->
HighLimit()));
857 vars.push_back(parVal);
double delta
Adaptation target acceptance statistic; 0 < delta < 1.
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 ¶meter_writer)
std::unique_ptr< TMatrixD > invMetric
virtual double HighLimit() const =0
T scalar_lub_constrain(const TL lb, const TU ub)
std::size_t NumSamples() const
How many samples do we have?
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.
virtual std::unique_ptr< SystShifts > Copy() const
void get_param_names(std::vector< std::string > &names) const
Return names of parameters. (Required by Stan interface.)
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
std::string DemangledTypeName(T *obj)
utility method to figure out exactly what kind of object I am
Simple record of shifts applied to systematic parameters.
double stepsize
Integrator step size used in each simulation. Typically this is adaptive and the default value is a f...
std::vector< T > & data_r()
virtual double LowLimit() const =0
double SampleValue(const IFitVar *var, std::size_t idx) const
Get the value of Var var for sample number .
::xsd::cxx::tree::exception< char > exception
const std::vector< const IFitVar * > & Vars() const
Which Vars are sampled in these samples?
void write_array(RNG &base_rng__, std::vector< double > ¶ms_real, std::vector< int > ¶ms_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.
Verbosity verbosity
How verbose do you want me to be?
void DecodePars(const std::vector< T > &pars, osc::IOscCalcAdjustableStan *calc) const
void scalar_lub_unconstrain(double lb, double ub, T &y)
bool save_warmup
Save the warmup steps?
stan::io::array_var_context BuildInitContext(osc::IOscCalcAdjustable *seed, const SystShifts &systSeed) const
Helper function to build the initial seed point context for Stan.
virtual std::vector< double > vals_r(const std::string &name) const =0
Encapsulate code to systematically shift a caf::SRProxy.
virtual bool contains_r(const std::string &name) const =0
int isnan(const stan::math::var &a)
void scalar_unconstrain(T &y)
T__ log_prob(std::vector< T__ > ¶ms_real, std::vector< int > ¶ms_i__, std::ostream *pstream__=0) const
std::vector< const ISyst * > fSysts
unsigned int term_buffer
Width of final fast adaptation interval.
double init_radius
Size of the range in unconstrained parameter space where the initial point for un-specified parameter...
StanFitter(const IExperiment *expt, std::vector< const IFitVar * > vars, std::vector< const ISyst * > systs={})
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 > ¶ms_int, std::vector< double > ¶ms_real, std::ostream *pstream__) const
Take the seed points and convert them into Stan's unconstrained internal space.
const XML_Char int const XML_Char * value
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
int refresh
Number of iterations between printout updates. Only relevant when fitting in 'verbose' mode...
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.
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.
T GetShift(const ISyst *syst) const
int num_warmup
Number of initial steps in the Markov chain (used to enter the typical set). These are usually discar...
void constrained_param_names(std::vector< std::string > ¶m_names, bool include_tparams=true, bool include_gqs=true) const
Return names of parameters . (Required by Stan interface.)
const IExperiment * fExpt
void validate_diag_inv_metric(const Eigen::VectorXd &inv_metric, callbacks::logger &logger)
void ReuseWarmup(MCMCSamples &&warmup)
const Hyperparameters & Hyperparams() const
const std::vector< const ana::ISyst * > & Systs() const
Which Systs are sampled in these samples?
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
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)
bool denseMassMx
Should the mass matrix used in HMC be diagonal (default) or dense?
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.
Eigen::MatrixXd EigenMatrixXdFromTMatrixD(const TMatrixD *mat)
void domain_error(const char *function, const char *name, const T &y, const char *msg1, const char *msg2)
virtual stan::math::var LogPrior(const stan::math::var &var, const osc::IOscCalcAdjustableStan *calc) const
double SampleLL(std::size_t idx) const
Get the LL for sample number .
virtual stan::math::var LogLikelihood(osc::IOscCalcAdjustableStan *osc, const SystShifts &syst=kNoShift) const
virtual void info(const std::string &message)
StanConfig fStanConfig
Configuration passed to Stan for fitting. See the StanConfig struct documentation for ideas...
std::unique_ptr< osc::IOscCalcAdjustableStan > fCalc
double stepsize_jitter
Randomly "jitter" the step size by this fraction of the step size. (Sometimes can help if the integra...
std::vector< const IFitVar * > fVars
int max_depth
Depth of the binary tree used by the HMC when taking its leapfrog steps. Note: the number of leapfrog...
const std::string & ShortName() const
double gamma
Adaptation regularization scale; gamma > 0.
std::unique_ptr< IFitSummary > FitHelperSeeded(osc::IOscCalcAdjustable *seed, SystShifts &systSeed, Verbosity verb) const override
Implement workhorse method from IFitter interface.
void CreateCalculator(osc::IOscCalcAdjustable *seed) const
Convert a 'normal' calculator into the Stan-aware variant used internally.
unsigned int init_buffer
Width of initial fast adaptation interval.
int num_thin
Number of Markov chain steps between saved samples when sampling.
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.
std::unique_ptr< MemoryTupleWriter > fValueWriter
void TestGradients(osc::IOscCalcAdjustable *seed, SystShifts &systSeed) const
Run Stan's test of its auto-differentiation (comparing to a finite-differences calculation) ...
std::unique_ptr< osc::IOscCalcAdjustable > fOscCalcCache
assert(nhit_max >=nhit_nbins)
Interface definition for fittable variables.
double kappa
Adaptation relaxation exponent; kappa > 0.
std::unique_ptr< SystShifts > fShifts
size_t num_params_r() const
std::vector< int > & data_i()
double t0
Adaptation iteration offset; t0 > 0.
const XML_Char XML_Content * model
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)
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.
boost::ecuyer1988 create_rng(unsigned int seed, unsigned int chain)
int num_samples
Number of steps in the Markov chain retained for analysis (after warmup).