KJB
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
mcmcda_proposer.h
Go to the documentation of this file.
1 /* =========================================================================== *
2  |
3  | Copyright (c) 1994-2011 by Kobus Barnard (author)
4  |
5  | Personal and educational use of this code is granted, provided that this
6  | header is kept intact, and that the authorship is not misrepresented, that
7  | its use is acknowledged in publications, and relevant papers are cited.
8  |
9  | For other use contact the author (kobus AT cs DOT arizona DOT edu).
10  |
11  | Please note that the code in this file has not necessarily been adequately
12  | tested. Naturally, there is no guarantee of performance, support, or fitness
13  | for any particular task. Nonetheless, I am interested in hearing about
14  | problems that you encounter.
15  |
16  | Author: Ernesto Brau
17  * =========================================================================== */
18 
19 #ifndef MCMCDA_PROPOSER_H_INCLUDED
20 #define MCMCDA_PROPOSER_H_INCLUDED
21 
24 #include <mcmcda_cpp/mcmcda_data.h>
25 #include <m_cpp/m_vector.h>
26 #include <m_cpp/m_matrix.h>
27 #include <n_cpp/n_eig.h>
29 #include <prob_cpp/prob_sample.h>
30 #include <prob_cpp/prob_pdf.h>
31 #include <prob_cpp/prob_util.h>
32 #include <prob_cpp/prob_stat.h>
33 #include <l_cpp/l_util.h>
34 #include <l_cpp/l_exception.h>
35 #include <l/l_sys_lib.h>
36 #include <set>
37 #include <algorithm>
38 #include <iterator>
39 #include <functional>
40 #include <list>
41 #include <utility>
42 #include <limits>
43 #include <map>
44 #include <vector>
45 #include <cmath>
46 #include <boost/foreach.hpp>
47 #include <boost/bimap.hpp>
48 #include <boost/bimap/set_of.hpp>
49 #include <boost/bimap/multiset_of.hpp>
50 
51 #ifdef KJB_HAVE_ERGO
52 #include <ergo/mh.h>
53 #else
54 // if we don't have libergo installed, this still compiles
55 namespace ergo {
56 
58 {
59  double fwd;
60  double rev;
61  std::string name;
62 
63  mh_proposal_result(double fp, double rp, const std::string& nm = "") :
64  fwd(fp), rev(rp), name(nm) {}
65 };
66 
67 }
68 #endif
69 
70 namespace kjb {
71 namespace mcmcda {
72 
73 #warning "[Code police] scoped enums are a c++0x feature."
74 enum Move : size_t
75 {
78  //MCMCDA_SPLIT,
79  //MCMCDA_MERGE,
86 };
87 
88 using namespace boost::bimaps;
89 
91 template <class Track>
92 class Proposer
93 {
94 private:
95  typedef Association<Track> Assoc;
96  typedef typename Track::Element Element;
97  typedef typename Assoc::const_iterator Assoc_const_iterator;
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;
105 
106 
107 private:
108  std::vector<Categorical_distribution<size_t> > m_p_move;
109  double m_v_bar;
110  int m_d_bar;
111  int m_b_bar;
112  double m_gamma;
113  typename Data<Element>::Convert m_convert;
114  Avg m_average;
115  Feature_prob m_feature_prob;
116  double m_noise_sigma;
117  Track def_track;
118  bool use_feature;
119 
120 public:
121  /* @brief C-tor. */
122  Proposer
123  (
124  const Categorical_distribution<size_t>& move_distribution,
125  double v_bar,
126  int d_bar,
127  int b_bar,
128  double gamma,
129  const typename Data<Element>::Convert& convert_to_vector,
130  const Avg& avg,
131  double detection_noise_variance,
132  const Track& empty_track
133  ) :
134  m_v_bar(v_bar),
135  m_d_bar(d_bar),
136  m_b_bar(b_bar),
137  m_gamma(gamma),
138  m_convert(convert_to_vector),
139  m_average(avg),
140  m_noise_sigma(detection_noise_variance),
141  def_track(empty_track),
142  use_feature(false)
143  {
144  make_p_move(move_distribution);
145  }
146 
147  Proposer
148  (
149  const Categorical_distribution<size_t>& move_distribution,
150  double v_bar,
151  int d_bar,
152  int b_bar,
153  double gamma,
154  const typename Data<Element>::Convert& convert_to_vector,
155  const Avg& avg,
156  const Feature_prob& feature_prob,
157  double detection_noise_variance,
158  const Track& empty_track
159  ) :
160  m_v_bar(v_bar),
161  m_d_bar(d_bar),
162  m_b_bar(b_bar),
163  m_gamma(gamma),
164  m_convert(convert_to_vector),
165  m_average(avg),
166  m_feature_prob(feature_prob),
167  m_noise_sigma(detection_noise_variance),
168  def_track(empty_track),
169  use_feature(true)
170  {
171  make_p_move(move_distribution);
172  }
173 
174  /* @brief Apply this proposer. */
175  ergo::mh_proposal_result operator()(const Assoc& w, Assoc& w_p) const;
176 
178  double v_bar() const { return m_v_bar; }
179 
181  int d_bar() const { return m_d_bar; }
182 
184  double gamma() const { return m_gamma; }
185 
187  //Move sample_move(const Assoc& w) const
188  size_t sample_move(const Assoc& w) const
189  {
190  size_t K = w.size() < 2 ? w.size() : 2;
191  return kjb::sample(m_p_move[K]);
192  }
193 
195  double move_log_pdf(Move m, const Assoc& w) const
196  {
197  size_t K = w.size() < 2 ? w.size() : 2;
198  return log_pdf(m_p_move[K], m);
199  }
200 
201 public:
203  double propose_birth(const Assoc& w, Assoc& w_p) const;
204 
206  double propose_death(const Assoc& w, Assoc& w_p) const;
207 
209  double propose_split(const Assoc& w, Assoc& w_p) const;
210 
212  double propose_merge(const Assoc& w, Assoc& w_p) const;
213 
215  double propose_extension(const Assoc& w, Assoc& w_p) const;
216 
218  double propose_reduction(const Assoc& w, Assoc& w_p) const;
219 
221  double propose_switch(const Assoc& w, Assoc& w_p) const;
222 
224  double propose_secretion(const Assoc& w, Assoc& w_p) const;
225 
227  double propose_absorption(const Assoc& w, Assoc& w_p) const;
228 
229 public:
231  double p_birth(const Assoc& w, const Assoc& w_p) const;
232 
234  double p_death(const Assoc& w, const Assoc& w_p) const;
235 
237  double p_split(const Assoc& w, const Assoc& w_p) const;
238 
240  double p_merge(const Assoc& w, const Assoc& w_p) const;
241 
243  double p_extension(const Assoc& w, const Assoc& w_p) const;
244 
246  double p_reduction(const Assoc& w, const Assoc& w_p) const;
247 
249  double p_switch(const Assoc& w, const Assoc& w_p) const;
250 
252  double p_secretion(const Assoc& w, const Assoc& w_p) const;
253 
255  double p_absorption(const Assoc& w, const Assoc& w_p) const;
256 
257 private:
259  bool is_valid_birth(const Assoc& w, const Assoc& w_p) const;
260 
262  bool is_valid_death(const Assoc& w, const Assoc& w_p) const;
263 
265  bool is_valid_split(const Assoc& w, const Assoc& w_p) const;
266 
268  bool is_valid_merge(const Assoc& w, const Assoc& w_p) const;
269 
271  bool is_valid_extension(const Assoc& w, const Assoc& w_p) const;
272 
274  bool is_valid_reduction(const Assoc& w, const Assoc& w_p) const;
275 
277  bool is_valid_switch(const Assoc& w, const Assoc& w_p) const;
278 
280  bool is_valid_secretion(const Assoc& w, const Assoc& w_p) const;
281 
283  bool is_valid_absorption(const Assoc& w, const Assoc& w_p) const;
284 
285 public:
287  void grow_track_forward(Track& track, const Assoc& w) const;
288 
290  void grow_track_backward(Track& track, const Assoc& w) const;
291 
293  double p_grow_track_forward
294  (
295  const Track& track,
296  const Assoc& w,
297  int t
298  ) const;
299 
301  double p_grow_track_backward
302  (
303  const Track& track,
304  const Assoc& w,
305  int t
306  ) const;
307 
309  double p_choose_nothing
310  (
311  const Probability_map& probs
312  ) const;
313 
315  void make_p_move(const Categorical_distribution<size_t>& move_distribution)
316  {
317  std::vector<double> ps(MCMCDA_NUM_MOVES, 0.0);
318 
319  // distribution for K = 0
320  ps[MCMCDA_BIRTH] = 1.0;
321  m_p_move.push_back(Categorical_distribution<size_t>(ps, 0));
322 
323  // distribution for K = 1
324  ps[MCMCDA_BIRTH] = 0.3;
325  ps[MCMCDA_DEATH] = 0.1;
326  ps[MCMCDA_EXTENSION] = 0.4;
327  ps[MCMCDA_REDUCTION] = 0.1;
328  ps[MCMCDA_SWITCH] = 0.0;
329  ps[MCMCDA_SECRETION] = 0.1;
330  ps[MCMCDA_ABSORPTION] = 0.0;
331  m_p_move.push_back(Categorical_distribution<size_t>(ps, 0));
332 
333  /*std::fill(ps.begin(), ps.end(), 1.0 / (MCMCDA_NUM_MOVES - 2));
334  //ps[MCMCDA_MERGE] = 0.0;
335  ps[MCMCDA_SWITCH] = 0.0;
336  ps[MCMCDA_ABSORPTION] = 0.0;
337  m_p_move.push_back(Categorical_distribution<size_t>(ps, 0));*/
338 
339  // distribution for K > 1
340  m_p_move.push_back(move_distribution);
341  }
342 
343 public:
345  Probability_map get_growing_probabilities
346  (
347  const Track& track,
348  int t,
349  std::set<const Element*>& candidates,
350  const Vector& x,
351  const Vector& vel,
352  int d
353  ) const;
354 
356  double get_velocity_noise_score(int d, const Vector& vel) const
357  {
358  d = std::abs(d);
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);
362  Normal_distribution N_u(0.0, sg_u);
363  Normal_distribution N_v(0.0, sg_v);
364  //double p_good = 0.5;
365  return pdf(N_u, sg_u) * pdf(N_v, sg_v);
366  }
367 
369  std::pair<Vector, Vector> track_future
370  (
371  const Track& track,
372  int t = -1
373  ) const;
374 
376  std::pair<Vector, Vector> track_past
377  (
378  const Track& track,
379  int t = -1
380  ) const;
381 
383  std::pair<Vector, Vector> track_future_ls
384  (
385  const Track& track,
386  int t = -1
387  ) const;
388 
390  std::pair<Vector, Vector> track_past_ls
391  (
392  const Track& track,
393  int t = -1
394  ) const;
395 
397  std::pair<Vector, Vector> track_future_tls
398  (
399  const Track& track,
400  int t = -1
401  ) const;
402 
404  std::pair<Vector, Vector> track_past_tls
405  (
406  const Track& track,
407  int t = -1
408  ) const;
409 
410 private:
412  size_t count_split_points(const Assoc& w) const;
413 
415  size_t count_merge_pairs(const Assoc& w) const;
416 
418  size_t count_secretion_tracks(const Assoc& w) const;
419 
421  size_t count_absorption_pairs(const Assoc& w) const;
422 
424  double negative_infinity() const
425  {
427  }
428 
429  // constants
430  static const int GROW_FORWARD;
431  static const int GROW_BACKWARD;
432 
433 private:
435  struct Birth_info
436  {
437  Assoc_const_iterator new_track_p;
438  };
439 
441  struct Death_info
442  {
443  size_t num_tracks;
444  };
445 
447  struct Split_info
448  {
449  size_t num_split_points;
450  };
451 
453  struct Merge_info
454  {
455  size_t num_merge_pairs;
456  };
457 
459  struct Extension_info
460  {
461  size_t num_tracks;
462  Assoc_const_iterator extended_track_p;
463  int direction;
464  int previous_end;
465  };
466 
468  struct Reduction_info
469  {
470  size_t num_tracks;
471  size_t reduced_track_size;
472  };
473 
475  struct Switch_info
476  {
477  size_t num_switch_points;
478  };
479 
481  struct Secretion_info
482  {
483  size_t num_valid_tracks;
484  size_t window_size;
485  };
486 
488  struct Absorption_info
489  {
490  size_t num_valid_track_pairs;
491  };
492 
493 private:
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;
503 };
504 
505 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
506 
507 template <class Track>
508 const int Proposer<Track>::GROW_FORWARD = 1;
509 
510 template <class Track>
511 const int Proposer<Track>::GROW_BACKWARD = 2;
512 
513 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
514 
515 template <class Track>
516 ergo::mh_proposal_result Proposer<Track>::operator()
517 (
518  const Assoc& w,
519  Assoc& w_p
520 ) const
521 {
522  double fwd, rev;
523  bool keep_going = true;
524  std::string move_name = "association-";
525 
526  while(keep_going)
527  {
528  //Move m = sample_move(w);
529  size_t m = sample_move(w);
530 
531  switch(m)
532  {
533  case MCMCDA_BIRTH:
534  {
535  fwd = propose_birth(w, w_p);
536  if(fwd == negative_infinity())
537  {
538  break;
539  }
540 
541  rev = p_death(w_p, w);
542  if(rev == negative_infinity())
543  {
544  break;
545  }
546 
547  fwd += move_log_pdf(MCMCDA_BIRTH, w);
548  rev += move_log_pdf(MCMCDA_DEATH, w_p);
549 
550  move_name += "birth";
551  break;
552  }
553 
554  case MCMCDA_DEATH:
555  {
556  fwd = propose_death(w, w_p);
557  if(fwd == negative_infinity())
558  {
559  break;
560  }
561 
562  rev = p_birth(w_p, w);
563  if(rev == negative_infinity())
564  {
565  break;
566  }
567 
568  fwd += move_log_pdf(MCMCDA_DEATH, w);
569  rev += move_log_pdf(MCMCDA_BIRTH, w_p);
570 
571  move_name += "death";
572  break;
573  }
574 
575  /*case MCMCDA_SPLIT:
576  {
577  fwd = propose_split(w, w_p);
578  if(fwd == negative_infinity())
579  {
580  break;
581  }
582 
583  rev = p_merge(w_p, w);
584  if(rev == negative_infinity())
585  {
586  break;
587  }
588 
589  fwd += move_log_pdf(MCMCDA_SPLIT, w);
590  rev += move_log_pdf(MCMCDA_MERGE, w_p);
591 
592  move_name += "split";
593  break;
594  }
595 
596  case MCMCDA_MERGE:
597  {
598  fwd = propose_merge(w, w_p);
599  if(fwd == negative_infinity())
600  {
601  break;
602  }
603 
604  rev = p_split(w_p, w);
605  if(rev == negative_infinity())
606  {
607  break;
608  }
609 
610  fwd += move_log_pdf(MCMCDA_MERGE, w);
611  rev += move_log_pdf(MCMCDA_SPLIT, w_p);
612 
613  move_name += "merge";
614  break;
615  }*/
616 
617  case MCMCDA_EXTENSION:
618  {
619  fwd = propose_extension(w, w_p);
620  if(fwd == negative_infinity())
621  {
622  break;
623  }
624 
625  rev = p_reduction(w_p, w);
626  if(rev == negative_infinity())
627  {
628  break;
629  }
630 
631  fwd += move_log_pdf(MCMCDA_EXTENSION, w);
632  rev += move_log_pdf(MCMCDA_REDUCTION, w_p);
633 
634  move_name += "extension";
635  break;
636  }
637 
638  case MCMCDA_REDUCTION:
639  {
640  fwd = propose_reduction(w, w_p);
641  if(fwd == negative_infinity())
642  {
643  break;
644  }
645 
646  rev = p_extension(w_p, w);
647  if(rev == negative_infinity())
648  {
649  break;
650  }
651 
652  fwd += move_log_pdf(MCMCDA_REDUCTION, w);
653  rev += move_log_pdf(MCMCDA_EXTENSION, w_p);
654 
655  move_name += "reduction";
656  break;
657  }
658 
659  case MCMCDA_SWITCH:
660  {
661  fwd = propose_switch(w, w_p);
662  if(fwd == negative_infinity())
663  {
664  break;
665  }
666 
667  fwd += move_log_pdf(MCMCDA_SWITCH, w);
668 
669  // recall that switch move is symmetric
670  rev = fwd;
671 
672  move_name += "switch";
673  break;
674  }
675 
676  case MCMCDA_SECRETION:
677  {
678  fwd = propose_secretion(w, w_p);
679  if(fwd == negative_infinity())
680  {
681  break;
682  }
683 
684  rev = p_absorption(w_p, w);
685  if(rev == negative_infinity())
686  {
687  break;
688  }
689 
690  fwd += move_log_pdf(MCMCDA_SECRETION, w);
691  rev += move_log_pdf(MCMCDA_ABSORPTION, w_p);
692  //fwd = rev = 0.0;
693 
694  move_name += "secretion";
695  break;
696  }
697 
698  case MCMCDA_ABSORPTION:
699  {
700  fwd = propose_absorption(w, w_p);
701  if(fwd == negative_infinity())
702  {
703  break;
704  }
705 
706  rev = p_secretion(w_p, w);
707  if(rev == negative_infinity())
708  {
709  break;
710  }
711 
712  fwd += move_log_pdf(MCMCDA_ABSORPTION, w);
713  rev += move_log_pdf(MCMCDA_SECRETION, w_p);
714  //fwd = rev = 0.0;
715 
716  move_name += "absorption";
717  break;
718  }
719 
720  default:
721  {
723  "MCMCDA: move %d does not exist.", (m));
724  break;
725  }
726  }
727 
728  if(fwd != negative_infinity() && rev != negative_infinity())
729  {
730  keep_going = false;
731  }
732  }
733 
734  return ergo::mh_proposal_result(fwd, rev, move_name);
735 }
736 
737 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
738 
739 template <class Track>
740 double Proposer<Track>::propose_birth(const Assoc& w, Assoc& w_p) const
741 {
742  //KJB(ASSERT(w.is_valid(m_v_bar, m_d_bar, m_noise_sigma, m_convert)));
743 
744  Track new_track = def_track;
745 
746  // pick first point in new track
747  size_t t_1 = kjb::sample(
748  Categorical_distribution<size_t>(1, w.get_data().size(), 1));
749 
750  std::set<const Element*> L_1 = w.get_dead_points_at_time(t_1);
751 
752  if(L_1.empty())
753  {
754  return negative_infinity();
755  }
756 
757  //const Element* x_1 = *element_uar(L_1.begin(), L_1.size());
758  std::vector<double> ps(L_1.size(), 1.0);
760  size_t n = kjb::sample(P);
761  typename std::set<const Element*>::const_iterator x_1_p = L_1.begin();
762 
763  std::advance(x_1_p, n);
764  new_track.insert(std::make_pair(t_1, *x_1_p));
765 
766  // grow track
767  grow_track_forward(new_track, w);
768  if(new_track.real_size() == 1)
769  {
770  return negative_infinity();
771  }
772 
773  // add new track to sampled association
774  w_p = w;
775  Assoc_const_iterator new_track_p = w_p.insert(new_track).first;
776 
777  // set changed flag
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());
780 
781  // compute forward probability
782  birth_info.new_track_p = new_track_p;
783  death_info.num_tracks = w_p.size();
784  return p_birth(w, w_p);
785 }
786 
787 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
788 
789 template <class Track>
790 double Proposer<Track>::propose_death(const Assoc& w, Assoc& w_p) const
791 {
792  //KJB(ASSERT(w.is_valid(m_v_bar, m_d_bar, m_noise_sigma, m_convert)));
793 
794  if(w.empty())
795  {
796  std::cout << "WARNING: cannot propose death for empty association.\n";
797  return std::log(0.0);
798  }
799 
800  // remove track UAR
801  w_p = w;
802  Assoc_const_iterator new_track_p = element_uar(w.begin(), w.size());
803  w_p.erase(*new_track_p);
804 
805  // compute forward probability
806  death_info.num_tracks = w.size();
807  birth_info.new_track_p = new_track_p;
808  return p_death(w, w_p);
809 }
810 
811 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
812 
813 template <class Track>
814 double Proposer<Track>::propose_split(const Assoc& w, Assoc& w_p) const
815 {
816  //KJB(ASSERT(w.is_valid(m_v_bar, m_d_bar, m_noise_sigma, m_convert)));
817 
818  if(w.empty())
819  {
820  std::cout << "WARNING: cannot propose split for empty association.\n";
821  return std::log(0.0);
822  }
823 
824  size_t nsp = count_split_points(w);
825  if(nsp == 0)
826  {
827  return negative_infinity();
828  }
829 
830  // pick random split point and compute it
831  size_t n = kjb::sample(Categorical_distribution<size_t>(1, nsp, 1));
832 
833  w_p = w;
834  Assoc_const_iterator track_p;
835  Track_const_iterator rp;
836  size_t count = 0;
837  for(track_p = w_p.begin(); track_p != w_p.end(); track_p++)
838  {
839  size_t tsz = track_p->real_size();
840  if(tsz >= 4)
841  {
842  if(count + tsz - 3 >= n)
843  {
844  int r = 0;
845  int t = 0;
846  Track_const_iterator pair_p;
847  for(pair_p = track_p->begin(); pair_p != track_p->end();
848  pair_p++)
849  {
850  if(pair_p->first == t)
851  {
852  continue;
853  }
854  t = pair_p->first;
855 
856  r++;
857  if(r == n - count + 1)
858  {
859  break;
860  }
861  }
862  rp = pair_p;
863  rp++;
864  break;
865  }
866  count += (tsz - 3);
867  }
868  }
869 
870  // split track into two tracks
871  Track beg = def_track;
872  Track fin = def_track;
873  beg.insert(track_p->begin(), rp);
874  fin.insert(rp, track_p->end());
875 
876  // set changed flags
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());
881 
882  // remove split track and add what it got split into
883  w_p.erase(track_p);
884  w_p.insert(beg);
885  w_p.insert(fin);
886 
887  // compute forward probability
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);
891 }
892 
893 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
894 
895 template <class Track>
896 double Proposer<Track>::propose_merge(const Assoc& w, Assoc& w_p) const
897 {
898  //KJB(ASSERT(w.is_valid(m_v_bar, m_d_bar, m_noise_sigma, m_convert)));
899 
900  if(w.size() < 2)
901  {
902  std::cout << "WARNING: cannot propose merge for single-track "
903  "association.\n";
904  return std::log(0.0);
905  }
906 
907  w_p = w;
908 
909  // compute merge candidates
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();
913  track_p1++)
914  {
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();
918  track_p2++)
919  {
920  int t_1 = track_p2->get_start_time();
921  if(track_p1 != track_p2 && t_f < t_1)
922  {
923  if(is_neighbor(track_p1->get_end_point(),
924  track_p2->get_start_point(),
925  t_1 - t_f, m_d_bar, m_v_bar,
926  m_noise_sigma, m_convert))
927  {
928  M.push_back(std::make_pair(track_p1, track_p2));
929  }
930  }
931  }
932  }
933 
934  if(M.empty())
935  {
936  return negative_infinity();
937  }
938 
939  // choose a candidate pair and merge its tracks into a single track
940  std::pair<Assoc_const_iterator, Assoc_const_iterator> ttm
941  = *element_uar(M.begin(), M.size());
942  Track new_track = *ttm.first;
943  new_track.insert(ttm.second->begin(), ttm.second->end());
944 
945  // set changed flags
946  new_track.set_changed_start(ttm.second->get_start_time());
947  new_track.set_changed_end(ttm.second->get_end_time());
948 
949  // remove merged tracks and add the result of merging them
950  w_p.erase(ttm.first);
951  w_p.erase(ttm.second);
952  w_p.insert(new_track);
953 
954  // compute forward probability
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);
958 }
959 
960 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
961 
962 template <class Track>
964 (
965  const Assoc& w,
966  Assoc& w_p
967 ) const
968 {
969  //KJB(ASSERT(w.is_valid(m_v_bar, m_d_bar, m_noise_sigma, m_convert)));
970 
971  if(w.empty())
972  {
973  std::cout << "WARNING: cannot propose extension for empty "
974  "association.\n";
975  return std::log(0.0);
976  }
977 
978  w_p = w;
979 
980  // choose track to extend
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();
984 
985  // grow track
986  int dir;
987  int prev_end;
988  double u = kjb::sample(Uniform_distribution());
989  if(u < 0.5)
990  {
991  dir = GROW_FORWARD;
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());
996  }
997  else
998  {
999  dir = GROW_BACKWARD;
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);
1004  }
1005 
1006 
1007  // if it did not grow...kill it!
1008  if(extended_track.size() == original_track_size)
1009  {
1010  return negative_infinity();
1011  }
1012 
1013  // replace track with extended version
1014  w_p.erase(track_p);
1015  Assoc_const_iterator extended_track_p = w_p.insert(extended_track).first;
1016 
1017  // compute forward probability
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);
1025 }
1026 
1027 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
1028 
1029 template <class Track>
1032  const Assoc& w,
1033  Assoc& w_p
1034 ) const
1035 {
1036  //KJB(ASSERT(w.is_valid(m_v_bar, m_d_bar, m_noise_sigma, m_convert)));
1037 
1038  if(w.empty())
1039  {
1040  std::cout << "WARNING: cannot propose reduction for empty "
1041  "association.\n";
1042  return std::log(0.0);
1043  }
1044 
1045  w_p = w;
1046 
1047  // choose track to reduce
1048  typename Assoc::iterator track_p = element_uar(w_p.begin(), w_p.size());
1049 
1050  // must be longer than two points
1051  if(track_p->real_size() <= 2)
1052  {
1053  return negative_infinity();
1054  }
1055 
1056  Track reduced_track = *track_p;
1057  size_t rtsz = reduced_track.real_size();
1058  size_t n = kjb::sample(Categorical_distribution<size_t>(2, rtsz - 1, 1));
1059  int t = reduced_track.get_nth_time(n);
1060 
1061  // choose direction and reduce
1062  int dir;
1063  if(kjb::sample(Uniform_distribution()) <= 0.5)
1064  {
1065  reduced_track.erase(reduced_track.upper_bound(t), reduced_track.end());
1066  dir = GROW_FORWARD;
1067  }
1068  else
1069  {
1070  reduced_track.erase(reduced_track.begin(), reduced_track.lower_bound(t));
1071  dir = GROW_BACKWARD;
1072  }
1073 
1074  // set changed flags
1075  reduced_track.set_changed_start(-1);
1076  reduced_track.set_changed_end(-1);
1077 
1078  // replace track with reduced version
1079  Assoc_const_iterator orig_track_p = w.find(*track_p);
1080  w_p.erase(track_p);
1081  w_p.insert(reduced_track);
1082 
1083  // compute forward probability
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);
1091 }
1092 
1093 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
1094 
1095 template <class Track>
1098  const Assoc& w,
1099  Assoc& w_p
1100 ) const
1101 {
1102  //KJB(ASSERT(w.is_valid(m_v_bar, m_d_bar, m_noise_sigma, m_convert)));
1103 
1104  using namespace boost;
1105 
1106  typedef tuple<const Track&, int, int, const Track&, int, int> Switch_point;
1107 
1108  std::list<Switch_point> switch_points;
1109  for(Assoc_const_iterator track_p1 = w.begin();
1110  track_p1 != w.end();
1111  track_p1++)
1112  {
1113  Assoc_const_iterator track_p2 = track_p1;
1114  for(track_p2++; track_p2 != w.end(); track_p2++)
1115  {
1116  // get first element in last frame of track1
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();
1120  pair_p1 != last_p1;
1121  pair_p1++)
1122  {
1123  // advance until last element of current frame
1124  int curt = pair_p1->first;
1125  while(pair_p1->first == curt)
1126  {
1127  pair_p1++;
1128  }
1129  pair_p1--;
1130 
1131  // get first element in last frame of track2
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();
1135  pair_p2 != last_p2;
1136  pair_p2++)
1137  {
1138  // advance until last element of current frame
1139  curt = pair_p2->first;
1140  while(pair_p2->first == curt)
1141  {
1142  pair_p2++;
1143  }
1144  pair_p2--;
1145 
1146  Track_const_iterator pair_q1 = pair_p1;
1147  Track_const_iterator pair_q2 = pair_p2;
1148 
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
1154  && is_neighbor(*pair_p1->second,
1155  *pair_q2->second, tq1 - t1, m_d_bar,
1156  m_v_bar, m_noise_sigma, m_convert)
1157  && is_neighbor(*pair_p2->second,
1158  *pair_q1->second, tp1 - t2, m_d_bar,
1159  m_v_bar, m_noise_sigma, m_convert))
1160  {
1161  switch_points.push_back(
1162  Switch_point(*track_p1, t1, tp1,
1163  *track_p2, t2, tq1));
1164  }
1165  }
1166  }
1167  }
1168  }
1169 
1170  if(switch_points.empty())
1171  {
1172  return negative_infinity();
1173  }
1174 
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);
1178  swap_tracks(track1_new, track2_new, get<1>(sp),
1179  get<2>(sp), get<4>(sp), get<5>(sp));
1180 
1181  // set changed flags
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());
1186 
1187  w_p = w;
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);
1192 
1193  // compute forward probability
1194  switch_info.num_switch_points = switch_points.size();
1195  return p_switch(w, w_p);
1196 }
1197 
1198 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
1199 
1200 template <class Track>
1203  const Assoc& w,
1204  Assoc& w_p
1205 ) const
1206 {
1207  //KJB(ASSERT(w.is_valid(m_v_bar, m_d_bar, m_noise_sigma, m_convert)));
1208 
1209  if(w.empty())
1210  {
1211  std::cout << "WARNING: cannot propose secretion for empty "
1212  "association.\n";
1213  return std::log(0.0);
1214  }
1215 
1216  w_p = w;
1217 
1218  typedef typename Assoc::iterator Assoc_iterator;
1219  //typedef typename Track::iterator Track_iterator;
1220 
1221  std::vector<Assoc_iterator> long_tracks;
1222  for(Assoc_iterator track_p = w_p.begin(); track_p != w_p.end(); track_p++)
1223  {
1224  if(track_p->real_size() > 4)
1225  {
1226  long_tracks.push_back(track_p);
1227  }
1228  else
1229  {
1230  if(track_p->count(track_p->get_start_time()) > 1
1231  && track_p->count(track_p->get_end_time()) > 1)
1232  {
1233  long_tracks.push_back(track_p);
1234  }
1235  }
1236  }
1237 
1238  if(long_tracks.empty())
1239  {
1240  return negative_infinity();
1241  }
1242 
1243  // choose random track
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;
1248 
1249  // choose separators
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);
1253 
1254  Track_const_iterator small_p, big_p;
1255  if(pair_p1->first < pair_p2->first)
1256  {
1257  small_p = pair_p1;
1258  big_p = pair_p2;
1259  }
1260  else if(pair_p1->first > pair_p2->first)
1261  {
1262  small_p = pair_p2;
1263  big_p = pair_p1;
1264  }
1265  else
1266  {
1267  if(std::distance(orig_p->begin(), pair_p1)
1268  < std::distance(orig_p->begin(), pair_p2))
1269  {
1270  small_p = pair_p1;
1271  big_p = pair_p2;
1272  }
1273  else
1274  {
1275  small_p = pair_p2;
1276  big_p = pair_p1;
1277  }
1278  }
1279 
1280  // do initial copying
1281  orig_track.insert(orig_p->begin(), ++small_p);
1282 
1284  if(kjb::sample(U) < 0.5)
1285  {
1286  orig_track.insert(++big_p, orig_p->end());
1287  }
1288  else
1289  {
1290  new_track.insert(++big_p, orig_p->end());
1291  }
1292 
1293  // coin flips in mid section
1294  size_t win_sz = 0;
1295  while(small_p != big_p)
1296  {
1297  if(kjb::sample(U) < 0.5)
1298  {
1299  orig_track.insert(*small_p);
1300  }
1301  else
1302  {
1303  new_track.insert(*small_p);
1304  }
1305 
1306  small_p++;
1307  win_sz++;
1308  }
1309 
1310  // make sure new tracks are valid
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))
1313  {
1314  return negative_infinity();
1315  }
1316 
1317  // set changed flags
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());
1322 
1323  // add new tracks / remove old one
1324  w_p.erase(orig_p);
1325  w_p.insert(orig_track);
1326  w_p.insert(new_track);
1327 
1328  // compute forward probability
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);
1333 }
1334 
1335 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
1336 
1337 template <class Track>
1340  const Assoc& w,
1341  Assoc& w_p
1342 ) const
1343 {
1344  //KJB(ASSERT(w.is_valid(m_v_bar, m_d_bar, m_noise_sigma, m_convert)));
1345 
1346  if(w.size() < 2)
1347  {
1348  std::cout << "WARNING: cannot propose absorption for single-track "
1349  "association.\n";
1350  return std::log(0.0);
1351  }
1352 
1353  w_p = w;
1354 
1355  // compute merge candidates
1356  using namespace boost;
1357  typedef tuple<Assoc_const_iterator, Assoc_const_iterator, Track> iit_tuple;
1358 
1359  std::vector<double> ps;
1360  ps.reserve((w.size()*(w.size()-1))/2);
1361  std::list<iit_tuple> M;
1362  Assoc w_e(w_p.get_data());
1363  for(Assoc_const_iterator track_p1 = w_p.begin();
1364  track_p1 != w_p.end();
1365  track_p1++)
1366  {
1367  Assoc_const_iterator track_p2 = track_p1;
1368  for(++track_p2; track_p2 != w_p.end(); ++track_p2)
1369  {
1370  // take longer track as "new" track
1371  Track new_track = def_track;
1372 
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();
1377 
1378  int l1 = ef1 - sf1 + 1;
1379  int l2 = ef2 - sf2 + 1;
1380 
1381  bool olap = (sf1 <= ef2 && sf2 <= ef1);
1382  size_t chst, ched;
1383  if(l1 > l2)
1384  {
1385  new_track = *track_p1;
1386  new_track.insert(track_p2->begin(), track_p2->end());
1387  if(olap)
1388  {
1389  chst = sf2;
1390  ched = ef2;
1391  }
1392  else
1393  {
1394  if(sf2 < sf1)
1395  {
1396  chst = sf2;
1397  ched = sf1;
1398  }
1399  else
1400  {
1401  chst = ef1;
1402  ched = ef2;
1403  }
1404  }
1405  }
1406  else
1407  {
1408  new_track = *track_p2;
1409  new_track.insert(track_p1->begin(), track_p1->end());
1410  if(olap)
1411  {
1412  chst = sf1;
1413  ched = ef1;
1414  }
1415  else
1416  {
1417  if(sf1 < sf2)
1418  {
1419  chst = sf1;
1420  ched = sf2;
1421  }
1422  else
1423  {
1424  chst = ef2;
1425  ched = ef1;
1426  }
1427  }
1428  }
1429 
1430  new_track.set_changed_start(chst);
1431  new_track.set_changed_end(ched);
1432 
1433  double pg = p_grow_track_forward(new_track, w_e,
1434  new_track.get_start_time());
1435  if(pg != negative_infinity())
1436  {
1437  M.push_back(boost::make_tuple(track_p1, track_p2, new_track));
1438  ps.push_back(pg);
1439  }
1440  }
1441  }
1442 
1443  if(ps.empty())
1444  {
1445  return negative_infinity();
1446  }
1447 
1448  // normalize
1449  typename std::list<iit_tuple>::const_iterator iit_p = M.begin();
1450  for(size_t i = 0; i < ps.size(); i++, iit_p++)
1451  {
1452  ps[i] /= get<2>(*iit_p).real_size();
1453  ps[i] *= -1.0;
1454  ps[i] = std::sqrt(ps[i]);
1455  ps[i] *= -1.0;
1456  }
1457 
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));
1463 
1464  // choose a candidate pair and merge its tracks into a single track
1466  typename std::list<iit_tuple>::iterator tuple_p = M.begin();
1467  std::advance(tuple_p, idx - 1);
1468  iit_tuple& ttm = *tuple_p;
1469 
1470  // figure out which track is which
1471  const Element* st_elem_p = get<0>(ttm)->begin()->second;
1472  if(get<0>(ttm)->begin()->second == get<2>(ttm).begin()->second)
1473  {
1474  st_elem_p = get<1>(ttm)->begin()->second;
1475  }
1476 
1477  const Element* en_elem_p = get<0>(ttm)->rbegin()->second;
1478  if(get<0>(ttm)->rbegin()->second == get<2>(ttm).rbegin()->second)
1479  {
1480  en_elem_p = get<1>(ttm)->rbegin()->second;
1481  }
1482 
1483  // compute length of overlap
1484  int last_start;
1485  int first_end;
1486  size_t i = 0;
1487  BOOST_FOREACH(const typename Track::value_type& pr, get<2>(ttm))
1488  {
1489  if(pr.second == st_elem_p)
1490  {
1491  last_start = i;
1492  }
1493 
1494  if(pr.second == en_elem_p)
1495  {
1496  first_end = i;
1497  }
1498 
1499  i++;
1500  }
1501 
1502  // remove merged tracks and add the result of merging them
1503  w_p.erase(get<0>(ttm));
1504  w_p.erase(get<1>(ttm));
1505  w_p.insert(get<2>(ttm));
1506 
1507  // compute forward probability
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);
1512 }
1513 
1514 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
1515 
1516 template <class Track>
1517 double Proposer<Track>::p_birth(const Assoc& w, const Assoc& w_p) const
1518 {
1519  //KJB(ASSERT(w.is_valid(m_v_bar, m_d_bar, m_noise_sigma, m_convert)));
1520  //KJB(ASSERT(w_p.is_valid(m_v_bar, m_d_bar, m_noise_sigma, m_convert)));
1521 
1522 #ifdef TEST
1523  if(!is_valid_birth(w, w_p))
1524  {
1525  std::cout << "WARNING: invalid birth." << std::endl;
1526  return std::log(0.0);
1527  }
1528 #endif
1529 
1530  const Track& extra_track = *birth_info.new_track_p;
1531  double p = 0.0;
1532 
1533  // compute first point and probability of choosing that time first
1534  Track_const_iterator pair_p = extra_track.begin();
1535  int t_1 = pair_p->first;
1536  p -= std::log(w.get_data().size());
1537 
1538  // compute probability of choosing first point
1539  int num_valid = w.count_dead_points_at_time(t_1);
1540 
1541  if(num_valid == 0)
1542  {
1543  return negative_infinity();
1544  }
1545 
1546  p -= std::log(num_valid);
1547 
1548  // probability of growing track to what it is
1549  p += p_grow_track_forward(extra_track, w, t_1);
1550 
1551  return p;
1552 }
1553 
1554 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
1555 
1556 template <class Track>
1557 double Proposer<Track>::p_death(const Assoc& w, const Assoc& w_p) const
1558 {
1559  //KJB(ASSERT(w.is_valid(m_v_bar, m_d_bar, m_noise_sigma, m_convert)));
1560  //KJB(ASSERT(w_p.is_valid(m_v_bar, m_d_bar, m_noise_sigma, m_convert)));
1561 
1562 #ifdef TEST
1563  if(!is_valid_death(w, w_p))
1564  {
1565  std::cout << "WARNING: invalid death." << std::endl;
1566  return std::log(0.0);
1567  }
1568 #endif
1569 
1570  return -std::log(death_info.num_tracks);
1571 }
1572 
1573 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
1574 
1575 template <class Track>
1576 double Proposer<Track>::p_split(const Assoc& w, const Assoc& w_p) const
1577 {
1578  //KJB(ASSERT(w.is_valid(m_v_bar, m_d_bar, m_noise_sigma, m_convert)));
1579  //KJB(ASSERT(w_p.is_valid(m_v_bar, m_d_bar, m_noise_sigma, m_convert)));
1580 
1581 #ifdef TEST
1582  if(!is_valid_split(w, w_p))
1583  {
1584  std::cout << "WARNING: invalid split." << std::endl;
1585  return std::log(0.0);
1586  }
1587 #endif
1588 
1589  size_t nsp = split_info.num_split_points;
1590  if(nsp == 0)
1591  {
1592  return negative_infinity();
1593  }
1594 
1595  return -std::log(nsp);
1596 }
1597 
1598 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
1599 
1600 template <class Track>
1601 double Proposer<Track>::p_merge(const Assoc& w, const Assoc& w_p) const
1602 {
1603  //KJB(ASSERT(w.is_valid(m_v_bar, m_d_bar, m_noise_sigma, m_convert)));
1604  //KJB(ASSERT(w_p.is_valid(m_v_bar, m_d_bar, m_noise_sigma, m_convert)));
1605 
1606 #ifdef TEST
1607  if(!is_valid_merge(w, w_p))
1608  {
1609  std::cout << "WARNING: invalid merge." << std::endl;
1610  return std::log(0.0);
1611  }
1612 #endif
1613 
1614  size_t nmp = merge_info.num_merge_pairs;
1615  if(nmp == 0)
1616  {
1617  return negative_infinity();
1618  }
1619 
1620  return -std::log(nmp);
1621 }
1622 
1623 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
1624 
1625 template <class Track>
1628  const Assoc& w,
1629  const Assoc& w_p
1630 ) const
1631 {
1632  //KJB(ASSERT(w.is_valid(m_v_bar, m_d_bar, m_noise_sigma, m_convert)));
1633  //KJB(ASSERT(w_p.is_valid(m_v_bar, m_d_bar, m_noise_sigma, m_convert)));
1634 
1635 #ifdef TEST
1636  if(!is_valid_extension(w, w_p))
1637  {
1638  std::cout << "WARNING: invalid extension." << std::endl;
1639  return std::log(0.0);
1640  }
1641 #endif
1642 
1643  double p = -std::log(2 * extension_info.num_tracks);
1644 
1645  // probability of growing to what it is
1646  if(extension_info.direction == GROW_FORWARD)
1647  {
1648  p += p_grow_track_forward(*extension_info.extended_track_p, w,
1649  extension_info.previous_end);
1650  }
1651  else
1652  {
1653  p += p_grow_track_backward(*extension_info.extended_track_p, w,
1654  extension_info.previous_end);
1655  }
1656 
1657  return p;
1658 }
1659 
1660 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
1661 
1662 template <class Track>
1665  const Assoc& w,
1666  const Assoc& w_p
1667 ) const
1668 {
1669  //KJB(ASSERT(w.is_valid(m_v_bar, m_d_bar, m_noise_sigma, m_convert)));
1670  //KJB(ASSERT(w_p.is_valid(m_v_bar, m_d_bar, m_noise_sigma, m_convert)));
1671 
1672 #ifdef TEST
1673  if(!is_valid_reduction(w, w_p))
1674  {
1675  std::cout << "WARNING: invalid reduction." << std::endl;
1676  return std::log(0.0);
1677  }
1678 #endif
1679 
1680  return std::log(0.5 / (reduction_info.num_tracks
1681  * (reduction_info.reduced_track_size - 2)));
1682 }
1683 
1684 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
1685 
1686 template <class Track>
1687 double Proposer<Track>::p_switch(const Assoc& w, const Assoc& w_p) const
1688 {
1689  //KJB(ASSERT(w.is_valid(m_v_bar, m_d_bar, m_noise_sigma, m_convert)));
1690  //KJB(ASSERT(w_p.is_valid(m_v_bar, m_d_bar, m_noise_sigma, m_convert)));
1691 
1692 #ifdef TEST
1693  if(!is_valid_switch(w, w_p))
1694  {
1695  std::cout << "WARNING: invalid switch." << std::endl;
1696  return std::log(0.0);
1697  }
1698 #endif
1699 
1700  size_t nsp = switch_info.num_switch_points;
1701  if(nsp == 0)
1702  {
1703  return negative_infinity();
1704  }
1705 
1706  return -std::log(nsp);
1707 }
1708 
1709 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
1710 
1711 template <class Track>
1714  const Assoc& w,
1715  const Assoc& w_p
1716 ) const
1717 {
1718  //KJB(ASSERT(w.is_valid(m_v_bar, m_d_bar, m_noise_sigma, m_convert)));
1719  //KJB(ASSERT(w_p.is_valid(m_v_bar, m_d_bar, m_noise_sigma, m_convert)));
1720 
1721 #ifdef TEST
1722  if(!is_valid_secretion(w, w_p))
1723  {
1724  std::cout << "WARNING: invalid secretion." << std::endl;
1725  return std::log(0.0);
1726  }
1727 #endif
1728 
1729  if(secretion_info.num_valid_tracks == 0)
1730  {
1731  return negative_infinity();
1732  }
1733 
1734  double p = -std::log(secretion_info.num_valid_tracks);
1735 
1736  p += (secretion_info.window_size)*std::log(0.5);
1737 
1738  return p;
1739 }
1740 
1741 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
1742 
1743 template <class Track>
1746  const Assoc& w,
1747  const Assoc& w_p
1748 ) const
1749 {
1750  //KJB(ASSERT(w.is_valid(m_v_bar, m_d_bar, m_noise_sigma, m_convert)));
1751  //KJB(ASSERT(w_p.is_valid(m_v_bar, m_d_bar, m_noise_sigma, m_convert)));
1752 
1753 #ifdef TEST
1754  if(!is_valid_absorption(w, w_p))
1755  {
1756  std::cout << "WARNING: invalid absorption." << std::endl;
1757  return std::log(0.0);
1758  }
1759 #endif
1760 
1761  return -std::log(absorption_info.num_valid_track_pairs);
1762 }
1763 
1764 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
1765 
1766 template <class Track>
1767 inline
1769 (
1770  const Assoc& w,
1771  const Assoc& w_p
1772 ) const
1773 {
1774  if(w_p.empty() || w_p.size() - w.size() != 1)
1775  {
1776  return false;
1777  }
1778 
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())
1782  {
1783  return true;
1784  }
1785 
1786  diff_p.second++;
1787  diff_p = std::mismatch(diff_p.first, w.end(), diff_p.second);
1788  if(diff_p.first == w.end())
1789  {
1790  return true;
1791  }
1792 
1793  return false;
1794 }
1795 
1796 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
1797 
1798 template <class Track>
1799 inline
1800 bool Proposer<Track>::is_valid_death
1801 (
1802  const Assoc& w,
1803  const Assoc& w_p
1804 ) const
1805 {
1806  return is_valid_birth(w_p, w);
1807 }
1808 
1809 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
1810 
1811 template <class Track>
1812 inline
1813 bool Proposer<Track>::is_valid_split
1814 (
1815  const Assoc& w,
1816  const Assoc& w_p
1817 ) const
1818 {
1819  if(w.empty() || w_p.size() - w.size() != 1)
1820  {
1821  return false;
1822  }
1823 
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>());
1828 
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>());
1833 
1834  if(w_minus_w_p.size() != 1 || w_p_minus_w.size() != 2)
1835  {
1836  return false;
1837  }
1838 
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());
1842 
1843  if(*w_minus_w_p.begin() != merged_track)
1844  {
1845  return false;
1846  }
1847 
1848  return true;
1849 }
1850 
1851 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
1852 
1853 template <class Track>
1854 inline
1855 bool Proposer<Track>::is_valid_merge
1856 (
1857  const Assoc& w,
1858  const Assoc& w_p
1859 ) const
1860 {
1861  return is_valid_split(w_p, w);
1862 }
1863 
1864 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
1865 
1866 template <class Track>
1867 bool Proposer<Track>::is_valid_extension
1868 (
1869  const Assoc& w,
1870  const Assoc& w_p
1871 ) const
1872 {
1873  if(w.empty() || w.size() != w_p.size())
1874  {
1875  return false;
1876  }
1877 
1878  // Find original track
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)
1883  {
1884  return false;
1885  }
1886 
1887  const Track& original_track = *w.find(*diff.begin());
1888 
1889  // Find extended track
1890  diff.clear();
1891  std::set_difference(w_p.begin(), w_p.end(), w.begin(), w.end(),
1892  std::inserter(diff, diff.begin()), Compare_tracks<Track>());
1893 
1894  if(diff.size() != 1)
1895  {
1896  return false;
1897  }
1898 
1899  Assoc_const_iterator extended_track_p = w_p.find(*diff.begin());
1900 
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())
1905  {
1906  return true;
1907  }
1908 
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())
1914  {
1915  return true;
1916  }
1917 
1918  return false;
1919 }
1920 
1921 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
1922 
1923 template <class Track>
1924 inline
1925 bool Proposer<Track>::is_valid_reduction
1926 (
1927  const Assoc& w,
1928  const Assoc& w_p
1929 ) const
1930 {
1931  return is_valid_extension(w_p, w);
1932 }
1933 
1934 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
1935 
1936 template <class Track>
1937 bool Proposer<Track>::is_valid_switch
1938 (
1939  const Assoc& w,
1940  const Assoc& w_p
1941 ) const
1942 {
1943  if(w.size() < 2 || w.size() != w_p.size())
1944  {
1945  return false;
1946  }
1947 
1948  // Find original tracks
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)
1953  {
1954  return false;
1955  }
1956 
1957  const Track& original_track_1 = *w.find(*diff.begin());
1958  const Track& original_track_2 = *w.find(*diff.rbegin());
1959 
1960  // Find swapped tracks
1961  diff.clear();
1962  std::set_difference(w_p.begin(), w_p.end(), w.begin(), w.end(),
1963  std::inserter(diff, diff.begin()), Compare_tracks<Track>());
1964 
1965  if(diff.size() != 2)
1966  {
1967  return false;
1968  }
1969 
1970  const Track& swapped_track_1 = *w_p.find(*diff.begin());
1971  const Track& swapped_track_2 = *w_p.find(*diff.rbegin());
1972 
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());
1976 
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());
1980 
1981  if(track1_break.first == original_track_1.end())
1982  {
1983  return false;
1984  }
1985 
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))
1994  {
1995  return false;
1996  }
1997 
1998  return true;
1999 }
2000 
2001 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
2002 
2003 template <class Track>
2004 inline
2005 bool Proposer<Track>::is_valid_secretion
2006 (
2007  const Assoc& w,
2008  const Assoc& w_p
2009 ) const
2010 {
2011  if(w.empty() || w_p.size() - w.size() != 1)
2012  {
2013  return false;
2014  }
2015 
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>());
2020 
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>());
2025 
2026  if(w_minus_w_p.size() != 1 || w_p_minus_w.size() != 2)
2027  {
2028  return false;
2029  }
2030 
2031  /*Assoc_const_iterator orig_track_p = w_minus_w_p.begin();
2032  Assoc_const_iterator new_track_p1 = w_p_minus_w.begin();
2033  Assoc_const_iterator new_track_p2 = w_p_minus_w.end();
2034  new_track_p2--;
2035 
2036  Track merged_track = *new_track_p1;
2037  merged_track.insert(new_track_p2->begin(), new_track_p2->end());
2038 
2039  if(*orig_track_p != merged_track)
2040  {
2041  return false;
2042  }
2043 
2044  if(orig_track_p->begin()->second != new_track_p1->begin()->second
2045  && orig_track_p->begin()->second != new_track_p2->begin()->second)
2046  {
2047  return false;
2048  }*/
2049 
2050  return true;
2051 }
2052 
2053 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
2054 
2055 template <class Track>
2056 inline
2057 bool Proposer<Track>::is_valid_absorption
2058 (
2059  const Assoc& w,
2060  const Assoc& w_p
2061 ) const
2062 {
2063  return is_valid_secretion(w_p, w);
2064 }
2065 
2066 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
2067 
2068 template <class Track>
2071  Track& track,
2072  const Assoc& w
2073 ) const
2074 {
2075  int t = track.get_end_time() + 1;
2076 
2077  if(t >= w.get_data().size()) return;
2078 
2079  for(;;)
2080  {
2081  int prev_t = track.get_end_time();
2082  int d = t - prev_t;
2083  if(d > m_b_bar) break;
2084 
2085  std::pair<Vector, Vector> fut = track_future_ls(track);
2086  const Vector& x = fut.first;
2087  const Vector& vel = fut.second;
2088 
2089  // get candidates and their probabilities
2090  //std::set<const Element*> dpt = w.dead_neighbors(x, prev_t, d,
2091  // m_b_bar, m_v_bar, m_noise_sigma, m_convert);
2092  std::set<const Element*> dpt = w.get_dead_points_at_time(t);
2093  Probability_map probs =
2094  get_growing_probabilities(track, prev_t, dpt, x, vel, d);
2095 
2096  size_t num_added = 0;
2097  BOOST_FOREACH(const typename Probability_map::left_map::value_type& pr,
2098  probs.left)
2099  {
2100  double p_det = std::exp(pr.first);
2101  if(p_det > 0.15)
2102  {
2103  if(kjb::sample(Uniform_distribution()) < p_det)
2104  {
2105  track.insert(std::make_pair(t, pr.second));
2106  num_added++;
2107  }
2108  }
2109  }
2110 
2111  // if nothing sampled
2112  if(num_added == 0)
2113  {
2114  // stop with probability gamma
2115  if(kjb::sample(Uniform_distribution()) < m_gamma) break;
2116  }
2117 
2118  if(t == w.get_data().size()) break;
2119  t++;
2120  }
2121 }
2122 
2123 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
2124 
2125 template <class Track>
2128  Track& track,
2129  const Assoc& w
2130 ) const
2131 {
2132  int t = track.get_start_time() - 1;
2133 
2134  if(t <= 1) return;
2135 
2136  for(;;)
2137  {
2138  int prev_t = track.get_start_time();
2139  int d = t - prev_t;
2140  if(std::abs(d) > m_b_bar) break;
2141 
2142  std::pair<Vector, Vector> fut = track_past_ls(track);
2143  const Vector& x = fut.first;
2144  const Vector& vel = fut.second;
2145 
2146  // get candidates and their probabilities
2147  //std::set<const Element*> dpt = w.dead_neighbors(x, prev_t, d,
2148  // m_b_bar, m_v_bar, m_noise_sigma, m_convert);
2149  std::set<const Element*> dpt = w.get_dead_points_at_time(t);
2150  Probability_map probs =
2151  get_growing_probabilities(track, prev_t, dpt, x, vel, d);
2152 
2153  size_t num_added = 0;
2154  BOOST_FOREACH(const typename Probability_map::left_map::value_type& pr,
2155  probs.left)
2156  {
2157  double p_det = std::exp(pr.first);
2158 
2159  if(kjb::sample(Uniform_distribution()) < p_det)
2160  {
2161  track.insert(std::make_pair(t, pr.second));
2162  num_added++;
2163  }
2164  }
2165 
2166  // if nothing sampled
2167  if(num_added == 0)
2168  {
2169  // stop with probability gamma
2170  if(kjb::sample(Uniform_distribution()) < m_gamma) break;
2171  }
2172 
2173  if(t == 1) break;
2174  t--;
2175  }
2176 }
2177 
2178 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
2179 
2180 template <class Track>
2183  const Track& track,
2184  const Assoc& w,
2185  int t
2186 ) const
2187 {
2188  double p = 0.0;
2189 
2190  // make sure 't' exists
2191  if(!track.count(t))
2192  {
2193  return std::log(0);
2194  }
2195 
2196  // current end
2197  Track_const_iterator eot = --track.end();
2198  if(eot->first == t)
2199  {
2200  return negative_infinity();
2201  }
2202 
2203  Track_const_iterator prev_p = track.upper_bound(t);
2204  Track_const_iterator next_p = prev_p--;
2205  t++;
2206 
2207  while(prev_p != eot)
2208  {
2209  int prev_t = prev_p->first;
2210 
2211  if(t - prev_t > m_d_bar)
2212  {
2213  return negative_infinity();
2214  }
2215 
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;
2219 
2220  std::set<const Element*> dpt = w.get_dead_points_at_time(t);
2221  Probability_map probs = get_growing_probabilities(track, prev_t,
2222  dpt, prev_x, vel, t - prev_t);
2223 
2224  // start by assuming none got added...
2225  p += p_choose_nothing(probs);
2226  if(t == next_p->first)
2227  {
2228  // for each added correct by...
2229  while(next_p != track.end() && next_p->first == t)
2230  {
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())
2235  {
2236  std::cout << "WARNING: this should never happen!\n";
2237  }
2238  else
2239  {
2240  // ...multiply by prob of choosing and divide by
2241  // prob of not choosing
2242  if(p != negative_infinity())
2243  {
2244  p += (pr_p->second) - log(1.0 - exp(pr_p->second));
2245  }
2246  }
2247  next_p++;
2248  }
2249 
2250  prev_p = next_p;
2251  prev_p--;
2252  }
2253  else
2254  {
2255  p += std::log(1.0 - m_gamma);
2256  }
2257 
2258  t++;
2259  }
2260 
2261  if(t > w.get_data().size())
2262  {
2263  return p;
2264  }
2265 
2266  // compute probability of stopping; add probabliities of stopping
2267  // in different ways
2268  double p_last = 0.0;
2269 
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;
2274  //std::set<const Element*> dpt = w.dead_neighbors(px, pt, 1,
2275  // m_d_bar, m_v_bar, m_noise_sigma, m_convert);
2276  std::set<const Element*> dpt = w.get_dead_points_at_time(pt + 1);
2277  Probability_map probs =
2278  get_growing_probabilities(track, pt, dpt, px, vel, 1);
2279 
2280  double a_t_1 = std::exp(p_choose_nothing(probs)) * m_gamma;
2281  p_last = a_t_1;
2282 
2283  t++;
2284  int last_t = std::min(pt + m_d_bar, (int)w.get_data().size());
2285  for(; t <= last_t; t++)
2286  {
2287  //dpt = w.dead_neighbors(px, pt, t - pt,
2288  // m_d_bar, m_v_bar, m_noise_sigma, m_convert);
2289  dpt = w.get_dead_points_at_time(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));
2293  p_last += a_t;
2294  a_t_1 = a_t;
2295  }
2296 
2297  p += std::log(p_last);
2298 
2299  return p;
2300 }
2301 
2302 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
2303 
2304 template <class Track>
2307  const Track& track,
2308  const Assoc& w,
2309  int t
2310 ) const
2311 {
2312  double p = 0.0;
2313 
2314  // make sure 't' exists
2315  if(!track.count(t))
2316  {
2317  return std::log(0.0);
2318  }
2319 
2320  // current end
2321  Track_const_riterator eot = --track.rend();
2322  if(eot->first == t)
2323  {
2324  return negative_infinity();
2325  }
2326 
2327  Track_const_riterator prev_p(track.lower_bound(t));
2328  Track_const_riterator next_p = prev_p--;
2329  t--;
2330 
2331  while(prev_p != eot)
2332  {
2333  int prev_t = prev_p->first;
2334 
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;
2338 
2339  std::set<const Element*> dpt = w.get_dead_points_at_time(t);
2340  Probability_map probs = get_growing_probabilities(track, prev_t, dpt,
2341  prev_x, vel, t - prev_t);
2342 
2343  // start by assuming none got added
2344  p += p_choose_nothing(probs);
2345  if(t == next_p->first)
2346  {
2347  // for each added correct by...
2348  while(next_p != track.rend() && next_p->first == t)
2349  {
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())
2354  {
2355  std::cout << "WARNING: this should never happen!\n";
2356  }
2357  else
2358  {
2359  // ...multiply by prob of choosing and divide by
2360  // prob of not choosing
2361  if(p != negative_infinity())
2362  {
2363  p += (pr_p->second) - log(1.0 - exp(pr_p->second));
2364  }
2365  }
2366  next_p++;
2367  }
2368 
2369  prev_p = next_p;
2370  prev_p--;
2371  }
2372  else
2373  {
2374  p += std::log(1.0 - m_gamma);
2375  }
2376 
2377  t--;
2378  }
2379 
2380  if(t < 1)
2381  {
2382  return p;
2383  }
2384 
2385  // compute probability of stopping; add probabliities of stopping
2386  // in different ways
2387  double p_last = 0.0;
2388 
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;
2393  //std::set<const Element*> dpt = w.dead_neighbors(px, pt, -1,
2394  // m_d_bar, m_v_bar, m_noise_sigma, m_convert);
2395  std::set<const Element*> dpt = w.get_dead_points_at_time(pt - 1);
2396  Probability_map probs =
2397  get_growing_probabilities(track, pt, dpt, px, vel, -1);
2398 
2399  double a_t_1 = std::exp(p_choose_nothing(probs)) * m_gamma;
2400  p_last = a_t_1;
2401 
2402  t--;
2403  int last_t = std::max(pt - m_d_bar, 1);
2404  for(; t >= last_t; t--)
2405  {
2406  //dpt = w.dead_neighbors(px, pt, t - pt,
2407  // m_d_bar, m_v_bar, m_noise_sigma, m_convert);
2408  dpt = w.get_dead_points_at_time(t);
2409  probs =
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));
2413  p_last += a_t;
2414  a_t_1 = a_t;
2415  }
2416 
2417  p += std::log(p_last);
2418 
2419  return p;
2420 }
2421 
2422 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
2423 
2424 template <class Track>
2427  const Probability_map& probabilities
2428 ) const
2429 {
2430  double p = 0.0;
2431  BOOST_FOREACH(const typename Probability_map::left_map::value_type& pr,
2432  probabilities.left)
2433  {
2434  p += std::log(1.0 - std::exp(pr.first));
2435  }
2436 
2437  return p;
2438 }
2439 
2440 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
2441 
2442 template <class Track>
2443 typename Proposer<Track>::Probability_map
2446  const Track& track,
2447  int t,
2448  std::set<const Element*>& candidates,
2449  const Vector& x,
2450  const Vector& vel,
2451  int d
2452 ) const
2453 {
2454  Probability_map probabilities;
2455  int dist = std::abs(d);
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);
2459  Normal_distribution N_u(0.0, sg_u);
2460  Normal_distribution N_v(0.0, sg_v);
2461 
2462  Vector v_x = x;
2463  if(!vel.empty()) v_x += dist*vel;
2464 
2465  double noise_score = get_velocity_noise_score(d, vel);
2466 
2467  BOOST_FOREACH(const Element* e_p, candidates)
2468  {
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);
2472 
2473  // probabilities of features
2474  if(m_feature_prob)
2475  {
2476  std::vector<double> f_ps
2477  = m_feature_prob(&track, t, e_p, t + d, m_b_bar);
2478 
2479  // Take the geometric mean
2480  double pp = 1.0;
2481  BOOST_FOREACH(double f_p, f_ps)
2482  {
2483  pp *= f_p;
2484  }
2485 
2486  p = std::pow(p * pp, 1.0/(1 + f_ps.size()));
2487  }
2488 
2489  typedef typename Probability_map::left_map::value_type value_type;
2490  double prob = 0.0;
2491  if(p == 0.0)
2492  {
2493  prob = negative_infinity();
2494  }
2495  else
2496  {
2497  prob = std::log(p);
2498  }
2499 
2500  probabilities.left.insert(value_type(prob, e_p));
2501  }
2502 
2503  return probabilities;
2504 }
2505 
2506 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
2507 
2508 template <class Track>
2509 std::pair<Vector, Vector> Proposer<Track>::track_future
2511  const Track& track,
2512  int t
2513 ) const
2514 {
2515  // constants
2516  const size_t wsz = m_b_bar / 4;
2517  const double vel_thresh = 0.0;
2518 
2519  Track_const_riterator pair_p;
2520  if(t < 0)
2521  {
2522  pair_p = track.rbegin();
2523  t = pair_p->first;
2524  }
2525  else
2526  {
2527  //pair_p = Track_const_riterator(track.lower_bound(t));
2528  //pair_p--;
2529  pair_p = Track_const_riterator(track.upper_bound(t));
2530  IFT(pair_p->first == t, Illegal_argument,
2531  "The frame passed to track_future() must be valid.");
2532  }
2533 
2534  // average final points
2535  std::vector<const Element*> cur_elems;
2536  while(pair_p != track.rend() && pair_p->first == t)
2537  {
2538  cur_elems.push_back(pair_p->second);
2539  pair_p++;
2540  }
2541 
2542  Element x = m_average(cur_elems);
2543  if(pair_p == track.rend())
2544  {
2545  return std::make_pair(m_convert(x), Vector());
2546  }
2547 
2548  // compute direction in window
2549  size_t wnd = t - pair_p->first;
2550  while(wnd < wsz)
2551  {
2552  pair_p++;
2553  if(pair_p == track.rend()) break;
2554  wnd = t - pair_p->first;
2555  }
2556 
2557  Vector vel;
2558  if(wnd >= wsz)
2559  {
2560  std::vector<const Element*> cur_elems;
2561  while(pair_p != track.rend() && pair_p->first == t - wnd)
2562  {
2563  cur_elems.push_back(pair_p->second);
2564  pair_p++;
2565  }
2566 
2567  Element from = m_average(cur_elems);
2568 
2569  vel = m_convert(x) - m_convert(from);
2570  vel /= wnd;
2571 
2572  if(vel.magnitude() <= vel_thresh)
2573  {
2574  vel.set(0.0, 0.0);
2575  }
2576  }
2577 
2578  return std::make_pair(m_convert(x), vel);
2579 }
2580 
2581 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
2582 
2583 template <class Track>
2584 std::pair<Vector, Vector> Proposer<Track>::track_past
2586  const Track& track,
2587  int t
2588 ) const
2589 {
2590  // constants
2591  const size_t wsz = m_b_bar / 4;
2592  const double vel_thresh = 0.0;
2593 
2594  Track_const_iterator pair_p;
2595  if(t < 0)
2596  {
2597  pair_p = track.begin();
2598  t = pair_p->first;
2599  }
2600  else
2601  {
2602  //pair_p = track.upper_bound(t);
2603  //pair_p--;
2604  pair_p = track.lower_bound(t);
2605  IFT(pair_p->first == t, Illegal_argument,
2606  "The frame passed to track_past() must be valid.");
2607  }
2608 
2609  // average final points
2610  std::vector<const Element*> cur_elems;
2611  while(pair_p != track.end() && pair_p->first == t)
2612  {
2613  cur_elems.push_back(pair_p->second);
2614  pair_p++;
2615  }
2616 
2617  Element x = m_average(cur_elems);
2618  if(pair_p == track.end())
2619  {
2620  return std::make_pair(m_convert(x), Vector());
2621  }
2622 
2623  // compute direction in window
2624  size_t wnd = pair_p->first - t;
2625  while(wnd < wsz)
2626  {
2627  pair_p++;
2628  if(pair_p == track.end()) break;
2629  wnd = pair_p->first - t;
2630  }
2631 
2632  Vector vel;
2633  if(wnd >= wsz)
2634  {
2635  std::vector<const Element*> cur_elems;
2636  while(pair_p != track.end() && pair_p->first == t + wnd)
2637  {
2638  cur_elems.push_back(pair_p->second);
2639  pair_p++;
2640  }
2641 
2642  Element from = m_average(cur_elems);
2643 
2644  vel = m_convert(x) - m_convert(from);
2645  vel /= wnd;
2646 
2647  if(vel.magnitude() <= vel_thresh)
2648  {
2649  vel.set(0.0, 0.0);
2650  }
2651  }
2652 
2653  return std::make_pair(m_convert(x), vel);
2654 }
2655 
2656 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
2657 
2658 template <class Track>
2659 std::pair<Vector, Vector> Proposer<Track>::track_future_ls
2661  const Track& track,
2662  int t
2663 ) const
2664 {
2665  // constants
2666  const size_t wsz = m_b_bar;
2667  const double vel_thresh = 0.0;
2668 
2669  Track_const_riterator pair_p;
2670  if(t < 0)
2671  {
2672  pair_p = track.rbegin();
2673  t = pair_p->first;
2674  }
2675  else
2676  {
2677  //pair_p = Track_const_riterator(track.lower_bound(t));
2678  //pair_p--;
2679  pair_p = Track_const_riterator(track.upper_bound(t));
2680  IFT(pair_p->first == t, Illegal_argument,
2681  "The frame passed to track_future() must be valid.");
2682  }
2683 
2684  double x_N = m_convert(*pair_p->second)[0];
2685  double y_N = m_convert(*pair_p->second)[1];
2686  double sum_x = x_N;
2687  double sum_y = y_N;
2688  double sum_x_squared = x_N*x_N;
2689  double sum_xy = x_N*y_N;
2690  size_t N = 1;
2691  double x_1;
2692  double y_1;
2693 
2694  // average final points
2695  std::vector<const Element*> cur_elems;
2696  size_t wnd = t - pair_p->first;
2697  while(wnd < wsz)
2698  {
2699  if(pair_p->first == t)
2700  {
2701  cur_elems.push_back(pair_p->second);
2702  }
2703  pair_p++;
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;
2708  sum_x += x_i;
2709  sum_y += y_i;
2710  sum_x_squared += (x_i * x_i);
2711  sum_xy += (x_i * y_i);
2712  x_1 = x_i;
2713  y_1 = y_i;
2714  N++;
2715  }
2716 
2717  if(wnd >= wsz)
2718  {
2719  Vector x;
2720  Vector vel;
2721 
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);
2725  if(den != 0)
2726  {
2727  double m = (sum_xy - (sum_x*sum_y)/N) / den;
2728  double b = sum_y / N - m * sum_x / N;
2729 
2730  double f_1 = m * x_1 + b;
2731  double f_N = m * p_N[0] + b;
2732 
2733  x.set(p_N[0], f_N);
2734  vel.set(p_N[0] - x_1, f_N - f_1);
2735  vel /= wnd;
2736  if(vel.magnitude() <= vel_thresh)
2737  {
2738  vel.set(0.0, 0.0);
2739  }
2740  }
2741  else
2742  {
2743  x = p_N;
2744  vel.set(0.0, p_N[1] - y_1);
2745  }
2746 
2747  return std::make_pair(x, vel);
2748  }
2749 
2750  // if we can't fit a line we just
2751  // do it the bad way
2752  return track_future(track, t);
2753 }
2754 
2755 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
2756 
2757 template <class Track>
2758 std::pair<Vector, Vector> Proposer<Track>::track_past_ls
2760  const Track& track,
2761  int t
2762 ) const
2763 {
2764  // constants
2765  const size_t wsz = m_b_bar;
2766  const double vel_thresh = 0.0;
2767 
2768  Track_const_iterator pair_p;
2769  if(t < 0)
2770  {
2771  pair_p = track.begin();
2772  t = pair_p->first;
2773  }
2774  else
2775  {
2776  //pair_p = track.upper_bound(t);
2777  //pair_p--;
2778  pair_p = track.lower_bound(t);
2779  IFT(pair_p->first == t, Illegal_argument,
2780  "The frame passed to track_future() must be valid.");
2781  }
2782 
2783  double x_N = m_convert(*pair_p->second)[0];
2784  double y_N = m_convert(*pair_p->second)[1];
2785  double sum_x = x_N;
2786  double sum_y = y_N;
2787  double sum_x_squared = x_N*x_N;
2788  double sum_xy = x_N*y_N;
2789  size_t N = 1;
2790  double x_1;
2791  double y_1;
2792 
2793  // average the final points
2794  std::vector<const Element*> cur_elems;
2795 
2796  size_t wnd = pair_p->first - t;
2797  while(wnd < wsz)
2798  {
2799  if(pair_p->first == t)
2800  {
2801  cur_elems.push_back(pair_p->second);
2802  }
2803  pair_p++;
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;
2808  sum_x += x_i;
2809  sum_y += y_i;
2810  sum_x_squared += (x_i * x_i);
2811  sum_xy += (x_i * y_i);
2812  x_1 = x_i;
2813  y_1 = y_i;
2814  N++;
2815  }
2816 
2817  if(wnd >= wsz)
2818  {
2819  Vector x;
2820  Vector vel;
2821 
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);
2825  if(den != 0)
2826  {
2827  double m = (sum_xy - (sum_x*sum_y)/N) / den;
2828  double b = sum_y / N - m * sum_x / N;
2829 
2830  double f_1 = m * x_1 + b;
2831  double f_N = m * p_N[0] + b;
2832 
2833  x.set(p_N[0], f_N);
2834  vel.set(p_N[0] - x_1, f_N - f_1);
2835  vel /= wnd;
2836  if(vel.magnitude() <= vel_thresh)
2837  {
2838  vel.set(0.0, 0.0);
2839  }
2840  }
2841  else
2842  {
2843  x = p_N;
2844  vel.set(0.0, p_N[1] - y_1);
2845  }
2846 
2847  return std::make_pair(x, vel);
2848  }
2849 
2850  // if we can't fit a line we just
2851  // do it the bad way
2852  return track_future(track, t);
2853 }
2854 
2855 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
2856 
2857 template <class Track>
2858 std::pair<Vector, Vector> Proposer<Track>::track_future_tls
2860  const Track& track,
2861  int t
2862 ) const
2863 {
2864  // constants
2865  const size_t wsz = m_b_bar;
2866  const double vel_thresh = 0.0;
2867 
2868  Track_const_riterator pair_p;
2869  if(t < 0)
2870  {
2871  pair_p = track.rbegin();
2872  t = pair_p->first;
2873  }
2874  else
2875  {
2876  //pair_p = Track_const_riterator(track.lower_bound(t));
2877  //pair_p--;
2878  pair_p = Track_const_riterator(track.upper_bound(t));
2879  IFT(pair_p->first == t, Illegal_argument,
2880  "The frame passed to track_future() must be valid.");
2881  }
2882 
2883  Vector x = Vector().set(1, m_convert(*pair_p->second)[0]);;
2884  Vector y = Vector().set(1, m_convert(*pair_p->second)[1]);;
2885  size_t wnd = t - pair_p->first;
2886  while(wnd < wsz)
2887  {
2888  pair_p++;
2889  if(pair_p == track.rend()) break;
2890 
2891  wnd = t - pair_p->first;
2892 
2893  x.push_back(m_convert(*pair_p->second)[0]);
2894  y.push_back(m_convert(*pair_p->second)[1]);
2895  }
2896 
2897  if(wnd >= wsz)
2898  {
2899  double x_1 = x.back();
2900  double x_n = x.front();
2901  double x_bar = mean(x.begin(), x.end());
2902  double y_bar = mean(y.begin(), y.end());
2903  x -= Vector(x.get_length(), x_bar);
2904  y -= Vector(y.get_length(), y_bar);
2905 
2906  Matrix A(x.get_length(), 2);
2907  A.set_col(0, x);
2908  A.set_col(1, y);
2909  Matrix V;
2910  Vector d;
2911  diagonalize(matrix_transpose(A)*A, V, d, true);
2912 
2913  double a = V(0, 1);
2914  double b = V(1, 1);
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;
2918 
2919  Vector vel = Vector().set(x_n - x_1, y_n - y_1) / wnd;
2920  Vector x_s = Vector().set(x_n, y_n);
2921  if(vel.magnitude() <= vel_thresh)
2922  {
2923  vel.set(0.0, 0.0);
2924  }
2925 
2926  return std::make_pair(x_s, vel);
2927  }
2928 
2929  // if we can't fit a line we just
2930  // do it the bad way
2931  return track_future(track, t);
2932 }
2933 
2934 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
2935 
2936 template <class Track>
2937 std::pair<Vector, Vector> Proposer<Track>::track_past_tls
2939  const Track& track,
2940  int t
2941 ) const
2942 {
2943  // constants
2944  const size_t wsz = m_b_bar;
2945  const double vel_thresh = 0.0;
2946 
2947  Track_const_iterator pair_p;
2948  if(t < 0)
2949  {
2950  pair_p = track.begin();
2951  t = pair_p->first;
2952  }
2953  else
2954  {
2955  //pair_p = track.upper_bound(t);
2956  //pair_p--;
2957  pair_p = track.lower_bound(t);
2958  IFT(pair_p->first == t, Illegal_argument,
2959  "The frame passed to track_future() must be valid.");
2960  }
2961 
2962  Vector x = Vector().set(1, m_convert(*pair_p->second)[0]);;
2963  Vector y = Vector().set(1, m_convert(*pair_p->second)[1]);;
2964  size_t wnd = pair_p->first - t;
2965  while(wnd < wsz)
2966  {
2967  pair_p++;
2968  if(pair_p == track.end()) break;
2969 
2970  wnd = pair_p->first - t;
2971 
2972  x.push_back(m_convert(*pair_p->second)[0]);
2973  y.push_back(m_convert(*pair_p->second)[1]);
2974  }
2975 
2976  if(wnd >= wsz)
2977  {
2978  double x_1 = x.back();
2979  double x_n = x.front();
2980  double x_bar = mean(x.begin(), x.end());
2981  double y_bar = mean(y.begin(), y.end());
2982  x -= Vector(x.get_length(), x_bar);
2983  y -= Vector(y.get_length(), y_bar);
2984 
2985  Matrix A(x.get_length(), 2);
2986  A.set_col(0, x);
2987  A.set_col(1, y);
2988  Matrix V;
2989  Vector d;
2990  diagonalize(matrix_transpose(A)*A, V, d, true);
2991 
2992  double a = V(0, 1);
2993  double b = V(1, 1);
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;
2997 
2998  Vector vel = Vector().set(x_n - x_1, y_n - y_1) / wnd;
2999  Vector x_s = Vector().set(x_n, y_n);
3000  if(vel.magnitude() <= vel_thresh)
3001  {
3002  vel.set(0.0, 0.0);
3003  }
3004 
3005  return std::make_pair(x_s, vel);
3006  }
3007 
3008  // if we can't fit a line we just
3009  // do it the bad way
3010  return track_future(track, t);
3011 }
3012 
3013 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
3014 
3015 template <class Track>
3016 size_t Proposer<Track>::count_split_points(const Assoc& w) const
3017 {
3018  size_t nsp = 0;
3019  BOOST_FOREACH(const Track& track, w)
3020  {
3021  size_t track_sz = track.real_size();
3022  if(track_sz >= 4)
3023  {
3024  nsp += (track_sz - 3);
3025  }
3026  }
3027 
3028  return nsp;
3029 }
3030 
3031 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
3032 
3033 template <class Track>
3034 size_t Proposer<Track>::count_merge_pairs(const Assoc& w) const
3035 {
3036  size_t nmp = 0;
3037  for(Assoc_const_iterator track_p1 = w.begin();
3038  track_p1 != w.end();
3039  track_p1++)
3040  {
3041  int t_f = track_p1->get_end_time();
3042  for(Assoc_const_iterator track_p2 = w.begin();
3043  track_p2 != w.end();
3044  track_p2++)
3045  {
3046  int t_1 = track_p2->get_start_time();
3047  if(track_p1 != track_p2 && t_f < t_1)
3048  {
3049  if(is_neighbor(track_p1->get_end_point(),
3050  track_p2->get_start_point(),
3051  t_1 - t_f, m_d_bar, m_v_bar,
3052  m_noise_sigma, m_convert))
3053  {
3054  nmp++;
3055  }
3056  }
3057  }
3058  }
3059 
3060  return nmp;
3061 }
3062 
3063 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
3064 
3065 template <class Track>
3066 size_t Proposer<Track>::count_secretion_tracks(const Assoc& w) const
3067 {
3068  size_t nst = 0;
3069  BOOST_FOREACH(const Track& track, w)
3070  {
3071  if(track.real_size() >= 4)
3072  {
3073  nst++;
3074  }
3075  else
3076  {
3077  if(track.count(track.get_start_time()) > 1
3078  && track.count(track.get_end_time()) > 1)
3079  {
3080  nst++;
3081  }
3082  }
3083  }
3084 
3085  return nst;
3086 }
3087 
3088 /* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
3089 
3090 template <class Track>
3091 size_t Proposer<Track>::count_absorption_pairs(const Assoc& w) const
3092 {
3093  /*size_t nap = 0;
3094  for(Assoc_const_iterator track_p1 = w.begin();
3095  track_p1 != w.end();
3096  track_p1++)
3097  {
3098  Assoc_const_iterator track_p2 = track_p1;
3099  for(track_p2++; track_p2 != w.end(); track_p2++)
3100  {
3101  Track new_track = *track_p1;
3102  new_track.insert(track_p2->begin(), track_p2->end());
3103  if(new_track.is_valid(m_v_bar, m_d_bar, m_noise_sigma, m_convert))
3104  {
3105  nap++;
3106  }
3107  }
3108  }
3109 
3110  return nap;*/
3111  return (w.size() * (w.size() - 1)) / 2;
3112 }
3113 
3114 }} //namespace kjb::mcmcda
3115 
3116 #endif /*MCMCDA_PROPOSER_H_INCLUDED */
3117 
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