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