1 /* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
2 Copyright (C) 2002 Martijn Leisink [martijn@mbfys.kun.nl]
3 Radboud University Nijmegen, The Netherlands
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_factor_h
24 #define __defined_libdai_factor_h
30 #include <dai/varset.h>
31 #include <dai/index.h>
37 template<typename T
> class TFactor
;
38 typedef TFactor
<Real
> Factor
;
42 template<typename T
> Real
dist( const TFactor
<T
> & x
, const TFactor
<T
> & y
, Prob::DistType dt
);
43 template<typename T
> Real
KL_dist( const TFactor
<T
> & p
, const TFactor
<T
> & q
);
44 template<typename T
> Real
MutualInfo( const TFactor
<T
> & p
);
45 template<typename T
> TFactor
<T
> max( const TFactor
<T
> & P
, const TFactor
<T
> & Q
);
46 template<typename T
> TFactor
<T
> min( const TFactor
<T
> & P
, const TFactor
<T
> & Q
);
47 template<typename T
> std::ostream
& operator<< (std::ostream
& os
, const TFactor
<T
>& P
);
50 // T should be castable from and to double
51 template <typename T
> class TFactor
{
57 // Construct Factor with empty VarSet but nonempty _p
58 TFactor ( Real p
= 1.0 ) : _vs(), _p(1,p
) {}
60 // Construct Factor from VarSet
61 TFactor( const VarSet
& ns
) : _vs(ns
), _p(nrStates(_vs
)) {}
63 // Construct Factor from VarSet and initial value
64 TFactor( const VarSet
& ns
, Real p
) : _vs(ns
), _p(nrStates(_vs
),p
) {}
66 // Construct Factor from VarSet and initial array
67 TFactor( const VarSet
& ns
, const Real
* p
) : _vs(ns
), _p(nrStates(_vs
),p
) {}
69 // Construct Factor from VarSet and TProb<T>
70 TFactor( const VarSet
& ns
, const TProb
<T
>& p
) : _vs(ns
), _p(p
) {
72 assert( nrStates(_vs
) == _p
.size() );
76 // Construct Factor from Var
77 TFactor( const Var
& n
) : _vs(n
), _p(n
.states()) {}
80 TFactor( const TFactor
<T
> &x
) : _vs(x
._vs
), _p(x
._p
) {}
82 // Assignment operator
83 TFactor
<T
> & operator= (const TFactor
<T
> &x
) {
91 const TProb
<T
> & p() const { return _p
; }
92 TProb
<T
> & p() { return _p
; }
93 const VarSet
& vars() const { return _vs
; }
94 size_t states() const { return _p
.size(); }
96 T
operator[] (size_t i
) const { return _p
[i
]; }
97 T
& operator[] (size_t i
) { return _p
[i
]; }
98 TFactor
<T
> & fill (T p
)
99 { _p
.fill( p
); return(*this); }
100 TFactor
<T
> & randomize ()
101 { _p
.randomize(); return(*this); }
102 TFactor
<T
> operator* (T x
) const {
103 Factor result
= *this;
107 TFactor
<T
>& operator*= (T x
) {
111 TFactor
<T
> operator/ (T x
) const {
112 Factor result
= *this;
116 TFactor
<T
>& operator/= (T x
) {
120 TFactor
<T
> operator* (const TFactor
<T
>& Q
) const;
121 TFactor
<T
> operator/ (const TFactor
<T
>& Q
) const;
122 TFactor
<T
>& operator*= (const TFactor
<T
>& Q
) { return( *this = (*this * Q
) ); }
123 TFactor
<T
>& operator/= (const TFactor
<T
>& Q
) { return( *this = (*this / Q
) ); }
124 TFactor
<T
> operator+ (const TFactor
<T
>& Q
) const {
126 assert( Q
._vs
== _vs
);
128 TFactor
<T
> sum(*this);
132 TFactor
<T
> operator- (const TFactor
<T
>& Q
) const {
134 assert( Q
._vs
== _vs
);
136 TFactor
<T
> sum(*this);
140 TFactor
<T
>& operator+= (const TFactor
<T
>& Q
) {
142 assert( Q
._vs
== _vs
);
147 TFactor
<T
>& operator-= (const TFactor
<T
>& Q
) {
149 assert( Q
._vs
== _vs
);
154 TFactor
<T
>& operator+= (T q
) {
158 TFactor
<T
>& operator-= (T q
) {
162 TFactor
<T
> operator+ (T q
) const {
163 TFactor
<T
> result(*this);
167 TFactor
<T
> operator- (T q
) const {
168 TFactor
<T
> result(*this);
173 TFactor
<T
> operator^ (Real a
) const { TFactor
<T
> x
; x
._vs
= _vs
; x
._p
= _p
^a
; return x
; }
174 TFactor
<T
>& operator^= (Real a
) { _p
^= a
; return *this; }
176 TFactor
<T
>& makeZero( Real epsilon
) {
177 _p
.makeZero( epsilon
);
181 TFactor
<T
>& makePositive( Real epsilon
) {
182 _p
.makePositive( epsilon
);
186 TFactor
<T
> inverse() const {
189 inv
._p
= _p
.inverse(true); // FIXME
193 TFactor
<T
> divided_by( const TFactor
<T
>& denom
) const {
195 assert( denom
._vs
== _vs
);
197 TFactor
<T
> quot(*this);
202 TFactor
<T
>& divide( const TFactor
<T
>& denom
) {
204 assert( denom
._vs
== _vs
);
210 TFactor
<T
> exp() const {
217 TFactor
<T
> abs() const {
224 TFactor
<T
> log() const {
231 TFactor
<T
> log0() const {
238 T
normalize( typename
Prob::NormType norm
= Prob::NORMPROB
) { return _p
.normalize( norm
); }
239 TFactor
<T
> normalized( typename
Prob::NormType norm
= Prob::NORMPROB
) const {
242 result
._p
= _p
.normalized( norm
);
246 // returns slice of this factor where the subset ns is in state ns_state
247 Factor
slice( const VarSet
& ns
, size_t ns_state
) const {
249 VarSet nsrem
= _vs
/ ns
;
250 Factor
result( nsrem
, 0.0 );
253 IndexFor
i_ns (ns
, _vs
);
254 IndexFor
i_nsrem (nsrem
, _vs
);
255 for( size_t i
= 0; i
< states(); i
++, ++i_ns
, ++i_nsrem
)
256 if( (size_t)i_ns
== ns_state
)
257 result
._p
[i_nsrem
] = _p
[i
];
262 // returns unnormalized marginal; ns should be a subset of vars()
263 TFactor
<T
> partSum(const VarSet
& ns
) const;
264 // returns (normalized by default) marginal; ns should be a subset of vars()
265 TFactor
<T
> marginal(const VarSet
& ns
, bool normed
= true) const { if(normed
) return partSum(ns
).normalized(); else return partSum(ns
); }
266 // sums out all variables except those in ns
267 TFactor
<T
> notSum(const VarSet
& ns
) const { return partSum(vars() ^ ns
); }
269 // embeds this factor in larger varset ns
270 TFactor
<T
> embed(const VarSet
& ns
) const {
276 return (*this) * Factor(ns
/ vs
, 1.0);
279 bool hasNaNs() const { return _p
.hasNaNs(); }
280 bool hasNegatives() const { return _p
.hasNegatives(); }
281 T
totalSum() const { return _p
.totalSum(); }
282 T
maxAbs() const { return _p
.maxAbs(); }
283 T
maxVal() const { return _p
.maxVal(); }
284 T
minVal() const { return _p
.minVal(); }
285 Real
entropy() const { return _p
.entropy(); }
286 T
strength( const Var
&i
, const Var
&j
) const;
288 friend Real
dist( const TFactor
<T
> & x
, const TFactor
<T
> & y
, Prob::DistType dt
) {
289 if( x
._vs
.empty() || y
._vs
.empty() )
293 assert( x
._vs
== y
._vs
);
295 return dist( x
._p
, y
._p
, dt
);
298 friend Real KL_dist
<> (const TFactor
<T
> & p
, const TFactor
<T
> & q
);
299 friend Real MutualInfo
<> ( const TFactor
<T
> & P
);
300 template<class U
> friend std::ostream
& operator<< (std::ostream
& os
, const TFactor
<U
>& P
);
304 template<typename T
> TFactor
<T
> TFactor
<T
>::partSum(const VarSet
& ns
) const {
309 TFactor
<T
> res( ns
, 0.0 );
311 IndexFor
i_res( ns
, _vs
);
312 for( size_t i
= 0; i
< _p
.size(); i
++, ++i_res
)
313 res
._p
[i_res
] += _p
[i
];
319 template<typename T
> std::ostream
& operator<< (std::ostream
& os
, const TFactor
<T
>& P
) {
320 os
<< "(" << P
.vars() << " <";
321 for( size_t i
= 0; i
< P
._p
.size(); i
++ )
322 os
<< P
._p
[i
] << " ";
328 template<typename T
> TFactor
<T
> TFactor
<T
>::operator* (const TFactor
<T
>& Q
) const {
329 TFactor
<T
> prod( _vs
| Q
._vs
, 0.0 );
331 IndexFor
i1(_vs
, prod
._vs
);
332 IndexFor
i2(Q
._vs
, prod
._vs
);
334 for( size_t i
= 0; i
< prod
._p
.size(); i
++, ++i1
, ++i2
)
335 prod
._p
[i
] += _p
[i1
] * Q
._p
[i2
];
341 template<typename T
> TFactor
<T
> TFactor
<T
>::operator/ (const TFactor
<T
>& Q
) const {
342 TFactor
<T
> quot( _vs
+ Q
._vs
, 0.0 );
344 IndexFor
i1(_vs
, quot
._vs
);
345 IndexFor
i2(Q
._vs
, quot
._vs
);
347 for( size_t i
= 0; i
< quot
._p
.size(); i
++, ++i1
, ++i2
)
348 quot
._p
[i
] += _p
[i1
] / Q
._p
[i2
];
354 template<typename T
> Real
KL_dist(const TFactor
<T
> & P
, const TFactor
<T
> & Q
) {
355 if( P
._vs
.empty() || Q
._vs
.empty() )
359 assert( P
._vs
== Q
._vs
);
361 return KL_dist( P
._p
, Q
._p
);
366 // calculate mutual information of x_i and x_j where P.vars() = \{x_i,x_j\}
367 template<typename T
> Real
MutualInfo(const TFactor
<T
> & P
) {
368 assert( P
._vs
.size() == 2 );
369 VarSet::const_iterator it
= P
._vs
.begin();
370 Var i
= *it
; it
++; Var j
= *it
;
371 TFactor
<T
> projection
= P
.marginal(i
) * P
.marginal(j
);
372 return real( KL_dist( P
.normalized(), projection
) );
376 template<typename T
> TFactor
<T
> max( const TFactor
<T
> & P
, const TFactor
<T
> & Q
) {
377 assert( P
._vs
== Q
._vs
);
378 return TFactor
<T
>( P
._vs
, min( P
.p(), Q
.p() ) );
381 template<typename T
> TFactor
<T
> min( const TFactor
<T
> & P
, const TFactor
<T
> & Q
) {
382 assert( P
._vs
== Q
._vs
);
383 return TFactor
<T
>( P
._vs
, max( P
.p(), Q
.p() ) );
386 // calculate N(psi, i, j)
387 template<typename T
> T TFactor
<T
>::strength( const Var
&i
, const Var
&j
) const {
389 assert( _vs
.contains( i
) );
390 assert( _vs
.contains( j
) );
396 for( size_t alpha1
= 0; alpha1
< i
.states(); alpha1
++ )
397 for( size_t alpha2
= 0; alpha2
< i
.states(); alpha2
++ )
398 if( alpha2
!= alpha1
)
399 for( size_t beta1
= 0; beta1
< j
.states(); beta1
++ )
400 for( size_t beta2
= 0; beta2
< j
.states(); beta2
++ )
401 if( beta2
!= beta1
) {
402 size_t as
= 1, bs
= 1;
407 T f1
= slice( ij
, alpha1
* as
+ beta1
* bs
).p().divide( slice( ij
, alpha2
* as
+ beta1
* bs
).p() ).maxVal();
408 T f2
= slice( ij
, alpha2
* as
+ beta2
* bs
).p().divide( slice( ij
, alpha1
* as
+ beta2
* bs
).p() ).maxVal();
414 return std::tanh( 0.25 * std::log( max
) );
418 template<typename T
> TFactor
<T
> RemoveFirstOrderInteractions( const TFactor
<T
> & psi
) {
419 TFactor
<T
> result
= psi
;
421 VarSet vars
= psi
.vars();
422 for( size_t iter
= 0; iter
< 100; iter
++ ) {
423 for( VarSet::const_iterator n
= vars
.begin(); n
!= vars
.end(); n
++ )
424 result
= result
* result
.partSum(*n
).inverse();
432 } // end of namespace dai