ParameterHandlerBase.h

Go to the documentation of this file.

#pragma once
 
// MaCh3 includes
#include "Manager/Manager.h"
#include "Parameters/ParameterHandlerUtils.h"
#include "Parameters/AdaptiveMCMCHandler.h"
#include "Parameters/PCAHandler.h"
#include "Parameters/ParameterTunes.h"
 
class ParameterHandlerBase {
 public:
  ParameterHandlerBase(const std::vector<std::string>& YAMLFile, std::string name, double threshold = -1, int FirstPCAdpar = -999, int LastPCAdpar = -999);
  ParameterHandlerBase(std::string name, std::string file, double threshold = -1, int FirstPCAdpar = -999, int LastPCAdpar = -999);
 
  virtual ~ParameterHandlerBase();
  
 
 
  // ETA - maybe need to add checks to index on the setters? i.e. if( i > _fPropVal.size()){throw;}
  void SetCovMatrix(TMatrixDSym *cov);
  void SetName(const std::string& name) { matrixName = name; }
  void SetParName(const int i, const std::string& name) { _fNames.at(i) = name; }
  void SetSingleParameter(const int parNo, const double parVal);
  void SetPar(const int i, const double val);
  void SetParCurrProp(const int i, const double val);
  void SetParProp(const int i, const double val) {
    _fPropVal[i] = val;
    if (pca) PCAObj->TransferToPCA();
  }
  void SetParameters(const std::vector<double>& pars = {});
  void SetFlatPrior(const int i, const bool eL);
  
  void SetRandomThrow(const int i, const double rand) { randParams[i] = rand;}
  double GetRandomThrow(const int i) const { return randParams[i];}
 
  void SetBranches(TTree &tree, const bool SaveProposal = false);
  void SetStepScale(const double scale);
  void SetIndivStepScale(const int ParameterIndex, const double StepScale){ _fIndivStepScale.at(ParameterIndex) = StepScale; }
  void SetIndivStepScale(const std::vector<double>& stepscale);
  inline void SetPrintLength(const unsigned int PriLen) { PrintLength = PriLen; }
 
  void SaveUpdatedMatrixConfig();
 
  void ThrowParProp(const double mag = 1.);
 
  void ThrowParCurr(const double mag = 1.);
  void ThrowParameters();
  void RandomConfiguration();
  
  int CheckBounds() const _noexcept_;
  double CalcLikelihood() const _noexcept_;
  virtual double GetLikelihood();
 
  TMatrixDSym *GetCovMatrix() const { return covMatrix; }
  TMatrixDSym *GetInvCovMatrix() const { return invCovMatrix; }
  double GetInvCovMatrix(const int i, const int j) const { return InvertCovMatrix[i][j]; }
 
  double GetCorrThrows(const int i) const { return corr_throw[i]; }
 
  inline bool GetFlatPrior(const int i) const { return _fFlatPrior[i]; }
 
  std::string GetName() const { return matrixName; }
  std::string GetParName(const int i) const {return _fNames[i];}
 
  int GetParIndex(const std::string& name) const;
 
  std::string GetParFancyName(const int i) const {return _fFancyNames[i];}
  std::string GetInputFile() const { return inputFile; }
 
  inline double GetDiagonalError(const int i) const { return std::sqrt((*covMatrix)(i,i)); }
  inline double GetError(const int i) const {return _fError[i];}
 
  void ResetIndivStepScale();
 
  void InitialiseAdaption(const YAML::Node& adapt_manager);
  void SaveAdaptiveToFile(const std::string& outFileName, const std::string& systematicName) {
    AdaptiveHandler->SaveAdaptiveToFile(outFileName, systematicName); }
 
  bool GetDoAdaption() const {return use_adaptive;}
  void SetThrowMatrix(TMatrixDSym *cov);
  void UpdateThrowMatrix(TMatrixDSym *cov);
  inline void SetNumberOfSteps(const int nsteps) {
    AdaptiveHandler->SetTotalSteps(nsteps);
    if(AdaptiveHandler->AdaptionUpdate()) ResetIndivStepScale();
  }
 
  inline TMatrixDSym *GetThrowMatrix() const {return throwMatrix;}
  double GetThrowMatrix(const int i, const int j) const { return throwMatrixCholDecomp[i][j];}
  inline TMatrixD *GetThrowMatrix_CholDecomp() const {return throwMatrix_CholDecomp;}
 
  TH2D* GetCorrelationMatrix();
 
  inline const double* RetPointer(const int iParam) {return &(_fPropVal.data()[iParam]);}
 
  inline const std::vector<double> &GetParPropVec() {return _fPropVal;}
 
  //Some Getters
  inline int  GetNumParams() const {return _fNumPar;}
  std::vector<double> GetPreFitValues() const {return _fPreFitValue;}
  std::vector<double> GetProposed() const;
  inline double GetParProp(const int i) const { return _fPropVal[i]; }
  inline double GetParCurr(const int i) const { return _fCurrVal[i]; }
  inline double GetParInit(const int i) const { return _fPreFitValue[i]; }
  inline double GetUpperBound(const int i) const { return _fUpBound[i];}
  inline double GetLowerBound(const int i) const { return _fLowBound[i]; }
  inline double GetIndivStepScale(const int ParameterIndex) const {return _fIndivStepScale.at(ParameterIndex); }
  inline double GetGlobalStepScale() const {return _fGlobalStepScale; }
 
  inline int GetNParameters() const {
    if (pca) return PCAObj->GetNumberPCAedParameters();
    else return _fNumPar;
  }
 
  void PrintNominal() const;
  void PrintNominalCurrProp() const;
  void PrintParameters() const {PrintNominalCurrProp();};
  void PrintIndivStepScale() const;
 
  virtual void ProposeStep();
  void Randomize() _noexcept_;
  void CorrelateSteps() _noexcept_;
  void UpdateAdaptiveCovariance();
 
  void AcceptStep() _noexcept_;
 
  void ToggleFixAllParameters();
  void ToggleFixParameter(const int i);
  void ToggleFixParameter(const std::string& name);
 
  bool IsParameterFixed(const int i) const {
    if (_fError[i] < 0) { return true; }
    else                { return false; }
  }
  bool IsParameterFixed(const std::string& name) const;
 
  void ConstructPCA(const double eigen_threshold, int FirstPCAdpar, int LastPCAdpar);
 
  inline bool IsPCA() const { return pca; }
 
  YAML::Node GetConfig() const { return _fYAMLDoc; }
 
  inline adaptive_mcmc::AdaptiveMCMCHandler* GetAdaptiveHandler() const  {
    if (!use_adaptive) {
      MACH3LOG_ERROR("Am not running in Adaptive mode");
      throw MaCh3Exception(__FILE__ , __LINE__ );
    }
    return AdaptiveHandler.get();
  }
 
  void SetTune(const std::string& TuneName);
  
  inline PCAHandler* GetPCAHandler() const {
    if (!pca) {
      MACH3LOG_ERROR("Am not running in PCA mode");
      throw MaCh3Exception(__FILE__ , __LINE__ );
    }
    return PCAObj.get();
  }
 
protected:
  void Init(const std::string& name, const std::string& file);
  void Init(const std::vector<std::string>& YAMLFile);
  void ReserveMemory(const int size);
 
  void MakePosDef(TMatrixDSym *cov = nullptr);
 
  void MakeClosestPosDef(TMatrixDSym *cov);
 
  void SetThrowMatrixFromFile(const std::string& matrix_file_name, const std::string& matrix_name, const std::string& means_name);
 
  bool AppliesToSample(const int SystIndex, const std::string& SampleName) const;
 
  void FlipParameterValue(const int index, const double FlipPoint);
 
  void CircularParBounds(const int i, const double LowBound, const double UpBound);
 
  void EnableSpecialProposal(const YAML::Node& param, const int Index);
 
  void SpecialStepProposal();
 
  bool doSpecialStepProposal;
 
  const std::string inputFile;
 
  std::string matrixName;
  TMatrixDSym *covMatrix;
  TMatrixDSym *invCovMatrix;
  std::vector<std::vector<double>> InvertCovMatrix;
    
  std::vector<std::unique_ptr<TRandom3>> random_number;
 
  double* randParams;
  double* corr_throw;
  double _fGlobalStepScale;
 
  int PrintLength;
 
  std::vector<std::string> _fNames;
  std::vector<std::string> _fFancyNames;
  YAML::Node _fYAMLDoc;
  int _fNumPar;
  std::vector<double> _fPreFitValue;
  std::vector<double> _fCurrVal;
  std::vector<double> _fPropVal;
  std::vector<double> _fError;
  std::vector<double> _fLowBound;
  std::vector<double> _fUpBound;
  std::vector<double> _fIndivStepScale;
  std::vector<bool> _fFlatPrior;
  std::vector<std::vector<std::string>> _fSampleNames;
 
  TMatrixDSym* throwMatrix;
  TMatrixD* throwMatrix_CholDecomp;
  double** throwMatrixCholDecomp;
 
  bool pca;
  bool use_adaptive;
 
  std::unique_ptr<PCAHandler> PCAObj;
  std::unique_ptr<adaptive_mcmc::AdaptiveMCMCHandler> AdaptiveHandler;
  std::unique_ptr<ParameterTunes> Tunes;
 
  std::vector<int>    FlipParameterIndex;
  std::vector<double> FlipParameterPoint;
  std::vector<int>    CircularBoundsIndex;
  std::vector<std::pair<double,double>> CircularBoundsValues;
};