src/factorgraph.cpp

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