src/bp_dual.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_CBP
  11
  12
  13 #include <iostream>
  14 #include <sstream>
  15 #include <algorithm>
  16
  17 #include <dai/bp_dual.h>
  18 #include <dai/util.h>
  19 #include <dai/bipgraph.h>
  20
  21
  22 namespace dai {
  23
  24
  25 using namespace std;
  26
  27
  28 void BP_dual::init() {
  29     regenerateMessages();
  30     regenerateBeliefs();
  31     calcMessages();
  32     calcBeliefs();
  33 }
  34
  35
  36 void BP_dual::regenerateMessages() {
  37     size_t nv = fg().nrVars();
  38     _msgs.Zn.resize(nv);
  39     _msgs.Zm.resize(nv);
  40     _msgs.m.resize(nv);
  41     _msgs.n.resize(nv);
  42     for( size_t i = 0; i < nv; i++ ) {
  43         size_t nvf = fg().nbV(i).size();
  44         _msgs.Zn[i].resize(nvf, 1.0);
  45         _msgs.Zm[i].resize(nvf, 1.0);
  46         size_t states = fg().var(i).states();
  47         _msgs.n[i].resize(nvf, Prob(states));
  48         _msgs.m[i].resize(nvf, Prob(states));
  49     }
  50 }
  51
  52
  53 void BP_dual::regenerateBeliefs() {
  54     _beliefs.b1.clear();
  55     _beliefs.b1.reserve(fg().nrVars());
  56     _beliefs.Zb1.resize(fg().nrVars(), 1.0);
  57     _beliefs.b2.clear();
  58     _beliefs.b2.reserve(fg().nrFactors());
  59     _beliefs.Zb2.resize(fg().nrFactors(), 1.0);
  60
  61     for( size_t i = 0; i < fg().nrVars(); i++ )
  62         _beliefs.b1.push_back( Prob( fg().var(i).states() ) );
  63     for( size_t I = 0; I < fg().nrFactors(); I++ )
  64         _beliefs.b2.push_back( Prob( fg().factor(I).nrStates() ) );
  65 }
  66
  67
  68 void BP_dual::calcMessages() {
  69     // calculate 'n' messages from "factor marginal / factor"
  70     for( size_t I = 0; I < fg().nrFactors(); I++ ) {
  71         Factor f = _ia->beliefF(I) / fg().factor(I);
  72         bforeach( const Neighbor &i, fg().nbF(I) )
  73             msgN(i, i.dual) = f.marginal( fg().var(i) ).p();
  74     }
  75     // calculate 'm' messages and normalizers from 'n' messages
  76     for( size_t i = 0; i < fg().nrVars(); i++ )
  77         bforeach( const Neighbor &I, fg().nbV(i) )
  78             calcNewM( i, I.iter );
  79     // recalculate 'n' messages and normalizers from 'm' messages
  80     for( size_t i = 0; i < fg().nrVars(); i++ )
  81         bforeach( const Neighbor &I, fg().nbV(i) )
  82             calcNewN(i, I.iter);
  83 }
  84
  85
  86 void BP_dual::calcNewM( size_t i, size_t _I ) {
  87     // calculate updated message I->i
  88     const Neighbor &I = fg().nbV(i)[_I];
  89     Prob prod( fg().factor(I).p() );
  90     bforeach( const Neighbor &j, fg().nbF(I) )
  91         if( j != i ) { // for all j in I \ i
  92             Prob &n = msgN(j,j.dual);
  93             IndexFor ind( fg().var(j), fg().factor(I).vars() );
  94             for( size_t x = 0; ind.valid(); x++, ++ind )
  95                 prod.set( x, prod[x] * n[ind] );
  96         }
  97     // Marginalize onto i
  98     Prob marg( fg().var(i).states(), 0.0 );
  99     // ind is the precalculated Index(i,I) i.e. to x_I == k corresponds x_i == ind[k]
 100     IndexFor ind( fg().var(i), fg().factor(I).vars() );
 101     for( size_t x = 0; ind.valid(); x++, ++ind )
 102         marg.set( ind, marg[ind] + prod[x] );
 103
 104     _msgs.Zm[i][_I] = marg.normalize();
 105     _msgs.m[i][_I] = marg;
 106 }
 107
 108
 109 void BP_dual::calcNewN( size_t i, size_t _I ) {
 110     // calculate updated message i->I
 111     const Neighbor &I = fg().nbV(i)[_I];
 112     Prob prod( fg().var(i).states(), 1.0 );
 113     bforeach( const Neighbor &J, fg().nbV(i) )
 114         if( J.node != I.node ) // for all J in i \ I
 115             prod *= msgM(i,J.iter);
 116     _msgs.Zn[i][_I] = prod.normalize();
 117     _msgs.n[i][_I] = prod;
 118 }
 119
 120
 121 void BP_dual::calcBeliefs() {
 122     for( size_t i = 0; i < fg().nrVars(); i++ )
 123         calcBeliefV(i);  // calculate b_i
 124     for( size_t I = 0; I < fg().nrFactors(); I++ )
 125         calcBeliefF(I);  // calculate b_I
 126 }
 127
 128
 129 void BP_dual::calcBeliefV( size_t i ) {
 130     Prob prod( fg().var(i).states(), 1.0 );
 131     bforeach( const Neighbor &I, fg().nbV(i) )
 132         prod *= msgM(i,I.iter);
 133     _beliefs.Zb1[i] = prod.normalize();
 134     _beliefs.b1[i] = prod;
 135 }
 136
 137
 138 void BP_dual::calcBeliefF( size_t I ) {
 139     Prob prod( fg().factor(I).p() );
 140     bforeach( const Neighbor &j, fg().nbF(I) ) {
 141         IndexFor ind( fg().var(j), fg().factor(I).vars() );
 142         Prob n( msgN(j,j.dual) );
 143         for( size_t x = 0; ind.valid(); x++, ++ind )
 144             prod.set( x, prod[x] * n[ind] );
 145     }
 146     _beliefs.Zb2[I] = prod.normalize();
 147     _beliefs.b2[I] = prod;
 148 }
 149
 150
 151 } // end of namespace dai
 152
 153
 154 #endif