19 std::vector<double>* fProp_Val) {
27 const double eigen_thresh,
const int _fNumPar) {
34 if (CovMatrix == NULL) {
41 MACH3LOG_ERROR(
"FirstPCAdpar and LastPCAdpar are higher than the number of parameters");
46 MACH3LOG_ERROR(
"FirstPCAdpar and LastPCAdpar are less than 0 but not default -999");
59 TMatrixDSymEigen eigen(submat);
97 temp2.ResizeTo(CovMatrix->GetNrows(),
NumParPCA);
100 for(
int iRow = 0; iRow < CovMatrix->GetNrows(); iRow++){
101 for(
int iCol = 0; iCol <
NumParPCA; iCol++){
102 temp2[iRow][iCol] = 0;
108 temp2[iRow][iRow] = 1;
117 for(
int iRow = 0;iRow < (CovMatrix->GetNrows()-1)-
LastPCAdpar; iRow++){
130 DebugPCA(sum, temp, submat, CovMatrix->GetNrows());
156 constexpr double correlation_threshold = 1e-6;
158 bool found_significant_correlation =
false;
160 int N = CovMatrix->GetNrows();
162 for (
int j = 0; j < N; ++j) {
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);
175 if (found_significant_correlation) {
176 MACH3LOG_ERROR(
"There are correlations between undecomposed and decomposed part of matrices, this will not work");
187 #pragma omp parallel for
199 const double GlobalStepScale,
205 #pragma omp parallel for
207 for (
int i = 0; i < NumParPCA; ++i)
209 if (_fErrorPCA[i] > 0.)
211 double IndStepScale = 1.;
213 if(isDecomposedPCA[i] >= 0)
215 IndStepScale *= IndivStepScale[isDecomposedPCA[i]];
216 IndStepScale *= corr_throw[isDecomposedPCA[i]];
221 IndStepScale *= randParams[i];
223 IndStepScale *= IndivStepScale[FirstPCAdpar];
225 _fParPropPCA(i) = _fParCurrPCA(i)+GlobalStepScale*IndStepScale*eigen_values_master[i];
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];
270 #pragma omp parallel for
272 for(
int i = 0; i < static_cast<int>(
_pCurrVal->size()); ++i) {
273 (*_pPropVal)[i] = fParProp_vec(i);
274 (*_pCurrVal)[i] = fParCurr_vec(i);
315 tree.Branch(Form(
"%s_PCA", Names[i].c_str()), &
_fParCurrPCA.GetMatrixArray()[i], Form(
"%s_PCA/D", Names[i].c_str()));
321 tree.Branch(Form(
"%s_PCA_Prop", Names[i].c_str()), &
_fParPropPCA.GetMatrixArray()[i], Form(
"%s_PCA_Prop/D", Names[i].c_str()));
342 MACH3LOG_ERROR(
"Parameter {} is PCA decomposed can't fix this", Names[i]);
346 MACH3LOG_INFO(
"Setting un-decomposed {}(parameter {}/{} in PCA base) to fixed at {}", Names[i], i, isDecom, (*
_pCurrVal)[i]);
353 double** throwMatrixCholDecomp,
356 const std::vector<double>& fPreFitValue,
357 const std::vector<double>& fLowBound,
358 const std::vector<double>& fUpBound,
359 const int _fNumPar) {
368 (*_pPropVal)[i] = fPreFitValue[i] + corr_throw[i];
374 (*_pPropVal)[i] = fPreFitValue[i] + corr_throw_single;
379 MACH3LOG_WARN(
"Tried {} times to throw parameter {} but failed",
throws, i);
382 (*_pPropVal)[i] = fPreFitValue[i];
386 (*_pCurrVal)[i] = (*_pPropVal)[i];
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];
403#pragma GCC diagnostic ignored "-Wfloat-conversion"
406void PCAHandler::DebugPCA(
const double sum, TMatrixD temp, TMatrixDSym submat,
int NumPar) {
409 TFile *PCA_Debug =
new TFile(
"Debug_PCA.root",
"RECREATE");
412 bool PlotText =
true;
414 if(NumPar > 200) PlotText =
false;
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");
422 double Cumulative = 0;
426 heigen_cumulative->SetBinContent(i+1, (
eigen_values)(i)/sum + Cumulative);
430 heigen_values->Write(
"heigen_values");
432 heigen_cumulative->Write(
"heigen_values_cumulative");
433 heigen_frac->Write(
"heigen_values_frac");
436 heigen_vectors->GetXaxis()->SetTitle(
"Parameter in Normal Base");
437 heigen_vectors->GetYaxis()->SetTitle(
"Parameter in Decomposed Base");
438 heigen_vectors->Write(
"heigen_vectors");
441 TH2D* SubsetPCA =
new TH2D(temp);
442 SubsetPCA->GetXaxis()->SetTitle(
"Parameter in Normal Base");
443 SubsetPCA->GetYaxis()->SetTitle(
"Parameter in Decomposed Base");
445 SubsetPCA->Write(
"hSubsetPCA");
446 temp.Write(
"SubsetPCA");
448 hTransferMat->GetXaxis()->SetTitle(
"Parameter in Normal Base");
449 hTransferMat->GetYaxis()->SetTitle(
"Parameter in Decomposed Base");
452 hTransferMatT->GetXaxis()->SetTitle(
"Parameter in Decomposed Base");
453 hTransferMatT->GetYaxis()->SetTitle(
"Parameter in Normal Base");
455 hTransferMat->Write(
"hTransferMat");
457 hTransferMatT->Write(
"hTransferMatT");
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);
467 gStyle->SetPaintTextFormat(
"4.1f");
468 gStyle->SetOptFit(0);
469 gStyle->SetOptStat(0);
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);
483 c1->Print(
"Debug_PCA.pdf[");
485 EigenLine->SetLineColor(kPink);
486 EigenLine->SetLineWidth(2);
487 EigenLine->SetLineStyle(kSolid);
489 TText* text =
new TText(0.5, 0.5, Form(
"Threshold = %g",
eigen_threshold));
490 text->SetTextFont (43);
491 text->SetTextSize (40);
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);
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");
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");
514 leg->SetLineColor(0);
515 leg->SetLineStyle(0);
516 leg->SetFillColor(0);
517 leg->SetFillStyle(0);
520 c1->Print(
"Debug_PCA.pdf");
521 c1->SetRightMargin(0.15);
526 delete heigen_values;
528 delete heigen_cumulative;
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");
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");
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");
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");
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;
575 Eigen::MatrixXd Submat_Eigen(submat.GetNrows(), submat.GetNcols());
578 #pragma omp parallel for
580 for(
int i = 0; i < submat.GetNrows(); i++)
582 for(
int j = 0; j < submat.GetNcols(); j++)
584 Submat_Eigen(i,j) = (submat)(i,j);
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++)
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];
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); } );
603 for(
auto const vect : eigen_vectors_and_values)
605 eigen_val(index) = std::get<0>(vect);
606 eigen_vect.row(index) = std::get<1>(vect);
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");
615 double Cumulative_Eigen = 0;
616 for(
int i = 0; i < eigen_val.size(); i++)
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;
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);
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");
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");
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");
648 c1->Print(
"Debug_PCA.pdf");
650 delete heigen_values_Eigen;
651 delete heigen_cumulative_Eigen;
652 delete heigen_frac_Eigen;
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());
657 for(
int i = 0; i < eigen_val.size(); i++)
659 for(
int j = 0; j < eigen_val.size(); j++)
662 heigen_vectors_Eigen->SetBinContent(i+1,j+1, eigen_vect(i,j));
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));
672 if(PlotText) heigen_vectors_Eigen->Draw(
"COLZ TEXT");
673 else heigen_vectors_Eigen->Draw(
"COLZ");
674 c1->Print(
"Debug_PCA.pdf");
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;
687 delete heigen_vectors;
689 c1->Print(
"Debug_PCA.pdf]");
#define _noexcept_
KS: noexcept can help with performance but is terrible for debugging, this is meant to help easy way ...
#define _restrict_
KS: Using restrict limits the effects of pointer aliasing, aiding optimizations. While reading I foun...
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...
TMatrixD TransferMat
Matrix used to converting from PCA base to normal base.
std::vector< double > eigen_values_master
Eigen values which have dimension equal to _fNumParPCA, and can be used in CorrelateSteps.
TMatrixD eigen_vectors
Eigen vectors only of params which are being decomposed.
TVectorD _fParPropPCA
CW: Current parameter value in PCA base.
void ThrowParProp(const double mag, const double *_restrict_ randParams)
Throw the proposed parameter by mag sigma.
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...
int NumParPCA
Number of parameters in PCA base.
std::vector< double > * _pCurrVal
Pointer to current value of the parameter.
TMatrixD TransferMatT
Matrix used to converting from normal base to PCA base.
void AcceptStep() _noexcept_
Accepted this step.
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.
void SetInitialParameters(std::vector< double > &IndStepScale)
KS: Transfer the starting parameters to the PCA basis, you don't want to start with zero....
void TransferToPCA()
Transfer param values from normal base to PCA base.
virtual ~PCAHandler()
Destructor.
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 ...
std::vector< double > * _pPropVal
Pointer to proposed value of the parameter.
TVectorD eigen_values
Eigen value only of particles which are being decomposed.
void TransferToParam()
Transfer param values from PCA base to normal base.
int LastPCAdpar
Index of the last param that is being decomposed.
int nKeptPCApars
Total number that remained after applying PCA Threshold.
std::vector< double > _fErrorPCA
Tells if parameter is fixed in PCA base or not.
void Print()
KS: Print info about PCA parameters.
std::vector< double > _fPreFitValuePCA
Prefit value for PCA params.
void ThrowParCurr(const double mag, const double *_restrict_ randParams)
Helper function to throw the current parameter by mag sigma.
TVectorD _fParCurrPCA
CW: Proposed parameter value in PCA base.
void SanitisePCA(TMatrixDSym *CovMatrix)
@biref KS: Make sure decomposed matrix isn't correlated with undecomposed
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.
int FirstPCAdpar
Index of the first param that is being decomposed.
double eigen_threshold
CW: Threshold based on which we remove parameters in eigen base.
bool IsParameterFixedPCA(const int i) const
Is parameter fixed in PCA base or not.
double GetParPropPCA(const int i) const
Get current parameter value using PCA.
double GetPreFitValuePCA(const int i) const
Get current parameter value using PCA.
bool IsParameterDecomposed(const int i) const
Check if parameter in PCA base is decomposed or not.
void ToggleFixAllParameters()
fix parameters at prior values
void SetBranches(TTree &tree, bool SaveProposal, const std::vector< std::string > &Names)
set branches for output file
void SetParPropPCA(const int i, const double value)
Set proposed value for parameter in PCA base.
void ToggleFixParameter(const int i, const std::vector< std::string > &Names)
fix parameters at prior values
void SetParCurrPCA(const int i, const double value)
Set current value for parameter in PCA base.
int GetThreadIndex()
thread index inside parallel loop
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.