feat(SimulationHelpers): add STatistic monitoring

Author: francoisbonneau

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

@@ -22,13 +22,8 @@
  */
 
 #pragma once
+#include <geode/basic/range.hpp>
 #include <geode/stochastic/common.hpp>
-// #include <geode/stochastic/sampling/mcmc/metropolis_hasting_sampler.hpp>
-// #include <geode/stochastic/sampling/mcmc/models/energy_term_collection.hpp>
-
-#include <absl/strings/str_join.h>
-#include <filesystem>
-#include <fstream>
 
 namespace geode
 {
@@ -44,40 +39,50 @@ namespace geode
 
         StatisticsMonitor( const index_t nb_energy_terms )
         {
-            sum.resize( nb_energy_terms, 0.0 );
-            sum_squares.resize( nb_energy_terms, 0.0 );
-            means.resize( nb_energy_terms, 0.0 );
-            variances.resize( nb_energy_terms, 0.0 );
+            sum_.resize( nb_energy_terms, 0.0 );
+            sum_squares_.resize( nb_energy_terms, 0.0 );
+            means_.resize( nb_energy_terms, 0.0 );
+            variances_.resize( nb_energy_terms, 0.0 );
         }
 
         void add_realization( const std::vector< double >& values )
         {
-            for( const auto stat_id : Range{ values.size() } )
+            OPENGEODE_EXCEPTION( values.size() == sum_.size(),
+                "[StatisticsMonitor] - Mismatch between realization size and "
+                "expected number of statistics." );
+            ++count_;
+            for( size_t i = 0; i < values.size(); ++i )
             {
-                sum[stat_id] += values[stat_id];
-                sum_squares[stat_id] += values[stat_id] * values[stat_id];
+                double delta = values[i] - means_[i];
+                means_[i] += delta / count_;
+                if( count_ > 1 )
+                    variances_[i] = ( ( count_ - 2 ) * variances_[i]
+                                        + delta * ( values[i] - means_[i] ) )
+                                    / ( count_ - 1 );
             }
         }
 
-        void finalize( const index_t nb_realizations )
+        const index_t statiscal_count() const
         {
-            for( const auto stat_id : Range{ sum.size() } )
-            {
-                means[stat_id] = sum[stat_id] / nb_realizations;
-                double variance =
-                    ( sum_squares[stat_id]
-                        - ( sum[stat_id] * sum[stat_id] ) / nb_realizations )
-                    / ( nb_realizations - 1 );
-                variances[stat_id] = variance;
-                // stddevs[stat_id] =std::sqrt( std::max( variance, 0.0 ) );
-            }
+            return count_;
         }
 
-    public:
-        std::vector< double > sum;
-        std::vector< double > sum_squares;
-        std::vector< double > means;
-        std::vector< double > variances;
+        const std::vector< double >& means() const
+        {
+            return means_;
+        }
+
+        const std::vector< double >& variances() const
+        {
+            return variances_;
+        }
+
+    private:
+        std::vector< double > sum_;
+        std::vector< double > sum_squares_;
+        std::vector< double > means_;
+        std::vector< double > variances_;
+        index_t count_{ 0 };
     };
 
 } // namespace geode

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

@@ -22,9 +22,9 @@
  */
 
 #pragma once
+
 #include <geode/stochastic/common.hpp>
-// #include <geode/stochastic/sampling/mcmc/metropolis_hasting_sampler.hpp>
-// #include <geode/stochastic/sampling/mcmc/models/energy_term_collection.hpp>
+
 #include <geode/stochastic/sampling/mcmc/helpers/simulation_monitor.hpp>
 #include <geode/stochastic/spatial/object_sets.hpp>
 
@@ -46,7 +46,7 @@ namespace geode
         std::string realisations_prefix{ "pattern_" };
         index_t realisations_print_frequency{ 100 };
 
-        std::string output_folder{ "." };
+        std::string output_folder{ std::filesystem::current_path().string() };
     };
 
     class SimulationPrinter
@@ -59,7 +59,7 @@ namespace geode
 
         // Print statistics to the configured statistics file
         void print_statistics(
-            const std::vector< double >& stats, absl::string_view header )
+            const std::vector< double >& stats, absl::string_view header ) const
         {
             if( !config_.print_statistics )
                 return;
@@ -72,9 +72,10 @@ namespace geode
 
         template < typename ObjectType >
         void print_object_sets( const ObjectSets< ObjectType >& object_sets,
-            index_t realization_id )
+            index_t realization_id ) const
         {
-            if( !config_.print_realisations )
+            if( !config_.print_realisations
+                || realization_id % config_.realisations_print_frequency != 0 )
                 return;
 
             const auto filename =
@@ -96,7 +97,7 @@ namespace geode
             }
         }
         void print_statistics_summary( const StatisticsMonitor& monitor,
-            absl::string_view energy_term_names )
+            absl::string_view energy_term_names ) const
         {
             if( !config_.print_statistics_summary )
                 return;
@@ -109,16 +110,17 @@ namespace geode
             std::ofstream file = open_file_with_dirs( summary_path );
             file << "# Summary statistics\n";
             file << energy_term_names.data() << "\n";
-            file << absl::StrCat( "# Count:\n", monitor.count, "\n" );
             file << absl::StrCat(
-                "# Means:\n", absl::StrJoin( monitor.means, " ; " ), "\n" );
+                "# Count:\n", monitor.statiscal_count(), "\n" );
+            file << absl::StrCat(
+                "# Means:\n", absl::StrJoin( monitor.means(), " ; " ), "\n" );
             file << absl::StrCat( "# Variances:\n",
-                absl::StrJoin( monitor.variances, " ; " ), "\n" );
+                absl::StrJoin( monitor.variances(), " ; " ), "\n" );
         }
 
     private:
         void write_header_if_new(
-            absl::string_view filename, absl::string_view header )
+            absl::string_view filename, absl::string_view header ) const
         {
             namespace fs = std::filesystem;
             fs::path file_path{ std::string( filename ) };
@@ -148,7 +150,8 @@ namespace geode
 
             return file;
         }
-        const std::string& create_statistics_file( absl::string_view header )
+        const std::string& create_statistics_file(
+            absl::string_view header ) const
         {
             stats_file_path_ = ( std::filesystem::path( config_.output_folder )
                                  / config_.statistics_filename )
@@ -165,7 +168,7 @@ namespace geode
 
     private:
         SimulationPrinterConfigurator config_;
-        std::optional< std::string > stats_file_path_;
+        mutable std::optional< std::string > stats_file_path_;
     };
 
 } // namespace geode