fix straus test

Author: francoisbonneau

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

@@ -82,22 +82,6 @@ namespace geode
 
 namespace geode
 {
-    //    struct EnergyTermDescription
-    //    {
-    //        geode::uuid id;
-    //        std::string name;
-    //        std::string type;
-    //        double parameter_value;
-    //        std::optional< uuid > targeted_set_id{};
-    //    }
-    //
-    //    struct StatisticalDescription
-    //    {
-    //        std::string label;
-    //        double value;
-    //        std::optional< uuid > targeted_set_id{};
-    //    };
-
     template < typename ObjectType >
     class EnergyTerm : public Identifier
     {

tests/stochastic/CMakeLists.txt

@@ -129,14 +129,14 @@ add_geode_test(
          ${PROJECT_NAME}::stochastic
 )
 
-#add_geode_test(
-#    SOURCE  "sampling/mcmc/test-mh-strauss.cpp"
-#    DEPENDENCIES
-#        OpenGeode::basic
-#        OpenGeode::geometry
-#        ${PROJECT_NAME}::stochastic
-#    UNSTABLE
-#)
+add_geode_test(
+    SOURCE  "sampling/mcmc/test-mh-strauss.cpp"
+    DEPENDENCIES
+        OpenGeode::basic
+        OpenGeode::geometry
+        ${PROJECT_NAME}::stochastic
+    UNSTABLE
+)
 
 
 add_geode_test(

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

@@ -30,61 +30,80 @@
 #include <geode/stochastic/spatial/object_sets.hpp>
 namespace
 {
-    struct PoissonSetDescription
+    struct SetDescription
     {
         std::string name;
-        double density;
-        double target_count;
-
         double birth_ratio{ 1.0 };
         double death_ratio{ 1.0 };
         double change_ratio{ 1.0 };
     };
 
+    struct PoissonDensityDescription
+    {
+        std::string set_name;
+        double density;
+        double target_count;
+    };
+
     class PoissonSimulationRunner
         : public geode::SimulationRunner< geode::Point2D >
     {
     public:
-        PoissonSimulationRunner( const geode::BoundingBox2D& box )
-            : box_( box ){};
+        PoissonSimulationRunner( const geode::BoundingBox2D& box ) : box_( box )
+        {
+        }
 
-        void add_set_descriptor( const PoissonSetDescription& descriptor )
+        void add_set_descriptor( const SetDescription& descriptor )
         {
             set_descriptors_.push_back( descriptor );
         }
 
+        void add_density_descriptor(
+            const PoissonDensityDescription& descriptor )
+        {
+            density_descriptors_.push_back( descriptor );
+        }
+
         void initialize() override
         {
             auto proposal_kernel =
                 std::make_unique< geode::ProposalKernel< geode::Point2D > >();
 
-            for( const auto& descriptor : set_descriptors_ )
+            // Mapping set names -> UUID
+            std::unordered_map< std::string, geode::uuid > name_to_uuid;
+
+            // Step 1: create object sets and samplers
+            for( const auto& set_desc : set_descriptors_ )
             {
-                const auto set_id =
-                    this->object_sets_.add_set( descriptor.name );
+                const auto set_id = this->object_sets_.add_set( set_desc.name );
+                name_to_uuid[set_desc.name] = set_id;
 
                 this->set_samplers_.push_back(
                     std::make_unique< geode::UniformPointSetSampler< 2 > >(
                         box_ ) );
 
-                // density energy term
+                geode::add_birth_death_change_moves( this->set_samplers_.back(),
+                    *proposal_kernel, set_id, set_desc.birth_ratio,
+                    set_desc.death_ratio, set_desc.change_ratio );
+            }
+
+            // Step 2: create energy terms
+            for( const auto& energy_desc : density_descriptors_ )
+            {
+                const auto set_id = name_to_uuid.at( energy_desc.set_name );
+
                 this->ordered_energy_terms_.push_back(
                     this->energy_terms_collection_.add_energy_term(
                         std::make_unique<
                             geode::DensityTerm< geode::Point2D > >(
-                            absl::StrCat( descriptor.name, "_density" ),
-                            descriptor.density,
+                            absl::StrCat( energy_desc.set_name, "_density" ),
+                            energy_desc.density,
                             absl::flat_hash_set< geode::uuid >{ set_id } ) ) );
 
-                // target statistics
                 this->ordered_target_statistics_.push_back(
-                    descriptor.target_count );
-
-                // proposal moves
-                geode::add_birth_death_change_moves( this->set_samplers_.back(),
-                    *proposal_kernel, set_id, descriptor.birth_ratio,
-                    descriptor.death_ratio, descriptor.change_ratio );
+                    energy_desc.target_count );
             }
+
             this->mh_sampler_ =
                 std::make_unique< geode::MetropolisHastings< geode::Point2D > >(
                     this->energy_terms_collection_,
@@ -102,33 +121,37 @@ namespace
 
                 const auto expected_means =
                     this->ordered_target_statistics_[stat_id];
+
                 const auto target_vs_mean_error =
                     std::abs(
                         statistic_monitoring.means[stat_id] - expected_means )
                     / expected_means;
+
                 OPENGEODE_EXCEPTION( target_vs_mean_error < 0.02,
                     "[MH test] statistic value ",
                     statistic_monitoring.means[stat_id],
                     " for energy term: ", term.name().data(),
-                    " not close to enought to expected value ", expected_means,
+                    " not close enough to expected value ", expected_means,
                     " --> error : ", target_vs_mean_error );
 
                 const auto target_vs_variance_error =
                     std::abs( statistic_monitoring.variances[stat_id]
                               - expected_means )
                     / expected_means;
+
                 OPENGEODE_EXCEPTION( target_vs_variance_error < 0.15,
-                    "[MH test] The variance of the staistic value ",
+                    "[MH test] variance of statistic ",
                     statistic_monitoring.variances[stat_id],
                     " for energy term: ", term.name().data(),
-                    " is not close to enought to expected value ",
-                    expected_means, " --> error : ", target_vs_variance_error );
+                    " not close enough to expected value ", expected_means,
+                    " --> error : ", target_vs_variance_error );
             }
         }
 
     private:
         geode::BoundingBox2D box_;
-        std::vector< PoissonSetDescription > set_descriptors_;
+        std::vector< SetDescription > set_descriptors_;
+        std::vector< PoissonDensityDescription > density_descriptors_;
     };
 
     void test_single_type_poisson()
@@ -144,27 +167,36 @@ namespace
 
         std::array< double, 4 > birth_ratio{ 0.1, 0.5, 2., 4. };
         std::array< double, 4 > change_ratio{ 0., 1., 1., 0. };
+
         for( const auto config : geode::Range{ birth_ratio.size() } )
         {
-            PoissonSetDescription setA;
+            // --- Set description
+            SetDescription setA;
             setA.name = "A";
-            setA.density = 0.3;
-            setA.target_count = 30.0;
             setA.birth_ratio = birth_ratio[config];
             setA.death_ratio = 1.0;
             setA.change_ratio = change_ratio[config];
 
+            // --- Energy term description
+            PoissonDensityDescription densityA;
+            densityA.set_name = "A";
+            densityA.density = 0.3;
+            densityA.target_count = 30.0;
+
             PoissonSimulationRunner runner( box );
             runner.add_set_descriptor( setA );
+            runner.add_density_descriptor( densityA );
             runner.initialize();
 
             constexpr geode::index_t steps = 1000;
             constexpr geode::index_t nb_realizations = 1000;
+
             runner.run( engine, steps );
             auto statistic_monitoring = runner.run_print_and_monitor(
                 "poisson_statistics", engine, steps, nb_realizations );
             runner.check_statistics( statistic_monitoring );
         }
+
         geode::Logger::info( "--> SUCCESS!" );
     }
 
@@ -179,48 +211,34 @@ namespace
         box.add_point( geode::Point2D{ { 0.0, 0.0 } } );
         box.add_point( geode::Point2D{ { 10.0, 10.0 } } );
 
-        std::array< double, 4 > birth_ratio{ 0.1, 0.5, 2., 4. };
-        std::array< double, 4 > change_ratio{ 0., 1., 1., 0. };
-        for( const auto config : geode::Range{ birth_ratio.size() } )
-        {
-            PoissonSetDescription set01;
-            set01.name = "set01";
-            set01.density = 0.1;
-            set01.target_count = 10.0;
-            set01.birth_ratio = 1.;
-            set01.death_ratio = 3.;
-            set01.change_ratio = 1.;
-
-            PoissonSetDescription set02;
-            set02.name = "set02";
-            set02.density = 0.4;
-            set02.target_count = 40.0;
-            set02.birth_ratio = 3.;
-            set02.death_ratio = 0.5;
-            set02.change_ratio = 1.;
-
-            PoissonSetDescription set03;
-            set03.name = "set03";
-            set03.density = 0.3;
-            set03.target_count = 30.0;
-            set03.birth_ratio = 4.;
-            set03.death_ratio = 1.;
-            set03.change_ratio = 1.;
+        // --- Set descriptions
+        SetDescription set01{ "set01", 2.0, 3.0, 1.0 };
+        SetDescription set02{ "set02", 3.0, 0.5, 1.0 };
+        SetDescription set03{ "set03", 4.0, 1.0, 1.0 };
 
-            PoissonSimulationRunner runner( box );
-            runner.add_set_descriptor( set01 );
-            runner.add_set_descriptor( set02 );
-            runner.add_set_descriptor( set03 );
+        // --- Energy term descriptions
+        PoissonDensityDescription density01{ "set01", 0.1, 10.0 };
+        PoissonDensityDescription density02{ "set02", 0.4, 40.0 };
+        PoissonDensityDescription density03{ "set03", 0.3, 30.0 };
 
-            runner.initialize();
+        PoissonSimulationRunner runner( box );
+        runner.add_set_descriptor( set01 );
+        runner.add_set_descriptor( set02 );
+        runner.add_set_descriptor( set03 );
 
-            constexpr geode::index_t steps = 1000;
-            constexpr geode::index_t nb_realizations = 1000;
+        runner.add_density_descriptor( density01 );
+        runner.add_density_descriptor( density02 );
+        runner.add_density_descriptor( density03 );
+
+        runner.initialize();
+
+        constexpr geode::index_t steps = 1000;
+        constexpr geode::index_t nb_realizations = 1000;
+
+        auto statistic_monitoring = runner.run_print_and_monitor(
+            "poisson_statistics", engine, steps, nb_realizations );
+        runner.check_statistics( statistic_monitoring );
 
-            auto statistic_monitoring = runner.run_print_and_monitor(
-                "poisson_statistics", engine, steps, nb_realizations );
-            runner.check_statistics( statistic_monitoring );
-        }
         geode::Logger::info( "--> SUCCESS!" );
     }
 } // namespace