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>
26 #include <dai/exceptions.h>
35 const char *HAK::Name
= "HAK";
38 bool HAK::checkProperties() {
39 if( !HasProperty("tol") )
41 if (!HasProperty("maxiter") )
43 if (!HasProperty("verbose") )
45 if( !HasProperty("doubleloop") )
47 if( !HasProperty("clusters") )
50 ConvertPropertyTo
<double>("tol");
51 ConvertPropertyTo
<size_t>("maxiter");
52 ConvertPropertyTo
<size_t>("verbose");
53 ConvertPropertyTo
<bool>("doubleloop");
54 ConvertPropertyTo
<ClustersType
>("clusters");
56 if( HasProperty("loopdepth") )
57 ConvertPropertyTo
<size_t>("loopdepth");
58 else if( Clusters() == ClustersType::LOOP
)
65 void HAK::constructMessages() {
66 // Create outer beliefs
69 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ )
70 _Qa
.push_back( Factor( OR(alpha
).vars() ) );
72 // Create inner beliefs
75 for( size_t beta
= 0; beta
< nrIRs(); beta
++ )
76 _Qb
.push_back( Factor( IR(beta
) ) );
80 _muab
.reserve( nrORs() );
82 _muba
.reserve( nrORs() );
83 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ ) {
84 _muab
.push_back( vector
<Factor
>() );
85 _muba
.push_back( vector
<Factor
>() );
86 _muab
[alpha
].reserve( nbOR(alpha
).size() );
87 _muba
[alpha
].reserve( nbOR(alpha
).size() );
88 foreach( const Neighbor
&beta
, nbOR(alpha
) ) {
89 _muab
[alpha
].push_back( Factor( IR(beta
) ) );
90 _muba
[alpha
].push_back( Factor( IR(beta
) ) );
96 HAK::HAK(const RegionGraph
& rg
, const Properties
&opts
) : DAIAlgRG(rg
, opts
) {
97 assert( checkProperties() );
103 void HAK::findLoopClusters( const FactorGraph
& fg
, std::set
<VarSet
> &allcl
, VarSet newcl
, const Var
& root
, size_t length
, VarSet vars
) {
104 for( VarSet::const_iterator in
= vars
.begin(); in
!= vars
.end(); in
++ ) {
105 VarSet ind
= fg
.delta( fg
.findVar( *in
) );
106 if( (newcl
.size()) >= 2 && (ind
>> root
) ) {
107 allcl
.insert( newcl
| *in
);
109 else if( length
> 1 )
110 findLoopClusters( fg
, allcl
, newcl
| *in
, root
, length
- 1, ind
/ newcl
);
115 HAK::HAK(const FactorGraph
& fg
, const Properties
&opts
) : DAIAlgRG(opts
) {
116 assert( checkProperties() );
119 if( Clusters() == ClustersType::MIN
) {
121 } else if( Clusters() == ClustersType::DELTA
) {
122 for( size_t i
= 0; i
< fg
.nrVars(); i
++ )
123 cl
.push_back(fg
.Delta(i
));
124 } else if( Clusters() == ClustersType::LOOP
) {
127 for( size_t i0
= 0; i0
< fg
.nrVars(); i0
++ ) {
128 VarSet i0d
= fg
.delta(i0
);
129 if( LoopDepth() > 1 )
130 findLoopClusters( fg
, scl
, fg
.var(i0
), fg
.var(i0
), LoopDepth() - 1, fg
.delta(i0
) );
132 for( set
<VarSet
>::const_iterator c
= scl
.begin(); c
!= scl
.end(); c
++ )
134 if( Verbose() >= 3 ) {
135 cout
<< "HAK uses the following clusters: " << endl
;
136 for( vector
<VarSet
>::const_iterator cli
= cl
.begin(); cli
!= cl
.end(); cli
++ )
137 cout
<< *cli
<< endl
;
140 DAI_THROW(INTERNAL_ERROR
);
142 RegionGraph
rg(fg
,cl
);
143 RegionGraph::operator=(rg
);
147 cout
<< "HAK regiongraph: " << *this << endl
;
151 string
HAK::identify() const {
152 stringstream
result (stringstream::out
);
153 result
<< Name
<< GetProperties();
158 void HAK::init( const VarSet
&ns
) {
159 for( vector
<Factor
>::iterator alpha
= _Qa
.begin(); alpha
!= _Qa
.end(); alpha
++ )
160 if( alpha
->vars().intersects( ns
) )
161 alpha
->fill( 1.0 / alpha
->states() );
163 for( size_t beta
= 0; beta
< nrIRs(); beta
++ )
164 if( IR(beta
).intersects( ns
) ) {
165 _Qb
[beta
].fill( 1.0 );
166 foreach( const Neighbor
&alpha
, nbIR(beta
) ) {
167 size_t _beta
= alpha
.dual
;
168 muab( alpha
, _beta
).fill( 1.0 / IR(beta
).states() );
169 muba( alpha
, _beta
).fill( 1.0 / IR(beta
).states() );
176 assert( checkProperties() );
178 for( vector
<Factor
>::iterator alpha
= _Qa
.begin(); alpha
!= _Qa
.end(); alpha
++ )
179 alpha
->fill( 1.0 / alpha
->states() );
181 for( vector
<Factor
>::iterator beta
= _Qb
.begin(); beta
!= _Qb
.end(); beta
++ )
182 beta
->fill( 1.0 / beta
->states() );
184 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ )
185 foreach( const Neighbor
&beta
, nbOR(alpha
) ) {
186 size_t _beta
= beta
.iter
;
187 muab( alpha
, _beta
).fill( 1.0 / muab( alpha
, _beta
).states() );
188 muba( alpha
, _beta
).fill( 1.0 / muab( alpha
, _beta
).states() );
193 double HAK::doGBP() {
195 cout
<< "Starting " << identify() << "...";
201 // Check whether counting numbers won't lead to problems
202 for( size_t beta
= 0; beta
< nrIRs(); beta
++ )
203 assert( nbIR(beta
).size() + IR(beta
).c() != 0.0 );
205 // Keep old beliefs to check convergence
206 vector
<Factor
> old_beliefs
;
207 old_beliefs
.reserve( nrVars() );
208 for( size_t i
= 0; i
< nrVars(); i
++ )
209 old_beliefs
.push_back( belief( var(i
) ) );
211 // Differences in single node beliefs
212 Diffs
diffs(nrVars(), 1.0);
215 // do several passes over the network until maximum number of iterations has
216 // been reached or until the maximum belief difference is smaller than tolerance
217 for( iter
= 0; iter
< MaxIter() && diffs
.maxDiff() > Tol(); iter
++ ) {
218 for( size_t beta
= 0; beta
< nrIRs(); beta
++ ) {
219 foreach( const Neighbor
&alpha
, nbIR(beta
) ) {
220 size_t _beta
= alpha
.dual
;
221 muab( alpha
, _beta
) = _Qa
[alpha
].marginal(IR(beta
)).divided_by( muba(alpha
,_beta
) );
225 foreach( const Neighbor
&alpha
, nbIR(beta
) ) {
226 size_t _beta
= alpha
.dual
;
227 Qb_new
*= muab(alpha
,_beta
) ^ (1 / (nbIR(beta
).size() + IR(beta
).c()));
230 Qb_new
.normalize( _normtype
);
231 if( Qb_new
.hasNaNs() ) {
232 cout
<< "HAK::doGBP: Qb_new has NaNs!" << endl
;
235 // _Qb[beta] = Qb_new.makeZero(1e-100); // damping?
238 foreach( const Neighbor
&alpha
, nbIR(beta
) ) {
239 size_t _beta
= alpha
.dual
;
241 muba(alpha
,_beta
) = _Qb
[beta
].divided_by( muab(alpha
,_beta
) );
243 Factor Qa_new
= OR(alpha
);
244 foreach( const Neighbor
&gamma
, nbOR(alpha
) )
245 Qa_new
*= muba(alpha
,gamma
.iter
);
246 Qa_new
^= (1.0 / OR(alpha
).c());
247 Qa_new
.normalize( _normtype
);
248 if( Qa_new
.hasNaNs() ) {
249 cout
<< "HAK::doGBP: Qa_new has NaNs!" << endl
;
252 // _Qa[alpha] = Qa_new.makeZero(1e-100); // damping?
257 // Calculate new single variable beliefs and compare with old ones
258 for( size_t i
= 0; i
< nrVars(); i
++ ) {
259 Factor new_belief
= belief( var( i
) );
260 diffs
.push( dist( new_belief
, old_beliefs
[i
], Prob::DISTLINF
) );
261 old_beliefs
[i
] = new_belief
;
265 cout
<< "HAK::doGBP: maxdiff " << diffs
.maxDiff() << " after " << iter
+1 << " passes" << endl
;
268 updateMaxDiff( diffs
.maxDiff() );
270 if( Verbose() >= 1 ) {
271 if( diffs
.maxDiff() > Tol() ) {
274 cout
<< "HAK::doGBP: WARNING: not converged within " << MaxIter() << " passes (" << toc() - tic
<< " clocks)...final maxdiff:" << diffs
.maxDiff() << endl
;
277 cout
<< "HAK::doGBP: ";
278 cout
<< "converged in " << iter
<< " passes (" << toc() - tic
<< " clocks)." << endl
;
282 return diffs
.maxDiff();
286 double HAK::doDoubleLoop() {
288 cout
<< "Starting " << identify() << "...";
294 // Save original outer regions
295 vector
<FRegion
> org_ORs
= ORs
;
297 // Save original inner counting numbers and set negative counting numbers to zero
298 vector
<double> org_IR_cs( nrIRs(), 0.0 );
299 for( size_t beta
= 0; beta
< nrIRs(); beta
++ ) {
300 org_IR_cs
[beta
] = IR(beta
).c();
301 if( IR(beta
).c() < 0.0 )
305 // Keep old beliefs to check convergence
306 vector
<Factor
> old_beliefs
;
307 old_beliefs
.reserve( nrVars() );
308 for( size_t i
= 0; i
< nrVars(); i
++ )
309 old_beliefs
.push_back( belief( var(i
) ) );
311 // Differences in single node beliefs
312 Diffs
diffs(nrVars(), 1.0);
314 size_t outer_maxiter
= MaxIter();
315 double outer_tol
= Tol();
316 size_t outer_verbose
= Verbose();
317 double org_maxdiff
= MaxDiff();
319 // Set parameters for inner loop
321 Verbose( outer_verbose
? outer_verbose
- 1 : 0 );
323 size_t outer_iter
= 0;
324 for( outer_iter
= 0; outer_iter
< outer_maxiter
&& diffs
.maxDiff() > outer_tol
; outer_iter
++ ) {
325 // Calculate new outer regions
326 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ ) {
327 OR(alpha
) = org_ORs
[alpha
];
328 foreach( const Neighbor
&beta
, nbOR(alpha
) )
329 OR(alpha
) *= _Qb
[beta
] ^ ((IR(beta
).c() - org_IR_cs
[beta
]) / nbIR(beta
).size());
333 if( isnan( doGBP() ) )
336 // Calculate new single variable beliefs and compare with old ones
337 for( size_t i
= 0; i
< nrVars(); ++i
) {
338 Factor new_belief
= belief( var( i
) );
339 diffs
.push( dist( new_belief
, old_beliefs
[i
], Prob::DISTLINF
) );
340 old_beliefs
[i
] = new_belief
;
344 cout
<< "HAK::doDoubleLoop: maxdiff " << diffs
.maxDiff() << " after " << outer_iter
+1 << " passes" << endl
;
347 // restore _maxiter, _verbose and _maxdiff
348 MaxIter( outer_maxiter
);
349 Verbose( outer_verbose
);
350 MaxDiff( org_maxdiff
);
352 updateMaxDiff( diffs
.maxDiff() );
354 // Restore original outer regions
357 // Restore original inner counting numbers
358 for( size_t beta
= 0; beta
< nrIRs(); ++beta
)
359 IR(beta
).c() = org_IR_cs
[beta
];
361 if( Verbose() >= 1 ) {
362 if( diffs
.maxDiff() > Tol() ) {
365 cout
<< "HAK::doDoubleLoop: WARNING: not converged within " << outer_maxiter
<< " passes (" << toc() - tic
<< " clocks)...final maxdiff:" << diffs
.maxDiff() << endl
;
368 cout
<< "HAK::doDoubleLoop: ";
369 cout
<< "converged in " << outer_iter
<< " passes (" << toc() - tic
<< " clocks)." << endl
;
373 return diffs
.maxDiff();
379 return doDoubleLoop();
385 Factor
HAK::belief( const VarSet
&ns
) const {
386 vector
<Factor
>::const_iterator beta
;
387 for( beta
= _Qb
.begin(); beta
!= _Qb
.end(); beta
++ )
388 if( beta
->vars() >> ns
)
390 if( beta
!= _Qb
.end() )
391 return( beta
->marginal(ns
) );
393 vector
<Factor
>::const_iterator alpha
;
394 for( alpha
= _Qa
.begin(); alpha
!= _Qa
.end(); alpha
++ )
395 if( alpha
->vars() >> ns
)
397 assert( alpha
!= _Qa
.end() );
398 return( alpha
->marginal(ns
) );
403 Factor
HAK::belief( const Var
&n
) const {
404 return belief( (VarSet
)n
);
408 vector
<Factor
> HAK::beliefs() const {
409 vector
<Factor
> result
;
410 for( size_t beta
= 0; beta
< nrIRs(); beta
++ )
411 result
.push_back( Qb(beta
) );
412 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ )
413 result
.push_back( Qa(alpha
) );
418 Complex
HAK::logZ() const {
420 for( size_t beta
= 0; beta
< nrIRs(); beta
++ )
421 sum
+= Complex(IR(beta
).c()) * Qb(beta
).entropy();
422 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ ) {
423 sum
+= Complex(OR(alpha
).c()) * Qa(alpha
).entropy();
424 sum
+= (OR(alpha
).log0() * Qa(alpha
)).totalSum();
430 } // end of namespace dai