MaCh3 2.2.1
Reference Guide
Loading...
Searching...
No Matches
ParameterHandlerBase.cpp
Go to the documentation of this file.
2#include "regex"
3
4// ********************************************
5ParameterHandlerBase::ParameterHandlerBase(std::string name, std::string file, double threshold, int FirstPCA, int LastPCA)
6 : inputFile(file), pca(false) {
7// ********************************************
8 MACH3LOG_DEBUG("Constructing instance of ParameterHandler");
10 if (threshold < 0 || threshold >= 1) {
11 MACH3LOG_INFO("NOTE: {} {}", name, file);
12 MACH3LOG_INFO("Principal component analysis but given the threshold for the principal components to be less than 0, or greater than (or equal to) 1. This will not work");
13 MACH3LOG_INFO("Please specify a number between 0 and 1");
14 MACH3LOG_INFO("You specified: ");
15 MACH3LOG_INFO("Am instead calling the usual non-PCA constructor...");
16 pca = false;
17 }
18 Init(name, file);
19
20 // Call the innocent helper function
21 if (pca) ConstructPCA(threshold, FirstPCA, LastPCA);
22}
23// ********************************************
24ParameterHandlerBase::ParameterHandlerBase(const std::vector<std::string>& YAMLFile, std::string name, double threshold, int FirstPCA, int LastPCA)
25 : inputFile(YAMLFile[0].c_str()), matrixName(name), pca(true) {
26// ********************************************
27 MACH3LOG_INFO("Constructing instance of ParameterHandler using");
29 for(unsigned int i = 0; i < YAMLFile.size(); i++)
30 {
31 MACH3LOG_INFO("{}", YAMLFile[i]);
32 }
33 MACH3LOG_INFO("as an input");
34
35 if (threshold < 0 || threshold >= 1) {
36 MACH3LOG_INFO("Principal component analysis but given the threshold for the principal components to be less than 0, or greater than (or equal to) 1. This will not work");
37 MACH3LOG_INFO("Please specify a number between 0 and 1");
38 MACH3LOG_INFO("You specified: ");
39 MACH3LOG_INFO("Am instead calling the usual non-PCA constructor...");
40 pca = false;
41 }
42
43 Init(YAMLFile);
44 // Call the innocent helper function
45 if (pca) ConstructPCA(threshold, FirstPCA, LastPCA);
46}
47
48// ********************************************
49//Destructor
51// ********************************************
52 delete[] randParams;
53 delete[] corr_throw;
54
55 if (covMatrix != nullptr) delete covMatrix;
56 if (invCovMatrix != nullptr) delete invCovMatrix;
58 if (throwMatrix != nullptr) delete throwMatrix;
59 for(int i = 0; i < _fNumPar; i++) {
60 delete[] throwMatrixCholDecomp[i];
61 }
62 delete[] throwMatrixCholDecomp;
63}
64
65// ********************************************
66void ParameterHandlerBase::ConstructPCA(const double eigen_threshold, int FirstPCAdpar, int LastPCAdpar) {
67// ********************************************
68 if(AdaptiveHandler) {
69 MACH3LOG_ERROR("Adaption has been enabled and now trying to enable PCA. Right now both configuration don't work with each other");
70 throw MaCh3Exception(__FILE__ , __LINE__ );
71 }
72
73 PCAObj = std::make_unique<PCAHandler>();
74 //Check whether first and last pcadpar are set and if not just PCA everything
75 if(FirstPCAdpar == -999 || LastPCAdpar == -999){
76 if(FirstPCAdpar == -999 && LastPCAdpar == -999){
77 FirstPCAdpar = 0;
78 LastPCAdpar = covMatrix->GetNrows()-1;
79 }
80 else{
81 MACH3LOG_ERROR("You must either leave FirstPCAdpar and LastPCAdpar at -999 or set them both to something");
82 throw MaCh3Exception(__FILE__ , __LINE__ );
83 }
84 }
85
86 PCAObj->ConstructPCA(covMatrix, FirstPCAdpar, LastPCAdpar, eigen_threshold, _fNumPar);
87 PCAObj->SetupPointers(&_fCurrVal, &_fPropVal);
88 // Make a note that we have now done PCA
89 pca = true;
90}
91
92// ********************************************
93void ParameterHandlerBase::Init(const std::string& name, const std::string& file) {
94// ********************************************
95 // Set the covariance matrix from input ROOT file (e.g. flux, ND280, NIWG)
96 TFile *infile = new TFile(file.c_str(), "READ");
97 if (infile->IsZombie()) {
98 MACH3LOG_ERROR("Could not open input covariance ROOT file {} !!!", file);
99 MACH3LOG_ERROR("Was about to retrieve matrix with name {}", name);
100 throw MaCh3Exception(__FILE__ , __LINE__ );
101 }
102
103 TMatrixDSym *CovMat = static_cast<TMatrixDSym*>(infile->Get(name.c_str()));
104
105 if (!CovMat) {
106 MACH3LOG_ERROR("Could not find covariance matrix name {} in file {}", name, file);
107 MACH3LOG_ERROR("Are you really sure {} exists in the file?", name);
108 throw MaCh3Exception(__FILE__ , __LINE__ );
109 }
110
111 PrintLength = 35;
112
113 const int nThreads = M3::GetNThreads();
114 //KS: set Random numbers for each thread so each thread has different seed
115 //or for one thread if without MULTITHREAD
116 random_number.reserve(nThreads);
117 for (int iThread = 0; iThread < nThreads; iThread++) {
118 random_number.emplace_back(std::make_unique<TRandom3>(0));
119 }
120 // Not using adaptive by default
121 use_adaptive = false;
122 // Set the covariance matrix
123 _fNumPar = CovMat->GetNrows();
124
125 InvertCovMatrix.resize(_fNumPar, std::vector<double>(_fNumPar, 0.0));
126 throwMatrixCholDecomp = new double*[_fNumPar]();
127 // Set the defaults to true
128 for(int i = 0; i < _fNumPar; i++) {
129 throwMatrixCholDecomp[i] = new double[_fNumPar]();
130 for (int j = 0; j < _fNumPar; j++) {
131 throwMatrixCholDecomp[i][j] = 0.;
132 }
133 }
134 SetName(name);
135 MakePosDef(CovMat);
136 SetCovMatrix(CovMat);
137 if (_fNumPar <= 0) {
138 MACH3LOG_CRITICAL("Covariance matrix {} has {} entries!", GetName(), _fNumPar);
139 throw MaCh3Exception(__FILE__ , __LINE__ );
140 }
141
143
144 infile->Close();
145
146 MACH3LOG_INFO("Created covariance matrix named: {}", GetName());
147 MACH3LOG_INFO("from file: {}", file);
148 delete infile;
149}
150
151// ********************************************
152// ETA An init function for the YAML constructor
153// All you really need from the YAML file is the number of Systematics
154void ParameterHandlerBase::Init(const std::vector<std::string>& YAMLFile) {
155// ********************************************
156 _fYAMLDoc["Systematics"] = YAML::Node(YAML::NodeType::Sequence);
157 for(unsigned int i = 0; i < YAMLFile.size(); i++)
158 {
159 YAML::Node YAMLDocTemp = M3OpenConfig(YAMLFile[i]);
160 for (const auto& item : YAMLDocTemp["Systematics"]) {
161 _fYAMLDoc["Systematics"].push_back(item);
162 }
163 }
164
165 const int nThreads = M3::GetNThreads();
166 //KS: set Random numbers for each thread so each thread has different seed
167 //or for one thread if without MULTITHREAD
168 random_number.reserve(nThreads);
169 for (int iThread = 0; iThread < nThreads; iThread++) {
170 random_number.emplace_back(std::make_unique<TRandom3>(0));
171 }
172 PrintLength = 35;
173
174 // Set the covariance matrix
175 _fNumPar = int(_fYAMLDoc["Systematics"].size());
176
177 use_adaptive = false;
178
179 InvertCovMatrix.resize(_fNumPar, std::vector<double>(_fNumPar, 0.0));
180 throwMatrixCholDecomp = new double*[_fNumPar]();
181 for(int i = 0; i < _fNumPar; i++) {
182 throwMatrixCholDecomp[i] = new double[_fNumPar]();
183 for (int j = 0; j < _fNumPar; j++) {
184 throwMatrixCholDecomp[i][j] = 0.;
185 }
186 }
188
189 TMatrixDSym* _fCovMatrix = new TMatrixDSym(_fNumPar);
190 int i = 0;
191 std::vector<std::map<std::string,double>> Correlations(_fNumPar);
192 std::map<std::string, int> CorrNamesMap;
193
194 //ETA - read in the systematics. Would be good to add in some checks to make sure
195 //that there are the correct number of entries i.e. are the _fNumPar for Names,
196 //PreFitValues etc etc.
197 for (auto const &param : _fYAMLDoc["Systematics"])
198 {
199 _fFancyNames[i] = Get<std::string>(param["Systematic"]["Names"]["FancyName"], __FILE__ , __LINE__);
200 _fPreFitValue[i] = Get<double>(param["Systematic"]["ParameterValues"]["PreFitValue"], __FILE__ , __LINE__);
201 _fIndivStepScale[i] = Get<double>(param["Systematic"]["StepScale"]["MCMC"], __FILE__ , __LINE__);
202 _fError[i] = Get<double>(param["Systematic"]["Error"], __FILE__ , __LINE__);
203 _fSampleNames[i] = GetFromManager<std::vector<std::string>>(param["Systematic"]["SampleNames"], {}, __FILE__, __LINE__);
204 if(_fError[i] <= 0) {
205 MACH3LOG_ERROR("Error for param {}({}) is negative and equal to {}", _fFancyNames[i], i, _fError[i]);
206 throw MaCh3Exception(__FILE__ , __LINE__ );
207 }
208 //ETA - a bit of a fudge but works
209 auto TempBoundsVec = GetBounds(param["Systematic"]["ParameterBounds"]);
210 _fLowBound[i] = TempBoundsVec[0];
211 _fUpBound[i] = TempBoundsVec[1];
212
213 //ETA - now for parameters which are optional and have default values
214 _fFlatPrior[i] = GetFromManager<bool>(param["Systematic"]["FlatPrior"], false, __FILE__ , __LINE__);
215
216 // Allow to fix param, this setting should be used only for params which are permanently fixed like baseline, please use global config for fixing param more flexibly
217 if(GetFromManager<bool>(param["Systematic"]["FixParam"], false, __FILE__ , __LINE__)) {
219 }
220
221 if(param["Systematic"]["SpecialProposal"]) {
222 EnableSpecialProposal(param["Systematic"]["SpecialProposal"], i);
223 }
224
225 //Fill the map to get the correlations later as well
226 CorrNamesMap[param["Systematic"]["Names"]["FancyName"].as<std::string>()]=i;
227
228 //Also loop through the correlations
229 if(param["Systematic"]["Correlations"]) {
230 for(unsigned int Corr_i = 0; Corr_i < param["Systematic"]["Correlations"].size(); ++Corr_i){
231 for (YAML::const_iterator it = param["Systematic"]["Correlations"][Corr_i].begin(); it!=param["Systematic"]["Correlations"][Corr_i].end();++it) {
232 Correlations[i][it->first.as<std::string>()] = it->second.as<double>();
233 }
234 }
235 }
236 i++;
237 } // end loop over para
238 if(i != _fNumPar) {
239 MACH3LOG_CRITICAL("Inconsistent number of params in Yaml {} vs {}, this indicate wrong syntax", i, i, _fNumPar);
240 throw MaCh3Exception(__FILE__ , __LINE__ );
241 }
242 // ETA Now that we've been through all systematic let's fill the covmatrix
243 //This makes the root TCov from YAML
244 for(int j = 0; j < _fNumPar; j++) {
245 (*_fCovMatrix)(j, j) = _fError[j]*_fError[j];
246 //Get the map of parameter name to correlation from the Correlations object
247 for (auto const& pair : Correlations[j]) {
248 auto const& key = pair.first;
249 auto const& val = pair.second;
250 int index = -1;
251 //If you found the parameter name then get the index
252 if (CorrNamesMap.find(key) != CorrNamesMap.end()) {
253 index = CorrNamesMap[key];
254 } else {
255 MACH3LOG_ERROR("Parameter {} not in list! Check your spelling?", key);
256 throw MaCh3Exception(__FILE__ , __LINE__ );
257 }
258 double Corr1 = val;
259 double Corr2 = 0;
260 if(Correlations[index].find(_fFancyNames[j]) != Correlations[index].end()) {
261 Corr2 = Correlations[index][_fFancyNames[j]];
262 //Do they agree to better than float precision?
263 if(std::abs(Corr2 - Corr1) > FLT_EPSILON) {
264 MACH3LOG_ERROR("Correlations are not equal between {} and {}", _fFancyNames[j], key);
265 MACH3LOG_ERROR("Got : {} and {}", Corr2, Corr1);
266 throw MaCh3Exception(__FILE__ , __LINE__ );
267 }
268 } else {
269 MACH3LOG_ERROR("Correlation does not appear reciprocally between {} and {}", _fFancyNames[j], key);
270 throw MaCh3Exception(__FILE__ , __LINE__ );
271 }
272 (*_fCovMatrix)(j, index)= (*_fCovMatrix)(index, j) = Corr1*_fError[j]*_fError[index];
273 }
274 }
275
276 //Now make positive definite
277 MakePosDef(_fCovMatrix);
278 SetCovMatrix(_fCovMatrix);
279
280 if (_fNumPar <= 0) {
281 MACH3LOG_ERROR("ParameterHandler object has {} systematics!", _fNumPar);
282 throw MaCh3Exception(__FILE__ , __LINE__ );
283 }
284
285 Tunes = std::make_unique<ParameterTunes>(_fYAMLDoc["Systematics"]);
286
287 MACH3LOG_INFO("Created covariance matrix from files: ");
288 for(const auto &file : YAMLFile){
289 MACH3LOG_INFO("{} ", file);
290 }
291 MACH3LOG_INFO("----------------");
292 MACH3LOG_INFO("Found {} systematics parameters in total", _fNumPar);
293 MACH3LOG_INFO("----------------");
294}
295
296// ********************************************
297void ParameterHandlerBase::EnableSpecialProposal(const YAML::Node& param, const int Index){
298// ********************************************
300
301 bool CircEnabled = false;
302 bool FlipEnabled = false;
303
304 if (param["CircularBounds"]) {
305 CircEnabled = true;
306 }
307
308 if (param["FlipParameter"]) {
309 FlipEnabled = true;
310 }
311
312 if (!CircEnabled && !FlipEnabled) {
313 MACH3LOG_ERROR("None of Special Proposal were enabled even though param {}, has SpecialProposal entry in Yaml", GetParFancyName(Index));
314 throw MaCh3Exception(__FILE__, __LINE__);
315 }
316
317 if (CircEnabled) {
318 CircularBoundsIndex.push_back(Index);
319 CircularBoundsValues.push_back(Get<std::pair<double, double>>(param["CircularBounds"], __FILE__, __LINE__));
320 MACH3LOG_INFO("Enabling CircularBounds for parameter {} with range [{}, {}]",
321 GetParFancyName(Index),
322 CircularBoundsValues.back().first,
323 CircularBoundsValues.back().second);
324 // KS: Make sure circular bounds are within physical bounds. If we are outside of physics bound MCMC will never explore such phase space region
325 if (CircularBoundsValues.back().first < _fLowBound.at(Index) || CircularBoundsValues.back().second > _fUpBound.at(Index)) {
326 MACH3LOG_ERROR("Circular bounds [{}, {}] for parameter {} exceed physical bounds [{}, {}]",
327 CircularBoundsValues.back().first, CircularBoundsValues.back().second,
328 GetParFancyName(Index),
329 _fLowBound.at(Index), _fUpBound.at(Index));
330 throw MaCh3Exception(__FILE__, __LINE__);
331 }
332 }
333
334 if (FlipEnabled) {
335 FlipParameterIndex.push_back(Index);
336 FlipParameterPoint.push_back(Get<double>(param["FlipParameter"], __FILE__, __LINE__));
337 MACH3LOG_INFO("Enabling Flipping for parameter {} with value {}",
338 GetParFancyName(Index),
339 FlipParameterPoint.back());
340 }
341
342 if (CircEnabled && FlipEnabled) {
343 if (FlipParameterPoint.back() < CircularBoundsValues.back().first || FlipParameterPoint.back() > CircularBoundsValues.back().second) {
344 MACH3LOG_ERROR("FlipParameter value {} for parameter {} is outside the CircularBounds [{}, {}]",
345 FlipParameterPoint.back(), GetParFancyName(Index), CircularBoundsValues.back().first, CircularBoundsValues.back().second);
346 throw MaCh3Exception(__FILE__, __LINE__);
347 }
348
349 const double low = CircularBoundsValues.back().first;
350 const double high = CircularBoundsValues.back().second;
351
352 // Sanity check: ensure flipping any x in [low, high] keeps the result in [low, high]
353 const double flipped_low = 2 * FlipParameterPoint.back() - low;
354 const double flipped_high = 2 * FlipParameterPoint.back() - high;
355 const double min_flip = std::min(flipped_low, flipped_high);
356 const double max_flip = std::max(flipped_low, flipped_high);
357
358 if (min_flip < low || max_flip > high) {
359 MACH3LOG_ERROR("Flipping about point {} for parameter {} would leave circular bounds [{}, {}]",
360 FlipParameterPoint.back(), GetParFancyName(Index), low, high);
361 throw MaCh3Exception(__FILE__, __LINE__);
362 }
363 }
364}
365
366// ********************************************
367// Set the covariance matrix for this class
369// ********************************************
370 if (cov == nullptr) {
371 MACH3LOG_ERROR("Could not find covariance matrix you provided to {}", __func__ );
372 throw MaCh3Exception(__FILE__ , __LINE__ );
373 }
374 covMatrix = cov;
375
376 invCovMatrix = static_cast<TMatrixDSym *>(cov->Clone());
377 invCovMatrix->Invert();
378 //KS: ROOT has bad memory management, using standard double means we can decrease most operation by factor 2 simply due to cache hits
379 for (int i = 0; i < _fNumPar; i++)
380 {
381 for (int j = 0; j < _fNumPar; ++j)
382 {
383 InvertCovMatrix[i][j] = (*invCovMatrix)(i,j);
384 }
385 }
386
387 SetThrowMatrix(cov);
388}
389// ********************************************
390void ParameterHandlerBase::ReserveMemory(const int SizeVec) {
391// ********************************************
392 _fNames = std::vector<std::string>(SizeVec);
393 _fFancyNames = std::vector<std::string>(SizeVec);
394 _fPreFitValue = std::vector<double>(SizeVec);
395 _fError = std::vector<double>(SizeVec);
396 _fCurrVal = std::vector<double>(SizeVec);
397 _fPropVal = std::vector<double>(SizeVec);
398 _fLowBound = std::vector<double>(SizeVec);
399 _fUpBound = std::vector<double>(SizeVec);
400 _fFlatPrior = std::vector<bool>(SizeVec);
401 _fIndivStepScale = std::vector<double>(SizeVec);
402 _fSampleNames = std::vector<std::vector<std::string>>(_fNumPar);
403
404 corr_throw = new double[SizeVec];
405 // set random parameter vector (for correlated steps)
406 randParams = new double[SizeVec];
407
408 // Set the defaults to true
409 for(int i = 0; i < SizeVec; i++) {
410 _fPreFitValue.at(i) = 1.;
411 _fError.at(i) = 1.;
412 _fCurrVal.at(i) = 0.;
413 _fPropVal.at(i) = 0.;
414 _fLowBound.at(i) = -999.99;
415 _fUpBound.at(i) = 999.99;
416 _fFlatPrior.at(i) = false;
417 _fIndivStepScale.at(i) = 1.;
418 corr_throw[i] = 0.0;
419 randParams[i] = 0.0;
420 }
421
422 _fGlobalStepScale = 1.0;
423}
424
425// ********************************************
426// Set all the covariance matrix parameters to a user-defined value
427// Might want to split this
428void ParameterHandlerBase::SetPar(int i , double val) {
429// ********************************************
430 MACH3LOG_INFO("Over-riding {}: ", GetParName(i));
431 MACH3LOG_INFO("_fPropVal ({}), _fCurrVal ({}), _fPreFitValue ({}) to ({})", _fPropVal[i], _fCurrVal[i], _fPreFitValue[i], val);
432
433 _fPropVal[i] = val;
434 _fCurrVal[i] = val;
435 _fPreFitValue[i] = val;
436
437 // Transfer the parameter values to the PCA basis
438 if (pca) PCAObj->TransferToPCA();
439}
440
441// ********************************************
442std::vector<double> ParameterHandlerBase::GetProposed() const {
443// ********************************************
444 std::vector<double> props(_fNumPar);
445 for (int i = 0; i < _fNumPar; ++i) props[i] = _fPropVal[i];
446 return props;
447}
448
449// *************************************
450// Throw the parameters according to the covariance matrix
451// This shouldn't be used in MCMC code ase it can break Detailed Balance;
453// *************************************
454 // First draw new randParams
455 Randomize();
456
458
459 // KS: We use PCA very rarely on top PCA functionality isn't implemented for this function.
460 // Use __builtin_expect to give compiler a hint which option is more likely, which should help
461 // with better optimisation. This isn't critical but more to have example
462 if (__builtin_expect(!pca, 1)) {
463 #ifdef MULTITHREAD
464 #pragma omp parallel for
465 #endif
466 for (int i = 0; i < _fNumPar; ++i) {
467 // Check if parameter is fixed first: if so don't randomly throw
468 if (IsParameterFixed(i)) continue;
469
471 int throws = 0;
472 // Try again if we the initial parameter proposal falls outside of the range of the parameter
473 while (_fPropVal[i] > _fUpBound[i] || _fPropVal[i] < _fLowBound[i]) {
474 randParams[i] = random_number[M3::GetThreadIndex()]->Gaus(0, 1);
475 const double corr_throw_single = M3::MatrixVectorMultiSingle(throwMatrixCholDecomp, randParams, _fNumPar, i);
476 _fPropVal[i] = _fPreFitValue[i] + corr_throw_single;
477 if (throws > 10000)
478 {
479 //KS: Since we are multithreading there is danger that those messages
480 //will be all over the place, small price to pay for faster code
481 MACH3LOG_WARN("Tried {} times to throw parameter {} but failed", throws, i);
482 MACH3LOG_WARN("Matrix: {}", matrixName);
483 MACH3LOG_WARN("Param: {}", _fNames[i]);
484 MACH3LOG_WARN("Setting _fPropVal: {} to {}", _fPropVal[i], _fPreFitValue[i]);
485 MACH3LOG_WARN("I live at {}:{}", __FILE__, __LINE__);
486 _fPropVal[i] = _fPreFitValue[i];
487 //throw MaCh3Exception(__FILE__ , __LINE__ );
488 }
489 throws++;
490 }
491 _fCurrVal[i] = _fPropVal[i];
492 }
493 }
494 else
495 {
499 } // end if pca
500}
501
502// *************************************
503// Throw each parameter within their 1 sigma range
504// Used to start the chain in different states
506// *************************************
507 // Have the 1 sigma for each parameter in each covariance class, sweet!
508 // Don't want to change the prior array because that's what determines our likelihood
509 // Want to change the _fPropVal, _fCurrVal, _fPreFitValue
510 // _fPreFitValue and the others will already be set
511 for (int i = 0; i < _fNumPar; ++i) {
512 // Check if parameter is fixed first: if so don't randomly throw
513 if (IsParameterFixed(i)) continue;
514 // Check that the sigma range is larger than the parameter range
515 // If not, throw in the valid parameter range instead
516 const double paramrange = _fUpBound[i] - _fLowBound[i];
517 const double sigma = sqrt((*covMatrix)(i,i));
518 double throwrange = sigma;
519 if (paramrange < sigma) throwrange = paramrange;
520
521 _fPropVal[i] = _fPreFitValue[i] + random_number[0]->Gaus(0, 1)*throwrange;
522 // Try again if we the initial parameter proposal falls outside of the range of the parameter
523 int throws = 0;
524 while (_fPropVal[i] > _fUpBound[i] || _fPropVal[i] < _fLowBound[i]) {
525 if (throws > 1000) {
526 MACH3LOG_WARN("Tried {} times to throw parameter {} but failed", throws, i);
527 MACH3LOG_WARN("Matrix: {}", matrixName);
528 MACH3LOG_WARN("Param: {}", _fNames[i]);
529 throw MaCh3Exception(__FILE__ , __LINE__ );
530 }
531 _fPropVal[i] = _fPreFitValue[i] + random_number[0]->Gaus(0, 1)*throwrange;
532 throws++;
533 }
534 MACH3LOG_INFO("Setting current step in {} param {} = {} from {}", matrixName, i, _fPropVal[i], _fCurrVal[i]);
535 _fCurrVal[i] = _fPropVal[i];
536 }
537 if (pca) PCAObj->TransferToPCA();
538}
539
540// *************************************
541// Set a single parameter
542void ParameterHandlerBase::SetSingleParameter(const int parNo, const double parVal) {
543// *************************************
544 _fPropVal[parNo] = parVal;
545 _fCurrVal[parNo] = parVal;
546 MACH3LOG_DEBUG("Setting {} (parameter {}) to {})", GetParName(parNo), parNo, parVal);
547 if (pca) PCAObj->TransferToPCA();
548}
549
550// ********************************************
551void ParameterHandlerBase::SetParCurrProp(const int parNo, const double parVal) {
552// ********************************************
553 _fPropVal[parNo] = parVal;
554 _fCurrVal[parNo] = parVal;
555 MACH3LOG_DEBUG("Setting {} (parameter {}) to {})", GetParName(parNo), parNo, parVal);
556 if (pca) PCAObj->TransferToPCA();
557}
558
559// ************************************************
560// Propose a step for the set of systematics parameters this covariance class holds
562// ************************************************
563 // Make the random numbers for the step proposal
564 Randomize();
567
568 // KS: According to Dr Wallace we update using previous not proposed step
569 // this way we do special proposal after adaptive after.
570 // This way we can shortcut and skip rest of proposal
571 if(!doSpecialStepProposal) return;
572
574}
575
576// ************************************************
578// ************************************************
580
581 // HW It should now automatically set dcp to be with [-pi, pi]
582 for (size_t i = 0; i < CircularBoundsIndex.size(); ++i) {
584 }
585
586 // Okay now we've done the standard steps, we can add in our nice flips hierarchy flip first
587 for (size_t i = 0; i < FlipParameterIndex.size(); ++i) {
589 }
590}
591
592// ************************************************
593// "Randomize" the parameters in the covariance class for the proposed step
594// Used the proposal kernel and the current parameter value to set proposed step
595// Also get a new random number for the randParams
597// ************************************************
598 if (!pca) {
599 //KS: By multithreading here we gain at least factor 2 with 8 threads with ND only fit
600 #ifdef MULTITHREAD
601 #pragma omp parallel for
602 #endif
603 for (int i = 0; i < _fNumPar; ++i) {
604 // If parameter isn't fixed
605 if (!IsParameterFixed(i) > 0.0) {
606 randParams[i] = random_number[M3::GetThreadIndex()]->Gaus(0, 1);
607 // If parameter IS fixed
608 } else {
609 randParams[i] = 0.0;
610 }
611 } // end for
612 // If we're in the PCA basis we instead throw parameters there (only _fNumParPCA parameter)
613 } else {
614 // Scale the random parameters by the sqrt of eigen values for the throw
615 #ifdef MULTITHREAD
616 #pragma omp parallel for
617 #endif
618 for (int i = 0; i < PCAObj->GetNumberPCAedParameters(); ++i)
619 {
620 // If parameter IS fixed or out of bounds
621 if (PCAObj->IsParameterFixedPCA(i)) {
622 randParams[i] = 0.0;
623 } else {
625 }
626 }
627 }
628}
629
630// ************************************************
631// Correlate the steps by setting the proposed step of a parameter to its current value + some correlated throw
633// ************************************************
634 //KS: Using custom function compared to ROOT one with 8 threads we have almost factor 2 performance increase, by replacing TMatrix with just double we increase it even more
636
637 // If not doing PCA
638 if (!pca) {
639 #ifdef MULTITHREAD
640 #pragma omp parallel for
641 #endif
642 for (int i = 0; i < _fNumPar; ++i) {
643 if (!IsParameterFixed(i) > 0.) {
645 }
646 }
647 // If doing PCA throw uncorrelated in PCA basis (orthogonal basis by definition)
648 } else {
650 }
651}
652// ********************************************
653// Update so that current step becomes the previously proposed step
655// ********************************************
656 if (!pca) {
657 #ifdef MULTITHREAD
658 #pragma omp parallel for
659 #endif
660 for (int i = 0; i < _fNumPar; ++i) {
661 // Update state so that current state is proposed state
662 _fCurrVal[i] = _fPropVal[i];
663 }
664 } else {
665 PCAObj->AcceptStep();
666 }
667}
668
669// *************************************
670//HW: This method is a tad hacky but modular arithmetic gives me a headache.
671void ParameterHandlerBase::CircularParBounds(const int index, const double LowBound, const double UpBound) {
672// *************************************
673 if(_fPropVal[index] > UpBound) {
674 _fPropVal[index] = LowBound + std::fmod(_fPropVal[index] - UpBound, UpBound - LowBound);
675 } else if (_fPropVal[index] < LowBound) {
676 _fPropVal[index] = UpBound - std::fmod(LowBound - _fPropVal[index], UpBound - LowBound);
677 }
678}
679
680// *************************************
681void ParameterHandlerBase::FlipParameterValue(const int index, const double FlipPoint) {
682// *************************************
683 if(random_number[0]->Uniform() < 0.5) {
684 _fPropVal[index] = 2 * FlipPoint - _fPropVal[index];
685 }
686}
687
688// ********************************************
689// Throw the proposed parameter by mag sigma
690// Should really just have the user specify this throw by having argument double
692// ********************************************
693 Randomize();
694 if (!pca) {
695 // Make the correlated throw
697 // Number of sigmas we throw
698 for (int i = 0; i < _fNumPar; i++) {
699 if (!IsParameterFixed(i) > 0.)
700 _fPropVal[i] = _fCurrVal[i] + corr_throw[i]*mag;
701 }
702 } else {
703 PCAObj->ThrowParProp(mag, randParams);
704 }
705}
706// ********************************************
707// Helper function to throw the current parameter by mag sigmas
708// Can study bias in MCMC with this; put different starting parameters
710// ********************************************
711 Randomize();
712 if (!pca) {
713 // Get the correlated throw vector
715 // The number of sigmas to throw
716 // Should probably have this as a default parameter input to the function instead
717 for (int i = 0; i < _fNumPar; i++) {
718 if (!IsParameterFixed(i) > 0.){
719 _fCurrVal[i] = corr_throw[i]*mag;
720 }
721 }
722 } else {
723 PCAObj->ThrowParCurr(mag, randParams);
724 }
725}
726// ********************************************
727// Function to print the prior values
729// ********************************************
730 MACH3LOG_INFO("Prior values for {} ParameterHandler:", GetName());
731 for (int i = 0; i < _fNumPar; i++) {
732 MACH3LOG_INFO(" {} {} ", GetParFancyName(i), GetParInit(i));
733 }
734}
735
736// ********************************************
737// Function to print the prior, current and proposed values
739// ********************************************
740 MACH3LOG_INFO("Printing parameters for {}", GetName());
741 // Dump out the PCA parameters too
742 if (pca) {
743 PCAObj->Print();
744 }
745 MACH3LOG_INFO("{:<30} {:<10} {:<10} {:<10}", "Name", "Prior", "Current", "Proposed");
746 for (int i = 0; i < _fNumPar; ++i) {
747 MACH3LOG_INFO("{:<30} {:<10.2f} {:<10.2f} {:<10.2f}", GetParFancyName(i), _fPreFitValue[i], _fCurrVal[i], _fPropVal[i]);
748 }
749}
750
751// ********************************************
752// Get the likelihood in the case where we want to include priors on the parameters
753// fParEvalLikelihood stores if we want to evaluate the likelihood for the given parameter
754// true = evaluate likelihood (so run with a prior)
755// false = don't evaluate likelihood (so run without a prior)
757// ********************************************
758 double logL = 0.0;
759 #ifdef MULTITHREAD
760 #pragma omp parallel for reduction(+:logL)
761 #endif
762 for(int i = 0; i < _fNumPar; ++i){
763 if(_fFlatPrior[i]){
764 //HW: Flat prior, no need to calculate anything
765 continue;
766 }
767
768 #ifdef MULTITHREAD
769 #pragma omp simd
770 #endif
771 for (int j = 0; j <= i; ++j) {
772 if (!_fFlatPrior[j]) {
773 //KS: Since matrix is symmetric we can calculate non diagonal elements only once and multiply by 2, can bring up to factor speed decrease.
774 double scale = (i != j) ? 1. : 0.5;
775 logL += scale * (_fPropVal[i] - _fPreFitValue[i])*(_fPropVal[j] - _fPreFitValue[j])*InvertCovMatrix[i][j];
776 }
777 }
778 }
779 return logL;
780}
781
782// ********************************************
784// ********************************************
785 int NOutside = 0;
786 #ifdef MULTITHREAD
787 #pragma omp parallel for reduction(+:NOutside)
788 #endif
789 for (int i = 0; i < _fNumPar; ++i){
790 if(_fPropVal[i] > _fUpBound[i] || _fPropVal[i] < _fLowBound[i]){
791 NOutside++;
792 }
793 }
794 return NOutside;
795}
796
797// ********************************************
799// ********************************************
800 // Default behaviour is to reject negative values + do std llh calculation
801 const int NOutside = CheckBounds();
802
803 if(NOutside > 0) return NOutside*M3::_LARGE_LOGL_;
804
805 return CalcLikelihood();
806}
807
808// ********************************************
809// Sets the proposed parameters to the prior values
810void ParameterHandlerBase::SetParameters(const std::vector<double>& pars) {
811// ********************************************
812 // If empty, set the proposed to prior
813 if (pars.empty()) {
814 // For xsec this means setting to the prior (because prior is the prior)
815 for (int i = 0; i < _fNumPar; i++) {
816 _fPropVal[i] = _fPreFitValue[i];
817 }
818 // If not empty, set the parameters to the specified
819 } else {
820 if (pars.size() != size_t(_fNumPar)) {
821 MACH3LOG_ERROR("Parameter arrays of incompatible size! Not changing parameters! {} has size {} but was expecting {}", matrixName, pars.size(), _fNumPar);
822 throw MaCh3Exception(__FILE__ , __LINE__ );
823 }
824 int parsSize = int(pars.size());
825 for (int i = 0; i < parsSize; i++) {
826 //Make sure that you are actually passing a number to set the parameter to
827 if(std::isnan(pars[i])) {
828 MACH3LOG_ERROR("Trying to set parameter value to a nan for parameter {} in matrix {}. This will not go well!", GetParName(i), matrixName);
829 throw MaCh3Exception(__FILE__ , __LINE__ );
830 } else {
831 _fPropVal[i] = pars[i];
832 }
833 }
834 }
835 // And if pca make the transfer
836 if (pca) {
837 PCAObj->TransferToPCA();
838 PCAObj->TransferToParam();
839 }
840}
841
842// ********************************************
843void ParameterHandlerBase::SetBranches(TTree &tree, bool SaveProposal) {
844// ********************************************
845 // loop over parameters and set a branch
846 for (int i = 0; i < _fNumPar; ++i) {
847 tree.Branch(_fNames[i].c_str(), &_fCurrVal[i], Form("%s/D", _fNames[i].c_str()));
848 }
849 // When running PCA, also save PCA parameters
850 if (pca) {
851 PCAObj->SetBranches(tree, SaveProposal, _fNames);
852 }
853 if(SaveProposal)
854 {
855 // loop over parameters and set a branch
856 for (int i = 0; i < _fNumPar; ++i) {
857 tree.Branch(Form("%s_Prop", _fNames[i].c_str()), &_fPropVal[i], Form("%s_Prop/D", _fNames[i].c_str()));
858 }
859 }
860}
861
862// ********************************************
863void ParameterHandlerBase::SetStepScale(const double scale) {
864// ********************************************
865 if(scale <= 0) {
866 MACH3LOG_ERROR("You are trying so set StepScale to 0 or negative this will not work");
867 throw MaCh3Exception(__FILE__ , __LINE__ );
868 }
869 MACH3LOG_INFO("{} setStepScale() = {}", GetName(), scale);
870 const double SuggestedScale = 2.38*2.38/_fNumPar;
871 if(std::fabs(scale - SuggestedScale)/SuggestedScale > 1) {
872 MACH3LOG_WARN("Defined Global StepScale is {}, while suggested suggested {}", scale, SuggestedScale);
873 }
874
875 _fGlobalStepScale = scale;
876}
877
878// ********************************************
879int ParameterHandlerBase::GetParIndex(const std::string& name) const {
880// ********************************************
881 int Index = M3::_BAD_INT_;
882 for (int i = 0; i <_fNumPar; ++i) {
883 if(name == _fFancyNames[i]) {
884 Index = i;
885 break;
886 }
887 }
888 return Index;
889}
890
891// ********************************************
893// ********************************************
894 // fix or unfix all parameters by multiplying by -1
895 if(!pca) {
896 for (int i = 0; i < _fNumPar; i++) _fError[i] *= -1.0;
897 } else{
898 PCAObj->ToggleFixAllParameters();
899 }
900}
901
902// ********************************************
904// ********************************************
905 if(!pca) {
906 if (i > _fNumPar) {
907 MACH3LOG_ERROR("Can't {} for parameter {} because size of covariance ={}", __func__, i, _fNumPar);
908 MACH3LOG_ERROR("Fix this in your config file please!");
909 throw MaCh3Exception(__FILE__ , __LINE__ );
910 } else {
911 _fError[i] *= -1.0;
912 MACH3LOG_INFO("Setting {}(parameter {}) to fixed at {}", GetParFancyName(i), i, _fCurrVal[i]);
913 }
914 } else {
915 PCAObj->ToggleFixParameter(i, _fNames);
916 }
917}
918
919// ********************************************
920void ParameterHandlerBase::ToggleFixParameter(const std::string& name) {
921// ********************************************
922 const int Index = GetParIndex(name);
923 if(Index != M3::_BAD_INT_) {
924 ToggleFixParameter(Index);
925 return;
926 }
927
928 MACH3LOG_WARN("I couldn't find parameter with name {}, therefore will not fix it", name);
929}
930
931// ********************************************
932bool ParameterHandlerBase::IsParameterFixed(const std::string& name) const {
933// ********************************************
934 const int Index = GetParIndex(name);
935 if(Index != M3::_BAD_INT_) {
936 return IsParameterFixed(Index);
937 }
938
939 MACH3LOG_WARN("I couldn't find parameter with name {}, therefore don't know if it fixed", name);
940 return false;
941}
942
943// ********************************************
944void ParameterHandlerBase::SetFlatPrior(const int i, const bool eL) {
945// ********************************************
946 if (i > _fNumPar) {
947 MACH3LOG_INFO("Can't {} for Cov={}/Param={} because size of Covariance = {}", __func__, GetName(), i, _fNumPar);
948 MACH3LOG_ERROR("Fix this in your config file please!");
949 throw MaCh3Exception(__FILE__ , __LINE__ );
950 } else {
951 if(eL){
952 MACH3LOG_INFO("Setting {} (parameter {}) to flat prior", GetParName(i), i);
953 }
954 else{
955 // HW :: This is useful
956 MACH3LOG_INFO("Setting {} (parameter {}) to non-flat prior", GetParName(i), i);
957 }
958 _fFlatPrior[i] = eL;
959 }
960}
961
962// ********************************************
963void ParameterHandlerBase::SetIndivStepScale(const std::vector<double>& stepscale) {
964// ********************************************
965 if (int(stepscale.size()) != _fNumPar)
966 {
967 MACH3LOG_WARN("Stepscale vector not equal to number of parameters. Quitting..");
968 MACH3LOG_WARN("Size of argument vector: {}", stepscale.size());
969 MACH3LOG_WARN("Expected size: {}", _fNumPar);
970 return;
971 }
972
973 for (int iParam = 0 ; iParam < _fNumPar; iParam++) {
974 _fIndivStepScale[iParam] = stepscale[iParam];
975 }
977}
978
979// ********************************************
981// ********************************************
982 MACH3LOG_INFO("============================================================");
983 MACH3LOG_INFO("{:<{}} | {:<11}", "Parameter:", PrintLength, "Step scale:");
984 for (int iParam = 0; iParam < _fNumPar; iParam++) {
985 MACH3LOG_INFO("{:<{}} | {:<11}", _fFancyNames[iParam].c_str(), PrintLength, _fIndivStepScale[iParam]);
986 }
987 MACH3LOG_INFO("============================================================");
988}
989
990// ********************************************
991//Makes sure that matrix is positive-definite by adding a small number to on-diagonal elements
992void ParameterHandlerBase::MakePosDef(TMatrixDSym *cov) {
993// ********************************************
994 if(cov == nullptr){
995 cov = &*covMatrix;
996 MACH3LOG_WARN("Passed nullptr to cov matrix in {}", matrixName);
997 }
998
999 //DB Save original warning state and then increase it in this function to suppress 'matrix not positive definite' messages
1000 //Means we no longer need to overload
1001 int originalErrorWarning = gErrorIgnoreLevel;
1002 gErrorIgnoreLevel = kFatal;
1003
1004 //DB Loop 1000 times adding 1e-9 which tops out at 1e-6 shift on the diagonal before throwing error
1005 int MaxAttempts = 1e5;
1006 int iAttempt = 0;
1007 bool CanDecomp = false;
1008 TDecompChol chdcmp;
1009
1010 for (iAttempt = 0; iAttempt < MaxAttempts; iAttempt++) {
1011 chdcmp = TDecompChol(*cov);
1012 if (chdcmp.Decompose()) {
1013 CanDecomp = true;
1014 break;
1015 } else {
1016 #ifdef MULTITHREAD
1017 #pragma omp parallel for
1018 #endif
1019 for (int iVar = 0 ; iVar < _fNumPar; iVar++) {
1020 (*cov)(iVar,iVar) += pow(10, -9);
1021 }
1022 }
1023 }
1024
1025 if (!CanDecomp) {
1026 MACH3LOG_ERROR("Tried {} times to shift diagonal but still can not decompose the matrix", MaxAttempts);
1027 MACH3LOG_ERROR("This indicates that something is wrong with the input matrix");
1028 throw MaCh3Exception(__FILE__ , __LINE__ );
1029 }
1030 if(!use_adaptive || AdaptiveHandler->GetTotalSteps() < 2) {
1031 MACH3LOG_INFO("Had to shift diagonal {} time(s) to allow the covariance matrix to be decomposed", iAttempt);
1032 }
1033 //DB Resetting warning level
1034 gErrorIgnoreLevel = originalErrorWarning;
1035}
1036
1037// ********************************************
1039// ********************************************
1040 std::vector<double> stepScales(_fNumPar);
1041 for (int i = 0; i <_fNumPar; i++) {
1042 stepScales[i] = 1.;
1043 }
1044 _fGlobalStepScale = 1.0;
1045 SetIndivStepScale(stepScales);
1046}
1047
1048// ********************************************
1049// HW: Code for throwing from separate throw matrix, needs to be set after init to ensure pos-def
1051// ********************************************
1052 if (cov == nullptr) {
1053 MACH3LOG_ERROR("Could not find covariance matrix you provided to {}", __func__);
1054 throw MaCh3Exception(__FILE__ , __LINE__ );
1055 }
1056
1057 if (covMatrix->GetNrows() != cov->GetNrows()) {
1058 MACH3LOG_ERROR("Matrix given for throw Matrix is not the same size as the covariance matrix stored in object!");
1059 MACH3LOG_ERROR("Stored covariance matrix size: {}", covMatrix->GetNrows());
1060 MACH3LOG_ERROR("Given matrix size: {}", cov->GetNrows());
1061 throw MaCh3Exception(__FILE__ , __LINE__ );
1062 }
1063
1064 throwMatrix = static_cast<TMatrixDSym*>(cov->Clone());
1066 else MakePosDef(throwMatrix);
1067
1068 TDecompChol TDecompChol_throwMatrix(*throwMatrix);
1069
1070 if(!TDecompChol_throwMatrix.Decompose()) {
1071 MACH3LOG_ERROR("Cholesky decomposition failed for {} trying to make positive definite", matrixName);
1072 throw MaCh3Exception(__FILE__ , __LINE__ );
1073 }
1074
1075 throwMatrix_CholDecomp = new TMatrixD(TDecompChol_throwMatrix.GetU());
1077
1078 //KS: ROOT has bad memory management, using standard double means we can decrease most operation by factor 2 simply due to cache hits
1079 #ifdef MULTITHREAD
1080 #pragma omp parallel for collapse(2)
1081 #endif
1082 for (int i = 0; i < _fNumPar; ++i)
1083 {
1084 for (int j = 0; j < _fNumPar; ++j)
1085 {
1086 throwMatrixCholDecomp[i][j] = (*throwMatrix_CholDecomp)(i,j);
1087 }
1088 }
1089}
1090
1091// ********************************************
1093// ********************************************
1094 delete throwMatrix;
1095 throwMatrix = nullptr;
1097 throwMatrix_CholDecomp = nullptr;
1098 SetThrowMatrix(cov);
1099}
1100
1101// ********************************************
1102// HW : Here be adaption
1103void ParameterHandlerBase::InitialiseAdaption(const YAML::Node& adapt_manager){
1104// ********************************************
1105 if(PCAObj){
1106 MACH3LOG_ERROR("PCA has been enabled and now trying to enable Adaption. Right now both configuration don't work with each other");
1107 throw MaCh3Exception(__FILE__ , __LINE__ );
1108 }
1109 AdaptiveHandler = std::make_unique<adaptive_mcmc::AdaptiveMCMCHandler>();
1110 // Now we read the general settings [these SHOULD be common across all matrices!]
1111 bool success = AdaptiveHandler->InitFromConfig(adapt_manager, matrixName, &_fCurrVal, &_fError);
1112 if(!success) return;
1113 AdaptiveHandler->Print();
1114
1115 // Next let"s check for external matrices
1116 // We"re going to grab this info from the YAML manager
1117 if(!GetFromManager<bool>(adapt_manager["AdaptionOptions"]["Covariance"][matrixName]["UseExternalMatrix"], false, __FILE__ , __LINE__)) {
1118 MACH3LOG_WARN("Not using external matrix for {}, initialising adaption from scratch", matrixName);
1119 // If we don't have a covariance matrix to start from for adaptive tune we need to make one!
1120 use_adaptive = true;
1121 AdaptiveHandler->CreateNewAdaptiveCovariance();
1122 return;
1123 }
1124
1125 // Finally, we accept that we want to read the matrix from a file!
1126 auto external_file_name = GetFromManager<std::string>(adapt_manager["AdaptionOptions"]["Covariance"][matrixName]["ExternalMatrixFileName"], "", __FILE__ , __LINE__);
1127 auto external_matrix_name = GetFromManager<std::string>(adapt_manager["AdaptionOptions"]["Covariance"][matrixName]["ExternalMatrixName"], "", __FILE__ , __LINE__);
1128 auto external_mean_name = GetFromManager<std::string>(adapt_manager["AdaptionOptions"]["Covariance"][matrixName]["ExternalMeansName"], "", __FILE__ , __LINE__);
1129
1130 AdaptiveHandler->SetThrowMatrixFromFile(external_file_name, external_matrix_name, external_mean_name, use_adaptive);
1131 SetThrowMatrix(AdaptiveHandler->GetAdaptiveCovariance());
1132
1133 MACH3LOG_INFO("Successfully Set External Throw Matrix Stored in {}", external_file_name);
1134}
1135
1136// ********************************************
1137// Truely adaptive MCMC!
1139// ********************************************
1140 // Updates adaptive matrix
1141 // First we update the total means
1142
1143 // Skip this if we're at a large number of steps
1144 if(AdaptiveHandler->SkipAdaption()) {
1145 AdaptiveHandler->IncrementNSteps();
1146 return;
1147 }
1148
1149 // Call main adaption function
1150 AdaptiveHandler->UpdateAdaptiveCovariance();
1151
1152 //This is likely going to be the slow bit!
1153 if(AdaptiveHandler->IndivStepScaleAdapt()) {
1155 }
1156
1157 if(AdaptiveHandler->UpdateMatrixAdapt()) {
1158 TMatrixDSym* update_matrix = static_cast<TMatrixDSym*>(AdaptiveHandler->GetAdaptiveCovariance()->Clone());
1159 UpdateThrowMatrix(update_matrix); //Now we update and continue!
1160 //Also Save the adaptive to file
1161 AdaptiveHandler->SaveAdaptiveToFile(AdaptiveHandler->GetOutFileName(), GetName());
1162 }
1163
1164 AdaptiveHandler->IncrementNSteps();
1165}
1166
1167// ********************************************
1168//HW: Finds closest possible positive definite matrix in Frobenius Norm ||.||_frob
1169// Where ||X||_frob=sqrt[sum_ij(x_ij^2)] (basically just turns an n,n matrix into vector in n^2 space
1170// then does Euclidean norm)
1172// ********************************************
1173 // Want to get cov' = (cov_sym+cov_polar)/2
1174 // cov_sym=(cov+cov^T)/2
1175 // cov_polar-> SVD cov to cov=USV^T then cov_polar=VSV^T
1176
1177 //Get frob norm of cov
1178 // Double_t cov_norm=cov->E2Norm();
1179
1180 TMatrixDSym* cov_trans = cov;
1181 cov_trans->T();
1182 TMatrixDSym cov_sym = 0.5*(*cov+*cov_trans); //If cov is symmetric does nothing, otherwise just ensures symmetry
1183
1184 //Do SVD to get polar form
1185 TDecompSVD cov_sym_svd=TDecompSVD(cov_sym);
1186 if(!cov_sym_svd.Decompose()){
1187 MACH3LOG_ERROR("Cannot do SVD on input matrix!");
1188 throw MaCh3Exception(__FILE__ , __LINE__ );
1189 }
1190
1191 TMatrixD cov_sym_v = cov_sym_svd.GetV();
1192 TMatrixD cov_sym_vt = cov_sym_v;
1193 cov_sym_vt.T();
1194 //SVD returns as vector (grrr) so need to get into matrix form for multiplying!
1195 TVectorD cov_sym_sigvect = cov_sym_svd.GetSig();
1196
1197 const Int_t nCols = cov_sym_v.GetNcols(); //square so only need rows hence lack of cols
1198 TMatrixDSym cov_sym_sig(nCols);
1199 TMatrixDDiag cov_sym_sig_diag(cov_sym_sig);
1200 cov_sym_sig_diag=cov_sym_sigvect;
1201
1202 //Can finally get H=VSV
1203 TMatrixDSym cov_sym_polar = cov_sym_sig.SimilarityT(cov_sym_vt);//V*S*V^T (this took forver to find!)
1204
1205 //Now we can construct closest approximater Ahat=0.5*(B+H)
1206 TMatrixDSym cov_closest_approx = 0.5*(cov_sym+cov_sym_polar);//Not fully sure why this is even needed since symmetric B -> U=V
1207 //Get norm of transformed
1208 // Double_t approx_norm=cov_closest_approx.E2Norm();
1209 //MACH3LOG_INFO("Initial Norm: {:.6f} | Norm after transformation: {:.6f} | Ratio: {:.6f}", cov_norm, approx_norm, cov_norm / approx_norm);
1210
1211 *cov = cov_closest_approx;
1212 //Now can just add a makeposdef!
1213 MakePosDef(cov);
1214}
1215
1216// ********************************************
1217// KS: Convert covariance matrix to correlation matrix and return TH2D which can be used for fancy plotting
1219// ********************************************
1220 TH2D* hMatrix = new TH2D(GetName().c_str(), GetName().c_str(), _fNumPar, 0.0, _fNumPar, _fNumPar, 0.0, _fNumPar);
1221
1222 for(int i = 0; i < _fNumPar; i++)
1223 {
1224 hMatrix->SetBinContent(i+1, i+1, 1.);
1225 hMatrix->GetXaxis()->SetBinLabel(i+1, GetParFancyName(i).c_str());
1226 hMatrix->GetYaxis()->SetBinLabel(i+1, GetParFancyName(i).c_str());
1227 }
1228
1229 #ifdef MULTITHREAD
1230 #pragma omp parallel for
1231 #endif
1232 for(int i = 0; i < _fNumPar; i++)
1233 {
1234 for(int j = 0; j <= i; j++)
1235 {
1236 const double Corr = (*covMatrix)(i,j) / ( GetDiagonalError(i) * GetDiagonalError(j));
1237 hMatrix->SetBinContent(i+1, j+1, Corr);
1238 hMatrix->SetBinContent(j+1, i+1, Corr);
1239 }
1240 }
1241 return hMatrix;
1242}
1243
1244// ********************************************
1245// KS: After step scale, prefit etc. value were modified save this modified config.
1247// ********************************************
1248 if (!_fYAMLDoc)
1249 {
1250 MACH3LOG_CRITICAL("Yaml node hasn't been initialised for matrix {}, something is not right", matrixName);
1251 MACH3LOG_CRITICAL("I am not throwing error but should be investigated");
1252 return;
1253 }
1254
1255 YAML::Node copyNode = _fYAMLDoc;
1256 int i = 0;
1257
1258 for (YAML::Node param : copyNode["Systematics"])
1259 {
1260 //KS: Feel free to update it, if you need updated prefit value etc
1261 param["Systematic"]["StepScale"]["MCMC"] = MaCh3Utils::FormatDouble(_fIndivStepScale[i], 4);
1262 i++;
1263 }
1264 // Save the modified node to a file
1265 std::ofstream fout("Modified_Matrix.yaml");
1266 fout << copyNode;
1267 fout.close();
1268}
1269
1270// ********************************************
1271bool ParameterHandlerBase::AppliesToSample(const int SystIndex, const std::string& SampleName) const {
1272// ********************************************
1273 // Empty means apply to all
1274 if (_fSampleNames[SystIndex].size() == 0) return true;
1275
1276 // Make a copy and to lower case to not be case sensitive
1277 std::string SampleNameCopy = SampleName;
1278 std::transform(SampleNameCopy.begin(), SampleNameCopy.end(), SampleNameCopy.begin(), ::tolower);
1279
1280 // Check for unsupported wildcards in SampleNameCopy
1281 if (SampleNameCopy.find('*') != std::string::npos) {
1282 MACH3LOG_ERROR("Wildcards ('*') are not supported in sample name: '{}'", SampleName);
1283 throw MaCh3Exception(__FILE__ , __LINE__ );
1284 }
1285
1286 bool Applies = false;
1287
1288 for (size_t i = 0; i < _fSampleNames[SystIndex].size(); i++) {
1289 // Convert to low case to not be case sensitive
1290 std::string pattern = _fSampleNames[SystIndex][i];
1291 std::transform(pattern.begin(), pattern.end(), pattern.begin(), ::tolower);
1292
1293 // Replace '*' in the pattern with '.*' for regex matching
1294 std::string regexPattern = "^" + std::regex_replace(pattern, std::regex("\\*"), ".*") + "$";
1295 try {
1296 std::regex regex(regexPattern);
1297 if (std::regex_match(SampleNameCopy, regex)) {
1298 Applies = true;
1299 break;
1300 }
1301 } catch (const std::regex_error& e) {
1302 // Handle regex error (for invalid patterns)
1303 MACH3LOG_ERROR("Regex error: {}", e.what());
1304 }
1305 }
1306 return Applies;
1307}
1308
1309// ********************************************
1310// Set proposed parameter values vector to be base on tune values
1311void ParameterHandlerBase::SetTune(const std::string& TuneName) {
1312// ********************************************
1313 if(Tunes == nullptr) {
1314 MACH3LOG_ERROR("Tunes haven't been initialised, which are being loaded from YAML, have you used some deprecated constructor");
1315 throw MaCh3Exception(__FILE__, __LINE__);
1316 }
1317 auto Values = Tunes->GetTune(TuneName);
1318
1319 SetParameters(Values);
1320}
#define _noexcept_
KS: noexcept can help with performance but is terrible for debugging, this is meant to help easy way ...
Definition: Core.h:83
int size
#define MACH3LOG_CRITICAL
Definition: MaCh3Logger.h:26
#define MACH3LOG_DEBUG
Definition: MaCh3Logger.h:22
#define MACH3LOG_ERROR
Definition: MaCh3Logger.h:25
#define MACH3LOG_INFO
Definition: MaCh3Logger.h:23
#define MACH3LOG_WARN
Definition: MaCh3Logger.h:24
Type Get(const YAML::Node &node, const std::string File, const int Line)
Get content of config file.
Definition: YamlHelper.h:265
#define M3OpenConfig(filename)
Macro to simplify calling LoadYaml with file and line info.
Definition: YamlHelper.h:560
#define GetBounds(filename)
Definition: YamlHelper.h:561
Custom exception class for MaCh3 errors.
std::vector< int > CircularBoundsIndex
Indices of parameters with circular bounds.
void ToggleFixParameter(const int i)
fix parameter at prior values
bool use_adaptive
Are we using AMCMC?
virtual void ProposeStep()
Generate a new proposed state.
double * randParams
Random number taken from gaussian around prior error used for corr_throw.
double ** throwMatrixCholDecomp
Throw matrix that is being used in the fit, much faster as TMatrixDSym cache miss.
std::unique_ptr< adaptive_mcmc::AdaptiveMCMCHandler > AdaptiveHandler
Struct containing information about adaption.
virtual ~ParameterHandlerBase()
Destructor.
TMatrixDSym * invCovMatrix
The inverse covariance matrix.
std::vector< int > FlipParameterIndex
Indices of parameters with flip symmetry.
bool AppliesToSample(const int SystIndex, const std::string &SampleName) const
Check if parameter is affecting given sample name.
void AcceptStep() _noexcept_
Accepted this step.
void ThrowParProp(const double mag=1.)
Throw the proposed parameter by mag sigma. Should really just have the user specify this throw by hav...
std::unique_ptr< ParameterTunes > Tunes
Struct containing information about adaption.
void InitialiseAdaption(const YAML::Node &adapt_manager)
Initialise adaptive MCMC.
TMatrixDSym * throwMatrix
Matrix which we use for step proposal before Cholesky decomposition (not actually used for step propo...
void ThrowParCurr(const double mag=1.)
Helper function to throw the current parameter by mag sigma. Can study bias in MCMC with this; put di...
double _fGlobalStepScale
Global step scale applied to all params in this class.
void ReserveMemory(const int size)
Initialise vectors with parameters information.
std::vector< std::string > _fFancyNames
Fancy name for example rather than xsec_0 it is MAQE, useful for human reading.
bool doSpecialStepProposal
Check if any of special step proposal were enabled.
std::vector< bool > _fFlatPrior
Whether to apply flat prior or not.
void ThrowParameters()
Throw the parameters according to the covariance matrix. This shouldn't be used in MCMC code ase it c...
void Init(const std::string &name, const std::string &file)
Initialisation of the class using matrix from root file.
std::string matrixName
Name of cov matrix.
void MakeClosestPosDef(TMatrixDSym *cov)
HW: Finds closest possible positive definite matrix in Frobenius Norm ||.||_frob Where ||X||_frob=sqr...
void RandomConfiguration()
Randomly throw the parameters in their 1 sigma range.
TMatrixD * throwMatrix_CholDecomp
Matrix which we use for step proposal after Cholesky decomposition.
std::vector< double > FlipParameterPoint
Central points around which parameters are flipped.
void UpdateAdaptiveCovariance()
Method to update adaptive MCMC .
std::vector< double > _fError
Prior error on the parameter.
void PrintNominal() const
Print prior value for every parameter.
YAML::Node _fYAMLDoc
Stores config describing systematics.
bool IsParameterFixed(const int i) const
Is parameter fixed or not.
void UpdateThrowMatrix(TMatrixDSym *cov)
void CircularParBounds(const int i, const double LowBound, const double UpBound)
HW :: This method is a tad hacky but modular arithmetic gives me a headache.
void SaveUpdatedMatrixConfig()
KS: After step scale, prefit etc. value were modified save this modified config.
void ResetIndivStepScale()
Adaptive Step Tuning Stuff.
std::unique_ptr< PCAHandler > PCAObj
Struct containing information about PCA.
int _fNumPar
Number of systematic parameters.
double CalcLikelihood() const _noexcept_
Calc penalty term based on inverted covariance matrix.
void FlipParameterValue(const int index, const double FlipPoint)
KS: Flip parameter around given value, for example mass ordering around 0.
void ConstructPCA(const double eigen_threshold, int FirstPCAdpar, int LastPCAdpar)
CW: Calculate eigen values, prepare transition matrices and remove param based on defined threshold.
std::vector< double > _fLowBound
Lowest physical bound, parameter will not be able to go beyond it.
std::vector< double > _fCurrVal
Current value of the parameter.
void PrintNominalCurrProp() const
Print prior, current and proposed value for each parameter.
ParameterHandlerBase(const std::vector< std::string > &YAMLFile, std::string name, double threshold=-1, int FirstPCAdpar=-999, int LastPCAdpar=-999)
ETA - constructor for a YAML file.
std::vector< std::vector< double > > InvertCovMatrix
KS: Same as above but much faster as TMatrixDSym cache miss.
TMatrixDSym * covMatrix
The covariance matrix.
std::vector< double > _fPreFitValue
Parameter value dictated by the prior model. Based on it penalty term is calculated.
void Randomize() _noexcept_
"Randomize" the parameters in the covariance class for the proposed step. Used the proposal kernel an...
void PrintIndivStepScale() const
Print step scale for each parameter.
std::vector< std::string > _fNames
ETA _fNames is set automatically in the covariance class to be something like xsec_i,...
std::vector< std::pair< double, double > > CircularBoundsValues
Circular bounds for each parameter (lower, upper)
void EnableSpecialProposal(const YAML::Node &param, const int Index)
Enable special proposal.
double * corr_throw
Result of multiplication of Cholesky matrix and randParams.
std::vector< double > _fUpBound
Upper physical bound, parameter will not be able to go beyond it.
bool pca
perform PCA or not
std::vector< double > _fIndivStepScale
Individual step scale used by MCMC algorithm.
std::vector< std::unique_ptr< TRandom3 > > random_number
KS: Set Random numbers for each thread so each thread has different seed.
void MakePosDef(TMatrixDSym *cov=nullptr)
Make matrix positive definite by adding small values to diagonal, necessary for inverting matrix.
void SpecialStepProposal()
Perform Special Step Proposal.
int PrintLength
KS: This is used when printing parameters, sometimes we have super long parameters name,...
void CorrelateSteps() _noexcept_
Use Cholesky throw matrix for better step proposal.
std::vector< double > _fPropVal
Proposed value of the parameter.
std::vector< std::vector< std::string > > _fSampleNames
Tells to which samples object param should be applied.
int CheckBounds() const _noexcept_
Check if parameters were proposed outside physical boundary.
void ToggleFixAllParameters()
fix parameters at prior values
std::string GetParName(const int i) const
Get name of parameter.
std::vector< double > GetProposed() const
Get vector of all proposed parameter values.
virtual double GetLikelihood()
Return CalcLikelihood if some params were thrown out of boundary return LARGE_LOGL
int GetParIndex(const std::string &name) const
Get index based on name.
std::string GetParFancyName(const int i) const
Get fancy name of the Parameter.
TH2D * GetCorrelationMatrix()
KS: Convert covariance matrix to correlation matrix and return TH2D which can be used for fancy plott...
std::string GetName() const
Get name of covariance.
double GetParInit(const int i) const
Get prior parameter value.
double GetDiagonalError(const int i) const
Get diagonal error for ith parameter.
void SetName(const std::string &name)
Set matrix name.
void SetCovMatrix(TMatrixDSym *cov)
Set covariance matrix.
void SetFlatPrior(const int i, const bool eL)
Set if parameter should have flat prior or not.
void SetStepScale(const double scale)
Set global step scale for covariance object.
void SetThrowMatrix(TMatrixDSym *cov)
Use new throw matrix, used in adaptive MCMC.
void SetParameters(const std::vector< double > &pars={})
Set parameter values using vector, it has to have same size as covariance class.
void SetBranches(TTree &tree, const bool SaveProposal=false)
set branches for output file
void SetParCurrProp(const int i, const double val)
Set current parameter value.
void SetSingleParameter(const int parNo, const double parVal)
Set value of single param to a given value.
void SetPar(const int i, const double val)
Set all the covariance matrix parameters to a user-defined value.
void SetIndivStepScale(const int ParameterIndex, const double StepScale)
DB Function to set fIndivStepScale from a vector (Can be used from execs and inside covariance constr...
void SetTune(const std::string &TuneName)
KS: Set proposed parameter values vector to be base on tune values, for example set proposed values t...
void MatrixVectorMulti(double *_restrict_ VecMulti, double **_restrict_ matrix, const double *_restrict_ vector, const int n)
KS: Custom function to perform multiplication of matrix and vector with multithreading.
int GetThreadIndex()
thread index inside parallel loop
Definition: Monitor.h:82
double MatrixVectorMultiSingle(double **_restrict_ matrix, const double *_restrict_ vector, const int Length, const int i)
KS: Custom function to perform multiplication of matrix and single element which is thread safe.
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
int GetNThreads()
number of threads which we need for example for TRandom3
Definition: Monitor.cpp:322
static constexpr const int _BAD_INT_
Default value used for int initialisation.
Definition: Core.h:45
std::string FormatDouble(const double value, const int precision)
Convert double into string for precision, useful for playing with yaml if you don't want to have in c...