Merge pull request #30 from Geode-solutions/contrain_p21

feat(FractureModel): add p21 controle

Author: BotellaA

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 ) {

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)

include/geode/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp

@@ -29,6 +29,7 @@
 #include <geode/stochastic/sampling/direct/object_set_sampler/segment_set_sampler.hpp>
 #include <geode/stochastic/sampling/mcmc/helpers/simulation_runner.hpp>
 #include <geode/stochastic/sampling/mcmc/models/components/density_term.hpp>
+#include <geode/stochastic/sampling/mcmc/models/components/intensity_term.hpp>
 #include <geode/stochastic/sampling/mcmc/models/components/pairwise_term.hpp>
 #include <geode/stochastic/sampling/mcmc/proposal/classical_proposals.hpp>
 #include <geode/stochastic/spatial/pairwise_interactions.hpp>
@@ -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. };
     };

include/geode/stochastic/sampling/mcmc/models/components/density_term.hpp

@@ -24,79 +24,39 @@
 
 #include <geode/stochastic/spatial/object_sets.hpp>
 
-#include <geode/stochastic/sampling/mcmc/models/components/energy_term.hpp>
+#include <geode/stochastic/sampling/mcmc/models/components/single_object_term.hpp>
 
 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

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