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_factor_h
23 #define __defined_libdai_factor_h
29 #include <dai/varset.h>
30 #include <dai/index.h>
36 template<typename T
> class TFactor
;
37 typedef TFactor
<Real
> Factor
;
38 typedef TFactor
<Complex
> CFactor
;
42 template<typename T
> Real
dist( const TFactor
<T
> & x
, const TFactor
<T
> & y
, Prob::DistType dt
);
43 template<typename T
> Complex
KL_dist( const TFactor
<T
> & p
, const TFactor
<T
> & q
);
44 template<typename T
> std::ostream
& operator<< (std::ostream
& os
, const TFactor
<T
>& P
);
47 // T should be castable from and to double and to complex
48 template <typename T
> class TFactor
{
54 // Default constructor
55 TFactor () : _vs(), _p(1,1.0) {}
57 // Construct Factor from VarSet
58 TFactor( const VarSet
& ns
) : _vs(ns
), _p(_vs
.stateSpace()) {}
60 // Construct Factor from VarSet and initial value
61 TFactor( const VarSet
& ns
, Real p
) : _vs(ns
), _p(_vs
.stateSpace(),p
) {}
63 // Construct Factor from VarSet and initial array
64 TFactor( const VarSet
& ns
, const Real
* p
) : _vs(ns
), _p(_vs
.stateSpace(),p
) {}
66 // Construct Factor from VarSet and TProb<T>
67 TFactor( const VarSet
& ns
, const TProb
<T
> p
) : _vs(ns
), _p(p
) {
69 assert( _vs
.stateSpace() == _p
.size() );
73 // Construct Factor from Var
74 TFactor( const Var
& n
) : _vs(n
), _p(n
.states()) {}
77 TFactor( const TFactor
<T
> &x
) : _vs(x
._vs
), _p(x
._p
) {}
79 // Assignment operator
80 TFactor
<T
> & operator= (const TFactor
<T
> &x
) {
88 const TProb
<T
> & p() const { return _p
; }
89 TProb
<T
> & p() { return _p
; }
90 const VarSet
& vars() const { return _vs
; }
91 size_t stateSpace() const {
93 assert( _vs
.stateSpace() == _p
.size() );
98 T
operator[] (size_t i
) const { return _p
[i
]; }
99 T
& operator[] (size_t i
) { return _p
[i
]; }
100 TFactor
<T
> & fill (T p
)
101 { _p
.fill( p
); return(*this); }
102 TFactor
<T
> & randomize ()
103 { _p
.randomize(); return(*this); }
104 TFactor
<T
> operator* (T x
) const {
105 Factor result
= *this;
109 TFactor
<T
>& operator*= (T x
) {
113 TFactor
<T
> operator/ (T x
) const {
114 Factor result
= *this;
118 TFactor
<T
>& operator/= (T x
) {
122 TFactor
<T
> operator* (const TFactor
<T
>& Q
) const;
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
);
155 TFactor
<T
> operator^ (Real a
) const { TFactor
<T
> x
; x
._vs
= _vs
; x
._p
= _p
^a
; return x
; }
156 TFactor
<T
>& operator^= (Real a
) { _p
^= a
; return *this; }
158 TFactor
<T
>& makeZero( Real epsilon
) {
159 _p
.makeZero( epsilon
);
163 TFactor
<T
> inverse() const {
166 inv
._p
= _p
.inverse(true); // FIXME
170 TFactor
<T
> divided_by( const TFactor
<T
>& denom
) const {
172 assert( denom
._vs
== _vs
);
174 TFactor
<T
> quot(*this);
179 TFactor
<T
>& divide( const TFactor
<T
>& denom
) {
181 assert( denom
._vs
== _vs
);
187 TFactor
<T
> exp() const {
194 TFactor
<T
> log() const {
201 TFactor
<T
> log0() const {
208 CFactor
clog0() const {
215 T
normalize( typename
Prob::NormType norm
) { return _p
.normalize( norm
); }
216 TFactor
<T
> normalized( typename
Prob::NormType norm
) const {
219 result
._p
= _p
.normalized( norm
);
223 // returns slice of this factor where the subset ns is in state ns_state
224 Factor
slice( const VarSet
& ns
, size_t ns_state
) const {
226 VarSet nsrem
= _vs
/ ns
;
227 Factor
result( nsrem
, 0.0 );
230 Index
i_ns (ns
, _vs
);
231 Index
i_nsrem (nsrem
, _vs
);
232 for( size_t i
= 0; i
< stateSpace(); i
++, ++i_ns
, ++i_nsrem
)
233 if( (size_t)i_ns
== ns_state
)
234 result
._p
[i_nsrem
] = _p
[i
];
239 // returns unnormalized marginal
240 TFactor
<T
> part_sum(const VarSet
& ns
) const;
241 // returns normalized marginal
242 TFactor
<T
> marginal(const VarSet
& ns
) const { return part_sum(ns
).normalized( Prob::NORMPROB
); }
244 bool hasNaNs() const { return _p
.hasNaNs(); }
245 bool hasNegatives() const { return _p
.hasNegatives(); }
246 T
totalSum() const { return _p
.totalSum(); }
247 T
maxAbs() const { return _p
.maxAbs(); }
248 T
max() const { return _p
.max(); }
249 Complex
entropy() const { return _p
.entropy(); }
250 T
strength( const Var
&i
, const Var
&j
) const;
252 friend Real
dist( const TFactor
<T
> & x
, const TFactor
<T
> & y
, Prob::DistType dt
) {
253 if( x
._vs
.empty() || y
._vs
.empty() )
257 assert( x
._vs
== y
._vs
);
259 return dist( x
._p
, y
._p
, dt
);
262 friend Complex KL_dist
<> (const TFactor
<T
> & p
, const TFactor
<T
> & q
);
263 template<class U
> friend std::ostream
& operator<< (std::ostream
& os
, const TFactor
<U
>& P
);
267 template<typename T
> TFactor
<T
> TFactor
<T
>::part_sum(const VarSet
& ns
) const {
272 TFactor
<T
> res( ns
, 0.0 );
274 Index
i_res( ns
, _vs
);
275 for( size_t i
= 0; i
< _p
.size(); i
++, ++i_res
)
276 res
._p
[i_res
] += _p
[i
];
282 template<typename T
> std::ostream
& operator<< (std::ostream
& os
, const TFactor
<T
>& P
) {
283 os
<< "(" << P
.vars() << " <";
284 for( size_t i
= 0; i
< P
._p
.size(); i
++ )
285 os
<< P
._p
[i
] << " ";
291 template<typename T
> TFactor
<T
> TFactor
<T
>::operator* (const TFactor
<T
>& Q
) const {
292 TFactor
<T
> prod( _vs
| Q
._vs
, 0.0 );
294 Index
i1(_vs
, prod
._vs
);
295 Index
i2(Q
._vs
, prod
._vs
);
297 for( size_t i
= 0; i
< prod
._p
.size(); i
++, ++i1
, ++i2
)
298 prod
._p
[i
] += _p
[i1
] * Q
._p
[i2
];
304 template<typename T
> Complex
KL_dist(const TFactor
<T
> & P
, const TFactor
<T
> & Q
) {
305 if( P
._vs
.empty() || Q
._vs
.empty() )
309 assert( P
._vs
== Q
._vs
);
311 return KL_dist( P
._p
, Q
._p
);
316 // calculate N(psi, i, j)
317 template<typename T
> T TFactor
<T
>::strength( const Var
&i
, const Var
&j
) const {
326 for( size_t alpha1
= 0; alpha1
< i
.states(); alpha1
++ )
327 for( size_t alpha2
= 0; alpha2
< i
.states(); alpha2
++ )
328 if( alpha2
!= alpha1
)
329 for( size_t beta1
= 0; beta1
< j
.states(); beta1
++ )
330 for( size_t beta2
= 0; beta2
< j
.states(); beta2
++ )
331 if( beta2
!= beta1
) {
332 size_t as
= 1, bs
= 1;
337 T f1
= slice( ij
, alpha1
* as
+ beta1
* bs
).p().divide( slice( ij
, alpha2
* as
+ beta1
* bs
).p() ).max();
338 T f2
= slice( ij
, alpha2
* as
+ beta2
* bs
).p().divide( slice( ij
, alpha1
* as
+ beta2
* bs
).p() ).max();
344 return std::tanh( 0.25 * std::log( max
) );
348 template<typename T
> TFactor
<T
> RemoveFirstOrderInteractions( const TFactor
<T
> & psi
) {
349 TFactor
<T
> result
= psi
;
351 VarSet vars
= psi
.vars();
352 for( size_t iter
= 0; iter
< 100; iter
++ ) {
353 for( VarSet::const_iterator n
= vars
.begin(); n
!= vars
.end(); n
++ )
354 result
= result
* result
.part_sum(*n
).inverse();
355 result
.normalize( Prob::NORMPROB
);
362 } // end of namespace dai