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