977a563831d6d28ac2b3c499585ba257e325ede1
[libdai.git] / src / glc.cpp
1 /* This file is part of libDAI - http://www.libdai.org/
2 *
3 * Copyright (c) 2006-2012, The libDAI authors. All rights reserved.
4 *
5 * Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
6 */
7
8
9 #include <iostream>
10 #include <algorithm>
11 #include <map>
12 #include <set>
13 #include <dai/glc.h>
14 #include <dai/util.h>
15 #include <dai/alldai.h>
16
17
18 namespace dai {
19
20
21 using namespace std;
22
23
24 void GLC::setProperties( const PropertySet &opts ) {
25 DAI_ASSERT( opts.hasKey("tol") );
26 DAI_ASSERT( opts.hasKey("maxiter") );
27 DAI_ASSERT( opts.hasKey("cavity") );
28 DAI_ASSERT( opts.hasKey("updates") );
29 DAI_ASSERT( opts.hasKey("rgntype") );
30 DAI_ASSERT( opts.hasKey("cavainame") );
31 DAI_ASSERT( opts.hasKey("cavaiopts") );
32 DAI_ASSERT( opts.hasKey("inainame") );
33 DAI_ASSERT( opts.hasKey("inaiopts") );
34
35 if( opts.hasKey("maxtime") )
36 props.maxtime = opts.getStringAs<Real>("maxtime");
37 else
38 props.maxtime = INFINITY;
39
40 props.tol = opts.getStringAs<Real>("tol");
41 props.maxiter = opts.getStringAs<size_t>("maxiter");
42 props.cavity = opts.getStringAs<Properties::CavityType>("cavity");
43 props.updates = opts.getStringAs<Properties::UpdateType>("updates");
44 props.rgntype = opts.getStringAs<Properties::RegionType>("rgntype");
45 if(opts.hasKey("neighbors"))
46 props.neighbors = opts.getStringAs<CobwebGraph::NeighborType>("neighbors");
47 else
48 props.neighbors = CobwebGraph::NeighborType::CLOSEST;
49 if( opts.hasKey("verbose") )
50 props.verbose = opts.getStringAs<size_t>("verbose");
51 else
52 props.verbose = 0;
53 if( opts.hasKey("loopdepth") )
54 props.loopdepth = opts.getStringAs<size_t>("loopdepth");
55 else
56 props.loopdepth = 4;
57 if( opts.hasKey("cavainame") )
58 props.cavainame = opts.getStringAs<string>("cavainame");
59 if( opts.hasKey("cavaiopts") )
60 props.cavaiopts = opts.getStringAs<PropertySet>("cavaiopts");
61 if( opts.hasKey("inainame") )
62 props.inainame = opts.getStringAs<string>("inainame");
63 if( opts.hasKey("inaiopts") )
64 props.inaiopts = opts.getStringAs<PropertySet>("inaiopts");
65 if( opts.hasKey("reinit") )
66 props.reinit = opts.getStringAs<bool>("reinit");
67 }
68
69
70 PropertySet GLC::getProperties() const {
71 PropertySet opts;
72 opts.set( "tol", props.tol );
73 opts.set( "maxiter", props.maxiter );
74 opts.set( "verbose", props.verbose );
75 opts.set( "loopdepth", props.loopdepth );
76 opts.set( "cavity", props.cavity );
77 opts.set( "neighbors", props.neighbors );
78 opts.set( "updates", props.updates );
79 opts.set( "rgntype", props.rgntype );
80 opts.set( "cavainame", props.cavainame );
81 opts.set( "cavaiopts", props.cavaiopts );
82 opts.set( "inainame", props.inainame );
83 opts.set( "inaiopts", props.inaiopts );
84 opts.set( "reinit", props.reinit );
85 opts.set( "maxtime", props.maxtime );
86
87 return opts;
88 }
89
90
91 string GLC::printProperties() const {
92 stringstream s( stringstream::out );
93 s << "[";
94 s << "tol=" << props.tol << ",";
95 s << "maxiter=" << props.maxiter << ",";
96 s << "verbose=" << props.verbose << ",";
97 s << "loopdepth=" << props.loopdepth << ",";
98 s << "cavity=" << props.cavity << ",";
99 s << "updates=" << props.updates << ",";
100 s << "neighbors=" << props.neighbors << ",";
101 s << "rgntype=" << props.rgntype << ",";
102 s << "cavainame=" << props.cavainame << ",";
103 s << "cavaiopts=" << props.cavaiopts << ",";
104 s << "inainame=" << props.inainame << ",";
105 s << "inaiopts=" << props.inaiopts << ",";
106 s << "reinit=" << props.reinit << ",";
107 s << "maxtime=" << props.maxtime << "]";
108 return s.str();
109 }
110
111
112 GLC::GLC( const FactorGraph& fg, const PropertySet &opts ) : DAIAlgCG(fg), _CWs(),_outOffset(), _cavitydists(), _beliefs(), _factorBeliefs(), _maxdiff(0.0), _iters(0), props() {
113 setProperties( opts );
114 setRgn( calcRegions(), props.neighbors, false );
115 if( props.verbose >= 3 ){
116 cerr << "#regions: " << nrCWs() << endl;
117 for( size_t i = 0; i < nrCWs(); i++ )
118 cerr << INRs(i).elements() << endl;
119 }
120 if( props.verbose >= 3 )
121 cerr << "initializing cavity dists" << endl;
122 _cavitydists.resize( nrCWs() );
123 _maxdiff = InitCavityDists( props.cavainame, props.cavaiopts );
124 // build the vector of Cobweb Regions
125 setCWs( props.inainame, props.inaiopts );
126 if( props.verbose >= 3 )
127 cerr << "Regions are built" << endl;
128 // create beliefs
129 _beliefs.reserve( nrVars() );
130 for( size_t i = 0; i < nrVars(); i++ )
131 _beliefs.push_back( Factor(var(i)) );
132 }
133
134
135 vector<SmallSet<size_t> > GLC::calcRegions() {
136 // contains the variable indices for each region
137 vector<SmallSet<size_t> > regions;
138 // single variables
139 if( props.rgntype == Properties::RegionType::SINGLE ) {
140 regions.resize( nrVars() );
141 for( size_t i = 0; i < nrVars(); i++ )
142 regions[i] = SmallSet<size_t>(i);
143 // partitioning factors
144 } else if( props.rgntype == Properties::RegionType::FACTOR ) {
145 SmallSet<size_t> remainingVars;
146 for( size_t i = 0; i < nrVars(); i++ )
147 remainingVars.insert( i );
148 vector<SmallSet<size_t> > facVars( nrFactors() );
149 for( size_t I = 0; I < nrFactors(); I++ )
150 facVars[I] = findVars( factor(I).vars() );
151 bool canadd = true;
152 while( canadd ) {
153 canadd = false;
154 for( size_t I = 0; I < nrFactors(); I++ ) {
155 if( facVars[I] << remainingVars ) {
156 regions.push_back( facVars[I] );
157 remainingVars /= facVars[I];
158 canadd = true;
159 break;
160 }
161 }
162 }
163 // add the variables that are not covered as single variable regions
164 bforeach(size_t i, remainingVars)
165 regions.push_back( SmallSet<size_t>(i) );
166 // overlapping factors (non-partitioning)
167 } else if( props.rgntype == Properties::RegionType::OVFACTOR ) {
168 SmallSet<size_t> remainingVars;
169 for( size_t I = 0; I < nrFactors(); I++ )
170 regions.push_back( findVars( factor(I).vars() ) );
171 // partitioning: adds loops over variables that are not covered recursively
172 } else if( props.rgntype == Properties::RegionType::LOOP ) {
173 set<SmallSet<size_t> > scl; // to contain clusters
174 SmallSet<size_t> remaining; // variables not in a cluster yet
175 // build the set of all var inds
176 for( size_t v = 0; v < nrVars(); v++ )
177 remaining.insert( v );
178 // while there for a remaining var we can find a loop
179 bool changing = true;
180 while( changing ) {
181 changing = false;
182 for( size_t i0 = 0; i0 < nrVars(); i0++ ) {
183 if( !remaining.contains(i0) )
184 continue;
185 size_t currentsize = scl.size();
186 if( props.loopdepth > 1 )
187 findLoopClusters( remaining, scl, SmallSet<size_t>(i0),i0, props.loopdepth - 1, deltai(i0) & remaining );
188 if(scl.size() > currentsize)
189 changing = true;
190 }
191 }
192 copy( scl.begin(), scl.end(), back_inserter(regions) );
193 // add the vars that are not covered as single variable regions
194 bforeach( size_t i, remaining )
195 regions.push_back( SmallSet<size_t>(i) );
196 // adds overlapping loops of maximum size loopdepth recursively
197 } else if( props.rgntype == Properties::RegionType::OVLOOP ) {
198 for( size_t I = 0; I < nrFactors(); I++ )
199 regions.push_back( findVars( factor(I).vars() ) );
200
201 SmallSet<size_t> remaining; // variables not in a cluster yet
202 // build the set of all var inds
203 for( size_t v = 0; v < nrVars(); v++ )
204 remaining.insert( v );
205
206 set<SmallSet<size_t> > scl; // to contain clusters
207 // while there for a remaining var we can find a loop
208
209 for( size_t i0 = 0; i0 < nrVars(); i0++ )
210 findOVLoopClusters( remaining, scl, SmallSet<size_t>(i0),i0, props.loopdepth - 1, deltai(i0) );
211 copy( scl.begin(), scl.end(), back_inserter(regions) );
212 // non-partitioning: overlapping (variable+markov blanket)s
213 } else if( props.rgntype == Properties::RegionType::OVDELTA ) {
214 for( size_t I = 0; I < nrFactors(); I++ )
215 regions.push_back( findVars(factor(I).vars()) );
216 SmallSet<size_t> remaining; // variables not in a cluster yet
217 // build the set of all var inds
218 for( size_t v = 0; v < nrVars(); v++ )
219 remaining.insert( v );
220 set<SmallSet<size_t> > scl; // to contain clusters
221 // while there for a remaining var we can find a loop
222 for( size_t i0 = 0; i0 < nrVars(); i0++ )
223 findOVLoopClusters( remaining, scl, SmallSet<size_t>(i0),i0, props.loopdepth - 1, deltai(i0) );
224 copy( scl.begin(), scl.end(), back_inserter(regions) );
225 } else if( props.rgntype == Properties::RegionType::DELTA ) {
226 SmallSet<size_t> remaining; // variables not in a cluster yet
227 // build the set of all var inds
228 for( size_t v = 0; v < nrVars(); v++ )
229 remaining.insert( v );
230 while( remaining.size() > 0 ) {
231 size_t cv = remaining.elements()[ rand() % remaining.size() ];
232 regions.push_back( Deltai(cv) & remaining );
233 remaining /= Deltai(cv);
234 }
235 }
236 return regions;
237 }
238
239
240 void GLC::findOVLoopClusters( SmallSet<size_t>& remaining, set<SmallSet<size_t> > &allcl, SmallSet<size_t> newcl, const size_t& root, size_t length, SmallSet<size_t> vars ) {
241 for( SmallSet<size_t>::const_iterator in = vars.begin(); in != vars.end(); in++ ) {
242 SmallSet<size_t> ind = deltai( *in );
243 if( (newcl.size()) >= 2 && ind.contains( root ) ) {
244 allcl.insert( newcl | *in );
245 remaining /= (newcl | *in);
246 } else if( length > 1 )
247 findOVLoopClusters( remaining, allcl, newcl | *in, root, length - 1, ind / newcl );
248 }
249 }
250
251
252 void GLC::findLoopClusters( SmallSet<size_t>& remaining, set<SmallSet<size_t> > &allcl, SmallSet<size_t> newcl, const size_t& root, size_t length, SmallSet<size_t> vars ) {
253 for( SmallSet<size_t>::const_iterator in = vars.begin(); in != vars.end(); in++ ) {
254 SmallSet<size_t> ind = deltai( *in ) & remaining;
255 if( (newcl.size()) >= 2 && ind.contains( root ) ) {
256 allcl.insert( newcl | *in );
257 remaining /= (newcl | *in);
258 break;
259 } else if( length > 1 )
260 findLoopClusters( remaining, allcl, newcl | *in, root, length - 1, ind / newcl );
261 }
262 }
263
264
265 void GLC::CalcBelief( size_t i, bool isFinal ) {
266 if( !isFinal )
267 _beliefs[i] = CW( var2CW(i)[0] ).belief( var(i) );
268 else{ // take geometric average
269 _beliefs[i] = Factor(var(i), (Real)1);
270 for( size_t rind = 0; rind < var2CW(i).size(); rind++ )
271 _beliefs[i] *= CW( var2CW(i)[rind] ).belief( var(i) );
272 _beliefs[i] ^= (Real)1 / var2CW(i).size();
273 }
274 }
275
276
277 void GLC::CalcFactorBelief( size_t I ) {
278 VarSet ns = factor(I).vars();
279 for( size_t R = 0; R < nrCWs(); R++ )
280 if( ns << inds2vars( INRs(R).elements() ) )
281 _factorBeliefs[ns] = CW(R).belief(ns);
282 }
283
284
285 Factor GLC::belief( const VarSet &ns ) const {
286 if( ns.size() == 0 )
287 return Factor();
288 else if( ns.size() == 1 )
289 return beliefV( findVar( *(ns.begin()) ) );
290 else {
291 map<VarSet, Factor>::const_iterator it = _factorBeliefs.find(ns);
292 if( it != _factorBeliefs.end() )
293 return it->second;
294 }
295 DAI_THROW(BELIEF_NOT_AVAILABLE);
296 return Factor();
297 }
298
299
300 vector<Factor> GLC::beliefs() const {
301 vector<Factor> result = _beliefs;
302 for( map<VarSet, Factor>::const_iterator it = _factorBeliefs.begin(); it != _factorBeliefs.end(); it++ )
303 result.push_back( it->second );
304 return result;
305 }
306
307
308 Real GLC::CalcCavityDist( size_t R, const std::string &name, const PropertySet& opts ) {
309 vector<Factor> pairbeliefs;
310 Real maxdiff = 0;
311
312 if( props.cavity != Properties::CavityType::UNIFORM ) {
313 InfAlg *cav = newInfAlg( name, *this, opts );
314 cav->makeRegionCavity( Rfs(R).elements() );
315
316 if( props.cavity == Properties::CavityType::FULL )
317 pairbeliefs.push_back( calcMarginal( *cav, inds2vars(EXRs(R).elements()), props.reinit ) );
318 if( props.cavity == Properties::CavityType::PAIR )
319 pairbeliefs = calcPairBeliefs( *cav, inds2vars(EXRs(R).elements()), props.reinit, false );
320 else if( props.cavity == Properties::CavityType::PAIR2 )
321 pairbeliefs = calcPairBeliefs( *cav, inds2vars(EXRs(R).elements()), props.reinit, true );
322 maxdiff = cav->maxDiff();
323 delete cav;
324 }
325 if( props.verbose >= 3 )
326 cerr << "R:" << R << "cavity of size " << pairbeliefs.size() << endl;
327 _cavitydists[R] = pairbeliefs;
328 return maxdiff;
329 }
330
331
332 Real GLC::InitCavityDists( const std::string &name, const PropertySet &opts ) {
333 double tic = toc();
334
335 if( props.verbose >= 2 ) {
336 cerr << this->name() << "::InitCavityDists: ";
337 if( props.cavity == Properties::CavityType::UNIFORM )
338 cerr << "Using uniform initial cavity distributions" << endl;
339 else if( props.cavity == Properties::CavityType::FULL )
340 cerr << "Using full " << name << opts << "...";
341 else if( props.cavity == Properties::CavityType::PAIR )
342 cerr << "Using pairwise " << name << opts << "...";
343 else if( props.cavity == Properties::CavityType::PAIR2 )
344 cerr << "Using pairwise(new) " << name << opts << "...";
345 }
346
347 Real maxdiff = 0.0;
348 for( size_t R = 0; R < nrCWs(); R++ ) {
349 Real md = CalcCavityDist(R, name, opts);
350 if( md > maxdiff )
351 maxdiff = md;
352 }
353 if( props.verbose >= 2 )
354 cerr << this->name() << "::InitCavityDists used " << toc() - tic << " seconds." << endl;
355 return maxdiff;
356 }
357
358
359 void GLC::NewPancake( size_t R, size_t _R2 ) {
360 Factor newmsg;
361 // calculate the new msg
362 newmsg = CW(M(R, _R2).his).marginal( M(R, _R2).msg.vars(), M(R, _R2).fc) / ((CW(R).marginal( M(R, _R2).msg.vars(), M(R, _R2).fc))*M(R, _R2).msg.inverse());
363 newmsg.normalize();
364 // update the cobweb-graph with this new msg
365 M(R, _R2).msg = newmsg;
366 // update the corresponding factor in the cobweb region so to be used in future marginalizations (-_R2 - 1) is the index (see cobweb class documents)
367 CW(R).updateFactor( -_R2 - 1 ,newmsg);
368 }
369
370
371 void GLC::OVNewPancake( size_t R ) {
372 Factor newmsg, allmsgs, marg;
373 for( size_t _R2 = 0; _R2 < M(R).size(); _R2++ ) { // for all _R2 that send messages to R
374 // calculate the message
375 newmsg = (CW(M(R, _R2).his).marginal( M(R, _R2).msg.vars(), M(R, _R2).fc) / ((CW(R).marginal( M(R, _R2).msg.vars(), M(R, _R2).fc))))*M(R, _R2).msg;
376 newmsg.normalize();
377 // allmsgs takes care of the double-counts the way using the region-graph does in the paper. It is more compact but not as efficient
378 allmsgs = newmsg;
379 for( size_t sub = 0; sub < M(R,_R2).subregions.size(); sub++ ) { //for all descendants of rho of (\ominus r_{R, R2})
380 Real c = (Real)cn(R)[M(R, _R2).subregions[sub]].first; // counting number for rho
381 size_t nrParents = cn(R)[M(R, _R2).subregions[sub]].second.size(); // number of rho's parents
382 marg = newmsg.marginal(M(R, _R2).subregions[sub], false); // marginal of the message from R2 to R over rho
383 allmsgs *= (marg^(c / nrParents)); // discount the double counting of rho from the message
384 for( size_t pind = 0; pind < nrParents; pind++ ) // this operation corresponds to the upward pass
385 if( cn(R)[M(R, _R2).subregions[sub]].second[pind] != _R2 )
386 M(R, cn(R)[M(R, _R2).subregions[sub]].second[pind]).newmsg *= marg^(-c / (nrParents * (nrParents - (Real)1)));
387 }
388 // set the message after taking double-countings into account
389 M(R, _R2).msg = allmsgs.normalized();
390 }
391
392 for( size_t _R2 = 0; _R2 < M(R).size(); _R2++ ) {
393 M(R, _R2).newmsg *= M(R, _R2).msg;
394 M(R, _R2).newmsg.normalize();
395 // update the corresponding factor in the cobweb region: to be used in future marginalizations
396 CW(R).updateFactor( -_R2 - 1 , M(R, _R2).msg, false);
397 M(R, _R2).msg = M(R, _R2).newmsg;
398 M(R, _R2).newmsg.fill( (Real)1 );
399 }
400 }
401
402
403 Real GLC::run() {
404 if( props.verbose >= 2 )
405 cerr << "Starting " << identify() << "...";
406 if( props.verbose >= 2 )
407 cerr << endl;
408
409 double tic = toc();
410
411 if( props.verbose >= 3 )
412 cerr << "Initializing older beliefs" << endl;
413 vector<vector<Factor> > oldBeliefsM;
414 vector<Factor> oldBeliefsP;
415 if( isPartition ) {
416 oldBeliefsM.reserve( nrCWs() );
417 for( size_t R = 0; R < nrCWs(); R++ ) {
418 oldBeliefsM.push_back( vector<Factor>() );
419 oldBeliefsM[R].reserve( M(R).size() );
420 for( size_t m = 0; m < M(R).size(); m++ )
421 oldBeliefsM[R].push_back( M(R,m).msg );
422 }
423 } else {
424 oldBeliefsP.reserve( nrVars() );
425 for( size_t i = 0; i < nrVars(); i++ ) {
426 CalcBelief( i, true );
427 oldBeliefsP.push_back( _beliefs[i] );
428 }
429 }
430 if( props.verbose >= 3 )
431 cerr << "Setting the update sequence" <<endl;
432
433 size_t nredges = 0;
434 vector<Edge> update_seq;
435 update_seq.reserve( nredges );
436 for( size_t R = 0; R < nrCWs(); ++R )
437 for( size_t R2 = 0; R2 < M(R).size(); R2++ ) {
438 update_seq.push_back( Edge( R, R2 ) );
439 nredges++;
440 }
441
442 vector<size_t> update_seq_node;
443 update_seq_node.reserve( nrCWs() );
444 for( size_t R = 0; R < nrCWs(); ++R )
445 update_seq_node.push_back( R );
446 if( props.verbose >= 3 )
447 cerr << "Going into the main loop" << endl;
448
449 Real maxDiff = INFINITY;
450 for( _iters = 0; _iters < props.maxiter && maxDiff > props.tol && (toc() - tic) < props.maxtime; _iters++ ) {
451 if( props.updates == Properties::UpdateType::SEQRND ) {
452 random_shuffle( update_seq.begin(), update_seq.end() );
453 random_shuffle( update_seq_node.begin(), update_seq_node.end() );
454 }
455 if( isPartition ) {
456 for( size_t t = 0; t < nredges; t++ ) {
457 size_t R = update_seq[t].first;
458 size_t _R2 = update_seq[t].second;
459 NewPancake( R, _R2);
460 }
461 } else {
462 for( size_t t = 0; t < nrCWs(); t++ )
463 OVNewPancake( update_seq_node[t] );
464 }
465
466 if( !isPartition )
467 for( size_t i = 0; i < nrVars(); i++ )
468 CalcBelief( i , true );
469
470 maxDiff = -INFINITY;
471 if( isPartition ) {
472 for( size_t R = 0; R < nrCWs(); R++ ) {
473 for( size_t m = 0; m < M(R).size(); m++ ) {
474 maxDiff = std::max( maxDiff, dist( M(R, m).msg, oldBeliefsM[R][m], DISTLINF ) );
475 oldBeliefsM[R][m] = M(R,m).msg;
476 }
477 }
478 } else {
479 for( size_t i = 0; i < nrVars(); i++ ) {
480 maxDiff = std::max( maxDiff, dist( _beliefs[i], oldBeliefsP[i], DISTLINF ) );
481 oldBeliefsP[i] = _beliefs[i];
482 }
483 }
484
485 if( props.verbose >= 2 )
486 cerr << this->name() << "::run: maxdiff " << maxDiff << " after " << _iters+1 << " passes" << endl;
487 }
488
489 if( isPartition )
490 for( size_t i = 0; i < nrVars(); i++ )
491 CalcBelief( i, true );
492
493 if( maxDiff > _maxdiff )
494 _maxdiff = maxDiff;
495
496 if( props.verbose >= 2 ) {
497 if( maxDiff > props.tol ) {
498 if( props.verbose == 2 )
499 cerr << endl;
500 cerr << this->name() << "::run: WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " seconds)...final maxdiff:" << maxDiff << endl;
501 } else {
502 if( props.verbose >= 2 )
503 cerr << this->name() << "::run: ";
504 cerr << "converged in " << _iters << " passes (" << toc() - tic << " seconds)." << endl;
505 }
506 }
507
508 for( size_t I = 0; I < nrFactors(); I++ )
509 CalcFactorBelief( I );
510
511 return maxDiff;
512 }
513
514
515 void GLC::setCWs( const std::string &name, const PropertySet &opts ) {
516 // needs to generate a vector of relevant factors, and a dictionary for mapping between local and global factors
517 _CWs.clear();
518 _CWs.reserve( _INRs.size() );
519 _outOffset.clear();
520 _outOffset.reserve( _INRs.size() );
521 // main factors
522 if( props.verbose >= 3 )
523 cerr << "Setting CWs..." << endl;
524
525 for( size_t R = 0; R < nrCWs(); R++ ) { // for each region
526 map<int, size_t> g2l;
527 vector<Factor> Rfactors; // contains all factors
528 size_t potsize = Rfs(R).size() + M(R).size() + _cavitydists[R].size();
529 Rfactors.reserve( potsize );
530 size_t _I = 0;
531 // adding main factors
532 for( SmallSet<size_t>::const_iterator it = Rfs(R).begin(); it != Rfs(R).end(); it++ ) {
533 Rfactors.push_back( factor(*it) );
534 g2l[(int)(*it)] = _I;
535 _I++;
536 }
537
538 // adding factors for incoming messages
539 size_t m = 0;
540 size_t offset = Rfs(R).size();
541 for( vector<Connection>::const_iterator it = M(R).begin(); it != M(R).end(); it++ ) {
542 Rfactors.push_back( (*it).msg );
543 g2l[(int)-m-1] = m + offset;
544 m++;
545 }
546 // cavities
547 offset = Rfs(R).size() + M(R).size();
548 for( size_t c = 0; c < _cavitydists[R].size(); c++ ) {
549 Rfactors.push_back( _cavitydists[R][c] );
550 g2l[(int)-M(R).size() - c - 1] = offset + c;
551 }
552 // outgoing msgs
553 offset = Rfs(R).size() + M(R).size() + _cavitydists[R].size();
554 // _outOffset.push_back(-M(R).size() - _cavitydists[R].size() - 1);
555 size_t counter = 0;
556 if( props.verbose >= 3 )
557 cerr << "going to the loop!" << endl;
558 for( size_t c = 0; c < outM(R).size(); c++ ) {
559 bool isAvailable = false;
560 for( size_t I = 0; I < Rfactors.size(); I++ )
561 if( outM(R,c) << Rfactors[I].vars() ) {
562 isAvailable = true;
563 break;
564 }
565 if( !isAvailable ) {
566 Rfactors.push_back( Factor( outM(R)[c], (Real)1 ) );
567 g2l[(int)-M(R).size() - _cavitydists[R].size() - 1 - counter] = offset + counter;
568 counter++;
569 }
570 }
571
572 _CWs.push_back( Cobweb(g2l) );
573 InfAlg *alg = newInfAlg( name, FactorGraph(Rfactors), opts );
574 DAIAlgFG *dalg = static_cast<DAIAlgFG*> (alg);
575 _CWs[R].setInfAlg( dalg );
576 }
577 }
578
579
580 void GLC::initCWs(){
581 for( size_t R = 0; R < nrCWs(); R++ ) {
582 _CWs[R].initialize();
583 for( size_t m = 0; m < M(R).size(); m++ ) {
584 M(R)[m].msg.fill( (Real)1 );
585 M(R)[m].newmsg.fill( (Real)1 );
586 }
587 }
588 }
589
590
591 } // end of namespace dai