Skip to content

Commit aa47dbe

Browse files
francoisbonneaufrancoisbonneau
committed
fix(MH): add convergence test
1 parent 23c514e commit aa47dbe

5 files changed

Lines changed: 167 additions & 22 deletions

File tree

‎include/geode/stochastic/sampling/mcmc/proposal/moves.hpp‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -85,7 +85,8 @@ namespace geode
8585
{
8686
OPENGEODE_EXCEPTION( birth_ratio_ > 0. && birth_ratio_ < 1.,
8787
"[BirthDeathMove]-the ratio of birth over mover should be in "
88-
"]0,1[." );
88+
"]0,1[. (here = ",
89+
birth_ratio_, ")" );
8990
log_p_birth_ = std::log( this->p_move_ * birth_ratio );
9091
log_p_death_ = std::log( this->p_move_ * ( 1.0 - birth_ratio ) );
9192
}

‎include/geode/stochastic/sampling/mcmc/proposal/proposal_kernel.hpp‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -42,6 +42,10 @@ namespace geode
4242
{
4343
OPENGEODE_EXCEPTION( !moves_.empty(),
4444
"[MCMC Proposal Kernel] - no move are defined in the Kernel." );
45+
if( cumulative_probs_.size() == 1 )
46+
{
47+
return moves_[0]->propose_move( current, engine );
48+
}
4549
auto rnd = engine.sample_uniform( uniform_closed_double_ );
4650
for( const auto proba_id : Range{ cumulative_probs_.size() } )
4751
{
@@ -88,10 +92,6 @@ namespace geode
8892
} );
8993
cumulative_probs_.back() = 1.0;
9094
Logger::info( " PROPOSAL: " );
91-
for( const auto proba : cumulative_probs_ )
92-
{
93-
Logger::info( proba );
94-
}
9595
}
9696

9797
private:

‎tests/stochastic/CMakeLists.txt‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -104,6 +104,13 @@ add_geode_test(
104104
OpenGeode::geometry
105105
${PROJECT_NAME}::stochastic
106106
)
107+
add_geode_test(
108+
SOURCE "sampling/mcmc/test-mh-poisson-statistics-and-convergence.cpp"
109+
DEPENDENCIES
110+
OpenGeode::basic
111+
OpenGeode::geometry
112+
${PROJECT_NAME}::stochastic
113+
)
107114

108115
add_geode_test(
109116
SOURCE "sampling/test-random-engine.cpp"

‎tests/stochastic/sampling/mcmc/test-metropolis-hasting-sampler.cpp‎

Lines changed: 1 addition & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -71,23 +71,9 @@ namespace
7171
exception_thrown, "[MH test] negative beta did not throw." );
7272
}
7373

74-
geode::Configuration< geode::Point2D > create_configuration()
75-
{
76-
geode::Point2D p1{ { 0., 0. } };
77-
geode::MarkedObject< geode::Point2D > mp1{ std::move( p1 ) };
78-
geode::Point2D p2{ { 1., 1. } };
79-
geode::MarkedObject< geode::Point2D > mp2{ std::move( p2 ) };
80-
81-
geode::Configuration< geode::Point2D > pattern;
82-
pattern.add_object( std::move( mp1 ) );
83-
pattern.add_object( std::move( mp2 ) );
84-
85-
return pattern;
86-
}
87-
8874
void test_steps( const geode::MetropolisHastings< geode::Point2D >& mh )
8975
{
90-
auto state = create_configuration();
76+
geode::Configuration< geode::Point2D > state;
9177
geode::RandomEngine engine;
9278

9379
geode::index_t stat_sum{ 0 };
@@ -175,7 +161,6 @@ int main()
175161

176162
geode::UniformMarkedPointSampler< 2 > sampler( box, std::nullopt );
177163

178-
// Create classical birth-death-change kernel
179164
auto kernel = geode::create_birth_death_change_kernel< geode::Point2D >(
180165
sampler, 0.1, 0.1 );
181166

@@ -187,7 +172,6 @@ int main()
187172

188173
geode::MetropolisHastings< geode::Point2D > mh(
189174
poisson_energy, std::move( kernel ) );
190-
191175
test_steps( mh );
192176
test_beta_setter( mh );
193177
test_acceptance_prob_helper();

‎tests/stochastic/sampling/mcmc/test-mh-poisson-statistics-and-convergence.cpp‎

Lines changed: 153 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,153 @@
1+
/*
2+
* Copyright (c) 2019 - 2025 Geode-solutions
3+
*
4+
* Permission is hereby granted, free of charge, to any person obtaining a copy
5+
* of this software and associated documentation files (the "Software"), to deal
6+
* in the Software without restriction, including without limitation the rights
7+
* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
8+
* copies of the Software, and to permit persons to whom the Software is
9+
* furnished to do so, subject to the following conditions:
10+
*
11+
* The above copyright notice and this permission notice shall be included in
12+
* all copies or substantial portions of the Software.
13+
*
14+
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
15+
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
16+
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
17+
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
18+
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
19+
* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
20+
* SOFTWARE.
21+
*
22+
*/
23+
#include <geode/geometry/point.hpp>
24+
#include <geode/stochastic/sampling/mcmc/metropolis_hasting_sampler.hpp>
25+
#include <geode/stochastic/sampling/mcmc/models/components/intensity_term.hpp>
26+
#include <geode/stochastic/sampling/mcmc/models/gibbs_energy.hpp>
27+
#include <geode/stochastic/sampling/mcmc/proposal/classical_proposals.hpp>
28+
#include <geode/stochastic/sampling/mcmc/proposal/marked_object_sampler/uniform_marked_point_sampler.hpp>
29+
namespace
30+
{
31+
32+
// ------------------------------------------------------------
33+
// Convergence test: mean number of points ≈ λ × area
34+
// ------------------------------------------------------------
35+
void test_convergence( geode::MetropolisHastings< geode::Point2D >& mh,
36+
double expected_points,
37+
geode::RandomEngine& engine )
38+
{
39+
geode::Configuration< geode::Point2D > state;
40+
41+
constexpr geode::index_t burn_in{ 10000 };
42+
constexpr geode::index_t N{ 1000000 };
43+
44+
// Burn-in
45+
for( const auto i : geode::Range{ burn_in } )
46+
{
47+
geode_unused( i );
48+
mh.step( state, engine );
49+
}
50+
51+
// Sampling
52+
double sum_points = 0.0;
53+
double sum_sq = 0.0;
54+
for( const auto i : geode::Range{ N } )
55+
{
56+
mh.step( state, engine );
57+
auto n = static_cast< double >( state.size() );
58+
sum_points += n;
59+
sum_sq += n * n;
60+
}
61+
62+
double mean_points = sum_points / N;
63+
double var = ( sum_sq / N ) - ( mean_points * mean_points );
64+
65+
geode::Logger::info( "[MH test] mean points = ", mean_points,
66+
" and var = ", var, " (expected ", expected_points, ")" );
67+
68+
OPENGEODE_EXCEPTION(
69+
std::abs( mean_points - expected_points ) < 0.01 * expected_points,
70+
"[MH test] mean number of points not close to expected." );
71+
72+
// ------------------------------------------------------------
73+
// Variance test: Poisson => Var(N) ≈ E[N]
74+
// ------------------------------------------------------------
75+
OPENGEODE_EXCEPTION(
76+
std::abs( var - mean_points ) < 0.1 * expected_points,
77+
"[MH test] variance not close to Poisson expectation." );
78+
}
79+
80+
// ------------------------------------------------------------
81+
// Dynamic independence: different kernels → same stationary law
82+
// ------------------------------------------------------------
83+
void test_dynamic_independence(
84+
double domain_length, double poisson_density, double birth_ratio )
85+
{
86+
geode::RandomEngine engine;
87+
engine.set_seed( "@mh-test@" );
88+
89+
geode::Point2D min_point{ { 0., 0. } };
90+
geode::Point2D max_point{ { domain_length, domain_length } };
91+
92+
geode::BoundingBox2D box;
93+
box.add_point( min_point );
94+
box.add_point( max_point );
95+
96+
double area = domain_length * domain_length;
97+
98+
geode::UniformMarkedPointSampler< 2 > sampler( box, std::nullopt );
99+
100+
geode::GibbsEnergy< geode::Point2D > poisson_energy;
101+
poisson_energy.add_energy_term(
102+
std::make_unique< geode::IntensityTerm< geode::Point2D > >(
103+
poisson_density ) );
104+
105+
// Kernel with only birth/death
106+
auto kernel1 = geode::create_birth_death_kernel< geode::Point2D >(
107+
sampler, birth_ratio );
108+
geode::MetropolisHastings< geode::Point2D > mh1(
109+
poisson_energy, std::move( kernel1 ) );
110+
111+
// Kernel with birth/death/change
112+
auto kernel2 =
113+
geode::create_birth_death_change_kernel< geode::Point2D >(
114+
sampler, birth_ratio, 0.1 );
115+
geode::MetropolisHastings< geode::Point2D > mh2(
116+
poisson_energy, std::move( kernel2 ) );
117+
118+
geode::Logger::info(
119+
"[MH test] Testing kernel with birth/death only..." );
120+
test_convergence( mh1, poisson_density * area, engine );
121+
122+
geode::Logger::info(
123+
"[MH test] Testing kernel with birth/death/change..." );
124+
test_convergence( mh2, poisson_density * area, engine );
125+
}
126+
} // namespace
127+
128+
int main()
129+
{
130+
try
131+
{
132+
geode::StochasticLibrary::initialize();
133+
134+
geode::StochasticLibrary::initialize();
135+
double domain_length{ 10. };
136+
double poisson_density_lambda{ 0.5 };
137+
std::array< double, 7 > birth_ratio{ 0.1, 0.2, 0.3, 0.4, 0.5, 0.6,
138+
0.8 };
139+
for( const auto b : birth_ratio )
140+
{
141+
geode::Logger::info( b );
142+
test_dynamic_independence(
143+
domain_length, poisson_density_lambda, b );
144+
}
145+
146+
geode::Logger::info( "MH STATISTICS AND CONVERGENCE TEST SUCCESS" );
147+
return 0;
148+
}
149+
catch( ... )
150+
{
151+
return geode::geode_lippincott();
152+
}
153+
}

0 commit comments

Comments
 (0)