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