libDAI
prob.h
Go to the documentation of this file.
1 /* This file is part of libDAI - http://www.libdai.org/
2  *
3  * Copyright (c) 2006-2011, The libDAI authors. All rights reserved.
4  *
5  * Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
6  */
7 
8 
11 
12 
13 #ifndef __defined_libdai_prob_h
14 #define __defined_libdai_prob_h
15 
16 
17 #include <cmath>
18 #include <vector>
19 #include <ostream>
20 #include <algorithm>
21 #include <numeric>
22 #include <functional>
23 #include <dai/util.h>
24 #include <dai/exceptions.h>
25 
26 
27 namespace dai {
28 
29 
31 template<typename T> struct fo_id : public std::unary_function<T, T> {
33  T operator()( const T &x ) const {
34  return x;
35  }
36 };
37 
38 
40 template<typename T> struct fo_abs : public std::unary_function<T, T> {
42  T operator()( const T &x ) const {
43  if( x < (T)0 )
44  return -x;
45  else
46  return x;
47  }
48 };
49 
50 
52 template<typename T> struct fo_exp : public std::unary_function<T, T> {
54  T operator()( const T &x ) const {
55  return exp( x );
56  }
57 };
58 
59 
61 template<typename T> struct fo_log : public std::unary_function<T, T> {
63  T operator()( const T &x ) const {
64  return log( x );
65  }
66 };
67 
68 
70 template<typename T> struct fo_log0 : public std::unary_function<T, T> {
72  T operator()( const T &x ) const {
73  if( x )
74  return log( x );
75  else
76  return 0;
77  }
78 };
79 
80 
82 template<typename T> struct fo_inv : public std::unary_function<T, T> {
84  T operator()( const T &x ) const {
85  return 1 / x;
86  }
87 };
88 
89 
91 template<typename T> struct fo_inv0 : public std::unary_function<T, T> {
93  T operator()( const T &x ) const {
94  if( x )
95  return 1 / x;
96  else
97  return 0;
98  }
99 };
100 
101 
103 template<typename T> struct fo_plog0p : public std::unary_function<T, T> {
105  T operator()( const T &p ) const {
106  return p * dai::log0(p);
107  }
108 };
109 
110 
112 template<typename T> struct fo_divides0 : public std::binary_function<T, T, T> {
114  T operator()( const T &x, const T &y ) const {
115  if( y == (T)0 )
116  return (T)0;
117  else
118  return x / y;
119  }
120 };
121 
122 
124 template<typename T> struct fo_KL : public std::binary_function<T, T, T> {
126  T operator()( const T &p, const T &q ) const {
127  if( p == (T)0 )
128  return (T)0;
129  else
130  return p * (log(p) - log(q));
131  }
132 };
133 
134 
136 template<typename T> struct fo_Hellinger : public std::binary_function<T, T, T> {
138  T operator()( const T &p, const T &q ) const {
139  T x = sqrt(p) - sqrt(q);
140  return x * x;
141  }
142 };
143 
144 
146 template<typename T> struct fo_pow : public std::binary_function<T, T, T> {
148  T operator()( const T &x, const T &y ) const {
149  if( y != 1 )
150  return pow( x, y );
151  else
152  return x;
153  }
154 };
155 
156 
158 template<typename T> struct fo_max : public std::binary_function<T, T, T> {
160  T operator()( const T &x, const T &y ) const {
161  return (x > y) ? x : y;
162  }
163 };
164 
165 
167 template<typename T> struct fo_min : public std::binary_function<T, T, T> {
169  T operator()( const T &x, const T &y ) const {
170  return (x > y) ? y : x;
171  }
172 };
173 
174 
176 template<typename T> struct fo_absdiff : public std::binary_function<T, T, T> {
178  T operator()( const T &x, const T &y ) const {
179  return dai::abs( x - y );
180  }
181 };
182 
183 
185 
193 template <typename T>
194 class TProb {
195  public:
197  typedef std::vector<T> container_type;
198 
201 
202  private:
205 
206  public:
208 
209 
210  TProb() : _p() {}
211 
213  explicit TProb( size_t n ) : _p( n, (T)1 / n ) {}
214 
216  explicit TProb( size_t n, T p ) : _p( n, p ) {}
217 
219 
225  template <typename TIterator>
226  TProb( TIterator begin, TIterator end, size_t sizeHint ) : _p() {
227  _p.reserve( sizeHint );
228  _p.insert( _p.begin(), begin, end );
229  }
230 
232 
235  template <typename S>
236  TProb( const std::vector<S> &v ) : _p() {
237  _p.reserve( v.size() );
238  _p.insert( _p.begin(), v.begin(), v.end() );
239  }
241 
243  typedef typename container_type::const_iterator const_iterator;
245  typedef typename container_type::iterator iterator;
247  typedef typename container_type::const_reverse_iterator const_reverse_iterator;
249  typedef typename container_type::reverse_iterator reverse_iterator;
250 
252 
253 
254  iterator begin() { return _p.begin(); }
256  const_iterator begin() const { return _p.begin(); }
257 
259  iterator end() { return _p.end(); }
261  const_iterator end() const { return _p.end(); }
262 
264  reverse_iterator rbegin() { return _p.rbegin(); }
266  const_reverse_iterator rbegin() const { return _p.rbegin(); }
267 
269  reverse_iterator rend() { return _p.rend(); }
271  const_reverse_iterator rend() const { return _p.rend(); }
273 
275 
276  void resize( size_t sz ) {
277  _p.resize( sz );
278  }
280 
282 
283 
284  T get( size_t i ) const {
285 #ifdef DAI_DEBUG
286  return _p.at(i);
287 #else
288  return _p[i];
289 #endif
290  }
291 
293  void set( size_t i, T val ) {
294  DAI_DEBASSERT( i < _p.size() );
295  _p[i] = val;
296  }
298 
300 
301 
302  const container_type& p() const { return _p; }
303 
305  container_type& p() { return _p; }
306 
308  T operator[]( size_t i ) const { return get(i); }
309 
311  size_t size() const { return _p.size(); }
312 
314 
322  template<typename unOp> T accumulateSum( T init, unOp op ) const {
323  T t = op(init);
324  for( const_iterator it = begin(); it != end(); it++ )
325  t += op(*it);
326  return t;
327  }
328 
330 
338  template<typename unOp> T accumulateMax( T init, unOp op, bool minimize ) const {
339  T t = op(init);
340  if( minimize ) {
341  for( const_iterator it = begin(); it != end(); it++ )
342  t = std::min( t, op(*it) );
343  } else {
344  for( const_iterator it = begin(); it != end(); it++ )
345  t = std::max( t, op(*it) );
346  }
347  return t;
348  }
349 
351  T entropy() const { return -accumulateSum( (T)0, fo_plog0p<T>() ); }
352 
354  T max() const { return accumulateMax( (T)(-INFINITY), fo_id<T>(), false ); }
355 
357  T min() const { return accumulateMax( (T)INFINITY, fo_id<T>(), true ); }
358 
360  T sum() const { return accumulateSum( (T)0, fo_id<T>() ); }
361 
363  T sumAbs() const { return accumulateSum( (T)0, fo_abs<T>() ); }
364 
366  T maxAbs() const { return accumulateMax( (T)0, fo_abs<T>(), false ); }
367 
369  bool hasNaNs() const {
370  bool foundnan = false;
371  for( const_iterator x = _p.begin(); x != _p.end(); x++ )
372  if( dai::isnan( *x ) ) {
373  foundnan = true;
374  break;
375  }
376  return foundnan;
377  }
378 
380  bool hasNegatives() const {
381  return (std::find_if( _p.begin(), _p.end(), std::bind2nd( std::less<T>(), (T)0 ) ) != _p.end());
382  }
383 
385  std::pair<size_t,T> argmax() const {
386  T max = _p[0];
387  size_t arg = 0;
388  for( size_t i = 1; i < size(); i++ ) {
389  if( _p[i] > max ) {
390  max = _p[i];
391  arg = i;
392  }
393  }
394  return std::make_pair( arg, max );
395  }
396 
398  size_t draw() {
399  Real x = rnd_uniform() * sum();
400  T s = 0;
401  for( size_t i = 0; i < size(); i++ ) {
402  s += get(i);
403  if( s > x )
404  return i;
405  }
406  return( size() - 1 );
407  }
408 
410 
412  bool operator<( const this_type& q ) const {
413  DAI_DEBASSERT( size() == q.size() );
414  return lexicographical_compare( begin(), end(), q.begin(), q.end() );
415  }
416 
418  bool operator==( const this_type& q ) const {
419  if( size() != q.size() )
420  return false;
421  return p() == q.p();
422  }
424 
426 
427 
428  template<typename unaryOp> this_type pwUnaryTr( unaryOp op ) const {
429  this_type r;
430  r._p.reserve( size() );
431  std::transform( _p.begin(), _p.end(), back_inserter( r._p ), op );
432  return r;
433  }
434 
436  this_type operator- () const { return pwUnaryTr( std::negate<T>() ); }
437 
439  this_type abs() const { return pwUnaryTr( fo_abs<T>() ); }
440 
442  this_type exp() const { return pwUnaryTr( fo_exp<T>() ); }
443 
445 
447  this_type log(bool zero=false) const {
448  if( zero )
449  return pwUnaryTr( fo_log0<T>() );
450  else
451  return pwUnaryTr( fo_log<T>() );
452  }
453 
455 
457  this_type inverse(bool zero=true) const {
458  if( zero )
459  return pwUnaryTr( fo_inv0<T>() );
460  else
461  return pwUnaryTr( fo_inv<T>() );
462  }
463 
465 
467  this_type normalized( ProbNormType norm = dai::NORMPROB ) const {
468  T Z = 0;
469  if( norm == dai::NORMPROB )
470  Z = sum();
471  else if( norm == dai::NORMLINF )
472  Z = maxAbs();
473  if( Z == (T)0 ) {
474  DAI_THROW(NOT_NORMALIZABLE);
475  return *this;
476  } else
477  return pwUnaryTr( std::bind2nd( std::divides<T>(), Z ) );
478  }
480 
482 
483 
484  template<typename unaryOp> this_type& pwUnaryOp( unaryOp op ) {
485  std::transform( _p.begin(), _p.end(), _p.begin(), op );
486  return *this;
487  }
488 
491  std::generate( _p.begin(), _p.end(), rnd_uniform );
492  return *this;
493  }
494 
497  fill( (T)1 / size() );
498  return *this;
499  }
500 
502  this_type& takeAbs() { return pwUnaryOp( fo_abs<T>() ); }
503 
505  this_type& takeExp() { return pwUnaryOp( fo_exp<T>() ); }
506 
508 
510  this_type& takeLog(bool zero=false) {
511  if( zero ) {
512  return pwUnaryOp( fo_log0<T>() );
513  } else
514  return pwUnaryOp( fo_log<T>() );
515  }
516 
518 
520  T normalize( ProbNormType norm=dai::NORMPROB ) {
521  T Z = 0;
522  if( norm == dai::NORMPROB )
523  Z = sum();
524  else if( norm == dai::NORMLINF )
525  Z = maxAbs();
526  if( Z == (T)0 )
527  DAI_THROW(NOT_NORMALIZABLE);
528  else
529  *this /= Z;
530  return Z;
531  }
533 
535 
536 
537  this_type& fill( T x ) {
538  std::fill( _p.begin(), _p.end(), x );
539  return *this;
540  }
541 
544  if( x != 0 )
545  return pwUnaryOp( std::bind2nd( std::plus<T>(), x ) );
546  else
547  return *this;
548  }
549 
552  if( x != 0 )
553  return pwUnaryOp( std::bind2nd( std::minus<T>(), x ) );
554  else
555  return *this;
556  }
557 
560  if( x != 1 )
561  return pwUnaryOp( std::bind2nd( std::multiplies<T>(), x ) );
562  else
563  return *this;
564  }
565 
568  if( x != 1 )
569  return pwUnaryOp( std::bind2nd( fo_divides0<T>(), x ) );
570  else
571  return *this;
572  }
573 
576  if( x != (T)1 )
577  return pwUnaryOp( std::bind2nd( fo_pow<T>(), x) );
578  else
579  return *this;
580  }
582 
584 
585 
586  this_type operator+ (T x) const { return pwUnaryTr( std::bind2nd( std::plus<T>(), x ) ); }
587 
589  this_type operator- (T x) const { return pwUnaryTr( std::bind2nd( std::minus<T>(), x ) ); }
590 
592  this_type operator* (T x) const { return pwUnaryTr( std::bind2nd( std::multiplies<T>(), x ) ); }
593 
595  this_type operator/ (T x) const { return pwUnaryTr( std::bind2nd( fo_divides0<T>(), x ) ); }
596 
598  this_type operator^ (T x) const { return pwUnaryTr( std::bind2nd( fo_pow<T>(), x ) ); }
600 
602 
603 
604 
608  template<typename binaryOp> this_type& pwBinaryOp( const this_type &q, binaryOp op ) {
609  DAI_DEBASSERT( size() == q.size() );
610  std::transform( _p.begin(), _p.end(), q._p.begin(), _p.begin(), op );
611  return *this;
612  }
613 
615 
617  this_type& operator+= (const this_type & q) { return pwBinaryOp( q, std::plus<T>() ); }
618 
620 
622  this_type& operator-= (const this_type & q) { return pwBinaryOp( q, std::minus<T>() ); }
623 
625 
627  this_type& operator*= (const this_type & q) { return pwBinaryOp( q, std::multiplies<T>() ); }
628 
630 
633  this_type& operator/= (const this_type & q) { return pwBinaryOp( q, fo_divides0<T>() ); }
634 
636 
639  this_type& divide (const this_type & q) { return pwBinaryOp( q, std::divides<T>() ); }
640 
642 
644  this_type& operator^= (const this_type & q) { return pwBinaryOp( q, fo_pow<T>() ); }
646 
648 
649 
650 
654  template<typename binaryOp> this_type pwBinaryTr( const this_type &q, binaryOp op ) const {
655  DAI_DEBASSERT( size() == q.size() );
656  TProb<T> r;
657  r._p.reserve( size() );
658  std::transform( _p.begin(), _p.end(), q._p.begin(), back_inserter( r._p ), op );
659  return r;
660  }
661 
663 
665  this_type operator+ ( const this_type& q ) const { return pwBinaryTr( q, std::plus<T>() ); }
666 
668 
670  this_type operator- ( const this_type& q ) const { return pwBinaryTr( q, std::minus<T>() ); }
671 
673 
675  this_type operator* ( const this_type &q ) const { return pwBinaryTr( q, std::multiplies<T>() ); }
676 
678 
681  this_type operator/ ( const this_type &q ) const { return pwBinaryTr( q, fo_divides0<T>() ); }
682 
684 
687  this_type divided_by( const this_type &q ) const { return pwBinaryTr( q, std::divides<T>() ); }
688 
690 
692  this_type operator^ ( const this_type &q ) const { return pwBinaryTr( q, fo_pow<T>() ); }
694 
696 
698  template<typename binOp1, typename binOp2> T innerProduct( const this_type &q, T init, binOp1 binaryOp1, binOp2 binaryOp2 ) const {
699  DAI_DEBASSERT( size() == q.size() );
700  return std::inner_product( begin(), end(), q.begin(), init, binaryOp1, binaryOp2 );
701  }
702 };
703 
704 
706 
709 template<typename T> T dist( const TProb<T> &p, const TProb<T> &q, ProbDistType dt ) {
710  switch( dt ) {
711  case DISTL1:
712  return p.innerProduct( q, (T)0, std::plus<T>(), fo_absdiff<T>() );
713  case DISTLINF:
714  return p.innerProduct( q, (T)0, fo_max<T>(), fo_absdiff<T>() );
715  case DISTTV:
716  return p.innerProduct( q, (T)0, std::plus<T>(), fo_absdiff<T>() ) / 2;
717  case DISTKL:
718  return p.innerProduct( q, (T)0, std::plus<T>(), fo_KL<T>() );
719  case DISTHEL:
720  return p.innerProduct( q, (T)0, std::plus<T>(), fo_Hellinger<T>() ) / 2;
721  default:
722  DAI_THROW(UNKNOWN_ENUM_VALUE);
723  return INFINITY;
724  }
725 }
726 
727 
729 
731 template<typename T> std::ostream& operator<< (std::ostream& os, const TProb<T>& p) {
732  os << "(";
733  for( size_t i = 0; i < p.size(); i++ )
734  os << ((i != 0) ? ", " : "") << p.get(i);
735  os << ")";
736  return os;
737 }
738 
739 
741 
744 template<typename T> TProb<T> min( const TProb<T> &a, const TProb<T> &b ) {
745  return a.pwBinaryTr( b, fo_min<T>() );
746 }
747 
748 
750 
753 template<typename T> TProb<T> max( const TProb<T> &a, const TProb<T> &b ) {
754  return a.pwBinaryTr( b, fo_max<T>() );
755 }
756 
757 
760 
761 
762 } // end of namespace dai
763 
764 
765 #endif