fix merge conflcit

.github/workflows/CIPythonValidations.yml

@@ -16,7 +16,7 @@ jobs:
     runs-on: ubuntu-latest
 
     container:
-      image: ghcr.io/mach3-software/mach3:alma9latest
+      image: ghcr.io/mach3-software/mach3:alma9v1.3.0
 
     steps:
     - uses: actions/checkout@v4

.github/workflows/CIValidations.yml

@@ -42,7 +42,7 @@ jobs:
             test_3: empty
 
     container:
-      image: ghcr.io/mach3-software/mach3:alma9latest
+      image: ghcr.io/mach3-software/mach3:alma9v1.3.0
 
     steps:
     - uses: actions/checkout@v4

.github/workflows/CodeQL.yml

@@ -11,7 +11,7 @@ jobs:
     name: Analyze C++ Code with CodeQL
     runs-on: ubuntu-latest
     container:
-      image: ghcr.io/mach3-software/mach3:alma9latest
+      image: ghcr.io/mach3-software/mach3:alma9v1.3.0
     permissions:
       security-events: write
       packages: read

CIValidations/Benchmark/CMakeLists.txt

@@ -0,0 +1,8 @@
+foreach(TEST
+    MaCh3Benchmark
+    )
+
+    add_executable(${TEST} ${TEST}.cpp)
+    target_link_libraries(${TEST} PRIVATE Catch2::Catch2WithMain MaCh3::All MaCh3Tutorial::samplePDFTutorial)
+    catch_discover_tests(${TEST})
+endforeach()

CIValidations/Benchmark/MaCh3Benchmark.cpp

@@ -0,0 +1,38 @@
+// MaCh3 spline includes
+#include "mcmc/MaCh3Factory.h"
+#include "samplePDF/samplePDFTutorial.h"
+
+#include "catch2/catch_test_macros.hpp"
+#include "catch2/benchmark/catch_benchmark.hpp"
+
+TEST_CASE("Benchmark MaCh3") {
+  // Initialise manger responsible for config handling
+  auto FitManager = std::make_unique<manager>("Inputs/FitterConfig.yaml");
+
+  // Initialise covariance class reasonable for Systematics
+  covarianceXsec* xsec = MaCh3CovarianceFactory(FitManager.get(), "Xsec");
+  covarianceOsc*  osc  = MaCh3CovarianceFactory<covarianceOsc>(FitManager.get(), "Osc");
+
+  // Initialise samplePDF
+  auto SampleConfig = FitManager->raw()["General"]["TutorialSamples"].as<std::vector<std::string>>();
+  auto mySamples = MaCh3SamplePDFFactory<samplePDFTutorial>(SampleConfig, xsec, osc);
+
+  // Create MCMC Class
+  std::unique_ptr<FitterBase> MaCh3Fitter = MaCh3FitterFactory(FitManager.get());
+  // Add covariance to MCM
+  MaCh3Fitter->addSystObj(xsec);
+  MaCh3Fitter->addSystObj(osc);
+  for (size_t i = 0; i < SampleConfig.size(); ++i) {
+    MaCh3Fitter->addSamplePDF(mySamples[i]);
+  }
+  // Benchmark
+  BENCHMARK("MaCh3Fitter::DragRace") {
+    MaCh3Fitter->DragRace(1);
+  };
+
+  delete xsec;
+  delete osc;
+  for (size_t i = 0; i < SampleConfig.size(); ++i) {
+    delete mySamples[i];
+  }
+}

CIValidations/CMakeLists.txt

@@ -39,3 +39,7 @@ install(DIRECTORY
 if(MaCh3Tutorial_UNITTESTS_ENABLED)
   add_subdirectory(UnitTests)
 endif()
+
+if(MaCh3Tutorial_Benchmark_ENABLED)
+  add_subdirectory(Benchmark)
+endif()

CIValidations/CovarianceValidations.cpp

@@ -17,7 +17,7 @@ int main(int argc, char *argv[])
   std::vector<std::string> xsecCovMatrixFile = {"Inputs/SystematicModel.yaml"};
   auto xsec = std::make_unique<covarianceXsec>(xsecCovMatrixFile, "xsec_cov");
 
-  std::vector<double> ParProp = {1.05, 0.90, 1.10, 1.05, 0.25, 1.70, 3.20, -1.10, -1.70};
+  std::vector<double> ParProp = {1.05, 0.90, 1.10, 1.05, 1.05, 1.05, 1.05, 1.05, 1.70, 3.20, -1.10, -1.70};
   xsec->setParameters(ParProp);
   xsec->printNominalCurrProp();
 
@@ -40,18 +40,23 @@ int main(int argc, char *argv[])
   xsec->setParameters(ParProp);
   xsec->acceptStep();
   ///// Test Params from DetId /////
-  const int Det_Id = 1;
-  for (int s = 0; s < kSystTypes; ++s)
+  const std::vector<int> Det_Id = {1, 984, 985};
+  for (int id = 0; id < kSystTypes; ++id)
   {
-    int NParams = xsec->GetNumParamsFromDetID(Det_Id, static_cast<SystType>(s));
-    outFile << "Found " << NParams << " Params of type " << SystType_ToString(static_cast<SystType>(s)) << std::endl;
+    for (int s = 0; s < kSystTypes; ++s)
+    {
+      int NParams = xsec->GetNumParamsFromDetID(Det_Id[id], static_cast<SystType>(s));
+      outFile << "Found " << NParams << " for "<< Det_Id[id] << " Params of type " << SystType_ToString(static_cast<SystType>(s)) << std::endl;
+
+      // These return vector and I am too lazy to make check so here let just run them to so if things don't break
+      xsec->GetSystIndexFromDetID(Det_Id[id], static_cast<SystType>(s));
+      xsec->GetParsIndexFromDetID(Det_Id[id], static_cast<SystType>(s));
+      xsec->GetParsNamesFromDetID(Det_Id[id], static_cast<SystType>(s));
+    }
 
-    // These return vector and I am too lazy to make check so here let just run them to so if things don't break
-    xsec->GetSystIndexFromDetID(Det_Id, static_cast<SystType>(s));
-    xsec->GetParsIndexFromDetID(Det_Id, static_cast<SystType>(s));
-    xsec->GetParsNamesFromDetID(Det_Id, static_cast<SystType>(s));
+    auto Norm = xsec->GetNormParsFromDetID(Det_Id[id]);
+    outFile << "Found " << Norm.size()   << " for "<< Det_Id[id] << " From GetNormParsFromDetID" << std::endl;
   }
-
   for (int i = 0; i < xsec->GetNumParams(); ++i) {
     for (int j = 0; j < xsec->GetNumParams(); ++j) {
       outFile << "Inv Cov Matrix for " << i << " and " << j << " is equal to=" << xsec->GetInvCovMatrix(i, j) << std::endl;
@@ -143,9 +148,20 @@ std::string yamlContent = R"(
   //KS: Let's make Doctor Wallace proud
   Adapt->initialiseAdaption(AdaptSetting);
 
-  std::vector<double> ParAdapt = {1.05, 0.90, 1.10, 1.05, 0.25, 1.70, 3.20, -1.10, -1.70, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 };
+  std::vector<double> ParAdapt = {1.05, 0.90, 1.10, 1.05, 1.05, 1.05, 1.05, 1.05, 1.70, 3.20, -1.10, -1.70, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 };
   Adapt->setParameters(ParAdapt);
+  bool increase = true;
   for(int i = 0; i < 50000; ++i ) {
+
+    // Determine the direction of adjustment
+    if (i % 20 == 10) { // Switch direction every 10 iterations
+      increase = !increase;
+    }
+    // Adjust parameters
+    for (double& param : ParAdapt) {
+      param += (increase ? 0.001 : -0.001);
+    }
+    Adapt->setParameters(ParAdapt);
     Adapt->acceptStep();
   }
   Adapt->saveAdaptiveToFile("Wacky.root", "xsec");

CIValidations/LLHValidation.cpp

@@ -20,15 +20,15 @@ int main(int argc, char *argv[])
   }
 
   TFile* file = new TFile("LLH_Test.root", "UPDATE");
-  std::vector<std::string> Names = {"Norm_Param_0", "Norm_Param_1", "Norm_Param_2"};
+  std::vector<std::string> Names = {"Norm_Param_0", "Norm_Param_1", "Norm_Param_2", "BinnedSplineParam1", "BinnedSplineParam2", "BinnedSplineParam3", "BinnedSplineParam4", "BinnedSplineParam5"};
   std::ofstream outFile("LLH.txt");
   for (size_t i = 0; i < Names.size(); ++i) {
     std::string histPath = "Sample_LLH/" + Names[i] + "_sam";
     TH1D* hist = static_cast<TH1D*>(file->Get(histPath.c_str()));
 
     for (int j = 1; j <= hist->GetNbinsX(); ++j) {
       double binContent = hist->GetBinContent(j);
-      outFile << "Dial " << Names[i] << " bin = " << j << " content = " << binContent << std::endl;
+      outFile << "Dial " << Names[i] << " bin = " << j << " content = " << std::fixed << std::setprecision(6) << binContent << std::endl;
     }
   }
   outFile.close();

CIValidations/MaCh3CLI

@@ -1,11 +1,16 @@
 #!/bin/bash
+#set -x  # Enable command tracing
+set -e  # Exit on error
 
 # Check if no arguments are provided
 if [ "$#" -eq 0 ]; then
   echo "[ERROR]: No argument specified, please use --help"
   exit 1
 fi
 
+mach3tutorial-config --has-feature UNITTESTS || true
+HAVE_UNITTEST_SUPPORT=${?}
+
 # Display help message
 show_help() {
   echo "Arguments:"
@@ -34,13 +39,17 @@ run_fitter_validations() {
   ${MaCh3_ROOT}/bin/MatrixPlotter bin/TutorialDiagConfig.yaml MCMC_Test_drawCorr.root
   ${MaCh3_ROOT}/bin/GetPostfitParamPlots MCMC_Test_drawCorr.root
   ${MaCh3_ROOT}/bin/PlotLLH MCMC_Test.root
+  ${MaCh3_ROOT}/bin/PlotMCMCDiag MCMC_Test_MCMC_Diag.root MCMC_Test_MCMC_Diag.root
 }
 
 run_covariance_validations() {
   ./CIValidations/CovarianceValidations
   ./CIValidations/MaCh3ModeValidations
   ./CIValidations/pValueValidations
-  ./CIValidations/UnitTests/manager_tests
+  if [ "IF[HAVE_UNITTEST_SUPPORT]" == "0" ]; then
+    ./CIValidations/UnitTests/manager_tests
+    ./CIValidations/UnitTests/histogram_tests
+  fi
 }
 
 run_nu_mcmc_tools_validations() {

CIValidations/NuOscillatorInterfaceValidations.cpp

@@ -1,3 +1,4 @@
+#include <Constants/OscillatorConstants.h>
 #include "Utils/Comparison.h"
 
 #include "Oscillator/OscillatorFactory.h"
@@ -25,7 +26,7 @@ int main() {
   auto osc = OscillatorFactory().CreateOscillator(
       "NuOscillatorInterfaceValidations.yml");
 
-  std::vector<double> Energies;
+  std::vector<FLOAT_T> Energies;
   double start = -3;
   double end = 3;
   size_t nbins = 100;
@@ -37,7 +38,7 @@ int main() {
   osc->SetEnergyArrayInCalcer(Energies);
   osc->Setup();
 
-  std::vector<double> OscParProp = {0.3,      0.5, 0.020, 7.53e-5,
+  std::vector<FLOAT_T> OscParProp = {0.3,      0.5, 0.020, 7.53e-5,
                                     2.494e-3, 0.0, 1300,  2.6, 0.5};
   osc->CalculateProbabilities(OscParProp);
 

CIValidations/README.md

@@ -43,3 +43,22 @@
 * Dump and compare event weights
 * Dump and compare weight for each event
 * Compare LLH using different test statistics
+
+## UnitTests
+
+### manager_tests
+* Test for YamlHelpers
+
+### histogram_tests
+* Test fo Histogram Utils and TH2Poly related stuff
+
+## Python Tests
+
+### test_fitter_sample_pdf
+* Fit with sample pdf
+
+### test_fitter_stat_only
+* Fit with cov only
+
+### test_splines
+* Tesing spline interpolation

CIValidations/Scripts/CreateBinnedSplineFile.cpp

@@ -0,0 +1,79 @@
+#include "TRandom3.h"
+#include "TSpline.h"
+#include "TFile.h"
+#include "TGraph.h"
+#include "TString.h"
+#include "TH3F.h"
+
+#include <iostream>
+#include <algorithm>
+
+void CreateBinnedSplineFile(){
+
+  int nTrueEnergyBins = 5;
+  int nXBins = 10;
+  //To make splines in Etrue, x and y change this to be a larger number
+  int nYBins = 1;
+  double TrueEnergyBins[] = {0., 1.0, 2.0, 3.0, 4.0, 5.0};
+  double XBins[] = {0., 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0};
+  double YBins[] = {0.,10.0};
+
+  //Let's try out 5 systematics
+  //This gives us a total of 5*10*10*5 = 2500 splines
+  std::vector<int> SystematicKnots = {5, 7, 7, 2, 3};//, 7, 2, 7, 7, 7};
+  std::vector<std::vector<std::string>> SystematicModeNames = {{"ccqe"}, {"cc1pipm"}, {"2p2h"}, {"other"}, {"ccqe"}};
+  //Values to centre the gaussian response on
+  std::vector<double> SystematicMeanResponse = {1.0, 2.0, 1.0, 1.5, 1.0};
+  std::vector<int> SystematicNominalKnot = {3, 4, 4, 1, 2};
+  std::vector<std::string> SystematicNames = {"mysyst1", "mysyst2", "mysyst3","mysyst4","mysyst5"};
+
+  if(SystematicKnots.size() != SystematicNames.size() && SystematicNames.size() != SystematicNominalKnot.size() && SystematicKnots.size() != SystematicMeanResponse.size()) {
+    std::cerr << "Systematic vectors are not of equal length this might cause problems, please fix!" << std::endl;
+    throw;
+  }
+
+  TRandom3 Random = TRandom3(0);
+
+  auto OutputFile = std::unique_ptr<TFile>(TFile::Open("BinnedSplinesTutorialInputs.root", "RECREATE"));
+
+  //Firstly let's write the dev.tmp.0.0 file which sets the binning
+  TH3F* BinningHist = new TH3F("dev_tmp.0.0", "dev_tmp.0.0", nTrueEnergyBins, TrueEnergyBins, nXBins, XBins, nYBins, YBins);
+
+  for(auto iSyst = 0 ; iSyst < SystematicNames.size() ; ++iSyst){
+    for(auto SystematicModeName : SystematicModeNames[iSyst]){
+      for(auto SystematicKnotNumber : SystematicKnots) {
+        for(auto TrueEnergyBin_i = 0 ; TrueEnergyBin_i < nTrueEnergyBins ; ++TrueEnergyBin_i){
+          for(auto XBin_i = 0 ; XBin_i < nXBins ; ++XBin_i){
+            for(auto YBin_i = 0 ; YBin_i < nYBins ; ++YBin_i){
+              double knot_w = 1;
+              TGraph *graph = new TGraph(SystematicKnotNumber);
+
+              for(auto iKnot = 0 ; iKnot < SystematicKnotNumber ; ++iKnot) { 
+                //Check on if you are the nominal knot
+                if(iKnot != SystematicNominalKnot[iSyst]) {
+                  //Do a random throw from a gaussian 
+                  knot_w = Random.Gaus(SystematicMeanResponse[iSyst], 1.0);
+                  knot_w = std::max(0., knot_w);
+                  // std::cout << "Knot_w is " << knot_w << std::endl;
+                  //point number, x-val, y-val
+                }
+                graph->SetPoint(iKnot, iKnot, knot_w);
+              }
+              TSpline3 *Spline = new TSpline3(Form("dev.%s.%s.sp.%i.%i.%i", SystematicNames[iSyst].c_str(), SystematicModeName.c_str(), TrueEnergyBin_i, XBin_i, YBin_i),graph);
+              Spline->SetName(Form("dev.%s.%s.sp.%i.%i.%i", SystematicNames[iSyst].c_str(), SystematicModeName.c_str(), TrueEnergyBin_i, XBin_i, YBin_i));
+              //This makes things slow but removes many backup-cycles being saved to the file
+              Spline->Write(Spline->GetName(), TDirectoryFile::kOverwrite);
+              delete graph;
+              delete Spline;
+            }
+          }
+        }
+      }
+    }
+  }
+
+  OutputFile->Write();
+  OutputFile->Close();
+
+  return;
+}