src/regiongraph.cpp

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