feat(SpatialDomain): account for edge effects

Author: francoisbonneau

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,
@@ -159,7 +163,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_;
 

src/geode/stochastic/CMakeLists.txt

@@ -40,6 +40,7 @@ add_geode_library(
         "spatial/object_neighborhood.hpp"
         "spatial/object_helpers.hpp"
         "spatial/pairwise_interactions.hpp"
+        "spatial/spatial_domain.hpp"
         "spatial/details/RTree.hpp"
         "sampling/direct/object_set_sampler/object_set_sampler.hpp"
         "sampling/direct/object_set_sampler/point_set_sampler.hpp"

tests/stochastic/sampling/mcmc/proposal/test-proposal-kernel.cpp

@@ -56,8 +56,9 @@ namespace
         geode::BoundingBox2D box;
         box.add_point( min_point );
         box.add_point( max_point );
+        geode::SpatialDomain domain( box, 0. );
 
-        geode::UniformPointSetSampler< 2 > sampler( box );
+        geode::UniformPointSetSampler< 2 > sampler( domain );
 
         // Create classical birth-death-change kernel
         auto kernel = geode::create_birth_death_change_kernel< geode::Point2D >(

tests/stochastic/sampling/mcmc/test-metropolis-hasting-sampler.cpp

@@ -159,10 +159,11 @@ int main()
         geode::BoundingBox2D box;
         box.add_point( min_point );
         box.add_point( max_point );
+        geode::SpatialDomain domain( box, 0. );
 
         geode::ObjectSets< geode::Point2D > state;
         const auto set_id = state.add_set( "default_name" );
-        geode::UniformPointSetSampler< 2 > sampler( box );
+        geode::UniformPointSetSampler< 2 > sampler( domain );
 
         double birth_prob = 0.3;
         double death_prob = 0.1;

tests/stochastic/sampling/mcmc/test-mh-fractures.cpp

@@ -36,7 +36,8 @@ namespace
         geode::BoundingBox2D box;
         box.add_point( geode::Point2D{ { 0.0, 0.0 } } );
         box.add_point( geode::Point2D{ { 100.0, 100.0 } } );
-
+        // todo change
+        geode::SpatialDomain domain( box, 0. );
         // --- Object set
         geode::FractureSetDescription setA;
         setA.name = "A";
@@ -57,7 +58,7 @@ namespace
         setA.p20 = 0.05;
         setA.minimal_spacing = 1.;
 
-        geode::FractureSimulationRunner runner( box );
+        geode::FractureSimulationRunner runner( domain );
         runner.add_fracture_set_descriptor( setA );
 
         runner.initialize();
@@ -90,6 +91,8 @@ namespace
         geode::BoundingBox2D box;
         box.add_point( geode::Point2D{ { 0.0, 0.0 } } );
         box.add_point( geode::Point2D{ { 100.0, 100.0 } } );
+        // todo change
+        geode::SpatialDomain domain( box, 0. );
 
         // --- Object set
         geode::FractureSetDescription setA;
@@ -133,7 +136,7 @@ namespace
         setB.p20 = 0.05;
         setB.minimal_spacing = 2.;
 
-        geode::FractureSimulationRunner runner( box );
+        geode::FractureSimulationRunner runner( domain );
         runner.add_fracture_set_descriptor( setA );
         runner.add_fracture_set_descriptor( setB );
 

tests/stochastic/sampling/mcmc/test-mh-poisson.cpp

@@ -49,8 +49,8 @@ namespace
         : public geode::SimulationRunner< geode::Point2D >
     {
     public:
-        PoissonSimulationRunner( const geode::BoundingBox2D& box )
-            : box_( box ) {};
+        PoissonSimulationRunner( const geode::SpatialDomain< 2 >& domain )
+            : geode::SimulationRunner< geode::Point2D >( domain ){};
 
         void add_set_descriptor( const SetDescription& descriptor )
         {
@@ -79,7 +79,7 @@ namespace
 
                 this->set_samplers_.push_back(
                     std::make_unique< geode::UniformPointSetSampler< 2 > >(
-                        box_ ) );
+                        domain_ ) );
 
                 geode::add_birth_death_change_moves( this->set_samplers_.back(),
                     *proposal_kernel, set_id, set_desc.birth_ratio,
@@ -163,6 +163,7 @@ namespace
         geode::BoundingBox2D box;
         box.add_point( geode::Point2D{ { 0.0, 0.0 } } );
         box.add_point( geode::Point2D{ { 10.0, 10.0 } } );
+        geode::SpatialDomain domain( box, 0. );
 
         std::array< double, 4 > birth_ratio{ 0.1, 0.5, 2., 4. };
         std::array< double, 4 > change_ratio{ 0., 1., 1., 0. };
@@ -182,7 +183,7 @@ namespace
             densityA.density = 0.3;
             densityA.target_count = 30.0;
 
-            PoissonSimulationRunner runner( box );
+            PoissonSimulationRunner runner( domain );
             runner.add_set_descriptor( setA );
             runner.add_density_descriptor( densityA );
             runner.initialize();
@@ -216,6 +217,7 @@ namespace
         geode::BoundingBox2D box;
         box.add_point( geode::Point2D{ { 0.0, 0.0 } } );
         box.add_point( geode::Point2D{ { 10.0, 10.0 } } );
+        geode::SpatialDomain domain( box, 0. );
 
         // --- Set descriptions
         SetDescription set01{ "set01", 2.0, 3.0, 1.0 };
@@ -227,7 +229,7 @@ namespace
         PoissonDensityDescription density02{ "set02", 0.4, 40.0 };
         PoissonDensityDescription density03{ "set03", 0.3, 30.0 };
 
-        PoissonSimulationRunner runner( box );
+        PoissonSimulationRunner runner( domain );
         runner.add_set_descriptor( set01 );
         runner.add_set_descriptor( set02 );
         runner.add_set_descriptor( set03 );