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