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
;
39 typedef TFactor
<Complex
> CFactor
;
43 template<typename T
> Real
dist( const TFactor
<T
> & x
, const TFactor
<T
> & y
, Prob::DistType dt
);
44 template<typename T
> Complex
KL_dist( const TFactor
<T
> & p
, const TFactor
<T
> & q
);
45 template<typename T
> std::ostream
& operator<< (std::ostream
& os
, const TFactor
<T
>& P
);
48 // T should be castable from and to double and to complex
49 template <typename T
> class TFactor
{
55 // Default constructor
56 TFactor () : _vs(), _p(1,1.0) {}
58 // Construct Factor from VarSet
59 TFactor( const VarSet
& ns
) : _vs(ns
), _p(_vs
.states()) {}
61 // Construct Factor from VarSet and initial value
62 TFactor( const VarSet
& ns
, Real p
) : _vs(ns
), _p(_vs
.states(),p
) {}
64 // Construct Factor from VarSet and initial array
65 TFactor( const VarSet
& ns
, const Real
* p
) : _vs(ns
), _p(_vs
.states(),p
) {}
67 // Construct Factor from VarSet and TProb<T>
68 TFactor( const VarSet
& ns
, const TProb
<T
> p
) : _vs(ns
), _p(p
) {
70 assert( _vs
.states() == _p
.size() );
74 // Construct Factor from Var
75 TFactor( const Var
& n
) : _vs(n
), _p(n
.states()) {}
78 TFactor( const TFactor
<T
> &x
) : _vs(x
._vs
), _p(x
._p
) {}
80 // Assignment operator
81 TFactor
<T
> & operator= (const TFactor
<T
> &x
) {
89 const TProb
<T
> & p() const { return _p
; }
90 TProb
<T
> & p() { return _p
; }
91 const VarSet
& vars() const { return _vs
; }
92 size_t states() const {
94 assert( _vs
.states() == _p
.size() );
99 T
operator[] (size_t i
) const { return _p
[i
]; }
100 T
& operator[] (size_t i
) { return _p
[i
]; }
101 TFactor
<T
> & fill (T p
)
102 { _p
.fill( p
); return(*this); }
103 TFactor
<T
> & randomize ()
104 { _p
.randomize(); return(*this); }
105 TFactor
<T
> operator* (T x
) const {
106 Factor result
= *this;
110 TFactor
<T
>& operator*= (T x
) {
114 TFactor
<T
> operator/ (T x
) const {
115 Factor result
= *this;
119 TFactor
<T
>& operator/= (T x
) {
123 TFactor
<T
> operator* (const TFactor
<T
>& Q
) const;
124 TFactor
<T
>& operator*= (const TFactor
<T
>& Q
) { return( *this = (*this * Q
) ); }
125 TFactor
<T
> operator+ (const TFactor
<T
>& Q
) const {
127 assert( Q
._vs
== _vs
);
129 TFactor
<T
> sum(*this);
133 TFactor
<T
> operator- (const TFactor
<T
>& Q
) const {
135 assert( Q
._vs
== _vs
);
137 TFactor
<T
> sum(*this);
141 TFactor
<T
>& operator+= (const TFactor
<T
>& Q
) {
143 assert( Q
._vs
== _vs
);
148 TFactor
<T
>& operator-= (const TFactor
<T
>& Q
) {
150 assert( Q
._vs
== _vs
);
156 TFactor
<T
> operator^ (Real a
) const { TFactor
<T
> x
; x
._vs
= _vs
; x
._p
= _p
^a
; return x
; }
157 TFactor
<T
>& operator^= (Real a
) { _p
^= a
; return *this; }
159 TFactor
<T
>& makeZero( Real epsilon
) {
160 _p
.makeZero( epsilon
);
164 TFactor
<T
> inverse() const {
167 inv
._p
= _p
.inverse(true); // FIXME
171 TFactor
<T
> divided_by( const TFactor
<T
>& denom
) const {
173 assert( denom
._vs
== _vs
);
175 TFactor
<T
> quot(*this);
180 TFactor
<T
>& divide( const TFactor
<T
>& denom
) {
182 assert( denom
._vs
== _vs
);
188 TFactor
<T
> exp() const {
195 TFactor
<T
> log() const {
202 TFactor
<T
> log0() const {
209 CFactor
clog0() const {
216 T
normalize( typename
Prob::NormType norm
) { return _p
.normalize( norm
); }
217 TFactor
<T
> normalized( typename
Prob::NormType norm
) const {
220 result
._p
= _p
.normalized( norm
);
224 // returns slice of this factor where the subset ns is in state ns_state
225 Factor
slice( const VarSet
& ns
, size_t ns_state
) const {
227 VarSet nsrem
= _vs
/ ns
;
228 Factor
result( nsrem
, 0.0 );
231 IndexFor
i_ns (ns
, _vs
);
232 IndexFor
i_nsrem (nsrem
, _vs
);
233 for( size_t i
= 0; i
< states(); i
++, ++i_ns
, ++i_nsrem
)
234 if( (size_t)i_ns
== ns_state
)
235 result
._p
[i_nsrem
] = _p
[i
];
240 // returns unnormalized marginal
241 TFactor
<T
> part_sum(const VarSet
& ns
) const;
242 // returns normalized marginal
243 TFactor
<T
> marginal(const VarSet
& ns
) const { return part_sum(ns
).normalized( Prob::NORMPROB
); }
245 bool hasNaNs() const { return _p
.hasNaNs(); }
246 bool hasNegatives() const { return _p
.hasNegatives(); }
247 T
totalSum() const { return _p
.totalSum(); }
248 T
maxAbs() const { return _p
.maxAbs(); }
249 T
max() const { return _p
.max(); }
250 Complex
entropy() const { return _p
.entropy(); }
251 T
strength( const Var
&i
, const Var
&j
) const;
253 friend Real
dist( const TFactor
<T
> & x
, const TFactor
<T
> & y
, Prob::DistType dt
) {
254 if( x
._vs
.empty() || y
._vs
.empty() )
258 assert( x
._vs
== y
._vs
);
260 return dist( x
._p
, y
._p
, dt
);
263 friend Complex KL_dist
<> (const TFactor
<T
> & p
, const TFactor
<T
> & q
);
264 template<class U
> friend std::ostream
& operator<< (std::ostream
& os
, const TFactor
<U
>& P
);
268 template<typename T
> TFactor
<T
> TFactor
<T
>::part_sum(const VarSet
& ns
) const {
273 TFactor
<T
> res( ns
, 0.0 );
275 IndexFor
i_res( ns
, _vs
);
276 for( size_t i
= 0; i
< _p
.size(); i
++, ++i_res
)
277 res
._p
[i_res
] += _p
[i
];
283 template<typename T
> std::ostream
& operator<< (std::ostream
& os
, const TFactor
<T
>& P
) {
284 os
<< "(" << P
.vars() << " <";
285 for( size_t i
= 0; i
< P
._p
.size(); i
++ )
286 os
<< P
._p
[i
] << " ";
292 template<typename T
> TFactor
<T
> TFactor
<T
>::operator* (const TFactor
<T
>& Q
) const {
293 TFactor
<T
> prod( _vs
| Q
._vs
, 0.0 );
295 IndexFor
i1(_vs
, prod
._vs
);
296 IndexFor
i2(Q
._vs
, prod
._vs
);
298 for( size_t i
= 0; i
< prod
._p
.size(); i
++, ++i1
, ++i2
)
299 prod
._p
[i
] += _p
[i1
] * Q
._p
[i2
];
305 template<typename T
> Complex
KL_dist(const TFactor
<T
> & P
, const TFactor
<T
> & Q
) {
306 if( P
._vs
.empty() || Q
._vs
.empty() )
310 assert( P
._vs
== Q
._vs
);
312 return KL_dist( P
._p
, Q
._p
);
317 // calculate N(psi, i, j)
318 template<typename T
> T TFactor
<T
>::strength( const Var
&i
, const Var
&j
) const {
327 for( size_t alpha1
= 0; alpha1
< i
.states(); alpha1
++ )
328 for( size_t alpha2
= 0; alpha2
< i
.states(); alpha2
++ )
329 if( alpha2
!= alpha1
)
330 for( size_t beta1
= 0; beta1
< j
.states(); beta1
++ )
331 for( size_t beta2
= 0; beta2
< j
.states(); beta2
++ )
332 if( beta2
!= beta1
) {
333 size_t as
= 1, bs
= 1;
338 T f1
= slice( ij
, alpha1
* as
+ beta1
* bs
).p().divide( slice( ij
, alpha2
* as
+ beta1
* bs
).p() ).max();
339 T f2
= slice( ij
, alpha2
* as
+ beta2
* bs
).p().divide( slice( ij
, alpha1
* as
+ beta2
* bs
).p() ).max();
345 return std::tanh( 0.25 * std::log( max
) );
349 template<typename T
> TFactor
<T
> RemoveFirstOrderInteractions( const TFactor
<T
> & psi
) {
350 TFactor
<T
> result
= psi
;
352 VarSet vars
= psi
.vars();
353 for( size_t iter
= 0; iter
< 100; iter
++ ) {
354 for( VarSet::const_iterator n
= vars
.begin(); n
!= vars
.end(); n
++ )
355 result
= result
* result
.part_sum(*n
).inverse();
356 result
.normalize( Prob::NORMPROB
);
363 } // end of namespace dai