random.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// Copyright (C) 2021 CERN & University of Surrey for the benefit of the
// BioDynaMo collaboration. All Rights Reserved.
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
//
// See the LICENSE file distributed with this work for details.
// See the NOTICE file distributed with this work for additional information
// regarding copyright ownership.
//
// -----------------------------------------------------------------------------
 
#ifndef CORE_UTIL_RANDOM_H_
#define CORE_UTIL_RANDOM_H_
 
#include <unordered_map>
#include "core/container/fixed_size_vector.h"
#include "core/container/math_array.h"
#include "core/util/root.h"
 
class TRandom;
class TF1;
class TF2;
class TF3;
 
namespace bdm {
 
// -----------------------------------------------------------------------------
template <typename TSample>
class DistributionRng {
 public:
  DistributionRng() = default;
  DistributionRng(TRootIOCtor*) {}
  virtual ~DistributionRng() = default;
  TSample Sample();
  MathArray<TSample, 2> Sample2();
  MathArray<TSample, 3> Sample3();
 
  template <uint64_t N>
  MathArray<TSample, N> SampleArray() {
    MathArray<TSample, N> ret;
    for (uint64_t i = 0; i < N; i++) {
      ret[i] = Sample();
    }
    return ret;
  }
 
 protected:
  virtual TSample SampleImpl(TRandom* rng) = 0;
  virtual MathArray<TSample, 2> Sample2Impl(TRandom* rng);
  virtual MathArray<TSample, 3> Sample3Impl(TRandom* rng);
  BDM_CLASS_DEF(DistributionRng, 1);
};
 
// -----------------------------------------------------------------------------
class UniformRng : public DistributionRng<real_t> {
 public:
  UniformRng(real_t min, real_t max);
  ~UniformRng() override;
 
 private:
  real_t min_, max_;
  real_t SampleImpl(TRandom* rng) override;
  BDM_CLASS_DEF_OVERRIDE(UniformRng, 1);
};
 
// -----------------------------------------------------------------------------
class GausRng : public DistributionRng<real_t> {
 public:
  GausRng(real_t mean, real_t sigma);
  virtual ~GausRng() override;
 
 private:
  real_t mean_, sigma_;
  real_t SampleImpl(TRandom* rng) override;
  BDM_CLASS_DEF_OVERRIDE(GausRng, 1);
};
 
// -----------------------------------------------------------------------------
class ExpRng : public DistributionRng<real_t> {
 public:
  ExpRng(real_t tau);
  ~ExpRng() override;
 
 private:
  real_t tau_;
  real_t SampleImpl(TRandom* rng) override;
  BDM_CLASS_DEF_OVERRIDE(ExpRng, 1);
};
 
// -----------------------------------------------------------------------------
class LandauRng : public DistributionRng<real_t> {
 public:
  LandauRng(real_t mean, real_t sigma);
  ~LandauRng() override;
 
 private:
  real_t mean_, sigma_;
  real_t SampleImpl(TRandom* rng) override;
  BDM_CLASS_DEF_OVERRIDE(LandauRng, 1);
};
 
// -----------------------------------------------------------------------------
class PoissonDRng : public DistributionRng<real_t> {
 public:
  PoissonDRng(real_t mean);
  ~PoissonDRng() override;
 
 private:
  real_t mean_;
  real_t SampleImpl(TRandom* rng) override;
  BDM_CLASS_DEF_OVERRIDE(PoissonDRng, 1);
};
 
// -----------------------------------------------------------------------------
class BreitWignerRng : public DistributionRng<real_t> {
 public:
  BreitWignerRng(real_t mean, real_t gamma);
  ~BreitWignerRng() override;
 
 private:
  real_t mean_, gamma_;
  real_t SampleImpl(TRandom* rng) override;
  BDM_CLASS_DEF_OVERRIDE(BreitWignerRng, 1);
};
 
// -----------------------------------------------------------------------------
class UserDefinedDistRng1D : public DistributionRng<real_t> {
 public:
  UserDefinedDistRng1D(TRootIOCtor* ioctor) : DistributionRng<real_t>(ioctor) {}
  UserDefinedDistRng1D(TF1* function, const char* option);
  ~UserDefinedDistRng1D() override;
  void Draw(const char* option = "");
  TF1* GetTF1();
 
 private:
  real_t SampleImpl(TRandom* rng) override;
  // TODO use shared_ptr once ROOT supports IO of them
  TF1* function_ = nullptr;
  const char* option_ = nullptr;
  BDM_CLASS_DEF_OVERRIDE(UserDefinedDistRng1D, 1);
};
 
// -----------------------------------------------------------------------------
class UserDefinedDistRng2D : public DistributionRng<real_t> {
 public:
  UserDefinedDistRng2D(TRootIOCtor* ioctor) : DistributionRng<real_t>(ioctor) {}
  UserDefinedDistRng2D(TF2* function, const char* option);
  ~UserDefinedDistRng2D() override;
  void Draw(const char* option = "");
  TF2* GetTF2();
 
 private:
  real_t SampleImpl(TRandom* rng) override;
  MathArray<real_t, 2> Sample2Impl(TRandom* rng) override;
  // TODO use shared_ptr once ROOT supports IO of them
  TF2* function_ = nullptr;
  const char* option_ = nullptr;
  BDM_CLASS_DEF_OVERRIDE(UserDefinedDistRng2D, 1);
};
 
// -----------------------------------------------------------------------------
class UserDefinedDistRng3D : public DistributionRng<real_t> {
 public:
  UserDefinedDistRng3D(TRootIOCtor* ioctor) : DistributionRng<real_t>(ioctor) {}
  UserDefinedDistRng3D(TF3* function, const char* option);
  ~UserDefinedDistRng3D() override;
  void Draw(const char* option = "");
  TF3* GetTF3();
 
 private:
  real_t SampleImpl(TRandom* rng) override;
  MathArray<real_t, 2> Sample2Impl(TRandom* rng) override;
  MathArray<real_t, 3> Sample3Impl(TRandom* rng) override;
  // TODO use shared_ptr once ROOT supports IO of them
  TF3* function_ = nullptr;
  const char* option_ = nullptr;
  BDM_CLASS_DEF_OVERRIDE(UserDefinedDistRng3D, 1);
};
 
// -----------------------------------------------------------------------------
struct UserDefinedDist {
  double (*ud_function)(const double*, const double*) = nullptr;
  FixedSizeVector<real_t, 10> parameters;
  real_t xmin = 0;
  real_t xmax = 0;
  real_t ymin = 0;
  real_t ymax = 0;
  real_t zmin = 0;
  real_t zmax = 0;
  bool operator==(const UserDefinedDist& other) const {
    return (ud_function == other.ud_function) &&
           (parameters == other.parameters) && (xmin == other.xmin) &&
           (xmax == other.xmax) && (ymin == other.ymin) &&
           (ymax == other.ymax) && (zmin == other.zmin) && (zmax == other.zmax);
  }
};
 
}  // namespace bdm
namespace std {
 
template <>
struct hash<bdm::UserDefinedDist> {
  std::size_t operator()(const bdm::UserDefinedDist& uddist) const noexcept {
    std::size_t h = 7;
    h *= 31 * reinterpret_cast<std::size_t>(uddist.ud_function);
    for (auto el : uddist.parameters) {
      h *= 31 * std::size_t(el * 1000);
    }
    h *= 31 * std::size_t(uddist.xmin * 1000);
    h *= 31 * std::size_t(uddist.xmax * 1000);
    h *= 31 * std::size_t(uddist.ymin * 1000);
    h *= 31 * std::size_t(uddist.ymax * 1000);
    h *= 31 * std::size_t(uddist.zmin * 1000);
    h *= 31 * std::size_t(uddist.zmax * 1000);
    return h;
  }
};
 
}  // namespace std
namespace bdm {
 
// -----------------------------------------------------------------------------
class BinomialRng : public DistributionRng<int> {
 public:
  BinomialRng(int ntot, real_t prob);
  ~BinomialRng() override;
 
 private:
  int ntot_;
  real_t prob_;
  int SampleImpl(TRandom* rng) override;
  BDM_CLASS_DEF_OVERRIDE(BinomialRng, 1);
};
 
// -----------------------------------------------------------------------------
class PoissonRng : public DistributionRng<int> {
 public:
  PoissonRng(real_t mean);
  ~PoissonRng() override;
 
 private:
  real_t mean_;
  int SampleImpl(TRandom* rng) override;
  BDM_CLASS_DEF_OVERRIDE(PoissonRng, 1);
};
 
// -----------------------------------------------------------------------------
class Random {
 public:
  Random();
  explicit Random(TRootIOCtor*);
  Random(const Random& other);
  ~Random();
  Random& operator=(const Random& other);
 
  real_t Uniform(real_t max = 1.0);
  real_t Uniform(real_t min, real_t max);
 
  template <uint64_t N>
  MathArray<real_t, N> UniformArray(real_t max = 1.0) {
    MathArray<real_t, N> ret;
    for (uint64_t i = 0; i < N; i++) {
      ret[i] = Uniform(max);
    }
    return ret;
  }
 
  template <uint64_t N>
  MathArray<real_t, N> UniformArray(real_t min, real_t max) {
    MathArray<real_t, N> ret;
    for (uint64_t i = 0; i < N; i++) {
      ret[i] = Uniform(min, max);
    }
    return ret;
  }
 
  real_t Gaus(real_t mean = 0.0, real_t sigma = 1.0);
  real_t Exp(real_t tau);
  real_t Landau(real_t mean = 0, real_t sigma = 1);
  real_t PoissonD(real_t mean);
  real_t BreitWigner(real_t mean = 0, real_t gamma = 1);
 
  unsigned Integer(int max);
  int Binomial(int ntot, real_t prob);
  int Poisson(real_t mean);
 
  MathArray<real_t, 2> Circle(real_t radius);
  MathArray<real_t, 3> Sphere(real_t radius);
 
  void SetSeed(uint64_t seed);
 
  uint64_t GetSeed() const;
 
  void SetGenerator(TRandom* new_rng);
 
  UniformRng GetUniformRng(real_t min = 0, real_t max = 1) const;
  GausRng GetGausRng(real_t mean = 0, real_t sigma = 1) const;
  ExpRng GetExpRng(real_t tau) const;
  LandauRng GetLandauRng(real_t mean = 0, real_t sigma = 1) const;
  PoissonDRng GetPoissonDRng(real_t mean) const;
  BreitWignerRng GetBreitWignerRng(real_t mean = 0, real_t gamma = 1) const;
 
  UserDefinedDistRng1D GetUserDefinedDistRng1D(
      double (*f)(const double*, const double*),
      const FixedSizeVector<real_t, 10>& params, real_t min, real_t max,
      const char* option = nullptr);
 
  UserDefinedDistRng2D GetUserDefinedDistRng2D(
      double (*f)(const double*, const double*),
      const FixedSizeVector<real_t, 10>& params, real_t xmin, real_t xmax,
      real_t ymin, real_t ymax, const char* option = nullptr);
 
  UserDefinedDistRng3D GetUserDefinedDistRng3D(
      double (*f)(const double*, const double*),
      const FixedSizeVector<real_t, 10>& params, real_t xmin, real_t xmax,
      real_t ymin, real_t ymax, real_t zmin, real_t zmax,
      const char* option = nullptr);
 
  BinomialRng GetBinomialRng(int ntot, real_t prob) const;
  PoissonRng GetPoissonRng(real_t mean) const;
 
 private:
  friend class DistributionRng<real_t>;
  friend class DistributionRng<int>;
 
  TRandom* generator_ = nullptr;
  std::unordered_map<UserDefinedDist, TF1*> udd_tf1_map_;  
  std::unordered_map<UserDefinedDist, TF2*> udd_tf2_map_;  
  std::unordered_map<UserDefinedDist, TF3*> udd_tf3_map_;  
  BDM_CLASS_DEF_NV(Random, 3);
};
 
}  // namespace bdm
 
#endif  // CORE_UTIL_RANDOM_H_