Removed deprecated interfaces
[libdai.git] / utils / uai2fg.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) 2008-2010 Joris Mooij [joris dot mooij at libdai dot org]
8 */
9
10
11 #include <iostream>
12 #include <fstream>
13 #include <dai/alldai.h>
14 #include <dai/util.h>
15 #include <dai/index.h>
16 #include <dai/jtree.h>
17
18
19 using namespace std;
20 using namespace dai;
21
22
23 /// Reads "evidence" (a mapping from observed variable labels to the observed values) from a UAI evidence file
24 map<size_t, size_t> ReadUAIEvidenceFile( char* filename ) {
25 map<size_t, size_t> evid;
26
27 // open file
28 ifstream is;
29 is.open( filename );
30 if( is.is_open() ) {
31 // read number of observed variables
32 size_t nr_evid;
33 is >> nr_evid;
34 if( is.fail() )
35 DAI_THROWE(INVALID_EVIDENCE_FILE,"Cannot read number of observed variables");
36
37 // for each observation, read the variable label and the observed value
38 for( size_t i = 0; i < nr_evid; i++ ) {
39 size_t label, val;
40 is >> label;
41 if( is.fail() )
42 DAI_THROWE(INVALID_EVIDENCE_FILE,"Cannot read label for " + toString(i) + "'th observed variable");
43 is >> val;
44 if( is.fail() )
45 DAI_THROWE(INVALID_EVIDENCE_FILE,"Cannot read value of " + toString(i) + "'th observed variable");
46 evid[label] = val;
47 }
48
49 // close file
50 is.close();
51 } else
52 DAI_THROWE(CANNOT_READ_FILE,"Cannot read from file " + std::string(filename));
53
54 return evid;
55 }
56
57
58 /// Reads factor graph (as a pair of a variable vector and factor vector) from a UAI factor graph file
59 void ReadUAIFGFile( const char *filename, size_t verbose, vector<Var>& vars, vector<Factor>& factors, vector<Permute>& permutations ) {
60 vars.clear();
61 factors.clear();
62 permutations.clear();
63
64 // open file
65 ifstream is;
66 is.open( filename );
67 if( is.is_open() ) {
68 size_t nrFacs, nrVars;
69 string line;
70
71 // read header line
72 getline(is,line);
73 if( is.fail() || (line != "BAYES" && line != "MARKOV" && line != "BAYES\r" && line != "MARKOV\r") )
74 DAI_THROWE(INVALID_FACTORGRAPH_FILE,"UAI factor graph file should start with \"BAYES\" or \"MARKOV\"");
75 if( verbose >= 2 )
76 cout << "Reading " << line << " network..." << endl;
77
78 // read number of variables
79 is >> nrVars;
80 if( is.fail() )
81 DAI_THROWE(INVALID_FACTORGRAPH_FILE,"Cannot read number of variables");
82 if( verbose >= 2 )
83 cout << "Reading " << nrVars << " variables..." << endl;
84
85 // for each variable, read its number of states
86 vars.reserve( nrVars );
87 for( size_t i = 0; i < nrVars; i++ ) {
88 size_t dim;
89 is >> dim;
90 if( is.fail() )
91 DAI_THROWE(INVALID_FACTORGRAPH_FILE,"Cannot read number of states of " + toString(i) + "'th variable");
92 vars.push_back( Var( i, dim ) );
93 }
94
95 // read number of factors
96 is >> nrFacs;
97 if( is.fail() )
98 DAI_THROWE(INVALID_FACTORGRAPH_FILE,"Cannot read number of factors");
99 if( verbose >= 2 )
100 cout << "Reading " << nrFacs << " factors..." << endl;
101
102 // for each factor, read the variables on which it depends
103 vector<vector<Var> > factorVars;
104 factors.reserve( nrFacs );
105 factorVars.reserve( nrFacs );
106 for( size_t I = 0; I < nrFacs; I++ ) {
107 if( verbose >= 3 )
108 cout << "Reading factor " << I << "..." << endl;
109
110 // read number of variables for factor I
111 size_t I_nrVars;
112 is >> I_nrVars;
113 if( is.fail() )
114 DAI_THROWE(INVALID_FACTORGRAPH_FILE,"Cannot read number of variables for " + toString(I) + "'th factor");
115 if( verbose >= 3 )
116 cout << " which depends on " << I_nrVars << " variables" << endl;
117
118 // read the variable labels
119 vector<long> I_labels;
120 vector<size_t> I_dims;
121 I_labels.reserve( I_nrVars );
122 I_dims.reserve( I_nrVars );
123 factorVars[I].reserve( I_nrVars );
124 for( size_t _i = 0; _i < I_nrVars; _i++ ) {
125 long label;
126 is >> label;
127 if( is.fail() )
128 DAI_THROWE(INVALID_FACTORGRAPH_FILE,"Cannot read variable labels for " + toString(I) + "'th factor");
129 I_labels.push_back( label );
130 I_dims.push_back( vars[label].states() );
131 factorVars[I].push_back( vars[label] );
132 }
133 if( verbose >= 3 )
134 cout << " labels: " << I_labels << ", dimensions " << I_dims << endl;
135
136 // add the factor and the labels
137 factors.push_back( Factor( VarSet( factorVars[I].begin(), factorVars[I].end(), factorVars[I].size() ), (Real)0 ) );
138 }
139
140 // for each factor, read its values
141 permutations.reserve( nrFacs );
142 for( size_t I = 0; I < nrFacs; I++ ) {
143 if( verbose >= 3 )
144 cout << "Reading factor " << I << "..." << endl;
145
146 // calculate permutation object, reversing the indexing in factorVars[I] first
147 Permute permindex( factorVars[I], true );
148 permutations.push_back( permindex );
149
150 // read factor values
151 size_t nrNonZeros;
152 is >> nrNonZeros;
153 if( is.fail() )
154 DAI_THROWE(INVALID_FACTORGRAPH_FILE,"Cannot read number of nonzero factor values for " + toString(I) + "'th factor");
155 if( verbose >= 3 )
156 cout << " number of nonzero values: " << nrNonZeros << endl;
157 DAI_ASSERT( nrNonZeros == factors[I].nrStates() );
158 for( size_t li = 0; li < nrNonZeros; li++ ) {
159 Real val;
160 is >> val;
161 if( is.fail() )
162 DAI_THROWE(INVALID_FACTORGRAPH_FILE,"Cannot read factor values of " + toString(I) + "'th factor");
163 // assign value after calculating its linear index corresponding to the permutation
164 if( verbose >= 4 )
165 cout << " " << li << "'th value " << val << " corresponds with index " << permindex.convertLinearIndex(li) << endl;
166 factors[I].set( permindex.convertLinearIndex( li ), val );
167 }
168 }
169 if( verbose >= 3 )
170 cout << "variables:" << vars << endl;
171 if( verbose >= 3 )
172 cout << "factors:" << factors << endl;
173
174 // close file
175 is.close();
176 } else
177 DAI_THROWE(CANNOT_READ_FILE,"Cannot read from file " + std::string(filename));
178 }
179
180
181 int main( int argc, char *argv[] ) {
182 if ( argc != 7 ) {
183 cout << "This program is part of libDAI - http://www.libdai.org/" << endl << endl;
184 cout << "Usage: ./uai2fg <filename.uai> <filename.uai.evid> <filename.fg> <type> <run_jtree> <verbose>" << endl << endl;
185 cout << "Converts input files in the UAI 2008 approximate inference evaluation format" << endl;
186 cout << "(see http://graphmod.ics.uci.edu/uai08/) to the libDAI factor graph format." << endl;
187 cout << "Reads factor graph <filename.uai> and evidence <filename.uai.evid>" << endl;
188 cout << "and writes the resulting clamped factor graph to <filename.fg>." << endl;
189 cout << "If type==0, uses surgery (recommended), otherwise, uses just adds delta factors." << endl;
190 cout << "If run_jtree!=0, runs a junction tree and reports the results in the UAI 2008 results file format." << endl;
191 return 1;
192 } else {
193 long verbose = atoi( argv[6] );
194 long type = atoi( argv[4] );
195 bool run_jtree = atoi( argv[5] );
196
197 // read factor graph
198 vector<Var> vars;
199 vector<Factor> facs;
200 vector<Permute> permutations;
201 ReadUAIFGFile( argv[1], verbose, vars, facs, permutations );
202
203 // read evidence
204 map<size_t,size_t> evid = ReadUAIEvidenceFile( argv[2] );
205
206 // construct unclamped factor graph
207 FactorGraph fg0( facs.begin(), facs.end(), vars.begin(), vars.end(), facs.size(), vars.size() );
208
209 // change factor graph to reflect observed evidence
210 if( type == 0 ) {
211 // replace factors with clamped variables with slices
212 for( size_t I = 0; I < facs.size(); I++ ) {
213 for( map<size_t,size_t>::const_iterator e = evid.begin(); e != evid.end(); e++ ) {
214 if( facs[I].vars() >> vars[e->first] ) {
215 if( verbose >= 2 )
216 cout << "Clamping " << e->first << " to value " << e->second << " in factor " << I << " = " << facs[I].vars() << endl;
217 facs[I] = facs[I].slice( vars[e->first], e->second );
218 if( verbose >= 2 )
219 cout << "...remaining vars: " << facs[I].vars() << endl;
220 }
221 }
222 }
223 // remove empty factors
224 double logZcorr = 0.0;
225 for( vector<Factor>::iterator I = facs.begin(); I != facs.end(); )
226 if( I->vars().size() == 0 ) {
227 logZcorr += std::log( (Real)(*I)[0] );
228 I = facs.erase( I );
229 } else
230 I++;
231 // multiply with logZcorr constant
232 if( facs.size() == 0 )
233 facs.push_back( Factor( VarSet(), std::exp(logZcorr) ) );
234 else
235 facs.front() *= std::exp(logZcorr);
236 }
237 // add delta factors corresponding to observed variable values
238 for( map<size_t,size_t>::const_iterator e = evid.begin(); e != evid.end(); e++ )
239 facs.push_back( createFactorDelta( vars[e->first], e->second ) );
240
241 // construct clamped factor graph
242 FactorGraph fg( facs.begin(), facs.end(), vars.begin(), vars.end(), facs.size(), vars.size() );
243
244 // write it to a file
245 fg.WriteToFile( argv[3] );
246
247 // if requested, perform various inference tasks
248 if( run_jtree ) {
249 // construct junction tree on unclamped factor graph
250 JTree jt0( fg0, PropertySet()("updates",string("HUGIN")) );
251 jt0.init();
252 jt0.run();
253
254 // construct junction tree on clamped factor graph
255 JTree jt( fg, PropertySet()("updates",string("HUGIN")) );
256 jt.init();
257 jt.run();
258
259 // output probability of evidence
260 cout.precision( 8 );
261 if( evid.size() )
262 cout << "z " << (jt.logZ() - jt0.logZ()) / dai::log((Real)10.0) << endl;
263 else
264 cout << "z " << jt.logZ() / dai::log((Real)10.0) << endl;
265
266 // output variable marginals
267 cout << "m " << jt.nrVars() << " ";
268 for( size_t i = 0; i < jt.nrVars(); i++ ) {
269 cout << jt.var(i).states() << " ";
270 for( size_t s = 0; s < jt.var(i).states(); s++ )
271 cout << jt.beliefV(i)[s] << " ";
272 }
273 cout << endl;
274
275 // calculate MAP state
276 jt.props.inference = JTree::Properties::InfType::MAXPROD;
277 jt.init();
278 jt.run();
279 vector<size_t> MAP = jt.findMaximum();
280 map<Var, size_t> state;
281 for( size_t i = 0; i < MAP.size(); i++ )
282 state[jt.var(i)] = MAP[i];
283 double log_MAP_prob = 0.0;
284 for( size_t I = 0; I < jt.nrFactors(); I++ )
285 log_MAP_prob += dai::log( jt.factor(I)[calcLinearState( jt.factor(I).vars(), state )] );
286
287 // output MAP state
288 cout << "s ";
289 if( evid.size() )
290 cout << (log_MAP_prob - jt0.logZ()) / dai::log((Real)10.0) << " ";
291 else
292 cout << log_MAP_prob / dai::log((Real)10.0) << " ";
293 cout << jt.nrVars() << " ";
294 for( size_t i = 0; i < jt.nrVars(); i++ )
295 cout << MAP[i] << " ";
296 cout << endl;
297 }
298 }
299
300 return 0;
301 }