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..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,6 +36,7 @@ namespace geode .def_readwrite( "length", &FractureSetDescription::length ) .def_readwrite( "azimuth", &FractureSetDescription::azimuth ) .def_readwrite( "p20", &FractureSetDescription::p20 ) + .def_readwrite( "p21", &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) 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..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 @@ -48,7 +49,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 @@ -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,52 +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_ ) ) ); - // 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_, @@ -229,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 7b39ee7..d316097 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,39 @@ #include -#include +#include 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 - { - 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/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/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..f7bc3af 100644 --- a/src/geode/stochastic/CMakeLists.txt +++ b/src/geode/stochastic/CMakeLists.txt @@ -54,9 +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/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/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..3a8eab3 --- /dev/null +++ b/tests/stochastic/sampling/mcmc/models/components/test-intensity-term.cpp @@ -0,0 +1,165 @@ +/* + * 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(); + } +} \ 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 );