src/regiongraph.cpp

   1 /*  This file is part of libDAI - http://www.libdai.org/
   2  *
   3  *  libDAI is licensed under the terms of the GNU General Public License version
   4  *  2, or (at your option) any later version. libDAI is distributed without any
   5  *  warranty. See the file COPYING for more details.
   6  *
   7  *  Copyright (C) 2006-2009  Joris Mooij  [joris dot mooij at libdai dot org]
   8  *  Copyright (C) 2006-2007  Radboud University Nijmegen, The Netherlands
   9  */
  10
  11
  12 #include <algorithm>
  13 #include <cmath>
  14 #include <boost/dynamic_bitset.hpp>
  15 #include <dai/regiongraph.h>
  16 #include <dai/factorgraph.h>
  17 #include <dai/clustergraph.h>
  18
  19
  20 namespace dai {
  21
  22
  23 using namespace std;
  24
  25
  26 RegionGraph::RegionGraph( const FactorGraph &fg, const std::vector<Region> &ors, const std::vector<Region> &irs, const std::vector<std::pair<size_t,size_t> > &edges) : FactorGraph(fg), G(), ORs(), IRs(irs), fac2OR() {
  27     // Copy outer regions (giving them counting number 1.0)
  28     ORs.reserve( ors.size() );
  29     for( vector<Region>::const_iterator alpha = ors.begin(); alpha != ors.end(); alpha++ )
  30         ORs.push_back( FRegion(Factor(*alpha, 1.0), 1.0) );
  31
  32     // For each factor, find an outer region that subsumes that factor.
  33     // Then, multiply the outer region with that factor.
  34     fac2OR.reserve( nrFactors() );
  35     for( size_t I = 0; I < nrFactors(); I++ ) {
  36         size_t alpha;
  37         for( alpha = 0; alpha < nrORs(); alpha++ )
  38             if( OR(alpha).vars() >> factor(I).vars() ) {
  39                 fac2OR.push_back( alpha );
  40                 break;
  41             }
  42         DAI_ASSERT( alpha != nrORs() );
  43     }
  44     RecomputeORs();
  45
  46     // Create bipartite graph
  47     G.construct( nrORs(), nrIRs(), edges.begin(), edges.end() );
  48
  49     // Check counting numbers
  50 #ifdef DAI_DEBUG
  51     checkCountingNumbers();
  52 #endif
  53 }
  54
  55
  56 // CVM style
  57 RegionGraph::RegionGraph( const FactorGraph &fg, const std::vector<VarSet> &cl ) : FactorGraph(fg), G(), ORs(), IRs(), fac2OR() {
  58     // Retain only maximal clusters
  59     ClusterGraph cg( cl );
  60     cg.eraseNonMaximal();
  61
  62     // Create outer regions, giving them counting number 1.0
  63     ORs.reserve( cg.size() );
  64     foreach( const VarSet &ns, cg.clusters )
  65         ORs.push_back( FRegion(Factor(ns, 1.0), 1.0) );
  66
  67     // For each factor, find an outer regions that subsumes that factor.
  68     // Then, multiply the outer region with that factor.
  69     fac2OR.reserve( nrFactors() );
  70     for( size_t I = 0; I < nrFactors(); I++ ) {
  71         size_t alpha;
  72         for( alpha = 0; alpha < nrORs(); alpha++ )
  73             if( OR(alpha).vars() >> factor(I).vars() ) {
  74                 fac2OR.push_back( alpha );
  75                 break;
  76             }
  77         DAI_ASSERT( alpha != nrORs() );
  78     }
  79     RecomputeORs();
  80
  81     // Create inner regions - first pass
  82     set<VarSet> betas;
  83     for( size_t alpha = 0; alpha < cg.clusters.size(); alpha++ )
  84         for( size_t alpha2 = alpha; (++alpha2) != cg.clusters.size(); ) {
  85             VarSet intersection = cg.clusters[alpha] & cg.clusters[alpha2];
  86             if( intersection.size() > 0 )
  87                 betas.insert( intersection );
  88         }
  89
  90     // Create inner regions - subsequent passes
  91     set<VarSet> new_betas;
  92     do {
  93         new_betas.clear();
  94         for( set<VarSet>::const_iterator gamma = betas.begin(); gamma != betas.end(); gamma++ )
  95             for( set<VarSet>::const_iterator gamma2 = gamma; (++gamma2) != betas.end(); ) {
  96                 VarSet intersection = (*gamma) & (*gamma2);
  97                 if( (intersection.size() > 0) && (betas.count(intersection) == 0) )
  98                     new_betas.insert( intersection );
  99             }
 100         betas.insert(new_betas.begin(), new_betas.end());
 101     } while( new_betas.size() );
 102
 103     // Create inner regions - store them in the bipartite graph
 104     IRs.reserve( betas.size() );
 105     for( set<VarSet>::const_iterator beta = betas.begin(); beta != betas.end(); beta++ )
 106         IRs.push_back( Region(*beta,0.0) );
 107
 108     // Create edges
 109     vector<pair<size_t,size_t> > edges;
 110     for( size_t beta = 0; beta < nrIRs(); beta++ ) {
 111         for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
 112             if( OR(alpha).vars() >> IR(beta) )
 113                 edges.push_back( pair<size_t,size_t>(alpha,beta) );
 114         }
 115     }
 116
 117     // Create bipartite graph
 118     G.construct( nrORs(), nrIRs(), edges.begin(), edges.end() );
 119
 120     // Calculate counting numbers
 121     calcCountingNumbers();
 122
 123     // Check counting numbers
 124 #ifdef DAI_DEBUG
 125     checkCountingNumbers();
 126 #endif
 127 }
 128
 129
 130 void RegionGraph::calcCountingNumbers() {
 131     // Calculates counting numbers of inner regions based upon counting numbers of outer regions
 132
 133     vector<vector<size_t> > ancestors(nrIRs());
 134     boost::dynamic_bitset<> assigned(nrIRs());
 135     for( size_t beta = 0; beta < nrIRs(); beta++ ) {
 136         IR(beta).c() = 0.0;
 137         for( size_t beta2 = 0; beta2 < nrIRs(); beta2++ )
 138             if( (beta2 != beta) && IR(beta2) >> IR(beta) )
 139                 ancestors[beta].push_back(beta2);
 140     }
 141
 142     bool new_counting;
 143     do {
 144         new_counting = false;
 145         for( size_t beta = 0; beta < nrIRs(); beta++ ) {
 146             if( !assigned[beta] ) {
 147                 bool has_unassigned_ancestor = false;
 148                 for( vector<size_t>::const_iterator beta2 = ancestors[beta].begin(); (beta2 != ancestors[beta].end()) && !has_unassigned_ancestor; beta2++ )
 149                     if( !assigned[*beta2] )
 150                         has_unassigned_ancestor = true;
 151                 if( !has_unassigned_ancestor ) {
 152                     Real c = 1.0;
 153                     foreach( const Neighbor &alpha, nbIR(beta) )
 154                         c -= OR(alpha).c();
 155                     for( vector<size_t>::const_iterator beta2 = ancestors[beta].begin(); beta2 != ancestors[beta].end(); beta2++ )
 156                         c -= IR(*beta2).c();
 157                     IR(beta).c() = c;
 158                     assigned.set(beta, true);
 159                     new_counting = true;
 160                 }
 161             }
 162         }
 163     } while( new_counting );
 164 }
 165
 166
 167 bool RegionGraph::checkCountingNumbers() const {
 168     // Checks whether the counting numbers satisfy the fundamental relation
 169
 170     bool all_valid = true;
 171     for( vector<Var>::const_iterator n = vars().begin(); n != vars().end(); n++ ) {
 172         Real c_n = 0.0;
 173         for( size_t alpha = 0; alpha < nrORs(); alpha++ )
 174             if( OR(alpha).vars().contains( *n ) )
 175                 c_n += OR(alpha).c();
 176         for( size_t beta = 0; beta < nrIRs(); beta++ )
 177             if( IR(beta).contains( *n ) )
 178                 c_n += IR(beta).c();
 179         if( fabs(c_n - 1.0) > 1e-15 ) {
 180             all_valid = false;
 181             cerr << "WARNING: counting numbers do not satisfy relation for " << *n << "(c_n = " << c_n << ")." << endl;
 182         }
 183     }
 184
 185     return all_valid;
 186 }
 187
 188
 189 void RegionGraph::RecomputeORs() {
 190     for( size_t alpha = 0; alpha < nrORs(); alpha++ )
 191         OR(alpha).fill( 1.0 );
 192     for( size_t I = 0; I < nrFactors(); I++ )
 193         if( fac2OR[I] != -1U )
 194             OR( fac2OR[I] ) *= factor( I );
 195 }
 196
 197
 198 void RegionGraph::RecomputeORs( const VarSet &ns ) {
 199     for( size_t alpha = 0; alpha < nrORs(); alpha++ )
 200         if( OR(alpha).vars().intersects( ns ) )
 201             OR(alpha).fill( 1.0 );
 202     for( size_t I = 0; I < nrFactors(); I++ )
 203         if( fac2OR[I] != -1U )
 204             if( OR( fac2OR[I] ).vars().intersects( ns ) )
 205                 OR( fac2OR[I] ) *= factor( I );
 206 }
 207
 208
 209 void RegionGraph::RecomputeOR( size_t I ) {
 210     DAI_ASSERT( I < nrFactors() );
 211     if( fac2OR[I] != -1U ) {
 212         size_t alpha = fac2OR[I];
 213         OR(alpha).fill( 1.0 );
 214         for( size_t J = 0; J < nrFactors(); J++ )
 215             if( fac2OR[J] == alpha )
 216                 OR(alpha) *= factor( J );
 217     }
 218 }
 219
 220
 221 /// Send RegionGraph to output stream
 222 ostream & operator << (ostream & os, const RegionGraph & rg) {
 223     os << "Outer regions" << endl;
 224     for( size_t alpha = 0; alpha < rg.nrORs(); alpha++ )
 225         os << alpha << ": " << rg.OR(alpha).vars() << ": c = " << rg.OR(alpha).c() << endl;
 226
 227     os << "Inner regions" << endl;
 228     for( size_t beta = 0; beta < rg.nrIRs(); beta++ )
 229         os << beta << ": " << (VarSet)rg.IR(beta) << ": c = " << rg.IR(beta).c() << endl;
 230
 231     os << "Edges" << endl;
 232     for( size_t alpha = 0; alpha < rg.nrORs(); alpha++ )
 233         foreach( const RegionGraph::Neighbor &beta, rg.nbOR(alpha) )
 234             os << alpha << "->" << beta << endl;
 235
 236     return(os);
 237 }
 238
 239
 240 } // end of namespace dai