Merge pull request #683 from gundam-organization/main

update

cmake/options.cmake

@@ -43,6 +43,7 @@ elseif( CMAKE_BUILD_TYPE STREQUAL "RELEASE" )
   cmessage( STATUS "Build type manually specified to: RELEASE." )
 elseif( CMAKE_BUILD_TYPE STREQUAL "DEBUG" )
   cmessage( STATUS "Build type manually specified to: DEBUG." )
+  add_definitions( -D DEBUG_BUILD )
 else()
   cmessage( WARNING "Build type not recognised: ${CMAKE_BUILD_TYPE}. Using default build type: RELEASE." )
   set( CMAKE_BUILD_TYPE "RELEASE" )

resources/doc/guides/installOnLXPLUS.md

@@ -1,72 +1,122 @@
-[![](./resources/guides/images/lxplusLogo.png)]
+![](./images/lxplusLogo.png)
 
-This is a rough guide to get going on lxplus. GUNDAM has pretty modern dependencies so we’ll need to pick these up from what we can find on cvmfs etc. It turns out the ROOT v.18/04 we often use (the one inside the LArSoft cvmfs) is too old but we can use Adrien’s 6.26/10. 
 
-## Environment Setup
-Go to the directory you cloned GUNDAM into and create the following script, filling in the correct location for REPO_DIR:
+This is a rough guide to get going on LXPlus. 
+GUNDAM has pretty modern dependencies so we’ll need to pick these up
+from what we can find on cvmfs etc.
+LXPlus has recently switched towards "Red Hat Enterprise Linux release 9.4",
+aka el9, which changed the dependencies paths.
+This guide has been updated with the current software libraries available
+on the platform.
+
+
+## Common software available
+
+Common software libraries can be found under `/cvmfs/sft.cern.ch/lcg/views/`.
+They are documented on https://lcginfo.cern.ch/. 
+To this date, the latest LCG release is 106 (Jun 19, 2024).
+We can set up this environment by sourcing the dedicated bash file:
+
+```bash
+source /cvmfs/sft.cern.ch/lcg/views/LCG_106/x86_64-el9-gcc13-opt/setup.sh
 ```
-export REPO_DIR="/path/to/dir/above/GUNDAM"
-export BUILD_DIR="$REPO_DIR/gundam/build"
-export INSTALL_DIR="$REPO_DIR/gundam/install"
 
-mkdir -p $BUILD_DIR
+For the rest of the guide we assume the software files are organised in 
+three directories: `$REPO_DIR`, `$BUILD_DIR` and `$INSTALL_DIR`.
+
+For instance in my configuration on LXPlus is:
+```bash
+export REPO_DIR=/eos/home-a/adblanch/software/repo
+export BUILD_DIR=/eos/home-a/adblanch/software/build
+export INSTALL_DIR=/eos/home-a/adblanch/software/install
+```
+
+You should create those directories before going any further in the guide:
+
+```bash
 mkdir -p $REPO_DIR
+mkdir -p $BUILD_DIR
 mkdir -p $INSTALL_DIR
+```
 
-# enabling a C++11 compatible compiler (i.e. GCC5+)
-source /cvmfs/sft.cern.ch/lcg/contrib/gcc/11/x86_64-centos7/setup.sh
 
-# ROOT
-export PATH="/afs/cern.ch/user/a/adblanch/public/software/install/root/6.26.10_x86_64_el7_gcc11_cxx17/bin:$PATH"
-export LD_LIBRARY_PATH="/afs/cern.ch/user/a/adblanch/public/software/install/root/6.26.10_x86_64_el7_gcc11_cxx17/lib:$LD_LIBRARY_PATH"
+## GUNDAM Installation
 
-# git with proper submodule support
-export PATH="/cvmfs/sft.cern.ch/lcg/contrib/git/2.28.0/x86_64-centos7/bin:$PATH"
-export LD_LIBRARY_PATH="/cvmfs/sft.cern.ch/lcg/contrib/git/2.28.0/x86_64-centos7/lib:$LD_LIBRARY_PATH"
+### Clone the source code
 
-# JSON (my own compiled version, couldn't find a suitable one)
-export PATH="/afs/cern.ch/user/a/adblanch/public/software/install/json/bin:$PATH"
-export LD_LIBRARY_PATH="/afs/cern.ch/user/a/adblanch/public/software/install/json/lib:$LD_LIBRARY_PATH"
+First we need to clone the repository:
 
-# yaml-cpp
-export PATH="/cvmfs/sft.cern.ch/lcg/releases/yamlcpp/0.6.3-d05b2/x86_64-centos7-gcc11-opt/bin:$PATH"
-export LD_LIBRARY_PATH="/cvmfs/sft.cern.ch/lcg/releases/yamlcpp/0.6.3-d05b2/x86_64-centos7-gcc11-opt/lib:$LD_LIBRARY_PATH"
+```bash
+cd $REPO_DIR
+git clone https://github.com/gundam-organization/gundam.git
+cd gundam
+```
 
-# zlib -> optional dependency of GUNDAM
-export PATH="/cvmfs/sft.cern.ch/lcg/releases/zlib/1.2.11-8af4c/x86_64-centos7-gcc11-opt//bin:$PATH"
-export LD_LIBRARY_PATH="/cvmfs/sft.cern.ch/lcg/releases/lzlib/1.2.11-8af4c/x86_64-centos7-gcc11-opt//lib:$LD_LIBRARY_PATH"
+Now need to check out a given version of GUNDAM and pull the associated submodules
+(sub git repositories). Here is an example on how to pull the latest tagged version:
 
-#setup the fitter (this won't work until we build it!):
-source $BUILD_DIR/Linux/setup.sh
+```bash
+./update.sh --latest
 ```
 
-Then source the script you just made.
+Then create the dedicated build and install directory for GUNDAM:
 
-## Building Gundam
+```bash
+mkdir $BUILD_DIR/gundam
+mkdir $INSTALL_DIR/gundam
+```
 
-To build, do the following:
+Now GUNDAM should be ready to build.
 
-```
-cd $BUILD_DIR
-cmake3 ../
-make -j4
+
+### Build the apps and libraries
+
+On LXPlus, the CUDA libraries are available which means we can be enabled the support for
+GPU.
+
+```bash
+cd $BUILD_DIR/gundam
+
+# remove the cache is you already built
+rm $BUILD_DIR/gundam/CMakeCache.txt
+
+CLUSTER_OPTIONS="-D USE_STATIC_LINKS=ON"
+CLUSTER_OPTIONS="$CLUSTER_OPTIONS -D WITH_CUDA_LIB=ON"
+CLUSTER_OPTIONS="$CLUSTER_OPTIONS -D CMAKE_CUDA_COMPILER=/cvmfs/sft.cern.ch/lcg/views/LCG_106_cuda/x86_64-el9-gcc11-opt/bin/nvcc"
+CLUSTER_OPTIONS="$CLUSTER_OPTIONS -D CMAKE_CUDA_ARCHITECTURES=all"
+
+cmake \
+      -D CMAKE_INSTALL_PREFIX:PATH=$INSTALL_DIR/gundam \
+      -D CMAKE_BUILD_TYPE=RELEASE \
+      $CLUSTER_OPTIONS \
+      $REPO_DIR/gundam/.
+
 make install
 ```
 
-Now let's re-source the setup script you just made or just run the last line yourself: 
 
+### Setup the environment
+
+Once the make command succeeded, we have to make the gundam bin and lib directories 
+referenced within the environment.
+Put these lines in your `.bash_profile`:
+
+```bash
+echo "Setting up env for GUNDAM"
+export PATH="$INSTALL_DIR/gundam/bin:$PATH"
+export LD_LIBRARY_PATH="$INSTALL_DIR/gundam/lib:$LD_LIBRARY_PATH"
 ```
-source $BUILD_DIR/Linux/setup.sh
+
+```bash
+source $HOME/.bash_profile
 ```
 
-Now you should be able to run your first gundam command:
+You should be able to run the gundam apps:
 
 ```sh
 gundamFitter
 ```
 
-The program should stop while complaining that you didn't provide a config
-file:
 
 ```
 11:59:42 [gundamFitter.cxx]: Usage:

src/DialDictionary/DialEngine/src/DialInputBuffer.cpp

@@ -5,6 +5,7 @@
 #include "DialInputBuffer.h"
 
 #include "Logger.h"
+#include "GundamBacktrace.h"
 
 
 #ifndef DISABLE_USER_HEADER
@@ -71,6 +72,16 @@ void DialInputBuffer::update(){
       // re-apply the offset
       tempBuffer += inputRef.mirrorEdges.minValue;
     }
+    if( std::isnan(tempBuffer) ){
+        // LogThrowIf is broken, but OK for real error traps, but this is
+        // checking user input it's critical that the error message is
+        // properly formated so print an error, a backtrace, and then exit.
+        LogError << "NaN while evaluating input buffer of "
+                 << inputRef.getParameter(_parSetListPtr_).getTitle()
+                 << std::endl;
+        LogError << GundamUtils::Backtrace << std::endl;
+        std::exit(EXIT_FAILURE);
+    }
 
     // has it been updated?
     if( _inputBuffer_[inputRef.bufferIndex] != tempBuffer ){

src/Fitter/Minimizer/src/MinimizerBase.cpp

@@ -162,13 +162,13 @@ double MinimizerBase::evalFit( const double* parArray_ ){
   if (_useNormalizedFitSpace_) {
     const double* v = parArray_;
     for( auto* par : _minimizerParameterPtrList_ ){
-      par->setParameterValue(ParameterSet::toRealParValue(*(v++),*par));
+      par->setParameterValue(ParameterSet::toRealParValue(*(v++),*par), true);
     }
   }
   else {
     const double* v = parArray_;
     for( auto* par : _minimizerParameterPtrList_ ){
-      par->setParameterValue(*(v++));
+      par->setParameterValue(*(v++), true);
     }
   }
 

src/ParametersManager/src/Parameter.cpp

@@ -66,25 +66,30 @@ void Parameter::setMaxMirror(double maxMirror) {
   _mirrorRange_.max = maxMirror;
 }
 void Parameter::setParameterValue(double parameterValue, bool force) {
+#ifdef DEBUG_BUILD
+  // those printouts will only show if the CMAKE_BUILD_TYPE is set to DEBUG
   if (not isInDomain(parameterValue, true)) {
     LogError << "New parameter value is not in domain: " << parameterValue
              << std::endl;
     LogError << GundamUtils::Backtrace;
     if (not force) std::exit(EXIT_FAILURE);
     else LogAlert << "Forced continuation with invalid parameter" << std::endl;
   }
+#endif
   if( _parameterValue_ != parameterValue ){
     _gotUpdated_ = true;
     _parameterValue_ = parameterValue;
   }
   else{ _gotUpdated_ = false; }
 }
 double Parameter::getParameterValue() const {
+#ifdef DEBUG_BUILD
   if ( isEnabled() and not isValueWithinBounds() ) {
     LogWarning << "Getting out of bounds parameter: "
                << getSummary() << std::endl;
     LogDebug << GundamUtils::Backtrace;
   }
+#endif
   return _parameterValue_;
 }
 void Parameter::setDialSetConfig(const JsonType &jsonConfig_) {
@@ -112,8 +117,7 @@ bool Parameter::isInDomain(double value_, bool verbose_) const {
   }
   if ( not std::isnan(_parameterLimits_.min) and value_ < _parameterLimits_.min ) {
     if (verbose_) {
-      LogError << "Value is below minimum: " << value_
-               << std::endl;
+      LogError << "Value is below minimum: " << value_ << std::endl;
       LogError << "Summary: " << getSummary() << std::endl;
     }
     return false;

src/StatisticalInference/JointProbability/include/BarlowBeestonBanff2022.h

@@ -30,6 +30,23 @@ namespace JointProbability{
     bool allowZeroMcWhenZeroData{true};
     bool usePoissonLikelihood{false};
     bool BBNoUpdateWeights{false}; // OA 2021 bug reimplementation
+    // OA2021 and BANFF fractional error limitation is only relevent with
+    // BBNoUpdateWeights is true, and is needed to reproduce bugs when it is
+    // true.  When BBNoUpdateWeights is false, the fractional error will
+    // naturally be limited to less than 100%.  Physically, the fractional
+    // uncertainty should be less than 100% since one MC event in a bin would
+    // have 100% fractional uncertainty [under the "Gaussian" approximation,
+    // so sqrt(1.0)/1.0].  The OA2021 behavior lets the fractional error grow,
+    // but the entire likelihood became discontinuous around a predicted value
+    // of 1E-16. Setting fractionalErrorLimit to 1E+19 matches OA2021 before
+    // the discontinuity.  The BANFF behavior has the likelihood failed around
+    // 1E-154.  Setting fractionalErrorLimit to 1E+150 matchs BANFF.  In both
+    // cases, the new likelihood behaves reasonably all the way to zero.
+    double fractionalErrorLimit{1.0E+150};
+    // OA2021 bug reimplmentation (set to numeric_limits::min() to reproduce
+    // the bug).
+    double expectedValueMinimum{-1.};
+
     mutable std::map<const Sample*, std::vector<double>> nomMcUncertList{}; // OA 2021 bug reimplementation
     mutable GenericToolbox::NoCopyWrapper<std::mutex> _mutex_{}; // for creating the nomMC
   };
@@ -39,9 +56,17 @@ namespace JointProbability{
     GenericToolbox::Json::fillValue(_config_, allowZeroMcWhenZeroData, "allowZeroMcWhenZeroData");
     GenericToolbox::Json::fillValue(_config_, usePoissonLikelihood, "usePoissonLikelihood");
     GenericToolbox::Json::fillValue(_config_, BBNoUpdateWeights, "BBNoUpdateWeights");
+    GenericToolbox::Json::fillValue(_config_, fractionalErrorLimit, "fractionalErrorLimit");
+    GenericToolbox::Json::fillValue(_config_, expectedValueMinimum, "expectedValueMinimum");
     GenericToolbox::Json::fillValue(_config_, verboseLevel, {{"verboseLevel"},{"isVerbose"}});
     GenericToolbox::Json::fillValue(_config_, throwIfInfLlh, "throwIfInfLlh");
 
+    // Place a hard limit on the fractional error to prevent numeric issues.
+    if( fractionalErrorLimit > 1.0E+152 ){
+      LogAlert << "Placing a hard limit on the fractional error to prevent numeric issues." << std::endl;
+      fractionalErrorLimit = 1.0E+152;
+    }
+
   }
   double BarlowBeestonBanff2022::eval(const SamplePair& samplePair_, int bin_) const {
     double dataVal = samplePair_.data->getHistogram().getBinContentList()[bin_].sumWeights;
@@ -98,6 +123,19 @@ namespace JointProbability{
       }
     }
 
+    // Add the (broken) expected value threshold that we saw in the old
+    // likelihood. This is here to help understand the old behavior.
+    if (predVal <= expectedValueMinimum) {
+      LogAlert << "Use incorrect expectation behavior --"
+               << " predVal: " << predVal
+               << " dataVal: " << dataVal
+               << std::endl;
+      if (predVal <= 0.0 and dataVal > 0.0) {
+        return std::numeric_limits<double>::infinity();
+      }
+      return 0.0;
+    }
+
     // Implementing Barlow-Beeston correction for LH calculation. The
     // following comments are inspired/copied from Clarence's comments in the
     // MaCh3 implementation of the same feature
@@ -115,9 +153,19 @@ namespace JointProbability{
 
     // Barlow-Beeston uses fractional uncertainty on MC, so sqrt(sum[w^2])/mc
     // -b/2a in quadratic equation
-    if (mcuncert > std::numeric_limits<double>::epsilon()
-        and predVal > std::numeric_limits<double>::epsilon()) {
-      double fractional = sqrt(mcuncert) / predVal;
+    if ( mcuncert > 0.0 ) {
+      // Physically, the
+      // fractional uncertainty should be less than 100% since one MC event in
+      // a bin would have 100% fractional uncertainty [under the "Gaussian"
+      // approximation, so sqrt(1.0)/1.0].  The OA2021 behavior lets the
+      // fractional error grow, but the entire likelihood became discontinuous
+      // around a predicted value of 1E-16. Setting fractionalLimit to 1E+19
+      // matches OA2021 before the discontinuity.  The BANFF behavior has the
+      // likelihood failed around 1E-154.  Setting fractionalLimit to 1E+150
+      // matchs BANFF.  In both cases, the new likelihood behaves reasonably
+      // all the way to zero.
+      double fractional = std::min( std::sqrt(mcuncert) / predVal, fractionalErrorLimit);
+
       double temp = predVal * fractional * fractional - 1;
       // b^2 - 4ac in quadratic equation
       double temp2 = temp * temp + 4 * dataVal * fractional * fractional;
@@ -233,6 +281,8 @@ namespace JointProbability{
       LogInfo << GET_VAR_NAME_VALUE(BBNoUpdateWeights) << std::endl;
       LogInfo << GET_VAR_NAME_VALUE(verboseLevel) << std::endl;
       LogInfo << GET_VAR_NAME_VALUE(throwIfInfLlh) << std::endl;
+      LogInfo << GET_VAR_NAME_VALUE(expectedValueMinimum) << std::endl;
+      LogInfo << GET_VAR_NAME_VALUE(fractionalErrorLimit) << std::endl;
     }
 
   }

src/StatisticalInference/Likelihood/include/LikelihoodInterface.h

@@ -74,7 +74,7 @@ class LikelihoodInterface : public JsonBaseClass  {
   [[nodiscard]] const Propagator& getDataPropagator() const{ return _dataPropagator_; }
   [[nodiscard]] const PlotGenerator& getPlotGenerator() const{ return _plotGenerator_; }
   [[nodiscard]] const JointProbability::JointProbabilityBase* getJointProbabilityPtr() const { return _jointProbabilityPtr_.get(); }
-  [[nodiscard]] const std::vector<DatasetDefinition>& getDatasetList() const { return _dataSetList_; }
+  [[nodiscard]] const std::vector<DatasetDefinition>& getDatasetList() const { return _datasetList_; }
   [[nodiscard]] const std::vector<SamplePair>& getSamplePairList() const { return _samplePairList_; }
   [[nodiscard]] const Buffer& getBuffer() const { return _buffer_; }
 
@@ -83,7 +83,7 @@ class LikelihoodInterface : public JsonBaseClass  {
   Propagator& getModelPropagator(){ return _modelPropagator_; }
   Propagator& getDataPropagator(){ return _dataPropagator_; }
   PlotGenerator& getPlotGenerator(){ return _plotGenerator_; }
-  std::vector<DatasetDefinition>& getDatasetList(){ return _dataSetList_; }
+  std::vector<DatasetDefinition>& getDatasetList(){ return _datasetList_; }
   std::vector<SamplePair>& getSamplePairList(){ return _samplePairList_; }
 
   // mutable core
@@ -135,7 +135,7 @@ class LikelihoodInterface : public JsonBaseClass  {
   GenericToolbox::ParallelWorker _threadPool_{};
 
   /// user defined datasets
-  std::vector<DatasetDefinition> _dataSetList_;
+  std::vector<DatasetDefinition> _datasetList_;
 
   /// this is where model and data are kept to be compared
   Propagator _modelPropagator_{};