Improved factor.h/cpp code and finished corresponding unit tests
[libdai.git] / include / dai / prob.h
1 /* This file is part of libDAI - http://www.libdai.org/
2 *
3 * libDAI is licensed under the terms of the GNU General Public License version
4 * 2, or (at your option) any later version. libDAI is distributed without any
5 * warranty. See the file COPYING for more details.
6 *
7 * Copyright (C) 2006-2009 Joris Mooij [joris dot mooij at libdai dot org]
8 * Copyright (C) 2006-2007 Radboud University Nijmegen, The Netherlands
9 */
10
11
12 /// \file
13 /// \brief Defines TProb<> and Prob classes which represent (probability) vectors (e.g., probability distributions of discrete random variables)
14
15
16 #ifndef __defined_libdai_prob_h
17 #define __defined_libdai_prob_h
18
19
20 #include <cmath>
21 #include <vector>
22 #include <ostream>
23 #include <algorithm>
24 #include <numeric>
25 #include <functional>
26 #include <dai/util.h>
27 #include <dai/exceptions.h>
28
29
30 namespace dai {
31
32
33 /// Function object that returns the value itself
34 template<typename T> struct fo_id : public std::unary_function<T, T> {
35 /// Returns \a x
36 T operator()( const T &x ) const {
37 return x;
38 }
39 };
40
41
42 /// Function object that takes the absolute value
43 template<typename T> struct fo_abs : public std::unary_function<T, T> {
44 /// Returns abs(\a x)
45 T operator()( const T &x ) const {
46 if( x < (T)0 )
47 return -x;
48 else
49 return x;
50 }
51 };
52
53
54 /// Function object that takes the exponent
55 template<typename T> struct fo_exp : public std::unary_function<T, T> {
56 /// Returns exp(\a x)
57 T operator()( const T &x ) const {
58 return exp( x );
59 }
60 };
61
62
63 /// Function object that takes the logarithm
64 template<typename T> struct fo_log : public std::unary_function<T, T> {
65 /// Returns log(\a x)
66 T operator()( const T &x ) const {
67 return log( x );
68 }
69 };
70
71
72 /// Function object that takes the logarithm, except that log(0) is defined to be 0
73 template<typename T> struct fo_log0 : public std::unary_function<T, T> {
74 /// Returns (\a x == 0 ? 0 : log(\a x))
75 T operator()( const T &x ) const {
76 if( x )
77 return log( x );
78 else
79 return 0;
80 }
81 };
82
83
84 /// Function object that takes the inverse
85 template<typename T> struct fo_inv : public std::unary_function<T, T> {
86 /// Returns 1 / \a x
87 T operator()( const T &x ) const {
88 return 1 / x;
89 }
90 };
91
92
93 /// Function object that takes the inverse, except that 1/0 is defined to be 0
94 template<typename T> struct fo_inv0 : public std::unary_function<T, T> {
95 /// Returns (\a x == 0 ? 0 : (1 / \a x))
96 T operator()( const T &x ) const {
97 if( x )
98 return 1 / x;
99 else
100 return 0;
101 }
102 };
103
104
105 /// Function object that returns p*log0(p)
106 template<typename T> struct fo_plog0p : public std::unary_function<T, T> {
107 /// Returns \a p * log0(\a p)
108 T operator()( const T &p ) const {
109 return p * dai::log0(p);
110 }
111 };
112
113
114 /// Function object similar to std::divides(), but different in that dividing by zero results in zero
115 template<typename T> struct fo_divides0 : public std::binary_function<T, T, T> {
116 /// Returns (\a y == 0 ? 0 : (\a x / \a y))
117 T operator()( const T &x, const T &y ) const {
118 if( y == (T)0 )
119 return (T)0;
120 else
121 return x / y;
122 }
123 };
124
125
126 /// Function object useful for calculating the KL distance
127 template<typename T> struct fo_KL : public std::binary_function<T, T, T> {
128 /// Returns (\a p == 0 ? 0 : (\a p * (log(\a p) - log(\a q))))
129 T operator()( const T &p, const T &q ) const {
130 if( p == (T)0 )
131 return (T)0;
132 else
133 return p * (log(p) - log(q));
134 }
135 };
136
137
138 /// Function object useful for calculating the Hellinger distance
139 template<typename T> struct fo_Hellinger : public std::binary_function<T, T, T> {
140 /// Returns (sqrt(\a p) - sqrt(\a q))^2
141 T operator()( const T &p, const T &q ) const {
142 T x = sqrt(p) - sqrt(q);
143 return x * x;
144 }
145 };
146
147
148 /// Function object that returns x to the power y
149 template<typename T> struct fo_pow : public std::binary_function<T, T, T> {
150 /// Returns (\a x ^ \a y)
151 T operator()( const T &x, const T &y ) const {
152 if( y != 1 )
153 return std::pow( x, y );
154 else
155 return x;
156 }
157 };
158
159
160 /// Function object that returns the maximum of two values
161 template<typename T> struct fo_max : public std::binary_function<T, T, T> {
162 /// Returns (\a x > y ? x : y)
163 T operator()( const T &x, const T &y ) const {
164 return (x > y) ? x : y;
165 }
166 };
167
168
169 /// Function object that returns the minimum of two values
170 template<typename T> struct fo_min : public std::binary_function<T, T, T> {
171 /// Returns (\a x > y ? y : x)
172 T operator()( const T &x, const T &y ) const {
173 return (x > y) ? y : x;
174 }
175 };
176
177
178 /// Function object that returns the absolute difference of x and y
179 template<typename T> struct fo_absdiff : public std::binary_function<T, T, T> {
180 /// Returns abs( \a x - \a y )
181 T operator()( const T &x, const T &y ) const {
182 return dai::abs( x - y );
183 }
184 };
185
186
187 /// Represents a vector with entries of type \a T.
188 /** It is simply a <tt>std::vector</tt><<em>T</em>> with an interface designed for dealing with probability mass functions.
189 *
190 * It is mainly used for representing measures on a finite outcome space, for example, the probability
191 * distribution of a discrete random variable. However, entries are not necessarily non-negative; it is also used to
192 * represent logarithms of probability mass functions.
193 *
194 * \tparam T Should be a scalar that is castable from and to dai::Real and should support elementary arithmetic operations.
195 */
196 template <typename T>
197 class TProb {
198 public:
199 /// Type of data structure used for storing the values
200 typedef std::vector<T> container_type;
201
202 private:
203 /// The data structure that stores the values
204 container_type _p;
205
206 public:
207 /// Enumerates different ways of normalizing a probability measure.
208 /**
209 * - NORMPROB means that the sum of all entries should be 1;
210 * - NORMLINF means that the maximum absolute value of all entries should be 1.
211 */
212 typedef enum { NORMPROB, NORMLINF } NormType;
213 /// Enumerates different distance measures between probability measures.
214 /**
215 * - DISTL1 is the \f$\ell_1\f$ distance (sum of absolute values of pointwise difference);
216 * - DISTLINF is the \f$\ell_\infty\f$ distance (maximum absolute value of pointwise difference);
217 * - DISTTV is the total variation distance (half of the \f$\ell_1\f$ distance);
218 * - DISTKL is the Kullback-Leibler distance (\f$\sum_i p_i (\log p_i - \log q_i)\f$).
219 * - DISTHEL is the Hellinger distance (\f$\frac{1}{2}\sum_i (\sqrt{p_i}-\sqrt{q_i})^2\f$).
220 */
221 typedef enum { DISTL1, DISTLINF, DISTTV, DISTKL, DISTHEL } DistType;
222
223 /// \name Constructors and destructors
224 //@{
225 /// Default constructor (constructs empty vector)
226 TProb() : _p() {}
227
228 /// Construct uniform probability distribution over \a n outcomes (i.e., a vector of length \a n with each entry set to \f$1/n\f$)
229 explicit TProb( size_t n ) : _p( n, (T)1 / n ) {}
230
231 /// Construct vector of length \a n with each entry set to \a p
232 explicit TProb( size_t n, T p ) : _p( n, p ) {}
233
234 /// Construct vector from a range
235 /** \tparam TIterator Iterates over instances that can be cast to \a T
236 * \param begin Points to first instance to be added.
237 * \param end Points just beyond last instance to be added.
238 * \param sizeHint For efficiency, the number of entries can be speficied by \a sizeHint.
239 */
240 template <typename TIterator>
241 TProb( TIterator begin, TIterator end, size_t sizeHint=0 ) : _p() {
242 _p.reserve( sizeHint );
243 _p.insert( _p.begin(), begin, end );
244 }
245
246 /// Construct vector from another vector
247 /** \tparam S type of elements in \a v (should be castable to type \a T)
248 * \param v vector used for initialization.
249 */
250 template <typename S>
251 TProb( const std::vector<S> &v ) : _p() {
252 _p.reserve( v.size() );
253 _p.insert( _p.begin(), v.begin(), v.end() );
254 }
255 //@}
256
257 /// Constant iterator over the elements
258 typedef typename container_type::const_iterator const_iterator;
259 /// Iterator over the elements
260 typedef typename container_type::iterator iterator;
261 /// Constant reverse iterator over the elements
262 typedef typename container_type::const_reverse_iterator const_reverse_iterator;
263 /// Reverse iterator over the elements
264 typedef typename container_type::reverse_iterator reverse_iterator;
265
266 /// \name Iterator interface
267 //@{
268 /// Returns iterator that points to the first element
269 iterator begin() { return _p.begin(); }
270 /// Returns constant iterator that points to the first element
271 const_iterator begin() const { return _p.begin(); }
272
273 /// Returns iterator that points beyond the last element
274 iterator end() { return _p.end(); }
275 /// Returns constant iterator that points beyond the last element
276 const_iterator end() const { return _p.end(); }
277
278 /// Returns reverse iterator that points to the last element
279 reverse_iterator rbegin() { return _p.rbegin(); }
280 /// Returns constant reverse iterator that points to the last element
281 const_reverse_iterator rbegin() const { return _p.rbegin(); }
282
283 /// Returns reverse iterator that points beyond the first element
284 reverse_iterator rend() { return _p.rend(); }
285 /// Returns constant reverse iterator that points beyond the first element
286 const_reverse_iterator rend() const { return _p.rend(); }
287 //@}
288
289 /// \name Miscellaneous operations
290 //@{
291 void resize( size_t sz, T c = T() ) {
292 _p.resize( sz, c );
293 }
294 //@}
295
296 /// \name Queries
297 //@{
298 /// Gets \a i 'th entry
299 T get( size_t i ) const {
300 #ifdef DAI_DEBUG
301 return _p.at(i);
302 #else
303 return _p[i];
304 #endif
305 }
306
307 /// Sets \a i 'th entry to \a val
308 void set( size_t i, T val ) {
309 DAI_DEBASSERT( i < _p.size() );
310 _p[i] = val;
311 }
312
313 /// Returns a const reference to the wrapped vector
314 const container_type& p() const { return _p; }
315
316 /// Returns a reference to the wrapped vector
317 container_type& p() { return _p; }
318
319 /// Returns a copy of the \a i 'th entry
320 T operator[]( size_t i ) const { return get(i); }
321
322 /// Returns reference to the \a i 'th entry
323 /** \deprecated Please use dai::TProb::set() instead
324 */
325 T& operator[]( size_t i ) { return _p[i]; }
326
327 /// Returns length of the vector (i.e., the number of entries)
328 size_t size() const { return _p.size(); }
329
330 /// Accumulate over all values, similar to std::accumulate
331 /** The following calculation is done:
332 * \code
333 * T t = op2(init);
334 * for( const_iterator it = begin(); it != end(); it++ )
335 * t = op1( t, op2(*it) );
336 * return t;
337 * \endcode
338 */
339 template<typename binOp, typename unOp> T accumulate( T init, binOp op1, unOp op2 ) const {
340 T t = op2(init);
341 for( const_iterator it = begin(); it != end(); it++ )
342 t = op1( t, op2(*it) );
343 return t;
344 }
345
346 /// Returns the Shannon entropy of \c *this, \f$-\sum_i p_i \log p_i\f$
347 T entropy() const { return -accumulate( (T)0, std::plus<T>(), fo_plog0p<T>() ); }
348
349 /// Returns maximum value of all entries
350 T max() const { return accumulate( (T)(-INFINITY), fo_max<T>(), fo_id<T>() ); }
351
352 /// Returns minimum value of all entries
353 T min() const { return accumulate( (T)INFINITY, fo_min<T>(), fo_id<T>() ); }
354
355 /// Returns sum of all entries
356 T sum() const { return accumulate( (T)0, std::plus<T>(), fo_id<T>() ); }
357
358 /// Return sum of absolute value of all entries
359 T sumAbs() const { return accumulate( (T)0, std::plus<T>(), fo_abs<T>() ); }
360
361 /// Returns maximum absolute value of all entries
362 T maxAbs() const { return accumulate( (T)0, fo_max<T>(), fo_abs<T>() ); }
363
364 /// Returns \c true if one or more entries are NaN
365 bool hasNaNs() const {
366 bool foundnan = false;
367 for( const_iterator x = _p.begin(); x != _p.end(); x++ )
368 if( isnan( *x ) ) {
369 foundnan = true;
370 break;
371 }
372 return foundnan;
373 }
374
375 /// Returns \c true if one or more entries are negative
376 bool hasNegatives() const {
377 return (std::find_if( _p.begin(), _p.end(), std::bind2nd( std::less<T>(), (T)0 ) ) != _p.end());
378 }
379
380 /// Returns a pair consisting of the index of the maximum value and the maximum value itself
381 std::pair<size_t,T> argmax() const {
382 T max = _p[0];
383 size_t arg = 0;
384 for( size_t i = 1; i < size(); i++ ) {
385 if( _p[i] > max ) {
386 max = _p[i];
387 arg = i;
388 }
389 }
390 return std::make_pair( arg, max );
391 }
392
393 /// Returns a random index, according to the (normalized) distribution described by *this
394 size_t draw() {
395 Real x = rnd_uniform() * sum();
396 T s = 0;
397 for( size_t i = 0; i < size(); i++ ) {
398 s += get(i);
399 if( s > x )
400 return i;
401 }
402 return( size() - 1 );
403 }
404
405 /// Lexicographical comparison
406 /** \pre <tt>this->size() == q.size()</tt>
407 */
408 bool operator<( const TProb<T>& q ) const {
409 DAI_DEBASSERT( size() == q.size() );
410 return lexicographical_compare( begin(), end(), q.begin(), q.end() );
411 }
412
413 /// Comparison
414 bool operator==( const TProb<T>& q ) const {
415 if( size() != q.size() )
416 return false;
417 return p() == q.p();
418 }
419 //@}
420
421 /// \name Unary transformations
422 //@{
423 /// Returns the result of applying operation \a op pointwise on \c *this
424 template<typename unaryOp> TProb<T> pwUnaryTr( unaryOp op ) const {
425 TProb<T> r;
426 r._p.reserve( size() );
427 std::transform( _p.begin(), _p.end(), back_inserter( r._p ), op );
428 return r;
429 }
430
431 /// Returns negative of \c *this
432 TProb<T> operator- () const { return pwUnaryTr( std::negate<T>() ); }
433
434 /// Returns pointwise absolute value
435 TProb<T> abs() const { return pwUnaryTr( fo_abs<T>() ); }
436
437 /// Returns pointwise exponent
438 TProb<T> exp() const { return pwUnaryTr( fo_exp<T>() ); }
439
440 /// Returns pointwise logarithm
441 /** If \a zero == \c true, uses <tt>log(0)==0</tt>; otherwise, <tt>log(0)==-Inf</tt>.
442 */
443 TProb<T> log(bool zero=false) const {
444 if( zero )
445 return pwUnaryTr( fo_log0<T>() );
446 else
447 return pwUnaryTr( fo_log<T>() );
448 }
449
450 /// Returns pointwise inverse
451 /** If \a zero == \c true, uses <tt>1/0==0</tt>; otherwise, <tt>1/0==Inf</tt>.
452 */
453 TProb<T> inverse(bool zero=true) const {
454 if( zero )
455 return pwUnaryTr( fo_inv0<T>() );
456 else
457 return pwUnaryTr( fo_inv<T>() );
458 }
459
460 /// Returns normalized copy of \c *this, using the specified norm
461 /** \throw NOT_NORMALIZABLE if the norm is zero
462 */
463 TProb<T> normalized( NormType norm = NORMPROB ) const {
464 T Z = 0;
465 if( norm == NORMPROB )
466 Z = sum();
467 else if( norm == NORMLINF )
468 Z = maxAbs();
469 if( Z == (T)0 ) {
470 DAI_THROW(NOT_NORMALIZABLE);
471 return *this;
472 } else
473 return pwUnaryTr( std::bind2nd( std::divides<T>(), Z ) );
474 }
475 //@}
476
477 /// \name Unary operations
478 //@{
479 /// Applies unary operation \a op pointwise
480 template<typename unaryOp> TProb<T>& pwUnaryOp( unaryOp op ) {
481 std::transform( _p.begin(), _p.end(), _p.begin(), op );
482 return *this;
483 }
484
485 /// Draws all entries i.i.d. from a uniform distribution on [0,1)
486 TProb<T>& randomize() {
487 std::generate( _p.begin(), _p.end(), rnd_uniform );
488 return *this;
489 }
490
491 /// Sets all entries to \f$1/n\f$ where \a n is the length of the vector
492 TProb<T>& setUniform () {
493 fill( (T)1 / size() );
494 return *this;
495 }
496
497 /// Applies absolute value pointwise
498 TProb<T>& takeAbs() { return pwUnaryOp( fo_abs<T>() ); }
499
500 /// Applies exponent pointwise
501 TProb<T>& takeExp() { return pwUnaryOp( fo_exp<T>() ); }
502
503 /// Applies logarithm pointwise
504 /** If \a zero == \c true, uses <tt>log(0)==0</tt>; otherwise, <tt>log(0)==-Inf</tt>.
505 */
506 TProb<T>& takeLog(bool zero=false) {
507 if( zero ) {
508 return pwUnaryOp( fo_log0<T>() );
509 } else
510 return pwUnaryOp( fo_log<T>() );
511 }
512
513 /// Normalizes vector using the specified norm
514 /** \throw NOT_NORMALIZABLE if the norm is zero
515 */
516 T normalize( NormType norm=NORMPROB ) {
517 T Z = 0;
518 if( norm == NORMPROB )
519 Z = sum();
520 else if( norm == NORMLINF )
521 Z = maxAbs();
522 if( Z == (T)0 )
523 DAI_THROW(NOT_NORMALIZABLE);
524 else
525 *this /= Z;
526 return Z;
527 }
528 //@}
529
530 /// \name Operations with scalars
531 //@{
532 /// Sets all entries to \a x
533 TProb<T> & fill( T x ) {
534 std::fill( _p.begin(), _p.end(), x );
535 return *this;
536 }
537
538 /// Adds scalar \a x to each entry
539 TProb<T>& operator+= (T x) {
540 if( x != 0 )
541 return pwUnaryOp( std::bind2nd( std::plus<T>(), x ) );
542 else
543 return *this;
544 }
545
546 /// Subtracts scalar \a x from each entry
547 TProb<T>& operator-= (T x) {
548 if( x != 0 )
549 return pwUnaryOp( std::bind2nd( std::minus<T>(), x ) );
550 else
551 return *this;
552 }
553
554 /// Multiplies each entry with scalar \a x
555 TProb<T>& operator*= (T x) {
556 if( x != 1 )
557 return pwUnaryOp( std::bind2nd( std::multiplies<T>(), x ) );
558 else
559 return *this;
560 }
561
562 /// Divides each entry by scalar \a x, where division by 0 yields 0
563 TProb<T>& operator/= (T x) {
564 if( x != 1 )
565 return pwUnaryOp( std::bind2nd( fo_divides0<T>(), x ) );
566 else
567 return *this;
568 }
569
570 /// Raises entries to the power \a x
571 TProb<T>& operator^= (T x) {
572 if( x != (T)1 )
573 return pwUnaryOp( std::bind2nd( fo_pow<T>(), x) );
574 else
575 return *this;
576 }
577 //@}
578
579 /// \name Transformations with scalars
580 //@{
581 /// Returns sum of \c *this and scalar \a x
582 TProb<T> operator+ (T x) const { return pwUnaryTr( std::bind2nd( std::plus<T>(), x ) ); }
583
584 /// Returns difference of \c *this and scalar \a x
585 TProb<T> operator- (T x) const { return pwUnaryTr( std::bind2nd( std::minus<T>(), x ) ); }
586
587 /// Returns product of \c *this with scalar \a x
588 TProb<T> operator* (T x) const { return pwUnaryTr( std::bind2nd( std::multiplies<T>(), x ) ); }
589
590 /// Returns quotient of \c *this and scalar \a x, where division by 0 yields 0
591 TProb<T> operator/ (T x) const { return pwUnaryTr( std::bind2nd( fo_divides0<T>(), x ) ); }
592
593 /// Returns \c *this raised to the power \a x
594 TProb<T> operator^ (T x) const { return pwUnaryTr( std::bind2nd( fo_pow<T>(), x ) ); }
595 //@}
596
597 /// \name Operations with other equally-sized vectors
598 //@{
599 /// Applies binary operation pointwise on two vectors
600 /** \tparam binaryOp Type of function object that accepts two arguments of type \a T and outputs a type \a T
601 * \param q Right operand
602 * \param op Operation of type \a binaryOp
603 */
604 template<typename binaryOp> TProb<T>& pwBinaryOp( const TProb<T> &q, binaryOp op ) {
605 DAI_DEBASSERT( size() == q.size() );
606 std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), op );
607 return *this;
608 }
609
610 /// Pointwise addition with \a q
611 /** \pre <tt>this->size() == q.size()</tt>
612 */
613 TProb<T>& operator+= (const TProb<T> & q) { return pwBinaryOp( q, std::plus<T>() ); }
614
615 /// Pointwise subtraction of \a q
616 /** \pre <tt>this->size() == q.size()</tt>
617 */
618 TProb<T>& operator-= (const TProb<T> & q) { return pwBinaryOp( q, std::minus<T>() ); }
619
620 /// Pointwise multiplication with \a q
621 /** \pre <tt>this->size() == q.size()</tt>
622 */
623 TProb<T>& operator*= (const TProb<T> & q) { return pwBinaryOp( q, std::multiplies<T>() ); }
624
625 /// Pointwise division by \a q, where division by 0 yields 0
626 /** \pre <tt>this->size() == q.size()</tt>
627 * \see divide(const TProb<T> &)
628 */
629 TProb<T>& operator/= (const TProb<T> & q) { return pwBinaryOp( q, fo_divides0<T>() ); }
630
631 /// Pointwise division by \a q, where division by 0 yields +Inf
632 /** \pre <tt>this->size() == q.size()</tt>
633 * \see operator/=(const TProb<T> &)
634 */
635 TProb<T>& divide (const TProb<T> & q) { return pwBinaryOp( q, std::divides<T>() ); }
636
637 /// Pointwise power
638 /** \pre <tt>this->size() == q.size()</tt>
639 */
640 TProb<T>& operator^= (const TProb<T> & q) { return pwBinaryOp( q, fo_pow<T>() ); }
641 //@}
642
643 /// \name Transformations with other equally-sized vectors
644 //@{
645 /// Returns the result of applying binary operation \a op pointwise on \c *this and \a q
646 /** \tparam binaryOp Type of function object that accepts two arguments of type \a T and outputs a type \a T
647 * \param q Right operand
648 * \param op Operation of type \a binaryOp
649 */
650 template<typename binaryOp> TProb<T> pwBinaryTr( const TProb<T> &q, binaryOp op ) const {
651 DAI_DEBASSERT( size() == q.size() );
652 TProb<T> r;
653 r._p.reserve( size() );
654 std::transform( _p.begin(), _p.end(), q._p.begin(), back_inserter( r._p ), op );
655 return r;
656 }
657
658 /// Returns sum of \c *this and \a q
659 /** \pre <tt>this->size() == q.size()</tt>
660 */
661 TProb<T> operator+ ( const TProb<T>& q ) const { return pwBinaryTr( q, std::plus<T>() ); }
662
663 /// Return \c *this minus \a q
664 /** \pre <tt>this->size() == q.size()</tt>
665 */
666 TProb<T> operator- ( const TProb<T>& q ) const { return pwBinaryTr( q, std::minus<T>() ); }
667
668 /// Return product of \c *this with \a q
669 /** \pre <tt>this->size() == q.size()</tt>
670 */
671 TProb<T> operator* ( const TProb<T> &q ) const { return pwBinaryTr( q, std::multiplies<T>() ); }
672
673 /// Returns quotient of \c *this with \a q, where division by 0 yields 0
674 /** \pre <tt>this->size() == q.size()</tt>
675 * \see divided_by(const TProb<T> &)
676 */
677 TProb<T> operator/ ( const TProb<T> &q ) const { return pwBinaryTr( q, fo_divides0<T>() ); }
678
679 /// Pointwise division by \a q, where division by 0 yields +Inf
680 /** \pre <tt>this->size() == q.size()</tt>
681 * \see operator/(const TProb<T> &)
682 */
683 TProb<T> divided_by( const TProb<T> &q ) const { return pwBinaryTr( q, std::divides<T>() ); }
684
685 /// Returns \c *this to the power \a q
686 /** \pre <tt>this->size() == q.size()</tt>
687 */
688 TProb<T> operator^ ( const TProb<T> &q ) const { return pwBinaryTr( q, fo_pow<T>() ); }
689 //@}
690
691 /// Performs a generalized inner product, similar to std::inner_product
692 /** \pre <tt>this->size() == q.size()</tt>
693 */
694 template<typename binOp1, typename binOp2> T innerProduct( const TProb<T> &q, T init, binOp1 binaryOp1, binOp2 binaryOp2 ) const {
695 DAI_DEBASSERT( size() == q.size() );
696 return std::inner_product( begin(), end(), q.begin(), init, binaryOp1, binaryOp2 );
697 }
698 };
699
700
701 /// Returns distance between \a p and \a q, measured using distance measure \a dt
702 /** \relates TProb
703 * \pre <tt>this->size() == q.size()</tt>
704 */
705 template<typename T> T dist( const TProb<T> &p, const TProb<T> &q, typename TProb<T>::DistType dt ) {
706 switch( dt ) {
707 case TProb<T>::DISTL1:
708 return p.innerProduct( q, (T)0, std::plus<T>(), fo_absdiff<T>() );
709 case TProb<T>::DISTLINF:
710 return p.innerProduct( q, (T)0, fo_max<T>(), fo_absdiff<T>() );
711 case TProb<T>::DISTTV:
712 return p.innerProduct( q, (T)0, std::plus<T>(), fo_absdiff<T>() ) / 2;
713 case TProb<T>::DISTKL:
714 return p.innerProduct( q, (T)0, std::plus<T>(), fo_KL<T>() );
715 case TProb<T>::DISTHEL:
716 return p.innerProduct( q, (T)0, std::plus<T>(), fo_Hellinger<T>() ) / 2;
717 default:
718 DAI_THROW(UNKNOWN_ENUM_VALUE);
719 return INFINITY;
720 }
721 }
722
723
724 /// Writes a TProb<T> to an output stream
725 /** \relates TProb
726 */
727 template<typename T> std::ostream& operator<< (std::ostream& os, const TProb<T>& p) {
728 os << "(";
729 for( typename TProb<T>::const_iterator it = p.begin(); it != p.end(); it++ )
730 os << (it != p.begin() ? ", " : "") << *it;
731 os << ")";
732 return os;
733 }
734
735
736 /// Returns the pointwise minimum of \a a and \a b
737 /** \relates TProb
738 * \pre <tt>this->size() == q.size()</tt>
739 */
740 template<typename T> TProb<T> min( const TProb<T> &a, const TProb<T> &b ) {
741 return a.pwBinaryTr( b, fo_min<T>() );
742 }
743
744
745 /// Returns the pointwise maximum of \a a and \a b
746 /** \relates TProb
747 * \pre <tt>this->size() == q.size()</tt>
748 */
749 template<typename T> TProb<T> max( const TProb<T> &a, const TProb<T> &b ) {
750 return a.pwBinaryTr( b, fo_max<T>() );
751 }
752
753
754 /// Represents a vector with entries of type dai::Real.
755 typedef TProb<Real> Prob;
756
757
758 } // end of namespace dai
759
760
761 #endif