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