diff --git a/include/geode/stochastic/sampling/mcmc/helpers/simulation_runner.hpp b/include/geode/stochastic/sampling/mcmc/helpers/simulation_runner.hpp new file mode 100644 index 0000000..4b9d338 --- /dev/null +++ b/include/geode/stochastic/sampling/mcmc/helpers/simulation_runner.hpp @@ -0,0 +1,232 @@ +/* + * 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 MonitoringStatistics + { + 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 ); + } + + 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; + }; + + template < typename ObjectType > + class SimulationRunner + { + public: + SimulationRunner() = default; + virtual ~SimulationRunner() = default; + + virtual void initialize() = 0; + + const ObjectSets< ObjectType >& run( + RandomEngine& engine, const index_t steps ) + { + mh_sampler_->walk( object_sets_, engine, steps ); + return object_sets_; + } + + void run_and_print( std::string_view filename, + RandomEngine& engine, + const index_t steps, + const index_t nb_realizations ) + { + const auto file_exist = + static_cast< bool >( std::ifstream( filename.data() ) ); + if( !file_exist ) + { + const auto header = statistics_header_file(); + print_to_file( filename, header ); + } + + for( const auto realization : Range{ nb_realizations } ) + { + run( engine, steps ); + const auto statistics = statistics_string(); + print_to_file( filename, statistics ); + } + } + + MonitoringStatistics run_print_and_monitor( std::string_view filename, + RandomEngine& engine, + const index_t steps, + const index_t nb_realizations ) + { + const auto file_exist = + static_cast< bool >( std::ifstream( filename.data() ) ); + if( !file_exist ) + { + const auto header = statistics_header_file(); + print_to_file( filename, header ); + } + MonitoringStatistics stat_monitoring( + energy_terms_collection_.size() ); + + for( const auto realization : Range{ nb_realizations } ) + { + run( engine, steps ); + const auto stats = statistics(); + print_to_file( filename, + absl::StrCat( absl::StrJoin( stats, " ; " ), "\n" ) ); + stat_monitoring.add_realization( stats ); + } + stat_monitoring.finalize( nb_realizations ); + return stat_monitoring; + } + + const ObjectSets< ObjectType >& current_pattern_realization() const + { + return object_sets_; + } + + std::vector< double > statistics() const + { + std::vector< double > statistic_values; + statistic_values.reserve( ordered_energy_terms_.size() ); + + for( const auto& energy_term_uuid : ordered_energy_terms_ ) + { + const auto& term = + energy_terms_collection_.get( energy_term_uuid ); + statistic_values.push_back( term.statistic( object_sets_ ) ); + } + + 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::vector< std::string > term_names; + term_names.reserve( ordered_energy_terms_.size() ); + + for( const auto& energy_term_uuid : ordered_energy_terms_ ) + { + const auto& term = + energy_terms_collection_.get( energy_term_uuid ); + term_names.push_back( term.name().data() ); + } + + 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_; + + std::vector< geode::uuid > ordered_energy_terms_; + std::vector< double > ordered_target_statistics_; + + EnergyTermCollection< ObjectType > energy_terms_collection_; + std::unique_ptr< geode::MetropolisHastings< ObjectType > > mh_sampler_; + + ObjectSets< ObjectType > object_sets_; + }; +} // namespace geode \ No newline at end of file diff --git a/include/geode/stochastic/sampling/mcmc/metropolis_hasting_sampler.hpp b/include/geode/stochastic/sampling/mcmc/metropolis_hasting_sampler.hpp index b31789c..8fa4214 100644 --- a/include/geode/stochastic/sampling/mcmc/metropolis_hasting_sampler.hpp +++ b/include/geode/stochastic/sampling/mcmc/metropolis_hasting_sampler.hpp @@ -49,9 +49,10 @@ namespace geode class MetropolisHastings { public: - MetropolisHastings( GibbsEnergy< ObjectType >& energy, + MetropolisHastings( + const EnergyTermCollection< ObjectType >& energy_term_collection, std::unique_ptr< ProposalKernel< ObjectType > > proposal_kernel ) - : energy_( energy ), + : gibbs_energy_{ energy_term_collection }, proposal_kernel_( std::move( proposal_kernel ) ) { OPENGEODE_ASSERT( @@ -188,7 +189,7 @@ namespace geode { const auto new_object = proposal.new_object(); const auto delta_log_energy = - energy_.delta_log_add( state, new_object ); + gibbs_energy_.delta_log_add( state, new_object ); return accept_or_reject( proposal, state, engine, delta_log_energy, []( auto& state, auto& proposal ) { state.add_object( @@ -203,7 +204,7 @@ namespace geode { const auto old_object_id = proposal.old_object_id(); const auto delta_log_energy = - energy_.delta_log_remove( state, old_object_id ); + gibbs_energy_.delta_log_remove( state, old_object_id ); return accept_or_reject( proposal, state, engine, delta_log_energy, []( auto& state, auto& proposal ) { state.remove_object( proposal.old_object_id() ); @@ -216,8 +217,8 @@ namespace geode { const auto new_object = proposal.new_object(); const auto old_object_id = proposal.old_object_id(); - const auto delta_log_energy = - energy_.delta_log_change( state, old_object_id, new_object ); + 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, @@ -229,7 +230,7 @@ namespace geode }; private: - const GibbsEnergy< ObjectType >& energy_; + GibbsEnergy< ObjectType > gibbs_energy_; std::unique_ptr< ProposalKernel< ObjectType > > proposal_kernel_; double beta_{ 1.0 }; }; diff --git a/include/geode/stochastic/sampling/mcmc/models/components/energy_term.hpp b/include/geode/stochastic/sampling/mcmc/models/components/energy_term.hpp index c327a5e..d3e4376 100644 --- a/include/geode/stochastic/sampling/mcmc/models/components/energy_term.hpp +++ b/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 { diff --git a/include/geode/stochastic/sampling/mcmc/models/components/pairwise_term.hpp b/include/geode/stochastic/sampling/mcmc/models/components/pairwise_term.hpp index 232893c..acf1835 100644 --- a/include/geode/stochastic/sampling/mcmc/models/components/pairwise_term.hpp +++ b/include/geode/stochastic/sampling/mcmc/models/components/pairwise_term.hpp @@ -155,7 +155,7 @@ namespace geode // continue; // } // if( neigh_obj_id.set_id == obj_id.set_id - // && neigh_obj_id.id <= obj_id.id ) + // && neigh_obj_id.index <= obj_id.index ) //{ // continue; // } diff --git a/include/geode/stochastic/sampling/mcmc/models/components/energy_term_collection.hpp b/include/geode/stochastic/sampling/mcmc/models/energy_term_collection.hpp similarity index 96% rename from include/geode/stochastic/sampling/mcmc/models/components/energy_term_collection.hpp rename to include/geode/stochastic/sampling/mcmc/models/energy_term_collection.hpp index 0a85fe6..f8ce986 100644 --- a/include/geode/stochastic/sampling/mcmc/models/components/energy_term_collection.hpp +++ b/include/geode/stochastic/sampling/mcmc/models/energy_term_collection.hpp @@ -63,14 +63,14 @@ namespace geode return energy_terms_.size(); } - [[nodiscard]] std::shared_ptr< const EnergyTerm< ObjectType > > get( + [[nodiscard]] const EnergyTerm< ObjectType >& get( const uuid& id ) const { auto it = energy_terms_.find( id ); OPENGEODE_EXCEPTION( it != energy_terms_.end(), absl::StrCat( "[EnergyTermCollection] Unknown energy term: ", id.string() ) ); - return it->second; + return *it->second; } [[nodiscard]] const absl::flat_hash_map< uuid, diff --git a/include/geode/stochastic/sampling/mcmc/models/gibbs_energy.hpp b/include/geode/stochastic/sampling/mcmc/models/gibbs_energy.hpp index 4a9a489..9bedbe9 100644 --- a/include/geode/stochastic/sampling/mcmc/models/gibbs_energy.hpp +++ b/include/geode/stochastic/sampling/mcmc/models/gibbs_energy.hpp @@ -22,7 +22,7 @@ */ #pragma once -#include +#include #include namespace geode @@ -65,8 +65,7 @@ namespace geode { double log_energy = 0.0; for( const auto& term : - energy_terms_collection_.terms_for_set( new_object - .set_id ) ) // energy_terms_collection_.all_terms() ) + energy_terms_collection_.terms_for_set( new_object.set_id ) ) { log_energy += term->delta_log_add( state, new_object ); } @@ -77,8 +76,8 @@ namespace geode const ObjectSets< ObjectType >& state, const ObjectId& id ) const { double log_energy = 0.0; - for( const auto& term : energy_terms_collection_.terms_for_set( - id.set_id ) ) // energy_terms_collection_.all_terms() ) + for( const auto& term : + energy_terms_collection_.terms_for_set( id.set_id ) ) { log_energy += term->delta_log_remove( state, id ); } @@ -90,8 +89,8 @@ namespace geode const ObjectRef< ObjectType >& new_object ) const { double log_energy = 0.0; - for( const auto& term : energy_terms_collection_.terms_for_set( - old_id.set_id ) ) // energy_terms_collection_.all_terms() ) + for( const auto& term : + energy_terms_collection_.terms_for_set( old_id.set_id ) ) { log_energy += term->delta_log_change( state, old_id, new_object ); diff --git a/include/geode/stochastic/sampling/mcmc/proposal/classical_proposals.hpp b/include/geode/stochastic/sampling/mcmc/proposal/classical_proposals.hpp index d98275e..c44c515 100644 --- a/include/geode/stochastic/sampling/mcmc/proposal/classical_proposals.hpp +++ b/include/geode/stochastic/sampling/mcmc/proposal/classical_proposals.hpp @@ -34,6 +34,7 @@ namespace geode void add_birth_death_change_moves( std::unique_ptr< geode::ObjectSetSampler< ObjectType > >& sampler, geode::ProposalKernel< ObjectType >& kernel, + const uuid& set_id, double birth_ratio, double death_ratio, double change_ratio ) @@ -44,14 +45,14 @@ namespace geode const auto p_birth = birth_ratio / total_ratio; kernel.add_move( - std::make_unique< geode::BirthDeathMove< ObjectType > >( - sampler, total_ratio, p_birth ) ); + set_id, std::make_unique< geode::BirthDeathMove< ObjectType > >( + *sampler, total_ratio, p_birth ) ); if( change_ratio > 0. ) { kernel.add_move( - std::make_unique< geode::ChangeMove< ObjectType > >( - sampler, change_ratio ) ); + set_id, std::make_unique< geode::ChangeMove< ObjectType > >( + *sampler, change_ratio ) ); } } diff --git a/include/geode/stochastic/spatial/object_set.hpp b/include/geode/stochastic/spatial/object_set.hpp index 4d191e3..52ebdba 100644 --- a/include/geode/stochastic/spatial/object_set.hpp +++ b/include/geode/stochastic/spatial/object_set.hpp @@ -36,6 +36,7 @@ namespace geode ObjectSet( ObjectSet&& ) noexcept = default; ObjectSet& operator=( ObjectSet&& ) noexcept = default; + void set_name( std::string_view name ); const Type& get_object( index_t index ) const; index_t add_object( Type&& object ); diff --git a/include/geode/stochastic/spatial/object_sets.hpp b/include/geode/stochastic/spatial/object_sets.hpp index 0ca68ff..a43f1e9 100644 --- a/include/geode/stochastic/spatial/object_sets.hpp +++ b/include/geode/stochastic/spatial/object_sets.hpp @@ -62,7 +62,7 @@ namespace geode index_t nb_objects_in_set( const uuid& set_id ) const; index_t nb_objects() const; - uuid add_set(); + uuid add_set( std::string_view name ); ObjectId add_object( Type&& object, const uuid& set_id ); void update_object( const ObjectId& object_id, Type&& object ); void remove_object( const ObjectId& object_id ); diff --git a/src/geode/stochastic/CMakeLists.txt b/src/geode/stochastic/CMakeLists.txt index 6a63ddf..5eb0368 100644 --- a/src/geode/stochastic/CMakeLists.txt +++ b/src/geode/stochastic/CMakeLists.txt @@ -48,10 +48,11 @@ add_geode_library( "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/models/components/pairwise_term.hpp" "sampling/mcmc/models/components/density_term.hpp" - "sampling/mcmc/models/components/energy_term_collection.hpp" "sampling/mcmc/models/components/energy_term.hpp" + "sampling/mcmc/models/energy_term_collection.hpp" "sampling/mcmc/models/gibbs_energy.hpp" "sampling/mcmc/proposal/classical_proposals.hpp" "sampling/mcmc/proposal/moves.hpp" diff --git a/src/geode/stochastic/spatial/object_set.cpp b/src/geode/stochastic/spatial/object_set.cpp index 6a726bf..c2a25a3 100644 --- a/src/geode/stochastic/spatial/object_set.cpp +++ b/src/geode/stochastic/spatial/object_set.cpp @@ -31,6 +31,12 @@ namespace geode { + template < typename Type > + void ObjectSet< Type >::set_name( std::string_view name ) + { + IdentifierBuilder builder( *this ); + builder.set_name( name ); + } template < typename Type > const Type& ObjectSet< Type >::get_object( index_t index ) const { diff --git a/src/geode/stochastic/spatial/object_sets.cpp b/src/geode/stochastic/spatial/object_sets.cpp index c090997..28dc776 100644 --- a/src/geode/stochastic/spatial/object_sets.cpp +++ b/src/geode/stochastic/spatial/object_sets.cpp @@ -90,9 +90,10 @@ namespace geode } template < typename Type > - uuid ObjectSets< Type >::add_set() + uuid ObjectSets< Type >::add_set( std::string_view name ) { ObjectSet< Type > new_set; + new_set.set_name( name ); const auto new_set_id = new_set.id(); auto [it, inserted] = sets_.emplace( new_set_id, std::move( new_set ) ); OPENGEODE_EXCEPTION( inserted, "[ObjectSet]- group (", diff --git a/tests/stochastic/CMakeLists.txt b/tests/stochastic/CMakeLists.txt index f02b78e..7705ba6 100644 --- a/tests/stochastic/CMakeLists.txt +++ b/tests/stochastic/CMakeLists.txt @@ -113,30 +113,29 @@ add_geode_test( ${PROJECT_NAME}::stochastic ) -# add_geode_test( -# SOURCE "sampling/mcmc/test-metropolis-hasting-sampler.cpp" -# DEPENDENCIES -# OpenGeode::basic -# OpenGeode::geometry -# ${PROJECT_NAME}::stochastic -# ) - -#add_geode_test( -# SOURCE "sampling/mcmc/test-mh-poisson.cpp" -# DEPENDENCIES -# OpenGeode::basic -# OpenGeode::geometry -# ${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-metropolis-hasting-sampler.cpp" + DEPENDENCIES + OpenGeode::basic + OpenGeode::geometry + ${PROJECT_NAME}::stochastic +) + +add_geode_test( + SOURCE "sampling/mcmc/test-mh-poisson.cpp" + DEPENDENCIES + OpenGeode::basic + OpenGeode::geometry + ${PROJECT_NAME}::stochastic +) + +add_geode_test( + SOURCE "sampling/mcmc/test-mh-strauss.cpp" + DEPENDENCIES + OpenGeode::basic + OpenGeode::geometry + ${PROJECT_NAME}::stochastic +) add_geode_test( diff --git a/tests/stochastic/sampling/mcmc/models/components/test-intensity-term.cpp b/tests/stochastic/sampling/mcmc/models/components/test-intensity-term.cpp index 18e4b0a..34e3f71 100644 --- a/tests/stochastic/sampling/mcmc/models/components/test-intensity-term.cpp +++ b/tests/stochastic/sampling/mcmc/models/components/test-intensity-term.cpp @@ -30,7 +30,7 @@ geode::uuid init_object_set( geode::ObjectSets< geode::Point2D > &pattern ) geode::Point2D p1{ { 0., 0. } }; geode::Point2D p2{ { 1., 1. } }; - auto set_id = pattern.add_set(); + auto set_id = pattern.add_set( "default_name" ); pattern.add_object( std::move( p1 ), set_id ); pattern.add_object( std::move( p2 ), set_id ); @@ -97,7 +97,6 @@ int main() try { geode::StochasticLibrary::initialize(); - geode::ObjectSets< geode::Point2D > pattern; auto set_id = init_object_set( pattern ); 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 c88c86d..49e4325 100644 --- a/tests/stochastic/sampling/mcmc/models/components/test-pairwise-term.cpp +++ b/tests/stochastic/sampling/mcmc/models/components/test-pairwise-term.cpp @@ -33,7 +33,7 @@ geode::uuid init_object_set( geode::ObjectSets< geode::Point2D >& pattern ) geode::Point2D p1{ { 0., 0. } }; geode::Point2D p2{ { 1., 1. } }; - auto set_id = pattern.add_set(); + auto set_id = pattern.add_set( "default_name" ); pattern.add_object( std::move( p1 ), set_id ); pattern.add_object( std::move( p2 ), set_id ); diff --git a/tests/stochastic/sampling/mcmc/models/test-gibbs-energy.cpp b/tests/stochastic/sampling/mcmc/models/test-gibbs-energy.cpp index 589168e..b34013c 100644 --- a/tests/stochastic/sampling/mcmc/models/test-gibbs-energy.cpp +++ b/tests/stochastic/sampling/mcmc/models/test-gibbs-energy.cpp @@ -26,8 +26,8 @@ #include #include -#include #include +#include #include #include diff --git a/tests/stochastic/sampling/mcmc/proposal/test-proposal-kernel.cpp b/tests/stochastic/sampling/mcmc/proposal/test-proposal-kernel.cpp index e776b62..15cd1ad 100644 --- a/tests/stochastic/sampling/mcmc/proposal/test-proposal-kernel.cpp +++ b/tests/stochastic/sampling/mcmc/proposal/test-proposal-kernel.cpp @@ -38,7 +38,7 @@ namespace geode::Point2D p1{ { 0., 0. } }; geode::Point2D p2{ { 1., 1. } }; - auto set_id = pattern.add_set(); + auto set_id = pattern.add_set( "default_name" ); pattern.add_object( std::move( p1 ), set_id ); pattern.add_object( std::move( p2 ), set_id ); @@ -156,7 +156,7 @@ namespace // --- Edge case: empty object_set --- geode::ObjectSets< geode::Point2D > empty_config; - const auto empty_set_id = empty_config.add_set(); + const auto empty_set_id = empty_config.add_set( "default_name" ); auto empty_kernel = geode::create_birth_death_change_kernel< geode::Point2D >( empty_set_id, sampler, 0.4, 0.4 ); diff --git a/tests/stochastic/sampling/mcmc/test-metropolis-hasting-sampler.cpp b/tests/stochastic/sampling/mcmc/test-metropolis-hasting-sampler.cpp index 8f352e6..8357cca 100644 --- a/tests/stochastic/sampling/mcmc/test-metropolis-hasting-sampler.cpp +++ b/tests/stochastic/sampling/mcmc/test-metropolis-hasting-sampler.cpp @@ -24,7 +24,7 @@ #include #include #include -#include +#include #include #include namespace @@ -73,21 +73,10 @@ namespace } void test_steps( const geode::MetropolisHastings< geode::Point2D >& mh, - const geode::uuid& set_id ) + geode::ObjectSets< geode::Point2D >& state ) { geode::RandomEngine engine; - absl::flat_hash_map< geode::uuid, geode::index_t > targets = { { set_id, - 20 } }; - // geode::ObjectSets< geode::Point2D > state = - // mh.initialize_object_set_with_sampling( engine, targets ); - - geode::ObjectSets< geode::Point2D > state; - state.add_set( set_id ); - state.add_object( geode::Point2D{ { 1., 1. } }, set_id ); - state.add_object( geode::Point2D{ { 2., 2. } }, set_id ); - state.add_object( geode::Point2D{ { 3., 3. } }, set_id ); - geode::index_t stat_sum{ 0 }; constexpr geode::index_t N{ 100000 }; @@ -135,7 +124,7 @@ namespace } } // should be change... only pone group here - stat_sum += state.nb_objects_in_set( set_id ); + stat_sum += state.nb_objects(); if( count % 1000 == 0 ) { @@ -171,7 +160,8 @@ int main() box.add_point( min_point ); box.add_point( max_point ); - geode::uuid set_id; + geode::ObjectSets< geode::Point2D > state; + const auto set_id = state.add_set( "default_name" ); geode::UniformPointSetSampler< 2 > sampler( box ); double birth_prob = 0.3; @@ -187,12 +177,14 @@ int main() "intensity", 0.5, absl::flat_hash_set< geode::uuid >{ set_id } ) ); - geode::GibbsEnergy< geode::Point2D > poisson_energy( energy_terms ); - geode::MetropolisHastings< geode::Point2D > mh( - poisson_energy, std::move( kernel ) ); + energy_terms, std::move( kernel ) ); + + state.add_object( geode::Point2D{ { 1., 1. } }, set_id ); + state.add_object( geode::Point2D{ { 2., 2. } }, set_id ); + state.add_object( geode::Point2D{ { 3., 3. } }, set_id ); - test_steps( mh, set_id ); + test_steps( mh, state ); test_beta_setter( mh ); test_acceptance_prob_helper(); diff --git a/tests/stochastic/sampling/mcmc/test-mh-poisson.cpp b/tests/stochastic/sampling/mcmc/test-mh-poisson.cpp index 448821d..f95d319 100644 --- a/tests/stochastic/sampling/mcmc/test-mh-poisson.cpp +++ b/tests/stochastic/sampling/mcmc/test-mh-poisson.cpp @@ -22,6 +22,7 @@ */ #include #include +#include #include #include #include @@ -29,161 +30,128 @@ #include namespace { - struct PoissonDescription + struct SetDescription { std::string name; - double density; - - // mh dynamic - double death_birth_ratio{ 2. }; - double birth_ratio{ 0.5 }; - double change_ratio{ 1. }; + double birth_ratio{ 1.0 }; + double death_ratio{ 1.0 }; + double change_ratio{ 1.0 }; }; - struct MultitypePoissonDescription + struct PoissonDensityDescription { - // voi - geode::Point2D min_point; - geode::Point2D max_point; - - // object sets - std::vector< PoissonDescription > set_desc; - - // mh - geode::index_t nb_steps{ 1000 }; - geode::index_t nb_realizations{ 1000 }; + std::string set_name; + double density; + double target_count; }; - struct UserProblem + + class PoissonSimulationRunner + : public geode::SimulationRunner< geode::Point2D > { - geode::BoundingBox2D box; - geode::ObjectSets< geode::Point2D > object_set; - std::vector< geode::uuid > object_set_id; + public: + PoissonSimulationRunner( const geode::BoundingBox2D& box ) + : box_( box ) {}; - geode::GibbsEnergy< geode::Point2D > gibbs_energy; + void add_set_descriptor( const SetDescription& descriptor ) + { + set_descriptors_.push_back( descriptor ); + } + + void add_density_descriptor( + const PoissonDensityDescription& descriptor ) + { + density_descriptors_.push_back( descriptor ); + } - absl::flat_hash_map< geode::uuid, geode::uuid > set_energy_term_ids; + void initialize() override + { + auto proposal_kernel = + std::make_unique< geode::ProposalKernel< geode::Point2D > >(); - absl::flat_hash_map< geode::uuid, double > set_stats_targets; + // Mapping set names -> UUID + std::unordered_map< std::string, geode::uuid > name_to_uuid; - absl::flat_hash_map< geode::uuid, geode::UniformPointSetSampler< 2 > > - set_samplers; + // 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; - std::unique_ptr< geode::MetropolisHastings< geode::Point2D > > - mh_sampler; - }; + this->set_samplers_.push_back( + std::make_unique< geode::UniformPointSetSampler< 2 > >( + box_ ) ); - UserProblem create_problems( - const MultitypePoissonDescription& description ) - { - UserProblem problem; - problem.box.add_point( description.min_point ); - problem.box.add_point( description.max_point ); - double area = problem.box.n_volume(); + 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 ); + } - std::unique_ptr< geode::ProposalKernel< geode::Point2D > > - proposal_kernel = - std::make_unique< geode::ProposalKernel< geode::Point2D > >(); - for( const auto& points_desc : description.set_desc ) - { - auto set_id = problem.object_set.add_set(); - problem.object_set_id.push_back( set_id ); - - // this should be linked to the object subset - // flat_hash_map - problem.set_samplers.emplace( set_id, - geode::UniformPointSetSampler< 2 >{ problem.box, set_id } ); - OPENGEODE_EXCEPTION( points_desc.death_birth_ratio > 0., - "Object cannot be add or removed. Please set a BIRTH-DEATH " - "with a positive probability." ); - proposal_kernel->add_move( - std::make_unique< geode::BirthDeathMove< geode::Point2D > >( - problem.set_samplers.at( set_id ), - points_desc.death_birth_ratio, points_desc.birth_ratio ) ); - if( points_desc.change_ratio > 0. ) + // Step 2: create energy terms + for( const auto& energy_desc : density_descriptors_ ) { - proposal_kernel->add_move( - std::make_unique< geode::ChangeMove< geode::Point2D > >( - problem.set_samplers.at( set_id ), - points_desc.change_ratio ) ); + 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( energy_desc.set_name, "_density" ), + energy_desc.density, + absl::flat_hash_set< geode::uuid >{ set_id } ) ) ); + + this->ordered_target_statistics_.push_back( + energy_desc.target_count ); } - // energy terms here we define intra subset terms (can be - // several of them) need to add inter set energy terms - problem.set_energy_term_ids.emplace( set_id, - problem.gibbs_energy.add_energy_term( - std::make_unique< geode::DensityTerm< geode::Point2D > >( - points_desc.name, points_desc.density, set_id ) ) ); - - problem.set_stats_targets.emplace( - set_id, points_desc.density * area ); + this->mh_sampler_ = + std::make_unique< geode::MetropolisHastings< geode::Point2D > >( + this->energy_terms_collection_, + std::move( proposal_kernel ) ); } - problem.mh_sampler = - std::make_unique< geode::MetropolisHastings< geode::Point2D > >( - problem.gibbs_energy, std::move( proposal_kernel ) ); - return problem; - } - void test_convergence( - const MultitypePoissonDescription& problem_description, - geode::RandomEngine& engine ) - { - auto problem = create_problems( problem_description ); - - problem.mh_sampler->walk( problem.object_set, engine, 500 ); - - std::vector< double > sum_points( problem.object_set_id.size(), 0. ); - std::vector< double > sum_sq_points( problem.object_set_id.size(), 0. ); - auto N = problem_description.nb_realizations; - - for( const auto i : geode::Range{ N } ) + void check_statistics( + const geode::MonitoringStatistics& statistic_monitoring ) const { - problem.mh_sampler->walk( - problem.object_set, engine, problem_description.nb_steps ); - for( const auto set_id : - geode::Range{ problem.object_set_id.size() } ) + for( const auto stat_id : + geode::Range{ this->energy_terms_collection_.size() } ) { - const auto& energy_term_uuid = problem.set_energy_term_ids.at( - problem.object_set_id[set_id] ); - const auto nb_points = - problem.gibbs_energy.energy_term_statistic( - problem.object_set, energy_term_uuid ); - sum_points[set_id] += nb_points; - sum_sq_points[set_id] += nb_points * nb_points; + const auto& term = energy_terms_collection_.get( + ordered_energy_terms_[stat_id] ); + + 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.05, + "[MH test] statistic value ", + statistic_monitoring.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 ) + / expected_means; + + OPENGEODE_EXCEPTION( target_vs_variance_error < 0.15, + "[MH test] variance of statistic ", + statistic_monitoring.variances[stat_id], + " for energy term: ", term.name().data(), + " not close enough to expected value ", expected_means, + " --> error : ", target_vs_variance_error ); } } - std::transform( sum_points.begin(), sum_points.end(), - sum_points.begin(), [N]( double p ) { - return p / N; - } ); - std::transform( sum_sq_points.begin(), sum_sq_points.end(), - sum_sq_points.begin(), [N]( double p ) { - return p / N; - } ); - for( const auto set_id : geode::Range{ problem.object_set_id.size() } ) - { - const auto& set_id = problem.object_set_id[set_id]; - const auto expected_points = problem.set_stats_targets.at( set_id ); - - const auto variance = sum_sq_points[set_id] - - ( sum_points[set_id] * sum_points[set_id] ); - geode::Logger::info( "[MH test] mean points = ", sum_points[set_id], - " and var = ", variance, " (expected ", expected_points, ")" ); - const auto error_stat = - std::abs( sum_points[set_id] - expected_points ) - / expected_points; - OPENGEODE_EXCEPTION( error_stat < 0.02, - "[MH test] mean number of points not close to enought to " - "expected value --> error : ", - error_stat ); - // Variance test: Poisson => Var(N) ≈ E[N] - const auto error_var = - std::abs( variance - expected_points ) / expected_points; - OPENGEODE_EXCEPTION( error_var < 0.15, - "[MH test] Variance not close to enought to " - "expected value --> error : ", - error_var ); - } - } + + private: + geode::BoundingBox2D box_; + std::vector< SetDescription > set_descriptors_; + std::vector< PoissonDensityDescription > density_descriptors_; + }; void test_single_type_poisson() { @@ -192,25 +160,43 @@ namespace geode::RandomEngine engine; engine.set_seed( "@mh-test-single-POISSON@" ); - std::array< double, 4 > birth_ratio{ 0.1, 0.3, 0.7, 0.9 }; + geode::BoundingBox2D box; + 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() } ) { - MultitypePoissonDescription problem_description; - problem_description.min_point = geode::Point2D{ { 0., 0. } }; - problem_description.max_point = geode::Point2D{ { 10., 10. } }; - - std::vector< PoissonDescription > poisson_description{ { "set01", - 0.5, 1., birth_ratio[config], change_ratio[config] } }; - problem_description.set_desc = poisson_description; - - problem_description.nb_steps = 1000.; - problem_description.nb_realizations = 1000.; - - test_convergence( problem_description, engine ); + // --- Set description + SetDescription setA; + setA.name = "A"; + 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( "MH SINGLE TYPE POISSON -- SUCCESS!" ); + + geode::Logger::info( "--> SUCCESS!" ); } void test_multitype_poisson() @@ -220,21 +206,39 @@ namespace geode::RandomEngine engine; engine.set_seed( "@mh-test-POISSON-multi@" ); - MultitypePoissonDescription problem_description; - problem_description.min_point = geode::Point2D{ { 0., 0. } }; - problem_description.max_point = geode::Point2D{ { 10., 10. } }; + geode::BoundingBox2D box; + box.add_point( geode::Point2D{ { 0.0, 0.0 } } ); + box.add_point( geode::Point2D{ { 10.0, 10.0 } } ); + + // --- 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 }; + + // --- Energy term descriptions + PoissonDensityDescription density01{ "set01", 0.1, 10.0 }; + PoissonDensityDescription density02{ "set02", 0.4, 40.0 }; + PoissonDensityDescription density03{ "set03", 0.3, 30.0 }; + + PoissonSimulationRunner runner( box ); + runner.add_set_descriptor( set01 ); + runner.add_set_descriptor( set02 ); + runner.add_set_descriptor( set03 ); + + runner.add_density_descriptor( density01 ); + runner.add_density_descriptor( density02 ); + runner.add_density_descriptor( density03 ); + + runner.initialize(); - std::vector< PoissonDescription > poisson_description{ - { "set01", 0.1, 3., 0.25, 1. }, { "set02", 0.4, 1., 0.75, 1. }, - { "set03", 0.3, 5., 0.5, 1. } - }; - problem_description.set_desc = poisson_description; + constexpr geode::index_t steps = 1000; + constexpr geode::index_t nb_realizations = 1000; - problem_description.nb_steps = 1000.; - problem_description.nb_realizations = 1000.; + auto statistic_monitoring = runner.run_print_and_monitor( + "poisson_statistics", engine, steps, nb_realizations ); + runner.check_statistics( statistic_monitoring ); - test_convergence( problem_description, engine ); - geode::Logger::info( "MH MULTITYPE POISSON -- SUCCESS!" ); + geode::Logger::info( "--> SUCCESS!" ); } } // namespace @@ -243,6 +247,7 @@ int main() try { geode::StochasticLibrary::initialize(); + geode::Logger::set_level( geode::Logger::LEVEL::debug ); test_single_type_poisson(); test_multitype_poisson(); return 0; diff --git a/tests/stochastic/sampling/mcmc/test-mh-strauss.cpp b/tests/stochastic/sampling/mcmc/test-mh-strauss.cpp index 47155eb..621965a 100644 --- a/tests/stochastic/sampling/mcmc/test-mh-strauss.cpp +++ b/tests/stochastic/sampling/mcmc/test-mh-strauss.cpp @@ -22,407 +22,295 @@ */ #include #include +#include #include #include #include #include #include #include +#include namespace { - struct StraussDescription + struct SetDescription { std::string name; - double density; - double expected_number_of_objects; - - // interaction - geode::PairwiseInteraction< geode::Point2D >::SCOPE scope; - double distance_treshold; - double gamma; - double expected_number_of_intersections; - - // mh dynamic - double death_birth_ratio{ 1. }; - double birth_ratio{ 0.5 }; - double change_ratio{ 1. }; + double birth_ratio{ 1.0 }; + double death_ratio{ 1.0 }; + double change_ratio{ 1.0 }; }; - struct MultitypeStraussDescription + struct PoissonDensityDescription { - // voi - geode::Point2D min_point; - geode::Point2D max_point; - - // object sets - std::vector< StraussDescription > set_desc; + std::string set_name; + double density; + double target_count; + }; - // mh - geode::index_t nb_steps{ 10000 }; - geode::index_t nb_realizations{ 1000 }; + struct PairwiseInteractionDescription + { + std::vector< std::string > set_names; + double strength; + double distance_threshold; + // geode::PairwiseInteraction::SCOPE interaction_scope; + double target_interaction_count; }; - struct UserProblem + + class StraussSimulationRunner + : public geode::SimulationRunner< geode::Point2D > { - geode::BoundingBox2D box; - geode::ObjectSets< geode::Point2D > object_set; - std::vector< geode::uuid > object_set_id; + public: + StraussSimulationRunner( const geode::BoundingBox2D& box ) : box_( box ) + { + } - geode::GibbsEnergy< geode::Point2D > gibbs_energy; + void add_set_descriptor( const SetDescription& descriptor ) + { + set_descriptors_.push_back( descriptor ); + } - absl::flat_hash_map< geode::uuid, std::vector< geode::uuid > > - set_energy_term_ids; - absl::flat_hash_map< geode::uuid, std::vector< double > > - set_stats_targets; + void add_density_descriptor( + const PoissonDensityDescription& descriptor ) + { + density_descriptors_.push_back( descriptor ); + } - absl::flat_hash_map< geode::uuid, geode::UniformPointSetSampler< 2 > > - set_samplers; - std::unique_ptr< geode::MetropolisHastings< geode::Point2D > > - mh_sampler; + void add_interaction_descriptor( + const PairwiseInteractionDescription& descriptor ) + { + interaction_descriptors_.push_back( descriptor ); + } - std::string string() const + void initialize() override { - std::string message( "User Problem for stochastic simulation: " ); - absl::StrAppend( - &message, "\n\t - VOI: BoundingBox ", box.string() ); - absl::StrAppend( &message, "\n\t - ", object_set.string() ); - absl::StrAppend( &message, "\n\t - subset uuid list: " ); - for( const auto& uuid : object_set_id ) + auto proposal_kernel = + std::make_unique< geode::ProposalKernel< geode::Point2D > >(); + + // 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_ ) { - absl::StrAppend( - &message, "\n\t --> subset uuid: ", uuid.string() ); + 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_ ) ); + + 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 ); } - absl::StrAppend( &message, "\n\t - ", gibbs_energy.string() ); - absl::StrAppend( - &message, "\n\t - subset energy term uuid list: " ); - for( const auto& [uuid, et_uuids] : set_energy_term_ids ) + + // Step 2: create density energy terms + for( const auto& density_desc : density_descriptors_ ) { - absl::StrAppend( - &message, "\n\t --> subset uuid: ", uuid.string() ); - for( const auto& et_id : et_uuids ) - { - absl::StrAppend( - &message, "; energy term uuid: ", et_id.string() ); - } + const auto set_id = name_to_uuid.at( density_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( density_desc.set_name, "_density" ), + density_desc.density, + absl::flat_hash_set< geode::uuid >{ set_id } ) ) ); + + this->ordered_target_statistics_.push_back( + density_desc.target_count ); } - absl::StrAppend( &message, "\n\t - subset stat target list: " ); - for( const auto& [uuid, stats] : set_stats_targets ) + // Step 3: create pairwise interaction terms + for( const auto& interaction_desc : interaction_descriptors_ ) { - absl::StrAppend( - &message, "\n\t --> subset uuid: ", uuid.string() ); - for( const auto& stat : stats ) + absl::flat_hash_set< geode::uuid > set_ids; + for( const auto& name : interaction_desc.set_names ) { - absl::StrAppend( &message, "; stat: ", stat ); + set_ids.emplace( name_to_uuid.at( name ) ); } + + auto interaction = std::make_unique< + geode::EuclideanCutoffInteraction< geode::Point2D > >( + interaction_desc.distance_threshold + /*,interaction_desc.interaction_scope*/ ); + + this->ordered_energy_terms_.push_back( + this->energy_terms_collection_.add_energy_term( + std::make_unique< + geode::PairwiseTerm< geode::Point2D > >( + absl::StrCat( absl::StrJoin( + interaction_desc.set_names, "_" ), + "_interaction" ), + interaction_desc.strength, set_ids, + std::move( interaction ) ) ) ); + + this->ordered_target_statistics_.push_back( + interaction_desc.target_interaction_count ); } - absl::StrAppend( - &message, "\n\t - subset sampler: ", set_samplers.size() ); - for( const auto& [uuid, sampler] : set_samplers ) - { - absl::StrAppend( - &message, "\n\t --> subset uuid: ", uuid.string() ); - } - absl::StrAppend( &message, "\n\t END " ); - return message; - } - }; - UserProblem create_problems( - const MultitypeStraussDescription& description ) - { - UserProblem problem; - problem.box.add_point( description.min_point ); - problem.box.add_point( description.max_point ); - double area = problem.box.n_volume(); + this->mh_sampler_ = + std::make_unique< geode::MetropolisHastings< geode::Point2D > >( + this->energy_terms_collection_, + std::move( proposal_kernel ) ); + } - std::unique_ptr< geode::ProposalKernel< geode::Point2D > > - proposal_kernel = - std::make_unique< geode::ProposalKernel< geode::Point2D > >(); - for( const auto& points_desc : description.set_desc ) + void check_statistics( + const geode::MonitoringStatistics& statistic_monitoring ) const { - auto set_id = problem.object_set.add_set(); - problem.object_set_id.push_back( set_id ); - - problem.set_samplers.emplace( set_id, - geode::UniformPointSetSampler< 2 >{ problem.box, set_id } ); - OPENGEODE_EXCEPTION( points_desc.death_birth_ratio > 0., - "Object cannot be add or removed. Please set a BIRTH-DEATH " - "with a positive probability." ); - proposal_kernel->add_move( - std::make_unique< geode::BirthDeathMove< geode::Point2D > >( - problem.set_samplers.at( set_id ), - points_desc.death_birth_ratio, points_desc.birth_ratio ) ); - if( points_desc.change_ratio > 0. ) + for( const auto stat_id : + geode::Range{ this->energy_terms_collection_.size() } ) { - proposal_kernel->add_move( - std::make_unique< geode::ChangeMove< geode::Point2D > >( - problem.set_samplers.at( set_id ), - points_desc.change_ratio ) ); - } + const auto& term = energy_terms_collection_.get( + ordered_energy_terms_[stat_id] ); + + const auto expected_mean = + this->ordered_target_statistics_[stat_id]; + auto target_vs_mean_error = std::abs( + statistic_monitoring.means[stat_id] - expected_mean ); + if( expected_mean > 0 ) + { + target_vs_mean_error /= expected_mean; + } - // energy terms here we define intra subset terms (can be - // several of them) need to add inter set energy terms - std::vector< geode::uuid > energy_terms; - energy_terms.push_back( problem.gibbs_energy.add_energy_term( - std::make_unique< geode::DensityTerm< geode::Point2D > >( - points_desc.name, points_desc.density, set_id ) ) ); - - auto interaction = std::make_unique< - geode::EuclideanCutoffInteraction< geode::Point2D > >( - points_desc.distance_treshold, points_desc.scope ); - energy_terms.push_back( problem.gibbs_energy.add_energy_term( - std::make_unique< geode::PairwiseTerm< geode::Point2D > >( - "interaction", points_desc.gamma, std::move( interaction ), - set_id ) ) ); - problem.set_energy_term_ids.emplace( set_id, energy_terms ); - - std::vector< double > expected_statistics; - expected_statistics.push_back( - points_desc.expected_number_of_objects ); - expected_statistics.push_back( - points_desc.expected_number_of_intersections ); - problem.set_stats_targets.emplace( set_id, expected_statistics ); + OPENGEODE_EXCEPTION( target_vs_mean_error < 0.1, + "[MH test] Statistic value ", + statistic_monitoring.means[stat_id], + " for energy term: ", term.name().data(), + " not close enough to expected value ", expected_mean, + " --> error: ", target_vs_mean_error ); + } } - DEBUG( proposal_kernel->string() ); - problem.mh_sampler = - std::make_unique< geode::MetropolisHastings< geode::Point2D > >( - problem.gibbs_energy, std::move( proposal_kernel ) ); - return problem; - } - void test_convergence( - const MultitypeStraussDescription& problem_description, - geode::RandomEngine& engine ) - { - auto problem = create_problems( problem_description ); - SDEBUG( problem ); - - problem.mh_sampler->walk( problem.object_set, engine, 500 ); - - // std::vector< double > sum_points( - // problem.object_set_id.size(), 0. ); std::vector< double > - // sum_nb_interactions( - // problem.object_set_id.size(), 0. ); - // - // auto N = problem_description.nb_realizations; - - // for( const auto i : geode::Range{ N } ) - // { - // problem.mh_sampler->walk( - // problem.object_set, engine, - // problem_description.nb_steps ); - // for( const auto set_id : - // geode::Range{ problem.object_set_id.size() } ) - // { - // const auto& energy_term_uuids = - // problem.set_energy_term_ids.at( - // problem.object_set_id[set_id] ); - // sum_points[set_id] += - // problem.gibbs_energy.energy_term_statistic( - // problem.object_set, energy_term_uuids[0] ); - // sum_nb_interactions[set_id] += - // problem.gibbs_energy.energy_term_statistic( - // problem.object_set, energy_term_uuids[1] ); - // } - // } - // std::transform( sum_points.begin(), sum_points.end(), - // sum_points.begin(), [N]( double p ) { - // return p / N; - // } ); - // std::transform( sum_nb_interactions.begin(), - // sum_nb_interactions.end(), - // sum_nb_interactions.begin(), [N]( double p ) { - // return p / N; - // } ); - // for( const auto set_id : - // geode::Range{ problem.object_set_id.size() } ) - // { - // const auto& set_id = - // problem.object_set_id[set_id]; const auto& - // expected_stats = - // problem.set_stats_targets.at( set_id ); - // - // geode::Logger::info( "[MH test] mean points = ", - // sum_points[set_id], " (expected ", - // expected_stats[0], - // ")" - // " and mean interactions = ", - // sum_nb_interactions[set_id], " (expected ", - // expected_stats[1], ")" ); - // const auto error_stat = - // std::abs( sum_points[set_id] - expected_stats[0] ) - // / expected_stats[0]; - // // OPENGEODE_EXCEPTION( error_stat < 0.015, - // // "[MH test] mean number of points not - // close to - // // enought to " " expected value-- > error: - // ", - // // error_stat ); - // if( expected_stats[1] == 0 ) - // { - // OPENGEODE_EXCEPTION( - // sum_nb_interactions[set_id] < - // geode::GLOBAL_EPSILON, - // "[MH test] Number of interactions not close to - // enought to " " expected value-- > error : ", - // sum_nb_interactions[set_id] ); - // } - // else - // { - // const auto error_interactions = - // std::abs( - // sum_nb_interactions[set_id] - - // expected_stats[1] ) - // / expected_stats[1]; - // // OPENGEODE_EXCEPTION( - // error_interactions < 0.1, - // // "[MH test] Number of - // interactions not - // // close to enought to " " expected - // value-- > - // // error : ", error_interactions ); - // } - // } - SDEBUG( problem ); - } + private: + geode::BoundingBox2D box_; + std::vector< SetDescription > set_descriptors_; + std::vector< PoissonDensityDescription > density_descriptors_; + std::vector< PairwiseInteractionDescription > interaction_descriptors_; + }; - void test_single_type_poisson() + void test_single_type_strauss() { - geode::Logger::info( "TEST - MH SINGLE TYPE STRAUSS" ); + geode::Logger::info( + "TEST - MH SINGLE TYPE STRAUSS (with intra-set interactions)" ); geode::RandomEngine engine; - engine.set_seed( "@mh-test-Strauss-single@" ); + engine.set_seed( "@mh-test-single-STRAUSS@" ); + + geode::BoundingBox2D box; + box.add_point( geode::Point2D{ { 0.0, 0.0 } } ); + box.add_point( geode::Point2D{ { 10.0, 10.0 } } ); + std::array< double, 5 > gamma_values{ 0, 0.3, 0.5, 0.7, 1.0 }; std::array< double, 5 > nb_points{ 22.6, 27.4, 31.3, 36.1, 50. }; - std::array< double, 5 > nb_interactions{ 0, 4, 8, 13, 36 }; + std::array< double, 5 > nb_interactions{ 0, 4, 8, 14, 36 }; for( const auto config : geode::Range{ gamma_values.size() } ) { - MultitypeStraussDescription problem_description; - problem_description.min_point = geode::Point2D{ { 0., 0. } }; - problem_description.max_point = geode::Point2D{ { 10., 10. } }; - - StraussDescription description1; - description1.name = "set1"; - description1.density = 0.5; - description1.expected_number_of_objects = nb_points[config]; - - description1.scope = - geode::PairwiseInteraction< geode::Point2D >::SCOPE::all_set; - description1.distance_treshold = 1; - description1.gamma = gamma_values[config]; - description1.expected_number_of_intersections = - nb_interactions[config]; - - // mh dynamic - description1.death_birth_ratio = 1.; - description1.birth_ratio = 0.5; - description1.change_ratio = 0.; - - std::vector< StraussDescription > poisson_description{ - description1 - }; - problem_description.set_desc = poisson_description; - - problem_description.nb_steps = 1000.; - problem_description.nb_realizations = 1000.; - - test_convergence( problem_description, engine ); + // --- Object set + SetDescription setA; + setA.name = "A"; + setA.birth_ratio = 1.0; + setA.death_ratio = 1.0; + setA.change_ratio = 1.0; + + // --- Density term + PoissonDensityDescription densityA; + densityA.set_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.strength = gamma_values[config]; + intraA.distance_threshold = 1; + // intraA.interaction_scope = + // geode::PairwiseInteraction::SCOPE::INTRA; + intraA.target_interaction_count = nb_interactions[config]; + + StraussSimulationRunner runner( box ); + runner.add_set_descriptor( setA ); + runner.add_density_descriptor( densityA ); + runner.add_interaction_descriptor( intraA ); + + 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 ); } - geode::Logger::info( "TEST - MH SINGLE TYPE STRAUSS ... SUCCESS!" ); + + geode::Logger::info( "--> SUCCESS!" ); } - void test_multi_type_strauss() + void test_multitype_strauss() { - geode::Logger::info( "TEST - MH MULTI TYPE STRAUSS" ); + geode::Logger::info( + "TEST - MH MULTITYPE STRAUSS (with inter-set interactions)" ); geode::RandomEngine engine; - engine.set_seed( "@mh-test-Strauss-multi@" ); - std::array< double, 4 > gamma_values{ 0, 0.3, 0.7, 1.0 }; - std::array< double, 4 > nb_points{ 22.6, 27.4, 31.3, 36.1 }; - std::array< double, 4 > nb_interactions{ 0, 4, 8, 13 }; - - MultitypeStraussDescription problem_description; - problem_description.min_point = geode::Point2D{ { 0., 0. } }; - problem_description.max_point = geode::Point2D{ { 10., 10. } }; - - StraussDescription description1; - description1.name = "set1"; - description1.density = 0.5; - description1.expected_number_of_objects = nb_points[0]; - - description1.scope = - geode::PairwiseInteraction< geode::Point2D >::SCOPE::same_set; - description1.distance_treshold = 1; - description1.gamma = gamma_values[0]; - description1.expected_number_of_intersections = nb_interactions[0]; - - // mh dynamic - description1.death_birth_ratio = 1.; - description1.birth_ratio = 0.5; - description1.change_ratio = 0.; - - StraussDescription description2; - description2.name = "set2"; - description2.density = 0.5; - description2.expected_number_of_objects = nb_points[1]; - - description2.scope = - geode::PairwiseInteraction< geode::Point2D >::SCOPE::same_set; - description2.distance_treshold = 2; - description2.gamma = gamma_values[1]; - description2.expected_number_of_intersections = nb_interactions[1]; - - // mh dynamic - description2.death_birth_ratio = 1.; - description2.birth_ratio = 0.5; - description2.change_ratio = 0.; - - StraussDescription description3; - description3.name = "set3"; - description3.density = 0.3; - description3.expected_number_of_objects = nb_points[2]; - - description3.scope = - geode::PairwiseInteraction< geode::Point2D >::SCOPE::same_set; - description3.distance_treshold = 1; - description3.gamma = gamma_values[2]; - description3.expected_number_of_intersections = nb_interactions[2]; - - // mh dynamic - description3.death_birth_ratio = 1.; - description3.birth_ratio = 0.5; - description3.change_ratio = 0.; - - StraussDescription description4; - description4.name = "set4"; - description4.density = 0.5; - description4.expected_number_of_objects = nb_points[3]; - - description4.scope = - geode::PairwiseInteraction< geode::Point2D >::SCOPE::same_set; - description4.distance_treshold = 1; - description4.gamma = gamma_values[3]; - description4.expected_number_of_intersections = nb_interactions[3]; - - // mh dynamic - description4.death_birth_ratio = 1.; - description4.birth_ratio = 0.5; - description4.change_ratio = 0.; - - std::vector< StraussDescription > strauss_description{ // description1, - // description2, - description3, description4 - }; - problem_description.set_desc = strauss_description; - - problem_description.nb_steps = 10000.; - problem_description.nb_realizations = 1000.; - - test_convergence( problem_description, engine ); - - geode::Logger::info( "TEST - MH SINGLE TYPE STRAUSS ... SUCCESS!" ); + engine.set_seed( "@mh-test-multi-STRAUSS@" ); + + geode::BoundingBox2D box; + box.add_point( geode::Point2D{ { 0.0, 0.0 } } ); + box.add_point( geode::Point2D{ { 10.0, 10.0 } } ); + + std::array< double, 3 > gamma_values{ 0, 0.5, 1.0 }; + std::array< double, 3 > nb_points01{ 3.5, 5, 10.0 }; + std::array< double, 3 > nb_points02{ 14, 21, 40.0 }; + std::array< double, 3 > nb_points03{ 11, 16, 30. }; + std::array< double, 3 > nb_interactions01{ 0, 15, 95 }; + std::array< double, 3 > nb_interactions02{ 16, 41, 171 }; + + for( const auto config : geode::Range{ gamma_values.size() } ) + { + // --- Sets + SetDescription set01{ "set01", 1.0, 3.0, 1.0 }; + SetDescription set02{ "set02", 3.0, 0.5, 1.0 }; + SetDescription set03{ "set03", 4.0, 1.0, 1.0 }; + + // --- Density terms + PoissonDensityDescription d01{ "set01", 0.1, nb_points01[config] }; + PoissonDensityDescription d02{ "set02", 0.4, nb_points02[config] }; + PoissonDensityDescription d03{ "set03", 0.3, nb_points03[config] }; + + // --- Pairwise interactions + // 1. Intra-type (repulsion within same set) + PairwiseInteractionDescription intra01{ { "set01", "set02", + "set03" }, + gamma_values[config], 1., nb_interactions01[config] }; + PairwiseInteractionDescription intra02{ { "set02" }, 1., 2., + nb_interactions02[config] }; + + StraussSimulationRunner runner( box ); + runner.add_set_descriptor( set01 ); + runner.add_set_descriptor( set02 ); + runner.add_set_descriptor( set03 ); + + runner.add_density_descriptor( d01 ); + runner.add_density_descriptor( d02 ); + runner.add_density_descriptor( d03 ); + + runner.add_interaction_descriptor( intra01 ); + runner.add_interaction_descriptor( intra02 ); + + 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 ); + } + + geode::Logger::info( "--> SUCCESS!" ); } } // namespace @@ -431,9 +319,9 @@ int main() try { geode::StochasticLibrary::initialize(); - geode::Logger::set_level( geode::Logger::LEVEL::trace ); - // test_single_type_poisson(); - test_multi_type_strauss(); + geode::Logger::set_level( geode::Logger::LEVEL::debug ); + test_single_type_strauss(); + test_multitype_strauss(); return 0; } catch( ... ) diff --git a/tests/stochastic/spatial/test-object-sets.cpp b/tests/stochastic/spatial/test-object-sets.cpp index 3bbcc75..5fbef25 100644 --- a/tests/stochastic/spatial/test-object-sets.cpp +++ b/tests/stochastic/spatial/test-object-sets.cpp @@ -30,8 +30,8 @@ namespace void test_add_sets_and_objects() { ObjectSets< geode::Point2D > sets; - const auto set_id1 = sets.add_set(); - const auto set_id2 = sets.add_set(); + const auto set_id1 = sets.add_set( "default_name" ); + const auto set_id2 = sets.add_set( "default_name" ); OPENGEODE_EXCEPTION( sets.nb_sets() == 2, "[TestObjectSets] - Expected 2 sets" ); @@ -59,7 +59,7 @@ namespace void test_neighbors_by_object_id() { ObjectSets< geode::Point2D > sets; - const auto set_id = sets.add_set(); + const auto set_id = sets.add_set( "default_name" ); const auto obj0 = sets.add_object( geode::Point2D{ { 0.0, 0.0 } }, set_id ); @@ -84,7 +84,7 @@ namespace void test_neighbors_by_object_value() { ObjectSets< geode::Point2D > sets; - const auto set_id = sets.add_set(); + const auto set_id = sets.add_set( "default_name" ); sets.add_object( geode::Point2D{ { 0.0, 0.0 } }, set_id ); sets.add_object( geode::Point2D{ { 1.0, 0.0 } }, set_id ); @@ -102,7 +102,7 @@ namespace void test_string_representation() { ObjectSets< geode::Point2D > sets; - const auto set_id = sets.add_set(); + const auto set_id = sets.add_set( "default_name" ); sets.add_object( geode::Point2D{ { 0.0, 0.0 } }, set_id ); sets.add_object( geode::Point2D{ { 1.0, 1.0 } }, set_id );