feat(applications): implement fracture application

Author: francoisbonneau

include/geode/stochastic/applications/fractures.hpp

@@ -0,0 +1,136 @@
+#pragma once
+
+#include <geode/stochastic/inference/target_statistics.hpp>
+#include <geode/stochastic/sampling/direct/object_set_sampler/segment_set_sampler.hpp>
+#include <geode/stochastic/sampling/mcmc/helpers/simulation_context.hpp>
+#include <geode/stochastic/sampling/mcmc/simulation_runner.hpp>
+
+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 };
+
+        std::string density_name() const
+        {
+            return absl::StrCat( fset_name, "_p20" );
+        }
+        double p20{ 0. };
+        std::optional< double > expected_number;
+
+        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;
+
+        std::string spacing_name() const
+        {
+            return absl::StrCat( fset_name, "_spacing" );
+        }
+        double minimal_spacing{ 0. };
+
+        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 FractureNetworkDescription
+    {
+        std::string fnet_name;
+
+        SpatialDomainConfig< 2 > domain;
+
+        std::vector< FractureSetDescription > fracture_sets;
+
+        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;
+        }
+
+        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;
+        //        }
+    };
+
+    FractureSimulationContext build_fractures_simulation_context(
+        const FractureNetworkDescription& description );
+
+    std::vector< geode::TargetStatisticConfig > build_fractures_targeted_stat(
+        const FractureNetworkDescription& description );
+
+} // namespace geode

include/geode/stochastic/sampling/direct/double_sampler.hpp

@@ -48,6 +48,13 @@ namespace geode
 
         struct DistributionDescription
         {
+            DistributionDescription( absl::string_view name_in )
+                : name{ name_in }
+            {
+            }
+
+            DistributionDescription() = default;
+
             std::string name{ "default_distribution" };
             DistributionType distribution_type;
 

include/geode/stochastic/sampling/direct/object_set_sampler/segment_set_sampler.hpp

@@ -37,8 +37,9 @@ namespace geode
         double move_ratio{ 0.1 };
         // NOLINTEND(*-magic-numbers)
 
-        DoubleSampler::DistributionDescription length;
-        DoubleSampler::DistributionDescription azimuth;
+        DoubleSampler::DistributionDescription length{ "length" };
+
+        DoubleSampler::DistributionDescription azimuth{ "azimuth" };
     };
 
     class UniformSegmentSetSampler : public ObjectSetSampler< OwnerSegment2D >

include/geode/stochastic/sampling/mcmc/helpers/simulation_context.hpp

@@ -116,6 +116,11 @@ namespace geode
         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,

include/geode/stochastic/spatial/spatial_domain.hpp

@@ -152,6 +152,16 @@ namespace geode
         Point< dimension > min_point;
         Point< dimension > max_point;
         double buffer_size{ 0.0 };
+        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 >

src/geode/stochastic/CMakeLists.txt

@@ -22,6 +22,7 @@ add_geode_library(
     NAME stochastic
     FOLDER "geode/stochastic"
     SOURCES
+        "applications/fractures.cpp"
         "applications/poisson_process.cpp"
         "applications/strauss_process.cpp"
         "spatial/object_set.cpp"
@@ -41,6 +42,7 @@ add_geode_library(
         "sampling/random_engine.cpp"
         "common.cpp"
     PUBLIC_HEADERS
+        "applications/fractures.hpp"
         "applications/poisson_process.hpp"
         "applications/strauss_process.hpp"
         #"inference/abc_shadow.hpp"

src/geode/stochastic/applications/fractures.cpp

@@ -0,0 +1,114 @@
+/*
+ * 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 <geode/stochastic/applications/fractures.hpp>
+#include <geode/stochastic/sampling/direct/object_set_sampler/segment_set_sampler.hpp>
+namespace
+{
+    std::vector< geode::Fracture > build_observed_fractures(
+        const std::vector< std::array< geode::Point2D, 2 > >&
+            fracture_extremities )
+    {
+        std::vector< geode::Fracture > fractures;
+        for( const auto& extremities : fracture_extremities )
+        {
+            fractures.emplace_back(
+                geode::OwnerSegment2D{ 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;
+            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

tests/stochastic/applications/CMakeLists.txt

@@ -18,13 +18,13 @@
 # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 # SOFTWARE.
 
-#add_geode_test(
-#     SOURCE "sampling/mcmc/test-mh-fractures.cpp"
-#     DEPENDENCIES
-#         OpenGeode::basic
-#         OpenGeode::geometry
-#         ${PROJECT_NAME}::stochastic
-#)
+add_geode_test(
+     SOURCE "test-fractures.cpp"
+     DEPENDENCIES
+         OpenGeode::basic
+         OpenGeode::geometry
+         ${PROJECT_NAME}::stochastic
+)
 
 add_geode_test(
       SOURCE "test-poisson-process.cpp"