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 void HAK::setProperties( const PropertySet
&opts
) {
38 assert( opts
.hasKey("tol") );
39 assert( opts
.hasKey("maxiter") );
40 assert( opts
.hasKey("verbose") );
41 assert( opts
.hasKey("doubleloop") );
42 assert( opts
.hasKey("clusters") );
44 props
.tol
= opts
.getStringAs
<double>("tol");
45 props
.maxiter
= opts
.getStringAs
<size_t>("maxiter");
46 props
.verbose
= opts
.getStringAs
<size_t>("verbose");
47 props
.doubleloop
= opts
.getStringAs
<bool>("doubleloop");
48 props
.clusters
= opts
.getStringAs
<Properties::ClustersType
>("clusters");
50 if( opts
.hasKey("loopdepth") )
51 props
.loopdepth
= opts
.getStringAs
<size_t>("loopdepth");
53 assert( props
.clusters
!= Properties::ClustersType::LOOP
);
57 PropertySet
HAK::getProperties() const {
59 opts
.Set( "tol", props
.tol
);
60 opts
.Set( "maxiter", props
.maxiter
);
61 opts
.Set( "verbose", props
.verbose
);
62 opts
.Set( "doubleloop", props
.doubleloop
);
63 opts
.Set( "clusters", props
.clusters
);
64 opts
.Set( "loopdepth", props
.loopdepth
);
69 void HAK::constructMessages() {
70 // Create outer beliefs
73 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ )
74 _Qa
.push_back( Factor( OR(alpha
).vars() ) );
76 // Create inner beliefs
79 for( size_t beta
= 0; beta
< nrIRs(); beta
++ )
80 _Qb
.push_back( Factor( IR(beta
) ) );
84 _muab
.reserve( nrORs() );
86 _muba
.reserve( nrORs() );
87 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ ) {
88 _muab
.push_back( vector
<Factor
>() );
89 _muba
.push_back( vector
<Factor
>() );
90 _muab
[alpha
].reserve( nbOR(alpha
).size() );
91 _muba
[alpha
].reserve( nbOR(alpha
).size() );
92 foreach( const Neighbor
&beta
, nbOR(alpha
) ) {
93 _muab
[alpha
].push_back( Factor( IR(beta
) ) );
94 _muba
[alpha
].push_back( Factor( IR(beta
) ) );
100 HAK::HAK(const RegionGraph
& rg
, const PropertySet
&opts
) : DAIAlgRG(rg
) {
101 setProperties( opts
);
107 void HAK::findLoopClusters( const FactorGraph
& fg
, std::set
<VarSet
> &allcl
, VarSet newcl
, const Var
& root
, size_t length
, VarSet vars
) {
108 for( VarSet::const_iterator in
= vars
.begin(); in
!= vars
.end(); in
++ ) {
109 VarSet ind
= fg
.delta( fg
.findVar( *in
) );
110 if( (newcl
.size()) >= 2 && (ind
>> root
) ) {
111 allcl
.insert( newcl
| *in
);
113 else if( length
> 1 )
114 findLoopClusters( fg
, allcl
, newcl
| *in
, root
, length
- 1, ind
/ newcl
);
119 HAK::HAK(const FactorGraph
& fg
, const PropertySet
&opts
) : DAIAlgRG(), props(), maxdiff(0.0) {
120 setProperties( opts
);
123 if( props
.clusters
== Properties::ClustersType::MIN
) {
125 } else if( props
.clusters
== Properties::ClustersType::DELTA
) {
126 for( size_t i
= 0; i
< fg
.nrVars(); i
++ )
127 cl
.push_back(fg
.Delta(i
));
128 } else if( props
.clusters
== Properties::ClustersType::LOOP
) {
131 for( size_t i0
= 0; i0
< fg
.nrVars(); i0
++ ) {
132 VarSet i0d
= fg
.delta(i0
);
133 if( props
.loopdepth
> 1 )
134 findLoopClusters( fg
, scl
, fg
.var(i0
), fg
.var(i0
), props
.loopdepth
- 1, fg
.delta(i0
) );
136 for( set
<VarSet
>::const_iterator c
= scl
.begin(); c
!= scl
.end(); c
++ )
138 if( props
.verbose
>= 3 ) {
139 cout
<< "HAK uses the following clusters: " << endl
;
140 for( vector
<VarSet
>::const_iterator cli
= cl
.begin(); cli
!= cl
.end(); cli
++ )
141 cout
<< *cli
<< endl
;
144 throw "Invalid Clusters type";
146 RegionGraph
rg(fg
,cl
);
147 RegionGraph::operator=(rg
);
150 if( props
.verbose
>= 3 )
151 cout
<< "HAK regiongraph: " << *this << endl
;
155 string
HAK::identify() const {
156 stringstream
result (stringstream::out
);
157 result
<< Name
<< getProperties();
162 void HAK::init( const VarSet
&ns
) {
163 for( vector
<Factor
>::iterator alpha
= _Qa
.begin(); alpha
!= _Qa
.end(); alpha
++ )
164 if( alpha
->vars().intersects( ns
) )
165 alpha
->fill( 1.0 / alpha
->states() );
167 for( size_t beta
= 0; beta
< nrIRs(); beta
++ )
168 if( IR(beta
).intersects( ns
) ) {
169 _Qb
[beta
].fill( 1.0 );
170 foreach( const Neighbor
&alpha
, nbIR(beta
) ) {
171 size_t _beta
= alpha
.dual
;
172 muab( alpha
, _beta
).fill( 1.0 / IR(beta
).states() );
173 muba( alpha
, _beta
).fill( 1.0 / IR(beta
).states() );
180 for( vector
<Factor
>::iterator alpha
= _Qa
.begin(); alpha
!= _Qa
.end(); alpha
++ )
181 alpha
->fill( 1.0 / alpha
->states() );
183 for( vector
<Factor
>::iterator beta
= _Qb
.begin(); beta
!= _Qb
.end(); beta
++ )
184 beta
->fill( 1.0 / beta
->states() );
186 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ )
187 foreach( const Neighbor
&beta
, nbOR(alpha
) ) {
188 size_t _beta
= beta
.iter
;
189 muab( alpha
, _beta
).fill( 1.0 / muab( alpha
, _beta
).states() );
190 muba( alpha
, _beta
).fill( 1.0 / muab( alpha
, _beta
).states() );
195 double HAK::doGBP() {
196 if( props
.verbose
>= 1 )
197 cout
<< "Starting " << identify() << "...";
198 if( props
.verbose
>= 3)
203 // Check whether counting numbers won't lead to problems
204 for( size_t beta
= 0; beta
< nrIRs(); beta
++ )
205 assert( nbIR(beta
).size() + IR(beta
).c() != 0.0 );
207 // Keep old beliefs to check convergence
208 vector
<Factor
> old_beliefs
;
209 old_beliefs
.reserve( nrVars() );
210 for( size_t i
= 0; i
< nrVars(); i
++ )
211 old_beliefs
.push_back( belief( var(i
) ) );
213 // Differences in single node beliefs
214 Diffs
diffs(nrVars(), 1.0);
217 // do several passes over the network until maximum number of iterations has
218 // been reached or until the maximum belief difference is smaller than tolerance
219 for( iter
= 0; iter
< props
.maxiter
&& diffs
.maxDiff() > props
.tol
; iter
++ ) {
220 for( size_t beta
= 0; beta
< nrIRs(); beta
++ ) {
221 foreach( const Neighbor
&alpha
, nbIR(beta
) ) {
222 size_t _beta
= alpha
.dual
;
223 muab( alpha
, _beta
) = _Qa
[alpha
].marginal(IR(beta
)).divided_by( muba(alpha
,_beta
) );
227 foreach( const Neighbor
&alpha
, nbIR(beta
) ) {
228 size_t _beta
= alpha
.dual
;
229 Qb_new
*= muab(alpha
,_beta
) ^ (1 / (nbIR(beta
).size() + IR(beta
).c()));
232 Qb_new
.normalize( _normtype
);
233 if( Qb_new
.hasNaNs() ) {
234 cout
<< "HAK::doGBP: Qb_new has NaNs!" << endl
;
237 // _Qb[beta] = Qb_new.makeZero(1e-100); // damping?
240 foreach( const Neighbor
&alpha
, nbIR(beta
) ) {
241 size_t _beta
= alpha
.dual
;
243 muba(alpha
,_beta
) = _Qb
[beta
].divided_by( muab(alpha
,_beta
) );
245 Factor Qa_new
= OR(alpha
);
246 foreach( const Neighbor
&gamma
, nbOR(alpha
) )
247 Qa_new
*= muba(alpha
,gamma
.iter
);
248 Qa_new
^= (1.0 / OR(alpha
).c());
249 Qa_new
.normalize( _normtype
);
250 if( Qa_new
.hasNaNs() ) {
251 cout
<< "HAK::doGBP: Qa_new has NaNs!" << endl
;
254 // _Qa[alpha] = Qa_new.makeZero(1e-100); // damping?
259 // Calculate new single variable beliefs and compare with old ones
260 for( size_t i
= 0; i
< nrVars(); i
++ ) {
261 Factor new_belief
= belief( var( i
) );
262 diffs
.push( dist( new_belief
, old_beliefs
[i
], Prob::DISTLINF
) );
263 old_beliefs
[i
] = new_belief
;
266 if( props
.verbose
>= 3 )
267 cout
<< "HAK::doGBP: maxdiff " << diffs
.maxDiff() << " after " << iter
+1 << " passes" << endl
;
270 if( diffs
.maxDiff() > maxdiff
)
271 maxdiff
= diffs
.maxDiff();
273 if( props
.verbose
>= 1 ) {
274 if( diffs
.maxDiff() > props
.tol
) {
275 if( props
.verbose
== 1 )
277 cout
<< "HAK::doGBP: WARNING: not converged within " << props
.maxiter
<< " passes (" << toc() - tic
<< " clocks)...final maxdiff:" << diffs
.maxDiff() << endl
;
279 if( props
.verbose
>= 2 )
280 cout
<< "HAK::doGBP: ";
281 cout
<< "converged in " << iter
<< " passes (" << toc() - tic
<< " clocks)." << endl
;
285 return diffs
.maxDiff();
289 double HAK::doDoubleLoop() {
290 if( props
.verbose
>= 1 )
291 cout
<< "Starting " << identify() << "...";
292 if( props
.verbose
>= 3)
297 // Save original outer regions
298 vector
<FRegion
> org_ORs
= ORs
;
300 // Save original inner counting numbers and set negative counting numbers to zero
301 vector
<double> org_IR_cs( nrIRs(), 0.0 );
302 for( size_t beta
= 0; beta
< nrIRs(); beta
++ ) {
303 org_IR_cs
[beta
] = IR(beta
).c();
304 if( IR(beta
).c() < 0.0 )
308 // Keep old beliefs to check convergence
309 vector
<Factor
> old_beliefs
;
310 old_beliefs
.reserve( nrVars() );
311 for( size_t i
= 0; i
< nrVars(); i
++ )
312 old_beliefs
.push_back( belief( var(i
) ) );
314 // Differences in single node beliefs
315 Diffs
diffs(nrVars(), 1.0);
317 size_t outer_maxiter
= props
.maxiter
;
318 double outer_tol
= props
.tol
;
319 size_t outer_verbose
= props
.verbose
;
320 double org_maxdiff
= maxdiff
;
322 // Set parameters for inner loop
324 props
.verbose
= outer_verbose
? outer_verbose
- 1 : 0;
326 size_t outer_iter
= 0;
327 for( outer_iter
= 0; outer_iter
< outer_maxiter
&& diffs
.maxDiff() > outer_tol
; outer_iter
++ ) {
328 // Calculate new outer regions
329 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ ) {
330 OR(alpha
) = org_ORs
[alpha
];
331 foreach( const Neighbor
&beta
, nbOR(alpha
) )
332 OR(alpha
) *= _Qb
[beta
] ^ ((IR(beta
).c() - org_IR_cs
[beta
]) / nbIR(beta
).size());
336 if( isnan( doGBP() ) )
339 // Calculate new single variable beliefs and compare with old ones
340 for( size_t i
= 0; i
< nrVars(); ++i
) {
341 Factor new_belief
= belief( var( i
) );
342 diffs
.push( dist( new_belief
, old_beliefs
[i
], Prob::DISTLINF
) );
343 old_beliefs
[i
] = new_belief
;
346 if( props
.verbose
>= 3 )
347 cout
<< "HAK::doDoubleLoop: maxdiff " << diffs
.maxDiff() << " after " << outer_iter
+1 << " passes" << endl
;
350 // restore _maxiter, _verbose and _maxdiff
351 props
.maxiter
= outer_maxiter
;
352 props
.verbose
= outer_verbose
;
353 maxdiff
= org_maxdiff
;
355 if( diffs
.maxDiff() > maxdiff
)
356 maxdiff
= diffs
.maxDiff();
358 // Restore original outer regions
361 // Restore original inner counting numbers
362 for( size_t beta
= 0; beta
< nrIRs(); ++beta
)
363 IR(beta
).c() = org_IR_cs
[beta
];
365 if( props
.verbose
>= 1 ) {
366 if( diffs
.maxDiff() > props
.tol
) {
367 if( props
.verbose
== 1 )
369 cout
<< "HAK::doDoubleLoop: WARNING: not converged within " << outer_maxiter
<< " passes (" << toc() - tic
<< " clocks)...final maxdiff:" << diffs
.maxDiff() << endl
;
371 if( props
.verbose
>= 3 )
372 cout
<< "HAK::doDoubleLoop: ";
373 cout
<< "converged in " << outer_iter
<< " passes (" << toc() - tic
<< " clocks)." << endl
;
377 return diffs
.maxDiff();
382 if( props
.doubleloop
)
383 return doDoubleLoop();
389 Factor
HAK::belief( const VarSet
&ns
) const {
390 vector
<Factor
>::const_iterator beta
;
391 for( beta
= _Qb
.begin(); beta
!= _Qb
.end(); beta
++ )
392 if( beta
->vars() >> ns
)
394 if( beta
!= _Qb
.end() )
395 return( beta
->marginal(ns
) );
397 vector
<Factor
>::const_iterator alpha
;
398 for( alpha
= _Qa
.begin(); alpha
!= _Qa
.end(); alpha
++ )
399 if( alpha
->vars() >> ns
)
401 assert( alpha
!= _Qa
.end() );
402 return( alpha
->marginal(ns
) );
407 Factor
HAK::belief( const Var
&n
) const {
408 return belief( (VarSet
)n
);
412 vector
<Factor
> HAK::beliefs() const {
413 vector
<Factor
> result
;
414 for( size_t beta
= 0; beta
< nrIRs(); beta
++ )
415 result
.push_back( Qb(beta
) );
416 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ )
417 result
.push_back( Qa(alpha
) );
422 Complex
HAK::logZ() const {
424 for( size_t beta
= 0; beta
< nrIRs(); beta
++ )
425 sum
+= Complex(IR(beta
).c()) * Qb(beta
).entropy();
426 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ ) {
427 sum
+= Complex(OR(alpha
).c()) * Qa(alpha
).entropy();
428 sum
+= (OR(alpha
).log0() * Qa(alpha
)).totalSum();
434 } // end of namespace dai