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
27 #ifndef __defined_libdai_prob_h
28 #define __defined_libdai_prob_h
44 /// Real number (alias for double, could be changed to long double if necessary)
47 template<typename T
> class TProb
;
49 /// Represents a probability measure, with entries of type Real.
50 typedef TProb
<Real
> Prob
;
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.
58 template <typename T
> class TProb
{
60 /// The probability measure
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
;
69 /// Enumerates different ways of normalizing a probability measure.
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.
74 typedef enum { NORMPROB
, NORMLINF
} NormType
;
75 /// Enumerates different distance measures between probability measures.
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.
82 typedef enum { DISTL1
, DISTLINF
, DISTTV
, DISTKL
} DistType
;
84 /// Default constructor
87 /// Construct uniform distribution of given length
88 explicit TProb( size_t n
) : _p(std::vector
<T
>(n
, 1.0 / n
)) {}
90 /// Construct from given length and initial value
91 TProb( size_t n
, Real p
) : _p(n
, (T
)p
) {}
93 /// Construct from given length and initial array
94 TProb( size_t n
, const Real
* p
) : _p(p
, p
+ n
) {}
96 /// Returns a const reference to the probability vector
97 const std::vector
<T
> & p() const { return _p
; }
99 /// Returns a reference to the probability vector
100 std::vector
<T
> & p() { return _p
; }
102 /// Returns a copy of the i'th probability entry
103 T
operator[]( size_t i
) const {
111 /// Returns a reference to the i'th probability entry
112 T
& operator[]( size_t i
) { return _p
[i
]; }
114 /// Returns iterator pointing to first entry
115 iterator
begin() { return _p
.begin(); }
117 /// Returns const iterator pointing to first entry
118 const_iterator
begin() const { return _p
.begin(); }
120 /// Returns iterator pointing beyond last entry
121 iterator
end() { return _p
.end(); }
123 /// Returns const iterator pointing beyond last entry
124 const_iterator
end() const { return _p
.end(); }
126 /// Sets all elements to x
127 TProb
<T
> & fill(T x
) {
128 std::fill( _p
.begin(), _p
.end(), x
);
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
);
138 /// Returns number of elements
139 size_t size() const {
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
)
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
) )
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
) );
165 /// Returns product of *this with x
166 TProb
<T
> operator* (T x
) const {
167 TProb
<T
> prod( *this );
172 /// Divides each entry by x
173 TProb
<T
>& operator/= (T x
) {
177 std::transform( _p
.begin(), _p
.end(), _p
.begin(), std::bind2nd( std::divides
<T
>(), x
) );
181 /// Returns quotient of *this and x
182 TProb
<T
> operator/ (T x
) const {
183 TProb
<T
> quot( *this );
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
) );
194 /// Returns sum of *this and x
195 TProb
<T
> operator+ (T x
) const {
196 TProb
<T
> sum( *this );
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
) );
207 /// Returns difference of *this and x
208 TProb
<T
> operator- (T x
) const {
209 TProb
<T
> diff( *this );
214 /// Pointwise comparison
215 bool operator<= (const TProb
<T
> & q
) const {
217 assert( size() == q
.size() );
219 for( size_t i
= 0; i
< size(); i
++ )
220 if( !(_p
[i
] <= q
[i
]) )
225 /// Pointwise multiplication with q
226 TProb
<T
>& operator*= (const TProb
<T
> & q
) {
228 assert( size() == q
.size() );
230 std::transform( _p
.begin(), _p
.end(), q
._p
.begin(), _p
.begin(), std::multiplies
<T
>() );
234 /// Return product of *this with q
235 TProb
<T
> operator* (const TProb
<T
> & q
) const {
237 assert( size() == q
.size() );
239 TProb
<T
> prod( *this );
244 /// Pointwise addition with q
245 TProb
<T
>& operator+= (const TProb
<T
> & q
) {
247 assert( size() == q
.size() );
249 std::transform( _p
.begin(), _p
.end(), q
._p
.begin(), _p
.begin(), std::plus
<T
>() );
253 /// Return sum of *this and q
254 TProb
<T
> operator+ (const TProb
<T
> & q
) const {
256 assert( size() == q
.size() );
258 TProb
<T
> sum( *this );
263 /// Pointwise subtraction of q
264 TProb
<T
>& operator-= (const TProb
<T
> & q
) {
266 assert( size() == q
.size() );
268 std::transform( _p
.begin(), _p
.end(), q
._p
.begin(), _p
.begin(), std::minus
<T
>() );
272 /// Return *this minus q
273 TProb
<T
> operator- (const TProb
<T
> & q
) const {
275 assert( size() == q
.size() );
277 TProb
<T
> diff( *this );
282 /// Pointwise division by q, where division by zero yields zero
283 TProb
<T
>& operator/= (const TProb
<T
> & q
) {
285 assert( size() == q
.size() );
287 for( size_t i
= 0; i
< size(); i
++ ) {
296 /// Pointwise division by q, where division by zero yields infinity
297 TProb
<T
>& divide (const TProb
<T
> & q
) {
299 assert( size() == q
.size() );
301 std::transform( _p
.begin(), _p
.end(), q
._p
.begin(), _p
.begin(), std::divides
<T
>() );
305 /// Returns quotient of *this with q
306 TProb
<T
> operator/ (const TProb
<T
> & q
) const {
308 assert( size() == q
.size() );
310 TProb
<T
> quot( *this );
315 /// Returns pointwise inverse
316 TProb
<T
> inverse(bool zero
= false) const {
318 inv
._p
.reserve( size() );
320 for( size_t i
= 0; i
< size(); i
++ )
321 inv
._p
.push_back( _p
[i
] == 0.0 ? 0.0 : 1.0 / _p
[i
] );
323 for( size_t i
= 0; i
< size(); i
++ ) {
325 assert( _p
[i
] != 0.0 );
327 inv
._p
.push_back( 1.0 / _p
[i
] );
332 /// Raises elements to the power a
333 TProb
<T
>& operator^= (Real a
) {
335 std::transform( _p
.begin(), _p
.end(), _p
.begin(), std::bind2nd( std::ptr_fun
<T
, Real
, T
>(std::pow
), a
) );
339 /// Returns *this raised to the power a
340 TProb
<T
> operator^ (Real a
) const {
341 TProb
<T
> power(*this);
346 /// Returns pointwise signum
347 TProb
<T
> sgn() const {
349 x
._p
.reserve( size() );
350 for( size_t i
= 0; i
< size(); i
++ ) {
361 /// Returns pointwise absolute value
362 TProb
<T
> abs() const {
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
] );
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
) );
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
) );
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
] ) );
389 /// Returns pointwise exp
390 TProb
<T
> exp() const {
396 /// Returns pointwise log
397 TProb
<T
> log() const {
403 /// Returns pointwise log (defining log(0)=0)
404 TProb
<T
> log0() const {
410 /// Returns distance of p and q, measured using dt
411 friend Real
dist( const TProb
<T
> &p
, const TProb
<T
> &q
, DistType dt
) {
413 assert( p
.size() == q
.size() );
418 for( size_t i
= 0; i
< p
.size(); i
++ )
419 result
+= fabs((Real
)p
[i
] - (Real
)q
[i
]);
423 for( size_t i
= 0; i
< p
.size(); i
++ ) {
424 Real z
= fabs((Real
)p
[i
] - (Real
)q
[i
]);
431 for( size_t i
= 0; i
< p
.size(); i
++ )
432 result
+= fabs((Real
)p
[i
] - (Real
)q
[i
]);
437 for( size_t i
= 0; i
< p
.size(); i
++ ) {
439 result
+= p
[i
] * (std::log(p
[i
]) - std::log(q
[i
]));
445 /// Returns sum of all entries
447 T Z
= std::accumulate( _p
.begin(), _p
.end(), (T
)0 );
451 /// Returns maximum absolute value of entries
454 for( size_t i
= 0; i
< size(); i
++ ) {
455 Real mag
= fabs( (Real
) _p
[i
] );
462 /// Returns maximum value of entries
464 T Z
= *std::max_element( _p
.begin(), _p
.end() );
468 /// Returns minimum value of entries
470 T Z
= *std::min_element( _p
.begin(), _p
.end() );
474 /// Normalizes using the specified norm
475 T
normalize( NormType norm
= NORMPROB
) {
477 if( norm
== NORMPROB
)
479 else if( norm
== NORMLINF
)
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
);
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());
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());
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());
511 Real
entropy() const {
513 for( size_t i
= 0; i
< size(); i
++ )
518 /// Writes a TProb<T> to an output stream
519 friend std::ostream
& operator<< (std::ostream
& os
, const TProb
<T
>& P
) {
521 std::copy( P
._p
.begin(), P
._p
.end(), std::ostream_iterator
<T
>(os
, " ") );
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
)); }
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
++ )
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
++ )
558 } // end of namespace dai