Updated copyrights
[libdai.git] / src / factorgraph.cpp
1 /* Copyright (C) 2006-2008 Joris Mooij [joris dot mooij at tuebingen dot mpg dot de]
2 Radboud University Nijmegen, The Netherlands /
3 Max Planck Institute for Biological Cybernetics, Germany
4
5 This file is part of libDAI.
6
7 libDAI is free software; you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation; either version 2 of the License, or
10 (at your option) any later version.
11
12 libDAI is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with libDAI; if not, write to the Free Software
19 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
20 */
21
22
23 #include <iostream>
24 #include <iterator>
25 #include <map>
26 #include <set>
27 #include <fstream>
28 #include <string>
29 #include <algorithm>
30 #include <functional>
31 #include <dai/factorgraph.h>
32 #include <dai/util.h>
33 #include <dai/exceptions.h>
34
35
36 namespace dai {
37
38
39 using namespace std;
40
41
42 FactorGraph::FactorGraph( const std::vector<Factor> &P ) : G(), _backup() {
43 // add factors, obtain variables
44 set<Var> _vars;
45 _factors.reserve( P.size() );
46 size_t nrEdges = 0;
47 for( vector<Factor>::const_iterator p2 = P.begin(); p2 != P.end(); p2++ ) {
48 _factors.push_back( *p2 );
49 copy( p2->vars().begin(), p2->vars().end(), inserter( _vars, _vars.begin() ) );
50 nrEdges += p2->vars().size();
51 }
52
53 // add _vars
54 vars.reserve( _vars.size() );
55 for( set<Var>::const_iterator p1 = _vars.begin(); p1 != _vars.end(); p1++ )
56 vars.push_back( *p1 );
57
58 // create graph structure
59 constructGraph( nrEdges );
60 }
61
62
63 /// Part of constructors (creates edges, neighbours and adjacency matrix)
64 void FactorGraph::constructGraph( size_t nrEdges ) {
65 // create a mapping for indices
66 hash_map<size_t, size_t> hashmap;
67
68 for( size_t i = 0; i < vars.size(); i++ )
69 hashmap[var(i).label()] = i;
70
71 // create edge list
72 vector<Edge> edges;
73 edges.reserve( nrEdges );
74 for( size_t i2 = 0; i2 < nrFactors(); i2++ ) {
75 const VarSet& ns = factor(i2).vars();
76 for( VarSet::const_iterator q = ns.begin(); q != ns.end(); q++ )
77 edges.push_back( Edge(hashmap[q->label()], i2) );
78 }
79
80 // create bipartite graph
81 G.construct( nrVars(), nrFactors(), edges.begin(), edges.end() );
82 }
83
84
85 ostream& operator << (ostream& os, const FactorGraph& fg) {
86 os << fg.nrFactors() << endl;
87
88 for( size_t I = 0; I < fg.nrFactors(); I++ ) {
89 os << endl;
90 os << fg.factor(I).vars().size() << endl;
91 for( VarSet::const_iterator i = fg.factor(I).vars().begin(); i != fg.factor(I).vars().end(); i++ )
92 os << i->label() << " ";
93 os << endl;
94 for( VarSet::const_iterator i = fg.factor(I).vars().begin(); i != fg.factor(I).vars().end(); i++ )
95 os << i->states() << " ";
96 os << endl;
97 size_t nr_nonzeros = 0;
98 for( size_t k = 0; k < fg.factor(I).states(); k++ )
99 if( fg.factor(I)[k] != 0.0 )
100 nr_nonzeros++;
101 os << nr_nonzeros << endl;
102 for( size_t k = 0; k < fg.factor(I).states(); k++ )
103 if( fg.factor(I)[k] != 0.0 ) {
104 char buf[20];
105 sprintf(buf,"%18.14g", fg.factor(I)[k]);
106 os << k << " " << buf << endl;
107 }
108 }
109
110 return(os);
111 }
112
113
114 istream& operator >> (istream& is, FactorGraph& fg) {
115 long verbose = 0;
116
117 try {
118 vector<Factor> facs;
119 size_t nr_Factors;
120 string line;
121
122 while( (is.peek()) == '#' )
123 getline(is,line);
124 is >> nr_Factors;
125 if( is.fail() )
126 DAI_THROW(INVALID_FACTORGRAPH_FILE);
127 if( verbose >= 2 )
128 cout << "Reading " << nr_Factors << " factors..." << endl;
129
130 getline (is,line);
131 if( is.fail() )
132 DAI_THROW(INVALID_FACTORGRAPH_FILE);
133
134 map<long,size_t> vardims;
135 for( size_t I = 0; I < nr_Factors; I++ ) {
136 if( verbose >= 3 )
137 cout << "Reading factor " << I << "..." << endl;
138 size_t nr_members;
139 while( (is.peek()) == '#' )
140 getline(is,line);
141 is >> nr_members;
142 if( verbose >= 3 )
143 cout << " nr_members: " << nr_members << endl;
144
145 vector<long> labels;
146 for( size_t mi = 0; mi < nr_members; mi++ ) {
147 long mi_label;
148 while( (is.peek()) == '#' )
149 getline(is,line);
150 is >> mi_label;
151 labels.push_back(mi_label);
152 }
153 if( verbose >= 3 )
154 cout << " labels: " << labels << endl;
155
156 vector<size_t> dims;
157 for( size_t mi = 0; mi < nr_members; mi++ ) {
158 size_t mi_dim;
159 while( (is.peek()) == '#' )
160 getline(is,line);
161 is >> mi_dim;
162 dims.push_back(mi_dim);
163 }
164 if( verbose >= 3 )
165 cout << " dimensions: " << dims << endl;
166
167 // add the Factor
168 VarSet I_vars;
169 for( size_t mi = 0; mi < nr_members; mi++ ) {
170 map<long,size_t>::iterator vdi = vardims.find( labels[mi] );
171 if( vdi != vardims.end() ) {
172 // check whether dimensions are consistent
173 if( vdi->second != dims[mi] )
174 DAI_THROW(INVALID_FACTORGRAPH_FILE);
175 } else
176 vardims[labels[mi]] = dims[mi];
177 I_vars |= Var(labels[mi], dims[mi]);
178 }
179 facs.push_back( Factor( I_vars, 0.0 ) );
180
181 // calculate permutation sigma (internally, members are sorted)
182 vector<size_t> sigma(nr_members,0);
183 VarSet::iterator j = I_vars.begin();
184 for( size_t mi = 0; mi < nr_members; mi++,j++ ) {
185 long search_for = j->label();
186 vector<long>::iterator j_loc = find(labels.begin(),labels.end(),search_for);
187 sigma[mi] = j_loc - labels.begin();
188 }
189 if( verbose >= 3 )
190 cout << " sigma: " << sigma << endl;
191
192 // calculate multindices
193 Permute permindex( dims, sigma );
194
195 // read values
196 size_t nr_nonzeros;
197 while( (is.peek()) == '#' )
198 getline(is,line);
199 is >> nr_nonzeros;
200 if( verbose >= 3 )
201 cout << " nonzeroes: " << nr_nonzeros << endl;
202 for( size_t k = 0; k < nr_nonzeros; k++ ) {
203 size_t li;
204 double val;
205 while( (is.peek()) == '#' )
206 getline(is,line);
207 is >> li;
208 while( (is.peek()) == '#' )
209 getline(is,line);
210 is >> val;
211
212 // store value, but permute indices first according
213 // to internal representation
214 facs.back()[permindex.convert_linear_index( li )] = val;
215 }
216 }
217
218 if( verbose >= 3 )
219 cout << "factors:" << facs << endl;
220
221 fg = FactorGraph(facs);
222 } catch (char *e) {
223 cout << e << endl;
224 }
225
226 return is;
227 }
228
229
230 VarSet FactorGraph::delta( unsigned i ) const {
231 return( Delta(i) / var(i) );
232 }
233
234
235 VarSet FactorGraph::Delta( unsigned i ) const {
236 // calculate Markov Blanket
237 VarSet Del;
238 foreach( const Neighbor &I, nbV(i) ) // for all neighboring factors I of i
239 foreach( const Neighbor &j, nbF(I) ) // for all neighboring variables j of I
240 Del |= var(j);
241
242 return Del;
243 }
244
245
246 VarSet FactorGraph::Delta( const VarSet &ns ) const {
247 VarSet result;
248 for( VarSet::const_iterator n = ns.begin(); n != ns.end(); n++ )
249 result |= Delta(findVar(*n));
250 return result;
251 }
252
253
254 void FactorGraph::makeCavity( unsigned i, bool backup ) {
255 // fills all Factors that include var(i) with ones
256 map<size_t,Factor> newFacs;
257 foreach( const Neighbor &I, nbV(i) ) // for all neighboring factors I of i
258 newFacs[I] = Factor(factor(I).vars(), 1.0);
259 setFactors( newFacs, backup );
260 }
261
262
263 void FactorGraph::ReadFromFile( const char *filename ) {
264 ifstream infile;
265 infile.open( filename );
266 if( infile.is_open() ) {
267 infile >> *this;
268 infile.close();
269 } else
270 DAI_THROW(CANNOT_READ_FILE);
271 }
272
273
274 void FactorGraph::WriteToFile( const char *filename ) const {
275 ofstream outfile;
276 outfile.open( filename );
277 if( outfile.is_open() ) {
278 outfile << *this;
279 outfile.close();
280 } else
281 DAI_THROW(CANNOT_WRITE_FILE);
282 }
283
284
285 void FactorGraph::printDot( std::ostream &os ) const {
286 os << "graph G {" << endl;
287 os << "node[shape=circle,width=0.4,fixedsize=true];" << endl;
288 for( size_t i = 0; i < nrVars(); i++ )
289 os << "\tv" << var(i).label() << ";" << endl;
290 os << "node[shape=box,width=0.3,height=0.3,fixedsize=true];" << endl;
291 for( size_t I = 0; I < nrFactors(); I++ )
292 os << "\tf" << I << ";" << endl;
293 for( size_t i = 0; i < nrVars(); i++ )
294 foreach( const Neighbor &I, nbV(i) ) // for all neighboring factors I of i
295 os << "\tv" << var(i).label() << " -- f" << I << ";" << endl;
296 os << "}" << endl;
297 }
298
299
300 vector<VarSet> FactorGraph::Cliques() const {
301 vector<VarSet> result;
302
303 for( size_t I = 0; I < nrFactors(); I++ ) {
304 bool maximal = true;
305 for( size_t J = 0; (J < nrFactors()) && maximal; J++ )
306 if( (factor(J).vars() >> factor(I).vars()) && (factor(J).vars() != factor(I).vars()) )
307 maximal = false;
308
309 if( maximal )
310 result.push_back( factor(I).vars() );
311 }
312
313 return result;
314 }
315
316
317 void FactorGraph::clamp( const Var & n, size_t i, bool backup ) {
318 assert( i <= n.states() );
319
320 // Multiply each factor that contains the variable with a delta function
321
322 Factor delta_n_i(n,0.0);
323 delta_n_i[i] = 1.0;
324
325 map<size_t, Factor> newFacs;
326 // For all factors that contain n
327 for( size_t I = 0; I < nrFactors(); I++ )
328 if( factor(I).vars().contains( n ) )
329 // Multiply it with a delta function
330 newFacs[I] = factor(I) * delta_n_i;
331 setFactors( newFacs, backup );
332
333 return;
334 }
335
336
337 void FactorGraph::backupFactor( size_t I ) {
338 map<size_t,Factor>::iterator it = _backup.find( I );
339 if( it != _backup.end() )
340 DAI_THROW( MULTIPLE_UNDO );
341 _backup[I] = factor(I);
342 }
343
344
345 void FactorGraph::restoreFactor( size_t I ) {
346 map<size_t,Factor>::iterator it = _backup.find( I );
347 if( it != _backup.end() ) {
348 setFactor(I, it->second);
349 _backup.erase(it);
350 }
351 }
352
353
354 void FactorGraph::backupFactors( const VarSet &ns ) {
355 for( size_t I = 0; I < nrFactors(); I++ )
356 if( factor(I).vars().intersects( ns ) )
357 backupFactor( I );
358 }
359
360
361 void FactorGraph::restoreFactors( const VarSet &ns ) {
362 map<size_t,Factor> facs;
363 for( map<size_t,Factor>::iterator uI = _backup.begin(); uI != _backup.end(); ) {
364 if( factor(uI->first).vars().intersects( ns ) ) {
365 facs.insert( *uI );
366 _backup.erase(uI++);
367 } else
368 uI++;
369 }
370 setFactors( facs );
371 }
372
373
374 void FactorGraph::restoreFactors() {
375 setFactors( _backup );
376 _backup.clear();
377 }
378
379 void FactorGraph::backupFactors( const std::set<size_t> & facs ) {
380 for( std::set<size_t>::const_iterator fac = facs.begin(); fac != facs.end(); fac++ )
381 backupFactor( *fac );
382 }
383
384
385 bool FactorGraph::isPairwise() const {
386 bool pairwise = true;
387 for( size_t I = 0; I < nrFactors() && pairwise; I++ )
388 if( factor(I).vars().size() > 2 )
389 pairwise = false;
390 return pairwise;
391 }
392
393
394 bool FactorGraph::isBinary() const {
395 bool binary = true;
396 for( size_t i = 0; i < nrVars() && binary; i++ )
397 if( var(i).states() > 2 )
398 binary = false;
399 return binary;
400 }
401
402
403 FactorGraph FactorGraph::clamped( const Var & v_i, size_t state ) const {
404 Real zeroth_order = 1.0;
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_i ) )
410 new_factor = factor(I).slice( v_i, 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