Opened by mistake, please disregard

CMakeLists.txt

@@ -153,6 +153,7 @@ option(valgrind  "Enable valgrind tests and make build compatible with valgrind
 option(rpath     "Link libraries with built-in RPATH (run-time search path)." OFF)
 option(real_t    "Define data type for real numbers. Currently supported: float, double" double)
 option(boost     "Enable Boost-based behaviors (e.g. RegulatoryNetwork)." ON)
+option(BDM_ANY_PYTHON "Accept any Python >= 3.9 (needed when ROOT does not come from the pinned prebuilt bundle)" OFF)
 
 if(APPLE)
   # ParaView on Apple devices
@@ -244,7 +245,11 @@ else()
       # Same Python ABI as the reused Ubuntu 24.04 archives; supports OpenSSL 3.
       set(pythonvers "3.9.25")
    endif()
-   find_package(Python ${pythonvers} EXACT COMPONENTS Interpreter Development)
+   if(BDM_ANY_PYTHON)
+     find_package(Python 3.9 COMPONENTS Interpreter Development)
+   else()
+     find_package(Python ${pythonvers} EXACT COMPONENTS Interpreter Development)
+   endif()
    if (NOT Python_FOUND)
       message(FATAL_ERROR "We did not find the Python ${pythonvers} interpreter from PyEnv. "
          "Make sure you executed \"pyenv shell ${pythonvers}\" before running cmake. "

src/core/model_initializer.h

@@ -15,7 +15,6 @@
 #ifndef CORE_MODEL_INITIALIZER_H_
 #define CORE_MODEL_INITIALIZER_H_
 
-#include <Math/DistFunc.h>
 #include <omp.h>
 #include <ctime>
 #include <string>

src/core/multi_simulation/experiment.h

@@ -18,7 +18,7 @@
 #include <functional>
 #include <vector>
 
-#include "TMath.h"
+#include <numeric>
 
 #include "core/analysis/time_series.h"
 #include "core/functor.h"
@@ -63,7 +63,8 @@ inline real_t Experiment(
                     [](const std::vector<real_t> all_y_values, real_t* y,
                        real_t* eh, real_t* el) {
                       *y =
-                          TMath::Mean(all_y_values.begin(), all_y_values.end());
+                          std::accumulate(all_y_values.begin(), all_y_values.end(), 0.0) /
+                              static_cast<double>(all_y_values.size());
                     });
 
   // Execute post-simulation lambda (e.g. plotting or exporting of simulated

src/core/param/command_line_options.cc

@@ -146,19 +146,19 @@ void CommandLineOptions::HandleCoreOptions() {
   // Handle "verbose" argument
   // If set in etc/bdm.rootrc use that value, command line argument will
   // override it
-  Int_t ll = kWarning;
+  int ll = kLogWarning;
   TString slevel = "Warning";
   TEnvRec* rec = gEnv->Lookup("Root.ErrorIgnoreLevel");
   if (rec) {
     if (rec->GetLevel() == kEnvUser)
       slevel = rec->GetValue();
   }
   if (!slevel.CompareTo("Print", TString::kIgnoreCase))
-    ll = kPrint;
+    ll = kLogPrint;
   else if (!slevel.CompareTo("Info", TString::kIgnoreCase))
-    ll = kInfo;
+    ll = kLogInfo;
   else if (!slevel.CompareTo("Warning", TString::kIgnoreCase))
-    ll = kWarning;
+    ll = kLogWarning;
   else if (!slevel.CompareTo("Error", TString::kIgnoreCase))
     ll = kError;
 
@@ -168,21 +168,22 @@ void CommandLineOptions::HandleCoreOptions() {
     switch (verbosity) {
       // case 0 can never occur; we wouldn't go into this if statement
       case 1:
-        ll = kWarning;
+        ll = kLogWarning;
         break;
       case 2:
-        ll = kInfo;
+        ll = kLogInfo;
         break;
       case 3:
-        ll = kPrint;
+        ll = kLogPrint;
         break;
       default:
-        ll = kPrint;
+        ll = kLogPrint;
         break;
     }
   }
   // Global variable of ROOT that determines verbosity of logging functions
-  gErrorIgnoreLevel = ll;
+  Log::SetLevel(ll);
+  gErrorIgnoreLevel = ll;  // ROOT's own messages, until ROOT is gone
 
   if (parser_->count("output-default-json")) {
     Param param;

src/core/substance_initializers.h

@@ -20,6 +20,7 @@
 
 #include "Math/DistFunc.h"
 
+#include "core/util/math.h"
 #include "core/diffusion/diffusion_grid.h"
 
 namespace bdm {
@@ -108,11 +109,11 @@ class GaussianBand {
   real_t operator()(real_t x, real_t y, real_t z) {
     switch (axis_) {
       case Axis::kXAxis:
-        return scaling_ * ROOT::Math::normal_pdf(x, sigma_, mean_);
+        return scaling_ * NormalPdf(x, sigma_, mean_);
       case Axis::kYAxis:
-        return scaling_ * ROOT::Math::normal_pdf(y, sigma_, mean_);
+        return scaling_ * NormalPdf(y, sigma_, mean_);
       case Axis::kZAxis: {
-        return scaling_ * ROOT::Math::normal_pdf(z, sigma_, mean_);
+        return scaling_ * NormalPdf(z, sigma_, mean_);
       }
       default:
         throw std::logic_error("You have chosen an non-existing axis!");
@@ -127,7 +128,7 @@ class GaussianBand {
 };
 
 /// An initializer that follows a Poisson (normal) distribution along one axis
-/// The function ROOT::Math::poisson_pdfd(X, lambda) follows the normal
+/// The function PoissonPdf(X, lambda) follows the normal
 /// probability density function:
 /// {e^( - lambda ) * lambda ^x )} / x!
 class PoissonBand {
@@ -153,11 +154,11 @@ class PoissonBand {
   real_t operator()(real_t x, real_t y, real_t z) {
     switch (axis_) {
       case Axis::kXAxis:
-        return ROOT::Math::poisson_pdf(x, lambda_);
+        return PoissonPdf(x, lambda_);
       case Axis::kYAxis:
-        return ROOT::Math::poisson_pdf(y, lambda_);
+        return PoissonPdf(y, lambda_);
       case Axis::kZAxis:
-        return ROOT::Math::poisson_pdf(z, lambda_);
+        return PoissonPdf(z, lambda_);
       default:
         throw std::logic_error("You have chosen an non-existing axis!");
     }

src/core/util/log.h

@@ -15,121 +15,110 @@
 #ifndef CORE_UTIL_LOG_H_
 #define CORE_UTIL_LOG_H_
 
-#include <TError.h>
-
+#include <atomic>
+#include <cerrno>
 #include <cstdio>
+#include <cstring>
 #include <cstdlib>
 #include <functional>
 #include <iostream>
+#include <mutex>
 #include <string>
 
 #include "core/util/string.h"
 
 namespace bdm {
 
-/// @brief Wrapper class over ROOT logging module
+/// Log levels. The numbers match ROOT's (kPrint, kInfo, ...) so existing
+/// verbosity settings keep their meaning; messages below the current level are dropped.
+enum LogLevel : int {
+  kLogPrint = 0,
+  kLogInfo = 1000,
+  kLogWarning = 2000,
+  kLogError = 3000,
+  kLogBreak = 4000,
+  kLogSysError = 5000,
+  kLogFatal = 6000
+};
+
+/// @brief Logging without ROOT. Output keeps ROOT's format: "<Level> in <location>: message"
 ///
 class Log {
  public:
-  /// @brief      Prints debug message
-  ///
-  /// @param[in]  location  The location of the message
-  /// @param[in]  parts     objects that compose the entire message
-  ///
+  /// Messages with a level below this are not printed (ROOT: gErrorIgnoreLevel).
+  static int GetLevel() { return LevelStorage().load(std::memory_order_relaxed); }
+  static void SetLevel(int level) { LevelStorage().store(level, std::memory_order_relaxed); }
+
   template <typename... Args>
   static void Debug(const std::string& location, const Args&... parts) {
-    // kPrint has the highest level of verbosity
-    if (gErrorIgnoreLevel <= kPrint) {
-      std::string message = Concat(parts...);
-      // Mimic ROOT logging output
-      fprintf(stderr, "Debug in <%s>: %s\n", location.c_str(), message.c_str());
+    // kLogPrint has the highest level of verbosity
+    if (GetLevel() <= kLogPrint) {
+      Print("Debug", location, Concat(parts...));
     }
   }
 
-  /// @brief      Prints information message
-  ///
-  /// @param[in]  location  The location of the message
-  /// @param[in]  parts     objects that compose the entire message
-  ///
   template <typename... Args>
   static void Info(const std::string& location, const Args&... parts) {
-    std::string message = Concat(parts...);
-    // ROOT function
-    ::Info(location.c_str(), "%s", message.c_str());
+    Emit(kLogInfo, "Info", location, Concat(parts...));
   }
 
-  /// @brief      Prints warning message
-  ///
-  /// @param[in]  location  The location of the message
-  /// @param[in]  parts     objects that compose the entire message
-  ///
   template <typename... Args>
   static void Warning(const std::string& location, const Args&... parts) {
-    std::string message = Concat(parts...);
-    // ROOT function
-    ::Warning(location.c_str(), "%s", message.c_str());
+    Emit(kLogWarning, "Warning", location, Concat(parts...));
   }
 
-  /// @brief Prints error message
-  ///
-  /// @param[in]  location  The location of the message
-  /// @param[in]  parts     objects that compose the entire message
-  ///
   template <typename... Args>
   static void Error(const std::string& location, const Args&... parts) {
-    std::string message = Concat(parts...);
-    // ROOT function
-    ::Error(location.c_str(), "%s", message.c_str());
+    Emit(kLogError, "Error", location, Concat(parts...));
   }
 
-  /// @brief Prints break message
-  ///
-  /// @param[in]  location  The location of the message
-  /// @param[in]  parts     objects that compose the entire message
-  ///
   template <typename... Args>
   static void Break(const std::string& location, const Args&... parts) {
-    std::string message = Concat(parts...);
-    // ROOT function
-    ::Break(location.c_str(), "%s", message.c_str());
+    Emit(kLogBreak, "Break", location, Concat(parts...));
   }
 
-  /// @brief Prints system error message
-  ///
-  /// @param[in]  location  The location of the message
-  /// @param[in]  parts     objects that compose the entire message
-  ///
   template <typename... Args>
   static void SysError(const std::string& location, const Args&... parts) {
-    std::string message = Concat(parts...);
-    // ROOT function
-    ::SysError(location.c_str(), "%s", message.c_str());
+    const int err = errno;
+    Emit(kLogSysError, "SysError", location,
+         Concat(parts...) + " (" + std::strerror(err) + ")");
   }
 
-  /// @brief Prints fatal error message
-  ///
-  /// @param[in]  location  The location of the message
-  /// @param[in]  parts     objects that compose the entire message
-  ///
   template <typename... Args>
   static void Fatal(const std::string& location, const Args&... parts) {
-    std::string message = Concat(parts...);
-    // ROOT function
-    ::Error(location.c_str(), "%s", message.c_str());
+    Emit(kLogError, "Error", location, Concat(parts...));
     // std::runtime_error() will fail the DeathTests
     exit(1);
   }
 
   template <typename... Args>
   static void Condition(const std::function<bool()>& lambda,
                         const Args&... parts) {
-    // kPrint has the highest level of verbosity
     if (lambda()) {
       std::string message = Concat(parts...);
-      // Mimic ROOT logging output
       fprintf(stdout, "%s\n", message.c_str());
     }
   }
+
+ private:
+  static std::atomic<int>& LevelStorage() {
+    static std::atomic<int> level{kLogPrint};
+    return level;
+  }
+  static void Emit(int level, const char* tag, const std::string& location,
+                   const std::string& message) {
+    if (level >= GetLevel()) {
+      Print(tag, location, message);
+    }
+  }
+  static void Print(const char* tag, const std::string& location,
+                    const std::string& message) {
+    static std::mutex mutex;
+    std::lock_guard<std::mutex> lock(mutex);
+    fflush(stdout);
+    fprintf(stderr, "%s in <%s>: %s\n", tag, location.c_str(), message.c_str());
+    fflush(stderr);
+  }
 };
 }  // namespace bdm
 

src/core/util/math.h

@@ -20,7 +20,7 @@
 #include <numeric>
 #include <vector>
 
-#include <TMath.h>
+#include <cmath>
 
 #include "core/container/math_array.h"
 #include "core/util/log.h"
@@ -30,7 +30,7 @@ namespace bdm {
 
 struct Math {
   /// value of pi
-  static constexpr real_t kPi = static_cast<real_t>(TMath::Pi());
+  static constexpr real_t kPi = static_cast<real_t>(3.14159265358979323846);
   /// Helpful constant to identify 'infinity'
   static constexpr real_t kInfinity = 1e20;
 
@@ -133,6 +133,27 @@ struct Math {
   }
 };
 
+/// Normal probability density, same definition as ROOT::Math::normal_pdf(x, sigma, mean).
+inline real_t NormalPdf(real_t x, real_t sigma, real_t mean) {
+  const real_t z = (x - mean) / sigma;
+  return std::exp(static_cast<real_t>(-0.5) * z * z) /
+         (std::abs(sigma) * std::sqrt(static_cast<real_t>(2) * Math::kPi));
+}
+
+/// Poisson probability of floor(x) events, as ROOT::Math::poisson_pdf(unsigned n, mu)
+/// for x >= 0. Returns 0 for x < 0, where the ROOT call converted a negative value
+/// to unsigned (undefined behaviour).
+inline real_t PoissonPdf(real_t x, real_t mu) {
+  if (x < 0) {
+    return 0;
+  }
+  const real_t n = std::floor(x);
+  if (n == 0) {
+    return std::exp(-mu);
+  }
+  return std::exp(n * std::log(mu) - std::lgamma(n + 1) - mu);
+}
+
 }  // namespace bdm
 
 #endif  // CORE_UTIL_MATH_H_