Significant improvement of documentation
[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
26
27 #ifndef __defined_libdai_prob_h
28 #define __defined_libdai_prob_h
29
30
31 #include <cmath>
32 #include <vector>
33 #include <ostream>
34 #include <cassert>
35 #include <algorithm>
36 #include <numeric>
37 #include <functional>
38 #include <dai/util.h>
39
40
41 namespace dai {
42
43
44 /// Real number (alias for double, could be changed to long double if necessary)
45 typedef double Real;
46
47 template<typename T> class TProb;
48
49 /// Represents a probability measure, with entries of type Real.
50 typedef TProb<Real> Prob;
51
52
53 /// Represents a probability measure on a finite outcome space (i.e., corresponding to a discrete random variable).
54 /** It is implemented as a std::vector<T> but adds a convenient interface.
55 * It is not necessarily normalized at all times.
56 * \tparam T Should be castable from and to double.
57 */
58 template <typename T> class TProb {
59 private:
60 /// The probability measure
61 std::vector<T> _p;
62
63 public:
64 /// Enumerates different ways of normalizing a probability measure.
65 /**
66 * - NORMPROB means that the sum of all entries should be 1;
67 * - NORMLINF means that the maximum absolute value of all entries should be 1.
68 */
69 typedef enum { NORMPROB, NORMLINF } NormType;
70 /// Enumerates different distance measures between probability measures.
71 /**
72 * - DISTL1 is the L-1 distance (sum of absolute values of pointwise difference);
73 * - DISTLINF is the L-inf distance (maximum absolute value of pointwise difference);
74 * - DISTTV is the Total Variation distance;
75 * - DISTKL is the Kullback-Leibler distance.
76 */
77 typedef enum { DISTL1, DISTLINF, DISTTV, DISTKL } DistType;
78
79 /// Default constructor
80 TProb() : _p() {}
81
82 /// Construct uniform distribution of given length
83 explicit TProb( size_t n ) : _p(std::vector<T>(n, 1.0 / n)) {}
84
85 /// Construct from given length and initial value
86 TProb( size_t n, Real p ) : _p(n, (T)p) {}
87
88 /// Construct from given length and initial array
89 TProb( size_t n, const Real* p ) : _p(p, p + n ) {}
90
91 /// Returns a const reference to the probability vector
92 const std::vector<T> & p() const { return _p; }
93
94 /// Returns a reference to the probability vector
95 std::vector<T> & p() { return _p; }
96
97 /// Returns a copy of the i'th probability entry
98 T operator[]( size_t i ) const {
99 #ifdef DAI_DEBUG
100 return _p.at(i);
101 #else
102 return _p[i];
103 #endif
104 }
105
106 /// Returns a reference to the i'th probability entry
107 T& operator[]( size_t i ) { return _p[i]; }
108
109 /// Sets all elements to x
110 TProb<T> & fill(T x) {
111 std::fill( _p.begin(), _p.end(), x );
112 return *this;
113 }
114
115 /// Sets all elements to i.i.d. random numbers from a uniform[0,1) distribution
116 TProb<T> & randomize() {
117 std::generate(_p.begin(), _p.end(), rnd_uniform);
118 return *this;
119 }
120
121 /// Returns number of elements
122 size_t size() const {
123 return _p.size();
124 }
125
126 /// Sets entries that are smaller than epsilon to zero
127 TProb<T>& makeZero( Real epsilon ) {
128 for( size_t i = 0; i < size(); i++ )
129 if( fabs(_p[i]) < epsilon )
130 _p[i] = 0;
131 return *this;
132 }
133
134 /// Sets entries that are smaller than epsilon to epsilon
135 TProb<T>& makePositive (Real epsilon) {
136 for( size_t i = 0; i < size(); i++ )
137 if( (0 < (Real)_p[i]) && ((Real)_p[i] < epsilon) )
138 _p[i] = epsilon;
139 return *this;
140 }
141
142 /// Multiplies each entry with x
143 TProb<T>& operator*= (T x) {
144 std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::multiplies<T>(), x) );
145 return *this;
146 }
147
148 /// Returns product of *this with x
149 TProb<T> operator* (T x) const {
150 TProb<T> prod( *this );
151 prod *= x;
152 return prod;
153 }
154
155 /// Divides each entry by x
156 TProb<T>& operator/= (T x) {
157 #ifdef DAI_DEBUG
158 assert( x != 0.0 );
159 #endif
160 std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::divides<T>(), x ) );
161 return *this;
162 }
163
164 /// Returns quotient of *this and x
165 TProb<T> operator/ (T x) const {
166 TProb<T> quot( *this );
167 quot /= x;
168 return quot;
169 }
170
171 /// Adds x to each entry
172 TProb<T>& operator+= (T x) {
173 std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::plus<T>(), x ) );
174 return *this;
175 }
176
177 /// Returns sum of *this and x
178 TProb<T> operator+ (T x) const {
179 TProb<T> sum( *this );
180 sum += x;
181 return sum;
182 }
183
184 /// Subtracts x from each entry
185 TProb<T>& operator-= (T x) {
186 std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::minus<T>(), x ) );
187 return *this;
188 }
189
190 /// Returns difference of *this and x
191 TProb<T> operator- (T x) const {
192 TProb<T> diff( *this );
193 diff -= x;
194 return diff;
195 }
196
197 /// Pointwise comparison
198 bool operator<= (const TProb<T> & q) const {
199 #ifdef DAI_DEBUG
200 assert( size() == q.size() );
201 #endif
202 for( size_t i = 0; i < size(); i++ )
203 if( !(_p[i] <= q[i]) )
204 return false;
205 return true;
206 }
207
208 /// Pointwise multiplication with q
209 TProb<T>& operator*= (const TProb<T> & q) {
210 #ifdef DAI_DEBUG
211 assert( size() == q.size() );
212 #endif
213 std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::multiplies<T>() );
214 return *this;
215 }
216
217 /// Return product of *this with q
218 TProb<T> operator* (const TProb<T> & q) const {
219 #ifdef DAI_DEBUG
220 assert( size() == q.size() );
221 #endif
222 TProb<T> prod( *this );
223 prod *= q;
224 return prod;
225 }
226
227 /// Pointwise addition with q
228 TProb<T>& operator+= (const TProb<T> & q) {
229 #ifdef DAI_DEBUG
230 assert( size() == q.size() );
231 #endif
232 std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::plus<T>() );
233 return *this;
234 }
235
236 /// Return sum of *this and q
237 TProb<T> operator+ (const TProb<T> & q) const {
238 #ifdef DAI_DEBUG
239 assert( size() == q.size() );
240 #endif
241 TProb<T> sum( *this );
242 sum += q;
243 return sum;
244 }
245
246 /// Pointwise subtraction of q
247 TProb<T>& operator-= (const TProb<T> & q) {
248 #ifdef DAI_DEBUG
249 assert( size() == q.size() );
250 #endif
251 std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::minus<T>() );
252 return *this;
253 }
254
255 /// Return *this minus q
256 TProb<T> operator- (const TProb<T> & q) const {
257 #ifdef DAI_DEBUG
258 assert( size() == q.size() );
259 #endif
260 TProb<T> diff( *this );
261 diff -= q;
262 return diff;
263 }
264
265 /// Pointwise division by q, where division by zero yields zero
266 TProb<T>& operator/= (const TProb<T> & q) {
267 #ifdef DAI_DEBUG
268 assert( size() == q.size() );
269 #endif
270 for( size_t i = 0; i < size(); i++ ) {
271 if( q[i] == 0.0 )
272 _p[i] = 0.0;
273 else
274 _p[i] /= q[i];
275 }
276 return *this;
277 }
278
279 /// Pointwise division by q, where division by zero yields infinity
280 TProb<T>& divide (const TProb<T> & q) {
281 #ifdef DAI_DEBUG
282 assert( size() == q.size() );
283 #endif
284 std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::divides<T>() );
285 return *this;
286 }
287
288 /// Returns quotient of *this with q
289 TProb<T> operator/ (const TProb<T> & q) const {
290 #ifdef DAI_DEBUG
291 assert( size() == q.size() );
292 #endif
293 TProb<T> quot( *this );
294 quot /= q;
295 return quot;
296 }
297
298 /// Returns pointwise inverse
299 TProb<T> inverse(bool zero = false) 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] == 0.0 ? 0.0 : 1.0 / _p[i] );
305 else
306 for( size_t i = 0; i < size(); i++ ) {
307 #ifdef DAI_DEBUG
308 assert( _p[i] != 0.0 );
309 #endif
310 inv._p.push_back( 1.0 / _p[i] );
311 }
312 return inv;
313 }
314
315 /// Raises elements 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 const TProb<T>& takeLog() {
361 std::transform( _p.begin(), _p.end(), _p.begin(), std::ptr_fun<T, T>(std::log) );
362 return *this;
363 }
364
365 /// Applies log pointwise (defining log(0)=0)
366 const TProb<T>& takeLog0() {
367 for( size_t i = 0; i < size(); i++ )
368 _p[i] = ( (_p[i] == 0.0) ? 0.0 : std::log( _p[i] ) );
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 TProb<T> log() const {
381 TProb<T> l(*this);
382 l.takeLog();
383 return l;
384 }
385
386 /// Returns pointwise log (defining log(0)=0)
387 TProb<T> log0() const {
388 TProb<T> l0(*this);
389 l0.takeLog0();
390 return l0;
391 }
392
393 /// Returns distance of p and q, measured using dt
394 friend Real dist( const TProb<T> &p, const TProb<T> &q, DistType dt ) {
395 #ifdef DAI_DEBUG
396 assert( p.size() == q.size() );
397 #endif
398 Real result = 0.0;
399 switch( dt ) {
400 case DISTL1:
401 for( size_t i = 0; i < p.size(); i++ )
402 result += fabs((Real)p[i] - (Real)q[i]);
403 break;
404
405 case DISTLINF:
406 for( size_t i = 0; i < p.size(); i++ ) {
407 Real z = fabs((Real)p[i] - (Real)q[i]);
408 if( z > result )
409 result = z;
410 }
411 break;
412
413 case DISTTV:
414 for( size_t i = 0; i < p.size(); i++ )
415 result += fabs((Real)p[i] - (Real)q[i]);
416 result *= 0.5;
417 break;
418
419 case DISTKL:
420 for( size_t i = 0; i < p.size(); i++ ) {
421 if( p[i] != 0.0 )
422 result += p[i] * (std::log(p[i]) - std::log(q[i]));
423 }
424 }
425 return result;
426 }
427
428 /// Returns sum of all entries
429 T totalSum() const {
430 T Z = std::accumulate( _p.begin(), _p.end(), (T)0 );
431 return Z;
432 }
433
434 /// Returns maximum absolute value of entries
435 T maxAbs() const {
436 T Z = 0;
437 for( size_t i = 0; i < size(); i++ ) {
438 Real mag = fabs( (Real) _p[i] );
439 if( mag > Z )
440 Z = mag;
441 }
442 return Z;
443 }
444
445 /// Returns maximum value of entries
446 T maxVal() const {
447 T Z = *std::max_element( _p.begin(), _p.end() );
448 return Z;
449 }
450
451 /// Returns minimum value of entries
452 T minVal() const {
453 T Z = *std::min_element( _p.begin(), _p.end() );
454 return Z;
455 }
456
457 /// Normalizes using the specified norm
458 T normalize( NormType norm = NORMPROB ) {
459 T Z = 0.0;
460 if( norm == NORMPROB )
461 Z = totalSum();
462 else if( norm == NORMLINF )
463 Z = maxAbs();
464 #ifdef DAI_DEBUG
465 assert( Z != 0.0 );
466 #endif
467 *this /= Z;
468 return Z;
469 }
470
471 /// Returns normalized copy of *this, using the specified norm
472 TProb<T> normalized( NormType norm = NORMPROB ) const {
473 TProb<T> result(*this);
474 result.normalize( norm );
475 return result;
476 }
477
478 /// Returns true if one or more entries are NaN
479 bool hasNaNs() const {
480 return (std::find_if( _p.begin(), _p.end(), isnan ) != _p.end());
481 }
482
483 /// Returns true if one or more entries are negative
484 bool hasNegatives() const {
485 return (std::find_if( _p.begin(), _p.end(), std::bind2nd( std::less<Real>(), 0.0 ) ) != _p.end());
486 }
487
488 /// Returns true if one or more entries are non-positive (causes problems with logscale)
489 bool hasNonPositives() const {
490 return (std::find_if( _p.begin(), _p.end(), std::bind2nd( std::less_equal<Real>(), 0.0 ) ) != _p.end());
491 }
492
493 /// Returns entropy
494 Real entropy() const {
495 Real S = 0.0;
496 for( size_t i = 0; i < size(); i++ )
497 S -= xlogx(_p[i]);
498 return S;
499 }
500
501 /// Writes a TProb<T> to an output stream
502 friend std::ostream& operator<< (std::ostream& os, const TProb<T>& P) {
503 os << "[";
504 std::copy( P._p.begin(), P._p.end(), std::ostream_iterator<T>(os, " ") );
505 os << "]";
506 return os;
507 }
508
509 private:
510 /// Returns x*log(x), or 0 if x == 0
511 Real xlogx( Real x ) const { return( x == 0.0 ? 0.0 : x * std::log(x)); }
512 };
513
514
515 /// Returns TProb<T> containing the pointwise minimum of a and b (which should have equal size)
516 template<typename T> TProb<T> min( const TProb<T> &a, const TProb<T> &b ) {
517 assert( a.size() == b.size() );
518 TProb<T> result( a.size() );
519 for( size_t i = 0; i < a.size(); i++ )
520 if( a[i] < b[i] )
521 result[i] = a[i];
522 else
523 result[i] = b[i];
524 return result;
525 }
526
527
528 /// Returns TProb<T> containing the pointwise maximum of a and b (which should have equal size)
529 template<typename T> TProb<T> max( const TProb<T> &a, const TProb<T> &b ) {
530 assert( a.size() == b.size() );
531 TProb<T> result( a.size() );
532 for( size_t i = 0; i < a.size(); i++ )
533 if( a[i] > b[i] )
534 result[i] = a[i];
535 else
536 result[i] = b[i];
537 return result;
538 }
539
540
541 } // end of namespace dai
542
543
544 #endif