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