From a68ffce199aada54cb8f7e0bcccc2c952be7625d Mon Sep 17 00:00:00 2001 From: Francois Bonneau Date: Mon, 22 Dec 2025 12:49:15 +0100 Subject: [PATCH 1/4] feat(FractureModel): add p21 controle --- .../helpers/fracture_simulation_runner.hpp | 13 +- .../mcmc/models/components/density_term.hpp | 160 +++++++++------ .../mcmc/models/components/energy_term.hpp | 18 +- .../mcmc/models/components/pairwise_term.hpp | 77 +++++--- .../models/components/single_object_term.hpp | 61 +++--- src/geode/stochastic/CMakeLists.txt | 1 + tests/stochastic/CMakeLists.txt | 8 + .../models/components/test-intensity-term.cpp | 182 ++++++++++++++++++ .../sampling/mcmc/test-mh-fractures.cpp | 1 + 9 files changed, 393 insertions(+), 128 deletions(-) create mode 100644 tests/stochastic/sampling/mcmc/models/components/test-intensity-term.cpp 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 f885a50..3d47683 100644 --- a/include/geode/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp +++ b/include/geode/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp @@ -48,7 +48,7 @@ namespace geode // positionning double p20; std::vector< std::array< geode::Point2D, 2 > > observed_fractures; - // double p21; + double p21{ 1. }; double minimal_spacing{ 0. }; // mh dynamique @@ -147,6 +147,17 @@ namespace geode absl::StrCat( set_desc.name, "_density" ), set_desc.p20, std::vector< uuid >{ set_id }, this->domain_ ) ) ); + // intentisty + if( std::fabs( set_desc.p21 - 1. ) > GLOBAL_EPSILON ) + { + this->ordered_energy_terms_.push_back( + this->energy_terms_collection_.add_energy_term( + std::make_unique< IntensityTerm >( + absl::StrCat( set_desc.name, "_intensity" ), + set_desc.p21, std::vector< uuid >{ set_id }, + 0.5 * this->domain_.smallest_length(), + this->domain_ ) ) ); + } // spacing if( set_desc.minimal_spacing > GLOBAL_EPSILON ) { diff --git a/include/geode/stochastic/sampling/mcmc/models/components/density_term.hpp b/include/geode/stochastic/sampling/mcmc/models/components/density_term.hpp index 7b39ee7..53c0a86 100644 --- a/include/geode/stochastic/sampling/mcmc/models/components/density_term.hpp +++ b/include/geode/stochastic/sampling/mcmc/models/components/density_term.hpp @@ -24,79 +24,121 @@ #include -#include +#include + +double length_inside_box( const geode::Point2D& p0, + const geode::Point2D& p1, + const geode::BoundingBox2D& box ) +{ + double dx = p1.value( 0 ) - p0.value( 0 ); + double dy = p1.value( 1 ) - p0.value( 1 ); + + double t_min = 0.0; + double t_max = 1.0; + + auto update_interval = [&]( double p, double q_min, double q_max ) -> bool { + if( p == 0.0 ) + { + // Segment parallel to axis + return ( q_min <= 0.0 && 0.0 <= q_max ); + } + + double t1 = q_min / p; + double t2 = q_max / p; + if( t1 > t2 ) + std::swap( t1, t2 ); + + t_min = std::max( t_min, t1 ); + t_max = std::min( t_max, t2 ); + + return t_min <= t_max; + }; + + // X axis + if( !update_interval( dx, box.min().value( 0 ) - p0.value( 0 ), + box.max().value( 0 ) - p0.value( 0 ) ) ) + { + return 0.0; + } + + // Y axis + if( !update_interval( dy, box.min().value( 1 ) - p0.value( 1 ), + box.max().value( 1 ) - p0.value( 1 ) ) ) + { + return 0.0; + } + + if( t_max <= t_min ) + { + return 0.0; + } + + double clipped_dx = dx * ( t_max - t_min ); + double clipped_dy = dy * ( t_max - t_min ); + + return std::sqrt( clipped_dx * clipped_dx + clipped_dy * clipped_dy ); +} namespace geode { - /// ObjectCountTerm template < typename ObjectType > - class DensityTerm : public EnergyTerm< ObjectType > + class DensityTerm : public SingleObjectTerm< ObjectType, + std::function< double( const ObjectType&, + const SpatialDomain< ObjectType::dim >& ) > > { public: explicit DensityTerm( std::string_view name, double lambda, std::vector< uuid > targeted_set_ids, const SpatialDomain< ObjectType::dim >& domain ) - : EnergyTerm< ObjectType >( - name, lambda, std::move( targeted_set_ids ), domain ) + : SingleObjectTerm< ObjectType, + std::function< double( const ObjectType&, + const SpatialDomain< ObjectType::dim >& ) > >( + name, + lambda, + std::move( targeted_set_ids ), + 1.0, // scale by domain area to get density per unit + []( const ObjectType& obj, + const SpatialDomain< ObjectType::dim >& domain ) { + if( SpatialDomainChecker< + ObjectType >::is_anchored_in_domain( domain, + obj ) ) + { + return 1.0; + } + return 0.0; + }, // contribution = 1 anchoredin domain + domain ) { } + }; - double total_log( const ObjectSets< ObjectType >& state ) const override - { - const auto n = this->statistic( state ); - return this->contribution( n ); - } - - double delta_log_add( const ObjectSets< ObjectType >& /*state*/, - const ObjectRef< ObjectType >& new_object ) const override - { - if( !this->is_targeted_set( new_object.set_id ) - || !this->is_anchored_in_domain( new_object.object ) ) - { - return 0.0; - } - return this->contribution( 1.0 ); - } - - double delta_log_remove( const ObjectSets< ObjectType >& state, - const ObjectId& object_id ) const override - { - if( !this->is_targeted_set( object_id.set_id ) - || !this->is_anchored_in_domain( - state.get_object( object_id ) ) ) - { - return 0.0; - } - return this->contribution( -1.0 ); - } - - double delta_log_change( const ObjectSets< ObjectType >& state, - const ObjectId& old_object_id, - const ObjectRef< ObjectType >& new_object ) const override - { - auto old_in = this->is_anchored_in_domain( - state.get_object( old_object_id ) ); - auto new_in = this->is_anchored_in_domain( new_object.object ); - if( old_in == new_in ) - { - return 0.0; - } - return this->contribution( new_in ? 1.0 : -1.0 ); - } - - double statistic( const ObjectSets< ObjectType >& state ) const override + class IntensityTerm + : public SingleObjectTerm< OwnerSegment2D, + std::function< double( const OwnerSegment2D&, + const SpatialDomain< OwnerSegment2D::dim >& ) > > + { + public: + explicit IntensityTerm( std::string_view name, + double lambda, + std::vector< uuid > targeted_set_ids, + double caracteristic_length, + const SpatialDomain< OwnerSegment2D::dim >& domain ) + : SingleObjectTerm< OwnerSegment2D, + std::function< double( const OwnerSegment2D&, + const SpatialDomain< OwnerSegment2D::dim >& ) > >( + name, + lambda, + std::move( targeted_set_ids ), + 1.0 / caracteristic_length, + []( const OwnerSegment2D& segment, + const SpatialDomain< OwnerSegment2D::dim >& domain ) { + auto seg_extremities = segment.vertices(); + return length_inside_box( seg_extremities[0], + seg_extremities[1], domain.box() ); + }, + domain ) { - double sum = 0.0; - this->for_each_targeted_object( - state, [&]( const ObjectId& obj_id ) { - if( this->is_anchored_in_domain( - state.get_object( obj_id ) ) ) - { - sum += 1.0; - } - } ); - return sum; } }; } // namespace geode \ No newline at end of file diff --git a/include/geode/stochastic/sampling/mcmc/models/components/energy_term.hpp b/include/geode/stochastic/sampling/mcmc/models/components/energy_term.hpp index d07be4a..416ad9c 100644 --- a/include/geode/stochastic/sampling/mcmc/models/components/energy_term.hpp +++ b/include/geode/stochastic/sampling/mcmc/models/components/energy_term.hpp @@ -156,17 +156,17 @@ namespace geode targeted_set_ids_.begin(), targeted_set_ids_.end(), set_id ); } - bool is_anchored_in_domain( const ObjectType& obj ) const + const SpatialDomain< ObjectType::dim >& domain() const { - return SpatialDomainChecker< ObjectType >::is_anchored_in_domain( - domain_, obj ); + return domain_; } - bool intersects_domain( const ObjectType& obj ) const - { - return SpatialDomainChecker< ObjectType >::intersects_domain( - domain_, obj ); - } + // bool intersects_domain( const ObjectType& obj ) const + // { + // return SpatialDomainChecker< ObjectType + // >::intersects_domain( + // domain_, obj ); + // } template < typename Func > void for_each_targeted_object( @@ -185,6 +185,6 @@ namespace geode private: detail::EnergyScale energy_scale_; std::vector< uuid > targeted_set_ids_; - const SpatialDomain< ObjectType::dim > domain_; + const SpatialDomain< ObjectType::dim >& domain_; }; } // namespace geode \ No newline at end of file diff --git a/include/geode/stochastic/sampling/mcmc/models/components/pairwise_term.hpp b/include/geode/stochastic/sampling/mcmc/models/components/pairwise_term.hpp index f84162d..082c30f 100644 --- a/include/geode/stochastic/sampling/mcmc/models/components/pairwise_term.hpp +++ b/include/geode/stochastic/sampling/mcmc/models/components/pairwise_term.hpp @@ -68,8 +68,13 @@ namespace geode geode::ObjectRef< ObjectType > neigh_object{ state.get_object( neigh_id ), neigh_id.set_id }; - if( this->is_anchored_in_domain( new_object.object ) - || this->is_anchored_in_domain( neigh_object.object ) ) + // is that reallythe good test? + // intersect? + if( SpatialDomainChecker< ObjectType >::is_anchored_in_domain( + this->domain(), new_object.object ) + || SpatialDomainChecker< + ObjectType >::is_anchored_in_domain( this->domain(), + neigh_object.object ) ) { delta += interaction_->evaluate( new_object, neigh_object ); } @@ -96,8 +101,11 @@ namespace geode ObjectRef< ObjectType > neigh_object{ state.get_object( neigh_id ), neigh_id.set_id }; - if( this->is_anchored_in_domain( object_to_remove.object ) - || this->is_anchored_in_domain( neigh_object.object ) ) + if( SpatialDomainChecker< ObjectType >::is_anchored_in_domain( + this->domain(), object_to_remove.object ) + || SpatialDomainChecker< + ObjectType >::is_anchored_in_domain( this->domain(), + neigh_object.object ) ) { delta += interaction_->evaluate( object_to_remove, neigh_object ); @@ -129,8 +137,11 @@ namespace geode ObjectRef< ObjectType > neigh_object{ state.get_object( neigh_id ), neigh_id.set_id }; - if( this->is_anchored_in_domain( object_to_remove.object ) - || this->is_anchored_in_domain( neigh_object.object ) ) + if( SpatialDomainChecker< ObjectType >::is_anchored_in_domain( + this->domain(), object_to_remove.object ) + || SpatialDomainChecker< + ObjectType >::is_anchored_in_domain( this->domain(), + neigh_object.object ) ) { delta -= interaction_->evaluate( object_to_remove, neigh_object ); @@ -150,8 +161,11 @@ namespace geode ObjectRef< ObjectType > neigh_object{ state.get_object( neigh_id ), neigh_id.set_id }; - if( this->is_anchored_in_domain( new_object.object ) - || this->is_anchored_in_domain( neigh_object.object ) ) + if( SpatialDomainChecker< ObjectType >::is_anchored_in_domain( + this->domain(), new_object.object ) + || SpatialDomainChecker< + ObjectType >::is_anchored_in_domain( this->domain(), + neigh_object.object ) ) { delta += interaction_->evaluate( new_object, neigh_object ); } @@ -162,31 +176,31 @@ namespace geode double statistic( const ObjectSets< ObjectType >& state ) const override { double sum = 0.0; - this->for_each_targeted_object( - state, [&]( const ObjectId& obj_id ) { - const auto& cur_obj = state.get_object( obj_id ); - if( !this->is_anchored_in_domain( cur_obj ) ) + this->for_each_targeted_object( state, [&]( const ObjectId& + obj_id ) { + const auto& cur_obj = state.get_object( obj_id ); + if( !SpatialDomainChecker< ObjectType >::is_anchored_in_domain( + this->domain(), cur_obj ) ) + { + return; + } + ObjectRef< ObjectType > object{ cur_obj, obj_id.set_id }; + const auto neighbors = + state.neighbors( obj_id, this->targeted_set_ids(), + interaction_->neighborhood_searching_distance() ); + for( const auto& neigh_obj_id : neighbors ) + { + if( !is_new_pair( cur_obj, obj_id, neigh_obj_id ) ) { - return; + continue; } - ObjectRef< ObjectType > object{ cur_obj, obj_id.set_id }; - const auto neighbors = - state.neighbors( obj_id, this->targeted_set_ids(), - interaction_->neighborhood_searching_distance() ); - for( const auto& neigh_obj_id : neighbors ) - { - if( !is_new_pair( cur_obj, obj_id, neigh_obj_id ) ) - { - continue; - } - ObjectRef< ObjectType > neigh_object{ - state.get_object( neigh_obj_id ), - neigh_obj_id.set_id - }; + ObjectRef< ObjectType > neigh_object{ + state.get_object( neigh_obj_id ), neigh_obj_id.set_id + }; - sum += interaction_->evaluate( object, neigh_object ); - } - } ); + sum += interaction_->evaluate( object, neigh_object ); + } + } ); return sum; } @@ -195,7 +209,8 @@ namespace geode const ObjectId& obj_id, const ObjectId& neigh_obj_id ) const { - if( !this->is_anchored_in_domain( obj ) ) + if( !SpatialDomainChecker< ObjectType >::is_anchored_in_domain( + this->domain(), obj ) ) { return true; } diff --git a/include/geode/stochastic/sampling/mcmc/models/components/single_object_term.hpp b/include/geode/stochastic/sampling/mcmc/models/components/single_object_term.hpp index ac415c0..e3b0eea 100644 --- a/include/geode/stochastic/sampling/mcmc/models/components/single_object_term.hpp +++ b/include/geode/stochastic/sampling/mcmc/models/components/single_object_term.hpp @@ -25,6 +25,7 @@ #include #include +#include namespace geode { @@ -34,67 +35,71 @@ namespace geode public: explicit SingleObjectTerm( std::string_view name, double lambda, - absl::flat_hash_set< uuid > targeted_set_ids, - ObjectContributionFunc contribution_func ) - : EnergyTerm< ObjectType >( name, lambda, targeted_set_ids ), + std::vector< uuid > targeted_set_ids, + double scale, + ObjectContributionFunc contribution_func, + const SpatialDomain< ObjectType::dim >& domain ) + : EnergyTerm< ObjectType >( + name, lambda, std::move( targeted_set_ids ), domain ), + scale_( scale ), contribution_func_( std::move( contribution_func ) ) { } double total_log( const ObjectSets< ObjectType >& state ) const override { - auto total = this->statistic( state ); - return this->contribution( total ); + return this->contribution( scale_ * statistic( state ) ); } double delta_log_add( const ObjectSets< ObjectType >& /*state*/, - const ObjectType& new_object, - const uuid& new_object_set_id ) const override + const ObjectRef< ObjectType >& new_object ) const override { - if( !this->is_targeted_set( new_object_set_id ) ) + if( !this->is_targeted_set( new_object.set_id ) ) { return 0.0; } - return this->contribution( contribution_func_( new_object ) ); + return this->contribution( + scale_ + * contribution_func_( new_object.object, this->domain() ) ); } double delta_log_remove( const ObjectSets< ObjectType >& state, - ObjectId object_id ) const override + const ObjectId& object_id ) const override { if( !this->is_targeted_set( object_id.set_id ) ) { return 0.0; } return this->contribution( - -contribution_func_( state.get_object( object_id ) ) ); + -scale_ + * contribution_func_( + state.get_object( object_id ), this->domain() ) ); } double delta_log_change( const ObjectSets< ObjectType >& state, - ObjectId old_object_id, - const ObjectType& new_object, - const uuid& new_object_set_id ) const override + const ObjectId& old_object_id, + const ObjectRef< ObjectType >& new_object ) const override { - if( this->is_targeted_set( old_object_id.set_id ) - && this->is_targeted_set( new_object_set_id ) ) - { - return 0.0; - } - auto delta = - contribution_func_( new_object ) - - contribution_func_( state.get_object( old_object_id ) ); - return this->contribution( delta ); + double delta = + contribution_func_( new_object.object, this->domain() ) + - contribution_func_( + state.get_object( old_object_id ), this->domain() ); + return this->contribution( scale_ * delta ); } double statistic( const ObjectSets< ObjectType >& state ) const override { - double total = 0.0; - this->for_each_targeted_object( state, [&]( const ObjectId& id ) { - total += contribution_func_( state.get_object( id ) ); - } ); - return total; + double sum = 0.0; + this->for_each_targeted_object( + state, [&]( const ObjectId& obj_id ) { + const auto& obj = state.get_object( obj_id ); + sum += contribution_func_( obj, this->domain() ); + } ); + return sum; } private: + double scale_{ 1. }; ObjectContributionFunc contribution_func_; }; } // namespace geode \ No newline at end of file diff --git a/src/geode/stochastic/CMakeLists.txt b/src/geode/stochastic/CMakeLists.txt index a1baa80..3749d7a 100644 --- a/src/geode/stochastic/CMakeLists.txt +++ b/src/geode/stochastic/CMakeLists.txt @@ -56,6 +56,7 @@ add_geode_library( "sampling/mcmc/helpers/simulation_monitor.hpp" "sampling/mcmc/models/components/pairwise_term.hpp" "sampling/mcmc/models/components/density_term.hpp" + "sampling/mcmc/models/components/single_object_term.hpp" "sampling/mcmc/models/components/energy_term.hpp" "sampling/mcmc/models/energy_term_collection.hpp" "sampling/mcmc/models/gibbs_energy.hpp" diff --git a/tests/stochastic/CMakeLists.txt b/tests/stochastic/CMakeLists.txt index d003af1..7e9a61a 100644 --- a/tests/stochastic/CMakeLists.txt +++ b/tests/stochastic/CMakeLists.txt @@ -113,6 +113,14 @@ add_geode_test( ${PROJECT_NAME}::stochastic ) +add_geode_test( + SOURCE "sampling/mcmc/models/components/test-intensity-term.cpp" + DEPENDENCIES + OpenGeode::basic + OpenGeode::geometry + ${PROJECT_NAME}::stochastic +) + add_geode_test( SOURCE "sampling/mcmc/models/components/test-pairwise-term.cpp" DEPENDENCIES diff --git a/tests/stochastic/sampling/mcmc/models/components/test-intensity-term.cpp b/tests/stochastic/sampling/mcmc/models/components/test-intensity-term.cpp new file mode 100644 index 0000000..4db911c --- /dev/null +++ b/tests/stochastic/sampling/mcmc/models/components/test-intensity-term.cpp @@ -0,0 +1,182 @@ +/* + * Copyright (c) 2019 - 2025 Geode-solutions + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + * + */ +#include + +#include +#include +#include + +geode::uuid init_segment_set( + geode::ObjectSets< geode::OwnerSegment2D >& pattern ) +{ + using Segment = geode::OwnerSegment2D; + + // Fully inside domain + Segment s1{ geode::Point2D{ { -0.5, 0.0 } }, + geode::Point2D{ { 0.5, 0.0 } } }; // length = 1.0 + + // Partially inside domain + Segment s2{ geode::Point2D{ { -2.0, 0.0 } }, + geode::Point2D{ { 0.0, 0.0 } } }; // clipped length = 1.0 + + // Fully outside domain (buffer) + Segment s_buffer{ geode::Point2D{ { 2.0, 2.0 } }, + geode::Point2D{ { 3.0, 3.0 } } }; // clipped = 0.0 + + auto set_id = pattern.add_set( "segments" ); + pattern.add_object( std::move( s1 ), set_id, false ); + pattern.add_object( std::move( s2 ), set_id, false ); + pattern.add_object( std::move( s_buffer ), set_id, false ); + + return set_id; +} + +geode::SpatialDomain< 2 > init_domain() +{ + geode::BoundingBox2D box; + box.add_point( geode::Point2D{ { -1.0, -1.0 } } ); + box.add_point( geode::Point2D{ { 1.0, 1.0 } } ); + return geode::SpatialDomain< 2 >{ box, 0.5 }; +} + +void run_intensity_test( double lambda, + double characteristic_length, + const geode::ObjectSets< geode::OwnerSegment2D >& pattern, + const geode::uuid& set_id, + const geode::SpatialDomain< 2 >& domain ) +{ + geode::IntensityTerm term( + "intensity", lambda, { set_id }, characteristic_length, domain ); + + const double neg_log_lambda = + ( lambda > 0. ? -std::log( lambda ) + : std::numeric_limits< double >::infinity() ); + + // Total clipped length inside domain: + // s1: 1.0 + // s2: 1.0 + // s_buffer: 0.0 + const double total_length = 2.0; + const double scaled_total = total_length / characteristic_length; + + const double expected_total = + ( lambda > 0. ? neg_log_lambda * scaled_total + : std::numeric_limits< double >::infinity() ); + + // --- Total log + double total = term.total_log( pattern ); + OPENGEODE_EXCEPTION( + total == expected_total, "[IntensityTerm] total_log wrong" ); + + // --- Delta add (segment fully inside) + geode::OwnerSegment2D s_inside{ geode::Point2D{ { 0.0, -0.5 } }, + geode::Point2D{ { 0.0, 0.5 } } }; // length 1 + geode::ObjectRef< geode::OwnerSegment2D > ref_inside{ s_inside, set_id }; + + double expected_add = + ( lambda > 0. ? neg_log_lambda * ( 1.0 / characteristic_length ) + : std::numeric_limits< double >::infinity() ); + + double delta = term.delta_log_add( pattern, ref_inside ); + OPENGEODE_EXCEPTION( + delta == expected_add, "[IntensityTerm] delta_log_add inside wrong" ); + + // --- Delta add outside domain + geode::OwnerSegment2D s_outside{ geode::Point2D{ { 2.0, 2.0 } }, + geode::Point2D{ { 3.0, 3.0 } } }; + geode::ObjectRef< geode::OwnerSegment2D > ref_out{ s_outside, set_id }; + + delta = term.delta_log_add( pattern, ref_out ); + OPENGEODE_EXCEPTION( + delta == 0.0, "[IntensityTerm] delta_log_add outside wrong" ); + + // --- Delta remove (first segment) + geode::ObjectId obj_id{ 0, false, set_id }; + double expected_remove = ( lambda > 0. ? -expected_add : 0.0 ); + + delta = term.delta_log_remove( pattern, obj_id ); + OPENGEODE_EXCEPTION( + delta == expected_remove, "[IntensityTerm] delta_log_remove wrong" ); + + // --- Delta change: inside → outside + delta = term.delta_log_change( pattern, obj_id, ref_out ); + OPENGEODE_EXCEPTION( delta == expected_remove, + "[IntensityTerm] delta_log_change inside→outside wrong" ); + + // --- Delta change: outside → inside + geode::ObjectId buffer_id{ 2, false, set_id }; + delta = term.delta_log_change( pattern, buffer_id, ref_inside ); + OPENGEODE_EXCEPTION( delta == expected_add, + "[IntensityTerm] delta_log_change outside→inside wrong" ); + + // --- Delta change: inside → inside (same length) + geode::OwnerSegment2D s_same{ geode::Point2D{ { -0.2, 0.0 } }, + geode::Point2D{ { 0.8, 0.0 } } }; // length still 1 + geode::ObjectRef< geode::OwnerSegment2D > ref_same{ s_same, set_id }; + + delta = term.delta_log_change( pattern, obj_id, ref_same ); + OPENGEODE_EXCEPTION( + delta == 0.0, "[IntensityTerm] delta_log_change inside→inside wrong" ); +} + +int main() +{ + try + { + geode::StochasticLibrary::initialize(); + geode::Logger::set_level( geode::Logger::LEVEL::debug ); + + geode::ObjectSets< geode::OwnerSegment2D > pattern; + auto set_id = init_segment_set( pattern ); + auto domain = init_domain(); + + const double L0 = 1.0; + + run_intensity_test( 0.5, L0, pattern, set_id, domain ); + run_intensity_test( + geode::GLOBAL_EPSILON, L0, pattern, set_id, domain ); + run_intensity_test( 50.0, L0, pattern, set_id, domain ); + run_intensity_test( 0.0, L0, pattern, set_id, domain ); + } + catch( ... ) + { + return geode::geode_lippincott(); + } + + try + { + geode::StochasticLibrary::initialize(); + geode::uuid set_id; + auto domain = init_domain(); + + geode::DensityTerm< geode::Point2D > term( + "zero", -geode::GLOBAL_EPSILON, { set_id }, domain ); + geode::Logger::info( "TEST FAILED" ); + return 1; + } + catch( geode::OpenGeodeException& expt ) + { + geode::Logger::info( "TEST SUCCESS" ); + return 0; + } +} \ No newline at end of file diff --git a/tests/stochastic/sampling/mcmc/test-mh-fractures.cpp b/tests/stochastic/sampling/mcmc/test-mh-fractures.cpp index 97fc99c..50f79c9 100644 --- a/tests/stochastic/sampling/mcmc/test-mh-fractures.cpp +++ b/tests/stochastic/sampling/mcmc/test-mh-fractures.cpp @@ -59,6 +59,7 @@ namespace // positionning setA.p20 = 0.05; + setA.p21 = 200; setA.minimal_spacing = 1.; geode::FractureSimulationRunner runner( domain ); From f4279185c5496bbcb5d2ecc4ef4d29a230fce222 Mon Sep 17 00:00:00 2001 From: Francois Bonneau Date: Mon, 22 Dec 2025 12:53:06 +0100 Subject: [PATCH 2/4] add binding --- .../sampling/mcmc/helpers/fracture_simulation_runner.hpp | 1 + bindings/python/tests/stochastic/test-py-mh-fractures.py | 1 + 2 files changed, 2 insertions(+) diff --git a/bindings/python/src/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp b/bindings/python/src/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp index 760369a..84cdbc2 100644 --- a/bindings/python/src/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp +++ b/bindings/python/src/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp @@ -36,6 +36,7 @@ namespace geode .def_readwrite( "length", &FractureSetDescription::length ) .def_readwrite( "azimuth", &FractureSetDescription::azimuth ) .def_readwrite( "p20", &FractureSetDescription::p20 ) + .def_readwrite( "p20", &FractureSetDescription::p21 ) .def( "add_observed_fracture", []( FractureSetDescription& self, const geode::Point2D& a, const geode::Point2D& b ) { diff --git a/bindings/python/tests/stochastic/test-py-mh-fractures.py b/bindings/python/tests/stochastic/test-py-mh-fractures.py index 7450563..4ee9a80 100644 --- a/bindings/python/tests/stochastic/test-py-mh-fractures.py +++ b/bindings/python/tests/stochastic/test-py-mh-fractures.py @@ -60,6 +60,7 @@ def test_fracture_simulator(): # positioning setA.p20 = 0.06 + setA.p21 = 10 setA.minimal_spacing = 1.0 runner = stochastic.FractureSimulationRunner(domain) From 9d53c3b01e178fe4f18b16051fc8eac05e5c94b8 Mon Sep 17 00:00:00 2001 From: Francois Bonneau Date: Mon, 22 Dec 2025 13:01:04 +0100 Subject: [PATCH 3/4] add binding --- .../sampling/mcmc/helpers/fracture_simulation_runner.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/bindings/python/src/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp b/bindings/python/src/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp index 84cdbc2..34abc84 100644 --- a/bindings/python/src/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp +++ b/bindings/python/src/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp @@ -36,7 +36,7 @@ namespace geode .def_readwrite( "length", &FractureSetDescription::length ) .def_readwrite( "azimuth", &FractureSetDescription::azimuth ) .def_readwrite( "p20", &FractureSetDescription::p20 ) - .def_readwrite( "p20", &FractureSetDescription::p21 ) + .def_readwrite( "p21", &FractureSetDescription::p21 ) .def( "add_observed_fracture", []( FractureSetDescription& self, const geode::Point2D& a, const geode::Point2D& b ) { From ab52b5744edff21c8583bb5aec012cfe91fe9440 Mon Sep 17 00:00:00 2001 From: Francois Bonneau Date: Mon, 5 Jan 2026 10:41:49 +0100 Subject: [PATCH 4/4] review arnaud --- .../helpers/fracture_simulation_runner.hpp | 137 ++++++++++-------- .../mcmc/models/components/density_term.hpp | 82 ----------- .../mcmc/models/components/intensity_term.hpp | 125 ++++++++++++++++ src/geode/stochastic/CMakeLists.txt | 5 +- .../models/components/test-intensity-term.cpp | 19 +-- 5 files changed, 209 insertions(+), 159 deletions(-) create mode 100644 include/geode/stochastic/sampling/mcmc/models/components/intensity_term.hpp 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 3d47683..dc67046 100644 --- a/include/geode/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp +++ b/include/geode/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp @@ -29,6 +29,7 @@ #include #include #include +#include #include #include #include @@ -107,7 +108,8 @@ namespace geode std::make_unique< ProposalKernel< OwnerSegment2D > >(); // Mapping set names -> UUID - std::unordered_map< std::string, uuid > name_to_uuid; + // std::unordered_map< std::string, uuid > + // set_name_to_uuid_; // Step 1: create object sets and samplers for( const auto& set_desc : set_descriptors_ ) @@ -120,7 +122,7 @@ namespace geode fixed_object[0], fixed_object[1] }, set_id, true ); } - name_to_uuid[set_desc.name] = set_id; + set_name_to_uuid_[set_desc.name] = set_id; auto length_distribution = DoubleSampler::create_distribution( set_desc.length ); @@ -139,63 +141,12 @@ namespace geode // Step 2: create density energy terms for( const auto& set_desc : set_descriptors_ ) { - const auto set_id = name_to_uuid.at( set_desc.name ); - // p20 - this->ordered_energy_terms_.push_back( - this->energy_terms_collection_.add_energy_term( - std::make_unique< DensityTerm< OwnerSegment2D > >( - absl::StrCat( set_desc.name, "_density" ), - set_desc.p20, std::vector< uuid >{ set_id }, - this->domain_ ) ) ); - // intentisty - if( std::fabs( set_desc.p21 - 1. ) > GLOBAL_EPSILON ) - { - this->ordered_energy_terms_.push_back( - this->energy_terms_collection_.add_energy_term( - std::make_unique< IntensityTerm >( - absl::StrCat( set_desc.name, "_intensity" ), - set_desc.p21, std::vector< uuid >{ set_id }, - 0.5 * this->domain_.smallest_length(), - this->domain_ ) ) ); - } - // spacing - if( set_desc.minimal_spacing > GLOBAL_EPSILON ) - { - auto interaction = std::make_unique< - EuclideanCutoffInteraction< OwnerSegment2D > >( - set_desc.minimal_spacing, - PairwiseInteraction< - OwnerSegment2D >::SCOPE::same_set ); - - this->ordered_energy_terms_.push_back( - this->energy_terms_collection_.add_energy_term( - std::make_unique< PairwiseTerm< OwnerSegment2D > >( - absl::StrCat( set_desc.name, "_min_spacing" ), - 0., std::vector< uuid >{ set_id }, - std::move( interaction ), this->domain_ ) ) ); - } + set_density_term( set_desc ); + set_intensity_term( set_desc ); + set_minimal_spacing_term( set_desc ); } - // x node monitoring - if( std::fabs( beta_x_node_ - 1. ) > GLOBAL_EPSILON ) - { - std::vector< uuid > set_uuids; - set_uuids.reserve( name_to_uuid.size() ); - for( const auto& [name, id] : name_to_uuid ) - { - set_uuids.push_back( id ); - } - auto interaction = std::make_unique< - EuclideanCutoffInteraction< OwnerSegment2D > >( - 0., PairwiseInteraction< - OwnerSegment2D >::SCOPE::different_set ); + set_x_intersection_term(); - this->ordered_energy_terms_.push_back( - this->energy_terms_collection_.add_energy_term( - std::make_unique< PairwiseTerm< OwnerSegment2D > >( - absl::StrCat( "inter_set_x_nodes" ), beta_x_node_, - set_uuids, std::move( interaction ), - this->domain_ ) ) ); - } this->mh_sampler_ = std::make_unique< MetropolisHastings< OwnerSegment2D > >( this->energy_terms_collection_, @@ -240,8 +191,80 @@ namespace geode return message; } + private: + void set_density_term( const FractureSetDescription& set_desc ) + { + const auto set_id = set_name_to_uuid_.at( set_desc.name ); + this->ordered_energy_terms_.push_back( + this->energy_terms_collection_.add_energy_term( + std::make_unique< DensityTerm< OwnerSegment2D > >( + absl::StrCat( set_desc.name, "_density" ), set_desc.p20, + std::vector< uuid >{ set_id }, this->domain_ ) ) ); + } + + void set_intensity_term( const FractureSetDescription& set_desc ) + { + if( std::fabs( set_desc.p21 - 1. ) < GLOBAL_EPSILON ) + { + return; + } + const auto set_id = set_name_to_uuid_.at( set_desc.name ); + this->ordered_energy_terms_.push_back( + this->energy_terms_collection_.add_energy_term( + std::make_unique< IntensityTerm >( + absl::StrCat( set_desc.name, "_intensity" ), + set_desc.p21, std::vector< uuid >{ set_id }, + 0.5 * this->domain_.smallest_length(), + this->domain_ ) ) ); + } + + void set_minimal_spacing_term( const FractureSetDescription& set_desc ) + { + if( set_desc.minimal_spacing < GLOBAL_EPSILON ) + { + return; + } + const auto set_id = set_name_to_uuid_.at( set_desc.name ); + auto interaction = std::make_unique< + EuclideanCutoffInteraction< OwnerSegment2D > >( + set_desc.minimal_spacing, + PairwiseInteraction< OwnerSegment2D >::SCOPE::same_set ); + + this->ordered_energy_terms_.push_back( + this->energy_terms_collection_.add_energy_term( + std::make_unique< PairwiseTerm< OwnerSegment2D > >( + absl::StrCat( set_desc.name, "_min_spacing" ), 0., + std::vector< uuid >{ set_id }, std::move( interaction ), + this->domain_ ) ) ); + } + + void set_x_intersection_term() + { + if( std::fabs( beta_x_node_ - 1. ) > GLOBAL_EPSILON ) + { + std::vector< uuid > set_uuids; + set_uuids.reserve( set_name_to_uuid_.size() ); + for( const auto& [name, id] : set_name_to_uuid_ ) + { + set_uuids.push_back( id ); + } + auto interaction = std::make_unique< + EuclideanCutoffInteraction< OwnerSegment2D > >( + 0., PairwiseInteraction< + OwnerSegment2D >::SCOPE::different_set ); + + this->ordered_energy_terms_.push_back( + this->energy_terms_collection_.add_energy_term( + std::make_unique< PairwiseTerm< OwnerSegment2D > >( + absl::StrCat( "inter_set_x_nodes" ), beta_x_node_, + set_uuids, std::move( interaction ), + this->domain_ ) ) ); + } + } + private: std::vector< FractureSetDescription > set_descriptors_; + std::unordered_map< std::string, uuid > set_name_to_uuid_; // x node monitoring double beta_x_node_{ 1. }; }; diff --git a/include/geode/stochastic/sampling/mcmc/models/components/density_term.hpp b/include/geode/stochastic/sampling/mcmc/models/components/density_term.hpp index 53c0a86..d316097 100644 --- a/include/geode/stochastic/sampling/mcmc/models/components/density_term.hpp +++ b/include/geode/stochastic/sampling/mcmc/models/components/density_term.hpp @@ -26,59 +26,6 @@ #include -double length_inside_box( const geode::Point2D& p0, - const geode::Point2D& p1, - const geode::BoundingBox2D& box ) -{ - double dx = p1.value( 0 ) - p0.value( 0 ); - double dy = p1.value( 1 ) - p0.value( 1 ); - - double t_min = 0.0; - double t_max = 1.0; - - auto update_interval = [&]( double p, double q_min, double q_max ) -> bool { - if( p == 0.0 ) - { - // Segment parallel to axis - return ( q_min <= 0.0 && 0.0 <= q_max ); - } - - double t1 = q_min / p; - double t2 = q_max / p; - if( t1 > t2 ) - std::swap( t1, t2 ); - - t_min = std::max( t_min, t1 ); - t_max = std::min( t_max, t2 ); - - return t_min <= t_max; - }; - - // X axis - if( !update_interval( dx, box.min().value( 0 ) - p0.value( 0 ), - box.max().value( 0 ) - p0.value( 0 ) ) ) - { - return 0.0; - } - - // Y axis - if( !update_interval( dy, box.min().value( 1 ) - p0.value( 1 ), - box.max().value( 1 ) - p0.value( 1 ) ) ) - { - return 0.0; - } - - if( t_max <= t_min ) - { - return 0.0; - } - - double clipped_dx = dx * ( t_max - t_min ); - double clipped_dy = dy * ( t_max - t_min ); - - return std::sqrt( clipped_dx * clipped_dx + clipped_dy * clipped_dy ); -} - namespace geode { template < typename ObjectType > @@ -112,33 +59,4 @@ namespace geode { } }; - - class IntensityTerm - : public SingleObjectTerm< OwnerSegment2D, - std::function< double( const OwnerSegment2D&, - const SpatialDomain< OwnerSegment2D::dim >& ) > > - { - public: - explicit IntensityTerm( std::string_view name, - double lambda, - std::vector< uuid > targeted_set_ids, - double caracteristic_length, - const SpatialDomain< OwnerSegment2D::dim >& domain ) - : SingleObjectTerm< OwnerSegment2D, - std::function< double( const OwnerSegment2D&, - const SpatialDomain< OwnerSegment2D::dim >& ) > >( - name, - lambda, - std::move( targeted_set_ids ), - 1.0 / caracteristic_length, - []( const OwnerSegment2D& segment, - const SpatialDomain< OwnerSegment2D::dim >& domain ) { - auto seg_extremities = segment.vertices(); - return length_inside_box( seg_extremities[0], - seg_extremities[1], domain.box() ); - }, - domain ) - { - } - }; } // namespace geode \ No newline at end of file diff --git a/include/geode/stochastic/sampling/mcmc/models/components/intensity_term.hpp b/include/geode/stochastic/sampling/mcmc/models/components/intensity_term.hpp new file mode 100644 index 0000000..c31b546 --- /dev/null +++ b/include/geode/stochastic/sampling/mcmc/models/components/intensity_term.hpp @@ -0,0 +1,125 @@ +/* + * Copyright (c) 2019 - 2025 Geode-solutions + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + * + */ +#pragma once + +#include + +#include + +double length_inside_box( const geode::Point2D& segment_start, + const geode::Point2D& segment_end, + const geode::BoundingBox2D& box ) +{ + // Segment direction + const double dir_x = segment_end.value( 0 ) - segment_start.value( 0 ); + const double dir_y = segment_end.value( 1 ) - segment_start.value( 1 ); + + // Parameter interval where the segment is inside the box + double t_enter = 0.0; // start of valid interval + double t_exit = 1.0; // end of valid interval + + // Clips the parametric segment against one axis interval + auto clip_against_axis_interval = [&t_enter, &t_exit]( double direction, + double min_distance, + double max_distance ) -> bool { + // Segment is parallel to this axis + if( std::fabs( direction ) < geode::GLOBAL_EPSILON ) + { + // Outside the interval → no intersection + return min_distance <= 0.0 && 0.0 <= max_distance; + } + + double t0 = min_distance / direction; + double t1 = max_distance / direction; + + if( t0 > t1 ) + { + std::swap( t0, t1 ); + } + + t_enter = std::max( t_enter, t0 ); + t_exit = std::min( t_exit, t1 ); + + return t_enter <= t_exit; + }; + + // Clip against X interval of the box + if( !clip_against_axis_interval( dir_x, + box.min().value( 0 ) - segment_start.value( 0 ), + box.max().value( 0 ) - segment_start.value( 0 ) ) ) + { + return 0.0; + } + + // Clip against Y interval of the box + if( !clip_against_axis_interval( dir_y, + box.min().value( 1 ) - segment_start.value( 1 ), + box.max().value( 1 ) - segment_start.value( 1 ) ) ) + { + return 0.0; + } + + // No portion of the segment is inside the box + if( t_exit <= t_enter ) + { + return 0.0; + } + + // Length of the clipped segment + const double clipped_dx = dir_x * ( t_exit - t_enter ); + const double clipped_dy = dir_y * ( t_exit - t_enter ); + + return std::sqrt( clipped_dx * clipped_dx + clipped_dy * clipped_dy ); +} + +namespace geode +{ + class IntensityTerm + : public SingleObjectTerm< OwnerSegment2D, + std::function< double( const OwnerSegment2D&, + const SpatialDomain< OwnerSegment2D::dim >& ) > > + { + public: + explicit IntensityTerm( std::string_view name, + double lambda, + std::vector< uuid > targeted_set_ids, + double caracteristic_length, + const SpatialDomain< OwnerSegment2D::dim >& domain ) + : SingleObjectTerm< OwnerSegment2D, + std::function< double( const OwnerSegment2D&, + const SpatialDomain< OwnerSegment2D::dim >& ) > >( + name, + lambda, + std::move( targeted_set_ids ), + 1.0 / caracteristic_length, + []( const OwnerSegment2D& segment, + const SpatialDomain< OwnerSegment2D::dim >& domain ) { + auto seg_extremities = segment.vertices(); + return length_inside_box( seg_extremities[0], + seg_extremities[1], domain.box() ); + }, + domain ) + { + } + }; +} // namespace geode \ No newline at end of file diff --git a/src/geode/stochastic/CMakeLists.txt b/src/geode/stochastic/CMakeLists.txt index 3749d7a..f7bc3af 100644 --- a/src/geode/stochastic/CMakeLists.txt +++ b/src/geode/stochastic/CMakeLists.txt @@ -54,10 +54,11 @@ add_geode_library( "sampling/mcmc/helpers/simulation_runner.hpp" "sampling/mcmc/helpers/simulation_printer.hpp" "sampling/mcmc/helpers/simulation_monitor.hpp" - "sampling/mcmc/models/components/pairwise_term.hpp" "sampling/mcmc/models/components/density_term.hpp" - "sampling/mcmc/models/components/single_object_term.hpp" "sampling/mcmc/models/components/energy_term.hpp" + "sampling/mcmc/models/components/intensity_term.hpp" + "sampling/mcmc/models/components/pairwise_term.hpp" + "sampling/mcmc/models/components/single_object_term.hpp" "sampling/mcmc/models/energy_term_collection.hpp" "sampling/mcmc/models/gibbs_energy.hpp" "sampling/mcmc/proposal/classical_proposals.hpp" diff --git a/tests/stochastic/sampling/mcmc/models/components/test-intensity-term.cpp b/tests/stochastic/sampling/mcmc/models/components/test-intensity-term.cpp index 4db911c..3a8eab3 100644 --- a/tests/stochastic/sampling/mcmc/models/components/test-intensity-term.cpp +++ b/tests/stochastic/sampling/mcmc/models/components/test-intensity-term.cpp @@ -20,7 +20,7 @@ * SOFTWARE. * */ -#include +#include #include #include @@ -162,21 +162,4 @@ int main() { return geode::geode_lippincott(); } - - try - { - geode::StochasticLibrary::initialize(); - geode::uuid set_id; - auto domain = init_domain(); - - geode::DensityTerm< geode::Point2D > term( - "zero", -geode::GLOBAL_EPSILON, { set_id }, domain ); - geode::Logger::info( "TEST FAILED" ); - return 1; - } - catch( geode::OpenGeodeException& expt ) - { - geode::Logger::info( "TEST SUCCESS" ); - return 0; - } } \ No newline at end of file