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
5 This file is part of libDAI.
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.
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.
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
24 /// \brief Defines TProb<T> and Prob classes
25 /// \todo Rename to Vector<T>
28 #ifndef __defined_libdai_prob_h
29 #define __defined_libdai_prob_h
40 #include <dai/exceptions.h>
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.
52 template <typename T
> class TProb
{
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
;
63 /// Enumerates different ways of normalizing a probability measure.
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.
68 typedef enum { NORMPROB
, NORMLINF
} NormType
;
69 /// Enumerates different distance measures between probability measures.
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.
76 typedef enum { DISTL1
, DISTLINF
, DISTTV
, DISTKL
} DistType
;
78 /// Default constructor
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
)) {}
84 /// Construct vector of length n with each entry set to p
85 explicit TProb( size_t n
, Real p
) : _p(n
, (T
)p
) {}
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.
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
);
99 /// Returns a const reference to the vector
100 const std::vector
<T
> & p() const { return _p
; }
102 /// Returns a reference to the vector
103 std::vector
<T
> & p() { return _p
; }
105 /// Returns a copy of the i'th entry
106 T
operator[]( size_t i
) const {
114 /// Returns reference to the i'th entry
115 T
& operator[]( size_t i
) { return _p
[i
]; }
117 /// Returns iterator pointing to first entry
118 iterator
begin() { return _p
.begin(); }
120 /// Returns const iterator pointing to first entry
121 const_iterator
begin() const { return _p
.begin(); }
123 /// Returns iterator pointing beyond last entry
124 iterator
end() { return _p
.end(); }
126 /// Returns const iterator pointing beyond last entry
127 const_iterator
end() const { return _p
.end(); }
129 /// Sets all entries to x
130 TProb
<T
> & fill(T x
) {
131 std::fill( _p
.begin(), _p
.end(), x
);
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
);
141 /// Returns length of the vector, i.e., the number of entries
142 size_t size() const {
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
)
154 /// Sets entries that are smaller than epsilon to epsilon
155 TProb
<T
>& makePositive( Real epsilon
) {
156 for( size_t i
= 0; i
< size(); i
++ )
157 if( (0 < (Real
)_p
[i
]) && ((Real
)_p
[i
] < epsilon
) )
162 /// Multiplies each entry with scalar x
163 TProb
<T
>& operator*= (T x
) {
164 std::transform( _p
.begin(), _p
.end(), _p
.begin(), std::bind2nd( std::multiplies
<T
>(), x
) );
168 /// Returns product of *this with scalar x
169 TProb
<T
> operator* (T x
) const {
170 TProb
<T
> prod( *this );
175 /// Divides each entry by scalar x
176 TProb
<T
>& operator/= (T x
) {
180 std::transform( _p
.begin(), _p
.end(), _p
.begin(), std::bind2nd( std::divides
<T
>(), x
) );
184 /// Returns quotient of *this and scalar x
185 TProb
<T
> operator/ (T x
) const {
186 TProb
<T
> quot( *this );
191 /// Adds scalar x to each entry
192 TProb
<T
>& operator+= (T x
) {
193 std::transform( _p
.begin(), _p
.end(), _p
.begin(), std::bind2nd( std::plus
<T
>(), x
) );
197 /// Returns sum of *this and scalar x
198 TProb
<T
> operator+ (T x
) const {
199 TProb
<T
> sum( *this );
204 /// Subtracts scalar x from each entry
205 TProb
<T
>& operator-= (T x
) {
206 std::transform( _p
.begin(), _p
.end(), _p
.begin(), std::bind2nd( std::minus
<T
>(), x
) );
210 /// Returns difference of *this and scalar x
211 TProb
<T
> operator- (T x
) const {
212 TProb
<T
> diff( *this );
217 /// Lexicographical comparison (sizes should be identical)
218 bool operator<= (const TProb
<T
> & q
) const {
220 assert( size() == q
.size() );
222 for( size_t i
= 0; i
< size(); i
++ )
223 if( !(_p
[i
] <= q
[i
]) )
228 /// Pointwise multiplication with q (sizes should be identical)
229 TProb
<T
>& operator*= (const TProb
<T
> & q
) {
231 assert( size() == q
.size() );
233 std::transform( _p
.begin(), _p
.end(), q
._p
.begin(), _p
.begin(), std::multiplies
<T
>() );
237 /// Return product of *this with q (sizes should be identical)
238 TProb
<T
> operator* (const TProb
<T
> & q
) const {
240 assert( size() == q
.size() );
242 TProb
<T
> prod( *this );
247 /// Pointwise addition with q (sizes should be identical)
248 TProb
<T
>& operator+= (const TProb
<T
> & q
) {
250 assert( size() == q
.size() );
252 std::transform( _p
.begin(), _p
.end(), q
._p
.begin(), _p
.begin(), std::plus
<T
>() );
256 /// Returns sum of *this and q (sizes should be identical)
257 TProb
<T
> operator+ (const TProb
<T
> & q
) const {
259 assert( size() == q
.size() );
261 TProb
<T
> sum( *this );
266 /// Pointwise subtraction of q (sizes should be identical)
267 TProb
<T
>& operator-= (const TProb
<T
> & q
) {
269 assert( size() == q
.size() );
271 std::transform( _p
.begin(), _p
.end(), q
._p
.begin(), _p
.begin(), std::minus
<T
>() );
275 /// Return *this minus q (sizes should be identical)
276 TProb
<T
> operator- (const TProb
<T
> & q
) const {
278 assert( size() == q
.size() );
280 TProb
<T
> diff( *this );
285 /// Pointwise division by q, where division by 0 yields 0 (sizes should be identical)
286 TProb
<T
>& operator/= (const TProb
<T
> & q
) {
288 assert( size() == q
.size() );
290 for( size_t i
= 0; i
< size(); i
++ ) {
299 /// Pointwise division by q, where division by 0 yields +Inf (sizes should be identical)
300 TProb
<T
>& divide (const TProb
<T
> & q
) {
302 assert( size() == q
.size() );
304 std::transform( _p
.begin(), _p
.end(), q
._p
.begin(), _p
.begin(), std::divides
<T
>() );
308 /// Returns quotient of *this with q (sizes should be identical)
309 TProb
<T
> operator/ (const TProb
<T
> & q
) const {
311 assert( size() == q
.size() );
313 TProb
<T
> quot( *this );
318 /// Returns pointwise inverse
319 /** If zero==true; uses 1/0==0, otherwise 1/0==Inf.
321 TProb
<T
> inverse(bool zero
=true) const {
323 inv
._p
.reserve( size() );
325 for( size_t i
= 0; i
< size(); i
++ )
326 inv
._p
.push_back( _p
[i
] == 0.0 ? 0.0 : 1.0 / _p
[i
] );
328 for( size_t i
= 0; i
< size(); i
++ )
329 inv
._p
.push_back( 1.0 / _p
[i
] );
333 /// Raises entries to the power a
334 TProb
<T
>& operator^= (Real a
) {
336 std::transform( _p
.begin(), _p
.end(), _p
.begin(), std::bind2nd( std::ptr_fun
<T
, Real
, T
>(std::pow
), a
) );
340 /// Returns *this raised to the power a
341 TProb
<T
> operator^ (Real a
) const {
342 TProb
<T
> power(*this);
347 /// Returns pointwise signum
348 TProb
<T
> sgn() const {
350 x
._p
.reserve( size() );
351 for( size_t i
= 0; i
< size(); i
++ ) {
362 /// Returns pointwise absolute value
363 TProb
<T
> abs() const {
365 x
._p
.reserve( size() );
366 for( size_t i
= 0; i
< size(); i
++ )
367 x
._p
.push_back( _p
[i
] < 0 ? (-p
[i
]) : p
[i
] );
371 /// Applies exp pointwise
372 const TProb
<T
>& takeExp() {
373 std::transform( _p
.begin(), _p
.end(), _p
.begin(), std::ptr_fun
<T
, T
>(std::exp
) );
377 /// Applies log pointwise
378 /** If zero==true, uses log(0)==0; otherwise, log(0)==-Inf.
380 const TProb
<T
>& takeLog(bool zero
=false) {
382 for( size_t i
= 0; i
< size(); i
++ )
383 _p
[i
] = ( (_p
[i
] == 0.0) ? 0.0 : std::log( _p
[i
] ) );
385 std::transform( _p
.begin(), _p
.end(), _p
.begin(), std::ptr_fun
<T
, T
>(std::log
) );
389 /// Returns pointwise exp
390 TProb
<T
> exp() const {
396 /// Returns pointwise log
397 /** If zero==true, uses log(0)==0; otherwise, log(0)==-Inf.
399 TProb
<T
> log(bool zero
=false) const {
405 /// Returns sum of all entries
407 T Z
= std::accumulate( _p
.begin(), _p
.end(), (T
)0 );
411 /// Returns maximum absolute value of all entries
414 for( size_t i
= 0; i
< size(); i
++ ) {
415 Real mag
= fabs( (Real
) _p
[i
] );
422 /// Returns maximum value of all entries
424 T Z
= *std::max_element( _p
.begin(), _p
.end() );
428 /// Returns minimum value of all entries
430 T Z
= *std::min_element( _p
.begin(), _p
.end() );
434 /// Normalizes vector using the specified norm
435 T
normalize( NormType norm
=NORMPROB
) {
437 if( norm
== NORMPROB
)
439 else if( norm
== NORMLINF
)
442 DAI_THROW(NOT_NORMALIZABLE
);
448 /// Returns normalized copy of *this, using the specified norm
449 TProb
<T
> normalized( NormType norm
= NORMPROB
) const {
450 TProb
<T
> result(*this);
451 result
.normalize( norm
);
455 /// Returns true if one or more entries are NaN
456 bool hasNaNs() const {
457 bool foundnan
= false;
458 for( typename
std::vector
<T
>::const_iterator x
= _p
.begin(); x
!= _p
.end(); x
++ )
466 /// Returns true if one or more entries are negative
467 bool hasNegatives() const {
468 return (std::find_if( _p
.begin(), _p
.end(), std::bind2nd( std::less
<Real
>(), 0.0 ) ) != _p
.end());
471 /// Returns entropy of *this
472 Real
entropy() const {
474 for( size_t i
= 0; i
< size(); i
++ )
475 S
-= (_p
[i
] == 0 ? 0 : _p
[i
] * std::log(_p
[i
]));
479 /// Returns a random index, according to the (normalized) distribution described by *this
481 double x
= rnd_uniform() * totalSum();
483 for( size_t i
= 0; i
< size(); i
++ ) {
488 return( size() - 1 );
493 /// Returns distance of p and q (sizes should be identical), measured using distance measure dt
496 template<typename T
> Real
dist( const TProb
<T
> &p
, const TProb
<T
> &q
, typename TProb
<T
>::DistType dt
) {
498 assert( p
.size() == q
.size() );
502 case TProb
<T
>::DISTL1
:
503 for( size_t i
= 0; i
< p
.size(); i
++ )
504 result
+= fabs((Real
)p
[i
] - (Real
)q
[i
]);
507 case TProb
<T
>::DISTLINF
:
508 for( size_t i
= 0; i
< p
.size(); i
++ ) {
509 Real z
= fabs((Real
)p
[i
] - (Real
)q
[i
]);
515 case TProb
<T
>::DISTTV
:
516 for( size_t i
= 0; i
< p
.size(); i
++ )
517 result
+= fabs((Real
)p
[i
] - (Real
)q
[i
]);
521 case TProb
<T
>::DISTKL
:
522 for( size_t i
= 0; i
< p
.size(); i
++ ) {
524 result
+= p
[i
] * (std::log(p
[i
]) - std::log(q
[i
]));
531 /// Writes a TProb<T> to an output stream
534 template<typename T
> std::ostream
& operator<< (std::ostream
& os
, const TProb
<T
>& P
) {
536 std::copy( P
.p().begin(), P
.p().end(), std::ostream_iterator
<T
>(os
, " ") );
542 /// Returns the TProb<T> containing the pointwise minimum of a and b (which should have equal size)
545 template<typename T
> TProb
<T
> min( const TProb
<T
> &a
, const TProb
<T
> &b
) {
546 assert( a
.size() == b
.size() );
547 TProb
<T
> result( a
.size() );
548 for( size_t i
= 0; i
< a
.size(); i
++ )
557 /// Returns the TProb<T> containing the pointwise maximum of a and b (which should have equal size)
560 template<typename T
> TProb
<T
> max( const TProb
<T
> &a
, const TProb
<T
> &b
) {
561 assert( a
.size() == b
.size() );
562 TProb
<T
> result( a
.size() );
563 for( size_t i
= 0; i
< a
.size(); i
++ )
572 /// Represents a vector with entries of type Real.
573 typedef TProb
<Real
> Prob
;
576 } // end of namespace dai