src/mr.cpp

   1 /*  This file is part of libDAI - http://www.libdai.org/
   2  *
   3  *  Copyright (c) 2006-2011, 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_MR
  11
  12
  13 #include <cstdio>
  14 #include <ctime>
  15 #include <cmath>
  16 #include <cstdlib>
  17 #include <dai/mr.h>
  18 #include <dai/bp.h>
  19 #include <dai/jtree.h>
  20 #include <dai/util.h>
  21 #include <dai/bbp.h>
  22
  23
  24 namespace dai {
  25
  26
  27 using namespace std;
  28
  29
  30 void MR::setProperties( const PropertySet &opts ) {
  31     DAI_ASSERT( opts.hasKey("tol") );
  32     DAI_ASSERT( opts.hasKey("updates") );
  33     DAI_ASSERT( opts.hasKey("inits") );
  34
  35     props.tol = opts.getStringAs<Real>("tol");
  36     props.updates = opts.getStringAs<Properties::UpdateType>("updates");
  37     props.inits = opts.getStringAs<Properties::InitType>("inits");
  38     if( opts.hasKey("verbose") )
  39         props.verbose = opts.getStringAs<size_t>("verbose");
  40     else
  41         props.verbose = 0;
  42 }
  43
  44
  45 PropertySet MR::getProperties() const {
  46     PropertySet opts;
  47     opts.set( "tol", props.tol );
  48     opts.set( "verbose", props.verbose );
  49     opts.set( "updates", props.updates );
  50     opts.set( "inits", props.inits );
  51     return opts;
  52 }
  53
  54
  55 string MR::printProperties() const {
  56     stringstream s( stringstream::out );
  57     s << "[";
  58     s << "tol=" << props.tol << ",";
  59     s << "verbose=" << props.verbose << ",";
  60     s << "updates=" << props.updates << ",";
  61     s << "inits=" << props.inits << "]";
  62     return s.str();
  63 }
  64
  65
  66 Real MR::T(size_t i, sub_nb A) {
  67     sub_nb _nbi_min_A(G.nb(i).size());
  68     _nbi_min_A.set();
  69     _nbi_min_A &= ~A;
  70
  71     Real res = theta[i];
  72     for( size_t _j = 0; _j < _nbi_min_A.size(); _j++ )
  73         if( _nbi_min_A.test(_j) )
  74             res += atanh(tJ[i][_j] * M[i][_j]);
  75     return tanh(res);
  76 }
  77
  78
  79 Real MR::T(size_t i, size_t _j) {
  80     sub_nb j(G.nb(i).size());
  81     j.set(_j);
  82     return T(i,j);
  83 }
  84
  85
  86 Real MR::Omega(size_t i, size_t _j, size_t _l) {
  87     sub_nb jl(G.nb(i).size());
  88     jl.set(_j);
  89     jl.set(_l);
  90     Real Tijl = T(i,jl);
  91     return Tijl / (1.0 + tJ[i][_l] * M[i][_l] * Tijl);
  92 }
  93
  94
  95 Real MR::Gamma(size_t i, size_t _j, size_t _l1, size_t _l2) {
  96     sub_nb jll(G.nb(i).size());
  97     jll.set(_j);
  98     Real Tij = T(i,jll);
  99     jll.set(_l1);
 100     jll.set(_l2);
 101     Real Tijll = T(i,jll);
 102
 103     return (Tijll - Tij) / (1.0 + tJ[i][_l1] * tJ[i][_l2] * M[i][_l1] * M[i][_l2] + tJ[i][_l1] * M[i][_l1] * Tijll + tJ[i][_l2] * M[i][_l2] * Tijll);
 104 }
 105
 106
 107 Real MR::Gamma(size_t i, size_t _l1, size_t _l2) {
 108     sub_nb ll(G.nb(i).size());
 109     Real Ti = T(i,ll);
 110     ll.set(_l1);
 111     ll.set(_l2);
 112     Real Till = T(i,ll);
 113
 114     return (Till - Ti) / (1.0 + tJ[i][_l1] * tJ[i][_l2] * M[i][_l1] * M[i][_l2] + tJ[i][_l1] * M[i][_l1] * Till + tJ[i][_l2] * M[i][_l2] * Till);
 115 }
 116
 117
 118 Real MR::_tJ(size_t i, sub_nb A) {
 119     sub_nb::size_type _j = A.find_first();
 120     if( _j == sub_nb::npos )
 121         return 1.0;
 122     else
 123         return tJ[i][_j] * _tJ(i, A.reset(_j));
 124 }
 125
 126
 127 Real MR::appM(size_t i, sub_nb A) {
 128     sub_nb::size_type _j = A.find_first();
 129     if( _j == sub_nb::npos )
 130         return 1.0;
 131     else {
 132         sub_nb A_j(A); A_j.reset(_j);
 133
 134         Real result = M[i][_j] * appM(i, A_j);
 135         for( size_t _k = 0; _k < A_j.size(); _k++ )
 136             if( A_j.test(_k) ) {
 137                 sub_nb A_jk(A_j); A_jk.reset(_k);
 138                 result += cors[i][_j][_k] * appM(i,A_jk);
 139             }
 140
 141         return result;
 142     }
 143 }
 144
 145
 146 void MR::sum_subs(size_t j, sub_nb A, Real *sum_even, Real *sum_odd) {
 147     *sum_even = 0.0;
 148     *sum_odd = 0.0;
 149
 150     sub_nb B(A.size());
 151     do {
 152         if( B.count() % 2 )
 153             *sum_odd += _tJ(j,B) * appM(j,B);
 154         else
 155             *sum_even += _tJ(j,B) * appM(j,B);
 156
 157         // calc next subset B
 158         size_t bit = 0;
 159         for( ; bit < A.size(); bit++ )
 160             if( A.test(bit) ) {
 161                 if( B.test(bit) )
 162                     B.reset(bit);
 163                 else {
 164                     B.set(bit);
 165                     break;
 166                 }
 167             }
 168     } while (!B.none());
 169 }
 170
 171
 172 void MR::propagateCavityFields() {
 173     Real sum_even, sum_odd;
 174     Real maxdev;
 175     size_t maxruns = 1000;
 176
 177     for( size_t i = 0; i < G.nrNodes(); i++ )
 178         bforeach( const Neighbor &j, G.nb(i) )
 179             M[i][j.iter] = 0.1;
 180
 181     size_t run=0;
 182     do {
 183         maxdev=0.0;
 184         run++;
 185         for( size_t i = 0; i < G.nrNodes(); i++ ) {
 186             bforeach( const Neighbor &j, G.nb(i) ) {
 187                 size_t _j = j.iter;
 188                 size_t _i = G.findNb(j,i);
 189                 DAI_ASSERT( G.nb(j,_i) == i );
 190
 191                 Real newM = 0.0;
 192                 if( props.updates == Properties::UpdateType::FULL ) {
 193                     // find indices in nb(j) that do not correspond with i
 194                     sub_nb _nbj_min_i(G.nb(j).size());
 195                     _nbj_min_i.set();
 196                     _nbj_min_i.reset(_i);
 197
 198                     // find indices in nb(i) that do not correspond with j
 199                     sub_nb _nbi_min_j(G.nb(i).size());
 200                     _nbi_min_j.set();
 201                     _nbi_min_j.reset(_j);
 202
 203                     sum_subs(j, _nbj_min_i, &sum_even, &sum_odd);
 204                     newM = (tanh(theta[j]) * sum_even + sum_odd) / (sum_even + tanh(theta[j]) * sum_odd);
 205
 206                     sum_subs(i, _nbi_min_j, &sum_even, &sum_odd);
 207                     Real denom = sum_even + tanh(theta[i]) * sum_odd;
 208                     Real numer = 0.0;
 209                     for(size_t _k=0; _k < G.nb(i).size(); _k++) if(_k != _j) {
 210                         sub_nb _nbi_min_jk(_nbi_min_j);
 211                         _nbi_min_jk.reset(_k);
 212                         sum_subs(i, _nbi_min_jk, &sum_even, &sum_odd);
 213                         numer += tJ[i][_k] * cors[i][_j][_k] * (tanh(theta[i]) * sum_even + sum_odd);
 214                     }
 215                     newM -= numer / denom;
 216                 } else if( props.updates == Properties::UpdateType::LINEAR ) {
 217                     newM = T(j,_i);
 218                     for(size_t _l=0; _l<G.nb(i).size(); _l++) if( _l != _j )
 219                         newM -= Omega(i,_j,_l) * tJ[i][_l] * cors[i][_j][_l];
 220                     for(size_t _l1=0; _l1<G.nb(j).size(); _l1++) if( _l1 != _i )
 221                         for( size_t _l2=_l1+1; _l2<G.nb(j).size(); _l2++) if( _l2 != _i)
 222                             newM += Gamma(j,_i,_l1,_l2) * tJ[j][_l1] * tJ[j][_l2] * cors[j][_l1][_l2];
 223                 }
 224
 225                 Real dev = newM - M[i][_j];
 226 //              dev *= 0.02;
 227                 if( abs(dev) >= maxdev )
 228                     maxdev = abs(dev);
 229
 230                 newM = M[i][_j] + dev;
 231                 if( abs(newM) > 1.0 )
 232                     newM = (newM > 0.0) ? 1.0 : -1.0;
 233                 M[i][_j] = newM;
 234             }
 235         }
 236     } while((maxdev>props.tol)&&(run<maxruns));
 237
 238     _iters = run;
 239     if( maxdev > _maxdiff )
 240         _maxdiff = maxdev;
 241
 242     if(run==maxruns){
 243         if( props.verbose >= 1 )
 244             cerr << "MR::propagateCavityFields: Convergence not reached (maxdev=" << maxdev << ")..." << endl;
 245     }
 246 }
 247
 248
 249 void MR::calcMagnetizations() {
 250     for( size_t i = 0; i < G.nrNodes(); i++ ) {
 251         if( props.updates == Properties::UpdateType::FULL ) {
 252             // find indices in nb(i)
 253             sub_nb _nbi( G.nb(i).size() );
 254             _nbi.set();
 255
 256             // calc numerator1 and denominator1
 257             Real sum_even, sum_odd;
 258             sum_subs(i, _nbi, &sum_even, &sum_odd);
 259
 260             Mag[i] = (tanh(theta[i]) * sum_even + sum_odd) / (sum_even + tanh(theta[i]) * sum_odd);
 261
 262         } else if( props.updates == Properties::UpdateType::LINEAR ) {
 263             sub_nb empty( G.nb(i).size() );
 264             Mag[i] = T(i,empty);
 265
 266             for( size_t _l1 = 0; _l1 < G.nb(i).size(); _l1++ )
 267                 for( size_t _l2 = _l1 + 1; _l2 < G.nb(i).size(); _l2++ )
 268                     Mag[i] += Gamma(i,_l1,_l2) * tJ[i][_l1] * tJ[i][_l2] * cors[i][_l1][_l2];
 269         }
 270         if( abs( Mag[i] ) > 1.0 )
 271             Mag[i] = (Mag[i] > 0.0) ? 1.0 : -1.0;
 272     }
 273 }
 274
 275
 276 Real MR::calcCavityCorrelations() {
 277     Real md = 0.0;
 278     for( size_t i = 0; i < nrVars(); i++ ) {
 279         vector<Factor> pairq;
 280         if( props.inits == Properties::InitType::EXACT ) {
 281             JTree jtcav(*this, PropertySet()("updates", string("HUGIN"))("verbose", (size_t)0) );
 282             jtcav.makeCavity( i );
 283             pairq = calcPairBeliefs( jtcav, delta(i), false, true );
 284         } else if( props.inits == Properties::InitType::CLAMPING ) {
 285             BP bpcav(*this, PropertySet()("updates", string("SEQMAX"))("tol", (Real)1.0e-9)("maxiter", (size_t)10000)("verbose", (size_t)0)("logdomain", false));
 286             bpcav.makeCavity( i );
 287
 288             pairq = calcPairBeliefs( bpcav, delta(i), false, true );
 289             md = std::max( md, bpcav.maxDiff() );
 290         } else if( props.inits == Properties::InitType::RESPPROP ) {
 291             BP bpcav(*this, PropertySet()("updates", string("SEQMAX"))("tol", (Real)1.0e-9)("maxiter", (size_t)10000)("verbose", (size_t)0)("logdomain", false));
 292             bpcav.makeCavity( i );
 293             bpcav.makeCavity( i );
 294             bpcav.init();
 295             bpcav.run();
 296
 297             BBP bbp( &bpcav, PropertySet()("verbose",(size_t)0)("tol",(Real)1.0e-9)("maxiter",(size_t)10000)("damping",(Real)0.0)("updates",string("SEQ_MAX")) );
 298             bforeach( const Neighbor &j, G.nb(i) ) {
 299                 // Create weights for magnetization of some spin
 300                 Prob p( 2, 0.0 );
 301                 p.set( 0, -1.0 );
 302                 p.set( 1, 1.0 );
 303
 304                 // BBP cost function would be the magnetization of spin j
 305                 vector<Prob> b1_adj;
 306                 b1_adj.reserve( nrVars() );
 307                 for( size_t l = 0; l < nrVars(); l++ )
 308                     if( l == j )
 309                         b1_adj.push_back( p );
 310                     else
 311                         b1_adj.push_back( Prob( 2, 0.0 ) );
 312                 bbp.init_V( b1_adj );
 313
 314                 // run BBP to estimate adjoints
 315                 bbp.run();
 316
 317                 bforeach( const Neighbor &k, G.nb(i) ) {
 318                     if( k != j )
 319                         cors[i][j.iter][k.iter] = (bbp.adj_psi_V(k)[1] - bbp.adj_psi_V(k)[0]);
 320                     else
 321                         cors[i][j.iter][k.iter] = 0.0;
 322                 }
 323             }
 324         }
 325
 326         if( props.inits != Properties::InitType::RESPPROP ) {
 327             for( size_t jk = 0; jk < pairq.size(); jk++ ) {
 328                 VarSet::const_iterator kit = pairq[jk].vars().begin();
 329                 size_t j = findVar( *(kit) );
 330                 size_t k = findVar( *(++kit) );
 331                 pairq[jk].normalize();
 332                 Real cor = (pairq[jk][3] - pairq[jk][2] - pairq[jk][1] + pairq[jk][0]) - (pairq[jk][3] + pairq[jk][2] - pairq[jk][1] - pairq[jk][0]) * (pairq[jk][3] - pairq[jk][2] + pairq[jk][1] - pairq[jk][0]);
 333
 334                 size_t _j = G.findNb(i,j);
 335                 size_t _k = G.findNb(i,k);
 336                 cors[i][_j][_k] = cor;
 337                 cors[i][_k][_j] = cor;
 338             }
 339         }
 340
 341     }
 342     return md;
 343 }
 344
 345
 346 Real MR::run() {
 347     if( supported ) {
 348         if( props.verbose >= 1 )
 349             cerr << "Starting " << identify() << "...";
 350
 351         double tic = toc();
 352
 353         // approximate correlations of cavity spins
 354         Real md = calcCavityCorrelations();
 355         if( md > _maxdiff )
 356             _maxdiff = md;
 357
 358         // solve messages
 359         propagateCavityFields();
 360
 361         // calculate magnetizations
 362         calcMagnetizations();
 363
 364         if( props.verbose >= 1 )
 365             cerr << name() << " needed " << toc() - tic << " seconds." << endl;
 366
 367         return _maxdiff;
 368     } else
 369         return 1.0;
 370 }
 371
 372
 373 Factor MR::beliefV( size_t i ) const {
 374     if( supported ) {
 375         Real x[2];
 376         x[0] = 0.5 - Mag[i] / 2.0;
 377         x[1] = 0.5 + Mag[i] / 2.0;
 378
 379         return Factor( var(i), x );
 380     } else
 381         return Factor();
 382 }
 383
 384
 385 Factor MR::belief (const VarSet &ns) const {
 386     if( ns.size() == 0 )
 387         return Factor();
 388     else if( ns.size() == 1 )
 389         return beliefV( findVar( *(ns.begin()) ) );
 390     else {
 391         DAI_THROW(BELIEF_NOT_AVAILABLE);
 392         return Factor();
 393     }
 394 }
 395
 396
 397 vector<Factor> MR::beliefs() const {
 398     vector<Factor> result;
 399     for( size_t i = 0; i < nrVars(); i++ )
 400         result.push_back( beliefV( i ) );
 401     return result;
 402 }
 403
 404
 405 MR::MR( const FactorGraph &fg, const PropertySet &opts ) : DAIAlgFG(fg), supported(true), _maxdiff(0.0), _iters(0) {
 406     setProperties( opts );
 407
 408     size_t N = fg.nrVars();
 409
 410     // check whether all vars in fg are binary
 411     for( size_t i = 0; i < N; i++ )
 412         if( (fg.var(i).states() > 2) ) {
 413             supported = false;
 414             break;
 415         }
 416     if( !supported )
 417         DAI_THROWE(NOT_IMPLEMENTED,"MR only supports binary variables");
 418
 419     // check whether all interactions are pairwise or single
 420     // and construct Markov graph
 421     G = GraphAL(N);
 422     for( size_t I = 0; I < fg.nrFactors(); I++ ) {
 423         const Factor &psi = fg.factor(I);
 424         if( psi.vars().size() > 2 ) {
 425             supported = false;
 426             break;
 427         } else if( psi.vars().size() == 2 ) {
 428             VarSet::const_iterator jit = psi.vars().begin();
 429             size_t i = fg.findVar( *(jit) );
 430             size_t j = fg.findVar( *(++jit) );
 431             G.addEdge( i, j, false );
 432         }
 433     }
 434     if( !supported )
 435         DAI_THROWE(NOT_IMPLEMENTED,"MR does not support higher order interactions (only single and pairwise are supported)");
 436
 437     // construct theta
 438     theta.clear();
 439     theta.resize( N, 0.0 );
 440
 441     // construct tJ
 442     tJ.resize( N );
 443     for( size_t i = 0; i < N; i++ )
 444         tJ[i].resize( G.nb(i).size(), 0.0 );
 445
 446     // initialize theta and tJ
 447     for( size_t I = 0; I < fg.nrFactors(); I++ ) {
 448         const Factor &psi = fg.factor(I);
 449         if( psi.vars().size() == 1 ) {
 450             size_t i = fg.findVar( *(psi.vars().begin()) );
 451             theta[i] += 0.5 * log(psi[1] / psi[0]);
 452         } else if( psi.vars().size() == 2 ) {
 453             VarSet::const_iterator jit = psi.vars().begin();
 454             size_t i = fg.findVar( *(jit) );
 455             size_t j = fg.findVar( *(++jit) );
 456
 457             Real w_ij = 0.25 * log(psi[3] * psi[0] / (psi[2] * psi[1]));
 458             tJ[i][G.findNb(i,j)] += w_ij;
 459             tJ[j][G.findNb(j,i)] += w_ij;
 460
 461             theta[i] += 0.25 * log(psi[3] / psi[2] * psi[1] / psi[0]);
 462             theta[j] += 0.25 * log(psi[3] / psi[1] * psi[2] / psi[0]);
 463         }
 464     }
 465     for( size_t i = 0; i < N; i++ )
 466         bforeach( const Neighbor &j, G.nb(i) )
 467             tJ[i][j.iter] = tanh( tJ[i][j.iter] );
 468
 469     // construct M
 470     M.resize( N );
 471     for( size_t i = 0; i < N; i++ )
 472         M[i].resize( G.nb(i).size() );
 473
 474     // construct cors
 475     cors.resize( N );
 476     for( size_t i = 0; i < N; i++ )
 477         cors[i].resize( G.nb(i).size() );
 478     for( size_t i = 0; i < N; i++ )
 479         for( size_t _j = 0; _j < cors[i].size(); _j++ )
 480             cors[i][_j].resize( G.nb(i).size() );
 481
 482     // construct Mag
 483     Mag.resize( N );
 484 }
 485
 486
 487 } // end of namespace dai
 488
 489
 490 #endif