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
25 #include <dai/diffs.h>
34 const char *HAK::Name
= "HAK";
37 bool HAK::checkProperties() {
38 if( !HasProperty("tol") )
40 if (!HasProperty("maxiter") )
42 if (!HasProperty("verbose") )
44 if( !HasProperty("doubleloop") )
46 if( !HasProperty("clusters") )
49 ConvertPropertyTo
<double>("tol");
50 ConvertPropertyTo
<size_t>("maxiter");
51 ConvertPropertyTo
<size_t>("verbose");
52 ConvertPropertyTo
<bool>("doubleloop");
53 ConvertPropertyTo
<ClustersType
>("clusters");
55 if( HasProperty("loopdepth") )
56 ConvertPropertyTo
<size_t>("loopdepth");
57 else if( Clusters() == ClustersType::LOOP
)
64 void HAK::constructMessages() {
65 // Create outer beliefs
67 _Qa
.reserve(nr_ORs());
68 for( size_t alpha
= 0; alpha
< nr_ORs(); alpha
++ )
69 _Qa
.push_back( Factor( OR(alpha
).vars() ) );
71 // Create inner beliefs
73 _Qb
.reserve(nr_IRs());
74 for( size_t beta
= 0; beta
< nr_IRs(); beta
++ )
75 _Qb
.push_back( Factor( IR(beta
) ) );
79 _muab
.reserve(nr_Redges());
81 _muba
.reserve(nr_Redges());
82 for( vector
<R_edge_t
>::const_iterator ab
= Redges().begin(); ab
!= Redges().end(); ab
++ ) {
83 _muab
.push_back( Factor( IR(ab
->second
) ) );
84 _muba
.push_back( Factor( IR(ab
->second
) ) );
89 HAK::HAK(const RegionGraph
& rg
, const Properties
&opts
) : DAIAlgRG(rg
, opts
) {
90 assert( checkProperties() );
96 void HAK::findLoopClusters( const FactorGraph
& fg
, set
<VarSet
> &allcl
, VarSet newcl
, const Var
& root
, size_t length
, VarSet vars
) {
97 for( VarSet::const_iterator in
= vars
.begin(); in
!= vars
.end(); in
++ ) {
98 VarSet ind
= fg
.delta( *in
);
99 if( (newcl
.size()) >= 2 && (ind
>> root
) ) {
100 allcl
.insert( newcl
| *in
);
102 else if( length
> 1 )
103 findLoopClusters( fg
, allcl
, newcl
| *in
, root
, length
- 1, ind
/ newcl
);
108 HAK::HAK(const FactorGraph
& fg
, const Properties
&opts
) : DAIAlgRG(opts
) {
109 assert( checkProperties() );
112 if( Clusters() == ClustersType::MIN
) {
114 } else if( Clusters() == ClustersType::DELTA
) {
115 for( size_t i
= 0; i
< fg
.nrVars(); i
++ )
116 cl
.push_back(fg
.Delta(fg
.var(i
)));
117 } else if( Clusters() == ClustersType::LOOP
) {
120 for( vector
<Var
>::const_iterator i0
= fg
.vars().begin(); i0
!= fg
.vars().end(); i0
++ ) {
121 VarSet i0d
= fg
.delta(*i0
);
122 if( LoopDepth() > 1 )
123 findLoopClusters( fg
, scl
, *i0
, *i0
, LoopDepth() - 1, fg
.delta(*i0
) );
125 for( set
<VarSet
>::const_iterator c
= scl
.begin(); c
!= scl
.end(); c
++ )
127 if( Verbose() >= 3 ) {
128 cout
<< "HAK uses the following clusters: " << endl
;
129 for( vector
<VarSet
>::const_iterator cli
= cl
.begin(); cli
!= cl
.end(); cli
++ )
130 cout
<< *cli
<< endl
;
133 throw "Invalid Clusters type";
135 RegionGraph
rg(fg
,cl
);
136 RegionGraph::operator=(rg
);
140 cout
<< "HAK regiongraph: " << *this << endl
;
144 string
HAK::identify() const {
145 stringstream
result (stringstream::out
);
146 result
<< Name
<< GetProperties();
151 void HAK::init( const VarSet
&ns
) {
152 for( vector
<Factor
>::iterator alpha
= _Qa
.begin(); alpha
!= _Qa
.end(); alpha
++ )
153 if( alpha
->vars() && ns
)
154 alpha
->fill( 1.0 / alpha
->stateSpace() );
156 for( size_t beta
= 0; beta
< nr_IRs(); beta
++ )
157 if( IR(beta
) && ns
) {
158 _Qb
[beta
].fill( 1.0 / IR(beta
).stateSpace() );
159 for( R_nb_cit alpha
= nbIR(beta
).begin(); alpha
!= nbIR(beta
).end(); alpha
++ ) {
160 muab(*alpha
,beta
).fill( 1.0 / IR(beta
).stateSpace() );
161 muba(beta
,*alpha
).fill( 1.0 / IR(beta
).stateSpace() );
168 assert( checkProperties() );
170 for( vector
<Factor
>::iterator alpha
= _Qa
.begin(); alpha
!= _Qa
.end(); alpha
++ )
171 alpha
->fill( 1.0 / alpha
->stateSpace() );
173 for( vector
<Factor
>::iterator beta
= _Qb
.begin(); beta
!= _Qb
.end(); beta
++ )
174 beta
->fill( 1.0 / beta
->stateSpace() );
176 for( size_t ab
= 0; ab
< nr_Redges(); ab
++ ) {
177 _muab
[ab
].fill( 1.0 / _muab
[ab
].stateSpace() );
178 _muba
[ab
].fill( 1.0 / _muba
[ab
].stateSpace() );
183 double HAK::doGBP() {
185 cout
<< "Starting " << identify() << "...";
191 // Check whether counting numbers won't lead to problems
192 for( size_t beta
= 0; beta
< nr_IRs(); beta
++ )
193 assert( nbIR(beta
).size() + IR(beta
).c() != 0.0 );
195 // Keep old beliefs to check convergence
196 vector
<Factor
> old_beliefs
;
197 old_beliefs
.reserve( nrVars() );
198 for( size_t i
= 0; i
< nrVars(); i
++ )
199 old_beliefs
.push_back( belief( var(i
) ) );
201 // Differences in single node beliefs
202 Diffs
diffs(nrVars(), 1.0);
205 // do several passes over the network until maximum number of iterations has
206 // been reached or until the maximum belief difference is smaller than tolerance
207 for( iter
= 0; iter
< MaxIter() && diffs
.max() > Tol(); iter
++ ) {
208 for( size_t beta
= 0; beta
< nr_IRs(); beta
++ ) {
209 for( R_nb_cit alpha
= nbIR(beta
).begin(); alpha
!= nbIR(beta
).end(); alpha
++ )
210 muab(*alpha
,beta
) = _Qa
[*alpha
].marginal(IR(beta
)).divided_by( muba(beta
,*alpha
) );
213 for( R_nb_cit alpha
= nbIR(beta
).begin(); alpha
!= nbIR(beta
).end(); alpha
++ )
214 Qb_new
*= muab(*alpha
,beta
) ^ (1 / (nbIR(beta
).size() + IR(beta
).c()));
215 Qb_new
.normalize( _normtype
);
216 if( Qb_new
.hasNaNs() ) {
217 cout
<< "HAK::doGBP: Qb_new has NaNs!" << endl
;
220 // _Qb[beta] = Qb_new.makeZero(1e-100); // damping?
223 for( R_nb_cit alpha
= nbIR(beta
).begin(); alpha
!= nbIR(beta
).end(); alpha
++ ) {
224 muba(beta
,*alpha
) = _Qb
[beta
].divided_by( muab(*alpha
,beta
) );
226 Factor Qa_new
= OR(*alpha
);
227 for( R_nb_cit gamma
= nbOR(*alpha
).begin(); gamma
!= nbOR(*alpha
).end(); gamma
++ )
228 Qa_new
*= muba(*gamma
,*alpha
);
229 Qa_new
^= (1.0 / OR(*alpha
).c());
230 Qa_new
.normalize( _normtype
);
231 if( Qa_new
.hasNaNs() ) {
232 cout
<< "HAK::doGBP: Qa_new has NaNs!" << endl
;
235 // _Qa[*alpha] = Qa_new.makeZero(1e-100); // damping?
236 _Qa
[*alpha
] = Qa_new
;
240 // Calculate new single variable beliefs and compare with old ones
241 for( size_t i
= 0; i
< nrVars(); i
++ ) {
242 Factor new_belief
= belief( var( i
) );
243 diffs
.push( dist( new_belief
, old_beliefs
[i
], Prob::DISTLINF
) );
244 old_beliefs
[i
] = new_belief
;
248 cout
<< "HAK::doGBP: maxdiff " << diffs
.max() << " after " << iter
+1 << " passes" << endl
;
251 updateMaxDiff( diffs
.max() );
253 if( Verbose() >= 1 ) {
254 if( diffs
.max() > Tol() ) {
257 cout
<< "HAK::doGBP: WARNING: not converged within " << MaxIter() << " passes (" << toc() - tic
<< " clocks)...final maxdiff:" << diffs
.max() << endl
;
260 cout
<< "HAK::doGBP: ";
261 cout
<< "converged in " << iter
<< " passes (" << toc() - tic
<< " clocks)." << endl
;
269 double HAK::doDoubleLoop() {
271 cout
<< "Starting " << identify() << "...";
277 // Save original outer regions
278 vector
<FRegion
> org_ORs
= ORs();
280 // Save original inner counting numbers and set negative counting numbers to zero
281 vector
<double> org_IR_cs( nr_IRs(), 0.0 );
282 for( size_t beta
= 0; beta
< nr_IRs(); beta
++ ) {
283 org_IR_cs
[beta
] = IR(beta
).c();
284 if( IR(beta
).c() < 0.0 )
288 // Keep old beliefs to check convergence
289 vector
<Factor
> old_beliefs
;
290 old_beliefs
.reserve( nrVars() );
291 for( size_t i
= 0; i
< nrVars(); i
++ )
292 old_beliefs
.push_back( belief( var(i
) ) );
294 // Differences in single node beliefs
295 Diffs
diffs(nrVars(), 1.0);
297 size_t outer_maxiter
= MaxIter();
298 double outer_tol
= Tol();
299 size_t outer_verbose
= Verbose();
300 double org_maxdiff
= MaxDiff();
302 // Set parameters for inner loop
304 Verbose( outer_verbose
? outer_verbose
- 1 : 0 );
306 size_t outer_iter
= 0;
307 for( outer_iter
= 0; outer_iter
< outer_maxiter
&& diffs
.max() > outer_tol
; outer_iter
++ ) {
308 // Calculate new outer regions
309 for( size_t alpha
= 0; alpha
< nr_ORs(); alpha
++ ) {
310 OR(alpha
) = org_ORs
[alpha
];
311 for( R_nb_cit beta
= nbOR(alpha
).begin(); beta
!= nbOR(alpha
).end(); beta
++ )
312 OR(alpha
) *= _Qb
[*beta
] ^ ((IR(*beta
).c() - org_IR_cs
[*beta
]) / nbIR(*beta
).size());
316 if( isnan( doGBP() ) )
319 // Calculate new single variable beliefs and compare with old ones
320 for( size_t i
= 0; i
< nrVars(); ++i
) {
321 Factor new_belief
= belief( var( i
) );
322 diffs
.push( dist( new_belief
, old_beliefs
[i
], Prob::DISTLINF
) );
323 old_beliefs
[i
] = new_belief
;
327 cout
<< "HAK::doDoubleLoop: maxdiff " << diffs
.max() << " after " << outer_iter
+1 << " passes" << endl
;
330 // restore _maxiter, _verbose and _maxdiff
331 MaxIter( outer_maxiter
);
332 Verbose( outer_verbose
);
333 MaxDiff( org_maxdiff
);
335 updateMaxDiff( diffs
.max() );
337 // Restore original outer regions
340 // Restore original inner counting numbers
341 for( size_t beta
= 0; beta
< nr_IRs(); ++beta
)
342 IR(beta
).c() = org_IR_cs
[beta
];
344 if( Verbose() >= 1 ) {
345 if( diffs
.max() > Tol() ) {
348 cout
<< "HAK::doDoubleLoop: WARNING: not converged within " << outer_maxiter
<< " passes (" << toc() - tic
<< " clocks)...final maxdiff:" << diffs
.max() << endl
;
351 cout
<< "HAK::doDoubleLoop: ";
352 cout
<< "converged in " << outer_iter
<< " passes (" << toc() - tic
<< " clocks)." << endl
;
362 return doDoubleLoop();
368 Factor
HAK::belief( const VarSet
&ns
) const {
369 vector
<Factor
>::const_iterator beta
;
370 for( beta
= _Qb
.begin(); beta
!= _Qb
.end(); beta
++ )
371 if( beta
->vars() >> ns
)
373 if( beta
!= _Qb
.end() )
374 return( beta
->marginal(ns
) );
376 vector
<Factor
>::const_iterator alpha
;
377 for( alpha
= _Qa
.begin(); alpha
!= _Qa
.end(); alpha
++ )
378 if( alpha
->vars() >> ns
)
380 assert( alpha
!= _Qa
.end() );
381 return( alpha
->marginal(ns
) );
386 Factor
HAK::belief( const Var
&n
) const {
387 return belief( (VarSet
)n
);
391 vector
<Factor
> HAK::beliefs() const {
392 vector
<Factor
> result
;
393 for( size_t beta
= 0; beta
< nr_IRs(); beta
++ )
394 result
.push_back( Qb(beta
) );
395 for( size_t alpha
= 0; alpha
< nr_ORs(); alpha
++ )
396 result
.push_back( Qa(alpha
) );
401 Complex
HAK::logZ() const {
403 for( size_t beta
= 0; beta
< nr_IRs(); beta
++ )
404 sum
+= Complex(IR(beta
).c()) * Qb(beta
).entropy();
405 for( size_t alpha
= 0; alpha
< nr_ORs(); alpha
++ ) {
406 sum
+= Complex(OR(alpha
).c()) * Qa(alpha
).entropy();
407 sum
+= (OR(alpha
).log0() * Qa(alpha
)).totalSum();
413 } // end of namespace dai