src/daialg.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 <vector>
  23 #include <dai/daialg.h>
  24
  25
  26 namespace dai {
  27
  28
  29 using namespace std;
  30
  31
  32 /// Calculate the marginal of obj on ns by clamping
  33 /// all variables in ns and calculating logZ for each joined state
  34 Factor calcMarginal( const InfAlg & obj, const VarSet & ns, bool reInit ) {
  35     Factor Pns (ns);
  36
  37     InfAlg *clamped = obj.clone();
  38     if( !reInit )
  39         clamped->init();
  40
  41     Real logZ0 = 0.0;
  42     for( State s(ns); s.valid(); s++ ) {
  43         // save unclamped factors connected to ns
  44         clamped->saveProbs( ns );
  45
  46         // set clamping Factors to delta functions
  47         for( VarSet::const_iterator n = ns.begin(); n != ns.end(); n++ )
  48             clamped->clamp( *n, s(*n) );
  49
  50         // run DAIAlg, calc logZ, store in Pns
  51         if( reInit )
  52             clamped->init();
  53         clamped->run();
  54
  55         Real Z;
  56         if( s == 0 ) {
  57             logZ0 = clamped->logZ();
  58             Z = 1.0;
  59         } else {
  60             // subtract logZ0 to avoid very large numbers
  61             Z = exp(clamped->logZ() - logZ0);
  62         }
  63
  64         Pns[s] = Z;
  65
  66         // restore clamped factors
  67         clamped->undoProbs( ns );
  68     }
  69
  70     delete clamped;
  71
  72     return( Pns.normalized(Prob::NORMPROB) );
  73 }
  74
  75
  76 vector<Factor> calcPairBeliefs( const InfAlg & obj, const VarSet& ns, bool reInit ) {
  77     // convert ns to vector<VarSet>
  78     size_t N = ns.size();
  79     vector<Var> vns;
  80     vns.reserve( N );
  81     for( VarSet::const_iterator n = ns.begin(); n != ns.end(); n++ )
  82         vns.push_back( *n );
  83
  84     vector<Factor> pairbeliefs;
  85     pairbeliefs.reserve( N * N );
  86     for( size_t j = 0; j < N; j++ )
  87         for( size_t k = 0; k < N; k++ )
  88             if( j == k )
  89                 pairbeliefs.push_back(Factor());
  90             else
  91                 pairbeliefs.push_back(Factor(vns[j] | vns[k]));
  92
  93     InfAlg *clamped = obj.clone();
  94     if( !reInit )
  95         clamped->init();
  96
  97     Real logZ0 = 0.0;
  98     for( size_t j = 0; j < N; j++ ) {
  99         // clamp Var j to its possible values
 100         for( size_t j_val = 0; j_val < vns[j].states(); j_val++ ) {
 101             // save unclamped factors connected to ns
 102             clamped->saveProbs( ns );
 103
 104             clamped->clamp( vns[j], j_val );
 105             if( reInit )
 106                 clamped->init();
 107             clamped->run();
 108
 109             //if( j == 0 )
 110             //  logZ0 = obj.logZ();
 111             double Z_xj = 1.0;
 112             if( j == 0 && j_val == 0 ) {
 113                 logZ0 = clamped->logZ();
 114             } else {
 115                 // subtract logZ0 to avoid very large numbers
 116                 Z_xj = exp(clamped->logZ() - logZ0);
 117             }
 118
 119             for( size_t k = 0; k < N; k++ )
 120                 if( k != j ) {
 121                     Factor b_k = clamped->belief(vns[k]);
 122                     for( size_t k_val = 0; k_val < vns[k].states(); k_val++ )
 123                         if( vns[j].label() < vns[k].label() )
 124                             pairbeliefs[j * N + k][j_val + (k_val * vns[j].states())] = Z_xj * b_k[k_val];
 125                         else
 126                             pairbeliefs[j * N + k][k_val + (j_val * vns[k].states())] = Z_xj * b_k[k_val];
 127                 }
 128
 129             // restore clamped factors
 130             clamped->undoProbs( ns );
 131         }
 132     }
 133
 134     delete clamped;
 135
 136     // Calculate result by taking the geometric average
 137     vector<Factor> result;
 138     result.reserve( N * (N - 1) / 2 );
 139     for( size_t j = 0; j < N; j++ )
 140         for( size_t k = j+1; k < N; k++ )
 141             result.push_back( (pairbeliefs[j * N + k] * pairbeliefs[k * N + j]) ^ 0.5 );
 142
 143     return result;
 144 }
 145
 146
 147 Factor calcMarginal2ndO( const InfAlg & obj, const VarSet& ns, bool reInit ) {
 148     // returns a a probability distribution whose 1st order interactions
 149     // are unspecified, whose 2nd order interactions approximate those of
 150     // the marginal on ns, and whose higher order interactions are absent.
 151
 152     vector<Factor> pairbeliefs = calcPairBeliefs( obj, ns, reInit );
 153
 154     Factor Pns (ns);
 155     for( size_t ij = 0; ij < pairbeliefs.size(); ij++ )
 156         Pns *= pairbeliefs[ij];
 157
 158     return( Pns.normalized(Prob::NORMPROB) );
 159 }
 160
 161
 162 vector<Factor> calcPairBeliefsNew( const InfAlg & obj, const VarSet& ns, bool reInit ) {
 163     vector<Factor> result;
 164     result.reserve( ns.size() * (ns.size() - 1) / 2 );
 165
 166     InfAlg *clamped = obj.clone();
 167     if( !reInit )
 168         clamped->init();
 169
 170     Real logZ0 = 0.0;
 171     VarSet::const_iterator nj = ns.begin();
 172     for( long j = 0; j < (long)ns.size() - 1; j++, nj++ ) {
 173         size_t k = 0;
 174         for( VarSet::const_iterator nk = nj; (++nk) != ns.end(); k++ ) {
 175             Factor pairbelief( *nj | *nk );
 176
 177             // clamp Vars j and k to their possible values
 178             for( size_t j_val = 0; j_val < nj->states(); j_val++ )
 179                 for( size_t k_val = 0; k_val < nk->states(); k_val++ ) {
 180                     // save unclamped factors connected to ns
 181                     clamped->saveProbs( ns );
 182
 183                     clamped->clamp( *nj, j_val );
 184                     clamped->clamp( *nk, k_val );
 185                     if( reInit )
 186                         clamped->init();
 187                     clamped->run();
 188
 189                     double Z_xj = 1.0;
 190                     if( j_val == 0 && k_val == 0 ) {
 191                         logZ0 = clamped->logZ();
 192                     } else {
 193                         // subtract logZ0 to avoid very large numbers
 194                         Z_xj = exp(clamped->logZ() - logZ0);
 195                     }
 196
 197                     // we assume that j.label() < k.label()
 198                     // i.e. we make an assumption here about the indexing
 199                     pairbelief[j_val + (k_val * nj->states())] = Z_xj;
 200
 201                     // restore clamped factors
 202                     clamped->undoProbs( ns );
 203                 }
 204
 205             result.push_back( pairbelief );
 206         }
 207     }
 208
 209     delete clamped;
 210
 211     assert( result.size() == (ns.size() * (ns.size() - 1) / 2) );
 212
 213     return result;
 214 }
 215
 216
 217 } // end of namespace dai