Merged tests/*, matlab/*, utils/* from SVN head...
[libdai.git] / utils / fginfo.cpp
1 /* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
2 Radboud University Nijmegen, The Netherlands
3
4 This file is part of libDAI.
5
6 libDAI is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2 of the License, or
9 (at your option) any later version.
10
11 libDAI is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with libDAI; if not, write to the Free Software
18 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
19 */
20
21
22 #include <iostream>
23 #include <cstdlib>
24 #include <dai/factorgraph.h>
25 #include <dai/jtree.h>
26
27
28 using namespace std;
29 using namespace dai;
30
31
32 void findLoopClusters( const FactorGraph & fg, std::set<VarSet> &allcl, VarSet newcl, const Var & root, size_t length, VarSet vars ) {
33 for( VarSet::const_iterator in = vars.begin(); in != vars.end(); in++ ) {
34 VarSet ind = fg.delta( *in );
35 if( (newcl.size()) >= 2 && ind.contains( root ) ) {
36 allcl.insert( newcl | *in );
37 }
38 else if( length > 1 )
39 findLoopClusters( fg, allcl, newcl | *in, root, length - 1, ind / newcl );
40 }
41 }
42
43
44 size_t countLoops( const FactorGraph & fg, size_t loopdepth ) {
45 set<VarSet> loops;
46 for( vector<Var>::const_iterator i0 = fg.vars.begin(); i0 != fg.vars.end(); i0++ ) {
47 VarSet i0d = fg.delta(*i0);
48 if( loopdepth > 1 )
49 findLoopClusters( fg, loops, *i0, *i0, loopdepth - 1, i0d );
50 }
51 return loops.size();
52 }
53
54
55 bool hasShortLoops( const std::vector<Factor> &P ) {
56 bool found = false;
57 vector<Factor>::const_iterator I, J;
58 for( I = P.begin(); I != P.end(); I++ ) {
59 J = I;
60 J++;
61 for( ; J != P.end(); J++ )
62 if( (I->vars() & J->vars()).size() >= 2 ) {
63 found = true;
64 break;
65 }
66 if( found )
67 break;
68 }
69 return found;
70 }
71
72
73 bool hasNegatives( const std::vector<Factor> &P ) {
74 bool found = false;
75 for( size_t I = 0; I < P.size(); I++ )
76 if( P[I].hasNegatives() ) {
77 found = true;
78 break;
79 }
80 return found;
81 }
82
83
84 int main( int argc, char *argv[] ) {
85 if( argc != 3 ) {
86 cout << "Usage: " << argv[0] << " <in.fg> <tw>" << endl << endl;
87 cout << "Reports some characteristics of the .fg network." << endl;
88 cout << "Also calculates treewidth (which may take some time) unless <tw> == 0." << endl;
89 return 1;
90 } else {
91 // Read factorgraph
92 FactorGraph fg;
93 char *infile = argv[1];
94 int calc_tw = atoi(argv[2]);
95 fg.ReadFromFile( infile );
96
97 cout << "Number of variables: " << fg.nrVars() << endl;
98 cout << "Number of factors: " << fg.nrFactors() << endl;
99 cout << "Connected: " << fg.isConnected() << endl;
100 cout << "Tree: " << fg.isTree() << endl;
101 cout << "Has short loops: " << hasShortLoops(fg.factors()) << endl;
102 cout << "Has negatives: " << hasNegatives(fg.factors()) << endl;
103 cout << "Binary variables? " << fg.isBinary() << endl;
104 cout << "Pairwise interactions? " << fg.isPairwise() << endl;
105 if( calc_tw ) {
106 std::pair<size_t,size_t> tw = treewidth(fg);
107 cout << "Treewidth: " << tw.first << endl;
108 cout << "Largest cluster for JTree has " << tw.second << " states " << endl;
109 }
110 double stsp = 1.0;
111 for( size_t i = 0; i < fg.nrVars(); i++ )
112 stsp *= fg.var(i).states();
113 cout << "Total state space: " << stsp << endl;
114
115 double cavsum_lcbp = 0.0;
116 double cavsum_lcbp2 = 0.0;
117 size_t max_Delta_size = 0;
118 map<size_t,size_t> cavsizes;
119 for( size_t i = 0; i < fg.nrVars(); i++ ) {
120 VarSet di = fg.delta(i);
121 if( cavsizes.count(di.size()) )
122 cavsizes[di.size()]++;
123 else
124 cavsizes[di.size()] = 1;
125 size_t Ds = fg.Delta(i).states();
126 if( Ds > max_Delta_size )
127 max_Delta_size = Ds;
128 cavsum_lcbp += di.states();
129 for( VarSet::const_iterator j = di.begin(); j != di.end(); j++ )
130 cavsum_lcbp2 += j->states();
131 }
132 cout << "Maximum pancake has " << max_Delta_size << " states" << endl;
133 cout << "LCBP with full cavities needs " << cavsum_lcbp << " BP runs" << endl;
134 cout << "LCBP with only pairinteractions needs " << cavsum_lcbp2 << " BP runs" << endl;
135 cout << "Cavity sizes: ";
136 for( map<size_t,size_t>::const_iterator it = cavsizes.begin(); it != cavsizes.end(); it++ )
137 cout << it->first << "(" << it->second << ") ";
138 cout << endl;
139
140 cout << "Type: " << (fg.isPairwise() ? "pairwise" : "higher order") << " interactions, " << (fg.isBinary() ? "binary" : "nonbinary") << " variables" << endl;
141
142 if( fg.isPairwise() ) {
143 bool girth_reached = false;
144 size_t loopdepth;
145 for( loopdepth = 2; loopdepth <= fg.nrVars() && !girth_reached; loopdepth++ ) {
146 size_t nr_loops = countLoops( fg, loopdepth );
147 cout << "Loops up to " << loopdepth << " variables: " << nr_loops << endl;
148 if( nr_loops > 0 )
149 girth_reached = true;
150 }
151 if( girth_reached )
152 cout << "Girth: " << loopdepth-1 << endl;
153 else
154 cout << "Girth: infinity" << endl;
155 }
156
157 return 0;
158 }
159 }