Improved documentation of include/dai/enum.h and include/dai/factor.h
[libdai.git] / include / dai / prob.h
1 /* Copyright (C) 2006-2008 Joris Mooij [joris dot mooij at tuebingen dot mpg dot de]
2 Radboud University Nijmegen, The Netherlands /
3 Max Planck Institute for Biological Cybernetics, Germany
4
5 This file is part of libDAI.
6
7 libDAI is free software; you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation; either version 2 of the License, or
10 (at your option) any later version.
11
12 libDAI is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with libDAI; if not, write to the Free Software
19 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
20 */
21
22
23 /// \file
24 /// \brief Defines TProb<T> and Prob classes
25 /// \todo Rename to Vector<T>
26
27
28 #ifndef __defined_libdai_prob_h
29 #define __defined_libdai_prob_h
30
31
32 #include <cmath>
33 #include <vector>
34 #include <ostream>
35 #include <cassert>
36 #include <algorithm>
37 #include <numeric>
38 #include <functional>
39 #include <dai/util.h>
40
41
42 namespace dai {
43
44
45 /// Represents a vector with entries of type \a T.
46 /** A TProb<T> is a std::vector<T> with an interface designed for dealing with probability mass functions.
47 * It is mainly used for representing measures on a finite outcome space, e.g., the probability
48 * distribution of a discrete random variable.
49 * \tparam T Should be a scalar that is castable from and to double and should support elementary arithmetic operations.
50 */
51 template <typename T> class TProb {
52 private:
53 /// The vector
54 std::vector<T> _p;
55
56 public:
57 /// Enumerates different ways of normalizing a probability measure.
58 /**
59 * - NORMPROB means that the sum of all entries should be 1;
60 * - NORMLINF means that the maximum absolute value of all entries should be 1.
61 */
62 typedef enum { NORMPROB, NORMLINF } NormType;
63 /// Enumerates different distance measures between probability measures.
64 /**
65 * - DISTL1 is the L-1 distance (sum of absolute values of pointwise difference);
66 * - DISTLINF is the L-inf distance (maximum absolute value of pointwise difference);
67 * - DISTTV is the Total Variation distance;
68 * - DISTKL is the Kullback-Leibler distance.
69 */
70 typedef enum { DISTL1, DISTLINF, DISTTV, DISTKL } DistType;
71
72 /// Default constructor
73 TProb() : _p() {}
74
75 /// Construct uniform distribution over n outcomes, i.e., a vector of length n with each entry set to 1/n
76 explicit TProb( size_t n ) : _p(std::vector<T>(n, 1.0 / n)) {}
77
78 /// Construct vector of length n with each entry set to p
79 explicit TProb( size_t n, Real p ) : _p(n, (T)p) {}
80
81 /// Construct vector of length n by copying the elements between p and p+n
82 TProb( size_t n, const Real* p ) : _p(p, p + n ) {}
83
84 /// Returns a const reference to the vector
85 const std::vector<T> & p() const { return _p; }
86
87 /// Returns a reference to the vector
88 std::vector<T> & p() { return _p; }
89
90 /// Returns a copy of the i'th entry
91 T operator[]( size_t i ) const {
92 #ifdef DAI_DEBUG
93 return _p.at(i);
94 #else
95 return _p[i];
96 #endif
97 }
98
99 /// Returns reference to the i'th entry
100 T& operator[]( size_t i ) { return _p[i]; }
101
102 /// Sets all entries to x
103 TProb<T> & fill(T x) {
104 std::fill( _p.begin(), _p.end(), x );
105 return *this;
106 }
107
108 /// Draws all entries i.i.d. from a uniform distribution on [0,1)
109 TProb<T> & randomize() {
110 std::generate(_p.begin(), _p.end(), rnd_uniform);
111 return *this;
112 }
113
114 /// Returns length of the vector, i.e., the number of entries
115 size_t size() const {
116 return _p.size();
117 }
118
119 /// Sets entries that are smaller than epsilon to 0
120 TProb<T>& makeZero( Real epsilon ) {
121 for( size_t i = 0; i < size(); i++ )
122 if( fabs(_p[i]) < epsilon )
123 _p[i] = 0;
124 return *this;
125 }
126
127 /// Sets entries that are smaller than epsilon to epsilon
128 TProb<T>& makePositive( Real epsilon ) {
129 for( size_t i = 0; i < size(); i++ )
130 if( (0 < (Real)_p[i]) && ((Real)_p[i] < epsilon) )
131 _p[i] = epsilon;
132 return *this;
133 }
134
135 /// Multiplies each entry with scalar x
136 TProb<T>& operator*= (T x) {
137 std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::multiplies<T>(), x) );
138 return *this;
139 }
140
141 /// Returns product of *this with scalar x
142 TProb<T> operator* (T x) const {
143 TProb<T> prod( *this );
144 prod *= x;
145 return prod;
146 }
147
148 /// Divides each entry by scalar x
149 TProb<T>& operator/= (T x) {
150 #ifdef DAI_DEBUG
151 assert( x != 0.0 );
152 #endif
153 std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::divides<T>(), x ) );
154 return *this;
155 }
156
157 /// Returns quotient of *this and scalar x
158 TProb<T> operator/ (T x) const {
159 TProb<T> quot( *this );
160 quot /= x;
161 return quot;
162 }
163
164 /// Adds scalar x to each entry
165 TProb<T>& operator+= (T x) {
166 std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::plus<T>(), x ) );
167 return *this;
168 }
169
170 /// Returns sum of *this and scalar x
171 TProb<T> operator+ (T x) const {
172 TProb<T> sum( *this );
173 sum += x;
174 return sum;
175 }
176
177 /// Subtracts scalar x from each entry
178 TProb<T>& operator-= (T x) {
179 std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::minus<T>(), x ) );
180 return *this;
181 }
182
183 /// Returns difference of *this and scalar x
184 TProb<T> operator- (T x) const {
185 TProb<T> diff( *this );
186 diff -= x;
187 return diff;
188 }
189
190 /// Lexicographical comparison (sizes should be identical)
191 bool operator<= (const TProb<T> & q) const {
192 #ifdef DAI_DEBUG
193 assert( size() == q.size() );
194 #endif
195 for( size_t i = 0; i < size(); i++ )
196 if( !(_p[i] <= q[i]) )
197 return false;
198 return true;
199 }
200
201 /// Pointwise multiplication with q (sizes should be identical)
202 TProb<T>& operator*= (const TProb<T> & q) {
203 #ifdef DAI_DEBUG
204 assert( size() == q.size() );
205 #endif
206 std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::multiplies<T>() );
207 return *this;
208 }
209
210 /// Return product of *this with q (sizes should be identical)
211 TProb<T> operator* (const TProb<T> & q) const {
212 #ifdef DAI_DEBUG
213 assert( size() == q.size() );
214 #endif
215 TProb<T> prod( *this );
216 prod *= q;
217 return prod;
218 }
219
220 /// Pointwise addition with q (sizes should be identical)
221 TProb<T>& operator+= (const TProb<T> & q) {
222 #ifdef DAI_DEBUG
223 assert( size() == q.size() );
224 #endif
225 std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::plus<T>() );
226 return *this;
227 }
228
229 /// Returns sum of *this and q (sizes should be identical)
230 TProb<T> operator+ (const TProb<T> & q) const {
231 #ifdef DAI_DEBUG
232 assert( size() == q.size() );
233 #endif
234 TProb<T> sum( *this );
235 sum += q;
236 return sum;
237 }
238
239 /// Pointwise subtraction of q (sizes should be identical)
240 TProb<T>& operator-= (const TProb<T> & q) {
241 #ifdef DAI_DEBUG
242 assert( size() == q.size() );
243 #endif
244 std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::minus<T>() );
245 return *this;
246 }
247
248 /// Return *this minus q (sizes should be identical)
249 TProb<T> operator- (const TProb<T> & q) const {
250 #ifdef DAI_DEBUG
251 assert( size() == q.size() );
252 #endif
253 TProb<T> diff( *this );
254 diff -= q;
255 return diff;
256 }
257
258 /// Pointwise division by q, where division by 0 yields 0 (sizes should be identical)
259 TProb<T>& operator/= (const TProb<T> & q) {
260 #ifdef DAI_DEBUG
261 assert( size() == q.size() );
262 #endif
263 for( size_t i = 0; i < size(); i++ ) {
264 if( q[i] == 0.0 )
265 _p[i] = 0.0;
266 else
267 _p[i] /= q[i];
268 }
269 return *this;
270 }
271
272 /// Pointwise division by q, where division by 0 yields +Inf (sizes should be identical)
273 TProb<T>& divide (const TProb<T> & q) {
274 #ifdef DAI_DEBUG
275 assert( size() == q.size() );
276 #endif
277 std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::divides<T>() );
278 return *this;
279 }
280
281 /// Returns quotient of *this with q (sizes should be identical)
282 TProb<T> operator/ (const TProb<T> & q) const {
283 #ifdef DAI_DEBUG
284 assert( size() == q.size() );
285 #endif
286 TProb<T> quot( *this );
287 quot /= q;
288 return quot;
289 }
290
291 /// Returns pointwise inverse
292 /** If zero==true; uses 1/0==0, otherwise 1/0==Inf.
293 */
294 TProb<T> inverse(bool zero=true) const {
295 TProb<T> inv;
296 inv._p.reserve( size() );
297 if( zero )
298 for( size_t i = 0; i < size(); i++ )
299 inv._p.push_back( _p[i] == 0.0 ? 0.0 : 1.0 / _p[i] );
300 else
301 for( size_t i = 0; i < size(); i++ )
302 inv._p.push_back( 1.0 / _p[i] );
303 return inv;
304 }
305
306 /// Raises entries to the power a
307 TProb<T>& operator^= (Real a) {
308 if( a != 1.0 )
309 std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::ptr_fun<T, Real, T>(std::pow), a) );
310 return *this;
311 }
312
313 /// Returns *this raised to the power a
314 TProb<T> operator^ (Real a) const {
315 TProb<T> power(*this);
316 power ^= a;
317 return power;
318 }
319
320 /// Returns pointwise signum
321 TProb<T> sgn() const {
322 TProb<T> x;
323 x._p.reserve( size() );
324 for( size_t i = 0; i < size(); i++ ) {
325 T s = 0;
326 if( _p[i] > 0 )
327 s = 1;
328 else if( _p[i] < 0 )
329 s = -1;
330 x._p.push_back( s );
331 }
332 return x;
333 }
334
335 /// Returns pointwise absolute value
336 TProb<T> abs() const {
337 TProb<T> x;
338 x._p.reserve( size() );
339 for( size_t i = 0; i < size(); i++ )
340 x._p.push_back( _p[i] < 0 ? (-p[i]) : p[i] );
341 return x;
342 }
343
344 /// Applies exp pointwise
345 const TProb<T>& takeExp() {
346 std::transform( _p.begin(), _p.end(), _p.begin(), std::ptr_fun<T, T>(std::exp) );
347 return *this;
348 }
349
350 /// Applies log pointwise
351 /** If zero==true, uses log(0)==0; otherwise, log(0)==-Inf.
352 */
353 const TProb<T>& takeLog(bool zero=false) {
354 if( zero ) {
355 for( size_t i = 0; i < size(); i++ )
356 _p[i] = ( (_p[i] == 0.0) ? 0.0 : std::log( _p[i] ) );
357 } else
358 std::transform( _p.begin(), _p.end(), _p.begin(), std::ptr_fun<T, T>(std::log) );
359 return *this;
360 }
361
362 /// Returns pointwise exp
363 TProb<T> exp() const {
364 TProb<T> e(*this);
365 e.takeExp();
366 return e;
367 }
368
369 /// Returns pointwise log
370 /** If zero==true, uses log(0)==0; otherwise, log(0)==-Inf.
371 */
372 TProb<T> log(bool zero=false) const {
373 TProb<T> l(*this);
374 l.takeLog(zero);
375 return l;
376 }
377
378 /// Returns sum of all entries
379 T totalSum() const {
380 T Z = std::accumulate( _p.begin(), _p.end(), (T)0 );
381 return Z;
382 }
383
384 /// Returns maximum absolute value of all entries
385 T maxAbs() const {
386 T Z = 0;
387 for( size_t i = 0; i < size(); i++ ) {
388 Real mag = fabs( (Real) _p[i] );
389 if( mag > Z )
390 Z = mag;
391 }
392 return Z;
393 }
394
395 /// Returns maximum value of all entries
396 T maxVal() const {
397 T Z = *std::max_element( _p.begin(), _p.end() );
398 return Z;
399 }
400
401 /// Returns minimum value of all entries
402 T minVal() const {
403 T Z = *std::min_element( _p.begin(), _p.end() );
404 return Z;
405 }
406
407 /// Normalizes vector using the specified norm
408 T normalize( NormType norm=NORMPROB ) {
409 T Z = 0.0;
410 if( norm == NORMPROB )
411 Z = totalSum();
412 else if( norm == NORMLINF )
413 Z = maxAbs();
414 #ifdef DAI_DEBUG
415 assert( Z != 0.0 );
416 #endif
417 *this /= Z;
418 return Z;
419 }
420
421 /// Returns normalized copy of *this, using the specified norm
422 TProb<T> normalized( NormType norm = NORMPROB ) const {
423 TProb<T> result(*this);
424 result.normalize( norm );
425 return result;
426 }
427
428 /// Returns true if one or more entries are NaN
429 bool hasNaNs() const {
430 return (std::find_if( _p.begin(), _p.end(), isnan ) != _p.end());
431 }
432
433 /// Returns true if one or more entries are negative
434 bool hasNegatives() const {
435 return (std::find_if( _p.begin(), _p.end(), std::bind2nd( std::less<Real>(), 0.0 ) ) != _p.end());
436 }
437
438 /// Returns entropy of *this
439 Real entropy() const {
440 Real S = 0.0;
441 for( size_t i = 0; i < size(); i++ )
442 S -= (_p[i] == 0 ? 0 : _p[i] * std::log(_p[i]));
443 return S;
444 }
445 };
446
447
448 /// Returns distance of p and q (sizes should be identical), measured using distance measure dt
449 /** \relates TProb
450 */
451 template<typename T> Real dist( const TProb<T> &p, const TProb<T> &q, typename TProb<T>::DistType dt ) {
452 #ifdef DAI_DEBUG
453 assert( p.size() == q.size() );
454 #endif
455 Real result = 0.0;
456 switch( dt ) {
457 case TProb<T>::DISTL1:
458 for( size_t i = 0; i < p.size(); i++ )
459 result += fabs((Real)p[i] - (Real)q[i]);
460 break;
461
462 case TProb<T>::DISTLINF:
463 for( size_t i = 0; i < p.size(); i++ ) {
464 Real z = fabs((Real)p[i] - (Real)q[i]);
465 if( z > result )
466 result = z;
467 }
468 break;
469
470 case TProb<T>::DISTTV:
471 for( size_t i = 0; i < p.size(); i++ )
472 result += fabs((Real)p[i] - (Real)q[i]);
473 result *= 0.5;
474 break;
475
476 case TProb<T>::DISTKL:
477 for( size_t i = 0; i < p.size(); i++ ) {
478 if( p[i] != 0.0 )
479 result += p[i] * (std::log(p[i]) - std::log(q[i]));
480 }
481 }
482 return result;
483 }
484
485
486 /// Writes a TProb<T> to an output stream
487 /** \relates TProb
488 */
489 template<typename T> std::ostream& operator<< (std::ostream& os, const TProb<T>& P) {
490 os << "[";
491 std::copy( P.p().begin(), P.p().end(), std::ostream_iterator<T>(os, " ") );
492 os << "]";
493 return os;
494 }
495
496
497 /// Returns the TProb<T> containing the pointwise minimum of a and b (which should have equal size)
498 /** \relates TProb
499 */
500 template<typename T> TProb<T> min( const TProb<T> &a, const TProb<T> &b ) {
501 assert( a.size() == b.size() );
502 TProb<T> result( a.size() );
503 for( size_t i = 0; i < a.size(); i++ )
504 if( a[i] < b[i] )
505 result[i] = a[i];
506 else
507 result[i] = b[i];
508 return result;
509 }
510
511
512 /// Returns the TProb<T> containing the pointwise maximum of a and b (which should have equal size)
513 /** \relates TProb
514 */
515 template<typename T> TProb<T> max( const TProb<T> &a, const TProb<T> &b ) {
516 assert( a.size() == b.size() );
517 TProb<T> result( a.size() );
518 for( size_t i = 0; i < a.size(); i++ )
519 if( a[i] > b[i] )
520 result[i] = a[i];
521 else
522 result[i] = b[i];
523 return result;
524 }
525
526
527 /// Represents a vector with entries of type Real.
528 typedef TProb<Real> Prob;
529
530
531 } // end of namespace dai
532
533
534 #endif