diff --git a/include/geode/stochastic/sampling/direct/double_sampler.hpp b/include/geode/stochastic/sampling/direct/double_sampler.hpp index 6bd25bc..ac962f2 100644 --- a/include/geode/stochastic/sampling/direct/double_sampler.hpp +++ b/include/geode/stochastic/sampling/direct/double_sampler.hpp @@ -43,6 +43,20 @@ namespace geode Gaussian, TruncatedGaussian >; + struct DistributionDescription + { + std::string name{ "default_distribution" }; + DistributionType distribution_type; + + std::optional< double > min_value; + std::optional< double > max_value; + std::optional< double > mean; + std::optional< double > standard_deviation; + }; + + static Distribution create_distribution( + const DistributionDescription& desc ); + static double sample( RandomEngine& engine, const Distribution& dist ); }; diff --git a/include/geode/stochastic/sampling/direct/object_set_sampler/point_set_sampler.hpp b/include/geode/stochastic/sampling/direct/object_set_sampler/point_set_sampler.hpp index 32d0364..391e332 100644 --- a/include/geode/stochastic/sampling/direct/object_set_sampler/point_set_sampler.hpp +++ b/include/geode/stochastic/sampling/direct/object_set_sampler/point_set_sampler.hpp @@ -81,7 +81,7 @@ namespace geode } private: - const BoundingBox< dimension >& box_; + BoundingBox< dimension > box_; }; } // namespace geode \ No newline at end of file diff --git a/include/geode/stochastic/sampling/direct/object_set_sampler/segment_set_sampler.hpp b/include/geode/stochastic/sampling/direct/object_set_sampler/segment_set_sampler.hpp new file mode 100644 index 0000000..115b306 --- /dev/null +++ b/include/geode/stochastic/sampling/direct/object_set_sampler/segment_set_sampler.hpp @@ -0,0 +1,103 @@ +/* + * Copyright (c) 2019 - 2025 Geode-solutions + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + * + */ + +#pragma once + +#include +#include + +#include +#include +#include + +namespace geode +{ + class UniformSegmentSetSampler : public ObjectSetSampler< OwnerSegment2D > + { + public: + UniformSegmentSetSampler( const BoundingBox< 2 >& box, + const DoubleSampler::Distribution& length, + const DoubleSampler::Distribution& azimuth ) + : ObjectSetSampler< OwnerSegment2D >{}, + box_{ box }, + length_{ length }, + azimuth_{ azimuth } + { + auto volume = box_.n_volume(); + OPENGEODE_EXCEPTION( volume != 0., + "[SegmentSetSampler] - Undefined Bounding Box (volume == 0)." ); + this->log_pdf_ = -std::log( volume ); + } + + OwnerSegment2D sample( RandomEngine& engine ) const override + { + auto seg = SegmentUniformSampler::sample( + engine, box_, length_, azimuth_ ); + return seg; + } + + OwnerSegment2D change( + const OwnerSegment2D& obj, RandomEngine& engine ) const override + { + double ratio = 0.1; + const auto& extremities = obj.vertices(); + const auto current = + static_cast< local_index_t >( engine.sample_bernoulli( 0.5 ) ); + + geode::Sphere< 2 > ball{ extremities[current], + ratio * obj.length() }; + + auto new_point = PointUniformSampler::sample< 2 >( engine, ball ); + constexpr index_t max_try{ 100 }; + for( const auto try_id : geode::Range{ max_try } ) + { + if( box_.contains( new_point ) ) + { + OwnerSegment2D new_segment{ obj }; + new_segment.set_point( current, new_point ); + return new_segment; + } + new_point = PointUniformSampler::sample< 2 >( engine, ball ); + } + throw OpenGeodeException( absl::StrCat( + "[PointSampler] - Cannot find a point in the box: ", + box_.string() ) ); + return obj; + } + + private: + bool is_valid_object( const OwnerSegment2D& obj ) const override + { + const auto& extremities = obj.vertices(); + + return box_.contains( extremities[0] ) + && box_.contains( extremities[1] ); + } + + private: + BoundingBox2D box_; + DoubleSampler::Distribution length_; + DoubleSampler::Distribution azimuth_; + }; + +} // namespace geode \ No newline at end of file diff --git a/include/geode/stochastic/sampling/distributions.hpp b/include/geode/stochastic/sampling/distributions.hpp index 3818eda..d100911 100644 --- a/include/geode/stochastic/sampling/distributions.hpp +++ b/include/geode/stochastic/sampling/distributions.hpp @@ -27,8 +27,16 @@ #include +#include + namespace geode { + struct DistributionTag + { + }; + + using DistributionType = NamedType< std::string, DistributionTag >; + template < typename Type > struct UniformClosed { @@ -37,6 +45,16 @@ namespace geode Type min_value{ static_cast< Type >( 0 ) }; Type max_value{ static_cast< Type >( 1 ) }; + + [[nodiscard]] static DistributionType distribution_type_static() + { + return DistributionType{ "UniformClosed" }; + } + + [[nodiscard]] DistributionType distribution_type() const + { + return distribution_type_static(); + } }; template < typename Type > @@ -47,6 +65,16 @@ namespace geode Type min_value{ static_cast< Type >( 0 ) }; Type max_value{ static_cast< Type >( 1 ) }; + + [[nodiscard]] static DistributionType distribution_type_static() + { + return DistributionType{ "UniformClosedOpen" }; + } + + [[nodiscard]] DistributionType distribution_type() const + { + return distribution_type_static(); + } }; struct opengeode_stochastic_stochastic_api Gaussian @@ -56,6 +84,16 @@ namespace geode double mean{ 0. }; double standard_deviation{ 1. }; + + [[nodiscard]] static DistributionType distribution_type_static() + { + return DistributionType{ "Gaussian" }; + } + + [[nodiscard]] DistributionType distribution_type() const + { + return distribution_type_static(); + } }; struct opengeode_stochastic_stochastic_api TruncatedGaussian @@ -69,6 +107,16 @@ namespace geode std::optional< double > min_value; std::optional< double > max_value; + + [[nodiscard]] static DistributionType distribution_type_static() + { + return DistributionType{ "TruncatedGaussian" }; + } + + [[nodiscard]] DistributionType distribution_type() const + { + return distribution_type_static(); + } }; } // namespace geode \ No newline at end of file diff --git a/include/geode/stochastic/sampling/mcmc/helpers/simulation_manager.hpp b/include/geode/stochastic/sampling/mcmc/helpers/simulation_manager.hpp new file mode 100644 index 0000000..e69de29 diff --git a/include/geode/stochastic/sampling/mcmc/helpers/simulation_monitor.hpp b/include/geode/stochastic/sampling/mcmc/helpers/simulation_monitor.hpp new file mode 100644 index 0000000..57923eb --- /dev/null +++ b/include/geode/stochastic/sampling/mcmc/helpers/simulation_monitor.hpp @@ -0,0 +1,88 @@ +/* + * Copyright (c) 2019 - 2025 Geode-solutions + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + * + */ + +#pragma once +#include +#include + +namespace geode +{ + + class StatisticsMonitor + { + public: + StatisticsMonitor( StatisticsMonitor&& ) = default; + StatisticsMonitor( const StatisticsMonitor& ) = default; + StatisticsMonitor& operator=( StatisticsMonitor&& ) noexcept = default; + StatisticsMonitor& operator=( + const StatisticsMonitor& ) noexcept = default; + + 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 ); + } + + void add_realization( const std::vector< double >& values ) + { + 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 ) + { + 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 ); + } + } + + const index_t statiscal_count() const + { + return count_; + } + + 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 \ No newline at end of file diff --git a/include/geode/stochastic/sampling/mcmc/helpers/simulation_printer.hpp b/include/geode/stochastic/sampling/mcmc/helpers/simulation_printer.hpp new file mode 100644 index 0000000..93dc7fd --- /dev/null +++ b/include/geode/stochastic/sampling/mcmc/helpers/simulation_printer.hpp @@ -0,0 +1,174 @@ +/* + * Copyright (c) 2019 - 2025 Geode-solutions + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + * + */ + +#pragma once + +#include + +#include +#include + +#include +#include +#include + +namespace geode +{ + struct SimulationPrinterConfigurator + { + bool print_statistics{ true }; + std::string statistics_filename{ "statistics.txt" }; + + bool print_statistics_summary{ true }; + std::string statistics_summary_filename{ "statistics_summary.txt" }; + + bool print_realisations{ true }; + std::string realisations_prefix{ "pattern_" }; + index_t realisations_print_frequency{ 100 }; + + std::string output_folder{ std::filesystem::current_path().string() }; + }; + + class SimulationPrinter + { + public: + SimulationPrinter( const SimulationPrinterConfigurator& config ) + : config_( config ) + { + } + + // Print statistics to the configured statistics file + void print_statistics( + const std::vector< double >& stats, absl::string_view header ) const + { + if( !config_.print_statistics ) + return; + const auto stats_file_path = + stats_file_path_.value_or( create_statistics_file( header ) ); + std::ofstream file = + open_file_with_dirs( stats_file_path, std::ios::app ); + file << absl::StrJoin( stats, " ; " ) << "\n"; + } + + template < typename ObjectType > + void print_object_sets( const ObjectSets< ObjectType >& object_sets, + index_t realization_id ) const + { + if( !config_.print_realisations + || realization_id % config_.realisations_print_frequency != 0 ) + return; + + const auto filename = + ( std::filesystem::path( config_.output_folder ) + / absl::StrCat( + config_.realisations_prefix, realization_id, ".txt" ) ) + .string(); + + std::ofstream file = open_file_with_dirs( filename ); + + const auto all_objects = object_sets.get_all_object(); + file << "#nb_objects\t" << all_objects.size() << "\n"; + + for( const auto& object_id : all_objects ) + { + const auto& object = object_sets.get_object( object_id ); + file << object.string() << "\t" << object_id.set_id.string() + << "\n"; + } + } + void print_statistics_summary( const StatisticsMonitor& monitor, + absl::string_view energy_term_names ) const + { + if( !config_.print_statistics_summary ) + return; + + const auto summary_path = + ( std::filesystem::path( config_.output_folder ) + / config_.statistics_summary_filename ) + .string(); + + 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.statiscal_count(), "\n" ); + file << absl::StrCat( + "# Means:\n", absl::StrJoin( monitor.means(), " ; " ), "\n" ); + file << absl::StrCat( "# Variances:\n", + absl::StrJoin( monitor.variances(), " ; " ), "\n" ); + } + + private: + void write_header_if_new( + absl::string_view filename, absl::string_view header ) const + { + namespace fs = std::filesystem; + fs::path file_path{ std::string( filename ) }; + if( !fs::exists( file_path ) ) + { + std::ofstream file = open_file_with_dirs( filename ); + file << header; + } + } + + std::ofstream open_file_with_dirs( absl::string_view path_filename, + std::ios::openmode mode = std::ofstream::out ) const + { + namespace fs = std::filesystem; + fs::path file_path{ std::string( path_filename ) }; + + if( !file_path.has_parent_path() ) + file_path = fs::current_path() / file_path; + + if( file_path.has_parent_path() ) + fs::create_directories( file_path.parent_path() ); + + std::ofstream file( file_path, mode ); + if( !file.is_open() ) + throw geode::OpenGeodeException( + "Cannot open file: " + file_path.string() ); + + return file; + } + const std::string& create_statistics_file( + absl::string_view header ) const + { + stats_file_path_ = ( std::filesystem::path( config_.output_folder ) + / config_.statistics_filename ) + .string(); + + if( config_.print_statistics ) + { + write_header_if_new( *stats_file_path_, + absl::StrCat( + "# Simulation Statistics\n", header.data(), "\n" ) ); + } + return *stats_file_path_; + } + + private: + SimulationPrinterConfigurator config_; + mutable std::optional< std::string > stats_file_path_; + }; + +} // namespace geode \ No newline at end of file diff --git a/include/geode/stochastic/sampling/mcmc/helpers/simulation_runner.hpp b/include/geode/stochastic/sampling/mcmc/helpers/simulation_runner.hpp index 4b9d338..a46de63 100644 --- a/include/geode/stochastic/sampling/mcmc/helpers/simulation_runner.hpp +++ b/include/geode/stochastic/sampling/mcmc/helpers/simulation_runner.hpp @@ -23,60 +23,22 @@ #pragma once #include +#include +#include #include #include #include -#include namespace geode { - class MonitoringStatistics + struct SimulationConfigurator { - public: - MonitoringStatistics( MonitoringStatistics&& ) = default; - MonitoringStatistics( const MonitoringStatistics& ) = default; - MonitoringStatistics& operator=( - MonitoringStatistics&& ) noexcept = default; - MonitoringStatistics& operator=( - const MonitoringStatistics& ) noexcept = default; - - MonitoringStatistics( 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 ); - } + index_t realizations{ 1000 }; + index_t metropolis_hasting_steps{ 1000 }; + index_t burn_in_steps{ 1000 }; - void add_realization( const std::vector< double >& values ) - { - for( const auto stat_id : Range{ values.size() } ) - { - sum[stat_id] += values[stat_id]; - sum_squares[stat_id] += values[stat_id] * values[stat_id]; - } - } - - void finalize( const index_t nb_realizations ) - { - 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 ) ); - } - } - - public: - std::vector< double > sum; - std::vector< double > sum_squares; - std::vector< double > means; - std::vector< double > variances; + std::optional< SimulationPrinter > printer{ std::nullopt }; }; template < typename ObjectType > @@ -95,60 +57,63 @@ namespace geode return object_sets_; } - void run_and_print( std::string_view filename, - RandomEngine& engine, - const index_t steps, - const index_t nb_realizations ) + void run( + RandomEngine& engine, const SimulationConfigurator& configurator ) { - const auto file_exist = - static_cast< bool >( std::ifstream( filename.data() ) ); - if( !file_exist ) + if( configurator.burn_in_steps > 0 ) { - const auto header = statistics_header_file(); - print_to_file( filename, header ); + run( engine, configurator.burn_in_steps ); } - - for( const auto realization : Range{ nb_realizations } ) + for( const auto realization : Range{ configurator.realizations } ) { - run( engine, steps ); - const auto statistics = statistics_string(); - print_to_file( filename, statistics ); + run( engine, configurator.metropolis_hasting_steps ); + + if( configurator.printer.has_value() ) + { + configurator.printer->print_statistics( + state_statistics(), model_energy_term_names() ); + configurator.printer->print_object_sets( + object_sets_, realization ); + } } } - MonitoringStatistics run_print_and_monitor( std::string_view filename, - RandomEngine& engine, - const index_t steps, - const index_t nb_realizations ) + StatisticsMonitor run_and_monitor( + RandomEngine& engine, const SimulationConfigurator& configurator ) { - const auto file_exist = - static_cast< bool >( std::ifstream( filename.data() ) ); - if( !file_exist ) + if( configurator.burn_in_steps > 0 ) { - const auto header = statistics_header_file(); - print_to_file( filename, header ); + run( engine, configurator.burn_in_steps ); } - MonitoringStatistics stat_monitoring( + StatisticsMonitor stat_monitoring( energy_terms_collection_.size() ); - - for( const auto realization : Range{ nb_realizations } ) + for( const auto realization : Range{ configurator.realizations } ) { - run( engine, steps ); - const auto stats = statistics(); - print_to_file( filename, - absl::StrCat( absl::StrJoin( stats, " ; " ), "\n" ) ); + run( engine, configurator.metropolis_hasting_steps ); + const auto stats = state_statistics(); stat_monitoring.add_realization( stats ); + if( configurator.printer.has_value() ) + { + configurator.printer->print_statistics( + state_statistics(), model_energy_term_names() ); + configurator.printer->print_object_sets( + object_sets_, realization ); + } + } + if( configurator.printer.has_value() ) + { + configurator.printer->print_statistics_summary( + stat_monitoring, model_energy_term_names() ); } - stat_monitoring.finalize( nb_realizations ); return stat_monitoring; } - const ObjectSets< ObjectType >& current_pattern_realization() const + const ObjectSets< ObjectType >& state_realization() const { return object_sets_; } - std::vector< double > statistics() const + std::vector< double > state_statistics() const { std::vector< double > statistic_values; statistic_values.reserve( ordered_energy_terms_.size() ); @@ -163,24 +128,7 @@ namespace geode return statistic_values; } - std::string statistics_log_info() const - { - std::string message( "Pattern statistics: " ); - for( const auto term_id : - geode::Range{ ordered_energy_terms_.size() } ) - { - const auto& energy_term = energy_terms_collection_.get( - ordered_energy_terms_[term_id] ); - const double value = energy_term.statistic( object_sets_ ); - absl::StrAppend( &message, " \t Term(", energy_term.name(), - ") --> value/traget: ", value, " / ", - ordered_target_statistics_[term_id] ); - } - return message; - } - - protected: - std::string energy_term_names() const + std::string model_energy_term_names() const { std::vector< std::string > term_names; term_names.reserve( ordered_energy_terms_.size() ); @@ -195,28 +143,6 @@ namespace geode return absl::StrCat( absl::StrJoin( term_names, " ; " ), "\n" ); } - std::string statistics_string() const - { - return absl::StrCat( absl::StrJoin( statistics(), " ; " ), "\n" ); - } - - std::string statistics_header_file() - { - std::string message( "Sufficient statistics mcmc iterations:\n" ); - absl::StrAppend( &message, energy_term_names() ); - return message; - } - - void print_to_file( - absl::string_view filename, absl::string_view message ) - { - std::ofstream file( - filename.data(), std::ofstream::out | std::ofstream::app ); - file << message; - file.close(); - return; - } - protected: std::vector< std::unique_ptr< geode::ObjectSetSampler< ObjectType > > > set_samplers_; diff --git a/include/geode/stochastic/sampling/mcmc/metropolis_hasting_sampler.hpp b/include/geode/stochastic/sampling/mcmc/metropolis_hasting_sampler.hpp index 8fa4214..933b96a 100644 --- a/include/geode/stochastic/sampling/mcmc/metropolis_hasting_sampler.hpp +++ b/include/geode/stochastic/sampling/mcmc/metropolis_hasting_sampler.hpp @@ -219,8 +219,6 @@ namespace geode const auto old_object_id = proposal.old_object_id(); const auto delta_log_energy = gibbs_energy_.delta_log_change( state, old_object_id, new_object ); - // should we test that objects are in the same group? - // should be ensured by the dynamic return accept_or_reject( proposal, state, engine, delta_log_energy, []( auto& state, auto& proposal ) { state.update_object( proposal.old_object_id(), diff --git a/include/geode/stochastic/spatial/object_sets.hpp b/include/geode/stochastic/spatial/object_sets.hpp index a43f1e9..1d1882a 100644 --- a/include/geode/stochastic/spatial/object_sets.hpp +++ b/include/geode/stochastic/spatial/object_sets.hpp @@ -26,9 +26,11 @@ #include #include + #include #include #include + #include #include diff --git a/include/geode/stochastic/spatial/pairwise_interactions.hpp b/include/geode/stochastic/spatial/pairwise_interactions.hpp index 3bda5bc..6915304 100644 --- a/include/geode/stochastic/spatial/pairwise_interactions.hpp +++ b/include/geode/stochastic/spatial/pairwise_interactions.hpp @@ -25,6 +25,7 @@ #include +#include #include #include @@ -74,38 +75,26 @@ namespace geode SCOPE scope_{ SCOPE::all_set }; }; + /*! + * EuclideanCutoffInteraction + * A pairwise interaction that returns 1 if the Euclidean distance + * between objects is within a cutoff radius, otherwise 0. + */ template < typename Type > class EuclideanCutoffInteraction : public PairwiseInteraction< Type > { public: - EuclideanCutoffInteraction( double cutoff_distance ) - : PairwiseInteraction< Type >(), cutoff_distance_( cutoff_distance ) - { - } + explicit EuclideanCutoffInteraction( double cutoff_distance ); EuclideanCutoffInteraction( double cutoff_distance, - typename PairwiseInteraction< Type >::SCOPE scope ) - : PairwiseInteraction< Type >( scope ), - cutoff_distance_( cutoff_distance ) - { - } + typename PairwiseInteraction< Type >::SCOPE scope ); - double neighborhood_searching_distance() const override - { - return cutoff_distance_; - } + double neighborhood_searching_distance() const override; protected: double compute( const ObjectRef< Type >& object_a, - const ObjectRef< Type >& object_b ) const override - { - auto dist = geode::point_point_distance( - object_barycenter( object_a.object ), - object_barycenter( object_b.object ) ); - return dist <= cutoff_distance_ ? 1.0 : 0.0; - } + const ObjectRef< Type >& object_b ) const override; private: double cutoff_distance_{ GLOBAL_EPSILON }; }; - } // namespace geode \ No newline at end of file diff --git a/src/geode/stochastic/CMakeLists.txt b/src/geode/stochastic/CMakeLists.txt index 5eb0368..ebd5d89 100644 --- a/src/geode/stochastic/CMakeLists.txt +++ b/src/geode/stochastic/CMakeLists.txt @@ -25,6 +25,7 @@ add_geode_library( "spatial/object_set.cpp" "spatial/object_sets.cpp" "spatial/object_neighborhood.cpp" + "spatial/pairwise_interactions.cpp" "sampling/direct/ball_sampler.cpp" "sampling/direct/bounding_box_sampler.cpp" "sampling/direct/double_sampler.cpp" @@ -43,12 +44,16 @@ add_geode_library( "spatial/details/RTree.hpp" "sampling/direct/object_set_sampler/object_set_sampler.hpp" "sampling/direct/object_set_sampler/point_set_sampler.hpp" + "sampling/direct/object_set_sampler/segment_set_sampler.hpp" "sampling/direct/ball_sampler.hpp" "sampling/direct/bounding_box_sampler.hpp" "sampling/direct/double_sampler.hpp" "sampling/direct/point_uniform_sampler.hpp" "sampling/direct/segment_uniform_sampler.hpp" "sampling/mcmc/helpers/simulation_runner.hpp" + "sampling/mcmc/helpers/simulation_manager.hpp" + "sampling/mcmc/helpers/simulation_printer.hpp" + "sampling/mcmc/helpers/simulation_monitor.hpp" "sampling/mcmc/models/components/pairwise_term.hpp" "sampling/mcmc/models/components/density_term.hpp" "sampling/mcmc/models/components/energy_term.hpp" diff --git a/src/geode/stochastic/sampling/direct/double_sampler.cpp b/src/geode/stochastic/sampling/direct/double_sampler.cpp index 3736028..5bc6604 100644 --- a/src/geode/stochastic/sampling/direct/double_sampler.cpp +++ b/src/geode/stochastic/sampling/direct/double_sampler.cpp @@ -25,9 +25,110 @@ #include #include +namespace +{ + struct DistributionParams + { + double min; + double max; + double mean; + double std; + }; + + DistributionParams complete_distribution_params( + const geode::DoubleSampler::DistributionDescription& d ) + { + DistributionParams p{}; + if( d.min_value && d.max_value ) + { + p.min = *d.min_value; + p.max = *d.max_value; + p.mean = ( p.min + p.max ) / 2.0; + p.std = ( p.max - p.min ) / std::sqrt( 12.0 ); + return p; + } + + if( d.mean && d.standard_deviation ) + { + p.mean = *d.mean; + p.std = *d.standard_deviation; + p.min = p.mean - std::sqrt( 3.0 ) * p.std; + p.max = p.mean + std::sqrt( 3.0 ) * p.std; + return p; + } + + throw geode::OpenGeodeException( + "[DistributionDescripption] Incomplete distribution description: " + "need at least (min,max) or (mean,std)." ); + return p; + } +} // namespace namespace geode { + struct DistributionTypeHasher + { + std::size_t operator()( const DistributionType& d ) const noexcept + { + // Use the underlying string from NamedType + return absl::Hash< std::string >{}( d.get() ); + } + }; + + using DistributionFactory = std::function< DoubleSampler::Distribution( + const DoubleSampler::DistributionDescription& ) >; + + static absl::flat_hash_map< DistributionType, // key type + DistributionFactory, // value type + DistributionTypeHasher // custom hasher + > + distribution_registry = { + { UniformClosed< double >::distribution_type_static(), + []( const DoubleSampler::DistributionDescription& d ) { + auto p = complete_distribution_params( d ); + UniformClosed< double > dist; + dist.min_value = p.min; + dist.max_value = p.max; + return dist; + } }, + { UniformClosedOpen< double >::distribution_type_static(), + []( const DoubleSampler::DistributionDescription& d ) { + auto p = complete_distribution_params( d ); + UniformClosedOpen< double > dist; + dist.min_value = p.min; + dist.max_value = p.max; + return dist; + } }, + { Gaussian::distribution_type_static(), + []( const DoubleSampler::DistributionDescription& d ) { + auto p = complete_distribution_params( d ); + Gaussian dist; + dist.mean = p.mean; + dist.standard_deviation = p.std; + return dist; + } }, + { TruncatedGaussian::distribution_type_static(), + []( const DoubleSampler::DistributionDescription& d ) { + auto p = complete_distribution_params( d ); + TruncatedGaussian dist; + dist.mean = p.mean; + dist.standard_deviation = p.std; + dist.min_value = p.min; + dist.max_value = p.max; + return dist; + } }, + }; + + DoubleSampler::Distribution DoubleSampler::create_distribution( + const DistributionDescription& desc ) + { + auto it = distribution_registry.find( desc.distribution_type ); + if( it == distribution_registry.end() ) + throw geode::OpenGeodeException( absl::StrCat( + "Unknown distribution type: ", desc.distribution_type.get() ) ); + return it->second( desc ); + } + double DoubleSampler::sample( RandomEngine& engine, const Distribution& dist ) { diff --git a/src/geode/stochastic/spatial/object_set.cpp b/src/geode/stochastic/spatial/object_set.cpp index c2a25a3..8e10128 100644 --- a/src/geode/stochastic/spatial/object_set.cpp +++ b/src/geode/stochastic/spatial/object_set.cpp @@ -87,8 +87,8 @@ namespace geode template < typename Type > std::string ObjectSet< Type >::string() const { - return absl::StrCat( - "ObjectSet ", this->id().string(), " with ", size(), " objects" ); + return absl::StrCat( "ObjectSet ", this->name(), " (", + this->id().string(), ") contains ", size(), " objects" ); } // Explicit template instantiation diff --git a/src/geode/stochastic/spatial/object_sets.cpp b/src/geode/stochastic/spatial/object_sets.cpp index 28dc776..23f7806 100644 --- a/src/geode/stochastic/spatial/object_sets.cpp +++ b/src/geode/stochastic/spatial/object_sets.cpp @@ -4,31 +4,7 @@ #include #include -namespace -{ - template < geode::index_t dimension > - geode::Vector< dimension > enlarge_vector( double distance ); - - template <> - geode::Vector2D enlarge_vector< 2 >( double distance ) - { - return geode::Vector2D{ { distance, distance } }; - } - template <> - geode::Vector3D enlarge_vector< 3 >( double distance ) - { - return geode::Vector3D{ { distance, distance, distance } }; - } - template < geode::index_t dimension > - void enlarge_box( - geode::BoundingBox< dimension >& box, double search_distance ) - { - auto vec = enlarge_vector< dimension >( search_distance ); - box.add_point( box.min() - vec ); - box.add_point( box.max() + vec ); - } -} // namespace namespace geode { template < typename Type > @@ -150,7 +126,7 @@ namespace geode const ObjectId& object_id, double searching_distance ) const { auto box = object_bounding_box( get_object( object_id ) ); - enlarge_box< Type::dim >( box, searching_distance ); + box.extends( searching_distance * 2. ); return neighborhood_.get_all_neighbor_ids( box, object_id ); } @@ -159,7 +135,7 @@ namespace geode const Type& object, double searching_distance ) const { auto box = object_bounding_box( object ); - enlarge_box< Type::dim >( box, searching_distance ); + box.extends( searching_distance * 2. ); return neighborhood_.get_all_neighbor_ids( box, std::nullopt ); } @@ -173,12 +149,11 @@ namespace geode template < typename Type > std::string ObjectSets< Type >::string() const { - auto message = - absl::StrCat( "ObjectSets with ", nb_sets(), " ObjectSet:" ); + auto message = absl::StrCat( "ObjectSets with ", nb_objects(), + " objects in, ", nb_sets(), " sets" ); for( const auto& [set_id, objs] : sets_ ) { - absl::StrAppend( &message, "\n\t - set uuid: ", set_id.string(), - "--> ", objs.string() ); + absl::StrAppend( &message, "\n\t --> ", objs.string() ); } return message; } diff --git a/src/geode/stochastic/spatial/pairwise_interactions.cpp b/src/geode/stochastic/spatial/pairwise_interactions.cpp new file mode 100644 index 0000000..77c5d65 --- /dev/null +++ b/src/geode/stochastic/spatial/pairwise_interactions.cpp @@ -0,0 +1,86 @@ +/* + * Copyright (c) 2019 - 2025 Geode-solutions + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + * + */ + +#include + +#include +#include + +namespace +{ + template < geode::index_t dimension > + double compute_euclidean_distance( const geode::Point< dimension >& point1, + const geode::Point< dimension >& point2 ) + { + return geode::point_point_distance( point1, point2 ); + } + template < geode::index_t dimension > + double compute_euclidean_distance( const geode::Segment< dimension >& seg1, + const geode::Segment< dimension >& seg2 ) + { + return std::get< 0 >( geode::segment_segment_distance( seg1, seg2 ) ); + } +} // namespace +namespace geode +{ + template < typename Type > + EuclideanCutoffInteraction< Type >::EuclideanCutoffInteraction( + double cutoff_distance ) + : PairwiseInteraction< Type >(), cutoff_distance_( cutoff_distance ) + { + } + + template < typename Type > + EuclideanCutoffInteraction< Type >::EuclideanCutoffInteraction( + double cutoff_distance, + typename PairwiseInteraction< Type >::SCOPE scope ) + : PairwiseInteraction< Type >( scope ), + cutoff_distance_( cutoff_distance ) + { + } + + template < typename Type > + double EuclideanCutoffInteraction< Type >::neighborhood_searching_distance() + const + { + return cutoff_distance_; + } + + template < typename Type > + double EuclideanCutoffInteraction< Type >::compute( + const ObjectRef< Type >& object_a, + const ObjectRef< Type >& object_b ) const + { + auto dist = compute_euclidean_distance< Type::dim >( + object_a.object, object_b.object ); + return dist <= cutoff_distance_ ? 1.0 : 0.0; + } + + template class opengeode_stochastic_stochastic_api + EuclideanCutoffInteraction< Point< 2 > >; + template class opengeode_stochastic_stochastic_api + EuclideanCutoffInteraction< Point< 3 > >; + + template class opengeode_stochastic_stochastic_api + EuclideanCutoffInteraction< OwnerSegment< 2 > >; +} // namespace geode \ No newline at end of file diff --git a/tests/stochastic/CMakeLists.txt b/tests/stochastic/CMakeLists.txt index 7705ba6..fbe5d60 100644 --- a/tests/stochastic/CMakeLists.txt +++ b/tests/stochastic/CMakeLists.txt @@ -81,6 +81,22 @@ add_geode_test( ${PROJECT_NAME}::stochastic ) +add_geode_test( + SOURCE "sampling/mcmc/helpers/test-simulation_printer.cpp" + DEPENDENCIES + OpenGeode::basic + OpenGeode::geometry + ${PROJECT_NAME}::stochastic +) + +add_geode_test( + SOURCE "sampling/mcmc/helpers/test-simulation_monitor.cpp" + DEPENDENCIES + OpenGeode::basic + OpenGeode::geometry + ${PROJECT_NAME}::stochastic +) + add_geode_test( SOURCE "sampling/mcmc/models/components/test-intensity-term.cpp" DEPENDENCIES @@ -121,6 +137,14 @@ add_geode_test( ${PROJECT_NAME}::stochastic ) +add_geode_test( + SOURCE "sampling/mcmc/test-mh-fractures.cpp" + DEPENDENCIES + OpenGeode::basic + OpenGeode::geometry + ${PROJECT_NAME}::stochastic +) + add_geode_test( SOURCE "sampling/mcmc/test-mh-poisson.cpp" DEPENDENCIES diff --git a/tests/stochastic/sampling/direct/test-double-sampler.cpp b/tests/stochastic/sampling/direct/test-double-sampler.cpp index c8e5ead..04569ae 100644 --- a/tests/stochastic/sampling/direct/test-double-sampler.cpp +++ b/tests/stochastic/sampling/direct/test-double-sampler.cpp @@ -51,7 +51,7 @@ void test_double_sampler( geode::RandomEngine& engine ) geode::TruncatedGaussian t_gaussian_double; t_gaussian_double.mean = 0.; - t_gaussian_double.mean = 1.; + t_gaussian_double.standard_deviation = 1.; t_gaussian_double.max_value = 1.; t_gaussian_double.min_value = -1.; value = 1000.; @@ -69,6 +69,43 @@ void test_double_sampler( geode::RandomEngine& engine ) value != 1000., "[Gaussian] - value did not changed." ); } +// Test function: create a uniform closed distribution from description +void test_create_uniform_closed( geode::RandomEngine& engine ) +{ + geode::DoubleSampler::DistributionDescription desc; + desc.name = "uniform_closed"; + desc.distribution_type = + geode::UniformClosed< double >::distribution_type_static(); + desc.min_value = 2.; + desc.max_value = 5.; + + auto dist = geode::DoubleSampler::create_distribution( desc ); + double value = geode::DoubleSampler::sample( engine, dist ); + OPENGEODE_EXCEPTION( + value >= 2. && value <= 5., "[UniformClosed] - value out of bounds" ); +} + +// Test function: create a truncated Gaussian distribution from description +void test_create_truncated_gaussian( geode::RandomEngine& engine ) +{ + geode::DoubleSampler::DistributionDescription desc; + desc.name = "truncated_gaussian"; + desc.distribution_type = + geode::TruncatedGaussian::distribution_type_static(); + desc.min_value = -2.; + desc.max_value = 2.; + desc.mean = 0.0; + desc.standard_deviation = 1.0; + + auto dist = geode::DoubleSampler::create_distribution( desc ); + for( int i = 0; i < 100; ++i ) + { + double value = geode::DoubleSampler::sample( engine, dist ); + OPENGEODE_EXCEPTION( value >= -2. && value <= 2., + "[TruncatedGaussian] - value out of bounds" ); + } +} + int main() { try @@ -78,6 +115,9 @@ int main() test_double_sampler( random_engine ); + test_create_uniform_closed( random_engine ); + test_create_truncated_gaussian( random_engine ); + geode::Logger::info( "TEST SUCCESS" ); return 0; } diff --git a/tests/stochastic/sampling/mcmc/helpers/test-simulation_monitor.cpp b/tests/stochastic/sampling/mcmc/helpers/test-simulation_monitor.cpp new file mode 100644 index 0000000..dd985af --- /dev/null +++ b/tests/stochastic/sampling/mcmc/helpers/test-simulation_monitor.cpp @@ -0,0 +1,83 @@ +/* + * Copyright (c) 2019 - 2025 Geode-solutions + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + * + */ + +#include +#include +#include + +#include +#include + +namespace +{ + void test_statistics_monitor_basic() + { + geode::Logger::info( "TEST - StatisticsMonitor basic functionality" ); + + // --- Create monitor with 3 energy terms + geode::StatisticsMonitor monitor( 3 ); + + // --- Add 2 realizations + monitor.add_realization( { 1.0, 2.0, 3.0 } ); + monitor.add_realization( { 2.0, 3.0, 4.0 } ); + + // --- Check count + OPENGEODE_EXCEPTION( monitor.statiscal_count() == 2, "Count mismatch" ); + + const auto& means = monitor.means(); + const auto& variances = monitor.variances(); + + // --- Check means + OPENGEODE_EXCEPTION( + std::fabs( means[0] - 1.5 ) < 1e-12, "Mean[0] incorrect" ); + OPENGEODE_EXCEPTION( + std::fabs( means[1] - 2.5 ) < 1e-12, "Mean[1] incorrect" ); + OPENGEODE_EXCEPTION( + std::fabs( means[2] - 3.5 ) < 1e-12, "Mean[2] incorrect" ); + + // --- Check variances + OPENGEODE_EXCEPTION( + std::fabs( variances[0] - 0.5 ) < 1e-12, "Variance[0] incorrect" ); + OPENGEODE_EXCEPTION( + std::fabs( variances[1] - 0.5 ) < 1e-12, "Variance[1] incorrect" ); + OPENGEODE_EXCEPTION( + std::fabs( variances[2] - 0.5 ) < 1e-12, "Variance[2] incorrect" ); + + geode::Logger::info( "--> SUCCESS!" ); + } +} // namespace + +int main() +{ + try + { + geode::StochasticLibrary::initialize(); + geode::Logger::set_level( geode::Logger::LEVEL::debug ); + test_statistics_monitor_basic(); + return 0; + } + catch( ... ) + { + return geode::geode_lippincott(); + } +} \ No newline at end of file diff --git a/tests/stochastic/sampling/mcmc/helpers/test-simulation_printer.cpp b/tests/stochastic/sampling/mcmc/helpers/test-simulation_printer.cpp new file mode 100644 index 0000000..7806b9e --- /dev/null +++ b/tests/stochastic/sampling/mcmc/helpers/test-simulation_printer.cpp @@ -0,0 +1,154 @@ +/* + * Copyright (c) 2019 - 2025 Geode-solutions + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + * + */ + +#include +#include +#include + +#include +#include + +namespace +{ + void test_print_statistics( + const geode::SimulationPrinterConfigurator& config ) + { + geode::Logger::info( "[TEST] SimulationPrinter print statistics" ); + + geode::SimulationPrinter printer( config ); + + // --- Test print_statistics + std::vector< double > stats = { 1.0, 2.5, 3.7 }; + printer.print_statistics( + stats, "EnergyTerm1;EnergyTerm2;EnergyTerm3" ); + + const std::filesystem::path temp_folder = config.output_folder; + const auto stats_path = temp_folder / config.statistics_filename; + OPENGEODE_EXCEPTION( std::filesystem::exists( stats_path ), + "Statistics file not created" ); + + std::ifstream stats_file( stats_path ); + std::string line; + std::getline( stats_file, line ); + OPENGEODE_EXCEPTION( + line == "# Simulation Statistics", "Header not correctly written" ); + std::getline( stats_file, line ); + OPENGEODE_EXCEPTION( line == "EnergyTerm1;EnergyTerm2;EnergyTerm3", + "Header not correctly written" ); + std::getline( stats_file, line ); + OPENGEODE_EXCEPTION( + line == "1 ; 2.5 ; 3.7", "Statistics line not correctly written" ); + + geode::Logger::info( "--> Success!" ); + } + + void test_statistics_summary( + const geode::SimulationPrinterConfigurator& config ) + { + geode::Logger::info( + "[TEST] SimulationPrinter print statistics summary" ); + + geode::SimulationPrinter printer( config ); + geode::StatisticsMonitor monitor( 2 ); + monitor.add_realization( { 1, 2 } ); + + printer.print_statistics_summary( monitor, "EnergyTerm1;EnergyTerm2" ); + + const std::filesystem::path temp_folder = config.output_folder; + const auto summary_path = + temp_folder / config.statistics_summary_filename; + + OPENGEODE_EXCEPTION( std::filesystem::exists( summary_path ), + "Summary file not created" ); + + std::ifstream summary_file( summary_path ); + std::string content( + ( std::istreambuf_iterator< char >( summary_file ) ), + std::istreambuf_iterator< char >() ); + OPENGEODE_EXCEPTION( + content.find( "EnergyTerm1;EnergyTerm2" ) != std::string::npos, + "Energy term names missing" ); + OPENGEODE_EXCEPTION( + content.find( "2" ) != std::string::npos, "Count missing" ); + OPENGEODE_EXCEPTION( + content.find( "1 ; 2" ) != std::string::npos, "Means missing" ); + OPENGEODE_EXCEPTION( + content.find( "0 ; 0" ) != std::string::npos, "Variances missing" ); + + geode::Logger::info( "--> Success!" ); + } + + void test_print_objects( + const geode::SimulationPrinterConfigurator& config ) + { + geode::Logger::info( "[TEST] SimulationPrinter print object" ); + + geode::ObjectSets< geode::Point2D > object_sets; + const auto set_id = object_sets.add_set( "default_name" ); + + object_sets.add_object( geode::Point2D{ { 0.0, 0.0 } }, set_id ); + object_sets.add_object( geode::Point2D{ { 1.0, 0.0 } }, set_id ); + object_sets.add_object( geode::Point2D{ { 3.0, 0.0 } }, set_id ); + + geode::SimulationPrinter printer( config ); + printer.print_object_sets( object_sets, 0 ); + + const std::filesystem::path temp_folder = config.output_folder; + + const auto obj_path = temp_folder / "pattern_0.txt"; + OPENGEODE_EXCEPTION( std::filesystem::exists( obj_path ), + "Object sets file not created" ); + + geode::Logger::info( "--> Success!" ); + } + +} // namespace + +int main() +{ + try + { + geode::StochasticLibrary::initialize(); + geode::Logger::set_level( geode::Logger::LEVEL::debug ); + + const std::filesystem::path temp_folder = + std::filesystem::current_path() / "simprinter_test"; + SDEBUG( temp_folder ); + std::filesystem::create_directories( temp_folder ); + + geode::SimulationPrinterConfigurator config; + config.output_folder = temp_folder.string(); + config.statistics_filename = "stats.txt"; + config.statistics_summary_filename = "summary.txt"; + config.realisations_prefix = "pattern_"; + + test_print_statistics( config ); + test_statistics_summary( config ); + test_print_objects( config ); + return 0; + } + catch( ... ) + { + return geode::geode_lippincott(); + } +} \ No newline at end of file diff --git a/tests/stochastic/sampling/mcmc/models/components/test-pairwise-term.cpp b/tests/stochastic/sampling/mcmc/models/components/test-pairwise-term.cpp index 49e4325..4dce460 100644 --- a/tests/stochastic/sampling/mcmc/models/components/test-pairwise-term.cpp +++ b/tests/stochastic/sampling/mcmc/models/components/test-pairwise-term.cpp @@ -92,7 +92,7 @@ void test_zero_pairwise( double gamma, auto total = term.total_log( pattern ); OPENGEODE_EXCEPTION( - std::isinf( total ), "[test zero pairwise] - log_total wrong value." ); + std::isinf( total ), "[test zero pairwise] - total_log wrong value." ); geode::Point2D p3{ { 0.5, 0.5 } }; geode::ObjectRef< geode::Point2D > p_ref{ p3, set_id }; @@ -114,6 +114,7 @@ int main() try { geode::StochasticLibrary::initialize(); + geode::Logger::set_level( geode::Logger::LEVEL::debug ); geode::ObjectSets< geode::Point2D > pattern; auto set_id = init_object_set( pattern ); diff --git a/tests/stochastic/sampling/mcmc/test-mh-fractures.cpp b/tests/stochastic/sampling/mcmc/test-mh-fractures.cpp new file mode 100644 index 0000000..2eccec7 --- /dev/null +++ b/tests/stochastic/sampling/mcmc/test-mh-fractures.cpp @@ -0,0 +1,301 @@ +/* + * Copyright (c) 2019 - 2025 Geode-solutions + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + * + */ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace +{ + struct FractureSetDescription + { + std::string name; + + geode::DoubleSampler::DistributionDescription length; + geode::DoubleSampler::DistributionDescription azimuth; + + // positionning + double p20; + // double p21; + double minimal_spacing; + + // mh dynamique + double birth_ratio{ 1.0 }; + double death_ratio{ 1.0 }; + double change_ratio{ 1.0 }; + }; + + class FractureSimulationRunner + : public geode::SimulationRunner< geode::OwnerSegment2D > + { + public: + FractureSimulationRunner( const geode::BoundingBox2D& box ) + : box_( box ) + { + } + + void add_fracture_set_descriptor( + const FractureSetDescription& descriptor ) + { + set_descriptors_.push_back( descriptor ); + } + + void initialize() override + { + auto proposal_kernel = std::make_unique< + geode::ProposalKernel< geode::OwnerSegment2D > >(); + + // 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( set_desc.name ); + name_to_uuid[set_desc.name] = set_id; + + auto length_distribution = + geode::DoubleSampler::create_distribution( + set_desc.length ); + auto azimuth_distribution = + geode::DoubleSampler::create_distribution( + set_desc.azimuth ); + this->set_samplers_.push_back( + std::make_unique< geode::UniformSegmentSetSampler >( + box_, length_distribution, azimuth_distribution ) ); + + 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 density energy terms + for( const auto& set_desc : set_descriptors_ ) + { + const auto set_id = name_to_uuid.at( set_desc.name ); + // p20 + this->ordered_energy_terms_.push_back( + this->energy_terms_collection_.add_energy_term( + std::make_unique< + geode::DensityTerm< geode::OwnerSegment2D > >( + absl::StrCat( set_desc.name, "_density" ), + set_desc.p20, + absl::flat_hash_set< geode::uuid >{ set_id } ) ) ); + // spacing + auto interaction = + std::make_unique< geode::EuclideanCutoffInteraction< + geode::OwnerSegment2D > >( set_desc.minimal_spacing, + geode::PairwiseInteraction< + geode::OwnerSegment2D >::SCOPE::same_set ); + + this->ordered_energy_terms_.push_back( + this->energy_terms_collection_.add_energy_term( + std::make_unique< + geode::PairwiseTerm< geode::OwnerSegment2D > >( + absl::StrCat( set_desc.name, "_min_spacing" ), 0., + absl::flat_hash_set< geode::uuid >{ set_id }, + std::move( interaction ) ) ) ); + } + + this->mh_sampler_ = std::make_unique< + geode::MetropolisHastings< geode::OwnerSegment2D > >( + this->energy_terms_collection_, std::move( proposal_kernel ) ); + } + + void check_statistics( + const geode::StatisticsMonitor& statistic_monitoring ) const + { + const auto& computed_means = statistic_monitoring.means(); + + for( const auto stat_id : + geode::Range{ this->energy_terms_collection_.size() } ) + { + const auto& term = energy_terms_collection_.get( + ordered_energy_terms_[stat_id] ); + geode::Logger::info( "[MH test] Statistic value ", + computed_means[stat_id], + " for energy term: ", term.name().data() ); + } + } + + private: + geode::BoundingBox2D box_; + std::vector< FractureSetDescription > set_descriptors_; + }; + + void test_fracture_simulator() + { + geode::Logger::info( "TEST - MH SINGLE SET FRACTURE SIMULATOR (with " + "intra-set interactions)" ); + + geode::RandomEngine engine; + engine.set_seed( "@mh-test-single-Fracture-set@" ); + + geode::BoundingBox2D box; + box.add_point( geode::Point2D{ { 0.0, 0.0 } } ); + box.add_point( geode::Point2D{ { 100.0, 100.0 } } ); + + // --- Object set + FractureSetDescription setA; + setA.name = "A"; + + // length + setA.length.distribution_type = + geode::UniformClosed< double >::distribution_type_static(); + setA.length.min_value = 1; + setA.length.max_value = 10.; + + // azimuth + setA.azimuth.distribution_type = + geode::UniformClosed< double >::distribution_type_static(); + setA.azimuth.min_value = 1; + setA.azimuth.max_value = 10.; + + // positionning + setA.p20 = 0.1; + setA.minimal_spacing = 1.; + + FractureSimulationRunner runner( box ); + runner.add_fracture_set_descriptor( setA ); + + runner.initialize(); + + // run simulation + geode::SimulationPrinterConfigurator printer_config; + printer_config.output_folder = absl::StrCat( + printer_config.output_folder, "/single_fracture_set" ); + + geode::SimulationConfigurator sim_config; + sim_config.realizations = 1000; + sim_config.metropolis_hasting_steps = 1000; + sim_config.burn_in_steps = 1000; + sim_config.printer = printer_config; + + // runner.run( engine, sim_config ); + auto statistic_monitoring = + runner.run_and_monitor( engine, sim_config ); + runner.check_statistics( statistic_monitoring ); + + geode::Logger::info( "--> SUCCESS!" ); + } + + void test_two_fracture_sets_simulator() + { + geode::Logger::info( "TEST - MH TWO SET FRACTURE SIMULATOR (with " + "intra-set interactions)" ); + + geode::RandomEngine engine; + engine.set_seed( "@mh-test-single-Fracture-set@" ); + + geode::BoundingBox2D box; + box.add_point( geode::Point2D{ { 0.0, 0.0 } } ); + box.add_point( geode::Point2D{ { 100.0, 100.0 } } ); + + // --- Object set + FractureSetDescription setA; + setA.name = "A"; + + // length + setA.length.distribution_type = + geode::UniformClosed< double >::distribution_type_static(); + setA.length.min_value = 1; + setA.length.max_value = 10.; + + // azimuth + setA.azimuth.distribution_type = + geode::UniformClosed< double >::distribution_type_static(); + setA.azimuth.min_value = 1; + setA.azimuth.max_value = 10.; + + // positionning + setA.p20 = 0.1; + setA.minimal_spacing = 1.; + + // --- Object set + FractureSetDescription setB; + setB.name = "B"; + + // length + setB.length.distribution_type = + geode::UniformClosed< double >::distribution_type_static(); + setB.length.min_value = 1; + setB.length.max_value = 10.; + + // azimuth + setB.azimuth.distribution_type = + geode::UniformClosed< double >::distribution_type_static(); + setB.azimuth.min_value = 90.; + setB.azimuth.max_value = 100.; + + // positionning + setB.p20 = 0.1; + setB.minimal_spacing = 2.; + + FractureSimulationRunner runner( box ); + runner.add_fracture_set_descriptor( setA ); + runner.add_fracture_set_descriptor( setB ); + + runner.initialize(); + + // run simulation + geode::SimulationPrinterConfigurator printer_config; + printer_config.output_folder = + absl::StrCat( printer_config.output_folder, "/two_fracture_sets" ); + + geode::SimulationConfigurator sim_config; + sim_config.realizations = 1000; + sim_config.metropolis_hasting_steps = 1000; + sim_config.burn_in_steps = 1000; + sim_config.printer = printer_config; + + // runner.run( engine, sim_config ); + auto statistic_monitoring = + runner.run_and_monitor( engine, sim_config ); + runner.check_statistics( statistic_monitoring ); + + geode::Logger::info( "--> SUCCESS!" ); + } +} // namespace + +int main() +{ + try + { + geode::StochasticLibrary::initialize(); + geode::Logger::set_level( geode::Logger::LEVEL::debug ); + test_fracture_simulator(); + test_two_fracture_sets_simulator(); + return 0; + } + catch( ... ) + { + return geode::geode_lippincott(); + } +} \ No newline at end of file diff --git a/tests/stochastic/sampling/mcmc/test-mh-poisson.cpp b/tests/stochastic/sampling/mcmc/test-mh-poisson.cpp index f95d319..97caaf5 100644 --- a/tests/stochastic/sampling/mcmc/test-mh-poisson.cpp +++ b/tests/stochastic/sampling/mcmc/test-mh-poisson.cpp @@ -40,7 +40,7 @@ namespace struct PoissonDensityDescription { - std::string set_name; + std::string name; double density; double target_count; }; @@ -89,13 +89,13 @@ namespace // Step 2: create energy terms for( const auto& energy_desc : density_descriptors_ ) { - const auto set_id = name_to_uuid.at( energy_desc.set_name ); + const auto set_id = name_to_uuid.at( energy_desc.name ); this->ordered_energy_terms_.push_back( this->energy_terms_collection_.add_energy_term( std::make_unique< geode::DensityTerm< geode::Point2D > >( - absl::StrCat( energy_desc.set_name, "_density" ), + absl::StrCat( energy_desc.name, "_density" ), energy_desc.density, absl::flat_hash_set< geode::uuid >{ set_id } ) ) ); @@ -110,8 +110,11 @@ namespace } void check_statistics( - const geode::MonitoringStatistics& statistic_monitoring ) const + const geode::StatisticsMonitor& statistic_monitoring ) const { + const auto& computed_means = statistic_monitoring.means(); + const auto& computed_variances = statistic_monitoring.variances(); + for( const auto stat_id : geode::Range{ this->energy_terms_collection_.size() } ) { @@ -122,25 +125,22 @@ namespace this->ordered_target_statistics_[stat_id]; const auto target_vs_mean_error = - std::abs( - statistic_monitoring.means[stat_id] - expected_means ) + std::abs( computed_means[stat_id] - expected_means ) / expected_means; OPENGEODE_EXCEPTION( target_vs_mean_error < 0.05, - "[MH test] statistic value ", - statistic_monitoring.means[stat_id], + "[MH test] statistic value ", computed_means[stat_id], " for energy term: ", term.name().data(), " 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 ) + std::abs( computed_variances[stat_id] - expected_means ) / expected_means; OPENGEODE_EXCEPTION( target_vs_variance_error < 0.15, "[MH test] variance of statistic ", - statistic_monitoring.variances[stat_id], + computed_variances[stat_id], " for energy term: ", term.name().data(), " not close enough to expected value ", expected_means, " --> error : ", target_vs_variance_error ); @@ -178,7 +178,7 @@ namespace // --- Energy term description PoissonDensityDescription densityA; - densityA.set_name = "A"; + densityA.name = "A"; densityA.density = 0.3; densityA.target_count = 30.0; @@ -187,12 +187,21 @@ namespace 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 ); + // run simulation + geode::SimulationPrinterConfigurator printer_config; + printer_config.output_folder = + absl::StrCat( printer_config.output_folder, + "/sim_point_poisson_test_", config ); + + geode::SimulationConfigurator sim_config; + sim_config.realizations = 1000; + sim_config.metropolis_hasting_steps = 1000; + sim_config.burn_in_steps = 1000; + sim_config.printer = printer_config; + + // runner.run( engine, sim_config ); + auto statistic_monitoring = + runner.run_and_monitor( engine, sim_config ); runner.check_statistics( statistic_monitoring ); } @@ -231,11 +240,20 @@ namespace runner.initialize(); - constexpr geode::index_t steps = 1000; - constexpr geode::index_t nb_realizations = 1000; + // run simulation + geode::SimulationPrinterConfigurator printer_config; + printer_config.output_folder = absl::StrCat( + printer_config.output_folder, "/sim_point_multitype_poisson_test" ); + + geode::SimulationConfigurator sim_config; + sim_config.realizations = 1000; + sim_config.metropolis_hasting_steps = 1000; + sim_config.burn_in_steps = 1000; + sim_config.printer = printer_config; - auto statistic_monitoring = runner.run_print_and_monitor( - "poisson_statistics", engine, steps, nb_realizations ); + // runner.run( engine, sim_config ); + auto statistic_monitoring = + runner.run_and_monitor( engine, sim_config ); runner.check_statistics( statistic_monitoring ); geode::Logger::info( "--> SUCCESS!" ); diff --git a/tests/stochastic/sampling/mcmc/test-mh-strauss-statistics-and-convergence.cpp b/tests/stochastic/sampling/mcmc/test-mh-strauss-statistics-and-convergence.cpp deleted file mode 100644 index 6117a14..0000000 --- a/tests/stochastic/sampling/mcmc/test-mh-strauss-statistics-and-convergence.cpp +++ /dev/null @@ -1,143 +0,0 @@ -/* - * Copyright (c) 2019 - 2025 Geode-solutions - * - * Permission is hereby granted, free of charge, to any person obtaining a copy - * of this software and associated documentation files (the "Software"), to deal - * in the Software without restriction, including without limitation the rights - * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell - * copies of the Software, and to permit persons to whom the Software is - * furnished to do so, subject to the following conditions: - * - * The above copyright notice and this permission notice shall be included in - * all copies or substantial portions of the Software. - * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR - * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, - * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE - * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER - * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, - * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE - * SOFTWARE. - * - */ -#include -#include -#include -#include -#include -#include -#include -namespace -{ - - void test_convergence( double domain_length, - double poisson_density, - double gamma, - double nb_points, - double nb_paires ) - { - geode::Point2D min_point{ { 0., 0. } }; - geode::Point2D max_point{ { domain_length, domain_length } }; - - geode::BoundingBox2D box; - box.add_point( min_point ); - box.add_point( max_point ); - - double area = domain_length * domain_length; - - geode::uuid set_id; - geode::UniformPointSetSampler< 2 > sampler( box, set_id ); - - geode::GibbsEnergy< geode::Point2D > energy; - energy.add_energy_term( - std::make_unique< geode::DensityTerm< geode::Point2D > >( - "intensity", poisson_density, set_id ) ); - - // auto interaction_fn = - // []( const geode::Point2D& a, const geode::uuid& a_uuid, - // const geode::Point2D& b, const geode::uuid& b_uuid ) { - // geode_unused( a_uuid ); - // geode_unused( b_uuid ); - // auto dx = a.value( 0 ) - b.value( 0 ); - // auto dy = a.value( 1 ) - b.value( 1 ); - // auto dist_sq = dx * dx + dy * dy; - // return dist_sq < 2; - // }; - auto interaction = std::make_unique< - geode::EuclideanCutoffInteraction< geode::Point2D > >( - std::sqrt( 2 ) ); - - energy.add_energy_term( - std::make_unique< geode::PairwiseTerm< geode::Point2D > >( - "interaction", gamma, std::move( interaction ) ) ); - - auto kernel1 = - geode::create_birth_death_change_kernel< geode::Point2D >( - sampler, 0.33, 0.33 ); - geode::MetropolisHastings< geode::Point2D > mh( - energy, std::move( kernel1 ) ); - geode::RandomEngine engine; - engine.set_seed( "@mh-test@" ); - absl::flat_hash_map< geode::uuid, geode::index_t > targets = { { set_id, - 0 } }; - geode::ObjectSets< geode::Point2D > state = - mh.initialize_object_set_with_sampling( engine, targets ); - mh.walk( state, engine, 1000 ); - constexpr geode::index_t N{ 5000 }; - - // Sampling - double sum_points = 0.0; - double sum_paires = 0.0; - for( const auto i : geode::Range{ N } ) - { - mh.walk( state, engine, 1000 ); - auto stats = energy.ordered_energy_term_statistics( state ); - sum_points += stats[1]; - sum_paires += stats[0]; - } - - double mean_points = sum_points / N; - double mean_paires = sum_paires / N; - - geode::Logger::info( "[MH test] mean points = ", mean_points, - " and mean paires = ", mean_paires, " (expected ", nb_points, " ; ", - nb_paires, ") " ); - - OPENGEODE_EXCEPTION( std::abs( mean_points - nb_points ) < 3, - "[MH test] unexpected nb points." ); - - if( nb_paires == 0 ) - { - OPENGEODE_EXCEPTION( - mean_paires == 0, "[MH test] unexpected nb paires." ); - } - else - { - OPENGEODE_EXCEPTION( std::abs( mean_paires - nb_paires ) < 5, - "[MH test] unexpected nb paires." ); - } - } -} // namespace - -int main() -{ - try - { - geode::StochasticLibrary::initialize(); - - double domain_length{ 10. }; - double poisson_density_lambda{ 0.5 }; - std::array< double, 3 > gamma{ 0, 0.5, 1. }; - - test_convergence( domain_length, poisson_density_lambda, 0., 15, 0 ); - test_convergence( domain_length, poisson_density_lambda, 1., 50, 70 ); - test_convergence( domain_length, poisson_density_lambda, 0.5, 24, 10 ); - - geode::Logger::info( "MH STATISTICS AND CONVERGENCE TEST SUCCESS" ); - return 0; - } - catch( ... ) - { - return geode::geode_lippincott(); - } -} \ No newline at end of file diff --git a/tests/stochastic/sampling/mcmc/test-mh-strauss.cpp b/tests/stochastic/sampling/mcmc/test-mh-strauss.cpp index 621965a..13242c2 100644 --- a/tests/stochastic/sampling/mcmc/test-mh-strauss.cpp +++ b/tests/stochastic/sampling/mcmc/test-mh-strauss.cpp @@ -43,14 +43,14 @@ namespace struct PoissonDensityDescription { - std::string set_name; + std::string name; double density; double target_count; }; struct PairwiseInteractionDescription { - std::vector< std::string > set_names; + std::vector< std::string > names; double strength; double distance_threshold; // geode::PairwiseInteraction::SCOPE interaction_scope; @@ -108,13 +108,13 @@ namespace // Step 2: create density energy terms for( const auto& density_desc : density_descriptors_ ) { - const auto set_id = name_to_uuid.at( density_desc.set_name ); + const auto set_id = name_to_uuid.at( density_desc.name ); this->ordered_energy_terms_.push_back( this->energy_terms_collection_.add_energy_term( std::make_unique< geode::DensityTerm< geode::Point2D > >( - absl::StrCat( density_desc.set_name, "_density" ), + absl::StrCat( density_desc.name, "_density" ), density_desc.density, absl::flat_hash_set< geode::uuid >{ set_id } ) ) ); @@ -126,7 +126,7 @@ namespace for( const auto& interaction_desc : interaction_descriptors_ ) { absl::flat_hash_set< geode::uuid > set_ids; - for( const auto& name : interaction_desc.set_names ) + for( const auto& name : interaction_desc.names ) { set_ids.emplace( name_to_uuid.at( name ) ); } @@ -140,8 +140,8 @@ namespace this->energy_terms_collection_.add_energy_term( std::make_unique< geode::PairwiseTerm< geode::Point2D > >( - absl::StrCat( absl::StrJoin( - interaction_desc.set_names, "_" ), + absl::StrCat( + absl::StrJoin( interaction_desc.names, "_" ), "_interaction" ), interaction_desc.strength, set_ids, std::move( interaction ) ) ) ); @@ -157,8 +157,10 @@ namespace } void check_statistics( - const geode::MonitoringStatistics& statistic_monitoring ) const + const geode::StatisticsMonitor& statistic_monitoring ) const { + const auto& computed_means = statistic_monitoring.means(); + for( const auto stat_id : geode::Range{ this->energy_terms_collection_.size() } ) { @@ -167,16 +169,15 @@ namespace const auto expected_mean = this->ordered_target_statistics_[stat_id]; - auto target_vs_mean_error = std::abs( - statistic_monitoring.means[stat_id] - expected_mean ); + auto target_vs_mean_error = + std::abs( computed_means[stat_id] - expected_mean ); if( expected_mean > 0 ) { target_vs_mean_error /= expected_mean; } OPENGEODE_EXCEPTION( target_vs_mean_error < 0.1, - "[MH test] Statistic value ", - statistic_monitoring.means[stat_id], + "[MH test] Statistic value ", computed_means[stat_id], " for energy term: ", term.name().data(), " not close enough to expected value ", expected_mean, " --> error: ", target_vs_mean_error ); @@ -217,13 +218,13 @@ namespace // --- Density term PoissonDensityDescription densityA; - densityA.set_name = "A"; + densityA.name = "A"; densityA.density = 0.5; densityA.target_count = nb_points[config]; // --- Intra-set pairwise interaction (Strauss process) PairwiseInteractionDescription intraA; - intraA.set_names = { "A" }; // same set + intraA.names = { "A" }; // same set intraA.strength = gamma_values[config]; intraA.distance_threshold = 1; // intraA.interaction_scope = @@ -237,13 +238,22 @@ namespace runner.initialize(); - constexpr geode::index_t steps = 1000; - constexpr geode::index_t nb_realizations = 750; - - runner.run( engine, steps ); - auto stats = runner.run_print_and_monitor( - "single_strauss_stats", engine, steps, nb_realizations ); - runner.check_statistics( stats ); + // run simulation + geode::SimulationPrinterConfigurator printer_config; + printer_config.output_folder = + absl::StrCat( printer_config.output_folder, + "/sim_point_strauss_test_", config ); + + geode::SimulationConfigurator sim_config; + sim_config.realizations = 1000; + sim_config.metropolis_hasting_steps = 1000; + sim_config.burn_in_steps = 1000; + sim_config.printer = printer_config; + + // runner.run( engine, sim_config ); + auto statistic_monitoring = + runner.run_and_monitor( engine, sim_config ); + runner.check_statistics( statistic_monitoring ); } geode::Logger::info( "--> SUCCESS!" ); @@ -302,12 +312,22 @@ namespace runner.initialize(); - constexpr geode::index_t steps = 1000; - constexpr geode::index_t nb_realizations = 500; - - auto stats = runner.run_print_and_monitor( - "multi_strauss_stats", engine, steps, nb_realizations ); - runner.check_statistics( stats ); + // run simulation + geode::SimulationPrinterConfigurator printer_config; + printer_config.output_folder = + absl::StrCat( printer_config.output_folder, + "/sim_point_multitype_strauss_test" ); + + geode::SimulationConfigurator sim_config; + sim_config.realizations = 750; + sim_config.metropolis_hasting_steps = 1000; + sim_config.burn_in_steps = 1000; + sim_config.printer = printer_config; + + // runner.run( engine, sim_config ); + auto statistic_monitoring = + runner.run_and_monitor( engine, sim_config ); + runner.check_statistics( statistic_monitoring ); } geode::Logger::info( "--> SUCCESS!" );