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