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