Miscellaneuos changes thanks to Dan Preston
[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 #include <dai/exceptions.h>
41
42
43 namespace dai {
44
45
46 /// Represents a vector with entries of type \a T.
47 /** A TProb<T> is a std::vector<T> with an interface designed for dealing with probability mass functions.
48 * It is mainly used for representing measures on a finite outcome space, e.g., the probability
49 * distribution of a discrete random variable.
50 * \tparam T Should be a scalar that is castable from and to double and should support elementary arithmetic operations.
51 */
52 template <typename T> class TProb {
53 private:
54 /// The vector
55 std::vector<T> _p;
56
57 public:
58 /// Enumerates different ways of normalizing a probability measure.
59 /**
60 * - NORMPROB means that the sum of all entries should be 1;
61 * - NORMLINF means that the maximum absolute value of all entries should be 1.
62 */
63 typedef enum { NORMPROB, NORMLINF } NormType;
64 /// Enumerates different distance measures between probability measures.
65 /**
66 * - DISTL1 is the L-1 distance (sum of absolute values of pointwise difference);
67 * - DISTLINF is the L-inf distance (maximum absolute value of pointwise difference);
68 * - DISTTV is the Total Variation distance;
69 * - DISTKL is the Kullback-Leibler distance.
70 */
71 typedef enum { DISTL1, DISTLINF, DISTTV, DISTKL } DistType;
72
73 /// Default constructor
74 TProb() : _p() {}
75
76 /// Construct uniform distribution over n outcomes, i.e., a vector of length n with each entry set to 1/n
77 explicit TProb( size_t n ) : _p(std::vector<T>(n, 1.0 / n)) {}
78
79 /// Construct vector of length n with each entry set to p
80 explicit TProb( size_t n, Real p ) : _p(n, (T)p) {}
81
82 /// Construct vector from a range
83 /** \tparam Iterator Iterates over instances that can be cast to T.
84 * \param begin Points to first instance to be added.
85 * \param end Points just beyond last instance to be added.
86 * \param sizeHint For efficiency, the number of entries can be speficied by sizeHint.
87 */
88 template <typename Iterator>
89 TProb( Iterator begin, Iterator end, size_t sizeHint=0 ) : _p() {
90 _p.reserve( sizeHint );
91 _p.insert( _p.begin(), begin, end );
92 }
93
94 /// Returns a const reference to the vector
95 const std::vector<T> & p() const { return _p; }
96
97 /// Returns a reference to the vector
98 std::vector<T> & p() { return _p; }
99
100 /// Returns a copy of the i'th entry
101 T operator[]( size_t i ) const {
102 #ifdef DAI_DEBUG
103 return _p.at(i);
104 #else
105 return _p[i];
106 #endif
107 }
108
109 /// Returns reference to the i'th entry
110 T& operator[]( size_t i ) { return _p[i]; }
111
112 /// Sets all entries to x
113 TProb<T> & fill(T x) {
114 std::fill( _p.begin(), _p.end(), x );
115 return *this;
116 }
117
118 /// Draws all entries i.i.d. from a uniform distribution on [0,1)
119 TProb<T> & randomize() {
120 std::generate(_p.begin(), _p.end(), rnd_uniform);
121 return *this;
122 }
123
124 /// Returns length of the vector, i.e., the number of entries
125 size_t size() const {
126 return _p.size();
127 }
128
129 /// Sets entries that are smaller than epsilon to 0
130 TProb<T>& makeZero( Real epsilon ) {
131 for( size_t i = 0; i < size(); i++ )
132 if( fabs(_p[i]) < epsilon )
133 _p[i] = 0;
134 return *this;
135 }
136
137 /// Sets entries that are smaller than epsilon to epsilon
138 TProb<T>& makePositive( Real epsilon ) {
139 for( size_t i = 0; i < size(); i++ )
140 if( (0 < (Real)_p[i]) && ((Real)_p[i] < epsilon) )
141 _p[i] = epsilon;
142 return *this;
143 }
144
145 /// Multiplies each entry with scalar x
146 TProb<T>& operator*= (T x) {
147 std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::multiplies<T>(), x) );
148 return *this;
149 }
150
151 /// Returns product of *this with scalar x
152 TProb<T> operator* (T x) const {
153 TProb<T> prod( *this );
154 prod *= x;
155 return prod;
156 }
157
158 /// Divides each entry by scalar x
159 TProb<T>& operator/= (T x) {
160 #ifdef DAI_DEBUG
161 assert( x != 0.0 );
162 #endif
163 std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::divides<T>(), x ) );
164 return *this;
165 }
166
167 /// Returns quotient of *this and scalar x
168 TProb<T> operator/ (T x) const {
169 TProb<T> quot( *this );
170 quot /= x;
171 return quot;
172 }
173
174 /// Adds scalar x to each entry
175 TProb<T>& operator+= (T x) {
176 std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::plus<T>(), x ) );
177 return *this;
178 }
179
180 /// Returns sum of *this and scalar x
181 TProb<T> operator+ (T x) const {
182 TProb<T> sum( *this );
183 sum += x;
184 return sum;
185 }
186
187 /// Subtracts scalar x from each entry
188 TProb<T>& operator-= (T x) {
189 std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::minus<T>(), x ) );
190 return *this;
191 }
192
193 /// Returns difference of *this and scalar x
194 TProb<T> operator- (T x) const {
195 TProb<T> diff( *this );
196 diff -= x;
197 return diff;
198 }
199
200 /// Lexicographical comparison (sizes should be identical)
201 bool operator<= (const TProb<T> & q) const {
202 #ifdef DAI_DEBUG
203 assert( size() == q.size() );
204 #endif
205 for( size_t i = 0; i < size(); i++ )
206 if( !(_p[i] <= q[i]) )
207 return false;
208 return true;
209 }
210
211 /// Pointwise multiplication with q (sizes should be identical)
212 TProb<T>& operator*= (const TProb<T> & q) {
213 #ifdef DAI_DEBUG
214 assert( size() == q.size() );
215 #endif
216 std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::multiplies<T>() );
217 return *this;
218 }
219
220 /// Return product of *this with q (sizes should be identical)
221 TProb<T> operator* (const TProb<T> & q) const {
222 #ifdef DAI_DEBUG
223 assert( size() == q.size() );
224 #endif
225 TProb<T> prod( *this );
226 prod *= q;
227 return prod;
228 }
229
230 /// Pointwise addition with q (sizes should be identical)
231 TProb<T>& operator+= (const TProb<T> & q) {
232 #ifdef DAI_DEBUG
233 assert( size() == q.size() );
234 #endif
235 std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::plus<T>() );
236 return *this;
237 }
238
239 /// Returns sum of *this and q (sizes should be identical)
240 TProb<T> operator+ (const TProb<T> & q) const {
241 #ifdef DAI_DEBUG
242 assert( size() == q.size() );
243 #endif
244 TProb<T> sum( *this );
245 sum += q;
246 return sum;
247 }
248
249 /// Pointwise subtraction of q (sizes should be identical)
250 TProb<T>& operator-= (const TProb<T> & q) {
251 #ifdef DAI_DEBUG
252 assert( size() == q.size() );
253 #endif
254 std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::minus<T>() );
255 return *this;
256 }
257
258 /// Return *this minus q (sizes should be identical)
259 TProb<T> operator- (const TProb<T> & q) const {
260 #ifdef DAI_DEBUG
261 assert( size() == q.size() );
262 #endif
263 TProb<T> diff( *this );
264 diff -= q;
265 return diff;
266 }
267
268 /// Pointwise division by q, where division by 0 yields 0 (sizes should be identical)
269 TProb<T>& operator/= (const TProb<T> & q) {
270 #ifdef DAI_DEBUG
271 assert( size() == q.size() );
272 #endif
273 for( size_t i = 0; i < size(); i++ ) {
274 if( q[i] == 0.0 )
275 _p[i] = 0.0;
276 else
277 _p[i] /= q[i];
278 }
279 return *this;
280 }
281
282 /// Pointwise division by q, where division by 0 yields +Inf (sizes should be identical)
283 TProb<T>& divide (const TProb<T> & q) {
284 #ifdef DAI_DEBUG
285 assert( size() == q.size() );
286 #endif
287 std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::divides<T>() );
288 return *this;
289 }
290
291 /// Returns quotient of *this with q (sizes should be identical)
292 TProb<T> operator/ (const TProb<T> & q) const {
293 #ifdef DAI_DEBUG
294 assert( size() == q.size() );
295 #endif
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 totalSum() const {
390 T Z = std::accumulate( _p.begin(), _p.end(), (T)0 );
391 return Z;
392 }
393
394 /// Returns maximum absolute value of all entries
395 T maxAbs() const {
396 T Z = 0;
397 for( size_t i = 0; i < size(); i++ ) {
398 Real mag = fabs( (Real) _p[i] );
399 if( mag > Z )
400 Z = mag;
401 }
402 return Z;
403 }
404
405 /// Returns maximum value of all entries
406 T maxVal() const {
407 T Z = *std::max_element( _p.begin(), _p.end() );
408 return Z;
409 }
410
411 /// Returns minimum value of all entries
412 T minVal() const {
413 T Z = *std::min_element( _p.begin(), _p.end() );
414 return Z;
415 }
416
417 /// Normalizes vector using the specified norm
418 T normalize( NormType norm=NORMPROB ) {
419 T Z = 0.0;
420 if( norm == NORMPROB )
421 Z = totalSum();
422 else if( norm == NORMLINF )
423 Z = maxAbs();
424 if( Z == 0.0 )
425 DAI_THROW(NOT_NORMALIZABLE);
426 else
427 *this /= Z;
428 return Z;
429 }
430
431 /// Returns normalized copy of *this, using the specified norm
432 TProb<T> normalized( NormType norm = NORMPROB ) const {
433 TProb<T> result(*this);
434 result.normalize( norm );
435 return result;
436 }
437
438 /// Returns true if one or more entries are NaN
439 bool hasNaNs() const {
440 bool foundnan = false;
441 for( typename std::vector<T>::const_iterator x = _p.begin(); x != _p.end(); x++ )
442 if( isnan( *x ) ) {
443 foundnan = true;
444 break;
445 }
446 return foundnan;
447 }
448
449 /// Returns true if one or more entries are negative
450 bool hasNegatives() const {
451 return (std::find_if( _p.begin(), _p.end(), std::bind2nd( std::less<Real>(), 0.0 ) ) != _p.end());
452 }
453
454 /// Returns entropy of *this
455 Real entropy() const {
456 Real S = 0.0;
457 for( size_t i = 0; i < size(); i++ )
458 S -= (_p[i] == 0 ? 0 : _p[i] * std::log(_p[i]));
459 return S;
460 }
461
462 /// Returns a random index, according to the (normalized) distribution described by *this
463 size_t draw() {
464 double x = rnd_uniform() * totalSum();
465 T s = 0;
466 for( size_t i = 0; i < size(); i++ ) {
467 s += _p[i];
468 if( s > x )
469 return i;
470 }
471 return( size() - 1 );
472 }
473 };
474
475
476 /// Returns distance of p and q (sizes should be identical), measured using distance measure dt
477 /** \relates TProb
478 */
479 template<typename T> Real dist( const TProb<T> &p, const TProb<T> &q, typename TProb<T>::DistType dt ) {
480 #ifdef DAI_DEBUG
481 assert( p.size() == q.size() );
482 #endif
483 Real result = 0.0;
484 switch( dt ) {
485 case TProb<T>::DISTL1:
486 for( size_t i = 0; i < p.size(); i++ )
487 result += fabs((Real)p[i] - (Real)q[i]);
488 break;
489
490 case TProb<T>::DISTLINF:
491 for( size_t i = 0; i < p.size(); i++ ) {
492 Real z = fabs((Real)p[i] - (Real)q[i]);
493 if( z > result )
494 result = z;
495 }
496 break;
497
498 case TProb<T>::DISTTV:
499 for( size_t i = 0; i < p.size(); i++ )
500 result += fabs((Real)p[i] - (Real)q[i]);
501 result *= 0.5;
502 break;
503
504 case TProb<T>::DISTKL:
505 for( size_t i = 0; i < p.size(); i++ ) {
506 if( p[i] != 0.0 )
507 result += p[i] * (std::log(p[i]) - std::log(q[i]));
508 }
509 }
510 return result;
511 }
512
513
514 /// Writes a TProb<T> to an output stream
515 /** \relates TProb
516 */
517 template<typename T> std::ostream& operator<< (std::ostream& os, const TProb<T>& P) {
518 os << "[";
519 std::copy( P.p().begin(), P.p().end(), std::ostream_iterator<T>(os, " ") );
520 os << "]";
521 return os;
522 }
523
524
525 /// Returns the TProb<T> containing the pointwise minimum of a and b (which should have equal size)
526 /** \relates TProb
527 */
528 template<typename T> TProb<T> min( const TProb<T> &a, const TProb<T> &b ) {
529 assert( a.size() == b.size() );
530 TProb<T> result( a.size() );
531 for( size_t i = 0; i < a.size(); i++ )
532 if( a[i] < b[i] )
533 result[i] = a[i];
534 else
535 result[i] = b[i];
536 return result;
537 }
538
539
540 /// Returns the TProb<T> containing the pointwise maximum of a and b (which should have equal size)
541 /** \relates TProb
542 */
543 template<typename T> TProb<T> max( const TProb<T> &a, const TProb<T> &b ) {
544 assert( a.size() == b.size() );
545 TProb<T> result( a.size() );
546 for( size_t i = 0; i < a.size(); i++ )
547 if( a[i] > b[i] )
548 result[i] = a[i];
549 else
550 result[i] = b[i];
551 return result;
552 }
553
554
555 /// Represents a vector with entries of type Real.
556 typedef TProb<Real> Prob;
557
558
559 } // end of namespace dai
560
561
562 #endif