diff --git a/include/geode/stochastic/inference/statistics_tools.hpp b/include/geode/stochastic/inference/statistics_tools.hpp index ae0abf2..4a426be 100644 --- a/include/geode/stochastic/inference/statistics_tools.hpp +++ b/include/geode/stochastic/inference/statistics_tools.hpp @@ -27,50 +27,47 @@ #include #include -namespace geode +namespace geode::statistics { - namespace statistics + template < typename ObjectType > + void validate( const StatisticsTracker< ObjectType >& tracker, + const TargetStatistics< ObjectType >& targets ) { - template < typename ObjectType > - void validate( const StatisticsTracker< ObjectType >& tracker, - const TargetStatistics< ObjectType >& targets ) - { - const auto& model = targets.model(); + const auto& model = targets.model(); - for( const auto& term_uuid : targets.active_terms() ) - { - const auto mean = tracker.mean( term_uuid ); - const auto target = targets.target( term_uuid ); + for( const auto& term_uuid : targets.active_terms() ) + { + const auto mean = tracker.mean( term_uuid ); + const auto target = targets.target( term_uuid ); - const auto rel_error = - std::fabs( mean - target ) - / ( std::fabs( target ) + geode::GLOBAL_EPSILON ); + const auto rel_error = + std::fabs( mean - target ) + / ( std::fabs( target ) + geode::GLOBAL_EPSILON ); - OpenGeodeStochasticStochasticException::check_exception( - rel_error < targets.tolerance( term_uuid ), nullptr, - OpenGeodeException::TYPE::result, - "[StatisticsValidator] Failure for term ", - model.term_name( term_uuid ), "\n mean = ", mean, - "\n target = ", target, "\n error = ", rel_error, - "\n tol = ", targets.tolerance( term_uuid ) ); - } + OpenGeodeStochasticStochasticException::check_exception( + rel_error < targets.tolerance( term_uuid ), nullptr, + OpenGeodeException::TYPE::result, + "[StatisticsValidator] Failure \n --> term ", + model.term_name( term_uuid ), "\n mean = ", mean, + "\n target = ", target, "\n error = ", rel_error, + "\n tol = ", targets.tolerance( term_uuid ) ); } + } - template < typename ObjectType > - double quadratic_loss( const StatisticsTracker< ObjectType >& tracker, - const TargetStatistics< ObjectType >& targets ) - { - double loss = 0.0; - - for( const auto& term_uuid : targets.active_terms() ) - { - const auto diff = - tracker.mean( term_uuid ) - targets.value( term_uuid ); + template < typename ObjectType > + double quadratic_loss( const StatisticsTracker< ObjectType >& tracker, + const TargetStatistics< ObjectType >& targets ) + { + double loss = 0.0; - loss += diff * diff; - } + for( const auto& term_uuid : targets.active_terms() ) + { + const auto diff = + tracker.mean( term_uuid ) - targets.value( term_uuid ); - return loss; + loss += diff * diff; } - } // namespace statistics -} // namespace geode \ No newline at end of file + + return loss; + } +} // namespace geode::statistics diff --git a/include/geode/stochastic/models/energy_terms/energy_term_builder.hpp b/include/geode/stochastic/models/energy_terms/energy_term_builder.hpp index c3b1503..119f01c 100644 --- a/include/geode/stochastic/models/energy_terms/energy_term_builder.hpp +++ b/include/geode/stochastic/models/energy_terms/energy_term_builder.hpp @@ -30,6 +30,8 @@ #include #include +#include + #include #include #include @@ -108,9 +110,9 @@ namespace geode template < typename ObjectType, typename NewEnergyTypeConfig > std::unique_ptr< EnergyTerm< ObjectType > > build_energy_term_impl( - const NewEnergyTypeConfig&, - const ObjectSets< ObjectType >&, - const SpatialDomain< ObjectType::dim >& ) + const NewEnergyTypeConfig& /*unused*/, + const ObjectSets< ObjectType >& /*unused*/, + const SpatialDomain< ObjectType::dim >& /*unused*/ ) { static_assert( sizeof( NewEnergyTypeConfig ) == 0, "Unsupported EnergyTermConfig type" ); @@ -119,9 +121,9 @@ namespace geode template < typename ObjectType > std::unique_ptr< EnergyTerm< ObjectType > > build_energy_term_impl( - const std::monostate&, - const ObjectSets< ObjectType >&, - const SpatialDomain< ObjectType::dim >& ) + const std::monostate& /*unused*/, + const ObjectSets< ObjectType >& /*unused*/, + const SpatialDomain< ObjectType::dim >& /*unused*/ ) { throw OpenGeodeStochasticStochasticException{ nullptr, OpenGeodeException::TYPE::data, @@ -135,7 +137,7 @@ namespace geode const SpatialDomain< ObjectType::dim >& domain ) { return std::visit( - [&]( const auto& term_cfg ) + [&object_sets, &domain]( const auto& term_cfg ) -> std::unique_ptr< EnergyTerm< ObjectType > > { return build_energy_term_impl< ObjectType >( term_cfg, object_sets, domain ); diff --git a/include/geode/stochastic/models/energy_terms/pairwise_term.hpp b/include/geode/stochastic/models/energy_terms/pairwise_term.hpp index 747ab99..c4feabb 100644 --- a/include/geode/stochastic/models/energy_terms/pairwise_term.hpp +++ b/include/geode/stochastic/models/energy_terms/pairwise_term.hpp @@ -106,8 +106,8 @@ namespace geode { double sum = 0.0; this->for_each_object_in_sets( state, this->impacted_set_ids(), - [&]( const ObjectId& cur_obj_id ) { - sum += accumulate_interactions_with_neighbors( + [&sum, &state, this]( const ObjectId& cur_obj_id ) { + sum += this->accumulate_interactions_with_neighbors( cur_obj_id, state ); } ); return sum; @@ -119,10 +119,16 @@ namespace geode std::optional< ObjectId > exclude_id, const ObjectSets< ObjectType >& state ) const { - double sum = 0.0; + const auto impacted_set_it = + objectset_adjacency_map_.find( object_ref.set_id ); + if( impacted_set_it == objectset_adjacency_map_.end() ) + { + return 0.; + } const auto neighbors = state.neighbors( object_ref.object, - objectset_adjacency_map_.at( object_ref.set_id ), + impacted_set_it->second, interaction_->neighborhood_searching_distance(), exclude_id ); + double sum = 0.0; for( const auto& neigh_id : neighbors ) { ObjectRef< ObjectType > neigh_object{ @@ -147,11 +153,17 @@ namespace geode { return 0.; } - double sum = 0.0; + const auto impacted_set_it = + objectset_adjacency_map_.find( object_id.set_id ); + if( impacted_set_it == objectset_adjacency_map_.end() ) + { + return 0.; + } const auto neighbors = state.neighbors( cur_obj, - objectset_adjacency_map_.at( object_id.set_id ), + impacted_set_it->second, interaction_->neighborhood_searching_distance(), object_id ); ObjectRef< ObjectType > object_ref{ cur_obj, object_id.set_id }; + double sum = 0.0; for( const auto& neigh_id : neighbors ) { const auto& neigh_obj = state.get_object( neigh_id ); diff --git a/include/geode/stochastic/models/energy_terms/single_object_term.hpp b/include/geode/stochastic/models/energy_terms/single_object_term.hpp index ae00091..fab3286 100644 --- a/include/geode/stochastic/models/energy_terms/single_object_term.hpp +++ b/include/geode/stochastic/models/energy_terms/single_object_term.hpp @@ -95,10 +95,10 @@ namespace geode const ObjectSets< ObjectType >& state ) const override { double sum = 0.0; - this->for_each_object_in_sets( - state, this->impacted_set_ids(), [&]( const ObjectId& obj_id ) { + this->for_each_object_in_sets( state, this->impacted_set_ids(), + [&state, &sum, this]( const ObjectId& obj_id ) { const auto& obj = state.get_object( obj_id ); - sum += feature_->evaluate( obj, this->domain() ); + sum += this->feature_->evaluate( obj, this->domain() ); } ); return sum; } diff --git a/include/geode/stochastic/spatial/pairwise_interactions/pairwise_interactions_builder.hpp b/include/geode/stochastic/spatial/pairwise_interactions/pairwise_interactions_builder.hpp index 04a6fbf..4421e0f 100644 --- a/include/geode/stochastic/spatial/pairwise_interactions/pairwise_interactions_builder.hpp +++ b/include/geode/stochastic/spatial/pairwise_interactions/pairwise_interactions_builder.hpp @@ -21,6 +21,8 @@ * */ +#include + #include #include @@ -47,7 +49,8 @@ namespace geode template < typename ObjectType, typename NewInteractionConfig > std::unique_ptr< PairwiseInteraction< ObjectType > > - build_pairwise_interaction_impl( const NewInteractionConfig& ) + build_pairwise_interaction_impl( + const NewInteractionConfig& /*unused*/ ) { static_assert( sizeof( NewInteractionConfig ) == 0, "Unsupported PairwiseInteractionConfig type" ); @@ -56,7 +59,7 @@ namespace geode template < typename ObjectType > std::unique_ptr< PairwiseInteraction< ObjectType > > - build_pairwise_interaction_impl( const std::monostate& ) + build_pairwise_interaction_impl( const std::monostate& /*unused*/ ) { throw OpenGeodeStochasticStochasticException{ nullptr, OpenGeodeException::TYPE::data, @@ -68,7 +71,7 @@ namespace geode build_pairwise_interaction( const PairwiseInteractionConfig& cfg ) { return std::visit( - [&]( auto&& interaction_cfg ) + []( auto&& interaction_cfg ) -> std::unique_ptr< PairwiseInteraction< ObjectType > > { return build_pairwise_interaction_impl< ObjectType >( interaction_cfg ); diff --git a/include/geode/stochastic/spatial/single_object_features/single_object_feature_builder.hpp b/include/geode/stochastic/spatial/single_object_features/single_object_feature_builder.hpp index edd3c88..97c954d 100644 --- a/include/geode/stochastic/spatial/single_object_features/single_object_feature_builder.hpp +++ b/include/geode/stochastic/spatial/single_object_features/single_object_feature_builder.hpp @@ -24,6 +24,8 @@ #include +#include + #include #include @@ -52,22 +54,25 @@ namespace geode { throw OpenGeodeStochasticStochasticException{ nullptr, OpenGeodeException::TYPE::data, - "SegmentLengthInsideBoxFeature not valid for this ObjectType" }; + "[SingleObjectFeatureBuilder] SegmentLengthInsideBoxFeature " + "not valid for this ObjectType" }; } } template < typename ObjectType, typename NewObjectFeatureConfig > std::unique_ptr< SingleObjectFeature< ObjectType > > - build_single_object_feature_impl( const NewObjectFeatureConfig& ) + build_single_object_feature_impl( + const NewObjectFeatureConfig& /*unused*/ ) { static_assert( sizeof( NewObjectFeatureConfig ) == 0, - "Unsupported SingleObjectFeatureConfig type" ); + "[SingleObjectFeatureBuilder] Unsupported " + "SingleObjectFeatureConfig type" ); return nullptr; } template < typename ObjectType > std::unique_ptr< SingleObjectFeature< ObjectType > > - build_single_object_feature_impl( const std::monostate& ) + build_single_object_feature_impl( const std::monostate& /*unused*/ ) { throw OpenGeodeStochasticStochasticException{ nullptr, OpenGeodeException::TYPE::data, @@ -80,7 +85,7 @@ namespace geode build_single_object_feature( const SingleObjectFeatureConfig& cfg ) { return std::visit( - [&]( auto&& interaction_cfg ) + []( auto&& interaction_cfg ) -> std::unique_ptr< SingleObjectFeature< ObjectType > > { return build_single_object_feature_impl< ObjectType >( interaction_cfg ); diff --git a/src/geode/stochastic/models/energy_terms/energy_term_builder.cpp b/src/geode/stochastic/models/energy_terms/energy_term_builder.cpp index ba865cd..4f9e157 100644 --- a/src/geode/stochastic/models/energy_terms/energy_term_builder.cpp +++ b/src/geode/stochastic/models/energy_terms/energy_term_builder.cpp @@ -20,7 +20,6 @@ * SOFTWARE. * */ -#pragma once #include @@ -63,12 +62,12 @@ namespace namespace geode { - std::pair< std::vector< geode::uuid >, + std::pair< std::vector< uuid >, absl::flat_hash_map< uuid, std::vector< uuid > > > pairwise_builder_initialize_interactions_helper( const std::vector< std::pair< uuid, uuid > >& interacting_sets ) { - std::vector< geode::uuid > interacting_set_ids; + std::vector< uuid > interacting_set_ids; absl::flat_hash_map< uuid, std::vector< uuid > > objectset_adjacency_map; diff --git a/tests/stochastic/CMakeLists.txt b/tests/stochastic/CMakeLists.txt index b18c043..10cbb1e 100644 --- a/tests/stochastic/CMakeLists.txt +++ b/tests/stochastic/CMakeLists.txt @@ -176,13 +176,13 @@ add_geode_test( ${PROJECT_NAME}::stochastic ) -#add_geode_test( -# SOURCE "sampling/mcmc/test-mh-strauss.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/test-mh-poisson.cpp b/tests/stochastic/sampling/mcmc/test-mh-poisson.cpp index eabfe37..bdd5932 100644 --- a/tests/stochastic/sampling/mcmc/test-mh-poisson.cpp +++ b/tests/stochastic/sampling/mcmc/test-mh-poisson.cpp @@ -48,7 +48,8 @@ namespace : public geode::SimulationRunner< geode::Point2D > { public: - PoissonSimulationRunner( const geode::SpatialDomain< 2 >& domain ) + explicit PoissonSimulationRunner( + const geode::SpatialDomain< 2 >& domain ) : geode::SimulationRunner< geode::Point2D >( domain ) {}; void add_set_descriptor( const SetDescription& descriptor ) @@ -94,7 +95,7 @@ namespace geode::ModelConfig config; for( const auto& energy_desc : density_descriptors_ ) { - config.terms.push_back( energy_desc ); + config.terms.emplace_back( energy_desc ); } model_ = std::move( geode::build_model< geode::Point2D >( @@ -122,52 +123,7 @@ namespace create_target_statistics(); } - void check_statistics( const geode::StatisticsTracker< geode::Point2D >& - 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() } ) - // { - // 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::fabs( computed_means[stat_id] - - // expected_means ) / expected_means; - // - // geode::OpenGeodeStochasticStochasticException::test( - // target_vs_mean_error < 0.05, "[MH test] - // statistic value ", computed_means[stat_id], " - // for energy term: ", term.name().value_or( - // term.id().string() ), " not close enough to - // expected value ", expected_means, " --> error - // : ", target_vs_mean_error ); - // - // const auto target_vs_variance_error = - // std::fabs( computed_variances[stat_id] - - // expected_means ) / expected_means; - // - // geode::OpenGeodeStochasticStochasticException::test( - // target_vs_variance_error < 0.15, - // "[MH test] variance of statistic ", - // computed_variances[stat_id], " for energy - // term: ", term.name().value_or( - // term.id().string() ), " not close enough to - // expected value ", expected_means, " --> error - // : ", target_vs_variance_error ); - // } - } - private: - geode::BoundingBox2D box_; std::vector< SetDescription > set_descriptors_; std::vector< PoissonDensityDescription > density_descriptors_; std::vector< geode::TargetStatisticConfig > @@ -181,6 +137,7 @@ namespace geode::RandomEngine engine; engine.set_seed( "@mh-test-single-POISSON@" ); + // NOLINTBEGIN(*-magic-numbers) geode::BoundingBox2D box; box.add_point( geode::Point2D{ { 0.0, 0.0 } } ); box.add_point( geode::Point2D{ { 10.0, 10.0 } } ); @@ -188,29 +145,28 @@ namespace std::array< double, 4 > birth_ratio{ 0.1, 0.5, 2., 4. }; std::array< double, 4 > change_ratio{ 0., 1., 1., 0. }; - for( const auto config : geode::Range{ birth_ratio.size() } ) { // --- Set description - SetDescription setA; - setA.name = "A"; - setA.birth_ratio = birth_ratio[config]; - setA.death_ratio = 1.0; - setA.change_ratio = change_ratio[config]; + SetDescription set_a; + set_a.name = "A"; + set_a.birth_ratio = birth_ratio[config]; + set_a.death_ratio = 1.0; + set_a.change_ratio = change_ratio[config]; // --- Energy term description - PoissonDensityDescription densityA; - densityA.term_name = "density"; - densityA.object_set_names = { "A" }; - densityA.lambda = 0.3; - densityA.object_feature = geode::ObjectInDomainFeatureConfig{}; + PoissonDensityDescription density_a; + density_a.term_name = "density"; + density_a.object_set_names = { "A" }; + density_a.lambda = 0.3; + density_a.object_feature = geode::ObjectInDomainFeatureConfig{}; - geode::TargetStatisticConfig statA{ "density", 30.0, 0.15 }; + geode::TargetStatisticConfig stat_a{ "density", 30.0, 0.15 }; PoissonSimulationRunner runner( domain ); - runner.add_set_descriptor( setA ); - runner.add_density_descriptor( densityA ); - runner.add_target_statistics( statA ); + runner.add_set_descriptor( set_a ); + runner.add_density_descriptor( density_a ); + runner.add_target_statistics( stat_a ); runner.initialize(); // run simulation @@ -229,6 +185,7 @@ namespace geode::statistics::validate( statistic_tracker, runner.target_statistics() ); } + // NOLINTEND(*-magic-numbers) geode::Logger::info( "--> SUCCESS!" ); } @@ -239,6 +196,7 @@ namespace geode::RandomEngine engine; engine.set_seed( "@mh-test-POISSON-multi@" ); + // NOLINTBEGIN(*-magic-numbers) geode::BoundingBox2D box; box.add_point( geode::Point2D{ { 0.0, 0.0 } } ); @@ -300,6 +258,7 @@ namespace sim_config.metropolis_hasting_steps = 1000; sim_config.burn_in_steps = 3000; sim_config.printer = printer_config; + // NOLINTEND(*-magic-numbers) auto statistic_tracker = runner.run( engine, sim_config ); geode::statistics::validate( diff --git a/tests/stochastic/sampling/mcmc/test-mh-strauss.cpp b/tests/stochastic/sampling/mcmc/test-mh-strauss.cpp index 61cc11d..5f09044 100644 --- a/tests/stochastic/sampling/mcmc/test-mh-strauss.cpp +++ b/tests/stochastic/sampling/mcmc/test-mh-strauss.cpp @@ -21,15 +21,16 @@ * */ #include -#include -#include +#include +#include #include -#include -#include -#include #include +#include +#include +#include #include -#include +// #include +#include namespace { @@ -41,27 +42,15 @@ namespace double change_ratio{ 1.0 }; }; - struct PoissonDensityDescription - { - std::string name; - double density; - double target_count; - }; - - struct PairwiseInteractionDescription - { - std::vector< std::string > names; - double strength; - double distance_threshold; - // geode::PairwiseInteraction::SCOPE interaction_scope; - double target_interaction_count; - }; + using PoissonDensityDescription = geode::SingleObjectTermConfig; + using PairwiseInteractionDescription = geode::PairwiseTermConfig; class StraussSimulationRunner : public geode::SimulationRunner< geode::Point2D > { public: - StraussSimulationRunner( const geode::SpatialDomain< 2 >& domain ) + explicit StraussSimulationRunner( + const geode::SpatialDomain< 2 >& domain ) : geode::SimulationRunner< geode::Point2D >( domain ) { } @@ -83,113 +72,75 @@ namespace interaction_descriptors_.push_back( descriptor ); } - void initialize() override + void add_target_statistics( + const geode::TargetStatisticConfig& statistic_descriptor ) + { + targeted_statistics_descriptors_.push_back( statistic_descriptor ); + } + std::unique_ptr< geode::ProposalKernel< geode::Point2D > > + create_sets_and_set_samplers() { 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_ ) { 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 > >( - this->domain_ ) ); + domain_ ) ); 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 ); } + return proposal_kernel; + } - // Step 2: create density energy terms - for( const auto& density_desc : density_descriptors_ ) + void create_model() + { + geode::ModelConfig config; + for( const auto& energy_desc : density_descriptors_ ) { - 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.name, "_density" ), - density_desc.density, - std::vector< geode::uuid >{ set_id }, - this->domain_ ) ) ); - - this->ordered_target_statistics_.push_back( - density_desc.target_count ); + config.terms.emplace_back( energy_desc ); } - - // Step 3: create pairwise interaction terms for( const auto& interaction_desc : interaction_descriptors_ ) { - std::vector< geode::uuid > set_ids; - for( const auto& name : interaction_desc.names ) - { - set_ids.emplace_back( name_to_uuid.at( name ) ); - } - - auto interaction = std::make_unique< - geode::MinimalDistanceCutoff< 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.names, "_" ), - "_interaction" ), - interaction_desc.strength, set_ids, - std::move( interaction ), this->domain_ ) ) ); - - this->ordered_target_statistics_.push_back( - interaction_desc.target_interaction_count ); + config.terms.emplace_back( interaction_desc ); } - this->mh_sampler_ = - std::make_unique< geode::MetropolisHastings< geode::Point2D > >( - this->energy_terms_collection_, - std::move( proposal_kernel ) ); + model_ = std::move( geode::build_model< geode::Point2D >( + config, object_sets_, domain_ ) ); + create_target_statistics(); } - void check_statistics( - const geode::StatisticsTracker& statistic_monitoring ) const + void create_target_statistics() { - const auto& computed_means = statistic_monitoring.means(); - - for( const auto stat_id : - geode::Range{ this->energy_terms_collection_.size() } ) + target_statistics_.emplace( *model_ ); + for( const auto& target_stat : targeted_statistics_descriptors_ ) { - 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::fabs( computed_means[stat_id] - expected_mean ); - if( expected_mean > 0 ) - { - target_vs_mean_error /= expected_mean; - } - - geode::OpenGeodeStochasticStochasticException::test( - target_vs_mean_error < 0.1, "[MH test] Statistic value ", - computed_means[stat_id], " for energy term: ", - term.name().value_or( term.id().string() ), - " not close enough to expected value ", expected_mean, - " --> error: ", target_vs_mean_error ); + target_statistics_->set_target( target_stat ); } } + void initialize() override + { + auto proposal_kernel = create_sets_and_set_samplers(); + create_model(); + + this->mh_sampler_ = + std::make_unique< geode::MetropolisHastings< geode::Point2D > >( + *model_, std::move( proposal_kernel ) ); + create_target_statistics(); + } + private: std::vector< SetDescription > set_descriptors_; std::vector< PoissonDensityDescription > density_descriptors_; std::vector< PairwiseInteractionDescription > interaction_descriptors_; + + std::vector< geode::TargetStatisticConfig > + targeted_statistics_descriptors_; }; void test_single_type_strauss() @@ -199,46 +150,54 @@ namespace geode::RandomEngine engine; engine.set_seed( "@mh-test-single-STRAUSS@" ); - + // NOLINTBEGIN(*-magic-numbers) geode::BoundingBox2D box; box.add_point( geode::Point2D{ { 0.0, 0.0 } } ); box.add_point( geode::Point2D{ { 10.0, 10.0 } } ); - // todo change!! geode::SpatialDomain domain( box, 1. ); std::array< double, 5 > gamma_values{ 0, 0.3, 0.5, 0.7, 1.0 }; std::array< double, 5 > nb_points{ 19.5, 24.4, 31.3, 36.1, 50. }; - std::array< double, 5 > nb_interactions{ 0, 4, 8, 15, 43 }; - + std::array< double, 5 > nb_interactions{ 0, 4.7, 9.8, 18.7, 50.3 }; for( const auto config : geode::Range{ gamma_values.size() } ) { // --- 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.name = "A"; - densityA.density = 0.5; - densityA.target_count = nb_points[config]; + SetDescription set_a; + set_a.name = "A"; + set_a.birth_ratio = 1.0; + set_a.death_ratio = 1.0; + set_a.change_ratio = 1.0; + + // --- Energy term description + PoissonDensityDescription density_a; + density_a.term_name = "density_a"; + density_a.object_set_names = { "A" }; + density_a.lambda = 0.5; + density_a.object_feature = geode::ObjectInDomainFeatureConfig{}; + + geode::TargetStatisticConfig stat_a{ "density_a", nb_points[config], + 0.1 }; // --- Intra-set pairwise interaction (Strauss process) - PairwiseInteractionDescription intraA; - intraA.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]; + PairwiseInteractionDescription interaction_a; + interaction_a.term_name = "interactionA"; + interaction_a.object_set_names_interactions = { { "A", "A" } }; + interaction_a.gamma = gamma_values[config]; - StraussSimulationRunner runner( domain ); - runner.add_set_descriptor( setA ); - runner.add_density_descriptor( densityA ); - runner.add_interaction_descriptor( intraA ); + interaction_a.interaction_config = + geode::MinimalDistanceCutoffConfig{ + 1. + }; /// ici cela devrait etre un paramètre utilisateur + + geode::TargetStatisticConfig stat_intra_a{ "interactionA", + nb_interactions[config], 0.1 }; + StraussSimulationRunner runner( domain ); + runner.add_set_descriptor( set_a ); + runner.add_density_descriptor( density_a ); + runner.add_interaction_descriptor( interaction_a ); + runner.add_target_statistics( stat_a ); + runner.add_target_statistics( stat_intra_a ); runner.initialize(); // run simulation @@ -253,10 +212,11 @@ namespace sim_config.burn_in_steps = 1000; sim_config.printer = printer_config; - auto statistic_monitoring = runner.run( engine, sim_config ); - runner.check_statistics( statistic_monitoring ); + auto statistic_tracker = runner.run( engine, sim_config ); + geode::statistics::validate( + statistic_tracker, runner.target_statistics() ); } - + // NOLINTEND(*-magic-numbers) geode::Logger::info( "--> SUCCESS!" ); } @@ -267,20 +227,18 @@ namespace geode::RandomEngine engine; engine.set_seed( "@mh-test-multi-STRAUSS@" ); - + // NOLINTBEGIN(*-magic-numbers) geode::BoundingBox2D box; box.add_point( geode::Point2D{ { 0.0, 0.0 } } ); box.add_point( geode::Point2D{ { 10.0, 10.0 } } ); - // todo change!! - geode::SpatialDomain domain( box, 0. ); + geode::SpatialDomain domain( box, 2. ); 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{ 8, 20, 85 }; - + std::array< double, 3 > nb_points01{ 6.7, 8, 10.0 }; + std::array< double, 3 > nb_points02{ 17.5, 24.6, 40.0 }; + std::array< double, 3 > nb_points03{ 14.6, 19.4, 30. }; + std::array< double, 3 > nb_interactions01{ 0, 15, 59.8 }; + std::array< double, 3 > nb_interactions02{ 37.2, 70, 174 }; for( const auto config : geode::Range{ gamma_values.size() } ) { // --- Sets @@ -289,17 +247,59 @@ namespace 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] }; + PoissonDensityDescription d01; + d01.term_name = "density_set01"; + d01.object_set_names = { "set01" }; + d01.lambda = 0.1; + d01.object_feature = geode::ObjectInDomainFeatureConfig{}; + + geode::TargetStatisticConfig stat01{ "density_set01", + nb_points01[config], 0.1 }; + + PoissonDensityDescription d02; + d02.term_name = "density_set02"; + d02.object_set_names = { "set02" }; + d02.lambda = 0.4; + d02.object_feature = geode::ObjectInDomainFeatureConfig{}; + + geode::TargetStatisticConfig stat02{ "density_set02", + nb_points02[config], 0.1 }; + + PoissonDensityDescription d03; + d03.term_name = "density_set03"; + d03.object_set_names = { "set03" }; + d03.lambda = 0.3; + d03.object_feature = geode::ObjectInDomainFeatureConfig{}; + + geode::TargetStatisticConfig stat03{ "density_set03", + nb_points03[config], 0.1 }; // --- 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] }; + PairwiseInteractionDescription intra01; + intra01.term_name = "interaction01"; + intra01.object_set_names_interactions = { { "set01", "set01" }, + { "set02", "set02" }, { "set03", "set03" } }; + intra01.gamma = gamma_values[config]; + + intra01.interaction_config = geode::MinimalDistanceCutoffConfig{ + 1. + }; /// ici cela devrait etre un paramètre utilisateur + + geode::TargetStatisticConfig stat_intra_01{ "interaction01", + nb_interactions01[config], 0.1 }; + + PairwiseInteractionDescription intra02; + intra02.term_name = "interaction02"; + intra02.object_set_names_interactions = { { "set02", "set02" } }; + intra02.gamma = 1.; // gamma_values[config]; + + intra02.interaction_config = geode::MinimalDistanceCutoffConfig{ + 2. + }; /// ici cela devrait etre un paramètre utilisateur + + geode::TargetStatisticConfig stat_intra_02{ "interaction02", + nb_interactions02[config], 0.1 }; StraussSimulationRunner runner( domain ); runner.add_set_descriptor( set01 ); @@ -313,6 +313,12 @@ namespace runner.add_interaction_descriptor( intra01 ); runner.add_interaction_descriptor( intra02 ); + runner.add_target_statistics( stat01 ); + runner.add_target_statistics( stat02 ); + runner.add_target_statistics( stat03 ); + runner.add_target_statistics( stat_intra_01 ); + runner.add_target_statistics( stat_intra_02 ); + runner.initialize(); // run simulation @@ -322,14 +328,16 @@ namespace "/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.realizations = 600; + sim_config.metropolis_hasting_steps = 750; + sim_config.burn_in_steps = 5000; sim_config.printer = printer_config; - auto statistic_monitoring = runner.run( engine, sim_config ); - runner.check_statistics( statistic_monitoring ); + auto statistic_tracker = runner.run( engine, sim_config ); + geode::statistics::validate( + statistic_tracker, runner.target_statistics() ); } + // NOLINTEND(*-magic-numbers) geode::Logger::info( "--> SUCCESS!" ); }