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
38 #include <dai/varset.h>
39 #include <dai/index.h>
46 // Function object similar to std::divides(), but different in that dividing by zero results in zero
47 template<typename T
> struct divides0
: public std::binary_function
<T
, T
, T
> {
48 // Returns (j == 0 ? 0 : (i/j))
49 T
operator()( const T
&i
, const T
&j
) const {
58 /// Represents a (probability) factor.
59 /** Mathematically, a \e factor is a function mapping joint states of some
60 * variables to the nonnegative real numbers.
61 * More formally, denoting a discrete variable with label \f$l\f$ by
62 * \f$x_l\f$ and its state space by \f$X_l = \{0,1,\dots,S_l-1\}\f$,
63 * then a factor depending on the variables \f$\{x_l\}_{l\in L}\f$ is
64 * a function \f$f_L : \prod_{l\in L} X_l \to [0,\infty)\f$.
66 * In libDAI, a factor is represented by a TFactor<\a T> object, which has two
68 * \arg a VarSet, corresponding with the set of variables \f$\{x_l\}_{l\in L}\f$
69 * that the factor depends on;
70 * \arg a TProb<\a T>, a vector containing the value of the factor for each possible
71 * joint state of the variables.
73 * The factor values are stored in the entries of the TProb<\a T> in a particular
74 * ordering, which is defined by the one-to-one correspondence of a joint state
75 * in \f$\prod_{l\in L} X_l\f$ with a linear index in
76 * \f$\{0,1,\dots,\prod_{l\in L} S_l-1\}\f$ according to the mapping \f$\sigma\f$
77 * induced by VarSet::calcState(const std::map<Var,size_t> &).
79 * \tparam T Should be a scalar that is castable from and to double and should support elementary arithmetic operations.
80 * \todo Define a better fileformat for .fg files (maybe using XML)?
81 * \todo Add support for sparse factors.
83 template <typename T
> class TFactor
{
89 /// Iterator over factor entries
90 typedef typename TProb
<T
>::iterator iterator
;
92 /// Const iterator over factor entries
93 typedef typename TProb
<T
>::const_iterator const_iterator
;
95 /// Constructs TFactor depending on no variables, with value p
96 TFactor ( Real p
= 1.0 ) : _vs(), _p(1,p
) {}
98 /// Constructs TFactor depending on variables in vars, with uniform distribution
99 TFactor( const VarSet
& vars
) : _vs(vars
), _p(_vs
.nrStates()) {}
101 /// Constructs TFactor depending on variables in vars, with all values set to p
102 TFactor( const VarSet
& vars
, Real p
) : _vs(vars
), _p(_vs
.nrStates(),p
) {}
104 /// Constructs TFactor depending on variables in vars, copying the values from the range starting at begin
105 /** \param vars contains the variables that the new TFactor should depend on.
106 * \tparam Iterator Iterates over instances of type T; should support addition of size_t.
107 * \param begin Points to first element to be added.
109 template<typename TIterator
>
110 TFactor( const VarSet
& vars
, TIterator begin
) : _vs(vars
), _p(begin
, begin
+ _vs
.nrStates(), _vs
.nrStates()) {}
112 /// Constructs TFactor depending on variables in vars, with values set to the TProb p
113 TFactor( const VarSet
& vars
, const TProb
<T
> &p
) : _vs(vars
), _p(p
) {
114 DAI_DEBASSERT( _vs
.nrStates() == _p
.size() );
117 /// Constructs TFactor depending on variables in vars, permuting the values given in TProb p
118 TFactor( const std::vector
<Var
> &vars
, const std::vector
<T
> &p
) : _vs(vars
.begin(), vars
.end(), vars
.size()), _p(p
.size()) {
119 Permute
permindex(vars
);
120 for( size_t li
= 0; li
< p
.size(); ++li
)
121 _p
[permindex
.convert_linear_index(li
)] = p
[li
];
124 /// Constructs TFactor depending on the variable v, with uniform distribution
125 TFactor( const Var
&v
) : _vs(v
), _p(v
.states()) {}
127 /// Returns const reference to value vector
128 const TProb
<T
>& p() const { return _p
; }
129 /// Returns reference to value vector
130 TProb
<T
>& p() { return _p
; }
132 /// Returns const reference to variable set
133 const VarSet
& vars() const { return _vs
; }
135 /// Returns the number of possible joint states of the variables
136 /** \note This is equal to the length of the value vector.
138 size_t states() const { return _p
.size(); }
140 /// Returns a copy of the i'th entry of the value vector
141 T
operator[] (size_t i
) const { return _p
[i
]; }
143 /// Returns a reference to the i'th entry of the value vector
144 T
& operator[] (size_t i
) { return _p
[i
]; }
146 /// Returns iterator pointing to first entry
147 iterator
begin() { return _p
.begin(); }
148 /// Returns const iterator pointing to first entry
149 const_iterator
begin() const { return _p
.begin(); }
150 /// Returns iterator pointing beyond last entry
151 iterator
end() { return _p
.end(); }
152 /// Returns const iterator pointing beyond last entry
153 const_iterator
end() const { return _p
.end(); }
155 /// Sets all values to p
156 TFactor
<T
> & fill (T p
) { _p
.fill( p
); return(*this); }
158 /// Draws all values i.i.d. from a uniform distribution on [0,1)
159 TFactor
<T
> & randomize () { _p
.randomize(); return(*this); }
162 /// Multiplies *this with scalar t
163 TFactor
<T
>& operator*= (T t
) {
168 /// Divides *this by scalar t
169 TFactor
<T
>& operator/= (T t
) {
174 /// Adds scalar t to *this
175 TFactor
<T
>& operator+= (T t
) {
180 /// Subtracts scalar t from *this
181 TFactor
<T
>& operator-= (T t
) {
186 /// Raises *this to the power a
187 TFactor
<T
>& operator^= (Real a
) { _p
^= a
; return *this; }
190 /// Returns product of *this with scalar t
191 TFactor
<T
> operator* (T t
) const {
192 TFactor
<T
> result
= *this;
197 /// Returns quotient of *this with scalar t
198 TFactor
<T
> operator/ (T t
) const {
199 TFactor
<T
> result
= *this;
204 /// Returns sum of *this and scalar t
205 TFactor
<T
> operator+ (T t
) const {
206 TFactor
<T
> result(*this);
211 /// Returns *this minus scalar t
212 TFactor
<T
> operator- (T t
) const {
213 TFactor
<T
> result(*this);
218 /// Returns *this raised to the power a
219 TFactor
<T
> operator^ (Real a
) const {
226 /// Multiplies *this with the TFactor f
227 TFactor
<T
>& operator*= (const TFactor
<T
>& f
) {
228 if( f
._vs
== _vs
) // optimize special case
235 /// Divides *this by the TFactor f
236 TFactor
<T
>& operator/= (const TFactor
<T
>& f
) {
237 if( f
._vs
== _vs
) // optimize special case
244 /// Adds the TFactor f to *this
245 TFactor
<T
>& operator+= (const TFactor
<T
>& f
) {
246 if( f
._vs
== _vs
) // optimize special case
253 /// Subtracts the TFactor f from *this
254 TFactor
<T
>& operator-= (const TFactor
<T
>& f
) {
255 if( f
._vs
== _vs
) // optimize special case
262 /// Returns product of *this with the TFactor f
263 /** The product of two factors is defined as follows: if
264 * \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
265 * \f[fg : \prod_{l\in L\cup M} X_l \to [0,\infty) : x \mapsto f(x_L) g(x_M).\f]
267 TFactor
<T
> operator* (const TFactor
<T
>& f
) const {
268 return pointwiseOp(*this,f
,std::multiplies
<T
>());
271 /// Returns quotient of *this by the TFactor f
272 /** The quotient of two factors is defined as follows: if
273 * \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
274 * \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]
276 TFactor
<T
> operator/ (const TFactor
<T
>& f
) const {
277 return pointwiseOp(*this,f
,divides0
<T
>());
280 /// Returns sum of *this and the TFactor f
281 /** The sum of two factors is defined as follows: if
282 * \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
283 * \f[f+g : \prod_{l\in L\cup M} X_l \to [0,\infty) : x \mapsto f(x_L) + g(x_M).\f]
285 TFactor
<T
> operator+ (const TFactor
<T
>& f
) const {
286 return pointwiseOp(*this,f
,std::plus
<T
>());
289 /// Returns *this minus the TFactor f
290 /** The difference of two factors is defined as follows: if
291 * \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
292 * \f[f-g : \prod_{l\in L\cup M} X_l \to [0,\infty) : x \mapsto f(x_L) - g(x_M).\f]
294 TFactor
<T
> operator- (const TFactor
<T
>& f
) const {
295 return pointwiseOp(*this,f
,std::minus
<T
>());
299 /// Sets all values that are smaller than epsilon to 0
300 TFactor
<T
>& makeZero( T epsilon
) {
301 _p
.makeZero( epsilon
);
305 /// Sets all values that are smaller than epsilon to epsilon
306 TFactor
<T
>& makePositive( T epsilon
) {
307 _p
.makePositive( epsilon
);
311 /// Returns pointwise inverse of *this.
312 /** If zero == true, uses 1 / 0 == 0; otherwise 1 / 0 == Inf.
314 TFactor
<T
> inverse(bool zero
=true) const {
317 inv
._p
= _p
.inverse(zero
);
321 /// Returns pointwise exp of *this
322 TFactor
<T
> exp() const {
329 /// Returns pointwise logarithm of *this
330 /** If zero==true, uses log(0)==0; otherwise, log(0)=-Inf.
332 TFactor
<T
> log(bool zero
=false) const {
339 /// Returns pointwise absolute value of *this
340 TFactor
<T
> abs() const {
347 /// Normalizes *this TFactor according to the specified norm
348 T
normalize( typename
Prob::NormType norm
=Prob::NORMPROB
) { return _p
.normalize( norm
); }
350 /// Returns a normalized copy of *this, according to the specified norm
351 TFactor
<T
> normalized( typename
Prob::NormType norm
=Prob::NORMPROB
) const {
354 result
._p
= _p
.normalized( norm
);
358 /// Returns a slice of this TFactor, where the subset ns is in state nsState
359 /** \pre \a ns sould be a subset of vars()
360 * \pre \a nsState < ns.states()
362 * The result is a TFactor that depends on the variables in this->vars() except those in \a ns,
363 * obtained by setting the variables in \a ns to the joint state specified by the linear index
364 * \a nsState. Formally, if *this corresponds with the factor \f$f : \prod_{l\in L} X_l \to [0,\infty)\f$,
365 * \f$M \subset L\f$ corresponds with \a ns and \a nsState corresponds with a mapping \f$s\f$ that
366 * 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
367 * returned corresponds with the factor \f$g : \prod_{l \in L \setminus M} X_l \to [0,\infty)\f$
368 * 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$.
370 TFactor
<T
> slice( const VarSet
& ns
, size_t nsState
) const {
372 VarSet nsrem
= _vs
/ ns
;
373 TFactor
<T
> result( nsrem
, T(0) );
376 IndexFor
i_ns (ns
, _vs
);
377 IndexFor
i_nsrem (nsrem
, _vs
);
378 for( size_t i
= 0; i
< states(); i
++, ++i_ns
, ++i_nsrem
)
379 if( (size_t)i_ns
== nsState
)
380 result
._p
[i_nsrem
] = _p
[i
];
385 /// Returns marginal on ns, obtained by summing out all variables except those in ns, and normalizing the result if normed==true
386 TFactor
<T
> marginal(const VarSet
& ns
, bool normed
=true) const;
388 /// Returns max-marginal on ns, obtained by maximizing all variables except those in ns, and normalizing the result if normed==true
389 TFactor
<T
> maxMarginal(const VarSet
& ns
, bool normed
=true) const;
391 /// Embeds this factor in a larger VarSet
392 /** \pre vars() should be a subset of ns
394 * If *this corresponds with \f$f : \prod_{l\in L} X_l \to [0,\infty)\f$ and \f$L \subset M\f$, then
395 * the embedded factor corresponds with \f$g : \prod_{m\in M} X_m \to [0,\infty) : x \mapsto f(x_L)\f$.
397 TFactor
<T
> embed(const VarSet
& ns
) const {
402 return (*this) * TFactor
<T
>(ns
/ _vs
, (T
)1);
405 /// Returns true if *this has NaN values
406 bool hasNaNs() const { return _p
.hasNaNs(); }
408 /// Returns true if *this has negative values
409 bool hasNegatives() const { return _p
.hasNegatives(); }
411 /// Returns total sum of values
412 T
sum() const { return _p
.sum(); }
414 /// Returns maximum absolute value
415 T
maxAbs() const { return _p
.maxAbs(); }
417 /// Returns maximum value
418 T
max() const { return _p
.max(); }
420 /// Returns minimum value
421 T
min() const { return _p
.min(); }
423 /// Returns entropy of *this TFactor
424 Real
entropy() const { return _p
.entropy(); }
426 /// Returns strength of *this TFactor (between variables i and j), as defined in eq. (52) of [\ref MoK07b]
427 T
strength( const Var
&i
, const Var
&j
) const;
431 template<typename T
> TFactor
<T
> TFactor
<T
>::marginal(const VarSet
& ns
, bool normed
) const {
432 VarSet res_ns
= ns
& _vs
;
434 TFactor
<T
> res( res_ns
, 0.0 );
436 IndexFor
i_res( res_ns
, _vs
);
437 for( size_t i
= 0; i
< _p
.size(); i
++, ++i_res
)
438 res
._p
[i_res
] += _p
[i
];
441 res
.normalize( Prob::NORMPROB
);
447 template<typename T
> TFactor
<T
> TFactor
<T
>::maxMarginal(const VarSet
& ns
, bool normed
) const {
448 VarSet res_ns
= ns
& _vs
;
450 TFactor
<T
> res( res_ns
, 0.0 );
452 IndexFor
i_res( res_ns
, _vs
);
453 for( size_t i
= 0; i
< _p
.size(); i
++, ++i_res
)
454 if( _p
[i
] > res
._p
[i_res
] )
455 res
._p
[i_res
] = _p
[i
];
458 res
.normalize( Prob::NORMPROB
);
464 /// Apply binary operator pointwise on two factors
465 template<typename T
, typename binaryOp
> TFactor
<T
> pointwiseOp( const TFactor
<T
> &f
, const TFactor
<T
> &g
, binaryOp op
) {
466 if( f
.vars() == g
.vars() ) { // optimizate special case
467 TFactor
<T
> result(f
);
468 for( size_t i
= 0; i
< result
.states(); i
++ )
469 result
[i
] = op( result
[i
], g
[i
] );
472 TFactor
<T
> result( f
.vars() | g
.vars(), 0.0 );
474 IndexFor
i1(f
.vars(), result
.vars());
475 IndexFor
i2(g
.vars(), result
.vars());
477 for( size_t i
= 0; i
< result
.states(); i
++, ++i1
, ++i2
)
478 result
[i
] = op( f
[i1
], g
[i2
] );
485 template<typename T
> T TFactor
<T
>::strength( const Var
&i
, const Var
&j
) const {
486 DAI_DEBASSERT( _vs
.contains( i
) );
487 DAI_DEBASSERT( _vs
.contains( j
) );
488 DAI_DEBASSERT( i
!= j
);
492 for( size_t alpha1
= 0; alpha1
< i
.states(); alpha1
++ )
493 for( size_t alpha2
= 0; alpha2
< i
.states(); alpha2
++ )
494 if( alpha2
!= alpha1
)
495 for( size_t beta1
= 0; beta1
< j
.states(); beta1
++ )
496 for( size_t beta2
= 0; beta2
< j
.states(); beta2
++ )
497 if( beta2
!= beta1
) {
498 size_t as
= 1, bs
= 1;
503 T f1
= slice( ij
, alpha1
* as
+ beta1
* bs
).p().divide( slice( ij
, alpha2
* as
+ beta1
* bs
).p() ).max();
504 T f2
= slice( ij
, alpha2
* as
+ beta2
* bs
).p().divide( slice( ij
, alpha1
* as
+ beta2
* bs
).p() ).max();
510 return std::tanh( 0.25 * std::log( max
) );
514 /// Writes a TFactor to an output stream
517 template<typename T
> std::ostream
& operator<< (std::ostream
& os
, const TFactor
<T
>& P
) {
518 os
<< "(" << P
.vars() << ", (";
519 for( size_t i
= 0; i
< P
.states(); i
++ )
520 os
<< (i
== 0 ? "" : ", ") << P
[i
];
526 /// Returns distance between two TFactors f and g, according to the distance measure dt
528 * \pre f.vars() == g.vars()
530 template<typename T
> Real
dist( const TFactor
<T
> &f
, const TFactor
<T
> &g
, Prob::DistType dt
) {
531 if( f
.vars().empty() || g
.vars().empty() )
534 DAI_DEBASSERT( f
.vars() == g
.vars() );
535 return dist( f
.p(), g
.p(), dt
);
540 /// Returns the pointwise maximum of two TFactors
542 * \pre f.vars() == g.vars()
544 template<typename T
> TFactor
<T
> max( const TFactor
<T
> &f
, const TFactor
<T
> &g
) {
545 assert( f
._vs
== g
._vs
);
546 return TFactor
<T
>( f
._vs
, min( f
.p(), g
.p() ) );
550 /// Returns the pointwise minimum of two TFactors
552 * \pre f.vars() == g.vars()
554 template<typename T
> TFactor
<T
> min( const TFactor
<T
> &f
, const TFactor
<T
> &g
) {
555 assert( f
._vs
== g
._vs
);
556 return TFactor
<T
>( f
._vs
, max( f
.p(), g
.p() ) );
560 /// Calculates the mutual information between the two variables that f depends on, under the distribution given by f
562 * \pre f.vars().size() == 2
564 template<typename T
> Real
MutualInfo(const TFactor
<T
> &f
) {
565 assert( f
.vars().size() == 2 );
566 VarSet::const_iterator it
= f
.vars().begin();
567 Var i
= *it
; it
++; Var j
= *it
;
568 TFactor
<T
> projection
= f
.marginal(i
) * f
.marginal(j
);
569 return dist( f
.normalized(), projection
, Prob::DISTKL
);
573 /// Represents a factor with values of type Real.
574 typedef TFactor
<Real
> Factor
;
577 } // end of namespace dai