Replaced doubles by Reals, fixed two bugs
[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
14 /// \todo Rename to Vector<>
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 /// Represents a vector with entries of type \a T.
35 /** It is simply a <tt>std::vector</tt><<em>T</em>> with an interface designed for dealing with probability mass functions.
36 *
37 * It is mainly used for representing measures on a finite outcome space, for example, the probability
38 * distribution of a discrete random variable. However, entries are not necessarily non-negative; it is also used to
39 * represent logarithms of probability mass functions.
40 *
41 * \tparam T Should be a scalar that is castable from and to dai::Real and should support elementary arithmetic operations.
42 */
43 template <typename T> class TProb {
44 private:
45 /// The vector
46 std::vector<T> _p;
47
48 public:
49 /// Enumerates different ways of normalizing a probability measure.
50 /**
51 * - NORMPROB means that the sum of all entries should be 1;
52 * - NORMLINF means that the maximum absolute value of all entries should be 1.
53 */
54 typedef enum { NORMPROB, NORMLINF } NormType;
55 /// Enumerates different distance measures between probability measures.
56 /**
57 * - DISTL1 is the \f$\ell_1\f$ distance (sum of absolute values of pointwise difference);
58 * - DISTLINF is the \f$\ell_\infty\f$ distance (maximum absolute value of pointwise difference);
59 * - DISTTV is the total variation distance (half of the \f$\ell_1\f$ distance);
60 * - DISTKL is the Kullback-Leibler distance (\f$\sum_i p_i (\log p_i - \log q_i)\f$).
61 */
62 typedef enum { DISTL1, DISTLINF, DISTTV, DISTKL } DistType;
63
64 /// \name Constructors and destructors
65 //@{
66 /// Default constructor (constructs empty vector)
67 TProb() : _p() {}
68
69 /// 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$)
70 explicit TProb( size_t n ) : _p(std::vector<T>(n, (T)1 / n)) {}
71
72 /// Construct vector of length \a n with each entry set to \a p
73 explicit TProb( size_t n, T p ) : _p(n, p) {}
74
75 /// Construct vector from a range
76 /** \tparam TIterator Iterates over instances that can be cast to \a T
77 * \param begin Points to first instance to be added.
78 * \param end Points just beyond last instance to be added.
79 * \param sizeHint For efficiency, the number of entries can be speficied by \a sizeHint.
80 */
81 template <typename TIterator>
82 TProb( TIterator begin, TIterator end, size_t sizeHint=0 ) : _p() {
83 _p.reserve( sizeHint );
84 _p.insert( _p.begin(), begin, end );
85 }
86
87 /// Construct vector from another vector
88 /** \tparam S type of elements in \a v (should be castable to type \a T)
89 * \param v vector used for initialization
90 */
91 template <typename S>
92 TProb( const std::vector<S> &v ) : _p() {
93 _p.reserve( v.size() );
94 _p.insert( _p.begin(), v.begin(), v.end() );
95 }
96 //@}
97
98 /// Constant iterator over the elements
99 typedef typename std::vector<T>::const_iterator const_iterator;
100 /// Iterator over the elements
101 typedef typename std::vector<T>::iterator iterator;
102 /// Constant reverse iterator over the elements
103 typedef typename std::vector<T>::const_reverse_iterator const_reverse_iterator;
104 /// Reverse iterator over the elements
105 typedef typename std::vector<T>::reverse_iterator reverse_iterator;
106
107 /// @name Iterator interface
108 //@{
109 /// Returns iterator that points to the first element
110 iterator begin() { return _p.begin(); }
111 /// Returns constant iterator that points to the first element
112 const_iterator begin() const { return _p.begin(); }
113
114 /// Returns iterator that points beyond the last element
115 iterator end() { return _p.end(); }
116 /// Returns constant iterator that points beyond the last element
117 const_iterator end() const { return _p.end(); }
118
119 /// Returns reverse iterator that points to the last element
120 reverse_iterator rbegin() { return _p.rbegin(); }
121 /// Returns constant reverse iterator that points to the last element
122 const_reverse_iterator rbegin() const { return _p.rbegin(); }
123
124 /// Returns reverse iterator that points beyond the first element
125 reverse_iterator rend() { return _p.rend(); }
126 /// Returns constant reverse iterator that points beyond the first element
127 const_reverse_iterator rend() const { return _p.rend(); }
128 //@}
129
130 /// \name Queries
131 //@{
132 /// Returns a const reference to the wrapped vector
133 const std::vector<T> & p() const { return _p; }
134
135 /// Returns a reference to the wrapped vector
136 std::vector<T> & p() { return _p; }
137
138 /// Returns a copy of the \a i 'th entry
139 T operator[]( size_t i ) const {
140 #ifdef DAI_DEBUG
141 return _p.at(i);
142 #else
143 return _p[i];
144 #endif
145 }
146
147 /// Returns reference to the \a i 'th entry
148 T& operator[]( size_t i ) { return _p[i]; }
149
150 /// Returns length of the vector (i.e., the number of entries)
151 size_t size() const { return _p.size(); }
152
153 /// Returns the Shannon entropy of \c *this, \f$-\sum_i p_i \log p_i\f$
154 T entropy() const {
155 T S = 0;
156 for( size_t i = 0; i < size(); i++ )
157 S -= (_p[i] == 0 ? 0 : _p[i] * dai::log(_p[i]));
158 return S;
159 }
160
161 /// Returns maximum value of all entries
162 T max() const {
163 T Z = *std::max_element( _p.begin(), _p.end() );
164 return Z;
165 }
166
167 /// Returns minimum value of all entries
168 T min() const {
169 T Z = *std::min_element( _p.begin(), _p.end() );
170 return Z;
171 }
172
173 /// Returns a pair consisting of the index of the maximum value and the maximum value itself
174 std::pair<size_t,T> argmax() const {
175 T max = _p[0];
176 size_t arg = 0;
177 for( size_t i = 1; i < size(); i++ ) {
178 if( _p[i] > max ) {
179 max = _p[i];
180 arg = i;
181 }
182 }
183 return std::make_pair(arg,max);
184 }
185
186 /// Returns sum of all entries
187 T sum() const {
188 T Z = std::accumulate( _p.begin(), _p.end(), (T)0 );
189 return Z;
190 }
191
192 /// Return sum of absolute value of all entries
193 T sumAbs() const {
194 T s = 0;
195 for( size_t i = 0; i < size(); i++ )
196 s += dai::abs(_p[i]);
197 return s;
198 }
199
200 /// Returns maximum absolute value of all entries
201 T maxAbs() const {
202 T Z = 0;
203 for( size_t i = 0; i < size(); i++ ) {
204 T mag = dai::abs(_p[i]);
205 if( mag > Z )
206 Z = mag;
207 }
208 return Z;
209 }
210
211 /// Returns \c true if one or more entries are NaN
212 bool hasNaNs() const {
213 bool foundnan = false;
214 for( typename std::vector<T>::const_iterator x = _p.begin(); x != _p.end(); x++ )
215 if( isnan( *x ) ) {
216 foundnan = true;
217 break;
218 }
219 return foundnan;
220 }
221
222 /// Returns \c true if one or more entries are negative
223 bool hasNegatives() const {
224 return (std::find_if( _p.begin(), _p.end(), std::bind2nd( std::less<T>(), (T)0 ) ) != _p.end());
225 }
226
227 /// Returns a random index, according to the (normalized) distribution described by *this
228 size_t draw() {
229 Real x = rnd_uniform() * sum();
230 T s = 0;
231 for( size_t i = 0; i < size(); i++ ) {
232 s += _p[i];
233 if( s > x )
234 return i;
235 }
236 return( size() - 1 );
237 }
238
239 /// Lexicographical comparison
240 /** \pre <tt>this->size() == q.size()</tt>
241 */
242 bool operator<= (const TProb<T> & q) const {
243 DAI_DEBASSERT( size() == q.size() );
244 for( size_t i = 0; i < size(); i++ )
245 if( !(_p[i] <= q[i]) )
246 return false;
247 return true;
248 }
249 //@}
250
251 /// \name Unary transformations
252 //@{
253 // OBSOLETE
254 /// Returns pointwise signum
255 /** \note Obsolete, to be removed soon
256 */
257 TProb<T> sgn() const {
258 TProb<T> x;
259 x._p.reserve( size() );
260 for( size_t i = 0; i < size(); i++ ) {
261 T s = 0;
262 if( _p[i] > 0 )
263 s = 1;
264 else if( _p[i] < 0 )
265 s = -1;
266 x._p.push_back( s );
267 }
268 return x;
269 }
270
271 /// Returns pointwise absolute value
272 TProb<T> abs() const {
273 TProb<T> x;
274 x._p.reserve( size() );
275 for( size_t i = 0; i < size(); i++ )
276 x._p.push_back( abs(_p[i]) );
277 return x;
278 }
279
280 /// Returns pointwise exponent
281 TProb<T> exp() const {
282 TProb<T> e(*this);
283 e.takeExp();
284 return e;
285 }
286
287 /// Returns pointwise logarithm
288 /** If \a zero == \c true, uses <tt>log(0)==0</tt>; otherwise, <tt>log(0)==-Inf</tt>.
289 */
290 TProb<T> log(bool zero=false) const {
291 TProb<T> l(*this);
292 l.takeLog(zero);
293 return l;
294 }
295
296 /// Returns pointwise inverse
297 /** If \a zero == \c true, uses <tt>1/0==0</tt>; otherwise, <tt>1/0==Inf</tt>.
298 */
299 TProb<T> inverse(bool zero=true) const {
300 TProb<T> inv;
301 inv._p.reserve( size() );
302 if( zero )
303 for( size_t i = 0; i < size(); i++ )
304 inv._p.push_back( _p[i] == (T)0 ? (T)0 : (T)1 / _p[i] );
305 else
306 for( size_t i = 0; i < size(); i++ )
307 inv._p.push_back( (T)1 / _p[i] );
308 return inv;
309 }
310
311 /// Returns normalized copy of \c *this, using the specified norm
312 TProb<T> normalized( NormType norm = NORMPROB ) const {
313 TProb<T> result(*this);
314 result.normalize( norm );
315 return result;
316 }
317 //@}
318
319 /// \name Unary operations
320 //@{
321 /// Draws all entries i.i.d. from a uniform distribution on [0,1)
322 TProb<T>& randomize() {
323 std::generate( _p.begin(), _p.end(), rnd_uniform );
324 return *this;
325 }
326
327 /// Sets all entries to \f$1/n\f$ where \a n is the length of the vector
328 TProb<T>& setUniform () {
329 fill( (T)1 / size() );
330 return *this;
331 }
332
333 /// Applies exponent pointwise
334 const TProb<T>& takeExp() {
335 for( size_t i = 0; i < size(); i++ )
336 _p[i] = dai::exp(_p[i]);
337 return *this;
338 }
339
340 /// Applies logarithm pointwise
341 /** If \a zero == \c true, uses <tt>log(0)==0</tt>; otherwise, <tt>log(0)==-Inf</tt>.
342 */
343 const TProb<T>& takeLog(bool zero=false) {
344 if( zero ) {
345 for( size_t i = 0; i < size(); i++ )
346 _p[i] = ( (_p[i] == 0.0) ? 0.0 : dai::log(_p[i]) );
347 } else
348 for( size_t i = 0; i < size(); i++ )
349 _p[i] = dai::log(_p[i]);
350 return *this;
351 }
352
353 /// Normalizes vector using the specified norm
354 T normalize( NormType norm=NORMPROB ) {
355 T Z = 0.0;
356 if( norm == NORMPROB )
357 Z = sum();
358 else if( norm == NORMLINF )
359 Z = maxAbs();
360 if( Z == 0.0 )
361 DAI_THROW(NOT_NORMALIZABLE);
362 else
363 *this /= Z;
364 return Z;
365 }
366 //@}
367
368 /// \name Operations with scalars
369 //@{
370 /// Sets all entries to \a x
371 TProb<T> & fill(T x) {
372 std::fill( _p.begin(), _p.end(), x );
373 return *this;
374 }
375
376 // OBSOLETE
377 /// Sets entries that are smaller (in absolute value) than \a epsilon to 0
378 /** \note Obsolete, to be removed soon
379 */
380 TProb<T>& makeZero( T epsilon ) {
381 for( size_t i = 0; i < size(); i++ )
382 if( (_p[i] < epsilon) && (_p[i] > -epsilon) )
383 _p[i] = 0;
384 return *this;
385 }
386
387 // OBSOLETE
388 /// Sets entries that are smaller than \a epsilon to \a epsilon
389 /** \note Obsolete, to be removed soon
390 */
391 TProb<T>& makePositive( T epsilon ) {
392 for( size_t i = 0; i < size(); i++ )
393 if( (0 < _p[i]) && (_p[i] < epsilon) )
394 _p[i] = epsilon;
395 return *this;
396 }
397
398 /// Adds scalar \a x to each entry
399 TProb<T>& operator+= (T x) {
400 std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::plus<T>(), x ) );
401 return *this;
402 }
403
404 /// Subtracts scalar \a x from each entry
405 TProb<T>& operator-= (T x) {
406 std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::minus<T>(), x ) );
407 return *this;
408 }
409
410 /// Multiplies each entry with scalar \a x
411 TProb<T>& operator*= (T x) {
412 std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::multiplies<T>(), x) );
413 return *this;
414 }
415
416 /// Divides each entry by scalar \a x
417 TProb<T>& operator/= (T x) {
418 DAI_DEBASSERT( x != 0 );
419 std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::divides<T>(), x ) );
420 return *this;
421 }
422
423 /// Raises entries to the power \a x
424 TProb<T>& operator^= (T x) {
425 if( x != (T)1 )
426 std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::ptr_fun<T, T, T>(std::pow), x) );
427 return *this;
428 }
429 //@}
430
431 /// \name Transformations with scalars
432 //@{
433 /// Returns sum of \c *this and scalar \a x
434 TProb<T> operator+ (T x) const {
435 TProb<T> sum( *this );
436 sum += x;
437 return sum;
438 }
439
440 /// Returns difference of \c *this and scalar \a x
441 TProb<T> operator- (T x) const {
442 TProb<T> diff( *this );
443 diff -= x;
444 return diff;
445 }
446
447 /// Returns product of \c *this with scalar \a x
448 TProb<T> operator* (T x) const {
449 TProb<T> prod( *this );
450 prod *= x;
451 return prod;
452 }
453
454 /// Returns quotient of \c *this and scalar \a x
455 TProb<T> operator/ (T x) const {
456 TProb<T> quot( *this );
457 quot /= x;
458 return quot;
459 }
460
461 /// Returns \c *this raised to the power \a x
462 TProb<T> operator^ (T x) const {
463 TProb<T> power(*this);
464 power ^= x;
465 return power;
466 }
467 //@}
468
469 /// \name Operations with other equally-sized vectors
470 //@{
471 /// Pointwise addition with \a q
472 /** \pre <tt>this->size() == q.size()</tt>
473 */
474 TProb<T>& operator+= (const TProb<T> & q) {
475 DAI_DEBASSERT( size() == q.size() );
476 std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::plus<T>() );
477 return *this;
478 }
479
480 /// Pointwise subtraction of \a q
481 /** \pre <tt>this->size() == q.size()</tt>
482 */
483 TProb<T>& operator-= (const TProb<T> & q) {
484 DAI_DEBASSERT( size() == q.size() );
485 std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::minus<T>() );
486 return *this;
487 }
488
489 /// Pointwise multiplication with \a q
490 /** \pre <tt>this->size() == q.size()</tt>
491 */
492 TProb<T>& operator*= (const TProb<T> & q) {
493 DAI_DEBASSERT( size() == q.size() );
494 std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::multiplies<T>() );
495 return *this;
496 }
497
498 /// Pointwise division by \a q, where division by 0 yields 0
499 /** \pre <tt>this->size() == q.size()</tt>
500 * \see divide(const TProb<T> &)
501 */
502 TProb<T>& operator/= (const TProb<T> & q) {
503 DAI_DEBASSERT( size() == q.size() );
504 for( size_t i = 0; i < size(); i++ ) {
505 if( q[i] == 0.0 )
506 _p[i] = 0.0;
507 else
508 _p[i] /= q[i];
509 }
510 return *this;
511 }
512
513 /// Pointwise division by \a q, where division by 0 yields +Inf
514 /** \pre <tt>this->size() == q.size()</tt>
515 * \see operator/=(const TProb<T> &)
516 */
517 TProb<T>& divide (const TProb<T> & q) {
518 DAI_DEBASSERT( size() == q.size() );
519 std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::divides<T>() );
520 return *this;
521 }
522 //@}
523
524 /// \name Transformations with other equally-sized vectors
525 //@{
526 /// Returns sum of \c *this and \a q
527 /** \pre <tt>this->size() == q.size()</tt>
528 */
529 TProb<T> operator+ (const TProb<T> & q) const {
530 DAI_DEBASSERT( size() == q.size() );
531 TProb<T> sum( *this );
532 sum += q;
533 return sum;
534 }
535
536 /// Return \c *this minus \a q
537 /** \pre <tt>this->size() == q.size()</tt>
538 */
539 TProb<T> operator- (const TProb<T> & q) const {
540 DAI_DEBASSERT( size() == q.size() );
541 TProb<T> diff( *this );
542 diff -= q;
543 return diff;
544 }
545
546 /// Return product of \c *this with \a q
547 /** \pre <tt>this->size() == q.size()</tt>
548 */
549 TProb<T> operator* (const TProb<T> & q) const {
550 DAI_DEBASSERT( size() == q.size() );
551 TProb<T> prod( *this );
552 prod *= q;
553 return prod;
554 }
555
556 /// Returns quotient of \c *this with \a q, where division by 0 yields +Inf
557 /** \pre <tt>this->size() == q.size()</tt>
558 * \see divided_by(const TProb<T> &)
559 */
560 TProb<T> operator/ (const TProb<T> & q) const {
561 DAI_DEBASSERT( size() == q.size() );
562 TProb<T> quot( *this );
563 quot /= q;
564 return quot;
565 }
566
567 /// Pointwise division by \a q, where division by 0 yields 0
568 /** \pre <tt>this->size() == q.size()</tt>
569 * \see operator/(const TProb<T> &)
570 */
571 TProb<T> divided_by (const TProb<T> & q) const {
572 DAI_DEBASSERT( size() == q.size() );
573 TProb<T> quot( *this );
574 quot.divide(q);
575 return quot;
576 }
577 //@}
578 };
579
580
581 /// Returns distance between \a p and \a q, measured using distance measure \a dt
582 /** \relates TProb
583 * \pre <tt>this->size() == q.size()</tt>
584 */
585 template<typename T> T dist( const TProb<T> &p, const TProb<T> &q, typename TProb<T>::DistType dt ) {
586 DAI_DEBASSERT( p.size() == q.size() );
587 T result = 0;
588 switch( dt ) {
589 case TProb<T>::DISTL1:
590 for( size_t i = 0; i < p.size(); i++ )
591 result += abs(p[i] - q[i]);
592 break;
593
594 case TProb<T>::DISTLINF:
595 for( size_t i = 0; i < p.size(); i++ ) {
596 T z = abs(p[i] - q[i]);
597 if( z > result )
598 result = z;
599 }
600 break;
601
602 case TProb<T>::DISTTV:
603 for( size_t i = 0; i < p.size(); i++ )
604 result += abs(p[i] - q[i]);
605 result /= 2;
606 break;
607
608 case TProb<T>::DISTKL:
609 for( size_t i = 0; i < p.size(); i++ ) {
610 if( p[i] != 0.0 )
611 result += p[i] * (dai::log(p[i]) - dai::log(q[i]));
612 }
613 }
614 return result;
615 }
616
617
618 /// Writes a TProb<T> to an output stream
619 /** \relates TProb
620 */
621 template<typename T> std::ostream& operator<< (std::ostream& os, const TProb<T>& p) {
622 os << "[";
623 std::copy( p.p().begin(), p.p().end(), std::ostream_iterator<T>(os, " ") );
624 os << "]";
625 return os;
626 }
627
628
629 /// Returns the pointwise minimum of \a a and \a b
630 /** \relates TProb
631 * \pre <tt>this->size() == q.size()</tt>
632 */
633 template<typename T> TProb<T> min( const TProb<T> &a, const TProb<T> &b ) {
634 DAI_ASSERT( a.size() == b.size() );
635 TProb<T> result( a.size() );
636 for( size_t i = 0; i < a.size(); i++ )
637 if( a[i] < b[i] )
638 result[i] = a[i];
639 else
640 result[i] = b[i];
641 return result;
642 }
643
644
645 /// Returns the pointwise maximum of \a a and \a b
646 /** \relates TProb
647 * \pre <tt>this->size() == q.size()</tt>
648 */
649 template<typename T> TProb<T> max( const TProb<T> &a, const TProb<T> &b ) {
650 DAI_ASSERT( a.size() == b.size() );
651 TProb<T> result( a.size() );
652 for( size_t i = 0; i < a.size(); i++ )
653 if( a[i] > b[i] )
654 result[i] = a[i];
655 else
656 result[i] = b[i];
657 return result;
658 }
659
660
661 /// Represents a vector with entries of type dai::Real.
662 typedef TProb<Real> Prob;
663
664
665 } // end of namespace dai
666
667
668 #endif