1+ /*
2+ * Copyright (c) 2019 - 2025 Geode-solutions
3+ *
4+ * Permission is hereby granted, free of charge, to any person obtaining a copy
5+ * of this software and associated documentation files (the "Software"), to deal
6+ * in the Software without restriction, including without limitation the rights
7+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
8+ * copies of the Software, and to permit persons to whom the Software is
9+ * furnished to do so, subject to the following conditions:
10+ *
11+ * The above copyright notice and this permission notice shall be included in
12+ * all copies or substantial portions of the Software.
13+ *
14+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
15+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
16+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
17+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
18+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
19+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
20+ * SOFTWARE.
21+ *
22+ */
23+
24+ #pragma once
25+ #include < geode/stochastic/sampling/mcmc/metropolis_hasting_sampler.hpp>
26+ #include < geode/stochastic/sampling/mcmc/models/energy_term_collection.hpp>
27+
28+ #include " absl/strings/str_join.h"
29+ #include < fstream>
30+
31+ namespace geode
32+ {
33+ class opengeode_stochastic_stochastic_api MonitoringStatistics
34+ {
35+ public:
36+ std::vector< double > sum;
37+ std::vector< double > sum_squares;
38+ std::vector< double > means;
39+ std::vector< double > variances;
40+
41+ MonitoringStatistics ( const index_t nb_energy_terms )
42+ {
43+ sum.resize ( nb_energy_terms, 0.0 );
44+ sum_squares.resize ( nb_energy_terms, 0.0 );
45+ means.resize ( nb_energy_terms, 0.0 );
46+ variances.resize ( nb_energy_terms, 0.0 );
47+ }
48+
49+ void add_realization ( const std::vector< double >& values )
50+ {
51+ for ( const auto stat_id : Range{ values.size () } )
52+ {
53+ sum[stat_id] += values[stat_id];
54+ sum_squares[stat_id] += values[stat_id] * values[stat_id];
55+ }
56+ }
57+
58+ void finalize ( const index_t nb_realizations )
59+ {
60+ for ( const auto stat_id : Range{ sum.size () } )
61+ {
62+ means[stat_id] = sum[stat_id] / nb_realizations;
63+ double variance =
64+ ( sum_squares[stat_id]
65+ - ( sum[stat_id] * sum[stat_id] ) / nb_realizations )
66+ / ( nb_realizations - 1 );
67+ variances[stat_id] = variance;
68+ // stddevs[stat_id] =std::sqrt( std::max( variance, 0.0 ) );
69+ }
70+ }
71+ };
72+
73+ template < typename ObjectType >
74+ class SimulationRunner
75+ {
76+ public:
77+ SimulationRunner () = default ;
78+ virtual ~SimulationRunner () = default ;
79+
80+ virtual void initialize () = 0;
81+
82+ const ObjectSets< ObjectType >& run (
83+ RandomEngine& engine, const index_t steps )
84+ {
85+ mh_sampler_->walk ( object_sets_, engine, steps );
86+ return object_sets_;
87+ }
88+
89+ void run_and_print ( std::string_view filename,
90+ RandomEngine& engine,
91+ const index_t steps,
92+ const index_t nb_realizations )
93+ {
94+ const auto file_exist =
95+ static_cast < bool >( std::ifstream ( filename.data () ) );
96+ if ( !file_exist )
97+ {
98+ const auto header = statistics_header_file ();
99+ print_to_file ( filename, header );
100+ }
101+
102+ for ( const auto realization : Range{ nb_realizations } )
103+ {
104+ run ( engine, steps );
105+ const auto statistics = statistics_string ();
106+ print_to_file ( filename, statistics );
107+ }
108+ }
109+
110+ MonitoringStatistics run_print_and_monitor ( std::string_view filename,
111+ RandomEngine& engine,
112+ const index_t steps,
113+ const index_t nb_realizations )
114+ {
115+ const auto file_exist =
116+ static_cast < bool >( std::ifstream ( filename.data () ) );
117+ if ( !file_exist )
118+ {
119+ const auto header = statistics_header_file ();
120+ print_to_file ( filename, header );
121+ }
122+ MonitoringStatistics stat_monitoring (
123+ energy_terms_collection_.size () );
124+
125+ for ( const auto realization : Range{ nb_realizations } )
126+ {
127+ run ( engine, steps );
128+ const auto stats = statistics ();
129+ print_to_file ( filename,
130+ absl::StrCat ( absl::StrJoin ( stats, " ; " " \n "
131+ stat_monitoring.add_realization ( stats );
132+ }
133+ stat_monitoring.finalize ( nb_realizations );
134+ return stat_monitoring;
135+ }
136+
137+ const ObjectSets< ObjectType >& current_pattern_realization () const
138+ {
139+ return object_sets_;
140+ }
141+
142+ std::vector< double > statistics () const
143+ {
144+ std::vector< double > statistic_values;
145+ statistic_values.reserve ( ordered_energy_terms_.size () );
146+
147+ for ( const auto & energy_term_uuid : ordered_energy_terms_ )
148+ {
149+ const auto & term =
150+ energy_terms_collection_.get ( energy_term_uuid );
151+ statistic_values.push_back ( term.statistic ( object_sets_ ) );
152+ }
153+
154+ return statistic_values;
155+ }
156+
157+ std::string statistics_log_info () const
158+ {
159+ std::string message ( " Pattern statistics: "
160+ for ( const auto term_id :
161+ geode::Range{ ordered_energy_terms_.size () } )
162+ {
163+ const auto & energy_term = energy_terms_collection_.get (
164+ ordered_energy_terms_[term_id] );
165+ const double value = energy_term.statistic ( object_sets_ );
166+ absl::StrAppend ( &message, " \t Term(" name (),
167+ " ) --> value/traget: " " / "
168+ ordered_target_statistics_[term_id] );
169+ }
170+ return message;
171+ }
172+
173+ protected:
174+ std::string energy_term_names () const
175+ {
176+ std::vector< std::string > term_names;
177+ term_names.reserve ( ordered_energy_terms_.size () );
178+
179+ for ( const auto & energy_term_uuid : ordered_energy_terms_ )
180+ {
181+ const auto & term =
182+ energy_terms_collection_.get ( energy_term_uuid );
183+ term_names.push_back ( term.name ().data () );
184+ }
185+
186+ return absl::StrCat ( absl::StrJoin ( term_names, " ; " " \n "
187+ }
188+
189+ std::string statistics_string () const
190+ {
191+ return absl::StrCat ( absl::StrJoin ( statistics (), " ; " " \n "
192+ }
193+
194+ std::string statistics_header_file ()
195+ {
196+ std::string message ( " Sufficient statistics mcmc iterations:\n "
197+ absl::StrAppend ( &message, energy_term_names () );
198+ return message;
199+ }
200+
201+ void print_to_file (
202+ absl::string_view filename, absl::string_view message )
203+ {
204+ std::ofstream file (
205+ filename.data (), std::ofstream::out | std::ofstream::app );
206+ file << message;
207+ file.close ();
208+ return ;
209+ }
210+
211+ protected:
212+ std::vector< std::unique_ptr< geode::ObjectSetSampler< ObjectType > > >
213+ set_samplers_;
214+
215+ std::vector< geode::uuid > ordered_energy_terms_;
216+ std::vector< double > ordered_target_statistics_;
217+
218+ EnergyTermCollection< ObjectType > energy_terms_collection_;
219+ std::unique_ptr< geode::MetropolisHastings< ObjectType > > mh_sampler_;
220+
221+ ObjectSets< ObjectType > object_sets_;
222+ };
223+ } // namespace geode
0 commit comments