diff --git a/.github/pull_request_template.md b/.github/pull_request_template.md new file mode 100644 index 000000000..359b0b96a --- /dev/null +++ b/.github/pull_request_template.md @@ -0,0 +1,7 @@ +# Pull request description: + + +## Changes or fixes: + + +## Examples: diff --git a/.github/workflows/CDImage.yml b/.github/workflows/CDImage.yml new file mode 100644 index 000000000..1957a2f23 --- /dev/null +++ b/.github/workflows/CDImage.yml @@ -0,0 +1,52 @@ +name: Image CD + +# The events that trigger the workflow +on: + push: + branches: [ develop ] + release: + types: [ published ] + +permissions: + contents: read + packages: write + +jobs: + build: + runs-on: ubuntu-latest + strategy: + matrix: + include: + - os: Ubuntu22.04 + file: Doc/MaCh3DockerFiles/Ubuntu22.04/Dockerfile + tag_latest: ubuntu22.04latest + tag_release: ubuntu22.04v1.1.0 + - os: Alma9 + file: Doc/MaCh3DockerFiles/Alma9/Dockerfile + tag_latest: alma9latest + tag_release: alma9v1.1.0 + + name: Image CD ${{ matrix.os }} + + steps: + - name: Checkout repository + uses: actions/checkout@v3 + + - name: Log in to GitHub Container Registry + run: echo "${{ secrets.GITHUB_TOKEN }}" | docker login ghcr.io -u ${{ github.actor }} --password-stdin + + - name: Build Docker image + run: | + if [ "${{ github.event_name }}" == 'release' ]; then + docker build . --file ${{ matrix.file }} --tag ghcr.io/${{ github.repository_owner }}/mach3:${{ matrix.tag_release }} --build-arg MACH3_VERSION=develop + else + docker build . --file ${{ matrix.file }} --tag ghcr.io/${{ github.repository_owner }}/mach3:${{ matrix.tag_latest }} --build-arg MACH3_VERSION=develop + fi + + - name: Push Docker image + run: | + if [ "${{ github.event_name }}" == 'release' ]; then + docker push ghcr.io/${{ github.repository_owner }}/mach3:${{ matrix.tag_release }} + else + docker push ghcr.io/${{ github.repository_owner }}/mach3:${{ matrix.tag_latest }} + fi diff --git a/.github/workflows/CIBuild.yml b/.github/workflows/CIBuild.yml new file mode 100644 index 000000000..beb19f526 --- /dev/null +++ b/.github/workflows/CIBuild.yml @@ -0,0 +1,34 @@ +name: Build CI + +# The events that trigger the workflow +on: + pull_request: + branches: [ develop ] + +permissions: + contents: read + packages: write + +jobs: + build: + runs-on: ubuntu-latest + strategy: + matrix: + include: + - os: Ubuntu22.04 + file: Doc/MaCh3DockerFiles/Ubuntu22.04/Dockerfile + tag: ubuntu22.04latest + - os: Alma9 + file: Doc/MaCh3DockerFiles/Alma9/Dockerfile + tag: alma9latest + + name: Build CI ${{ matrix.os }} + + steps: + - uses: actions/checkout@v3 + + - name: Log in to GitHub Container Registry + run: echo "${{ secrets.GITHUB_TOKEN }}" | docker login ghcr.io -u ${{ github.actor }} --password-stdin + + - name: Build the Docker image + run: docker build . --file ${{ matrix.file }} --tag ghcr.io/${{ github.repository_owner }}/mach3:${{ matrix.tag }} --build-arg MACH3_VERSION=${{ github.head_ref }} diff --git a/.github/workflows/test-result-comment.yml b/.github/workflows/test-result-comment.yml new file mode 100644 index 000000000..dbb61f326 --- /dev/null +++ b/.github/workflows/test-result-comment.yml @@ -0,0 +1,48 @@ +#Bases on https://github.com/root-project/root/blob/master/.github/workflows/test-result-comment.yml + +name: Test Summary PR comment + +on: + workflow_run: + # do NOT use quotes: https://stackoverflow.com/a/72551795/17876693 + workflows: [MaCh3 CI] + types: + - completed +permissions: {} + +jobs: + comment-test-results: + name: Publish Test Results + + if: (github.event.workflow_run.event == 'pull_request' || github.event.workflow_run.event == 'schedule') && github.event.workflow_run.conclusion != 'skipped' + + runs-on: ubuntu-latest + + permissions: + checks: write + pull-requests: write + actions: read + + steps: + - name: Download and Extract Artifacts + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + run: | + mkdir -p artifacts && cd artifacts + + artifacts_url=${{ github.event.workflow_run.artifacts_url }} + + gh api --paginate "$artifacts_url" -q '.artifacts[] | [.name, .archive_download_url] | @tsv' | while read artifact + do + IFS=$'\t' read name url <<< "$artifact" + gh api $url > "$name.zip" + unzip -d "$name" "$name.zip" + done + + - name: Publish Test Results + uses: EnricoMi/publish-unit-test-result-action@v2 + with: + commit: ${{ github.event.workflow_run.head_sha }} + event_file: artifacts/Event File/event.json + event_name: ${{ github.event.workflow_run.event }} + files: "artifacts/**/*.xml" diff --git a/.gitignore b/.gitignore index 736556653..af795af04 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,6 @@ #ignore any object files *.o +*.sif #ignore some stuff in libconfig *.lo diff --git a/CMakeLists.txt b/CMakeLists.txt index 45f6730d8..b373bfa76 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -67,7 +67,7 @@ if(NOT TARGET ROOT::ROOT) endif() if(ROOT_VERSION VERSION_LESS 6.18.0) - message(FATAL_ERROR "Using ROOT version smaller than 6.18.0, this may lead to unexpected results") + cmessage(FATAL_ERROR "Using ROOT version smaller than 6.18.0, this may lead to unexpected results") endif() #YAML for reading in config files @@ -108,7 +108,7 @@ if(NOT DEFINED USE_PROB3) SET(USE_PROB3 FALSE) endif() -# Oscillation calcualtion +# Oscillation calculation # In the future which osc calc we use might be set with a flag SET(MaCh3_Oscillator_ENABLED "") if (USE_PROB3) @@ -244,6 +244,12 @@ endif() #Include logger include(${CMAKE_CURRENT_LIST_DIR}/cmake/Modules/Logger.cmake) ################################# Features ################################## +#KS: Retrieve and print the compile options and link libraries for maximal verbose +get_target_property(compile_options MaCh3CompilerOptions INTERFACE_COMPILE_OPTIONS) +get_target_property(link_libraries MaCh3CompilerOptions INTERFACE_LINK_LIBRARIES) + +cmessage(STATUS "Compile options for MaCh3CompilerOptions: ${compile_options}") +cmessage(STATUS "Link libraries for MaCh3CompilerOptions: ${link_libraries}") LIST(APPEND ALL_FEATURES DEBUG @@ -286,21 +292,12 @@ add_library(MaCh3::All ALIAS MaCh3) configure_file(cmake/Templates/setup.MaCh3.sh.in "${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/setup.MaCh3.sh" @ONLY) install(FILES "${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/setup.MaCh3.sh" DESTINATION ${CMAKE_INSTALL_PREFIX}/bin) - -set(export_destinations - ${CMAKE_INSTALL_PREFIX}/lib/cmake/ - ${CMAKE_INSTALL_PREFIX}/ +install(EXPORT MaCh3-targets + FILE MaCh3Targets.cmake + NAMESPACE MaCh3:: + DESTINATION ${CMAKE_INSTALL_PREFIX}/ ) -foreach(dest ${export_destinations}) - install(EXPORT MaCh3-targets - FILE MaCh3Targets.cmake - NAMESPACE MaCh3:: - DESTINATION ${dest} - ) -endforeach() - - #KS: Options to print dependency graph if(NOT DEFINED MaCh3_DependancyGraph) set(MaCh3_DependancyGraph FALSE) diff --git a/Diagnostics/CombineMaCh3Chains.cpp b/Diagnostics/CombineMaCh3Chains.cpp index aa734782c..120282e98 100644 --- a/Diagnostics/CombineMaCh3Chains.cpp +++ b/Diagnostics/CombineMaCh3Chains.cpp @@ -106,7 +106,7 @@ void CombineChain() bool openedFile = fileMerger->OutputFile(OutFileName.c_str(), outFileOption.c_str(), targetCompression); if (!openedFile){ MACH3LOG_ERROR("Failed to create output file."); - throw; + throw MaCh3Exception(__FILE__ , __LINE__ ); } TFile *prevFile = nullptr; @@ -120,10 +120,10 @@ void CombineChain() TFile *file = new TFile(fileName.c_str()); if(file->Get("posteriors")->GetEntries() == 0){ - MACH3LOG_WARN("Hmmm, file {} Doesn't seem to have any entries", fileName.c_str()); - MACH3LOG_WARN("That's weird but I guess there's no rule that says a file can't be empty"); - MACH3LOG_WARN("I'll skip it but maybe double check that this doesn't indicate some deeper problem"); - continue; + MACH3LOG_WARN("Hmmm, file {} Doesn't seem to have any entries", fileName.c_str()); + MACH3LOG_WARN("That's weird but I guess there's no rule that says a file can't be empty"); + MACH3LOG_WARN("I'll skip it but maybe double check that this doesn't indicate some deeper problem"); + continue; } // EM: need to set this in the initial case @@ -143,14 +143,14 @@ void CombineChain() weirdFile = true; if(weirdFile){ - MACH3LOG_ERROR(""); - MACH3LOG_ERROR("====================================================================================="); - MACH3LOG_ERROR("This is not a great idea and could lead to weird outputs and cause some big headaches"); - MACH3LOG_ERROR("further down the road. But if you reeeeally wanna do it and you know what you're"); - MACH3LOG_ERROR("doing you can come here and remove the 'throw'"); - MACH3LOG_ERROR("{}:{}", __FILE__, __LINE__ + 2); - MACH3LOG_ERROR("====================================================================================="); - throw; + MACH3LOG_ERROR(""); + MACH3LOG_ERROR("====================================================================================="); + MACH3LOG_ERROR("This is not a great idea and could lead to weird outputs and cause some big headaches"); + MACH3LOG_ERROR("further down the road. But if you reeeeally wanna do it and you know what you're"); + MACH3LOG_ERROR("doing you can come here and remove the 'throw'"); + MACH3LOG_ERROR("{}:{}", __FILE__, __LINE__ + 2); + MACH3LOG_ERROR("====================================================================================="); + throw MaCh3Exception(__FILE__ , __LINE__ ); } // EM: file seems good, we'll add the trees to the lists fileMerger->AddFile(file); @@ -195,75 +195,75 @@ void CombineChain() } void usage(){ - std::cout << "Combine MaCh3 Chains files, very similar to hadd, but will compare embedded version info in the files to avoid accidentally combining files made with different software versions. Also avoids having a hige dump of separate version files in the output that happens with hadd." << std::endl; - std::cout << "Cmd line syntax should be:" << std::endl; - std::cout << "combineND280Splines [-h] [-c [0-9]] [-f] [-o ] input1.root [input2.root, input3.root ...] " << std::endl << std::endl; - std::cout << "inputX.root : names of individual spline files to combine, can specify any number, need at least one" << std::endl; - std::cout << "output file : name of combined spline file. optional: if not specified, the app will just use the first input file as the output, the same as hadd'" << std::endl; - std::cout << "-c : target compression level for the combined file, default is 1, in line with hadd" << std::endl; - std::cout << "-f : force overwtite the output file if it exists already" << std::endl; - std::cout << "-h : print this message and exit" << std::endl; + MACH3LOG_INFO("Combine MaCh3 Chains files, very similar to hadd, but will compare embedded version info in the files to avoid accidentally combining files made with different software versions. Also avoids having a hige dump of separate version files in the output that happens with hadd."); + MACH3LOG_INFO("Cmd line syntax should be:"); + MACH3LOG_INFO("combineND280Splines [-h] [-c [0-9]] [-f] [-o ] input1.root [input2.root, input3.root ...]"); + MACH3LOG_INFO("inputX.root : names of individual spline files to combine, can specify any number, need at least one"); + MACH3LOG_INFO("output file : name of combined spline file. optional: if not specified, the app will just use the first input file as the output, the same as hadd'"); + MACH3LOG_INFO("-c : target compression level for the combined file, default is 1, in line with hadd"); + MACH3LOG_INFO("-f : force overwrite the output file if it exists already"); + MACH3LOG_INFO("-h : print this message and exit"); } void ParseArg(int argc, char *argv[]){ - if(argc < 2){ - MACH3LOG_ERROR("Too few arguments!!"); - MACH3LOG_ERROR("USAGE:"); - usage(); - throw; + if(argc < 2){ + MACH3LOG_ERROR("Too few arguments!!"); + MACH3LOG_ERROR("USAGE:"); + usage(); + throw MaCh3Exception(__FILE__ , __LINE__ ); + } + + int c; + for(;;) { + c = getopt(argc, argv, "o:c:hf"); + if (c == -1){ // loop over the remaining arguments + while (optind < argc){ + // any non option input is assumed to be a root file + std::string fName = std::string(argv[optind]); + MACH3LOG_DEBUG("adding {} to file list", fName.c_str()); + inpFileList.push_back(fName); + optind ++; + } + break; } - - int c; - for(;;) { - c = getopt(argc, argv, "o:c:hf"); - if (c == -1){ // loop over the remaining arguments - while (optind < argc){ - // any non option input is assumed to be a root file - std::string fName = std::string(argv[optind]); - MACH3LOG_DEBUG("adding {} to file list", fName.c_str()); - inpFileList.push_back(fName); - optind ++; - } - break; - } - else{ - switch (c) { - case 'o': { - OutFileName = optarg; - break; - } - case 'f': { - forceOverwrite = true; - break; - } - case 'c': { - targetCompression = atoi(optarg); - break; - } - case 'h': { - usage(); - exit(0); - } - default: { - MACH3LOG_ERROR("Un recognised option"); - usage(); - exit(1); - } - } + else{ + switch (c) { + case 'o': { + OutFileName = optarg; + break; + } + case 'f': { + forceOverwrite = true; + break; + } + case 'c': { + targetCompression = atoi(optarg); + break; + } + case 'h': { + usage(); + exit(0); } + default: { + MACH3LOG_ERROR("Un recognised option"); + usage(); + exit(1); + } + } } + } - if(OutFileName == ""){ - MACH3LOG_INFO("Using first file in list as output: ", inpFileList[0].c_str()); - OutFileName = inpFileList[0]; - inpFileList.erase(inpFileList.begin()); - } + if(OutFileName == ""){ + MACH3LOG_INFO("Using first file in list as output: ", inpFileList[0].c_str()); + OutFileName = inpFileList[0]; + inpFileList.erase(inpFileList.begin()); + } - if(forceOverwrite){ - MACH3LOG_INFO("Will overwrite {} if it exists already", OutFileName.c_str()); - } + if(forceOverwrite){ + MACH3LOG_INFO("Will overwrite {} if it exists already", OutFileName.c_str()); + } - MACH3LOG_INFO("Combining a total of {} files into {}", inpFileList.size(), OutFileName.c_str()); + MACH3LOG_INFO("Combining a total of {} files into {}", inpFileList.size(), OutFileName.c_str()); } int main(int argc, char *argv[]) diff --git a/Diagnostics/DiagMCMC.cpp b/Diagnostics/DiagMCMC.cpp index 2b874df44..382cbe2c3 100644 --- a/Diagnostics/DiagMCMC.cpp +++ b/Diagnostics/DiagMCMC.cpp @@ -6,14 +6,13 @@ void DiagMCMC(std::string inputFile, std::string config); int main(int argc, char *argv[]) { - + SetMaCh3LoggerFormat(); if (argc != 3) { - std::cerr << "How to use: DiagMCMC MCMC_Output.root config" << std::endl; - exit(-1); + MACH3LOG_ERROR("How to use: DiagMCMC MCMC_Output.root config"); + throw MaCh3Exception(__FILE__ , __LINE__ ); } - - std::cout << "Producing single fit output" << std::endl; + MACH3LOG_INFO("Producing single fit output"); std::string filename = argv[1]; std::string config = argv[2]; DiagMCMC(filename, config); @@ -24,28 +23,28 @@ int main(int argc, char *argv[]) { void DiagMCMC(std::string inputFile, std::string config) { - std::cout << "File for study: " << inputFile << std::endl; - - YAML::Node Settings = YAML::LoadFile(config); + MACH3LOG_INFO("File for study: {}", inputFile); - // Make the processor - MCMCProcessor* Processor = new MCMCProcessor(inputFile, false); - - Processor->SetOutputSuffix("_MCMC_Diag"); - //KS:Turn off plotting detector and some other setting - Processor->SetExcludedTypes(GetFromManager>(Settings["DiagMCMC"]["ExcludedTypes"], {""})); - Processor->SetExcludedNames(GetFromManager>(Settings["DiagMCMC"]["ExcludedNames"], {""})); - Processor->SetPlotRelativeToPrior(GetFromManager(Settings["DiagMCMC"]["PlotRelativeToPrior"], false)); - //KS: Use 20 batches for batched means - Processor->SetnBatches(GetFromManager(Settings["DiagMCMC"]["nBatches"], 20)); - Processor->SetnLags(GetFromManager(Settings["DiagMCMC"]["nLags"], 25000)); - Processor->SetPrintToPDF(GetFromManager(Settings["PrintToPDF"], true)); - Processor->Initialise(); - - //KS: finally call main method - Processor->DiagMCMC(); - - delete Processor; + YAML::Node Settings = YAML::LoadFile(config); + + // Make the processor + MCMCProcessor* Processor = new MCMCProcessor(inputFile, false); + + Processor->SetOutputSuffix("_MCMC_Diag"); + //KS:Turn off plotting detector and some other setting + Processor->SetExcludedTypes(GetFromManager>(Settings["DiagMCMC"]["ExcludedTypes"], {""})); + Processor->SetExcludedNames(GetFromManager>(Settings["DiagMCMC"]["ExcludedNames"], {""})); + Processor->SetPlotRelativeToPrior(GetFromManager(Settings["DiagMCMC"]["PlotRelativeToPrior"], false)); + //KS: Use 20 batches for batched means + Processor->SetnBatches(GetFromManager(Settings["DiagMCMC"]["nBatches"], 20)); + Processor->SetnLags(GetFromManager(Settings["DiagMCMC"]["nLags"], 25000)); + Processor->SetPrintToPDF(GetFromManager(Settings["PrintToPDF"], true)); + Processor->Initialise(); + + //KS: finally call main method + Processor->DiagMCMC(); + + delete Processor; } diff --git a/Diagnostics/GetPenaltyTerm.cpp b/Diagnostics/GetPenaltyTerm.cpp index bc0ec1890..61efda29e 100644 --- a/Diagnostics/GetPenaltyTerm.cpp +++ b/Diagnostics/GetPenaltyTerm.cpp @@ -47,9 +47,9 @@ int main(int argc, char *argv[]) SetMaCh3LoggerFormat(); if (argc != 3 ) { - std::cerr<< " Something went wrong " << std::endl; - std::cerr << "./GetPenaltyTerm root_file_to_analyse.root " << std::endl; - exit(-1); + MACH3LOG_WARN("Something went wrong "); + MACH3LOG_WARN("./GetPenaltyTerm root_file_to_analyse.root "); + throw MaCh3Exception(__FILE__ , __LINE__ ); } std::string filename = argv[1]; std::string config = argv[2]; @@ -140,7 +140,7 @@ void GetPenaltyTerm(std::string inputFile, std::string configFile) { isRelevantParam[i].resize(size); int counter = 0; - //Loop over paramters in the Covariance object + //Loop over parameters in the Covariance object for (int j = 0; j < size; j++) { isRelevantParam[i][j] = false; @@ -176,7 +176,7 @@ void GetPenaltyTerm(std::string inputFile, std::string configFile) } } } - std::cout<<" Found "<GetEntries(); @@ -188,9 +188,9 @@ void GetPenaltyTerm(std::string inputFile, std::string configFile) hLogL[i]->SetLineColor(kBlue); } double* logL = new double[NSets](); - for(int n = 0; n < AllEvents; n++) + for(int n = 0; n < AllEvents; n++) { - if(n%10000 == 0) std::cout<GetEntry(n); @@ -317,7 +317,7 @@ void ReadXSecFile(std::string inputFile) if (Config == nullptr) { MACH3LOG_ERROR("Didn't find MaCh3_Config tree in MCMC file! {}", inputFile); TempFile->ls(); - throw; + throw MaCh3Exception(__FILE__ , __LINE__ ); } YAML::Node Settings = TMacroToYAML(*Config); @@ -327,8 +327,8 @@ void ReadXSecFile(std::string inputFile) if(XsecCovPos.back() == "none") { MACH3LOG_WARN("Couldn't find XsecCov branch in output"); - std::cout< ", argv[0]); + throw MaCh3Exception(__FILE__ , __LINE__ ); + } + + if (argc == 3) + { + MACH3LOG_INFO("Producing single fit output"); + config = argv[1]; + std::string filename = argv[2]; + ProcessMCMC(filename); + } + // If we want to compare two or more fits (e.g. binning changes or introducing new params/priors) + else if (argc == 6 || argc == 8) + { + MACH3LOG_INFO("Producing two fit comparison"); + config = argv[1]; + + FileNames.push_back(argv[2]); + TitleNames.push_back(argv[3]); + + FileNames.push_back(argv[4]); + TitleNames.push_back(argv[5]); + //KS: If there is third file add it + if(argc == 8) { - std::cerr << "How to use: "<< argv[0] <<" " << std::endl; - exit(-1); + FileNames.push_back(argv[6]); + TitleNames.push_back(argv[7]); } - - if (argc == 3) - { - std::cout << "Producing single fit output" << std::endl; - config = argv[1]; - std::string filename = argv[2]; - ProcessMCMC(filename); - } - // If we want to compare two or more fits (e.g. binning changes or introducing new params/priors) - else if (argc == 6 || argc == 8) - { - std::cout << "Producing two fit comparison" << std::endl; - config = argv[1]; - - FileNames.push_back(argv[2]); - TitleNames.push_back(argv[3]); - - FileNames.push_back(argv[4]); - TitleNames.push_back(argv[5]); - //KS: If there is third file add it - if(argc == 8) - { - FileNames.push_back(argv[6]); - TitleNames.push_back(argv[7]); - } - MultipleProcessMCMC(); - } - + MultipleProcessMCMC(); + } + return 0; } @@ -149,7 +150,7 @@ void MultipleProcessMCMC() Processor = new MCMCProcessor*[nFiles]; for (int ik = 0; ik < nFiles; ik++) { - std::cout << "File for study: " << FileNames[ik] << std::endl; + MACH3LOG_INFO("File for study: {}", FileNames[ik]); // Make the processor Processor[ik] = new MCMCProcessor(FileNames[ik], false); Processor[ik]->SetOutputSuffix(("_" + std::to_string(ik)).c_str()); @@ -274,7 +275,7 @@ void MultipleProcessMCMC() hpd[ik]->SetLineStyle(kSolid); } - // Find the maximum value to nicley resize hist + // Find the maximum value to nicely resize hist double maximum = 0; for (int ik = 0; ik < nFiles; ik++) maximum = std::max(maximum, hpost[ik]->GetMaximum()); for (int ik = 0; ik < nFiles; ik++) hpost[ik]->SetMaximum(1.3*maximum); @@ -296,7 +297,7 @@ void MultipleProcessMCMC() } delete[] hpost; delete[] hpd; - }//End loop over paramters + }//End loop over parameters // Finally draw the parameter plot onto the PDF // Close the .pdf file with all the posteriors @@ -315,7 +316,7 @@ void MultipleProcessMCMC() delete[] Processor; } -// KS: Calculate Bayes factor for a given hiphothesis, most onformative are those related to osc params. However, it make realtive easy interpreation for switch dials +// KS: Calculate Bayes factor for a given hypothesis, most informative are those related to osc params. However, it make relative easy interpretation for switch dials void CalcBayesFactor(MCMCProcessor* Processor) { YAML::Node card_yaml = YAML::LoadFile(config.c_str()); @@ -472,7 +473,7 @@ void DiagnoseCovarianceMatrix(MCMCProcessor* Processor, std::string inputFile) TH2D *CovarianceDiff = (TH2D*)CovarianceHist->Clone("Covariance_Ratio"); TH2D *CorrelationDiff = (TH2D*)CorrelationHist->Clone("Correlation_Ratio"); - //KS: Bit messy but quite often covariance is 0 is divided by 0 is problemiatic so + //KS: Bit messy but quite often covariance is 0 is divided by 0 is problematic so #ifdef MULTITHREAD #pragma omp parallel for #endif @@ -674,7 +675,6 @@ void KolmogorovSmirnovTest(MCMCProcessor** Processor, TCanvas* Posterior, TStrin for (int j = 1; j < NumberOfBins+1; ++j) { Cumulative += hpost[ik]->GetBinContent(j)/Integral; - CumulativeDistribution[ik]->SetBinContent(j, Cumulative); } //KS: Set overflow to 1 just in case @@ -695,7 +695,7 @@ void KolmogorovSmirnovTest(MCMCProcessor** Processor, TCanvas* Posterior, TStrin { double BinValue = CumulativeDistribution[0]->GetBinCenter(j); int BinNumber = CumulativeDistribution[ik]->FindBin(BinValue); - //KS: Calculate D statistic for this bin, only save it if it's bigger than previosly found value + //KS: Calculate D statistic for this bin, only save it if it's bigger than previously found value double TempDstat = std::fabs(CumulativeDistribution[0]->GetBinContent(j) - CumulativeDistribution[ik]->GetBinContent(BinNumber)); if(TempDstat > TestStatD[ik]) { @@ -737,9 +737,9 @@ void KolmogorovSmirnovTest(MCMCProcessor** Processor, TCanvas* Posterior, TStrin delete leg; for (int ik = 0; ik < nFiles; ik++) { - delete hpost[ik]; - delete CumulativeDistribution[ik]; - delete LineD[ik]; + delete hpost[ik]; + delete CumulativeDistribution[ik]; + delete LineD[ik]; } delete[] hpost; delete[] CumulativeDistribution; diff --git a/Diagnostics/README.md b/Diagnostics/README.md new file mode 100644 index 000000000..3f933a3e6 --- /dev/null +++ b/Diagnostics/README.md @@ -0,0 +1,51 @@ +### Plotting and Diagnostic +**ProcessMCMC** - The main app you want to use for analysing the ND280 chain. It prints posterior distribution after burn-in the cut. Moreover, you can compare two/three different chains. There are a few options you can modify easily inside the app like selection, burn-in cut, and whether to plot xse+flux or only flux. Other functionality +
    +
  1. Produce a covariance matrix with multithreading (which requires lots of RAM due to caching)
    2. +
  2. Violin plots
    3. +
  3. Credible intervals and regions
    4. +
  4. Calculate Bayes factor and give significance based on Jeffreys scale
    5. +
  5. Produce triangle plots
    6. +
  6. Study covariance matrix stability
    7. +
+ +**GetPostfitParamPlots** - This will plot output from ProcessMCMC for nice plots which can go to TN. Bits are hardcoded to make plots nicer users should be careful when using the non-conventional xsec model. If you used `ProcessMCMC` with `PlotDet` you will also get an overlay of detector parameters (ND or ND+FD depending on chain type). If Violin plot was produced in `ProcessMCMC` you will get fancy plots here as well. + +**GetPenaltyTerm** - Since xsec and flux and ND spline systematic are treated as the same systematic object we cannot just take log_xsec, hence we need this script, use `GetFluxPenaltyTerm MCMCChain.root config`. Parameters of relevance are loaded via config, thus you can study any combination you want. Time needed increases with number of sets :( + +**DiagMCMC** - Perform MCMC diagnostic like autocorrelation or trace plots. + +**RHat** - Performs RHat diagnostic to study if all used chains converged to the same stationary distribution. +``` +./RHat Ntoys MCMCchain_1.root MCMCchain_2.root MCMCchain_3.root ... [how many you like] +``` + +**PlotLLH** - Plot LLH scans, flexible and configurable in command line. can take any number of LLH scans as input, will use the first one as a baseline when making e.g. ratio plots. The first file must be a MaCh3 scan. +options: + + -r overlay ratio plots + + -s also make plots split by sample contribution, to use this option, the LLH scan must have been run with the option `LLH_SCAN_BY_SAMPLE = true` in the config file + + -g draw a grid on the plots + + -l a string specifying the labels for each scan to plot in the legent. this should be a string of labels separated by semi colons, e.g.: -`l "label1;label2;label3"` + + -o the name of the output pdf + + -d a string specifying additional drawing options to pass to the histogram draw calls, e.g. `-d "C"` will plot smooth curves through the histogram bins. See https://root.cern/doc/master/classTHistPainter.html#HP01a for possible options. + + +**CombineMaCh3Chains** - will combine chains files produced by **MCMC**, enforcing the condition that all the files to combine were made using the exact same software versions and config files +``` +CombineMaCh3Chains [-h] [-c [0-9]] [-f] [-o ] file1.root [file2.root, file3.root ...] +``` +*fileX.root* are the individual spline files to combine, can specify any number, need at least one + +-c target compression level for the combined file, default is 1, in line with hadd + +-f force overwrite of the combined file if it exists already + +-h print usage message and exit + +*Output file* (optional) name of the combined file. If not specified, will just use *file1.root*, the first in the list of files, same as *hadd*. diff --git a/Diagnostics/RHat.cpp b/Diagnostics/RHat.cpp index 5d4cc77ec..9c33c7229 100644 --- a/Diagnostics/RHat.cpp +++ b/Diagnostics/RHat.cpp @@ -37,7 +37,6 @@ int Nchains; int nDraw; -bool VERBOSE; std::vector BranchNames; std::vector MCMCFile; std::vector ValidPar; @@ -86,40 +85,40 @@ int main(int argc, char *argv[]) { // ******************* SetMaCh3LoggerFormat(); + MaCh3Utils::MaCh3Welcome(); - Draws = NULL; - Mean = NULL; - StandardDeviation = NULL; + Draws = nullptr; + Mean = nullptr; + StandardDeviation = nullptr; - MeanGlobal = NULL; - StandardDeviationGlobal = NULL; + MeanGlobal = nullptr; + StandardDeviationGlobal = nullptr; - BetweenChainVariance = NULL; - MarginalPosteriorVariance = NULL; - RHat = NULL; - EffectiveSampleSize = NULL; + BetweenChainVariance = nullptr; + MarginalPosteriorVariance = nullptr; + RHat = nullptr; + EffectiveSampleSize = nullptr; - DrawsFolded = NULL; - MedianArr = NULL; - MeanFolded = NULL; - StandardDeviationFolded = NULL; + DrawsFolded = nullptr; + MedianArr = nullptr; + MeanFolded = nullptr; + StandardDeviationFolded = nullptr; - MeanGlobalFolded = NULL; - StandardDeviationGlobalFolded = NULL; + MeanGlobalFolded = nullptr; + StandardDeviationGlobalFolded = nullptr; - BetweenChainVarianceFolded = NULL; - MarginalPosteriorVarianceFolded = NULL; - RHatFolded = NULL; - EffectiveSampleSizeFolded = NULL; + BetweenChainVarianceFolded = nullptr; + MarginalPosteriorVarianceFolded = nullptr; + RHatFolded = nullptr; + EffectiveSampleSizeFolded = nullptr; - VERBOSE = false; Nchains = 0; - if (argc == 1 && argc == 2) + if (argc == 1 || argc == 2) { - std::cerr << "Wrong arguments" << std::endl; - std::cerr << "./RHat Ntoys MCMCchain_1.root MCMCchain_2.root MCMCchain_3.root ... [how many you like]" << std::endl; - exit(-1); + MACH3LOG_ERROR("Wrong arguments"); + MACH3LOG_ERROR("./RHat Ntoys MCMCchain_1.root MCMCchain_2.root MCMCchain_3.root ... [how many you like]"); + throw MaCh3Exception(__FILE__ , __LINE__ ); } Ntoys = atoi(argv[1]); @@ -128,22 +127,22 @@ int main(int argc, char *argv[]) { for (int i = 2; i < argc; i++) { MCMCFile.push_back(std::string(argv[i])); - std::cout<<"Adding file: "<Add(MCMCFile[m].c_str()); - + MACH3LOG_INFO("On file: {}", MCMCFile[m].c_str()); nEntries[m] = Chain->GetEntries(); // Set the step cut to be 20% @@ -237,7 +237,7 @@ void PrepareChains() { } } Chain->SetBranchStatus(bname, true); - if (VERBOSE) std::cout<= nEntries[m]) { - std::cerr<<"You are running on a chain shorter than BurnIn cut"< 10 && i % (Ntoys/10) == 0) { - std::cout << "On toy " << i << "/" << Ntoys << " (" << int(double(i)/double(Ntoys)*100.0) << "%)" << std::endl; - std::cout << " Getting random entry " << entry << std::endl; + MaCh3Utils::PrintProgressBar(i+m*Ntoys, Ntoys*Nchains); + MACH3LOG_DEBUG("Getting random entry {}", entry); } // Set the branch addresses for params @@ -356,7 +354,7 @@ void PrepareChains() { delete[] nBranches; clock.Stop(); - std::cout << "Finsihed producing Toys, it took " << clock.RealTime() << "s to finish" << std::endl; + MACH3LOG_INFO("Finished calculating Toys, it took {:.2f}s to finish", clock.RealTime()); } // ******************* @@ -364,8 +362,7 @@ void PrepareChains() { void InitialiseArrays() { // ******************* - std::cout << "Initialising arrays " << std::endl; - + MACH3LOG_INFO("Initialising arrays"); Mean = new double*[Nchains](); StandardDeviation = new double*[Nchains](); @@ -577,7 +574,7 @@ void CalcRhat() { #endif clock.Stop(); - std::cout << "Finsihed calcaulting RHat, it took " << clock.RealTime() << "s to finish" << std::endl; + MACH3LOG_INFO("Finished calculating RHat, it took {:.2f}s to finish", clock.RealTime()); } @@ -585,7 +582,7 @@ void CalcRhat() { void SaveResults() { // ******************* std::string NameTemp = ""; - //KS: If we run over many many chains there is danger that name will be so absurdly long we run over systmat limit and job will be killed :( + //KS: If we run over many many chains there is danger that name will be so absurdly long we run over system limit and job will be killed :( if(Nchains < 5) { for (int i = 0; i < Nchains; i++) @@ -599,9 +596,8 @@ void SaveResults() { NameTemp = NameTemp + temp + "_"; } } - else - { - NameTemp = std::to_string(Nchains) + "Chains" + "_"; + else { + NameTemp = std::to_string(Nchains) + "Chains" + "_"; } NameTemp += "diag.root"; @@ -653,13 +649,12 @@ void SaveResults() { } } //KS: We set criterium of 1.1 based on Gelman et al. (2003) Bayesian Data Analysis - std::cout<<" Number of parameters which has R hat greater than 1.1 is "< 1.1 || RHatFolded[j] > 1.1) && ValidPar[j]) { - std::cout<<"Parameter "<Write(); @@ -709,7 +704,7 @@ void SaveResults() { Legend->Draw("same"); TempCanvas->Write("Rhat"); delete Legend; - Legend = NULL; + Legend = nullptr; //Now R hat for log L RhatLogPlot->GetXaxis()->SetTitle("R hat for LogL"); @@ -734,7 +729,7 @@ void SaveResults() { Legend->Draw("same"); TempCanvas->Write("RhatLog"); delete Legend; - Legend = NULL; + Legend = nullptr; //Now canvas for effective sample size EffectiveSampleSizePlot->GetXaxis()->SetTitle("S_{eff, BDA2}"); @@ -785,7 +780,7 @@ void SaveResults() { DiagFile->Close(); delete DiagFile; - std::cout<<"Finished and wrote results to "< MaCh3 @@ -110,10 +109,21 @@ Most of external libraries are being handled through CPM. The only external libr Based on several test here are recommended version: ``` - GCC: ... + GCC: >= 8.5 [lower versions may work] CMake: >= 3.14 ROOT: >= 6.18 ``` +### Supported operational systems +| Name | Status | +|-------------|--------| +| Alma9 | ✅ | +| Ubuntu22.04 | ✅ | +| CentOS7 | ❔ | +| Windows | ❌ | + +✅ - Part of CI/CD
+❔ - Not part of CI/CD but used by some users/developers so it might work
+❌ - Not supported and no plans right now
# How To Use This is an example how your executable can look like using MaCh3: @@ -143,55 +153,4 @@ This is an example how your executable can look like using MaCh3: ### 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. Currently available utils include: - -**ProcessMCMC** - The main app you want to use for analysing the ND280 chain. It prints posterior distribution after burn-in the cut. Moreover, you can compare two/three different chains. There are a few options you can modify easily inside the app like selection, burn-in cut, and whether to plot xse+flux or only flux. Other functionality -
    -
  1. Produce a covariance matrix with multithreading (which requires lots of RAM due to caching)
    2. -
  2. Violin plots
    3. -
  3. Credible intervals and regions
    4. -
  4. Calculate Bayes factor and give significance based on Jeffreys scale
    5. -
  5. Produce triangle plots
    6. -
  6. Study covariance matrix stability
    7. -
- -**GetPostfitParamPlots** - This will plot output from ProcessMCMC for nice plots which can go to TN. Bits are hardcoded to make plots nicer users should be careful when using the non-conventional xsec model. If you used `ProcessMCMC` with `PlotDet` you will also get an overlay of detector parameters (ND or ND+FD depending on chain type). If Violin plot was produced in `ProcessMCMC` you will get fancy plots here as well. - -**GetPenaltyTerm** - Since xsec and flux and ND spline systematic are treated as the same systematic object we cannot just take log_xsec, hence we need this script, use `GetFluxPenaltyTerm MCMCChain.root config`. Parameters of relevance are loaded via config, thus you can study any combination you want. Time needed increases with number of sets :( - -**DiagMCMC** - Perform MCMC diagnostic like autocorrelation or trace plots. - -**RHat** - Performs RHat diagnostic to study if all used chains converged to the same stationary distribution. -``` -./RHat Ntoys MCMCchain_1.root MCMCchain_2.root MCMCchain_3.root ... [how many you like] -``` - -**PlotLLH** - Plot LLH scans, flexible and configurable in command line. can take any number of LLH scans as input, will use the first one as a baseline when making e.g. ratio plots. The first file must be a MaCh3 scan. -options: - - -r overlay ratio plots - - -s also make plots split by sample contribution, to use this option, the LLH scan must have been run with the option `LLH_SCAN_BY_SAMPLE = true` in the config file - - -g draw a grid on the plots - - -l a string specifying the labels for each scan to plot in the legent. this should be a string of labels separated by semi colons, e.g.: -`l "label1;label2;label3"` - - -o the name of the output pdf - - -d a string specifying additional drawing options to pass to the histogram draw calls, e.g. `-d "C"` will plot smooth curves through the histogram bins. See https://root.cern/doc/master/classTHistPainter.html#HP01a for possible options. - - -**CombineMaCh3Chains** - will combine chains files produced by **MCMC**, enforcing the condition that all the files to combine were made using the exact same software versions and config files -``` -CombineMaCh3Chains [-h] [-c [0-9]] [-f] [-o ] file1.root [file2.root, file3.root ...] -``` -*fileX.root* are the individual spline files to combine, can specify any number, need at least one - --c target compression level for the combined file, default is 1, in line with hadd - --f force overwrite of the combined file if it exists already - --h print usage message and exit - -*Output file* (optional) name of the combined file. If not specified, will just use *file1.root*, the first in the list of files, same as *hadd*. +Example of chain diagnostic utils can be found [here](https://github.com/mach3-software/MaCh3/tree/develop/Diagnostics) with example of config. diff --git a/cmake/Modules/CUDASamples.cmake b/cmake/Modules/CUDASamples.cmake new file mode 100644 index 000000000..389236180 --- /dev/null +++ b/cmake/Modules/CUDASamples.cmake @@ -0,0 +1,33 @@ +# Variable to hold the paths once found +set(CMAKE_CUDA_SAMPLES_PATH "") +set(CUDASAMPLES_FOUND FALSE) + +#KS: List of directories to search for helper_functions.h and cuda_runtime.h +# Please expand if needed +set(CUDA_SAMPLE_SEARCH_PATHS + "${CUDAToolkit_TARGET_DIR}/samples/common/inc/" +) + +# Loop over the directories and append to CMAKE_CUDA_SAMPLES_PATH if headers are found +foreach(dir ${CUDA_SAMPLE_SEARCH_PATHS}) + if(EXISTS "${dir}/helper_functions.h") + list(APPEND CMAKE_CUDA_SAMPLES_PATH ${dir}) + set(CUDASAMPLES_FOUND TRUE) + endif() +endforeach() + +#HW: If we're doing GPU stuff, we need the CUDA helper module +if(NOT CUDASAMPLES_FOUND) + cmessage(WARNING "Adding Built-in CUDA Sample Library, it make take some time") + CPMAddPackage( + NAME cuda-samples + GITHUB_REPOSITORY "NVIDIA/cuda-samples" + GIT_TAG v12.3 + SOURCE_SUBDIR Common + DOWNLOAD_ONLY YES + ) + list(APPEND CMAKE_CUDA_SAMPLES_PATH ${cuda-samples_SOURCE_DIR}/Common) +endif() + +cmessage(STATUS "Using the following CUDA samples paths: ${CMAKE_CUDA_SAMPLES_PATH}") +target_include_directories(MaCh3CompilerOptions INTERFACE ${CMAKE_CUDA_SAMPLES_PATH}) diff --git a/cmake/Modules/CUDASetup.cmake b/cmake/Modules/CUDASetup.cmake index fca720298..093678767 100644 --- a/cmake/Modules/CUDASetup.cmake +++ b/cmake/Modules/CUDASetup.cmake @@ -53,7 +53,7 @@ if(NOT MaCh3_DEBUG_ENABLED) "$<$:-Xcompiler=-fpic;-Xcompiler=-O3;-Xcompiler=-Wall;-Xcompiler=-Wextra;-Xcompiler=-Werror;-Xcompiler=-Wno-error=unused-parameter>" ) else() -#CWret: -g and -G for debug flags to use cuda-gdb; slows stuff A LOT +#CW: -g and -G for debug flags to use cuda-gdb; slows stuff A LOT #-pxtas-options=-v, -maxregcount=N target_compile_options(MaCh3CompilerOptions INTERFACE "$<$:-prec-sqrt=false;-use_fast_math;-Werror;cross-execution-space-call;-w>" @@ -63,42 +63,8 @@ else() target_compile_definitions(MaCh3CompilerOptions INTERFACE "$<$:CUDA_ERROR_CHECK>") endif() - -target_link_options(MaCh3CompilerOptions INTERFACE -I$ENV{CUDAPATH}/lib64 -I$ENV{CUDAPATH}/include -I$ENV{CUDAPATH}/common/inc -I$ENV{CUDAPATH}/samples/common/inc) - -if(NOT DEFINED ENV{CUDAPATH}) - cmessage(FATAL_ERROR "CUDAPATH environment variable is not defined. Please set it to the root directory of your CUDA installation.") -endif() - -#if(NOT DEFINED CUDA_SAMPLES) -# cmessage(FATAL_ERROR "When using CUDA, CUDA_SAMPLES must be defined to point to the CUDAToolkit samples directory (should contain common/helper_functions.h).") -#endif() - -# Set the CUDA samples path -set(CMAKE_CUDA_SAMPLES_PATH - $ENV{CUDAPATH}/include - $ENV{CUDAPATH}/samples/common/inc -) - -cmessage(STATUS "Using the following CUDA samples paths: ${CMAKE_CUDA_SAMPLES_PATH}") -target_include_directories(MaCh3CompilerOptions INTERFACE ${CMAKE_CUDA_SAMPLES_PATH}) - -#HW: If we're doing GPU stuff, we need the CUDA helper module -if(BUILTIN_CUDASAMPLES_ENABLED ) - cmessage("Adding Builtin CUDA Sample Library") - CPMAddPackage( - NAME cuda-samples - GITHUB_REPOSITORY "NVIDIA/cuda-samples" - GIT_TAG v12.3 - SOURCE_SUBDIR Common - DOWNLOAD_ONLY YES - ) -# Now we add the library - if(cuda-samples_ADDED) - add_library(cuda-samples INTERFACE IMPORTED) - target_include_directories(cuda-samples INTERFACE ${cuda-samples_SOURCE_DIR}) - endif() -endif() +target_include_directories(MaCh3CompilerOptions INTERFACE ${CUDAToolkit_INCLUDE_DIRS}) +include(${CMAKE_CURRENT_LIST_DIR}/CUDASamples.cmake) #KS: Keep this for backward compatibility @@ -106,7 +72,7 @@ endif() #-Xcompiler "-fpic" -c #-prec-sqrt=false -use_fast_math -#CWRET comment: could change arch here. also, prec-div is only used in (seg+khig)/2; could replace by *0.5 +#CW comment: could change arch here. also, prec-div is only used in (seg+khig)/2; could replace by *0.5 # -prec-sqrt is not needed (no sqrts in program code!), -use_fast_math forces intrinsic # tried this on GTX 660 and no speed increase, buu and P6 6c # don't use fastmath! diff --git a/manager/MaCh3Logger.h b/manager/MaCh3Logger.h index 5d2f6bbef..6045b8dc5 100644 --- a/manager/MaCh3Logger.h +++ b/manager/MaCh3Logger.h @@ -1,6 +1,10 @@ #pragma once #include "spdlog/spdlog.h" +#include +#include +#include +#include //KS: Based on this https://github.com/gabime/spdlog/blob/a2b4262090fd3f005c2315dcb5be2f0f1774a005/include/spdlog/spdlog.h#L284 @@ -23,3 +27,48 @@ inline void SetMaCh3LoggerFormat() spdlog::set_pattern("[%H:%M:%S][%s][%^%l%$] %v"); #endif } + +/// @brief KS: This is bit convoluted but this is to allow redirecting cout and errors from external library into MaCh3 logger format +/// @tparam Func The type of the function to be called. +/// @tparam Args The types of the arguments to be passed to the function. +/// @param LibName The name of the library to be included in the log output. +/// @param func The function to be called whose output needs to be captured. +/// @param args The arguments to be passed to the function. +/// @code +/// void ExampleFunc(int x, const std::string& str) { +/// std::cout << "Output from exampleFunc: " << x << ", " << str << "\n"; +/// std::cerr << "Error from exampleFunc: " << x << ", " << str << "\n"; +/// } +/// +/// int main() { +/// LoggerPrint("ExampleLib", static_cast(blarb), 666, "Number of the BEAST"); +/// return 0; +/// } +/// @endcode +template +void LoggerPrint(const std::string& LibName, Func&& func, Args&&... args) +{ + // Create a stringstream to capture the output + std::stringstream sss; + + // Save the original stream buffers + std::streambuf* coutBuf = std::cout.rdbuf(); + std::streambuf* cerrBuf = std::cerr.rdbuf(); + + // Redirect std::cout and std::cerr to the stringstream buffer + std::cout.rdbuf(sss.rdbuf()); + std::cerr.rdbuf(sss.rdbuf()); + + // Call the provided function + func(std::forward(args)...); + + // Restore the original stream buffers + std::cout.rdbuf(coutBuf); + std::cerr.rdbuf(cerrBuf); + + std::string line; + while (std::getline(sss, line)) + { + MACH3LOG_INFO("[{}] {}", LibName, line); + } +} diff --git a/mcmc/CMakeLists.txt b/mcmc/CMakeLists.txt index a2d2bb890..56e3cfb84 100644 --- a/mcmc/CMakeLists.txt +++ b/mcmc/CMakeLists.txt @@ -2,7 +2,7 @@ set(HEADERS FitterBase.h mcmc.h LikelihoodFit.h - $<$>:MinuitFit.h> + MinuitFit.h PSO.h MCMCProcessor.h SampleSummary.h @@ -13,7 +13,7 @@ add_library(MCMC SHARED FitterBase.cpp mcmc.cpp LikelihoodFit.cpp - $<$>:MinuitFit.cpp> + $<$:MinuitFit.cpp> PSO.cpp MCMCProcessor.cpp SampleSummary.cpp