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 Copyright (C) 2002 Martijn Leisink [martijn@mbfys.kun.nl]
6 Radboud University Nijmegen, The Netherlands
8 This file is part of libDAI.
10 libDAI is free software; you can redistribute it and/or modify
11 it under the terms of the GNU General Public License as published by
12 the Free Software Foundation; either version 2 of the License, or
13 (at your option) any later version.
15 libDAI is distributed in the hope that it will be useful,
16 but WITHOUT ANY WARRANTY; without even the implied warranty of
17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 GNU General Public License for more details.
20 You should have received a copy of the GNU General Public License
21 along with libDAI; if not, write to the Free Software
22 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
27 /// \brief Defines TFactor<T> and Factor classes
30 #ifndef __defined_libdai_factor_h
31 #define __defined_libdai_factor_h
37 #include <dai/varset.h>
38 #include <dai/index.h>
44 /// Represents a (probability) factor.
45 /** Mathematically, a \e factor is a function mapping joint states of some
46 * variables to the nonnegative real numbers.
47 * More formally, denoting a discrete variable with label \f$l\f$ by
48 * \f$x_l\f$ and its state space by \f$X_l = \{0,1,\dots,S_l-1\}\f$,
49 * then a factor depending on the variables \f$\{x_l\}_{l\in L}\f$ is
50 * a function \f$f_L : \prod_{l\in L} X_l \to [0,\infty)\f$.
52 * In libDAI, a factor is represented by a TFactor<\a T> object, which has two
54 * \arg a VarSet, corresponding with the set of variables \f$\{x_l\}_{l\in L}\f$
55 * that the factor depends on;
56 * \arg a TProb<\a T>, a vector containing the value of the factor for each possible
57 * joint state of the variables.
59 * The factor values are stored in the entries of the TProb<\a T> in a particular
60 * ordering, which is defined by the one-to-one correspondence of a joint state
61 * in \f$\prod_{l\in L} X_l\f$ with a linear index in
62 * \f$\{0,1,\dots,\prod_{l\in L} S_l-1\}\f$ according to the mapping \f$\sigma\f$
63 * induced by VarSet::calcState(const std::map<Var,size_t> &).
65 * \tparam T Should be a scalar that is castable from and to double and should support elementary arithmetic operations.
66 * \todo Define a better fileformat for .fg files (maybe using XML)?
67 * \todo Add support for sparse factors.
69 template <typename T
> class TFactor
{
75 /// Constructs TFactor depending on no variables, with value p
76 TFactor ( Real p
= 1.0 ) : _vs(), _p(1,p
) {}
78 /// Constructs TFactor depending on variables in ns, with uniform distribution
79 TFactor( const VarSet
& ns
) : _vs(ns
), _p(_vs
.nrStates()) {}
81 /// Constructs TFactor depending on variables in ns, with all values set to p
82 TFactor( const VarSet
& ns
, Real p
) : _vs(ns
), _p(_vs
.nrStates(),p
) {}
84 /// Constructs TFactor depending on variables in ns, copying the values from the range starting at begin
85 /** \param ns contains the variables that the new TFactor should depend on.
86 * \tparam Iterator Iterates over instances of type T; should support addition of size_t.
87 * \param begin Points to first element to be added.
89 template<typename TIterator
>
90 TFactor( const VarSet
& ns
, TIterator begin
) : _vs(ns
), _p(begin
, begin
+ _vs
.nrStates(), _vs
.nrStates()) {}
92 /// Constructs TFactor depending on variables in ns, with values set to the TProb p
93 TFactor( const VarSet
& ns
, const TProb
<T
>& p
) : _vs(ns
), _p(p
) {
95 assert( _vs
.nrStates() == _p
.size() );
99 /// Constructs TFactor depending on the variable n, with uniform distribution
100 TFactor( const Var
& n
) : _vs(n
), _p(n
.states()) {}
103 TFactor( const TFactor
<T
> &x
) : _vs(x
._vs
), _p(x
._p
) {}
105 /// Assignment operator
106 TFactor
<T
> & operator= (const TFactor
<T
> &x
) {
114 /// Returns const reference to value vector
115 const TProb
<T
> & p() const { return _p
; }
116 /// Returns reference to value vector
117 TProb
<T
> & p() { return _p
; }
119 /// Returns const reference to variable set
120 const VarSet
& vars() const { return _vs
; }
122 /// Returns the number of possible joint states of the variables
123 /** \note This is equal to the length of the value vector.
125 size_t states() const { return _p
.size(); }
127 /// Returns a copy of the i'th entry of the value vector
128 T
operator[] (size_t i
) const { return _p
[i
]; }
130 /// Returns a reference to the i'th entry of the value vector
131 T
& operator[] (size_t i
) { return _p
[i
]; }
133 /// Sets all values to p
134 TFactor
<T
> & fill (T p
) { _p
.fill( p
); return(*this); }
136 /// Draws all values i.i.d. from a uniform distribution on [0,1)
137 TFactor
<T
> & randomize () { _p
.randomize(); return(*this); }
140 /// Multiplies *this with scalar t
141 TFactor
<T
>& operator*= (T t
) {
146 /// Divides *this by scalar t
147 TFactor
<T
>& operator/= (T t
) {
152 /// Adds scalar t to *this
153 TFactor
<T
>& operator+= (T t
) {
158 /// Subtracts scalar t from *this
159 TFactor
<T
>& operator-= (T t
) {
164 /// Raises *this to the power a
165 TFactor
<T
>& operator^= (Real a
) { _p
^= a
; return *this; }
168 /// Returns product of *this with scalar t
169 TFactor
<T
> operator* (T t
) const {
170 TFactor
<T
> result
= *this;
175 /// Returns quotient of *this with scalar t
176 TFactor
<T
> operator/ (T t
) const {
177 TFactor
<T
> result
= *this;
182 /// Returns sum of *this and scalar t
183 TFactor
<T
> operator+ (T t
) const {
184 TFactor
<T
> result(*this);
189 /// Returns *this minus scalar t
190 TFactor
<T
> operator- (T t
) const {
191 TFactor
<T
> result(*this);
196 /// Returns *this raised to the power a
197 TFactor
<T
> operator^ (Real a
) const {
204 /// Multiplies *this with the TFactor f
205 TFactor
<T
>& operator*= (const TFactor
<T
>& f
) {
206 if( f
._vs
== _vs
) // optimize special case
213 /// Divides *this by the TFactor f
214 TFactor
<T
>& operator/= (const TFactor
<T
>& f
) {
215 if( f
._vs
== _vs
) // optimize special case
222 /// Returns product of *this with the TFactor f
223 /** The product of two factors is defined as follows: if
224 * \f$f : \prod_{l\in L} X_l \to [0,\infty)\f$ and \f$g : \prod_{m\in M} X_m \to [0,\infty)\f$, then
225 * \f[fg : \prod_{l\in L\cup M} X_l \to [0,\infty) : x \mapsto f(x_L) g(x_M).\f]
227 TFactor
<T
> operator* (const TFactor
<T
>& f
) const;
229 /// Returns quotient of *this by the TFactor f
230 /** The quotient of two factors is defined as follows: if
231 * \f$f : \prod_{l\in L} X_l \to [0,\infty)\f$ and \f$g : \prod_{m\in M} X_m \to [0,\infty)\f$, then
232 * \f[\frac{f}{g} : \prod_{l\in L\cup M} X_l \to [0,\infty) : x \mapsto \frac{f(x_L)}{g(x_M)}.\f]
234 TFactor
<T
> operator/ (const TFactor
<T
>& f
) const;
236 /// Adds the TFactor f to *this
237 /** \pre this->vars() == f.vars()
239 TFactor
<T
>& operator+= (const TFactor
<T
>& f
) {
241 assert( f
._vs
== _vs
);
247 /// Subtracts the TFactor f from *this
248 /** \pre this->vars() == f.vars()
250 TFactor
<T
>& operator-= (const TFactor
<T
>& f
) {
252 assert( f
._vs
== _vs
);
258 /// Returns sum of *this and the TFactor f
259 /** \pre this->vars() == f.vars()
261 TFactor
<T
> operator+ (const TFactor
<T
>& f
) const {
263 assert( f
._vs
== _vs
);
265 TFactor
<T
> sum(*this);
270 /// Returns *this minus the TFactor f
271 /** \pre this->vars() == f.vars()
273 TFactor
<T
> operator- (const TFactor
<T
>& f
) const {
275 assert( f
._vs
== _vs
);
277 TFactor
<T
> sum(*this);
283 /// Sets all values that are smaller than epsilon to 0
284 TFactor
<T
>& makeZero( T epsilon
) {
285 _p
.makeZero( epsilon
);
289 /// Sets all values that are smaller than epsilon to epsilon
290 TFactor
<T
>& makePositive( T epsilon
) {
291 _p
.makePositive( epsilon
);
295 /// Returns pointwise inverse of *this.
296 /** If zero == true, uses 1 / 0 == 0; otherwise 1 / 0 == Inf.
298 TFactor
<T
> inverse(bool zero
=true) const {
301 inv
._p
= _p
.inverse(zero
);
305 /// Returns pointwise exp of *this
306 TFactor
<T
> exp() const {
313 /// Returns pointwise logarithm of *this
314 /** If zero==true, uses log(0)==0; otherwise, log(0)=-Inf.
316 TFactor
<T
> log(bool zero
=false) const {
323 /// Returns pointwise absolute value of *this
324 TFactor
<T
> abs() const {
331 /// Normalizes *this TFactor according to the specified norm
332 T
normalize( typename
Prob::NormType norm
=Prob::NORMPROB
) { return _p
.normalize( norm
); }
334 /// Returns a normalized copy of *this, according to the specified norm
335 TFactor
<T
> normalized( typename
Prob::NormType norm
=Prob::NORMPROB
) const {
338 result
._p
= _p
.normalized( norm
);
342 /// Returns a slice of this TFactor, where the subset ns is in state nsState
343 /** \pre \a ns sould be a subset of vars()
344 * \pre \a nsState < ns.states()
346 * The result is a TFactor that depends on the variables in this->vars() except those in \a ns,
347 * obtained by setting the variables in \a ns to the joint state specified by the linear index
348 * \a nsState. Formally, if *this corresponds with the factor \f$f : \prod_{l\in L} X_l \to [0,\infty)\f$,
349 * \f$M \subset L\f$ corresponds with \a ns and \a nsState corresponds with a mapping \f$s\f$ that
350 * maps a variable \f$x_m\f$ with \f$m\in M\f$ to its state \f$s(x_m) \in X_m\f$, then the slice
351 * returned corresponds with the factor \f$g : \prod_{l \in L \setminus M} X_l \to [0,\infty)\f$
352 * defined by \f$g(\{x_l\}_{l\in L \setminus M}) = f(\{x_l\}_{l\in L \setminus M}, \{s(x_m)\}_{m\in M})\f$.
354 TFactor
<T
> slice( const VarSet
& ns
, size_t nsState
) const {
356 VarSet nsrem
= _vs
/ ns
;
357 TFactor
<T
> result( nsrem
, T(0) );
360 IndexFor
i_ns (ns
, _vs
);
361 IndexFor
i_nsrem (nsrem
, _vs
);
362 for( size_t i
= 0; i
< states(); i
++, ++i_ns
, ++i_nsrem
)
363 if( (size_t)i_ns
== nsState
)
364 result
._p
[i_nsrem
] = _p
[i
];
369 /// Returns marginal on ns, obtained by summing out all variables except those in ns, and normalizing the result if normed==true
370 TFactor
<T
> marginal(const VarSet
& ns
, bool normed
=true) const;
372 /// Embeds this factor in a larger VarSet
373 /** \pre vars() should be a subset of ns
375 * If *this corresponds with \f$f : \prod_{l\in L} X_l \to [0,\infty)\f$ and \f$L \subset M\f$, then
376 * the embedded factor corresponds with \f$g : \prod_{m\in M} X_m \to [0,\infty) : x \mapsto f(x_L)\f$.
378 TFactor
<T
> embed(const VarSet
& ns
) const {
383 return (*this) * TFactor
<T
>(ns
/ _vs
, 1);
386 /// Returns true if *this has NaN values
387 bool hasNaNs() const { return _p
.hasNaNs(); }
389 /// Returns true if *this has negative values
390 bool hasNegatives() const { return _p
.hasNegatives(); }
392 /// Returns total sum of values
393 T
totalSum() const { return _p
.totalSum(); }
395 /// Returns maximum absolute value
396 T
maxAbs() const { return _p
.maxAbs(); }
398 /// Returns maximum value
399 T
maxVal() const { return _p
.maxVal(); }
401 /// Returns minimum value
402 T
minVal() const { return _p
.minVal(); }
404 /// Returns entropy of *this TFactor
405 Real
entropy() const { return _p
.entropy(); }
407 /// Returns strength of *this TFactor (between variables i and j), as defined in eq. (52) of [\ref MoK07b]
408 T
strength( const Var
&i
, const Var
&j
) const;
412 template<typename T
> TFactor
<T
> TFactor
<T
>::marginal(const VarSet
& ns
, bool normed
) const {
413 VarSet res_ns
= ns
& _vs
;
415 TFactor
<T
> res( res_ns
, 0.0 );
417 IndexFor
i_res( res_ns
, _vs
);
418 for( size_t i
= 0; i
< _p
.size(); i
++, ++i_res
)
419 res
._p
[i_res
] += _p
[i
];
422 res
.normalize( Prob::NORMPROB
);
428 template<typename T
> TFactor
<T
> TFactor
<T
>::operator* (const TFactor
<T
>& f
) const {
429 if( f
._vs
== _vs
) { // optimizate special case
430 TFactor
<T
> prod(*this);
434 TFactor
<T
> prod( _vs
| f
._vs
, 0.0 );
436 IndexFor
i1(_vs
, prod
._vs
);
437 IndexFor
i2(f
._vs
, prod
._vs
);
439 for( size_t i
= 0; i
< prod
._p
.size(); i
++, ++i1
, ++i2
)
440 prod
._p
[i
] += _p
[i1
] * f
._p
[i2
];
447 template<typename T
> TFactor
<T
> TFactor
<T
>::operator/ (const TFactor
<T
>& f
) const {
448 if( f
._vs
== _vs
) { // optimizate special case
449 TFactor
<T
> quot(*this);
453 TFactor
<T
> quot( _vs
| f
._vs
, 0.0 );
455 IndexFor
i1(_vs
, quot
._vs
);
456 IndexFor
i2(f
._vs
, quot
._vs
);
458 for( size_t i
= 0; i
< quot
._p
.size(); i
++, ++i1
, ++i2
)
459 quot
._p
[i
] += _p
[i1
] / f
._p
[i2
];
466 template<typename T
> T TFactor
<T
>::strength( const Var
&i
, const Var
&j
) const {
468 assert( _vs
.contains( i
) );
469 assert( _vs
.contains( j
) );
475 for( size_t alpha1
= 0; alpha1
< i
.states(); alpha1
++ )
476 for( size_t alpha2
= 0; alpha2
< i
.states(); alpha2
++ )
477 if( alpha2
!= alpha1
)
478 for( size_t beta1
= 0; beta1
< j
.states(); beta1
++ )
479 for( size_t beta2
= 0; beta2
< j
.states(); beta2
++ )
480 if( beta2
!= beta1
) {
481 size_t as
= 1, bs
= 1;
486 T f1
= slice( ij
, alpha1
* as
+ beta1
* bs
).p().divide( slice( ij
, alpha2
* as
+ beta1
* bs
).p() ).maxVal();
487 T f2
= slice( ij
, alpha2
* as
+ beta2
* bs
).p().divide( slice( ij
, alpha1
* as
+ beta2
* bs
).p() ).maxVal();
493 return std::tanh( 0.25 * std::log( max
) );
497 /// Writes a TFactor to an output stream
500 template<typename T
> std::ostream
& operator<< (std::ostream
& os
, const TFactor
<T
>& P
) {
501 os
<< "(" << P
.vars() << ", (";
502 for( size_t i
= 0; i
< P
.states(); i
++ )
503 os
<< (i
== 0 ? "" : ", ") << P
[i
];
509 /// Returns distance between two TFactors f and g, according to the distance measure dt
511 * \pre f.vars() == g.vars()
513 template<typename T
> Real
dist( const TFactor
<T
> &f
, const TFactor
<T
> &g
, Prob::DistType dt
) {
514 if( f
.vars().empty() || g
.vars().empty() )
518 assert( f
.vars() == g
.vars() );
520 return dist( f
.p(), g
.p(), dt
);
525 /// Returns the pointwise maximum of two TFactors
527 * \pre f.vars() == g.vars()
529 template<typename T
> TFactor
<T
> max( const TFactor
<T
> &f
, const TFactor
<T
> &g
) {
530 assert( f
._vs
== g
._vs
);
531 return TFactor
<T
>( f
._vs
, min( f
.p(), g
.p() ) );
535 /// Returns the pointwise minimum of two TFactors
537 * \pre f.vars() == g.vars()
539 template<typename T
> TFactor
<T
> min( const TFactor
<T
> &f
, const TFactor
<T
> &g
) {
540 assert( f
._vs
== g
._vs
);
541 return TFactor
<T
>( f
._vs
, max( f
.p(), g
.p() ) );
545 /// Calculates the mutual information between the two variables that f depends on, under the distribution given by f
547 * \pre f.vars().size() == 2
549 template<typename T
> Real
MutualInfo(const TFactor
<T
> &f
) {
550 assert( f
.vars().size() == 2 );
551 VarSet::const_iterator it
= f
.vars().begin();
552 Var i
= *it
; it
++; Var j
= *it
;
553 TFactor
<T
> projection
= f
.marginal(i
) * f
.marginal(j
);
554 return real( dist( f
.normalized(), projection
, Prob::DISTKL
) );
558 /// Represents a factor with values of type Real.
559 typedef TFactor
<Real
> Factor
;
562 } // end of namespace dai