1 /* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
2 Radboud University Nijmegen, The Netherlands
4 This file is part of libDAI.
6 libDAI is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2 of the License, or
9 (at your option) any later version.
11 libDAI is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
16 You should have received a copy of the GNU General Public License
17 along with libDAI; if not, write to the Free Software
18 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
22 #ifndef __defined_libdai_prob_h
23 #define __defined_libdai_prob_h
41 typedef std::complex<double> Complex
;
43 template<typename T
> class TProb
;
44 typedef TProb
<Real
> Prob
;
45 typedef TProb
<Complex
> CProb
;
48 /// TProb<T> implements a probability vector of type T.
49 /// T should be castable from and to double and to complex.
50 template <typename T
> class TProb
{
56 /// Calculate x times log(x), or 0 if x == 0
57 Complex
xlogx( Real x
) const { return( x
== 0.0 ? 0.0 : Complex(x
) * std::log(Complex(x
))); }
60 /// NORMPROB means that the sum of all entries should be 1
61 /// NORMLINF means that the maximum absolute value of all entries should be 1
62 typedef enum { NORMPROB
, NORMLINF
} NormType
;
63 /// DISTL1 is the L-1 distance (sum of absolute values of pointwise difference)
64 /// DISTLINF is the L-inf distance (maximum absolute value of pointwise difference)
65 /// DISTTV is the Total Variation distance
66 typedef enum { DISTL1
, DISTLINF
, DISTTV
} DistType
;
68 /// Default constructor
71 /// Construct uniform distribution of given length
72 explicit TProb( size_t n
) : _p(std::vector
<T
>(n
, 1.0 / n
)) {}
74 /// Construct with given length and initial value
75 TProb( size_t n
, Real p
) : _p(n
, (T
)p
) {}
77 /// Construct with given length and initial array
78 TProb( size_t n
, const Real
* p
) : _p(p
, p
+ n
) {}
80 /// Provide read access to _p
81 const std::vector
<T
> & p() const { return _p
; }
83 /// Provide full access to _p
84 std::vector
<T
> & p() { return _p
; }
86 /// Provide read access to ith element of _p
87 T
operator[]( size_t i
) const { return _p
[i
]; }
89 /// Provide full access to ith element of _p
90 T
& operator[]( size_t i
) { return _p
[i
]; }
92 /// Set all elements to x
93 TProb
<T
> & fill(T x
) {
94 std::fill( _p
.begin(), _p
.end(), x
);
98 /// Set all elements to iid random numbers from uniform(0,1) distribution
99 TProb
<T
> & randomize() {
100 std::generate(_p
.begin(), _p
.end(), rnd_uniform
);
105 size_t size() const {
109 /// Make entries zero if (Real) absolute value smaller than epsilon
110 TProb
<T
>& makeZero (Real epsilon
) {
111 for( size_t i
= 0; i
< size(); i
++ )
112 if( fabs((Real
)_p
[i
]) < epsilon
)
114 // std::replace_if( _p.begin(), _p.end(), fabs((Real)boost::lambda::_1) < epsilon, 0.0 );
118 /// Multiplication with T x
119 TProb
<T
>& operator*= (T x
) {
120 std::transform( _p
.begin(), _p
.end(), _p
.begin(), std::bind2nd( std::multiplies
<T
>(), x
) );
124 /// Return product of *this with T x
125 TProb
<T
> operator* (T x
) const {
126 TProb
<T
> prod( *this );
132 TProb
<T
>& operator/= (T x
) {
136 std::transform( _p
.begin(), _p
.end(), _p
.begin(), std::bind2nd( std::divides
<T
>(), x
) );
140 /// Return quotient of *this and T x
141 TProb
<T
> operator/ (T x
) const {
142 TProb
<T
> quot( *this );
148 TProb
<T
>& operator+= (T x
) {
149 std::transform( _p
.begin(), _p
.end(), _p
.begin(), std::bind2nd( std::plus
<T
>(), x
) );
153 /// Return sum of *this with T x
154 TProb
<T
> operator+ (T x
) const {
155 TProb
<T
> sum( *this );
160 /// Difference by T x
161 TProb
<T
>& operator-= (T x
) {
162 std::transform( _p
.begin(), _p
.end(), _p
.begin(), std::bind2nd( std::minus
<T
>(), x
) );
166 /// Return difference of *this and T x
167 TProb
<T
> operator- (T x
) const {
168 TProb
<T
> diff( *this );
173 /// Pointwise multiplication with q
174 TProb
<T
>& operator*= (const TProb
<T
> & q
) {
176 assert( size() == q
.size() );
178 std::transform( _p
.begin(), _p
.end(), q
._p
.begin(), _p
.begin(), std::multiplies
<T
>() );
182 /// Return product of *this with q
183 TProb
<T
> operator* (const TProb
<T
> & q
) const {
185 assert( size() == q
.size() );
187 TProb
<T
> prod( *this );
192 /// Pointwise addition with q
193 TProb
<T
>& operator+= (const TProb
<T
> & q
) {
195 assert( size() == q
.size() );
197 std::transform( _p
.begin(), _p
.end(), q
._p
.begin(), _p
.begin(), std::plus
<T
>() );
201 /// Pointwise subtraction of 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::minus
<T
>() );
210 /// Return sum of *this and q
211 TProb
<T
> operator+ (const TProb
<T
> & q
) const {
213 assert( size() == q
.size() );
215 TProb
<T
> sum( *this );
220 /// Return *this minus q
221 TProb
<T
> operator- (const TProb
<T
> & q
) const {
223 assert( size() == q
.size() );
225 TProb
<T
> diff( *this );
230 /// Pointwise division by q (division by zero yields zero)
231 TProb
<T
>& operator/= (const TProb
<T
> & q
) {
233 assert( size() == q
.size() );
235 for( size_t i
= 0; i
< size(); i
++ ) {
237 // assert( q[i] != 0.0 );
247 /// Pointwise division by q (division by zero yields infinity)
248 TProb
<T
>& divide (const TProb
<T
> & q
) {
250 assert( size() == q
.size() );
252 std::transform( _p
.begin(), _p
.end(), q
._p
.begin(), _p
.begin(), std::divides
<T
>() );
256 /// Return quotient of *this with q
257 TProb
<T
> operator/ (const TProb
<T
> & q
) const {
259 assert( size() == q
.size() );
261 TProb
<T
> quot( *this );
266 /// Return pointwise inverse
267 TProb
<T
> inverse(bool zero
= false) const {
269 inv
._p
.reserve( size() );
271 for( size_t i
= 0; i
< size(); i
++ )
272 inv
._p
.push_back( _p
[i
] == 0.0 ? 0.0 : 1.0 / _p
[i
] );
274 for( size_t i
= 0; i
< size(); i
++ ) {
276 assert( _p
[i
] != 0.0 );
278 inv
._p
.push_back( 1.0 / _p
[i
] );
283 /// Return *this to the power of a (pointwise)
284 TProb
<T
>& operator^= (Real a
) {
286 std::transform( _p
.begin(), _p
.end(), _p
.begin(), std::bind2nd( std::ptr_fun
<T
, Real
, T
>(std::pow
), a
) );
290 /// Pointwise power of a
291 TProb
<T
> operator^ (Real a
) const {
292 TProb
<T
> power(*this);
298 const TProb
<T
>& takeExp() {
299 std::transform( _p
.begin(), _p
.end(), _p
.begin(), std::ptr_fun
<T
, T
>(std::exp
) );
304 const TProb
<T
>& takeLog() {
305 std::transform( _p
.begin(), _p
.end(), _p
.begin(), std::ptr_fun
<T
, T
>(std::log
) );
309 /// Pointwise log (or 0 if == 0)
310 const TProb
<T
>& takeLog0() {
311 for( size_t i
= 0; i
< size(); i
++ )
312 _p
[i
] = ( (_p
[i
] == 0.0) ? 0.0 : std::log( _p
[i
] ) );
317 TProb
<T
> exp() const {
324 TProb
<T
> log() const {
330 /// Pointwise log (or 0 if == 0)
331 TProb
<T
> log0() const {
337 /// Pointwise (complex) log (or 0 if == 0)
338 /* CProb clog0() const {
340 l0._p.reserve( size() );
341 for( size_t i = 0; i < size(); i++ )
342 l0._p.push_back( (_p[i] == 0.0) ? 0.0 : std::log( Complex( _p[i] ) ) );
346 /// Return distance of p and q
347 friend Real
dist( const TProb
<T
> & p
, const TProb
<T
> & q
, DistType dt
) {
349 assert( p
.size() == q
.size() );
354 for( size_t i
= 0; i
< p
.size(); i
++ )
355 result
+= fabs((Real
)p
[i
] - (Real
)q
[i
]);
359 for( size_t i
= 0; i
< p
.size(); i
++ ) {
360 Real z
= fabs((Real
)p
[i
] - (Real
)q
[i
]);
367 for( size_t i
= 0; i
< p
.size(); i
++ )
368 result
+= fabs((Real
)p
[i
] - (Real
)q
[i
]);
375 /// Return (complex) Kullback-Leibler distance with q
376 friend Complex
KL_dist( const TProb
<T
> & p
, const TProb
<T
> & q
) {
378 assert( p
.size() == q
.size() );
380 Complex result
= 0.0;
381 for( size_t i
= 0; i
< p
.size(); i
++ ) {
382 if( (Real
) p
[i
] != 0.0 ) {
385 result
+= p_i
* (std::log(p_i
) - std::log(q_i
));
391 /// Return sum of all entries
393 T Z
= std::accumulate( _p
.begin(), _p
.end(), (T
)0 );
397 /// Converts entries to Real and returns maximum absolute value
400 for( size_t i
= 0; i
< size(); i
++ ) {
401 Real mag
= fabs( (Real
) _p
[i
] );
408 /// Returns maximum value
410 T Z
= *std::max_element( _p
.begin(), _p
.end() );
414 /// Normalize, using the specified norm
415 T
normalize( NormType norm
) {
417 if( norm
== NORMPROB
)
419 else if( norm
== NORMLINF
)
428 /// Return normalized copy of *this, using the specified norm
429 TProb
<T
> normalized( NormType norm
) const {
430 TProb
<T
> result(*this);
431 result
.normalize( norm
);
435 /// Returns true if one or more entries are NaN
436 bool hasNaNs() const {
437 return (std::find_if( _p
.begin(), _p
.end(), isnan
) != _p
.end());
440 /// Returns true if one or more entries are negative
441 bool hasNegatives() const {
442 return (std::find_if( _p
.begin(), _p
.end(), std::bind2nd( std::less
<Real
>(), 0.0 ) ) != _p
.end());
445 /// Returns true if one or more entries are non-positive (causes problems with logscale)
446 bool hasNonPositives() const {
447 return (std::find_if( _p
.begin(), _p
.end(), std::bind2nd( std::less_equal
<Real
>(), 0.0 ) ) != _p
.end());
450 /// Returns (complex) entropy
451 Complex
entropy() const {
453 for( size_t i
= 0; i
< size(); i
++ )
458 friend std::ostream
& operator<< (std::ostream
& os
, const TProb
<T
>& P
) {
459 std::copy( P
._p
.begin(), P
._p
.end(), std::ostream_iterator
<T
>(os
, " ") );
466 } // end of namespace dai