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     // calc iI
  75     for( size_t i=0; i < nrVars(); i++ ) {
  76         for( _nb_cit I = nb1(i).begin(); I != nb1(i).end(); I++ ) {
  77                 _iI_type _iI_entry;
  78                 _iI_entry.i = i;
  79                 _iI_entry.I = *I;
  80
  81                 _iI.push_back(_iI_entry);
  82             }
  83     }
  84
  85     // create pancakes
  86     _pancakes.resize(nrVars());
  87
  88     // create cavitydists
  89     for( size_t i=0; i < nrVars(); i++ )
  90         _cavitydists.push_back(Factor(delta(var(i))));
  91
  92     // create phis
  93     _phis.reserve(nr_edges());
  94     for( size_t iI = 0; iI < nr_edges(); iI++ ) {
  95         size_t i = edge(iI).first;
  96         size_t I = edge(iI).second;
  97         _phis.push_back( Factor( factor(I).vars() / var(i) ) );
  98     }
  99
 100     // create beliefs
 101     for( size_t i=0; i < nrVars(); i++ )
 102         _beliefs.push_back(Factor(var(i)));
 103 }
 104
 105
 106 string LC::identify() const {
 107     stringstream result (stringstream::out);
 108     result << Name << GetProperties();
 109     return result.str();
 110 }
 111
 112
 113 void LC::CalcBelief (size_t i) {
 114     _beliefs[i] = _pancakes[i].marginal(var(i));
 115 }
 116
 117
 118 double LC::CalcCavityDist (size_t i, const string &name, const Properties &opts) {
 119     Factor Bi;
 120     double maxdiff = 0;
 121
 122     if( Verbose() >= 2 )
 123         cout << "Initing cavity " << var(i) << "(" << delta(var(i)).size() << " vars, " << delta(var(i)).stateSpace() << " states)" << endl;
 124
 125     if( Cavity() == CavityType::UNIFORM )
 126         Bi = Factor(delta(var(i)));
 127     else {
 128         InfAlg *cav = newInfAlg( name, *this, opts );
 129         cav->makeCavity( var(i) );
 130
 131         if( Cavity() == CavityType::FULL )
 132             Bi = calcMarginal( *cav, cav->delta(var(i)), reInit() );
 133         else if( Cavity() == CavityType::PAIR )
 134             Bi = calcMarginal2ndO( *cav, cav->delta(var(i)), reInit() );
 135         else if( Cavity() == CavityType::PAIR2 ) {
 136             vector<Factor> pairbeliefs = calcPairBeliefsNew( *cav, cav->delta(var(i)), reInit() );
 137             for( size_t ij = 0; ij < pairbeliefs.size(); ij++ )
 138                 Bi *= pairbeliefs[ij];
 139         } else if( Cavity() == CavityType::PAIRINT ) {
 140             Bi = calcMarginal( *cav, cav->delta(var(i)), reInit() );
 141
 142             // Set interactions of order > 2 to zero
 143             size_t N = delta(var(i)).size();
 144             Real *p = &(*Bi.p().p().begin());
 145             x2x::p2logp (N, p);
 146             x2x::logp2w (N, p);
 147             x2x::fill (N, p, 2, 0.0);
 148             x2x::w2logp (N, p);
 149 //            x2x::logpnorm (N, p);
 150             x2x::logp2p (N, p);
 151         } else if( Cavity() == CavityType::PAIRCUM ) {
 152             Bi = calcMarginal( *cav, cav->delta(var(i)), reInit() );
 153
 154             // Set cumulants of order > 2 to zero
 155             size_t N = delta(var(i)).size();
 156             Real *p = &(*Bi.p().p().begin());
 157             x2x::p2m (N, p);
 158             x2x::m2c (N, p, N);
 159             x2x::fill (N, p, 2, 0.0);
 160             x2x::c2m (N, p, N);
 161             x2x::m2p (N, p);
 162         }
 163         maxdiff = cav->MaxDiff();
 164         delete cav;
 165     }
 166     Bi.normalize( _normtype );
 167     _cavitydists[i] = Bi;
 168
 169     return maxdiff;
 170 }
 171
 172
 173 double LC::InitCavityDists (const string &name, const Properties &opts) {
 174     clock_t tic = toc();
 175
 176     if( Verbose() >= 1 ) {
 177         cout << "LC::InitCavityDists:  ";
 178         if( Cavity() == CavityType::UNIFORM )
 179             cout << "Using uniform initial cavity distributions" << endl;
 180         else if( Cavity() == CavityType::FULL )
 181             cout << "Using full " << name << opts << "...";
 182         else if( Cavity() == CavityType::PAIR )
 183             cout << "Using pairwise " << name << opts << "...";
 184         else if( Cavity() == CavityType::PAIR2 )
 185             cout << "Using pairwise(new) " << name << opts << "...";
 186     }
 187
 188     double maxdiff = 0.0;
 189     for( size_t i = 0; i < nrVars(); i++ ) {
 190         double md = CalcCavityDist(i, name, opts);
 191         if( md > maxdiff )
 192             maxdiff = md;
 193     }
 194     init();
 195
 196     if( Verbose() >= 1 ) {
 197         cout << "used " << toc() - tic << " clocks." << endl;
 198     }
 199
 200     return maxdiff;
 201 }
 202
 203
 204 long LC::SetCavityDists (vector<Factor> &Q) {
 205     if( Verbose() >= 1 )
 206         cout << "LC::SetCavityDists:  Setting initial cavity distributions" << endl;
 207     if( Q.size() != nrVars() )
 208         return -1;
 209     for( size_t i = 0; i < nrVars(); i++ ) {
 210         if( _cavitydists[i].vars() != Q[i].vars() ) {
 211             return i+1;
 212         } else
 213             _cavitydists[i] = Q[i];
 214     }
 215     init();
 216     return 0;
 217 }
 218
 219
 220 void LC::init() {
 221     for( size_t iI = 0; iI < nr_edges(); iI++ ) {
 222         if( Updates() == UpdateType::SEQRND )
 223             _phis[iI].randomize();
 224         else
 225             _phis[iI].fill(1.0);
 226     }
 227     for( size_t i = 0; i < nrVars(); i++ ) {
 228         _pancakes[i] = _cavitydists[i];
 229
 230         for( _nb_cit I = nb1(i).begin(); I != nb1(i).end(); I++ ) {
 231             _pancakes[i] *= factor(*I);
 232             if( Updates() == UpdateType::SEQRND )
 233               _pancakes[i] *= _phis[VV2E(i,*I)];
 234         }
 235
 236         _pancakes[i].normalize( _normtype );
 237
 238         CalcBelief(i);
 239     }
 240 }
 241
 242
 243 Factor LC::NewPancake (size_t iI, bool & hasNaNs) {
 244     size_t i = _iI[iI].i;
 245     size_t I = _iI[iI].I;
 246     iI = VV2E(i, I);
 247
 248     Factor piet = _pancakes[i];
 249
 250     // recalculate _pancake[i]
 251     VarSet Ivars = factor(I).vars();
 252     Factor A_I;
 253     for( VarSet::const_iterator k = Ivars.begin(); k != Ivars.end(); k++ )
 254         if( var(i) != *k )
 255             A_I *= (_pancakes[findVar(*k)] * factor(I).inverse()).part_sum( Ivars / var(i) );
 256     if( Ivars.size() > 1 )
 257         A_I ^= (1.0 / (Ivars.size() - 1));
 258     Factor A_Ii = (_pancakes[i] * factor(I).inverse() * _phis[iI].inverse()).part_sum( Ivars / var(i) );
 259     Factor quot = A_I.divided_by(A_Ii);
 260
 261     piet *= quot.divided_by( _phis[iI] ).normalized( _normtype );
 262     _phis[iI] = quot.normalized( _normtype );
 263
 264     piet.normalize( _normtype );
 265
 266     if( piet.hasNaNs() ) {
 267         cout << "LC::NewPancake(" << iI << "):  has NaNs!" << endl;
 268         hasNaNs = true;
 269     }
 270
 271     return piet;
 272 }
 273
 274
 275 double LC::run() {
 276     if( Verbose() >= 1 )
 277         cout << "Starting " << identify() << "...";
 278     if( Verbose() >= 2 )
 279         cout << endl;
 280
 281     clock_t tic = toc();
 282     Diffs diffs(nrVars(), 1.0);
 283
 284     updateMaxDiff( InitCavityDists(GetPropertyAs<string>("cavainame"), GetPropertyAs<Properties>("cavaiopts")) );
 285
 286     vector<Factor> old_beliefs;
 287     for(size_t i=0; i < nrVars(); i++ )
 288         old_beliefs.push_back(belief(i));
 289
 290     bool hasNaNs = false;
 291     for( size_t i=0; i < nrVars(); i++ )
 292         if( _pancakes[i].hasNaNs() ) {
 293             hasNaNs = true;
 294             break;
 295         }
 296     if( hasNaNs ) {
 297         cout << "LC::run:  initial _pancakes has NaNs!" << endl;
 298         return NAN;
 299     }
 300
 301     vector<long> update_seq(nr_iI(),0);
 302     for( size_t k=0; k < nr_iI(); k++ )
 303         update_seq[k] = k;
 304
 305     size_t iter=0;
 306
 307     // do several passes over the network until maximum number of iterations has
 308     // been reached or until the maximum belief difference is smaller than tolerance
 309     for( iter=0; iter < MaxIter() && diffs.max() > Tol(); iter++ ) {
 310         // Sequential updates
 311         if( Updates() == UpdateType::SEQRND )
 312             random_shuffle( update_seq.begin(), update_seq.end() );
 313
 314         for( size_t t=0; t < nr_iI(); t++ ) {
 315             long iI = update_seq[t];
 316             long i = _iI[iI].i;
 317             _pancakes[i] = NewPancake(iI, hasNaNs);
 318             if( hasNaNs )
 319                 return NAN;
 320             CalcBelief(i);
 321         }
 322
 323         // compare new beliefs with old ones
 324         for(size_t i=0; i < nrVars(); i++ ) {
 325             diffs.push( dist( belief(i), old_beliefs[i], Prob::DISTLINF ) );
 326             old_beliefs[i] = belief(i);
 327         }
 328
 329         if( Verbose() >= 3 )
 330             cout << "LC::run:  maxdiff " << diffs.max() << " after " << iter+1 << " passes" << endl;
 331     }
 332
 333     updateMaxDiff( diffs.max() );
 334
 335     if( Verbose() >= 1 ) {
 336         if( diffs.max() > Tol() ) {
 337             if( Verbose() == 1 )
 338                 cout << endl;
 339                 cout << "LC::run:  WARNING: not converged within " << MaxIter() << " passes (" << toc() - tic << " clocks)...final maxdiff:" << diffs.max() << endl;
 340         } else {
 341             if( Verbose() >= 2 )
 342                 cout << "LC::run:  ";
 343                 cout << "converged in " << iter << " passes (" << toc() - tic << " clocks)." << endl;
 344         }
 345     }
 346
 347     return diffs.max();
 348 }
 349
 350
 351 } // end of namespace dai