Now compiles also with Visual Studio 2008 under Windows (still buggy!)
[libdai.git] / utils / createfg.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 <fstream>
24 #include <vector>
25 #include <iterator>
26 #include <algorithm>
27 #include <boost/program_options.hpp>
28 #include <boost/numeric/ublas/matrix_sparse.hpp>
29 #include <boost/numeric/ublas/io.hpp>
30 #include <dai/factorgraph.h>
31 #include <dai/weightedgraph.h>
32 #include <dai/util.h>
33
34
35 using namespace std;
36 using namespace dai;
37 namespace po = boost::program_options;
38 typedef boost::numeric::ublas::compressed_matrix<double> matrix;
39 typedef matrix::value_array_type::const_iterator matrix_vcit;
40 typedef matrix::index_array_type::const_iterator matrix_icit;
41
42
43 Factor BinaryFactor( const Var &n, double field ) {
44 assert( n.states() == 2 );
45 double buf[2];
46 buf[0] = exp(-field);
47 buf[1] = exp(field);
48 return Factor(n, &buf[0]);
49 }
50
51
52 Factor BinaryFactor( const Var &n1, const Var &n2, double coupling ) {
53 assert( n1.states() == 2 );
54 assert( n2.states() == 2 );
55 assert( n1 != n2 );
56 double buf[4];
57 buf[0] = (buf[3] = exp(coupling));
58 buf[1] = (buf[2] = exp(-coupling));
59 return Factor(n1 | n2, &buf[0]);
60 }
61
62
63 Factor RandomFactor( const VarSet &ns, double beta ) {
64 Factor fac( ns );
65 for( size_t t = 0; t < fac.states(); t++ )
66 fac[t] = exp(rnd_stdnormal() * beta);
67 return fac;
68 }
69
70
71 Factor PottsFactor( const Var &n1, const Var &n2, double beta ) {
72 Factor fac( n1 | n2, 1.0 );
73 assert( n1.states() == n2.states() );
74 for( size_t s = 0; s < n1.states(); s++ )
75 fac[ s * (n1.states() + 1) ] = exp(beta);
76 return fac;
77 }
78
79
80 void MakeHOIFG( size_t N, size_t M, size_t k, double sigma, FactorGraph &fg ) {
81 vector<Var> vars;
82 vector<Factor> factors;
83
84 vars.reserve(N);
85 for( size_t i = 0; i < N; i++ )
86 vars.push_back(Var(i,2));
87
88 for( size_t I = 0; I < M; I++ ) {
89 VarSet vars;
90 while( vars.size() < k ) {
91 do {
92 size_t newind = (size_t)(N * rnd_uniform());
93 Var newvar = Var(newind, 2);
94 if( !vars.contains( newvar ) ) {
95 vars |= newvar;
96 break;
97 }
98 } while( 1 );
99 }
100 factors.push_back( RandomFactor( vars, sigma ) );
101 }
102
103 fg = FactorGraph( factors.begin(), factors.end(), vars.begin(), vars.end(), factors.size(), vars.size() );
104 }
105
106
107 // w should be upper triangular or lower triangular
108 void WTh2FG( const matrix &w, const vector<double> &th, FactorGraph &fg ) {
109 vector<Var> vars;
110 vector<Factor> factors;
111
112 size_t N = th.size();
113 assert( (w.size1() == N) && (w.size2() == N) );
114
115 vars.reserve(N);
116 for( size_t i = 0; i < N; i++ )
117 vars.push_back(Var(i,2));
118
119 factors.reserve( w.nnz() + N );
120 // walk through the sparse array structure
121 // this is similar to matlab sparse arrays
122 // index2 gives the column index (similar to ir in matlab)
123 // index1 gives the starting indices for each row (similar to jc in matlab)
124 size_t i = 0;
125 for( size_t pos = 0; pos < w.nnz(); pos++ ) {
126 while( pos == w.index1_data()[i+1] )
127 i++;
128 size_t j = w.index2_data()[pos];
129 double w_ij = w.value_data()[pos];
130 factors.push_back( BinaryFactor( vars[i], vars[j], w_ij ) );
131 }
132 for( size_t i = 0; i < N; i++ )
133 factors.push_back( BinaryFactor( vars[i], th[i] ) );
134
135 fg = FactorGraph( factors.begin(), factors.end(), vars.begin(), vars.end(), factors.size(), vars.size() );
136 }
137
138
139 void MakeFullFG( size_t N, double mean_w, double mean_th, double sigma_w, double sigma_th, FactorGraph &fg ) {
140 matrix w(N,N,N*(N-1)/2);
141 vector<double> th(N,0.0);
142
143 for( size_t i = 0; i < N; i++ ) {
144 for( size_t j = i+1; j < N; j++ )
145 w(i,j) = rnd_stdnormal() * sigma_w + mean_w;
146 th[i] = rnd_stdnormal() * sigma_th + mean_th;
147 }
148
149 WTh2FG( w, th, fg );
150 }
151
152
153 void Make3DPotts( size_t n1, size_t n2, size_t n3, size_t states, double beta, FactorGraph &fg ) {
154 vector<Var> vars;
155 vector<Factor> factors;
156
157 for( size_t i1 = 0; i1 < n1; i1++ )
158 for( size_t i2 = 0; i2 < n2; i2++ )
159 for( size_t i3 = 0; i3 < n3; i3++ ) {
160 vars.push_back( Var( i1*n2*n3 + i2*n3 + i3, states ) );
161 if( i1 )
162 factors.push_back( PottsFactor( vars.back(), vars[ (i1-1)*n2*n3 + i2*n3 + i3 ], beta ) );
163 if( i2 )
164 factors.push_back( PottsFactor( vars.back(), vars[ i1*n2*n3 + (i2-1)*n3 + i3 ], beta ) );
165 if( i3 )
166 factors.push_back( PottsFactor( vars.back(), vars[ i1*n2*n3 + i2*n3 + (i3-1) ], beta ) );
167 }
168
169 fg = FactorGraph( factors.begin(), factors.end(), vars.begin(), vars.end(), factors.size(), vars.size() );
170 }
171
172
173 void MakeGridFG( long periodic, size_t n, double mean_w, double mean_th, double sigma_w, double sigma_th, FactorGraph &fg ) {
174 size_t N = n*n;
175
176 matrix w(N,N,2*N);
177 vector<double> th(N,0.0);
178
179 for( size_t i = 0; i < n; i++ )
180 for( size_t j = 0; j < n; j++ ) {
181 if( i+1 < n || periodic )
182 w(i*n+j, ((i+1)%n)*n+j) = rnd_stdnormal() * sigma_w + mean_w;
183 if( j+1 < n || periodic )
184 w(i*n+j, i*n+((j+1)%n)) = rnd_stdnormal() * sigma_w + mean_w;
185 th[i*n+j] = rnd_stdnormal() * sigma_th + mean_th;
186 }
187
188 WTh2FG( w, th, fg );
189 }
190
191
192 void MakeGridNonbinaryFG( bool periodic, size_t n, size_t states, double beta, FactorGraph &fg ) {
193 size_t N = n*n;
194
195 vector<Var> vars;
196 vector<Factor> factors;
197
198 vars.reserve(N);
199 for( size_t i = 0; i < N; i++ )
200 vars.push_back(Var(i, states));
201
202 factors.reserve( 2 * N );
203 for( size_t i = 0; i < n; i++ ) {
204 for( size_t j = 0; j < n; j++ ) {
205 if( i+1 < n || periodic )
206 factors.push_back( RandomFactor( vars[i*n+j] | vars[((i+1)%n)*n+j], beta ) );
207 if( j+1 < n || periodic )
208 factors.push_back( RandomFactor( vars[i*n+j] | vars[i*n+((j+1)%n)], beta ) );
209 }
210 }
211
212 fg = FactorGraph( factors.begin(), factors.end(), vars.begin(), vars.end(), factors.size(), vars.size() );
213 }
214
215
216 void MakeLoopFG( size_t N, double mean_w, double mean_th, double sigma_w, double sigma_th, FactorGraph &fg ) {
217 matrix w(N,N,N);
218 vector<double> th(N,0.0);
219
220 for( size_t i = 0; i < N; i++ ) {
221 w(i, (i+1)%N) = rnd_stdnormal() * sigma_w + mean_w;
222 th[i] = rnd_stdnormal() * sigma_th + mean_th;
223 }
224
225 WTh2FG( w, th, fg );
226 }
227
228
229 void MakeLoopNonbinaryFG( size_t N, size_t states, double beta, FactorGraph &fg ) {
230 vector<Var> vars;
231 vector<Factor> factors;
232
233 vars.reserve(N);
234 for( size_t i = 0; i < N; i++ )
235 vars.push_back(Var(i, states));
236
237 factors.reserve( N );
238 for( size_t i = 0; i < N; i++ ) {
239 factors.push_back( RandomFactor( vars[i] | vars[(i+1)%N], beta ) );
240 }
241
242 fg = FactorGraph( factors.begin(), factors.end(), vars.begin(), vars.end(), factors.size(), vars.size() );
243 }
244
245
246 void MakeTreeFG( size_t N, double mean_w, double mean_th, double sigma_w, double sigma_th, FactorGraph &fg ) {
247 matrix w(N,N,N-1);
248 vector<double> th(N,0.0);
249
250 for( size_t i = 0; i < N; i++ ) {
251 th[i] = rnd_stdnormal() * sigma_th + mean_th;
252 if( i > 0 ) {
253 size_t j = rnd_int( 0, i-1 );
254 w(i,j) = rnd_stdnormal() * sigma_w + mean_w;
255 }
256 }
257
258 WTh2FG( w, th, fg );
259 }
260
261
262 void MakeDRegFG( size_t N, size_t d, double mean_w, double mean_th, double sigma_w, double sigma_th, FactorGraph &fg ) {
263 matrix w(N,N,(d*N)/2);
264 vector<double> th(N,0.0);
265
266 UEdgeVec g = RandomDRegularGraph( N, d );
267 for( size_t i = 0; i < g.size(); i++ )
268 w(g[i].n1, g[i].n2) = rnd_stdnormal() * sigma_w + mean_w;
269
270 for( size_t i = 0; i < N; i++ )
271 th[i] = rnd_stdnormal() * sigma_th + mean_th;
272
273 WTh2FG( w, th, fg );
274 }
275
276
277 // N = number of variables
278 // n = size of variable neighborhoods
279 // K = number of factors
280 // k = size of factor neighborhoods
281 // asserts: N * n == K * k
282 BipartiteGraph CreateRandomBipartiteGraph( size_t N, size_t K, size_t n, size_t k ) {
283 BipartiteGraph G;
284
285 assert( N * n == K * k );
286
287 // build lists of degree-repeated vertex numbers
288 std::vector<size_t> stubs1(N*n,0);
289 for( size_t i = 0; i < N; i++ )
290 for( size_t t = 0; t < n; t++ )
291 stubs1[i*n + t] = i;
292
293 // build lists of degree-repeated vertex numbers
294 std::vector<size_t> stubs2(K*k,0);
295 for( size_t I = 0; I < K; I++ )
296 for( size_t t = 0; t < k; t++ )
297 stubs2[I*k + t] = I;
298
299 // shuffle lists
300 random_shuffle( stubs1.begin(), stubs1.end() );
301 random_shuffle( stubs2.begin(), stubs2.end() );
302
303 // add edges
304 vector<BipartiteGraph::Edge> edges;
305 edges.reserve( N*n );
306 for( size_t e = 0; e < N*n; e++ )
307 edges.push_back( BipartiteGraph::Edge(stubs1[e], stubs2[e]) );
308
309 // finish construction
310 G.create( N, K, edges.begin(), edges.end() );
311
312 return G;
313 }
314
315
316 // Returns x**n % p, assuming p is prime
317 int powmod (int x, int n, int p) {
318 int y = 1;
319 for( int m = 0; m < n; m++ )
320 y = (x * y) % p;
321 return y;
322 }
323
324
325 // Returns order of x in GF(p) with p prime
326 size_t order (int x, int p) {
327 x = x % p;
328 assert( x != 0 );
329 size_t n = 0;
330 size_t prod = 1;
331 do {
332 prod = (prod * x) % p;
333 n++;
334 } while( prod != 1 );
335 return n;
336 }
337
338
339 // Returns whether n is a prime number
340 bool isPrime (size_t n) {
341 bool result = true;
342 for( size_t k = 2; (k < n) && result; k++ )
343 if( n % k == 0 )
344 result = false;
345 return result;
346 }
347
348
349 // Make a regular LDPC graph with N=6, j=2, K=4, k=3
350 BipartiteGraph CreateSmallLDPCGraph() {
351 BipartiteGraph G;
352 size_t N=4, j=3, K=4; // k=3;
353
354 typedef BipartiteGraph::Edge Edge;
355 vector<Edge> edges;
356 edges.reserve( N*j );
357 edges.push_back( Edge(0,0) ); edges.push_back( Edge(1,0) ); edges.push_back( Edge(2,0) );
358 edges.push_back( Edge(0,1) ); edges.push_back( Edge(1,1) ); edges.push_back( Edge(3,1) );
359 edges.push_back( Edge(0,2) ); edges.push_back( Edge(2,2) ); edges.push_back( Edge(3,2) );
360 edges.push_back( Edge(1,3) ); edges.push_back( Edge(2,3) ); edges.push_back( Edge(3,3) );
361
362 // finish construction
363 G.create( N, K, edges.begin(), edges.end() );
364
365 return G;
366 }
367
368
369 // Use construction described in "A Class of Group-Structured LDPC Codes"
370 // by R. M. Tanner, D. Sridhara and T. Fuja
371 // Proceedings of ICSTA, 2001
372 //
373 // Example parameters: (p,j,k) = (31,3,5)
374 // j and k must be divisors of p-1
375 BipartiteGraph CreateGroupStructuredLDPCGraph( size_t p, size_t j, size_t k ) {
376 BipartiteGraph G;
377
378 size_t n = j;
379 size_t N = p * k;
380 size_t K = p * j;
381
382 size_t a, b;
383 for( a = 2; a < p; a++ )
384 if( order(a,p) == k )
385 break;
386 assert( a != p );
387 for( b = 2; b < p; b++ )
388 if( order(b,p) == j )
389 break;
390 assert( b != p );
391 // cout << "# order(a=" << a << ") = " << order(a,p) << endl;
392 // cout << "# order(b=" << b << ") = " << order(b,p) << endl;
393
394 assert( N * n == K * k );
395
396 typedef BipartiteGraph::Edge Edge;
397 vector<Edge> edges;
398 edges.reserve( N * n );
399
400 for( size_t s = 0; s < j; s++ )
401 for( size_t t = 0; t < k; t++ ) {
402 size_t P = (powmod(b,s,p) * powmod(a,t,p)) % p;
403 for( size_t m = 0; m < p; m++ )
404 edges.push_back( Edge(t*p + m, s*p + ((m + P) % p)) );
405 }
406
407 // finish construction
408 G.create( N, K, edges.begin(), edges.end() );
409
410 return G;
411 }
412
413
414 // Make parity check table
415 void MakeParityCheck( double *result, size_t n, double eps ) {
416 size_t N = 1 << n;
417 for( size_t i = 0; i < N; i++ ) {
418 size_t c = 0;
419 for( size_t t = 0; t < n; t++ )
420 if( i & (1 << t) )
421 c ^= 1;
422 if( c )
423 result[i] = eps;
424 else
425 result[i] = 1.0 - eps;
426 }
427 return;
428 }
429
430
431 const char *FULL_TYPE = "full";
432 const char *GRID_TYPE = "grid";
433 const char *GRID_TORUS_TYPE = "grid_torus";
434 const char *DREG_TYPE = "dreg";
435 const char *HOI_TYPE = "hoi";
436 const char *LDPC_RANDOM_TYPE = "ldpc_random";
437 const char *LDPC_GROUP_TYPE = "ldpc_group";
438 const char *LDPC_SMALL_TYPE = "ldpc_small";
439 const char *LOOP_TYPE = "loop";
440 const char *TREE_TYPE = "tree";
441 const char *POTTS3D_TYPE = "potts3d";
442
443
444 int main( int argc, char *argv[] ) {
445 try {
446 size_t N, K, k, d, j, n1, n2, n3;
447 size_t prime;
448 size_t seed;
449 double beta, sigma_w, sigma_th, noise, mean_w, mean_th;
450 string type;
451 size_t states = 2;
452
453 // Declare the supported options.
454 po::options_description desc("Allowed options");
455 desc.add_options()
456 ("help", "produce help message")
457 ("type", po::value<string>(&type), "factor graph type:\n\t'full', 'grid', 'grid_torus', 'dreg', 'loop', 'tree', 'hoi', 'ldpc_random', 'ldpc_group', 'ldpc_small', 'potts3d'")
458 ("seed", po::value<size_t>(&seed), "random number seed (tries to read from /dev/urandom if not specified)")
459 ("N", po::value<size_t>(&N), "number of variables (not for type=='ldpc_small')")
460 ("n1", po::value<size_t>(&n1), "width of 3D grid (only for type=='potts3d')")
461 ("n2", po::value<size_t>(&n2), "height of 3D grid (only for type=='potts3d')")
462 ("n3", po::value<size_t>(&n3), "length of 3D grid (only for type=='potts3d')")
463 ("K", po::value<size_t>(&K), "number of factors\n\t(only for type=='hoi' and 'type=='ldpc_{random,group}')")
464 ("k", po::value<size_t>(&k), "number of variables per factor\n\t(only for type=='hoi' and type=='ldpc_{random,group}')")
465 ("d", po::value<size_t>(&d), "variable connectivity\n\t(only for type=='dreg')")
466 ("j", po::value<size_t>(&j), "number of parity checks per bit\n\t(only for type=='ldpc_{random,group}')")
467 ("prime", po::value<size_t>(&prime), "prime number for construction of LDPC code\n\t(only for type=='ldpc_group')")
468 ("beta", po::value<double>(&beta), "stddev of log-factor entries\n\t(only for type=='hoi', 'potts3d', 'grid' if states>2)")
469 ("mean_w", po::value<double>(&mean_w), "mean of pairwise interactions w_{ij}\n\t(not for type=='hoi', 'ldpc_*', 'potts3d')")
470 ("mean_th", po::value<double>(&mean_th), "mean of singleton interactions th_i\n\t(not for type=='hoi', 'ldpc_*', 'potts3d')")
471 ("sigma_w", po::value<double>(&sigma_w), "stddev of pairwise interactions w_{ij}\n\t(not for type=='hoi', 'ldpc_*', 'potts3d')")
472 ("sigma_th", po::value<double>(&sigma_th), "stddev of singleton interactions th_i\n\t(not for type=='hoi', 'ldpc_*', 'potts3d'")
473 ("noise", po::value<double>(&noise), "bitflip probability for binary symmetric channel (only for type=='ldpc')")
474 ("states", po::value<size_t>(&states), "number of states of each variable (should be 2 for all but type=='grid', 'grid_torus', 'loop', 'potts3d')")
475 ;
476
477 po::variables_map vm;
478 po::store(po::parse_command_line(argc, argv, desc), vm);
479 po::notify(vm);
480
481 if( vm.count("help") || !vm.count("type") ) {
482 if( vm.count("type") ) {
483 if( type == FULL_TYPE ) {
484 cout << "Creates fully connected pairwise graphical model of <N> binary variables;" << endl;
485 } else if( type == GRID_TYPE ) {
486 cout << "Creates (non-periodic) 2D Ising grid of (approx.) <N> variables (which need not be binary);" << endl;
487 } else if( type == GRID_TORUS_TYPE ) {
488 cout << "Creates periodic 2D Ising grid of (approx.) <N> variables (which need not be binary);" << endl;
489 } else if( type == DREG_TYPE ) {
490 cout << "Creates random d-regular graph of <N> binary variables with uniform degree <d>" << endl;
491 cout << "(where <d><N> should be even);" << endl;
492 } else if( type == LOOP_TYPE ) {
493 cout << "Creates a pairwise graphical model consisting of a single loop of" << endl;
494 cout << "<N> variables (which need not be binary);" << endl;
495 } else if( type == TREE_TYPE ) {
496 cout << "Creates a pairwise, connected graphical model without cycles (i.e., a tree)" << endl;
497 cout << "of <N> binary variables;" << endl;
498 } else if( type == HOI_TYPE ) {
499 cout << "Creates a random factor graph of <N> binary variables and" << endl;
500 cout << "<K> factors, each factor being an interaction of <k> variables." << endl;
501 cout << "The entries of the factors are exponentials of i.i.d. Gaussian" << endl;
502 cout << "variables with mean 0 and standard deviation <beta>." << endl;
503 } else if( type == LDPC_RANDOM_TYPE ) {
504 cout << "Simulates LDPC decoding problem, using a LDPC code of <N> bits and <K> parity" << endl;
505 cout << "checks, with <k> bits per check and <j> checks per bit, transmitted on a binary" << endl;
506 cout << "symmetric channel with probability <noise> of flipping a bit. The transmitted" << endl;
507 cout << "codeword has all bits set to zero. The LDPC code is randomly generated." << endl;
508 } else if( type == LDPC_GROUP_TYPE ) {
509 cout << "Simulates LDPC decoding problem, using a LDPC code of <N> bits and <K> parity" << endl;
510 cout << "checks, with <k> bits per check and <j> checks per bit, transmitted on a binary" << endl;
511 cout << "symmetric channel with probability <noise> of flipping a bit. The transmitted" << endl;
512 cout << "codeword has all bits set to zero. The LDPC code is constructed (using group" << endl;
513 cout << "theory) using a parameter <prime>; <j> and <k> should both be divisors of <prime>-1." << endl;
514 } else if( type == LDPC_SMALL_TYPE ) {
515 cout << "Simulates LDPC decoding problem, using a LDPC code of 4 bits and 4 parity" << endl;
516 cout << "checks, with 3 bits per check and 3 checks per bit, transmitted on a binary" << endl;
517 cout << "symmetric channel with probability <noise> of flipping a bit. The transmitted" << endl;
518 cout << "codeword has all bits set to zero. The LDPC code is fixed." << endl;
519 } else if( type == POTTS3D_TYPE ) {
520 cout << "Builds 3D Potts model of size <n1>x<n2>x<n3> with nearest-neighbour Potts" << endl;
521 cout << "interactions with <states> states and inverse temperature <beta>." << endl;
522 } else
523 cerr << "Unknown type (should be one of 'full', 'grid', 'grid_torus', 'dreg', 'loop', 'tree', 'hoi', 'ldpc_random', 'ldpc_group', 'ldpc_small', 'potts3d')" << endl;
524
525 if( type == FULL_TYPE || type == GRID_TYPE || type == GRID_TORUS_TYPE || type == DREG_TYPE || type == LOOP_TYPE || type == TREE_TYPE ) {
526 if( type == GRID_TYPE || type == GRID_TORUS_TYPE || type == LOOP_TYPE ) {
527 cout << "if <states> > 2: factor entries are exponents of Gaussians with mean 0 and standard deviation beta; otherwise," << endl;
528 }
529 cout << "singleton interactions are Gaussian with mean <mean_th> and standard" << endl;
530 cout << "deviation <sigma_th>; pairwise interactions are Gaussian with mean" << endl;
531 cout << "<mean_w> and standard deviation <sigma_w>." << endl;
532 }
533 }
534 cout << endl << desc << endl;
535 return 1;
536 }
537
538 if( !vm.count("states") )
539 states = 2;
540
541 if( !vm.count("seed") ) {
542 ifstream infile;
543 bool success;
544 infile.open( "/dev/urandom" );
545 success = infile.is_open();
546 if( success ) {
547 infile.read( (char *)&seed, sizeof(size_t) / sizeof(char) );
548 success = infile.good();
549 infile.close();
550 }
551 if( !success )
552 throw "Please specify random number seed.";
553 }
554 rnd_seed( seed );
555
556 FactorGraph fg;
557
558 cout << "# Factor graph made by " << argv[0] << endl;
559 cout << "# type = " << type << endl;
560
561 if( type == FULL_TYPE ) {
562 if( !vm.count("N") || !vm.count("mean_w") || !vm.count("mean_th") || !vm.count("sigma_w") || !vm.count("sigma_th") )
563 throw "Please specify all required arguments";
564 MakeFullFG( N, mean_w, mean_th, sigma_w, sigma_th, fg );
565
566 cout << "# N = " << N << endl;
567 cout << "# mean_w = " << mean_w << endl;
568 cout << "# mean_th = " << mean_th << endl;
569 cout << "# sigma_w = " << sigma_w << endl;
570 cout << "# sigma_th = " << sigma_th << endl;
571 } else if( type == GRID_TYPE || type == GRID_TORUS_TYPE ) {
572 if( states > 2 ) {
573 if( !vm.count("N") || !vm.count("beta") )
574 throw "Please specify all required arguments";
575 } else {
576 if( !vm.count("N") || !vm.count("mean_w") || !vm.count("mean_th") || !vm.count("sigma_w") || !vm.count("sigma_th") )
577 throw "Please specify all required arguments";
578 }
579
580 size_t n = (size_t)sqrt((long double)N);
581 N = n * n;
582
583 bool periodic = false;
584 if( type == GRID_TYPE )
585 periodic = false;
586 else
587 periodic = true;
588
589 if( states > 2 )
590 MakeGridNonbinaryFG( periodic, n, states, beta, fg );
591 else
592 MakeGridFG( periodic, n, mean_w, mean_th, sigma_w, sigma_th, fg );
593
594 cout << "# n = " << n << endl;
595 cout << "# N = " << N << endl;
596
597 if( states > 2 )
598 cout << "# beta = " << beta << endl;
599 else {
600 cout << "# mean_w = " << mean_w << endl;
601 cout << "# mean_th = " << mean_th << endl;
602 cout << "# sigma_w = " << sigma_w << endl;
603 cout << "# sigma_th = " << sigma_th << endl;
604 }
605 } else if( type == DREG_TYPE ) {
606 if( !vm.count("N") || !vm.count("mean_w") || !vm.count("mean_th") || !vm.count("sigma_w") || !vm.count("sigma_th") || !vm.count("d") )
607 throw "Please specify all required arguments";
608
609 MakeDRegFG( N, d, mean_w, mean_th, sigma_w, sigma_th, fg );
610
611 cout << "# N = " << N << endl;
612 cout << "# d = " << d << endl;
613 cout << "# mean_w = " << mean_w << endl;
614 cout << "# mean_th = " << mean_th << endl;
615 cout << "# sigma_w = " << sigma_w << endl;
616 cout << "# sigma_th = " << sigma_th << endl;
617 } else if( type == LOOP_TYPE ) {
618 if( states > 2 ) {
619 if( !vm.count("N") || !vm.count("beta") )
620 throw "Please specify all required arguments";
621 } else {
622 if( !vm.count("N") || !vm.count("mean_w") || !vm.count("mean_th") || !vm.count("sigma_w") || !vm.count("sigma_th") )
623 throw "Please specify all required arguments";
624 }
625 if( states > 2 )
626 MakeLoopNonbinaryFG( N, states, beta, fg );
627 else
628 MakeLoopFG( N, mean_w, mean_th, sigma_w, sigma_th, fg );
629
630 cout << "# N = " << N << endl;
631
632 if( states > 2 )
633 cout << "# beta = " << beta << endl;
634 else {
635 cout << "# mean_w = " << mean_w << endl;
636 cout << "# mean_th = " << mean_th << endl;
637 cout << "# sigma_w = " << sigma_w << endl;
638 cout << "# sigma_th = " << sigma_th << endl;
639 }
640 } else if( type == TREE_TYPE ) {
641 if( !vm.count("N") || !vm.count("mean_w") || !vm.count("mean_th") || !vm.count("sigma_w") || !vm.count("sigma_th") )
642 throw "Please specify all required arguments";
643 MakeTreeFG( N, mean_w, mean_th, sigma_w, sigma_th, fg );
644
645 cout << "# N = " << N << endl;
646 cout << "# mean_w = " << mean_w << endl;
647 cout << "# mean_th = " << mean_th << endl;
648 cout << "# sigma_w = " << sigma_w << endl;
649 cout << "# sigma_th = " << sigma_th << endl;
650 } else if( type == HOI_TYPE ) {
651 if( !vm.count("N") || !vm.count("K") || !vm.count("k") || !vm.count("beta") )
652 throw "Please specify all required arguments";
653 do {
654 MakeHOIFG( N, K, k, beta, fg );
655 } while( !fg.G.isConnected() );
656
657 cout << "# N = " << N << endl;
658 cout << "# K = " << K << endl;
659 cout << "# k = " << k << endl;
660 cout << "# beta = " << beta << endl;
661 } else if( type == LDPC_RANDOM_TYPE || type == LDPC_GROUP_TYPE || type == LDPC_SMALL_TYPE ) {
662 if( !vm.count("noise") )
663 throw "Please specify all required arguments";
664
665 if( type == LDPC_RANDOM_TYPE ) {
666 if( !vm.count("N") || !vm.count("K") || !vm.count("j") || !vm.count("k") )
667 throw "Please specify all required arguments";
668
669 if( N * j != K * k )
670 throw "Parameters should satisfy N * j == K * k";
671 } else if( type == LDPC_GROUP_TYPE ) {
672 if( !vm.count("prime") || !vm.count("j") || !vm.count("k") )
673 throw "Please specify all required arguments";
674
675 if( !isPrime(prime) )
676 throw "Parameter <prime> should be prime";
677 if( !((prime-1) % j == 0 ) )
678 throw "Parameters should satisfy (prime-1) % j == 0";
679 if( !((prime-1) % k == 0 ) )
680 throw "Parameters should satisfy (prime-1) % k == 0";
681
682 N = prime * k;
683 K = prime * j;
684 } else if( type == LDPC_SMALL_TYPE ) {
685 N = 4;
686 K = 4;
687 j = 3;
688 k = 3;
689 }
690
691 cout << "# N = " << N << endl;
692 cout << "# K = " << K << endl;
693 cout << "# j = " << j << endl;
694 cout << "# k = " << k << endl;
695 if( type == LDPC_GROUP_TYPE )
696 cout << "# prime = " << prime << endl;
697 cout << "# noise = " << noise << endl;
698
699 // p = 31, j = 3, k = 5
700 // p = 37, j = 3, k = 4
701 // p = 7 , j = 2, k = 3
702 // p = 29, j = 2, k = 4
703
704 // Construct likelihood and paritycheck factors
705 double likelihood[4] = {1.0 - noise, noise, noise, 1.0 - noise};
706 double *paritycheck = new double[1 << k];
707 MakeParityCheck(paritycheck, k, 0.0);
708
709 // Create LDPC structure
710 BipartiteGraph ldpcG;
711 bool regular;
712 do {
713 if( type == LDPC_GROUP_TYPE )
714 ldpcG = CreateGroupStructuredLDPCGraph( prime, j, k );
715 else if( type == LDPC_RANDOM_TYPE )
716 ldpcG = CreateRandomBipartiteGraph( N, K, j, k );
717 else if( type == LDPC_SMALL_TYPE )
718 ldpcG = CreateSmallLDPCGraph();
719
720 regular = true;
721 for( size_t i = 0; i < N; i++ )
722 if( ldpcG.nb1(i).size() != j )
723 regular = false;
724 for( size_t I = 0; I < K; I++ )
725 if( ldpcG.nb2(I).size() != k )
726 regular = false;
727 } while( !regular && !ldpcG.isConnected() );
728
729 // Convert to FactorGraph
730 vector<Factor> factors;
731 for( size_t I = 0; I < K; I++ ) {
732 VarSet vs;
733 for( size_t _i = 0; _i < k; _i++ ) {
734 size_t i = ldpcG.nb2(I)[_i];
735 vs |= Var( i, 2 );
736 }
737 factors.push_back( Factor( vs, paritycheck ) );
738 }
739 delete paritycheck;
740
741 // Generate noise vector
742 vector<char> noisebits(N,0);
743 size_t bitflips = 0;
744 for( size_t i = 0; i < N; i++ ) {
745 if( rnd_uniform() < noise ) {
746 noisebits[i] = 1;
747 bitflips++;
748 }
749 }
750 cout << "# bitflips = " << bitflips << endl;
751
752 // Simulate transmission of all-zero codeword
753 vector<char> input(N,0);
754 vector<char> output(N,0);
755 for( size_t i = 0; i < N; i++ )
756 output[i] = (input[i] + noisebits[i]) & 1;
757
758 // Add likelihoods
759 for( size_t i = 0; i < N; i++ )
760 factors.push_back( Factor(Var(i,2), likelihood + output[i]*2) );
761
762 // Construct Factor Graph
763 fg = FactorGraph( factors );
764 } else if( type == POTTS3D_TYPE ) {
765 if( !vm.count("n1") || !vm.count("n2") || !vm.count("n3") || !vm.count("beta") || !vm.count("states") )
766 throw "Please specify all required arguments";
767 Make3DPotts( n1, n2, n3, states, beta, fg );
768
769 cout << "# N = " << n1*n2*n3 << endl;
770 cout << "# n1 = " << n1 << endl;
771 cout << "# n2 = " << n2 << endl;
772 cout << "# n3 = " << n3 << endl;
773 cout << "# beta = " << beta << endl;
774 cout << "# states = " << states << endl;
775 } else {
776 throw "Invalid type";
777 }
778
779 cout << "# seed = " << seed << endl;
780 cout << fg;
781 }
782 catch(exception& e) {
783 cerr << "Error: " << e.what() << endl;
784 return 1;
785 }
786 catch(const char * e) {
787 cerr << "Error: " << e << endl;
788 return 1;
789 }
790 catch(...) {
791 cerr << "Exception of unknown type!" << endl;
792 }
793
794 return 0;
795 }