Merge pull request #27 from Geode-solutions/spatial_window

Spatial_window

Author: BotellaA

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/direct/object_set_sampler/point_set_sampler.hpp

@@ -28,60 +28,72 @@
 
 #include <geode/stochastic/sampling/direct/object_set_sampler/object_set_sampler.hpp>
 #include <geode/stochastic/sampling/direct/point_uniform_sampler.hpp>
+#include <geode/stochastic/spatial/spatial_domain.hpp>
 
 namespace geode
 {
     template < index_t dimension >
     class UniformPointSetSampler : public ObjectSetSampler< Point< dimension > >
     {
     public:
-        UniformPointSetSampler( const BoundingBox< dimension >& box )
-            : ObjectSetSampler< Point< dimension > >{}, box_( box )
+        UniformPointSetSampler( const SpatialDomain< dimension >& domain )
+            : ObjectSetSampler< Point< dimension > >{}, domain_( domain )
         {
-            auto volume = box_.n_volume();
+            auto volume = domain_.extended_n_volume();
             OPENGEODE_EXCEPTION( volume != 0.,
-                "[PointSetSampler] - Undefined Bounding Box (volume ==0)." );
+                "[PointSetSampler] - Undefined Extended Bounding "
+                "Box (volume ==0)." );
             this->log_pdf_ = -std::log( volume );
+            step_move_ = define_step_for_move();
+            OPENGEODE_EXCEPTION( step_move_ > 0.,
+                "[PointSetSampler] - Undefined step length for move (value == ",
+                step_move_, ")." );
         }
 
         Point< dimension > sample( RandomEngine& engine ) const override
         {
-            return PointUniformSampler::sample< dimension >( engine, box_ );
+            return PointUniformSampler::sample< dimension >(
+                engine, domain_.extended_box() );
         }
 
         Point< dimension > change(
             const Point< dimension >& obj, RandomEngine& engine ) const override
         {
-            double ratio = 0.1;
-            geode::Sphere< dimension > ball{ obj,
-                ratio * std::get< 1 >( box_.smallest_length() ) };
+            geode::Sphere< dimension > ball{ obj, step_move_ };
 
             auto new_point =
                 PointUniformSampler::sample< dimension >( engine, ball );
             constexpr index_t max_try{ 100 };
             for( const auto try_id : geode::Range{ max_try } )
             {
-                if( box_.contains( new_point ) )
+                if( domain_.extended_contains( new_point ) )
                 {
                     return new_point;
                 }
                 new_point =
                     PointUniformSampler::sample< dimension >( engine, ball );
             }
-            throw OpenGeodeException( absl::StrCat(
-                "[PointSampler] - Cannot find a point in the box: ",
-                box_.string() ) );
+            throw OpenGeodeException(
+                absl::StrCat( "[PointSampler] - Cannot find a point in the "
+                              "extended domain" ) );
             return obj;
         }
 
     private:
+        double define_step_for_move()
+        {
+            double ratio = 0.1;
+            return ratio * domain_.smallest_length();
+        }
+
         bool is_valid_object( const Point< dimension >& obj ) const override
         {
-            return box_.contains( obj );
+            return domain_.extended_contains( obj );
         }
 
     private:
-        BoundingBox< dimension > box_;
+        const SpatialDomain< dimension >& domain_;
+        double step_move_{ 0. };
     };
 
 } // namespace geode

include/geode/stochastic/sampling/direct/object_set_sampler/segment_set_sampler.hpp

@@ -29,30 +29,32 @@
 #include <geode/stochastic/sampling/direct/object_set_sampler/object_set_sampler.hpp>
 #include <geode/stochastic/sampling/direct/point_uniform_sampler.hpp>
 #include <geode/stochastic/sampling/direct/segment_uniform_sampler.hpp>
+#include <geode/stochastic/spatial/spatial_domain.hpp>
 
 namespace geode
 {
     class UniformSegmentSetSampler : public ObjectSetSampler< OwnerSegment2D >
     {
     public:
-        UniformSegmentSetSampler( const BoundingBox< 2 >& box,
+        UniformSegmentSetSampler( const SpatialDomain< 2 >& domain,
             const DoubleSampler::Distribution& length,
             const DoubleSampler::Distribution& azimuth )
             : ObjectSetSampler< OwnerSegment2D >{},
-              box_{ box },
+              domain_{ domain },
               length_{ length },
               azimuth_{ azimuth }
         {
-            auto volume = box_.n_volume();
+            auto volume = domain_.extended_n_volume();
             OPENGEODE_EXCEPTION( volume != 0.,
-                "[SegmentSetSampler] - Undefined Bounding Box (volume == 0)." );
+                "[SegmentSetSampler] - Undefined Extended Bounding "
+                "Box (volume ==0)." );
             this->log_pdf_ = -std::log( volume );
         }
 
         OwnerSegment2D sample( RandomEngine& engine ) const override
         {
             auto seg = SegmentUniformSampler::sample(
-                engine, box_, length_, azimuth_ );
+                engine, domain_.extended_box(), length_, azimuth_ );
             return seg;
         }
 
@@ -72,8 +74,8 @@ namespace geode
             constexpr index_t max_try{ 100 };
             for( const auto try_id : geode::Range{ max_try } )
             {
-                if( box_.contains( new_point )
-                    || box_.contains( extremities[other] ) )
+                if( domain_.extended_contains( new_point )
+                    || domain_.extended_contains( extremities[other] ) )
                 {
                     OwnerSegment2D new_segment{ obj };
                     new_segment.set_point( current, new_point );
@@ -82,22 +84,20 @@ namespace geode
                 new_point = PointUniformSampler::sample< 2 >( engine, ball );
             }
             throw OpenGeodeException( absl::StrCat(
-                "[PointSampler] - Cannot find a point in the box: ",
-                box_.string() ) );
+                "[SegmentSetSampler] - Cannot find a point in the box" ) );
             return obj;
         }
 
     private:
         bool is_valid_object( const OwnerSegment2D& obj ) const override
         {
             const auto& extremities = obj.vertices();
-
-            return box_.contains( extremities[0] )
-                   || box_.contains( extremities[1] );
+            return domain_.extended_contains( extremities[0] )
+                   || domain_.extended_contains( extremities[1] );
         }
 
     private:
-        BoundingBox2D box_;
+        const SpatialDomain< 2 >& domain_;
         DoubleSampler::Distribution length_;
         DoubleSampler::Distribution azimuth_;
     };

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

@@ -32,6 +32,7 @@
 #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>
+#include <geode/stochastic/spatial/spatial_domain.hpp>
 
 namespace geode
 {
@@ -72,7 +73,10 @@ namespace geode
     class FractureSimulationRunner : public SimulationRunner< OwnerSegment2D >
     {
     public:
-        FractureSimulationRunner( const BoundingBox2D& box ) : box_( box ) {}
+        FractureSimulationRunner( const SpatialDomain< 2 >& domain )
+            : SimulationRunner< OwnerSegment2D >( domain )
+        {
+        }
 
         void add_fracture_set_descriptor(
             const FractureSetDescription& descriptor )
@@ -101,7 +105,7 @@ namespace geode
                         set_desc.azimuth );
                 this->set_samplers_.push_back(
                     std::make_unique< UniformSegmentSetSampler >(
-                        box_, length_distribution, azimuth_distribution ) );
+                        domain_, length_distribution, azimuth_distribution ) );
 
                 add_birth_death_change_moves( this->set_samplers_.back(),
                     *proposal_kernel, set_id, set_desc.birth_ratio,
@@ -117,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 )
                 {
@@ -133,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_ =
@@ -159,7 +164,6 @@ namespace geode
         }
 
     private:
-        BoundingBox2D box_;
         std::vector< FractureSetDescription > set_descriptors_;
     };
 

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

@@ -27,6 +27,7 @@
 #include <geode/stochastic/sampling/mcmc/helpers/simulation_printer.hpp>
 #include <geode/stochastic/sampling/mcmc/metropolis_hasting_sampler.hpp>
 #include <geode/stochastic/sampling/mcmc/models/energy_term_collection.hpp>
+#include <geode/stochastic/spatial/spatial_domain.hpp>
 
 #include <absl/strings/str_join.h>
 
@@ -62,7 +63,8 @@ namespace geode
     class SimulationRunner
     {
     public:
-        SimulationRunner() = default;
+        SimulationRunner( const SpatialDomain< ObjectType::dim >& domain )
+            : domain_( domain ) {};
         virtual ~SimulationRunner() = default;
 
         virtual void initialize() = 0;
@@ -153,6 +155,7 @@ namespace geode
         }
 
     protected:
+        SpatialDomain< ObjectType::dim > domain_;
         std::vector< std::unique_ptr< geode::ObjectSetSampler< ObjectType > > >
             set_samplers_;