Use standard library math

doc/user_guide/fp_precision.md

@@ -13,7 +13,7 @@ BioDynaMo can be compiled with different floating point precisions.
 Currently, single-precision (`float`) and double-precision (`double`) are supported.
 By default, BioDynaMo is compiled with double-precision.
 
-Reduced floating-point precision reduces the required main memory, the file size of simulation backups, and might reduce the simulation runtime.
+Reduced floating-point precision reduces the required main memory and stored numerical data, and can reduce the simulation runtime.
 
 You can print the used precision with:
 

doc/user_guide/substance_initializers.md

@@ -114,9 +114,9 @@ struct GaussianBand {
 
   real_t operator()(real_t x, real_t y, real_t z) {
     switch(axis_) {
-      case Axis::kXAxis: return ROOT::Math::normal_pdf(x, sigma_, mean_);
-      case Axis::kYAxis: return ROOT::Math::normal_pdf(y, sigma_, mean_);
-      case Axis::kZAxis: return ROOT::Math::normal_pdf(z, sigma_, mean_);
+      case Axis::kXAxis: return Math::NormalPdf(x, sigma_, mean_);
+      case Axis::kYAxis: return Math::NormalPdf(y, sigma_, mean_);
+      case Axis::kZAxis: return Math::NormalPdf(z, sigma_, mean_);
       default: throw std::logic_error("You have chosen an non-existing axis!");
     }
   }
@@ -140,7 +140,7 @@ the following lambda:
 
 ```cpp
 auto gaussian_band = [](real_t x, real_t y, real_t z) {
-  return ROOT::Math::normal_pdf(x, 5, 0);
+  return Math::NormalPdf(x, 5, 0);
 };
 ```
 

src/core/container/math_array.h

@@ -16,6 +16,7 @@
 #define CORE_CONTAINER_MATH_ARRAY_H_
 
 #include <algorithm>
+#include <array>
 #include <cassert>
 #include <cmath>
 #include <numeric>
@@ -25,7 +26,6 @@
 
 #include "core/real_t.h"
 #include "core/util/log.h"
-#include "core/util/root.h"
 
 namespace bdm {
 
@@ -36,38 +36,31 @@ template <class T, std::size_t N>
 class MathArray {  // NOLINT
  public:
   /// Default constructor
-  MathArray() {
-#pragma omp simd
-    for (size_t i = 0; i < N; i++) {
-      data_[i] = T();
-    }
-  }
+  MathArray() = default;
 
   /// Constructor which accepts an std::initializer_list to set
   /// the array's content.
   /// \param l an initializer list
-  constexpr MathArray(std::initializer_list<T> l) {
-    assert(l.size() <= N);
-    auto it = l.begin();
-    for (uint64_t i = 0; i < N; i++) {
-      data_[i] = *(it++);
-    }
-    for (uint64_t i = l.size(); i < N; i++) {
-      data_[i] = T();
+  MathArray(std::initializer_list<T> l) {
+    if (l.size() > N) {
+      throw std::length_error("MathArray initializer exceeds its capacity");
     }
+    std::copy(l.begin(), l.end(), data_.begin());
   }
 
   /// Return a pointer to the underlying data.
   /// \return cont T pointer to the first entry of the array.
-  inline const T* data() const { return &data_[0]; }  // NOLINT
+  inline T* data() { return data_.data(); }  // NOLINT
+
+  inline const T* data() const { return data_.data(); }  // NOLINT
 
   /// Return the size of the array.
   /// \return integer denoting the array's size.
-  inline const size_t size() const { return N; }  // NOLINT
+  inline size_t size() const { return data_.size(); }  // NOLINT
 
   /// Check if the array is empty.
   /// \return true if size() == 0, false otherwise.
-  inline const bool empty() const { return N == 0; }  // NOLINT
+  inline bool empty() const { return data_.empty(); }  // NOLINT
 
   /// Overloaded array subscript operator. It does not perform
   /// any boundary checks.
@@ -85,58 +78,36 @@ class MathArray {  // NOLINT
   /// of the array's boundaries.
   /// \param idx the index of the element.
   /// \return the requested element.
-  T& at(size_t idx) noexcept(false) {  // NOLINT
-    if (idx > size() || idx < 0) {
-      throw std::out_of_range("The index is out of range");
-    }
-    return data_[idx];
-  }
+  T& at(size_t idx) { return data_.at(idx); }  // NOLINT
+
+  const T& at(size_t idx) const { return data_.at(idx); }  // NOLINT
 
-  const T* begin() const { return &(data_[0]); }  // NOLINT
+  const T* begin() const { return data_.begin(); }  // NOLINT
 
-  const T* end() const { return &(data_[N]); }  // NOLINT
+  const T* end() const { return data_.end(); }  // NOLINT
 
-  T* begin() { return &(data_[0]); }  // NOLINT
+  T* begin() { return data_.begin(); }  // NOLINT
 
-  T* end() { return &(data_[N]); }  // NOLINT
+  T* end() { return data_.end(); }  // NOLINT
 
   /// Returns the element at the beginning of the array.
   /// \return first element.
-  T& front() { return *(this->begin()); }  // NOLINT
+  T& front() { return data_.front(); }  // NOLINT
+
+  const T& front() const { return data_.front(); }  // NOLINT
 
   /// Return the element at the end of the array.
   /// \return last element.
-  T& back() {  // NOLINT
-    auto tmp = this->end();
-    tmp--;
-    return *tmp;
-  }
+  T& back() { return data_.back(); }  // NOLINT
 
-  /// Assignment operator.
-  /// \param other the other MathArray instance.
-  /// \return the current MathArray.
-  MathArray& operator=(const MathArray& other) {
-    if (this != &other) {
-      assert(other.size() == N);
-      std::copy(other.data_, other.data_ + other.size(), data_);
-    }
-    return *this;
-  }
+  const T& back() const { return data_.back(); }  // NOLINT
+
+  MathArray& operator=(const MathArray& other) = default;
 
   /// Equality operator.
   /// \param other a MathArray instance.
   /// \return true if they have the same content, false otherwise.
-  bool operator==(const MathArray& other) const {
-    if (other.size() != N) {
-      return false;
-    }
-    for (size_t i = 0; i < N; i++) {
-      if (other[i] != data_[i]) {
-        return false;
-      }
-    }
-    return true;
-  }
+  bool operator==(const MathArray& other) const { return data_ == other.data_; }
 
   bool operator!=(const MathArray& other) const { return !operator==(other); }
 
@@ -329,13 +300,13 @@ class MathArray {  // NOLINT
   /// \param k the constant value
   /// \return the array
   MathArray& fill(const T& k) {  // NOLINT
-    std::fill(std::begin(data_), std::end(data_), k);
+    data_.fill(k);
     return *this;
   }
 
   /// Return the sum of all the array's elements.
   /// \return sum of the array's content.
-  T Sum() const { return std::accumulate(begin(), end(), 0); }
+  T Sum() const { return std::accumulate(begin(), end(), T{}); }
 
   /// Checks if vector is a zero vector, e.g. if all entries are zero.
   bool IsZero() const {
@@ -404,8 +375,7 @@ class MathArray {  // NOLINT
   }
 
  private:
-  T data_[N];
-  BDM_CLASS_DEF_NV(MathArray, 1);  // NOLINT
+  std::array<T, N> data_{};
 };
 
 template <class T, std::size_t N>

src/core/randomized_rm.h

@@ -32,7 +32,6 @@ namespace bdm {
 template <typename TBaseRm>
 class RandomizedRm : public TBaseRm {
  public:
-  explicit RandomizedRm(TRootIOCtor* r) {}
   RandomizedRm(bool auto_randomize = true);
   virtual ~RandomizedRm();
 
@@ -43,7 +42,6 @@ class RandomizedRm : public TBaseRm {
   // Automatically randomize the agent order at the end of each iteration
   // If false, you can call RandomizeAgentsOrder() manually
   bool auto_randomize_ = true;
-  BDM_CLASS_DEF_NV(RandomizedRm, 1);
 };
 
 // -----------------------------------------------------------------------------
@@ -55,19 +53,6 @@ RandomizedRm<TBaseRm>::RandomizedRm(bool auto_randomize)
 template <typename TBaseRm>
 RandomizedRm<TBaseRm>::~RandomizedRm() = default;
 
-struct Ubrng {
-  using result_type = uint32_t;
-  Random* random;
-  Ubrng(Random* random) : random(random) {}
-  static constexpr result_type min() { return 0; }
-  static constexpr result_type max() {
-    return std::numeric_limits<result_type>::max();
-  }
-  result_type operator()() {
-    return random->Integer(std::numeric_limits<result_type>::max());
-  }
-};
-
 template <typename TBaseRm>
 void RandomizedRm<TBaseRm>::RandomizeAgentsOrder() {
   // shuffle
@@ -78,8 +63,7 @@ void RandomizedRm<TBaseRm>::RandomizeAgentsOrder() {
                                    this->agents_[n].end());
 #else
     auto* random = Simulation::GetActive()->GetRandom();
-    std::shuffle(this->agents_[n].begin(), this->agents_[n].end(),
-                 Ubrng(random));
+    std::shuffle(this->agents_[n].begin(), this->agents_[n].end(), *random);
 #endif  // LINUX
   }
 

src/core/substance_initializers.h

@@ -18,9 +18,8 @@
 #include <stdexcept>
 #include <vector>
 
-#include "Math/DistFunc.h"
-
 #include "core/diffusion/diffusion_grid.h"
+#include "core/util/math.h"
 
 namespace bdm {
 
@@ -77,8 +76,7 @@ class Uniform {
 };
 
 /// An initializer that follows a Gaussian (normal) distribution along one axis
-/// We use ROOT's built-in statistics function `normal_pdf(X, sigma, mean)`,
-/// that follows the normal probability density function:
+/// This initializer follows the normal probability density function:
 /// ( 1/( sigma * sqrt(2*pi) ))*e^( (-(x - mean )^2) / (2*sigma^2))
 class GaussianBand {
  public:
@@ -108,11 +106,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_ * Math::NormalPdf(x, sigma_, mean_);
       case Axis::kYAxis:
-        return scaling_ * ROOT::Math::normal_pdf(y, sigma_, mean_);
+        return scaling_ * Math::NormalPdf(y, sigma_, mean_);
       case Axis::kZAxis: {
-        return scaling_ * ROOT::Math::normal_pdf(z, sigma_, mean_);
+        return scaling_ * Math::NormalPdf(z, sigma_, mean_);
       }
       default:
         throw std::logic_error("You have chosen an non-existing axis!");
@@ -126,9 +124,8 @@ class GaussianBand {
   uint8_t axis_;
 };
 
-/// An initializer that follows a Poisson (normal) distribution along one axis
-/// The function ROOT::Math::poisson_pdfd(X, lambda) follows the normal
-/// probability density function:
+/// An initializer that follows a Poisson distribution along one axis.
+/// The Poisson probability mass function is:
 /// {e^( - lambda ) * lambda ^x )} / x!
 class PoissonBand {
  public:
@@ -153,11 +150,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 Math::PoissonPmf(static_cast<uint64_t>(x), lambda_);
       case Axis::kYAxis:
-        return ROOT::Math::poisson_pdf(y, lambda_);
+        return Math::PoissonPmf(static_cast<uint64_t>(y), lambda_);
       case Axis::kZAxis:
-        return ROOT::Math::poisson_pdf(z, lambda_);
+        return Math::PoissonPmf(static_cast<uint64_t>(z), lambda_);
       default:
         throw std::logic_error("You have chosen an non-existing axis!");
     }

src/core/util/math.h

@@ -17,26 +17,55 @@
 
 #include <array>
 #include <cmath>
+#include <limits>
 #include <numeric>
+#include <stdexcept>
 #include <vector>
 
-#include <TMath.h>
-
 #include "core/container/math_array.h"
 #include "core/util/log.h"
 #include "core/util/random.h"
 
 namespace bdm {
 
 struct Math {
-  /// value of pi
-  static constexpr real_t kPi = static_cast<real_t>(TMath::Pi());
+  // This long literal preserves pi when real_t is configured as double.
+  static constexpr real_t kPi =
+      static_cast<real_t>(3.141592653589793238462643383279502884L);
   /// Helpful constant to identify 'infinity'
-  static constexpr real_t kInfinity = 1e20;
+  static constexpr real_t kInfinity = std::numeric_limits<real_t>::infinity();
 
   static real_t ToDegree(real_t rad) { return rad * (180 / kPi); }
   static real_t ToRadian(real_t deg) { return deg * (kPi / 180); }
 
+  static real_t NormalPdf(real_t value, real_t sigma, real_t mean) {
+    if (sigma <= 0) {
+      throw std::invalid_argument("normal distribution sigma must be positive");
+    }
+    const auto normalized = (value - mean) / sigma;
+    return std::exp(-0.5 * normalized * normalized) /
+           (sigma * std::sqrt(2 * kPi));
+  }
+
+  static real_t NormalCdf(real_t value, real_t sigma, real_t mean) {
+    if (sigma <= 0) {
+      throw std::invalid_argument("normal distribution sigma must be positive");
+    }
+    return 0.5 *
+           (1 + std::erf((value - mean) / (sigma * std::sqrt(real_t{2}))));
+  }
+
+  static real_t PoissonPmf(uint64_t value, real_t mean) {
+    if (mean < 0) {
+      throw std::invalid_argument(
+          "Poisson distribution mean must be nonnegative");
+    }
+    if (mean == 0) {
+      return value == 0 ? 1 : 0;
+    }
+    return std::exp(value * std::log(mean) - mean - std::lgamma(value + 1));
+  }
+
   // Helper function that returns distance (L2 norm) between two positions in 3D
   static real_t GetL2Distance(const Real3& pos1, const Real3& pos2) {
     Real3 dist_array;