1 /* This file is part of libDAI - http://www.libdai.org/
3 * libDAI is licensed under the terms of the GNU General Public License version
4 * 2, or (at your option) any later version. libDAI is distributed without any
5 * warranty. See the file COPYING for more details.
7 * Copyright (C) 2006-2009 Joris Mooij [joris dot mooij at libdai dot org]
8 * Copyright (C) 2006-2007 Radboud University Nijmegen, The Netherlands
13 /// \brief Defines TProb<T> and Prob classes
14 /// \todo Rename to Vector<T>
17 #ifndef __defined_libdai_prob_h
18 #define __defined_libdai_prob_h
29 #include <dai/exceptions.h>
35 /// Represents a vector with entries of type \a T.
36 /** A TProb<T> is a std::vector<T> with an interface designed for dealing with probability mass functions.
37 * It is mainly used for representing measures on a finite outcome space, e.g., the probability
38 * distribution of a discrete random variable.
39 * \tparam T Should be a scalar that is castable from and to double and should support elementary arithmetic operations.
41 template <typename T
> class TProb
{
47 /// Iterator over entries
48 typedef typename
std::vector
<T
>::iterator iterator
;
49 /// Const iterator over entries
50 typedef typename
std::vector
<T
>::const_iterator const_iterator
;
52 /// Enumerates different ways of normalizing a probability measure.
54 * - NORMPROB means that the sum of all entries should be 1;
55 * - NORMLINF means that the maximum absolute value of all entries should be 1.
57 typedef enum { NORMPROB
, NORMLINF
} NormType
;
58 /// Enumerates different distance measures between probability measures.
60 * - DISTL1 is the L-1 distance (sum of absolute values of pointwise difference);
61 * - DISTLINF is the L-inf distance (maximum absolute value of pointwise difference);
62 * - DISTTV is the Total Variation distance;
63 * - DISTKL is the Kullback-Leibler distance.
65 typedef enum { DISTL1
, DISTLINF
, DISTTV
, DISTKL
} DistType
;
67 /// Default constructor
70 /// Construct uniform distribution over n outcomes, i.e., a vector of length n with each entry set to 1/n
71 explicit TProb( size_t n
) : _p(std::vector
<T
>(n
, 1.0 / n
)) {}
73 /// Construct vector of length n with each entry set to p
74 explicit TProb( size_t n
, Real p
) : _p(n
, (T
)p
) {}
76 /// Construct vector from a range
77 /** \tparam Iterator Iterates over instances that can be cast to T.
78 * \param begin Points to first instance to be added.
79 * \param end Points just beyond last instance to be added.
80 * \param sizeHint For efficiency, the number of entries can be speficied by sizeHint.
82 template <typename Iterator
>
83 TProb( Iterator begin
, Iterator end
, size_t sizeHint
=0 ) : _p() {
84 _p
.reserve( sizeHint
);
85 _p
.insert( _p
.begin(), begin
, end
);
88 /// Construct vector from a vector
89 /** \tparam S type of elements in v
90 * \param v vector used for initialization
93 TProb( const std::vector
<S
> &v
) : _p() {
94 _p
.reserve( v
.size() );
95 _p
.insert( _p
.begin(), v
.begin(), v
.end() );
98 /// Returns a const reference to the vector
99 const std::vector
<T
> & p() const { return _p
; }
101 /// Returns a reference to the vector
102 std::vector
<T
> & p() { return _p
; }
104 /// Returns a copy of the i'th entry
105 T
operator[]( size_t i
) const {
113 /// Returns reference to the i'th entry
114 T
& operator[]( size_t i
) { return _p
[i
]; }
116 /// Returns iterator pointing to first entry
117 iterator
begin() { return _p
.begin(); }
119 /// Returns const iterator pointing to first entry
120 const_iterator
begin() const { return _p
.begin(); }
122 /// Returns iterator pointing beyond last entry
123 iterator
end() { return _p
.end(); }
125 /// Returns const iterator pointing beyond last entry
126 const_iterator
end() const { return _p
.end(); }
128 /// Sets all entries to x
129 TProb
<T
> & fill(T x
) {
130 std::fill( _p
.begin(), _p
.end(), x
);
134 /// Draws all entries i.i.d. from a uniform distribution on [0,1)
135 TProb
<T
> & randomize() {
136 std::generate(_p
.begin(), _p
.end(), rnd_uniform
);
140 /// Returns length of the vector, i.e., the number of entries
141 size_t size() const {
145 /// Sets entries that are smaller than epsilon to 0
146 TProb
<T
>& makeZero( Real epsilon
) {
147 for( size_t i
= 0; i
< size(); i
++ )
148 if( fabs(_p
[i
]) < epsilon
)
153 /// Set all entries to 1.0/size()
154 TProb
<T
>& setUniform () {
159 /// Sets entries that are smaller than epsilon to epsilon
160 TProb
<T
>& makePositive( Real epsilon
) {
161 for( size_t i
= 0; i
< size(); i
++ )
162 if( (0 < (Real
)_p
[i
]) && ((Real
)_p
[i
] < epsilon
) )
167 /// Multiplies each entry with scalar x
168 TProb
<T
>& operator*= (T x
) {
169 std::transform( _p
.begin(), _p
.end(), _p
.begin(), std::bind2nd( std::multiplies
<T
>(), x
) );
173 /// Returns product of *this with scalar x
174 TProb
<T
> operator* (T x
) const {
175 TProb
<T
> prod( *this );
180 /// Divides each entry by scalar x
181 TProb
<T
>& operator/= (T x
) {
182 DAI_DEBASSERT( x
!= 0.0 );
183 std::transform( _p
.begin(), _p
.end(), _p
.begin(), std::bind2nd( std::divides
<T
>(), x
) );
187 /// Returns quotient of *this and scalar x
188 TProb
<T
> operator/ (T x
) const {
189 TProb
<T
> quot( *this );
194 /// Adds scalar x to each entry
195 TProb
<T
>& operator+= (T x
) {
196 std::transform( _p
.begin(), _p
.end(), _p
.begin(), std::bind2nd( std::plus
<T
>(), x
) );
200 /// Returns sum of *this and scalar x
201 TProb
<T
> operator+ (T x
) const {
202 TProb
<T
> sum( *this );
207 /// Subtracts scalar x from each entry
208 TProb
<T
>& operator-= (T x
) {
209 std::transform( _p
.begin(), _p
.end(), _p
.begin(), std::bind2nd( std::minus
<T
>(), x
) );
213 /// Returns difference of *this and scalar x
214 TProb
<T
> operator- (T x
) const {
215 TProb
<T
> diff( *this );
220 /// Lexicographical comparison (sizes should be identical)
221 bool operator<= (const TProb
<T
> & q
) const {
222 DAI_DEBASSERT( size() == q
.size() );
223 for( size_t i
= 0; i
< size(); i
++ )
224 if( !(_p
[i
] <= q
[i
]) )
229 /// Pointwise multiplication with q (sizes should be identical)
230 TProb
<T
>& operator*= (const TProb
<T
> & q
) {
231 DAI_DEBASSERT( size() == q
.size() );
232 std::transform( _p
.begin(), _p
.end(), q
._p
.begin(), _p
.begin(), std::multiplies
<T
>() );
236 /// Return product of *this with q (sizes should be identical)
237 TProb
<T
> operator* (const TProb
<T
> & q
) const {
238 DAI_DEBASSERT( size() == q
.size() );
239 TProb
<T
> prod( *this );
244 /// Pointwise addition with q (sizes should be identical)
245 TProb
<T
>& operator+= (const TProb
<T
> & q
) {
246 DAI_DEBASSERT( size() == q
.size() );
247 std::transform( _p
.begin(), _p
.end(), q
._p
.begin(), _p
.begin(), std::plus
<T
>() );
251 /// Returns sum of *this and q (sizes should be identical)
252 TProb
<T
> operator+ (const TProb
<T
> & q
) const {
253 DAI_DEBASSERT( size() == q
.size() );
254 TProb
<T
> sum( *this );
259 /// Pointwise subtraction of q (sizes should be identical)
260 TProb
<T
>& operator-= (const TProb
<T
> & q
) {
261 DAI_DEBASSERT( size() == q
.size() );
262 std::transform( _p
.begin(), _p
.end(), q
._p
.begin(), _p
.begin(), std::minus
<T
>() );
266 /// Return *this minus q (sizes should be identical)
267 TProb
<T
> operator- (const TProb
<T
> & q
) const {
268 DAI_DEBASSERT( size() == q
.size() );
269 TProb
<T
> diff( *this );
274 /// Pointwise division by q, where division by 0 yields 0 (sizes should be identical)
275 TProb
<T
>& operator/= (const TProb
<T
> & q
) {
276 DAI_DEBASSERT( size() == q
.size() );
277 for( size_t i
= 0; i
< size(); i
++ ) {
286 /// Pointwise division by q, where division by 0 yields +Inf (sizes should be identical)
287 TProb
<T
>& divide (const TProb
<T
> & q
) {
288 DAI_DEBASSERT( size() == q
.size() );
289 std::transform( _p
.begin(), _p
.end(), q
._p
.begin(), _p
.begin(), std::divides
<T
>() );
293 /// Returns quotient of *this with q (sizes should be identical)
294 TProb
<T
> operator/ (const TProb
<T
> & q
) const {
295 DAI_DEBASSERT( size() == q
.size() );
296 TProb
<T
> quot( *this );
301 /// Returns pointwise inverse
302 /** If zero==true; uses 1/0==0, otherwise 1/0==Inf.
304 TProb
<T
> inverse(bool zero
=true) const {
306 inv
._p
.reserve( size() );
308 for( size_t i
= 0; i
< size(); i
++ )
309 inv
._p
.push_back( _p
[i
] == 0.0 ? 0.0 : 1.0 / _p
[i
] );
311 for( size_t i
= 0; i
< size(); i
++ )
312 inv
._p
.push_back( 1.0 / _p
[i
] );
316 /// Raises entries to the power a
317 TProb
<T
>& operator^= (Real a
) {
319 std::transform( _p
.begin(), _p
.end(), _p
.begin(), std::bind2nd( std::ptr_fun
<T
, Real
, T
>(std::pow
), a
) );
323 /// Returns *this raised to the power a
324 TProb
<T
> operator^ (Real a
) const {
325 TProb
<T
> power(*this);
330 /// Returns pointwise signum
331 TProb
<T
> sgn() const {
333 x
._p
.reserve( size() );
334 for( size_t i
= 0; i
< size(); i
++ ) {
345 /// Returns pointwise absolute value
346 TProb
<T
> abs() const {
348 x
._p
.reserve( size() );
349 for( size_t i
= 0; i
< size(); i
++ )
350 x
._p
.push_back( _p
[i
] < 0 ? (-_p
[i
]) : _p
[i
] );
354 /// Applies exp pointwise
355 const TProb
<T
>& takeExp() {
356 std::transform( _p
.begin(), _p
.end(), _p
.begin(), std::ptr_fun
<T
, T
>(std::exp
) );
360 /// Applies log pointwise
361 /** If zero==true, uses log(0)==0; otherwise, log(0)==-Inf.
363 const TProb
<T
>& takeLog(bool zero
=false) {
365 for( size_t i
= 0; i
< size(); i
++ )
366 _p
[i
] = ( (_p
[i
] == 0.0) ? 0.0 : std::log( _p
[i
] ) );
368 std::transform( _p
.begin(), _p
.end(), _p
.begin(), std::ptr_fun
<T
, T
>(std::log
) );
372 /// Returns pointwise exp
373 TProb
<T
> exp() const {
379 /// Returns pointwise log
380 /** If zero==true, uses log(0)==0; otherwise, log(0)==-Inf.
382 TProb
<T
> log(bool zero
=false) const {
388 /// Returns sum of all entries
390 T Z
= std::accumulate( _p
.begin(), _p
.end(), (T
)0 );
394 /// Return sum of absolute value of all entries
397 for( size_t i
= 0; i
< size(); i
++ )
398 s
+= fabs( (Real
) _p
[i
] );
402 /// Returns maximum absolute value of all entries
405 for( size_t i
= 0; i
< size(); i
++ ) {
406 Real mag
= fabs( (Real
) _p
[i
] );
413 /// Returns maximum value of all entries
415 T Z
= *std::max_element( _p
.begin(), _p
.end() );
419 /// Returns minimum value of all entries
421 T Z
= *std::min_element( _p
.begin(), _p
.end() );
425 /// Returns {arg,}maximum value
426 std::pair
<size_t,T
> argmax() const {
429 for( size_t i
= 1; i
< size(); i
++ ) {
435 return std::make_pair(arg
,max
);
438 /// Normalizes vector using the specified norm
439 T
normalize( NormType norm
=NORMPROB
) {
441 if( norm
== NORMPROB
)
443 else if( norm
== NORMLINF
)
446 DAI_THROW(NOT_NORMALIZABLE
);
452 /// Returns normalized copy of *this, using the specified norm
453 TProb
<T
> normalized( NormType norm
= NORMPROB
) const {
454 TProb
<T
> result(*this);
455 result
.normalize( norm
);
459 /// Returns true if one or more entries are NaN
460 bool hasNaNs() const {
461 bool foundnan
= false;
462 for( typename
std::vector
<T
>::const_iterator x
= _p
.begin(); x
!= _p
.end(); x
++ )
470 /// Returns true if one or more entries are negative
471 bool hasNegatives() const {
472 return (std::find_if( _p
.begin(), _p
.end(), std::bind2nd( std::less
<Real
>(), 0.0 ) ) != _p
.end());
475 /// Returns entropy of *this
476 Real
entropy() const {
478 for( size_t i
= 0; i
< size(); i
++ )
479 S
-= (_p
[i
] == 0 ? 0 : _p
[i
] * std::log(_p
[i
]));
483 /// Returns a random index, according to the (normalized) distribution described by *this
485 double x
= rnd_uniform() * sum();
487 for( size_t i
= 0; i
< size(); i
++ ) {
492 return( size() - 1 );
497 /// Returns distance of p and q (sizes should be identical), measured using distance measure dt
500 template<typename T
> Real
dist( const TProb
<T
> &p
, const TProb
<T
> &q
, typename TProb
<T
>::DistType dt
) {
501 DAI_DEBASSERT( p
.size() == q
.size() );
504 case TProb
<T
>::DISTL1
:
505 for( size_t i
= 0; i
< p
.size(); i
++ )
506 result
+= fabs((Real
)p
[i
] - (Real
)q
[i
]);
509 case TProb
<T
>::DISTLINF
:
510 for( size_t i
= 0; i
< p
.size(); i
++ ) {
511 Real z
= fabs((Real
)p
[i
] - (Real
)q
[i
]);
517 case TProb
<T
>::DISTTV
:
518 for( size_t i
= 0; i
< p
.size(); i
++ )
519 result
+= fabs((Real
)p
[i
] - (Real
)q
[i
]);
523 case TProb
<T
>::DISTKL
:
524 for( size_t i
= 0; i
< p
.size(); i
++ ) {
526 result
+= p
[i
] * (std::log(p
[i
]) - std::log(q
[i
]));
533 /// Writes a TProb<T> to an output stream
536 template<typename T
> std::ostream
& operator<< (std::ostream
& os
, const TProb
<T
>& P
) {
538 std::copy( P
.p().begin(), P
.p().end(), std::ostream_iterator
<T
>(os
, " ") );
544 /// Returns the TProb<T> containing the pointwise minimum of a and b (which should have equal size)
547 template<typename T
> TProb
<T
> min( const TProb
<T
> &a
, const TProb
<T
> &b
) {
548 assert( a
.size() == b
.size() );
549 TProb
<T
> result( a
.size() );
550 for( size_t i
= 0; i
< a
.size(); i
++ )
559 /// Returns the TProb<T> containing the pointwise maximum of a and b (which should have equal size)
562 template<typename T
> TProb
<T
> max( const TProb
<T
> &a
, const TProb
<T
> &b
) {
563 assert( a
.size() == b
.size() );
564 TProb
<T
> result( a
.size() );
565 for( size_t i
= 0; i
< a
.size(); i
++ )
574 /// Represents a vector with entries of type Real.
575 typedef TProb
<Real
> Prob
;
578 } // end of namespace dai