From 25e5ecea9b52fc002e2771e1042499510d2c3af2 Mon Sep 17 00:00:00 2001 From: Francois Bonneau Date: Fri, 14 Nov 2025 20:18:09 +0100 Subject: [PATCH 1/4] feat(Distribution): add power law and log normal distributions --- .../sampling/direct/double_sampler.hpp | 8 + .../src/stochastic/sampling/distributions.hpp | 30 +++ .../src/stochastic/sampling/random_engine.hpp | 12 +- .../tests/stochastic/test-py-mh-fractures.py | 8 +- .../sampling/direct/double_sampler.hpp | 20 +- .../stochastic/sampling/distributions.hpp | 130 +++++++++ .../stochastic/sampling/random_engine.hpp | 4 + .../sampling/direct/double_sampler.cpp | 146 +++++----- .../stochastic/sampling/random_engine.cpp | 253 +++++++++++++++--- .../sampling/mcmc/test-mh-fractures.cpp | 9 +- 10 files changed, 499 insertions(+), 121 deletions(-) diff --git a/bindings/python/src/stochastic/sampling/direct/double_sampler.hpp b/bindings/python/src/stochastic/sampling/direct/double_sampler.hpp index 8cddc09..c4461a9 100644 --- a/bindings/python/src/stochastic/sampling/direct/double_sampler.hpp +++ b/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__", diff --git a/bindings/python/src/stochastic/sampling/distributions.hpp b/bindings/python/src/stochastic/sampling/distributions.hpp index 54f979a..d1acdf8 100644 --- a/bindings/python/src/stochastic/sampling/distributions.hpp +++ b/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 diff --git a/bindings/python/src/stochastic/sampling/random_engine.hpp b/bindings/python/src/stochastic/sampling/random_engine.hpp index 546decf..c91e7b9 100644 --- a/bindings/python/src/stochastic/sampling/random_engine.hpp +++ b/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, diff --git a/bindings/python/tests/stochastic/test-py-mh-fractures.py b/bindings/python/tests/stochastic/test-py-mh-fractures.py index 8f979b4..f7d466a 100644 --- a/bindings/python/tests/stochastic/test-py-mh-fractures.py +++ b/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 diff --git a/include/geode/stochastic/sampling/direct/double_sampler.hpp b/include/geode/stochastic/sampling/direct/double_sampler.hpp index 5541fab..31f5c52 100644 --- a/include/geode/stochastic/sampling/direct/double_sampler.hpp +++ b/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; } }; diff --git a/include/geode/stochastic/sampling/distributions.hpp b/include/geode/stochastic/sampling/distributions.hpp index c2e3b56..3c808f1 100644 --- a/include/geode/stochastic/sampling/distributions.hpp +++ b/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 \ No newline at end of file diff --git a/include/geode/stochastic/sampling/random_engine.hpp b/include/geode/stochastic/sampling/random_engine.hpp index f0825eb..a6c0b5b 100644 --- a/include/geode/stochastic/sampling/random_engine.hpp +++ b/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 ); diff --git a/src/geode/stochastic/sampling/direct/double_sampler.cpp b/src/geode/stochastic/sampling/direct/double_sampler.cpp index fc8ac6c..36766f8 100644 --- a/src/geode/stochastic/sampling/direct/double_sampler.cpp +++ b/src/geode/stochastic/sampling/direct/double_sampler.cpp @@ -28,59 +28,6 @@ #include #include -namespace -{ - struct DistributionParams - { - double min; - double max; - double mean; - double std; - }; - - DistributionParams complete_distribution_params( - const geode::DoubleSampler::DistributionDescription& d ) - { - DistributionParams p{}; - - if( d.min_value && d.max_value ) - { - p.min = *d.min_value; - p.max = *d.max_value; - p.mean = ( p.min + p.max ) / 2.0; - p.std = ( p.max - p.min ) / std::sqrt( 12.0 ); - return p; - } - - if( d.mean && d.standard_deviation ) - { - p.mean = *d.mean; - p.std = *d.standard_deviation; - p.min = p.mean - std::sqrt( 3.0 ) * p.std; - p.max = p.mean + std::sqrt( 3.0 ) * p.std; - return p; - } - - throw geode::OpenGeodeException( - "[DistributionDescripption] Incomplete distribution description: " - "need at least (min,max) or (mean,std)." ); - return p; - } - - double compute_kappa_von_mises( double standard_deviation_rad ) - { - geode::Logger::info( - "[VonMises] approximate concentation (kappa) from with " - "1/(std*std)." ); - OPENGEODE_EXCEPTION( standard_deviation_rad > geode::GLOBAL_EPSILON, - "Cannot evaluate the VonMises concentration since standard " - "deviation is equal to ", - standard_deviation_rad, " (in rad)." ); - - return 1.0 / ( standard_deviation_rad * standard_deviation_rad ); - } - -} // namespace namespace geode { struct DistributionTypeHasher @@ -101,48 +48,87 @@ namespace geode > distribution_registry = { { UniformClosed< double >::distribution_type_static(), - []( const DoubleSampler::DistributionDescription& d ) { - auto p = complete_distribution_params( d ); + []( const DoubleSampler::DistributionDescription& desc ) { + OPENGEODE_EXCEPTION( desc.min_value && desc.max_value, + "[DoubleSampler] - Incomplete description for " + "Uniform distribution need at least min " + "and max values" ); UniformClosed< double > dist; - dist.min_value = p.min; - dist.max_value = p.max; + dist.min_value = desc.min_value.value(); + dist.max_value = desc.max_value.value(); return dist; } }, { UniformClosedOpen< double >::distribution_type_static(), - []( const DoubleSampler::DistributionDescription& d ) { - auto p = complete_distribution_params( d ); + []( const DoubleSampler::DistributionDescription& desc ) { + OPENGEODE_EXCEPTION( desc.min_value && desc.max_value, + "[DoubleSampler] - Incomplete description for " + "Uniform distribution need at least min " + "and max values" ); UniformClosedOpen< double > dist; - dist.min_value = p.min; - dist.max_value = p.max; + dist.min_value = desc.min_value.value(); + dist.max_value = desc.max_value.value(); return dist; } }, { Gaussian::distribution_type_static(), - []( const DoubleSampler::DistributionDescription& d ) { - auto p = complete_distribution_params( d ); + []( const DoubleSampler::DistributionDescription& desc ) { + OPENGEODE_EXCEPTION( desc.mean && desc.standard_deviation, + "[DoubleSampler] - Incomplete description for " + "Gaussian distribution need at least mean " + "and standard deviation values" ); Gaussian dist; - dist.mean = p.mean; - dist.standard_deviation = p.std; + dist.mean = desc.mean.value(); + dist.standard_deviation = desc.standard_deviation.value(); return dist; } }, { TruncatedGaussian::distribution_type_static(), - []( const DoubleSampler::DistributionDescription& d ) { - auto p = complete_distribution_params( d ); + []( const DoubleSampler::DistributionDescription& desc ) { + OPENGEODE_EXCEPTION( desc.mean && desc.standard_deviation, + "[DoubleSampler] - Incomplete description for " + "Truncated Gaussian distribution need at least mean " + "and standard deviation values" ); TruncatedGaussian dist; - dist.mean = p.mean; - dist.standard_deviation = p.std; - dist.min_value = p.min; - dist.max_value = p.max; + dist.mean = desc.mean.value(); + dist.standard_deviation = desc.standard_deviation.value(); + dist.min_value = desc.min_value; + dist.max_value = desc.max_value; + return dist; } }, { VonMises::distribution_type_static(), - []( const DoubleSampler::DistributionDescription& d ) { - OPENGEODE_EXCEPTION( d.mean && d.standard_deviation, - "[DoubleSampler] - Provide mean and Standard deviation " - "to set up a VonMises Distribution." ); + []( const DoubleSampler::DistributionDescription& desc ) { + OPENGEODE_EXCEPTION( desc.mean && desc.kappa, + "[DoubleSampler] - Incomplete description for " + "Von Mises distribution need at least mean " + "and concentration (kappa) values" ); VonMises dist; - dist.mean = *d.mean; - dist.concentration = - compute_kappa_von_mises( *d.standard_deviation ); + dist.mean = desc.mean.value(); + dist.concentration = desc.kappa.value(); + return dist; + } }, + { TruncatedLogNormal::distribution_type_static(), + []( const DoubleSampler::DistributionDescription& desc ) { + OPENGEODE_EXCEPTION( desc.mean && desc.standard_deviation, + "[DoubleSampler] - Incomplete description for " + "TruncatedLogNormal distribution need mean " + "and standard deviation values of the underlying " + "normal distribution." ); + TruncatedLogNormal dist; + dist.mean = desc.mean.value(); + dist.standard_deviation = desc.standard_deviation.value(); + dist.min_value = desc.min_value; + dist.max_value = desc.max_value; + return dist; + } }, + { TruncatedPowerLaw::distribution_type_static(), + []( const DoubleSampler::DistributionDescription& desc ) { + OPENGEODE_EXCEPTION( desc.alpha, + "[DoubleSampler] - Incomplete description for " + "TruncatedPowerLaw distribution need power law " + "exponent (alpha)." ); + TruncatedPowerLaw dist; + dist.alpha = desc.alpha.value(); + dist.min_value = desc.min_value; + dist.max_value = desc.max_value; return dist; } }, }; @@ -208,6 +194,10 @@ namespace geode return engine.sample_truncated_gaussian( d ); if constexpr( std::is_same_v< D, VonMises > ) return engine.sample_von_mises( d ); + if constexpr( std::is_same_v< D, TruncatedLogNormal > ) + return engine.sample_truncated_lognormal( d ); + if constexpr( std::is_same_v< D, TruncatedPowerLaw > ) + return engine.sample_truncated_powerlaw( d ); throw OpenGeodeException( "DoubleSampler - Unsupported " "distribution for double" ); }, diff --git a/src/geode/stochastic/sampling/random_engine.cpp b/src/geode/stochastic/sampling/random_engine.cpp index c7bba3d..8175291 100644 --- a/src/geode/stochastic/sampling/random_engine.cpp +++ b/src/geode/stochastic/sampling/random_engine.cpp @@ -41,6 +41,80 @@ namespace { + double normal_cdf( double x ) + { + return 0.5 * std::erfc( -x / std::sqrt( 2.0 ) ); + } + + // NormalQuantile based on Peter J. Acklam's inverse normal CDF + // approximation. Original algorithm released to the public domain. + // Reference: + // https://web.archive.org/web/20150910005011/http://home.online.no/~pjacklam + double normal_quantile( double p ) + { + OPENGEODE_EXCEPTION( p >= 0.0 && p <= 1.0, + "[normal_quantile] - p must be in (0,1). Check the consistencies " + "between min,mean and max values." ); + + static const double a1 = -3.969683028665376e+01; + static const double a2 = 2.209460984245205e+02; + static const double a3 = -2.759285104469687e+02; + static const double a4 = 1.383577518672690e+02; + static const double a5 = -3.066479806614716e+01; + static const double a6 = 2.506628277459239e+00; + + static const double b1 = -5.447609879822406e+01; + static const double b2 = 1.615858368580409e+02; + static const double b3 = -1.556989798598866e+02; + static const double b4 = 6.680131188771972e+01; + static const double b5 = -1.328068155288572e+01; + + static const double c1 = -7.784894002430293e-03; + static const double c2 = -3.223964580411365e-01; + static const double c3 = -2.400758277161838e+00; + static const double c4 = -2.549732539343734e+00; + static const double c5 = 4.374664141464968e+00; + static const double c6 = 2.938163982698783e+00; + + static const double d1 = 7.784695709041462e-03; + static const double d2 = 3.224671290700398e-01; + static const double d3 = 2.445134137142996e+00; + static const double d4 = 3.754408661907416e+00; + + const double plow = 0.02425; + const double phigh = 1 - plow; + + double q, r; + if( p < plow ) + { + // lower tail + q = std::sqrt( -2 * std::log( p ) ); + return ( ( ( ( ( c1 * q + c2 ) * q + c3 ) * q + c4 ) * q + c5 ) * q + + c6 ) + / ( ( ( ( d1 * q + d2 ) * q + d3 ) * q + d4 ) * q + 1 ); + } + else if( p > phigh ) + { + // upper tail + q = std::sqrt( -2 * std::log( 1 - p ) ); + return -( ( ( ( ( c1 * q + c2 ) * q + c3 ) * q + c4 ) * q + c5 ) * q + + c6 ) + / ( ( ( ( d1 * q + d2 ) * q + d3 ) * q + d4 ) * q + 1 ); + } + else + { + // central region + q = p - 0.5; + r = q * q; + return ( ( ( ( ( a1 * r + a2 ) * r + a3 ) * r + a4 ) * r + a5 ) * r + + a6 ) + * q + / ( ( ( ( ( b1 * r + b2 ) * r + b3 ) * r + b4 ) * r + b5 ) + * r + + 1 ); + } + } + std::seed_seq create_seed_seq( uint64_t seed ) { std::vector< uint32_t > seed_data = { static_cast< uint32_t >( @@ -108,38 +182,41 @@ namespace geode OPENGEODE_ASSERT( law.standard_deviation > 0 && std::isfinite( law.standard_deviation ) && std::isfinite( law.mean ), - "[Truncated Gaussian sampling] - Infinite " - "parameters or negative standard deviation N(", - law.mean, law.standard_deviation, ")." ); - const auto max = law.max_value.value_or( + "[Gaussian sampling] - Infinite parameters or " + "negative standard deviation N(", + law.mean, "," law.standard_deviation, ")." ); + + const double max = law.max_value.value_or( std::numeric_limits< double >::infinity() ); - const auto min = law.min_value.value_or( + const double min = law.min_value.value_or( -std::numeric_limits< double >::infinity() ); OPENGEODE_ASSERT( min < max, - "[Truncated Gaussian sampling] - Wrong truncation range ", min, + "[Gaussian sampling] - Wrong truncation range ", min, " is not < than ", max, "." ); - OPENGEODE_ASSERT( - min < law.mean + 6.0 * law.standard_deviation - && max > law.mean - 6.0 * law.standard_deviation, - "[Truncated Gaussian sampling] - Truncation range is too far " - "from mean." ); - constexpr index_t MAX_ATTEMPTS = 10000; - for( const auto attempt : geode::Range( MAX_ATTEMPTS ) ) - { - geode_unused( attempt ); - auto value = absl::gaussian_distribution< double >( - law.mean, law.standard_deviation )( rand_gen_ ); - if( value >= min && value <= max ) - { - return value; - } - } - OPENGEODE_ASSERT_NOT_REACHED( - "[Truncated Gaussian sampling] - number of attemps > ", - MAX_ATTEMPTS ); - return 0.; + // Standardize bounds + const double alpha = ( min - law.mean ) / law.standard_deviation; + const double beta = ( max - law.mean ) / law.standard_deviation; + + // Compute CDF of bounds, handling infinite alpha/beta + double F_min = std::isfinite( alpha ) ? normal_cdf( alpha ) : 0.0; + double F_max = std::isfinite( beta ) ? normal_cdf( beta ) : 1.0; + + // Clamp to avoid exact 0 or 1 (normal_quantile cannot handle them) + F_min = std::max( F_min, geode::GLOBAL_EPSILON ); + F_max = std::min( F_max, 1.0 - geode::GLOBAL_EPSILON ); + + OPENGEODE_EXCEPTION( F_min < F_max, + "[Gaussian sampling] - truncation " + "range is extreme please check inputs" ); + + // Sample uniform in [F_min, F_max] + std::uniform_real_distribution< double > uniform( F_min, F_max ); + const double u = uniform( rand_gen_ ); + + // Map through inverse CDF + return law.mean + law.standard_deviation * normal_quantile( u ); } double sample_von_mises( const VonMises& law ) @@ -149,6 +226,8 @@ namespace geode "[VonMises sampling] - Invalid parameters: mean=", law.mean, ", concentration=", law.concentration, "." ); + // Uniform approximation for very small concentration (nearly + // uniform) if( law.concentration < geode::GLOBAL_EPSILON ) { UniformClosedOpen< double > uniform_dist; @@ -156,7 +235,21 @@ namespace geode return sample_uniform( uniform_dist ); } - // Best & Fisher (1979) algorithm for von Mises sampling + // Normal approximation for very large concentration + const double LARGE_KAPPA = 1e3; // threshold, can be tuned + if( law.concentration > LARGE_KAPPA ) + { + // Variance of approximate normal around mean + const double stddev = 1.0 / std::sqrt( law.concentration ); + std::normal_distribution< double > normal_dist( + law.mean, stddev ); + double theta = normal_dist( rand_gen_ ); + // Wrap to [0, 2π) + return std::fmod( theta + 2.0 * M_PI, 2.0 * M_PI ); + } + + // Best & Fisher (1979) rejection algorithm for moderate + // concentration const double a = 1.0 + std::sqrt( @@ -179,22 +272,106 @@ namespace geode { theta = std::acos( f ); if( sample_bernoulli( 0.5 ) ) - { theta = -theta; - } break; } } - // Shift by mean and normalize to [0, 2π) + // Shift by mean and wrap to [0, 2π) theta += law.mean; - theta = std::fmod( theta, 2.0 * M_PI ); - if( theta < 0.0 ) - theta += 2.0 * M_PI; + theta = std::fmod( theta + 2.0 * M_PI, 2.0 * M_PI ); return theta; } + double sample_truncated_lognormal( const TruncatedLogNormal& law ) + { + // Basic sanity checks + OPENGEODE_ASSERT( law.standard_deviation > 0 + && std::isfinite( law.standard_deviation ) + && std::isfinite( law.mean ), + "[Truncated LogNormal sampling] - Infinite parameters or " + "negative standard deviation N(", + law.mean, ", ", law.standard_deviation, ")." ); + + // Determine min and max, respecting optional values + const double min_val = law.min_value.value_or( 0.0 ); + const double max_val = law.max_value.value_or( + std::numeric_limits< double >::infinity() ); + + OPENGEODE_ASSERT( min_val < max_val, + "[Truncated LogNormal sampling] - Wrong truncation range ", + min_val, " is not < than ", max_val, "." ); + + // Transform to standard normal space + const double zmin = + ( std::log( min_val ) - law.mean ) / law.standard_deviation; + const double zmax = + ( std::isfinite( max_val ) + ? ( std::log( max_val ) - law.mean ) + / law.standard_deviation + : std::numeric_limits< double >::infinity() ); + + // Compute CDF bounds, handling infinite zmin/zmax + double F_min = + std::max( normal_cdf( zmin ), geode::GLOBAL_EPSILON ); + double F_max = + std::min( normal_cdf( zmax ), 1.0 - geode::GLOBAL_EPSILON ); + + OPENGEODE_EXCEPTION( F_min < F_max, + "[Truncated LogNormal sampling] - truncation " + "range is extreme please chack inputs" ); + + // Sample uniform in [Fmin, Fmax] + std::uniform_real_distribution< double > uniform( F_min, F_max ); + const double u = uniform( rand_gen_ ); + + // Map through inverse CDF and exponentiate + const double z = normal_quantile( u ); + return std::exp( law.mean + law.standard_deviation * z ); + } + + double sample_truncated_powerlaw( const TruncatedPowerLaw& law ) + { + OPENGEODE_ASSERT( + alpha > 0, "Power-law exponent alpha must be > 0" ); + + // Set bounds + const double xmin = law.min_value.value_or( + geode::GLOBAL_EPSILON ); // default 1.0 if unspecified + const double xmax = law.max_value.value_or( + std::numeric_limits< double >::infinity() ); + + OPENGEODE_ASSERT( xmin < xmax, "Truncated power-law: min >= max" ); + + // Sample uniform + std::uniform_real_distribution< double > uniform( 0.0, 1.0 ); + const double u = uniform( rand_gen_ ); + + // Inverse CDF + if( std::abs( law.alpha - 1.0 ) > geode::GLOBAL_EPSILON ) + { + double xmin_pow = std::pow( xmin, 1.0 - law.alpha ); + double xmax_pow = + std::isfinite( xmax ) + ? std::pow( xmax, 1.0 - law.alpha ) + : std::numeric_limits< double >::infinity(); + if( !std::isfinite( xmax_pow ) ) + { + return std::pow( xmin_pow + u, 1.0 / ( 1.0 - law.alpha ) ); + } + return std::pow( u * ( xmax_pow - xmin_pow ) + xmin_pow, + 1.0 / ( 1.0 - law.alpha ) ); + } + else + { + // alpha == 1 + const double xmax_eff = + std::isfinite( xmax ) ? xmax : xmin * geode::GLOBAL_EPSILON; + return xmin * std::pow( xmax_eff / xmin, u ); + } + } + double sample_log() { return std::log( @@ -277,6 +454,16 @@ namespace geode { return impl_->sample_von_mises( law ); } + double RandomEngine::sample_truncated_lognormal( + const TruncatedLogNormal& law ) + { + return impl_->sample_truncated_lognormal( law ); + } + double RandomEngine::sample_truncated_powerlaw( + const TruncatedPowerLaw& law ) + { + return impl_->sample_truncated_powerlaw( law ); + } double RandomEngine::sample_log() { return impl_->sample_log(); diff --git a/tests/stochastic/sampling/mcmc/test-mh-fractures.cpp b/tests/stochastic/sampling/mcmc/test-mh-fractures.cpp index 0ae37ce..1859276 100644 --- a/tests/stochastic/sampling/mcmc/test-mh-fractures.cpp +++ b/tests/stochastic/sampling/mcmc/test-mh-fractures.cpp @@ -117,15 +117,16 @@ namespace // length setB.length.distribution_type = - geode::UniformClosed< double >::distribution_type_static(); + geode::TruncatedLogNormal::distribution_type_static(); setB.length.min_value = 1; - setB.length.max_value = 10.; - + setB.length.max_value = 50.; + setB.length.mean = 1; + setB.length.standard_deviation = 1.; // azimuth setB.azimuth.distribution_type = geode::VonMises::distribution_type_static(); setB.azimuth.mean = 60.; - setB.azimuth.standard_deviation = 15.; + setB.azimuth.kappa = 1.; // positionning setB.p20 = 0.05; From 8533097937d3310ee4329b43218fbe92cb933f6b Mon Sep 17 00:00:00 2001 From: Francois Bonneau Date: Fri, 14 Nov 2025 20:22:19 +0100 Subject: [PATCH 2/4] coma --- src/geode/stochastic/sampling/random_engine.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/geode/stochastic/sampling/random_engine.cpp b/src/geode/stochastic/sampling/random_engine.cpp index 8175291..843589f 100644 --- a/src/geode/stochastic/sampling/random_engine.cpp +++ b/src/geode/stochastic/sampling/random_engine.cpp @@ -184,7 +184,7 @@ namespace geode && std::isfinite( law.mean ), "[Gaussian sampling] - Infinite parameters or " "negative standard deviation N(", - law.mean, "," law.standard_deviation, ")." ); + law.mean, ",", law.standard_deviation, ")." ); const double max = law.max_value.value_or( std::numeric_limits< double >::infinity() ); From 33ce02eb0ef00e03911b011e95dff8b2caf1e951 Mon Sep 17 00:00:00 2001 From: Francois Bonneau Date: Fri, 14 Nov 2025 20:31:15 +0100 Subject: [PATCH 3/4] compile --- src/geode/stochastic/sampling/random_engine.cpp | 2 +- tests/stochastic/sampling/mcmc/test-mh-fractures.cpp | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/geode/stochastic/sampling/random_engine.cpp b/src/geode/stochastic/sampling/random_engine.cpp index 843589f..f28f575 100644 --- a/src/geode/stochastic/sampling/random_engine.cpp +++ b/src/geode/stochastic/sampling/random_engine.cpp @@ -334,7 +334,7 @@ namespace geode double sample_truncated_powerlaw( const TruncatedPowerLaw& law ) { OPENGEODE_ASSERT( - alpha > 0, "Power-law exponent alpha must be > 0" ); + law.alpha > 0, "Power-law exponent alpha must be > 0" ); // Set bounds const double xmin = law.min_value.value_or( diff --git a/tests/stochastic/sampling/mcmc/test-mh-fractures.cpp b/tests/stochastic/sampling/mcmc/test-mh-fractures.cpp index 1859276..9d85210 100644 --- a/tests/stochastic/sampling/mcmc/test-mh-fractures.cpp +++ b/tests/stochastic/sampling/mcmc/test-mh-fractures.cpp @@ -97,7 +97,8 @@ namespace // length setA.length.distribution_type = - geode::UniformClosed< double >::distribution_type_static(); + geode::TruncatedPowerLaw::distribution_type_static(); + setA.length.alpha = 2.; setA.length.min_value = 1; setA.length.max_value = 10.; From d32c1ef8cf8788d13f6d8005c21eb49a093b0bf8 Mon Sep 17 00:00:00 2001 From: Francois Bonneau Date: Fri, 14 Nov 2025 20:46:44 +0100 Subject: [PATCH 4/4] fix --- tests/stochastic/sampling/mcmc/test-mh-fractures.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/stochastic/sampling/mcmc/test-mh-fractures.cpp b/tests/stochastic/sampling/mcmc/test-mh-fractures.cpp index 9d85210..a8f9d14 100644 --- a/tests/stochastic/sampling/mcmc/test-mh-fractures.cpp +++ b/tests/stochastic/sampling/mcmc/test-mh-fractures.cpp @@ -68,7 +68,7 @@ namespace printer_config.output_folder, "/single_fracture_set" ); geode::SimulationConfigurator sim_config; - sim_config.realizations = 500; + sim_config.realizations = 300; sim_config.metropolis_hasting_steps = 1000; sim_config.burn_in_steps = 1000; sim_config.printer = printer_config; @@ -145,7 +145,7 @@ namespace absl::StrCat( printer_config.output_folder, "/two_fracture_sets" ); geode::SimulationConfigurator sim_config; - sim_config.realizations = 500; + sim_config.realizations = 300; sim_config.metropolis_hasting_steps = 1000; sim_config.burn_in_steps = 1000; sim_config.printer = printer_config;