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
23 #ifndef __defined_libdai_prob_h
24 #define __defined_libdai_prob_h
42 template<typename T
> class TProb
;
43 typedef TProb
<Real
> Prob
;
47 template<typename T
> TProb
<T
> min( const TProb
<T
> &a
, const TProb
<T
> &b
);
48 template<typename T
> TProb
<T
> max( const TProb
<T
> &a
, const TProb
<T
> &b
);
51 /// TProb<T> implements a probability vector of type T.
52 /// T should be castable from and to double.
53 template <typename T
> class TProb
{
59 /// Calculate x times log(x), or 0 if x == 0
60 Real
xlogx( Real x
) const { return( x
== 0.0 ? 0.0 : x
* std::log(x
)); }
63 /// NORMPROB means that the sum of all entries should be 1
64 /// NORMLINF means that the maximum absolute value of all entries should be 1
65 typedef enum { NORMPROB
, NORMLINF
} NormType
;
66 /// DISTL1 is the L-1 distance (sum of absolute values of pointwise difference)
67 /// DISTLINF is the L-inf distance (maximum absolute value of pointwise difference)
68 /// DISTTV is the Total Variation distance
69 typedef enum { DISTL1
, DISTLINF
, DISTTV
} DistType
;
71 /// Default constructor
74 /// Construct uniform distribution of given length
75 explicit TProb( size_t n
) : _p(std::vector
<T
>(n
, 1.0 / n
)) {}
77 /// Construct with given length and initial value
78 TProb( size_t n
, Real p
) : _p(n
, (T
)p
) {}
80 /// Construct with given length and initial array
81 TProb( size_t n
, const Real
* p
) : _p(p
, p
+ n
) {}
83 /// Provide read access to _p
84 const std::vector
<T
> & p() const { return _p
; }
86 /// Provide full access to _p
87 std::vector
<T
> & p() { return _p
; }
89 /// Provide read access to ith element of _p
90 T
operator[]( size_t i
) const {
98 /// Provide full access to ith element of _p
99 T
& operator[]( size_t i
) { return _p
[i
]; }
101 /// Set all elements to x
102 TProb
<T
> & fill(T x
) {
103 std::fill( _p
.begin(), _p
.end(), x
);
107 /// Set all elements to iid random numbers from uniform(0,1) distribution
108 TProb
<T
> & randomize() {
109 std::generate(_p
.begin(), _p
.end(), rnd_uniform
);
114 size_t size() const {
118 /// Make entries zero if (Real) absolute value smaller than epsilon
119 TProb
<T
>& makeZero (Real epsilon
) {
120 for( size_t i
= 0; i
< size(); i
++ )
121 if( fabs((Real
)_p
[i
]) < epsilon
)
123 // std::replace_if( _p.begin(), _p.end(), fabs((Real)boost::lambda::_1) < epsilon, 0.0 );
127 /// Make entries epsilon if they are smaller than 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
) )
135 /// Multiplication with T x
136 TProb
<T
>& operator*= (T x
) {
137 std::transform( _p
.begin(), _p
.end(), _p
.begin(), std::bind2nd( std::multiplies
<T
>(), x
) );
141 /// Return product of *this with T x
142 TProb
<T
> operator* (T x
) const {
143 TProb
<T
> prod( *this );
149 TProb
<T
>& operator/= (T x
) {
153 std::transform( _p
.begin(), _p
.end(), _p
.begin(), std::bind2nd( std::divides
<T
>(), x
) );
157 /// Return quotient of *this and T x
158 TProb
<T
> operator/ (T x
) const {
159 TProb
<T
> quot( *this );
165 TProb
<T
>& operator+= (T x
) {
166 std::transform( _p
.begin(), _p
.end(), _p
.begin(), std::bind2nd( std::plus
<T
>(), x
) );
170 /// Return sum of *this with T x
171 TProb
<T
> operator+ (T x
) const {
172 TProb
<T
> sum( *this );
177 /// Difference by T x
178 TProb
<T
>& operator-= (T x
) {
179 std::transform( _p
.begin(), _p
.end(), _p
.begin(), std::bind2nd( std::minus
<T
>(), x
) );
183 /// Return difference of *this and T x
184 TProb
<T
> operator- (T x
) const {
185 TProb
<T
> diff( *this );
190 /// Pointwise comparison
191 bool operator<= (const TProb
<T
> & q
) const {
193 assert( size() == q
.size() );
195 for( size_t i
= 0; i
< size(); i
++ )
196 if( !(_p
[i
] <= q
[i
]) )
201 /// Pointwise multiplication with q
202 TProb
<T
>& operator*= (const TProb
<T
> & q
) {
204 assert( size() == q
.size() );
206 std::transform( _p
.begin(), _p
.end(), q
._p
.begin(), _p
.begin(), std::multiplies
<T
>() );
210 /// Return product of *this with q
211 TProb
<T
> operator* (const TProb
<T
> & q
) const {
213 assert( size() == q
.size() );
215 TProb
<T
> prod( *this );
220 /// Pointwise addition with q
221 TProb
<T
>& operator+= (const TProb
<T
> & q
) {
223 assert( size() == q
.size() );
225 std::transform( _p
.begin(), _p
.end(), q
._p
.begin(), _p
.begin(), std::plus
<T
>() );
229 /// Pointwise subtraction of q
230 TProb
<T
>& operator-= (const TProb
<T
> & q
) {
232 assert( size() == q
.size() );
234 std::transform( _p
.begin(), _p
.end(), q
._p
.begin(), _p
.begin(), std::minus
<T
>() );
238 /// Return sum of *this and q
239 TProb
<T
> operator+ (const TProb
<T
> & q
) const {
241 assert( size() == q
.size() );
243 TProb
<T
> sum( *this );
248 /// Return *this minus q
249 TProb
<T
> operator- (const TProb
<T
> & q
) const {
251 assert( size() == q
.size() );
253 TProb
<T
> diff( *this );
258 /// Pointwise division by q (division by zero yields zero)
259 TProb
<T
>& operator/= (const TProb
<T
> & q
) {
261 assert( size() == q
.size() );
263 for( size_t i
= 0; i
< size(); i
++ ) {
272 /// Pointwise division by q (division by zero yields infinity)
273 TProb
<T
>& divide (const TProb
<T
> & q
) {
275 assert( size() == q
.size() );
277 std::transform( _p
.begin(), _p
.end(), q
._p
.begin(), _p
.begin(), std::divides
<T
>() );
281 /// Return quotient of *this with q
282 TProb
<T
> operator/ (const TProb
<T
> & q
) const {
284 assert( size() == q
.size() );
286 TProb
<T
> quot( *this );
291 /// Return pointwise inverse
292 TProb
<T
> inverse(bool zero
= false) const {
294 inv
._p
.reserve( size() );
296 for( size_t i
= 0; i
< size(); i
++ )
297 inv
._p
.push_back( _p
[i
] == 0.0 ? 0.0 : 1.0 / _p
[i
] );
299 for( size_t i
= 0; i
< size(); i
++ ) {
301 assert( _p
[i
] != 0.0 );
303 inv
._p
.push_back( 1.0 / _p
[i
] );
308 /// Return *this to the power of a (pointwise)
309 TProb
<T
>& operator^= (Real a
) {
311 std::transform( _p
.begin(), _p
.end(), _p
.begin(), std::bind2nd( std::ptr_fun
<T
, Real
, T
>(std::pow
), a
) );
315 /// Pointwise power of a
316 TProb
<T
> operator^ (Real a
) const {
317 TProb
<T
> power(*this);
323 TProb
<T
> sgn() const {
325 x
._p
.reserve( size() );
326 for( size_t i
= 0; i
< size(); i
++ ) {
337 /// Pointwise absolute value
338 TProb
<T
> abs() const {
340 x
._p
.reserve( size() );
341 for( size_t i
= 0; i
< size(); i
++ )
342 x
._p
.push_back( _p
[i
] < 0 ? (-p
[i
]) : p
[i
] );
347 const TProb
<T
>& takeExp() {
348 std::transform( _p
.begin(), _p
.end(), _p
.begin(), std::ptr_fun
<T
, T
>(std::exp
) );
353 const TProb
<T
>& takeLog() {
354 std::transform( _p
.begin(), _p
.end(), _p
.begin(), std::ptr_fun
<T
, T
>(std::log
) );
358 /// Pointwise log (or 0 if == 0)
359 const TProb
<T
>& takeLog0() {
360 for( size_t i
= 0; i
< size(); i
++ )
361 _p
[i
] = ( (_p
[i
] == 0.0) ? 0.0 : std::log( _p
[i
] ) );
366 TProb
<T
> exp() const {
373 TProb
<T
> log() const {
379 /// Pointwise log (or 0 if == 0)
380 TProb
<T
> log0() const {
386 /// Return distance of p and q
387 friend Real
dist( const TProb
<T
> & p
, const TProb
<T
> & q
, DistType dt
) {
389 assert( p
.size() == q
.size() );
394 for( size_t i
= 0; i
< p
.size(); i
++ )
395 result
+= fabs((Real
)p
[i
] - (Real
)q
[i
]);
399 for( size_t i
= 0; i
< p
.size(); i
++ ) {
400 Real z
= fabs((Real
)p
[i
] - (Real
)q
[i
]);
407 for( size_t i
= 0; i
< p
.size(); i
++ )
408 result
+= fabs((Real
)p
[i
] - (Real
)q
[i
]);
415 /// Return Kullback-Leibler distance with q
416 friend Real
KL_dist( const TProb
<T
> & p
, const TProb
<T
> & q
) {
418 assert( p
.size() == q
.size() );
421 for( size_t i
= 0; i
< p
.size(); i
++ ) {
422 if( (Real
) p
[i
] != 0.0 ) {
425 result
+= p_i
* (std::log(p_i
) - std::log(q_i
));
431 /// Return sum of all entries
433 T Z
= std::accumulate( _p
.begin(), _p
.end(), (T
)0 );
437 /// Converts entries to Real and returns maximum absolute value
440 for( size_t i
= 0; i
< size(); i
++ ) {
441 Real mag
= fabs( (Real
) _p
[i
] );
448 /// Returns maximum value
450 T Z
= *std::max_element( _p
.begin(), _p
.end() );
454 /// Returns minimum value
456 T Z
= *std::min_element( _p
.begin(), _p
.end() );
460 /// Normalize, using the specified norm
461 T
normalize( NormType norm
= NORMPROB
) {
463 if( norm
== NORMPROB
)
465 else if( norm
== NORMLINF
)
474 /// Return normalized copy of *this, using the specified norm
475 TProb
<T
> normalized( NormType norm
= NORMPROB
) const {
476 TProb
<T
> result(*this);
477 result
.normalize( norm
);
481 /// Returns true if one or more entries are NaN
482 bool hasNaNs() const {
483 return (std::find_if( _p
.begin(), _p
.end(), isnan
) != _p
.end());
486 /// Returns true if one or more entries are negative
487 bool hasNegatives() const {
488 return (std::find_if( _p
.begin(), _p
.end(), std::bind2nd( std::less
<Real
>(), 0.0 ) ) != _p
.end());
491 /// Returns true if one or more entries are non-positive (causes problems with logscale)
492 bool hasNonPositives() const {
493 return (std::find_if( _p
.begin(), _p
.end(), std::bind2nd( std::less_equal
<Real
>(), 0.0 ) ) != _p
.end());
497 Real
entropy() const {
499 for( size_t i
= 0; i
< size(); i
++ )
504 /// Returns TProb<T> containing the pointwise minimum of a and b (which should have equal size)
505 friend TProb
<T
> min
<> ( const TProb
<T
> &a
, const TProb
<T
> &b
);
507 /// Returns TProb<T> containing the pointwise maximum of a and b (which should have equal size)
508 friend TProb
<T
> max
<> ( const TProb
<T
> &a
, const TProb
<T
> &b
);
510 friend std::ostream
& operator<< (std::ostream
& os
, const TProb
<T
>& P
) {
512 std::copy( P
._p
.begin(), P
._p
.end(), std::ostream_iterator
<T
>(os
, " ") );
519 template<typename T
> TProb
<T
> min( const TProb
<T
> &a
, const TProb
<T
> &b
) {
520 assert( a
.size() == b
.size() );
521 TProb
<T
> result( a
.size() );
522 for( size_t i
= 0; i
< a
.size(); i
++ )
531 template<typename T
> TProb
<T
> max( const TProb
<T
> &a
, const TProb
<T
> &b
) {
532 assert( a
.size() == b
.size() );
533 TProb
<T
> result( a
.size() );
534 for( size_t i
= 0; i
< a
.size(); i
++ )
543 } // end of namespace dai