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