Fixed GLC::CalcFactorBelief (now uses geometric averaging) and GLC::belief
[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 Factor tmp( ns, 1.0 );
280 size_t counter = 0;
281 for( size_t R = 0; R < nrCWs(); R++ )
282 if( ns << inds2vars( INRs(R).elements() ) ) {
283 tmp *= CW(R).belief(ns);
284 counter++;
285 }
286 if( counter > 0 )
287 _factorBeliefs[ns] = tmp ^ (1.0 / counter);
288 }
289
290
291 Factor GLC::belief( const VarSet &ns ) const {
292 if( ns.size() == 0 )
293 return Factor();
294 else if( ns.size() == 1 )
295 return beliefV( findVar( *(ns.begin()) ) );
296 else {
297 map<VarSet, Factor>::const_iterator it = _factorBeliefs.find(ns);
298 if( it != _factorBeliefs.end() )
299 return it->second;
300 for( map<VarSet, Factor>::const_iterator it = _factorBeliefs.begin(); it != _factorBeliefs.end(); it++ )
301 if( ns << it->second.vars() )
302 return it->second.marginal( ns );
303 }
304 DAI_THROW(BELIEF_NOT_AVAILABLE);
305 return Factor();
306 }
307
308
309 vector<Factor> GLC::beliefs() const {
310 vector<Factor> result = _beliefs;
311 for( map<VarSet, Factor>::const_iterator it = _factorBeliefs.begin(); it != _factorBeliefs.end(); it++ )
312 result.push_back( it->second );
313 return result;
314 }
315
316
317 Real GLC::CalcCavityDist( size_t R, const std::string &name, const PropertySet& opts ) {
318 vector<Factor> pairbeliefs;
319 Real maxdiff = 0;
320
321 if( props.cavity != Properties::CavityType::UNIFORM ) {
322 InfAlg *cav = newInfAlg( name, *this, opts );
323 cav->makeRegionCavity( Rfs(R).elements() );
324
325 if( props.cavity == Properties::CavityType::FULL )
326 pairbeliefs.push_back( calcMarginal( *cav, inds2vars(EXRs(R).elements()), props.reinit ) );
327 if( props.cavity == Properties::CavityType::PAIR )
328 pairbeliefs = calcPairBeliefs( *cav, inds2vars(EXRs(R).elements()), props.reinit, false );
329 else if( props.cavity == Properties::CavityType::PAIR2 )
330 pairbeliefs = calcPairBeliefs( *cav, inds2vars(EXRs(R).elements()), props.reinit, true );
331 maxdiff = cav->maxDiff();
332 delete cav;
333 }
334 if( props.verbose >= 3 )
335 cerr << "R:" << R << "cavity of size " << pairbeliefs.size() << endl;
336 _cavitydists[R] = pairbeliefs;
337 return maxdiff;
338 }
339
340
341 Real GLC::InitCavityDists( const std::string &name, const PropertySet &opts ) {
342 double tic = toc();
343
344 if( props.verbose >= 2 ) {
345 cerr << this->name() << "::InitCavityDists: ";
346 if( props.cavity == Properties::CavityType::UNIFORM )
347 cerr << "Using uniform initial cavity distributions" << endl;
348 else if( props.cavity == Properties::CavityType::FULL )
349 cerr << "Using full " << name << opts << "...";
350 else if( props.cavity == Properties::CavityType::PAIR )
351 cerr << "Using pairwise " << name << opts << "...";
352 else if( props.cavity == Properties::CavityType::PAIR2 )
353 cerr << "Using pairwise(new) " << name << opts << "...";
354 }
355
356 Real maxdiff = 0.0;
357 for( size_t R = 0; R < nrCWs(); R++ ) {
358 Real md = CalcCavityDist(R, name, opts);
359 if( md > maxdiff )
360 maxdiff = md;
361 }
362 if( props.verbose >= 2 )
363 cerr << this->name() << "::InitCavityDists used " << toc() - tic << " seconds." << endl;
364 return maxdiff;
365 }
366
367
368 void GLC::NewPancake( size_t R, size_t _R2 ) {
369 Factor newmsg;
370 // calculate the new msg
371 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());
372 newmsg.normalize();
373 // update the cobweb-graph with this new msg
374 M(R, _R2).msg = newmsg;
375 // update the corresponding factor in the cobweb region so to be used in future marginalizations (-_R2 - 1) is the index (see cobweb class documents)
376 CW(R).updateFactor( -_R2 - 1 ,newmsg);
377 }
378
379
380 void GLC::OVNewPancake( size_t R ) {
381 Factor newmsg, allmsgs, marg;
382 for( size_t _R2 = 0; _R2 < M(R).size(); _R2++ ) { // for all _R2 that send messages to R
383 // calculate the message
384 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;
385 newmsg.normalize();
386 // 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
387 allmsgs = newmsg;
388 for( size_t sub = 0; sub < M(R,_R2).subregions.size(); sub++ ) { //for all descendants of rho of (\ominus r_{R, R2})
389 Real c = (Real)cn(R)[M(R, _R2).subregions[sub]].first; // counting number for rho
390 size_t nrParents = cn(R)[M(R, _R2).subregions[sub]].second.size(); // number of rho's parents
391 marg = newmsg.marginal(M(R, _R2).subregions[sub], false); // marginal of the message from R2 to R over rho
392 allmsgs *= (marg^(c / nrParents)); // discount the double counting of rho from the message
393 for( size_t pind = 0; pind < nrParents; pind++ ) // this operation corresponds to the upward pass
394 if( cn(R)[M(R, _R2).subregions[sub]].second[pind] != _R2 )
395 M(R, cn(R)[M(R, _R2).subregions[sub]].second[pind]).newmsg *= marg^(-c / (nrParents * (nrParents - (Real)1)));
396 }
397 // set the message after taking double-countings into account
398 M(R, _R2).msg = allmsgs.normalized();
399 }
400
401 for( size_t _R2 = 0; _R2 < M(R).size(); _R2++ ) {
402 M(R, _R2).newmsg *= M(R, _R2).msg;
403 M(R, _R2).newmsg.normalize();
404 // update the corresponding factor in the cobweb region: to be used in future marginalizations
405 CW(R).updateFactor( -_R2 - 1 , M(R, _R2).msg, false);
406 M(R, _R2).msg = M(R, _R2).newmsg;
407 M(R, _R2).newmsg.fill( (Real)1 );
408 }
409 }
410
411
412 Real GLC::run() {
413 if( props.verbose >= 2 )
414 cerr << "Starting " << identify() << "...";
415 if( props.verbose >= 2 )
416 cerr << endl;
417
418 double tic = toc();
419
420 if( props.verbose >= 3 )
421 cerr << "Initializing older beliefs" << endl;
422 vector<vector<Factor> > oldBeliefsM;
423 vector<Factor> oldBeliefsP;
424 if( isPartition ) {
425 oldBeliefsM.reserve( nrCWs() );
426 for( size_t R = 0; R < nrCWs(); R++ ) {
427 oldBeliefsM.push_back( vector<Factor>() );
428 oldBeliefsM[R].reserve( M(R).size() );
429 for( size_t m = 0; m < M(R).size(); m++ )
430 oldBeliefsM[R].push_back( M(R,m).msg );
431 }
432 } else {
433 oldBeliefsP.reserve( nrVars() );
434 for( size_t i = 0; i < nrVars(); i++ ) {
435 CalcBelief( i, true );
436 oldBeliefsP.push_back( _beliefs[i] );
437 }
438 }
439 if( props.verbose >= 3 )
440 cerr << "Setting the update sequence" <<endl;
441
442 size_t nredges = 0;
443 vector<Edge> update_seq;
444 update_seq.reserve( nredges );
445 for( size_t R = 0; R < nrCWs(); ++R )
446 for( size_t R2 = 0; R2 < M(R).size(); R2++ ) {
447 update_seq.push_back( Edge( R, R2 ) );
448 nredges++;
449 }
450
451 vector<size_t> update_seq_node;
452 update_seq_node.reserve( nrCWs() );
453 for( size_t R = 0; R < nrCWs(); ++R )
454 update_seq_node.push_back( R );
455 if( props.verbose >= 3 )
456 cerr << "Going into the main loop" << endl;
457
458 Real maxDiff = INFINITY;
459 for( _iters = 0; _iters < props.maxiter && maxDiff > props.tol && (toc() - tic) < props.maxtime; _iters++ ) {
460 if( props.updates == Properties::UpdateType::SEQRND ) {
461 random_shuffle( update_seq.begin(), update_seq.end() );
462 random_shuffle( update_seq_node.begin(), update_seq_node.end() );
463 }
464 if( isPartition ) {
465 for( size_t t = 0; t < nredges; t++ ) {
466 size_t R = update_seq[t].first;
467 size_t _R2 = update_seq[t].second;
468 NewPancake( R, _R2);
469 }
470 } else {
471 for( size_t t = 0; t < nrCWs(); t++ )
472 OVNewPancake( update_seq_node[t] );
473 }
474
475 if( !isPartition )
476 for( size_t i = 0; i < nrVars(); i++ )
477 CalcBelief( i , true );
478
479 maxDiff = -INFINITY;
480 if( isPartition ) {
481 for( size_t R = 0; R < nrCWs(); R++ ) {
482 for( size_t m = 0; m < M(R).size(); m++ ) {
483 maxDiff = std::max( maxDiff, dist( M(R, m).msg, oldBeliefsM[R][m], DISTLINF ) );
484 oldBeliefsM[R][m] = M(R,m).msg;
485 }
486 }
487 } else {
488 for( size_t i = 0; i < nrVars(); i++ ) {
489 maxDiff = std::max( maxDiff, dist( _beliefs[i], oldBeliefsP[i], DISTLINF ) );
490 oldBeliefsP[i] = _beliefs[i];
491 }
492 }
493
494 if( props.verbose >= 2 )
495 cerr << this->name() << "::run: maxdiff " << maxDiff << " after " << _iters+1 << " passes" << endl;
496 }
497
498 if( isPartition )
499 for( size_t i = 0; i < nrVars(); i++ )
500 CalcBelief( i, true );
501
502 if( maxDiff > _maxdiff )
503 _maxdiff = maxDiff;
504
505 if( props.verbose >= 2 ) {
506 if( maxDiff > props.tol ) {
507 if( props.verbose == 2 )
508 cerr << endl;
509 cerr << this->name() << "::run: WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " seconds)...final maxdiff:" << maxDiff << endl;
510 } else {
511 if( props.verbose >= 2 )
512 cerr << this->name() << "::run: ";
513 cerr << "converged in " << _iters << " passes (" << toc() - tic << " seconds)." << endl;
514 }
515 }
516
517 for( size_t I = 0; I < nrFactors(); I++ )
518 CalcFactorBelief( I );
519
520 return maxDiff;
521 }
522
523
524 void GLC::setCWs( const std::string &name, const PropertySet &opts ) {
525 // needs to generate a vector of relevant factors, and a dictionary for mapping between local and global factors
526 _CWs.clear();
527 _CWs.reserve( _INRs.size() );
528 _outOffset.clear();
529 _outOffset.reserve( _INRs.size() );
530 // main factors
531 if( props.verbose >= 3 )
532 cerr << "Setting CWs..." << endl;
533
534 for( size_t R = 0; R < nrCWs(); R++ ) { // for each region
535 map<int, size_t> g2l;
536 vector<Factor> Rfactors; // contains all factors
537 size_t potsize = Rfs(R).size() + M(R).size() + _cavitydists[R].size();
538 Rfactors.reserve( potsize );
539 size_t _I = 0;
540 // adding main factors
541 for( SmallSet<size_t>::const_iterator it = Rfs(R).begin(); it != Rfs(R).end(); it++ ) {
542 Rfactors.push_back( factor(*it) );
543 g2l[(int)(*it)] = _I;
544 _I++;
545 }
546
547 // adding factors for incoming messages
548 size_t m = 0;
549 size_t offset = Rfs(R).size();
550 for( vector<Connection>::const_iterator it = M(R).begin(); it != M(R).end(); it++ ) {
551 Rfactors.push_back( (*it).msg );
552 g2l[(int)-m-1] = m + offset;
553 m++;
554 }
555 // cavities
556 offset = Rfs(R).size() + M(R).size();
557 for( size_t c = 0; c < _cavitydists[R].size(); c++ ) {
558 Rfactors.push_back( _cavitydists[R][c] );
559 g2l[(int)-M(R).size() - c - 1] = offset + c;
560 }
561 // outgoing msgs
562 offset = Rfs(R).size() + M(R).size() + _cavitydists[R].size();
563 // _outOffset.push_back(-M(R).size() - _cavitydists[R].size() - 1);
564 size_t counter = 0;
565 if( props.verbose >= 3 )
566 cerr << "going to the loop!" << endl;
567 for( size_t c = 0; c < outM(R).size(); c++ ) {
568 bool isAvailable = false;
569 for( size_t I = 0; I < Rfactors.size(); I++ )
570 if( outM(R,c) << Rfactors[I].vars() ) {
571 isAvailable = true;
572 break;
573 }
574 if( !isAvailable ) {
575 Rfactors.push_back( Factor( outM(R)[c], (Real)1 ) );
576 g2l[(int)-M(R).size() - _cavitydists[R].size() - 1 - counter] = offset + counter;
577 counter++;
578 }
579 }
580
581 _CWs.push_back( Cobweb(g2l) );
582 InfAlg *alg = newInfAlg( name, FactorGraph(Rfactors), opts );
583 DAIAlgFG *dalg = static_cast<DAIAlgFG*> (alg);
584 _CWs[R].setInfAlg( dalg );
585 }
586 }
587
588
589 void GLC::initCWs(){
590 for( size_t R = 0; R < nrCWs(); R++ ) {
591 _CWs[R].initialize();
592 for( size_t m = 0; m < M(R).size(); m++ ) {
593 M(R)[m].msg.fill( (Real)1 );
594 M(R)[m].newmsg.fill( (Real)1 );
595 }
596 }
597 }
598
599
600 } // end of namespace dai