fix(Neighborhood): use rtree search in pairwise interactions.

commitlint.config.js

@@ -1,4 +1,4 @@
-export default {
+const Configuration = {
   extends: ["@commitlint/config-angular"],
   rules: {
     "scope-empty": [2, "never"],
@@ -12,5 +12,8 @@ export default {
     "subject-full-stop": [0],
     "type-case": [0],
     "type-empty": [0],
+    "type-enum": [2, "always", ["feat", "fix", "perf"]],
   },
 }
+
+export default Configuration

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

@@ -117,8 +117,7 @@ namespace geode
                     this->energy_terms_collection_.add_energy_term(
                         std::make_unique< DensityTerm< OwnerSegment2D > >(
                             absl::StrCat( set_desc.name, "_density" ),
-                            set_desc.p20,
-                            absl::flat_hash_set< uuid >{ set_id } ) ) );
+                            set_desc.p20, std::vector< uuid >{ set_id } ) ) );
                 // spacing
                 if( set_desc.minimal_spacing < GLOBAL_EPSILON )
                 {
@@ -133,7 +132,7 @@ namespace geode
                     this->energy_terms_collection_.add_energy_term(
                         std::make_unique< PairwiseTerm< OwnerSegment2D > >(
                             absl::StrCat( set_desc.name, "_min_spacing" ), 0.,
-                            absl::flat_hash_set< uuid >{ set_id },
+                            std::vector< uuid >{ set_id },
                             std::move( interaction ) ) ) );
             }
 

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

@@ -35,7 +35,7 @@ namespace geode
     public:
         explicit DensityTerm( std::string_view name,
             double lambda,
-            absl::flat_hash_set< uuid > targeted_set_ids )
+            std::vector< uuid > targeted_set_ids )
             : EnergyTerm< ObjectType >(
                   name, lambda, std::move( targeted_set_ids ) )
         {

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

@@ -88,10 +88,11 @@ namespace geode
     public:
         explicit EnergyTerm( std::string_view name,
             double param,
-            absl::flat_hash_set< uuid >&& targeted_set_ids )
+            std::vector< uuid >&& targeted_set_ids )
             : energy_scale_{ param },
               targeted_set_ids_{ std::move( targeted_set_ids ) }
         {
+            std::sort( targeted_set_ids_.begin(), targeted_set_ids_.end() );
             IdentifierBuilder builder( *this );
             builder.set_name( name );
         }
@@ -103,7 +104,7 @@ namespace geode
             return energy_scale_.parameter();
         }
 
-        const absl::flat_hash_set< uuid >& targeted_set_ids() const
+        const std::vector< uuid >& targeted_set_ids() const
         {
             return targeted_set_ids_;
         }
@@ -147,7 +148,8 @@ namespace geode
     protected:
         bool is_targeted_set( const uuid& set_id ) const
         {
-            return targeted_set_ids_.find( set_id ) != targeted_set_ids_.end();
+            return std::binary_search(
+                targeted_set_ids_.begin(), targeted_set_ids_.end(), set_id );
         }
 
         template < typename Func >
@@ -166,6 +168,6 @@ namespace geode
 
     private:
         detail::EnergyScale energy_scale_;
-        absl::flat_hash_set< uuid > targeted_set_ids_;
+        std::vector< uuid > targeted_set_ids_;
     };
 } // namespace geode

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

@@ -37,7 +37,7 @@ namespace geode
     public:
         explicit PairwiseTerm( std::string_view name,
             double gamma,
-            absl::flat_hash_set< uuid > targeted_set_ids,
+            std::vector< uuid > targeted_set_ids,
             std::unique_ptr< PairwiseInteraction< ObjectType > > interaction )
             : EnergyTerm< ObjectType >(
                   name, gamma, std::move( targeted_set_ids ) ),
@@ -59,8 +59,9 @@ namespace geode
                 return 0.0;
             }
             double delta = 0.0;
-            const auto neighbors = state.neighbors( new_object.object,
-                interaction_->neighborhood_searching_distance() );
+            const auto neighbors =
+                state.neighbors( new_object.object, this->targeted_set_ids(),
+                    interaction_->neighborhood_searching_distance() );
             for( const auto& neigh_id : neighbors )
             {
                 geode::ObjectRef< ObjectType > neigh_object{
@@ -83,8 +84,9 @@ namespace geode
             ObjectRef< ObjectType > object_to_remove{
                 state.get_object( object_id ), object_id.set_id
             };
-            const auto neighbors = state.neighbors(
-                object_id, interaction_->neighborhood_searching_distance() );
+            const auto neighbors =
+                state.neighbors( object_id, this->targeted_set_ids(),
+                    interaction_->neighborhood_searching_distance() );
             for( auto neigh_id : neighbors )
             {
                 ObjectRef< ObjectType > neigh_object{
@@ -111,8 +113,9 @@ namespace geode
             ObjectRef< ObjectType > object_to_remove{
                 state.get_object( old_object_id ), old_object_id.set_id
             };
-            const auto old_neighbors = state.neighbors( old_object_id,
-                interaction_->neighborhood_searching_distance() );
+            const auto old_neighbors =
+                state.neighbors( old_object_id, this->targeted_set_ids(),
+                    interaction_->neighborhood_searching_distance() );
             for( auto neigh_id : old_neighbors )
             {
                 ObjectRef< ObjectType > neigh_object{
@@ -123,8 +126,9 @@ namespace geode
             }
 
             // Add new object's interactions
-            const auto new_neighbors = state.neighbors( new_object.object,
-                interaction_->neighborhood_searching_distance() );
+            const auto new_neighbors =
+                state.neighbors( new_object.object, this->targeted_set_ids(),
+                    interaction_->neighborhood_searching_distance() );
             for( auto neigh_id : new_neighbors )
             {
                 if( old_object_id == neigh_id )
@@ -143,36 +147,33 @@ namespace geode
         double statistic( const ObjectSets< ObjectType >& state ) const override
         {
             double sum = 0.0;
-            this->for_each_targeted_object( state, [&]( const ObjectId&
-                                                           obj_id ) {
-                const auto& cur_obj = state.get_object( obj_id );
-                const auto neighbors =
-                    state.get_all_object(); // state.neighbors( obj_id, 1.1 );
-                for( const auto& neigh_obj_id : neighbors )
-                {
-                    // if( neigh_obj_id.set_id < obj_id.set_id )
-                    //{
-                    //     continue;
-                    // }
-                    // if( neigh_obj_id.set_id == obj_id.set_id
-                    //     && neigh_obj_id.index <= obj_id.index )
-                    //{
-                    //     continue;
-                    // }
-
-                    if( neigh_obj_id == obj_id )
-                    {
-                        continue;
-                    }
+            this->for_each_targeted_object(
+                state, [&]( const ObjectId& obj_id ) {
+                    const auto& cur_obj = state.get_object( obj_id );
                     ObjectRef< ObjectType > object{ cur_obj, obj_id.set_id };
-                    ObjectRef< ObjectType > neigh_object{
-                        state.get_object( neigh_obj_id ), neigh_obj_id.set_id
-                    };
+                    const auto neighbors =
+                        state.neighbors( obj_id, this->targeted_set_ids(),
+                            interaction_->neighborhood_searching_distance() );
+                    for( const auto& neigh_obj_id : neighbors )
+                    {
+                        if( neigh_obj_id.set_id < obj_id.set_id )
+                        {
+                            continue;
+                        }
+                        if( neigh_obj_id.set_id == obj_id.set_id
+                            && neigh_obj_id.index <= obj_id.index )
+                        {
+                            continue;
+                        }
+                        ObjectRef< ObjectType > neigh_object{
+                            state.get_object( neigh_obj_id ),
+                            neigh_obj_id.set_id
+                        };
 
-                    sum += interaction_->evaluate( object, neigh_object );
-                }
-            } );
-            return sum / 2.;
+                        sum += interaction_->evaluate( object, neigh_object );
+                    }
+                } );
+            return sum;
         }
 
     private:

include/geode/stochastic/spatial/object_neighborhood.hpp

@@ -64,6 +64,7 @@ namespace geode
 
         std::vector< ObjectId > get_all_neighbor_ids(
             const BoundingBox< dimension >& box,
+            const std::vector< uuid >& targeted_set_ids,
             std::optional< ObjectId > exclude_self_id ) const;
 
     private:

include/geode/stochastic/spatial/object_sets.hpp

@@ -70,12 +70,14 @@ namespace geode
         void remove_object( const ObjectId& object_id );
 
         // Object neighbor search by ObjectId (always excludes self)
-        std::vector< ObjectId > neighbors(
-            const ObjectId& object_id, double searching_distance ) const;
+        std::vector< ObjectId > neighbors( const ObjectId& object_id,
+            const std::vector< uuid >& targeted_set_ids,
+            double searching_distance ) const;
         // Object neighbor search by arbitrary object (return self if in the
         // object_set)
-        std::vector< ObjectId > neighbors(
-            const Type& object, double searching_distance ) const;
+        std::vector< ObjectId > neighbors( const Type& object,
+            const std::vector< uuid >& targeted_set_ids,
+            double searching_distance ) const;
 
         std::string string() const;
 

src/geode/stochastic/spatial/object_neighborhood.cpp

@@ -90,12 +90,19 @@ namespace geode
     std::vector< ObjectId >
         ObjectNeighborhood< dimension >::get_all_neighbor_ids(
             const BoundingBox< dimension >& box,
+            const std::vector< uuid >& targeted_set_ids,
             std::optional< ObjectId > exclude_self_id ) const
     {
         auto [min_bound, max_bound] = bounding_box_bounds< dimension >( box );
         std::vector< ObjectId > res;
         tree_.Search( min_bound.data(), max_bound.data(),
-            [&res, &exclude_self_id]( const ObjectId& cur_id ) -> bool {
+            [&res, &exclude_self_id, &targeted_set_ids](
+                const ObjectId& cur_id ) -> bool {
+                if( !std::binary_search( targeted_set_ids.begin(),
+                        targeted_set_ids.end(), cur_id.set_id ) )
+                {
+                    return true;
+                }
                 if( exclude_self_id && exclude_self_id.value() == cur_id )
                 {
                     return true;

src/geode/stochastic/spatial/object_sets.cpp

@@ -123,20 +123,25 @@ namespace geode
 
     template < typename Type >
     std::vector< ObjectId > ObjectSets< Type >::neighbors(
-        const ObjectId& object_id, double searching_distance ) const
+        const ObjectId& object_id,
+        const std::vector< uuid >& targeted_set_ids,
+        double searching_distance ) const
     {
         auto box = object_bounding_box( get_object( object_id ) );
         box.extends( searching_distance * 2. );
-        return neighborhood_.get_all_neighbor_ids( box, object_id );
+        return neighborhood_.get_all_neighbor_ids(
+            box, targeted_set_ids, object_id );
     }
 
     template < typename Type >
-    std::vector< ObjectId > ObjectSets< Type >::neighbors(
-        const Type& object, double searching_distance ) const
+    std::vector< ObjectId > ObjectSets< Type >::neighbors( const Type& object,
+        const std::vector< uuid >& targeted_set_ids,
+        double searching_distance ) const
     {
         auto box = object_bounding_box( object );
         box.extends( searching_distance * 2. );
-        return neighborhood_.get_all_neighbor_ids( box, std::nullopt );
+        return neighborhood_.get_all_neighbor_ids(
+            box, targeted_set_ids, std::nullopt );
     }
 
     template < typename Type >

tests/stochastic/sampling/mcmc/models/test-gibbs-energy.cpp

@@ -31,30 +31,21 @@
 #include <geode/stochastic/sampling/mcmc/models/gibbs_energy.hpp>
 #include <geode/stochastic/spatial/object_sets.hpp>
 
-namespace
+void test_gibbs_energy()
 {
-    geode::ObjectSets< geode::Point2D > create_object_set(
-        const geode::uuid& set_id )
-    {
-        geode::Point2D p1{ { 0., 0. } };
-        geode::Point2D p2{ { 1., 1. } };
-
-        geode::ObjectSets< geode::Point2D > pattern;
-        pattern.add_object( std::move( p1 ), set_id );
-        pattern.add_object( std::move( p2 ), set_id );
-
-        return pattern;
-    }
-} // namespace
+    geode::ObjectSets< geode::Point2D > pattern;
+    const auto set_id = pattern.add_set( "default_name" );
+    geode::Point2D p1{ { 0., 0. } };
+    geode::Point2D p2{ { 1., 1. } };
+    pattern.add_object( std::move( p1 ), set_id );
+    pattern.add_object( std::move( p2 ), set_id );
 
-void test_gibbs_energy( const geode::uuid& set_id )
-{
     geode::EnergyTermCollection< geode::Point2D > energy_terms;
 
     // Add intensity term
     energy_terms.add_energy_term(
         std::make_unique< geode::DensityTerm< geode::Point2D > >(
-            "intensity", 0.5, absl::flat_hash_set< geode::uuid >{ set_id } ) );
+            "intensity", 0.5, std::vector< geode::uuid >{ set_id } ) );
 
     // Add pairwise term with trivial interaction: always counts 1 for each pair
     auto interaction =
@@ -63,16 +54,14 @@ void test_gibbs_energy( const geode::uuid& set_id )
 
     energy_terms.add_energy_term(
         std::make_unique< geode::PairwiseTerm< geode::Point2D > >(
-            "interaction", 0.8, absl::flat_hash_set< geode::uuid >{ set_id },
+            "interaction", 0.8, std::vector< geode::uuid >{ set_id },
             std::move( interaction ) ) );
 
     OPENGEODE_EXCEPTION( energy_terms.size() == 2,
         "[test gibbs] Wrong number of components after adding terms." );
 
     geode::GibbsEnergy< geode::Point2D > gibbs_energy( energy_terms );
 
-    auto pattern = create_object_set( set_id );
-
     // Check total log-energy is finite
     double total_energy = gibbs_energy.total_log_energy( pattern );
     OPENGEODE_EXCEPTION( std::isfinite( total_energy ),
@@ -99,7 +88,7 @@ void test_gibbs_energy( const geode::uuid& set_id )
 
     // Clear components and verify
     energy_terms.clear();
-    OPENGEODE_EXCEPTION( energy_terms.size() == 2,
+    OPENGEODE_EXCEPTION( energy_terms.size() == 0,
         "[test gibbs] Components not cleared properly." );
 }
 
@@ -108,9 +97,9 @@ int main()
     try
     {
         geode::StochasticLibrary::initialize();
-        geode::uuid set_id;
-
-        // test_gibbs_energy( set_id );
+        test_gibbs_energy();
+        geode::Logger::info( "MH TEST GIBBS ENERGY SUCCESS" );
+        return 0;
     }
     catch( ... )
     {

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

@@ -174,8 +174,7 @@ int main()
         // Add intensity term
         energy_terms.add_energy_term(
             std::make_unique< geode::DensityTerm< geode::Point2D > >(
-                "intensity", 0.5,
-                absl::flat_hash_set< geode::uuid >{ set_id } ) );
+                "intensity", 0.5, std::vector< geode::uuid >{ set_id } ) );
 
         geode::MetropolisHastings< geode::Point2D > mh(
             energy_terms, std::move( kernel ) );

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

@@ -97,7 +97,7 @@ namespace
                             geode::DensityTerm< geode::Point2D > >(
                             absl::StrCat( energy_desc.name, "_density" ),
                             energy_desc.density,
-                            absl::flat_hash_set< geode::uuid >{ set_id } ) ) );
+                            std::vector< geode::uuid >{ set_id } ) ) );
 
                 this->ordered_target_statistics_.push_back(
                     energy_desc.target_count );