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
38 typedef std::complex<double> Complex
;
40 template<typename T
> class TProb
;
41 typedef TProb
<Real
> Prob
;
42 typedef TProb
<Complex
> CProb
;
45 /// TProb<T> implements a probability vector of type T.
46 /// T should be castable from and to double and to complex.
47 template <typename T
> class TProb
{
53 /// Calculate x times log(x), or 0 if x == 0
54 Complex
xlogx( Real x
) const { return( x
== 0.0 ? 0.0 : Complex(x
) * std::log(Complex(x
))); }
57 /// NORMPROB means that the sum of all entries should be 1
58 /// NORMLINF means that the maximum absolute value of all entries should be 1
59 typedef enum { NORMPROB
, NORMLINF
} NormType
;
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 typedef enum { DISTL1
, DISTLINF
, DISTTV
} DistType
;
65 /// Default constructor
68 /// Construct uniform distribution of given length
69 TProb( size_t n
) : _p(std::vector
<T
>(n
, 1.0 / n
)) {}
71 /// Construct with given length and initial value
72 TProb( size_t n
, Real p
) : _p(std::vector
<T
>(n
,(T
)p
)) {}
74 /// Construct with given length and initial array
75 TProb( size_t n
, const Real
* p
) {
76 // Reserve-push_back is faster than resize-copy
78 for( size_t i
= 0; i
< n
; i
++ )
83 TProb( const TProb
<T
> & x
) : _p(x
._p
) {}
85 /// Assignment operator
86 TProb
<T
> & operator=( const TProb
<T
> &x
) {
93 /// Provide read access to _p
94 const std::vector
<T
> & p() const { return _p
; }
96 /// Provide full access to _p
97 std::vector
<T
> & p() { return _p
; }
99 /// Provide read access to ith element of _p
100 T
operator[]( size_t i
) const { return _p
[i
]; }
102 /// Provide full access to ith element of _p
103 T
& operator[]( size_t i
) { return _p
[i
]; }
105 /// Set all elements to x
106 TProb
<T
> & fill(T x
) {
107 for( size_t i
= 0; i
< size(); i
++ )
112 /// Set all elements to iid random numbers from uniform(0,1) distribution
113 TProb
<T
> & randomize() {
114 for( size_t i
= 0; i
< size(); i
++ )
115 _p
[i
] = rnd_uniform();
120 size_t size() const {
124 /// Make entries zero if (Real) absolute value smaller than epsilon
125 TProb
<T
>& makeZero (Real epsilon
) {
126 for( size_t i
= 0; i
< size(); i
++ )
127 if( fabs((Real
)_p
[i
]) < epsilon
)
132 /// Multiplication with T x
133 TProb
<T
>& operator*= (T x
) {
134 for( size_t i
= 0; i
< size(); i
++ )
139 /// Return product of *this with T x
140 TProb
<T
> operator* (T x
) const {
141 TProb
<T
> prod( *this );
147 TProb
<T
>& operator/= (T x
) {
151 for( size_t i
= 0; i
< size(); i
++ )
156 /// Return quotient of *this and T x
157 TProb
<T
> operator/ (T x
) const {
158 TProb
<T
> prod( *this );
163 /// Pointwise multiplication with q
164 TProb
<T
>& operator*= (const TProb
<T
> & q
) {
166 assert( size() == q
.size() );
168 for( size_t i
= 0; i
< size(); i
++ )
173 /// Return product of *this with q
174 TProb
<T
> operator* (const TProb
<T
> & q
) const {
176 assert( size() == q
.size() );
178 TProb
<T
> prod( *this );
183 /// Pointwise addition with q
184 TProb
<T
>& operator+= (const TProb
<T
> & q
) {
186 assert( size() == q
.size() );
188 for( size_t i
= 0; i
< size(); i
++ )
193 /// Pointwise subtraction of q
194 TProb
<T
>& operator-= (const TProb
<T
> & q
) {
196 assert( size() == q
.size() );
198 for( size_t i
= 0; i
< size(); i
++ )
203 /// Return sum of *this and q
204 TProb
<T
> operator+ (const TProb
<T
> & q
) const {
206 assert( size() == q
.size() );
208 TProb
<T
> sum( *this );
213 /// Return *this minus q
214 TProb
<T
> operator- (const TProb
<T
> & q
) const {
216 assert( size() == q
.size() );
218 TProb
<T
> sum( *this );
223 /// Pointwise division by q
224 TProb
<T
>& operator/= (const TProb
<T
> & q
) {
226 assert( size() == q
.size() );
228 for( size_t i
= 0; i
< size(); i
++ ) {
230 // assert( q[i] != 0.0 );
232 if( q
[i
] == 0.0 ) // FIXME
240 /// Pointwise division by q, division by zero yields infinity
241 TProb
<T
>& divide (const TProb
<T
> & q
) {
243 assert( size() == q
.size() );
245 for( size_t i
= 0; i
< size(); i
++ )
250 /// Return quotient of *this with q
251 TProb
<T
> operator/ (const TProb
<T
> & q
) const {
253 assert( size() == q
.size() );
255 TProb
<T
> quot( *this );
260 /// Return pointwise inverse
261 TProb
<T
> inverse(bool zero
= false) const {
263 inv
._p
.reserve( size() );
265 for( size_t i
= 0; i
< size(); i
++ )
266 inv
._p
.push_back( _p
[i
] == 0.0 ? 0.0 : 1.0 / _p
[i
] );
268 for( size_t i
= 0; i
< size(); i
++ ) {
270 assert( _p
[i
] != 0.0 );
272 inv
._p
.push_back( 1.0 / _p
[i
] );
277 /// Return *this to the power of a (pointwise)
278 TProb
<T
>& operator^= (Real a
) {
280 for( size_t i
= 0; i
< size(); i
++ )
281 _p
[i
] = std::pow( _p
[i
], a
);
286 /// Pointwise power of a
287 TProb
<T
> operator^ (Real a
) const {
290 power
._p
.reserve( size() );
291 for( size_t i
= 0; i
< size(); i
++ )
292 power
._p
.push_back( std::pow( _p
[i
], a
) );
299 TProb
<T
> exp() const {
301 e
._p
.reserve( size() );
302 for( size_t i
= 0; i
< size(); i
++ )
303 e
._p
.push_back( std::exp( _p
[i
] ) );
308 TProb
<T
> log() const {
310 l
._p
.reserve( size() );
311 for( size_t i
= 0; i
< size(); i
++ )
312 l
._p
.push_back( std::log( _p
[i
] ) );
316 /// Pointwise log (or 0 if == 0)
317 TProb
<T
> log0() const {
319 l0
._p
.reserve( size() );
320 for( size_t i
= 0; i
< size(); i
++ )
321 l0
._p
.push_back( (_p
[i
] == 0.0) ? 0.0 : std::log( _p
[i
] ) );
325 /// Pointwise (complex) log (or 0 if == 0)
326 /* CProb clog0() const {
328 l0._p.reserve( size() );
329 for( size_t i = 0; i < size(); i++ )
330 l0._p.push_back( (_p[i] == 0.0) ? 0.0 : std::log( Complex( _p[i] ) ) );
334 /// Return distance of p and q
335 friend Real
dist( const TProb
<T
> & p
, const TProb
<T
> & q
, DistType dt
) {
337 assert( p
.size() == q
.size() );
342 for( size_t i
= 0; i
< p
.size(); i
++ )
343 result
+= fabs((Real
)p
[i
] - (Real
)q
[i
]);
347 for( size_t i
= 0; i
< p
.size(); i
++ ) {
348 Real z
= fabs((Real
)p
[i
] - (Real
)q
[i
]);
355 for( size_t i
= 0; i
< p
.size(); i
++ )
356 result
+= fabs((Real
)p
[i
] - (Real
)q
[i
]);
363 /// Return (complex) Kullback-Leibler distance with q
364 friend Complex
KL_dist( const TProb
<T
> & p
, const TProb
<T
> & q
) {
366 assert( p
.size() == q
.size() );
368 Complex result
= 0.0;
369 for( size_t i
= 0; i
< p
.size(); i
++ ) {
370 if( (Real
) p
[i
] != 0.0 ) {
373 result
+= p_i
* (std::log(p_i
) - std::log(q_i
));
379 /// Return sum of all entries
382 for( size_t i
= 0; i
< size(); i
++ )
387 /// Converts entries to Real and returns maximum absolute value
390 for( size_t i
= 0; i
< size(); i
++ ) {
391 Real mag
= fabs( (Real
) _p
[i
] );
398 /// Returns maximum value
401 for( size_t i
= 0; i
< size(); i
++ ) {
408 /// Normalize, using the specified norm
409 T
normalize( NormType norm
) {
411 if( norm
== NORMPROB
)
413 else if( norm
== NORMLINF
)
419 for( size_t i
= 0; i
< size(); i
++ )
424 /// Return normalized copy of *this, using the specified norm
425 TProb
<T
> normalized( NormType norm
) const {
427 if( norm
== NORMPROB
)
429 else if( norm
== NORMLINF
)
437 result
._p
.reserve( size() );
438 for( size_t i
= 0; i
< size(); i
++ )
439 result
._p
.push_back( _p
[i
] * Z
);
443 /// Returns true if one or more entries are NaN
444 bool hasNaNs() const {
446 for( size_t i
= 0; i
< size() && !NaNs
; i
++ )
452 /// Returns true if one or more entries are negative
453 bool hasNegatives() const {
454 bool Negatives
= false;
455 for( size_t i
= 0; i
< size() && !Negatives
; i
++ )
461 /// Returns (complex) entropy
462 Complex
entropy() const {
464 for( size_t i
= 0; i
< size(); i
++ )
469 friend std::ostream
& operator<< (std::ostream
& os
, const TProb
<T
>& P
) {
470 for( size_t i
= 0; i
< P
.size(); i
++ )
471 os
<< P
._p
[i
] << " ";
478 } // end of namespace dai