19 #ifndef MCMCDA_PROPOSER_H_INCLUDED
20 #define MCMCDA_PROPOSER_H_INCLUDED
35 #include <l/l_sys_lib.h>
46 #include <boost/foreach.hpp>
47 #include <boost/bimap.hpp>
48 #include <boost/bimap/set_of.hpp>
49 #include <boost/bimap/multiset_of.hpp>
73 #warning "[Code police] scoped enums are a c++0x feature."
88 using namespace boost::bimaps;
91 template <
class Track>
96 typedef typename Track::Element Element;
98 typedef typename Track::const_iterator Track_const_iterator;
99 typedef typename Track::const_reverse_iterator Track_const_riterator;
100 typedef boost::function1<Element, const std::vector<const Element*>&> Avg;
101 typedef boost::bimap<multiset_of<double, std::greater<double> >,
102 set_of<const Element*> > Probability_map;
103 typedef boost::function5<std::vector<double>,
const Track*, int,
104 const Element*, int,
size_t> Feature_prob;
108 std::vector<Categorical_distribution<size_t> > m_p_move;
115 Feature_prob m_feature_prob;
116 double m_noise_sigma;
131 double detection_noise_variance,
132 const Track& empty_track
138 m_convert(convert_to_vector),
140 m_noise_sigma(detection_noise_variance),
141 def_track(empty_track),
144 make_p_move(move_distribution);
156 const Feature_prob& feature_prob,
157 double detection_noise_variance,
158 const Track& empty_track
164 m_convert(convert_to_vector),
166 m_feature_prob(feature_prob),
167 m_noise_sigma(detection_noise_variance),
168 def_track(empty_track),
171 make_p_move(move_distribution);
178 double v_bar()
const {
return m_v_bar; }
181 int d_bar()
const {
return m_d_bar; }
184 double gamma()
const {
return m_gamma; }
190 size_t K = w.size() < 2 ? w.size() : 2;
197 size_t K = w.size() < 2 ? w.size() : 2;
198 return log_pdf(m_p_move[K], m);
203 double propose_birth(
const Assoc& w, Assoc& w_p)
const;
206 double propose_death(
const Assoc& w, Assoc& w_p)
const;
209 double propose_split(
const Assoc& w, Assoc& w_p)
const;
212 double propose_merge(
const Assoc& w, Assoc& w_p)
const;
215 double propose_extension(
const Assoc& w, Assoc& w_p)
const;
218 double propose_reduction(
const Assoc& w, Assoc& w_p)
const;
221 double propose_switch(
const Assoc& w, Assoc& w_p)
const;
224 double propose_secretion(
const Assoc& w, Assoc& w_p)
const;
227 double propose_absorption(
const Assoc& w, Assoc& w_p)
const;
231 double p_birth(
const Assoc& w,
const Assoc& w_p)
const;
234 double p_death(
const Assoc& w,
const Assoc& w_p)
const;
237 double p_split(
const Assoc& w,
const Assoc& w_p)
const;
240 double p_merge(
const Assoc& w,
const Assoc& w_p)
const;
243 double p_extension(
const Assoc& w,
const Assoc& w_p)
const;
246 double p_reduction(
const Assoc& w,
const Assoc& w_p)
const;
249 double p_switch(
const Assoc& w,
const Assoc& w_p)
const;
252 double p_secretion(
const Assoc& w,
const Assoc& w_p)
const;
255 double p_absorption(
const Assoc& w,
const Assoc& w_p)
const;
259 bool is_valid_birth(
const Assoc& w,
const Assoc& w_p)
const;
262 bool is_valid_death(
const Assoc& w,
const Assoc& w_p)
const;
265 bool is_valid_split(
const Assoc& w,
const Assoc& w_p)
const;
268 bool is_valid_merge(
const Assoc& w,
const Assoc& w_p)
const;
271 bool is_valid_extension(
const Assoc& w,
const Assoc& w_p)
const;
274 bool is_valid_reduction(
const Assoc& w,
const Assoc& w_p)
const;
277 bool is_valid_switch(
const Assoc& w,
const Assoc& w_p)
const;
280 bool is_valid_secretion(
const Assoc& w,
const Assoc& w_p)
const;
283 bool is_valid_absorption(
const Assoc& w,
const Assoc& w_p)
const;
287 void grow_track_forward(Track& track,
const Assoc& w)
const;
290 void grow_track_backward(Track& track,
const Assoc& w)
const;
293 double p_grow_track_forward
301 double p_grow_track_backward
309 double p_choose_nothing
311 const Probability_map& probs
340 m_p_move.push_back(move_distribution);
345 Probability_map get_growing_probabilities
349 std::set<const Element*>& candidates,
359 double vel_sigma = vel.
empty() ? d*d*m_v_bar*m_v_bar : 0.0;
360 double sg_u = sqrt(m_noise_sigma + vel_sigma);
361 double sg_v = sqrt(m_noise_sigma + vel_sigma);
365 return pdf(N_u, sg_u) *
pdf(N_v, sg_v);
369 std::pair<Vector, Vector> track_future
376 std::pair<Vector, Vector> track_past
383 std::pair<Vector, Vector> track_future_ls
390 std::pair<Vector, Vector> track_past_ls
397 std::pair<Vector, Vector> track_future_tls
404 std::pair<Vector, Vector> track_past_tls
412 size_t count_split_points(
const Assoc& w)
const;
415 size_t count_merge_pairs(
const Assoc& w)
const;
418 size_t count_secretion_tracks(
const Assoc& w)
const;
421 size_t count_absorption_pairs(
const Assoc& w)
const;
424 double negative_infinity()
const
430 static const int GROW_FORWARD;
431 static const int GROW_BACKWARD;
437 Assoc_const_iterator new_track_p;
449 size_t num_split_points;
455 size_t num_merge_pairs;
459 struct Extension_info
462 Assoc_const_iterator extended_track_p;
468 struct Reduction_info
471 size_t reduced_track_size;
477 size_t num_switch_points;
481 struct Secretion_info
483 size_t num_valid_tracks;
488 struct Absorption_info
490 size_t num_valid_track_pairs;
494 mutable Birth_info birth_info;
495 mutable Death_info death_info;
496 mutable Split_info split_info;
497 mutable Merge_info merge_info;
498 mutable Extension_info extension_info;
499 mutable Reduction_info reduction_info;
500 mutable Switch_info switch_info;
501 mutable Secretion_info secretion_info;
502 mutable Absorption_info absorption_info;
507 template <
class Track>
508 const int Proposer<Track>::GROW_FORWARD = 1;
510 template <
class Track>
511 const int Proposer<Track>::GROW_BACKWARD = 2;
515 template <
class Track>
523 bool keep_going =
true;
524 std::string move_name =
"association-";
529 size_t m = sample_move(w);
535 fwd = propose_birth(w, w_p);
536 if(fwd == negative_infinity())
541 rev = p_death(w_p, w);
542 if(rev == negative_infinity())
550 move_name +=
"birth";
556 fwd = propose_death(w, w_p);
557 if(fwd == negative_infinity())
562 rev = p_birth(w_p, w);
563 if(rev == negative_infinity())
571 move_name +=
"death";
619 fwd = propose_extension(w, w_p);
620 if(fwd == negative_infinity())
625 rev = p_reduction(w_p, w);
626 if(rev == negative_infinity())
634 move_name +=
"extension";
640 fwd = propose_reduction(w, w_p);
641 if(fwd == negative_infinity())
646 rev = p_extension(w_p, w);
647 if(rev == negative_infinity())
655 move_name +=
"reduction";
661 fwd = propose_switch(w, w_p);
662 if(fwd == negative_infinity())
672 move_name +=
"switch";
678 fwd = propose_secretion(w, w_p);
679 if(fwd == negative_infinity())
684 rev = p_absorption(w_p, w);
685 if(rev == negative_infinity())
694 move_name +=
"secretion";
700 fwd = propose_absorption(w, w_p);
701 if(fwd == negative_infinity())
706 rev = p_secretion(w_p, w);
707 if(rev == negative_infinity())
716 move_name +=
"absorption";
723 "MCMCDA: move %d does not exist.", (m));
728 if(fwd != negative_infinity() && rev != negative_infinity())
739 template <
class Track>
744 Track new_track = def_track;
754 return negative_infinity();
758 std::vector<double> ps(L_1.size(), 1.0);
761 typename std::set<const Element*>::const_iterator x_1_p = L_1.begin();
763 std::advance(x_1_p, n);
764 new_track.insert(std::make_pair(t_1, *x_1_p));
767 grow_track_forward(new_track, w);
768 if(new_track.real_size() == 1)
770 return negative_infinity();
775 Assoc_const_iterator new_track_p = w_p.insert(new_track).first;
778 new_track_p->set_changed_start(new_track_p->get_start_time());
779 new_track_p->set_changed_end(new_track_p->get_end_time());
782 birth_info.new_track_p = new_track_p;
783 death_info.num_tracks = w_p.size();
784 return p_birth(w, w_p);
789 template <
class Track>
796 std::cout <<
"WARNING: cannot propose death for empty association.\n";
797 return std::log(0.0);
802 Assoc_const_iterator new_track_p =
element_uar(w.begin(), w.size());
803 w_p.erase(*new_track_p);
806 death_info.num_tracks = w.size();
807 birth_info.new_track_p = new_track_p;
808 return p_death(w, w_p);
813 template <
class Track>
820 std::cout <<
"WARNING: cannot propose split for empty association.\n";
821 return std::log(0.0);
824 size_t nsp = count_split_points(w);
827 return negative_infinity();
834 Assoc_const_iterator track_p;
835 Track_const_iterator rp;
837 for(track_p = w_p.begin(); track_p != w_p.end(); track_p++)
839 size_t tsz = track_p->real_size();
842 if(count + tsz - 3 >= n)
846 Track_const_iterator pair_p;
847 for(pair_p = track_p->begin(); pair_p != track_p->end();
850 if(pair_p->first == t)
857 if(r == n - count + 1)
871 Track beg = def_track;
872 Track fin = def_track;
873 beg.insert(track_p->begin(), rp);
874 fin.insert(rp, track_p->end());
877 beg.set_changed_start(beg.get_start_time());
878 beg.set_changed_end(beg.get_end_time());
879 fin.set_changed_start(fin.get_start_time());
880 fin.set_changed_end(fin.get_end_time());
888 split_info.num_split_points = nsp;
889 merge_info.num_merge_pairs = count_merge_pairs(w_p);
890 return p_split(w, w_p);
895 template <
class Track>
902 std::cout <<
"WARNING: cannot propose merge for single-track "
904 return std::log(0.0);
910 std::list<std::pair<Assoc_const_iterator, Assoc_const_iterator> > M;
911 for(Assoc_const_iterator track_p1 = w_p.begin();
912 track_p1 != w_p.end();
915 int t_f = track_p1->get_end_time();
916 for(Assoc_const_iterator track_p2 = w_p.begin();
917 track_p2 != w_p.end();
920 int t_1 = track_p2->get_start_time();
921 if(track_p1 != track_p2 && t_f < t_1)
924 track_p2->get_start_point(),
925 t_1 - t_f, m_d_bar, m_v_bar,
926 m_noise_sigma, m_convert))
928 M.push_back(std::make_pair(track_p1, track_p2));
936 return negative_infinity();
940 std::pair<Assoc_const_iterator, Assoc_const_iterator> ttm
942 Track new_track = *ttm.first;
943 new_track.insert(ttm.second->begin(), ttm.second->end());
946 new_track.set_changed_start(ttm.second->get_start_time());
947 new_track.set_changed_end(ttm.second->get_end_time());
950 w_p.erase(ttm.first);
951 w_p.erase(ttm.second);
952 w_p.insert(new_track);
955 merge_info.num_merge_pairs = M.size();
956 split_info.num_split_points = count_split_points(w_p);
957 return p_merge(w, w_p);
962 template <
class Track>
973 std::cout <<
"WARNING: cannot propose extension for empty "
975 return std::log(0.0);
981 Assoc_const_iterator track_p =
element_uar(w_p.begin(), w_p.size());
982 Track extended_track = *track_p;
983 size_t original_track_size = extended_track.size();
992 prev_end = extended_track.get_end_time();
993 grow_track_forward(extended_track, w_p);
994 extended_track.set_changed_start(prev_end);
995 extended_track.set_changed_end(extended_track.get_end_time());
1000 prev_end = extended_track.get_start_time();
1001 grow_track_backward(extended_track, w_p);
1002 extended_track.set_changed_start(extended_track.get_start_time());
1003 extended_track.set_changed_end(prev_end);
1008 if(extended_track.size() == original_track_size)
1010 return negative_infinity();
1015 Assoc_const_iterator extended_track_p = w_p.insert(extended_track).first;
1018 extension_info.num_tracks = w.size();
1019 extension_info.extended_track_p = extended_track_p;
1020 extension_info.direction = dir;
1021 extension_info.previous_end = prev_end;
1022 reduction_info.num_tracks = w_p.size();
1023 reduction_info.reduced_track_size = extended_track.real_size();
1024 return p_extension(w, w_p);
1029 template <
class Track>
1040 std::cout <<
"WARNING: cannot propose reduction for empty "
1042 return std::log(0.0);
1051 if(track_p->real_size() <= 2)
1053 return negative_infinity();
1056 Track reduced_track = *track_p;
1057 size_t rtsz = reduced_track.real_size();
1059 int t = reduced_track.get_nth_time(n);
1065 reduced_track.erase(reduced_track.upper_bound(t), reduced_track.end());
1070 reduced_track.erase(reduced_track.begin(), reduced_track.lower_bound(t));
1071 dir = GROW_BACKWARD;
1075 reduced_track.set_changed_start(-1);
1076 reduced_track.set_changed_end(-1);
1079 Assoc_const_iterator orig_track_p = w.find(*track_p);
1081 w_p.insert(reduced_track);
1084 reduction_info.num_tracks = w.size();
1085 reduction_info.reduced_track_size = rtsz;
1086 extension_info.num_tracks = w_p.size();
1087 extension_info.extended_track_p = orig_track_p;
1088 extension_info.direction = dir;
1089 extension_info.previous_end = t;
1090 return p_reduction(w, w_p);
1095 template <
class Track>
1104 using namespace boost;
1106 typedef tuple<const Track&, int, int, const Track&, int, int> Switch_point;
1108 std::list<Switch_point> switch_points;
1109 for(Assoc_const_iterator track_p1 = w.begin();
1110 track_p1 != w.end();
1113 Assoc_const_iterator track_p2 = track_p1;
1114 for(track_p2++; track_p2 != w.end(); track_p2++)
1117 Track_const_iterator last_p1
1118 = track_p1->lower_bound(track_p1->rbegin()->first);
1119 for(Track_const_iterator pair_p1 = track_p1->begin();
1124 int curt = pair_p1->first;
1125 while(pair_p1->first == curt)
1132 Track_const_iterator last_p2
1133 = track_p2->lower_bound(track_p2->rbegin()->first);
1134 for(Track_const_iterator pair_p2 = track_p2->begin();
1139 curt = pair_p2->first;
1140 while(pair_p2->first == curt)
1146 Track_const_iterator pair_q1 = pair_p1;
1147 Track_const_iterator pair_q2 = pair_p2;
1149 int t1 = pair_p1->first;
1150 int t2 = pair_p2->first;
1151 int tp1 = (++pair_q1)->first;
1152 int tq1 = (++pair_q2)->first;
1153 if(tp1 > t2 && tq1 > t1
1155 *pair_q2->second, tq1 - t1, m_d_bar,
1156 m_v_bar, m_noise_sigma, m_convert)
1158 *pair_q1->second, tp1 - t2, m_d_bar,
1159 m_v_bar, m_noise_sigma, m_convert))
1161 switch_points.push_back(
1162 Switch_point(*track_p1, t1, tp1,
1163 *track_p2, t2, tq1));
1170 if(switch_points.empty())
1172 return negative_infinity();
1175 Switch_point sp = *
element_uar(switch_points.begin(), switch_points.size());
1176 Track track1_new = get<0>(sp);
1177 Track track2_new = get<3>(sp);
1179 get<2>(sp), get<4>(sp), get<5>(sp));
1182 track1_new.set_changed_start(get<1>(sp));
1183 track1_new.set_changed_end(track1_new.get_end_time());
1184 track2_new.set_changed_start(get<4>(sp));
1185 track2_new.set_changed_end(track2_new.get_end_time());
1188 w_p.erase(get<0>(sp));
1189 w_p.erase(get<3>(sp));
1190 w_p.insert(track1_new);
1191 w_p.insert(track2_new);
1194 switch_info.num_switch_points = switch_points.size();
1195 return p_switch(w, w_p);
1200 template <
class Track>
1211 std::cout <<
"WARNING: cannot propose secretion for empty "
1213 return std::log(0.0);
1221 std::vector<Assoc_iterator> long_tracks;
1222 for(Assoc_iterator track_p = w_p.begin(); track_p != w_p.end(); track_p++)
1224 if(track_p->real_size() > 4)
1226 long_tracks.push_back(track_p);
1230 if(track_p->count(track_p->get_start_time()) > 1
1231 && track_p->count(track_p->get_end_time()) > 1)
1233 long_tracks.push_back(track_p);
1238 if(long_tracks.empty())
1240 return negative_infinity();
1244 Assoc_iterator orig_p = *
element_uar(long_tracks.begin(),
1245 long_tracks.size());
1246 Track orig_track = def_track;
1247 Track new_track = def_track;
1250 size_t sz = orig_p->size();
1251 Track_const_iterator pair_p1 =
element_uar(orig_p->begin(), sz);
1252 Track_const_iterator pair_p2 =
element_uar(orig_p->begin(), sz);
1254 Track_const_iterator small_p, big_p;
1255 if(pair_p1->first < pair_p2->first)
1260 else if(pair_p1->first > pair_p2->first)
1267 if(std::distance(orig_p->begin(), pair_p1)
1268 < std::distance(orig_p->begin(), pair_p2))
1281 orig_track.insert(orig_p->begin(), ++small_p);
1286 orig_track.insert(++big_p, orig_p->end());
1290 new_track.insert(++big_p, orig_p->end());
1295 while(small_p != big_p)
1299 orig_track.insert(*small_p);
1303 new_track.insert(*small_p);
1311 if(!orig_track.is_valid(m_v_bar, m_d_bar, m_noise_sigma, m_convert)
1312 || !new_track.is_valid(m_v_bar, m_d_bar, m_noise_sigma, m_convert))
1314 return negative_infinity();
1318 orig_track.set_changed_start(orig_track.get_start_time());
1319 orig_track.set_changed_end(orig_track.get_end_time());
1320 new_track.set_changed_start(new_track.get_start_time());
1321 new_track.set_changed_end(new_track.get_end_time());
1325 w_p.insert(orig_track);
1326 w_p.insert(new_track);
1329 secretion_info.num_valid_tracks = long_tracks.size();
1330 secretion_info.window_size = win_sz;
1331 absorption_info.num_valid_track_pairs = count_absorption_pairs(w_p);
1332 return p_secretion(w, w_p);
1337 template <
class Track>
1348 std::cout <<
"WARNING: cannot propose absorption for single-track "
1350 return std::log(0.0);
1356 using namespace boost;
1357 typedef tuple<Assoc_const_iterator, Assoc_const_iterator, Track> iit_tuple;
1359 std::vector<double> ps;
1360 ps.reserve((w.size()*(w.size()-1))/2);
1361 std::list<iit_tuple> M;
1363 for(Assoc_const_iterator track_p1 = w_p.begin();
1364 track_p1 != w_p.end();
1367 Assoc_const_iterator track_p2 = track_p1;
1368 for(++track_p2; track_p2 != w_p.end(); ++track_p2)
1371 Track new_track = def_track;
1373 size_t sf1 = track_p1->get_start_time();
1374 size_t ef1 = track_p1->get_end_time();
1375 size_t sf2 = track_p2->get_start_time();
1376 size_t ef2 = track_p2->get_end_time();
1378 int l1 = ef1 - sf1 + 1;
1379 int l2 = ef2 - sf2 + 1;
1381 bool olap = (sf1 <= ef2 && sf2 <= ef1);
1385 new_track = *track_p1;
1386 new_track.insert(track_p2->begin(), track_p2->end());
1408 new_track = *track_p2;
1409 new_track.insert(track_p1->begin(), track_p1->end());
1430 new_track.set_changed_start(chst);
1431 new_track.set_changed_end(ched);
1433 double pg = p_grow_track_forward(new_track, w_e,
1434 new_track.get_start_time());
1435 if(pg != negative_infinity())
1437 M.push_back(boost::make_tuple(track_p1, track_p2, new_track));
1445 return negative_infinity();
1449 typename std::list<iit_tuple>::const_iterator iit_p = M.begin();
1450 for(
size_t i = 0;
i < ps.size();
i++, iit_p++)
1452 ps[
i] /= get<2>(*iit_p).real_size();
1454 ps[
i] = std::sqrt(ps[
i]);
1458 double lsum =
log_sum(ps.begin(), ps.end());
1459 std::transform(ps.begin(), ps.end(), ps.begin(),
1460 std::bind2nd(std::minus<double>(), lsum));
1461 std::transform(ps.begin(), ps.end(), ps.begin(),
1462 std::ptr_fun<double, double>(std::exp));
1466 typename std::list<iit_tuple>::iterator tuple_p = M.begin();
1467 std::advance(tuple_p, idx - 1);
1468 iit_tuple& ttm = *tuple_p;
1471 const Element* st_elem_p = get<0>(ttm)->begin()->second;
1472 if(get<0>(ttm)->begin()->second == get<2>(ttm).begin()->second)
1474 st_elem_p = get<1>(ttm)->begin()->second;
1477 const Element* en_elem_p = get<0>(ttm)->rbegin()->second;
1478 if(get<0>(ttm)->rbegin()->second == get<2>(ttm).rbegin()->second)
1480 en_elem_p = get<1>(ttm)->rbegin()->second;
1487 BOOST_FOREACH(
const typename Track::value_type& pr, get<2>(ttm))
1489 if(pr.second == st_elem_p)
1494 if(pr.second == en_elem_p)
1503 w_p.erase(get<0>(ttm));
1504 w_p.erase(get<1>(ttm));
1505 w_p.insert(get<2>(ttm));
1508 absorption_info.num_valid_track_pairs = M.size();
1509 secretion_info.num_valid_tracks = count_secretion_tracks(w_p);
1510 secretion_info.window_size =
std::abs(last_start - first_end + 1);
1511 return p_absorption(w, w_p);
1516 template <
class Track>
1523 if(!is_valid_birth(w, w_p))
1525 std::cout <<
"WARNING: invalid birth." << std::endl;
1526 return std::log(0.0);
1530 const Track& extra_track = *birth_info.new_track_p;
1534 Track_const_iterator pair_p = extra_track.begin();
1535 int t_1 = pair_p->first;
1536 p -= std::log(w.
get_data().size());
1543 return negative_infinity();
1546 p -= std::log(num_valid);
1549 p += p_grow_track_forward(extra_track, w, t_1);
1556 template <
class Track>
1563 if(!is_valid_death(w, w_p))
1565 std::cout <<
"WARNING: invalid death." << std::endl;
1566 return std::log(0.0);
1570 return -std::log(death_info.num_tracks);
1575 template <
class Track>
1582 if(!is_valid_split(w, w_p))
1584 std::cout <<
"WARNING: invalid split." << std::endl;
1585 return std::log(0.0);
1589 size_t nsp = split_info.num_split_points;
1592 return negative_infinity();
1595 return -std::log(nsp);
1600 template <
class Track>
1607 if(!is_valid_merge(w, w_p))
1609 std::cout <<
"WARNING: invalid merge." << std::endl;
1610 return std::log(0.0);
1614 size_t nmp = merge_info.num_merge_pairs;
1617 return negative_infinity();
1620 return -std::log(nmp);
1625 template <
class Track>
1636 if(!is_valid_extension(w, w_p))
1638 std::cout <<
"WARNING: invalid extension." << std::endl;
1639 return std::log(0.0);
1643 double p = -std::log(2 * extension_info.num_tracks);
1646 if(extension_info.direction == GROW_FORWARD)
1648 p += p_grow_track_forward(*extension_info.extended_track_p, w,
1649 extension_info.previous_end);
1653 p += p_grow_track_backward(*extension_info.extended_track_p, w,
1654 extension_info.previous_end);
1662 template <
class Track>
1673 if(!is_valid_reduction(w, w_p))
1675 std::cout <<
"WARNING: invalid reduction." << std::endl;
1676 return std::log(0.0);
1680 return std::log(0.5 / (reduction_info.num_tracks
1681 * (reduction_info.reduced_track_size - 2)));
1686 template <
class Track>
1693 if(!is_valid_switch(w, w_p))
1695 std::cout <<
"WARNING: invalid switch." << std::endl;
1696 return std::log(0.0);
1700 size_t nsp = switch_info.num_switch_points;
1703 return negative_infinity();
1706 return -std::log(nsp);
1711 template <
class Track>
1722 if(!is_valid_secretion(w, w_p))
1724 std::cout <<
"WARNING: invalid secretion." << std::endl;
1725 return std::log(0.0);
1729 if(secretion_info.num_valid_tracks == 0)
1731 return negative_infinity();
1734 double p = -std::log(secretion_info.num_valid_tracks);
1736 p += (secretion_info.window_size)*std::log(0.5);
1743 template <
class Track>
1754 if(!is_valid_absorption(w, w_p))
1756 std::cout <<
"WARNING: invalid absorption." << std::endl;
1757 return std::log(0.0);
1761 return -std::log(absorption_info.num_valid_track_pairs);
1766 template <
class Track>
1774 if(w_p.empty() || w_p.size() - w.size() != 1)
1779 std::pair<Assoc_const_iterator, Assoc_const_iterator> diff_p;
1780 diff_p = std::mismatch(w.begin(), w.end(), w_p.begin());
1781 if(diff_p.first == w.end())
1787 diff_p = std::mismatch(diff_p.first, w.end(), diff_p.second);
1788 if(diff_p.first == w.end())
1798 template <
class Track>
1800 bool Proposer<Track>::is_valid_death
1806 return is_valid_birth(w_p, w);
1811 template <
class Track>
1813 bool Proposer<Track>::is_valid_split
1819 if(w.empty() || w_p.size() - w.size() != 1)
1824 std::set<Track, Compare_tracks<Track> > w_minus_w_p;
1825 std::set_difference(w.begin(), w.end(), w_p.begin(), w_p.end(),
1826 std::inserter(w_minus_w_p, w_minus_w_p.begin()),
1827 Compare_tracks<Track>());
1829 std::set<Track, Compare_tracks<Track> > w_p_minus_w;
1830 std::set_difference(w_p.begin(), w_p.end(), w.begin(), w.end(),
1831 std::inserter(w_p_minus_w, w_p_minus_w.begin()),
1832 Compare_tracks<Track>());
1834 if(w_minus_w_p.size() != 1 || w_p_minus_w.size() != 2)
1839 Track merged_track = *w_p_minus_w.begin();
1840 merged_track.insert((++w_p_minus_w.begin())->begin(),
1841 (++w_p_minus_w.begin())->end());
1843 if(*w_minus_w_p.begin() != merged_track)
1853 template <
class Track>
1855 bool Proposer<Track>::is_valid_merge
1861 return is_valid_split(w_p, w);
1866 template <
class Track>
1867 bool Proposer<Track>::is_valid_extension
1873 if(w.empty() || w.size() != w_p.size())
1879 std::set<Track, Compare_tracks<Track> > diff;
1880 std::set_difference(w.begin(), w.end(), w_p.begin(), w_p.end(),
1881 std::inserter(diff, diff.begin()), Compare_tracks<Track>());
1882 if(diff.size() != 1)
1887 const Track& original_track = *w.find(*diff.begin());
1891 std::set_difference(w_p.begin(), w_p.end(), w.begin(), w.end(),
1892 std::inserter(diff, diff.begin()), Compare_tracks<Track>());
1894 if(diff.size() != 1)
1899 Assoc_const_iterator extended_track_p = w_p.find(*diff.begin());
1901 std::pair<Track_const_iterator, Track_const_iterator> diff_ps
1902 = std::mismatch(original_track.begin(), original_track.end(),
1903 extended_track_p->begin());
1904 if(diff_ps.first == original_track.end())
1909 std::pair<
typename Track::const_reverse_iterator,
1910 typename Track::const_reverse_iterator>
1911 r_diff_ps = std::mismatch(original_track.rbegin(), original_track.rend(),
1912 extended_track_p->rbegin());
1913 if(r_diff_ps.first == original_track.rend())
1923 template <
class Track>
1925 bool Proposer<Track>::is_valid_reduction
1931 return is_valid_extension(w_p, w);
1936 template <
class Track>
1937 bool Proposer<Track>::is_valid_switch
1943 if(w.size() < 2 || w.size() != w_p.size())
1949 std::set<Track, Compare_tracks<Track> > diff;
1950 std::set_difference(w.begin(), w.end(), w_p.begin(), w_p.end(),
1951 std::inserter(diff, diff.begin()), Compare_tracks<Track>());
1952 if(diff.size() != 2)
1957 const Track& original_track_1 = *w.find(*diff.begin());
1958 const Track& original_track_2 = *w.find(*diff.rbegin());
1962 std::set_difference(w_p.begin(), w_p.end(), w.begin(), w.end(),
1963 std::inserter(diff, diff.begin()), Compare_tracks<Track>());
1965 if(diff.size() != 2)
1970 const Track& swapped_track_1 = *w_p.find(*diff.begin());
1971 const Track& swapped_track_2 = *w_p.find(*diff.rbegin());
1973 std::pair<Track_const_iterator, Track_const_iterator> track1_break
1974 = std::mismatch(original_track_1.begin(), original_track_1.end(),
1975 swapped_track_1.begin());
1977 std::pair<Track_const_iterator, Track_const_iterator> track2_break
1978 = std::mismatch(original_track_2.begin(), original_track_2.end(),
1979 swapped_track_2.begin());
1981 if(track1_break.first == original_track_1.end())
1986 if(!std::equal(original_track_1.begin(), track1_break.first,
1987 swapped_track_1.begin()) ||
1988 !std::equal(track1_break.first, original_track_1.end(),
1989 track2_break.second) ||
1990 !std::equal(original_track_2.begin(), track2_break.first,
1991 swapped_track_2.begin()) ||
1992 !std::equal(track2_break.first, original_track_2.end(),
1993 track1_break.second))
2003 template <
class Track>
2005 bool Proposer<Track>::is_valid_secretion
2011 if(w.empty() || w_p.size() - w.size() != 1)
2016 std::set<Track, Compare_tracks<Track> > w_minus_w_p;
2017 std::set_difference(w.begin(), w.end(), w_p.begin(), w_p.end(),
2018 std::inserter(w_minus_w_p, w_minus_w_p.begin()),
2019 Compare_tracks<Track>());
2021 std::set<Track, Compare_tracks<Track> > w_p_minus_w;
2022 std::set_difference(w_p.begin(), w_p.end(), w.begin(), w.end(),
2023 std::inserter(w_p_minus_w, w_p_minus_w.begin()),
2024 Compare_tracks<Track>());
2026 if(w_minus_w_p.size() != 1 || w_p_minus_w.size() != 2)
2055 template <
class Track>
2057 bool Proposer<Track>::is_valid_absorption
2063 return is_valid_secretion(w_p, w);
2068 template <
class Track>
2075 int t = track.get_end_time() + 1;
2077 if(t >= w.
get_data().size())
return;
2081 int prev_t = track.get_end_time();
2083 if(d > m_b_bar)
break;
2085 std::pair<Vector, Vector> fut = track_future_ls(track);
2087 const Vector& vel = fut.second;
2093 Probability_map probs =
2094 get_growing_probabilities(track, prev_t, dpt, x, vel, d);
2096 size_t num_added = 0;
2097 BOOST_FOREACH(
const typename Probability_map::left_map::value_type& pr,
2100 double p_det = std::exp(pr.first);
2105 track.insert(std::make_pair(t, pr.second));
2118 if(t == w.
get_data().size())
break;
2125 template <
class Track>
2132 int t = track.get_start_time() - 1;
2138 int prev_t = track.get_start_time();
2142 std::pair<Vector, Vector> fut = track_past_ls(track);
2144 const Vector& vel = fut.second;
2150 Probability_map probs =
2151 get_growing_probabilities(track, prev_t, dpt, x, vel, d);
2153 size_t num_added = 0;
2154 BOOST_FOREACH(
const typename Probability_map::left_map::value_type& pr,
2157 double p_det = std::exp(pr.first);
2161 track.insert(std::make_pair(t, pr.second));
2180 template <
class Track>
2197 Track_const_iterator eot = --track.end();
2200 return negative_infinity();
2203 Track_const_iterator prev_p = track.upper_bound(t);
2204 Track_const_iterator next_p = prev_p--;
2207 while(prev_p != eot)
2209 int prev_t = prev_p->first;
2211 if(t - prev_t > m_d_bar)
2213 return negative_infinity();
2216 std::pair<Vector, Vector> fut = track_future_ls(track, prev_t);
2217 const Vector& prev_x = fut.first;
2218 const Vector& vel = fut.second;
2221 Probability_map probs = get_growing_probabilities(track, prev_t,
2222 dpt, prev_x, vel, t - prev_t);
2225 p += p_choose_nothing(probs);
2226 if(t == next_p->first)
2229 while(next_p != track.end() && next_p->first == t)
2231 const Element&
x = *next_p->second;
2232 typename Probability_map::right_map::const_iterator
2233 pr_p = probs.right.find(&x);
2234 if(pr_p == probs.right.end())
2236 std::cout <<
"WARNING: this should never happen!\n";
2242 if(p != negative_infinity())
2244 p += (pr_p->second) - log(1.0 - exp(pr_p->second));
2255 p += std::log(1.0 - m_gamma);
2268 double p_last = 0.0;
2270 int pt = prev_p->first;
2271 std::pair<Vector, Vector> fut = track_future_ls(track, pt);
2272 const Vector& px = fut.first;
2273 const Vector& vel = fut.second;
2277 Probability_map probs =
2278 get_growing_probabilities(track, pt, dpt, px, vel, 1);
2280 double a_t_1 = std::exp(p_choose_nothing(probs)) * m_gamma;
2285 for(; t <= last_t; t++)
2290 probs = get_growing_probabilities(track, pt, dpt, px, vel, t - pt);
2291 double a_t = a_t_1 * (1 - m_gamma)
2292 * std::exp(p_choose_nothing(probs));
2297 p += std::log(p_last);
2304 template <
class Track>
2317 return std::log(0.0);
2321 Track_const_riterator eot = --track.rend();
2324 return negative_infinity();
2327 Track_const_riterator prev_p(track.lower_bound(t));
2328 Track_const_riterator next_p = prev_p--;
2331 while(prev_p != eot)
2333 int prev_t = prev_p->first;
2335 std::pair<Vector, Vector> fut = track_past_ls(track, prev_t);
2336 const Vector& prev_x = fut.first;
2337 const Vector& vel = fut.second;
2340 Probability_map probs = get_growing_probabilities(track, prev_t, dpt,
2341 prev_x, vel, t - prev_t);
2344 p += p_choose_nothing(probs);
2345 if(t == next_p->first)
2348 while(next_p != track.rend() && next_p->first == t)
2350 const Element&
x = *next_p->second;
2351 typename Probability_map::right_map::const_iterator
2352 pr_p = probs.right.find(&x);
2353 if(pr_p == probs.right.end())
2355 std::cout <<
"WARNING: this should never happen!\n";
2361 if(p != negative_infinity())
2363 p += (pr_p->second) - log(1.0 - exp(pr_p->second));
2374 p += std::log(1.0 - m_gamma);
2387 double p_last = 0.0;
2389 int pt = prev_p->first;
2390 std::pair<Vector, Vector> fut = track_past_ls(track, pt);
2391 const Vector& px = fut.first;
2392 const Vector& vel = fut.second;
2396 Probability_map probs =
2397 get_growing_probabilities(track, pt, dpt, px, vel, -1);
2399 double a_t_1 = std::exp(p_choose_nothing(probs)) * m_gamma;
2403 int last_t =
std::max(pt - m_d_bar, 1);
2404 for(; t >= last_t; t--)
2410 get_growing_probabilities(track, pt, dpt, px, vel, t - pt);
2411 double a_t = a_t_1 * (1 - m_gamma)
2412 * std::exp(p_choose_nothing(probs));
2417 p += std::log(p_last);
2424 template <
class Track>
2427 const Probability_map& probabilities
2431 BOOST_FOREACH(
const typename Probability_map::left_map::value_type& pr,
2434 p += std::log(1.0 - std::exp(pr.first));
2442 template <
class Track>
2443 typename Proposer<Track>::Probability_map
2448 std::set<const Element*>& candidates,
2454 Probability_map probabilities;
2456 double vel_sigma = vel.
empty() ? dist*dist*m_v_bar*m_v_bar : 0.0;
2457 double sg_u = sqrt(m_noise_sigma + vel_sigma);
2458 double sg_v = sqrt(m_noise_sigma + vel_sigma);
2463 if(!vel.
empty()) v_x += dist*vel;
2465 double noise_score = get_velocity_noise_score(d, vel);
2467 BOOST_FOREACH(
const Element* e_p, candidates)
2469 Vector v = m_convert(*e_p);
2470 double score =
log_pdf(N_u, v_x[0]-v[0]) +
log_pdf(N_v, v_x[1]-v[1]);
2471 double p = exp(score) / (exp(score) + noise_score);
2476 std::vector<double> f_ps
2477 = m_feature_prob(&track, t, e_p, t + d, m_b_bar);
2481 BOOST_FOREACH(
double f_p, f_ps)
2486 p = std::pow(p * pp, 1.0/(1 + f_ps.size()));
2489 typedef typename Probability_map::left_map::value_type value_type;
2493 prob = negative_infinity();
2500 probabilities.left.insert(value_type(prob, e_p));
2503 return probabilities;
2508 template <
class Track>
2516 const size_t wsz = m_b_bar / 4;
2517 const double vel_thresh = 0.0;
2519 Track_const_riterator pair_p;
2522 pair_p = track.rbegin();
2529 pair_p = Track_const_riterator(track.upper_bound(t));
2531 "The frame passed to track_future() must be valid.");
2535 std::vector<const Element*> cur_elems;
2536 while(pair_p != track.rend() && pair_p->first == t)
2538 cur_elems.push_back(pair_p->second);
2542 Element x = m_average(cur_elems);
2543 if(pair_p == track.rend())
2545 return std::make_pair(m_convert(x),
Vector());
2549 size_t wnd = t - pair_p->first;
2553 if(pair_p == track.rend())
break;
2554 wnd = t - pair_p->first;
2560 std::vector<const Element*> cur_elems;
2561 while(pair_p != track.rend() && pair_p->first == t - wnd)
2567 Element from = m_average(cur_elems);
2569 vel = m_convert(x) - m_convert(from);
2578 return std::make_pair(m_convert(x), vel);
2583 template <
class Track>
2591 const size_t wsz = m_b_bar / 4;
2592 const double vel_thresh = 0.0;
2594 Track_const_iterator pair_p;
2597 pair_p = track.begin();
2604 pair_p = track.lower_bound(t);
2606 "The frame passed to track_past() must be valid.");
2610 std::vector<const Element*> cur_elems;
2611 while(pair_p != track.end() && pair_p->first == t)
2613 cur_elems.push_back(pair_p->second);
2617 Element x = m_average(cur_elems);
2618 if(pair_p == track.end())
2620 return std::make_pair(m_convert(x),
Vector());
2624 size_t wnd = pair_p->first - t;
2628 if(pair_p == track.end())
break;
2629 wnd = pair_p->first - t;
2635 std::vector<const Element*> cur_elems;
2636 while(pair_p != track.end() && pair_p->first == t + wnd)
2642 Element from = m_average(cur_elems);
2644 vel = m_convert(x) - m_convert(from);
2653 return std::make_pair(m_convert(x), vel);
2658 template <
class Track>
2666 const size_t wsz = m_b_bar;
2667 const double vel_thresh = 0.0;
2669 Track_const_riterator pair_p;
2672 pair_p = track.rbegin();
2679 pair_p = Track_const_riterator(track.upper_bound(t));
2681 "The frame passed to track_future() must be valid.");
2684 double x_N = m_convert(*pair_p->second)[0];
2685 double y_N = m_convert(*pair_p->second)[1];
2688 double sum_x_squared = x_N*x_N;
2689 double sum_xy = x_N*y_N;
2695 std::vector<const Element*> cur_elems;
2696 size_t wnd = t - pair_p->first;
2699 if(pair_p->first == t)
2701 cur_elems.push_back(pair_p->second);
2704 if(pair_p == track.rend())
break;
2705 double x_i = m_convert(*pair_p->second)[0];
2706 double y_i = m_convert(*pair_p->second)[1];
2707 wnd = t - pair_p->first;
2710 sum_x_squared += (x_i * x_i);
2711 sum_xy += (x_i * y_i);
2722 Element e_N = m_average(cur_elems);
2723 Vector p_N = m_convert(e_N);
2724 double den = (sum_x_squared - (sum_x*sum_x)/N);
2727 double m = (sum_xy - (sum_x*sum_y)/N) / den;
2728 double b = sum_y / N - m * sum_x / N;
2730 double f_1 = m * x_1 + b;
2731 double f_N = m * p_N[0] + b;
2734 vel.
set(p_N[0] - x_1, f_N - f_1);
2744 vel.
set(0.0, p_N[1] - y_1);
2747 return std::make_pair(x, vel);
2752 return track_future(track, t);
2757 template <
class Track>
2765 const size_t wsz = m_b_bar;
2766 const double vel_thresh = 0.0;
2768 Track_const_iterator pair_p;
2771 pair_p = track.begin();
2778 pair_p = track.lower_bound(t);
2780 "The frame passed to track_future() must be valid.");
2783 double x_N = m_convert(*pair_p->second)[0];
2784 double y_N = m_convert(*pair_p->second)[1];
2787 double sum_x_squared = x_N*x_N;
2788 double sum_xy = x_N*y_N;
2794 std::vector<const Element*> cur_elems;
2796 size_t wnd = pair_p->first - t;
2799 if(pair_p->first == t)
2801 cur_elems.push_back(pair_p->second);
2804 if(pair_p == track.end())
break;
2805 double x_i = m_convert(*pair_p->second)[0];
2806 double y_i = m_convert(*pair_p->second)[1];
2807 wnd = pair_p->first - t;
2810 sum_x_squared += (x_i * x_i);
2811 sum_xy += (x_i * y_i);
2822 Element e_N = m_average(cur_elems);
2823 Vector p_N = m_convert(e_N);
2824 double den = (sum_x_squared - (sum_x*sum_x)/N);
2827 double m = (sum_xy - (sum_x*sum_y)/N) / den;
2828 double b = sum_y / N - m * sum_x / N;
2830 double f_1 = m * x_1 + b;
2831 double f_N = m * p_N[0] + b;
2834 vel.
set(p_N[0] - x_1, f_N - f_1);
2844 vel.
set(0.0, p_N[1] - y_1);
2847 return std::make_pair(x, vel);
2852 return track_future(track, t);
2857 template <
class Track>
2865 const size_t wsz = m_b_bar;
2866 const double vel_thresh = 0.0;
2868 Track_const_riterator pair_p;
2871 pair_p = track.rbegin();
2878 pair_p = Track_const_riterator(track.upper_bound(t));
2880 "The frame passed to track_future() must be valid.");
2885 size_t wnd = t - pair_p->first;
2889 if(pair_p == track.rend())
break;
2891 wnd = t - pair_p->first;
2893 x.
push_back(m_convert(*pair_p->second)[0]);
2894 y.
push_back(m_convert(*pair_p->second)[1]);
2899 double x_1 = x.
back();
2900 double x_n = x.
front();
2903 x -=
Vector(x.get_length(), x_bar);
2906 Matrix A(x.get_length(), 2);
2915 double r = a*x_bar + b*y_bar;
2916 double y_1 = (r - a*x_1)/b;
2917 double y_n = (r - a*x_n)/b;
2926 return std::make_pair(x_s, vel);
2931 return track_future(track, t);
2936 template <
class Track>
2944 const size_t wsz = m_b_bar;
2945 const double vel_thresh = 0.0;
2947 Track_const_iterator pair_p;
2950 pair_p = track.begin();
2957 pair_p = track.lower_bound(t);
2959 "The frame passed to track_future() must be valid.");
2964 size_t wnd = pair_p->first - t;
2968 if(pair_p == track.end())
break;
2970 wnd = pair_p->first - t;
2972 x.
push_back(m_convert(*pair_p->second)[0]);
2973 y.
push_back(m_convert(*pair_p->second)[1]);
2978 double x_1 = x.
back();
2979 double x_n = x.
front();
2982 x -=
Vector(x.get_length(), x_bar);
2985 Matrix A(x.get_length(), 2);
2994 double r = a*x_bar + b*y_bar;
2995 double y_1 = (r - a*x_1)/b;
2996 double y_n = (r - a*x_n)/b;
3005 return std::make_pair(x_s, vel);
3010 return track_future(track, t);
3015 template <
class Track>
3019 BOOST_FOREACH(
const Track& track, w)
3021 size_t track_sz = track.real_size();
3024 nsp += (track_sz - 3);
3033 template <
class Track>
3034 size_t Proposer<Track>::count_merge_pairs(
const Assoc& w)
const
3037 for(Assoc_const_iterator track_p1 = w.begin();
3038 track_p1 != w.end();
3041 int t_f = track_p1->get_end_time();
3042 for(Assoc_const_iterator track_p2 = w.begin();
3043 track_p2 != w.end();
3046 int t_1 = track_p2->get_start_time();
3047 if(track_p1 != track_p2 && t_f < t_1)
3050 track_p2->get_start_point(),
3051 t_1 - t_f, m_d_bar, m_v_bar,
3052 m_noise_sigma, m_convert))
3065 template <
class Track>
3066 size_t Proposer<Track>::count_secretion_tracks(
const Assoc& w)
const
3069 BOOST_FOREACH(
const Track& track, w)
3071 if(track.real_size() >= 4)
3077 if(track.count(track.get_start_time()) > 1
3078 && track.count(track.get_end_time()) > 1)
3090 template <
class Track>
3091 size_t Proposer<Track>::count_absorption_pairs(
const Assoc& w)
const
3111 return (w.size() * (w.size() - 1)) / 2;
double gamma() const
Return probability of stopping a growing track.
Definition: mcmcda_proposer.h:184
void grow_track_forward(Track &track, const Assoc &w) const
Grows a track forward according to MCMCDA.
Definition: mcmcda_proposer.h:2070
Definition: mcmcda_proposer.h:57
double fwd
Definition: mcmcda_proposer.h:59
double p_birth(const Assoc &w, const Assoc &w_p) const
Probability of birth.
Definition: mcmcda_proposer.h:1517
double get_velocity_noise_score(int d, const Vector &vel) const
Helper function that gives score of a noise thingy.
Definition: mcmcda_proposer.h:356
Definition for the Matrix class, a thin wrapper on the KJB Matrix struct and its related functionalit...
Int_matrix::Value_type max(const Int_matrix &mat)
Return the maximum value in this matrix.
Definition: l_int_matrix.h:1397
reference back()
Returns the last element of the vector.
Definition: m_vector.h:1107
bool empty() const
Returns true iff size is zero. Required to comply with stl Container concept.
Definition: m_vector.h:526
void diagonalize(const Matrix &M, Matrix &eig_vectors, Vector &eig_values, bool symmetric=false)
KJB c-style syntax for eigenvalue decomposition.
Definition: n_eig.h:39
Elementp_set get_dead_points_at_time(int t) const
Computes points at time t that do not belong to any track.
Definition: mcmcda_association.h:170
Definition of various standard probability distributions.
void set_col(int col, const Generic_vector &v)
Replace a column of the matrix with the given vector. "vector" can be any collection of values conver...
Definition: m_matrix.h:1726
Parent::const_iterator const_iterator
Definition: mcmcda_association.h:122
int get_length() const
Return the length of the vector.
Definition: m_vector.h:501
Definition: mcmcda_proposer.h:82
Parent::iterator iterator
Definition: mcmcda_association.h:121
double p_extension(const Assoc &w, const Assoc &w_p) const
Probability of extension.
Definition: mcmcda_proposer.h:1627
void push_back(Value_type x)
inserts an element at the back of the vector in amortized constant time.
Definition: m_vector.h:1128
r
Definition: APPgetLargeConnectedEdges.m:127
Probability_map get_growing_probabilities(const Track &track, int t, std::set< const Element * > &candidates, const Vector &x, const Vector &vel, int d) const
Helper function that gives probabilities of growth candidates.
Definition: mcmcda_proposer.h:2445
double propose_absorption(const Assoc &w, Assoc &w_p) const
Propose absorption.
Definition: mcmcda_proposer.h:1339
double propose_extension(const Assoc &w, Assoc &w_p) const
Propose exteision.
Definition: mcmcda_proposer.h:964
Definition: mcmcda_proposer.h:76
This class implements vectors, in the linear-algebra sense, with real-valued elements.
Definition: m_vector.h:87
A class that represents a MCMCDA association.
Definition: mcmcda_association.h:66
double pdf(const MV_gaussian_distribution &P, const Vector &x)
Computes the joint PDF of a multivariate Gaussian at x.
Definition: prob_pdf.cpp:32
double p_switch(const Assoc &w, const Assoc &w_p) const
Probability of switch.
Definition: mcmcda_proposer.h:1687
std::pair< Vector, Vector > track_past_tls(const Track &track, int t=-1) const
Estimate the velocity of the detections going backward.
Definition: mcmcda_proposer.h:2938
Vector sample(const MV_gaussian_distribution &dist)
Sample from a multivariate normal distribution.
Definition: prob_sample.cpp:42
reference front()
Returns the first element of the vector.
Definition: m_vector.h:1086
double propose_reduction(const Assoc &w, Assoc &w_p) const
Propose reduction.
Definition: mcmcda_proposer.h:1031
Value_type mean(Iterator begin, Iterator end, const Value_type &)
Definition: prob_stat.h:56
size_t sample_move(const Assoc &w) const
Sample a move type.
Definition: mcmcda_proposer.h:188
Definition: mcmcda_proposer.h:83
#define IFT(a, ex, msg)
Definition: l_exception.h:101
double dist(const pt a, const pt b)
compute approx. Great Circle distance between two UTM points
Definition: layer.cpp:45
double p_grow_track_forward(const Track &track, const Assoc &w, int t) const
Probability of growing a track forward to what it is.
Definition: mcmcda_proposer.h:2182
double p_death(const Assoc &w, const Assoc &w_p) const
Probability of death.
Definition: mcmcda_proposer.h:1557
mh_proposal_result(double fp, double rp, const std::string &nm="")
Definition: mcmcda_proposer.h:63
Int_matrix matrix_transpose(const Int_matrix &op1)
Test for any difference between two matrices.
Definition: l_int_matrix.h:1331
iterator end()
Definition: m_vector.h:557
Definition: mcmcda_proposer.h:81
int d_bar() const
Return maximum number of missed frames.
Definition: mcmcda_proposer.h:181
kjb_c::Pixel abs(const kjb_c::Pixel &p)
Take the channel-wise absolute value of a kjb_c::Pixel.
Definition: i_pixel.h:354
iterator begin()
Definition: m_vector.h:537
void swap_tracks(Generic_track< Element > &track1, Generic_track< Element > &track2, int t1, int tp1, int t2, int tq1)
Mergest two tracks.
Definition: mcmcda_track.h:295
PDFs and CDFs for the different distributions defined in "prob_distribution.h".
double rev
Definition: mcmcda_proposer.h:60
x
Definition: APPgetLargeConnectedEdges.m:100
Vector & set(Value_type val)
Clone of zero_out(int)
Definition: m_vector.h:707
boost::function1< Vector, const Element & > Convert
Definition: mcmcda_data.h:53
double log_pdf(const MV_gaussian_distribution &P, const Vector &x)
Computes the log PDF a multivariate normal distribution at x.
Definition: prob_pdf.cpp:64
double propose_secretion(const Assoc &w, Assoc &w_p) const
Propose secretion.
Definition: mcmcda_proposer.h:1202
void grow_track_backward(Track &track, const Assoc &w) const
Grows a track backward according to MCMCDA.
Definition: mcmcda_proposer.h:2127
Value_type magnitude() const
Return this vector's magnitude.
Definition: m_vector.h:1490
const Data< Element > & get_data() const
Returns data set const-ref.
Definition: mcmcda_association.h:126
std::string name
Definition: mcmcda_proposer.h:61
double propose_death(const Assoc &w, Assoc &w_p) const
Propose death.
Definition: mcmcda_proposer.h:790
double propose_split(const Assoc &w, Assoc &w_p) const
Propose split.
Definition: mcmcda_proposer.h:814
count
Definition: APPgetLargeConnectedEdges.m:71
#define KJB_THROW_3(ex, fmt, params)
Definition: l_exception.h:56
Definition: mcmcda_proposer.h:84
Move
Definition: mcmcda_proposer.h:74
Various utility functions for probabiliity-related stuff.
double p_secretion(const Assoc &w, const Assoc &w_p) const
Probability of secretion.
Definition: mcmcda_proposer.h:1713
std::pair< Vector, Vector > track_past(const Track &track, int t=-1) const
Estimate the velocity of the detections going backward.
Definition: mcmcda_proposer.h:2585
int count_dead_points_at_time(int t) const
Counts points at time t that do not belong to any track. This is usually faster than computing the po...
Definition: mcmcda_association.h:191
double propose_switch(const Assoc &w, Assoc &w_p) const
Propose switch.
Definition: mcmcda_proposer.h:1097
Int_matrix::Value_type min(const Int_matrix &mat)
Return the minimum value in this matrix.
Definition: l_int_matrix.h:1385
boost::math::normal Normal_distribution
Definition: prob_distribution.h:68
double p_absorption(const Assoc &w, const Assoc &w_p) const
Probability of absorption.
Definition: mcmcda_proposer.h:1745
Definition: mcmcda_proposer.h:77
Sampling functionality for the different distributions defined in "prob_distributions.h".
double propose_birth(const Assoc &w, Assoc &w_p) const
Propose birth.
Definition: mcmcda_proposer.h:740
std::pair< Vector, Vector > track_future_ls(const Track &track, int t=-1) const
Estimate the velocity of the detections going forward.
Definition: mcmcda_proposer.h:2660
Iterator element_uar(Iterator first, distance_type size)
Pick an element uniformly at random (UAR) from a sequence, represented by a beginning iterator and a ...
Definition: prob_sample.h:228
std::pair< Vector, Vector > track_future_tls(const Track &track, int t=-1) const
Estimate the velocity of the detections going forward.
Definition: mcmcda_proposer.h:2859
double propose_merge(const Assoc &w, Assoc &w_p) const
Propose merge.
Definition: mcmcda_proposer.h:896
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
Object thrown when an argument to a function is not acceptable.
Definition: l_exception.h:377
Definition: mcmcda_proposer.h:92
void make_p_move(const Categorical_distribution< size_t > &move_distribution)
Creates move distribution.
Definition: mcmcda_proposer.h:315
double v_bar() const
Return maximum velocity of targets.
Definition: mcmcda_proposer.h:178
double p_merge(const Assoc &w, const Assoc &w_p) const
Probability of merge.
Definition: mcmcda_proposer.h:1601
get the indices of edges in each direction for i
Definition: APPgetLargeConnectedEdges.m:48
This class implements matrices, in the linear-algebra sense, with real-valued elements.
Definition: m_matrix.h:94
boost::math::uniform Uniform_distribution
Definition: prob_distribution.h:70
for m
Definition: APPgetLargeConnectedEdges.m:64
double log_sum(double n1, double n2)
Definition: prob_util.h:66
Support for error handling exception classes in libKJB.
double p_split(const Assoc &w, const Assoc &w_p) const
Probability of split.
Definition: mcmcda_proposer.h:1576
Definition: mcmcda_proposer.h:80
std::pair< Vector, Vector > track_future(const Track &track, int t=-1) const
Estimate the velocity of the detections going forward.
Definition: mcmcda_proposer.h:2510
Definition for the Vector class, a thin wrapper on the KJB Vector struct and its related functionalit...
double p_grow_track_backward(const Track &track, const Assoc &w, int t) const
Probability of growing a track backward to what it is.
Definition: mcmcda_proposer.h:2306
double p_reduction(const Assoc &w, const Assoc &w_p) const
Probability of reduction.
Definition: mcmcda_proposer.h:1664
Object thrown when computation fails somehow during execution.
Definition: l_exception.h:321
Definition: mcmcda_proposer.h:85
std::pair< Vector, Vector > track_past_ls(const Track &track, int t=-1) const
Estimate the velocity of the detections going backward.
Definition: mcmcda_proposer.h:2759
double p_choose_nothing(const Probability_map &probs) const
Probability of not choosing any detections.
Definition: mcmcda_proposer.h:2426
double move_log_pdf(Move m, const Assoc &w) const
Compute the log pdf of a move type.
Definition: mcmcda_proposer.h:195