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