MaCh3 2.2.1
Reference Guide
Loading...
Searching...
No Matches
PCAHandler.cpp
Go to the documentation of this file.
2
3
4// ********************************************
6// ********************************************
7 _pCurrVal = nullptr;
8 _pPropVal = nullptr;
9}
10
11// ********************************************
13// ********************************************
14
15}
16
17// ********************************************
18void PCAHandler::SetupPointers(std::vector<double>* fCurr_Val,
19 std::vector<double>* fProp_Val) {
20// ********************************************
21 _pCurrVal = fCurr_Val;
22 _pPropVal = fProp_Val;
23}
24
25// ********************************************
26void PCAHandler::ConstructPCA(TMatrixDSym* CovMatrix, const int firstPCAd, const int lastPCAd,
27 const double eigen_thresh, const int _fNumPar) {
28// ********************************************
29 FirstPCAdpar = firstPCAd;
30 LastPCAdpar = lastPCAd;
31 eigen_threshold = eigen_thresh;
32
33 // Check that covariance matrix exists
34 if (CovMatrix == NULL) {
35 MACH3LOG_ERROR("Covariance matrix for has not yet been set");
36 MACH3LOG_ERROR("Can not construct PCA until it is set");
37 throw MaCh3Exception(__FILE__ , __LINE__ );
38 }
39
40 if(FirstPCAdpar > CovMatrix->GetNrows()-1 || LastPCAdpar>CovMatrix->GetNrows()-1) {
41 MACH3LOG_ERROR("FirstPCAdpar and LastPCAdpar are higher than the number of parameters");
42 MACH3LOG_ERROR("first: {} last: {}, params: {}", FirstPCAdpar, LastPCAdpar, CovMatrix->GetNrows()-1);
43 throw MaCh3Exception(__FILE__ , __LINE__ );
44 }
45 if(FirstPCAdpar < 0 || LastPCAdpar < 0){
46 MACH3LOG_ERROR("FirstPCAdpar and LastPCAdpar are less than 0 but not default -999");
47 MACH3LOG_ERROR("first: {} last: {}", FirstPCAdpar, LastPCAdpar);
48 throw MaCh3Exception(__FILE__ , __LINE__ );
49 }
50 MACH3LOG_INFO("PCAing parameters {} through {} inclusive", FirstPCAdpar, LastPCAdpar);
51 int numunpcadpars = CovMatrix->GetNrows()-(LastPCAdpar-FirstPCAdpar+1);
52
53 // KS: Make sure we are not doing anything silly with PCA
54 SanitisePCA(CovMatrix);
55
56 TMatrixDSym submat(CovMatrix->GetSub(FirstPCAdpar,LastPCAdpar,FirstPCAdpar,LastPCAdpar));
57
58 //CW: Calculate how many eigen values this threshold corresponds to
59 TMatrixDSymEigen eigen(submat);
60 eigen_values.ResizeTo(eigen.GetEigenValues());
61 eigen_vectors.ResizeTo(eigen.GetEigenVectors());
62 eigen_values = eigen.GetEigenValues();
63 eigen_vectors = eigen.GetEigenVectors();
64 double sum = 0;
65 // Loop over eigen values and sum them up
66 for (int i = 0; i < eigen_values.GetNrows(); ++i) {
67 sum += eigen_values(i);
68 }
69 nKeptPCApars = eigen_values.GetNrows();
70 //CW: Now go through again and see how many eigen values correspond to threshold
71 for (int i = 0; i < eigen_values.GetNrows(); ++i) {
72 // Get the relative size of the eigen value
73 double sig = eigen_values(i)/sum;
74 // Check against the threshold
75 if (sig < eigen_threshold) {
76 nKeptPCApars = i;
77 break;
78 }
79 }
80 NumParPCA = numunpcadpars+nKeptPCApars;
81 MACH3LOG_INFO("Threshold of {} on eigen values relative sum of eigen value ({}) generates {} eigen vectors, plus we have {} unpcad pars, for a total of {}", eigen_threshold, sum, nKeptPCApars, numunpcadpars, NumParPCA);
82
83 //DB Create array of correct size so eigen_values can be used in CorrelateSteps
84 eigen_values_master = std::vector<double>(NumParPCA, 1.0);
86
87 // Now construct the transfer matrices
88 //These matrices will be as big as number of unPCAd pars plus number of eigenvalues kept
89 TransferMat.ResizeTo(CovMatrix->GetNrows(), NumParPCA);
90 TransferMatT.ResizeTo(CovMatrix->GetNrows(), NumParPCA);
91
92 // Get a subset of the eigen vector matrix
93 TMatrixD temp(eigen_vectors.GetSub(0, eigen_vectors.GetNrows()-1, 0, nKeptPCApars-1));
94
95 //Make transfer matrix which is two blocks of identity with a block of the PCA transfer matrix in between
96 TMatrixD temp2;
97 temp2.ResizeTo(CovMatrix->GetNrows(), NumParPCA);
98
99 //First set the whole thing to 0
100 for(int iRow = 0; iRow < CovMatrix->GetNrows(); iRow++){
101 for(int iCol = 0; iCol < NumParPCA; iCol++){
102 temp2[iRow][iCol] = 0;
103 }
104 }
105 //Set the first identity block
106 if(FirstPCAdpar != 0){
107 for(int iRow = 0; iRow < FirstPCAdpar; iRow++){
108 temp2[iRow][iRow] = 1;
109 }
110 }
111
112 //Set the transfer matrix block for the PCAd pars
113 temp2.SetSub(FirstPCAdpar,FirstPCAdpar,temp);
114
115 //Set the second identity block
116 if(LastPCAdpar != CovMatrix->GetNrows()-1){
117 for(int iRow = 0;iRow < (CovMatrix->GetNrows()-1)-LastPCAdpar; iRow++){
118 temp2[LastPCAdpar+1+iRow][FirstPCAdpar+nKeptPCApars+iRow] = 1;
119 }
120 }
121
122 TransferMat = temp2;
123 // Copy the contents
125 // And then transpose
126 TransferMatT.T();
127
128 #ifdef DEBUG_PCA
129 //KS: Let's dump all useful matrices to properly validate PCA
130 DebugPCA(sum, temp, submat, CovMatrix->GetNrows());
131 #endif
132
133 // Make the PCA parameter arrays
134 _fParCurrPCA.ResizeTo(NumParPCA);
135 _fParPropPCA.ResizeTo(NumParPCA);
137
138 //KS: make easy map so we could easily find un-decomposed parameters
140 _fErrorPCA.resize(NumParPCA);
141 for (int i = 0; i < NumParPCA; ++i)
142 {
143 _fErrorPCA[i] = 1;
144 isDecomposedPCA[i] = -1;
145 }
146 for (int i = 0; i < FirstPCAdpar; ++i) isDecomposedPCA[i] = i;
147
148 for (int i = FirstPCAdpar+nKeptPCApars+1; i < NumParPCA; ++i) isDecomposedPCA[i] = i+(_fNumPar-NumParPCA);
149}
150
151
152// ********************************************
153// Make sure decomposed matrix isn't correlated with undecomposed
154void PCAHandler::SanitisePCA(TMatrixDSym* CovMatrix) {
155// ********************************************
156 constexpr double correlation_threshold = 1e-6;
157
158 bool found_significant_correlation = false;
159
160 int N = CovMatrix->GetNrows();
161 for (int i = FirstPCAdpar; i <= LastPCAdpar; ++i) {
162 for (int j = 0; j < N; ++j) {
163 // Skip if j is inside the decomposed range (we only want cross-correlations)
164 if (j >= FirstPCAdpar && j <= LastPCAdpar) continue;
165
166 double corr_val = (*CovMatrix)(i, j);
167 if (std::fabs(corr_val) > correlation_threshold) {
168 found_significant_correlation = true;
169 MACH3LOG_ERROR("Significant correlation detected between decomposed parameter '{}' "
170 "and undecomposed parameter '{}': {:.6e}", i, j, corr_val);
171 }
172 }
173 }
174
175 if (found_significant_correlation) {
176 MACH3LOG_ERROR("There are correlations between undecomposed and decomposed part of matrices, this will not work");
177 throw MaCh3Exception(__FILE__ , __LINE__);
178 }
179}
180
181// ********************************************
182// Update so that current step becomes the previously proposed step
184// ********************************************
185 // Update the book-keeping for the output
186 #ifdef MULTITHREAD
187 #pragma omp parallel for
188 #endif
189 for (int i = 0; i < NumParPCA; ++i) {
191 }
192 // Then update the parameter basis
194}
195
196// ************************************************
197// Correlate the steps by setting the proposed step of a parameter to its current value + some correlated throw
198void PCAHandler::CorrelateSteps(const std::vector<double>& IndivStepScale,
199 const double GlobalStepScale,
200 const double* _restrict_ randParams,
201 const double* _restrict_ corr_throw) _noexcept_ {
202// ************************************************
203 // Throw around the current step
204 #ifdef MULTITHREAD
205 #pragma omp parallel for
206 #endif
207 for (int i = 0; i < NumParPCA; ++i)
208 {
209 if (_fErrorPCA[i] > 0.)
210 {
211 double IndStepScale = 1.;
212 //KS: If undecomposed parameter apply individual step scale and Cholesky for better acceptance rate
213 if(isDecomposedPCA[i] >= 0)
214 {
215 IndStepScale *= IndivStepScale[isDecomposedPCA[i]];
216 IndStepScale *= corr_throw[isDecomposedPCA[i]];
217 }
218 //If decomposed apply only random number
219 else
220 {
221 IndStepScale *= randParams[i];
222 //KS: All PCA-ed parameters have the same step scale
223 IndStepScale *= IndivStepScale[FirstPCAdpar];
224 }
225 _fParPropPCA(i) = _fParCurrPCA(i)+GlobalStepScale*IndStepScale*eigen_values_master[i];
226 }
227 }
228 // Then update the parameter basis
229 TransferToParam();
230}
231
232// ********************************************
233// Transfer a parameter variation in the parameter basis to the eigen basis
235// ********************************************
236 // Make the temporary vectors
237 TVectorD fParCurr_vec(static_cast<Int_t>(_pCurrVal->size()));
238 TVectorD fParProp_vec(static_cast<Int_t>(_pCurrVal->size()));
239 for(int i = 0; i < static_cast<int>(_pCurrVal->size()); ++i) {
240 fParCurr_vec(i) = (*_pCurrVal)[i];
241 fParProp_vec(i) = (*_pPropVal)[i];
242 }
243
244 _fParCurrPCA = TransferMatT*fParCurr_vec;
245 _fParPropPCA = TransferMatT*fParProp_vec;
246}
247
248// ********************************************
249void PCAHandler::SetInitialParameters(std::vector<double>& IndStepScale) {
250// ********************************************
252 for (int i = 0; i < NumParPCA; ++i) {
254 }
255 //DB Set Individual Step scale for PCA parameters to the LastPCAdpar fIndivStepScale because the step scale for those parameters is set by 'eigen_values[i]' but needs an overall step scale
256 // However, individual step scale for non-PCA parameters needs to be set correctly
257 for (int i = FirstPCAdpar; i <= LastPCAdpar; i++) {
258 IndStepScale[i] = IndStepScale[LastPCAdpar-1];
259 }
260}
261
262// ********************************************
263// Transfer a parameter variation in the eigen basis to the parameter basis
265// ********************************************
266 // Make the temporary vectors
267 TVectorD fParProp_vec = TransferMat*_fParPropPCA;
268 TVectorD fParCurr_vec = TransferMat*_fParCurrPCA;
269 #ifdef MULTITHREAD
270 #pragma omp parallel for
271 #endif
272 for(int i = 0; i < static_cast<int>(_pCurrVal->size()); ++i) {
273 (*_pPropVal)[i] = fParProp_vec(i);
274 (*_pCurrVal)[i] = fParCurr_vec(i);
275 }
276}
277
278// ********************************************
279// Throw the proposed parameter by mag sigma.
280void PCAHandler::ThrowParProp(const double mag, const double* _restrict_ randParams) {
281// ********************************************
282 for (int i = 0; i < NumParPCA; i++) {
283 if (_fErrorPCA[i] > 0.) {
284 _fParPropPCA(i) = _fParCurrPCA(i)+mag*randParams[i];
285 }
286 }
288}
289
290// ********************************************
291// Throw the proposed parameter by mag sigma.
292void PCAHandler::ThrowParCurr(const double mag, const double* _restrict_ randParams) {
293// ********************************************
294 for (int i = 0; i < NumParPCA; i++) {
295 if (_fErrorPCA[i] > 0.) {
296 _fParPropPCA(i) = mag*randParams[i];
297 }
298 }
300}
301
302// ********************************************
304// ********************************************
305 MACH3LOG_INFO("PCA:");
306 for (int i = 0; i < NumParPCA; ++i) {
307 MACH3LOG_INFO("PCA {:<2} Current: {:<10.2f} Proposed: {:<10.2f}", i, _fParCurrPCA(i), _fParPropPCA(i));
308 }
309}
310
311// ********************************************
312void PCAHandler::SetBranches(TTree &tree, bool SaveProposal, const std::vector<std::string>& Names) {
313// ********************************************
314 for (int i = 0; i < NumParPCA; ++i) {
315 tree.Branch(Form("%s_PCA", Names[i].c_str()), &_fParCurrPCA.GetMatrixArray()[i], Form("%s_PCA/D", Names[i].c_str()));
316 }
317
318 if(SaveProposal)
319 {
320 for (int i = 0; i < NumParPCA; ++i) {
321 tree.Branch(Form("%s_PCA_Prop", Names[i].c_str()), &_fParPropPCA.GetMatrixArray()[i], Form("%s_PCA_Prop/D", Names[i].c_str()));
322 }
323 }
324}
325
326// ********************************************
328// ********************************************
329 for (int i = 0; i < NumParPCA; i++) {
330 _fErrorPCA[i] *= -1.0;
331 }
332}
333
334// ********************************************
335void PCAHandler::ToggleFixParameter(const int i, const std::vector<std::string>& Names) {
336// ********************************************
337 int isDecom = -1;
338 for (int im = 0; im < NumParPCA; ++im) {
339 if(isDecomposedPCA[im] == i) {isDecom = im;}
340 }
341 if(isDecom < 0) {
342 MACH3LOG_ERROR("Parameter {} is PCA decomposed can't fix this", Names[i]);
343 //throw MaCh3Exception(__FILE__ , __LINE__ );
344 } else {
345 _fErrorPCA[isDecom] *= -1.0;
346 MACH3LOG_INFO("Setting un-decomposed {}(parameter {}/{} in PCA base) to fixed at {}", Names[i], i, isDecom, (*_pCurrVal)[i]);
347 }
348}
349
350
351// ********************************************
352void PCAHandler::ThrowParameters(const std::vector<std::unique_ptr<TRandom3>>& random_number,
353 double** throwMatrixCholDecomp,
354 double* randParams,
355 double* corr_throw,
356 const std::vector<double>& fPreFitValue,
357 const std::vector<double>& fLowBound,
358 const std::vector<double>& fUpBound,
359 const int _fNumPar) {
360// ********************************************
361 //KS: Do not multithread!
362 for (int i = 0; i < NumParPCA; ++i) {
363 // Check if parameter is fixed first: if so don't randomly throw
364 if (IsParameterFixedPCA(i)) continue;
365
367 {
368 (*_pPropVal)[i] = fPreFitValue[i] + corr_throw[i];
369 int throws = 0;
370 // Try again if we the initial parameter proposal falls outside of the range of the parameter
371 while ((*_pPropVal)[i] > fUpBound[i] || (*_pPropVal)[i] < fLowBound[i]) {
372 randParams[i] = random_number[M3::GetThreadIndex()]->Gaus(0, 1);
373 const double corr_throw_single = M3::MatrixVectorMultiSingle(throwMatrixCholDecomp, randParams, _fNumPar, i);
374 (*_pPropVal)[i] = fPreFitValue[i] + corr_throw_single;
375 if (throws > 10000)
376 {
377 //KS: Since we are multithreading there is danger that those messages
378 //will be all over the place, small price to pay for faster code
379 MACH3LOG_WARN("Tried {} times to throw parameter {} but failed", throws, i);
380 MACH3LOG_WARN("Setting _fPropVal: {} to {}", (*_pPropVal)[i], fPreFitValue[i]);
381 MACH3LOG_WARN("I live at {}:{}", __FILE__, __LINE__);
382 (*_pPropVal)[i] = fPreFitValue[i];
383 }
384 throws++;
385 }
386 (*_pCurrVal)[i] = (*_pPropVal)[i];
387
388 } else {
389 // KS: We have to multiply by number of parameters in PCA base
392 }
393 } // end of parameter loop
394
396 for (int i = 0; i < _fNumPar; ++i) {
397 (*_pPropVal)[i] = std::max(fLowBound[i], std::min((*_pPropVal)[i], fUpBound[i]));
398 (*_pCurrVal)[i] = (*_pPropVal)[i];
399 }
400}
401
402#ifdef DEBUG_PCA
403#pragma GCC diagnostic ignored "-Wfloat-conversion"
404// ********************************************
405//KS: Let's dump all useful matrices to properly validate PCA
406void PCAHandler::DebugPCA(const double sum, TMatrixD temp, TMatrixDSym submat, int NumPar) {
407// ********************************************
408 (void)submat;//This is used if DEBUG_PCA==2, this hack is to avoid compiler warnings
409 TFile *PCA_Debug = new TFile("Debug_PCA.root", "RECREATE");
410 PCA_Debug->cd();
411
412 bool PlotText = true;
413 //KS: If we have more than 200 plot becomes unreadable :(
414 if(NumPar > 200) PlotText = false;
415
416 TH1D* heigen_values = new TH1D("eigen_values", "Eigen Values", eigen_values.GetNrows(), 0.0, eigen_values.GetNrows());
417 TH1D* heigen_cumulative = new TH1D("heigen_cumulative", "heigen_cumulative", eigen_values.GetNrows(), 0.0, eigen_values.GetNrows());
418 TH1D* heigen_frac = new TH1D("heigen_fractional", "heigen_fractional", eigen_values.GetNrows(), 0.0, eigen_values.GetNrows());
419 heigen_values->GetXaxis()->SetTitle("Eigen Vector");
420 heigen_values->GetYaxis()->SetTitle("Eigen Value");
421
422 double Cumulative = 0;
423 for(int i = 0; i < eigen_values.GetNrows(); i++)
424 {
425 heigen_values->SetBinContent(i+1, (eigen_values)(i));
426 heigen_cumulative->SetBinContent(i+1, (eigen_values)(i)/sum + Cumulative);
427 heigen_frac->SetBinContent(i+1, (eigen_values)(i)/sum);
428 Cumulative += (eigen_values)(i)/sum;
429 }
430 heigen_values->Write("heigen_values");
431 eigen_values.Write("eigen_values");
432 heigen_cumulative->Write("heigen_values_cumulative");
433 heigen_frac->Write("heigen_values_frac");
434
435 TH2D* heigen_vectors = new TH2D(eigen_vectors);
436 heigen_vectors->GetXaxis()->SetTitle("Parameter in Normal Base");
437 heigen_vectors->GetYaxis()->SetTitle("Parameter in Decomposed Base");
438 heigen_vectors->Write("heigen_vectors");
439 eigen_vectors.Write("eigen_vectors");
440
441 TH2D* SubsetPCA = new TH2D(temp);
442 SubsetPCA->GetXaxis()->SetTitle("Parameter in Normal Base");
443 SubsetPCA->GetYaxis()->SetTitle("Parameter in Decomposed Base");
444
445 SubsetPCA->Write("hSubsetPCA");
446 temp.Write("SubsetPCA");
447 TH2D* hTransferMat = new TH2D(TransferMat);
448 hTransferMat->GetXaxis()->SetTitle("Parameter in Normal Base");
449 hTransferMat->GetYaxis()->SetTitle("Parameter in Decomposed Base");
450 TH2D* hTransferMatT = new TH2D(TransferMatT);
451
452 hTransferMatT->GetXaxis()->SetTitle("Parameter in Decomposed Base");
453 hTransferMatT->GetYaxis()->SetTitle("Parameter in Normal Base");
454
455 hTransferMat->Write("hTransferMat");
456 TransferMat.Write("TransferMat");
457 hTransferMatT->Write("hTransferMatT");
458 TransferMatT.Write("TransferMatT");
459
460 TCanvas *c1 = new TCanvas("c1"," ", 0, 0, 1024, 1024);
461 c1->SetBottomMargin(0.1);
462 c1->SetTopMargin(0.05);
463 c1->SetRightMargin(0.05);
464 c1->SetLeftMargin(0.12);
465 c1->SetGrid();
466
467 gStyle->SetPaintTextFormat("4.1f");
468 gStyle->SetOptFit(0);
469 gStyle->SetOptStat(0);
470 // Make pretty correlation colors (red to blue)
471 const int NRGBs = 5;
472 TColor::InitializeColors();
473 Double_t stops[NRGBs] = { 0.00, 0.25, 0.50, 0.75, 1.00 };
474 Double_t red[NRGBs] = { 0.00, 0.25, 1.00, 1.00, 0.50 };
475 Double_t green[NRGBs] = { 0.00, 0.25, 1.00, 0.25, 0.00 };
476 Double_t blue[NRGBs] = { 0.50, 1.00, 1.00, 0.25, 0.00 };
477 TColor::CreateGradientColorTable(5, stops, red, green, blue, 255);
478 gStyle->SetNumberContours(255);
479
480 double maxz = 0;
481 double minz = 0;
482
483 c1->Print("Debug_PCA.pdf[");
484 TLine *EigenLine = new TLine(nKeptPCApars, 0, nKeptPCApars, heigen_cumulative->GetMaximum());
485 EigenLine->SetLineColor(kPink);
486 EigenLine->SetLineWidth(2);
487 EigenLine->SetLineStyle(kSolid);
488
489 TText* text = new TText(0.5, 0.5, Form("Threshold = %g", eigen_threshold));
490 text->SetTextFont (43);
491 text->SetTextSize (40);
492
493 heigen_values->SetLineColor(kRed);
494 heigen_values->SetLineWidth(2);
495 heigen_cumulative->SetLineColor(kGreen);
496 heigen_cumulative->SetLineWidth(2);
497 heigen_frac->SetLineColor(kBlue);
498 heigen_frac->SetLineWidth(2);
499
500 c1->SetLogy();
501 heigen_values->SetMaximum(heigen_cumulative->GetMaximum()+heigen_cumulative->GetMaximum()*0.4);
502 heigen_values->Draw();
503 heigen_frac->Draw("SAME");
504 heigen_cumulative->Draw("SAME");
505 EigenLine->Draw("Same");
506 text->DrawTextNDC(0.42, 0.84,Form("Threshold = %g", eigen_threshold));
507
508 TLegend *leg = new TLegend(0.2, 0.2, 0.6, 0.5);
509 leg->SetTextSize(0.04);
510 leg->AddEntry(heigen_values, "Absolute", "l");
511 leg->AddEntry(heigen_frac, "Fractional", "l");
512 leg->AddEntry(heigen_cumulative, "Cumulative", "l");
513
514 leg->SetLineColor(0);
515 leg->SetLineStyle(0);
516 leg->SetFillColor(0);
517 leg->SetFillStyle(0);
518 leg->Draw("Same");
519
520 c1->Print("Debug_PCA.pdf");
521 c1->SetRightMargin(0.15);
522 c1->SetLogy(0);
523 delete EigenLine;
524 delete leg;
525 delete text;
526 delete heigen_values;
527 delete heigen_frac;
528 delete heigen_cumulative;
529
530 heigen_vectors->SetMarkerSize(0.2);
531 minz = heigen_vectors->GetMinimum();
532 if (fabs(0-maxz)>fabs(0-minz)) heigen_vectors->GetZaxis()->SetRangeUser(0-fabs(0-maxz),0+fabs(0-maxz));
533 else heigen_vectors->GetZaxis()->SetRangeUser(0-fabs(0-minz),0+fabs(0-minz));
534 if(PlotText) heigen_vectors->Draw("COLZ TEXT");
535 else heigen_vectors->Draw("COLZ");
536
537 TLine *Eigen_Line = new TLine(0, nKeptPCApars, LastPCAdpar-FirstPCAdpar, nKeptPCApars);
538 Eigen_Line->SetLineColor(kGreen);
539 Eigen_Line->SetLineWidth(2);
540 Eigen_Line->SetLineStyle(kDotted);
541 Eigen_Line->Draw("SAME");
542 c1->Print("Debug_PCA.pdf");
543 delete Eigen_Line;
544
545 SubsetPCA->SetMarkerSize(0.2);
546 minz = SubsetPCA->GetMinimum();
547 if (fabs(0-maxz)>fabs(0-minz)) SubsetPCA->GetZaxis()->SetRangeUser(0-fabs(0-maxz),0+fabs(0-maxz));
548 else SubsetPCA->GetZaxis()->SetRangeUser(0-fabs(0-minz),0+fabs(0-minz));
549 if(PlotText) SubsetPCA->Draw("COLZ TEXT");
550 else SubsetPCA->Draw("COLZ");
551 c1->Print("Debug_PCA.pdf");
552 delete SubsetPCA;
553
554 hTransferMat->SetMarkerSize(0.15);
555 minz = hTransferMat->GetMinimum();
556 if (fabs(0-maxz)>fabs(0-minz)) hTransferMat->GetZaxis()->SetRangeUser(0-fabs(0-maxz),0+fabs(0-maxz));
557 else hTransferMat->GetZaxis()->SetRangeUser(0-fabs(0-minz),0+fabs(0-minz));
558 if(PlotText) hTransferMat->Draw("COLZ TEXT");
559 else hTransferMat->Draw("COLZ");
560 c1->Print("Debug_PCA.pdf");
561 delete hTransferMat;
562
563 hTransferMatT->SetMarkerSize(0.15);
564 minz = hTransferMatT->GetMinimum();
565 if (fabs(0-maxz)>fabs(0-minz)) hTransferMatT->GetZaxis()->SetRangeUser(0-fabs(0-maxz),0+fabs(0-maxz));
566 else hTransferMatT->GetZaxis()->SetRangeUser(0-fabs(0-minz),0+fabs(0-minz));
567 if(PlotText) hTransferMatT->Draw("COLZ TEXT");
568 else hTransferMatT->Draw("COLZ");
569 c1->Print( "Debug_PCA.pdf");
570 delete hTransferMatT;
571
572
573 //KS: Crosscheck against Eigen library
574 #if DEBUG_PCA == 2
575 Eigen::MatrixXd Submat_Eigen(submat.GetNrows(), submat.GetNcols());
576
577 #ifdef MULTITHREAD
578 #pragma omp parallel for
579 #endif
580 for(int i = 0; i < submat.GetNrows(); i++)
581 {
582 for(int j = 0; j < submat.GetNcols(); j++)
583 {
584 Submat_Eigen(i,j) = (submat)(i,j);
585 }
586 }
587 Eigen::EigenSolver<Eigen::MatrixXd> EigenSolver;
588 EigenSolver.compute(Submat_Eigen);
589 Eigen::VectorXd eigen_val = EigenSolver.eigenvalues().real();
590 Eigen::MatrixXd eigen_vect = EigenSolver.eigenvectors().real();
591 std::vector<std::tuple<double, Eigen::VectorXd>> eigen_vectors_and_values;
592 double Sum_Eigen = 0;
593 for(int i = 0; i < eigen_val.size(); i++)
594 {
595 std::tuple<double, Eigen::VectorXd> vec_and_val(eigen_val[i], eigen_vect.row(i));
596 eigen_vectors_and_values.push_back(vec_and_val);
597 Sum_Eigen += eigen_val[i];
598 }
599 std::sort(eigen_vectors_and_values.begin(), eigen_vectors_and_values.end(),
600 [&](const std::tuple<double, Eigen::VectorXd>& a, const std::tuple<double, Eigen::VectorXd>& b) -> bool
601 { return std::get<0>(a) > std::get<0>(b); } );
602 int index = 0;
603 for(auto const vect : eigen_vectors_and_values)
604 {
605 eigen_val(index) = std::get<0>(vect);
606 eigen_vect.row(index) = std::get<1>(vect);
607 index++;
608 }
609 TH1D* heigen_values_Eigen = new TH1D("eig_values", "Eigen Values", eigen_val.size(), 0.0, eigen_val.size());
610 TH1D* heigen_cumulative_Eigen = new TH1D("eig_cumulative", "heigen_cumulative", eigen_val.size(), 0.0, eigen_val.size());
611 TH1D* heigen_frac_Eigen = new TH1D("eig_fractional", "heigen_fractional", eigen_val.size(), 0.0, eigen_val.size());
612 heigen_values_Eigen->GetXaxis()->SetTitle("Eigen Vector");
613 heigen_values_Eigen->GetYaxis()->SetTitle("Eigen Value");
614
615 double Cumulative_Eigen = 0;
616 for(int i = 0; i < eigen_val.size(); i++)
617 {
618 heigen_values_Eigen->SetBinContent(i+1, eigen_val(i));
619 heigen_cumulative_Eigen->SetBinContent(i+1, eigen_val(i)/sum + Cumulative_Eigen);
620 heigen_frac_Eigen->SetBinContent(i+1, eigen_val(i)/sum);
621 Cumulative_Eigen += eigen_val(i)/sum;
622 }
623 heigen_values_Eigen->SetLineColor(kRed);
624 heigen_values_Eigen->SetLineWidth(2);
625 heigen_cumulative_Eigen->SetLineColor(kGreen);
626 heigen_cumulative_Eigen->SetLineWidth(2);
627 heigen_frac_Eigen->SetLineColor(kBlue);
628 heigen_frac_Eigen->SetLineWidth(2);
629
630 c1->SetLogy();
631 heigen_values_Eigen->SetMaximum(heigen_cumulative_Eigen->GetMaximum()+heigen_cumulative_Eigen->GetMaximum()*0.4);
632 heigen_values_Eigen->Draw();
633 heigen_cumulative_Eigen->Draw("SAME");
634 heigen_frac_Eigen->Draw("SAME");
635
636 TLegend *leg_Eigen = new TLegend(0.2, 0.2, 0.6, 0.5);
637 leg_Eigen->SetTextSize(0.04);
638 leg_Eigen->AddEntry(heigen_values_Eigen, "Absolute", "l");
639 leg_Eigen->AddEntry(heigen_frac_Eigen, "Fractional", "l");
640 leg_Eigen->AddEntry(heigen_cumulative_Eigen, "Cumulative", "l");
641
642 leg_Eigen->SetLineColor(0);
643 leg_Eigen->SetLineStyle(0);
644 leg_Eigen->SetFillColor(0);
645 leg_Eigen->SetFillStyle(0);
646 leg_Eigen->Draw("Same");
647
648 c1->Print("Debug_PCA.pdf");
649 c1->SetLogy(0);
650 delete heigen_values_Eigen;
651 delete heigen_cumulative_Eigen;
652 delete heigen_frac_Eigen;
653 delete leg_Eigen;
654
655 TH2D* heigen_vectors_Eigen = new TH2D("Eigen_Vectors", "Eigen_Vectors", eigen_val.size(), 0.0, eigen_val.size(), eigen_val.size(), 0.0, eigen_val.size());
656
657 for(int i = 0; i < eigen_val.size(); i++)
658 {
659 for(int j = 0; j < eigen_val.size(); j++)
660 {
661 //KS: +1 because there is offset in histogram relative to TMatrix
662 heigen_vectors_Eigen->SetBinContent(i+1,j+1, eigen_vect(i,j));
663 }
664 }
665 heigen_vectors_Eigen->GetXaxis()->SetTitle("Parameter in Normal Base");
666 heigen_vectors_Eigen->GetYaxis()->SetTitle("Parameter in Decomposed Base");
667 heigen_vectors_Eigen->SetMarkerSize(0.15);
668 minz = heigen_vectors_Eigen->GetMinimum();
669 if (fabs(0-maxz)>fabs(0-minz)) heigen_vectors_Eigen->GetZaxis()->SetRangeUser(0-fabs(0-maxz),0+fabs(0-maxz));
670 else heigen_vectors_Eigen->GetZaxis()->SetRangeUser(0-fabs(0-minz),0+fabs(0-minz));
671
672 if(PlotText) heigen_vectors_Eigen->Draw("COLZ TEXT");
673 else heigen_vectors_Eigen->Draw("COLZ");
674 c1->Print( "Debug_PCA.pdf");
675
676 heigen_vectors->SetTitle("ROOT minus Eigen");
677 heigen_vectors->Add(heigen_vectors_Eigen, -1.);
678 minz = heigen_vectors->GetMinimum();
679 if (fabs(0-maxz)>fabs(0-minz)) heigen_vectors->GetZaxis()->SetRangeUser(0-fabs(0-maxz),0+fabs(0-maxz));
680 else heigen_vectors->GetZaxis()->SetRangeUser(0-fabs(0-minz),0+fabs(0-minz));
681 if(PlotText) heigen_vectors->Draw("COLZ TEXT");
682 else heigen_vectors->Draw("COLZ");
683 c1->Print("Debug_PCA.pdf");
684 delete heigen_vectors_Eigen;
685
686 #endif // end if Eigen enabled
687 delete heigen_vectors;
688
689 c1->Print("Debug_PCA.pdf]");
690 delete c1;
691 PCA_Debug->Close();
692 delete PCA_Debug;
693}
694#endif
#define _noexcept_
KS: noexcept can help with performance but is terrible for debugging, this is meant to help easy way ...
Definition: Core.h:83
#define _restrict_
KS: Using restrict limits the effects of pointer aliasing, aiding optimizations. While reading I foun...
Definition: Core.h:90
#define MACH3LOG_ERROR
Definition: MaCh3Logger.h:25
#define MACH3LOG_INFO
Definition: MaCh3Logger.h:23
#define MACH3LOG_WARN
Definition: MaCh3Logger.h:24
Custom exception class for MaCh3 errors.
std::vector< int > isDecomposedPCA
If param is decomposed this will return -1, if not this will return enumerator to param in normal bas...
Definition: PCAHandler.h:211
TMatrixD TransferMat
Matrix used to converting from PCA base to normal base.
Definition: PCAHandler.h:216
std::vector< double > eigen_values_master
Eigen values which have dimension equal to _fNumParPCA, and can be used in CorrelateSteps.
Definition: PCAHandler.h:224
TMatrixD eigen_vectors
Eigen vectors only of params which are being decomposed.
Definition: PCAHandler.h:222
TVectorD _fParPropPCA
CW: Current parameter value in PCA base.
Definition: PCAHandler.h:205
void ThrowParProp(const double mag, const double *_restrict_ randParams)
Throw the proposed parameter by mag sigma.
Definition: PCAHandler.cpp:280
void ThrowParameters(const std::vector< std::unique_ptr< TRandom3 > > &random_number, double **throwMatrixCholDecomp, double *randParams, double *corr_throw, const std::vector< double > &fPreFitValue, const std::vector< double > &fLowBound, const std::vector< double > &fUpBound, int _fNumPar)
Throw the parameters according to the covariance matrix. This shouldn't be used in MCMC code ase it c...
Definition: PCAHandler.cpp:352
int NumParPCA
Number of parameters in PCA base.
Definition: PCAHandler.h:213
std::vector< double > * _pCurrVal
Pointer to current value of the parameter.
Definition: PCAHandler.h:236
TMatrixD TransferMatT
Matrix used to converting from normal base to PCA base.
Definition: PCAHandler.h:218
void AcceptStep() _noexcept_
Accepted this step.
Definition: PCAHandler.cpp:183
void CorrelateSteps(const std::vector< double > &IndivStepScale, const double GlobalStepScale, const double *_restrict_ randParams, const double *_restrict_ corr_throw) _noexcept_
Use Cholesky throw matrix for better step proposal.
Definition: PCAHandler.cpp:198
void SetInitialParameters(std::vector< double > &IndStepScale)
KS: Transfer the starting parameters to the PCA basis, you don't want to start with zero....
Definition: PCAHandler.cpp:249
void TransferToPCA()
Transfer param values from normal base to PCA base.
Definition: PCAHandler.cpp:234
virtual ~PCAHandler()
Destructor.
Definition: PCAHandler.cpp:12
void SetupPointers(std::vector< double > *fCurr_Val, std::vector< double > *fProp_Val)
KS: Setup pointers to current and proposed parameter value which we need to convert them to PCA base ...
Definition: PCAHandler.cpp:18
std::vector< double > * _pPropVal
Pointer to proposed value of the parameter.
Definition: PCAHandler.h:238
TVectorD eigen_values
Eigen value only of particles which are being decomposed.
Definition: PCAHandler.h:220
void TransferToParam()
Transfer param values from PCA base to normal base.
Definition: PCAHandler.cpp:264
int LastPCAdpar
Index of the last param that is being decomposed.
Definition: PCAHandler.h:231
int nKeptPCApars
Total number that remained after applying PCA Threshold.
Definition: PCAHandler.h:227
std::vector< double > _fErrorPCA
Tells if parameter is fixed in PCA base or not.
Definition: PCAHandler.h:209
PCAHandler()
Constructor.
Definition: PCAHandler.cpp:5
void Print()
KS: Print info about PCA parameters.
Definition: PCAHandler.cpp:303
std::vector< double > _fPreFitValuePCA
Prefit value for PCA params.
Definition: PCAHandler.h:203
void ThrowParCurr(const double mag, const double *_restrict_ randParams)
Helper function to throw the current parameter by mag sigma.
Definition: PCAHandler.cpp:292
TVectorD _fParCurrPCA
CW: Proposed parameter value in PCA base.
Definition: PCAHandler.h:207
void SanitisePCA(TMatrixDSym *CovMatrix)
@biref KS: Make sure decomposed matrix isn't correlated with undecomposed
Definition: PCAHandler.cpp:154
void ConstructPCA(TMatrixDSym *CovMatrix, const int firstPCAd, const int lastPCAd, const double eigen_thresh, const int _fNumPar)
CW: Calculate eigen values, prepare transition matrices and remove param based on defined threshold.
Definition: PCAHandler.cpp:26
int FirstPCAdpar
Index of the first param that is being decomposed.
Definition: PCAHandler.h:229
double eigen_threshold
CW: Threshold based on which we remove parameters in eigen base.
Definition: PCAHandler.h:233
bool IsParameterFixedPCA(const int i) const
Is parameter fixed in PCA base or not.
Definition: PCAHandler.h:117
double GetParPropPCA(const int i) const
Get current parameter value using PCA.
Definition: PCAHandler.h:150
double GetPreFitValuePCA(const int i) const
Get current parameter value using PCA.
Definition: PCAHandler.h:157
bool IsParameterDecomposed(const int i) const
Check if parameter in PCA base is decomposed or not.
Definition: PCAHandler.h:188
void ToggleFixAllParameters()
fix parameters at prior values
Definition: PCAHandler.cpp:327
void SetBranches(TTree &tree, bool SaveProposal, const std::vector< std::string > &Names)
set branches for output file
Definition: PCAHandler.cpp:312
void SetParPropPCA(const int i, const double value)
Set proposed value for parameter in PCA base.
Definition: PCAHandler.h:132
void ToggleFixParameter(const int i, const std::vector< std::string > &Names)
fix parameters at prior values
Definition: PCAHandler.cpp:335
void SetParCurrPCA(const int i, const double value)
Set current value for parameter in PCA base.
Definition: PCAHandler.h:141
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.