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