Merge pull request #24 from Geode-solutions/fracture_orientation

BotellaA · web-flow · commit 4e64a688eec0 · 2025-11-13T21:33:52.000+01:00
feat(SampleAzimuth): modify the segment sampler.
diff --git a/bindings/python/src/stochastic/sampling/direct/double_sampler.hpp b/bindings/python/src/stochastic/sampling/direct/double_sampler.hpp
@@ -60,6 +60,11 @@ namespace geode
             .def_static( "create_distribution",
                 &DoubleSampler::create_distribution, pybind11::arg( "desc" ),
                 "Create a distribution from a description" )
+            .def_static( "create_rad_angle_distribution_from_degree",
+                &DoubleSampler::create_rad_angle_distribution_from_degree,
+                pybind11::arg( "desc" ),
+                "Create a angle distribution in radian from a description "
+                "provided in degree" )
             .def_static( "sample", &DoubleSampler::sample,
                 pybind11::arg( "engine" ), pybind11::arg( "dist" ),
                 "Sample a value from a distribution using a RandomEngine" );
diff --git a/bindings/python/src/stochastic/sampling/distributions.hpp b/bindings/python/src/stochastic/sampling/distributions.hpp
@@ -87,6 +87,17 @@ namespace geode
                 &TruncatedGaussian::distribution_type_static )
             .def( "distribution_type", &TruncatedGaussian::distribution_type )
             .def( "string", &TruncatedGaussian::string );
+
+        // VonMises
+        pybind11::class_< VonMises >( module, "VonMises" )
+            .def( pybind11::init<>() )
+            .def_readwrite( "mean", &VonMises::mean )
+            .def_readwrite( "concentration", &VonMises::concentration )
+            .def( "is_valid", &VonMises::is_valid )
+            .def_static( "distribution_type_static",
+                &VonMises::distribution_type_static )
+            .def( "distribution_type", &VonMises::distribution_type )
+            .def( "string", &VonMises::string );
     }
 
 } // namespace geode
diff --git a/bindings/python/src/stochastic/sampling/random_engine.hpp b/bindings/python/src/stochastic/sampling/random_engine.hpp
@@ -64,6 +64,9 @@ namespace geode
                 &RandomEngine::sample_truncated_gaussian,
                 pybind11::arg( "law" ),
                 "Sample a value from a truncated Gaussian" )
+            .def( "sample_von_mises", &RandomEngine::sample_von_mises,
+                pybind11::arg( "law" ),
+                "Sample a value from a Von Mises-Fisher" )
 
             // Other distributions
             .def( "sample_log", &RandomEngine::sample_log,
diff --git a/bindings/python/tests/stochastic/test-py-mh-fractures.py b/bindings/python/tests/stochastic/test-py-mh-fractures.py
@@ -57,7 +57,7 @@ def test_fracture_simulator():
     setA.azimuth.max_value = 10.0
 
     # positioning
-    setA.p20 = 0.1
+    setA.p20 = 0.06
     setA.minimal_spacing = 1.0
 
     runner = stochastic.FractureSimulationRunner(box)
@@ -69,7 +69,7 @@ def test_fracture_simulator():
     printer_config.output_folder = os.path.join(printer_config.output_folder , "py_single_fracture_set")
 
     sim_config = stochastic.SimulationConfigurator()
-    sim_config.realizations = 1000
+    sim_config.realizations = 500
     sim_config.metropolis_hasting_steps = 1000
     sim_config.burn_in_steps = 1000
     sim_config.printer = printer_config
@@ -97,10 +97,10 @@ def test_two_fracture_sets_simulator():
     setA.length.distribution_type = stochastic.DistributionType("UniformClosed")
     setA.length.min_value = 1.0
     setA.length.max_value = 10.0
-    setA.azimuth.distribution_type = stochastic.DistributionType("UniformClosed")
-    setA.azimuth.min_value = 1.0
-    setA.azimuth.max_value = 10.0
-    setA.p20 = 0.1
+    setA.azimuth.distribution_type = stochastic.DistributionType("VonMises")
+    setA.azimuth.mean = 45
+    setA.azimuth.standard_deviation = 10.0
+    setA.p20 = 0.05
     setA.minimal_spacing = 1.0
 
     # --- Object set B
@@ -112,7 +112,7 @@ def test_two_fracture_sets_simulator():
     setB.azimuth.distribution_type =stochastic.DistributionType("UniformClosed")
     setB.azimuth.min_value = 90.0
     setB.azimuth.max_value = 100.0
-    setB.p20 = 0.1
+    setB.p20 = 0.03
     setB.minimal_spacing = 2.0
 
     runner = stochastic.FractureSimulationRunner(box)
@@ -124,7 +124,7 @@ def test_two_fracture_sets_simulator():
     printer_config.output_folder = os.path.join(printer_config.output_folder ,printer_config.output_folder, "py_two_fracture_sets")
 
     sim_config = stochastic.SimulationConfigurator()
-    sim_config.realizations = 1000
+    sim_config.realizations = 500
     sim_config.metropolis_hasting_steps = 1000
     sim_config.burn_in_steps = 1000
     sim_config.printer = printer_config
diff --git a/include/geode/stochastic/sampling/direct/double_sampler.hpp b/include/geode/stochastic/sampling/direct/double_sampler.hpp
@@ -41,7 +41,8 @@ namespace geode
         using Distribution = std::variant< UniformClosed< double >,
             UniformClosedOpen< double >,
             Gaussian,
-            TruncatedGaussian >;
+            TruncatedGaussian,
+            VonMises >;
 
         struct DistributionDescription
         {
@@ -84,6 +85,8 @@ namespace geode
 
         static Distribution create_distribution(
             const DistributionDescription& desc );
+        static Distribution create_rad_angle_distribution_from_degree(
+            const DistributionDescription& angle_desc_deg );
 
         static double sample( RandomEngine& engine, const Distribution& dist );
     };
diff --git a/include/geode/stochastic/sampling/direct/object_set_sampler/object_set_sampler.hpp b/include/geode/stochastic/sampling/direct/object_set_sampler/object_set_sampler.hpp
@@ -41,6 +41,7 @@ namespace geode
         {
             if( !is_valid_object( obj ) )
             {
+                Logger::warn( "[ObjectSetSampler] - invalid object proposed." );
                 return LOG_PROB_INVALID;
             }
             return log_pdf_;
diff --git a/include/geode/stochastic/sampling/direct/object_set_sampler/segment_set_sampler.hpp b/include/geode/stochastic/sampling/direct/object_set_sampler/segment_set_sampler.hpp
@@ -63,6 +63,7 @@ namespace geode
             const auto& extremities = obj.vertices();
             const auto current =
                 static_cast< local_index_t >( engine.sample_bernoulli( 0.5 ) );
+            const auto other = 1 - current;
 
             geode::Sphere< 2 > ball{ extremities[current],
                 ratio * obj.length() };
@@ -71,7 +72,8 @@ namespace geode
             constexpr index_t max_try{ 100 };
             for( const auto try_id : geode::Range{ max_try } )
             {
-                if( box_.contains( new_point ) )
+                if( box_.contains( new_point )
+                    || box_.contains( extremities[other] ) )
                 {
                     OwnerSegment2D new_segment{ obj };
                     new_segment.set_point( current, new_point );
@@ -91,7 +93,7 @@ namespace geode
             const auto& extremities = obj.vertices();
 
             return box_.contains( extremities[0] )
-                   && box_.contains( extremities[1] );
+                   || box_.contains( extremities[1] );
         }
 
     private:
diff --git a/include/geode/stochastic/sampling/distributions.hpp b/include/geode/stochastic/sampling/distributions.hpp
@@ -176,4 +176,34 @@ namespace geode
         }
     };
 
+    struct opengeode_stochastic_stochastic_api VonMises
+    {
+        VonMises() = default;
+
+        // 2D von Mises parameters
+        double mean{ 0.0 }; // mean direction in radians
+        double concentration{ 1.0 }; // kappa (concentration parameter)
+
+        [[nodiscard]] static DistributionType distribution_type_static()
+        {
+            return DistributionType{ "VonMises" };
+        }
+
+        [[nodiscard]] DistributionType distribution_type() const
+        {
+            return distribution_type_static();
+        }
+
+        [[nodiscard]] bool is_valid() const
+        {
+            return concentration >= 0.0 && std::isfinite( mean );
+        }
+
+        std::string string() const
+        {
+            return absl::StrCat( distribution_type().get(), "(mean=", mean,
+                ", kappa=", concentration, ")" );
+        }
+    };
+
 } // namespace geode
diff --git a/include/geode/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp b/include/geode/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp
@@ -97,7 +97,8 @@ namespace geode
                 auto length_distribution =
                     DoubleSampler::create_distribution( set_desc.length );
                 auto azimuth_distribution =
-                    DoubleSampler::create_distribution( set_desc.azimuth );
+                    DoubleSampler::create_rad_angle_distribution_from_degree(
+                        set_desc.azimuth );
                 this->set_samplers_.push_back(
                     std::make_unique< UniformSegmentSetSampler >(
                         box_, length_distribution, azimuth_distribution ) );
diff --git a/include/geode/stochastic/sampling/random_engine.hpp b/include/geode/stochastic/sampling/random_engine.hpp
@@ -50,6 +50,7 @@ namespace geode
 
         double sample_gaussian( const Gaussian& law );
         double sample_truncated_gaussian( const TruncatedGaussian& law );
+        double sample_von_mises( const VonMises& law );
 
         double sample_log();
 
diff --git a/src/geode/stochastic/sampling/direct/double_sampler.cpp b/src/geode/stochastic/sampling/direct/double_sampler.cpp
@@ -23,8 +23,11 @@
 
 #include <geode/stochastic/sampling/direct/double_sampler.hpp>
 
+#include <geode/geometry/angle.hpp>
+
 #include <geode/stochastic/common.hpp>
 #include <geode/stochastic/sampling/random_engine.hpp>
+
 namespace
 {
     struct DistributionParams
@@ -63,6 +66,20 @@ namespace
             "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
 {
@@ -117,6 +134,17 @@ namespace geode
                     dist.max_value = p.max;
                     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." );
+                    VonMises dist;
+                    dist.mean = *d.mean;
+                    dist.concentration =
+                        compute_kappa_von_mises( *d.standard_deviation );
+                    return dist;
+                } },
         };
 
     DoubleSampler::Distribution DoubleSampler::create_distribution(
@@ -129,6 +157,41 @@ namespace geode
         return it->second( desc );
     }
 
+    DoubleSampler::Distribution
+        DoubleSampler::create_rad_angle_distribution_from_degree(
+            const DistributionDescription& desc_deg )
+    {
+        DistributionDescription desc_rad = desc_deg;
+        if( desc_rad.mean )
+        {
+            auto mean_angle = Angle::create_from_degrees( *( desc_rad.mean ) );
+            desc_rad.mean = mean_angle.normalized_between_0_and_2pi().radians();
+        }
+        if( desc_rad.standard_deviation )
+        {
+            auto std_angle =
+                Angle::create_from_degrees( *( desc_rad.standard_deviation ) );
+            desc_rad.standard_deviation =
+                std_angle.normalized_between_0_and_2pi().radians();
+        }
+        if( desc_rad.min_value )
+        {
+            auto min_angle =
+                Angle::create_from_degrees( *( desc_rad.min_value ) );
+            desc_rad.min_value =
+                min_angle.normalized_between_0_and_2pi().radians();
+        }
+        if( desc_rad.max_value )
+        {
+            auto max_angle =
+                Angle::create_from_degrees( *( desc_rad.max_value ) );
+            desc_rad.max_value =
+                max_angle.normalized_between_0_and_2pi().radians();
+        }
+        // Call the general create_distribution
+        return create_distribution( desc_rad );
+    }
+
     double DoubleSampler::sample(
         RandomEngine& engine, const Distribution& dist )
     {
@@ -143,6 +206,8 @@ namespace geode
                     return engine.sample_gaussian( d );
                 if constexpr( std::is_same_v< D, TruncatedGaussian > )
                     return engine.sample_truncated_gaussian( d );
+                if constexpr( std::is_same_v< D, VonMises > )
+                    return engine.sample_von_mises( d );
                 throw OpenGeodeException( "DoubleSampler - Unsupported "
                                           "distribution for double" );
             },
diff --git a/src/geode/stochastic/sampling/direct/segment_uniform_sampler.cpp b/src/geode/stochastic/sampling/direct/segment_uniform_sampler.cpp
@@ -26,20 +26,34 @@
 #include <geode/stochastic/common.hpp>
 #include <geode/stochastic/sampling/random_engine.hpp>
 
+#include <geode/geometry/angle.hpp>
 #include <geode/geometry/basic_objects/segment.hpp>
+#include <geode/geometry/vector.hpp>
 
+namespace
+{
+    geode::Vector2D direction_from_azimuth_angle( const geode::Angle& azimuth )
+    {
+        return geode::Vector2D( { azimuth.sin(), azimuth.cos() } );
+    }
+} // namespace
 namespace geode
 {
 
     OwnerSegment2D SegmentUniformSampler::sample( RandomEngine& engine,
         const PointUniformSampler::Object< 2 >& object,
         const DoubleSampler::Distribution& length,
-        const DoubleSampler::Distribution& azimuth )
+        const DoubleSampler::Distribution& azimuth_rad )
     {
-        geode_unused( length );
-        geode_unused( azimuth );
         auto point1 = PointUniformSampler::sample( engine, object );
-        auto point2 = PointUniformSampler::sample( engine, object );
+
+        auto segment_length = DoubleSampler::sample( engine, length );
+        auto segment_azimuth = Angle::create_from_radians(
+            DoubleSampler::sample( engine, azimuth_rad ) );
+
+        auto point2 =
+            point1
+            + direction_from_azimuth_angle( segment_azimuth ) * segment_length;
 
         return OwnerSegment2D{ point1, point2 };
     }
diff --git a/src/geode/stochastic/sampling/random_engine.cpp b/src/geode/stochastic/sampling/random_engine.cpp
diff --git a/tests/stochastic/sampling/direct/test-segment-uniform-sampler.cpp b/tests/stochastic/sampling/direct/test-segment-uniform-sampler.cpp
diff --git a/tests/stochastic/sampling/mcmc/test-mh-fractures.cpp b/tests/stochastic/sampling/mcmc/test-mh-fractures.cpp

Original file line number	Diff line number	Diff line change
`@@ -41,6 +41,7 @@ namespace geode`
`41`	`41`	`{`
`42`	`42`	`if( !is_valid_object( obj ) )`
`43`	`43`	`{`
	`44`	`+ Logger::warn( "[ObjectSetSampler] - invalid object proposed." );`
`44`	`45`	`return LOG_PROB_INVALID;`
`45`	`46`	`}`
`46`	`47`	`return log_pdf_;`