Updated copyright headers
[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
16
17 #ifndef __defined_libdai_prob_h
18 #define __defined_libdai_prob_h
19
20
21 #include <cmath>
22 #include <vector>
23 #include <ostream>
24 #include <cassert>
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 /// Returns a const reference to the vector
89 const std::vector<T> & p() const { return _p; }
90
91 /// Returns a reference to the vector
92 std::vector<T> & p() { return _p; }
93
94 /// Returns a copy of the i'th entry
95 T operator[]( size_t i ) const {
96 #ifdef DAI_DEBUG
97 return _p.at(i);
98 #else
99 return _p[i];
100 #endif
101 }
102
103 /// Returns reference to the i'th entry
104 T& operator[]( size_t i ) { return _p[i]; }
105
106 /// Returns iterator pointing to first entry
107 iterator begin() { return _p.begin(); }
108
109 /// Returns const iterator pointing to first entry
110 const_iterator begin() const { return _p.begin(); }
111
112 /// Returns iterator pointing beyond last entry
113 iterator end() { return _p.end(); }
114
115 /// Returns const iterator pointing beyond last entry
116 const_iterator end() const { return _p.end(); }
117
118 /// Sets all entries to x
119 TProb<T> & fill(T x) {
120 std::fill( _p.begin(), _p.end(), x );
121 return *this;
122 }
123
124 /// Draws all entries i.i.d. from a uniform distribution on [0,1)
125 TProb<T> & randomize() {
126 std::generate(_p.begin(), _p.end(), rnd_uniform);
127 return *this;
128 }
129
130 /// Returns length of the vector, i.e., the number of entries
131 size_t size() const {
132 return _p.size();
133 }
134
135 /// Sets entries that are smaller than epsilon to 0
136 TProb<T>& makeZero( Real epsilon ) {
137 for( size_t i = 0; i < size(); i++ )
138 if( fabs(_p[i]) < epsilon )
139 _p[i] = 0;
140 return *this;
141 }
142
143 /// Set all entries to 1.0/size()
144 TProb<T>& setUniform () {
145 fill(1.0/size());
146 return *this;
147 }
148
149 /// Sets entries that are smaller than epsilon to epsilon
150 TProb<T>& makePositive( Real epsilon ) {
151 for( size_t i = 0; i < size(); i++ )
152 if( (0 < (Real)_p[i]) && ((Real)_p[i] < epsilon) )
153 _p[i] = epsilon;
154 return *this;
155 }
156
157 /// Multiplies each entry with scalar x
158 TProb<T>& operator*= (T x) {
159 std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::multiplies<T>(), x) );
160 return *this;
161 }
162
163 /// Returns product of *this with scalar x
164 TProb<T> operator* (T x) const {
165 TProb<T> prod( *this );
166 prod *= x;
167 return prod;
168 }
169
170 /// Divides each entry by scalar x
171 TProb<T>& operator/= (T x) {
172 DAI_DEBASSERT( x != 0.0 );
173 std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::divides<T>(), x ) );
174 return *this;
175 }
176
177 /// Returns quotient of *this and scalar x
178 TProb<T> operator/ (T x) const {
179 TProb<T> quot( *this );
180 quot /= x;
181 return quot;
182 }
183
184 /// Adds scalar x to each entry
185 TProb<T>& operator+= (T x) {
186 std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::plus<T>(), x ) );
187 return *this;
188 }
189
190 /// Returns sum of *this and scalar x
191 TProb<T> operator+ (T x) const {
192 TProb<T> sum( *this );
193 sum += x;
194 return sum;
195 }
196
197 /// Subtracts scalar x from each entry
198 TProb<T>& operator-= (T x) {
199 std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::minus<T>(), x ) );
200 return *this;
201 }
202
203 /// Returns difference of *this and scalar x
204 TProb<T> operator- (T x) const {
205 TProb<T> diff( *this );
206 diff -= x;
207 return diff;
208 }
209
210 /// Lexicographical comparison (sizes should be identical)
211 bool operator<= (const TProb<T> & q) const {
212 DAI_DEBASSERT( size() == q.size() );
213 for( size_t i = 0; i < size(); i++ )
214 if( !(_p[i] <= q[i]) )
215 return false;
216 return true;
217 }
218
219 /// Pointwise multiplication with q (sizes should be identical)
220 TProb<T>& operator*= (const TProb<T> & q) {
221 DAI_DEBASSERT( size() == q.size() );
222 std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::multiplies<T>() );
223 return *this;
224 }
225
226 /// Return product of *this with q (sizes should be identical)
227 TProb<T> operator* (const TProb<T> & q) const {
228 DAI_DEBASSERT( size() == q.size() );
229 TProb<T> prod( *this );
230 prod *= q;
231 return prod;
232 }
233
234 /// Pointwise addition with q (sizes should be identical)
235 TProb<T>& operator+= (const TProb<T> & q) {
236 DAI_DEBASSERT( size() == q.size() );
237 std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::plus<T>() );
238 return *this;
239 }
240
241 /// Returns sum of *this and q (sizes should be identical)
242 TProb<T> operator+ (const TProb<T> & q) const {
243 DAI_DEBASSERT( size() == q.size() );
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 DAI_DEBASSERT( size() == q.size() );
252 std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::minus<T>() );
253 return *this;
254 }
255
256 /// Return *this minus q (sizes should be identical)
257 TProb<T> operator- (const TProb<T> & q) const {
258 DAI_DEBASSERT( size() == q.size() );
259 TProb<T> diff( *this );
260 diff -= q;
261 return diff;
262 }
263
264 /// Pointwise division by q, where division by 0 yields 0 (sizes should be identical)
265 TProb<T>& operator/= (const TProb<T> & q) {
266 DAI_DEBASSERT( size() == q.size() );
267 for( size_t i = 0; i < size(); i++ ) {
268 if( q[i] == 0.0 )
269 _p[i] = 0.0;
270 else
271 _p[i] /= q[i];
272 }
273 return *this;
274 }
275
276 /// Pointwise division by q, where division by 0 yields +Inf (sizes should be identical)
277 TProb<T>& divide (const TProb<T> & q) {
278 DAI_DEBASSERT( size() == q.size() );
279 std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), std::divides<T>() );
280 return *this;
281 }
282
283 /// Returns quotient of *this with q (sizes should be identical)
284 TProb<T> operator/ (const TProb<T> & q) const {
285 DAI_DEBASSERT( size() == q.size() );
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 sum() const {
380 T Z = std::accumulate( _p.begin(), _p.end(), (T)0 );
381 return Z;
382 }
383
384 /// Return sum of absolute value of all entries
385 T sumAbs() const {
386 T s = 0;
387 for( size_t i = 0; i < size(); i++ )
388 s += fabs( (Real) _p[i] );
389 return s;
390 }
391
392 /// Returns maximum absolute value of all entries
393 T maxAbs() const {
394 T Z = 0;
395 for( size_t i = 0; i < size(); i++ ) {
396 Real mag = fabs( (Real) _p[i] );
397 if( mag > Z )
398 Z = mag;
399 }
400 return Z;
401 }
402
403 /// Returns maximum value of all entries
404 T max() const {
405 T Z = *std::max_element( _p.begin(), _p.end() );
406 return Z;
407 }
408
409 /// Returns minimum value of all entries
410 T min() const {
411 T Z = *std::min_element( _p.begin(), _p.end() );
412 return Z;
413 }
414
415 /// Returns {arg,}maximum value
416 std::pair<size_t,T> argmax() const {
417 T max = _p[0];
418 size_t arg = 0;
419 for( size_t i = 1; i < size(); i++ ) {
420 if( _p[i] > max ) {
421 max = _p[i];
422 arg = i;
423 }
424 }
425 return std::make_pair(arg,max);
426 }
427
428 /// Normalizes vector using the specified norm
429 T normalize( NormType norm=NORMPROB ) {
430 T Z = 0.0;
431 if( norm == NORMPROB )
432 Z = sum();
433 else if( norm == NORMLINF )
434 Z = maxAbs();
435 if( Z == 0.0 )
436 DAI_THROW(NOT_NORMALIZABLE);
437 else
438 *this /= Z;
439 return Z;
440 }
441
442 /// Returns normalized copy of *this, using the specified norm
443 TProb<T> normalized( NormType norm = NORMPROB ) const {
444 TProb<T> result(*this);
445 result.normalize( norm );
446 return result;
447 }
448
449 /// Returns true if one or more entries are NaN
450 bool hasNaNs() const {
451 bool foundnan = false;
452 for( typename std::vector<T>::const_iterator x = _p.begin(); x != _p.end(); x++ )
453 if( isnan( *x ) ) {
454 foundnan = true;
455 break;
456 }
457 return foundnan;
458 }
459
460 /// Returns true if one or more entries are negative
461 bool hasNegatives() const {
462 return (std::find_if( _p.begin(), _p.end(), std::bind2nd( std::less<Real>(), 0.0 ) ) != _p.end());
463 }
464
465 /// Returns entropy of *this
466 Real entropy() const {
467 Real S = 0.0;
468 for( size_t i = 0; i < size(); i++ )
469 S -= (_p[i] == 0 ? 0 : _p[i] * std::log(_p[i]));
470 return S;
471 }
472
473 /// Returns a random index, according to the (normalized) distribution described by *this
474 size_t draw() {
475 double x = rnd_uniform() * sum();
476 T s = 0;
477 for( size_t i = 0; i < size(); i++ ) {
478 s += _p[i];
479 if( s > x )
480 return i;
481 }
482 return( size() - 1 );
483 }
484 };
485
486
487 /// Returns distance of p and q (sizes should be identical), measured using distance measure dt
488 /** \relates TProb
489 */
490 template<typename T> Real dist( const TProb<T> &p, const TProb<T> &q, typename TProb<T>::DistType dt ) {
491 DAI_DEBASSERT( p.size() == q.size() );
492 Real result = 0.0;
493 switch( dt ) {
494 case TProb<T>::DISTL1:
495 for( size_t i = 0; i < p.size(); i++ )
496 result += fabs((Real)p[i] - (Real)q[i]);
497 break;
498
499 case TProb<T>::DISTLINF:
500 for( size_t i = 0; i < p.size(); i++ ) {
501 Real z = fabs((Real)p[i] - (Real)q[i]);
502 if( z > result )
503 result = z;
504 }
505 break;
506
507 case TProb<T>::DISTTV:
508 for( size_t i = 0; i < p.size(); i++ )
509 result += fabs((Real)p[i] - (Real)q[i]);
510 result *= 0.5;
511 break;
512
513 case TProb<T>::DISTKL:
514 for( size_t i = 0; i < p.size(); i++ ) {
515 if( p[i] != 0.0 )
516 result += p[i] * (std::log(p[i]) - std::log(q[i]));
517 }
518 }
519 return result;
520 }
521
522
523 /// Writes a TProb<T> to an output stream
524 /** \relates TProb
525 */
526 template<typename T> std::ostream& operator<< (std::ostream& os, const TProb<T>& P) {
527 os << "[";
528 std::copy( P.p().begin(), P.p().end(), std::ostream_iterator<T>(os, " ") );
529 os << "]";
530 return os;
531 }
532
533
534 /// Returns the TProb<T> containing the pointwise minimum of a and b (which should have equal size)
535 /** \relates TProb
536 */
537 template<typename T> TProb<T> min( const TProb<T> &a, const TProb<T> &b ) {
538 assert( a.size() == b.size() );
539 TProb<T> result( a.size() );
540 for( size_t i = 0; i < a.size(); i++ )
541 if( a[i] < b[i] )
542 result[i] = a[i];
543 else
544 result[i] = b[i];
545 return result;
546 }
547
548
549 /// Returns the TProb<T> containing the pointwise maximum of a and b (which should have equal size)
550 /** \relates TProb
551 */
552 template<typename T> TProb<T> max( const TProb<T> &a, const TProb<T> &b ) {
553 assert( a.size() == b.size() );
554 TProb<T> result( a.size() );
555 for( size_t i = 0; i < a.size(); i++ )
556 if( a[i] > b[i] )
557 result[i] = a[i];
558 else
559 result[i] = b[i];
560 return result;
561 }
562
563
564 /// Represents a vector with entries of type Real.
565 typedef TProb<Real> Prob;
566
567
568 } // end of namespace dai
569
570
571 #endif