311bcb8b1151c511c5d16e90c16d38f6533764ee
[libdai.git] / utils / fginfo.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-2010 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 <cstdlib>
14 #include <dai/factorgraph.h>
15 #include <dai/jtree.h>
16
17
18 using namespace std;
19 using namespace dai;
20
21
22 void findLoopClusters( const FactorGraph & fg, std::set<VarSet> &allcl, VarSet newcl, const Var & root, size_t length, VarSet vars ) {
23 for( VarSet::const_iterator in = vars.begin(); in != vars.end(); in++ ) {
24 VarSet ind = fg.delta( *in );
25 if( (newcl.size()) >= 2 && ind.contains( root ) ) {
26 allcl.insert( newcl | *in );
27 }
28 else if( length > 1 )
29 findLoopClusters( fg, allcl, newcl | *in, root, length - 1, ind / newcl );
30 }
31 }
32
33
34 size_t countLoops( const FactorGraph & fg, size_t loopdepth ) {
35 set<VarSet> loops;
36 for( vector<Var>::const_iterator i0 = fg.vars().begin(); i0 != fg.vars().end(); i0++ ) {
37 VarSet i0d = fg.delta(*i0);
38 if( loopdepth > 1 )
39 findLoopClusters( fg, loops, *i0, *i0, loopdepth - 1, i0d );
40 }
41 return loops.size();
42 }
43
44
45 bool hasShortLoops( const std::vector<Factor> &P ) {
46 bool found = false;
47 vector<Factor>::const_iterator I, J;
48 for( I = P.begin(); I != P.end(); I++ ) {
49 J = I;
50 J++;
51 for( ; J != P.end(); J++ )
52 if( (I->vars() & J->vars()).size() >= 2 ) {
53 found = true;
54 break;
55 }
56 if( found )
57 break;
58 }
59 return found;
60 }
61
62
63 bool hasNegatives( const std::vector<Factor> &P ) {
64 bool found = false;
65 for( size_t I = 0; I < P.size(); I++ )
66 if( P[I].hasNegatives() ) {
67 found = true;
68 break;
69 }
70 return found;
71 }
72
73
74 int main( int argc, char *argv[] ) {
75 if( argc != 3 ) {
76 // Display help message if number of command line arguments is incorrect
77 cout << "This program is part of libDAI - http://www.libdai.org/" << endl << endl;
78 cout << "Usage: ./fginfo <in.fg> <maxstates>" << endl << endl;
79 cout << "Reports some detailed information about the factor graph <in.fg>." << endl;
80 cout << "Also calculates treewidth, with maximum total number of states" << endl;
81 cout << "given by <maxstates>, where 0 means unlimited." << endl << endl;
82 return 1;
83 } else {
84 // Read factorgraph
85 FactorGraph fg;
86 char *infile = argv[1];
87 size_t maxstates = fromString<size_t>( argv[2] );
88 fg.ReadFromFile( infile );
89
90 // Output various statistics
91 cout << "Number of variables: " << fg.nrVars() << endl;
92 cout << "Number of factors: " << fg.nrFactors() << endl;
93 cout << "Connected: " << fg.isConnected() << endl;
94 cout << "Tree: " << fg.isTree() << endl;
95 cout << "Has short loops: " << hasShortLoops(fg.factors()) << endl;
96 cout << "Has negatives: " << hasNegatives(fg.factors()) << endl;
97 cout << "Binary variables? " << fg.isBinary() << endl;
98 cout << "Pairwise interactions? " << fg.isPairwise() << endl;
99
100 // Calculate treewidth using various heuristics
101 std::pair<size_t,double> tw;
102 cout << "Treewidth (MinNeighbors): ";
103 try {
104 tw = boundTreewidth(fg, &eliminationCost_MinNeighbors, maxstates );
105 cout << tw.first << " (" << tw.second << " states)" << endl;
106 } catch( Exception &e ) {
107 if( e.code() == Exception::OUT_OF_MEMORY )
108 cout << "> " << maxstates << endl;
109 else
110 cout << "an exception occurred" << endl;
111 }
112
113 cout << "Treewidth (MinWeight): ";
114 try {
115 tw = boundTreewidth(fg, &eliminationCost_MinWeight, maxstates );
116 cout << tw.first << " (" << tw.second << " states)" << endl;
117 } catch( Exception &e ) {
118 if( e.code() == Exception::OUT_OF_MEMORY )
119 cout << "> " << maxstates << endl;
120 else
121 cout << "an exception occurred" << endl;
122 }
123
124 cout << "Treewidth (MinFill): ";
125 try {
126 tw = boundTreewidth(fg, &eliminationCost_MinFill, maxstates );
127 cout << tw.first << " (" << tw.second << " states)" << endl;
128 } catch( Exception &e ) {
129 if( e.code() == Exception::OUT_OF_MEMORY )
130 cout << "> " << maxstates << endl;
131 else
132 cout << "an exception occurred" << endl;
133 }
134
135 cout << "Treewidth (WeightedMinFill): ";
136 try {
137 tw = boundTreewidth(fg, &eliminationCost_WeightedMinFill, maxstates );
138 cout << tw.first << " (" << tw.second << " states)" << endl;
139 } catch( Exception &e ) {
140 if( e.code() == Exception::OUT_OF_MEMORY )
141 cout << "> " << maxstates << endl;
142 else
143 cout << "an exception occurred" << endl;
144 }
145
146 // Calculate total state space
147 long double stsp = 1.0;
148 for( size_t i = 0; i < fg.nrVars(); i++ )
149 stsp *= fg.var(i).states();
150 cout << "Total state space: " << stsp << endl;
151 // Output type of factor graph
152 cout << "Type: " << (fg.isPairwise() ? "pairwise" : "higher order") << " interactions, " << (fg.isBinary() ? "binary" : "nonbinary") << " variables" << endl;
153
154 // Calculate complexity for LCBP
155 long double cavsum_lcbp = 0.0;
156 long double cavsum_lcbp2 = 0.0;
157 size_t max_Delta_size = 0;
158 map<size_t,size_t> cavsizes;
159 for( size_t i = 0; i < fg.nrVars(); i++ ) {
160 VarSet di = fg.delta(i);
161 if( cavsizes.count(di.size()) )
162 cavsizes[di.size()]++;
163 else
164 cavsizes[di.size()] = 1;
165 size_t Ds = fg.Delta(i).nrStates();
166 if( Ds > max_Delta_size )
167 max_Delta_size = Ds;
168 cavsum_lcbp += di.nrStates();
169 for( VarSet::const_iterator j = di.begin(); j != di.end(); j++ )
170 cavsum_lcbp2 += j->states();
171 }
172 cout << "Maximum pancake has " << max_Delta_size << " states" << endl;
173 cout << "LCBP with full cavities needs " << cavsum_lcbp << " BP runs" << endl;
174 cout << "LCBP with only pairinteractions needs " << cavsum_lcbp2 << " BP runs" << endl;
175 cout << "Cavity sizes: ";
176 for( map<size_t,size_t>::const_iterator it = cavsizes.begin(); it != cavsizes.end(); it++ )
177 cout << it->first << "(" << it->second << ") ";
178 cout << endl;
179
180 // Calculate girth and length of loops
181 if( fg.isPairwise() ) {
182 bool girth_reached = false;
183 size_t loopdepth;
184 for( loopdepth = 2; loopdepth <= fg.nrVars() && !girth_reached; loopdepth++ ) {
185 size_t nr_loops = countLoops( fg, loopdepth );
186 cout << "Loops up to " << loopdepth << " variables: " << nr_loops << endl;
187 if( nr_loops > 0 )
188 girth_reached = true;
189 }
190 if( girth_reached )
191 cout << "Girth: " << loopdepth-1 << endl;
192 else
193 cout << "Girth: infinity" << endl;
194 }
195
196 // Output factor state spaces
197 map<size_t,size_t> facsizes;
198 for( size_t I = 0; I < fg.nrFactors(); I++ ) {
199 size_t Isize = fg.factor(I).vars().size();
200 if( facsizes.count( Isize ) )
201 facsizes[Isize]++;
202 else
203 facsizes[Isize] = 1;
204 }
205 cout << "Factor sizes: ";
206 for( map<size_t,size_t>::const_iterator it = facsizes.begin(); it != facsizes.end(); it++ )
207 cout << it->first << "(" << it->second << ") ";
208 cout << endl;
209
210 return 0;
211 }
212 }