19 #ifndef MCMCDA_GIBBS_PROPOSER_H_INCLUDED
20 #define MCMCDA_GIBBS_PROPOSER_H_INCLUDED
36 #include <boost/function.hpp>
45 template<
class Track,
class Lhood>
49 typedef typename Track::Element Element;
50 typedef std::pair<int, typename std::set<Element>::const_iterator>
53 const Lhood& m_likelihood;
56 mutable double previous_prior;
57 mutable double previous_likelihood;
66 const Lhood& likelihood,
67 int neighborhood_size = -1
70 m_likelihood(likelihood),
71 nsize(neighborhood_size),
72 previous_prior(-std::numeric_limits<double>::infinity())
81 const Lhood& likelihood,
82 const Track& default_track,
83 int neighborhood_size = -1
86 m_likelihood(likelihood),
87 def_track(default_track),
88 nsize(neighborhood_size),
89 previous_prior(-std::numeric_limits<double>::infinity())
96 boost::optional<double> operator()
108 const Location_pair& cur_location
122 template<
class Track>
127 std::mem_fun_ref(&std::set<typename Track::Element>::size));
136 template<
class Track,
class Lhood>
137 boost::optional<double> Gibbs_proposer<Track, Lhood>::operator()
143 using std::make_pair;
145 typedef typename Track::iterator Tr_it;
146 typedef typename Track::const_iterator Tr_cit;
155 Location_pair tap = get_time_and_place(var, tr_data);
157 typename std::set<Element>::const_iterator point_p = tap.second;
160 for(owner_p = w.begin(); owner_p != w.end(); owner_p++)
162 std::pair<Tr_cit, Tr_cit> pair_ps = owner_p->equal_range(t);
164 for(Tr_cit pair_p = pair_ps.first; pair_p != pair_ps.second; pair_p++)
166 if(pair_p->second == &(*point_p))
180 if(owner_p != w.end())
184 std::remove_copy(w.begin(), w.end(),
185 std::insert_iterator<Association<Track> >(w_p, w_p.begin()),
188 Track track = *owner_p;
189 std::pair<Tr_it, Tr_it> pair_ps = track.equal_range(t);
190 for(Tr_it pair_p = pair_ps.first; pair_p != pair_ps.second; pair_p++)
192 if(pair_p->second == &(*point_p))
212 std::vector<Track> affected_tracks;
213 std::remove_copy_if(w_p.begin(), w_p.end(),
214 std::back_inserter(affected_tracks),
221 m_likelihood.set_limits(t - nsize/2, t + nsize/2);
226 double likelihood_noise;
227 if(previous_prior != -std::numeric_limits<double>::infinity()
228 && was_noise && var != 0)
230 prior_noise = previous_prior;
231 likelihood_noise = previous_likelihood;
235 prior_noise = m_prior(w_p);
236 likelihood_noise = m_likelihood(w_p);
240 Track modified_track = def_track;
242 modified_track.insert(make_pair(t, &(*point_p)));
244 modified_p = w_p.insert(modified_track).first;
245 double prior_new = m_prior(w_p);
246 double likelihood_new = likelihood_noise
247 - m_likelihood.at_noise_point(*point_p)
248 + m_likelihood.at_track(modified_track);
250 w_p.erase(modified_p);
272 std::vector<double> densities(affected_tracks.size() + 2);
273 std::vector<double> priors(affected_tracks.size() + 2);
274 std::vector<double> likelihoods(affected_tracks.size() + 2);
276 priors[0] = prior_noise;
277 priors[1] = prior_new;
278 likelihoods[0] = likelihood_noise;
279 likelihoods[1] = likelihood_new;
280 densities[0] = prior_noise + likelihood_noise;
281 densities[1] = prior_new + likelihood_new;
282 for(
size_t i = 0;
i < affected_tracks.size();
i++)
284 modified_track = affected_tracks[
i];
285 Track original_track = modified_track;
287 modified_track.insert(make_pair(t, &(*point_p)));
288 likelihoods[
i + 2] = likelihood_noise
289 - m_likelihood.at_track(affected_tracks[
i])
290 - m_likelihood.at_noise_point(*point_p)
291 + m_likelihood.at_track(modified_track);
293 w_p.erase(original_track);
294 modified_p = w_p.insert(modified_track).first;
295 priors[i + 2] = m_prior(w_p);
296 w_p.erase(modified_p);
297 w_p.insert(original_track);
299 densities[i + 2] = priors[i + 2] + likelihoods[i + 2];
303 m_likelihood.reset_limits();
305 std::vector<double> probabilities(densities.begin(), densities.end());
306 if(affected_tracks.empty())
308 probabilities[0] += log(point_p->prob_noise);
309 probabilities[1] += log(1.0 - point_p->prob_noise);
312 double mx = *std::max_element(probabilities.begin(), probabilities.end());
313 std::transform(probabilities.begin(), probabilities.end(),
314 probabilities.begin(),
315 std::bind2nd(std::minus<double>(), mx));
316 std::transform(probabilities.begin(), probabilities.end(),
317 probabilities.begin(),
318 std::ptr_fun<double, double>(std::exp));
321 probabilities.end(), 0.0);
322 std::transform(probabilities.begin(), probabilities.end(),
323 probabilities.begin(),
324 std::bind2nd(std::divides<double>(), total));
347 modified_track.clear();
349 modified_track.insert(make_pair(t, &(*point_p)));
350 w_p.insert(modified_track);
354 w_p.erase(affected_tracks[idx - 2]);
356 affected_tracks[idx - 2].insert(make_pair(t, &(*point_p)));
357 w_p.insert(affected_tracks[idx - 2]);
360 previous_prior = priors[idx];
361 previous_likelihood = likelihoods[idx];
363 return densities[idx];
368 template<
class Track,
class Lhood>
372 const Location_pair& cur_location
375 int t = cur_location.first;
376 typename std::set<Element>::const_iterator point_p = cur_location.second;
377 const double v_bar = m_prior.get_v_bar();
378 const int d_bar = m_prior.get_d_bar();
381 typename Track::const_iterator lb_pair_p = track.lower_bound(t);
383 if(lb_pair_p == track.end())
385 typename Track::const_reverse_iterator last_pair_p = track.rbegin();
386 d = t - last_pair_p->first;
388 return is_neighbor(*last_pair_p->second, *point_p, d,
389 d_bar, v_bar, m_likelihood.get_noise_sigma(),
390 m_likelihood.get_convert());
393 if(lb_pair_p->first == t)
395 return is_neighbor(*point_p, *lb_pair_p->second, 0,
396 d_bar, v_bar, m_likelihood.get_noise_sigma(),
397 m_likelihood.get_convert());
400 if(lb_pair_p == track.begin())
402 d = lb_pair_p->first - t;
404 return is_neighbor(*point_p, *lb_pair_p->second, d,
405 d_bar, v_bar, m_likelihood.get_noise_sigma(),
406 m_likelihood.get_convert());
409 typename Track::const_iterator next_pair_p = lb_pair_p--;
410 d = next_pair_p->first - t;
411 if(!
is_neighbor(*next_pair_p->second, *point_p, d, d_bar, v_bar,
412 m_likelihood.get_noise_sigma(),
413 m_likelihood.get_convert()))
418 d = t - lb_pair_p->first;
419 if(!
is_neighbor(*point_p, *lb_pair_p->second, d, d_bar, v_bar,
420 m_likelihood.get_noise_sigma(), m_likelihood.get_convert()))
430 template<
class Track,
class Lhood>
431 typename Gibbs_proposer<Track, Lhood>::Location_pair
442 while(cur_var + time_p->size() <= var)
445 cur_var += time_p->size();
449 typename std::set<Element>::const_iterator point_p = time_p->begin();
450 std::advance(point_p, var - cur_var);
452 return std::make_pair(t, point_p);
double accumulate(const Matrix_d< R, C, T > &mat, double init)
Definition: m_matrix_d.impl.h:432
Parent::const_iterator const_iterator
Definition: mcmcda_association.h:122
A class that holds data for the tracking problem.
Definition: mcmcda_data.h:50
void swap(Perspective_camera &cam1, Perspective_camera &cam2)
Swap two cameras.
Definition: perspective_camera.h:599
Parent::iterator iterator
Definition: mcmcda_association.h:121
Gibbs_proposer(const Prior< Track > &prior, const Lhood &likelihood, int neighborhood_size=-1)
Ctor with default-constructible track.
Definition: mcmcda_gibbs_proposer.h:64
Computes the prior of an MCMCDA association.
Definition: mcmcda_prior.h:38
A class that represents a MCMCDA association.
Definition: mcmcda_association.h:66
size_t get_association_dimension(const Association< Track > &w)
Returns the dimension of an association.
Definition: mcmcda_gibbs_proposer.h:123
size_t dims(const Scene &scene, bool respect_changed=true, bool infer_head=true)
Computes the number of variables in this scene.
Definition: pt_scene.cpp:106
Location_pair get_time_and_place(size_t var, const Data< Element > &data) const
Converts a variable index into a time and point;.
Definition: mcmcda_gibbs_proposer.h:433
const Data< Element > & get_data() const
Returns data set const-ref.
Definition: mcmcda_association.h:126
std::pair< std::vector< Track >, std::vector< size_t > > sample(const Prior< Track > &prior, size_t T, const Track &def)
Sample an association from the prior.
Definition: mcmcda_prior.h:193
void swap(kjb::Gsl_Multimin_fdf &m1, kjb::Gsl_Multimin_fdf &m2)
Swap two wrapped multimin objects.
Definition: gsl_multimin.h:693
Parent::const_iterator const_iterator
Definition: mcmcda_data.h:83
Gibbs proposal mechanism for tracking. Complies with Gibbs proposer concept.
Definition: mcmcda_gibbs_proposer.h:46
bool is_neighbor(const Element &y, const Element &y_p, int d, int d_bar, double v_bar, double sg, const typename Data< Element >::Convert &to_vector)
Returns true if given point is a neighbor of second given point.
Definition: mcmcda_data.h:182
get the indices of edges in each direction for i
Definition: APPgetLargeConnectedEdges.m:48
bool is_track_affected(const Track &track, const Location_pair &cur_location) const
Tests whether a track is affected by current variable.
Definition: mcmcda_gibbs_proposer.h:370
Definition for the Vector class, a thin wrapper on the KJB Vector struct and its related functionalit...