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