diff --git a/include/geode/stochastic/applications/fractures.hpp b/include/geode/stochastic/applications/fractures.hpp index b10a8e2..930cef6 100644 --- a/include/geode/stochastic/applications/fractures.hpp +++ b/include/geode/stochastic/applications/fractures.hpp @@ -94,21 +94,6 @@ namespace geode } }; - // 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; @@ -125,6 +110,23 @@ namespace geode return fracture_set; } + void add_x_node_monitoring( double beta ) + { + OpenGeodeStochasticStochasticException::check_exception( + beta <= 1.0 && beta >= 0., nullptr, + OpenGeodeException::TYPE::data, + "[FractureSimulationRunner] x node should be inhibitated, " + "please provise a value in [0., 1.]." ); + beta_x_node = beta; + } + + [[nodiscard]] std::string x_node_interaction_name() const + { + return absl::StrCat( fnet_name, "_x_node" ); + } + double beta_x_node{ 1. }; + std::optional< double > expected_x_node; + [[nodiscard]] std::string string() const { auto message = @@ -134,6 +136,8 @@ namespace geode { absl::StrAppend( &message, "\n\t --> ", fset_desc.string() ); } + absl::StrAppend( &message, "\n\t --> ", x_node_interaction_name(), + ": ", beta_x_node ); return message; } // std::vector< StraussInteractionDescription< ObjectType > > diff --git a/include/geode/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp b/include/geode/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp deleted file mode 100644 index 160aba0..0000000 --- a/include/geode/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp +++ /dev/null @@ -1,274 +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 -#include -#include -#include -#include -#include -#include -#include - -#include - -namespace geode -{ - struct FractureSetDescription - { - std::string name; - - DoubleSampler::DistributionDescription length; - DoubleSampler::DistributionDescription azimuth; - - // positionning - double p20; - std::vector< std::array< geode::Point2D, 2 > > observed_fractures; - double p21{ 1. }; - double minimal_spacing{ 0. }; - - // mh dynamique - double birth_ratio{ 1.0 }; - double death_ratio{ 1.0 }; - double change_ratio{ 1.0 }; - - std::string string() const - { - auto message = absl::StrCat( "FractureSetDescription: ", 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: ", length.string() ); - absl::StrAppend( - &message, "\n\t --> azimuth distribution: ", azimuth.string() ); - absl::StrAppend( &message, "\n\t --> targeted p20: ", p20 ); - absl::StrAppend( - &message, "\n\t --> minimal spacing: ", minimal_spacing ); - absl::StrAppend( &message, - "\n\t --> MH move ratio - birth/death/change (", birth_ratio, - " / ", death_ratio, " / ", change_ratio, ")" ); - return message; - } - }; - - class FractureSimulationRunner : public SimulationRunner< OwnerSegment2D > - { - public: - FractureSimulationRunner( const SpatialDomain< 2 >& domain ) - : SimulationRunner< OwnerSegment2D >( domain ) - { - } - - 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; - } - void add_fracture_set_descriptor( - const FractureSetDescription& descriptor ) - { - set_descriptors_.push_back( descriptor ); - } - - void initialize() override - { - auto proposal_kernel = - std::make_unique< ProposalKernel< OwnerSegment2D > >(); - - // Mapping set names -> UUID - // std::unordered_map< std::string, uuid > - // set_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 ); - for( const auto& fixed_object : set_desc.observed_fractures ) - { - this->object_sets_.add_object( - geode::OwnerSegment2D{ - fixed_object[0], fixed_object[1] }, - set_id, true ); - } - set_name_to_uuid_[set_desc.name] = set_id; - - auto length_distribution = - DoubleSampler::create_distribution( set_desc.length ); - auto azimuth_distribution = - DoubleSampler::create_rad_angle_distribution_from_degree( - set_desc.azimuth ); - this->set_samplers_.push_back( - std::make_unique< UniformSegmentSetSampler >( - domain_, length_distribution, azimuth_distribution ) ); - - add_birth_death_change_moves( this->set_samplers_.back(), - *proposal_kernel, set_id, set_desc.birth_ratio, - set_desc.death_ratio, set_desc.change_ratio ); - } - - // Step 2: create density energy terms - for( const auto& set_desc : set_descriptors_ ) - { - set_density_term( set_desc ); - set_intensity_term( set_desc ); - set_minimal_spacing_term( set_desc ); - } - set_x_intersection_term(); - - this->mh_sampler_ = - std::make_unique< MetropolisHastings< OwnerSegment2D > >( - this->energy_terms_collection_, - std::move( proposal_kernel ) ); - } - - void check_statistics( - const StatisticsMonitor& statistic_monitoring ) const - { - const auto& computed_means = statistic_monitoring.means(); - - for( const auto stat_id : - Range{ this->energy_terms_collection_.size() } ) - { - const auto& term = energy_terms_collection_.get( - ordered_energy_terms_[stat_id] ); - Logger::info( "[MH test] Statistic value ", - computed_means[stat_id], " for energy term: ", - term.name().value_or( term.id().string() ) ); - } - } - - std::string string() const - { - auto message = - absl::StrCat( "Fracture Simulation Runner description" ); - for( const auto& desc : set_descriptors_ ) - { - absl::StrAppend( &message, "\n\t ", desc.string() ); - } - if( std::fabs( beta_x_node_ - 1. ) > GLOBAL_EPSILON ) - { - absl::StrAppend( &message, - "\n\t --> x node monitioring (beta inhibition value): ", - beta_x_node_ ); - } - else - { - absl::StrAppend( - &message, "\n\t --> x node monitioring : no inhibition." ); - } - return message; - } - - private: - void set_density_term( const FractureSetDescription& set_desc ) - { - const auto set_id = set_name_to_uuid_.at( set_desc.name ); - this->ordered_energy_terms_.push_back( - this->energy_terms_collection_.add_energy_term( - std::make_unique< DensityTerm< OwnerSegment2D > >( - absl::StrCat( set_desc.name, "_density" ), set_desc.p20, - std::vector< uuid >{ set_id }, this->domain_ ) ) ); - } - - void set_intensity_term( const FractureSetDescription& set_desc ) - { - if( std::fabs( set_desc.p21 - 1. ) < GLOBAL_EPSILON ) - { - return; - } - const auto set_id = set_name_to_uuid_.at( set_desc.name ); - this->ordered_energy_terms_.push_back( - this->energy_terms_collection_.add_energy_term( - std::make_unique< IntensityTerm >( - absl::StrCat( set_desc.name, "_intensity" ), - set_desc.p21, std::vector< uuid >{ set_id }, - 0.5 * this->domain_.smallest_length(), - this->domain_ ) ) ); - } - - void set_minimal_spacing_term( const FractureSetDescription& set_desc ) - { - if( set_desc.minimal_spacing < GLOBAL_EPSILON ) - { - return; - } - const auto set_id = set_name_to_uuid_.at( set_desc.name ); - auto interaction = std::make_unique< - EuclideanCutoffInteraction< OwnerSegment2D > >( - set_desc.minimal_spacing, - PairwiseInteraction< OwnerSegment2D >::SCOPE::same_set ); - - this->ordered_energy_terms_.push_back( - this->energy_terms_collection_.add_energy_term( - std::make_unique< PairwiseTerm< OwnerSegment2D > >( - absl::StrCat( set_desc.name, "_min_spacing" ), 0., - std::vector< uuid >{ set_id }, std::move( interaction ), - this->domain_ ) ) ); - } - - void set_x_intersection_term() - { - if( std::fabs( beta_x_node_ - 1. ) > GLOBAL_EPSILON ) - { - std::vector< uuid > set_uuids; - set_uuids.reserve( set_name_to_uuid_.size() ); - for( const auto& [name, id] : set_name_to_uuid_ ) - { - set_uuids.push_back( id ); - } - auto interaction = std::make_unique< - EuclideanCutoffInteraction< OwnerSegment2D > >( - 0., PairwiseInteraction< - OwnerSegment2D >::SCOPE::different_set ); - - this->ordered_energy_terms_.push_back( - this->energy_terms_collection_.add_energy_term( - std::make_unique< PairwiseTerm< OwnerSegment2D > >( - absl::StrCat( "inter_set_x_nodes" ), beta_x_node_, - set_uuids, std::move( interaction ), - this->domain_ ) ) ); - } - } - - private: - std::vector< FractureSetDescription > set_descriptors_; - std::unordered_map< std::string, uuid > set_name_to_uuid_; - // x node monitoring - double beta_x_node_{ 1. }; - }; - -} // namespace geode \ No newline at end of file diff --git a/src/geode/stochastic/CMakeLists.txt b/src/geode/stochastic/CMakeLists.txt index 81d5776..1790309 100644 --- a/src/geode/stochastic/CMakeLists.txt +++ b/src/geode/stochastic/CMakeLists.txt @@ -71,7 +71,6 @@ 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/simulation_context.hpp" "sampling/mcmc/helpers/simulation_printer.hpp" "models/energy_terms/energy_term.hpp" diff --git a/src/geode/stochastic/applications/fractures.cpp b/src/geode/stochastic/applications/fractures.cpp index 9d2c7fb..f592b92 100644 --- a/src/geode/stochastic/applications/fractures.cpp +++ b/src/geode/stochastic/applications/fractures.cpp @@ -37,6 +37,20 @@ namespace } return fractures; } + + std::vector< std::pair< std::string, std::string > > inter_set_interactions( + const std::vector< std::string >& set_names ) + { + std::vector< std::pair< std::string, std::string > > interactions; + 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] ); + } + } + return interactions; + } } // namespace namespace geode @@ -45,6 +59,8 @@ namespace geode using FractureIntensityDescription = SingleObjectTermConfig; using FractureSpacingDescription = PairwiseTermConfig; + using XNodeInteractionDescription = geode::PairwiseTermConfig; + using FractureSimulationConfig = SimulationContextConfig< Fracture >; FractureSimulationContext build_fractures_simulation_context( @@ -53,10 +69,14 @@ namespace geode FractureSimulationConfig simulation_config; simulation_config.domain = fnet_desc.domain; + std::vector< std::string > set_names; + set_names.reserve( fnet_desc.fracture_sets.size() ); for( const auto& fset_desc : fnet_desc.fracture_sets ) { auto& fset = simulation_config.add_set( fset_desc.fset_name ); + set_names.emplace_back( 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; @@ -91,6 +111,18 @@ namespace geode geode::MinimalDistanceCutoffConfig{ fset_desc.minimal_spacing }; simulation_config.model.terms.emplace_back( std::move( spacing ) ); } + if( set_names.size() > 1 ) + { + XNodeInteractionDescription interaction; + interaction.term_name = fnet_desc.x_node_interaction_name(); + interaction.object_set_names_interactions = + inter_set_interactions( set_names ); + interaction.gamma = fnet_desc.beta_x_node; + interaction.interaction_config = + geode::MinimalDistanceCutoffConfig{ 0. }; + simulation_config.model.terms.emplace_back( + std::move( interaction ) ); + } return build_simulation_context( simulation_config ); } diff --git a/tests/stochastic/applications/test-fractures.cpp b/tests/stochastic/applications/test-fractures.cpp index b09e725..e7e2255 100644 --- a/tests/stochastic/applications/test-fractures.cpp +++ b/tests/stochastic/applications/test-fractures.cpp @@ -44,14 +44,6 @@ namespace // --- 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; @@ -79,7 +71,7 @@ namespace geode::FractureSimulationRunner runner{ std::move( context ) }; // run simulation geode::SimulationConfigurator sim_config; - sim_config.realizations = 2000; + sim_config.realizations = 500; sim_config.metropolis_hasting_steps = 100; sim_config.burn_in_steps = 1000; @@ -109,86 +101,88 @@ namespace 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!" ); - // } + 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-two-Fracture-set@" ); + + // NOLINTBEGIN(*-magic-numbers) + geode::FractureNetworkDescription fnet_desc; + fnet_desc.fnet_name = "Two_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_01 = fnet_desc.add_fracture_set( "fset_01" ); + fset_01.sampler.length.distribution_type = + geode::TruncatedPowerLaw::distribution_type_static(); + fset_01.sampler.length.alpha = 2.; + fset_01.sampler.length.min_value = 1; + fset_01.sampler.length.max_value = 10.; + + fset_01.sampler.azimuth.distribution_type = + geode::UniformClosed< double >::distribution_type_static(); + fset_01.sampler.azimuth.min_value = 1; + fset_01.sampler.azimuth.max_value = 10.; + + fset_01.p20 = 0.05; + fset_01.p21 = 200; + fset_01.minimal_spacing = 1.; + + auto& fset_02 = fnet_desc.add_fracture_set( "fset_02" ); + fset_02.sampler.length.distribution_type = + geode::TruncatedLogNormal::distribution_type_static(); + fset_02.sampler.length.mean = 1; + fset_02.sampler.length.standard_deviation = 1.; + fset_02.sampler.length.min_value = 1; + fset_02.sampler.length.max_value = 50.; + + fset_02.sampler.azimuth.distribution_type = + geode::VonMises::distribution_type_static(); + fset_02.sampler.azimuth.mean = 60.; + fset_02.sampler.azimuth.kappa = 1.; + + fset_02.p20 = 0.05; + fset_02.p21 = 200; + fset_02.minimal_spacing = 2.; + + fnet_desc.add_x_node_monitoring( 0.3 ); + + // runner + auto context = build_fractures_simulation_context( fnet_desc ); + geode::FractureSimulationRunner runner{ std::move( context ) }; + // run simulation + geode::SimulationConfigurator sim_config; + sim_config.realizations = 500; + 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::Logger::info( "--> SUCCESS!" ); + } } // namespace int main() @@ -198,7 +192,7 @@ int main() geode::OpenGeodeStochasticStochasticLibrary::initialize(); geode::Logger::set_level( geode::Logger::LEVEL::debug ); test_fracture_simulator(); - // test_two_fracture_sets_simulator(); + test_two_fracture_sets_simulator(); return 0; } catch( ... )