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