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