KJB
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
mcmcda_prior.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_PRIOR_H_INCLUDED
20 #define MCMCDA_PRIOR_H_INCLUDED
21 
24 #include <prob_cpp/prob_sample.h>
25 #include <cmath>
26 #include <vector>
27 #include <utility>
28 
29 namespace kjb {
30 namespace mcmcda {
31 
37 template<class Track>
38 class Prior
39 {
40 private:
41  double kappa_;
42  double theta_;
43  double lambda_N_;
44  double lambda_A_;
45  double lambda_O_;
46 
47 public:
58  Prior(double kappa, double theta, double lambda_N, double lambda_A) :
59  kappa_(kappa),
60  theta_(theta),
61  lambda_N_(lambda_N),
62  lambda_A_(lambda_A),
63  lambda_O_(0.1)
64  {}
65 
72  double operator()(const Association<Track>& w) const;
73 
74  /* ==================================================== *
75  * GETTERS *
76  * -----------------------------------------------------*/
77 
79  double kappa() const { return kappa_; }
80 
82  double theta() const { return theta_; }
83 
85  double lambda_N() const { return lambda_N_; }
86 
88  double lambda_A() const { return lambda_A_; }
89 
91  double lambda_O() const { return lambda_O_; }
92 
93  /* ==================================================== *
94  * SETTERS *
95  * -----------------------------------------------------*/
96 
98  void set_kappa(double kappa) { kappa_ = kappa; }
99 
101  void set_theta(double theta) { theta_ = theta; }
102 
104  void set_lambda_N(double lambda_N) { lambda_N_ = lambda_N; }
105 
107  void set_lambda_A(double lambda_A) { lambda_A_ = lambda_A; }
108 
110  void set_lambda_O(double lambda_O) { lambda_O_ = lambda_O; }
111 };
112 
113 /*============================================================================*
114  * MEMBER FUNCTION DEFINITIONS *
115  *----------------------------------------------------------------------------*/
116 
117 template<class Track>
119 {
120  using std::log;
121 
122  size_t T = w.get_data().size();
123  size_t m;
124  size_t e;
125  size_t d;
126  size_t a;
127  size_t n;
128  size_t l;
129 
130  get_association_totals(w, m, e, d, a, n, l);
131 
132  double N_t_fac = 0.0;
133  double e_t_fac = 0.0;
134  double n_t_fac = 0.0;
135  double a_it_fac = 0.0;
136  for(size_t t = 1; t <= T; t++)
137  {
138  size_t N_t = w.get_data()[t - 1].size();
139  size_t e_t = 0;
140  size_t a_t = 0;
141  BOOST_FOREACH(const Track& track, w)
142  {
143  if(track.get_start_time() == t) e_t++;
144 
145  size_t a_it = track.count(t);
146  a_it_fac += lgamma(a_it + 1);
147 
148  a_t += a_it;
149  }
150 
151  size_t n_t = N_t - a_t;
152 
153  N_t_fac += lgamma(N_t + 1);
154  e_t_fac += lgamma(e_t + 1);
155  n_t_fac += lgamma(n_t + 1);
156  }
157 
158  // numerator
159  double p = m * (log(kappa_) - lambda_A_);
160  p += (e + d) * log(theta_);
161  p += n * log(lambda_N_);
162  p += a * log(lambda_A_);
163  p += -(kappa_ + (T - 1)*kappa_*theta_ + l*theta_ + T*lambda_N_);
164 
165  // denominator
166  p -= N_t_fac;
167  p -= e_t_fac;
168  p -= n_t_fac;
169  p -= a_it_fac;
170 
171  return p;
172 }
173 
174 /*============================================================================*
175  * NON-MEMBER FUNCTION *
176  *----------------------------------------------------------------------------*/
177 
191 template<class Track>
192 std::pair<std::vector<Track>, std::vector<size_t> > sample
193 (
194  const Prior<Track>& prior,
195  size_t T,
196  const Track& def
197 )
198 {
199  const typename Track::Element* null = 0;
200  IFT(T != 0, Illegal_argument, "Cannot sample from prior; T = 0.");
201 
202  Poisson_distribution P_e1(prior.kappa());
203  Poisson_distribution P_e(prior.kappa() * prior.theta());
204  Poisson_distribution P_l(prior.theta());
205  Poisson_distribution P_a(prior.lambda_A());
206  Poisson_distribution P_N(prior.lambda_N());
207 
208  std::vector<Track> tracks;
209  std::vector<size_t> false_alarms(T);
210  for(size_t t = 1; t <= T; t++)
211  {
212  // new tracks
213  size_t e_t = t == 1 ? kjb::sample(P_e1) : kjb::sample(P_e);
214  for(size_t r = 1; r <= e_t; r++)
215  {
216  // get track length
217  size_t l_r = kjb::sample(P_l);
218 
219  // create track
220  Track track = def;
221  track.insert(std::make_pair(t, null));
222  track.insert(std::make_pair(std::min(t + l_r, T), null));
223 
224  // add it to list of tracks
225  tracks.push_back(track);
226  }
227 
228  // add detections
229  BOOST_FOREACH(Track& track, tracks)
230  {
231  size_t a_rt = kjb::sample(P_a);
232 
233  // all tracks already have 1 detection in first and last frames
234  if(track.get_start_time() == t || track.get_end_time() == t)
235  {
236  if(a_rt != 0) a_rt--;
237  }
238 
239  for(size_t i = 1; i <= a_rt; i++)
240  {
241  track.insert(std::make_pair(t, null));
242  }
243  }
244 
245  // false alarms
246  size_t n_t = kjb::sample(P_N);
247  false_alarms[t - 1] = n_t;
248  }
249 
250  return std::make_pair(tracks, false_alarms);
251 }
252 
262 template<class Track>
263 inline
264 std::pair<std::vector<Track>, std::vector<size_t> > sample
265 (
266  const Prior<Track>& prior,
267  size_t T
268 )
269 {
270  Track def_track;
271  return sample(prior, T, def_track);
272 }
273 
274 /*============================================================================*
275  * OLD PRIORS (KEPT JUST IN CASE) *
276  *----------------------------------------------------------------------------*/
277 
278 /*template<class Track>
279 double Prior<Track>::operator()(const Association<Track>& w) const
280 {
281  int T = w.get_data().size();
282  int M = w.size();
283  int D = 0;
284  int B = 0;
285  int L = 0;
286  int N = 0;
287  int m_1;
288  int n_1;
289  double nblsum = 0.0;
290  double outf = 0.0;
291  for(int t = 1; t <= T; t++)
292  {
293  int b_t = 0;
294  BOOST_FOREACH(const Track& track, w)
295  {
296  if(t != T)
297  {
298  if(track.get_end_time() == t)
299  {
300  D++;
301  }
302  }
303 
304  if(track.get_start_time() == t)
305  {
306  b_t++;
307  }
308 
309  if(track.get_start_time() <= t && track.get_end_time() >= t)
310  {
311  L++;
312  double v_it = 0.9;
313  int a_it = track.count(t);
314  int fa_it = tgamma(a_it + 1);
315  double outf_ti
316  = ((exp(-m_lambda_A)*pow(m_lambda_A, a_it)*v_it)
317  + (exp(-m_lambda_O)*pow(m_lambda_O, a_it)*(1.0 - v_it)))
318  / fa_it;
319 
320  outf += log(outf_ti);
321  }
322  }
323 
324  int m_t = w.count_live_points_at_time(t);
325  if(t == 1)
326  {
327  m_1 = m_t;
328  }
329 
330  int n_t = w.get_data()[t - 1].size() - m_t;
331  N += n_t;
332  if(t == 1)
333  {
334  n_1 = n_t;
335  }
336 
337  if(t >= 2)
338  {
339  B += b_t;
340  nblsum += lgamma(n_t + 1) + lgamma(b_t + 1);
341  }
342  }
343 
344  double num = M*log(m_mu) + (B + D)*log(m_tau) + N*log(m_lambda_N)
345  + -(m_mu + (T - 1)*m_mu*m_tau + L*m_tau + T*m_lambda_N);
346 
347  double den = lgamma(m_1 + 1) + lgamma(n_1 + 1) + nblsum;
348 
349  const double f = 4.0;
350  double reward = 0.0;
351  if(M > 0)
352  {
353  reward = ((double)L/M) * std::log(f);
354  }
355 
356  return num - den + outf + reward;
357 }*/
358 
359 /*template<class Track>
360 double Prior<Track>::operator()(const Association<Track>& w) const
361 {
362  double log_prior = 0.0;
363  int e_t_1 = 0;
364 
365  for(int t = 1; t <= static_cast<int>(w.get_data().size()); t++)
366  {
367  int z_t = 0;
368  int a_t = 0;
369  for(typename Association<Track>::const_iterator track_p = w.begin();
370  track_p != w.end();
371  track_p++)
372  {
373  if(track_p->get_end_time() == t)
374  {
375  z_t++;
376  }
377 
378  if(track_p->get_start_time() == t)
379  {
380  a_t++;
381  }
382  }
383  int c_t = std::max(0, e_t_1 - z_t);
384  int d_t = w.count_live_points_at_time(t);
385  int g_t = c_t + a_t - d_t;
386  int f_t = w.get_data()[t - 1].size() - d_t;
387 
388  using std::log;
389  log_prior += ((z_t * log(m_p_z)) + (c_t * log(1 - m_p_z))
390  + (d_t * log(m_p_d)) + (g_t * log(1 - m_p_d))
391  + (a_t * log(m_lambda_b)) + (f_t * log(m_lambda_f))
392  - lgamma(a_t + 1) - lgamma(f_t + 1));
393 
394  e_t_1 = d_t;
395  }
396 
397  return log_prior;
398 }*/
399 
400 
401 }} // namespace kjb::mcmcda
402 
403 #endif /*MCMCDA_PRIOR_H_INCLUDED */
404 
Definition of various standard probability distributions.
boost::math::poisson Poisson_distribution
Definition: prob_distribution.h:69
r
Definition: APPgetLargeConnectedEdges.m:127
Computes the prior of an MCMCDA association.
Definition: mcmcda_prior.h:38
double kappa() const
Get average number of before start of video.
Definition: mcmcda_prior.h:79
A class that represents a MCMCDA association.
Definition: mcmcda_association.h:66
void set_kappa(double kappa)
Set average number of tracks per time.
Definition: mcmcda_prior.h:98
double lambda_N() const
Get false detection rate.
Definition: mcmcda_prior.h:85
Vector sample(const MV_gaussian_distribution &dist)
Sample from a multivariate normal distribution.
Definition: prob_sample.cpp:42
#define IFT(a, ex, msg)
Definition: l_exception.h:101
Prior(double kappa, double theta, double lambda_N, double lambda_A)
Constructs a MCMCDA prior functor.
Definition: mcmcda_prior.h:58
void set_lambda_A(double lambda_A)
Set detection rate.
Definition: mcmcda_prior.h:107
const Data< Element > & get_data() const
Returns data set const-ref.
Definition: mcmcda_association.h:126
std::pair< std::vector< Track >, std::vector< size_t > > sample(const Prior< Track > &prior, size_t T, const Track &def)
Sample an association from the prior.
Definition: mcmcda_prior.h:193
double operator()(const Association< Track > &w) const
Apply this functor to an association.
Definition: mcmcda_prior.h:118
double theta() const
Get average track length.
Definition: mcmcda_prior.h:82
Int_matrix::Value_type min(const Int_matrix &mat)
Return the minimum value in this matrix.
Definition: l_int_matrix.h:1385
void get_association_totals(const Association< Track > &w, size_t &m, size_t &e, size_t &d, size_t &a, size_t &n, size_t &l)
Get total tracks, entrances, exits, true detections, noisy detections, and track lengths of a given a...
Definition: mcmcda_association.h:605
Sampling functionality for the different distributions defined in "prob_distributions.h".
void set_lambda_N(double lambda_N)
Set false detection rate.
Definition: mcmcda_prior.h:104
Object thrown when an argument to a function is not acceptable.
Definition: l_exception.h:377
double lambda_O() const
Get occluded detection rate.
Definition: mcmcda_prior.h:91
get the indices of edges in each direction for i
Definition: APPgetLargeConnectedEdges.m:48
for m
Definition: APPgetLargeConnectedEdges.m:64
void set_lambda_O(double lambda_O)
Set occluded detection rate.
Definition: mcmcda_prior.h:110
double lambda_A() const
Get detection rate.
Definition: mcmcda_prior.h:88
void set_theta(double theta)
Set average track length.
Definition: mcmcda_prior.h:101