Merge branch 'develop' into feature_CUDAbot

mach3-software · Nov 7, 2024 · 5b2a4e1 · 5b2a4e1
2 parents 9b97005 + 66a293b
commit 5b2a4e1
Show file tree

Hide file tree

Showing 16 changed files with 385 additions and 323 deletions.
diff --git a/.github/workflows/Newsletter.yml b/.github/workflows/Newsletter.yml
@@ -27,10 +27,17 @@ jobs:
           latest_tag=$(git describe --tags --abbrev=0)
           echo "latest_tag=${latest_tag}" >> "$GITHUB_ENV"
 
-      - name: Get the current date
+      - name: Get the most active author in the last week for MaCh3
         run: |
-          current_date=$(date +"%d/%m/%Y")
-          echo "current_date=${current_date}" >> "$GITHUB_ENV"
+          most_active_author=$(git log --since='1 week ago' --pretty=format:"%an" | sort | uniq -c | sort -nr | head -n 1)
+          echo "most_active_author=${most_active_author}" >> "$GITHUB_ENV"
+
+      - name: Get the author with most commits (count and name)
+        run: |
+          COUNT=$(echo "$most_active_author" | awk '{print $1}')
+          NAME=$(echo "$most_active_author" | cut -d ' ' -f2-)
+          echo "maCh3_most_active_author_count=${COUNT}" >> "$GITHUB_ENV"
+          echo "maCh3_most_active_author_name=${NAME}" >> "$GITHUB_ENV"
 
       - name: Checkout the MaCh3Tutorial repository
         uses: actions/checkout@v4
@@ -52,6 +59,18 @@ jobs:
           tutorial_latest_tag=$(git -C MaCh3Tutorial describe --tags --abbrev=0)
           echo "tutorial_latest_tag=${tutorial_latest_tag}" >> "$GITHUB_ENV"
 
+      - name: Get the most active author in the last week for MaCh3Tutorial
+        run: |
+          most_active_author=$(git -C MaCh3Tutorial log --since='1 week ago' --pretty=format:"%an" | sort | uniq -c | sort -nr | head -n 1)
+          echo "tutorial_most_active_author=${most_active_author}" >> "$GITHUB_ENV"
+
+      - name: Get the author with most commits for MaCh3Tutorial (count and name)
+        run: |
+          COUNT=$(echo "$tutorial_most_active_author" | awk '{print $1}')
+          NAME=$(echo "$tutorial_most_active_author" | cut -d ' ' -f2-)
+          echo "maCh3Tutorial_most_active_author_count=${COUNT}" >> "$GITHUB_ENV"
+          echo "maCh3Tutorial_most_active_author_name=${NAME}" >> "$GITHUB_ENV"
+
       - name: Checkout the MaCh3-PythonUtils repository
         uses: actions/checkout@v4
         with:
@@ -72,6 +91,18 @@ jobs:
           python_utils_latest_tag=$(git -C MaCh3-PythonUtils describe --tags --abbrev=0)
           echo "python_utils_latest_tag=${python_utils_latest_tag}" >> "$GITHUB_ENV"
 
+      - name: Get the most active author in the last week for MaCh3-PythonUtils
+        run: |
+          most_active_author=$(git -C MaCh3-PythonUtils log --since='1 week ago' --pretty=format:"%an" | sort | uniq -c | sort -nr | head -n 1)
+          echo "python_utils_most_active_author=${most_active_author}" >> "$GITHUB_ENV"
+
+      - name: Get the author with most commits for MaCh3-PythonUtils (count and name)
+        run: |
+          COUNT=$(echo "$python_utils_most_active_author" | awk '{print $1}')
+          NAME=$(echo "$python_utils_most_active_author" | cut -d ' ' -f2-)
+          echo "maCh3PythonUtils_most_active_author_count=${COUNT}" >> "$GITHUB_ENV"
+          echo "maCh3PythonUtils_most_active_author_name=${NAME}" >> "$GITHUB_ENV"
+
       - name: Create a new GitHub Discussion
         uses: abirismyname/[email protected]
         env:
@@ -87,16 +118,19 @@ jobs:
             - Total number of commits: ${{ env.commit_count }}
             - Total number of commits in the last week: ${{ env.commits_last_week }}
             - Most recent tag: ${{ env.latest_tag }}
+            - Most active author in the last week: ${{ env.maCh3_most_active_author_name }} (Commits: ${{ env.maCh3_most_active_author_count }})
 
             **MaCh3Tutorial Repository**
             - Total number of commits: ${{ env.tutorial_commit_count }}
             - Total number of commits in the last week: ${{ env.tutorial_commits_last_week }}
             - Most recent tag: ${{ env.tutorial_latest_tag }}
+            - Most active author in the last week: ${{ env.maCh3Tutorial_most_active_author_name }} (Commits: ${{ env.maCh3Tutorial_most_active_author_count }})
 
             **MaCh3-PythonUtils Repository**
             - Total number of commits: ${{ env.python_utils_commit_count }}
             - Total number of commits in the last week: ${{ env.python_utils_commits_last_week }}
             - Most recent tag: ${{ env.python_utils_latest_tag }}
+            - Most active author in the last week: ${{ env.maCh3PythonUtils_most_active_author_name }} (Commits: ${{ env.maCh3PythonUtils_most_active_author_count }})
 
             Cheers,
             MaCh3-bot

diff --git a/Diagnostics/ProcessMCMC.cpp b/Diagnostics/ProcessMCMC.cpp
@@ -93,7 +93,11 @@ void ProcessMCMC(const std::string& inputFile)
   Processor->SetPost2DPlotThreshold(GetFromManager<double>(Settings["Post2DPlotThreshold"], 0.2));
 
   Processor->Initialise();
-
+  if(Settings["Thinning"])
+  {
+    bool DoThin = Settings["Thinning"][0].as<bool>();
+    if(DoThin) Processor->ThinMCMC(Settings["Thinning"][1].as<int>());
+  }
   // Make the postfit
   Processor->MakePostfit();
   Processor->DrawPostfit();

diff --git a/Doc/bibliography.bib b/Doc/bibliography.bib
@@ -168,3 +168,15 @@ @article{gelman2019
     doi          = "10.48550/arXiv.1903.08008",
     url          = "https://doi.org/10.48550/arXiv.1903.08008"
 }
+
+@article{2011ThinningMCMC,
+  author       = {William A. Link and Mitchell J. Eaton},
+  title        = {On thinning of chains in MCMC},
+  journal      = {Methods in Ecology and Evolution},
+  volume       = {2},
+  number       = {3},
+  pages        = {305-310},
+  year         = {2011},
+  doi          = {10.1111/j.2041-210X.2011.00131.x},
+  url          = {https://doi.org/10.1111/j.2041-210X.2011.00131.x},
+}
diff --git a/manager/gpuUtils.cuh b/manager/gpuUtils.cuh
@@ -8,7 +8,7 @@
 // CUDA specifics
 // Because CUDA is cuda, need to make sure we don't check C-style floats...
 #pragma GCC diagnostic push
-#pragma GCC diagnostics ignored "-Wold-style-cast"
+#pragma GCC diagnostic ignored "-Wold-style-cast"
 #include <cuda_runtime.h>
 #pragma GCC diagnostic pop
 

diff --git a/mcmc/MCMCProcessor.h b/mcmc/MCMCProcessor.h
@@ -159,6 +159,10 @@ class MCMCProcessor {
     void ParameterEvolution(const std::vector<std::string>& Names,
                             const std::vector<int>& NIntervals);
 
+    /// @brief Thin MCMC Chain, to save space and maintain low autocorrelations.
+    /// @param ThinningCut every which entry you want to thin
+    inline void ThinMCMC(const int ThinningCut) {ThinningMCMC(MCMCFile+".root", ThinningCut);};
+
     /// @brief KS: Perform MCMC diagnostic including Autocorrelation, Trace etc.
     void DiagMCMC();
 

diff --git a/mcmc/PSO.cpp b/mcmc/PSO.cpp
@@ -2,8 +2,9 @@
 
 #include <cmath>
 
+// ***************
 PSO::PSO(manager *man) : LikelihoodFit(man) {
-
+// ***************
   fConstriction = fitMan->raw()["General"]["PSO"]["Constriction"].as<double>();
   fInertia = fitMan->raw()["General"]["PSO"]["Inertia"].as<double>()*fConstriction;
   fOne = fitMan->raw()["General"]["PSO"]["One"].as<double>()*fConstriction;
@@ -20,8 +21,9 @@ PSO::PSO(manager *man) : LikelihoodFit(man) {
   }
 }
 
+// ***************
 void PSO::runMCMC(){
-
+// ***************
   PrepareFit();
 
   if(fTestLikelihood){
@@ -41,17 +43,16 @@ void PSO::runMCMC(){
   return;
 }
 
-
+// *************************
 void PSO::init(){
-
+// *************************
   fBestValue = 1234567890.0;
 
   //KS: For none PCA this will be eqaul to normal parameters
   //const int NparsPSOFull = NPars;
   //const int NparsPSO = NParsPCA;
 
-  std::cout << "Preparing PSO" << std::endl;
-
+  MACH3LOG_INFO("Preparing PSO");
   // Initialise bounds on parameters
   if(fTestLikelihood){
     for (int i = 0; i < fDim; i++){
@@ -105,9 +106,9 @@ void PSO::init(){
     }
   }
 
-  std::cout << "Printing Minimums and Maximums of Variables to be minimised" << std::endl;
-  for (int i =0; i<fDim; i++){
-    std::cout << "Variable "<< i<<" :" << ranges_min[i] << ", "<< ranges_max[i] << std::endl;
+  MACH3LOG_INFO("Printing Minimums and Maximums of Variables to be minimized");
+  for (int i = 0; i < fDim; i++){
+    MACH3LOG_INFO("Variable {} : {:.2f}, {:.2f}", i, ranges_min[i], ranges_max[i]);
   }
 
   // Initialise particle positions
@@ -143,10 +144,12 @@ void PSO::init(){
   }
 }
 
+// *************************
 std::vector<std::vector<double> > PSO::bisection(std::vector<double>position,double minimum, double range, double precision){
+// *************************
   std::vector<std::vector<double>> uncertainties_list;
-  for (unsigned int i = 0; i< position.size(); ++i){
-    std::cout << i << std::endl;
+  for (unsigned int i = 0; i< position.size(); ++i) {
+    MACH3LOG_INFO("{}", i);
     //std::vector<double> uncertainties;
     std::vector<double> new_position = position; new_position[i] = position[i]-range;
     double val_1 = CalcChi(new_position)-minimum-1.0;
@@ -201,7 +204,8 @@ std::vector<std::vector<double> > PSO::bisection(std::vector<double>position,dou
         position_list_p[1] = new_bisect_position_p;
         res = std::abs(position[2]-position[0]);
         res_p = std::abs(position_list_p[1][i]-position_list_p[0][i]);
-        //std::cout << "Pos midpoint is " << position_list_p[1][i] << std::endl;
+        MACH3LOG_TRACE("Pos midpoint is {:.2f}", position_list_p[1][i]);
+
       }
       else{
         std::vector<double> new_bisect_position_p = position_list_p[1];new_bisect_position_p[i] += (position_list_p[2][i]-position_list_p[1][i])/2;
@@ -211,17 +215,20 @@ std::vector<std::vector<double> > PSO::bisection(std::vector<double>position,dou
         value_list_p[1] = new_val_p;
         position_list_p[1] = new_bisect_position_p;
         res_p = std::abs(position_list_p[2][i]-position_list_p[1][i]);
-        //std::cout << "Pos midpoint is " << position_list_p[1][i] << std::endl;
+        MACH3LOG_TRACE("Pos midpoint is {:.2f}", position_list_p[1][i]);
       }
     }
     uncertainties_list.push_back({std::abs(position[i]-position_list[1][i]),std::abs(position[i]-position_list_p[1][i])});
-    std::cout << "Uncertainty finished for d = "<< i << std::endl;
-    std::cout << std::setprecision(10)<< "LLR values for ± positive and negative uncertainties are " << CalcChi(position_list[1]) << " and " << CalcChi(position_list_p[1]) << std::endl;
+    MACH3LOG_INFO("Uncertainty finished for d = {}", i);
+    MACH3LOG_INFO("LLR values for ± positive and negative uncertainties are {:<10.2f} and {:<10.2f}",
+                  CalcChi(position_list[1]), CalcChi(position_list_p[1]));
   }
   return uncertainties_list;
 }
 
-std::vector<std::vector<double>> PSO::calc_uncertainty(std::vector<double>position,double minimum) {
+// *************************
+std::vector<std::vector<double>> PSO::calc_uncertainty(std::vector<double>position, double minimum) {
+// *************************
   std::vector<double> pos_uncertainty(position.size());
   std::vector<double> neg_uncertainty(position.size());
   int num = 200;
@@ -255,7 +262,7 @@ std::vector<std::vector<double>> PSO::calc_uncertainty(std::vector<double>positi
       }
     }
     neg_uncertainty[i] = x[closest_index];
-    std::cout << "Neg" << std::endl;
+    MACH3LOG_INFO("Neg");
     x.assign(num, 0);
     y.assign(num, 0);
     StepPoint = (pos_stop-start) / (num - 1);
@@ -283,11 +290,13 @@ std::vector<std::vector<double>> PSO::calc_uncertainty(std::vector<double>positi
   return res;
 }
 
+// *************************
 void PSO::uncertainty_check(std::vector<double> previous_pos){
+// *************************
   std::vector<std::vector<double >> x_list;
   std::vector<std::vector<double >> y_list;
   std::vector<double> position = previous_pos;
-  int num = 5000;
+  constexpr int num = 5000;
   for (unsigned int i = 0;i<previous_pos.size();++i){
     double curr_ival = position[i];
     double start = previous_pos[i] - 1e-1;
@@ -296,30 +305,33 @@ void PSO::uncertainty_check(std::vector<double> previous_pos){
     std::vector<double> y(num);
     double StepPoint = (stop - start) / (num - 1);
     double value = start;
-    // std::cout << "result for fDim " << 1 << std::endl;
-    for (int j =0;j< num; ++j){
+    MACH3LOG_TRACE("result for fDim: {}", 1);
+    for (int j = 0;j < num; ++j) {
       position[i] = value;
       double LLR = CalcChi(position);
       x[j] = value;
       y[j] = LLR;
       value += StepPoint;
     }
     position[i] = curr_ival;
-    std::cout << " " << std::endl;
-    std::cout << "For fDim" << i+1  << " x list is " ;
-    for (unsigned int k= 0;k<x.size(); ++k){
-      std::cout << x[k] << " , " ;
-    } std::cout << " " << std::endl;
-    std::cout << "  " << std::endl;
-    std::cout << "For fDim" << i+1 << " y list is " ;
-    for (unsigned int k= 0;k<x.size(); ++k){
-      std::cout << y[k] << " , " ;
-    } std::cout <<  " " << std::endl;
-    std::cout << " " <<std::endl;
+    MACH3LOG_INFO("");
+    MACH3LOG_INFO("For fDim{} x list is", i+1);
+    for (unsigned int k = 0;k<x.size(); ++k){
+      MACH3LOG_INFO("  {}", x[k]);
+    }
+    MACH3LOG_INFO("");
+    MACH3LOG_INFO("For fDim{} y list is", i+1);
+    for (unsigned int k = 0; k < x.size(); ++k){
+      MACH3LOG_INFO("  {}", y[k]);
+    }
+    MACH3LOG_INFO("");
   }
 }
 
+// *************************
 double PSO::swarmIterate(){
+// *************************
+
   std::vector<double> total_pos(fDim,0.0);
 
   for (int i = 0; i < fParticles; ++i) {
@@ -377,8 +389,9 @@ double PSO::swarmIterate(){
   return mean_dist_sq;
 }
 
+// *************************
 void PSO::run() {
-
+// *************************
   double mean_dist_sq = 0;
 
   int iter = 0;
@@ -398,11 +411,11 @@ void PSO::run() {
     accCount++;
 
     if (i%100 == 0){
-      std::cout << "Mean Dist Sq = " << mean_dist_sq <<std::endl;
-      std::cout << "Current LLR = " << fBestValue << std::endl;
-      std::cout << "Position = " <<std::endl;
-      for (int j = 0; j< fDim; ++j){
-        std::cout << "    Dim " << j << " = " << std::setprecision(10) << get_best_particle()->get_personal_best_position()[j] << std::endl;
+      MACH3LOG_INFO("Mean Dist Sq = {:.2f}", mean_dist_sq);
+      MACH3LOG_INFO("Current LLR = {:.2f}", fBestValue);
+      MACH3LOG_INFO("Position:");
+      for (int j = 0; j < fDim; ++j){
+        MACH3LOG_INFO("    Dim {} = {:<10.2f}", j, get_best_particle()->get_personal_best_position()[j]);
       }
     }
     if(fConvergence > 0.0){
@@ -411,18 +424,20 @@ void PSO::run() {
       }
     }
   }
-  std::cout << "Finished after " << iter <<" runs out of "<< fIterations << std::endl;
-  std::cout << "Mean Dist " << mean_dist_sq <<std::endl;
-  std::cout << "Best LLR " << get_best_particle()->get_personal_best_value() << std::endl;
-
+  MACH3LOG_INFO("Finished after {} runs out of {}", iter, fIterations);
+  MACH3LOG_INFO("Mean Dist: {:.2f}", mean_dist_sq);
+  MACH3LOG_INFO("Best LLR: {:.2f}", get_best_particle()->get_personal_best_value());
   uncertainties = bisection(get_best_particle()->get_personal_best_position(),get_best_particle()->get_personal_best_value(),0.5,0.005);
-  std::cout << "Position for Global Minimum = "<<std::endl;
+  MACH3LOG_INFO("Position for Global Minimum = ");
   for (int i = 0; i< fDim; ++i){
-    std::cout << "    Dim " << i << " = " << std::setprecision(10) << get_best_particle()->get_personal_best_position()[i]  << " +" << uncertainties[i][1] << ", -" <<  uncertainties[i][0] << std::endl;
+    MACH3LOG_INFO("    Dim {} = {:<10.2f} +{:.2f}, -{:.2f}", i, get_best_particle()->get_personal_best_position()[i], uncertainties[i][1], uncertainties[i][0]);
   }
 }
 
+// *************************
 void PSO::WriteOutput(){
+// *************************
+
   outputFile->cd();
 
   TVectorD* PSOParValue = new TVectorD(fDim);

diff --git a/mcmc/PSO.h b/mcmc/PSO.h
@@ -1,9 +1,5 @@
 #pragma once
-//
-// Created by Emily Ip on 24/2/2023.
-//
-// Created by Emily Ip on 26/1/2023.
-//
+
 // C++ includes
 #include <algorithm>
 #include <functional>
@@ -15,6 +11,7 @@
 
 /// @brief Class particle - stores the position, velocity and personal best
 /// With functions which move particle and update velocity
+/// @note Created by Emily Ip on 24/2/2023.
 class particle{
   public:
       particle(){};
@@ -72,6 +69,7 @@ class particle{
  /// @brief Class PSO, consist of a vector of object Class Particle and global best
  /// Takes in the size (number of particle) and number of iteration
  /// functions includes: finding global best, updating velocity, actual minimisation function
+/// @note Created by Emily Ip on 24/2/2023.
 class PSO : public LikelihoodFit {
   public:
     /// @brief constructor