feat(EdgeEffect): manage Spatial Domain in stochastic simulation

Author: francoisbonneau

bindings/python/src/stochastic/CMakeLists.txt

@@ -28,6 +28,7 @@ add_geode_python_binding(
         "sampling/direct/double_sampler.hpp"
         "sampling/random_engine.hpp"
         "sampling/distributions.hpp"
+        "spatial/spatial_domain.hpp"
         "stochastic.cpp"
     DEPENDENCIES
         ${PROJECT_NAME}::stochastic

bindings/python/src/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp

@@ -53,7 +53,7 @@ namespace geode
         pybind11::class_< FractureSimulationRunner,
             std::shared_ptr< FractureSimulationRunner > >(
             module, "FractureSimulationRunner" )
-            .def( pybind11::init< const BoundingBox2D& >(),
+            .def( pybind11::init< const SpatialDomain< 2 >& >(),
                 pybind11::arg( "box" ) )
             .def( "add_fracture_set_descriptor",
                 &FractureSimulationRunner::add_fracture_set_descriptor,

bindings/python/src/stochastic/spatial/spatial_domain.hpp

@@ -0,0 +1,60 @@
+/*
+ * 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 <geode/stochastic/spatial/spatial_domain.hpp>
+
+namespace
+{
+    template < geode::index_t dimension >
+    void declare_spatial_domain( pybind11::module& module )
+    {
+        using Domain = geode::SpatialDomain< dimension >;
+        using BBox = geode::BoundingBox< dimension >;
+
+        const auto pyclass_name =
+            "SpatialDomain" + std::to_string( dimension ) + "D";
+
+        pybind11::class_< Domain >( module, pyclass_name.c_str() )
+            .def( pybind11::init< const BBox&, double >(),
+                pybind11::arg( "domain" ), pybind11::arg( "buffer_size" ),
+                R"doc(
+                Create a spatial domain composed of:
+
+                - a core domain (the VOI)
+                - a buffer zone
+                - the extended domain (domain + buffer)
+
+                Arguments:
+                    domain (BoundingBox): main domain / VOI
+                    buffer_size (float): buffer thickness
+            )doc" );
+    }
+} // namespace
+namespace geode
+{
+    void define_spatial_domain( pybind11::module& module )
+    {
+        declare_spatial_domain< 2 >( module );
+        declare_spatial_domain< 3 >( module );
+    }
+} // namespace geode

bindings/python/src/stochastic/stochastic.cpp

@@ -33,13 +33,16 @@
 
 #include "sampling/distributions.hpp"
 #include "sampling/random_engine.hpp"
+#include "spatial/spatial_domain.hpp"
 
 PYBIND11_MODULE( opengeode_stochastic_py_stochastic, module )
 {
     module.doc() = "OpenGeode-Stochastic Python binding";
     pybind11::class_< geode::StochasticLibrary >( module, "StochasticLibrary" )
         .def( "initialize", &geode::StochasticLibrary::initialize );
 
+    geode::define_spatial_domain( module );
+
     geode::define_distributions( module );
     geode::define_random_engine( module );
     geode::define_double_sampler( module );

bindings/python/tests/stochastic/test-py-mh-fractures.py

@@ -41,6 +41,7 @@ def test_fracture_simulator():
     box = og.BoundingBox2D()
     box.add_point(og.Point2D([0.0, 0.0]))
     box.add_point(og.Point2D([100.0, 100.0]))
+    domain = stochastic.SpatialDomain2D(box,10.)
 
     # --- Object set
     setA = stochastic.FractureSetDescription()
@@ -60,7 +61,7 @@ def test_fracture_simulator():
     setA.p20 = 0.06
     setA.minimal_spacing = 1.0
 
-    runner = stochastic.FractureSimulationRunner(box)
+    runner = stochastic.FractureSimulationRunner(domain)
     runner.add_fracture_set_descriptor(setA)
     runner.initialize()
 
@@ -90,7 +91,7 @@ def test_two_fracture_sets_simulator():
     box = og.BoundingBox2D()
     box.add_point(og.Point2D([0.0, 0.0]))
     box.add_point(og.Point2D([100.0, 100.0]))
-
+    domain = stochastic.SpatialDomain2D(box,10.)
     # --- Object set A
     setA = stochastic.FractureSetDescription()
     setA.name = "A"
@@ -117,7 +118,7 @@ def test_two_fracture_sets_simulator():
     setB.p20 = 0.03
     setB.minimal_spacing = 2.0
 
-    runner = stochastic.FractureSimulationRunner(box)
+    runner = stochastic.FractureSimulationRunner(domain)
     runner.add_fracture_set_descriptor(setA)
     runner.add_fracture_set_descriptor(setB)
     runner.initialize()

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

@@ -121,7 +121,8 @@ namespace geode
                     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 } ) ) );
+                            set_desc.p20, std::vector< uuid >{ set_id },
+                            this->domain_ ) ) );
                 // spacing
                 if( set_desc.minimal_spacing < GLOBAL_EPSILON )
                 {
@@ -137,7 +138,7 @@ namespace geode
                         std::make_unique< PairwiseTerm< OwnerSegment2D > >(
                             absl::StrCat( set_desc.name, "_min_spacing" ), 0.,
                             std::vector< uuid >{ set_id },
-                            std::move( interaction ) ) ) );
+                            std::move( interaction ), this->domain_ ) ) );
             }
 
             this->mh_sampler_ =

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

@@ -35,9 +35,10 @@ namespace geode
     public:
         explicit DensityTerm( std::string_view name,
             double lambda,
-            std::vector< uuid > targeted_set_ids )
+            std::vector< uuid > targeted_set_ids,
+            const SpatialDomain< ObjectType::dim >& domain )
             : EnergyTerm< ObjectType >(
-                  name, lambda, std::move( targeted_set_ids ) )
+                  name, lambda, std::move( targeted_set_ids ), domain )
         {
         }
 
@@ -50,38 +51,52 @@ namespace geode
         double delta_log_add( const ObjectSets< ObjectType >& /*state*/,
             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 )
+                || !this->is_anchored_in_domain( new_object.object ) )
             {
                 return 0.0;
             }
             return this->contribution( 1.0 );
         }
 
-        double delta_log_remove( const ObjectSets< ObjectType >& /*state*/,
+        double delta_log_remove( const ObjectSets< ObjectType >& state,
             const ObjectId& object_id ) const override
         {
-            if( !this->is_targeted_set( object_id.set_id ) )
+            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
+        double delta_log_change( const ObjectSets< ObjectType >& state,
+            const ObjectId& old_object_id,
+            const ObjectRef< ObjectType >& new_object ) const override
         {
-            return 0.0;
+            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
         {
-            index_t count{ 0 };
-            for( const auto& set_id : this->targeted_set_ids() )
-            {
-                count += state.nb_objects_in_set( set_id );
-            }
-            return static_cast< double >( count );
+            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

@@ -30,6 +30,8 @@
 #include <geode/stochastic/common.hpp>
 #include <geode/stochastic/spatial/object_sets.hpp>
 
+#include <geode/stochastic/spatial/spatial_domain.hpp>
+
 #include <optional>
 
 namespace geode
@@ -88,9 +90,11 @@ namespace geode
     public:
         explicit EnergyTerm( std::string_view name,
             double param,
-            std::vector< uuid >&& targeted_set_ids )
+            std::vector< uuid >&& targeted_set_ids,
+            const SpatialDomain< ObjectType::dim >& domain )
             : energy_scale_{ param },
-              targeted_set_ids_{ std::move( targeted_set_ids ) }
+              targeted_set_ids_{ std::move( targeted_set_ids ) },
+              domain_( domain )
         {
             std::sort( targeted_set_ids_.begin(), targeted_set_ids_.end() );
             IdentifierBuilder builder( *this );
@@ -152,6 +156,18 @@ namespace geode
                 targeted_set_ids_.begin(), targeted_set_ids_.end(), set_id );
         }
 
+        bool is_anchored_in_domain( const ObjectType& obj ) const
+        {
+            return SpatialDomainChecker< ObjectType >::is_anchored_in_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(
             const ObjectSets< ObjectType >& state, Func&& do_apply ) const
@@ -169,5 +185,6 @@ namespace geode
     private:
         detail::EnergyScale energy_scale_;
         std::vector< uuid > targeted_set_ids_;
+        const SpatialDomain< ObjectType::dim > domain_;
     };
 } // namespace geode