Merge pull request #71 from mach3-software/feature_FixParamBasedOnName

Fix param using NAME

CMakeLists.txt

@@ -136,6 +136,7 @@ if (NOT DEFINED CMAKE_CXX_STANDARD OR "${CMAKE_CXX_STANDARD} " STREQUAL " ")
   SET(CMAKE_CXX_STANDARD 11)
 endif()
 
+# KS: If C++ standard is lower than C++ standar used for ROOT compilation things will go terribly wrong
 if(DEFINED ROOT_CXX_STANDARD AND ROOT_CXX_STANDARD GREATER CMAKE_CXX_STANDARD)
   set(CMAKE_CXX_STANDARD ${ROOT_CXX_STANDARD})
 endif()
@@ -167,7 +168,8 @@ target_compile_options(MaCh3CompilerOptions INTERFACE
                                                     -Wshadow
                                                     -Wuninitialized
                                                     -Wnon-virtual-dtor
-                                                    -Woverloaded-virtual)
+                                                    -Woverloaded-virtual
+                                                    )
 #KS: If Debug is not defined disable it by default
 if(NOT DEFINED MaCh3_DEBUG_ENABLED)
   SET(MaCh3_DEBUG_ENABLED FALSE)
@@ -190,9 +192,16 @@ if(MaCh3_DEBUG_ENABLED)
 else()
   #KS: Consider in future __attribute__((always_inline)) see https://indico.cern.ch/event/386232/sessions/159923/attachments/771039/1057534/always_inline_performance.pdf
   #https://gcc.gnu.org/onlinedocs/gcc-3.3.6/gcc/Optimize-Options.html
-  target_compile_options(MaCh3CompilerOptions INTERFACE -O3 -funroll-loops --param=max-vartrack-size=100000000 -finline-limit=100000000)
+  target_compile_options(MaCh3CompilerOptions INTERFACE
+                                                        -O3
+                                                        -funroll-loops
+                                                        --param=max-vartrack-size=100000000
+                                                        -finline-limit=100000000
+                                                        #-flto # FIXME need more testing
+                                                        )
   #KS: addomg Link-Time Optimization (LTO)
-  target_link_options(MaCh3CompilerOptions INTERFACE -flto)
+  # FIXME previously it wasn't used correctly but would need more testing
+  #target_link_libraries(MaCh3CompilerOptions INTERFACE -flto)
 endif()
 
 #KS: If multithreading is not defined enable it by default

Diagnostics/GetPenaltyTerm.cpp

@@ -21,6 +21,7 @@
 #include "TH1.h"
 #include "TColor.h"
 #include "TObjString.h"
+#include "TROOT.h"
 
 #ifdef MULTITHREAD
 #include "omp.h"
@@ -76,6 +77,11 @@ void GetPenaltyTerm(std::string inputFile, std::string configFile)
 
   ReadXSecFile(inputFile);
 
+  // KS: This can reduce time necessary for caching even by half
+  #ifdef MULTITHREAD
+  //ROOT::EnableImplicitMT();
+  #endif
+
   // Open the Chain
   TChain* Chain = new TChain("posteriors","");
   Chain->Add(inputFile.c_str()); 
@@ -145,7 +151,7 @@ void GetPenaltyTerm(std::string inputFile, std::string configFile)
     {
       isRelevantParam[i][j] = false;
 
-      //KS: Here we loop over all names and if parameters wasn't matched then we set it is relevant. For xsec it is easier to do it liek this
+      //KS: Here we loop over all names and if parameters wasn't matched then we set it is relevant. For xsec it is easier to do it like this
       if(Exclude[i])
       {
         bool found = false;
@@ -219,7 +225,7 @@ void GetPenaltyTerm(std::string inputFile, std::string configFile)
               //Check if parameter is relevant for this set
               if (isRelevantParam[k][i] && isRelevantParam[k][j])
               {
-                //KS: Since matrix is symmetric we can calcaute non diagonal elements only once and multiply by 2, can bring up to factor speed decrease.
+                //KS: Since matrix is symmetric we can calculate non diagonal elements only once and multiply by 2, can bring up to factor speed decrease.
                 int scale = 1;
                 if(i != j) scale = 2;
                 logL_private[k] += scale * 0.5*(fParProp[i] - nominal[i])*(fParProp[j] - nominal[j])*invCovMatrix[i][j];
@@ -234,7 +240,7 @@ void GetPenaltyTerm(std::string inputFile, std::string configFile)
         #pragma omp atomic
         logL[k] += logL_private[k];
       }
-      //Delete private arrayss
+      //Delete private arrays
       delete[] logL_private;
     }//End omp range
 
@@ -251,7 +257,7 @@ void GetPenaltyTerm(std::string inputFile, std::string configFile)
             //Check if parameter is relevant for this set
             if (isRelevantParam[k][i] && isRelevantParam[k][j])
             {
-              //KS: Since matrix is symetric we can calcaute non daigonal elements only once and multiply by 2, can bring up to factor speed decrease.
+              //KS: Since matrix is symmetric we can calculate non diagonal elements only once and multiply by 2, can bring up to factor speed decrease.
               int scale = 1;
               if(i != j) scale = 2;
               logL[k] += scale * 0.5*(fParProp[i] - nominal[i])*(fParProp[j] - nominal[j])*invCovMatrix[i][j];
@@ -348,7 +354,6 @@ void ReadXSecFile(std::string inputFile)
     }
   }
 
-
   size = XSecMatrix->GetNrows();
 
   auto systematics = XSecFile["Systematics"];

Diagnostics/RHat.cpp

@@ -21,13 +21,15 @@
 #include "TStopwatch.h"
 #include "TColor.h"
 #include "TStyle.h"
+#include "TROOT.h"
 
 #ifdef MULTITHREAD
 #include "omp.h"
 #endif
 
 #include "manager/manager.h"
 
+
 //KS: This exe is meant to calculate R hat estimator. For a well converged this distribution should be centred at one. 
 //Based on Gelman et. al. arXiv:1903.08008v5
 
@@ -179,6 +181,11 @@ void PrepareChains() {
   Draws = new double**[Nchains]();
   DrawsFolded = new double**[Nchains]();
 
+  // KS: This can reduce time necessary for caching even by half
+  #ifdef MULTITHREAD
+  //ROOT::EnableImplicitMT();
+  #endif
+
   // Open the Chain
   //It is tempting to multithread here but unfortunately, ROOT files are not thread safe :(
   for (int m = 0; m < Nchains; m++)

Doc/mainpage.md

@@ -4,3 +4,6 @@
 Welcome to %MaCh3!
 
 This is the Reference Guide of the MaCh3 software.
+
+You can find additional documentation on our [Wiki](https://github.com/mach3-software/MaCh3/wiki)
+

README.md

@@ -150,7 +150,7 @@ This is an example how your executable can look like using MaCh3:
 - [Wiki](https://github.com/mach3-software/MaCh3/wiki)
 - [Mailing lists](https://www.jiscmail.ac.uk/cgi-bin/webadmin?A0=MACH3)
 - [Slack](https://t2k-experiment.slack.com/archives/C06EM0C6D7W/p1705599931356889)
-
+- [Discussions](https://github.com/mach3-software/MaCh3/discussions)
 
 ### Plotting and Diagnostic
 Example of chain diagnostic utils can be found [here](https://github.com/mach3-software/MaCh3/tree/develop/Diagnostics) with example of config.

covariance/covarianceBase.cpp

@@ -10,7 +10,7 @@ covarianceBase::covarianceBase(const char *name, const char *file) : inputFile(s
   LastPCAdpar = -999;
 }
 // ********************************************
-covarianceBase::covarianceBase(std::vector<std::string> YAMLFile, const char *name, double threshold, int FirstPCA, int LastPCA) : inputFile(YAMLFile[0].c_str()), matrixName(name), pca(true), eigen_threshold(threshold), FirstPCAdpar(FirstPCA), LastPCAdpar(LastPCA) {
+covarianceBase::covarianceBase(const std::vector<std::string>& YAMLFile, const char *name, double threshold, int FirstPCA, int LastPCA) : inputFile(YAMLFile[0].c_str()), matrixName(name), pca(true), eigen_threshold(threshold), FirstPCAdpar(FirstPCA), LastPCAdpar(LastPCA) {
 // ********************************************
 
   MACH3LOG_INFO("Constructing instance of covarianceBase using ");
@@ -50,8 +50,8 @@ covarianceBase::covarianceBase(const char *name, const char *file, int seed) : i
   init(name, file);
   FirstPCAdpar = -999;
   LastPCAdpar = -999;
-
 }
+
 // ********************************************
 covarianceBase::covarianceBase(const char *name, const char *file, int seed, double threshold, int firstpcapar, int lastpcapar) : inputFile(std::string(file)), pca(true), eigen_threshold(threshold), FirstPCAdpar(firstpcapar), LastPCAdpar(lastpcapar) {
 // ********************************************
@@ -325,7 +325,7 @@ void covarianceBase::init(const char *name, const char *file)
 // An init function for the YAML constructor
 // All you really need from the YAML file is the number of Systematics
 // Then get all the info from the YAML file in the covarianceXsec::ParseYAML function
-void covarianceBase::init(std::vector<std::string> YAMLFile)
+void covarianceBase::init(const std::vector<std::string>& YAMLFile)
 {
   _fYAMLDoc["Systematics"] = YAML::Node(YAML::NodeType::Sequence);
   for(unsigned int i = 0; i < YAMLFile.size(); i++)
@@ -417,7 +417,7 @@ void covarianceBase::init(std::vector<std::string> YAMLFile)
   //This makes the root TCov from YAML
   for(int j = 0; j < _fNumPar;j++) {
     (*_fCovMatrix)(j, j)=_fError[j]*_fError[j];
-    //Get the map of parameter name to correlation fomr the Correlations object
+    //Get the map of parameter name to correlation from the Correlations object
     for (auto const& [key, val] : Correlations[j]) {
       int index = -1;
 
@@ -555,7 +555,6 @@ void covarianceBase::ReserveMemory(const int SizeVec) {
 // Set all the covariance matrix parameters to a user-defined value
 // Might want to split this
 void covarianceBase::setPar(int i , double val) {
-
   MACH3LOG_INFO("Over-riding {}: ", GetParName(i));
   MACH3LOG_INFO("_fPropVal ({}), _fCurrVal ({}), _fPreFitValue ({}) to ({})", _fPropVal[i], _fCurrVal[i], _fPreFitValue[i], val);
 
@@ -614,7 +613,6 @@ const std::vector<double> covarianceBase::getProposed() const {
 
 // Throw nominal values
 void covarianceBase::throwNominal(bool nomValues, int seed) {
-
   TVectorD* vec = new TVectorD(_fNumPar);
   for (int i = 0; i < _fNumPar; i++) {
     (*vec)(i) = 1.0;
@@ -654,7 +652,6 @@ void covarianceBase::throwNominal(bool nomValues, int seed) {
       }
     }
   } else {
-
     // If we want nominal values, set all entries to 1 (defined as nominal in MaCh3)
     for (int i = 0; i < int(nominal.size()); i++) {
       nominal[i] = 1.0;
@@ -758,7 +755,6 @@ void covarianceBase::RandomConfiguration() {
     _fCurrVal[i] = _fPropVal[i];
   }
   if (pca) TransferToPCA();
-
 }
 
 // *************************************
@@ -1047,9 +1043,9 @@ double covarianceBase::GetLikelihood() {
   return CalcLikelihood();
 }
 
-
+// ********************************************
 void covarianceBase::printPars() {
-
+// ********************************************
   MACH3LOG_INFO("Number of pars: {}", _fNumPar);
   MACH3LOG_INFO("Current {} parameters:", matrixName);
   for(int i = 0; i < _fNumPar; i++) {
@@ -1059,9 +1055,10 @@ void covarianceBase::printPars() {
   return;
 }
 
+// ********************************************
 // Sets the proposed parameters to the nominal values
 void covarianceBase::setParameters(std::vector<double> pars) {
-
+// ********************************************
   // If empty, set the proposed to nominal
   if (pars.empty()) {
     // For xsec this means setting to the prior (because nominal is the prior)
@@ -1070,7 +1067,6 @@ void covarianceBase::setParameters(std::vector<double> pars) {
     }
     // If not empty, set the parameters to the specified
   } else {
-
 	if (pars.size() != size_t(_fNumPar)) {
       MACH3LOG_ERROR("Warning: parameter arrays of incompatible size! Not changing parameters! {} has size {} but was expecting {}", matrixName, pars.size(), _fNumPar);
       throw;
@@ -1154,30 +1150,55 @@ void covarianceBase::toggleFixAllParameters() {
 void covarianceBase::toggleFixParameter(const int i) {
 // ********************************************
   if(!pca) {
-	if (i > _fNumPar) {
+    if (i > _fNumPar) {
       MACH3LOG_ERROR("Can't toggleFixParameter for parameter {} because size of covariance ={}", i, _fNumPar);
       MACH3LOG_ERROR("Fix this in your config file please!");
-	  throw;
-	} else {
-	  _fError[i] *= -1.0;
-      MACH3LOG_INFO("Setting {}(parameter {}) to fixed at {}", GetParName(i), i, _fCurrVal[i]);
-	} 
+      throw;
+    } else {
+      _fError[i] *= -1.0;
+      MACH3LOG_INFO("Setting {}(parameter {}) to fixed at {}", GetParFancyName(i), i, _fCurrVal[i]);
+    }
   } else {
-	int isDecom = -1;
-	for (int im = 0; im < _fNumParPCA; ++im) {
-	  if(isDecomposed_PCA[im] == i) {isDecom = im;}
-	}
-	if(isDecom < 0) {
-      MACH3LOG_ERROR("Setting {}(parameter {}) to fixed at {}", GetParName(i), i, _fCurrVal[i]);
-	  //throw; 
-	} else {
-	  fParSigma_PCA[isDecom] *= -1.0;
+    int isDecom = -1;
+    for (int im = 0; im < _fNumParPCA; ++im) {
+      if(isDecomposed_PCA[im] == i) {isDecom = im;}
+    }
+    if(isDecom < 0) {
+      MACH3LOG_ERROR("Parameter {} is PCA decomposed can't fix this", GetParName(i));
+      //throw;
+    } else {
+      fParSigma_PCA[isDecom] *= -1.0;
       MACH3LOG_INFO("Setting un-decomposed {}(parameter {}/{} in PCA base) to fixed at {}", GetParName(i), i, isDecom, _fCurrVal[i]);
-	}
+    }
   }
 
   return;
 }
+
+// ********************************************
+void covarianceBase::toggleFixParameter(const std::string& name) {
+// ********************************************
+  for (int i = 0; i <_fNumPar; ++i) {
+    if(name == _fFancyNames[i]) {
+      toggleFixParameter(i);
+      return;
+    }
+  }
+  MACH3LOG_WARN("I couldn't find parameter with name {}, therefore will not fix it", name);
+}
+
+// ********************************************
+bool covarianceBase::isParameterFixed(const std::string& name) {
+// ********************************************
+  for (int i = 0; i <_fNumPar; ++i) {
+    if(name == _fFancyNames[i]) {
+      return isParameterFixed(i);
+    }
+  }
+  MACH3LOG_WARN("I couldn't find parameter with name {}, therefore will not fix it", name);
+  return false;
+}
+
 // ********************************************
 void covarianceBase::setFlatPrior(const int i, const bool eL) {
 // ********************************************
@@ -1186,7 +1207,6 @@ void covarianceBase::setFlatPrior(const int i, const bool eL) {
     MACH3LOG_ERROR("Fix this in your config file please!");
     throw;
   } else {
-
     if(eL){
       MACH3LOG_INFO("Setting {} (parameter {}) to flat prior", GetParName(i), i);
     }
@@ -1202,7 +1222,6 @@ void covarianceBase::setFlatPrior(const int i, const bool eL) {
 //KS: Custom function to perform multiplication of matrix and vector with multithreading
 void covarianceBase::MatrixVectorMulti(double* _restrict_ VecMulti, double** _restrict_ matrix, const double* _restrict_ vector, const int n) {
 // ********************************************
-
   #ifdef MULTITHREAD
   #pragma omp parallel for
   #endif