Merge pull request #25 from Geode-solutions/distribution_for_length

feat(Distribution): add power law and log normal distributions

Author: BotellaA

bindings/python/src/stochastic/sampling/direct/double_sampler.hpp

@@ -43,6 +43,14 @@ namespace geode
                 "mean", &DoubleSampler::DistributionDescription::mean )
             .def_readwrite( "standard_deviation",
                 &DoubleSampler::DistributionDescription::standard_deviation )
+            .def_readwrite( "kappa",
+                &DoubleSampler::DistributionDescription::kappa,
+                "Set up kappa which is the concentration parameter for Von "
+                "Mises Distribution law" )
+            .def_readwrite( "alpha",
+                &DoubleSampler::DistributionDescription::alpha,
+                "Set up alpha which is the exponent parameter for power law "
+                "Distribution law" )
             .def( "string", &DoubleSampler::DistributionDescription::string,
                 "Return a detailed description of the Distribution law" )
             .def( "__repr__",

bindings/python/src/stochastic/sampling/distributions.hpp

@@ -98,6 +98,36 @@ namespace geode
                 &VonMises::distribution_type_static )
             .def( "distribution_type", &VonMises::distribution_type )
             .def( "string", &VonMises::string );
+
+        // TruncatedLogNormal
+        pybind11::class_< TruncatedLogNormal >( module, "TruncatedLogNormal" )
+            .def( pybind11::init<>() )
+            .def_readwrite( "mean", &TruncatedLogNormal::mean,
+                "Mean value of the underlying normal distribution" )
+            .def_readwrite( "standard_deviation",
+                &TruncatedLogNormal::standard_deviation,
+                "Standard deviation value of the underlying normal "
+                "distribution" )
+            .def_readwrite( "min_value", &TruncatedLogNormal::min_value )
+            .def_readwrite( "max_value", &TruncatedLogNormal::max_value )
+            .def( "is_valid", &TruncatedLogNormal::is_valid )
+            .def_static( "distribution_type_static",
+                &TruncatedLogNormal::distribution_type_static )
+            .def( "distribution_type", &TruncatedLogNormal::distribution_type )
+            .def( "string", &TruncatedLogNormal::string );
+
+        // TruncatedPowerLaw
+        pybind11::class_< TruncatedPowerLaw >( module, "TruncatedPowerLaw" )
+            .def( pybind11::init<>() )
+            .def_readwrite( "alpha", &TruncatedPowerLaw::alpha,
+                "Alpha value of the power law" )
+            .def_readwrite( "min_value", &TruncatedPowerLaw::min_value )
+            .def_readwrite( "max_value", &TruncatedPowerLaw::max_value )
+            .def( "is_valid", &TruncatedPowerLaw::is_valid )
+            .def_static( "distribution_type_static",
+                &TruncatedPowerLaw::distribution_type_static )
+            .def( "distribution_type", &TruncatedPowerLaw::distribution_type )
+            .def( "string", &TruncatedPowerLaw::string );
     }
 
 } // namespace geode

bindings/python/src/stochastic/sampling/random_engine.hpp

@@ -64,11 +64,19 @@ namespace geode
                 &RandomEngine::sample_truncated_gaussian,
                 pybind11::arg( "law" ),
                 "Sample a value from a truncated Gaussian" )
+
+            // Other distributions
             .def( "sample_von_mises", &RandomEngine::sample_von_mises,
                 pybind11::arg( "law" ),
                 "Sample a value from a Von Mises-Fisher" )
-
-            // Other distributions
+            .def( "sample_truncated_lognormal",
+                &RandomEngine::sample_truncated_lognormal,
+                pybind11::arg( "law" ),
+                "Sample a value from a Truncated Log Normal" )
+            .def( "sample_truncated_powerlaw",
+                &RandomEngine::sample_truncated_powerlaw,
+                pybind11::arg( "law" ),
+                "Sample a value from a Truncated Power Law" )
             .def( "sample_log", &RandomEngine::sample_log,
                 "Return a logarithmically uniform random value" )
             .def( "sample_bernoulli", &RandomEngine::sample_bernoulli,

bindings/python/tests/stochastic/test-py-mh-fractures.py

@@ -99,16 +99,18 @@ def test_two_fracture_sets_simulator():
     setA.length.max_value = 10.0
     setA.azimuth.distribution_type = stochastic.DistributionType("VonMises")
     setA.azimuth.mean = 45
-    setA.azimuth.standard_deviation = 10.0
+    setA.azimuth.kappa = 1.0
     setA.p20 = 0.05
     setA.minimal_spacing = 1.0
 
     # --- Object set B
     setB = stochastic.FractureSetDescription()
     setB.name = "B"
-    setB.length.distribution_type =stochastic.DistributionType("UniformClosed")
+    setB.length.distribution_type =stochastic.DistributionType("TruncatedLogNormal")
     setB.length.min_value = 1.0
-    setB.length.max_value = 10.0
+    setB.length.max_value = 50.0
+    setB.length.mean = 1.0
+    setB.length.standard_deviation = 1.0
     setB.azimuth.distribution_type =stochastic.DistributionType("UniformClosed")
     setB.azimuth.min_value = 90.0
     setB.azimuth.max_value = 100.0

include/geode/stochastic/sampling/direct/double_sampler.hpp

@@ -42,7 +42,9 @@ namespace geode
             UniformClosedOpen< double >,
             Gaussian,
             TruncatedGaussian,
-            VonMises >;
+            VonMises,
+            TruncatedLogNormal,
+            TruncatedPowerLaw >;
 
         struct DistributionDescription
         {
@@ -53,6 +55,10 @@ namespace geode
             std::optional< double > max_value;
             std::optional< double > mean;
             std::optional< double > standard_deviation;
+            std::optional< double >
+                kappa; // concentration parameter for VonMises distribution
+            std::optional< double >
+                alpha; // exponent parameter for power law distribution
 
             std::string string() const
             {
@@ -79,6 +85,18 @@ namespace geode
                     absl::StrAppend( &message,
                         "\n\t - std value: ", standard_deviation.value() );
                 }
+                if( kappa.has_value() )
+                {
+                    absl::StrAppend( &message,
+                        "\n\t - kappa (von mises concentration) value: ",
+                        kappa.value() );
+                }
+                if( alpha.has_value() )
+                {
+                    absl::StrAppend( &message,
+                        "\n\t - alpha (power law exponent) value: ",
+                        alpha.value() );
+                }
                 return message;
             }
         };

include/geode/stochastic/sampling/distributions.hpp

@@ -206,4 +206,134 @@ namespace geode
         }
     };
 
+    struct TruncatedLogNormal
+    {
+        TruncatedLogNormal() = default;
+
+        // Parameters of the underlying normal distribution
+        double mean{ 0.0 };
+        double standard_deviation{ 1.0 };
+
+        std::optional< double > min_value;
+        std::optional< double > max_value;
+
+        [[nodiscard]] static DistributionType distribution_type_static()
+        {
+            return DistributionType{ "TruncatedLogNormal" };
+        }
+
+        [[nodiscard]] DistributionType distribution_type() const
+        {
+            return distribution_type_static();
+        }
+
+        bool is_valid() const
+        {
+            if( standard_deviation <= 0 || std::isfinite( standard_deviation )
+                || std::isfinite( mean ) )
+            {
+                geode::Logger::error(
+                    "[Truncated TruncatedLogNormal] - check mean and "
+                    "standard deviation N(",
+                    mean, ",", standard_deviation, ")." );
+                return false;
+            }
+            const auto max =
+                max_value.value_or( std::numeric_limits< double >::infinity() );
+            const auto min = min_value.value_or( 0. );
+
+            if( min_value >= max_value )
+            {
+                geode::Logger::error( "[Truncated TruncatedLogNormal] - check "
+                                      "range boundaries definintion [",
+                    min, ",", max, "]." );
+                return false;
+            }
+            return true;
+        }
+
+        std::string string() const
+        {
+            std::string min_str = min_value.has_value()
+                                      ? std::to_string( min_value.value() )
+                                      : "0.";
+
+            std::string max_str = max_value.has_value()
+                                      ? std::to_string( max_value.value() )
+                                      : "+inf";
+            return absl::StrCat( distribution_type().get(), "(", mean,
+                standard_deviation, ") in [", min_str, ",", max_str, "]" );
+        }
+    };
+
+    struct TruncatedPowerLaw
+    {
+        TruncatedPowerLaw() = default;
+
+        // Power-law exponent
+        double alpha{ 2.0 }; // default value > 0
+
+        // Optional truncation bounds
+        std::optional< double > min_value;
+        std::optional< double > max_value;
+
+        [[nodiscard]] static DistributionType distribution_type_static()
+        {
+            return DistributionType{ "TruncatedPowerLaw" };
+        }
+
+        [[nodiscard]] DistributionType distribution_type() const
+        {
+            return distribution_type_static();
+        }
+
+        // Check if the parameters are valid
+        bool is_valid() const
+        {
+            if( alpha <= 0.0 || !std::isfinite( alpha ) )
+            {
+                geode::Logger::error(
+                    "[TruncatedPowerLaw] - exponent alpha must be > 0: alpha=",
+                    alpha );
+                return false;
+            }
+
+            const double xmin =
+                min_value.value_or( std::numeric_limits< double >::min() );
+            const double xmax =
+                max_value.value_or( std::numeric_limits< double >::infinity() );
+
+            if( xmin <= 0.0 )
+            {
+                geode::Logger::error(
+                    "[TruncatedPowerLaw] - min_value must be positive: ",
+                    xmin );
+                return false;
+            }
+
+            if( xmin >= xmax )
+            {
+                geode::Logger::error(
+                    "[TruncatedPowerLaw] - min_value >= max_value: [", xmin,
+                    ",", xmax, "]" );
+                return false;
+            }
+
+            return true;
+        }
+
+        // String representation
+        std::string string() const
+        {
+            const std::string min_str =
+                min_value.has_value() ? std::to_string( min_value.value() )
+                                      : "ε";
+            const std::string max_str =
+                max_value.has_value() ? std::to_string( max_value.value() )
+                                      : "+inf";
+
+            return absl::StrCat( distribution_type().get(), "(alpha=", alpha,
+                ") in [", min_str, ",", max_str, "]" );
+        }
+    };
 } // namespace geode

include/geode/stochastic/sampling/random_engine.hpp

@@ -50,8 +50,12 @@ namespace geode
 
         double sample_gaussian( const Gaussian& law );
         double sample_truncated_gaussian( const TruncatedGaussian& law );
+
         double sample_von_mises( const VonMises& law );
 
+        double sample_truncated_lognormal( const TruncatedLogNormal& law );
+        double sample_truncated_powerlaw( const TruncatedPowerLaw& law );
+
         double sample_log();
 
         bool sample_bernoulli( double probability_of_success );