Skip to content

Commit 8a925a1

Browse files
francoisbonneaufrancoisbonneau
committed
feat(Proposal): make move agnostic from the subset id
1 parent 335d6f7 commit 8a925a1

9 files changed

Lines changed: 399 additions & 313 deletions

File tree

‎include/geode/stochastic/common.hpp‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -31,4 +31,7 @@
3131
namespace geode
3232
{
3333
OPENGEODE_LIBRARY( opengeode_stochastic_stochastic_api, Stochastic );
34+
35+
static constexpr double LOG_PROB_INVALID =
36+
-std::numeric_limits< double >::infinity();
3437
} // namespace geode

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

Lines changed: 10 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -32,34 +32,24 @@ namespace geode
3232
class ObjectSetSampler
3333
{
3434
public:
35-
ObjectSetSampler( const uuid& subset_id ) : subset_id_( subset_id ) {}
36-
virtual ~ObjectSetSampler() = default;
35+
[[nodiscard]] virtual Type sample( RandomEngine& engine ) const = 0;
3736

38-
virtual std::pair< Type, uuid > sample(
39-
RandomEngine& engine ) const = 0;
37+
[[nodiscard]] virtual Type change(
38+
const Type& object, RandomEngine& engine ) const = 0;
4039

41-
std::optional< ObjectId > sample_id(
42-
const ObjectSet< Type >& config, RandomEngine& engine ) const
40+
[[nodiscard]] double log_pdf( const Type& obj ) const
4341
{
44-
const auto max_obj_id = config.nb_objects_in_subset( subset_id_ );
45-
if( max_obj_id == 0 )
42+
if( !is_valid_object( obj ) )
4643
{
47-
return std::nullopt;
44+
return LOG_PROB_INVALID;
4845
}
49-
geode::UniformClosed< index_t > uniform_closed_index_t;
50-
uniform_closed_index_t.min_value = 0;
51-
uniform_closed_index_t.max_value = max_obj_id - 1;
52-
ObjectId result{ engine.sample_uniform( uniform_closed_index_t ),
53-
subset_id_ };
54-
return result;
46+
return log_pdf_;
5547
}
5648

57-
virtual std::pair< Type, uuid > change(
58-
const Type& object, RandomEngine& engine ) const = 0;
59-
60-
virtual double log_pdf( const Type& obj ) const = 0;
49+
protected:
50+
virtual bool is_valid_object( const Type& obj ) const = 0;
6151

6252
protected:
63-
const uuid& subset_id_;
53+
double log_pdf_{ LOG_PROB_INVALID };
6454
};
6555
} // namespace geode

‎include/geode/stochastic/sampling/direct/object_set_sampler/point_set_sampler.hpp‎

Lines changed: 16 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -35,25 +35,21 @@ namespace geode
3535
class UniformPointSetSampler : public ObjectSetSampler< Point< dimension > >
3636
{
3737
public:
38-
UniformPointSetSampler(
39-
const BoundingBox< dimension >& box, const uuid& subset_id )
40-
: ObjectSetSampler< Point< dimension > >( subset_id ), box_( box )
38+
UniformPointSetSampler( const BoundingBox< dimension >& box )
39+
: ObjectSetSampler< Point< dimension > >{}, box_( box )
4140
{
4241
auto volume = box_.n_volume();
43-
if( volume != 0. )
44-
{
45-
log_pdf_ = -std::log( volume );
46-
}
42+
OPENGEODE_EXCEPTION( volume != 0.,
43+
"[PointSetSampler] - Undefined Bounding Box (volume ==0)." );
44+
this->log_pdf_ = -std::log( volume );
4745
}
4846

49-
std::pair< Point< dimension >, uuid > sample(
50-
RandomEngine& engine ) const override
47+
Point< dimension > sample( RandomEngine& engine ) const override
5148
{
52-
return { PointUniformSampler::sample< dimension >( engine, box_ ),
53-
this->subset_id_ };
49+
return PointUniformSampler::sample< dimension >( engine, box_ );
5450
}
5551

56-
std::pair< Point< dimension >, uuid > change(
52+
Point< dimension > change(
5753
const Point< dimension >& obj, RandomEngine& engine ) const override
5854
{
5955
double ratio = 0.1;
@@ -67,30 +63,25 @@ namespace geode
6763
{
6864
if( box_.contains( new_point ) )
6965
{
70-
return { new_point, this->subset_id_ };
66+
return new_point;
7167
}
7268
new_point =
7369
PointUniformSampler::sample< dimension >( engine, ball );
7470
}
75-
OPENGEODE_EXCEPTION( true == false,
76-
absl::StrCat(
77-
"[PointSampler] - Cannot find a point in the box: ",
78-
box_.string() ) );
79-
return { obj, this->subset_id_ };
71+
throw OpenGeodeException( absl::StrCat(
72+
"[PointSampler] - Cannot find a point in the box: ",
73+
box_.string() ) );
74+
return obj;
8075
}
8176

82-
double log_pdf( const Point< dimension >& obj ) const override
77+
private:
78+
bool is_valid_object( const Point< dimension >& obj ) const override
8379
{
84-
if( !box_.contains( obj ) )
85-
{
86-
return -std::numeric_limits< double >::infinity();
87-
}
88-
return log_pdf_;
80+
return box_.contains( obj );
8981
}
9082

9183
private:
9284
const BoundingBox< dimension >& box_;
93-
double log_pdf_{ -std::numeric_limits< double >::infinity() };
9485
};
9586

9687
} // namespace geode

‎include/geode/stochastic/sampling/mcmc/metropolis_hasting_sampler.hpp‎

Lines changed: 52 additions & 104 deletions
Original file line numberDiff line numberDiff line change
@@ -36,69 +36,50 @@ namespace geode
3636
Undecided
3737
};
3838

39-
template < typename Type >
39+
template < typename ObjectType >
4040
struct StepResult
4141
{
4242
MHDecision decision{ MHDecision::Undecided };
43-
typename Proposal< Type >::Move move_type{
44-
Proposal< Type >::Move::Invalid
45-
};
43+
MoveType move_type{ MoveType::Invalid };
4644
double log_accept{ -std::numeric_limits< double >::infinity() };
4745
double delta_log_energy{ 0.0 };
4846
};
4947

50-
template < typename Type >
48+
template < typename ObjectType >
5149
class MetropolisHastings
5250
{
5351
public:
54-
MetropolisHastings( GibbsEnergy< Type >& energy,
55-
std::unique_ptr< ProposalKernel< Type > > proposal_kernel )
52+
MetropolisHastings( GibbsEnergy< ObjectType >& energy,
53+
std::unique_ptr< ProposalKernel< ObjectType > > proposal_kernel )
5654
: energy_( energy ),
5755
proposal_kernel_( std::move( proposal_kernel ) )
5856
{
5957
OPENGEODE_ASSERT(
6058
proposal_kernel_ != nullptr, "[MH] null proposal kernel" );
6159
}
6260

63-
ObjectSet< Type > initialize_object_set_with_sampling(
64-
RandomEngine& engine,
65-
const absl::flat_hash_map< uuid, index_t >& group_targets ) const
66-
{
67-
ObjectSet< Type > config;
68-
for( const auto& [subset_id, target] : group_targets )
69-
{
70-
config.add_subset( subset_id );
71-
while( config.nb_objects_in_subset( subset_id ) < target )
72-
{
73-
OPENGEODE_EXCEPTION( try_birth( config, subset_id, engine ),
74-
"[MH] Birth move need to be more probable for group: ",
75-
subset_id.string() );
76-
}
77-
}
78-
return config;
79-
}
80-
81-
StepResult< Type > step(
82-
ObjectSet< Type >& state, RandomEngine& engine ) const
61+
StepResult< ObjectType > step(
62+
ObjectSet< ObjectType >& state, RandomEngine& engine ) const
8363
{
84-
Proposal< Type > proposal =
64+
Proposal< ObjectType > proposal =
8565
proposal_kernel_->propose( state, engine );
86-
if( proposal.type == Proposal< Type >::Move::Birth )
66+
const auto& move_type = proposal.proposed_move.type;
67+
if( move_type == MoveType::Birth )
8768
{
8869
return birth_step( proposal, state, engine );
8970
}
90-
if( proposal.type == Proposal< Type >::Move::Death )
71+
if( move_type == MoveType::Death )
9172
{
9273
return death_step( proposal, state, engine );
9374
}
94-
if( proposal.type == Proposal< Type >::Move::Change )
75+
if( move_type == MoveType::Change )
9576
{
9677
return change_step( proposal, state, engine );
9778
}
98-
return StepResult< Type >{};
79+
return StepResult< ObjectType >{};
9980
}
10081

101-
void walk( ObjectSet< Type >& state,
82+
void walk( ObjectSet< ObjectType >& state,
10283
RandomEngine& engine,
10384
index_t nb_steps ) const
10485
{
@@ -109,7 +90,7 @@ namespace geode
10990
}
11091
}
11192

112-
ObjectSet< Type > walk_copy( ObjectSet< Type > initial,
93+
ObjectSet< ObjectType > walk_copy( ObjectSet< ObjectType > initial,
11394
RandomEngine& engine,
11495
index_t nb_steps ) const
11596
{
@@ -162,57 +143,34 @@ namespace geode
162143
return 0.0;
163144
if( log_accept >= 0.0 )
164145
return 1.0;
165-
// prevent expoential overflow
146+
// prevent exponential overflow
166147
constexpr double LOG_MIN = -745.0;
167148
if( log_accept < LOG_MIN )
168149
return 0.0;
169150
return std::exp( log_accept );
170151
}
171152

172153
private:
173-
bool try_birth( ObjectSet< Type >& config,
174-
const uuid& subset_id,
175-
RandomEngine& engine ) const
154+
const double compute_log_accept( const double deltaU,
155+
const ProposalProbabilities& proposal_probas ) const
176156
{
177-
// pbirth_should be > 0.01
178-
for( const auto attempt : geode::Range{ 100 } )
179-
{
180-
geode_unused( attempt );
181-
auto proposal = proposal_kernel_->propose( config, engine );
182-
if( proposal.type == Proposal< Type >::Move::Birth
183-
&& proposal.new_object->second == subset_id )
184-
{
185-
config.add_object( std::move( proposal.new_object->first ),
186-
proposal.new_object->second );
187-
return true;
188-
}
189-
}
190-
return false;
191-
}
192-
193-
const double compute_log_accept(
194-
const double deltaU, const Proposal< Type >& proposal ) const
195-
{
196-
OPENGEODE_ASSERT(
197-
std::isfinite( proposal.log_forward_prob )
198-
&& std::isfinite( proposal.log_backward_prob ),
199-
"[MH] Non-finite proposal log-probabilities" );
200-
return -beta_ * deltaU + proposal.log_backward_prob
201-
- proposal.log_forward_prob;
157+
return -beta_ * deltaU + proposal_probas.transition_probability();
202158
}
203159

204160
template < typename ApplyMove >
205-
StepResult< Type > accept_or_reject( Proposal< Type >& proposal,
206-
ObjectSet< Type >& state,
161+
StepResult< ObjectType > accept_or_reject(
162+
Proposal< ObjectType >& proposal,
163+
ObjectSet< ObjectType >& state,
207164
RandomEngine& engine,
208165
const double delta_log_energy,
209166
ApplyMove&& apply_move ) const
210167
{
211-
StepResult< Type > step_result;
212-
step_result.move_type = proposal.type;
168+
const auto& proposed_move = proposal.proposed_move;
169+
StepResult< ObjectType > step_result;
170+
step_result.move_type = proposed_move.type;
213171
step_result.delta_log_energy = delta_log_energy;
214-
step_result.log_accept =
215-
compute_log_accept( delta_log_energy, proposal );
172+
step_result.log_accept = compute_log_accept(
173+
delta_log_energy, proposed_move.proposal_probabilities );
216174

217175
double log_u = engine.sample_log();
218176
step_result.decision = ( log_u < step_result.log_accept )
@@ -224,65 +182,55 @@ namespace geode
224182
return step_result;
225183
}
226184

227-
StepResult< Type > birth_step( Proposal< Type >& proposal,
228-
ObjectSet< Type >& state,
185+
StepResult< ObjectType > birth_step( Proposal< ObjectType >& proposal,
186+
ObjectSet< ObjectType >& state,
229187
RandomEngine& engine ) const
230188
{
231-
OPENGEODE_ASSERT( proposal.new_object.has_value(),
232-
"[MH] Birth proposal has no new_object" );
233-
geode::ObjectRef< Type > new_object{
234-
proposal.new_object.value().first,
235-
proposal.new_object.value().second
236-
};
189+
const auto new_object = proposal.new_object();
237190
const auto delta_log_energy =
238191
energy_.delta_log_energy_add( state, new_object );
239192
return accept_or_reject( proposal, state, engine, delta_log_energy,
240-
[]( auto& s, auto& p ) {
241-
s.add_object( std::move( p.new_object.value().first ),
242-
p.new_object.value().second );
193+
[]( auto& state, auto& proposal ) {
194+
state.add_object(
195+
std::move( proposal.proposed_move.new_object.value() ),
196+
proposal.subset_id );
243197
} );
244198
};
245199

246-
StepResult< Type > death_step( Proposal< Type >& proposal,
247-
ObjectSet< Type >& state,
200+
StepResult< ObjectType > death_step( Proposal< ObjectType >& proposal,
201+
ObjectSet< ObjectType >& state,
248202
RandomEngine& engine ) const
249203
{
250-
OPENGEODE_ASSERT( proposal.old_object_id.has_value(),
251-
"[MH] Death proposal has no index" );
252-
const auto delta_log_energy = energy_.delta_log_energy_remove(
253-
state, proposal.old_object_id.value() );
204+
const auto old_object_id = proposal.old_object_id();
205+
const auto delta_log_energy =
206+
energy_.delta_log_energy_remove( state, old_object_id );
254207
return accept_or_reject( proposal, state, engine, delta_log_energy,
255-
[]( auto& s, auto& p ) {
256-
s.remove_object( p.old_object_id.value() );
208+
[]( auto& state, auto& proposal ) {
209+
state.remove_object( proposal.old_object_id() );
257210
} );
258211
};
259212

260-
StepResult< Type > change_step( Proposal< Type >& proposal,
261-
ObjectSet< Type >& state,
213+
StepResult< ObjectType > change_step( Proposal< ObjectType >& proposal,
214+
ObjectSet< ObjectType >& state,
262215
RandomEngine& engine ) const
263216
{
264-
OPENGEODE_ASSERT( proposal.new_object.has_value(),
265-
"[MH] Change proposal has no new_object" );
266-
OPENGEODE_ASSERT( proposal.old_object_id.has_value(),
267-
"[MH] Change proposal has no index" );
268-
geode::ObjectRef< Type > new_object{
269-
proposal.new_object.value().first,
270-
proposal.new_object.value().second
271-
};
217+
const auto new_object = proposal.new_object();
218+
const auto old_object_id = proposal.old_object_id();
272219
const auto delta_log_energy = energy_.delta_log_energy_change(
273-
state, proposal.old_object_id.value(), new_object );
220+
state, old_object_id, new_object );
274221
// should we test that objects are in the same group?
275222
// should be ensured by the dynamic
276223
return accept_or_reject( proposal, state, engine, delta_log_energy,
277-
[]( auto& s, auto& p ) {
278-
s.update_object( p.old_object_id.value(),
279-
std::move( p.new_object.value().first ) );
224+
[]( auto& state, auto& proposal ) {
225+
state.update_object( proposal.old_object_id(),
226+
std::move(
227+
proposal.proposed_move.new_object.value() ) );
280228
} );
281229
};
282230

283231
private:
284-
const GibbsEnergy< Type >& energy_;
285-
std::unique_ptr< ProposalKernel< Type > > proposal_kernel_;
232+
const GibbsEnergy< ObjectType >& energy_;
233+
std::unique_ptr< ProposalKernel< ObjectType > > proposal_kernel_;
286234
double beta_{ 1.0 };
287235
};
288236
} // namespace geode

0 commit comments

Comments
 (0)