diff --git a/bindings/python/src/stochastic/sampling/direct/double_sampler.hpp b/bindings/python/src/stochastic/sampling/direct/double_sampler.hpp index f23e7a0..8cddc09 100644 --- 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 index a749bc4..54f979a 100644 --- 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 index 8034f1f..546decf 100644 --- 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 index 7f14fa7..8f979b4 100644 --- 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 index 277f9e2..5541fab 100644 --- 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 index 8ab882a..d499bc6 100644 --- 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 index 115b306..45fedbf 100644 --- 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 index 1c4a600..c2e3b56 100644 --- 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 \ No newline at end of file diff --git a/include/geode/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp b/include/geode/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp index 20c5848..29cb619 100644 --- 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 index b7eb7b8..f0825eb 100644 --- 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 index 5bc6604..fc8ac6c 100644 --- a/src/geode/stochastic/sampling/direct/double_sampler.cpp +++ b/src/geode/stochastic/sampling/direct/double_sampler.cpp @@ -23,8 +23,11 @@ #include +#include + #include #include + 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 index c34059d..158a173 100644 --- 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 #include +#include #include +#include +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 index 1582ad0..c7bba3d 100644 --- a/src/geode/stochastic/sampling/random_engine.cpp +++ b/src/geode/stochastic/sampling/random_engine.cpp @@ -142,6 +142,59 @@ namespace geode return 0.; } + double sample_von_mises( const VonMises& law ) + { + OPENGEODE_ASSERT( + law.concentration >= 0.0 && std::isfinite( law.mean ), + "[VonMises sampling] - Invalid parameters: mean=", law.mean, + ", concentration=", law.concentration, "." ); + + if( law.concentration < geode::GLOBAL_EPSILON ) + { + UniformClosedOpen< double > uniform_dist; + uniform_dist.max_value = 2.0 * M_PI; + return sample_uniform( uniform_dist ); + } + + // Best & Fisher (1979) algorithm for von Mises sampling + const double a = + 1.0 + + std::sqrt( + 1.0 + 4.0 * law.concentration * law.concentration ); + const double b = + ( a - std::sqrt( 2.0 * a ) ) / ( 2.0 * law.concentration ); + const double r = ( 1.0 + b * b ) / ( 2.0 * b ); + + double theta; + UniformClosed< double > uniform_dist; + while( true ) + { + double u1 = sample_uniform( uniform_dist ); + double z = std::cos( M_PI * u1 ); + double f = ( 1.0 + r * std::abs( z ) ) / ( r + std::abs( z ) ); + double c = law.concentration * ( r - f ); + double u2 = sample_uniform( uniform_dist ); + + if( u2 < c * ( 2.0 - c ) || u2 <= c * std::exp( 1.0 - c ) ) + { + theta = std::acos( f ); + if( sample_bernoulli( 0.5 ) ) + { + theta = -theta; + } + break; + } + } + + // Shift by mean and normalize to [0, 2π) + theta += law.mean; + theta = std::fmod( theta, 2.0 * M_PI ); + if( theta < 0.0 ) + theta += 2.0 * M_PI; + + return theta; + } + double sample_log() { return std::log( @@ -220,6 +273,10 @@ namespace geode { return impl_->sample_truncated_gaussian( law ); } + double RandomEngine::sample_von_mises( const VonMises& law ) + { + return impl_->sample_von_mises( law ); + } double RandomEngine::sample_log() { return impl_->sample_log(); diff --git a/tests/stochastic/sampling/direct/test-segment-uniform-sampler.cpp b/tests/stochastic/sampling/direct/test-segment-uniform-sampler.cpp index 348ee9c..19d74d9 100644 --- a/tests/stochastic/sampling/direct/test-segment-uniform-sampler.cpp +++ b/tests/stochastic/sampling/direct/test-segment-uniform-sampler.cpp @@ -43,14 +43,16 @@ void test_sample_segment( { geode::UniformClosed< double > length; geode::UniformClosed< double > az; + auto box_enlarged = box; + box_enlarged.extends( 1.01 ); for( const auto i : geode::Range{ NUMBER_OF_SAMPLES } ) { auto value = geode::SegmentUniformSampler::sample( engine, box, length, az ); - OPENGEODE_EXCEPTION( box.contains( value.vertices()[0] ), + OPENGEODE_EXCEPTION( box_enlarged.contains( value.vertices()[0] ), "[SegmentUniformSampler] - segment out of box." ); - OPENGEODE_EXCEPTION( box.contains( value.vertices()[1] ), + OPENGEODE_EXCEPTION( box_enlarged.contains( value.vertices()[1] ), "[SegmentUniformSampler] - segment out of box." ); } } diff --git a/tests/stochastic/sampling/mcmc/test-mh-fractures.cpp b/tests/stochastic/sampling/mcmc/test-mh-fractures.cpp index ae70366..0ae37ce 100644 --- a/tests/stochastic/sampling/mcmc/test-mh-fractures.cpp +++ b/tests/stochastic/sampling/mcmc/test-mh-fractures.cpp @@ -54,7 +54,7 @@ namespace setA.azimuth.max_value = 10.; // positionning - setA.p20 = 0.1; + setA.p20 = 0.05; setA.minimal_spacing = 1.; geode::FractureSimulationRunner runner( box ); @@ -68,7 +68,7 @@ namespace printer_config.output_folder, "/single_fracture_set" ); geode::SimulationConfigurator sim_config; - 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; @@ -108,7 +108,7 @@ namespace setA.azimuth.max_value = 10.; // positionning - setA.p20 = 0.1; + setA.p20 = 0.05; setA.minimal_spacing = 1.; // --- Object set @@ -123,12 +123,12 @@ namespace // azimuth setB.azimuth.distribution_type = - geode::UniformClosed< double >::distribution_type_static(); - setB.azimuth.min_value = 90.; - setB.azimuth.max_value = 100.; + geode::VonMises::distribution_type_static(); + setB.azimuth.mean = 60.; + setB.azimuth.standard_deviation = 15.; // positionning - setB.p20 = 0.1; + setB.p20 = 0.05; setB.minimal_spacing = 2.; geode::FractureSimulationRunner runner( box ); @@ -143,7 +143,7 @@ namespace absl::StrCat( printer_config.output_folder, "/two_fracture_sets" ); geode::SimulationConfigurator sim_config; - 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;