src/lc.cpp

   1 /*  Copyright (C) 2006-2008  Joris Mooij  [j dot mooij at science dot ru dot nl]
   2     Radboud University Nijmegen, The Netherlands
   3
   4     This file is part of libDAI.
   5
   6     libDAI is free software; you can redistribute it and/or modify
   7     it under the terms of the GNU General Public License as published by
   8     the Free Software Foundation; either version 2 of the License, or
   9     (at your option) any later version.
  10
  11     libDAI is distributed in the hope that it will be useful,
  12     but WITHOUT ANY WARRANTY; without even the implied warranty of
  13     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14     GNU General Public License for more details.
  15
  16     You should have received a copy of the GNU General Public License
  17     along with libDAI; if not, write to the Free Software
  18     Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  19 */
  20
  21
  22 #include <iostream>
  23 #include <algorithm>
  24 #include <map>
  25 #include <set>
  26 #include <dai/lc.h>
  27 #include <dai/diffs.h>
  28 #include <dai/util.h>
  29 #include <dai/alldai.h>
  30 #include <dai/x2x.h>
  31
  32
  33 namespace dai {
  34
  35
  36 using namespace std;
  37
  38
  39 const char *LC::Name = "LC";
  40
  41
  42 bool LC::checkProperties() {
  43     if( !HasProperty("cavity") )
  44         return false;
  45     if( !HasProperty("updates") )
  46         return false;
  47     if( !HasProperty("tol") )
  48         return false;
  49     if (!HasProperty("maxiter") )
  50         return false;
  51     if (!HasProperty("verbose") )
  52         return false;
  53
  54     ConvertPropertyTo<CavityType>("cavity");
  55     ConvertPropertyTo<UpdateType>("updates");
  56     ConvertPropertyTo<double>("tol");
  57     ConvertPropertyTo<size_t>("maxiter");
  58     ConvertPropertyTo<size_t>("verbose");
  59
  60     if (HasProperty("cavainame") )
  61         ConvertPropertyTo<string>("cavainame");
  62     if (HasProperty("cavaiopts") )
  63         ConvertPropertyTo<Properties>("cavaiopts");
  64     if( HasProperty("reinit") )
  65         ConvertPropertyTo<bool>("reinit");
  66
  67     return true;
  68 }
  69
  70
  71 LC::LC(const FactorGraph & fg, const Properties &opts) : DAIAlgFG(fg, opts) {
  72     assert( checkProperties() );
  73
  74     // create pancakes
  75     _pancakes.resize(nrVars());
  76
  77     // create cavitydists
  78     for( size_t i=0; i < nrVars(); i++ )
  79         _cavitydists.push_back(Factor(delta(i)));
  80
  81     // create phis
  82     _phis.reserve( nrVars() );
  83     for( size_t i = 0; i < nrVars(); i++ ) {
  84         _phis.push_back( vector<Factor>() );
  85         _phis[i].reserve( nbV(i).size() );
  86         foreach( const Neighbor &I, nbV(i) )
  87             _phis[i].push_back( Factor( factor(I).vars() / var(i) ) );
  88     }
  89
  90     // create beliefs
  91     for( size_t i=0; i < nrVars(); i++ )
  92         _beliefs.push_back(Factor(var(i)));
  93 }
  94
  95
  96 string LC::identify() const {
  97     stringstream result (stringstream::out);
  98     result << Name << GetProperties();
  99     return result.str();
 100 }
 101
 102
 103 void LC::CalcBelief (size_t i) {
 104     _beliefs[i] = _pancakes[i].marginal(var(i));
 105 }
 106
 107
 108 double LC::CalcCavityDist (size_t i, const std::string &name, const Properties &opts) {
 109     Factor Bi;
 110     double maxdiff = 0;
 111
 112     if( Verbose() >= 2 )
 113         cout << "Initing cavity " << var(i) << "(" << delta(i).size() << " vars, " << delta(i).states() << " states)" << endl;
 114
 115     if( Cavity() == CavityType::UNIFORM )
 116         Bi = Factor(delta(i));
 117     else {
 118         InfAlg *cav = newInfAlg( name, *this, opts );
 119         cav->makeCavity( i );
 120
 121         if( Cavity() == CavityType::FULL )
 122             Bi = calcMarginal( *cav, cav->fg().delta(i), reInit() );
 123         else if( Cavity() == CavityType::PAIR )
 124             Bi = calcMarginal2ndO( *cav, cav->fg().delta(i), reInit() );
 125         else if( Cavity() == CavityType::PAIR2 ) {
 126             vector<Factor> pairbeliefs = calcPairBeliefsNew( *cav, cav->fg().delta(i), reInit() );
 127             for( size_t ij = 0; ij < pairbeliefs.size(); ij++ )
 128                 Bi *= pairbeliefs[ij];
 129         } else if( Cavity() == CavityType::PAIRINT ) {
 130             Bi = calcMarginal( *cav, cav->fg().delta(i), reInit() );
 131
 132             // Set interactions of order > 2 to zero
 133             size_t N = delta(i).size();
 134             Real *p = &(*Bi.p().p().begin());
 135             x2x::p2logp (N, p);
 136             x2x::logp2w (N, p);
 137             x2x::fill (N, p, 2, 0.0);
 138             x2x::w2logp (N, p);
 139 //            x2x::logpnorm (N, p);
 140             x2x::logp2p (N, p);
 141         } else if( Cavity() == CavityType::PAIRCUM ) {
 142             Bi = calcMarginal( *cav, cav->fg().delta(i), reInit() );
 143
 144             // Set cumulants of order > 2 to zero
 145             size_t N = delta(i).size();
 146             Real *p = &(*Bi.p().p().begin());
 147             x2x::p2m (N, p);
 148             x2x::m2c (N, p, N);
 149             x2x::fill (N, p, 2, 0.0);
 150             x2x::c2m (N, p, N);
 151             x2x::m2p (N, p);
 152         }
 153         maxdiff = cav->MaxDiff();
 154         delete cav;
 155     }
 156     Bi.normalize( _normtype );
 157     _cavitydists[i] = Bi;
 158
 159     return maxdiff;
 160 }
 161
 162
 163 double LC::InitCavityDists (const std::string &name, const Properties &opts) {
 164     double tic = toc();
 165
 166     if( Verbose() >= 1 ) {
 167         cout << "LC::InitCavityDists:  ";
 168         if( Cavity() == CavityType::UNIFORM )
 169             cout << "Using uniform initial cavity distributions" << endl;
 170         else if( Cavity() == CavityType::FULL )
 171             cout << "Using full " << name << opts << "...";
 172         else if( Cavity() == CavityType::PAIR )
 173             cout << "Using pairwise " << name << opts << "...";
 174         else if( Cavity() == CavityType::PAIR2 )
 175             cout << "Using pairwise(new) " << name << opts << "...";
 176     }
 177
 178     double maxdiff = 0.0;
 179     for( size_t i = 0; i < nrVars(); i++ ) {
 180         double md = CalcCavityDist(i, name, opts);
 181         if( md > maxdiff )
 182             maxdiff = md;
 183     }
 184     init();
 185
 186     if( Verbose() >= 1 ) {
 187         cout << "used " << toc() - tic << " clocks." << endl;
 188     }
 189
 190     return maxdiff;
 191 }
 192
 193
 194 long LC::SetCavityDists( std::vector<Factor> &Q ) {
 195     if( Verbose() >= 1 )
 196         cout << "LC::SetCavityDists:  Setting initial cavity distributions" << endl;
 197     if( Q.size() != nrVars() )
 198         return -1;
 199     for( size_t i = 0; i < nrVars(); i++ ) {
 200         if( _cavitydists[i].vars() != Q[i].vars() ) {
 201             return i+1;
 202         } else
 203             _cavitydists[i] = Q[i];
 204     }
 205     init();
 206     return 0;
 207 }
 208
 209
 210 void LC::init() {
 211     for( size_t i = 0; i < nrVars(); ++i )
 212         foreach( const Neighbor &I, nbV(i) )
 213             if( Updates() == UpdateType::SEQRND )
 214                 _phis[i][I.iter].randomize();
 215             else
 216                 _phis[i][I.iter].fill(1.0);
 217     for( size_t i = 0; i < nrVars(); i++ ) {
 218         _pancakes[i] = _cavitydists[i];
 219
 220         foreach( const Neighbor &I, nbV(i) ) {
 221             _pancakes[i] *= factor(I);
 222             if( Updates() == UpdateType::SEQRND )
 223               _pancakes[i] *= _phis[i][I.iter];
 224         }
 225
 226         _pancakes[i].normalize( _normtype );
 227
 228         CalcBelief(i);
 229     }
 230 }
 231
 232
 233 Factor LC::NewPancake (size_t i, size_t _I, bool & hasNaNs) {
 234     size_t I = nbV(i)[_I];
 235     Factor piet = _pancakes[i];
 236
 237     // recalculate _pancake[i]
 238     VarSet Ivars = factor(I).vars();
 239     Factor A_I;
 240     for( VarSet::const_iterator k = Ivars.begin(); k != Ivars.end(); k++ )
 241         if( var(i) != *k )
 242             A_I *= (_pancakes[findVar(*k)] * factor(I).inverse()).part_sum( Ivars / var(i) );
 243     if( Ivars.size() > 1 )
 244         A_I ^= (1.0 / (Ivars.size() - 1));
 245     Factor A_Ii = (_pancakes[i] * factor(I).inverse() * _phis[i][_I].inverse()).part_sum( Ivars / var(i) );
 246     Factor quot = A_I.divided_by(A_Ii);
 247
 248     piet *= quot.divided_by( _phis[i][_I] ).normalized( _normtype );
 249     _phis[i][_I] = quot.normalized( _normtype );
 250
 251     piet.normalize( _normtype );
 252
 253     if( piet.hasNaNs() ) {
 254         cout << "LC::NewPancake(" << i << ", " << _I << "):  has NaNs!" << endl;
 255         hasNaNs = true;
 256     }
 257
 258     return piet;
 259 }
 260
 261
 262 double LC::run() {
 263     if( Verbose() >= 1 )
 264         cout << "Starting " << identify() << "...";
 265     if( Verbose() >= 2 )
 266         cout << endl;
 267
 268     double tic = toc();
 269     Diffs diffs(nrVars(), 1.0);
 270
 271     updateMaxDiff( InitCavityDists(GetPropertyAs<string>("cavainame"), GetPropertyAs<Properties>("cavaiopts")) );
 272
 273     vector<Factor> old_beliefs;
 274     for(size_t i=0; i < nrVars(); i++ )
 275         old_beliefs.push_back(belief(i));
 276
 277     bool hasNaNs = false;
 278     for( size_t i=0; i < nrVars(); i++ )
 279         if( _pancakes[i].hasNaNs() ) {
 280             hasNaNs = true;
 281             break;
 282         }
 283     if( hasNaNs ) {
 284         cout << "LC::run:  initial _pancakes has NaNs!" << endl;
 285         return -1.0;
 286     }
 287
 288     size_t nredges = nrEdges();
 289     vector<Edge> update_seq;
 290     update_seq.reserve( nredges );
 291     for( size_t i = 0; i < nrVars(); ++i )
 292         foreach( const Neighbor &I, nbV(i) )
 293             update_seq.push_back( Edge( i, I.iter ) );
 294
 295     size_t iter = 0;
 296
 297     // do several passes over the network until maximum number of iterations has
 298     // been reached or until the maximum belief difference is smaller than tolerance
 299     for( iter=0; iter < MaxIter() && diffs.maxDiff() > Tol(); iter++ ) {
 300         // Sequential updates
 301         if( Updates() == UpdateType::SEQRND )
 302             random_shuffle( update_seq.begin(), update_seq.end() );
 303
 304         for( size_t t=0; t < nredges; t++ ) {
 305             size_t i = update_seq[t].first;
 306             size_t _I = update_seq[t].second;
 307             _pancakes[i] = NewPancake( i, _I, hasNaNs);
 308             if( hasNaNs )
 309                 return -1.0;
 310             CalcBelief( i );
 311         }
 312
 313         // compare new beliefs with old ones
 314         for(size_t i=0; i < nrVars(); i++ ) {
 315             diffs.push( dist( belief(i), old_beliefs[i], Prob::DISTLINF ) );
 316             old_beliefs[i] = belief(i);
 317         }
 318
 319         if( Verbose() >= 3 )
 320             cout << "LC::run:  maxdiff " << diffs.maxDiff() << " after " << iter+1 << " passes" << endl;
 321     }
 322
 323     updateMaxDiff( diffs.maxDiff() );
 324
 325     if( Verbose() >= 1 ) {
 326         if( diffs.maxDiff() > Tol() ) {
 327             if( Verbose() == 1 )
 328                 cout << endl;
 329                 cout << "LC::run:  WARNING: not converged within " << MaxIter() << " passes (" << toc() - tic << " clocks)...final maxdiff:" << diffs.maxDiff() << endl;
 330         } else {
 331             if( Verbose() >= 2 )
 332                 cout << "LC::run:  ";
 333                 cout << "converged in " << iter << " passes (" << toc() - tic << " clocks)." << endl;
 334         }
 335     }
 336
 337     return diffs.maxDiff();
 338 }
 339
 340
 341 } // end of namespace dai