1 /* This file is part of libDAI - http://www.libdai.org/
3 * Copyright (c) 2006-2011, The libDAI authors. All rights reserved.
5 * Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
18 #include <dai/factorgraph.h>
20 #include <dai/exceptions.h>
21 #include <boost/lexical_cast.hpp>
30 FactorGraph::FactorGraph( const std::vector
<Factor
> &P
) : _G(), _backup() {
31 // add factors, obtain variables
33 _factors
.reserve( P
.size() );
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();
42 _vars
.reserve( varset
.size() );
43 for( set
<Var
>::const_iterator p1
= varset
.begin(); p1
!= varset
.end(); p1
++ )
44 _vars
.push_back( *p1
);
46 // create graph structure
47 constructGraph( nrEdges
);
51 void FactorGraph::constructGraph( size_t nrEdges
) {
52 // create a mapping for indices
53 hash_map
<size_t, size_t> hashmap
;
55 for( size_t i
= 0; i
< vars().size(); i
++ )
56 hashmap
[var(i
).label()] = i
;
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
) );
67 // create bipartite graph
68 _G
.construct( nrVars(), nrFactors(), edges
.begin(), edges
.end() );
72 /// Writes a FactorGraph to an output stream
73 std::ostream
& operator<< ( std::ostream
&os
, const FactorGraph
&fg
) {
74 os
<< fg
.nrFactors() << endl
;
76 for( size_t I
= 0; I
< fg
.nrFactors(); I
++ ) {
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() << " ";
82 for( VarSet::const_iterator i
= fg
.factor(I
).vars().begin(); i
!= fg
.factor(I
).vars().end(); i
++ )
83 os
<< i
->states() << " ";
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 )
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
;
99 /// Reads a FactorGraph from an input stream
100 std::istream
& operator>> ( std::istream
& is
, FactorGraph
&fg
) {
107 while( (is
.peek()) == '#' )
111 DAI_THROWE(INVALID_FACTORGRAPH_FILE
,"Cannot read number of factors");
113 cerr
<< "Reading " << nr_Factors
<< " factors..." << endl
;
116 if( is
.fail() || line
.size() > 0 )
117 DAI_THROWE(INVALID_FACTORGRAPH_FILE
,"Expecting empty line");
119 map
<long,size_t> vardims
;
120 for( size_t I
= 0; I
< nr_Factors
; I
++ ) {
122 cerr
<< "Reading factor " << I
<< "..." << endl
;
124 while( (is
.peek()) == '#' )
128 cerr
<< " nr_members: " << nr_members
<< endl
;
131 for( size_t mi
= 0; mi
< nr_members
; mi
++ ) {
133 while( (is
.peek()) == '#' )
136 labels
.push_back(mi_label
);
139 cerr
<< " labels: " << labels
<< endl
;
142 for( size_t mi
= 0; mi
< nr_members
; mi
++ ) {
144 while( (is
.peek()) == '#' )
147 dims
.push_back(mi_dim
);
150 cerr
<< " dimensions: " << dims
<< endl
;
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.");
162 vardims
[labels
[mi
]] = dims
[mi
];
163 Ivars
.push_back( Var(labels
[mi
], dims
[mi
]) );
165 facs
.push_back( Factor( VarSet( Ivars
.begin(), Ivars
.end(), Ivars
.size() ), (Real
)0 ) );
167 cerr
<< " vardims: " << vardims
<< endl
;
169 // calculate permutation object
170 Permute
permindex( Ivars
);
174 while( (is
.peek()) == '#' )
178 cerr
<< " nonzeroes: " << nr_nonzeros
<< endl
;
179 for( size_t k
= 0; k
< nr_nonzeros
; k
++ ) {
182 while( (is
.peek()) == '#' )
185 while( (is
.peek()) == '#' )
189 // store value, but permute indices first according to internal representation
190 facs
.back().set( permindex
.convertLinearIndex( li
), val
);
195 cerr
<< "factors:" << facs
<< endl
;
197 fg
= FactorGraph(facs
);
203 VarSet
FactorGraph::Delta( size_t i
) const {
204 // calculate Markov Blanket
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
214 VarSet
FactorGraph::Delta( const VarSet
&ns
) const {
216 for( VarSet::const_iterator n
= ns
.begin(); n
!= ns
.end(); n
++ )
217 result
|= Delta( findVar(*n
) );
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
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
);
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
);
250 void FactorGraph::ReadFromFile( const char *filename
) {
252 infile
.open( filename
);
253 if( infile
.is_open() ) {
257 DAI_THROWE(CANNOT_READ_FILE
,"Cannot read from file " + std::string(filename
));
261 void FactorGraph::WriteToFile( const char *filename
, size_t precision
) const {
263 outfile
.open( filename
);
264 if( outfile
.is_open() ) {
265 outfile
.precision( precision
);
269 DAI_THROWE(CANNOT_WRITE_FILE
,"Cannot write to file " + std::string(filename
));
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
;
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
) )
294 G
.addEdge( i
, j
, true );
299 bool FactorGraph::isMaximal( size_t I
) const {
300 const VarSet
& I_vars
= factor(I
).vars();
301 size_t I_size
= I_vars
.size();
304 for( size_t J
= 0; J
< nrFactors(); J
++ )
306 if( factor(J
).vars().size() > 0 )
310 bforeach( const Neighbor
& i
, nbF(I
) ) {
311 bforeach( const Neighbor
& J
, nbV(i
) ) {
313 if( (factor(J
).vars() >> I_vars
) && (factor(J
).vars().size() != I_size
) )
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();
327 for( size_t J
= 0; J
< nrFactors(); J
++ )
329 if( factor(J
).vars().size() > 0 )
330 return maximalFactor( J
);
333 bforeach( const Neighbor
& i
, nbF(I
) ) {
334 bforeach( const Neighbor
& J
, nbV(i
) ) {
336 if( (factor(J
).vars() >> I_vars
) && (factor(J
).vars().size() != I_size
) )
337 return maximalFactor( J
);
345 vector
<VarSet
> FactorGraph::maximalFactorDomains() const {
346 vector
<VarSet
> result
;
348 for( size_t I
= 0; I
< nrFactors(); I
++ )
350 result
.push_back( factor(I
).vars() );
352 if( result
.size() == 0 )
353 result
.push_back( VarSet() );
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
];
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
369 for( size_t I
= 0; I
< nrFactors(); I
++ )
370 lS
+= dai::log( factor(I
)[BigInt_size_t(S(factor(I
).vars()))] );
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 );
380 map
<size_t, Factor
> newFacs
;
381 bforeach( const Neighbor
&I
, nbV(i
) )
382 newFacs
[I
] = factor(I
) * mask
;
383 setFactors( newFacs
, backup
);
389 void FactorGraph::clampVar( size_t i
, const vector
<size_t> &is
, bool backup
) {
391 Factor
mask_n( n
, (Real
)0 );
393 bforeach( size_t i
, is
) {
394 DAI_ASSERT( i
<= n
.states() );
395 mask_n
.set( i
, (Real
)1 );
398 map
<size_t, Factor
> newFacs
;
399 bforeach( const Neighbor
&I
, nbV(i
) )
400 newFacs
[I
] = factor(I
) * mask_n
;
401 setFactors( newFacs
, backup
);
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 );
409 bforeach( size_t i
, is
) {
410 DAI_ASSERT( i
<= st
);
411 newF
.set( i
, factor(I
)[i
] );
414 setFactor( I
, newF
, backup
);
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
);
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
);
432 DAI_THROW(OBJECT_NOT_FOUND
);
436 void FactorGraph::backupFactors( const VarSet
&ns
) {
437 for( size_t I
= 0; I
< nrFactors(); I
++ )
438 if( factor(I
).vars().intersects( ns
) )
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
) ) {
456 void FactorGraph::restoreFactors() {
457 setFactors( _backup
);
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
);
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 )
477 bool FactorGraph::isBinary() const {
479 for( size_t i
= 0; i
< nrVars() && binary
; i
++ )
480 if( var(i
).states() > 2 )
486 FactorGraph
FactorGraph::clamped( size_t i
, size_t state
) const {
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();
494 if( v_I
.intersects( v
) )
495 new_factor
= factor(I
).slice( v
, state
);
497 new_factor
= factor(I
);
499 if( new_factor
.vars().size() != 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
;
507 // otherwise, push it back
508 if( J
== clamped_facs
.size() || clamped_facs
.size() == 0 )
509 clamped_facs
.push_back( new_factor
);
511 zeroth_order
*= new_factor
[0];
513 *(clamped_facs
.begin()) *= zeroth_order
;
514 return FactorGraph( clamped_facs
);
518 FactorGraph
FactorGraph::maximalFactors() const {
519 vector
<size_t> maxfac( nrFactors() );
520 map
<size_t,size_t> newindex
;
522 for( size_t I
= 0; I
< nrFactors(); 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
) ) {
529 maxfacvars
= factor(maxfac
[I
]).vars();
533 newindex
[I
] = nrmax
++;
536 vector
<Factor
> facs( nrmax
);
537 for( size_t I
= 0; I
< nrFactors(); I
++ )
538 facs
[newindex
[maxfac
[I
]]] *= factor(I
);
540 return FactorGraph( facs
.begin(), facs
.end(), vars().begin(), vars().end(), facs
.size(), nrVars() );
544 } // end of namespace dai