Some small documentation updates
[libdai.git] / 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).nrStates(); 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).nrStates(); 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() || line.size() > 0 )
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().set( 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 // calculate Markov Blanket
206 VarSet Del;
207 foreach( const Neighbor &I, nbV(i) ) // for all neighboring factors I of i
208 foreach( const Neighbor &j, nbF(I) ) // for all neighboring variables j of I
209 Del |= var(j);
210
211 return Del;
212 }
213
214
215 VarSet FactorGraph::Delta( const VarSet &ns ) const {
216 VarSet result;
217 for( VarSet::const_iterator n = ns.begin(); n != ns.end(); n++ )
218 result |= Delta( findVar(*n) );
219 return result;
220 }
221
222
223 void FactorGraph::makeCavity( size_t i, bool backup ) {
224 // fills all Factors that include var(i) with ones
225 map<size_t,Factor> newFacs;
226 foreach( const Neighbor &I, nbV(i) ) // for all neighboring factors I of i
227 newFacs[I] = Factor( factor(I).vars(), (Real)1 );
228 setFactors( newFacs, backup );
229 }
230
231
232 void FactorGraph::ReadFromFile( const char *filename ) {
233 ifstream infile;
234 infile.open( filename );
235 if( infile.is_open() ) {
236 infile >> *this;
237 infile.close();
238 } else
239 DAI_THROWE(CANNOT_READ_FILE,"Cannot read from file " + std::string(filename));
240 }
241
242
243 void FactorGraph::WriteToFile( const char *filename, size_t precision ) const {
244 ofstream outfile;
245 outfile.open( filename );
246 if( outfile.is_open() ) {
247 outfile.precision( precision );
248 outfile << *this;
249 outfile.close();
250 } else
251 DAI_THROWE(CANNOT_WRITE_FILE,"Cannot write to file " + std::string(filename));
252 }
253
254
255 void FactorGraph::printDot( std::ostream &os ) const {
256 os << "graph FactorGraph {" << endl;
257 os << "node[shape=circle,width=0.4,fixedsize=true];" << endl;
258 for( size_t i = 0; i < nrVars(); i++ )
259 os << "\tv" << var(i).label() << ";" << endl;
260 os << "node[shape=box,width=0.3,height=0.3,fixedsize=true];" << endl;
261 for( size_t I = 0; I < nrFactors(); I++ )
262 os << "\tf" << I << ";" << endl;
263 for( size_t i = 0; i < nrVars(); i++ )
264 foreach( const Neighbor &I, nbV(i) ) // for all neighboring factors I of i
265 os << "\tv" << var(i).label() << " -- f" << I << ";" << endl;
266 os << "}" << endl;
267 }
268
269
270 GraphAL FactorGraph::MarkovGraph() const {
271 GraphAL G( nrVars() );
272 for( size_t i = 0; i < nrVars(); i++ )
273 foreach( const Neighbor &I, nbV(i) )
274 foreach( const Neighbor &j, nbF(I) )
275 if( i < j )
276 G.addEdge( i, j, true );
277 return G;
278 }
279
280
281 bool FactorGraph::isMaximal( size_t I ) const {
282 const VarSet& I_vars = factor(I).vars();
283 size_t I_size = I_vars.size();
284
285 if( I_size == 0 ) {
286 for( size_t J = 0; J < nrFactors(); J++ )
287 if( J != I )
288 if( factor(J).vars().size() > 0 )
289 return false;
290 return true;
291 } else {
292 foreach( const Neighbor& i, nbF(I) ) {
293 foreach( const Neighbor& J, nbV(i) ) {
294 if( J != I )
295 if( (factor(J).vars() >> I_vars) && (factor(J).vars().size() != I_size) )
296 return false;
297 }
298 }
299 return true;
300 }
301 }
302
303
304 size_t FactorGraph::maximalFactor( size_t I ) const {
305 const VarSet& I_vars = factor(I).vars();
306 size_t I_size = I_vars.size();
307
308 if( I_size == 0 ) {
309 for( size_t J = 0; J < nrFactors(); J++ )
310 if( J != I )
311 if( factor(J).vars().size() > 0 )
312 return maximalFactor( J );
313 return I;
314 } else {
315 foreach( const Neighbor& i, nbF(I) ) {
316 foreach( const Neighbor& J, nbV(i) ) {
317 if( J != I )
318 if( (factor(J).vars() >> I_vars) && (factor(J).vars().size() != I_size) )
319 return maximalFactor( J );
320 }
321 }
322 return I;
323 }
324 }
325
326
327 vector<VarSet> FactorGraph::maximalFactorDomains() const {
328 vector<VarSet> result;
329
330 for( size_t I = 0; I < nrFactors(); I++ )
331 if( isMaximal( I ) )
332 result.push_back( factor(I).vars() );
333
334 if( result.size() == 0 )
335 result.push_back( VarSet() );
336 return result;
337 }
338
339
340 Real FactorGraph::logScore( const std::vector<size_t>& statevec ) {
341 // Construct a State object that represents statevec
342 // This decouples the representation of the joint state in statevec from the factor graph
343 map<Var, size_t> statemap;
344 for( size_t i = 0; i < statevec.size(); i++ )
345 statemap[var(i)] = statevec[i];
346 State S(statemap);
347
348 // Evaluate the log probability of the joint configuration in statevec
349 // by summing the log factor entries of the factors that correspond to this joint configuration
350 Real lS = 0.0;
351 for( size_t I = 0; I < nrFactors(); I++ )
352 lS += dai::log( factor(I)[S(factor(I).vars())] );
353 return lS;
354 }
355
356
357 void FactorGraph::clamp( size_t i, size_t x, bool backup ) {
358 DAI_ASSERT( x <= var(i).states() );
359 Factor mask( var(i), (Real)0 );
360 mask.set( x, (Real)1 );
361
362 map<size_t, Factor> newFacs;
363 foreach( const Neighbor &I, nbV(i) )
364 newFacs[I] = factor(I) * mask;
365 setFactors( newFacs, backup );
366
367 return;
368 }
369
370
371 void FactorGraph::clampVar( size_t i, const vector<size_t> &is, bool backup ) {
372 Var n = var(i);
373 Factor mask_n( n, (Real)0 );
374
375 foreach( size_t i, is ) {
376 DAI_ASSERT( i <= n.states() );
377 mask_n.set( i, (Real)1 );
378 }
379
380 map<size_t, Factor> newFacs;
381 foreach( const Neighbor &I, nbV(i) )
382 newFacs[I] = factor(I) * mask_n;
383 setFactors( newFacs, backup );
384 }
385
386
387 void FactorGraph::clampFactor( size_t I, const vector<size_t> &is, bool backup ) {
388 size_t st = factor(I).nrStates();
389 Factor newF( factor(I).vars(), (Real)0 );
390
391 foreach( size_t i, is ) {
392 DAI_ASSERT( i <= st );
393 newF.set( i, factor(I)[i] );
394 }
395
396 setFactor( I, newF, backup );
397 }
398
399
400 void FactorGraph::backupFactor( size_t I ) {
401 map<size_t,Factor>::iterator it = _backup.find( I );
402 if( it != _backup.end() )
403 DAI_THROW(MULTIPLE_UNDO);
404 _backup[I] = factor(I);
405 }
406
407
408 void FactorGraph::restoreFactor( size_t I ) {
409 map<size_t,Factor>::iterator it = _backup.find( I );
410 if( it != _backup.end() ) {
411 setFactor(I, it->second);
412 _backup.erase(it);
413 } else
414 DAI_THROW(OBJECT_NOT_FOUND);
415 }
416
417
418 void FactorGraph::backupFactors( const VarSet &ns ) {
419 for( size_t I = 0; I < nrFactors(); I++ )
420 if( factor(I).vars().intersects( ns ) )
421 backupFactor( I );
422 }
423
424
425 void FactorGraph::restoreFactors( const VarSet &ns ) {
426 map<size_t,Factor> facs;
427 for( map<size_t,Factor>::iterator uI = _backup.begin(); uI != _backup.end(); ) {
428 if( factor(uI->first).vars().intersects( ns ) ) {
429 facs.insert( *uI );
430 _backup.erase(uI++);
431 } else
432 uI++;
433 }
434 setFactors( facs );
435 }
436
437
438 void FactorGraph::restoreFactors() {
439 setFactors( _backup );
440 _backup.clear();
441 }
442
443
444 void FactorGraph::backupFactors( const std::set<size_t> & facs ) {
445 for( std::set<size_t>::const_iterator fac = facs.begin(); fac != facs.end(); fac++ )
446 backupFactor( *fac );
447 }
448
449
450 bool FactorGraph::isPairwise() const {
451 bool pairwise = true;
452 for( size_t I = 0; I < nrFactors() && pairwise; I++ )
453 if( factor(I).vars().size() > 2 )
454 pairwise = false;
455 return pairwise;
456 }
457
458
459 bool FactorGraph::isBinary() const {
460 bool binary = true;
461 for( size_t i = 0; i < nrVars() && binary; i++ )
462 if( var(i).states() > 2 )
463 binary = false;
464 return binary;
465 }
466
467
468 FactorGraph FactorGraph::clamped( size_t i, size_t state ) const {
469 Var v = var( i );
470 Real zeroth_order = (Real)1;
471 vector<Factor> clamped_facs;
472 clamped_facs.push_back( createFactorDelta( v, state ) );
473 for( size_t I = 0; I < nrFactors(); I++ ) {
474 VarSet v_I = factor(I).vars();
475 Factor new_factor;
476 if( v_I.intersects( v ) )
477 new_factor = factor(I).slice( v, state );
478 else
479 new_factor = factor(I);
480
481 if( new_factor.vars().size() != 0 ) {
482 size_t J = 0;
483 // if it can be merged with a previous one, do that
484 for( J = 0; J < clamped_facs.size(); J++ )
485 if( clamped_facs[J].vars() == new_factor.vars() ) {
486 clamped_facs[J] *= new_factor;
487 break;
488 }
489 // otherwise, push it back
490 if( J == clamped_facs.size() || clamped_facs.size() == 0 )
491 clamped_facs.push_back( new_factor );
492 } else
493 zeroth_order *= new_factor[0];
494 }
495 *(clamped_facs.begin()) *= zeroth_order;
496 return FactorGraph( clamped_facs );
497 }
498
499
500 FactorGraph FactorGraph::maximalFactors() const {
501 vector<size_t> maxfac( nrFactors() );
502 map<size_t,size_t> newindex;
503 size_t nrmax = 0;
504 for( size_t I = 0; I < nrFactors(); I++ ) {
505 maxfac[I] = I;
506 VarSet maxfacvars = factor(maxfac[I]).vars();
507 for( size_t J = 0; J < nrFactors(); J++ ) {
508 VarSet Jvars = factor(J).vars();
509 if( Jvars >> maxfacvars && (Jvars != maxfacvars) ) {
510 maxfac[I] = J;
511 maxfacvars = factor(maxfac[I]).vars();
512 }
513 }
514 if( maxfac[I] == I )
515 newindex[I] = nrmax++;
516 }
517
518 vector<Factor> facs( nrmax );
519 for( size_t I = 0; I < nrFactors(); I++ )
520 facs[newindex[maxfac[I]]] *= factor(I);
521
522 return FactorGraph( facs.begin(), facs.end(), vars().begin(), vars().end(), facs.size(), nrVars() );
523 }
524
525
526 } // end of namespace dai