MaCh3 2.2.1
Reference Guide
Loading...
Searching...
No Matches
mcmc.cpp
Go to the documentation of this file.
1#include "mcmc.h"
2
3// *************************
4// Initialise the manager and make it an object of mcmc class
5// Now we can dump manager settings to the output file
7// *************************
8 // Beginning step number
9 stepStart = 0;
10
11 // Starting parameters should be thrown
12 reject = false;
13 chainLength = Get<unsigned>(fitMan->raw()["General"]["MCMC"]["NSteps"], __FILE__, __LINE__);
14
15 AnnealTemp = GetFromManager<double>(fitMan->raw()["General"]["MCMC"]["AnnealTemp"], -999);
16 if(AnnealTemp < 0) anneal = false;
17 else
18 {
19 MACH3LOG_INFO("Enabling simulated annealing with T = {}", AnnealTemp);
20 anneal = true;
21 }
22}
23
24// *************************
25// Destructor: close the logger and output file
27// *************************
28
29}
30
31// **********************
32// Do we accept the proposed step for all the parameters?
34// **********************
35 bool accept = false;
36
37 // Set the acceptance probability to zero
38 accProb = 0.0;
39
40 // Calculate acceptance probability
41 if (anneal) accProb = std::min(1., std::exp( -(logLProp-logLCurr) / (std::exp(-step/AnnealTemp))));
42 else accProb = std::min(1., std::exp(logLCurr-logLProp));
43
44 // Get the random number
45 const double fRandom = random->Rndm();
46
47 // Do the accept/reject
48 if (fRandom <= accProb) {
49 accept = true;
50 ++accCount;
51 } else {
52 accept = false;
53 }
54
55 #ifdef DEBUG
56 if (debug) debugFile << " logLProp: " << logLProp << " logLCurr: " << logLCurr << " accProb: " << accProb << " fRandom: " << fRandom << std::endl;
57 #endif
58
59 // Update all the handlers to accept the step
60 if (accept && !reject) {
62
63 // Loop over systematics and accept
64 for (size_t s = 0; s < systematics.size(); ++s) {
65 systematics[s]->AcceptStep();
66 }
67 }
68
69 stepClock->Stop();
70 stepTime = stepClock->RealTime();
71
72 // Write step to output tree
73 outTree->Fill();
74}
75
76// *******************
77// Run the Markov chain with all the systematic objects added
79// *******************
80 // Save the settings into the output file
82
83 // Prepare the output branches
85
86 // Remove obsolete memory and make other checks before fit starts
88
89 // Reconfigure the samples, systematics and oscillation for first weight
90 // ProposeStep sets logLProp
92 // Set the current logL to the proposed logL for the 0th step
93 // Accept the first step to set logLCurr: this shouldn't affect the MCMC because we ignore the first N steps in burn-in
95
96 // Begin MCMC
97 const auto StepEnd = stepStart + chainLength;
98 for (step = stepStart; step < StepEnd; ++step)
99 {
100 stepClock->Start();
101 // Set the initial rejection to false
102 reject = false;
103
104 // Print 10 steps in total
105 if ( (step-stepStart) % (chainLength/10) == 0) {
107 }
108
109 // Propose current step variation and save the systematic likelihood that results in this step being taken
110 // Updates logLProp
111 ProposeStep();
112
113 // Does the MCMC accept this step?
114 CheckStep();
115
116 // Auto save the output
117 if (step % auto_save == 0) outTree->AutoSave();
118 }
119
120 // Save all the MCMC output
121 SaveOutput();
122
123 //Save the Adaptive output
124 for(const auto& syst : systematics){
125 if(syst->GetDoAdaption()){
126 syst->GetAdaptiveHandler()->SaveAdaptiveToFile(syst->GetAdaptiveHandler()->GetOutFileName(), syst->GetName(), true);
127 }
128 }
129
130 // Process MCMC
131 ProcessMCMC();
132}
133
134// *******************
135// Do the initial reconfigure of the MCMC
137// *******************
138 // Initial likelihood
139 double llh = 0.0;
140
141 // Initiate to false
142 reject = false;
143
144 // Loop over the systematics and propose the initial step
145 for (size_t s = 0; s < systematics.size(); ++s) {
146 // Could throw the initial value here to do MCMC stability studies
147 // Propose the steps for the systematics
148 systematics[s]->ProposeStep();
149
150 // Get the likelihood from the systematics
151 syst_llh[s] = systematics[s]->GetLikelihood();
152 llh += syst_llh[s];
153
154 #ifdef DEBUG
155 if (debug) debugFile << "LLH after " << systematics[s]->GetName() << " " << llh << std::endl;
156 #endif
157 }
158
159 // Check if we've hit a boundary in the systematics
160 // In this case we can save time by not having to reconfigure the simulation
161 if (llh >= M3::_LARGE_LOGL_) {
162 reject = true;
163 #ifdef DEBUG
164 if (debug) debugFile << "Rejecting based on boundary" << std::endl;
165 #endif
166 }
167
168 // Only reweight when we have a good parameter configuration
169 // This speeds things up considerably because for every bad parameter configuration we don't have to reweight the MC
170 if (!reject)
171 {
172 // Could multi-thread this
173 // But since sample reweight is multi-threaded it's probably better to do that
174 for (size_t i = 0; i < samples.size(); ++i)
175 {
176 samples[i]->Reweight();
177 }
178
179 //DB for atmospheric event by event sample migration, need to fully reweight all samples to allow event passing prior to likelihood evaluation
180 for (size_t i = 0; i < samples.size(); ++i) {
181 // Get the sample likelihoods and add them
182 sample_llh[i] = samples[i]->GetLikelihood();
183 llh += sample_llh[i];
184 #ifdef DEBUG
185 if (debug) debugFile << "LLH after sample " << i << " " << llh << std::endl;
186 #endif
187 }
188
189 // For when we don't have to reweight, set sample to madness
190 } else {
191 for (size_t i = 0; i < samples.size(); ++i) {
192 // Set the sample_llh[i] to be madly high also to signify a step out of bounds
194 #ifdef DEBUG
195 if (debug) debugFile << "LLH after REJECT sample " << i << " " << llh << std::endl;
196 #endif
197 }
198 }
199
200 // Save the proposed likelihood (class member)
201 logLProp = llh;
202}
203
204// *******************
205// Print the fit output progress
207// *******************
208 MACH3LOG_INFO("Step:\t{}/{}, current: {:.2f}, proposed: {:.2f}", step - stepStart, chainLength, logLCurr, logLProp);
209 MACH3LOG_INFO("Accepted/Total steps: {}/{} = {:.2f}", accCount, step - stepStart, static_cast<double>(accCount) / static_cast<double>(step - stepStart));
210
211 for (ParameterHandlerBase *cov : systematics) {
212 cov->PrintNominalCurrProp();
213 }
214 #ifdef DEBUG
215 if (debug) {
216 debugFile << "\n-------------------------------------------------------" << std::endl;
217 debugFile << "Step:\t" << step + 1 << "/" << chainLength << " | current: " << logLCurr << " proposed: " << logLProp << std::endl;
218 }
219 #endif
220}
221
222// *******************
223void mcmc::StartFromPreviousFit(const std::string& FitName) {
224// *******************
225 // Use base class
227
228 // For MCMC we also need to set stepStart
229 TFile *infile = new TFile(FitName.c_str(), "READ");
230 TTree *posts = infile->Get<TTree>("posteriors");
231 int step_val = -1;
232
233 posts->SetBranchAddress("step",&step_val);
234 posts->GetEntry(posts->GetEntries()-1);
235
236 stepStart = step_val + 1;
237 // KS: Also update number of steps if using adaption
238 for(unsigned int i = 0; i < systematics.size(); ++i){
239 if(systematics[i]->GetDoAdaption()){
240 systematics[i]->SetNumberOfSteps(step_val);
241 }
242 }
243 infile->Close();
244 delete infile;
245}
MaCh3Plotting::PlottingManager * man
#define MACH3LOG_INFO
Definition: MaCh3Logger.h:23
Base class for implementing fitting algorithms.
Definition: FitterBase.h:23
std::unique_ptr< TRandom3 > random
Random number.
Definition: FitterBase.h:128
std::vector< SampleHandlerBase * > samples
Sample holder.
Definition: FitterBase.h:113
int stepStart
step start if restarting
Definition: FitterBase.h:105
double logLProp
proposed likelihood
Definition: FitterBase.h:99
void ProcessMCMC()
Process MCMC output.
Definition: FitterBase.cpp:371
int accCount
counts accepted steps
Definition: FitterBase.h:103
void SaveOutput()
Save output and close files.
Definition: FitterBase.cpp:221
void SaveSettings()
Save the settings that the MCMC was run with.
Definition: FitterBase.cpp:82
manager * fitMan
The manager.
Definition: FitterBase.h:92
unsigned int step
current state
Definition: FitterBase.h:95
void PrepareOutput()
Prepare the output file.
Definition: FitterBase.cpp:148
double accProb
current acceptance prob
Definition: FitterBase.h:101
virtual void StartFromPreviousFit(const std::string &FitName)
Allow to start from previous fit/chain.
Definition: FitterBase.cpp:295
std::vector< double > sample_llh
store the llh breakdowns
Definition: FitterBase.h:108
double stepTime
Time of single step.
Definition: FitterBase.h:125
std::unique_ptr< TStopwatch > stepClock
tells how long single step/fit iteration took
Definition: FitterBase.h:123
double logLCurr
current likelihood
Definition: FitterBase.h:97
std::vector< double > syst_llh
systematic llh breakdowns
Definition: FitterBase.h:110
int auto_save
auto save every N steps
Definition: FitterBase.h:139
TTree * outTree
Output tree with posteriors.
Definition: FitterBase.h:137
void SanitiseInputs()
Remove obsolete memory and make other checks before fit starts.
Definition: FitterBase.cpp:213
std::vector< ParameterHandlerBase * > systematics
Systematic holder.
Definition: FitterBase.h:118
Base class responsible for handling of systematic error parameters. Capable of using PCA or using ada...
The manager class is responsible for managing configurations and settings.
Definition: Manager.h:16
YAML::Node const & raw()
Return config.
Definition: Manager.h:41
void RunMCMC() override
Actual implementation of MCMC fitting algorithm.
Definition: mcmc.cpp:78
unsigned int chainLength
number of steps in chain
Definition: mcmc.h:41
double AnnealTemp
simulated annealing temperature
Definition: mcmc.h:46
bool anneal
simulated annealing
Definition: mcmc.h:44
bool reject
Do we reject based on hitting boundaries in systs.
Definition: mcmc.h:39
mcmc(manager *const fitMan)
Constructor.
Definition: mcmc.cpp:6
virtual ~mcmc()
Destructor.
Definition: mcmc.cpp:26
void CheckStep()
Do we accept the step.
Definition: mcmc.cpp:33
void ProposeStep()
Propose a step.
Definition: mcmc.cpp:136
void StartFromPreviousFit(const std::string &FitName) override
Allow to start from previous fit/chain.
Definition: mcmc.cpp:223
void PrintProgress()
Print the progress.
Definition: mcmc.cpp:206
static constexpr const double _LARGE_LOGL_
Large Likelihood is used it parameter go out of physical boundary, this indicates in MCMC that such s...
Definition: Core.h:70