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_object(
+                        geode::OwnerSegment2D{
+                            fixed_object[0], fixed_object[1] },
+                        set_id, true );
+                }
                 name_to_uuid[set_desc.name] = set_id;
 
                 auto length_distribution =

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

@@ -194,7 +194,7 @@ namespace geode
                 []( auto& state, auto& proposal ) {
                     state.add_object(
                         std::move( proposal.proposed_move.new_object.value() ),
-                        proposal.set_id );
+                        proposal.set_id, false );
                 } );
         };
 
@@ -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/models/components/energy_term.hpp

@@ -175,9 +175,9 @@ namespace geode
             for( const auto& targeted_set_id : targeted_set_ids_ )
             {
                 for( const auto id :
-                    geode::Range{ state.nb_objects_in_set( targeted_set_id ) } )
+                    state.get_objects_in_set( targeted_set_id ) )
                 {
-                    do_apply( ObjectId{ id, targeted_set_id } );
+                    do_apply( id );
                 }
             }
         }

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

@@ -138,14 +138,14 @@ 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;
@@ -212,26 +212,27 @@ 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 ) );
+                + this->sampler_.log_pdf(
+                    set.get_free_object( cur_object_id ) );
             return death;
         }
 
@@ -256,14 +257,14 @@ 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;
             }
             change.type = MoveType::Change;
             const auto& object_to_change =
-                set.get_object( change.old_object_id.value() );
+                set.get_free_object( change.old_object_id.value() );
             change.new_object =
                 this->sampler_.change( object_to_change, engine );
             change.proposal_probabilities.log_forward_prob =

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

@@ -46,7 +46,8 @@ namespace geode
         {
             OPENGEODE_EXCEPTION( proposed_move.old_object_id.has_value(),
                 "[Proposal] Proposal has no old_object_id" );
-            return ObjectId{ proposed_move.old_object_id.value(), set_id };
+            return ObjectId{ proposed_move.old_object_id.value(), false,
+                set_id };
         };
 
         std::string string() const

include/geode/stochastic/spatial/object_neighborhood.hpp

@@ -35,14 +35,17 @@ namespace geode
     struct ObjectId
     {
         index_t index;
+        bool fixed;
         uuid set_id;
         bool operator==( const ObjectId& other ) const noexcept
         {
-            return index == other.index && set_id == other.set_id;
+            return index == other.index && fixed == other.fixed
+                   && set_id == other.set_id;
         }
         bool operator!=( const ObjectId& other ) const noexcept
         {
-            return index != other.index || set_id != other.set_id;
+            return index != other.index || fixed != other.fixed
+                   || set_id != other.set_id;
         }
     };
 

include/geode/stochastic/spatial/object_set.hpp

@@ -37,18 +37,25 @@ namespace geode
         ObjectSet& operator=( ObjectSet&& ) noexcept = default;
 
         void set_name( std::string_view name );
-        const Type& get_object( index_t index ) const;
+        const Type& get_fixed_object( index_t index ) const;
+        const Type& get_free_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 );
+
+        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;
 
         std::string string() const;
 
     private:
-        std::vector< Type > objects_;
+        std::vector< Type > fixed_objects_;
+        std::vector< Type > free_objects_;
     };
 } // namespace geode

include/geode/stochastic/spatial/object_sets.hpp

@@ -59,15 +59,16 @@ namespace geode
         const ObjectSet< Type >& get_set( const uuid& set_id ) const;
         const Type& get_object( const ObjectId& object_id ) const;
         std::vector< ObjectId > get_all_object() const;
+        std::vector< ObjectId > get_objects_in_set( const uuid& set_id ) const;
 
         index_t nb_sets() const;
         index_t nb_objects_in_set( const uuid& set_id ) const;
         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_object( Type&& object, const uuid& set_id, bool fixed );
+        void update_free_object( const ObjectId& object_id, Type&& object );
+        void remove_free_object( const ObjectId& object_id );
 
         // Object neighbor search by ObjectId (always excludes self)
         std::vector< ObjectId > neighbors( const ObjectId& object_id,
@@ -83,6 +84,8 @@ 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

@@ -37,58 +37,89 @@ namespace geode
         IdentifierBuilder builder( *this );
         builder.set_name( name );
     }
+
+    template < typename Type >
+    const Type& ObjectSet< Type >::get_fixed_object( index_t index ) const
+    {
+        OPENGEODE_EXCEPTION( index < fixed_objects_.size(),
+            "[ObjectSet] - index for fixed object out of range." );
+        return fixed_objects_[index];
+    }
+
+    template < typename Type >
+    const Type& ObjectSet< Type >::get_free_object( index_t index ) const
+    {
+        OPENGEODE_EXCEPTION( index < free_objects_.size(),
+            "[ObjectSet] - index for free object out of range." );
+        return free_objects_[index];
+    }
+
     template < typename Type >
-    const Type& ObjectSet< Type >::get_object( index_t index ) const
+    index_t ObjectSet< Type >::add_fixed_object( Type&& object )
     {
-        OPENGEODE_EXCEPTION( index < objects_.size(),
-            "[ObjectSet] - object index out of range." );
-        return objects_[index];
+        fixed_objects_.push_back( std::move( object ) );
+        return fixed_objects_.size() - 1;
     }
 
     template < typename Type >
-    index_t ObjectSet< Type >::add_object( Type&& object )
+    index_t ObjectSet< Type >::add_free_object( Type&& object )
     {
-        objects_.push_back( std::move( object ) );
-        return objects_.size() - 1;
+        free_objects_.push_back( std::move( object ) );
+        return free_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 < objects_.size(),
-            "[ObjectSet] - object index out of range." );
-        objects_[index] = std::move( object );
+        OPENGEODE_EXCEPTION( index < free_objects_.size(),
+            "[ObjectSet] - free object index out of range." );
+        free_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 )
+        const index_t last = free_objects_.size() - 1;
+        OPENGEODE_EXCEPTION(
+            index <= last, "[ObjectSet] - free object index out of range." );
+        if( index != last )
         {
-            objects_[index] = std::move( objects_.back() );
+            std::swap( free_objects_[index], free_objects_[last] );
         }
-        objects_.pop_back();
+        free_objects_.pop_back();
+    }
+
+    template < typename Type >
+    index_t ObjectSet< Type >::nb_objects() const
+    {
+        return free_objects_.size() + fixed_objects_.size();
+    }
+
+    template < typename Type >
+    index_t ObjectSet< Type >::nb_fixed_objects() const
+    {
+        return fixed_objects_.size();
     }
 
     template < typename Type >
-    index_t ObjectSet< Type >::size() const
+    index_t ObjectSet< Type >::nb_free_objects() const
     {
-        return objects_.size();
+        return free_objects_.size();
     }
 
     template < typename Type >
     bool ObjectSet< Type >::empty() const
     {
-        return objects_.empty();
+        return free_objects_.empty() && fixed_objects_.empty();
     }
 
     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