Class responsible for handling Principal Component Analysis (PCA) of covariance matrix. More...

#include <Parameters/PCAHandler.h>

Collaboration diagram for PCAHandler:

Public Member Functions
	PCAHandler ()
	Constructor. More...

virtual	~PCAHandler ()
	Destructor. More...

void	Print () const
	KS: Print info about PCA parameters. More...

int	GetNumberPCAedParameters () const
	Retrieve number of parameters in PCA base. More...

void	SetupPointers (std::vector< double > fCurr_Val, std::vector< M3::float_t > fProp_Val)
	KS: Setup pointers to current and proposed parameter value which we need to convert them to PCA base each step. More...

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. More...

void	TransferToPCA ()
	Transfer param values from normal base to PCA base. More...

void	TransferToParam ()
	Transfer param values from PCA base to normal base. More...

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 can break Detailed Balance;. More...

void	AcceptStep () _noexcept_
	Accepted this step. More...

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. More...

void	SetInitialParameters ()
	KS: Transfer the starting parameters to the PCA basis, you don't want to start with zero.. More...

void	SetBranches (TTree &tree, const bool SaveProposal, const std::vector< std::string > &Names)
	set branches for output file More...

void	ToggleFixAllParameters (const std::vector< std::string > &Names)
	fix parameters at prior values More...

void	ToggleFixParameter (const int i, const std::vector< std::string > &Names)
	fix parameters at prior values More...

void	SetParametersPCA (const std::vector< double > &pars)
	Set values for PCA parameters in PCA base. More...

bool	IsParameterFixedPCA (const int i) const
	Is parameter fixed in PCA base or not. More...

const TMatrixD	GetEigenVectors () const
	Get eigen vectors of covariance matrix, only works with PCA. More...

void	SetParPropPCA (const int i, const double value)
	Set proposed value for parameter in PCA base. More...

void	SetParCurrPCA (const int i, const double value)
	Set current value for parameter in PCA base. More...

double	GetParPropPCA (const int i) const
	Get current parameter value using PCA. More...

double	GetPreFitValuePCA (const int i) const
	Get current parameter value using PCA. More...

double	GetParCurrPCA (const int i) const
	Get current parameter value using PCA. More...

const TMatrixD	GetTransferMatrix () const
	Get transfer matrix allowing to go from PCA base to normal base. More...

const TVectorD	GetEigenValues () const
	Get eigen values for all parameters, if you want for decomposed only parameters use GetEigenValuesMaster. More...

const std::vector< double >	GetEigenValuesMaster () const
	Get eigen value of only decomposed parameters, if you want for all parameters use GetEigenValues. More...

bool	IsParameterDecomposed (const int i) const
	Check if parameter in PCA base is decomposed or not. More...

Private Member Functions
void	SanitisePCA (TMatrixDSym *CovMatrix)
	KS: Make sure decomposed matrix isn't correlated with undecomposed. More...

Private Attributes
std::vector< double >	_fPreFitValuePCA
	Prefit value for PCA params. More...

TVectorD	_fParPropPCA
	CW: Current parameter value in PCA base. More...

TVectorD	_fParCurrPCA
	CW: Proposed parameter value in PCA base. More...

std::vector< double >	_fErrorPCA
	Tells if parameter is fixed in PCA base or not. More...

std::vector< int >	isDecomposedPCA
	If param is decomposed this will return -1, if not this will return enumerator to param in normal base. This way we can map stuff like step scale etc between normal base and undecomposed param in eigen base. More...

int	NumParPCA
	Number of parameters in PCA base. More...

TMatrixD	TransferMat
	Matrix used to converting from PCA base to normal base. More...

TMatrixD	TransferMatT
	Matrix used to converting from normal base to PCA base. More...

TVectorD	eigen_values
	Eigen value only of particles which are being decomposed. More...

TMatrixD	eigen_vectors
	Eigen vectors only of params which are being decomposed. More...

std::vector< double >	eigen_values_master
	Eigen values which have dimension equal to _fNumParPCA, and can be used in CorrelateSteps. More...

int	nKeptPCApars
	Total number that remained after applying PCA Threshold. More...

int	FirstPCAdpar
	Index of the first param that is being decomposed. More...

int	LastPCAdpar
	Index of the last param that is being decomposed. More...

double	eigen_threshold
	CW: Threshold based on which we remove parameters in eigen base. More...

std::vector< double > *	_pCurrVal
	Pointer to current value of the parameter. More...

std::vector< M3::float_t > *	_pPropVal
	Pointer to proposed value of the parameter. More...

Detailed Description

Class responsible for handling Principal Component Analysis (PCA) of covariance matrix.

Author: Clarence Wret

Introduction

This class applies Principal Component Analysis (PCA) to covariance matrices using eigen-decomposition. The main benefit is that PCA rotates the parameter space into a basis where parameters are uncorrelated. In practice, this helps MCMC fits move more efficiently, since correlated directions in the original space often slow down convergence.

Dimensionality Reduction

PCA makes it possible to drop directions in parameter space with very small eigenvalues. These correspond to combinations of parameters that are poorly constrained or dominated by statistical noise. By applying a threshold (e.g. discarding eigenvalues smaller than 1e-4 of the largest one), the effective number of parameters can be reduced.

This is particularly useful in large MCMC fits, where many parameters may contribute little information. Removing them reduces dimensionality while keeping the dominant structure of the parameter space intact, often leading to faster and more stable fits.

Partial Decomposition

PCA does not need to be applied to the full covariance matrix. This class supports eigen-decomposition of only a submatrix. This is useful when only a subset of parameters is strongly correlated and benefits from PCA, while the rest remain in the original parameter basis.

Definition at line 38 of file PCAHandler.h.

Constructor & Destructor Documentation

◆ PCAHandler()

PCAHandler::PCAHandler ( )

Constructor.

Definition at line 17 of file PCAHandler.cpp.

                        {
 // ********************************************
   _pCurrVal = nullptr;
   _pPropVal = nullptr;
 }

◆ ~PCAHandler()

PCAHandler::~PCAHandler ( )

virtual

Destructor.

Definition at line 24 of file PCAHandler.cpp.

                         {
 // ********************************************
 }

Member Function Documentation

◆ AcceptStep()

void PCAHandler::AcceptStep ( )

Accepted this step.

Definition at line 191 of file PCAHandler.cpp.

                                        {
 // ********************************************
   // Update the book-keeping for the output
   #ifdef MULTITHREAD
   #pragma omp parallel for
   #endif
   for (int i = 0; i < NumParPCA; ++i) {
     _fParCurrPCA(i) = _fParPropPCA(i);
   }
   // Then update the parameter basis
   TransferToParam();
 }

◆ ConstructPCA()

void PCAHandler::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.

Parameters

CovMatrix	Symmetric covariance matrix used for eigen decomposition.
firstPCAd	Index of the first PCA component to include.
lastPCAd	Index of the last PCA component to include.
eigen_thresh	Threshold for eigenvalues below which parameters are discarded.
_fNumPar	Total number of parameters in the original (non-PCA) basis.

Definition at line 37 of file PCAHandler.cpp.

                                                                              {
 // ********************************************
   FirstPCAdpar = firstPCAd;
   LastPCAdpar = lastPCAd;
   eigen_threshold = eigen_thresh;
  
   // Check that covariance matrix exists
   if (CovMatrix == nullptr) {
     MACH3LOG_ERROR("Covariance matrix for has not yet been set");
     MACH3LOG_ERROR("Can not construct PCA until it is set");
     throw MaCh3Exception(__FILE__ , __LINE__ );
   }
  
   if(FirstPCAdpar > CovMatrix->GetNrows()-1 || LastPCAdpar>CovMatrix->GetNrows()-1) {
     MACH3LOG_ERROR("FirstPCAdpar and LastPCAdpar are higher than the number of parameters");
     MACH3LOG_ERROR("first: {} last: {}, params: {}", FirstPCAdpar, LastPCAdpar, CovMatrix->GetNrows()-1);
     throw MaCh3Exception(__FILE__ , __LINE__ );
   }
   if(FirstPCAdpar < 0 || LastPCAdpar < 0){
     MACH3LOG_ERROR("FirstPCAdpar and LastPCAdpar are less than 0 but not default -999");
     MACH3LOG_ERROR("first: {} last: {}", FirstPCAdpar, LastPCAdpar);
     throw MaCh3Exception(__FILE__ , __LINE__ );
   }
   MACH3LOG_INFO("PCAing parameters {} through {} inclusive", FirstPCAdpar, LastPCAdpar);
   int NumUnPCAdPars = CovMatrix->GetNrows()-(LastPCAdpar-FirstPCAdpar+1);
  
   // KS: Make sure we are not doing anything silly with PCA
   SanitisePCA(CovMatrix);
  
   TMatrixDSym submat(CovMatrix->GetSub(FirstPCAdpar,LastPCAdpar,FirstPCAdpar,LastPCAdpar));
  
   //CW: Calculate how many eigen values this threshold corresponds to
   TMatrixDSymEigen eigen(submat);
   eigen_values.ResizeTo(eigen.GetEigenValues());
   eigen_vectors.ResizeTo(eigen.GetEigenVectors());
   eigen_values = eigen.GetEigenValues();
   eigen_vectors = eigen.GetEigenVectors();
   double sum = 0;
   // Loop over eigen values and sum them up
   for (int i = 0; i < eigen_values.GetNrows(); ++i) {
     sum += eigen_values(i);
   }
   nKeptPCApars = eigen_values.GetNrows();
   //CW: Now go through again and see how many eigen values correspond to threshold
   for (int i = 0; i < eigen_values.GetNrows(); ++i) {
     // Get the relative size of the eigen value
     double sig = eigen_values(i)/sum;
     // Check against the threshold
     if (sig < eigen_threshold) {
       nKeptPCApars = i;
       break;
     }
   }
   NumParPCA = NumUnPCAdPars+nKeptPCApars;
   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);
   //DB Create array of correct size so eigen_values can be used in CorrelateSteps
   eigen_values_master = std::vector<double>(NumParPCA, 1.0);
   for (int i = FirstPCAdpar; i < FirstPCAdpar+nKeptPCApars; ++i) {eigen_values_master[i] = eigen_values(i-FirstPCAdpar);}
  
   // Now construct the transfer matrices
   //These matrices will be as big as number of unPCAd pars plus number of eigenvalues kept
   TransferMat.ResizeTo(CovMatrix->GetNrows(), NumParPCA);
   TransferMatT.ResizeTo(CovMatrix->GetNrows(), NumParPCA);
  
   // Get a subset of the eigen vector matrix
   TMatrixD temp(eigen_vectors.GetSub(0, eigen_vectors.GetNrows()-1, 0, nKeptPCApars-1));
  
   //Make transfer matrix which is two blocks of identity with a block of the PCA transfer matrix in between
   TMatrixD temp2;
   temp2.ResizeTo(CovMatrix->GetNrows(), NumParPCA);
  
   //First set the whole thing to 0
   temp2.Zero();
  
   //Set the first identity block for non-PCAed params before PCA block, before PCA XRows == YRows
   for(int iRow = 0; iRow < FirstPCAdpar; iRow++) {
     temp2(iRow, iRow) = 1.0;
   }
  
   //Set the transfer matrix block for the PCAd pars
   temp2.SetSub(FirstPCAdpar,FirstPCAdpar,temp);
  
   //Set the second identity block starting after PCA block, remember XRows != YRows.
   // XRows -> normal base, YRows, PCA base
   for(int iRow = 0;iRow < (CovMatrix->GetNrows()-1)-LastPCAdpar; iRow++) {
     const int OrigRow = LastPCAdpar + 1 + iRow;
     const int PCARow = FirstPCAdpar + nKeptPCApars + iRow;
     temp2(OrigRow, PCARow) = 1.;
   }
  
   TransferMat = temp2;
   // Copy the contents
   TransferMatT = TransferMat;
   // And then transpose
   TransferMatT.T();
  
   // Make the PCA parameter arrays
   _fParCurrPCA.ResizeTo(NumParPCA);
   _fParPropPCA.ResizeTo(NumParPCA);
   _fPreFitValuePCA.resize(NumParPCA);
  
   //KS: make easy map so we could easily find un-decomposed parameters
   _fErrorPCA.assign(NumParPCA, 1);
   isDecomposedPCA.assign(NumParPCA, -1);
   // First non PCA-ed block, since this is before PCA-ed block we don't need any mapping
   for (int i = 0; i < FirstPCAdpar; ++i) {
     isDecomposedPCA[i] = i;
   }
  
   // Second non-PCA-ed block, keep in mind this is in PCA base so we cannot use LastPCAdpar
   // we must shift them back into the original parameter index range.
   const int shift = _fNumPar - NumParPCA;
   for (int i = FirstPCAdpar + nKeptPCApars; i < NumParPCA; ++i) {
     isDecomposedPCA[i] = i + shift;
   }
  
   #ifdef MACH3_DEBUG_PCA
   //KS: Let's dump all useful matrices to properly validate PCA
   DebugPCA(sum, temp, CovMatrix->GetNrows());
   #endif
 }

◆ CorrelateSteps()

void PCAHandler::CorrelateSteps	(	const std::vector< double > &	IndivStepScale,
		const double	GlobalStepScale,
		const double *_restrict_	randParams,
		const double *_restrict_	corr_throw
	)

Use Cholesky throw matrix for better step proposal.

Definition at line 206 of file PCAHandler.cpp.

                                                                                 {
 // ************************************************
   // Throw around the current step
   #ifdef MULTITHREAD
   #pragma omp parallel for
   #endif
   for (int i = 0; i < NumParPCA; ++i)
   {
     if (_fErrorPCA[i] > 0.)
     {
       double IndStepScale = 1.;
       //KS: If undecomposed parameter apply individual step scale and Cholesky for better acceptance rate
       if(isDecomposedPCA[i] >= 0)
       {
         IndStepScale *= IndivStepScale[isDecomposedPCA[i]];
         IndStepScale *= corr_throw[isDecomposedPCA[i]];
       }
       //If decomposed apply only random number
       else
       {
         IndStepScale *= randParams[i];
         //KS: All PCA-ed parameters have the same step scale
         IndStepScale *= IndivStepScale[FirstPCAdpar];
       }
       _fParPropPCA(i) = _fParCurrPCA(i)+GlobalStepScale*IndStepScale*eigen_values_master[i];
     }
   }
   // Then update the parameter basis
   TransferToParam();
 }

◆ GetEigenValues()

const TVectorD PCAHandler::GetEigenValues ( ) const

inline

Get eigen values for all parameters, if you want for decomposed only parameters use GetEigenValuesMaster.

Definition at line 180 of file PCAHandler.h.

                                         {
     return eigen_values;
   }

◆ GetEigenValuesMaster()

const std::vector<double> PCAHandler::GetEigenValuesMaster ( ) const

inline

Get eigen value of only decomposed parameters, if you want for all parameters use GetEigenValues.

Definition at line 185 of file PCAHandler.h.

                                                        {
     return eigen_values_master;
   }

◆ GetEigenVectors()

const TMatrixD PCAHandler::GetEigenVectors ( ) const

inline

Get eigen vectors of covariance matrix, only works with PCA.

Definition at line 128 of file PCAHandler.h.

                                          {
     return eigen_vectors;
   }

◆ GetNumberPCAedParameters()

int PCAHandler::GetNumberPCAedParameters ( ) const

inline

Retrieve number of parameters in PCA base.

Definition at line 49 of file PCAHandler.h.

49 { return NumParPCA; }

◆ GetParCurrPCA()

double PCAHandler::GetParCurrPCA ( const int i ) const

inline

Get current parameter value using PCA.

Parameters

i	Parameter index

Definition at line 168 of file PCAHandler.h.

                                           {
     return _fParCurrPCA(i);
   }

◆ GetParPropPCA()

double PCAHandler::GetParPropPCA ( const int i ) const

inline

Get current parameter value using PCA.

Parameters

i	Parameter index

Definition at line 154 of file PCAHandler.h.

                                           {
     return _fParPropPCA(i);
   }

◆ GetPreFitValuePCA()

double PCAHandler::GetPreFitValuePCA ( const int i ) const

inline

Get current parameter value using PCA.

Parameters

i	Parameter index

Definition at line 161 of file PCAHandler.h.

                                               {
     return _fPreFitValuePCA[i];
   }

◆ GetTransferMatrix()

const TMatrixD PCAHandler::GetTransferMatrix ( ) const

inline

Get transfer matrix allowing to go from PCA base to normal base.

Definition at line 174 of file PCAHandler.h.

                                            {
     return TransferMat;
   }

◆ IsParameterDecomposed()

bool PCAHandler::IsParameterDecomposed ( const int i ) const

inline

Check if parameter in PCA base is decomposed or not.

Parameters

i	Parameter index

Definition at line 192 of file PCAHandler.h.

                                                 {
     if(isDecomposedPCA[i] >= 0) return false;
     else return true;
   }

◆ IsParameterFixedPCA()

bool PCAHandler::IsParameterFixedPCA ( const int i ) const

inline

Is parameter fixed in PCA base or not.

Parameters

i	Parameter index

Definition at line 121 of file PCAHandler.h.

                                               {
     if (_fErrorPCA[i] < 0) { return true;  }
     else                   { return false; }
   }

◆ Print()

void PCAHandler::Print ( ) const

KS: Print info about PCA parameters.

Definition at line 285 of file PCAHandler.cpp.

                              {
 // ********************************************
   MACH3LOG_INFO("PCA:");
   for (int i = 0; i < NumParPCA; ++i) {
     MACH3LOG_INFO("PCA {:<2} Current: {:<10.2f} Proposed: {:<10.2f}", i, _fParCurrPCA(i), _fParPropPCA(i));
   }
 }

◆ SanitisePCA()

void PCAHandler::SanitisePCA ( TMatrixDSym * CovMatrix )

private

KS: Make sure decomposed matrix isn't correlated with undecomposed.

Definition at line 162 of file PCAHandler.cpp.

                                                    {
 // ********************************************
   constexpr double correlation_threshold = 1e-6;
  
   bool found_significant_correlation = false;
  
   int N = CovMatrix->GetNrows();
   for (int i = FirstPCAdpar; i <= LastPCAdpar; ++i) {
     for (int j = 0; j < N; ++j) {
       // Skip if j is inside the decomposed range (we only want cross-correlations)
       if (j >= FirstPCAdpar && j <= LastPCAdpar) continue;
  
       double corr_val = (*CovMatrix)(i, j);
       if (std::fabs(corr_val) > correlation_threshold) {
         found_significant_correlation = true;
         MACH3LOG_ERROR("Significant correlation detected between decomposed parameter '{}' "
         "and undecomposed parameter '{}': {:.6e}", i, j, corr_val);
       }
     }
   }
  
   if (found_significant_correlation) {
     MACH3LOG_ERROR("There are correlations between undecomposed and decomposed part of matrices, this will not work");
     throw MaCh3Exception(__FILE__ , __LINE__);
   }
 }

◆ SetBranches()

void PCAHandler::SetBranches	(	TTree &	tree,
		const bool	SaveProposal,
		const std::vector< std::string > &	Names
	)

set branches for output file

Definition at line 294 of file PCAHandler.cpp.

                                                                                                     {
 // ********************************************
   for (int i = 0; i < NumParPCA; ++i) {
     tree.Branch(Form("%s_PCA", Names[i].c_str()), &_fParCurrPCA.GetMatrixArray()[i], Form("%s_PCA/D", Names[i].c_str()));
   }
  
   if(SaveProposal)
   {
     for (int i = 0; i < NumParPCA; ++i) {
       tree.Branch(Form("%s_PCA_Prop", Names[i].c_str()), &_fParPropPCA.GetMatrixArray()[i], Form("%s_PCA_Prop/D", Names[i].c_str()));
     }
   }
 }

◆ SetInitialParameters()

void PCAHandler::SetInitialParameters ( )

KS: Transfer the starting parameters to the PCA basis, you don't want to start with zero..

Parameters

IndStepScale Per parameter step scale, we set so each PCA param uses same step scale [eigen value takes care of size]

Definition at line 257 of file PCAHandler.cpp.

                                       {
 // ********************************************
   TransferToPCA();
   for (int i = 0; i < NumParPCA; ++i) {
     _fPreFitValuePCA[i] = _fParCurrPCA(i);
   }
 }

◆ SetParametersPCA()

void PCAHandler::SetParametersPCA ( const std::vector< double > & pars )

inline

Set values for PCA parameters in PCA base.

Parameters

pars	vector with new values of PCA params

Definition at line 105 of file PCAHandler.h.

                                                        {
     if (int(pars.size()) != NumParPCA) {
       MACH3LOG_ERROR("Parameter arrays of incompatible size! Not changing parameters! has size {} but was expecting {}", pars.size(), NumParPCA);
       throw MaCh3Exception(__FILE__ , __LINE__ );
     }
     int parsSize = int(pars.size());
     for (int i = 0; i < parsSize; i++) {
       _fParPropPCA(i) = pars[i];
     }
     //KS: Transfer to normal base
     TransferToParam();
   }

◆ SetParCurrPCA()

void PCAHandler::SetParCurrPCA	(	const int	i,
		const double	value
	)

inline

Set current value for parameter in PCA base.

Parameters

i	Parameter index
value	new value

Definition at line 145 of file PCAHandler.h.

                                                       {
     _fParCurrPCA(i) = value;
     // And then transfer back to the parameter basis
     TransferToParam();
   }

◆ SetParPropPCA()

void PCAHandler::SetParPropPCA	(	const int	i,
		const double	value
	)

inline

Set proposed value for parameter in PCA base.

Parameters

i	Parameter index
value	new value

Definition at line 136 of file PCAHandler.h.

                                                       {
     _fParPropPCA(i) = value;
     // And then transfer back to the parameter basis
     TransferToParam();
   }

◆ SetupPointers()

void PCAHandler::SetupPointers	(	std::vector< double > *	fCurr_Val,
		std::vector< M3::float_t > *	fProp_Val
	)

KS: Setup pointers to current and proposed parameter value which we need to convert them to PCA base each step.

Parameters

fCurr_Val	pointer to current position of parameter
fProp_Val	pointer to proposed position of parameter

Definition at line 29 of file PCAHandler.cpp.

                                                                 {
 // ********************************************
   _pCurrVal = fCurr_Val;
   _pPropVal = fProp_Val;
 }

◆ ThrowParameters()

void PCAHandler::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 can break Detailed Balance;.

Todo:: KS: We don't check if param is out of bounds. This is more problematic for PCA params.

Definition at line 333 of file PCAHandler.cpp.

                                                      {
 // ********************************************
   #pragma GCC diagnostic push
   #pragma GCC diagnostic ignored "-Wuseless-cast"
   //KS: Do not multithread!
   for (int i = 0; i < NumParPCA; ++i) {
     // Check if parameter is fixed first: if so don't randomly throw
     if (IsParameterFixedPCA(i)) continue;
  
     if(!IsParameterDecomposed(i))
     {
       (*_pPropVal)[i] = static_cast<M3::float_t>(fPreFitValue[i] + corr_throw[i]);
       int throws = 0;
       // Try again if we the initial parameter proposal falls outside of the range of the parameter
       while ((*_pPropVal)[i] > fUpBound[i] || (*_pPropVal)[i] < fLowBound[i]) {
         randParams[i] = random_number[M3::GetThreadIndex()]->Gaus(0, 1);
         const double corr_throw_single = M3::MatrixVectorMultiSingle(throwMatrixCholDecomp, randParams, _fNumPar, i);
         (*_pPropVal)[i] = static_cast<M3::float_t>(fPreFitValue[i] + corr_throw_single);
         if (throws > 10000)
         {
           //KS: Since we are multithreading there is danger that those messages
           //will be all over the place, small price to pay for faster code
           MACH3LOG_WARN("Tried {} times to throw parameter {} but failed", throws, i);
           MACH3LOG_WARN("Setting _fPropVal:  {} to {}", (*_pPropVal)[i], fPreFitValue[i]);
           MACH3LOG_WARN("I live at {}:{}", __FILE__, __LINE__);
          (*_pPropVal)[i] = static_cast<M3::float_t>(fPreFitValue[i]);
         }
         throws++;
       }
       (*_pCurrVal)[i] = (*_pPropVal)[i];
  
     } else {
       // KS: We have to multiply by number of parameters in PCA base
       SetParPropPCA(i, GetPreFitValuePCA(i) + randParams[i] * eigen_values_master[i] * (LastPCAdpar - FirstPCAdpar));
       SetParCurrPCA(i, GetParPropPCA(i));
     }
   } // end of parameter loop
  
   for (int i = 0; i < _fNumPar; ++i) {
     auto low = static_cast<M3::float_t>(fLowBound[i]);
     auto up  = static_cast<M3::float_t>(fUpBound[i]);
     (*_pPropVal)[i] = std::max(up, std::min((*_pPropVal)[i], low));
     (*_pCurrVal)[i] = (*_pPropVal)[i];
   }
   #pragma GCC diagnostic pop
 }

◆ ToggleFixAllParameters()

void PCAHandler::ToggleFixAllParameters ( const std::vector< std::string > & Names )

fix parameters at prior values

Definition at line 309 of file PCAHandler.cpp.

                                                                          {
 // ********************************************
   for (int i = 0; i < NumParPCA; i++) ToggleFixParameter(i, Names);
 }

◆ ToggleFixParameter()

void PCAHandler::ToggleFixParameter	(	const int	i,
		const std::vector< std::string > &	Names
	)

fix parameters at prior values

Definition at line 315 of file PCAHandler.cpp.

                                                                                   {
 // ********************************************
   int isDecom = -1;
   for (int im = 0; im < NumParPCA; ++im) {
     if(isDecomposedPCA[im] == i) {isDecom = im;}
   }
   if(isDecom < 0) {
     MACH3LOG_ERROR("Parameter {} is PCA decomposed can't fix this", Names[i]);
     //throw MaCh3Exception(__FILE__ , __LINE__ );
   } else {
     _fErrorPCA[isDecom] *= -1.0;
     if(IsParameterFixedPCA(i)) MACH3LOG_INFO("Setting un-decomposed {}(parameter {}/{} in PCA base) to fixed at {}", Names[i], i, isDecom, (*_pCurrVal)[i]);
     else MACH3LOG_INFO("Setting un-decomposed {}(parameter {}/{} in PCA base) free", Names[i], i, isDecom);
   }
 }

◆ TransferToParam()

void PCAHandler::TransferToParam ( )

Transfer param values from PCA base to normal base.

Definition at line 267 of file PCAHandler.cpp.

                                  {
 // ********************************************
   #pragma GCC diagnostic push
   #pragma GCC diagnostic ignored "-Wuseless-cast"
   // Make the temporary vectors
   TVectorD fParProp_vec = TransferMat*_fParPropPCA;
   TVectorD fParCurr_vec = TransferMat*_fParCurrPCA;
   #ifdef MULTITHREAD
   #pragma omp parallel for
   #endif
   for(int i = 0; i < static_cast<int>(_pCurrVal->size()); ++i) {
     (*_pPropVal)[i] = static_cast<M3::float_t>(fParProp_vec(i));
     (*_pCurrVal)[i] = fParCurr_vec(i);
   }
   #pragma GCC diagnostic pop
 }

◆ TransferToPCA()

void PCAHandler::TransferToPCA ( )

Transfer param values from normal base to PCA base.

Definition at line 242 of file PCAHandler.cpp.

                                {
 // ********************************************
   // Make the temporary vectors
   TVectorD fParCurr_vec(static_cast<Int_t>(_pCurrVal->size()));
   TVectorD fParProp_vec(static_cast<Int_t>(_pCurrVal->size()));
   for(int i = 0; i < static_cast<int>(_pCurrVal->size()); ++i) {
     fParCurr_vec(i) = (*_pCurrVal)[i];
     fParProp_vec(i) = (*_pPropVal)[i];
   }
  
   _fParCurrPCA = TransferMatT*fParCurr_vec;
   _fParPropPCA = TransferMatT*fParProp_vec;
 }

Member Data Documentation

◆ _fErrorPCA

std::vector<double> PCAHandler::_fErrorPCA

private

Tells if parameter is fixed in PCA base or not.

Definition at line 213 of file PCAHandler.h.

◆ _fParCurrPCA

TVectorD PCAHandler::_fParCurrPCA

private

CW: Proposed parameter value in PCA base.

Definition at line 211 of file PCAHandler.h.

◆ _fParPropPCA

TVectorD PCAHandler::_fParPropPCA

private

CW: Current parameter value in PCA base.

Definition at line 209 of file PCAHandler.h.

◆ _fPreFitValuePCA

std::vector<double> PCAHandler::_fPreFitValuePCA

private

Prefit value for PCA params.

Definition at line 207 of file PCAHandler.h.

◆ _pCurrVal

std::vector<double>* PCAHandler::_pCurrVal

private

Pointer to current value of the parameter.

Definition at line 240 of file PCAHandler.h.

◆ _pPropVal

std::vector<M3::float_t>* PCAHandler::_pPropVal

private

Pointer to proposed value of the parameter.

Definition at line 242 of file PCAHandler.h.

◆ eigen_threshold

double PCAHandler::eigen_threshold

private

CW: Threshold based on which we remove parameters in eigen base.

Definition at line 237 of file PCAHandler.h.

◆ eigen_values

TVectorD PCAHandler::eigen_values

private

Eigen value only of particles which are being decomposed.

Definition at line 224 of file PCAHandler.h.

◆ eigen_values_master

std::vector<double> PCAHandler::eigen_values_master

private

Eigen values which have dimension equal to _fNumParPCA, and can be used in CorrelateSteps.

Definition at line 228 of file PCAHandler.h.

◆ eigen_vectors

TMatrixD PCAHandler::eigen_vectors

private

Eigen vectors only of params which are being decomposed.

Definition at line 226 of file PCAHandler.h.

◆ FirstPCAdpar

int PCAHandler::FirstPCAdpar

private

Index of the first param that is being decomposed.

Definition at line 233 of file PCAHandler.h.

◆ isDecomposedPCA

std::vector<int> PCAHandler::isDecomposedPCA

private

If param is decomposed this will return -1, if not this will return enumerator to param in normal base. This way we can map stuff like step scale etc between normal base and undecomposed param in eigen base.

Definition at line 215 of file PCAHandler.h.

◆ LastPCAdpar

int PCAHandler::LastPCAdpar

private

Index of the last param that is being decomposed.

Definition at line 235 of file PCAHandler.h.

◆ nKeptPCApars

int PCAHandler::nKeptPCApars

private

Total number that remained after applying PCA Threshold.

Definition at line 231 of file PCAHandler.h.

◆ NumParPCA

int PCAHandler::NumParPCA

private

Number of parameters in PCA base.

Definition at line 217 of file PCAHandler.h.

◆ TransferMat

TMatrixD PCAHandler::TransferMat

private

Matrix used to converting from PCA base to normal base.

Definition at line 220 of file PCAHandler.h.

◆ TransferMatT

TMatrixD PCAHandler::TransferMatT

private

Matrix used to converting from normal base to PCA base.

Definition at line 222 of file PCAHandler.h.

The documentation for this class was generated from the following files:

Parameters/PCAHandler.h
Parameters/PCAHandler.cpp

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Introduction

Dimensionality Reduction

Partial Decomposition

Constructor & Destructor Documentation

◆ PCAHandler()

◆ ~PCAHandler()

Member Function Documentation

◆ AcceptStep()

◆ ConstructPCA()

◆ CorrelateSteps()

◆ GetEigenValues()

◆ GetEigenValuesMaster()

◆ GetEigenVectors()

◆ GetNumberPCAedParameters()

◆ GetParCurrPCA()

◆ GetParPropPCA()

◆ GetPreFitValuePCA()

◆ GetTransferMatrix()

◆ IsParameterDecomposed()

◆ IsParameterFixedPCA()

◆ Print()

◆ SanitisePCA()

◆ SetBranches()

◆ SetInitialParameters()

◆ SetParametersPCA()

◆ SetParCurrPCA()

◆ SetParPropPCA()

◆ SetupPointers()

◆ ThrowParameters()

◆ ToggleFixAllParameters()

◆ ToggleFixParameter()

◆ TransferToParam()

◆ TransferToPCA()

Member Data Documentation

◆ _fErrorPCA

◆ _fParCurrPCA

◆ _fParPropPCA

◆ _fPreFitValuePCA

◆ _pCurrVal

◆ _pPropVal

◆ eigen_threshold

◆ eigen_values

◆ eigen_values_master

◆ eigen_vectors

◆ FirstPCAdpar

◆ isDecomposedPCA

◆ LastPCAdpar

◆ nKeptPCApars

◆ NumParPCA

◆ TransferMat

◆ TransferMatT