feat(DeterministicObject): add the possibility to constrained simulation with determininstic object.

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

@@ -36,6 +36,11 @@ namespace geode
             .def_readwrite( "length", &FractureSetDescription::length )
             .def_readwrite( "azimuth", &FractureSetDescription::azimuth )
             .def_readwrite( "p20", &FractureSetDescription::p20 )
+            .def( "add_observed_fracture",
+                []( FractureSetDescription& self, const geode::Point2D& a,
+                    const geode::Point2D& b ) {
+                    self.observed_fractures.push_back( { a, b } );
+                } )
             .def_readwrite(
                 "minimal_spacing", &FractureSetDescription::minimal_spacing )
             .def_readwrite(

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

@@ -46,6 +46,7 @@ def test_fracture_simulator():
     # --- Object set
     setA = stochastic.FractureSetDescription()
     setA.name = "A"
+    setA.add_observed_fracture(og.Point2D([10.0, 10.0]), og.Point2D([20.0, 20.0]))
 
     # length
     setA.length.distribution_type =stochastic.DistributionType("UniformClosed")

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

@@ -34,6 +34,8 @@
 #include <geode/stochastic/spatial/pairwise_interactions.hpp>
 #include <geode/stochastic/spatial/spatial_domain.hpp>
 
+#include <geode/geometry/basic_objects/segment.hpp>
+
 namespace geode
 {
     struct FractureSetDescription
@@ -45,6 +47,7 @@ namespace geode
 
         // positionning
         double p20;
+        std::vector< std::array< geode::Point2D, 2 > > observed_fractures;
         // double p21;
         double minimal_spacing{ 0. };
 
@@ -96,6 +99,13 @@ namespace geode
             for( const auto& set_desc : set_descriptors_ )
             {
                 const auto set_id = this->object_sets_.add_set( set_desc.name );
+                for( const auto& fixed_object : set_desc.observed_fractures )
+                {
+                    this->object_sets_.add_fixed_object(
+                        geode::OwnerSegment2D{
+                            fixed_object[0], fixed_object[1] },
+                        set_id );
+                }
                 name_to_uuid[set_desc.name] = set_id;
 
                 auto length_distribution =

include/geode/stochastic/sampling/mcmc/metropolis_hasting_sampler.hpp

@@ -192,7 +192,7 @@ namespace geode
                 gibbs_energy_.delta_log_add( state, new_object );
             return accept_or_reject( proposal, state, engine, delta_log_energy,
                 []( auto& state, auto& proposal ) {
-                    state.add_object(
+                    state.add_free_object(
                         std::move( proposal.proposed_move.new_object.value() ),
                         proposal.set_id );
                 } );
@@ -207,7 +207,7 @@ namespace geode
                 gibbs_energy_.delta_log_remove( state, old_object_id );
             return accept_or_reject( proposal, state, engine, delta_log_energy,
                 []( auto& state, auto& proposal ) {
-                    state.remove_object( proposal.old_object_id() );
+                    state.remove_free_object( proposal.old_object_id() );
                 } );
         };
 
@@ -221,7 +221,7 @@ namespace geode
                 state, old_object_id, new_object );
             return accept_or_reject( proposal, state, engine, delta_log_energy,
                 []( auto& state, auto& proposal ) {
-                    state.update_object( proposal.old_object_id(),
+                    state.update_free_object( proposal.old_object_id(),
                         std::move(
                             proposal.proposed_move.new_object.value() ) );
                 } );

include/geode/stochastic/sampling/mcmc/proposal/moves.hpp

@@ -138,18 +138,19 @@ namespace geode
         virtual std::string string() const = 0;
 
     protected:
-        std::optional< geode::index_t > draw_a_sample_id(
+        std::optional< geode::index_t > draw_a_free_sample_id(
             const ObjectSet< ObjectType >& set, RandomEngine& engine ) const
         {
-            if( set.empty() )
+            const auto max_obj_id = set.nb_free_objects();
+            if( max_obj_id == 0 )
             {
                 return std::nullopt;
             }
-            const auto max_obj_id = set.size();
             geode::UniformClosed< index_t > uniform_closed_index_t;
             uniform_closed_index_t.min_value = 0;
             uniform_closed_index_t.max_value = max_obj_id - 1;
-            return engine.sample_uniform( uniform_closed_index_t );
+            return set.nb_fixed_objects()
+                   + engine.sample_uniform( uniform_closed_index_t );
         }
 
     protected:
@@ -212,23 +213,23 @@ namespace geode
             birth.proposal_probabilities.log_forward_prob =
                 log_p_birth_ + this->sampler_.log_pdf( new_obj );
             birth.proposal_probabilities.log_backward_prob =
-                log_p_death_ - std::log( set.size() + 1.0 );
+                log_p_death_ - std::log( set.nb_objects() + 1.0 );
             return birth;
         }
 
         MoveResult< ObjectType > propose_death_move(
             const ObjectSet< ObjectType >& set, RandomEngine& engine ) const
         {
             MoveResult< ObjectType > death;
-            death.old_object_id = this->draw_a_sample_id( set, engine );
+            death.old_object_id = this->draw_a_free_sample_id( set, engine );
             if( !death.old_object_id.has_value() )
             {
                 return death;
             }
             const auto cur_object_id = death.old_object_id.value();
             death.type = MoveType::Death;
             death.proposal_probabilities.log_forward_prob =
-                log_p_death_ - std::log( set.size() );
+                log_p_death_ - std::log( set.nb_free_objects() );
             death.proposal_probabilities.log_backward_prob =
                 log_p_birth_
                 + this->sampler_.log_pdf( set.get_object( cur_object_id ) );
@@ -256,7 +257,7 @@ namespace geode
             RandomEngine& engine ) const override
         {
             MoveResult< ObjectType > change;
-            change.old_object_id = this->draw_a_sample_id( set, engine );
+            change.old_object_id = this->draw_a_free_sample_id( set, engine );
             if( !change.old_object_id.has_value() )
             {
                 return change;

include/geode/stochastic/spatial/object_set.hpp

@@ -39,16 +39,27 @@ namespace geode
         void set_name( std::string_view name );
         const Type& get_object( index_t index ) const;
 
-        index_t add_object( Type&& object );
-        void update_object( index_t index, Type&& object );
-        void remove_object( index_t index );
+        index_t add_fixed_object( Type&& object );
+        void remove_fixed_object( index_t index );
+
+        index_t add_free_object( Type&& object );
+        void update_free_object( index_t index, Type&& object );
+        void remove_free_object( index_t index );
+
+        index_t nb_objects() const;
+        index_t nb_fixed_objects() const;
+        index_t nb_free_objects() const;
 
-        index_t size() const;
         bool empty() const;
+        bool is_fixed( index_t index ) const;
 
         std::string string() const;
 
+    private:
+        void do_remove_object( index_t index );
+
     private:
         std::vector< Type > objects_;
+        index_t first_free_object_{ 0 };
     };
 } // namespace geode

include/geode/stochastic/spatial/object_sets.hpp

@@ -65,9 +65,11 @@ namespace geode
         index_t nb_objects() const;
 
         uuid add_set( std::string_view name );
-        ObjectId add_object( Type&& object, const uuid& set_id );
-        void update_object( const ObjectId& object_id, Type&& object );
-        void remove_object( const ObjectId& object_id );
+        ObjectId add_fixed_object( Type&& object, const uuid& set_id );
+        ObjectId add_free_object( Type&& object, const uuid& set_id );
+        void update_free_object( const ObjectId& object_id, Type&& object );
+        void remove_free_object( const ObjectId& object_id );
+        void remove_fixed_object( const ObjectId& object_id );
 
         // Object neighbor search by ObjectId (always excludes self)
         std::vector< ObjectId > neighbors( const ObjectId& object_id,
@@ -83,6 +85,7 @@ namespace geode
 
     private:
         ObjectSet< Type >& get_set( const uuid& set_id );
+        void update_neighborhood_remove_context( const ObjectId& object_id );
 
     private:
         absl::flat_hash_map< uuid, ObjectSet< Type > > sets_;

src/geode/stochastic/spatial/object_set.cpp

@@ -46,49 +46,102 @@ namespace geode
     }
 
     template < typename Type >
-    index_t ObjectSet< Type >::add_object( Type&& object )
+    index_t ObjectSet< Type >::add_fixed_object( Type&& object )
+    {
+        auto new_fixed_object_id = first_free_object_;
+        first_free_object_++;
+        objects_.push_back( std::move( object ) );
+        if( new_fixed_object_id != objects_.size() - 1 )
+        {
+            std::swap( objects_[new_fixed_object_id], objects_.back() );
+        }
+        return new_fixed_object_id;
+    }
+
+    template < typename Type >
+    index_t ObjectSet< Type >::add_free_object( Type&& object )
     {
         objects_.push_back( std::move( object ) );
         return objects_.size() - 1;
     }
 
     template < typename Type >
-    void ObjectSet< Type >::update_object( index_t index, Type&& object )
+    void ObjectSet< Type >::update_free_object( index_t index, Type&& object )
     {
+        OPENGEODE_EXCEPTION( index >= first_free_object_,
+            "[ObjectSet] - cannot update fixed object." );
         OPENGEODE_EXCEPTION( index < objects_.size(),
             "[ObjectSet] - object index out of range." );
         objects_[index] = std::move( object );
     }
 
     template < typename Type >
-    void ObjectSet< Type >::remove_object( index_t index )
+    void ObjectSet< Type >::remove_free_object( index_t index )
     {
-        OPENGEODE_EXCEPTION( index < objects_.size(),
-            "[ObjectSet] - object index out of range." );
-        if( index != objects_.size() - 1 )
+        OPENGEODE_EXCEPTION( index >= first_free_object_,
+            "[ObjectSet] - cannot remove fixed object." );
+        do_remove_object( index );
+    }
+
+    template < typename Type >
+    void ObjectSet< Type >::remove_fixed_object( index_t index )
+    {
+        OPENGEODE_EXCEPTION( index < first_free_object_,
+            "[ObjectSet] - cannot remove free object." );
+        do_remove_object( index );
+        first_free_object_--;
+    }
+
+    template < typename Type >
+    void ObjectSet< Type >::do_remove_object( index_t index )
+    {
+        const index_t last = objects_.size() - 1;
+        OPENGEODE_EXCEPTION(
+            index <= last, "[ObjectSet] - object index out of range." );
+        if( index != last )
         {
-            objects_[index] = std::move( objects_.back() );
+            std::swap( objects_[index], objects_[last] );
         }
         objects_.pop_back();
     }
 
     template < typename Type >
-    index_t ObjectSet< Type >::size() const
+    index_t ObjectSet< Type >::nb_objects() const
     {
         return objects_.size();
     }
 
+    template < typename Type >
+    index_t ObjectSet< Type >::nb_fixed_objects() const
+    {
+        return first_free_object_;
+    }
+
+    template < typename Type >
+    index_t ObjectSet< Type >::nb_free_objects() const
+    {
+        return objects_.size() - first_free_object_;
+    }
+
     template < typename Type >
     bool ObjectSet< Type >::empty() const
     {
         return objects_.empty();
     }
 
+    template < typename Type >
+    bool ObjectSet< Type >::is_fixed( index_t index ) const
+    {
+        return index < first_free_object_;
+    }
+
     template < typename Type >
     std::string ObjectSet< Type >::string() const
     {
         return absl::StrCat( "ObjectSet ", this->name(), " (",
-            this->id().string(), ") contains ", size(), " objects" );
+            this->id().string(), ") contains ", nb_objects(),
+            " objects (fixed: ", nb_fixed_objects(),
+            " - free: ", nb_free_objects(), ")" );
     }
 
     // Explicit template instantiation