src/factorgraph.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 <iostream>
  13 #include <iomanip>
  14 #include <iterator>
  15 #include <map>
  16 #include <set>
  17 #include <fstream>
  18 #include <string>
  19 #include <algorithm>
  20 #include <functional>
  21 #include <dai/factorgraph.h>
  22 #include <dai/util.h>
  23 #include <dai/exceptions.h>
  24 #include <boost/lexical_cast.hpp>
  25
  26
  27 namespace dai {
  28
  29
  30 using namespace std;
  31
  32
  33 FactorGraph::FactorGraph( const std::vector<Factor> &P ) : G(), _backup() {
  34     // add factors, obtain variables
  35     set<Var> varset;
  36     _factors.reserve( P.size() );
  37     size_t nrEdges = 0;
  38     for( vector<Factor>::const_iterator p2 = P.begin(); p2 != P.end(); p2++ ) {
  39         _factors.push_back( *p2 );
  40         copy( p2->vars().begin(), p2->vars().end(), inserter( varset, varset.begin() ) );
  41         nrEdges += p2->vars().size();
  42     }
  43
  44     // add vars
  45     _vars.reserve( varset.size() );
  46     for( set<Var>::const_iterator p1 = varset.begin(); p1 != varset.end(); p1++ )
  47         _vars.push_back( *p1 );
  48
  49     // create graph structure
  50     constructGraph( nrEdges );
  51 }
  52
  53
  54 void FactorGraph::constructGraph( size_t nrEdges ) {
  55     // create a mapping for indices
  56     hash_map<size_t, size_t> hashmap;
  57
  58     for( size_t i = 0; i < vars().size(); i++ )
  59         hashmap[var(i).label()] = i;
  60
  61     // create edge list
  62     vector<Edge> edges;
  63     edges.reserve( nrEdges );
  64     for( size_t i2 = 0; i2 < nrFactors(); i2++ ) {
  65         const VarSet& ns = factor(i2).vars();
  66         for( VarSet::const_iterator q = ns.begin(); q != ns.end(); q++ )
  67             edges.push_back( Edge(hashmap[q->label()], i2) );
  68     }
  69
  70     // create bipartite graph
  71     G.construct( nrVars(), nrFactors(), edges.begin(), edges.end() );
  72 }
  73
  74
  75 /// Writes a FactorGraph to an output stream
  76 std::ostream& operator<< ( std::ostream &os, const FactorGraph &fg ) {
  77     os << fg.nrFactors() << endl;
  78
  79     for( size_t I = 0; I < fg.nrFactors(); I++ ) {
  80         os << endl;
  81         os << fg.factor(I).vars().size() << endl;
  82         for( VarSet::const_iterator i = fg.factor(I).vars().begin(); i != fg.factor(I).vars().end(); i++ )
  83             os << i->label() << " ";
  84         os << endl;
  85         for( VarSet::const_iterator i = fg.factor(I).vars().begin(); i != fg.factor(I).vars().end(); i++ )
  86             os << i->states() << " ";
  87         os << endl;
  88         size_t nr_nonzeros = 0;
  89         for( size_t k = 0; k < fg.factor(I).states(); k++ )
  90             if( fg.factor(I)[k] != (Real)0 )
  91                 nr_nonzeros++;
  92         os << nr_nonzeros << endl;
  93         for( size_t k = 0; k < fg.factor(I).states(); k++ )
  94             if( fg.factor(I)[k] != (Real)0 )
  95                 os << k << " " << setw(os.precision()+4) << fg.factor(I)[k] << endl;
  96     }
  97
  98     return(os);
  99 }
 100
 101
 102 /// Reads a FactorGraph from an input stream
 103 std::istream& operator>> ( std::istream& is, FactorGraph &fg ) {
 104     long verbose = 0;
 105
 106     vector<Factor> facs;
 107     size_t nr_Factors;
 108     string line;
 109
 110     while( (is.peek()) == '#' )
 111         getline(is,line);
 112     is >> nr_Factors;
 113     if( is.fail() )
 114         DAI_THROWE(INVALID_FACTORGRAPH_FILE,"Cannot read number of factors");
 115     if( verbose >= 2 )
 116         cerr << "Reading " << nr_Factors << " factors..." << endl;
 117
 118     getline (is,line);
 119     if( is.fail() )
 120         DAI_THROWE(INVALID_FACTORGRAPH_FILE,"Expecting empty line");
 121
 122     map<long,size_t> vardims;
 123     for( size_t I = 0; I < nr_Factors; I++ ) {
 124         if( verbose >= 3 )
 125             cerr << "Reading factor " << I << "..." << endl;
 126         size_t nr_members;
 127         while( (is.peek()) == '#' )
 128             getline(is,line);
 129         is >> nr_members;
 130         if( verbose >= 3 )
 131             cerr << "  nr_members: " << nr_members << endl;
 132
 133         vector<long> labels;
 134         for( size_t mi = 0; mi < nr_members; mi++ ) {
 135             long mi_label;
 136             while( (is.peek()) == '#' )
 137                 getline(is,line);
 138             is >> mi_label;
 139             labels.push_back(mi_label);
 140         }
 141         if( verbose >= 3 )
 142             cerr << "  labels: " << labels << endl;
 143
 144         vector<size_t> dims;
 145         for( size_t mi = 0; mi < nr_members; mi++ ) {
 146             size_t mi_dim;
 147             while( (is.peek()) == '#' )
 148                 getline(is,line);
 149             is >> mi_dim;
 150             dims.push_back(mi_dim);
 151         }
 152         if( verbose >= 3 )
 153             cerr << "  dimensions: " << dims << endl;
 154
 155         // add the Factor
 156         vector<Var> Ivars;
 157         Ivars.reserve( nr_members );
 158         for( size_t mi = 0; mi < nr_members; mi++ ) {
 159             map<long,size_t>::iterator vdi = vardims.find( labels[mi] );
 160             if( vdi != vardims.end() ) {
 161                 // check whether dimensions are consistent
 162                 if( vdi->second != dims[mi] )
 163                     DAI_THROWE(INVALID_FACTORGRAPH_FILE,"Variable with label " + boost::lexical_cast<string>(labels[mi]) + " has inconsistent dimensions.");
 164             } else
 165                 vardims[labels[mi]] = dims[mi];
 166             Ivars.push_back( Var(labels[mi], dims[mi]) );
 167         }
 168         facs.push_back( Factor( VarSet( Ivars.begin(), Ivars.end(), Ivars.size() ), (Real)0 ) );
 169
 170         // calculate permutation object
 171         Permute permindex( Ivars );
 172
 173         // read values
 174         size_t nr_nonzeros;
 175         while( (is.peek()) == '#' )
 176             getline(is,line);
 177         is >> nr_nonzeros;
 178         if( verbose >= 3 )
 179             cerr << "  nonzeroes: " << nr_nonzeros << endl;
 180         for( size_t k = 0; k < nr_nonzeros; k++ ) {
 181             size_t li;
 182             Real val;
 183             while( (is.peek()) == '#' )
 184                 getline(is,line);
 185             is >> li;
 186             while( (is.peek()) == '#' )
 187                 getline(is,line);
 188             is >> val;
 189
 190             // store value, but permute indices first according to internal representation
 191             facs.back()[permindex.convertLinearIndex( li )] = val;
 192         }
 193     }
 194
 195     if( verbose >= 3 )
 196         cerr << "factors:" << facs << endl;
 197
 198     fg = FactorGraph(facs);
 199
 200     return is;
 201 }
 202
 203
 204 VarSet FactorGraph::delta( size_t i ) const {
 205     return( Delta(i) / var(i) );
 206 }
 207
 208
 209 VarSet FactorGraph::Delta( size_t i ) const {
 210     // calculate Markov Blanket
 211     VarSet Del;
 212     foreach( const Neighbor &I, nbV(i) ) // for all neighboring factors I of i
 213         foreach( const Neighbor &j, nbF(I) ) // for all neighboring variables j of I
 214             Del |= var(j);
 215
 216     return Del;
 217 }
 218
 219
 220 VarSet FactorGraph::Delta( const VarSet &ns ) const {
 221     VarSet result;
 222     for( VarSet::const_iterator n = ns.begin(); n != ns.end(); n++ )
 223         result |= Delta(findVar(*n));
 224     return result;
 225 }
 226
 227
 228 void FactorGraph::makeCavity( size_t i, bool backup ) {
 229     // fills all Factors that include var(i) with ones
 230     map<size_t,Factor> newFacs;
 231     foreach( const Neighbor &I, nbV(i) ) // for all neighboring factors I of i
 232         newFacs[I] = Factor( factor(I).vars(), (Real)1 );
 233     setFactors( newFacs, backup );
 234 }
 235
 236
 237 void FactorGraph::ReadFromFile( const char *filename ) {
 238     ifstream infile;
 239     infile.open( filename );
 240     if( infile.is_open() ) {
 241         infile >> *this;
 242         infile.close();
 243     } else
 244         DAI_THROWE(CANNOT_READ_FILE,"Cannot read from file " + std::string(filename));
 245 }
 246
 247
 248 void FactorGraph::WriteToFile( const char *filename, size_t precision ) const {
 249     ofstream outfile;
 250     outfile.open( filename );
 251     if( outfile.is_open() ) {
 252         outfile.precision( precision );
 253         outfile << *this;
 254         outfile.close();
 255     } else
 256         DAI_THROWE(CANNOT_WRITE_FILE,"Cannot write to file " + std::string(filename));
 257 }
 258
 259
 260 void FactorGraph::printDot( std::ostream &os ) const {
 261     os << "graph G {" << endl;
 262     os << "node[shape=circle,width=0.4,fixedsize=true];" << endl;
 263     for( size_t i = 0; i < nrVars(); i++ )
 264         os << "\tv" << var(i).label() << ";" << endl;
 265     os << "node[shape=box,width=0.3,height=0.3,fixedsize=true];" << endl;
 266     for( size_t I = 0; I < nrFactors(); I++ )
 267         os << "\tf" << I << ";" << endl;
 268     for( size_t i = 0; i < nrVars(); i++ )
 269         foreach( const Neighbor &I, nbV(i) )  // for all neighboring factors I of i
 270             os << "\tv" << var(i).label() << " -- f" << I << ";" << endl;
 271     os << "}" << endl;
 272 }
 273
 274
 275 vector<VarSet> FactorGraph::Cliques() const {
 276     vector<VarSet> result;
 277
 278     for( size_t I = 0; I < nrFactors(); I++ ) {
 279         bool maximal = true;
 280         for( size_t J = 0; (J < nrFactors()) && maximal; J++ )
 281             if( (factor(J).vars() >> factor(I).vars()) && (factor(J).vars() != factor(I).vars()) )
 282                 maximal = false;
 283
 284         if( maximal )
 285             result.push_back( factor(I).vars() );
 286     }
 287
 288     return result;
 289 }
 290
 291
 292 void FactorGraph::clamp( size_t i, size_t x, bool backup ) {
 293     DAI_ASSERT( x <= var(i).states() );
 294     Factor mask( var(i), (Real)0 );
 295     mask[x] = (Real)1;
 296
 297     map<size_t, Factor> newFacs;
 298     foreach( const Neighbor &I, nbV(i) )
 299         newFacs[I] = factor(I) * mask;
 300     setFactors( newFacs, backup );
 301
 302     return;
 303 }
 304
 305
 306 void FactorGraph::clampVar( size_t i, const vector<size_t> &is, bool backup ) {
 307     Var n = var(i);
 308     Factor mask_n( n, (Real)0 );
 309
 310     foreach( size_t i, is ) {
 311         DAI_ASSERT( i <= n.states() );
 312         mask_n[i] = (Real)1;
 313     }
 314
 315     map<size_t, Factor> newFacs;
 316     foreach( const Neighbor &I, nbV(i) )
 317         newFacs[I] = factor(I) * mask_n;
 318     setFactors( newFacs, backup );
 319 }
 320
 321
 322 void FactorGraph::clampFactor( size_t I, const vector<size_t> &is, bool backup ) {
 323     size_t st = factor(I).states();
 324     Factor newF( factor(I).vars(), (Real)0 );
 325
 326     foreach( size_t i, is ) {
 327         DAI_ASSERT( i <= st );
 328         newF[i] = factor(I)[i];
 329     }
 330
 331     setFactor( I, newF, backup );
 332 }
 333
 334
 335 void FactorGraph::backupFactor( size_t I ) {
 336     map<size_t,Factor>::iterator it = _backup.find( I );
 337     if( it != _backup.end() )
 338         DAI_THROW(MULTIPLE_UNDO);
 339     _backup[I] = factor(I);
 340 }
 341
 342
 343 void FactorGraph::restoreFactor( size_t I ) {
 344     map<size_t,Factor>::iterator it = _backup.find( I );
 345     if( it != _backup.end() ) {
 346         setFactor(I, it->second);
 347         _backup.erase(it);
 348     }
 349 }
 350
 351
 352 void FactorGraph::backupFactors( const VarSet &ns ) {
 353     for( size_t I = 0; I < nrFactors(); I++ )
 354         if( factor(I).vars().intersects( ns ) )
 355             backupFactor( I );
 356 }
 357
 358
 359 void FactorGraph::restoreFactors( const VarSet &ns ) {
 360     map<size_t,Factor> facs;
 361     for( map<size_t,Factor>::iterator uI = _backup.begin(); uI != _backup.end(); ) {
 362         if( factor(uI->first).vars().intersects( ns ) ) {
 363             facs.insert( *uI );
 364             _backup.erase(uI++);
 365         } else
 366             uI++;
 367     }
 368     setFactors( facs );
 369 }
 370
 371
 372 void FactorGraph::restoreFactors() {
 373     setFactors( _backup );
 374     _backup.clear();
 375 }
 376
 377
 378 void FactorGraph::backupFactors( const std::set<size_t> & facs ) {
 379     for( std::set<size_t>::const_iterator fac = facs.begin(); fac != facs.end(); fac++ )
 380         backupFactor( *fac );
 381 }
 382
 383
 384 bool FactorGraph::isPairwise() const {
 385     bool pairwise = true;
 386     for( size_t I = 0; I < nrFactors() && pairwise; I++ )
 387         if( factor(I).vars().size() > 2 )
 388             pairwise = false;
 389     return pairwise;
 390 }
 391
 392
 393 bool FactorGraph::isBinary() const {
 394     bool binary = true;
 395     for( size_t i = 0; i < nrVars() && binary; i++ )
 396         if( var(i).states() > 2 )
 397             binary = false;
 398     return binary;
 399 }
 400
 401
 402 FactorGraph FactorGraph::clamped( size_t i, size_t state ) const {
 403     Var v = var( i );
 404     Real zeroth_order = (Real)1;
 405     vector<Factor> clamped_facs;
 406     for( size_t I = 0; I < nrFactors(); I++ ) {
 407         VarSet v_I = factor(I).vars();
 408         Factor new_factor;
 409         if( v_I.intersects( v ) )
 410             new_factor = factor(I).slice( v, state );
 411         else
 412             new_factor = factor(I);
 413
 414         if( new_factor.vars().size() != 0 ) {
 415             size_t J = 0;
 416             // if it can be merged with a previous one, do that
 417             for( J = 0; J < clamped_facs.size(); J++ )
 418                 if( clamped_facs[J].vars() == new_factor.vars() ) {
 419                     clamped_facs[J] *= new_factor;
 420                     break;
 421                 }
 422             // otherwise, push it back
 423             if( J == clamped_facs.size() || clamped_facs.size() == 0 )
 424                 clamped_facs.push_back( new_factor );
 425         } else
 426             zeroth_order *= new_factor[0];
 427     }
 428     *(clamped_facs.begin()) *= zeroth_order;
 429     return FactorGraph( clamped_facs );
 430 }
 431
 432
 433 FactorGraph FactorGraph::maximalFactors() const {
 434     vector<size_t> maxfac( nrFactors() );
 435     map<size_t,size_t> newindex;
 436     size_t nrmax = 0;
 437     for( size_t I = 0; I < nrFactors(); I++ ) {
 438         maxfac[I] = I;
 439         VarSet maxfacvars = factor(maxfac[I]).vars();
 440         for( size_t J = 0; J < nrFactors(); J++ ) {
 441             VarSet Jvars = factor(J).vars();
 442             if( Jvars >> maxfacvars && (Jvars != maxfacvars) ) {
 443                 maxfac[I] = J;
 444                 maxfacvars = factor(maxfac[I]).vars();
 445             }
 446         }
 447         if( maxfac[I] == I )
 448             newindex[I] = nrmax++;
 449     }
 450
 451     vector<Factor> facs( nrmax );
 452     for( size_t I = 0; I < nrFactors(); I++ )
 453         facs[newindex[maxfac[I]]] *= factor(I);
 454
 455     return FactorGraph( facs.begin(), facs.end(), vars().begin(), vars().end(), facs.size(), nrVars() );
 456 }
 457
 458
 459 } // end of namespace dai