3cff7eb2c942eae4e8d9b78041c73adb03e46c43
[libdai.git] / 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