Skip to content

feat(FractureModel): add p21 controle #30

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
BotellaA merged 4 commits into from
Jan 5, 2026
Merged

feat(FractureModel): add p21 controle #30

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
a68ffce
feat(FractureModel): add p21 controle
francoisbonneau Dec 22, 2025
f427918
add binding
francoisbonneau Dec 22, 2025
9d53c3b
add binding
francoisbonneau Dec 22, 2025
ab52b57
review arnaud
francoisbonneau Jan 5, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions bindings/python/src/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@ namespace geode
.def_readwrite( "length", &FractureSetDescription::length )
.def_readwrite( "azimuth", &FractureSetDescription::azimuth )
.def_readwrite( "p20", &FractureSetDescription::p20 )
.def_readwrite( "p21", &FractureSetDescription::p21 )
.def( "add_observed_fracture",
[]( FractureSetDescription& self, const geode::Point2D& a,
const geode::Point2D& b ) {
Expand Down
1 change: 1 addition & 0 deletions bindings/python/tests/stochastic/test-py-mh-fractures.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,7 @@ def test_fracture_simulator():

# positioning
setA.p20 = 0.06
setA.p21 = 10
setA.minimal_spacing = 1.0

runner = stochastic.FractureSimulationRunner(domain)
Expand Down
128 changes: 81 additions & 47 deletions include/geode/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@
#include <geode/stochastic/sampling/direct/object_set_sampler/segment_set_sampler.hpp>
#include <geode/stochastic/sampling/mcmc/helpers/simulation_runner.hpp>
#include <geode/stochastic/sampling/mcmc/models/components/density_term.hpp>
#include <geode/stochastic/sampling/mcmc/models/components/intensity_term.hpp>
#include <geode/stochastic/sampling/mcmc/models/components/pairwise_term.hpp>
#include <geode/stochastic/sampling/mcmc/proposal/classical_proposals.hpp>
#include <geode/stochastic/spatial/pairwise_interactions.hpp>
Expand All @@ -48,7 +49,7 @@ namespace geode
// positionning
double p20;
std::vector< std::array< geode::Point2D, 2 > > observed_fractures;
// double p21;
double p21{ 1. };
double minimal_spacing{ 0. };

// mh dynamique
Expand Down Expand Up @@ -107,7 +108,8 @@ namespace geode
std::make_unique< ProposalKernel< OwnerSegment2D > >();

// Mapping set names -> UUID
std::unordered_map< std::string, uuid > name_to_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_ )
Expand All @@ -120,7 +122,7 @@ namespace geode
fixed_object[0], fixed_object[1] },
set_id, true );
}
name_to_uuid[set_desc.name] = set_id;
set_name_to_uuid_[set_desc.name] = set_id;

auto length_distribution =
DoubleSampler::create_distribution( set_desc.length );
Expand All @@ -139,52 +141,12 @@ namespace geode
// Step 2: create density energy terms
for( const auto& set_desc : set_descriptors_ )
{
const auto set_id = name_to_uuid.at( set_desc.name );
// p20
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_ ) ) );
// spacing
if( set_desc.minimal_spacing > GLOBAL_EPSILON )
{
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_ ) ) );
}
set_density_term( set_desc );
set_intensity_term( set_desc );
set_minimal_spacing_term( set_desc );
}
// x node monitoring
if( std::fabs( beta_x_node_ - 1. ) > GLOBAL_EPSILON )
{
std::vector< uuid > set_uuids;
set_uuids.reserve( name_to_uuid.size() );
for( const auto& [name, id] : name_to_uuid )
{
set_uuids.push_back( id );
}
auto interaction = std::make_unique<
EuclideanCutoffInteraction< OwnerSegment2D > >(
0., PairwiseInteraction<
OwnerSegment2D >::SCOPE::different_set );
set_x_intersection_term();

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_ ) ) );
}
this->mh_sampler_ =
std::make_unique< MetropolisHastings< OwnerSegment2D > >(
this->energy_terms_collection_,
Expand Down Expand Up @@ -229,8 +191,80 @@ namespace geode
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. };
};
Expand Down
84 changes: 22 additions & 62 deletions include/geode/stochastic/sampling/mcmc/models/components/density_term.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,79 +24,39 @@

#include <geode/stochastic/spatial/object_sets.hpp>

#include <geode/stochastic/sampling/mcmc/models/components/energy_term.hpp>
#include <geode/stochastic/sampling/mcmc/models/components/single_object_term.hpp>

namespace geode
{
/// ObjectCountTerm
template < typename ObjectType >
class DensityTerm : public EnergyTerm< 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 )
: EnergyTerm< ObjectType >(
name, lambda, std::move( targeted_set_ids ), 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 >& domain ) {
if( SpatialDomainChecker<
ObjectType >::is_anchored_in_domain( domain,
obj ) )
{
return 1.0;
}
return 0.0;
}, // contribution = 1 anchoredin domain
domain )
{
}

double total_log( const ObjectSets< ObjectType >& state ) const override
{
const auto n = this->statistic( state );
return this->contribution( n );
}

double delta_log_add( const ObjectSets< ObjectType >& /*state*/,
const ObjectRef< ObjectType >& new_object ) const override
{
if( !this->is_targeted_set( new_object.set_id )
|| !this->is_anchored_in_domain( new_object.object ) )
{
return 0.0;
}
return this->contribution( 1.0 );
}

double delta_log_remove( const ObjectSets< ObjectType >& state,
const ObjectId& object_id ) const override
{
if( !this->is_targeted_set( object_id.set_id )
|| !this->is_anchored_in_domain(
state.get_object( object_id ) ) )
{
return 0.0;
}
return this->contribution( -1.0 );
}

double delta_log_change( const ObjectSets< ObjectType >& state,
const ObjectId& old_object_id,
const ObjectRef< ObjectType >& new_object ) const override
{
auto old_in = this->is_anchored_in_domain(
state.get_object( old_object_id ) );
auto new_in = this->is_anchored_in_domain( new_object.object );
if( old_in == new_in )
{
return 0.0;
}
return this->contribution( new_in ? 1.0 : -1.0 );
}

double statistic( const ObjectSets< ObjectType >& state ) const override
{
double sum = 0.0;
this->for_each_targeted_object(
state, [&]( const ObjectId& obj_id ) {
if( this->is_anchored_in_domain(
state.get_object( obj_id ) ) )
{
sum += 1.0;
}
} );
return sum;
}
};
} // namespace geode
18 changes: 9 additions & 9 deletions include/geode/stochastic/sampling/mcmc/models/components/energy_term.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -156,17 +156,17 @@ namespace geode
targeted_set_ids_.begin(), targeted_set_ids_.end(), set_id );
}

bool is_anchored_in_domain( const ObjectType& obj ) const
const SpatialDomain< ObjectType::dim >& domain() const
{
return SpatialDomainChecker< ObjectType >::is_anchored_in_domain(
domain_, obj );
return domain_;
}

bool intersects_domain( const ObjectType& obj ) const
{
return SpatialDomainChecker< ObjectType >::intersects_domain(
domain_, obj );
}
// bool intersects_domain( const ObjectType& obj ) const
// {
// return SpatialDomainChecker< ObjectType
// >::intersects_domain(
// domain_, obj );
// }

template < typename Func >
void for_each_targeted_object(
Expand All @@ -185,6 +185,6 @@ namespace geode
private:
detail::EnergyScale energy_scale_;
std::vector< uuid > targeted_set_ids_;
const SpatialDomain< ObjectType::dim > domain_;
const SpatialDomain< ObjectType::dim >& domain_;
};
} // namespace geode
Loading