diff --git a/bindings/python/src/stochastic/CMakeLists.txt b/bindings/python/src/stochastic/CMakeLists.txt index 48bd7a0..a79ce85 100644 --- a/bindings/python/src/stochastic/CMakeLists.txt +++ b/bindings/python/src/stochastic/CMakeLists.txt @@ -21,7 +21,7 @@ add_geode_python_binding( NAME "py_stochastic" SOURCES - "sampling/mcmc/helpers/fracture_simulation_runner.hpp" + # "sampling/mcmc/helpers/fracture_simulation_runner.hpp" "sampling/mcmc/helpers/simulation_monitor.hpp" "sampling/mcmc/helpers/simulation_printer.hpp" "sampling/mcmc/simulation_runner.hpp" diff --git a/bindings/python/src/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp b/bindings/python/src/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp index dbb6e98..9766162 100644 --- a/bindings/python/src/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp +++ b/bindings/python/src/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp @@ -88,7 +88,7 @@ namespace geode // pybind11::arg( "engine" ), pybind11::arg( "steps" // ), "Run simulation for a fixed number of steps." ) .def( "run", - static_cast< StatisticsMonitor ( FractureSimulationRunner::* )( + static_cast< StatisticsTracker ( FractureSimulationRunner::* )( RandomEngine&, const SimulationConfigurator& ) >( &FractureSimulationRunner::run ), pybind11::arg( "engine" ), pybind11::arg( "config" ), diff --git a/bindings/python/src/stochastic/sampling/mcmc/helpers/simulation_monitor.hpp b/bindings/python/src/stochastic/sampling/mcmc/helpers/simulation_monitor.hpp index 3a18b60..9b2d960 100644 --- a/bindings/python/src/stochastic/sampling/mcmc/helpers/simulation_monitor.hpp +++ b/bindings/python/src/stochastic/sampling/mcmc/helpers/simulation_monitor.hpp @@ -21,30 +21,30 @@ * */ -#include +#include namespace geode { void define_simulation_monitor( pybind11::module &module ) { - pybind11::class_< geode::StatisticsMonitor >( - module, "StatisticsMonitor" ) + pybind11::class_< geode::StatisticsTracker >( + module, "StatisticsTracker" ) .def( pybind11::init< geode::index_t >(), pybind11::arg( "nb_energy_terms" ), - "Create a StatisticsMonitor for a given number of energy " + "Create a StatisticsTracker for a given number of energy " "terms" ) - .def( "add_realization", &geode::StatisticsMonitor::add_realization, + .def( "add_realization", &geode::StatisticsTracker::add_realization, pybind11::arg( "values" ), "Add a realization (vector of doubles) to update statistics" ) - .def( "statiscal_count", &geode::StatisticsMonitor::statiscal_count, + .def( "statiscal_count", &geode::StatisticsTracker::statiscal_count, "Return the number of realizations added" ) - .def_property_readonly( "means", &geode::StatisticsMonitor::means, + .def_property_readonly( "means", &geode::StatisticsTracker::means, "Return the computed mean values for each energy term" ) .def_property_readonly( "variances", - &geode::StatisticsMonitor::variances, + &geode::StatisticsTracker::variances, "Return the computed variances for each energy term" ) - .def( "__repr__", []( const geode::StatisticsMonitor &self ) { - return ""; } ); } diff --git a/bindings/python/src/stochastic/sampling/mcmc/helpers/simulation_printer.hpp b/bindings/python/src/stochastic/sampling/mcmc/helpers/simulation_printer.hpp index a1d36d1..ca8c889 100644 --- a/bindings/python/src/stochastic/sampling/mcmc/helpers/simulation_printer.hpp +++ b/bindings/python/src/stochastic/sampling/mcmc/helpers/simulation_printer.hpp @@ -64,7 +64,7 @@ namespace geode // &SimulationPrinter::print_statistics_summary, // pybind11::arg( "monitor" ), // pybind11::arg( "energy_term_names" ) = "", - // "Print statistics summary from a StatisticsMonitor." + // "Print statistics summary from a StatisticsTracker." // ); } } // namespace geode \ No newline at end of file diff --git a/bindings/python/src/stochastic/stochastic.cpp b/bindings/python/src/stochastic/stochastic.cpp index c1a7516..374fa86 100644 --- a/bindings/python/src/stochastic/stochastic.cpp +++ b/bindings/python/src/stochastic/stochastic.cpp @@ -26,8 +26,8 @@ #include "sampling/direct/double_sampler.hpp" -#include "sampling/mcmc/helpers/fracture_simulation_runner.hpp" -#include "sampling/mcmc/helpers/simulation_monitor.hpp" +// #include "sampling/mcmc/helpers/fracture_simulation_runner.hpp" +// #include "sampling/mcmc/helpers/simulation_monitor.hpp" #include "sampling/mcmc/helpers/simulation_printer.hpp" #include "sampling/mcmc/simulation_runner.hpp" @@ -49,8 +49,8 @@ PYBIND11_MODULE( opengeode_stochastic_py_stochastic, module ) geode::define_random_engine( module ); geode::define_double_sampler( module ); - geode::define_simulation_monitor( module ); + // geode::define_simulation_monitor( module ); geode::define_simulation_printer( module ); geode::define_simulation_runner( module ); - geode::define_fracture_simulation( module ); + // geode::define_fracture_simulation( module ); } \ No newline at end of file diff --git a/bindings/python/tests/stochastic/CMakeLists.txt b/bindings/python/tests/stochastic/CMakeLists.txt index ded7457..6e00e74 100644 --- a/bindings/python/tests/stochastic/CMakeLists.txt +++ b/bindings/python/tests/stochastic/CMakeLists.txt @@ -18,8 +18,8 @@ # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE # SOFTWARE. -add_geode_python_test( - SOURCE "test-py-mh-fractures.py" - DEPENDENCIES - ${PROJECT_NAME}::py_stochastic -) \ No newline at end of file +#add_geode_python_test( +# SOURCE "test-py-mh-fractures.py" +# DEPENDENCIES +# ${PROJECT_NAME}::py_stochastic +#) \ No newline at end of file diff --git a/include/geode/stochastic/applications/fractures.hpp b/include/geode/stochastic/applications/fractures.hpp new file mode 100644 index 0000000..b10a8e2 --- /dev/null +++ b/include/geode/stochastic/applications/fractures.hpp @@ -0,0 +1,163 @@ +/* + * Copyright (c) 2019 - 2026 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 + +namespace geode +{ + using Fracture = OwnerSegment2D; + using FractureSamplerConfig = ObjectSamplerConfig< Fracture >; + using FractureSimulationContext = SimulationContext< Fracture >; + using FractureSimulationRunner = SimulationRunner< Fracture >; + + struct FractureSetDescription + { + std::string fset_name; + + FractureSamplerConfig sampler; + double birth_ratio{ 1.0 }; + double death_ratio{ 1.0 }; + double change_ratio{ 1.0 }; + + [[nodiscard]] std::string density_name() const + { + return absl::StrCat( fset_name, "_p20" ); + } + double p20{ 0. }; + std::optional< double > expected_number; + + [[nodiscard]] std::string intensity_name() const + { + return absl::StrCat( fset_name, "_p21" ); + } + double p21{ 0. }; + std::optional< double > expected_total_length; + + std::vector< std::array< geode::Point2D, 2 > > observed_fractures; + + [[nodiscard]] std::string spacing_name() const + { + return absl::StrCat( fset_name, "_spacing" ); + } + double minimal_spacing{ 0. }; + + [[nodiscard]] std::string string() const + { + auto message = + absl::StrCat( "FractureSetDescription: ", fset_name ); + for( const auto& fixed_object : observed_fractures ) + { + absl::StrAppend( &message, + "\n\t --> observation (x,y,z)start: ", + fixed_object[0].string(), + " (x,y,z)end: ", fixed_object[1].string() ); + } + absl::StrAppend( &message, + "\n\t --> length distribution: ", sampler.length.string() ); + absl::StrAppend( &message, + "\n\t --> azimuth distribution: ", sampler.azimuth.string() ); + absl::StrAppend( &message, "\n\t --> ", density_name(), ": ", p20 ); + absl::StrAppend( + &message, "\n\t --> ", intensity_name(), ": ", p21 ); + absl::StrAppend( + &message, "\n\t --> ", spacing_name(), ": ", minimal_spacing ); + + absl::StrAppend( &message, + "\n\t --> dynamic move ratio - birth/death/change (", + birth_ratio, " / ", death_ratio, " / ", change_ratio, ")" ); + return message; + } + }; + + // struct FractureInterSetDescription + // { + // std::string interaction_name; + // + // std::vector< std::string > set_names; + // + // double gamma{ 1. }; + // double distance{ 0. }; + // + // bool include_intra_set{ true }; + // bool include_inter_set{ false }; + // + // std::optional< double > expected_nb_interactions; + // }; + + struct opengeode_stochastic_stochastic_api FractureNetworkDescription + { + std::string fnet_name; + + SpatialDomainConfig< 2 > domain; + + std::vector< FractureSetDescription > fracture_sets; + + [[nodiscard]] FractureSetDescription& add_fracture_set( + absl::string_view fset_name ) + { + auto& fracture_set = fracture_sets.emplace_back(); + fracture_set.fset_name = fset_name; + return fracture_set; + } + + [[nodiscard]] std::string string() const + { + auto message = + absl::StrCat( "FractureNetworkDescription: ", fnet_name ); + absl::StrAppend( &message, "\n\t --> ", domain.string() ); + for( const auto& fset_desc : fracture_sets ) + { + absl::StrAppend( &message, "\n\t --> ", fset_desc.string() ); + } + return message; + } + // std::vector< StraussInteractionDescription< ObjectType > > + // inter_set_interactions; + // + // void add_x_node_monitoring( double beta_x_node ) + // { + // OpenGeodeStochasticStochasticException::check_exception( + // beta_x_node <= 1.0 && beta_x_node >= 0., nullptr, + // OpenGeodeException::TYPE::data, + // "[FractureSimulationRunner] x node should be + // inhibitated, " "please provise a value in [0., 1.]." + // ); + // beta_x_node_ = beta_x_node; + // } + }; + + opengeode_stochastic_stochastic_api FractureSimulationContext + build_fractures_simulation_context( + const FractureNetworkDescription& description ); + + opengeode_stochastic_stochastic_api + std::vector< geode::TargetStatisticConfig > + build_fractures_targeted_stat( + const FractureNetworkDescription& description ); + +} // namespace geode \ No newline at end of file diff --git a/include/geode/stochastic/applications/poisson_process.hpp b/include/geode/stochastic/applications/poisson_process.hpp new file mode 100644 index 0000000..acf38ec --- /dev/null +++ b/include/geode/stochastic/applications/poisson_process.hpp @@ -0,0 +1,73 @@ +/* + * Copyright (c) 2019 - 2026 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 +{ + template < typename ObjectType > + struct PoissonSetDescription + { + std::string set_name; + + ObjectSamplerConfig< ObjectType > sampler; + + std::string density_name; + double lambda{ 0. }; + std::optional< double > expected_nb_objects; + + double birth_ratio{ 1.0 }; + double death_ratio{ 1.0 }; + double change_ratio{ 1.0 }; + }; + + template < typename ObjectType > + struct PoissonProcessDescription + { + std::string process_name{ "Poisson" }; + SpatialDomainConfig< ObjectType::dim > domain; + + std::vector< PoissonSetDescription< ObjectType > > sets; + + PoissonSetDescription< ObjectType >& add_set( + absl::string_view set_name ) + { + auto& set = sets.emplace_back(); + set.set_name = set_name; + set.density_name = absl::StrCat( set_name, "_density" ); + return set; + } + }; + + template < typename ObjectType > + SimulationContext< ObjectType > build_poisson_process( + const PoissonProcessDescription< ObjectType >& description ); + + template < typename ObjectType > + std::vector< geode::TargetStatisticConfig > build_poisson_targeted_stat( + const PoissonProcessDescription< ObjectType >& description ); + +} // namespace geode \ No newline at end of file diff --git a/include/geode/stochastic/applications/strauss_process.hpp b/include/geode/stochastic/applications/strauss_process.hpp new file mode 100644 index 0000000..50c033a --- /dev/null +++ b/include/geode/stochastic/applications/strauss_process.hpp @@ -0,0 +1,84 @@ +/* + * Copyright (c) 2019 - 2026 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 + +namespace geode +{ + template < typename ObjectType > + struct StraussInteractionDescription + { + std::string interaction_name; + + std::vector< std::string > set_names; + + double gamma{ 1. }; + double distance{ 0. }; + + bool include_intra_set{ true }; + bool include_inter_set{ false }; + + std::optional< double > expected_nb_interactions; + }; + + template < typename ObjectType > + struct StraussProcessDescription + { + SpatialDomainConfig< ObjectType::dim > domain; + + std::vector< PoissonSetDescription< ObjectType > > sets; + + std::vector< StraussInteractionDescription< ObjectType > > interactions; + + PoissonSetDescription< ObjectType >& add_set( + absl::string_view set_name ) + { + auto& set = sets.emplace_back(); + set.set_name = set_name; + set.density_name = absl::StrCat( set_name, "_density" ); + + return set; + } + + StraussInteractionDescription< ObjectType >& add_interaction( + absl::string_view interaction_name ) + { + auto& interaction = interactions.emplace_back(); + interaction.interaction_name = interaction_name; + return interaction; + } + }; + + template < typename ObjectType > + SimulationContext< ObjectType > build_strauss_process( + const StraussProcessDescription< ObjectType >& description ); + + template < typename ObjectType > + std::vector< geode::TargetStatisticConfig > build_strauss_targeted_stat( + const StraussProcessDescription< ObjectType >& description ); + +} // namespace geode \ No newline at end of file diff --git a/include/geode/stochastic/inference/statistics_tools.hpp b/include/geode/stochastic/inference/statistics_tools.hpp new file mode 100644 index 0000000..4a426be --- /dev/null +++ b/include/geode/stochastic/inference/statistics_tools.hpp @@ -0,0 +1,73 @@ +/* + * Copyright (c) 2019 - 2026 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 + +namespace geode::statistics +{ + template < typename ObjectType > + void validate( const StatisticsTracker< ObjectType >& tracker, + const TargetStatistics< ObjectType >& targets ) + { + 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 ); + + 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 \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 ); + + loss += diff * diff; + } + + return loss; + } +} // namespace geode::statistics diff --git a/include/geode/stochastic/inference/statistics_tracker.hpp b/include/geode/stochastic/inference/statistics_tracker.hpp new file mode 100644 index 0000000..be9cb32 --- /dev/null +++ b/include/geode/stochastic/inference/statistics_tracker.hpp @@ -0,0 +1,108 @@ +/* + * Copyright (c) 2019 - 2026 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 + +namespace geode +{ + template < typename ObjectType > + class StatisticsTracker + { + public: + StatisticsTracker( const Model< ObjectType >& model ) : model_{ model } + { + means_.resize( model.nb_terms(), 0.0 ); + m2_.resize( model.nb_terms(), 0.0 ); + } + + [[nodiscard]] index_t statiscal_count() const + { + return count_; + } + + void add_realization( const std::vector< double >& values ) + { + ++count_; + for( const auto value_id : geode::Range{ values.size() } ) + { + auto& value = values[value_id]; + auto& mean = means_[value_id]; + auto& sum_of_squares = m2_[value_id]; + + const auto delta = value - mean; + mean += delta / count_; + const auto delta2 = value - mean; + sum_of_squares += delta * delta2; + } + } + + [[nodiscard]] double mean( const uuid& term_uuid ) const + { + return means_[model_.term_index( term_uuid )]; + } + + [[nodiscard]] const std::vector< double >& means() const + { + return means_; + } + + [[nodiscard]] double variance( const uuid& term_uuid ) const + { + return variance( model_.term_index( term_uuid ) ); + } + + [[nodiscard]] std::vector< double > variances() const + { + std::vector< double > variances; + variances.reserve( model_.nb_terms() ); + for( const auto variance_id : geode::Range{ model_.nb_terms() } ) + { + variances.emplace_back( this->variance( variance_id ) ); + } + return variances; + } + + private: + [[nodiscard]] double variance( index_t term_index ) const + { + if( count_ < 2 ) + { + return 0.0; + } + return m2_[term_index] / ( count_ - 1 ); + } + + private: + const Model< ObjectType >& model_; + + std::vector< double > means_; + std::vector< double > m2_; + index_t count_{ 0 }; + }; +} // namespace geode \ No newline at end of file diff --git a/include/geode/stochastic/inference/target_statistics.hpp b/include/geode/stochastic/inference/target_statistics.hpp new file mode 100644 index 0000000..703f28d --- /dev/null +++ b/include/geode/stochastic/inference/target_statistics.hpp @@ -0,0 +1,111 @@ +/* + * Copyright (c) 2019 - 2026 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 +{ + constexpr double STATISTIC_TOLERANCE_VALUE{ 0.1 }; + + struct TargetStatisticConfig + { + std::string term_name; + double value; + double tolerance{ STATISTIC_TOLERANCE_VALUE }; + }; + + template < typename ObjectType > + class TargetStatistics + { + public: + explicit TargetStatistics( const Model< ObjectType >& model, + const std::vector< TargetStatisticConfig >& statistic_targets ) + : model_( model ) + { + values_.resize( model.nb_terms(), 0.0 ); + tolerances_.resize( model.nb_terms(), 0.0 ); + active_.resize( model.nb_terms(), false ); + for( const auto& target : statistic_targets ) + { + set_target( target ); + } + } + + [[nodiscard]] const Model< ObjectType >& model() const + { + return model_; + } + + [[nodiscard]] bool has_target( const uuid& term_uuid ) const + { + return active_[model_.term_index( term_uuid )]; + } + + [[nodiscard]] double target( const uuid& term_uuid ) const + { + return values_[model_.term_index( term_uuid )]; + } + + [[nodiscard]] double tolerance( const uuid& term_uuid ) const + { + return tolerances_[model_.term_index( term_uuid )]; + } + + [[nodiscard]] std::vector< uuid > active_terms() const + { + std::vector< uuid > active_terms_uuid; + + for( const auto& term : model_.terms().energy_terms() ) + { + const auto& term_id = term->id(); + if( active_[model_.term_index( term_id )] ) + { + active_terms_uuid.push_back( term_id ); + } + } + return active_terms_uuid; + } + + private: + void set_target( const TargetStatisticConfig& statistic ) + { + const auto term_uuid = + model_.terms().get_term_uuid( statistic.term_name ); + const auto idx = model_.term_index( term_uuid ); + + values_[idx] = statistic.value; + tolerances_[idx] = statistic.tolerance; + active_[idx] = true; + } + + private: + const Model< ObjectType >& model_; + + std::vector< double > values_; + std::vector< double > tolerances_; + std::vector< bool > active_; + }; +} // namespace geode diff --git a/include/geode/stochastic/models/energy_term_collection.hpp b/include/geode/stochastic/models/energy_term_collection.hpp new file mode 100644 index 0000000..d333259 --- /dev/null +++ b/include/geode/stochastic/models/energy_term_collection.hpp @@ -0,0 +1,131 @@ +/* + * Copyright (c) 2019 - 2026 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 +{ + template < typename ObjectType > + class EnergyTermCollection + { + OPENGEODE_DISABLE_COPY( EnergyTermCollection ); + + public: + EnergyTermCollection() = default; + EnergyTermCollection( EnergyTermCollection&& ) noexcept = default; + + ~EnergyTermCollection() = default; + + EnergyTermCollection& operator=( EnergyTermCollection&& ) = default; + + uuid add_energy_term( + std::unique_ptr< EnergyTerm< ObjectType > >&& term ) + { + auto term_idx = energy_terms_.size(); + + const auto term_name = term->name(); + OpenGeodeStochasticStochasticException::check_exception( + term_name.has_value(), nullptr, OpenGeodeException::TYPE::data, + absl::StrCat( "[EnergyTermCollection] Energy Term name is not " + "defined." ) ); + const auto term_uuid = term->id(); + auto [it, inserted_uuid] = + name_to_uuid_.emplace( term_name.value(), term_uuid ); + OpenGeodeStochasticStochasticException::check_exception( + inserted_uuid, nullptr, OpenGeodeException::TYPE::data, + absl::StrCat( "[EnergyTermCollection] Energy Term named ", + term_name.value(), " already exists." ) ); + + auto [it2, inserted_index] = + uuid_to_index_.emplace( term_uuid, term_idx ); + OpenGeodeStochasticStochasticException::check_exception( + inserted_index, nullptr, OpenGeodeException::TYPE::data, + absl::StrCat( "[EnergyTermCollection] Energy Term ", + term_uuid.string(), " already exists." ) ); + + energy_terms_.emplace_back( std::move( term ) ); + return term_uuid; + } + + [[nodiscard]] index_t size() const + { + return energy_terms_.size(); + } + + [[nodiscard]] index_t get_term_index( const uuid& term_uuid ) const + { + auto term_it = uuid_to_index_.find( term_uuid ); + OpenGeodeStochasticStochasticException::check_exception( + term_it != uuid_to_index_.end(), nullptr, + OpenGeodeException::TYPE::data, + absl::StrCat( "[EnergyTermCollection] Unknown energy term: ", + term_uuid.string() ) ); + return term_it->second; + } + + [[nodiscard]] const EnergyTerm< ObjectType >& get( + const uuid& term_uuid ) const + { + return *energy_terms_[get_term_index( term_uuid )]; + } + + [[nodiscard]] uuid get_term_uuid( std::string_view name ) const + { + auto uuid_it = name_to_uuid_.find( name ); + OpenGeodeStochasticStochasticException::check_exception( + uuid_it != name_to_uuid_.end(), nullptr, + OpenGeodeException::TYPE::data, + absl::StrCat( + "[EnergyTermCollection] Unknown energy term: ", name ) ); + return uuid_it->second; + } + + [[nodiscard]] const std::vector< + std::unique_ptr< EnergyTerm< ObjectType > > >& + energy_terms() const + { + return energy_terms_; + } + + [[nodiscard]] std::string string() const + { + auto message = absl::StrCat( + "EnergyTermCollection: ", energy_terms_.size(), " terms:" ); + for( const auto& term : energy_terms_ ) + { + absl::StrAppend( &message, "\n\t --> ", term->string() ); + } + return message; + } + + private: + absl::flat_hash_map< std::string, uuid > name_to_uuid_; + absl::flat_hash_map< uuid, index_t > uuid_to_index_; + std::vector< std::unique_ptr< EnergyTerm< ObjectType > > > + energy_terms_; + }; + +} // namespace geode \ No newline at end of file diff --git a/include/geode/stochastic/sampling/mcmc/energy_terms/energy_term.hpp b/include/geode/stochastic/models/energy_terms/energy_term.hpp similarity index 80% rename from include/geode/stochastic/sampling/mcmc/energy_terms/energy_term.hpp rename to include/geode/stochastic/models/energy_terms/energy_term.hpp index 4df8069..4db0ecc 100644 --- a/include/geode/stochastic/sampling/mcmc/energy_terms/energy_term.hpp +++ b/include/geode/stochastic/models/energy_terms/energy_term.hpp @@ -42,7 +42,7 @@ namespace geode::detail { OpenGeodeStochasticStochasticException::check_exception( param >= 0., nullptr, OpenGeodeException::TYPE::data, - "[Gibbs energy term] - The model parameter cannot be " + "[EnergyTerm] Model parameter cannot be " "negative." ); if( param >= geode::GLOBAL_EPSILON ) @@ -90,13 +90,17 @@ namespace geode public: explicit EnergyTerm( std::string_view name, double param, - std::vector< uuid >&& targeted_set_ids, + std::vector< uuid >&& impacted_set_ids, const SpatialDomain< ObjectType::dim >& domain ) : energy_scale_{ param }, - targeted_set_ids_{ std::move( targeted_set_ids ) }, + impacted_set_ids_{ std::move( impacted_set_ids ) }, domain_( domain ) { - std::sort( targeted_set_ids_.begin(), targeted_set_ids_.end() ); + absl::c_sort( impacted_set_ids_ ); + impacted_set_ids_.erase( std::unique( impacted_set_ids_.begin(), + impacted_set_ids_.end() ), + impacted_set_ids_.end() ); + impacted_set_ids_.shrink_to_fit(); IdentifierBuilder builder( *this ); builder.set_name( name ); } @@ -108,9 +112,9 @@ namespace geode return energy_scale_.parameter(); } - [[nodiscard]] const std::vector< uuid >& targeted_set_ids() const + [[nodiscard]] const std::vector< uuid >& impacted_set_ids() const { - return targeted_set_ids_; + return impacted_set_ids_; } /// Energy contribution for a given statistic multiplier @@ -144,9 +148,9 @@ namespace geode absl::StrCat( "Term : ", name().value_or( id().string() ), "; uuid: ", id().string(), " parameter value: ", energy_scale_.parameter(), - " applyied on ", targeted_set_ids_.size(), + " applyied on ", impacted_set_ids_.size(), " object subsets -->" ); - for( const auto& set_id : targeted_set_ids_ ) + for( const auto& set_id : impacted_set_ids_ ) { absl::StrAppend( &message, "\t", set_id.string() ); } @@ -154,10 +158,10 @@ namespace geode } protected: - [[nodiscard]] bool is_targeted_set( const uuid& set_id ) const + [[nodiscard]] bool is_impacted_set( const uuid& set_id ) const { return std::binary_search( - targeted_set_ids_.begin(), targeted_set_ids_.end(), set_id ); + impacted_set_ids_.begin(), impacted_set_ids_.end(), set_id ); } [[nodiscard]] const SpatialDomain< ObjectType::dim >& domain() const @@ -165,30 +169,24 @@ namespace geode return domain_; } - // bool intersects_domain( const ObjectType& obj ) const - // { - // return SpatialDomainChecker< ObjectType - // >::intersects_domain( - // domain_, obj ); - // } - template < typename Func > - void for_each_targeted_object( - const ObjectSets< ObjectType >& state, Func&& do_apply ) const + void for_each_object_in_sets( const ObjectSets< ObjectType >& state, + const std::vector< uuid >& set_ids, + Func&& do_apply ) const { - for( const auto& targeted_set_id : targeted_set_ids_ ) + for( const auto& set_id : set_ids ) { - for( const auto set_id : - state.get_objects_in_set( targeted_set_id ) ) + for( const auto object_ids : + state.get_objects_in_set( set_id ) ) { - std::forward< Func >( do_apply )( set_id ); + std::forward< Func >( do_apply )( object_ids ); } } } private: detail::EnergyScale energy_scale_; - std::vector< uuid > targeted_set_ids_; + std::vector< uuid > impacted_set_ids_; const SpatialDomain< ObjectType::dim >& domain_; }; } // namespace geode \ No newline at end of file diff --git a/include/geode/stochastic/models/energy_terms/energy_term_builder.hpp b/include/geode/stochastic/models/energy_terms/energy_term_builder.hpp new file mode 100644 index 0000000..119f01c --- /dev/null +++ b/include/geode/stochastic/models/energy_terms/energy_term_builder.hpp @@ -0,0 +1,148 @@ +/* + * Copyright (c) 2019 - 2026 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 + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include + +namespace geode +{ + template < typename ObjectType > + std::unique_ptr< EnergyTerm< ObjectType > > build_single_term( + const SingleObjectTermConfig& cfg, + const ObjectSets< ObjectType >& object_sets, + const SpatialDomain< ObjectType::dim >& domain ) + { + auto set_ids = object_sets.get_set_uuids( cfg.object_set_names ); + + auto object_feature = + build_single_object_feature< ObjectType >( cfg.object_feature ); + + return std::make_unique< SingleObjectTerm< ObjectType > >( + cfg.term_name, cfg.lambda, std::move( set_ids ), domain, + std::move( object_feature ) ); + } + + std::pair< std::vector< uuid >, + absl::flat_hash_map< uuid, std::vector< uuid > > > + opengeode_stochastic_stochastic_api + pairwise_builder_initialize_interactions_helper( + const std::vector< std::pair< uuid, uuid > >& interacting_sets ); + + template < typename ObjectType > + std::unique_ptr< EnergyTerm< ObjectType > > build_pairwise_term( + const PairwiseTermConfig& cfg, + const ObjectSets< ObjectType >& object_sets, + const SpatialDomain< ObjectType::dim >& domain ) + { + auto set_id_interactions = + object_sets.get_set_uuid_pairs( cfg.object_set_names_interactions ); + auto [interacting_set_ids, adjacent_set_uuids] = + pairwise_builder_initialize_interactions_helper( + set_id_interactions ); + auto interaction = + build_pairwise_interaction< ObjectType >( cfg.interaction_config ); + + return std::make_unique< geode::PairwiseTerm< ObjectType > >( + cfg.term_name, cfg.gamma, std::move( interacting_set_ids ), domain, + std::move( adjacent_set_uuids ), std::move( interaction ) ); + } + + template < typename ObjectType > + std::unique_ptr< EnergyTerm< ObjectType > > build_energy_term_impl( + const SingleObjectTermConfig& cfg, + const ObjectSets< ObjectType >& object_sets, + const SpatialDomain< ObjectType::dim >& domain ) + { + return build_single_term< ObjectType >( cfg, object_sets, domain ); + } + + template < typename ObjectType > + std::unique_ptr< EnergyTerm< ObjectType > > build_energy_term_impl( + const PairwiseTermConfig& cfg, + const ObjectSets< ObjectType >& object_sets, + const SpatialDomain< ObjectType::dim >& domain ) + { + return build_pairwise_term< ObjectType >( cfg, object_sets, domain ); + } + + template < typename ObjectType, typename NewEnergyTypeConfig > + std::unique_ptr< EnergyTerm< ObjectType > > build_energy_term_impl( + const NewEnergyTypeConfig& /*unused*/, + const ObjectSets< ObjectType >& /*unused*/, + const SpatialDomain< ObjectType::dim >& /*unused*/ ) + { + static_assert( sizeof( NewEnergyTypeConfig ) == 0, + "Unsupported EnergyTermConfig type" ); + return nullptr; + } + + template < typename ObjectType > + std::unique_ptr< EnergyTerm< ObjectType > > build_energy_term_impl( + const std::monostate& /*unused*/, + const ObjectSets< ObjectType >& /*unused*/, + const SpatialDomain< ObjectType::dim >& /*unused*/ ) + { + throw OpenGeodeStochasticStochasticException{ nullptr, + OpenGeodeException::TYPE::data, + "[EnergyTermBuilder] Energy term config not initialized" }; + } + + template < typename ObjectType > + std::unique_ptr< EnergyTerm< ObjectType > > build_energy_term( + const EnergyTermConfig& cfg, + const ObjectSets< ObjectType >& object_sets, + const SpatialDomain< ObjectType::dim >& domain ) + { + return std::visit( + [&object_sets, &domain]( const auto& term_cfg ) + -> std::unique_ptr< EnergyTerm< ObjectType > > { + return build_energy_term_impl< ObjectType >( + term_cfg, object_sets, domain ); + }, + cfg ); + } + +} // namespace geode \ No newline at end of file diff --git a/include/geode/stochastic/models/energy_terms/energy_term_config.hpp b/include/geode/stochastic/models/energy_terms/energy_term_config.hpp new file mode 100644 index 0000000..a897c5d --- /dev/null +++ b/include/geode/stochastic/models/energy_terms/energy_term_config.hpp @@ -0,0 +1,56 @@ +/* + * Copyright (c) 2019 - 2026 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 +{ + struct SingleObjectTermConfig + { + std::string term_name; + std::vector< std::string > object_set_names; + double lambda; + + SingleObjectFeatureConfig object_feature; + }; + + struct PairwiseTermConfig + { + std::string term_name; + std::vector< std::pair< std::string, std::string > > + object_set_names_interactions; + double gamma; + + PairwiseInteractionConfig interaction_config; + }; + + using EnergyTermConfig = std:: + variant< std::monostate, SingleObjectTermConfig, PairwiseTermConfig >; + +} // namespace geode \ No newline at end of file diff --git a/include/geode/stochastic/models/energy_terms/pairwise_term.hpp b/include/geode/stochastic/models/energy_terms/pairwise_term.hpp new file mode 100644 index 0000000..c4feabb --- /dev/null +++ b/include/geode/stochastic/models/energy_terms/pairwise_term.hpp @@ -0,0 +1,198 @@ +/* + * Copyright (c) 2019 - 2026 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 +{ + template < typename ObjectType > + class PairwiseTerm : public EnergyTerm< ObjectType > + { + public: + explicit PairwiseTerm( std::string_view name, + double gamma, + std::vector< uuid >&& impacted_set_ids, + const SpatialDomain< ObjectType::dim >& domain, + absl::flat_hash_map< uuid, std::vector< uuid > >&& + objectset_adjacency_map, + std::unique_ptr< PairwiseInteraction< ObjectType > > interaction ) + : EnergyTerm< ObjectType >( + name, gamma, std::move( impacted_set_ids ), domain ), + objectset_adjacency_map_( std::move( objectset_adjacency_map ) ), + interaction_( std::move( interaction ) ) + { + } + + [[nodiscard]] double total_log( + const ObjectSets< ObjectType >& state ) const final + { + const auto interaction_weight = statistic( state ); + return this->contribution( interaction_weight ); + } + + [[nodiscard]] double delta_log_add( + const ObjectSets< ObjectType >& state, + const ObjectRef< ObjectType >& new_object ) const final + { + if( !this->is_impacted_set( new_object.set_id ) ) + { + return 0.0; + } + auto delta = compute_local_interactions_with_neighbors( + new_object, std::nullopt, state ); + return this->contribution( delta ); + } + + [[nodiscard]] double delta_log_remove( + const ObjectSets< ObjectType >& state, + const ObjectId& object_id ) const override + { + if( !this->is_impacted_set( object_id.set_id ) ) + { + return 0.0; + } + ObjectRef< ObjectType > old_object{ state.get_object( object_id ), + object_id.set_id }; + auto delta = compute_local_interactions_with_neighbors( + old_object, object_id, state ); + return this->contribution( -delta ); + } + + [[nodiscard]] double delta_log_change( + const ObjectSets< ObjectType >& state, + const ObjectId& old_object_id, + const ObjectRef< ObjectType >& new_object ) const override + { + auto delta_new = compute_local_interactions_with_neighbors( + new_object, old_object_id, state ); + ObjectRef< ObjectType > old_object{ + state.get_object( old_object_id ), old_object_id.set_id + }; + auto delta_old = compute_local_interactions_with_neighbors( + old_object, old_object_id, state ); + auto delta = delta_new - delta_old; + return this->contribution( delta ); + } + + [[nodiscard]] double statistic( + const ObjectSets< ObjectType >& state ) const override + { + double sum = 0.0; + this->for_each_object_in_sets( state, this->impacted_set_ids(), + [&sum, &state, this]( const ObjectId& cur_obj_id ) { + sum += this->accumulate_interactions_with_neighbors( + cur_obj_id, state ); + } ); + return sum; + } + + private: + double compute_local_interactions_with_neighbors( + const ObjectRef< ObjectType >& object_ref, + std::optional< ObjectId > exclude_id, + const ObjectSets< ObjectType >& state ) const + { + 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, + impacted_set_it->second, + interaction_->neighborhood_searching_distance(), exclude_id ); + double sum = 0.0; + for( const auto& neigh_id : neighbors ) + { + ObjectRef< ObjectType > neigh_object{ + state.get_object( neigh_id ), neigh_id.set_id + }; + if( !is_any_in_domain( + object_ref.object, neigh_object.object ) ) + { + continue; + } + sum += interaction_->evaluate( object_ref, neigh_object ); + } + return sum; + } + + double accumulate_interactions_with_neighbors( + const ObjectId& object_id, + const ObjectSets< ObjectType >& state ) const + { + const auto& cur_obj = state.get_object( object_id ); + if( !is_in_domain( cur_obj ) ) + { + return 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, + 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 ); + if( neigh_id < object_id && is_in_domain( neigh_obj ) ) + { + continue; + } + ObjectRef< ObjectType > neigh_object{ neigh_obj, + neigh_id.set_id }; + sum += interaction_->evaluate( object_ref, neigh_object ); + } + return sum; + } + + bool is_any_in_domain( + const ObjectType& object1, const ObjectType& object2 ) const + { + return is_in_domain( object1 ) || is_in_domain( object2 ); + } + + bool is_in_domain( const ObjectType& object ) const + { + return SpatialDomainChecker< ObjectType >::is_anchored_in_domain( + this->domain(), object ); + } + + private: + absl::flat_hash_map< uuid, std::vector< uuid > > + objectset_adjacency_map_; + std::unique_ptr< PairwiseInteraction< ObjectType > > interaction_; + }; +} // namespace geode \ No newline at end of file diff --git a/include/geode/stochastic/sampling/mcmc/energy_terms/single_object_term.hpp b/include/geode/stochastic/models/energy_terms/single_object_term.hpp similarity index 66% rename from include/geode/stochastic/sampling/mcmc/energy_terms/single_object_term.hpp rename to include/geode/stochastic/models/energy_terms/single_object_term.hpp index d98240c..fab3286 100644 --- a/include/geode/stochastic/sampling/mcmc/energy_terms/single_object_term.hpp +++ b/include/geode/stochastic/models/energy_terms/single_object_term.hpp @@ -24,59 +24,55 @@ #include -#include +#include +#include #include namespace geode { - template < typename ObjectType, typename ObjectContributionFunc > + template < typename ObjectType > class SingleObjectTerm : public EnergyTerm< ObjectType > { public: explicit SingleObjectTerm( std::string_view name, double lambda, - std::vector< uuid > targeted_set_ids, - double scale, - ObjectContributionFunc contribution_func, - const SpatialDomain< ObjectType::dim >& domain ) + std::vector< uuid >&& impacted_set_ids, + const SpatialDomain< ObjectType::dim >& domain, + std::unique_ptr< SingleObjectFeature< ObjectType > > feature ) : EnergyTerm< ObjectType >( - name, lambda, std::move( targeted_set_ids ), domain ), - scale_( scale ), - contribution_func_( std::move( contribution_func ) ) + name, lambda, std::move( impacted_set_ids ), domain ), + feature_( std::move( feature ) ) { } [[nodiscard]] double total_log( const ObjectSets< ObjectType >& state ) const override { - return this->contribution( scale_ * statistic( state ) ); + return this->contribution( statistic( state ) ); } [[nodiscard]] double delta_log_add( const ObjectSets< ObjectType >& /*state*/, const ObjectRef< ObjectType >& new_object ) const override { - if( !this->is_targeted_set( new_object.set_id ) ) + if( !this->is_impacted_set( new_object.set_id ) ) { return 0.0; } return this->contribution( - scale_ - * contribution_func_( new_object.object, this->domain() ) ); + feature_->evaluate( new_object.object, this->domain() ) ); } [[nodiscard]] double delta_log_remove( const ObjectSets< ObjectType >& state, const ObjectId& object_id ) const override { - if( !this->is_targeted_set( object_id.set_id ) ) + if( !this->is_impacted_set( object_id.set_id ) ) { return 0.0; } - return this->contribution( - -scale_ - * contribution_func_( - state.get_object( object_id ), this->domain() ) ); + return this->contribution( -feature_->evaluate( + state.get_object( object_id ), this->domain() ) ); } [[nodiscard]] double delta_log_change( @@ -84,27 +80,30 @@ namespace geode const ObjectId& old_object_id, const ObjectRef< ObjectType >& new_object ) const override { + if( !this->is_impacted_set( new_object.set_id ) ) + { + return 0.0; + } double delta = - contribution_func_( new_object.object, this->domain() ) - - contribution_func_( + feature_->evaluate( new_object.object, this->domain() ) + - feature_->evaluate( state.get_object( old_object_id ), this->domain() ); - return this->contribution( scale_ * delta ); + return this->contribution( delta ); } [[nodiscard]] double statistic( const ObjectSets< ObjectType >& state ) const override { double sum = 0.0; - this->for_each_targeted_object( - state, [&]( 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 += contribution_func_( obj, this->domain() ); + sum += this->feature_->evaluate( obj, this->domain() ); } ); return sum; } private: - double scale_{ 1. }; - ObjectContributionFunc contribution_func_; + std::unique_ptr< SingleObjectFeature< ObjectType > > feature_; }; } // namespace geode \ No newline at end of file diff --git a/include/geode/stochastic/sampling/mcmc/energy_terms/gibbs_energy.hpp b/include/geode/stochastic/models/gibbs_energy.hpp similarity index 76% rename from include/geode/stochastic/sampling/mcmc/energy_terms/gibbs_energy.hpp rename to include/geode/stochastic/models/gibbs_energy.hpp index 2718107..e626878 100644 --- a/include/geode/stochastic/sampling/mcmc/energy_terms/gibbs_energy.hpp +++ b/include/geode/stochastic/models/gibbs_energy.hpp @@ -22,7 +22,7 @@ */ #pragma once -#include +#include #include namespace geode @@ -41,21 +41,8 @@ namespace geode const ObjectSets< ObjectType >& state ) const { double log_energy = 0.0; - const auto& energy_terms = energy_terms_collection_.all_terms(); - for( auto& [term_id, term] : energy_terms ) - { - geode_unused( term_id ); - log_energy += term->total_log( state ); - } - return log_energy; - } - - [[nodiscard]] double total_log_energy_for_set( - const ObjectSets< ObjectType >& state, const uuid& set_id ) const - { - double log_energy = 0.0; - for( const auto term : - energy_terms_collection_.terms_for_set( set_id ) ) + const auto& energy_terms = energy_terms_collection_.energy_terms(); + for( const auto& term : energy_terms ) { log_energy += term->total_log( state ); } @@ -67,8 +54,7 @@ namespace geode const ObjectRef< ObjectType >& new_object ) const { double log_energy = 0.0; - for( const auto& term : - energy_terms_collection_.terms_for_set( new_object.set_id ) ) + for( const auto& term : energy_terms_collection_.energy_terms() ) { log_energy += term->delta_log_add( state, new_object ); } @@ -80,8 +66,7 @@ namespace geode const ObjectId& object_id ) const { double log_energy = 0.0; - for( const auto& term : - energy_terms_collection_.terms_for_set( object_id.set_id ) ) + for( const auto& term : energy_terms_collection_.energy_terms() ) { log_energy += term->delta_log_remove( state, object_id ); } @@ -94,8 +79,7 @@ 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 ) ) + for( const auto& term : energy_terms_collection_.energy_terms() ) { log_energy += term->delta_log_change( state, old_id, new_object ); diff --git a/include/geode/stochastic/models/model.hpp b/include/geode/stochastic/models/model.hpp new file mode 100644 index 0000000..532cec8 --- /dev/null +++ b/include/geode/stochastic/models/model.hpp @@ -0,0 +1,150 @@ +/* + * Copyright (c) 2019 - 2026 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 + +#include +#include +#include + +#include + +#include +#include +namespace geode +{ + + struct ModelConfig + { + std::vector< EnergyTermConfig > terms; + }; + + template < typename ObjectType > + class Model + { + OPENGEODE_DISABLE_COPY_AND_MOVE( Model ); + + public: + Model() = delete; + ~Model() = default; + + explicit Model( EnergyTermCollection< ObjectType >&& energy_terms ) + : terms_collection_( std::move( energy_terms ) ), + energy_{ terms_collection_ } + { + } + + [[nodiscard]] index_t nb_terms() const + { + return terms_collection_.size(); + } + + [[nodiscard]] const EnergyTermCollection< ObjectType >& terms() const + { + return terms_collection_; + } + + [[nodiscard]] index_t term_index( const uuid& term_uuid ) const + { + return terms_collection_.get_term_index( term_uuid ); + } + + [[nodiscard]] const GibbsEnergy< ObjectType >& energy() const + { + return energy_; + } + + [[nodiscard]] std::vector< double > compute_statistics( + const ObjectSets< ObjectType >& state ) const + { + std::vector< double > statistic_values; + statistic_values.reserve( terms_collection_.size() ); + for( const auto& term : terms_collection_.energy_terms() ) + { + statistic_values.emplace_back( term->statistic( state ) ); + } + return statistic_values; + } + + [[nodiscard]] double compute_statistic( + const ObjectSets< ObjectType >& state, const uuid& term_uuid ) const + { + const auto& term = terms_collection_.get( term_uuid ); + return term.statistic( state ); + } + + [[nodiscard]] std::string term_name( const uuid& term_uuid ) const + { + const auto& term = terms_collection_.get( term_uuid ); + return term.name().value_or( "unnamed" ); + } + + [[nodiscard]] std::vector< std::string > term_names() const + { + std::vector< std::string > names; + names.reserve( terms_collection_.size() ); + + for( const auto& term : terms_collection_.energy_terms() ) + { + names.emplace_back( term->name().value_or( "unnamed" ) ); + } + + return names; + } + + [[nodiscard]] std::string string() const + { + return terms_collection_.string(); + } + + private: + EnergyTermCollection< ObjectType > terms_collection_; + GibbsEnergy< ObjectType > energy_; + }; + + template < typename ObjectType > + std::unique_ptr< Model< ObjectType > > build_model( + const ModelConfig& config, + const ObjectSets< ObjectType >& object_sets, + const SpatialDomain< ObjectType::dim >& domain ) + { + EnergyTermCollection< ObjectType > collection; + + for( const auto& term_cfg : config.terms ) + { + collection.add_energy_term( build_energy_term< ObjectType >( + term_cfg, object_sets, domain ) ); + } + + return std::make_unique< Model< ObjectType > >( + std::move( collection ) ); + } +} // namespace geode \ No newline at end of file diff --git a/include/geode/stochastic/sampling/direct/double_sampler.hpp b/include/geode/stochastic/sampling/direct/double_sampler.hpp index c02e4b6..2e94317 100644 --- a/include/geode/stochastic/sampling/direct/double_sampler.hpp +++ b/include/geode/stochastic/sampling/direct/double_sampler.hpp @@ -48,6 +48,13 @@ namespace geode struct DistributionDescription { + explicit DistributionDescription( absl::string_view name_in ) + : name{ name_in } + { + } + + DistributionDescription() = default; + std::string name{ "default_distribution" }; DistributionType distribution_type; @@ -60,7 +67,7 @@ namespace geode std::optional< double > alpha; // exponent parameter for power law distribution - std::string string() const + [[nodiscard]] std::string string() const { auto message = absl::StrCat( "Distribution - ", name ); absl::StrAppend( &message, diff --git a/include/geode/stochastic/sampling/direct/object_set_sampler/object_set_sampler.hpp b/include/geode/stochastic/sampling/direct/object_set_sampler/object_set_sampler.hpp index e139fab..c978a66 100644 --- a/include/geode/stochastic/sampling/direct/object_set_sampler/object_set_sampler.hpp +++ b/include/geode/stochastic/sampling/direct/object_set_sampler/object_set_sampler.hpp @@ -24,14 +24,20 @@ #pragma once #include + #include +#include namespace geode { + template < typename Type > class ObjectSetSampler { + OPENGEODE_DISABLE_COPY_AND_MOVE( ObjectSetSampler ); + public: + ObjectSetSampler() = default; virtual ~ObjectSetSampler() = default; [[nodiscard]] virtual Type sample( RandomEngine& engine ) const = 0; @@ -50,9 +56,21 @@ namespace geode } protected: - virtual bool is_valid_object( const Type& obj ) const = 0; + [[nodiscard]] virtual bool is_valid_object( const Type& obj ) const = 0; + void set_log_pdf( double value ) + { + log_pdf_ = value; + } - protected: + private: double log_pdf_{ LOG_PROB_INVALID }; }; + + template < typename ObjectType > + struct ObjectSamplerConfig; + + template < typename ObjectType > + std::unique_ptr< ObjectSetSampler< ObjectType > > build_objectset_sampler( + const SpatialDomain< ObjectType::dim >& domain, + const ObjectSamplerConfig< ObjectType >& config ); } // namespace geode \ No newline at end of file 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 3f2145f..db7749c 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 @@ -23,75 +23,35 @@ #pragma once -#include -#include - #include -#include -#include namespace geode { + template < index_t dimension > + struct ObjectSamplerConfig< Point< dimension > > + { + // use to define the step for change move (move_ratio*domain volume) + // NOLINTNEXTLINE(*-magic-numbers) + double move_ratio{ 0.1 }; + }; + template < index_t dimension > class UniformPointSetSampler : public ObjectSetSampler< Point< dimension > > { public: - UniformPointSetSampler( const SpatialDomain< dimension >& domain ) - : ObjectSetSampler< Point< dimension > >{}, domain_( domain ) - { - auto volume = domain_.extended_n_volume(); - OpenGeodeStochasticStochasticException::check_exception( - volume != 0., nullptr, OpenGeodeException::TYPE::data, - "[PointSetSampler] - Undefined Extended Bounding " - "Box (volume ==0)." ); - this->log_pdf_ = -std::log( volume ); - step_move_ = define_step_for_move(); - OpenGeodeStochasticStochasticException::check_exception( - step_move_ > 0., nullptr, OpenGeodeException::TYPE::data, - "[PointSetSampler] - Undefined step length for move (value == ", - step_move_, ")." ); - } - - Point< dimension > sample( RandomEngine& engine ) const override - { - return PointUniformSampler::sample< dimension >( - engine, domain_.extended_box() ); - } + UniformPointSetSampler( const SpatialDomain< dimension >& domain, + const ObjectSamplerConfig< Point< dimension > >& config ); - Point< dimension > change( - const Point< dimension >& obj, RandomEngine& engine ) const override - { - geode::Sphere< dimension > ball{ obj, step_move_ }; + [[nodiscard]] Point< dimension > sample( + RandomEngine& engine ) const override; - auto new_point = - PointUniformSampler::sample< dimension >( engine, ball ); - constexpr index_t max_try{ 100 }; - for( const auto try_id : geode::Range{ max_try } ) - { - if( domain_.extended_contains( new_point ) ) - { - return new_point; - } - new_point = - PointUniformSampler::sample< dimension >( engine, ball ); - } - throw OpenGeodeStochasticStochasticException{ nullptr, - OpenGeodeException::TYPE::internal, - "[PointSampler] - Cannot find a point in the " - "extended domain" }; - } + [[nodiscard]] Point< dimension > change( const Point< dimension >& obj, + RandomEngine& engine ) const override; private: - double define_step_for_move() - { - double ratio = 0.1; - return ratio * domain_.smallest_length(); - } - - bool is_valid_object( const Point< dimension >& obj ) const override - { - return domain_.extended_contains( obj ); - } + [[nodiscard]] double define_step_for_move( double ratio ); + [[nodiscard]] bool is_valid_object( + const Point< dimension >& obj ) const override; private: const SpatialDomain< dimension >& domain_; 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 index 0e3e07d..8d44851 100644 --- 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 @@ -23,85 +23,43 @@ #pragma once -#include -#include - +#include #include -#include -#include -#include namespace geode { + template <> + struct ObjectSamplerConfig< OwnerSegment2D > + { + // NOLINTNEXTLINE(*-magic-numbers) + double move_ratio{ 0.1 }; + + DoubleSampler::DistributionDescription length{ "length" }; + + DoubleSampler::DistributionDescription azimuth{ "azimuth" }; + }; + class UniformSegmentSetSampler : public ObjectSetSampler< OwnerSegment2D > { public: UniformSegmentSetSampler( const SpatialDomain< 2 >& domain, - const DoubleSampler::Distribution& length, - const DoubleSampler::Distribution& azimuth ) - : ObjectSetSampler< OwnerSegment2D >{}, - domain_{ domain }, - length_{ length }, - azimuth_{ azimuth } - { - auto volume = domain_.extended_n_volume(); - OpenGeodeStochasticStochasticException::check_exception( - volume != 0., nullptr, OpenGeodeException::TYPE::data, - "[SegmentSetSampler] - Undefined Extended Bounding " - "Box (volume ==0)." ); - this->log_pdf_ = -std::log( volume ); - } - - OwnerSegment2D sample( RandomEngine& engine ) const override - { - auto seg = SegmentUniformSampler::sample( - engine, domain_.extended_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 ) ); - const auto other = 1 - current; + const ObjectSamplerConfig< OwnerSegment2D >& config ); - geode::Sphere< 2 > ball{ extremities[current], - ratio * obj.length() }; + [[nodiscard]] OwnerSegment2D sample( + RandomEngine& engine ) const override; - auto new_point = PointUniformSampler::sample< 2 >( engine, ball ); - constexpr index_t max_try{ 100 }; - for( const auto try_id : geode::Range{ max_try } ) - { - if( domain_.extended_contains( new_point ) - || domain_.extended_contains( extremities[other] ) ) - { - OwnerSegment2D new_segment{ obj }; - new_segment.set_point( current, new_point ); - return new_segment; - } - new_point = PointUniformSampler::sample< 2 >( engine, ball ); - } - throw OpenGeodeStochasticStochasticException{ nullptr, - OpenGeodeException::TYPE::internal, - "[SegmentSetSampler] - Cannot find a point in the box" }; - return obj; - } + [[nodiscard]] OwnerSegment2D change( + const OwnerSegment2D& obj, RandomEngine& engine ) const override; private: - bool is_valid_object( const OwnerSegment2D& obj ) const override - { - const auto& extremities = obj.vertices(); - return domain_.extended_contains( extremities[0] ) - || domain_.extended_contains( extremities[1] ); - } + [[nodiscard]] bool is_valid_object( + const OwnerSegment2D& obj ) const override; private: const SpatialDomain< 2 >& domain_; DoubleSampler::Distribution length_; DoubleSampler::Distribution azimuth_; + double move_ratio_; }; } // namespace geode \ No newline at end of file diff --git a/include/geode/stochastic/sampling/mcmc/energy_terms/density_term.hpp b/include/geode/stochastic/sampling/mcmc/energy_terms/density_term.hpp deleted file mode 100644 index 91e3858..0000000 --- a/include/geode/stochastic/sampling/mcmc/energy_terms/density_term.hpp +++ /dev/null @@ -1,61 +0,0 @@ -/* - * Copyright (c) 2019 - 2026 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 -{ - template < typename ObjectType > - class DensityTerm : public SingleObjectTerm< ObjectType, - std::function< double( const ObjectType&, - const SpatialDomain< ObjectType::dim >& ) > > - { - public: - explicit DensityTerm( std::string_view name, - double lambda, - std::vector< uuid > targeted_set_ids, - const SpatialDomain< ObjectType::dim >& domain ) - : SingleObjectTerm< ObjectType, - std::function< double( const ObjectType&, - const SpatialDomain< ObjectType::dim >& ) > >( - name, - lambda, - std::move( targeted_set_ids ), - 1.0, // scale by domain area to get density per unit - []( const ObjectType& obj, - const SpatialDomain< ObjectType::dim >& spatial_domain ) { - if( SpatialDomainChecker< ObjectType >:: - is_anchored_in_domain( spatial_domain, obj ) ) - { - return 1.0; - } - return 0.0; - }, // contribution = 1 anchoredin domain - domain ) - { - } - }; -} // namespace geode \ No newline at end of file diff --git a/include/geode/stochastic/sampling/mcmc/energy_terms/energy_term_collection.hpp b/include/geode/stochastic/sampling/mcmc/energy_terms/energy_term_collection.hpp deleted file mode 100644 index bbb76ea..0000000 --- a/include/geode/stochastic/sampling/mcmc/energy_terms/energy_term_collection.hpp +++ /dev/null @@ -1,128 +0,0 @@ -#pragma once - -#include -#include - -namespace geode -{ - template < typename ObjectType > - class EnergyTermCollection - { - OPENGEODE_DISABLE_COPY( EnergyTermCollection ); - - public: - EnergyTermCollection() = default; - EnergyTermCollection( EnergyTermCollection&& ) noexcept = default; - - ~EnergyTermCollection() = default; - - EnergyTermCollection& operator=( EnergyTermCollection&& ) = default; - - [[nodiscard]] uuid add_energy_term( - std::shared_ptr< EnergyTerm< ObjectType > > term ) - { - const uuid term_id = term->id(); - energy_terms_.emplace( term_id, term ); - for( const uuid& set_id : term->targeted_set_ids() ) - { - set_to_terms_[set_id].push_back( term ); - } - return term_id; - } - - [[nodiscard]] bool remove_energy_term( const uuid& term_id ) - { - auto term_it = energy_terms_.find( term_id ); - if( term_it == energy_terms_.end() ) - { - return false; - } - - auto term = term_it->second; - - for( const uuid& set_id : term->targeted_set_ids() ) - { - auto vec_it = set_to_terms_.find( set_id ); - if( vec_it == set_to_terms_.end() ) - { - continue; - } - - auto& vec = vec_it->second; - vec.erase( - std::remove( vec.begin(), vec.end(), term ), vec.end() ); - - if( vec.empty() ) - { - set_to_terms_.erase( vec_it ); - } - } - - energy_terms_.erase( term_it ); - return true; - } - - void clear() - { - energy_terms_.clear(); - set_to_terms_.clear(); - } - - [[nodiscard]] index_t size() const - { - return energy_terms_.size(); - } - - [[nodiscard]] const EnergyTerm< ObjectType >& get( - const uuid& term_id ) const - { - auto term_it = energy_terms_.find( term_id ); - OpenGeodeStochasticStochasticException::check_exception( - term_it != energy_terms_.end(), nullptr, - OpenGeodeException::TYPE::data, - "[EnergyTermCollection] Unknown energy term: ", - term_id.string() ); - return *term_it->second; - } - - [[nodiscard]] const absl::flat_hash_map< uuid, - std::shared_ptr< EnergyTerm< ObjectType > > >& - all_terms() const - { - return energy_terms_; - } - - [[nodiscard]] const std::vector< - std::shared_ptr< EnergyTerm< ObjectType > > >& - terms_for_set( const uuid& set_id ) const - { - const auto it = set_to_terms_.find( set_id ); - OpenGeodeStochasticStochasticException::check_exception( - it != set_to_terms_.end(), nullptr, - OpenGeodeException::TYPE::data, - "[EnergyTermCollection] - Object Subset (", set_id.string(), - ") does not have any energy term." ); - return it->second; - } - - [[nodiscard]] std::string string() const - { - auto message = absl::StrCat( - "EnergyTermCollection: ", energy_terms_.size(), " terms:" ); - for( const auto& [id, term] : energy_terms_ ) - { - absl::StrAppend( &message, "\n\t --> ", term->string() ); - } - return message; - } - - private: - absl::flat_hash_map< uuid, std::shared_ptr< EnergyTerm< ObjectType > > > - energy_terms_; - - absl::flat_hash_map< uuid, - std::vector< std::shared_ptr< EnergyTerm< ObjectType > > > > - set_to_terms_; - }; - -} // namespace geode \ No newline at end of file diff --git a/include/geode/stochastic/sampling/mcmc/energy_terms/pairwise_term.hpp b/include/geode/stochastic/sampling/mcmc/energy_terms/pairwise_term.hpp deleted file mode 100644 index 4f12ad1..0000000 --- a/include/geode/stochastic/sampling/mcmc/energy_terms/pairwise_term.hpp +++ /dev/null @@ -1,237 +0,0 @@ -/* - * Copyright (c) 2019 - 2026 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 - -namespace geode -{ - // rename PairwiseInteractionTerm - template < typename ObjectType > - class PairwiseTerm : public EnergyTerm< ObjectType > - { - public: - explicit PairwiseTerm( std::string_view name, - double gamma, - std::vector< uuid > targeted_set_ids, - std::unique_ptr< PairwiseInteraction< ObjectType > > interaction, - const SpatialDomain< ObjectType::dim >& domain ) - : EnergyTerm< ObjectType >( - name, gamma, std::move( targeted_set_ids ), domain ), - interaction_( std::move( interaction ) ) - { - } - - [[nodiscard]] double total_log( - const ObjectSets< ObjectType >& state ) const final - { - const auto interaction_weight = statistic( state ); - return this->contribution( interaction_weight ); - } - - [[nodiscard]] double delta_log_add( - const ObjectSets< ObjectType >& state, - const ObjectRef< ObjectType >& new_object ) const final - { - if( !this->is_targeted_set( new_object.set_id ) ) - { - return 0.0; - } - double delta = 0.0; - const auto neighbors = - state.neighbors( new_object.object, this->targeted_set_ids(), - interaction_->neighborhood_searching_distance() ); - for( const auto& neigh_id : neighbors ) - { - geode::ObjectRef< ObjectType > neigh_object{ - state.get_object( neigh_id ), neigh_id.set_id - }; - // is that really the good test? - // intersect? - if( SpatialDomainChecker< ObjectType >::is_anchored_in_domain( - this->domain(), new_object.object ) - || SpatialDomainChecker< - ObjectType >::is_anchored_in_domain( this->domain(), - neigh_object.object ) ) - { - delta += interaction_->evaluate( new_object, neigh_object ); - } - } - return this->contribution( delta ); - } - - [[nodiscard]] double delta_log_remove( - const ObjectSets< ObjectType >& state, - const ObjectId& object_id ) const override - { - if( !this->is_targeted_set( object_id.set_id ) ) - { - return 0.0; - } - double delta = 0.0; - ObjectRef< ObjectType > object_to_remove{ - state.get_object( object_id ), object_id.set_id - }; - const auto neighbors = - state.neighbors( object_id, this->targeted_set_ids(), - interaction_->neighborhood_searching_distance() ); - for( auto neigh_id : neighbors ) - { - ObjectRef< ObjectType > neigh_object{ - state.get_object( neigh_id ), neigh_id.set_id - }; - if( SpatialDomainChecker< ObjectType >::is_anchored_in_domain( - this->domain(), object_to_remove.object ) - || SpatialDomainChecker< - ObjectType >::is_anchored_in_domain( this->domain(), - neigh_object.object ) ) - { - delta += interaction_->evaluate( - object_to_remove, neigh_object ); - } - } - return this->contribution( -delta ); - } - - [[nodiscard]] double delta_log_change( - const ObjectSets< ObjectType >& state, - const ObjectId& old_object_id, - const ObjectRef< ObjectType >& new_object ) const override - { - if( !this->is_targeted_set( old_object_id.set_id ) - || !this->is_targeted_set( new_object.set_id ) ) - { - return 0.0; - } - double delta = 0.0; - - // Remove old object's interactions - ObjectRef< ObjectType > object_to_remove{ - state.get_object( old_object_id ), old_object_id.set_id - }; - const auto old_neighbors = - state.neighbors( old_object_id, this->targeted_set_ids(), - interaction_->neighborhood_searching_distance() ); - for( auto neigh_id : old_neighbors ) - { - ObjectRef< ObjectType > neigh_object{ - state.get_object( neigh_id ), neigh_id.set_id - }; - if( SpatialDomainChecker< ObjectType >::is_anchored_in_domain( - this->domain(), object_to_remove.object ) - || SpatialDomainChecker< - ObjectType >::is_anchored_in_domain( this->domain(), - neigh_object.object ) ) - { - delta -= interaction_->evaluate( - object_to_remove, neigh_object ); - } - } - - // Add new object's interactions - const auto new_neighbors = - state.neighbors( new_object.object, this->targeted_set_ids(), - interaction_->neighborhood_searching_distance() ); - for( auto neigh_id : new_neighbors ) - { - if( old_object_id == neigh_id ) - { - continue; // avoid double-counting - } - ObjectRef< ObjectType > neigh_object{ - state.get_object( neigh_id ), neigh_id.set_id - }; - if( SpatialDomainChecker< ObjectType >::is_anchored_in_domain( - this->domain(), new_object.object ) - || SpatialDomainChecker< - ObjectType >::is_anchored_in_domain( this->domain(), - neigh_object.object ) ) - { - delta += interaction_->evaluate( new_object, neigh_object ); - } - } - return this->contribution( delta ); - } - - [[nodiscard]] double statistic( - const ObjectSets< ObjectType >& state ) const override - { - double sum = 0.0; - this->for_each_targeted_object( state, [&]( const ObjectId& - obj_id ) { - const auto& cur_obj = state.get_object( obj_id ); - if( !SpatialDomainChecker< ObjectType >::is_anchored_in_domain( - this->domain(), cur_obj ) ) - { - return; - } - ObjectRef< ObjectType > object{ cur_obj, obj_id.set_id }; - const auto neighbors = - state.neighbors( obj_id, this->targeted_set_ids(), - interaction_->neighborhood_searching_distance() ); - for( const auto& neigh_obj_id : neighbors ) - { - if( !is_new_pair( cur_obj, obj_id, neigh_obj_id ) ) - { - continue; - } - ObjectRef< ObjectType > neigh_object{ - state.get_object( neigh_obj_id ), neigh_obj_id.set_id - }; - - sum += interaction_->evaluate( object, neigh_object ); - } - } ); - return sum; - } - - private: - [[nodiscard]] bool is_new_pair( const ObjectType& obj, - const ObjectId& obj_id, - const ObjectId& neigh_obj_id ) const - { - if( !SpatialDomainChecker< ObjectType >::is_anchored_in_domain( - this->domain(), obj ) ) - { - return true; - } - if( neigh_obj_id.set_id < obj_id.set_id ) - { - return false; - } - if( neigh_obj_id.set_id == obj_id.set_id - && neigh_obj_id.index <= obj_id.index ) - { - return false; - } - return true; - } - - private: - std::unique_ptr< PairwiseInteraction< ObjectType > > interaction_; - }; -} // namespace geode \ No newline at end of file diff --git a/include/geode/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp b/include/geode/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp index f4042e2..160aba0 100644 --- a/include/geode/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp +++ b/include/geode/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp @@ -25,14 +25,14 @@ #include #include +#include +#include +#include #include #include -#include -#include -#include #include #include -#include +#include #include #include diff --git a/include/geode/stochastic/sampling/mcmc/helpers/simulation_context.hpp b/include/geode/stochastic/sampling/mcmc/helpers/simulation_context.hpp new file mode 100644 index 0000000..1108bbe --- /dev/null +++ b/include/geode/stochastic/sampling/mcmc/helpers/simulation_context.hpp @@ -0,0 +1,148 @@ +/* + * Copyright (c) 2019 - 2026 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 +#include + +#include + +#include +#include + +namespace geode +{ + template < typename ObjectType > + struct SimulationContext + { + [[nodiscard]] std::string string() const + { + auto message = std::string{ "SimulationContext: " }; + absl::StrAppend( &message, "\n\t --> ", domain->string() ); + absl::StrAppend( &message, "\n\t --> ", object_sets->string() ); + absl::StrAppend( + &message, "\n\t --> ", set_samplers.size(), " Sets samplers " ); + absl::StrAppend( &message, "\n\t --> ", model->string() ); + // absl::StrAppend( &message, "\n\t --> ", mh_sampler_ > string() ); + + return message; + } + + std::unique_ptr< SpatialDomain< ObjectType::dim > > domain; + + std::unique_ptr< ObjectSets< ObjectType > > object_sets{ + std::make_unique< ObjectSets< ObjectType > >() + }; + std::vector< std::unique_ptr< geode::ObjectSetSampler< ObjectType > > > + set_samplers; + std::unique_ptr< Model< ObjectType > > model; + std::unique_ptr< geode::MetropolisHastings< ObjectType > > mh_sampler; + }; + + template < typename ObjectType > + struct ObjectSetDefinition + { + std::string name; + std::vector< ObjectType > fixed_objects; + + ObjectSamplerConfig< ObjectType > sampler; + ObjectSetDynamicsConfig dynamics; + }; + + template < typename ObjectType > + struct SimulationContextConfig + { + ObjectSetDefinition< ObjectType >& add_set( absl::string_view name ) + { + auto& set = sets.emplace_back(); + set.name = name; + return set; + } + + SpatialDomainConfig< ObjectType::dim > domain; + + std::vector< ObjectSetDefinition< ObjectType > > sets; + + geode::ModelConfig model; + }; + + template < typename ObjectType > + [[nodiscard]] geode::SimulationContext< ObjectType > + build_simulation_context( + const SimulationContextConfig< ObjectType >& config ) + { + geode::SimulationContext< ObjectType > context; + + // ------------------------- + // Domain + // ------------------------- + context.domain = geode::build_spatial_domain( config.domain ); + + // ------------------------- + // Sets + // ------------------------- + auto proposal_kernel = + std::make_unique< geode::ProposalKernel< ObjectType > >(); + for( const auto& set_def : config.sets ) + { + auto set_id = context.object_sets->add_set( set_def.name ); + for( auto fixed_object : set_def.fixed_objects ) + { + context.object_sets->add_object( + std::move( fixed_object ), set_id, true ); + } + auto sampler = + build_objectset_sampler( *context.domain, set_def.sampler ); + geode::add_birth_death_change_moves( *sampler, *proposal_kernel, + set_id, set_def.dynamics.birth_ratio, + set_def.dynamics.death_ratio, set_def.dynamics.change_ratio ); + context.set_samplers.emplace_back( std::move( sampler ) ); + } + + // ------------------------- + // Model + // ------------------------- + context.model = geode::build_model< ObjectType >( + config.model, *context.object_sets, *context.domain ); + + // ------------------------- + // MH sampler + // ------------------------- + context.mh_sampler = + std::make_unique< geode::MetropolisHastings< ObjectType > >( + *context.model, std::move( proposal_kernel ) ); + + return context; + } + +} // namespace geode \ No newline at end of file diff --git a/include/geode/stochastic/sampling/mcmc/helpers/simulation_monitor.hpp b/include/geode/stochastic/sampling/mcmc/helpers/simulation_monitor.hpp deleted file mode 100644 index 9f08b69..0000000 --- a/include/geode/stochastic/sampling/mcmc/helpers/simulation_monitor.hpp +++ /dev/null @@ -1,86 +0,0 @@ -/* - * Copyright (c) 2019 - 2026 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 ) - { - means_.resize( nb_energy_terms, 0.0 ); - variances_.resize( nb_energy_terms, 0.0 ); - } - - void add_realization( const std::vector< double >& values ) - { - OpenGeodeStochasticStochasticException::check_exception( - values.size() == means_.size(), nullptr, - OpenGeodeException::TYPE::data, - "[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 > 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 index 3ab9aa6..4811835 100644 --- a/include/geode/stochastic/sampling/mcmc/helpers/simulation_printer.hpp +++ b/include/geode/stochastic/sampling/mcmc/helpers/simulation_printer.hpp @@ -25,7 +25,7 @@ #include -#include +#include #include #include @@ -75,42 +75,84 @@ namespace geode } }; + template < typename ObjectType > class SimulationPrinter { public: - SimulationPrinter( const SimulationPrinterConfigurator& config ) - : config_( config ) + SimulationPrinter( const Model< ObjectType >& model, + const SimulationPrinterConfigurator& config ) + : model_( model ), config_( config ) { } // Print statistics to the configured statistics file - void print_statistics( - const std::vector< double >& stats, absl::string_view header ) const + void print_statistics( const std::vector< double >& stats ) const { if( !config_.print_statistics ) + { return; - const auto stats_file_path = - stats_file_path_.value_or( create_statistics_file( header ) ); + } + if( !stats_file_path_ ) + { + stats_file_path_ = + ( std::filesystem::path( config_.output_folder ) + / config_.statistics_filename ) + .string(); + create_file_and_write_header( *stats_file_path_, + absl::StrCat( "# Simulation Statistics\n", + energy_terms_name_header() ) ); + } std::ofstream file = - open_file_with_dirs( stats_file_path, std::ios::app ); + open_file_with_dirs( *stats_file_path_, std::ios::app ); file << absl::StrJoin( stats, " ; " ) << "\n"; } - template < typename ObjectType > + void print_statistics_summary( + const StatisticsTracker< ObjectType >& tracker ) const + { + if( !config_.print_statistics_summary ) + { + return; + } + + const auto summary_path = + ( std::filesystem::path( config_.output_folder ) + / config_.statistics_summary_filename ) + .string(); + create_file_and_write_header( + summary_path, absl::StrCat( "# Summary statistics\n", + energy_terms_name_header() ) ); + std::ofstream file = + open_file_with_dirs( summary_path, std::ios::app ); + file << absl::StrCat( + "# Count:\n", tracker.statiscal_count(), "\n" ); + file << absl::StrCat( + "# Means:\n", absl::StrJoin( tracker.means(), " ; " ), "\n" ); + file << absl::StrCat( "# Variances:\n", + absl::StrJoin( tracker.variances(), " ; " ), "\n" ); + } + 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(); + // const auto filename = + // ( std::filesystem::path( config_.output_folder ) + // / absl::StrCat( + // config_.realisations_prefix, realization_id, ".txt" ) + // ) + // .string(); + const auto filename = std::filesystem::path( config_.output_folder ) + / absl::StrCat( config_.realisations_prefix, + realization_id, ".txt" ); - std::ofstream file = open_file_with_dirs( filename ); + std::ofstream file = + open_file_with_dirs( filename, std::ofstream::out ); const auto all_objects = object_sets.get_all_object(); file << "#nb_objects\t" << all_objects.size() << "\n"; @@ -122,82 +164,51 @@ namespace geode << "\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 + void create_file_and_write_header( + const std::filesystem::path& 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 file = + open_file_with_dirs( filename, std::ofstream::out ); + + file << header; } - std::ofstream open_file_with_dirs( absl::string_view path_filename, - std::ios::openmode mode = std::ofstream::out ) const + std::ofstream open_file_with_dirs( + const std::filesystem::path& file_path, + std::ios::openmode mode ) const { - namespace fs = std::filesystem; - fs::path file_path{ std::string( path_filename ) }; - - if( !file_path.has_parent_path() ) + auto absolute_path = file_path; + if( !absolute_path.has_parent_path() ) { - file_path = fs::current_path() / file_path; + absolute_path = std::filesystem::current_path() / absolute_path; } - if( file_path.has_parent_path() ) + if( absolute_path.has_parent_path() ) { - fs::create_directories( file_path.parent_path() ); + std::filesystem::create_directories( + absolute_path.parent_path() ); } - std::ofstream file( file_path, mode ); + std::ofstream file( absolute_path, mode ); if( !file.is_open() ) { throw geode::OpenGeodeStochasticStochasticException{ nullptr, OpenGeodeException::TYPE::data, - "Cannot open file: " + file_path.string() }; + "Cannot open file: " + absolute_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_; + std::string energy_terms_name_header() const + { + return absl::StrCat( + absl::StrJoin( model_.term_names(), " ; " ), "\n" ); } private: - SimulationPrinterConfigurator config_; + const Model< ObjectType >& model_; + const SimulationPrinterConfigurator config_; mutable std::optional< std::string > stats_file_path_; }; diff --git a/include/geode/stochastic/sampling/mcmc/metropolis_hasting_sampler.hpp b/include/geode/stochastic/sampling/mcmc/metropolis_hasting_sampler.hpp index 785bdd0..2f96b79 100644 --- a/include/geode/stochastic/sampling/mcmc/metropolis_hasting_sampler.hpp +++ b/include/geode/stochastic/sampling/mcmc/metropolis_hasting_sampler.hpp @@ -24,22 +24,22 @@ #pragma once #include -#include +#include #include namespace geode { - enum struct MHDecision + enum struct MH_DECISION : local_index_t { - Accepted, - Rejected, - Undecided + accepted, + rejected, + undecided }; template < typename ObjectType > struct StepResult { - MHDecision decision{ MHDecision::Undecided }; + MH_DECISION decision{ MH_DECISION::undecided }; MoveType move_type{ MoveType::Invalid }; double log_accept{ -std::numeric_limits< double >::infinity() }; double delta_log_energy{ 0.0 }; @@ -49,15 +49,14 @@ namespace geode class MetropolisHastings { public: - MetropolisHastings( - const EnergyTermCollection< ObjectType >& energy_term_collection, + MetropolisHastings( const Model< ObjectType >& model, std::unique_ptr< ProposalKernel< ObjectType > > proposal_kernel ) - : gibbs_energy_{ energy_term_collection }, - proposal_kernel_( std::move( proposal_kernel ) ) + : model_{ model }, proposal_kernel_( std::move( proposal_kernel ) ) { OpenGeodeStochasticStochasticException::check_exception( proposal_kernel_ != nullptr, nullptr, - OpenGeodeException::TYPE::data, "[MH] null proposal kernel" ); + OpenGeodeException::TYPE::data, + "[MetropolisHastings] Proposal kernel is not defined." ); } StepResult< ObjectType > step( @@ -92,7 +91,8 @@ namespace geode } } - ObjectSets< ObjectType > walk_copy( ObjectSets< ObjectType > initial, + [[nodiscard]] ObjectSets< ObjectType > walk_copy( + ObjectSets< ObjectType > initial, RandomEngine& engine, index_t nb_steps ) const { @@ -100,62 +100,68 @@ namespace geode return initial; } - double beta() const + [[nodiscard]] double beta() const { return beta_; } - void set_beta( double b ) + void set_beta( double beta ) { - OpenGeodeStochasticStochasticException::check_exception( b >= 0.0, - nullptr, OpenGeodeException::TYPE::data, - "[MH] beta must be >= 0" ); - if( b == 0 ) + OpenGeodeStochasticStochasticException::check_exception( + beta >= 0.0, nullptr, OpenGeodeException::TYPE::data, + "[MetropolisHastings] The teperature (beta) must be >= 0" ); + if( beta == 0 ) { - Logger::info( - "[Metropolis Hastings] - beta == 0 all move will be " - "accepted - Uniform sampling." ); + Logger::info( "[MetropolisHastings] beta == 0 all move will be " + "accepted - Uniform sampling." ); } - if( b < 1 ) + if( beta < 1 ) { Logger::info( - "[Metropolis Hastings] - beta < 1 moves that increase " + "[MetropolisHastings] beta < 1 moves that increase " "energy are more likely to be accepted - Hot system " "introduce randomness for exploration." ); } - if( b == 1 ) + if( beta == 1 ) { - Logger::info( "[Metropolis Hastings] - beta == 1 default " + Logger::info( "[MetropolisHastings] beta == 1 default " "choice no temperature - only consider energy." ); } - if( b > 1 ) + if( beta > 1 ) { - Logger::info( "[Metropolis Hastings] - beta > 1 moves that " + Logger::info( "[MetropolisHastings] beta > 1 moves that " "increase energy are less likely to be accepted " "- Cold system to ensure convergence but may " "find local minimum randomness." ); } - beta_ = b; + beta_ = beta; } - static double acceptance_prob_helper( double log_accept ) + [[nodiscard]] static double acceptance_prob_helper( double log_accept ) { if( std::isnan( log_accept ) ) + { return 0.0; + } if( log_accept >= 0.0 ) + { return 1.0; + } // prevent exponential overflow constexpr double LOG_MIN = -745.0; if( log_accept < LOG_MIN ) + { return 0.0; + } return std::exp( log_accept ); } private: - const double compute_log_accept( const double deltaU, - const ProposalProbabilities& proposal_probas ) const + [[nodiscard]] double compute_log_accept( + double deltaU, const ProposalProbabilities& proposal_probas ) const { - return -beta_ * deltaU + proposal_probas.transition_probability(); + return ( -beta_ * deltaU ) + + proposal_probas.transition_probability(); } template < typename ApplyMove > @@ -164,7 +170,7 @@ namespace geode ObjectSets< ObjectType >& state, RandomEngine& engine, const double delta_log_energy, - ApplyMove&& apply_move ) const + const ApplyMove& apply_move ) const { const auto& proposed_move = proposal.proposed_move; StepResult< ObjectType > step_result; @@ -175,11 +181,12 @@ namespace geode double log_u = engine.sample_log(); step_result.decision = ( log_u < step_result.log_accept ) - ? MHDecision::Accepted - : MHDecision::Rejected; - if( step_result.decision == MHDecision::Accepted ) + ? MH_DECISION::accepted + : MH_DECISION::rejected; + if( step_result.decision == MH_DECISION::accepted ) + { apply_move( state, proposal ); - + } return step_result; } @@ -189,12 +196,13 @@ namespace geode { const auto new_object = proposal.new_object(); const auto delta_log_energy = - gibbs_energy_.delta_log_add( state, new_object ); + model_.energy().delta_log_add( state, new_object ); return accept_or_reject( proposal, state, engine, delta_log_energy, - []( auto& state, auto& proposal ) { - state.add_object( - std::move( proposal.proposed_move.new_object.value() ), - proposal.set_id, false ); + []( auto& cur_state, auto& accepted_proposal ) { + cur_state.add_object( + std::move( accepted_proposal.proposed_move.new_object + .value() ), + accepted_proposal.set_id, false ); } ); }; @@ -204,10 +212,11 @@ namespace geode { const auto old_object_id = proposal.old_object_id(); const auto delta_log_energy = - gibbs_energy_.delta_log_remove( state, old_object_id ); + model_.energy().delta_log_remove( state, old_object_id ); return accept_or_reject( proposal, state, engine, delta_log_energy, - []( auto& state, auto& proposal ) { - state.remove_free_object( proposal.old_object_id() ); + []( auto& cur_state, auto& accepted_proposal ) { + cur_state.remove_free_object( + accepted_proposal.old_object_id() ); } ); }; @@ -217,18 +226,19 @@ namespace geode { const auto new_object = proposal.new_object(); const auto old_object_id = proposal.old_object_id(); - const auto delta_log_energy = gibbs_energy_.delta_log_change( + const auto delta_log_energy = model_.energy().delta_log_change( state, old_object_id, new_object ); return accept_or_reject( proposal, state, engine, delta_log_energy, - []( auto& state, auto& proposal ) { - state.update_free_object( proposal.old_object_id(), - std::move( - proposal.proposed_move.new_object.value() ) ); + []( auto& cur_state, auto& accepted_proposal ) { + cur_state.update_free_object( + accepted_proposal.old_object_id(), + std::move( accepted_proposal.proposed_move.new_object + .value() ) ); } ); }; private: - GibbsEnergy< ObjectType > gibbs_energy_; + const Model< ObjectType >& model_; std::unique_ptr< ProposalKernel< ObjectType > > proposal_kernel_; double beta_{ 1.0 }; }; diff --git a/include/geode/stochastic/sampling/mcmc/proposal/classical_proposals.hpp b/include/geode/stochastic/sampling/mcmc/proposal/classical_proposals.hpp index c46b0b0..866b882 100644 --- a/include/geode/stochastic/sampling/mcmc/proposal/classical_proposals.hpp +++ b/include/geode/stochastic/sampling/mcmc/proposal/classical_proposals.hpp @@ -32,7 +32,7 @@ namespace geode { template < typename ObjectType > void add_birth_death_change_moves( - std::unique_ptr< geode::ObjectSetSampler< ObjectType > >& sampler, + const geode::ObjectSetSampler< ObjectType >& sampler, geode::ProposalKernel< ObjectType >& kernel, const uuid& set_id, double birth_ratio, @@ -42,18 +42,19 @@ namespace geode const auto total_ratio = birth_ratio + death_ratio; OpenGeodeStochasticStochasticException::check_exception( total_ratio > 0., nullptr, OpenGeodeException::TYPE::data, - "BIRTH+DEATH ratio must be positive" ); + "[add_birth_death_change_moves] Birth + Death ratio must be " + "positive." ); const auto p_birth = birth_ratio / total_ratio; kernel.add_move( set_id, std::make_unique< geode::BirthDeathMove< ObjectType > >( - *sampler, total_ratio, p_birth ) ); + sampler, total_ratio, p_birth ) ); if( change_ratio > 0. ) { kernel.add_move( set_id, std::make_unique< geode::ChangeMove< ObjectType > >( - *sampler, change_ratio ) ); + sampler, change_ratio ) ); } } @@ -80,7 +81,9 @@ namespace geode auto birth_death_prob = birth_prob + death_prob; OpenGeodeStochasticStochasticException::check_exception( birth_death_prob < 1., nullptr, OpenGeodeException::TYPE::data, - "[Proposal Kernel] - changes should be allowed." ); + "[create_birth_death_change_kernel] The probability of birth + " + "death = ", + birth_death_prob, " but should be < 1 to enable change moves." ); auto kernel = std::make_unique< ProposalKernel< ObjectType > >(); kernel->add_move( set_id, std::make_unique< BirthDeathMove< ObjectType > >( sampler, diff --git a/include/geode/stochastic/sampling/mcmc/proposal/moves.hpp b/include/geode/stochastic/sampling/mcmc/proposal/moves.hpp index 82cde0e..718e152 100644 --- a/include/geode/stochastic/sampling/mcmc/proposal/moves.hpp +++ b/include/geode/stochastic/sampling/mcmc/proposal/moves.hpp @@ -54,7 +54,7 @@ namespace geode std::isfinite( log_forward_prob ) && std::isfinite( log_backward_prob ), nullptr, OpenGeodeException::TYPE::data, - "[Proposal Probabilities] Non-finite proposal " + "[ProposalProbabilities] Non-finite proposal " "log-probabilities" ); return log_backward_prob - log_forward_prob; } @@ -119,8 +119,7 @@ namespace geode OpenGeodeStochasticStochasticException::check_exception( proportion_weight_ > 0., nullptr, OpenGeodeException::TYPE::data, - "[Move] - the weight factor for a move should be in higher " - "than 0. (here = ", + "[Move] The weight factor for a move should be > 0. (here = ", proportion_weight_, ")" ); initialize_probability( 1. ); } @@ -190,7 +189,7 @@ namespace geode OpenGeodeStochasticStochasticException::check_exception( birth_ratio_ > 0. && birth_ratio_ < 1., nullptr, OpenGeodeException::TYPE::data, - "[BirthDeathMove]-the ratio of birth over mover should be in " + "[BirthDeathMove] The ratio of birth moves should be in " "]0,1[. (here = ", birth_ratio_, ")" ); log_p_birth_ = std::log( probability * birth_ratio_ ); diff --git a/include/geode/stochastic/models/spatial_model.hpp b/include/geode/stochastic/sampling/mcmc/proposal/object_set_dynamic_config.hpp similarity index 80% rename from include/geode/stochastic/models/spatial_model.hpp rename to include/geode/stochastic/sampling/mcmc/proposal/object_set_dynamic_config.hpp index 1daba84..7d8edb0 100644 --- a/include/geode/stochastic/models/spatial_model.hpp +++ b/include/geode/stochastic/sampling/mcmc/proposal/object_set_dynamic_config.hpp @@ -19,4 +19,20 @@ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE * SOFTWARE. * - */ \ No newline at end of file + */ + +#pragma once + +#include + +namespace geode +{ + struct ObjectSetDynamicsConfig + { + std::string object_set_name; + + double birth_ratio = 1.0; + double death_ratio = 1.0; + double change_ratio = 1.0; + }; +} // namespace geode \ No newline at end of file diff --git a/include/geode/stochastic/sampling/mcmc/proposal/proposal_kernel.hpp b/include/geode/stochastic/sampling/mcmc/proposal/proposal_kernel.hpp index f539b0b..61d0be1 100644 --- a/include/geode/stochastic/sampling/mcmc/proposal/proposal_kernel.hpp +++ b/include/geode/stochastic/sampling/mcmc/proposal/proposal_kernel.hpp @@ -23,6 +23,9 @@ #pragma once +#include +#include + #include #include @@ -38,7 +41,7 @@ namespace geode { OpenGeodeStochasticStochasticException::check_exception( proposed_move.new_object.has_value(), nullptr, - OpenGeodeException::TYPE::data, + OpenGeodeException::TYPE::internal, "[Proposal] Proposal has no new_object" ); return ObjectRef< ObjectType >{ proposed_move.new_object.value(), set_id }; @@ -48,7 +51,7 @@ namespace geode { OpenGeodeStochasticStochasticException::check_exception( proposed_move.old_object_id.has_value(), nullptr, - OpenGeodeException::TYPE::data, + OpenGeodeException::TYPE::internal, "[Proposal] Proposal has no old_object_id" ); return ObjectId{ proposed_move.old_object_id.value(), false, set_id }; @@ -72,7 +75,7 @@ namespace geode { OpenGeodeStochasticStochasticException::check_exception( !set_moves_.empty(), nullptr, OpenGeodeException::TYPE::data, - "[MCMC Proposal Kernel] - no move are defined in the Kernel." ); + "[ProposalKernel] No move are defined in the Kernel." ); auto rnd = engine.sample_uniform( uniform_distribution_closed_ ); for( const auto proba_id : Range{ cumulative_probs_.size() } ) { @@ -86,7 +89,7 @@ namespace geode } throw OpenGeodeStochasticStochasticException{ nullptr, OpenGeodeException::TYPE::internal, - "[MCMC Proposal Kernel]: Should not be reached move pdf is " + "[ProposalKernel]: Should not be reached move pdf is not" "correctly set." }; } @@ -134,8 +137,7 @@ namespace geode OpenGeodeStochasticStochasticException::check_exception( total > GLOBAL_EPSILON, nullptr, OpenGeodeException::TYPE::internal, - "[MCMC Proposal Kernel] - Total " - "probability is zero in Kernel." ); + "[ProposalKernel] - The probability of every moves is zero." ); // Normalize std::transform( probabilities.begin(), probabilities.end(), diff --git a/include/geode/stochastic/sampling/mcmc/simulation_runner.hpp b/include/geode/stochastic/sampling/mcmc/simulation_runner.hpp index dd432db..df4fb8e 100644 --- a/include/geode/stochastic/sampling/mcmc/simulation_runner.hpp +++ b/include/geode/stochastic/sampling/mcmc/simulation_runner.hpp @@ -23,25 +23,31 @@ #pragma once #include -#include -#include #include -#include -#include -#include +#include +#include namespace geode { + namespace detail + { + // NOLINTBEGIN(*-magic-numbers) + constexpr index_t DEFAULT_REALIZATIONS{ 1000 }; + constexpr index_t DEFAULT_SIMULATION_STEPS{ 1000 }; + constexpr index_t DEFAULT_BURN_IN_STEPS{ 1000 }; + // NOLINTEND(*-magic-numbers) + } // namespace detail + struct SimulationConfigurator { - index_t realizations{ 1000 }; - index_t metropolis_hasting_steps{ 1000 }; - index_t burn_in_steps{ 1000 }; + index_t realizations{ detail::DEFAULT_REALIZATIONS }; + index_t metropolis_hasting_steps{ detail::DEFAULT_SIMULATION_STEPS }; + index_t burn_in_steps{ detail::DEFAULT_BURN_IN_STEPS }; std::optional< SimulationPrinterConfigurator > printer{ std::nullopt }; - std::string string() const + [[nodiscard]] std::string string() const { auto message = absl::StrCat( "SimulationConfigurator: " ); absl::StrAppend( &message, "\n\t --> ", realizations, @@ -62,110 +68,78 @@ namespace geode template < typename ObjectType > class SimulationRunner { + OPENGEODE_DISABLE_COPY_AND_MOVE( SimulationRunner ); + public: - SimulationRunner( const SpatialDomain< ObjectType::dim >& domain ) - : domain_( domain ) {}; + SimulationRunner() = delete; + explicit SimulationRunner( SimulationContext< ObjectType >&& context ) + : context_( std::move( context ) ) + { + } virtual ~SimulationRunner() = default; - virtual void initialize() = 0; - - const ObjectSets< ObjectType >& run( + [[nodiscard]] const ObjectSets< ObjectType >& run( RandomEngine& engine, const index_t steps ) { - mh_sampler_->walk( object_sets_, engine, steps ); - return object_sets_; + context_.mh_sampler->walk( *context_.object_sets, engine, steps ); + return *context_.object_sets; } - StatisticsMonitor run( + [[nodiscard]] StatisticsTracker< ObjectType > run( RandomEngine& engine, const SimulationConfigurator& config ) { if( config.burn_in_steps > 0 ) { - mh_sampler_->walk( object_sets_, engine, config.burn_in_steps ); + context_.mh_sampler->walk( + *context_.object_sets, engine, config.burn_in_steps ); } // Initialize monitoring - StatisticsMonitor stats_monitor( energy_terms_collection_.size() ); - std::unique_ptr< SimulationPrinter > printer; + StatisticsTracker< ObjectType > stats_monitor( *context_.model ); + std::unique_ptr< SimulationPrinter< ObjectType > > printer; if( config.printer.has_value() ) { - printer = std::make_unique< SimulationPrinter >( - config.printer.value() ); + printer = std::make_unique< SimulationPrinter< ObjectType > >( + *context_.model, config.printer.value() ); } for( const auto realization : Range{ config.realizations } ) { - mh_sampler_->walk( - object_sets_, engine, config.metropolis_hasting_steps ); + context_.mh_sampler->walk( *context_.object_sets, engine, + config.metropolis_hasting_steps ); - const auto stats = state_statistics(); + const auto stats = + context_.model->compute_statistics( *context_.object_sets ); stats_monitor.add_realization( stats ); if( printer ) { - printer->print_statistics( - stats, model_energy_term_names() ); - printer->print_object_sets( object_sets_, realization ); + printer->print_statistics( stats ); + printer->print_object_sets( + *context_.object_sets, realization ); } } if( printer ) { - printer->print_statistics_summary( - stats_monitor, model_energy_term_names() ); + printer->print_statistics_summary( stats_monitor ); } return stats_monitor; } - const ObjectSets< ObjectType >& state_realization() const - { - return object_sets_; - } - - std::vector< double > state_statistics() const + [[nodiscard]] const ObjectSets< ObjectType >& state_realization() 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; + return *context_.object_sets; } - std::string model_energy_term_names() const + [[nodiscard]] const Model< ObjectType >& model() 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().value_or( term.id().string() ) ); - } - - return absl::StrCat( absl::StrJoin( term_names, " ; " ), "\n" ); + return *context_.model; } - protected: - SpatialDomain< ObjectType::dim > domain_; - 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_; + private: + SimulationContext< ObjectType > context_; }; } // namespace geode \ No newline at end of file diff --git a/include/geode/stochastic/spatial/object_neighborhood.hpp b/include/geode/stochastic/spatial/object_neighborhood.hpp index a3d752f..a960a08 100644 --- a/include/geode/stochastic/spatial/object_neighborhood.hpp +++ b/include/geode/stochastic/spatial/object_neighborhood.hpp @@ -47,6 +47,11 @@ namespace geode return index != other.index || fixed != other.fixed || set_id != other.set_id; } + bool operator<( const ObjectId& other ) const noexcept + { + return index < other.index + && ( set_id < other.set_id || set_id == other.set_id ); + } }; template < index_t dimension > diff --git a/include/geode/stochastic/spatial/object_set.hpp b/include/geode/stochastic/spatial/object_set.hpp index 12c2de4..876888b 100644 --- a/include/geode/stochastic/spatial/object_set.hpp +++ b/include/geode/stochastic/spatial/object_set.hpp @@ -31,14 +31,17 @@ namespace geode template < typename Type > class ObjectSet : public Identifier { + OPENGEODE_DISABLE_COPY( ObjectSet ); + public: ObjectSet() noexcept = default; + ~ObjectSet() = default; ObjectSet( ObjectSet&& ) noexcept = default; ObjectSet& operator=( ObjectSet&& ) noexcept = default; void set_name( std::string_view name ); - const Type& get_fixed_object( index_t index ) const; - const Type& get_free_object( index_t index ) const; + [[nodiscard]] const Type& get_fixed_object( index_t index ) const; + [[nodiscard]] const Type& get_free_object( index_t index ) const; index_t add_fixed_object( Type&& object ); @@ -46,13 +49,13 @@ namespace geode void update_free_object( index_t index, Type&& object ); void remove_free_object( index_t index ); - index_t nb_objects() const; - index_t nb_fixed_objects() const; - index_t nb_free_objects() const; + [[nodiscard]] index_t nb_objects() const; + [[nodiscard]] index_t nb_fixed_objects() const; + [[nodiscard]] index_t nb_free_objects() const; - bool empty() const; + [[nodiscard]] bool empty() const; - std::string string() const; + [[nodiscard]] std::string string() const; private: std::vector< Type > fixed_objects_; diff --git a/include/geode/stochastic/spatial/object_sets.hpp b/include/geode/stochastic/spatial/object_sets.hpp index c6881d4..584ee8c 100644 --- a/include/geode/stochastic/spatial/object_sets.hpp +++ b/include/geode/stochastic/spatial/object_sets.hpp @@ -53,40 +53,49 @@ namespace geode public: ObjectSets() noexcept = default; + ~ObjectSets() noexcept = default; + ObjectSets( ObjectSets&& ) noexcept = default; ObjectSets& operator=( ObjectSets&& ) noexcept = default; - const ObjectSet< Type >& get_set( const uuid& set_id ) const; - const Type& get_object( const ObjectId& object_id ) const; - std::vector< ObjectId > get_all_object() const; - std::vector< ObjectId > get_objects_in_set( const uuid& set_id ) const; + [[nodiscard]] const ObjectSet< Type >& get_set( + const uuid& set_id ) const; + [[nodiscard]] const Type& get_object( const ObjectId& object_id ) const; + [[nodiscard]] std::vector< ObjectId > get_all_object() const; + [[nodiscard]] std::vector< ObjectId > get_objects_in_set( + const uuid& set_id ) const; - index_t nb_sets() const; - index_t nb_objects_in_set( const uuid& set_id ) const; - index_t nb_objects() const; + [[nodiscard]] index_t nb_sets() const; + [[nodiscard]] index_t nb_objects_in_set( const uuid& set_id ) const; + [[nodiscard]] index_t nb_objects() const; uuid add_set( std::string_view name ); ObjectId add_object( Type&& object, const uuid& set_id, bool fixed ); void update_free_object( const ObjectId& object_id, Type&& object ); void remove_free_object( const ObjectId& object_id ); - // Object neighbor search by ObjectId (always excludes self) - std::vector< ObjectId > neighbors( const ObjectId& object_id, - const std::vector< uuid >& targeted_set_ids, - double searching_distance ) const; - // Object neighbor search by arbitrary object (return self if in the - // object_set) - std::vector< ObjectId > neighbors( const Type& object, + [[nodiscard]] std::vector< ObjectId > neighbors( const Type& object, const std::vector< uuid >& targeted_set_ids, - double searching_distance ) const; + double searching_distance, + std::optional< ObjectId > exclude_id ) const; - std::string string() const; + [[nodiscard]] std::string string() const; + + [[nodiscard]] uuid get_set_uuid( std::string_view set_name ) const; + [[nodiscard]] std::vector< uuid > get_set_uuids( + const std::vector< std::string >& set_names ) const; + [[nodiscard]] std::vector< std::pair< uuid, uuid > > get_set_uuid_pairs( + const std::vector< std::pair< std::string, std::string > >& + set_name_pairs ) const; private: - ObjectSet< Type >& get_set( const uuid& set_id ); + [[nodiscard]] ObjectSet< Type >& get_set( const uuid& set_id ); private: - absl::flat_hash_map< uuid, ObjectSet< Type > > sets_; + absl::flat_hash_map< std::string, uuid > name_to_uuid_; + absl::flat_hash_map< uuid, index_t > uuid_to_index_; + std::vector< ObjectSet< Type > > sets_; + ObjectNeighborhood< Type::dim > neighborhood_; }; } // namespace geode diff --git a/include/geode/stochastic/spatial/pairwise_interactions/distance_cutoff.hpp b/include/geode/stochastic/spatial/pairwise_interactions/distance_cutoff.hpp new file mode 100644 index 0000000..684703f --- /dev/null +++ b/include/geode/stochastic/spatial/pairwise_interactions/distance_cutoff.hpp @@ -0,0 +1,83 @@ +/* + * Copyright (c) 2019 - 2026 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 +#include + +namespace geode +{ + /*! + * MinimalDistanceCutoff + * A pairwise interaction that returns 1 if the Euclidean distance + * between objects is within a cutoff radius, otherwise 0. + */ + template < typename Type > + class CenterEuclideanDistanceCutoff : public PairwiseInteraction< Type > + { + public: + explicit CenterEuclideanDistanceCutoff( double cutoff_distance ); + CenterEuclideanDistanceCutoff( double cutoff_distance, + typename PairwiseInteraction< Type >::SCOPE scope ); + + double neighborhood_searching_distance() const override; + + protected: + double compute( const ObjectRef< Type >& object_a, + const ObjectRef< Type >& object_b ) const override; + + private: + double cutoff_distance_{ GLOBAL_EPSILON }; + }; + + /*! + * MinimalDistanceCutoff + * A pairwise interaction that returns 1 if the Minimal distance + * between objects is within a cutoff radius, otherwise 0. + */ + template < typename Type > + class MinimalDistanceCutoff : public PairwiseInteraction< Type > + { + public: + explicit MinimalDistanceCutoff( double cutoff_distance ); + MinimalDistanceCutoff( double cutoff_distance, + typename PairwiseInteraction< Type >::SCOPE scope ); + + double neighborhood_searching_distance() const override; + + protected: + double compute( const ObjectRef< Type >& object_a, + const ObjectRef< Type >& object_b ) const override; + + private: + double cutoff_distance_{ GLOBAL_EPSILON }; + }; +} // namespace geode \ No newline at end of file diff --git a/include/geode/stochastic/spatial/pairwise_interactions.hpp b/include/geode/stochastic/spatial/pairwise_interactions/pairwise_interactions.hpp similarity index 77% rename from include/geode/stochastic/spatial/pairwise_interactions.hpp rename to include/geode/stochastic/spatial/pairwise_interactions/pairwise_interactions.hpp index 8034667..c10113b 100644 --- a/include/geode/stochastic/spatial/pairwise_interactions.hpp +++ b/include/geode/stochastic/spatial/pairwise_interactions/pairwise_interactions.hpp @@ -75,26 +75,4 @@ 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: - explicit EuclideanCutoffInteraction( double cutoff_distance ); - EuclideanCutoffInteraction( double cutoff_distance, - typename PairwiseInteraction< Type >::SCOPE scope ); - - double neighborhood_searching_distance() const override; - - protected: - double compute( const ObjectRef< Type >& object_a, - const ObjectRef< Type >& object_b ) const override; - - private: - double cutoff_distance_{ GLOBAL_EPSILON }; - }; } // namespace geode \ No newline at end of file diff --git a/include/geode/stochastic/spatial/pairwise_interactions/pairwise_interactions_builder.hpp b/include/geode/stochastic/spatial/pairwise_interactions/pairwise_interactions_builder.hpp new file mode 100644 index 0000000..72e1de2 --- /dev/null +++ b/include/geode/stochastic/spatial/pairwise_interactions/pairwise_interactions_builder.hpp @@ -0,0 +1,81 @@ +/* + * Copyright (c) 2019 - 2026 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 + +namespace geode +{ + template < typename ObjectType > + std::unique_ptr< PairwiseInteraction< ObjectType > > + build_pairwise_interaction_impl( + const EuclideanDistanceCutoffConfig& cfg ) + { + return std::make_unique< CenterEuclideanDistanceCutoff< ObjectType > >( + cfg.threshold ); + } + + template < typename ObjectType > + std::unique_ptr< PairwiseInteraction< ObjectType > > + build_pairwise_interaction_impl( + const MinimalDistanceCutoffConfig& cfg ) + { + return std::make_unique< MinimalDistanceCutoff< ObjectType > >( + cfg.threshold ); + } + + template < typename ObjectType, typename NewInteractionConfig > + std::unique_ptr< PairwiseInteraction< ObjectType > > + build_pairwise_interaction_impl( + const NewInteractionConfig& /*unused*/ ) + { + static_assert( sizeof( NewInteractionConfig ) == 0, + "Unsupported PairwiseInteractionConfig type" ); + return nullptr; + } + + template < typename ObjectType > + std::unique_ptr< PairwiseInteraction< ObjectType > > + build_pairwise_interaction_impl( const std::monostate& /*unused*/ ) + { + throw OpenGeodeStochasticStochasticException{ nullptr, + OpenGeodeException::TYPE::data, + "[PairWiseInteractionBuilder] interaction config not initialized" }; + } + + template < typename ObjectType > + std::unique_ptr< PairwiseInteraction< ObjectType > > + build_pairwise_interaction( const PairwiseInteractionConfig& cfg ) + { + return std::visit( + []( auto&& interaction_cfg ) + -> std::unique_ptr< PairwiseInteraction< ObjectType > > { + return build_pairwise_interaction_impl< ObjectType >( + interaction_cfg ); + }, + cfg ); + } +} // namespace geode \ No newline at end of file diff --git a/include/geode/stochastic/sampling/mcmc/energy_terms/intensity_term.hpp b/include/geode/stochastic/spatial/pairwise_interactions/pairwise_interactions_config.hpp similarity index 62% rename from include/geode/stochastic/sampling/mcmc/energy_terms/intensity_term.hpp rename to include/geode/stochastic/spatial/pairwise_interactions/pairwise_interactions_config.hpp index 7cd9dc1..3c3b628 100644 --- a/include/geode/stochastic/sampling/mcmc/energy_terms/intensity_term.hpp +++ b/include/geode/stochastic/spatial/pairwise_interactions/pairwise_interactions_config.hpp @@ -22,28 +22,24 @@ */ #pragma once -#include +#include -#include +#include namespace geode { - FORWARD_DECLARATION_DIMENSION_CLASS( OwnerSegment ); - ALIAS_2D( OwnerSegment ); -} // namespace geode + struct EuclideanDistanceCutoffConfig + { + double threshold; + }; -namespace geode -{ - class opengeode_stochastic_stochastic_api IntensityTerm - : public SingleObjectTerm< OwnerSegment2D, - std::function< double( const OwnerSegment2D&, - const SpatialDomain< OwnerSegment2D::dim >& ) > > + struct MinimalDistanceCutoffConfig { - public: - explicit IntensityTerm( std::string_view name, - double lambda, - std::vector< uuid > targeted_set_ids, - double caracteristic_length, - const SpatialDomain< OwnerSegment2D::dim >& domain ); + double threshold; }; + + using PairwiseInteractionConfig = std::variant< std::monostate, + EuclideanDistanceCutoffConfig, + MinimalDistanceCutoffConfig >; + } // namespace geode \ No newline at end of file diff --git a/include/geode/stochastic/spatial/single_object_features/segment_length_feature.hpp b/include/geode/stochastic/spatial/single_object_features/segment_length_feature.hpp new file mode 100644 index 0000000..d445dc1 --- /dev/null +++ b/include/geode/stochastic/spatial/single_object_features/segment_length_feature.hpp @@ -0,0 +1,42 @@ +/* + * Copyright (c) 2019 - 2026 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 + +namespace geode +{ + class opengeode_stochastic_stochastic_api SegmentLengthInsideBoxFeature + : public SingleObjectFeature< OwnerSegment2D > + { + public: + explicit SegmentLengthInsideBoxFeature( double characteristic_length ); + + [[nodiscard]] double evaluate( const OwnerSegment2D& segment, + const SpatialDomain< 2 >& domain ) const override; + + private: + double inv_length_; + }; +} // namespace geode \ No newline at end of file diff --git a/include/geode/stochastic/spatial/single_object_features/single_object_feature.hpp b/include/geode/stochastic/spatial/single_object_features/single_object_feature.hpp new file mode 100644 index 0000000..67fb5fc --- /dev/null +++ b/include/geode/stochastic/spatial/single_object_features/single_object_feature.hpp @@ -0,0 +1,54 @@ +/* + * Copyright (c) 2019 - 2026 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 + +namespace geode +{ + template < typename ObjectType > + class SingleObjectFeature + { + public: + virtual ~SingleObjectFeature() = default; + + virtual double evaluate( const ObjectType& obj, + const SpatialDomain< ObjectType::dim >& domain ) const = 0; + }; + + template < typename ObjectType > + class ObjectInDomainFeature : public SingleObjectFeature< ObjectType > + { + public: + double evaluate( const ObjectType& obj, + const SpatialDomain< ObjectType::dim >& domain ) const override + { + return SpatialDomainChecker< ObjectType >::is_anchored_in_domain( + domain, obj ) + ? 1.0 + : 0.0; + } + }; + +} // namespace geode \ No newline at end of file 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 new file mode 100644 index 0000000..97c954d --- /dev/null +++ b/include/geode/stochastic/spatial/single_object_features/single_object_feature_builder.hpp @@ -0,0 +1,96 @@ +/* + * Copyright (c) 2019 - 2026 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 + +namespace geode +{ + template < typename ObjectType > + std::unique_ptr< SingleObjectFeature< ObjectType > > + build_single_object_feature_impl( + const ObjectInDomainFeatureConfig& cfg ) + { + geode_unused( cfg ); + return std::make_unique< ObjectInDomainFeature< ObjectType > >(); + }; + + template < typename ObjectType > + std::unique_ptr< SingleObjectFeature< ObjectType > > + build_single_object_feature_impl( + const SegmentLengthInsideBoxFeatureConfig& cfg ) + { + if constexpr( std::is_same_v< ObjectType, OwnerSegment2D > ) + { + return std::make_unique< SegmentLengthInsideBoxFeature >( + cfg.characteristic_length ); + } + else + { + throw OpenGeodeStochasticStochasticException{ nullptr, + OpenGeodeException::TYPE::data, + "[SingleObjectFeatureBuilder] SegmentLengthInsideBoxFeature " + "not valid for this ObjectType" }; + } + } + + template < typename ObjectType, typename NewObjectFeatureConfig > + std::unique_ptr< SingleObjectFeature< ObjectType > > + build_single_object_feature_impl( + const NewObjectFeatureConfig& /*unused*/ ) + { + static_assert( sizeof( NewObjectFeatureConfig ) == 0, + "[SingleObjectFeatureBuilder] Unsupported " + "SingleObjectFeatureConfig type" ); + return nullptr; + } + + template < typename ObjectType > + std::unique_ptr< SingleObjectFeature< ObjectType > > + build_single_object_feature_impl( const std::monostate& /*unused*/ ) + { + throw OpenGeodeStochasticStochasticException{ nullptr, + OpenGeodeException::TYPE::data, + "[SingleObjectFeatureBuilder] object feature " + "config not initialized" }; + } + + template < typename ObjectType > + std::unique_ptr< SingleObjectFeature< ObjectType > > + build_single_object_feature( const SingleObjectFeatureConfig& cfg ) + { + return std::visit( + []( auto&& interaction_cfg ) + -> std::unique_ptr< SingleObjectFeature< ObjectType > > { + return build_single_object_feature_impl< ObjectType >( + interaction_cfg ); + }, + cfg ); + } + +} // namespace geode \ No newline at end of file diff --git a/include/geode/stochastic/spatial/single_object_features/single_object_feature_config.hpp b/include/geode/stochastic/spatial/single_object_features/single_object_feature_config.hpp new file mode 100644 index 0000000..8389f23 --- /dev/null +++ b/include/geode/stochastic/spatial/single_object_features/single_object_feature_config.hpp @@ -0,0 +1,42 @@ +/* + * Copyright (c) 2019 - 2026 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 + +namespace geode +{ + struct ObjectInDomainFeatureConfig + { + }; + + struct SegmentLengthInsideBoxFeatureConfig + { + double characteristic_length; + }; + + using SingleObjectFeatureConfig = std::variant< std::monostate, + ObjectInDomainFeatureConfig, + SegmentLengthInsideBoxFeatureConfig >; + +} // namespace geode \ No newline at end of file diff --git a/include/geode/stochastic/spatial/spatial_domain.hpp b/include/geode/stochastic/spatial/spatial_domain.hpp index 3d8d08d..d1bccef 100644 --- a/include/geode/stochastic/spatial/spatial_domain.hpp +++ b/include/geode/stochastic/spatial/spatial_domain.hpp @@ -23,6 +23,8 @@ #pragma once +#include + #include #include #include @@ -32,9 +34,12 @@ namespace geode template < index_t dimension > class SpatialDomain { + OPENGEODE_DISABLE_COPY_AND_MOVE( SpatialDomain ); + public: - SpatialDomain( - const BoundingBox< dimension >& domain, double buffer_size ) + ~SpatialDomain() = default; + + SpatialDomain( BoundingBox< dimension > domain, double buffer_size ) : domain_{ domain }, buffer_size_{ buffer_size }, extended_domain_{ domain } @@ -42,12 +47,12 @@ namespace geode auto volume = domain_.n_volume(); OpenGeodeStochasticStochasticException::check_exception( volume > 0., nullptr, OpenGeodeException::TYPE::data, - "[SpatialDomain] - Undefined Spatial Domain (volume == ", - volume, ")." ); + "[SpatialDomain] Undefined Spatial Domain (volume == ", volume, + ")." ); OpenGeodeStochasticStochasticException::check_exception( buffer_size_ >= 0.0, nullptr, OpenGeodeException::TYPE::data, - "[SpatialDomain] buffer size must be non-negative ( buffer " - "== ", + "[SpatialDomain] Buffer size must not be < 0 ( buffer " + "= ", buffer_size_, ")" ); if( buffer_size_ != 0. ) { @@ -55,41 +60,49 @@ namespace geode } } - const BoundingBox< dimension > box() const + [[nodiscard]] const BoundingBox< dimension >& box() const { return domain_; } - bool contains( const Point< dimension >& point ) const + [[nodiscard]] bool contains( const Point< dimension >& point ) const { return domain_.contains( point ); } - double n_volume() const + [[nodiscard]] double n_volume() const { return domain_.n_volume(); } - double smallest_length() const + [[nodiscard]] double smallest_length() const { return std::get< 1 >( domain_.smallest_length() ); } - bool extended_contains( const Point< dimension >& point ) const + [[nodiscard]] bool extended_contains( + const Point< dimension >& point ) const { return extended_domain_.contains( point ); } - double extended_n_volume() const + [[nodiscard]] double extended_n_volume() const { return extended_domain_.n_volume(); } - const BoundingBox< dimension > extended_box() const + [[nodiscard]] const BoundingBox< dimension >& extended_box() const { return extended_domain_; } + [[nodiscard]] std::string string() const + { + return absl::StrCat( "Spatial Domain --> center ", domain_.string(), + " extended: ", extended_domain_.string(), + " buffer: ", buffer_size_ ); + } + private: BoundingBox< dimension > domain_; @@ -132,4 +145,33 @@ namespace geode return domain.box().intersects( seg ); } }; + + template < index_t dimension > + struct SpatialDomainConfig + { + Point< dimension > min_point; + Point< dimension > max_point; + double buffer_size{ 0.0 }; + + [[nodiscard]] std::string string() const + { + std::string message = "Domain: "; + absl::StrAppend( + &message, "\n\t --> Min point: ", min_point.string() ); + absl::StrAppend( + &message, "\n\t --> Max point: ", max_point.string() ); + absl::StrAppend( &message, "\n\t --> Buffer size: ", buffer_size ); + return message; + } + }; + + template < index_t dimension > + std::unique_ptr< SpatialDomain< dimension > > build_spatial_domain( + const SpatialDomainConfig< dimension >& config ) + { + BoundingBox< dimension > box{ config.min_point, config.max_point }; + + return std::make_unique< SpatialDomain< dimension > >( + std::move( box ), config.buffer_size ); + } } // namespace geode \ No newline at end of file diff --git a/src/geode/stochastic/CMakeLists.txt b/src/geode/stochastic/CMakeLists.txt index 0394ccc..81d5776 100644 --- a/src/geode/stochastic/CMakeLists.txt +++ b/src/geode/stochastic/CMakeLists.txt @@ -22,27 +22,45 @@ add_geode_library( NAME stochastic FOLDER "geode/stochastic" SOURCES + "applications/fractures.cpp" + "applications/poisson_process.cpp" + "applications/strauss_process.cpp" "spatial/object_set.cpp" "spatial/object_sets.cpp" "spatial/object_neighborhood.cpp" - "spatial/pairwise_interactions.cpp" + "spatial/single_object_features/segment_length_feature.cpp" + "spatial/pairwise_interactions/distance_cutoff.cpp" + "sampling/direct/object_set_sampler/segment_set_sampler.cpp" + "sampling/direct/object_set_sampler/point_set_sampler.cpp" "sampling/direct/ball_sampler.cpp" "sampling/direct/bounding_box_sampler.cpp" "sampling/direct/double_sampler.cpp" "sampling/direct/point_uniform_sampler.cpp" "sampling/direct/segment_uniform_sampler.cpp" + "models/energy_terms/energy_term_builder.cpp" "sampling/distributions.cpp" "sampling/random_engine.cpp" - "sampling/mcmc/energy_terms/intensity_term.cpp" "common.cpp" PUBLIC_HEADERS - "inference/abc_shadow.hpp" - "models/spatial_model.hpp" + "applications/fractures.hpp" + "applications/poisson_process.hpp" + "applications/strauss_process.hpp" + #"inference/abc_shadow.hpp" + "inference/statistics_tools.hpp" + "inference/statistics_tracker.hpp" + "inference/target_statistics.hpp" "spatial/object_set.hpp" "spatial/object_sets.hpp" "spatial/object_neighborhood.hpp" "spatial/object_helpers.hpp" - "spatial/pairwise_interactions.hpp" + "spatial/single_object_features/single_object_feature_builder.hpp" + "spatial/single_object_features/single_object_feature_config.hpp" + "spatial/single_object_features/segment_length_feature.hpp" + "spatial/single_object_features/single_object_feature.hpp" + "spatial/pairwise_interactions/pairwise_interactions_builder.hpp" + "spatial/pairwise_interactions/pairwise_interactions_config.hpp" + "spatial/pairwise_interactions/pairwise_interactions.hpp" + "spatial/pairwise_interactions/distance_cutoff.hpp" "spatial/spatial_domain.hpp" "spatial/details/RTree.hpp" "sampling/direct/object_set_sampler/object_set_sampler.hpp" @@ -53,18 +71,20 @@ add_geode_library( "sampling/direct/double_sampler.hpp" "sampling/direct/point_uniform_sampler.hpp" "sampling/direct/segment_uniform_sampler.hpp" - "sampling/mcmc/helpers/fracture_simulation_runner.hpp" + #"sampling/mcmc/helpers/fracture_simulation_runner.hpp" + "sampling/mcmc/helpers/simulation_context.hpp" "sampling/mcmc/helpers/simulation_printer.hpp" - "sampling/mcmc/helpers/simulation_monitor.hpp" - "sampling/mcmc/energy_terms/density_term.hpp" - "sampling/mcmc/energy_terms/energy_term.hpp" - "sampling/mcmc/energy_terms/intensity_term.hpp" - "sampling/mcmc/energy_terms/pairwise_term.hpp" - "sampling/mcmc/energy_terms/single_object_term.hpp" - "sampling/mcmc/energy_terms/energy_term_collection.hpp" - "sampling/mcmc/energy_terms/gibbs_energy.hpp" + "models/energy_terms/energy_term.hpp" + "models/energy_terms/pairwise_term.hpp" + "models/energy_terms/single_object_term.hpp" + "models/energy_terms/energy_term_config.hpp" + "models/energy_terms/energy_term_builder.hpp" + "models/energy_term_collection.hpp" + "models/gibbs_energy.hpp" + "models/model.hpp" "sampling/mcmc/proposal/classical_proposals.hpp" "sampling/mcmc/proposal/moves.hpp" + "sampling/mcmc/proposal/object_set_dynamic_config.hpp" "sampling/mcmc/proposal/proposal_kernel.hpp" "sampling/mcmc/metropolis_hasting_sampler.hpp" "sampling/mcmc/simulation_runner.hpp" diff --git a/src/geode/stochastic/applications/fractures.cpp b/src/geode/stochastic/applications/fractures.cpp new file mode 100644 index 0000000..9d2c7fb --- /dev/null +++ b/src/geode/stochastic/applications/fractures.cpp @@ -0,0 +1,115 @@ +/* + * Copyright (c) 2019 - 2026 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 +namespace +{ + std::vector< geode::Fracture > build_observed_fractures( + const std::vector< std::array< geode::Point2D, 2 > >& + fracture_extremities ) + { + std::vector< geode::Fracture > fractures; + fractures.reserve( fracture_extremities.size() ); + for( const auto& extremities : fracture_extremities ) + { + fractures.emplace_back( extremities[0], extremities[1] ); + } + return fractures; + } +} // namespace + +namespace geode +{ + using FractureDensityDescription = SingleObjectTermConfig; + using FractureIntensityDescription = SingleObjectTermConfig; + using FractureSpacingDescription = PairwiseTermConfig; + + using FractureSimulationConfig = SimulationContextConfig< Fracture >; + + FractureSimulationContext build_fractures_simulation_context( + const FractureNetworkDescription& fnet_desc ) + { + FractureSimulationConfig simulation_config; + simulation_config.domain = fnet_desc.domain; + + for( const auto& fset_desc : fnet_desc.fracture_sets ) + { + auto& fset = simulation_config.add_set( fset_desc.fset_name ); + + fset.sampler = fset_desc.sampler; + fset.dynamics.birth_ratio = fset_desc.birth_ratio; + fset.dynamics.death_ratio = fset_desc.death_ratio; + fset.dynamics.change_ratio = fset_desc.change_ratio; + fset.fixed_objects = + build_observed_fractures( fset_desc.observed_fractures ); + + FractureDensityDescription density; + density.term_name = fset_desc.density_name(); + density.object_set_names = { fset_desc.fset_name }; + density.lambda = fset_desc.p20; + density.object_feature = ObjectInDomainFeatureConfig{}; + simulation_config.model.terms.emplace_back( std::move( density ) ); + + FractureIntensityDescription intensity; + intensity.term_name = fset_desc.intensity_name(); + intensity.object_set_names = { fset_desc.fset_name }; + intensity.lambda = fset_desc.p21; + constexpr double CARACTERISTIC_LENGTH = 1.0; // mean fracture + // length? + intensity.object_feature = + SegmentLengthInsideBoxFeatureConfig{ CARACTERISTIC_LENGTH }; + simulation_config.model.terms.emplace_back( + std::move( intensity ) ); + + FractureSpacingDescription spacing; + spacing.term_name = fset_desc.spacing_name(); + spacing.object_set_names_interactions = { { fset_desc.fset_name, + fset_desc.fset_name } }; + spacing.gamma = 0.; + spacing.interaction_config = + geode::MinimalDistanceCutoffConfig{ fset_desc.minimal_spacing }; + simulation_config.model.terms.emplace_back( std::move( spacing ) ); + } + + return build_simulation_context( simulation_config ); + } + + std::vector< geode::TargetStatisticConfig > build_fractures_targeted_stat( + const FractureNetworkDescription& description ) + { + std::vector< geode::TargetStatisticConfig > targets; + + for( const auto& set_desc : description.fracture_sets ) + { + if( set_desc.expected_number ) + { + targets.push_back( geode::TargetStatisticConfig{ + set_desc.density_name(), *set_desc.expected_number } ); + } + } + + return targets; + } + +} // namespace geode \ No newline at end of file diff --git a/src/geode/stochastic/applications/poisson_process.cpp b/src/geode/stochastic/applications/poisson_process.cpp new file mode 100644 index 0000000..788711b --- /dev/null +++ b/src/geode/stochastic/applications/poisson_process.cpp @@ -0,0 +1,95 @@ +/* + * Copyright (c) 2019 - 2026 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 + +namespace geode +{ + using PoissonDensityDescription = geode::SingleObjectTermConfig; + + template < typename ObjectType > + SimulationContext< ObjectType > build_poisson_process( + const PoissonProcessDescription< ObjectType >& desc ) + { + SimulationContextConfig< ObjectType > config; + + config.domain = desc.domain; + + for( const auto& set_desc : desc.sets ) + { + auto& set = config.add_set( set_desc.set_name ); + + set.sampler = set_desc.sampler; + + set.dynamics.birth_ratio = set_desc.birth_ratio; + set.dynamics.death_ratio = set_desc.death_ratio; + set.dynamics.change_ratio = set_desc.change_ratio; + + PoissonDensityDescription density; + + density.term_name = set_desc.density_name; + density.object_set_names = { set_desc.set_name }; + density.lambda = set_desc.lambda; + density.object_feature = ObjectInDomainFeatureConfig{}; + + config.model.terms.emplace_back( std::move( density ) ); + } + + return build_simulation_context( config ); + } + + template opengeode_stochastic_stochastic_api SimulationContext< Point2D > + build_poisson_process< Point2D >( + const PoissonProcessDescription< Point2D >& ); + template opengeode_stochastic_stochastic_api SimulationContext< Point3D > + build_poisson_process< Point3D >( + const PoissonProcessDescription< Point3D >& ); + + template < typename ObjectType > + std::vector< geode::TargetStatisticConfig > build_poisson_targeted_stat( + const PoissonProcessDescription< ObjectType >& description ) + { + std::vector< geode::TargetStatisticConfig > targets; + + for( const auto& set_desc : description.sets ) + { + if( !set_desc.expected_nb_objects ) + { + continue; + } + targets.push_back( geode::TargetStatisticConfig{ + set_desc.density_name, *set_desc.expected_nb_objects } ); + } + + return targets; + } + + template opengeode_stochastic_stochastic_api + std::vector< geode::TargetStatisticConfig > + build_poisson_targeted_stat< Point2D >( + const PoissonProcessDescription< Point2D >& ); + template opengeode_stochastic_stochastic_api + std::vector< geode::TargetStatisticConfig > + build_poisson_targeted_stat< Point3D >( + const PoissonProcessDescription< Point3D >& ); +} // namespace geode \ No newline at end of file diff --git a/src/geode/stochastic/applications/strauss_process.cpp b/src/geode/stochastic/applications/strauss_process.cpp new file mode 100644 index 0000000..69f5198 --- /dev/null +++ b/src/geode/stochastic/applications/strauss_process.cpp @@ -0,0 +1,168 @@ +/* + * Copyright (c) 2019 - 2026 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 + +namespace +{ + template < typename ObjectType > + void add_intra_set_interactions( + std::vector< std::pair< std::string, std::string > >& interactions, + const std::vector< std::string >& set_names ) + { + for( const auto name_id : geode::Range{ set_names.size() } ) + { + interactions.emplace_back( set_names[name_id], set_names[name_id] ); + } + } + + template < typename ObjectType > + void add_inter_set_interactions( + std::vector< std::pair< std::string, std::string > >& interactions, + const std::vector< std::string >& set_names ) + { + for( const auto id1 : geode::Range{ set_names.size() } ) + { + for( const auto id2 : geode::Range{ id1 + 1, set_names.size() } ) + { + interactions.emplace_back( set_names[id1], set_names[id2] ); + } + } + } + + template < typename ObjectType > + std::vector< std::pair< std::string, std::string > > + build_interaction_set_names( + const geode::StraussInteractionDescription< ObjectType >& + interaction_desc ) + { + std::vector< std::pair< std::string, std::string > > interactions; + const auto& set_names = interaction_desc.set_names; + if( set_names.empty() ) + { + return interactions; + } + if( interaction_desc.include_intra_set ) + { + add_intra_set_interactions< ObjectType >( interactions, set_names ); + } + if( interaction_desc.include_inter_set ) + { + add_inter_set_interactions< ObjectType >( interactions, set_names ); + } + return interactions; + } +} // namespace + +namespace geode +{ + using DensityDescription = geode::SingleObjectTermConfig; + using InteractionDescription = geode::PairwiseTermConfig; + + template < typename ObjectType > + SimulationContext< ObjectType > build_strauss_process( + const StraussProcessDescription< ObjectType >& desc ) + { + SimulationContextConfig< ObjectType > config; + + config.domain = desc.domain; + + for( const auto& set_desc : desc.sets ) + { + auto& set = config.add_set( set_desc.set_name ); + + set.sampler = set_desc.sampler; + + set.dynamics.birth_ratio = set_desc.birth_ratio; + set.dynamics.death_ratio = set_desc.death_ratio; + set.dynamics.change_ratio = set_desc.change_ratio; + + DensityDescription density; + + density.term_name = set_desc.density_name; + density.object_set_names = { set_desc.set_name }; + density.lambda = set_desc.lambda; + density.object_feature = ObjectInDomainFeatureConfig{}; + config.model.terms.emplace_back( std::move( density ) ); + } + + for( const auto& interaction_desc : desc.interactions ) + { + InteractionDescription interaction; + interaction.term_name = interaction_desc.interaction_name; + interaction.object_set_names_interactions = + std::move( build_interaction_set_names( interaction_desc ) ); + interaction.gamma = interaction_desc.gamma; + interaction.interaction_config = + geode::MinimalDistanceCutoffConfig{ interaction_desc.distance }; + config.model.terms.emplace_back( std::move( interaction ) ); + } + + return build_simulation_context( config ); + } + + template opengeode_stochastic_stochastic_api SimulationContext< Point2D > + build_strauss_process< Point2D >( + const StraussProcessDescription< Point2D >& ); + template opengeode_stochastic_stochastic_api SimulationContext< Point3D > + build_strauss_process< Point3D >( + const StraussProcessDescription< Point3D >& ); + + template < typename ObjectType > + std::vector< geode::TargetStatisticConfig > build_strauss_targeted_stat( + const StraussProcessDescription< ObjectType >& description ) + { + std::vector< geode::TargetStatisticConfig > targets; + + for( const auto& set_desc : description.sets ) + { + if( !set_desc.expected_nb_objects ) + { + continue; + } + targets.push_back( geode::TargetStatisticConfig{ + set_desc.density_name, *set_desc.expected_nb_objects } ); + } + for( const auto& inter_desc : description.interactions ) + { + if( !inter_desc.expected_nb_interactions ) + { + continue; + } + targets.push_back( + geode::TargetStatisticConfig{ inter_desc.interaction_name, + *inter_desc.expected_nb_interactions } ); + } + + return targets; + } + + template opengeode_stochastic_stochastic_api + std::vector< geode::TargetStatisticConfig > + build_strauss_targeted_stat< Point2D >( + const StraussProcessDescription< Point2D >& ); + template opengeode_stochastic_stochastic_api + std::vector< geode::TargetStatisticConfig > + build_strauss_targeted_stat< Point3D >( + const StraussProcessDescription< Point3D >& ); +} // namespace geode \ No newline at end of file diff --git a/src/geode/stochastic/models/energy_terms/energy_term_builder.cpp b/src/geode/stochastic/models/energy_terms/energy_term_builder.cpp new file mode 100644 index 0000000..4f9e157 --- /dev/null +++ b/src/geode/stochastic/models/energy_terms/energy_term_builder.cpp @@ -0,0 +1,89 @@ +/* + * Copyright (c) 2019 - 2026 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 + +namespace +{ + void register_interaction( + const std::pair< geode::uuid, geode::uuid >& interaction, + std::vector< geode::uuid >& interacting_set_ids, + absl::flat_hash_map< geode::uuid, std::vector< geode::uuid > >& + adjacency_map ) + { + const auto& [set1, set2] = interaction; + + adjacency_map[set1].push_back( set2 ); + adjacency_map[set2].push_back( set1 ); + + interacting_set_ids.push_back( set1 ); + interacting_set_ids.push_back( set2 ); + } + + void sort_unique_shrink( std::vector< geode::uuid >& values ) + { + absl::c_sort( values ); + values.erase( + std::unique( values.begin(), values.end() ), values.end() ); + values.shrink_to_fit(); + } + + void deduplicate_adjacency_map( + absl::flat_hash_map< geode::uuid, std::vector< geode::uuid > >& + adjacency_map ) + { + for( auto& [set_uuid, adjacent_set_uuids] : adjacency_map ) + { + geode_unused( set_uuid ); + sort_unique_shrink( adjacent_set_uuids ); + } + } +} // namespace + +namespace geode +{ + 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< uuid > interacting_set_ids; + absl::flat_hash_map< uuid, std::vector< uuid > > + objectset_adjacency_map; + + interacting_set_ids.reserve( 2 * interacting_sets.size() ); + objectset_adjacency_map.reserve( 2 * interacting_sets.size() ); + + for( const auto& interaction : interacting_sets ) + { + register_interaction( + interaction, interacting_set_ids, objectset_adjacency_map ); + } + + sort_unique_shrink( interacting_set_ids ); + deduplicate_adjacency_map( objectset_adjacency_map ); + + return { std::move( interacting_set_ids ), + std::move( objectset_adjacency_map ) }; + } +} // namespace geode diff --git a/src/geode/stochastic/sampling/direct/direction_sampler.cpp b/src/geode/stochastic/sampling/direct/direction_sampler.cpp deleted file mode 100644 index e69de29..0000000 diff --git a/src/geode/stochastic/sampling/direct/double_sampler.cpp b/src/geode/stochastic/sampling/direct/double_sampler.cpp index 10d611a..a23679a 100644 --- a/src/geode/stochastic/sampling/direct/double_sampler.cpp +++ b/src/geode/stochastic/sampling/direct/double_sampler.cpp @@ -20,6 +20,7 @@ * SOFTWARE. * */ +#include #include @@ -52,7 +53,7 @@ namespace geode OpenGeodeStochasticStochasticException::check_exception( desc.min_value && desc.max_value, nullptr, OpenGeodeException::TYPE::data, - "[DoubleSampler] - Incomplete description for " + "[DoubleSampler] Incomplete description for " "Uniform distribution need at least min " "and max values" ); UniformClosed< double > dist; @@ -65,7 +66,7 @@ namespace geode OpenGeodeStochasticStochasticException::check_exception( desc.min_value && desc.max_value, nullptr, OpenGeodeException::TYPE::data, - "[DoubleSampler] - Incomplete description for " + "[DoubleSampler] Incomplete description for " "Uniform distribution need at least min " "and max values" ); UniformClosedOpen< double > dist; @@ -78,7 +79,7 @@ namespace geode OpenGeodeStochasticStochasticException::check_exception( desc.mean && desc.standard_deviation, nullptr, OpenGeodeException::TYPE::data, - "[DoubleSampler] - Incomplete description for " + "[DoubleSampler] Incomplete description for " "Gaussian distribution need at least mean " "and standard deviation values" ); Gaussian dist; @@ -91,7 +92,7 @@ namespace geode OpenGeodeStochasticStochasticException::check_exception( desc.mean && desc.standard_deviation, nullptr, OpenGeodeException::TYPE::data, - "[DoubleSampler] - Incomplete description for " + "[DoubleSampler] Incomplete description for " "Truncated Gaussian distribution need at least mean " "and standard deviation values" ); TruncatedGaussian dist; @@ -107,7 +108,7 @@ namespace geode OpenGeodeStochasticStochasticException::check_exception( desc.mean && desc.kappa, nullptr, OpenGeodeException::TYPE::data, - "[DoubleSampler] - Incomplete description for " + "[DoubleSampler] Incomplete description for " "Von Mises distribution need at least mean " "and concentration (kappa) values" ); VonMises dist; @@ -120,7 +121,7 @@ namespace geode OpenGeodeStochasticStochasticException::check_exception( desc.mean && desc.standard_deviation, nullptr, OpenGeodeException::TYPE::data, - "[DoubleSampler] - Incomplete description for " + "[DoubleSampler] Incomplete description for " "TruncatedLogNormal distribution need mean " "and standard deviation values of the underlying " "normal distribution." ); @@ -136,7 +137,7 @@ namespace geode OpenGeodeStochasticStochasticException::check_exception( desc.alpha.has_value(), nullptr, OpenGeodeException::TYPE::data, - "[DoubleSampler] - Incomplete description for " + "[DoubleSampler] Incomplete description for " "TruncatedPowerLaw distribution need power law " "exponent (alpha)." ); TruncatedPowerLaw dist; @@ -155,7 +156,8 @@ namespace geode { throw OpenGeodeStochasticStochasticException{ nullptr, OpenGeodeException::TYPE::data, - "Unknown distribution type: ", desc.distribution_type.get() }; + "[DoubleSampler] Unknown distribution type: ", + desc.distribution_type.get() }; } return it->second( desc ); } diff --git a/src/geode/stochastic/sampling/direct/object_set_sampler/point_set_sampler.cpp b/src/geode/stochastic/sampling/direct/object_set_sampler/point_set_sampler.cpp new file mode 100644 index 0000000..fc14962 --- /dev/null +++ b/src/geode/stochastic/sampling/direct/object_set_sampler/point_set_sampler.cpp @@ -0,0 +1,126 @@ +/* + * Copyright (c) 2019 - 2026 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 + +namespace geode +{ + template < index_t dimension > + UniformPointSetSampler< dimension >::UniformPointSetSampler( + const SpatialDomain< dimension >& domain, + const ObjectSamplerConfig< Point< dimension > >& config ) + : ObjectSetSampler< Point< dimension > >{}, + domain_{ domain }, + step_move_( define_step_for_move( config.move_ratio ) ) + { + auto volume = domain_.extended_n_volume(); + OpenGeodeStochasticStochasticException::check_exception( volume != 0., + nullptr, OpenGeodeException::TYPE::data, + "[UniformPointSetSampler] Undefined Extended Bounding " + "Box (volume ==0)." ); + this->set_log_pdf( -std::log( volume ) ); + OpenGeodeStochasticStochasticException::check_exception( + step_move_ > 0., nullptr, OpenGeodeException::TYPE::data, + "[UniformPointSetSampler] Undefined step length for move " + "(value == ", + step_move_, ")." ); + } + template < index_t dimension > + Point< dimension > UniformPointSetSampler< dimension >::sample( + RandomEngine& engine ) const + { + return PointUniformSampler::sample< dimension >( + engine, domain_.extended_box() ); + } + + template < index_t dimension > + Point< dimension > UniformPointSetSampler< dimension >::change( + const Point< dimension >& obj, RandomEngine& engine ) const + { + geode::Sphere< dimension > ball{ obj, step_move_ }; + + auto new_point = + PointUniformSampler::sample< dimension >( engine, ball ); + constexpr index_t MAX_TRY{ 100 }; + for( const auto try_id : geode::Range{ MAX_TRY } ) + { + geode_unused( try_id ); + if( domain_.extended_contains( new_point ) ) + { + return new_point; + } + new_point = + PointUniformSampler::sample< dimension >( engine, ball ); + } + throw OpenGeodeStochasticStochasticException{ nullptr, + OpenGeodeException::TYPE::internal, + "[UniformPointSetSampler] Cannot find a point in the " + "extended domain" }; + } + + template < index_t dimension > + double UniformPointSetSampler< dimension >::define_step_for_move( + double ratio ) + { + return ratio * domain_.smallest_length(); + } + + template < index_t dimension > + bool UniformPointSetSampler< dimension >::is_valid_object( + const Point< dimension >& obj ) const + { + return domain_.extended_contains( obj ); + } + + template class opengeode_stochastic_stochastic_api + UniformPointSetSampler< 2 >; + template class opengeode_stochastic_stochastic_api + UniformPointSetSampler< 3 >; + + template <> + std::unique_ptr< ObjectSetSampler< Point2D > > + opengeode_stochastic_stochastic_api build_objectset_sampler< Point2D >( + const SpatialDomain< 2 >& domain, + const ObjectSamplerConfig< Point2D >& config ) + { + return std::make_unique< UniformPointSetSampler< 2 > >( + domain, config ); + } + + template <> + opengeode_stochastic_stochastic_api + std::unique_ptr< ObjectSetSampler< Point3D > > + build_objectset_sampler< Point3D >( const SpatialDomain< 3 >& domain, + const ObjectSamplerConfig< Point3D >& config ) + { + return std::make_unique< UniformPointSetSampler< 3 > >( + domain, config ); + } + +} // namespace geode \ No newline at end of file diff --git a/src/geode/stochastic/sampling/direct/object_set_sampler/segment_set_sampler.cpp b/src/geode/stochastic/sampling/direct/object_set_sampler/segment_set_sampler.cpp new file mode 100644 index 0000000..6d3f659 --- /dev/null +++ b/src/geode/stochastic/sampling/direct/object_set_sampler/segment_set_sampler.cpp @@ -0,0 +1,107 @@ +/* + * Copyright (c) 2019 - 2026 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 geode +{ + UniformSegmentSetSampler::UniformSegmentSetSampler( + const SpatialDomain< 2 >& domain, + const ObjectSamplerConfig< OwnerSegment2D >& config ) + : domain_{ domain }, + length_{ DoubleSampler::create_distribution( config.length ) }, + azimuth_{ DoubleSampler::create_distribution( config.azimuth ) }, + move_ratio_{ config.move_ratio } + { + auto volume = domain_.extended_n_volume(); + OpenGeodeStochasticStochasticException::check_exception( volume != 0., + nullptr, OpenGeodeException::TYPE::data, + "[UniformSegmentSetSampler] Undefined Extended Bounding " + "Box (volume ==0)." ); + this->set_log_pdf( -std::log( volume ) ); + } + + OwnerSegment2D UniformSegmentSetSampler::sample( + RandomEngine& engine ) const + { + auto seg = SegmentUniformSampler::sample( + engine, domain_.extended_box(), length_, azimuth_ ); + return seg; + } + + OwnerSegment2D UniformSegmentSetSampler::change( + const OwnerSegment2D& obj, RandomEngine& engine ) const + { + const auto& extremities = obj.vertices(); + const auto current = + static_cast< local_index_t >( engine.sample_bernoulli( 0.5 ) ); + const auto other = 1 - current; + + geode::Sphere< 2 > ball{ extremities[current], + move_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 } ) + { + geode_unused( try_id ); + if( domain_.extended_contains( new_point ) + || domain_.extended_contains( extremities[other] ) ) + { + OwnerSegment2D new_segment{ obj }; + new_segment.set_point( current, new_point ); + return new_segment; + } + new_point = PointUniformSampler::sample< 2 >( engine, ball ); + } + throw OpenGeodeStochasticStochasticException{ nullptr, + OpenGeodeException::TYPE::internal, + "[SegmentSetSampler] - Cannot find a point in the box" }; + return obj; + } + + bool UniformSegmentSetSampler::is_valid_object( + const OwnerSegment2D& obj ) const + { + const auto& extremities = obj.vertices(); + return domain_.extended_contains( extremities[0] ) + || domain_.extended_contains( extremities[1] ); + } + + template <> + std::unique_ptr< ObjectSetSampler< OwnerSegment2D > > + build_objectset_sampler( const SpatialDomain< 2 >& domain, + const ObjectSamplerConfig< OwnerSegment2D >& config ) + { + return std::make_unique< UniformSegmentSetSampler >( domain, config ); + } + +} // namespace geode \ No newline at end of file diff --git a/src/geode/stochastic/sampling/random_engine.cpp b/src/geode/stochastic/sampling/random_engine.cpp index 7220625..4511798 100644 --- a/src/geode/stochastic/sampling/random_engine.cpp +++ b/src/geode/stochastic/sampling/random_engine.cpp @@ -52,7 +52,7 @@ namespace geode::OpenGeodeStochasticStochasticException::check_exception( p >= 0.0 && p <= 1.0, nullptr, geode::OpenGeodeException::TYPE::data, - "[normal_quantile] - p must be in (0,1). Check the consistencies " + "[normal_quantile] p must be in (0,1). Check the consistencies " "between min,mean and max values." ); static const double a1 = -3.969683028665376e+01; @@ -150,8 +150,8 @@ namespace geode OpenGeodeStochasticStochasticException::check_exception( law.min_value <= law.max_value, nullptr, OpenGeodeException::TYPE::data, - "[Uniform sampling] - Wrong range ", law.min_value, - " is not <= than ", law.max_value, "." ); + "[RandomEngine] Uniform sampling with wrong range , ", + law.min_value, " is not <= than ", law.max_value, "." ); return absl::Uniform( absl::IntervalClosed, rand_gen_, law.min_value, law.max_value ); } @@ -162,8 +162,8 @@ namespace geode OpenGeodeStochasticStochasticException::check_exception( law.min_value < law.max_value, nullptr, OpenGeodeException::TYPE::data, - "[Uniform sampling] - Wrong range ", law.min_value, - " is not < than ", law.max_value, "." ); + "[RandomEngine] Uniform sampling with wrong range, ", + law.min_value, " is not < than ", law.max_value, "." ); return absl::Uniform( absl::IntervalClosedOpen, rand_gen_, law.min_value, law.max_value ); } @@ -175,7 +175,7 @@ namespace geode && std::isfinite( law.standard_deviation ) && std::isfinite( law.mean ), nullptr, OpenGeodeException::TYPE::data, - "[Gaussian sampling] - Infinite " + "[RandomEngine] Gaussian sampling with infinite " "parameters or negative standard deviation N(", law.mean, law.standard_deviation, ")." ); return absl::gaussian_distribution< double >( @@ -189,27 +189,28 @@ namespace geode && std::isfinite( law.standard_deviation ) && std::isfinite( law.mean ), nullptr, OpenGeodeException::TYPE::data, - "[Gaussian sampling] - Infinite parameters or " - "negative standard deviation N(", + "[RandomEngine] Gaussian sampling with infinite " + "parameters or negative standard deviation N(", law.mean, ",", law.standard_deviation, ")." ); - const double max = law.max_value.value_or( + const auto max = law.max_value.value_or( std::numeric_limits< double >::infinity() ); - const double min = law.min_value.value_or( + const auto min = law.min_value.value_or( -std::numeric_limits< double >::infinity() ); OpenGeodeStochasticStochasticException::check_exception( min < max, nullptr, OpenGeodeException::TYPE::data, - "[Gaussian sampling] - Wrong truncation range ", min, - " is not < than ", max, "." ); + "[RandomEngine] Gaussian sampling truncation with wrong " + "range, ", + min, " is not < than ", max, "." ); // Standardize bounds - const double alpha = ( min - law.mean ) / law.standard_deviation; - const double beta = ( max - law.mean ) / law.standard_deviation; + const auto alpha = ( min - law.mean ) / law.standard_deviation; + const auto beta = ( max - law.mean ) / law.standard_deviation; // Compute CDF of bounds, handling infinite alpha/beta - double F_min = std::isfinite( alpha ) ? normal_cdf( alpha ) : 0.0; - double F_max = std::isfinite( beta ) ? normal_cdf( beta ) : 1.0; + auto F_min = std::isfinite( alpha ) ? normal_cdf( alpha ) : 0.0; + auto F_max = std::isfinite( beta ) ? normal_cdf( beta ) : 1.0; // Clamp to avoid exact 0 or 1 (normal_quantile cannot handle them) F_min = std::max( F_min, geode::GLOBAL_EPSILON ); @@ -217,12 +218,12 @@ namespace geode OpenGeodeStochasticStochasticException::check_exception( F_min < F_max, nullptr, OpenGeodeException::TYPE::data, - "[Gaussian sampling] - truncation " + "[RandomEngine] Gaussian sampling truncation " "range is extreme please check inputs" ); // Sample uniform in [F_min, F_max] std::uniform_real_distribution< double > uniform( F_min, F_max ); - const double u = uniform( rand_gen_ ); + const auto u = uniform( rand_gen_ ); // Map through inverse CDF return law.mean + law.standard_deviation * normal_quantile( u ); @@ -233,8 +234,9 @@ namespace geode OpenGeodeStochasticStochasticException::check_exception( law.concentration >= 0.0 && std::isfinite( law.mean ), nullptr, OpenGeodeException::TYPE::data, - "[VonMises sampling] - Invalid parameters: mean=", law.mean, - ", concentration=", law.concentration, "." ); + "[RandomEngine] VonMises sampling with invalid parameters: " + "mean=", + law.mean, ", concentration=", law.concentration, "." ); // Uniform approximation for very small concentration (nearly // uniform) @@ -250,33 +252,33 @@ namespace geode if( law.concentration > LARGE_KAPPA ) { // Variance of approximate normal around mean - const double stddev = 1.0 / std::sqrt( law.concentration ); + const auto stddev = 1.0 / std::sqrt( law.concentration ); std::normal_distribution< double > normal_dist( law.mean, stddev ); - double theta = normal_dist( rand_gen_ ); + auto theta = normal_dist( rand_gen_ ); // Wrap to [0, 2π) return std::fmod( theta + 2.0 * M_PI, 2.0 * M_PI ); } // Best & Fisher (1979) rejection algorithm for moderate // concentration - const double a = + const auto a = 1.0 + std::sqrt( 1.0 + 4.0 * law.concentration * law.concentration ); - const double b = + const auto b = ( a - std::sqrt( 2.0 * a ) ) / ( 2.0 * law.concentration ); - const double r = ( 1.0 + b * b ) / ( 2.0 * b ); + const auto r = ( 1.0 + b * b ) / ( 2.0 * b ); double theta; UniformClosed< double > uniform_dist; while( true ) { - double u1 = sample_uniform( uniform_dist ); - double z = std::cos( M_PI * u1 ); - double f = ( 1.0 + r * std::abs( z ) ) / ( r + std::abs( z ) ); - double c = law.concentration * ( r - f ); - double u2 = sample_uniform( uniform_dist ); + auto u1 = sample_uniform( uniform_dist ); + auto z = std::cos( M_PI * u1 ); + auto f = ( 1.0 + r * std::fabs( z ) ) / ( r + std::fabs( z ) ); + auto c = law.concentration * ( r - f ); + auto u2 = sample_uniform( uniform_dist ); if( u2 < c * ( 2.0 - c ) || u2 <= c * std::exp( 1.0 - c ) ) { @@ -302,46 +304,47 @@ namespace geode && std::isfinite( law.standard_deviation ) && std::isfinite( law.mean ), nullptr, OpenGeodeException::TYPE::data, - "[Truncated LogNormal sampling] - Infinite parameters or " - "negative standard deviation N(", + "[RandomEngine] Truncated LogNormal sampling with infinite " + "parameters or negative standard deviation N(", law.mean, ", ", law.standard_deviation, ")." ); // Determine min and max, respecting optional values - const double min_val = law.min_value.value_or( 0.0 ); - const double max_val = law.max_value.value_or( + const auto min_val = law.min_value.value_or( 0.0 ); + const auto max_val = law.max_value.value_or( std::numeric_limits< double >::infinity() ); OpenGeodeStochasticStochasticException::check_exception( min_val < max_val, nullptr, OpenGeodeException::TYPE::data, - "[Truncated LogNormal sampling] - Wrong truncation range ", + "[RandomEngine] Truncated LogNormal sampling with wrong " + "truncation range ", min_val, " is not < than ", max_val, "." ); // Transform to standard normal space - const double zmin = + const auto zmin = ( std::log( min_val ) - law.mean ) / law.standard_deviation; - const double zmax = + const auto zmax = ( std::isfinite( max_val ) ? ( std::log( max_val ) - law.mean ) / law.standard_deviation : std::numeric_limits< double >::infinity() ); // Compute CDF bounds, handling infinite zmin/zmax - double F_min = - std::max( normal_cdf( zmin ), geode::GLOBAL_EPSILON ); - double F_max = + auto F_min = std::max( normal_cdf( zmin ), geode::GLOBAL_EPSILON ); + auto F_max = std::min( normal_cdf( zmax ), 1.0 - geode::GLOBAL_EPSILON ); OpenGeodeStochasticStochasticException::check_exception( F_min < F_max, nullptr, OpenGeodeException::TYPE::data, - "[Truncated LogNormal sampling] - truncation " - "range is extreme please chack inputs" ); + "[RandomEngine] Truncated LogNormal sampling with extreme " + "truncation " + "range, please chack inputs" ); // Sample uniform in [Fmin, Fmax] std::uniform_real_distribution< double > uniform( F_min, F_max ); - const double u = uniform( rand_gen_ ); + const auto u = uniform( rand_gen_ ); // Map through inverse CDF and exponentiate - const double z = normal_quantile( u ); + const auto z = normal_quantile( u ); return std::exp( law.mean + law.standard_deviation * z ); } @@ -349,30 +352,32 @@ namespace geode { OpenGeodeStochasticStochasticException::check_exception( law.alpha > 0, nullptr, OpenGeodeException::TYPE::data, - "Power-law exponent alpha must be > 0" ); + "[RandomEngine] Power-law sampling with wrong exponent value " + "(alpha should be > 0), gets ", + law.alpha ); // Set bounds - const double xmin = law.min_value.value_or( + const auto xmin = law.min_value.value_or( geode::GLOBAL_EPSILON ); // default 1.0 if unspecified - const double xmax = law.max_value.value_or( + const auto xmax = law.max_value.value_or( std::numeric_limits< double >::infinity() ); OpenGeodeStochasticStochasticException::check_exception( xmin < xmax, nullptr, OpenGeodeException::TYPE::data, - "Truncated power-law: min >= max" ); + "[RandomEngine] Wrong power-law trucation, ", xmin, + "si not < than", xmax ); // Sample uniform std::uniform_real_distribution< double > uniform( 0.0, 1.0 ); - const double u = uniform( rand_gen_ ); + const auto u = uniform( rand_gen_ ); // Inverse CDF - if( std::abs( law.alpha - 1.0 ) > geode::GLOBAL_EPSILON ) + if( std::fabs( law.alpha - 1.0 ) > geode::GLOBAL_EPSILON ) { - double xmin_pow = std::pow( xmin, 1.0 - law.alpha ); - double xmax_pow = - std::isfinite( xmax ) - ? std::pow( xmax, 1.0 - law.alpha ) - : std::numeric_limits< double >::infinity(); + auto xmin_pow = std::pow( xmin, 1.0 - law.alpha ); + auto xmax_pow = std::isfinite( xmax ) + ? std::pow( xmax, 1.0 - law.alpha ) + : std::numeric_limits< double >::infinity(); if( !std::isfinite( xmax_pow ) ) { return std::pow( xmin_pow + u, 1.0 / ( 1.0 - law.alpha ) ); @@ -383,7 +388,7 @@ namespace geode else { // alpha == 1 - const double xmax_eff = + const auto xmax_eff = std::isfinite( xmax ) ? xmax : xmin * geode::GLOBAL_EPSILON; return xmin * std::pow( xmax_eff / xmin, u ); } @@ -400,7 +405,7 @@ namespace geode OpenGeodeStochasticStochasticException::check_exception( probability_of_success >= 0. && probability_of_success <= 1.0, nullptr, OpenGeodeException::TYPE::data, - "Bernoulli sampling cannot be done since ", + "[RandomEngine] Bernoulli sampling cannot be done since ", probability_of_success, " is not in [0.,1.]." ); return absl::bernoulli_distribution( probability_of_success )( rand_gen_ ); diff --git a/src/geode/stochastic/spatial/object_set.cpp b/src/geode/stochastic/spatial/object_set.cpp index 21ce69a..e5afdaa 100644 --- a/src/geode/stochastic/spatial/object_set.cpp +++ b/src/geode/stochastic/spatial/object_set.cpp @@ -44,7 +44,7 @@ namespace geode OpenGeodeStochasticStochasticException::check_exception( index < fixed_objects_.size(), nullptr, OpenGeodeException::TYPE::data, - "[ObjectSet] - index for fixed object out of range." ); + "[ObjectSet] Index for fixed object out of range." ); return fixed_objects_[index]; } @@ -78,7 +78,7 @@ namespace geode OpenGeodeStochasticStochasticException::check_exception( index < free_objects_.size(), nullptr, OpenGeodeException::TYPE::data, - "[ObjectSet] - free object index out of range." ); + "[ObjectSet] Index for free object out of range." ); free_objects_[index] = std::move( object ); } @@ -88,7 +88,7 @@ namespace geode const index_t last = free_objects_.size() - 1; OpenGeodeStochasticStochasticException::check_exception( index <= last, nullptr, OpenGeodeException::TYPE::data, - "[ObjectSet] - free object index out of range." ); + "[ObjectSet] index for free object out of range." ); if( index != last ) { std::swap( free_objects_[index], free_objects_[last] ); diff --git a/src/geode/stochastic/spatial/object_sets.cpp b/src/geode/stochastic/spatial/object_sets.cpp index 4a540e0..c471e3b 100644 --- a/src/geode/stochastic/spatial/object_sets.cpp +++ b/src/geode/stochastic/spatial/object_sets.cpp @@ -1,3 +1,25 @@ +/* + * Copyright (c) 2019 - 2026 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 @@ -11,11 +33,11 @@ namespace geode const ObjectSet< Type >& ObjectSets< Type >::get_set( const uuid& set_id ) const { - auto it = sets_.find( set_id ); + auto it = uuid_to_index_.find( set_id ); OpenGeodeStochasticStochasticException::check_exception( - it != sets_.end(), nullptr, OpenGeodeException::TYPE::data, - "[ObjectSet] - group (", set_id.string(), ") is not defined." ); - return it->second; + it != uuid_to_index_.end(), nullptr, OpenGeodeException::TYPE::data, + "[ObjectSets] Group (", set_id.string(), ") is not defined." ); + return sets_[it->second]; } template < typename Type > @@ -34,8 +56,9 @@ namespace geode { std::vector< ObjectId > result; result.reserve( nb_objects() ); - for( const auto& [set_id, objs] : sets_ ) + for( const auto& objs : sets_ ) { + auto set_id = objs.id(); for( const auto obj_id : geode::Range{ objs.nb_fixed_objects() } ) { result.push_back( { obj_id, true, set_id } ); @@ -82,9 +105,8 @@ namespace geode index_t ObjectSets< Type >::nb_objects() const { index_t nb_objects{ 0 }; - for( const auto& [set_id, objs] : sets_ ) + for( const auto& objs : sets_ ) { - geode_unused( set_id ); nb_objects += objs.nb_objects(); } return nb_objects; @@ -93,14 +115,23 @@ namespace geode template < typename Type > uuid ObjectSets< Type >::add_set( std::string_view name ) { - ObjectSet< Type > new_set; + auto set_index = sets_.size(); + auto& new_set = sets_.emplace_back( ObjectSet< Type >{} ); + const auto set_uuid = new_set.id(); + 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 ) ); - OpenGeodeStochasticStochasticException::check_exception( inserted, - new_set_id, OpenGeodeException::TYPE::data, "[ObjectSet]- group (", - new_set_id.string(), ") already exists." ); - return new_set_id; + auto [it_set_name, set_uuid_inserted] = + name_to_uuid_.emplace( name, set_uuid ); + OpenGeodeStochasticStochasticException::check_exception( + set_uuid_inserted, nullptr, OpenGeodeException::TYPE::data, + "[ObjectSet]- Group named ", name, " already exists." ); + + auto [it_set_uuid, set_index_inserted] = + uuid_to_index_.emplace( set_uuid, set_index ); + OpenGeodeStochasticStochasticException::check_exception( + set_index_inserted, nullptr, OpenGeodeException::TYPE::data, + "[ObjectSets] Group (", set_uuid.string(), ") already exists." ); + return set_uuid; } template < typename Type > @@ -129,12 +160,12 @@ namespace geode { OpenGeodeStochasticStochasticException::check_exception( !old_object_id.fixed, nullptr, OpenGeodeException::TYPE::data, - "[ObjectSet]- cannot modify fixed object." ); + "[ObjectSets]- Cannot modify a fixed object." ); auto& set = get_set( old_object_id.set_id ); OpenGeodeStochasticStochasticException::check_exception( old_object_id.index < set.nb_objects(), nullptr, - OpenGeodeException::TYPE::data, - "[ObjectSet]- index of object to update out of range." ); + OpenGeodeException::TYPE::internal, + "[ObjectSets] Index of an object to change is out of range." ); auto old_box = object_bounding_box( get_object( old_object_id ) ); auto new_box = object_bounding_box( new_object ); neighborhood_.update( old_box, new_box, old_object_id ); @@ -146,8 +177,8 @@ namespace geode { auto& set = get_set( object_id.set_id ); OpenGeodeStochasticStochasticException::check_exception( - !object_id.fixed, nullptr, OpenGeodeException::TYPE::data, - "[ObjectSet]- Cannot remove fixed object." ); + !object_id.fixed, nullptr, OpenGeodeException::TYPE::internal, + "[ObjectSets] Cannot remove fixed object." ); const auto& obj_to_remove = get_object( object_id ); neighborhood_.remove( object_bounding_box( obj_to_remove ), object_id ); @@ -164,27 +195,16 @@ namespace geode set.remove_free_object( object_id.index ); } - template < typename Type > - std::vector< ObjectId > ObjectSets< Type >::neighbors( - const ObjectId& object_id, - const std::vector< uuid >& targeted_set_ids, - double searching_distance ) const - { - auto box = object_bounding_box( get_object( object_id ) ); - box.extends( searching_distance * 2. ); - return neighborhood_.get_all_neighbor_ids( - box, targeted_set_ids, object_id ); - } - template < typename Type > std::vector< ObjectId > ObjectSets< Type >::neighbors( const Type& object, const std::vector< uuid >& targeted_set_ids, - double searching_distance ) const + double searching_distance, + std::optional< ObjectId > excluded_id ) const { auto box = object_bounding_box( object ); box.extends( searching_distance * 2. ); return neighborhood_.get_all_neighbor_ids( - box, targeted_set_ids, std::nullopt ); + box, targeted_set_ids, excluded_id ); } template < typename Type > @@ -194,12 +214,56 @@ namespace geode static_cast< const ObjectSets* >( this )->get_set( set_id ) ); } + template < typename Type > + uuid ObjectSets< Type >::get_set_uuid( const std::string_view name ) const + { + if( auto set_uuid = name_to_uuid_.find( name ); + set_uuid != name_to_uuid_.end() ) + { + return set_uuid->second; + } + throw OpenGeodeStochasticStochasticException{ nullptr, + OpenGeodeException::TYPE::data, + "[ObjectSets] ObjectSet uuid accessor of group named ", name, + " does not exist." }; + } + + template < typename Type > + std::vector< uuid > ObjectSets< Type >::get_set_uuids( + const std::vector< std::string >& set_names ) const + { + std::vector< geode::uuid > uuids; + uuids.reserve( set_names.size() ); + + for( const auto& name : set_names ) + { + uuids.emplace_back( get_set_uuid( name ) ); + } + return uuids; + } + + template < typename Type > + std::vector< std::pair< uuid, uuid > > + ObjectSets< Type >::get_set_uuid_pairs( + const std::vector< std::pair< std::string, std::string > >& + set_name_pairs ) const + { + std::vector< std::pair< uuid, uuid > > result; + result.reserve( set_name_pairs.size() ); + + for( const auto& [name1, name2] : set_name_pairs ) + { + result.emplace_back( get_set_uuid( name1 ), get_set_uuid( name2 ) ); + } + return result; + } + template < typename Type > std::string ObjectSets< Type >::string() const { auto message = absl::StrCat( "ObjectSets with ", nb_objects(), " objects in, ", nb_sets(), " sets" ); - for( const auto& [set_id, objs] : sets_ ) + for( const auto& objs : sets_ ) { absl::StrAppend( &message, "\n\t --> ", objs.string() ); } diff --git a/src/geode/stochastic/spatial/pairwise_interactions.cpp b/src/geode/stochastic/spatial/pairwise_interactions.cpp deleted file mode 100644 index 1b76b20..0000000 --- a/src/geode/stochastic/spatial/pairwise_interactions.cpp +++ /dev/null @@ -1,86 +0,0 @@ -/* - * Copyright (c) 2019 - 2026 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/src/geode/stochastic/spatial/pairwise_interactions/distance_cutoff.cpp b/src/geode/stochastic/spatial/pairwise_interactions/distance_cutoff.cpp new file mode 100644 index 0000000..913bd25 --- /dev/null +++ b/src/geode/stochastic/spatial/pairwise_interactions/distance_cutoff.cpp @@ -0,0 +1,144 @@ +/* + * Copyright (c) 2019 - 2026 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 + +namespace +{ + template < geode::index_t dimension > + double compute_center_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_center_euclidean_distance( + const geode::Segment< dimension >& seg1, + const geode::Segment< dimension >& seg2 ) + { + return geode::point_point_distance( + seg1.barycenter(), seg2.barycenter() ); + } + template < geode::index_t dimension > + double compute_min_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_min_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 > + CenterEuclideanDistanceCutoff< Type >::CenterEuclideanDistanceCutoff( + double cutoff_distance ) + : PairwiseInteraction< Type >(), cutoff_distance_( cutoff_distance ) + { + } + + template < typename Type > + CenterEuclideanDistanceCutoff< Type >::CenterEuclideanDistanceCutoff( + double cutoff_distance, + typename PairwiseInteraction< Type >::SCOPE scope ) + : PairwiseInteraction< Type >( scope ), + cutoff_distance_( cutoff_distance ) + { + } + + template < typename Type > + double + CenterEuclideanDistanceCutoff< Type >::neighborhood_searching_distance() + const + { + return cutoff_distance_; + } + + template < typename Type > + double CenterEuclideanDistanceCutoff< Type >::compute( + const ObjectRef< Type >& object_a, + const ObjectRef< Type >& object_b ) const + { + auto dist = compute_center_euclidean_distance< Type::dim >( + object_a.object, object_b.object ); + return dist <= cutoff_distance_ ? 1.0 : 0.0; + } + + template class opengeode_stochastic_stochastic_api + CenterEuclideanDistanceCutoff< Point< 2 > >; + template class opengeode_stochastic_stochastic_api + CenterEuclideanDistanceCutoff< Point< 3 > >; + + template class opengeode_stochastic_stochastic_api + CenterEuclideanDistanceCutoff< OwnerSegment< 2 > >; + + template < typename Type > + MinimalDistanceCutoff< Type >::MinimalDistanceCutoff( + double cutoff_distance ) + : PairwiseInteraction< Type >(), cutoff_distance_( cutoff_distance ) + { + } + + template < typename Type > + MinimalDistanceCutoff< Type >::MinimalDistanceCutoff( + double cutoff_distance, + typename PairwiseInteraction< Type >::SCOPE scope ) + : PairwiseInteraction< Type >( scope ), + cutoff_distance_( cutoff_distance ) + { + } + + template < typename Type > + double + MinimalDistanceCutoff< Type >::neighborhood_searching_distance() const + { + return cutoff_distance_; + } + + template < typename Type > + double MinimalDistanceCutoff< Type >::compute( + const ObjectRef< Type >& object_a, + const ObjectRef< Type >& object_b ) const + { + auto dist = compute_min_distance< Type::dim >( + object_a.object, object_b.object ); + return dist <= cutoff_distance_ ? 1.0 : 0.0; + } + + template class opengeode_stochastic_stochastic_api + MinimalDistanceCutoff< Point< 2 > >; + template class opengeode_stochastic_stochastic_api + MinimalDistanceCutoff< Point< 3 > >; + + template class opengeode_stochastic_stochastic_api + MinimalDistanceCutoff< OwnerSegment< 2 > >; +} // namespace geode \ No newline at end of file diff --git a/src/geode/stochastic/sampling/mcmc/energy_terms/intensity_term.cpp b/src/geode/stochastic/spatial/single_object_features/segment_length_feature.cpp similarity index 79% rename from src/geode/stochastic/sampling/mcmc/energy_terms/intensity_term.cpp rename to src/geode/stochastic/spatial/single_object_features/segment_length_feature.cpp index a207110..53d1002 100644 --- a/src/geode/stochastic/sampling/mcmc/energy_terms/intensity_term.cpp +++ b/src/geode/stochastic/spatial/single_object_features/segment_length_feature.cpp @@ -22,7 +22,7 @@ * */ -#include +#include #include @@ -95,27 +95,23 @@ namespace ( clipped_dx * clipped_dx ) + ( clipped_dy * clipped_dy ) ); } } // namespace + namespace geode { - IntensityTerm::IntensityTerm( std::string_view name, - double lambda, - std::vector< uuid > targeted_set_ids, - double caracteristic_length, - const SpatialDomain< OwnerSegment2D::dim >& domain ) - : SingleObjectTerm< OwnerSegment2D, - std::function< double( const OwnerSegment2D&, - const SpatialDomain< OwnerSegment2D::dim >& ) > >( - name, - lambda, - std::move( targeted_set_ids ), - 1.0 / caracteristic_length, - []( const OwnerSegment2D& segment, - const SpatialDomain< OwnerSegment2D::dim >& spatial_domain ) { - auto seg_extremities = segment.vertices(); - return length_inside_box( seg_extremities[0], - seg_extremities[1], spatial_domain.box() ); - }, - domain ) + SegmentLengthInsideBoxFeature::SegmentLengthInsideBoxFeature( + double characteristic_length ) + : inv_length_( 1.0 / characteristic_length ) + { + } + + double SegmentLengthInsideBoxFeature::evaluate( + const OwnerSegment2D& segment, const SpatialDomain< 2 >& domain ) const { + auto seg_extremities = segment.vertices(); + + const auto length = length_inside_box( + seg_extremities[0], seg_extremities[1], domain.box() ); + + return inv_length_ * length; } } // namespace geode \ No newline at end of file diff --git a/tests/stochastic/CMakeLists.txt b/tests/stochastic/CMakeLists.txt index 0239d46..d835bab 100644 --- a/tests/stochastic/CMakeLists.txt +++ b/tests/stochastic/CMakeLists.txt @@ -18,169 +18,7 @@ # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE # SOFTWARE. -add_geode_test( - SOURCE "spatial/test-object-set.cpp" - DEPENDENCIES - OpenGeode::basic - OpenGeode::geometry - ${PROJECT_NAME}::stochastic -) - -add_geode_test( - SOURCE "spatial/test-object-sets.cpp" - DEPENDENCIES - OpenGeode::basic - OpenGeode::geometry - ${PROJECT_NAME}::stochastic -) - -add_geode_test( - SOURCE "spatial/test-pairwise-interactions.cpp" - DEPENDENCIES - OpenGeode::basic - OpenGeode::geometry - ${PROJECT_NAME}::stochastic -) - -add_geode_test( - SOURCE "spatial/test-spatial-domain.cpp" - DEPENDENCIES - OpenGeode::basic - OpenGeode::geometry - ${PROJECT_NAME}::stochastic -) - -add_geode_test( - SOURCE "sampling/direct/test-ball-sampler.cpp" - DEPENDENCIES - OpenGeode::basic - OpenGeode::geometry - ${PROJECT_NAME}::stochastic -) - -add_geode_test( - SOURCE "sampling/direct/test-bounding-box-sampler.cpp" - DEPENDENCIES - OpenGeode::basic - OpenGeode::geometry - ${PROJECT_NAME}::stochastic -) - -add_geode_test( - SOURCE "sampling/direct/test-double-sampler.cpp" - DEPENDENCIES - OpenGeode::basic - ${PROJECT_NAME}::stochastic -) - -add_geode_test( - SOURCE "sampling/direct/test-point-uniform-sampler.cpp" - DEPENDENCIES - OpenGeode::basic - OpenGeode::geometry - ${PROJECT_NAME}::stochastic -) - -add_geode_test( - SOURCE "sampling/direct/test-segment-uniform-sampler.cpp" - DEPENDENCIES - OpenGeode::basic - OpenGeode::geometry - ${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/energy_terms/test-density-term.cpp" - DEPENDENCIES - OpenGeode::basic - OpenGeode::geometry - ${PROJECT_NAME}::stochastic -) - -add_geode_test( - SOURCE "sampling/mcmc/energy_terms/test-intensity-term.cpp" - DEPENDENCIES - OpenGeode::basic - OpenGeode::geometry - ${PROJECT_NAME}::stochastic -) - -add_geode_test( - SOURCE "sampling/mcmc/energy_terms/test-pairwise-term.cpp" - DEPENDENCIES - OpenGeode::basic - OpenGeode::geometry - ${PROJECT_NAME}::stochastic -) - -add_geode_test( - SOURCE "sampling/mcmc/energy_terms/test-gibbs-energy.cpp" - DEPENDENCIES - OpenGeode::basic - OpenGeode::geometry - ${PROJECT_NAME}::stochastic -) - -add_geode_test( - SOURCE "sampling/mcmc/proposal/test-proposal-kernel.cpp" - DEPENDENCIES - OpenGeode::basic - OpenGeode::geometry - ${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-fractures.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( - SOURCE "sampling/test-random-engine.cpp" - DEPENDENCIES - OpenGeode::basic - ${PROJECT_NAME}::stochastic -) +add_subdirectory(applications) +add_subdirectory(models) +add_subdirectory(sampling) +add_subdirectory(spatial) \ No newline at end of file diff --git a/tests/stochastic/applications/CMakeLists.txt b/tests/stochastic/applications/CMakeLists.txt new file mode 100644 index 0000000..7e1fa98 --- /dev/null +++ b/tests/stochastic/applications/CMakeLists.txt @@ -0,0 +1,44 @@ +# Copyright (c) 2019 - 2026 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. + +add_geode_test( + SOURCE "test-fractures.cpp" + DEPENDENCIES + OpenGeode::basic + OpenGeode::geometry + ${PROJECT_NAME}::stochastic +) + +add_geode_test( + SOURCE "test-poisson-process.cpp" + DEPENDENCIES + OpenGeode::basic + OpenGeode::geometry + ${PROJECT_NAME}::stochastic + ) + +add_geode_test( + SOURCE "test-strauss-process.cpp" + DEPENDENCIES + OpenGeode::basic + OpenGeode::geometry + ${PROJECT_NAME}::stochastic +) + diff --git a/tests/stochastic/applications/test-fractures.cpp b/tests/stochastic/applications/test-fractures.cpp new file mode 100644 index 0000000..b09e725 --- /dev/null +++ b/tests/stochastic/applications/test-fractures.cpp @@ -0,0 +1,208 @@ +/* + * Copyright (c) 2019 - 2026 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 + +namespace +{ + 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@" ); + + // NOLINTBEGIN(*-magic-numbers) + geode::FractureNetworkDescription fnet_desc; + fnet_desc.fnet_name = "One_Set_FNet"; + constexpr double DOMAIN_BUFFER{ 10 }; + fnet_desc.domain = { geode::Point2D{ { 0, 0 } }, + geode::Point2D{ { 100.0, 100.0 } }, DOMAIN_BUFFER }; + + // --- Object set + auto& fset = fnet_desc.add_fracture_set( "fset_A" ); + // fractures sampler + // desc.name = "uniform_closed"; + // desc.distribution_type = + // geode::UniformClosed< double + // >::distribution_type_static(); + // desc.min_value = 2.; + // desc.max_value = 5.; + + fset.sampler.length.distribution_type = + geode::UniformClosed< double >::distribution_type_static(); + fset.sampler.length.min_value = 1; + fset.sampler.length.max_value = 10.; + + fset.sampler.azimuth.distribution_type = + geode::UniformClosed< double >::distribution_type_static(); + fset.sampler.azimuth.min_value = 1; + fset.sampler.azimuth.max_value = 10.; + + fset.p20 = 0.05; + fset.p21 = 200; + fset.minimal_spacing = 1.; + + // observed fractures + fset.observed_fractures.push_back( { geode::Point2D{ { 0.0, 15. } }, + geode::Point2D{ { 15., 15. } } } ); + fset.observed_fractures.push_back( { geode::Point2D{ { 1.0, 11. } }, + geode::Point2D{ { 11., 20. } } } ); + + geode::Logger::info( fnet_desc.string() ); + + // runner + auto context = build_fractures_simulation_context( fnet_desc ); + geode::FractureSimulationRunner runner{ std::move( context ) }; + // run simulation + geode::SimulationConfigurator sim_config; + sim_config.realizations = 2000; + sim_config.metropolis_hasting_steps = 100; + sim_config.burn_in_steps = 1000; + + geode::SimulationPrinterConfigurator printer_config; + printer_config.output_folder = absl::StrCat( + printer_config.output_folder, "/sim_one_fracture_set_test" ); + sim_config.printer = printer_config; + + auto statistic_tracker = runner.run( engine, sim_config ); + + // NOLINTEND(*-magic-numbers) + + // const auto targeted_statistics_descriptors = + // build_fractures_targeted_stat( fnet_desc ); + // geode::TargetStatistics target_stats{ runner.model(), + // targeted_statistics_descriptors }; + // geode::statistics::validate( statistic_tracker, target_stats + // ); + // + const auto& fset_state = runner.state_realization().get_set( + runner.state_realization().get_set_uuid( + fnet_desc.fracture_sets[0].fset_name ) ); + + geode::OpenGeodeStochasticStochasticException::test( + fset_state.nb_fixed_objects() == 2, + "nd fixed object = ", fset_state.nb_fixed_objects() ); + 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 } } ); + // // todo change + // geode::SpatialDomain domain( box, 0. ); + // + // // --- Object set + // geode::FractureSetDescription setA; + // setA.name = "A"; + // + // // length + // setA.length.distribution_type = + // geode::TruncatedPowerLaw::distribution_type_static(); + // setA.length.alpha = 2.; + // 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.05; + // setA.minimal_spacing = 1.; + // + // // --- Object set + // geode::FractureSetDescription setB; + // setB.name = "B"; + // + // // length + // setB.length.distribution_type = + // geode::TruncatedLogNormal::distribution_type_static(); + // setB.length.min_value = 1; + // setB.length.max_value = 50.; + // setB.length.mean = 1; + // setB.length.standard_deviation = 1.; + // // azimuth + // setB.azimuth.distribution_type = + // geode::VonMises::distribution_type_static(); + // setB.azimuth.mean = 60.; + // setB.azimuth.kappa = 1.; + // + // // positionning + // setB.p20 = 0.05; + // setB.minimal_spacing = 2.; + // + // geode::FractureSimulationRunner runner( domain ); + // runner.add_fracture_set_descriptor( setA ); + // runner.add_fracture_set_descriptor( setB ); + // runner.add_x_node_monitoring( 0.3 ); + // + // 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 = 300; + // sim_config.metropolis_hasting_steps = 1000; + // sim_config.burn_in_steps = 1000; + // sim_config.printer = printer_config; + // + // auto statistic_monitoring = runner.run( engine, sim_config ); + // runner.check_statistics( statistic_monitoring ); + // + // geode::Logger::info( "--> SUCCESS!" ); + // } +} // namespace + +int main() +{ + try + { + geode::OpenGeodeStochasticStochasticLibrary::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/applications/test-poisson-process.cpp b/tests/stochastic/applications/test-poisson-process.cpp new file mode 100644 index 0000000..a04cb7b --- /dev/null +++ b/tests/stochastic/applications/test-poisson-process.cpp @@ -0,0 +1,163 @@ +/* + * Copyright (c) 2019 - 2026 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 + +namespace +{ + + void test_single_type_poisson() + { + geode::Logger::info( "TEST - MH SINGLE TYPE POISSON" ); + + geode::RandomEngine engine; + engine.set_seed( "@mh-test-single-POISSON@" ); + + // NOLINTBEGIN(*-magic-numbers) + 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() } ) + { + geode::PoissonProcessDescription< geode::Point2D > poisson; + poisson.domain = { geode::Point2D{ { 0, 0 } }, + geode::Point2D{ { 10, 10 } }, 0. }; + auto& set_config = poisson.add_set( "set_A" ); + set_config.lambda = 0.3; + set_config.expected_nb_objects = 30; + set_config.birth_ratio = birth_ratio[config]; + set_config.death_ratio = 1.0; + set_config.change_ratio = change_ratio[config]; + + auto simulation_context = build_poisson_process( poisson ); + geode::SimulationRunner< geode::Point2D > runner{ std::move( + simulation_context ) }; + + // run simulation + + geode::SimulationConfigurator sim_config; + sim_config.realizations = 2000; + sim_config.metropolis_hasting_steps = 100; + sim_config.burn_in_steps = 1000; + + geode::SimulationPrinterConfigurator printer_config; + printer_config.output_folder = + absl::StrCat( printer_config.output_folder, + "/sim_point_poisson_test_", config ); + sim_config.printer = printer_config; + + auto statistic_tracker = runner.run( engine, sim_config ); + + const auto targeted_statistics_descriptors = + build_poisson_targeted_stat( poisson ); + geode::TargetStatistics target_stats{ runner.model(), + targeted_statistics_descriptors }; + + geode::statistics::validate( statistic_tracker, target_stats ); + } + // NOLINTEND(*-magic-numbers) + + geode::Logger::info( "--> SUCCESS!" ); + } + + void test_multitype_poisson() + { + geode::Logger::info( "TEST - MH MULTITYPE POISSON" ); + + geode::RandomEngine engine; + engine.set_seed( "@mh-test-POISSON-multi@" ); + + // NOLINTBEGIN(*-magic-numbers) + geode::PoissonProcessDescription< geode::Point2D > poisson; + poisson.domain = { geode::Point2D{ { 0, 0 } }, + geode::Point2D{ { 10, 10 } }, 0. }; + auto& set_config_01 = poisson.add_set( "set01" ); + set_config_01.lambda = 0.1; + set_config_01.expected_nb_objects = 10; + set_config_01.birth_ratio = 2.0; + set_config_01.death_ratio = 3.0; + set_config_01.change_ratio = 1.0; + + auto& set_config_02 = poisson.add_set( "set02" ); + set_config_02.lambda = 0.4; + set_config_01.expected_nb_objects = 40; + set_config_02.birth_ratio = 3.0; + set_config_02.death_ratio = 0.5; + set_config_02.change_ratio = 1.0; + + auto& set_config_03 = poisson.add_set( "set03" ); + set_config_03.lambda = 0.3; + set_config_01.expected_nb_objects = 30; + set_config_03.birth_ratio = 4.0; + set_config_03.death_ratio = 1.0; + set_config_03.change_ratio = 1.0; + + auto context = build_poisson_process( poisson ); + geode::SimulationRunner< geode::Point2D > runner{ std::move( + context ) }; + + // run simulation + geode::SimulationConfigurator sim_config; + sim_config.realizations = 2000; + sim_config.metropolis_hasting_steps = 100; + sim_config.burn_in_steps = 1000; + + geode::SimulationPrinterConfigurator printer_config; + printer_config.output_folder = absl::StrCat( + printer_config.output_folder, "/sim_point_multitype_poisson_test" ); + sim_config.printer = printer_config; + + auto statistic_tracker = runner.run( engine, sim_config ); + + // NOLINTEND(*-magic-numbers) + + const auto targeted_statistics_descriptors = + build_poisson_targeted_stat( poisson ); + geode::TargetStatistics target_stats{ runner.model(), + targeted_statistics_descriptors }; + geode::statistics::validate( statistic_tracker, target_stats ); + + geode::Logger::info( "--> SUCCESS!" ); + } +} // namespace + +int main() +{ + try + { + geode::OpenGeodeStochasticStochasticLibrary::initialize(); + geode::Logger::set_level( geode::Logger::LEVEL::debug ); + test_single_type_poisson(); + test_multitype_poisson(); + return 0; + } + catch( ... ) + { + return geode::geode_lippincott(); + } +} \ No newline at end of file diff --git a/tests/stochastic/applications/test-strauss-process.cpp b/tests/stochastic/applications/test-strauss-process.cpp new file mode 100644 index 0000000..0ae9bb3 --- /dev/null +++ b/tests/stochastic/applications/test-strauss-process.cpp @@ -0,0 +1,203 @@ +/* + * Copyright (c) 2019 - 2026 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 + +namespace +{ + using PoissonDensityDescription = geode::SingleObjectTermConfig; + using PairwiseInteractionDescription = geode::PairwiseTermConfig; + + void test_single_type_strauss() + { + geode::Logger::info( + "TEST - MH SINGLE TYPE STRAUSS (with intra-set interactions)" ); + + geode::RandomEngine engine; + engine.set_seed( "@mh-test-single-STRAUSS@" ); + + // NOLINTBEGIN(*-magic-numbers) + 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.7, 9.8, 18.7, 50.3 }; + for( const auto config : geode::Range{ gamma_values.size() } ) + { + geode::StraussProcessDescription< geode::Point2D > strauss; + + strauss.domain = { geode::Point2D{ { 0.0, 0.0 } }, + geode::Point2D{ { 10.0, 10.0 } }, 1. }; + + auto& set_config = strauss.add_set( "set_A" ); + set_config.lambda = 0.5; + set_config.expected_nb_objects = nb_points[config]; + + auto& interaction_config = + strauss.add_interaction( "interaction_A" ); + interaction_config.set_names = { "set_A" }; + interaction_config.gamma = gamma_values[config]; + interaction_config.distance = 1.0; + interaction_config.expected_nb_interactions = + nb_interactions[config]; + + auto simulation_context = geode::build_strauss_process( strauss ); + geode::SimulationRunner< geode::Point2D > runner{ std::move( + simulation_context ) }; + + // run simulation + geode::SimulationConfigurator sim_config; + sim_config.realizations = 2000; + sim_config.metropolis_hasting_steps = 100; + sim_config.burn_in_steps = 1000; + + geode::SimulationPrinterConfigurator printer_config; + printer_config.output_folder = + absl::StrCat( printer_config.output_folder, + "/sim_point_strauss_test_", config ); + sim_config.printer = printer_config; + + auto statistic_tracker = runner.run( engine, sim_config ); + + const auto targeted_statistics_descriptors = + geode::build_strauss_targeted_stat( strauss ); + + geode::TargetStatistics target_stats{ runner.model(), + targeted_statistics_descriptors }; + + geode::statistics::validate( statistic_tracker, target_stats ); + } + // NOLINTEND(*-magic-numbers) + geode::Logger::info( "--> SUCCESS!" ); + } + + void test_multitype_strauss() + { + geode::Logger::info( + "TEST - MH MULTITYPE STRAUSS (with inter-set interactions ) " ); + + geode::RandomEngine engine; + engine.set_seed( "@mh-test-multi-STRAUSS@" ); + + // NOLINTBEGIN(*-magic-numbers) + std::array< double, 3 > gamma_values{ 0, 0.5, 1.0 }; + 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() } ) + { + geode::StraussProcessDescription< geode::Point2D > strauss; + strauss.domain = { geode::Point2D{ { 0, 0 } }, + geode::Point2D{ { 10, 10 } }, 2. }; + + auto& set_config_01 = strauss.add_set( "set01" ); + set_config_01.lambda = 0.1; + set_config_01.expected_nb_objects = 10; + set_config_01.birth_ratio = 1.0; + set_config_01.death_ratio = 3.0; + set_config_01.change_ratio = 1.0; + + auto& set_config_02 = strauss.add_set( "set02" ); + set_config_02.lambda = 0.4; + set_config_01.expected_nb_objects = 40; + set_config_02.birth_ratio = 3.0; + set_config_02.death_ratio = 0.5; + set_config_02.change_ratio = 1.0; + + auto& set_config_03 = strauss.add_set( "set03" ); + set_config_03.lambda = 0.3; + set_config_01.expected_nb_objects = 30; + set_config_03.birth_ratio = 4.0; + set_config_03.death_ratio = 1.0; + set_config_03.change_ratio = 1.0; + + auto& interaction_config_01 = + strauss.add_interaction( "interaction_01" ); + interaction_config_01.set_names = { "set01", "set02", "set03" }; + interaction_config_01.gamma = gamma_values[config]; + interaction_config_01.distance = 1.0; + interaction_config_01.expected_nb_interactions = + nb_interactions01[config]; + + auto& interaction_config_02 = + strauss.add_interaction( "interaction_02" ); + interaction_config_02.set_names = { "set02" }; + interaction_config_02.gamma = 1.; + interaction_config_02.distance = 2.0; + interaction_config_02.expected_nb_interactions = + nb_interactions02[config]; + + // --- Pairwise interactions + // 1. Intra-type (repulsion within same set) + + // run simulation + auto context = build_strauss_process( strauss ); + geode::SimulationRunner< geode::Point2D > runner{ std::move( + context ) }; + + geode::SimulationConfigurator sim_config; + sim_config.realizations = 2000; + sim_config.metropolis_hasting_steps = 100; + sim_config.burn_in_steps = 1000; + + geode::SimulationPrinterConfigurator printer_config; + printer_config.output_folder = + absl::StrCat( printer_config.output_folder, + "/sim_point_multitype_strauss_test" ); + sim_config.printer = printer_config; + + auto statistic_tracker = runner.run( engine, sim_config ); + + const auto targeted_statistics_descriptors = + geode::build_strauss_targeted_stat( strauss ); + geode::TargetStatistics target_stats{ runner.model(), + targeted_statistics_descriptors }; + + geode::statistics::validate( statistic_tracker, target_stats ); + } + // NOLINTEND(*-magic-numbers) + + geode::Logger::info( "--> SUCCESS!" ); + } +} // namespace + +int main() +{ + try + { + geode::OpenGeodeStochasticStochasticLibrary::initialize(); + geode::Logger::set_level( geode::Logger::LEVEL::debug ); + test_single_type_strauss(); + // test_multitype_strauss(); + return 0; + } + catch( ... ) + { + return geode::geode_lippincott(); + } +} \ No newline at end of file diff --git a/tests/stochastic/models/CMakeLists.txt b/tests/stochastic/models/CMakeLists.txt new file mode 100644 index 0000000..1d9a142 --- /dev/null +++ b/tests/stochastic/models/CMakeLists.txt @@ -0,0 +1,59 @@ +# Copyright (c) 2019 - 2026 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. + +add_geode_test( + SOURCE "energy_terms/test-density-term.cpp" + DEPENDENCIES + OpenGeode::basic + OpenGeode::geometry + ${PROJECT_NAME}::stochastic +) + +add_geode_test( + SOURCE "energy_terms/test-intensity-term.cpp" + DEPENDENCIES + OpenGeode::basic + OpenGeode::geometry + ${PROJECT_NAME}::stochastic +) + +add_geode_test( + SOURCE "energy_terms/test-pairwise-term.cpp" + DEPENDENCIES + OpenGeode::basic + OpenGeode::geometry + ${PROJECT_NAME}::stochastic +) + +#add_geode_test( +# SOURCE "energy_terms/test-gibbs-energy.cpp" +# DEPENDENCIES +# OpenGeode::basic +# OpenGeode::geometry +# ${PROJECT_NAME}::stochastic +#) + +#add_geode_test( +# SOURCE "energy_term-collection.cpp" +# DEPENDENCIES +# OpenGeode::basic +# OpenGeode::geometry +# ${PROJECT_NAME}::stochastic +#) \ No newline at end of file diff --git a/tests/stochastic/sampling/mcmc/energy_terms/test-density-term.cpp b/tests/stochastic/models/energy_terms/test-density-term.cpp similarity index 69% rename from tests/stochastic/sampling/mcmc/energy_terms/test-density-term.cpp rename to tests/stochastic/models/energy_terms/test-density-term.cpp index db1b648..fedbd49 100644 --- a/tests/stochastic/sampling/mcmc/energy_terms/test-density-term.cpp +++ b/tests/stochastic/models/energy_terms/test-density-term.cpp @@ -20,19 +20,24 @@ * SOFTWARE. * */ -#include +#include +#include + +#include #include #include #include +const std::string set_name{ "segments" }; + geode::uuid init_object_set( geode::ObjectSets< geode::Point2D >& pattern ) { geode::Point2D p1{ { 0., 0. } }; geode::Point2D p2{ { 1., 1. } }; geode::Point2D p_buffer{ { 1.3, 0.1 } }; - auto set_id = pattern.add_set( "default_name" ); + auto set_id = pattern.add_set( set_name ); pattern.add_object( std::move( p1 ), set_id, false ); pattern.add_object( std::move( p2 ), set_id, false ); pattern.add_object( std::move( p_buffer ), set_id, false ); // buffer last @@ -48,13 +53,14 @@ geode::SpatialDomain< 2 > init_domain() return geode::SpatialDomain< 2 >{ box, 0.5 }; } -void run_density_test( double lambda, +void run_density_test( const geode::SingleObjectTermConfig& term_config, const geode::ObjectSets< geode::Point2D >& pattern, - const geode::uuid& set_id, const geode::SpatialDomain< 2 >& domain ) { - geode::DensityTerm< geode::Point2D > term( - "density", lambda, { set_id }, domain ); + auto term = geode::build_energy_term< geode::Point2D >( + term_config, pattern, domain ); + auto lambda = term_config.lambda; + const auto& set_id = pattern.get_set_uuid( set_name ); auto neg_log_lambda = -std::log( lambda ); double expected_add = @@ -66,33 +72,33 @@ void run_density_test( double lambda, double expected_total = ( lambda > 0. ? neg_log_lambda * 2. : std::numeric_limits< double >::infinity() ); - double total = term.total_log( pattern ); + double total = term->total_log( pattern ); geode::OpenGeodeStochasticStochasticException::test( total == expected_total, "[DensityTerm] total_log wrong" ); // --- Delta add inside VOI geode::Point2D p_inside{ { 0.5, 0.5 } }; geode::ObjectRef< geode::Point2D > ref_inside{ p_inside, set_id }; - double delta = term.delta_log_add( pattern, ref_inside ); + double delta = term->delta_log_add( pattern, ref_inside ); geode::OpenGeodeStochasticStochasticException::test( delta == expected_add, "[DensityTerm] delta_log_add inside VOI wrong" ); // --- Delta add in buffer (outside VOI) geode::Point2D p_buffer{ { 1.3, 0.0 } }; geode::ObjectRef< geode::Point2D > ref_buffer{ p_buffer, set_id }; - delta = term.delta_log_add( pattern, ref_buffer ); + delta = term->delta_log_add( pattern, ref_buffer ); geode::OpenGeodeStochasticStochasticException::test( delta == 0., "[DensityTerm] delta_log_add outside VOI wrong" ); // --- Delta remove anchored object geode::ObjectId obj_id{ 0, false, set_id }; - delta = term.delta_log_remove( pattern, obj_id ); + delta = term->delta_log_remove( pattern, obj_id ); geode::OpenGeodeStochasticStochasticException::test( delta == expected_remove, "[DensityTerm] delta_log_remove wrong" ); // --- Delta change anchored → buffer geode::ObjectRef< geode::Point2D > new_buffer{ p_buffer, set_id }; - delta = term.delta_log_change( pattern, obj_id, new_buffer ); + delta = term->delta_log_change( pattern, obj_id, new_buffer ); geode::OpenGeodeStochasticStochasticException::test( delta == expected_remove, "[DensityTerm] delta_log_change anchored→buffer wrong" ); @@ -100,13 +106,13 @@ void run_density_test( double lambda, // --- Delta change anchored → anchored geode::Point2D p_anchored{ { 0.1, 0.1 } }; geode::ObjectRef< geode::Point2D > new_anchored{ p_anchored, set_id }; - delta = term.delta_log_change( pattern, obj_id, new_anchored ); + delta = term->delta_log_change( pattern, obj_id, new_anchored ); geode::OpenGeodeStochasticStochasticException::test( delta == 0., "[DensityTerm] delta_log_change anchored→anchored wrong" ); // --- Delta change buffer → anchored geode::ObjectId buffer_id{ 2, false, set_id }; - delta = term.delta_log_change( pattern, buffer_id, ref_inside ); + delta = term->delta_log_change( pattern, buffer_id, ref_inside ); geode::OpenGeodeStochasticStochasticException::test( delta == expected_add, "[DensityTerm] delta_log_change buffer→anchored wrong" ); } @@ -121,12 +127,22 @@ int main() geode::ObjectSets< geode::Point2D > pattern; auto set_id = init_object_set( pattern ); auto domain = init_domain(); + geode::SingleObjectTermConfig term_config; + geode::ObjectInDomainFeatureConfig object_density; + + term_config.term_name = "density"; + term_config.object_set_names = { set_name }; + term_config.lambda = 0.5; + term_config.object_feature = object_density; // Test different lambda values including near-zero - run_density_test( 0.5, pattern, set_id, domain ); - run_density_test( geode::GLOBAL_EPSILON, pattern, set_id, domain ); - run_density_test( 100.0021165, pattern, set_id, domain ); - run_density_test( 0., pattern, set_id, domain ); // zero lambda + run_density_test( term_config, pattern, domain ); + term_config.lambda = geode::GLOBAL_EPSILON; + run_density_test( term_config, pattern, domain ); + term_config.lambda = 100.0021165; + run_density_test( term_config, pattern, domain ); + term_config.lambda = 0.; + run_density_test( term_config, pattern, domain ); // zero lambda } catch( ... ) { @@ -136,11 +152,20 @@ int main() try { geode::OpenGeodeStochasticStochasticLibrary::initialize(); - geode::uuid set_id; + geode::ObjectSets< geode::Point2D > pattern; + auto set_id = init_object_set( pattern ); auto domain = init_domain(); - geode::DensityTerm< geode::Point2D > term( - "zero", -geode::GLOBAL_EPSILON, { set_id }, domain ); + geode::SingleObjectTermConfig term_config; + + term_config.term_name = "zero"; + term_config.object_set_names = { set_name }; + term_config.lambda = -geode::GLOBAL_EPSILON; + term_config.object_feature = geode::ObjectInDomainFeatureConfig{}; + + auto term = geode::build_energy_term< geode::Point2D >( + term_config, pattern, domain ); + geode::Logger::info( "TEST FAILED" ); return 1; } diff --git a/tests/stochastic/sampling/mcmc/energy_terms/test-gibbs-energy.cpp b/tests/stochastic/models/energy_terms/test-gibbs-energy.cpp similarity index 90% rename from tests/stochastic/sampling/mcmc/energy_terms/test-gibbs-energy.cpp rename to tests/stochastic/models/energy_terms/test-gibbs-energy.cpp index dd95573..9412993 100644 --- a/tests/stochastic/sampling/mcmc/energy_terms/test-gibbs-energy.cpp +++ b/tests/stochastic/models/energy_terms/test-gibbs-energy.cpp @@ -25,10 +25,10 @@ #include -#include -#include -#include -#include +#include +#include +#include +#include #include void test_gibbs_energy() @@ -54,12 +54,14 @@ void test_gibbs_energy() geode_unused( term_id ); // Add pairwise term with trivial interaction: always counts 1 for each pair auto interaction = - std::make_unique< geode::EuclideanCutoffInteraction< geode::Point2D > >( + std::make_unique< geode::MinimalDistanceCutoff< geode::Point2D > >( 1000000 ); const auto pwterm_id = energy_terms.add_energy_term( std::make_unique< geode::PairwiseTerm< geode::Point2D > >( - "interaction", 0.8, std::vector< geode::uuid >{ set_id }, + "interaction", 0.8, + std::vector< std::pair< geode::uuid, geode::uuid > >{ + { set_id, set_id } }, std::move( interaction ), domain ) ); geode_unused( pwterm_id ); diff --git a/tests/stochastic/sampling/mcmc/energy_terms/test-intensity-term.cpp b/tests/stochastic/models/energy_terms/test-intensity-term.cpp similarity index 69% rename from tests/stochastic/sampling/mcmc/energy_terms/test-intensity-term.cpp rename to tests/stochastic/models/energy_terms/test-intensity-term.cpp index 6899c5c..c4cafb9 100644 --- a/tests/stochastic/sampling/mcmc/energy_terms/test-intensity-term.cpp +++ b/tests/stochastic/models/energy_terms/test-intensity-term.cpp @@ -20,12 +20,16 @@ * SOFTWARE. * */ -#include - #include +#include +#include + #include +#include #include +const std::string set_name{ "segments" }; + geode::uuid init_segment_set( geode::ObjectSets< geode::OwnerSegment2D >& pattern ) { @@ -43,7 +47,7 @@ geode::uuid init_segment_set( Segment s_buffer{ geode::Point2D{ { 2.0, 2.0 } }, geode::Point2D{ { 3.0, 3.0 } } }; // clipped = 0.0 - auto set_id = pattern.add_set( "segments" ); + auto set_id = pattern.add_set( set_name ); pattern.add_object( std::move( s1 ), set_id, false ); pattern.add_object( std::move( s2 ), set_id, false ); pattern.add_object( std::move( s_buffer ), set_id, false ); @@ -59,18 +63,16 @@ geode::SpatialDomain< 2 > init_domain() return geode::SpatialDomain< 2 >{ box, 0.5 }; } -void run_intensity_test( double lambda, - double characteristic_length, +void run_intensity_test( const geode::SingleObjectTermConfig& term_config, const geode::ObjectSets< geode::OwnerSegment2D >& pattern, - const geode::uuid& set_id, - const geode::SpatialDomain< 2 >& domain ) + const geode::SpatialDomain< 2 >& domain, + double characteristic_length ) { - geode::IntensityTerm term( - "intensity", lambda, { set_id }, characteristic_length, domain ); - + auto term = build_energy_term( term_config, pattern, domain ); + const auto& set_id = pattern.get_set_uuid( set_name ); const double neg_log_lambda = - ( lambda > 0. ? -std::log( lambda ) - : std::numeric_limits< double >::infinity() ); + ( term_config.lambda > 0. ? -std::log( term_config.lambda ) + : std::numeric_limits< double >::infinity() ); // Total clipped length inside domain: // s1: 1.0 @@ -80,11 +82,11 @@ void run_intensity_test( double lambda, const double scaled_total = total_length / characteristic_length; const double expected_total = - ( lambda > 0. ? neg_log_lambda * scaled_total - : std::numeric_limits< double >::infinity() ); + ( term_config.lambda > 0. ? neg_log_lambda * scaled_total + : std::numeric_limits< double >::infinity() ); // --- Total log - double total = term.total_log( pattern ); + double total = term->total_log( pattern ); geode::OpenGeodeStochasticStochasticException::test( total == expected_total, "[IntensityTerm] total_log wrong" ); @@ -94,10 +96,11 @@ void run_intensity_test( double lambda, geode::ObjectRef< geode::OwnerSegment2D > ref_inside{ s_inside, set_id }; double expected_add = - ( lambda > 0. ? neg_log_lambda * ( 1.0 / characteristic_length ) - : std::numeric_limits< double >::infinity() ); + ( term_config.lambda > 0. + ? neg_log_lambda * ( 1.0 / characteristic_length ) + : std::numeric_limits< double >::infinity() ); - double delta = term.delta_log_add( pattern, ref_inside ); + double delta = term->delta_log_add( pattern, ref_inside ); geode::OpenGeodeStochasticStochasticException::test( delta == expected_add, "[IntensityTerm] delta_log_add inside wrong" ); @@ -106,27 +109,27 @@ void run_intensity_test( double lambda, geode::Point2D{ { 3.0, 3.0 } } }; geode::ObjectRef< geode::OwnerSegment2D > ref_out{ s_outside, set_id }; - delta = term.delta_log_add( pattern, ref_out ); + delta = term->delta_log_add( pattern, ref_out ); geode::OpenGeodeStochasticStochasticException::test( delta == 0.0, "[IntensityTerm] delta_log_add outside wrong" ); // --- Delta remove (first segment) geode::ObjectId obj_id{ 0, false, set_id }; - double expected_remove = ( lambda > 0. ? -expected_add : 0.0 ); + double expected_remove = ( term_config.lambda > 0. ? -expected_add : 0.0 ); - delta = term.delta_log_remove( pattern, obj_id ); + delta = term->delta_log_remove( pattern, obj_id ); geode::OpenGeodeStochasticStochasticException::test( delta == expected_remove, "[IntensityTerm] delta_log_remove wrong" ); // --- Delta change: inside → outside - delta = term.delta_log_change( pattern, obj_id, ref_out ); + delta = term->delta_log_change( pattern, obj_id, ref_out ); geode::OpenGeodeStochasticStochasticException::test( delta == expected_remove, "[IntensityTerm] delta_log_change inside→outside wrong" ); // --- Delta change: outside → inside geode::ObjectId buffer_id{ 2, false, set_id }; - delta = term.delta_log_change( pattern, buffer_id, ref_inside ); + delta = term->delta_log_change( pattern, buffer_id, ref_inside ); geode::OpenGeodeStochasticStochasticException::test( delta == expected_add, "[IntensityTerm] delta_log_change outside→inside wrong" ); @@ -135,7 +138,7 @@ void run_intensity_test( double lambda, geode::Point2D{ { 0.8, 0.0 } } }; // length still 1 geode::ObjectRef< geode::OwnerSegment2D > ref_same{ s_same, set_id }; - delta = term.delta_log_change( pattern, obj_id, ref_same ); + delta = term->delta_log_change( pattern, obj_id, ref_same ); geode::OpenGeodeStochasticStochasticException::test( delta == 0.0, "[IntensityTerm] delta_log_change inside→inside wrong" ); } @@ -153,11 +156,21 @@ int main() const double L0 = 1.0; - run_intensity_test( 0.5, L0, pattern, set_id, domain ); - run_intensity_test( - geode::GLOBAL_EPSILON, L0, pattern, set_id, domain ); - run_intensity_test( 50.0, L0, pattern, set_id, domain ); - run_intensity_test( 0.0, L0, pattern, set_id, domain ); + geode::SingleObjectTermConfig term_config; + geode::SegmentLengthInsideBoxFeatureConfig object_length{ L0 }; + + term_config.term_name = "intensity"; + term_config.object_set_names = { set_name }; + term_config.lambda = 0.5; + term_config.object_feature = object_length; + + run_intensity_test( term_config, pattern, domain, L0 ); + term_config.lambda = geode::GLOBAL_EPSILON; + run_intensity_test( term_config, pattern, domain, L0 ); + term_config.lambda = 50.0; + run_intensity_test( term_config, pattern, domain, L0 ); + term_config.lambda = 0.0; + run_intensity_test( term_config, pattern, domain, L0 ); } catch( ... ) { diff --git a/tests/stochastic/sampling/mcmc/energy_terms/test-pairwise-term.cpp b/tests/stochastic/models/energy_terms/test-pairwise-term.cpp similarity index 65% rename from tests/stochastic/sampling/mcmc/energy_terms/test-pairwise-term.cpp rename to tests/stochastic/models/energy_terms/test-pairwise-term.cpp index fc48579..971378b 100644 --- a/tests/stochastic/sampling/mcmc/energy_terms/test-pairwise-term.cpp +++ b/tests/stochastic/models/energy_terms/test-pairwise-term.cpp @@ -22,23 +22,28 @@ */ #include -#include -#include +#include +#include +#include + +#include #include #include #include +const std::string object_set_name = "single_set"; + geode::uuid init_object_set( geode::ObjectSets< geode::Point2D >& pattern ) { geode::Point2D p1{ { 0.1, 0.1 } }; // VOI geode::Point2D p2{ { 0.9, 0.9 } }; // VOI geode::Point2D p3{ { 1.3, 1.3 } }; // buffer - auto set_id = pattern.add_set( "default_name" ); + auto set_id = pattern.add_set( object_set_name ); + pattern.add_object( std::move( p3 ), set_id, false ); pattern.add_object( std::move( p1 ), set_id, false ); pattern.add_object( std::move( p2 ), set_id, false ); - pattern.add_object( std::move( p3 ), set_id, false ); return set_id; } @@ -51,146 +56,142 @@ geode::SpatialDomain< 2 > init_domain() return geode::SpatialDomain< 2 >{ box, 0.5 }; } -void test_pairwise_term( double gamma, +void test_pairwise_term( const geode::PairwiseTermConfig& term_config, const geode::ObjectSets< geode::Point2D >& pattern, - const geode::uuid& set_id, const geode::SpatialDomain< 2 >& domain ) { - auto interaction = - std::make_unique< geode::EuclideanCutoffInteraction< geode::Point2D > >( - 1 ); - geode::PairwiseTerm< geode::Point2D > term( - "strauss", gamma, { set_id }, std::move( interaction ), domain ); + auto term = geode::build_energy_term< geode::Point2D >( + term_config, pattern, domain ); + auto set_id = pattern.get_set_uuid( object_set_name ); - double neg_log_gamma = -std::log( gamma ); + double neg_log_gamma = -std::log( term_config.gamma ); // --- total_log - double total = term.total_log( pattern ); + double total = term->total_log( pattern ); // Only VOI-anchored interactions counted in statistic: p1-p2 geode::OpenGeodeStochasticStochasticException::test( - total == term.contribution( 1 ), "[PairwiseTerm] total_log wrong" ); + total == term->contribution( 1 ), "[PairwiseTerm] total_log wrong" ); // --- delta_add VOI → VOI geode::Point2D p4{ { 0.5, 0.5 } }; geode::ObjectRef< geode::Point2D > p4_ref{ p4, set_id }; - double delta = term.delta_log_add( pattern, p4_ref ); + double delta = term->delta_log_add( pattern, p4_ref ); // interacts with p1 and p2 → 2 pairs geode::OpenGeodeStochasticStochasticException::test( - delta == term.contribution( 2 ), + delta == term->contribution( 2 ), "[PairwiseTerm] delta_log_add VOI wrong" ); // --- delta_add buffer → buffer (outside VOI, inside buffer) geode::Point2D p_buffer{ { 1.2, 1.2 } }; geode::ObjectRef< geode::Point2D > buffer_ref{ p_buffer, set_id }; - delta = term.delta_log_add( pattern, buffer_ref ); + delta = term->delta_log_add( pattern, buffer_ref ); // energy counts interactions: buffer interacts with p2,p3 (p3 is in // buffer) → 2 pairs geode::OpenGeodeStochasticStochasticException::test( - delta == term.contribution( 1 ), + delta == term->contribution( 1 ), "[PairwiseTerm] delta_log_add buffer wrong" ); // --- delta_change VOI → VOI - geode::ObjectId obj_id{ 0, false, set_id }; // p1 - delta = term.delta_log_change( pattern, obj_id, p4_ref ); + geode::ObjectId obj_id{ 1, false, set_id }; // p1 + delta = term->delta_log_change( pattern, obj_id, p4_ref ); // p1 replaced by p4: interacts with p2 → add 1 interaction geode::OpenGeodeStochasticStochasticException::test( - delta == term.contribution( 1 ), + delta == term->contribution( 1 ), "[PairwiseTerm] delta_log_change VOI->VOI wrong" ); // --- delta_change buffer → VOI - geode::ObjectId old_buffer{ 2, false, set_id }; // p3 - delta = term.delta_log_change( pattern, old_buffer, p4_ref ); + geode::ObjectId old_buffer{ 0, false, set_id }; // p3 + delta = term->delta_log_change( pattern, old_buffer, p4_ref ); // old buffer interactions removed (-p3 pairs (-1)), new VOI interactions // added // (+p4 pairs (+2)) net change = 1 - double expected_delta = term.contribution( 1 ); // adjust based on count + double expected_delta = term->contribution( 1 ); // adjust based on count geode::OpenGeodeStochasticStochasticException::test( delta == expected_delta, "[PairwiseTerm] delta_log_change buffer->VOI wrong" ); // --- delta_change VOI → buffer - delta = term.delta_log_change( pattern, obj_id, buffer_ref ); + delta = term->delta_log_change( pattern, obj_id, buffer_ref ); // p1 → buffer: p1 old interaction = 0 // p4 → p_buffer interaction with p2 +1 - expected_delta = term.contribution( 1 ); + expected_delta = term->contribution( 1 ); geode::OpenGeodeStochasticStochasticException::test( delta == expected_delta, "[PairwiseTerm] delta_log_change VOI->buffer wrong" ); // --- delta_remove VOI - delta = term.delta_log_remove( pattern, obj_id ); + delta = term->delta_log_remove( pattern, obj_id ); // p1 removed: no interaction with p2 removed geode::OpenGeodeStochasticStochasticException::test( - delta == term.contribution( 0 ), + delta == term->contribution( 0 ), "[PairwiseTerm] delta_log_remove VOI wrong" ); // --- statistic (only anchored objects in VOI) - double stat = term.statistic( pattern ); + double stat = term->statistic( pattern ); // p1,p2 anchored → 1 pair geode::OpenGeodeStochasticStochasticException::test( stat == 1., "[PairwiseTerm] statistic wrong" ); } -void test_pairwise_term_zero_gamma( double gamma, +void test_pairwise_term_zero_gamma( + const geode::PairwiseTermConfig& term_config, const geode::ObjectSets< geode::Point2D >& pattern, - const geode::uuid& set_id, const geode::SpatialDomain< 2 >& domain ) { - auto interaction = - std::make_unique< geode::EuclideanCutoffInteraction< geode::Point2D > >( - 1 ); - geode::PairwiseTerm< geode::Point2D > term( - "strauss", gamma, { set_id }, std::move( interaction ), domain ); + auto term = geode::build_energy_term< geode::Point2D >( + term_config, pattern, domain ); + auto set_id = pattern.get_set_uuid( object_set_name ); // --- total_log - double total = term.total_log( pattern ); + double total = term->total_log( pattern ); geode::OpenGeodeStochasticStochasticException::test( std::isinf( total ), "[PairwiseTerm] total_log with gamma p4_ref{ p4, set_id }; - double delta = term.delta_log_add( pattern, p4_ref ); + double delta = term->delta_log_add( pattern, p4_ref ); geode::OpenGeodeStochasticStochasticException::test( std::isinf( delta ), "[PairwiseTerm] delta_log_add with gamma buffer_ref{ p_buffer, set_id }; - delta = term.delta_log_add( pattern, buffer_ref ); + delta = term->delta_log_add( pattern, buffer_ref ); geode::OpenGeodeStochasticStochasticException::test( std::isinf( delta ), "[PairwiseTerm] delta_log_add buffer with " "gammadelta_log_change( pattern, obj_id, p4_ref ); geode::OpenGeodeStochasticStochasticException::test( std::isinf( delta ), "[PairwiseTerm] delta_log_change VOI->VOI with " "gammadelta_log_change( pattern, old_buffer, p4_ref ); geode::OpenGeodeStochasticStochasticException::test( std::isinf( delta ), "[PairwiseTerm] delta_log_change buffer->VOI with " - "gammadelta_log_change( pattern, obj_id, buffer_ref ); geode::OpenGeodeStochasticStochasticException::test( std::isinf( delta ), "[PairwiseTerm] delta_log_change VOI->buffer with " - "gammadelta_log_remove( pattern, obj_id ); geode::OpenGeodeStochasticStochasticException::test( delta == 0, "[PairwiseTerm] delta_log_remove VOI with " "gammastatistic( pattern ); // p1,p2 anchored → 1 pair geode::OpenGeodeStochasticStochasticException::test( stat == 1., "[PairwiseTerm] statistic wrong" ); @@ -207,21 +208,25 @@ int main() auto set_id = init_object_set( pattern ); auto domain = init_domain(); - test_pairwise_term( 0.5, pattern, set_id, domain ); - test_pairwise_term( geode::GLOBAL_EPSILON, pattern, set_id, domain ); - test_pairwise_term( 2.0, pattern, set_id, domain ); + geode::MinimalDistanceCutoffConfig interaction_cfg{ 1. }; + + geode::PairwiseTermConfig pw_interaction_cfg; + pw_interaction_cfg.term_name = "strauss"; + pw_interaction_cfg.gamma = 0.5; + pw_interaction_cfg.object_set_names_interactions = { { object_set_name, + object_set_name } }; + pw_interaction_cfg.interaction_config = interaction_cfg; + + test_pairwise_term( pw_interaction_cfg, pattern, domain ); - test_pairwise_term_zero_gamma( - 0.9999 * geode::GLOBAL_EPSILON, pattern, set_id, domain ); + pw_interaction_cfg.gamma = geode::GLOBAL_EPSILON; + test_pairwise_term( pw_interaction_cfg, pattern, domain ); + pw_interaction_cfg.gamma = 2.0; + test_pairwise_term( pw_interaction_cfg, pattern, domain ); - // test_normal_positive_pairwise( 0.5, pattern, set_id ); - // test_normal_positive_pairwise( - // geode::GLOBAL_EPSILON, pattern, set_id ); - // test_normal_positive_pairwise( 2.0, pattern, set_id ); + pw_interaction_cfg.gamma = 0.9999 * geode::GLOBAL_EPSILON; + test_pairwise_term_zero_gamma( pw_interaction_cfg, pattern, domain ); - // test_zero_pairwise( 0., pattern, set_id ); - // test_zero_pairwise( - // 0.9999 * geode::GLOBAL_EPSILON, pattern, set_id ); geode::Logger::info( "TEST SUCCESS" ); return 0; } diff --git a/include/geode/stochastic/sampling/direct/direction_sampler.hpp b/tests/stochastic/models/test-energy-term-collection.cpp similarity index 100% rename from include/geode/stochastic/sampling/direct/direction_sampler.hpp rename to tests/stochastic/models/test-energy-term-collection.cpp diff --git a/tests/stochastic/sampling/CMakeLists.txt b/tests/stochastic/sampling/CMakeLists.txt new file mode 100644 index 0000000..ccb4743 --- /dev/null +++ b/tests/stochastic/sampling/CMakeLists.txt @@ -0,0 +1,97 @@ +# Copyright (c) 2019 - 2026 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. + +add_geode_test( + SOURCE "direct/test-ball-sampler.cpp" + DEPENDENCIES + OpenGeode::basic + OpenGeode::geometry + ${PROJECT_NAME}::stochastic +) + +add_geode_test( + SOURCE "direct/test-bounding-box-sampler.cpp" + DEPENDENCIES + OpenGeode::basic + OpenGeode::geometry + ${PROJECT_NAME}::stochastic +) + +add_geode_test( + SOURCE "direct/test-double-sampler.cpp" + DEPENDENCIES + OpenGeode::basic + ${PROJECT_NAME}::stochastic +) + +add_geode_test( + SOURCE "direct/test-point-uniform-sampler.cpp" + DEPENDENCIES + OpenGeode::basic + OpenGeode::geometry + ${PROJECT_NAME}::stochastic +) + +add_geode_test( + SOURCE "direct/test-segment-uniform-sampler.cpp" + DEPENDENCIES + OpenGeode::basic + OpenGeode::geometry + ${PROJECT_NAME}::stochastic +) + +#add_geode_test( +# SOURCE "mcmc/helpers/test-simulation_printer.cpp" +# DEPENDENCIES +# OpenGeode::basic +# OpenGeode::geometry +# ${PROJECT_NAME}::stochastic +#) + +#add_geode_test( +# SOURCE "mcmc/helpers/test-simulation_monitor.cpp" +# DEPENDENCIES +# OpenGeode::basic +# OpenGeode::geometry +# ${PROJECT_NAME}::stochastic +#) + +add_geode_test( + SOURCE "mcmc/proposal/test-proposal-kernel.cpp" + DEPENDENCIES + OpenGeode::basic + OpenGeode::geometry + ${PROJECT_NAME}::stochastic +) + +#add_geode_test( +# SOURCE "mcmc/test-metropolis-hasting-sampler.cpp" +# DEPENDENCIES +# OpenGeode::basic +# OpenGeode::geometry +# ${PROJECT_NAME}::stochastic +#) + +add_geode_test( + SOURCE "test-random-engine.cpp" + DEPENDENCIES + OpenGeode::basic + ${PROJECT_NAME}::stochastic +) diff --git a/tests/stochastic/sampling/mcmc/helpers/test-simulation_monitor.cpp b/tests/stochastic/sampling/mcmc/helpers/test-simulation_monitor.cpp index ed4512c..c643b65 100644 --- a/tests/stochastic/sampling/mcmc/helpers/test-simulation_monitor.cpp +++ b/tests/stochastic/sampling/mcmc/helpers/test-simulation_monitor.cpp @@ -23,7 +23,7 @@ #include #include -#include +#include #include #include @@ -32,10 +32,10 @@ namespace { void test_statistics_monitor_basic() { - geode::Logger::info( "TEST - StatisticsMonitor basic functionality" ); + geode::Logger::info( "TEST - StatisticsTracker basic functionality" ); // --- Create monitor with 3 energy terms - geode::StatisticsMonitor monitor( 3 ); + geode::StatisticsTracker monitor; // --- Add 2 realizations monitor.add_realization( { 1.0, 2.0, 3.0 } ); diff --git a/tests/stochastic/sampling/mcmc/helpers/test-simulation_printer.cpp b/tests/stochastic/sampling/mcmc/helpers/test-simulation_printer.cpp index a69d147..9a060e9 100644 --- a/tests/stochastic/sampling/mcmc/helpers/test-simulation_printer.cpp +++ b/tests/stochastic/sampling/mcmc/helpers/test-simulation_printer.cpp @@ -22,7 +22,7 @@ */ #include -#include +#include #include #include @@ -71,7 +71,8 @@ namespace "[TEST] SimulationPrinter print statistics summary" ); geode::SimulationPrinter printer( config ); - geode::StatisticsMonitor monitor( 2 ); + + geode::StatisticsTracker tracker( 2 ); monitor.add_realization( { 1, 2 } ); printer.print_statistics_summary( monitor, "EnergyTerm1;EnergyTerm2" ); diff --git a/tests/stochastic/sampling/mcmc/proposal/test-proposal-kernel.cpp b/tests/stochastic/sampling/mcmc/proposal/test-proposal-kernel.cpp index 37d25f6..51e4eaf 100644 --- a/tests/stochastic/sampling/mcmc/proposal/test-proposal-kernel.cpp +++ b/tests/stochastic/sampling/mcmc/proposal/test-proposal-kernel.cpp @@ -35,12 +35,11 @@ namespace { geode::uuid init_object_set( geode::ObjectSets< geode::Point2D >& pattern ) { - geode::Point2D p1{ { 0., 0. } }; - geode::Point2D p2{ { 1., 1. } }; - + // NOLINTBEGIN(*-magic-numbers) auto set_id = pattern.add_set( "default_name" ); - pattern.add_object( std::move( p1 ), set_id, false ); - pattern.add_object( std::move( p2 ), set_id, false ); + pattern.add_object( geode::Point2D{ { 0., 0. } }, set_id, false ); + pattern.add_object( geode::Point2D{ { 1., 1. } }, set_id, false ); + // NOLINTEND(*-magic-numbers) return set_id; } @@ -50,15 +49,14 @@ namespace geode::ObjectSets< geode::Point2D > config; auto set_id = init_object_set( config ); - geode::Point2D min_point{ { 0., 0. } }; - geode::Point2D max_point{ { 10., 100. } }; - + // NOLINTBEGIN(*-magic-numbers) geode::BoundingBox2D box; - box.add_point( min_point ); - box.add_point( max_point ); + box.add_point( geode::Point2D{ { 0., 0. } } ); + box.add_point( geode::Point2D{ { 10., 100. } } ); geode::SpatialDomain domain( box, 0. ); - geode::UniformPointSetSampler< 2 > sampler( domain ); + geode::ObjectSamplerConfig< geode::Point2D > sampler_config; + geode::UniformPointSetSampler< 2 > sampler( domain, sampler_config ); // Create classical birth-death-change kernel auto kernel = geode::create_birth_death_change_kernel< geode::Point2D >( @@ -66,10 +64,13 @@ namespace geode::RandomEngine engine; - bool saw_birth = false, saw_death = false, saw_change = false; + bool saw_birth = false; + bool saw_death = false; + bool saw_change = false; - for( const auto i : geode::Range{ 400 } ) + for( const auto count : geode::Range{ 400 } ) { + geode_unused( count ); auto proposal = kernel->propose( config, engine ); const auto& proposed_move = proposal.proposed_move; @@ -91,7 +92,7 @@ namespace "0." ); geode::OpenGeodeStochasticStochasticException::test( - std::abs( + std::fabs( proposed_move.proposal_probabilities .log_backward_prob - ( std::log( 0.8 * 0.5 ) @@ -173,6 +174,7 @@ namespace proposed_move.type == geode::MoveType::Birth || proposed_move.type == geode::MoveType::Invalid, "[test proposal] On empty config, only Birth should be possible." ); + // NOLINTEND(*-magic-numbers) } } // namespace int main() diff --git a/tests/stochastic/sampling/mcmc/test-metropolis-hasting-sampler.cpp b/tests/stochastic/sampling/mcmc/test-metropolis-hasting-sampler.cpp index a307bf1..ef5ae0f 100644 --- a/tests/stochastic/sampling/mcmc/test-metropolis-hasting-sampler.cpp +++ b/tests/stochastic/sampling/mcmc/test-metropolis-hasting-sampler.cpp @@ -21,16 +21,18 @@ * */ #include +#include +#include +#include #include -#include -#include -#include #include #include namespace { void test_acceptance_prob_helper() { + // NOLINTBEGIN(*-magic-numbers) + // log_accept >= 0 → prob = 1 geode::OpenGeodeStochasticStochasticException::test( geode::MetropolisHastings< geode::Point2D >::acceptance_prob_helper( @@ -46,24 +48,26 @@ namespace "[MH test] acceptance_prob_helper wrong for extreme negative." ); // moderate negative → exp(log_accept) - double val = + auto val = geode::MetropolisHastings< geode::Point2D >::acceptance_prob_helper( -1.0 ); geode::OpenGeodeStochasticStochasticException::test( - std::abs( val - std::exp( -1.0 ) ) < 1e-12, + std::fabs( val - std::exp( -1.0 ) ) < geode::GLOBAL_EPSILON, "[MH test] acceptance_prob_helper wrong for -1.0." ); + // NOLINTEND(*-magic-numbers) } - void test_beta_setter( geode::MetropolisHastings< geode::Point2D >& mh ) + void test_beta_setter( geode::MetropolisHastings< geode::Point2D >& mh_eng ) { - mh.set_beta( 0.5 ); + // NOLINTBEGIN(*-magic-numbers) + mh_eng.set_beta( 0.5 ); geode::OpenGeodeStochasticStochasticException::test( - mh.beta() == 0.5, "[MH test] beta not set correctly." ); + mh_eng.beta() == 0.5, "[MH test] beta not set correctly." ); bool exception_thrown = false; try { - mh.set_beta( -1.0 ); + mh_eng.set_beta( -1.0 ); } catch( ... ) { @@ -71,44 +75,46 @@ namespace } geode::OpenGeodeStochasticStochasticException::test( exception_thrown, "[MH test] negative beta did not throw." ); + // NOLINTEND(*-magic-numbers) } - void test_steps( const geode::MetropolisHastings< geode::Point2D >& mh, + void test_steps( const geode::MetropolisHastings< geode::Point2D >& mh_eng, geode::ObjectSets< geode::Point2D >& state ) { + // NOLINTBEGIN(*-magic-numbers) geode::RandomEngine engine; geode::index_t stat_sum{ 0 }; - constexpr geode::index_t N{ 100000 }; + constexpr geode::index_t NUNBER_ITR{ 100000 }; geode::index_t accepted_birth{ 0 }; geode::index_t accepted_death{ 0 }; geode::index_t accepted_change{ 0 }; geode::index_t nb_accepted{ 0 }; - for( const auto count : geode::Range{ N } ) + for( const auto count : geode::Range{ NUNBER_ITR } ) { - auto result = mh.step( state, engine ); + auto result = mh_eng.step( state, engine ); // Invariant: fixed object must remain geode::OpenGeodeStochasticStochasticException::test( - result.decision == geode::MHDecision::Accepted - || result.decision == geode::MHDecision::Rejected, - "[MH test] decision should be Accepted or Rejected." ); + result.decision == geode::MH_DECISION::accepted + || result.decision == geode::MH_DECISION::rejected, + "[MH test] decision should be accepted or rejected." ); // Log each step (optional: comment out if too verbose) // geode::Logger::info( "Step: ", count, // " move_type= ", static_cast< int >( // result.move_type ), " decision= ", result.decision - // == geode::MHDecision::Accepted ? "Accepted" + // == geode::MH_DECISION::accepted ? "accepted" // : - // "Rejected", + // "rejected", // " delta_log_energy = ", result.delta_log_energy, // " log_accept = ", result.log_accept, // " state_size = ", state.size() ); // Keep track of accepted moves by type - if( result.decision == geode::MHDecision::Accepted ) + if( result.decision == geode::MH_DECISION::accepted ) { nb_accepted++; switch( result.move_type ) @@ -135,7 +141,7 @@ namespace " Mean objects = ", static_cast< double >( stat_sum ) / static_cast< double >( count ), - " nb accepted = ", nb_accepted, " Accepted(B/D/C) = ", + " nb accepted = ", nb_accepted, " accepted(B/D/C) = ", static_cast< double >( accepted_birth ) / static_cast< double >( nb_accepted ), " ", @@ -146,6 +152,7 @@ namespace / static_cast< double >( nb_accepted ) ); } } + // NOLINTEND(*-magic-numbers) } } // namespace @@ -155,7 +162,7 @@ int main() try { geode::OpenGeodeStochasticStochasticLibrary::initialize(); - + // NOLINTBEGIN(*-magic-numbers) geode::Point2D min_point{ { 0., 0. } }; geode::Point2D max_point{ { 10., 10. } }; @@ -182,23 +189,24 @@ int main() domain ) ); geode_unused( term_id ); - geode::MetropolisHastings< geode::Point2D > mh( + geode::MetropolisHastings< geode::Point2D > mh_eng( energy_terms, std::move( kernel ) ); state.add_object( geode::Point2D{ { 1., 1. } }, set_id, true ); state.add_object( geode::Point2D{ { 2., 2. } }, set_id, false ); state.add_object( geode::Point2D{ { 3., 3. } }, set_id, true ); - test_steps( mh, state ); + test_steps( mh_eng, state ); // geode::OpenGeodeStochasticStochasticException::test( state.get_set( // set_id ).nb_fixed_objects() == 2 // ); - test_beta_setter( mh ); + test_beta_setter( mh_eng ); test_acceptance_prob_helper(); geode::Logger::info( "MH TEST SUCCESS" ); return 0; + // NOLINTEND(*-magic-numbers) } catch( ... ) { diff --git a/tests/stochastic/sampling/mcmc/test-mh-fractures.cpp b/tests/stochastic/sampling/mcmc/test-mh-fractures.cpp deleted file mode 100644 index 656b3be..0000000 --- a/tests/stochastic/sampling/mcmc/test-mh-fractures.cpp +++ /dev/null @@ -1,186 +0,0 @@ -/* - * Copyright (c) 2019 - 2026 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 - -namespace -{ - 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 } } ); - // todo change - geode::SpatialDomain domain( box, 0. ); - // --- Object set - geode::FractureSetDescription setA; - setA.name = "A"; - setA.observed_fractures.push_back( { geode::Point2D{ { 0.0, 15. } }, - geode::Point2D{ { 15., 15. } } } ); - setA.observed_fractures.push_back( { geode::Point2D{ { 1.0, 11. } }, - geode::Point2D{ { 11., 20. } } } ); - // 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.05; - setA.p21 = 200; - setA.minimal_spacing = 1.; - - geode::FractureSimulationRunner runner( domain ); - runner.add_fracture_set_descriptor( setA ); - - runner.initialize(); - - // geode::OpenGeodeStochasticStochasticException::test( - // runner.state_realization().nb_fixed_objects() == 2 ); - // 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 = 300; - sim_config.metropolis_hasting_steps = 1000; - sim_config.burn_in_steps = 1000; - sim_config.printer = printer_config; - - auto statistic_monitoring = runner.run( engine, sim_config ); - runner.check_statistics( statistic_monitoring ); - - // geode::OpenGeodeStochasticStochasticException::test( - // runner.state_realization().nb_fixed_objects() == 2 ); - 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 } } ); - // todo change - geode::SpatialDomain domain( box, 0. ); - - // --- Object set - geode::FractureSetDescription setA; - setA.name = "A"; - - // length - setA.length.distribution_type = - geode::TruncatedPowerLaw::distribution_type_static(); - setA.length.alpha = 2.; - 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.05; - setA.minimal_spacing = 1.; - - // --- Object set - geode::FractureSetDescription setB; - setB.name = "B"; - - // length - setB.length.distribution_type = - geode::TruncatedLogNormal::distribution_type_static(); - setB.length.min_value = 1; - setB.length.max_value = 50.; - setB.length.mean = 1; - setB.length.standard_deviation = 1.; - // azimuth - setB.azimuth.distribution_type = - geode::VonMises::distribution_type_static(); - setB.azimuth.mean = 60.; - setB.azimuth.kappa = 1.; - - // positionning - setB.p20 = 0.05; - setB.minimal_spacing = 2.; - - geode::FractureSimulationRunner runner( domain ); - runner.add_fracture_set_descriptor( setA ); - runner.add_fracture_set_descriptor( setB ); - runner.add_x_node_monitoring( 0.3 ); - - 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 = 300; - sim_config.metropolis_hasting_steps = 1000; - sim_config.burn_in_steps = 1000; - sim_config.printer = printer_config; - - auto statistic_monitoring = runner.run( engine, sim_config ); - runner.check_statistics( statistic_monitoring ); - - geode::Logger::info( "--> SUCCESS!" ); - } -} // namespace - -int main() -{ - try - { - geode::OpenGeodeStochasticStochasticLibrary::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 deleted file mode 100644 index 944311b..0000000 --- a/tests/stochastic/sampling/mcmc/test-mh-poisson.cpp +++ /dev/null @@ -1,278 +0,0 @@ -/* - * Copyright (c) 2019 - 2026 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 -namespace -{ - struct SetDescription - { - std::string name; - double birth_ratio{ 1.0 }; - double death_ratio{ 1.0 }; - double change_ratio{ 1.0 }; - }; - - struct PoissonDensityDescription - { - std::string name; - double density; - double target_count; - }; - - class PoissonSimulationRunner - : public geode::SimulationRunner< geode::Point2D > - { - public: - PoissonSimulationRunner( const geode::SpatialDomain< 2 >& domain ) - : geode::SimulationRunner< geode::Point2D >( domain ) {}; - - void add_set_descriptor( const SetDescription& descriptor ) - { - set_descriptors_.push_back( descriptor ); - } - - void add_density_descriptor( - const PoissonDensityDescription& descriptor ) - { - density_descriptors_.push_back( descriptor ); - } - - void initialize() override - { - 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 > >( - 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 ); - } - - // Step 2: create energy terms - for( const auto& energy_desc : density_descriptors_ ) - { - 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.name, "_density" ), - energy_desc.density, - std::vector< geode::uuid >{ set_id }, - this->domain_ ) ) ); - - this->ordered_target_statistics_.push_back( - energy_desc.target_count ); - } - - this->mh_sampler_ = - std::make_unique< geode::MetropolisHastings< geode::Point2D > >( - this->energy_terms_collection_, - std::move( proposal_kernel ) ); - } - - void check_statistics( - 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() } ) - { - 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( 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::abs( 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_; - }; - - void test_single_type_poisson() - { - geode::Logger::info( "TEST - MH SINGLE TYPE POISSON" ); - - geode::RandomEngine engine; - engine.set_seed( "@mh-test-single-POISSON@" ); - - geode::BoundingBox2D box; - box.add_point( geode::Point2D{ { 0.0, 0.0 } } ); - box.add_point( geode::Point2D{ { 10.0, 10.0 } } ); - geode::SpatialDomain domain( box, 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() } ) - { - // --- 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.name = "A"; - densityA.density = 0.3; - densityA.target_count = 30.0; - - PoissonSimulationRunner runner( domain ); - runner.add_set_descriptor( setA ); - runner.add_density_descriptor( densityA ); - runner.initialize(); - - // 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; - - auto statistic_monitoring = runner.run( engine, sim_config ); - runner.check_statistics( statistic_monitoring ); - } - - geode::Logger::info( "--> SUCCESS!" ); - } - - void test_multitype_poisson() - { - geode::Logger::info( "TEST - MH MULTITYPE POISSON" ); - - geode::RandomEngine engine; - engine.set_seed( "@mh-test-POISSON-multi@" ); - - geode::BoundingBox2D box; - box.add_point( geode::Point2D{ { 0.0, 0.0 } } ); - box.add_point( geode::Point2D{ { 10.0, 10.0 } } ); - geode::SpatialDomain domain( box, 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( domain ); - 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(); - - // 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 = 1500; - sim_config.metropolis_hasting_steps = 1000; - sim_config.burn_in_steps = 3000; - sim_config.printer = printer_config; - - auto statistic_monitoring = runner.run( engine, sim_config ); - runner.check_statistics( statistic_monitoring ); - - geode::Logger::info( "--> SUCCESS!" ); - } -} // namespace - -int main() -{ - try - { - geode::OpenGeodeStochasticStochasticLibrary::initialize(); - geode::Logger::set_level( geode::Logger::LEVEL::debug ); - test_single_type_poisson(); - test_multitype_poisson(); - 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 deleted file mode 100644 index 14d11e6..0000000 --- a/tests/stochastic/sampling/mcmc/test-mh-strauss.cpp +++ /dev/null @@ -1,352 +0,0 @@ -/* - * Copyright (c) 2019 - 2026 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 - -namespace -{ - struct SetDescription - { - std::string name; - double birth_ratio{ 1.0 }; - double death_ratio{ 1.0 }; - 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; - }; - - class StraussSimulationRunner - : public geode::SimulationRunner< geode::Point2D > - { - public: - StraussSimulationRunner( const geode::SpatialDomain< 2 >& domain ) - : geode::SimulationRunner< geode::Point2D >( domain ) - { - } - - void add_set_descriptor( const SetDescription& descriptor ) - { - set_descriptors_.push_back( descriptor ); - } - - void add_density_descriptor( - const PoissonDensityDescription& descriptor ) - { - density_descriptors_.push_back( descriptor ); - } - - void add_interaction_descriptor( - const PairwiseInteractionDescription& descriptor ) - { - interaction_descriptors_.push_back( descriptor ); - } - - void initialize() override - { - 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_ ) ); - - 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& density_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 ); - } - - // 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::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.names, "_" ), - "_interaction" ), - interaction_desc.strength, set_ids, - std::move( interaction ), this->domain_ ) ) ); - - this->ordered_target_statistics_.push_back( - interaction_desc.target_interaction_count ); - } - - this->mh_sampler_ = - std::make_unique< geode::MetropolisHastings< geode::Point2D > >( - 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] ); - - const auto expected_mean = - this->ordered_target_statistics_[stat_id]; - auto target_vs_mean_error = - std::abs( 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 ); - } - } - - private: - std::vector< SetDescription > set_descriptors_; - std::vector< PoissonDensityDescription > density_descriptors_; - std::vector< PairwiseInteractionDescription > interaction_descriptors_; - }; - - void test_single_type_strauss() - { - geode::Logger::info( - "TEST - MH SINGLE TYPE STRAUSS (with intra-set interactions)" ); - - geode::RandomEngine engine; - 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 } } ); - // 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 }; - - 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]; - - // --- 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]; - - StraussSimulationRunner runner( domain ); - runner.add_set_descriptor( setA ); - runner.add_density_descriptor( densityA ); - runner.add_interaction_descriptor( intraA ); - - runner.initialize(); - - // 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; - - auto statistic_monitoring = runner.run( engine, sim_config ); - runner.check_statistics( statistic_monitoring ); - } - - geode::Logger::info( "--> SUCCESS!" ); - } - - void test_multitype_strauss() - { - geode::Logger::info( - "TEST - MH MULTITYPE STRAUSS (with inter-set interactions)" ); - - geode::RandomEngine engine; - 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 } } ); - // todo change!! - geode::SpatialDomain domain( box, 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{ 8, 20, 85 }; - - 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( domain ); - 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(); - - // 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; - - auto statistic_monitoring = runner.run( engine, sim_config ); - runner.check_statistics( statistic_monitoring ); - } - - geode::Logger::info( "--> SUCCESS!" ); - } -} // namespace - -int main() -{ - try - { - geode::OpenGeodeStochasticStochasticLibrary::initialize(); - geode::Logger::set_level( geode::Logger::LEVEL::debug ); - test_single_type_strauss(); - test_multitype_strauss(); - return 0; - } - catch( ... ) - { - return geode::geode_lippincott(); - } -} \ No newline at end of file diff --git a/tests/stochastic/sampling/test-random-engine.cpp b/tests/stochastic/sampling/test-random-engine.cpp index 49b0f53..d002ee5 100644 --- a/tests/stochastic/sampling/test-random-engine.cpp +++ b/tests/stochastic/sampling/test-random-engine.cpp @@ -28,204 +28,218 @@ #include #include - -const int NUMBER_OF_SAMPLES = 10000; - -void test_reproducibility() +// NOLINTBEGIN(*-magic-numbers) +namespace { - std::vector< std::string > seeds = { "same-seed", "geode-solutions", - "dfgzejhèçsodj", "&%;:!§" }; - for( const auto seed : seeds ) - { // Create first engine geode:: - geode::RandomEngine engine1; - engine1.set_seed( seed ); - - // Create second engine with same seed - geode::RandomEngine engine2; - engine2.set_seed( seed ); - - // Define a uniform distribution - geode::UniformClosed< int > dist; - dist.min_value = 1; - dist.max_value = 100; - - // Sample more values to check sequence reproducibility - for( const auto value : geode::Range{ 100 } ) - { - geode::OpenGeodeStochasticStochasticException::test( - engine1.sample_uniform( dist ) - == engine2.sample_uniform( dist ), - "[REPRODUCIBILITY] - Same seed should produce same output." ); + const int NUMBER_OF_SAMPLES = 100000; + + void test_reproducibility() + { + std::vector< std::string > seeds = { "same-seed", "geode-solutions", + "dfgzejhèçsodj", "&%;:!§" }; + for( const auto& seed : seeds ) + { // Create first engine geode:: + geode::RandomEngine engine1; + engine1.set_seed( seed ); + + // Create second engine with same seed + geode::RandomEngine engine2; + engine2.set_seed( seed ); + + // Define a uniform distribution + geode::UniformClosed< int > dist; + dist.min_value = 1; + dist.max_value = 100; + + // Sample more values to check sequence reproducibility + for( const auto value : geode::Range{ 100 } ) + { + geode_unused( value ); + geode::OpenGeodeStochasticStochasticException::test( + engine1.sample_uniform( dist ) + == engine2.sample_uniform( dist ), + "[REPRODUCIBILITY] - Same seed should produce same " + "output." ); + } + geode::Logger::info( "Test Reproducibility: ", seed, " SUCCESS " ); } - geode::Logger::info( "Test Reproducibility: ", seed, " SUCCESS " ); + geode::Logger::info( "Test Reproducibility: SUCCESS " ); } - geode::Logger::info( "Test Reproducibility: SUCCESS " ); -} - -template < typename T > -double compute_mean( const std::vector< T >& data ) -{ - double sum = std::accumulate( data.begin(), data.end(), 0.0 ); - return sum / data.size(); -} -template < typename T > -double compute_variance( const std::vector< T >& data, double mean ) -{ - double accum = 0.0; - for( auto& x : data ) + template < typename T > + double compute_mean( const std::vector< T >& data ) { - double diff = x - mean; - accum += diff * diff; + double sum = std::accumulate( data.begin(), data.end(), 0.0 ); + return sum / data.size(); } - return accum / ( data.size() - 1 ); -} -template < typename T > -void test_distribution_mean_and_variance( const std::vector< T >& data, - double expected_mean, - double expected_var, - double k = 3.0 ) -{ - const auto n = data.size(); - const double mean = compute_mean( data ); - const double variance = compute_variance( data, mean ); - geode::Logger::info( "Test Distribution ", - ": mean / expected mean = ", mean, "/", expected_mean, - " variance / expected variance = ", variance, "/", expected_var ); - - const double se_mean = std::sqrt( expected_var / n ); - const double se_var = - std::sqrt( 2.0 * expected_var * expected_var / ( n - 1 ) ); - - geode::OpenGeodeStochasticStochasticException::test( - std::abs( mean - expected_mean ) < k * se_mean, - "[Uniform] - Wrong expected mean." ); - geode::OpenGeodeStochasticStochasticException::test( - std::abs( variance - expected_var ) < k * se_var, - "[Uniform] - Wrong expected std." ); -} -template < typename T > -double compute_expected_variance( T min_value, T max_value ) -{ - if constexpr( std::is_integral_v< T > ) + template < typename T > + double compute_variance( const std::vector< T >& data, double mean ) { - const double distance = - static_cast< double >( max_value - min_value + 1 ); - return ( distance * distance - 1.0 ) / 12.0; + double accum = 0.0; + for( auto& x : data ) + { + double diff = x - mean; + accum += diff * diff; + } + return accum / ( data.size() - 1 ); } - else if constexpr( std::is_floating_point_v< T > ) + // NOLINTBEGING(*-magic-numbers) + template < typename T > + void test_distribution_mean_and_variance( const std::vector< T >& data, + double expected_mean, + double expected_var, + double k_coef = 5.0 ) { - const double distance = static_cast< double >( max_value - min_value ); - return ( distance * distance ) / 12.0; + const auto data_size = data.size(); + const double mean = compute_mean( data ); + const double variance = compute_variance( data, mean ); + geode::Logger::info( "Test Distribution ", + ": mean / expected mean = ", mean, "/", expected_mean, + " variance / expected variance = ", variance, "/", expected_var ); + + const double se_mean = std::sqrt( expected_var / data_size ); + const double se_var = + std::sqrt( 2.0 * expected_var * expected_var / ( data_size - 1 ) ); + + geode::OpenGeodeStochasticStochasticException::test( + std::fabs( mean - expected_mean ) < k_coef * se_mean, + "[Uniform] - Wrong expected mean." ); + geode::OpenGeodeStochasticStochasticException::test( + std::fabs( variance - expected_var ) < k_coef * se_var, + "[Uniform] - Wrong expected std." ); } - else + + template < typename T > + double compute_expected_variance( T min_value, T max_value ) { - static_assert( std::is_arithmetic_v< T >, - "Unsupported type to compute theoretical variance." ); + if constexpr( std::is_integral_v< T > ) + { + const auto distance = + static_cast< double >( max_value - min_value + 1 ); + return ( distance * distance - 1.0 ) / 12.0; + } + else if constexpr( std::is_floating_point_v< T > ) + { + const auto distance = + static_cast< double >( max_value - min_value ); + return ( distance * distance ) / 12.0; + } + else + { + static_assert( std::is_arithmetic_v< T >, + "Unsupported type to compute theoretical variance." ); + } } -} - -template < typename Type > -void test_uniform( - const Type min_value, const Type max_value, geode::RandomEngine& engine ) -{ - std::vector< Type > samples; - samples.reserve( NUMBER_OF_SAMPLES ); - for( const auto i : geode::Range{ NUMBER_OF_SAMPLES } ) + template < typename Type > + void test_uniform( const Type min_value, + const Type max_value, + geode::RandomEngine& engine ) { - geode::UniformClosed< Type > dist_closed; - dist_closed.min_value = min_value; - dist_closed.max_value = max_value; + std::vector< Type > samples; + samples.reserve( NUMBER_OF_SAMPLES ); - Type value = engine.sample_uniform( dist_closed ); - samples.emplace_back( value ); - geode::OpenGeodeStochasticStochasticException::test( - value >= min_value && value <= max_value, - "[Uniform] - value out of range." ); - } - - double expected_mean = ( min_value + max_value ) / 2.0; - double expected_var = compute_expected_variance( min_value, max_value ); + for( const auto count : geode::Range{ NUMBER_OF_SAMPLES } ) + { + geode_unused( count ); + geode::UniformClosed< Type > dist_closed; + dist_closed.min_value = min_value; + dist_closed.max_value = max_value; - test_distribution_mean_and_variance( samples, expected_mean, expected_var ); -} + Type value = engine.sample_uniform( dist_closed ); + samples.emplace_back( value ); + geode::OpenGeodeStochasticStochasticException::test( + value >= min_value && value <= max_value, + "[Uniform] - value out of range." ); + } -void test_gaussian( - double mean_value, double std_value, geode::RandomEngine& engine ) -{ - std::vector< double > samples; - samples.reserve( NUMBER_OF_SAMPLES ); - geode::Gaussian spec; - spec.mean = mean_value; - spec.standard_deviation = std_value; + double expected_mean = ( min_value + max_value ) / 2.0; + double expected_var = compute_expected_variance( min_value, max_value ); - for( const auto i : geode::Range{ NUMBER_OF_SAMPLES } ) - { - double value = engine.sample_gaussian( spec ); - samples.emplace_back( value ); + test_distribution_mean_and_variance( + samples, expected_mean, expected_var ); } - test_distribution_mean_and_variance( samples, mean_value, std_value ); -} - -void test_truncated_gaussian( double mean_value, - double std_value, - std::optional< double > min_value, - std::optional< double > max_value, - geode::RandomEngine& engine ) -{ - const auto max = - max_value.value_or( std::numeric_limits< double >::infinity() ); - const auto min = - min_value.value_or( -std::numeric_limits< double >::infinity() ); + void test_gaussian( + double mean_value, double std_value, geode::RandomEngine& engine ) + { + std::vector< double > samples; + samples.reserve( NUMBER_OF_SAMPLES ); + geode::Gaussian spec; + spec.mean = mean_value; + spec.standard_deviation = std_value; - std::vector< double > samples; - samples.reserve( NUMBER_OF_SAMPLES ); + for( const auto count : geode::Range{ NUMBER_OF_SAMPLES } ) + { + geode_unused( count ); + double value = engine.sample_gaussian( spec ); + samples.emplace_back( value ); + } - geode::TruncatedGaussian spec; - spec.mean = mean_value; - spec.standard_deviation = std_value; - spec.min_value = min_value; - spec.max_value = max_value; + test_distribution_mean_and_variance( samples, mean_value, std_value ); + } - for( const auto i : geode::Range{ NUMBER_OF_SAMPLES } ) + void test_truncated_gaussian( double mean_value, + double std_value, + std::optional< double > min_value, + std::optional< double > max_value, + geode::RandomEngine& engine ) { - double value = engine.sample_truncated_gaussian( spec ); - samples.emplace_back( value ); - geode::OpenGeodeStochasticStochasticException::test( - value >= min && value <= max, "[Gaussian] - value out of range." ); - } + const auto max = + max_value.value_or( std::numeric_limits< double >::infinity() ); + const auto min = + min_value.value_or( -std::numeric_limits< double >::infinity() ); - // Implementer un test de Chi2 ou Kolmogorov–Smirnov -} + std::vector< double > samples; + samples.reserve( NUMBER_OF_SAMPLES ); -void test_bernoulli( - double probability_of_success, geode::RandomEngine& engine ) -{ - int success_count = 0; + geode::TruncatedGaussian spec; + spec.mean = mean_value; + spec.standard_deviation = std_value; + spec.min_value = min_value; + spec.max_value = max_value; - for( const auto i : geode::Range{ NUMBER_OF_SAMPLES } ) - { - if( engine.sample_bernoulli( probability_of_success ) ) + for( const auto count : geode::Range{ NUMBER_OF_SAMPLES } ) { - ++success_count; + geode_unused( count ); + double value = engine.sample_truncated_gaussian( spec ); + samples.emplace_back( value ); + geode::OpenGeodeStochasticStochasticException::test( + value >= min && value <= max, + "[Gaussian] - value out of range." ); } + + // Implementer un test de Chi2 ou Kolmogorov–Smirnov } - const double empirical_probability = - static_cast< double >( success_count ) / NUMBER_OF_SAMPLES; + void test_bernoulli( + double probability_of_success, geode::RandomEngine& engine ) + { + int success_count = 0; - geode::Logger::info( "Test Bernoulli ", - ": empirical / expected = ", empirical_probability, " / ", - probability_of_success ); + for( const auto count : geode::Range{ NUMBER_OF_SAMPLES } ) + { + geode_unused( count ); + if( engine.sample_bernoulli( probability_of_success ) ) + { + ++success_count; + } + } - geode::OpenGeodeStochasticStochasticException::test( - abs( empirical_probability - probability_of_success ) < 0.05, - "[Bernoulli] - Empirical probability out of tolerance." ); -} + const double empirical_probability = + static_cast< double >( success_count ) / NUMBER_OF_SAMPLES; + geode::Logger::info( "Test Bernoulli ", + ": empirical / expected = ", empirical_probability, " / ", + probability_of_success ); + + // NOLINTBEGING(*-magic-numbers) + geode::OpenGeodeStochasticStochasticException::test( + abs( empirical_probability - probability_of_success ) < 0.05, + "[Bernoulli] - Empirical probability out of tolerance." ); + } +} // namespace int main() { try @@ -234,6 +248,8 @@ int main() geode::RandomEngine random_engine; test_reproducibility(); + random_engine.set_seed( "@-test-radom-engine@" ); + geode::Logger::info( "TEST UNIFORM SAMPLING" ); test_uniform< geode::index_t >( 1, 15, random_engine ); test_uniform< geode::local_index_t >( 1, 6, random_engine ); @@ -262,7 +278,6 @@ int main() test_bernoulli( 0.06, random_engine ); test_bernoulli( 0.96, random_engine ); geode::Logger::info( "TEST BERNOUILLI SAMPLING - SUCCESS" ); - return 0; } catch( ... ) @@ -270,3 +285,4 @@ int main() return geode::geode_lippincott(); } } +// NOLINTEND(*-magic-numbers) diff --git a/tests/stochastic/spatial/CMakeLists.txt b/tests/stochastic/spatial/CMakeLists.txt new file mode 100644 index 0000000..447e969 --- /dev/null +++ b/tests/stochastic/spatial/CMakeLists.txt @@ -0,0 +1,51 @@ +# Copyright (c) 2019 - 2026 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. + +add_geode_test( + SOURCE "test-object-set.cpp" + DEPENDENCIES + OpenGeode::basic + OpenGeode::geometry + ${PROJECT_NAME}::stochastic +) + +add_geode_test( + SOURCE "test-object-sets.cpp" + DEPENDENCIES + OpenGeode::basic + OpenGeode::geometry + ${PROJECT_NAME}::stochastic +) + +add_geode_test( + SOURCE "test-pairwise-interactions.cpp" + DEPENDENCIES + OpenGeode::basic + OpenGeode::geometry + ${PROJECT_NAME}::stochastic +) + +add_geode_test( + SOURCE "test-spatial-domain.cpp" + DEPENDENCIES + OpenGeode::basic + OpenGeode::geometry + ${PROJECT_NAME}::stochastic +) diff --git a/tests/stochastic/spatial/test-object-sets.cpp b/tests/stochastic/spatial/test-object-sets.cpp index 848bea8..79f931f 100644 --- a/tests/stochastic/spatial/test-object-sets.cpp +++ b/tests/stochastic/spatial/test-object-sets.cpp @@ -29,8 +29,8 @@ namespace void test_add_sets_and_objects() { ObjectSets< geode::Point2D > sets; - const auto set_id1 = sets.add_set( "default_name" ); - const auto set_id2 = sets.add_set( "default_name" ); + const auto set_id1 = sets.add_set( "default_name1" ); + const auto set_id2 = sets.add_set( "default_name2" ); geode::OpenGeodeStochasticStochasticException::test( sets.nb_sets() == 2, "[TestObjectSets] - Expected 2 sets" ); @@ -122,8 +122,9 @@ namespace sets.add_object( geode::Point2D{ { 5.0, 0.0 } }, set_id, false ); std::vector< uuid > targeted_set_ids{ set_id }; - const auto near_neighbors = - sets.neighbors( obj0, targeted_set_ids, 2.0 ); + + const auto near_neighbors = sets.neighbors( + sets.get_object( obj0 ), targeted_set_ids, 2.0, obj0 ); geode::OpenGeodeStochasticStochasticException::test( near_neighbors.size() == 1, @@ -145,7 +146,8 @@ namespace const geode::Point2D query{ { 0.5, 0.0 } }; std::vector< uuid > targeted_set_ids{ set_id }; - const auto nearby = sets.neighbors( query, targeted_set_ids, 1.0 ); + const auto nearby = + sets.neighbors( query, targeted_set_ids, 1.0, std::nullopt ); geode::OpenGeodeStochasticStochasticException::test( nearby.size() == 2, "[TestObjectSets] - Expected 2 neighbors near query object" ); diff --git a/tests/stochastic/spatial/test-pairwise-interactions.cpp b/tests/stochastic/spatial/test-pairwise-interactions.cpp index 8208e2d..0ede96b 100644 --- a/tests/stochastic/spatial/test-pairwise-interactions.cpp +++ b/tests/stochastic/spatial/test-pairwise-interactions.cpp @@ -24,12 +24,12 @@ #include #include -#include +#include namespace { void test_interaction() { - geode::EuclideanCutoffInteraction< geode::Point2D > interaction( 2.0 ); + geode::MinimalDistanceCutoff< geode::Point2D > interaction( 2.0 ); geode::Point2D p1{ { 0., 0. } }; geode::Point2D p2{ { 1., 1. } }; // distance = 1.44 < 2.0 @@ -39,11 +39,11 @@ namespace interaction.evaluate( geode::ObjectRef< geode::Point2D >{ p1, id }, geode::ObjectRef< geode::Point2D >{ p2, id } ); geode::OpenGeodeStochasticStochasticException::test( result == 1.0, - "[EuclideanCutoffInteraction] Failed for inside cutoff case." ); + "[MinimalDistanceCutoff] Failed for inside cutoff case." ); } void test_no_interaction() { - geode::EuclideanCutoffInteraction< geode::Point2D > interaction( 2.0 ); + geode::MinimalDistanceCutoff< geode::Point2D > interaction( 2.0 ); geode::Point2D p1{ { 0., 0. } }; geode::Point2D p2{ { 3., 0. } }; // distance = 3.0 > 2.0 @@ -53,11 +53,11 @@ namespace interaction.evaluate( geode::ObjectRef< geode::Point2D >{ p1, id }, geode::ObjectRef< geode::Point2D >{ p2, id } ); geode::OpenGeodeStochasticStochasticException::test( result == 0.0, - "[EuclideanCutoffInteraction] Failed for outside cutoff case." ); + "[MinimalDistanceCutoff] Failed for outside cutoff case." ); } void test_limit_interaction() { - geode::EuclideanCutoffInteraction< geode::Point2D > interaction( 2.0 ); + geode::MinimalDistanceCutoff< geode::Point2D > interaction( 2.0 ); geode::Point2D p1{ { 0., 0. } }; geode::Point2D p2{ { 2., 0. } }; // distance = 2.0 == cutoff @@ -67,7 +67,7 @@ namespace interaction.evaluate( geode::ObjectRef< geode::Point2D >{ p1, id }, geode::ObjectRef< geode::Point2D >{ p2, id } ); geode::OpenGeodeStochasticStochasticException::test( result == 1.0, - "[EuclideanCutoffInteraction] Failed for exact cutoff case." ); + "[MinimalDistanceCutoff] Failed for exact cutoff case." ); } } // namespace