Various changes:
[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-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 <fstream>
14 #include <vector>
15 #include <iterator>
16 #include <algorithm>
17 #include <boost/program_options.hpp>
18 #include <dai/factorgraph.h>
19 #include <dai/weightedgraph.h>
20 #include <dai/util.h>
21 #include <dai/graph.h>
22 #include <dai/enum.h>
23 #include <dai/properties.h>
24
25
26 using namespace std;
27 using namespace dai;
28 namespace po = boost::program_options;
29
30
31 /// Possible factor types
32 DAI_ENUM(FactorType,ISING,EXPGAUSS,POTTS);
33
34
35 /// Creates a factor graph from a pairwise interactions graph
36 /** \param G Graph describing interactions between variables
37 * \param ft Type of factors to use for interactions
38 * \param states Number of states of the variables
39 * \param props Additional properties for generating the interactions
40 */
41 FactorGraph createFG( const GraphAL &G, FactorType ft, size_t states, const PropertySet &props ) {
42 size_t N = G.nrNodes();
43
44 if( states > 2 )
45 DAI_ASSERT( ft != FactorType::ISING );
46
47 // Get inverse temperature
48 Real beta = 1.0;
49 if( ft != FactorType::ISING )
50 beta = props.GetAs<Real>("beta");
51
52 // Get properties for Ising factors
53 Real mean_h = 0.0;
54 Real sigma_h = 0.0;
55 Real mean_J = 0.0;
56 Real sigma_J = 0.0;
57 if( ft == FactorType::ISING ) {
58 mean_h = props.GetAs<Real>("mean_th");
59 sigma_h = props.GetAs<Real>("sigma_th");
60 mean_J = props.GetAs<Real>("mean_w");
61 sigma_J = props.GetAs<Real>("sigma_w");
62 }
63
64 // Create variables
65 vector<Var> vars;
66 vars.reserve( N );
67 for( size_t i = 0; i < N; i++ )
68 vars.push_back( Var( i, states ) );
69
70 // Create factors
71 vector<Factor> factors;
72 factors.reserve( G.nrEdges() + N );
73 // Pairwise factors
74 for( size_t i = 0; i < N; i++ )
75 foreach( const GraphAL::Neighbor &j, G.nb(i) )
76 if( i < j ) {
77 if( ft == FactorType::POTTS )
78 factors.push_back( createFactorPotts( vars[i], vars[j], beta ) );
79 else if( ft == FactorType::EXPGAUSS )
80 factors.push_back( createFactorExpGauss( VarSet( vars[i], vars[j] ), beta ) );
81 else if( ft == FactorType::ISING ) {
82 Real J = rnd_stdnormal() * sigma_J + mean_J;
83 factors.push_back( createFactorIsing( vars[i], vars[j], J ) );
84 }
85 }
86 // Unary factors
87 if( ft == FactorType::ISING )
88 for( size_t i = 0; i < N; i++ ) {
89 Real h = rnd_stdnormal() * sigma_h + mean_h;
90 factors.push_back( createFactorIsing( vars[i], h ) );
91 }
92
93 return FactorGraph( factors.begin(), factors.end(), vars.begin(), vars.end(), factors.size(), vars.size() );
94 }
95
96
97 /// Return a random factor graph with higher-order interactions
98 /** \param N number of variables
99 * \param M number of factors
100 * \param k number of variables that each factor depends on
101 * \param beta standard-deviation of Gaussian log-factor entries
102 */
103 FactorGraph createHOIFG( size_t N, size_t M, size_t k, Real beta ) {
104 vector<Var> vars;
105 vector<Factor> factors;
106
107 vars.reserve(N);
108 for( size_t i = 0; i < N; i++ )
109 vars.push_back(Var(i,2));
110
111 for( size_t I = 0; I < M; I++ ) {
112 VarSet vars;
113 while( vars.size() < k ) {
114 do {
115 size_t newind = (size_t)(N * rnd_uniform());
116 Var newvar = Var(newind, 2);
117 if( !vars.contains( newvar ) ) {
118 vars |= newvar;
119 break;
120 }
121 } while( 1 );
122 }
123 factors.push_back( createFactorExpGauss( vars, beta ) );
124 }
125
126 return FactorGraph( factors.begin(), factors.end(), vars.begin(), vars.end(), factors.size(), vars.size() );
127 }
128
129
130 /// Creates a regular random bipartite graph
131 /** \param N1 = number of nodes of type 1
132 * \param d1 = size of neighborhoods of nodes of type 1
133 * \param N2 = number of nodes of type 2
134 * \param d2 = size of neighborhoods of nodes of type 2
135 * \note asserts that N1 * d1 == N2 * d2
136 */
137 BipartiteGraph createRandomBipartiteGraph( size_t N1, size_t N2, size_t d1, size_t d2 ) {
138 BipartiteGraph G;
139
140 DAI_ASSERT( N1 * d1 == N2 * d2 );
141
142 // build lists of degree-repeated vertex numbers
143 std::vector<size_t> stubs1( N1*d1, 0 );
144 for( size_t n1 = 0; n1 < N1; n1++ )
145 for( size_t t = 0; t < d1; t++ )
146 stubs1[n1*d1 + t] = n1;
147
148 // build lists of degree-repeated vertex numbers
149 std::vector<size_t> stubs2( N2*d2, 0 );
150 for( size_t n2 = 0; n2 < N2; n2++ )
151 for( size_t t = 0; t < d2; t++ )
152 stubs2[n2*d2 + t] = n2;
153
154 // shuffle lists
155 random_shuffle( stubs1.begin(), stubs1.end() );
156 random_shuffle( stubs2.begin(), stubs2.end() );
157
158 // add edges
159 vector<BipartiteGraph::Edge> edges;
160 edges.reserve( N1*d1 );
161 for( size_t e = 0; e < N1*d1; e++ )
162 edges.push_back( BipartiteGraph::Edge(stubs1[e], stubs2[e]) );
163
164 // finish construction
165 G.construct( N1, N2, edges.begin(), edges.end() );
166
167 return G;
168 }
169
170
171 /// Returns x**n % p, assuming p is prime
172 int powmod (int x, int n, int p) {
173 int y = 1;
174 for( int m = 0; m < n; m++ )
175 y = (x * y) % p;
176 return y;
177 }
178
179
180 /// Returns order of x in GF(p) with p prime
181 size_t order (int x, int p) {
182 x = x % p;
183 DAI_ASSERT( x != 0 );
184 size_t n = 0;
185 size_t prod = 1;
186 do {
187 prod = (prod * x) % p;
188 n++;
189 } while( prod != 1 );
190 return n;
191 }
192
193
194 /// Returns whether n is a prime number
195 bool isPrime (size_t n) {
196 bool result = true;
197 for( size_t k = 2; (k < n) && result; k++ )
198 if( n % k == 0 )
199 result = false;
200 return result;
201 }
202
203
204 /// Constructs a regular LDPC graph with N=6, j=2, K=4, k=3
205 BipartiteGraph createSmallLDPCGraph() {
206 BipartiteGraph G;
207 size_t N=4, j=3, K=4; // k=3;
208
209 typedef BipartiteGraph::Edge Edge;
210 vector<Edge> edges;
211 edges.reserve( N*j );
212 edges.push_back( Edge(0,0) ); edges.push_back( Edge(1,0) ); edges.push_back( Edge(2,0) );
213 edges.push_back( Edge(0,1) ); edges.push_back( Edge(1,1) ); edges.push_back( Edge(3,1) );
214 edges.push_back( Edge(0,2) ); edges.push_back( Edge(2,2) ); edges.push_back( Edge(3,2) );
215 edges.push_back( Edge(1,3) ); edges.push_back( Edge(2,3) ); edges.push_back( Edge(3,3) );
216
217 // finish construction
218 G.construct( N, K, edges.begin(), edges.end() );
219
220 return G;
221 }
222
223
224 /// Creates group-structured LDPC code
225 /** Use construction described in "A Class of Group-Structured LDPC Codes"
226 * by R. M. Tanner, D. Sridhara and T. Fuja
227 * Proceedings of ICSTA, 2001
228 *
229 * Example parameters: (p,j,k) = (31,3,5)
230 * (p,j,k) = (37,3,4)
231 * (p,j,k) = (7,2,4)
232 * (p,j,k) = (29,2,4)
233 *
234 * j and k must be divisors of p-1
235 */
236 BipartiteGraph createGroupStructuredLDPCGraph( size_t p, size_t j, size_t k ) {
237 BipartiteGraph G;
238
239 size_t n = j;
240 size_t N = p * k;
241 size_t K = p * j;
242
243 size_t a, b;
244 for( a = 2; a < p; a++ )
245 if( order(a,p) == k )
246 break;
247 DAI_ASSERT( a != p );
248 for( b = 2; b < p; b++ )
249 if( order(b,p) == j )
250 break;
251 DAI_ASSERT( b != p );
252 // cout << "# order(a=" << a << ") = " << order(a,p) << endl;
253 // cout << "# order(b=" << b << ") = " << order(b,p) << endl;
254
255 DAI_ASSERT( N * n == K * k );
256
257 typedef BipartiteGraph::Edge Edge;
258 vector<Edge> edges;
259 edges.reserve( N * n );
260
261 for( size_t s = 0; s < j; s++ )
262 for( size_t t = 0; t < k; t++ ) {
263 size_t P = (powmod(b,s,p) * powmod(a,t,p)) % p;
264 for( size_t m = 0; m < p; m++ )
265 edges.push_back( Edge(t*p + m, s*p + ((m + P) % p)) );
266 }
267
268 // finish construction
269 G.construct( N, K, edges.begin(), edges.end() );
270
271 return G;
272 }
273
274
275 // Constructs a parity check table
276 void createParityCheck( Real *result, size_t n, Real eps ) {
277 size_t N = 1 << n;
278 for( size_t i = 0; i < N; i++ ) {
279 size_t c = 0;
280 for( size_t t = 0; t < n; t++ )
281 if( i & (1 << t) )
282 c ^= 1;
283 if( c )
284 result[i] = eps;
285 else
286 result[i] = 1.0 - eps;
287 }
288 return;
289 }
290
291
292 /// Predefined names of various factor graph types
293 const char *FULL_TYPE = "FULL";
294 const char *DREG_TYPE = "DREG";
295 const char *LOOP_TYPE = "LOOP";
296 const char *TREE_TYPE = "TREE";
297 const char *GRID_TYPE = "GRID";
298 const char *GRID3D_TYPE = "GRID3D";
299 const char *HOI_TYPE = "HOI";
300 const char *LDPC_TYPE = "LDPC";
301
302
303 /// Possible LDPC structures
304 DAI_ENUM(LDPCType,SMALL,GROUP,RANDOM);
305
306
307 /// Main function
308 int main( int argc, char *argv[] ) {
309 try {
310 // Variables for storing command line arguments
311 size_t seed;
312 size_t states = 2;
313 string type;
314 size_t d, N, K, k, j, n1, n2, n3, prime;
315 bool periodic = false;
316 FactorType ft;
317 LDPCType ldpc;
318 Real beta, sigma_w, sigma_th, mean_w, mean_th, noise;
319
320 // Declare the supported options.
321 po::options_description opts("General command line options");
322 opts.add_options()
323 ("help", "produce help message")
324 ("seed", po::value<size_t>(&seed), "random number seed (tries to read from /dev/urandom if not specified)")
325 ("states", po::value<size_t>(&states), "number of states of each variable (default=2 for binary variables)")
326 ;
327
328 // Graph structure options
329 po::options_description opts_graph("Options for specifying graph structure");
330 opts_graph.add_options()
331 ("type", po::value<string>(&type), "factor graph type (one of 'FULL', 'DREG', 'LOOP', 'TREE', 'GRID', 'GRID3D', 'HOI', 'LDPC')")
332 ("d", po::value<size_t>(&d), "variable connectivity (only for type=='DREG');\n\t<d><N> should be even")
333 ("N", po::value<size_t>(&N), "number of variables (not for type=='GRID','GRID3D')")
334 ("n1", po::value<size_t>(&n1), "width of grid (only for type=='GRID','GRID3D')")
335 ("n2", po::value<size_t>(&n2), "height of grid (only for type=='GRID','GRID3D')")
336 ("n3", po::value<size_t>(&n3), "length of grid (only for type=='GRID3D')")
337 ("periodic", po::value<bool>(&periodic), "periodic grid? (only for type=='GRID','GRID3D'; default=0)")
338 ("K", po::value<size_t>(&K), "number of factors (only for type=='HOI','LDPC')")
339 ("k", po::value<size_t>(&k), "number of variables per factor (only for type=='HOI','LDPC')")
340 ;
341
342 // Factor options
343 po::options_description opts_factors("Options for specifying factors");
344 opts_factors.add_options()
345 ("factors", po::value<FactorType>(&ft), "factor type (one of 'EXPGAUSS','POTTS','ISING')")
346 ("beta", po::value<Real>(&beta), "inverse temperature (ignored for factors=='ISING')")
347 ("mean_w", po::value<Real>(&mean_w), "mean of pairwise interactions w_{ij} (only for factors=='ISING')")
348 ("mean_th", po::value<Real>(&mean_th), "mean of unary interactions th_i (only for factors=='ISING')")
349 ("sigma_w", po::value<Real>(&sigma_w), "stddev of pairwise interactions w_{ij} (only for factors=='ISING')")
350 ("sigma_th", po::value<Real>(&sigma_th), "stddev of unary interactions th_i (only for factors=='ISING'")
351 ;
352
353 // LDPC options
354 po::options_description opts_ldpc("Options for specifying LDPC code factor graphs");
355 opts_ldpc.add_options()
356 ("ldpc", po::value<LDPCType>(&ldpc), "type of LDPC code (one of 'SMALL','GROUP','RANDOM')")
357 ("j", po::value<size_t>(&j), "number of parity checks per bit (only for type=='LDPC')")
358 ("noise", po::value<Real>(&noise), "bitflip probability for binary symmetric channel (only for type=='LDPC')")
359 ("prime", po::value<size_t>(&prime), "prime number for construction of LDPC code (only for type=='LDPC' with ldpc='GROUP'))")
360 ;
361
362 // All options
363 opts.add(opts_graph).add(opts_factors).add(opts_ldpc);
364
365 // Parse command line arguments
366 po::variables_map vm;
367 po::store(po::parse_command_line(argc, argv, opts), vm);
368 po::notify(vm);
369
370 // Display help message if necessary
371 if( vm.count("help") || !vm.count("type") ) {
372 cout << "This program is part of libDAI - http://www.libdai.org/" << endl << endl;
373 cout << "Usage: ./createfg [options]" << endl << endl;
374 cout << "Creates a factor graph according to the specified options." << endl << endl;
375
376 cout << endl << opts << endl;
377
378 cout << "The following factor graph types with pairwise interactions can be created:" << endl;
379 cout << "\t'FULL': fully connected graph of <N> variables" << endl;
380 cout << "\t'DREG': random regular graph of <N> variables where each variable is connected with <d> others" << endl;
381 cout << "\t'LOOP': a single loop of <N> variables" << endl;
382 cout << "\t'TREE': random tree-structured (acyclic, connected) graph of <N> variables" << endl;
383 cout << "\t'GRID': 2D grid of <n1>x<n2> variables" << endl;
384 cout << "\t'GRID3D': 3D grid of <n1>x<n2>x<n3> variables" << endl;
385 cout << "The following higher-order interactions factor graphs can be created:" << endl;
386 cout << "\t'HOI': random factor graph consisting of <N> variables and <K> factors," << endl;
387 cout << "\t each factor being an interaction of <k> variables." << endl;
388 cout << "The following LDPC code factor graphs can be created:" << endl;
389 cout << "\t'LDPC': simulates LDPC decoding problem, using an LDPC code of <N> bits and <K>" << endl;
390 cout << "\t parity checks, with <k> bits per check and <j> checks per bit, transmitted" << endl;
391 cout << "\t on a binary symmetric channel with probability <noise> of flipping a bit." << endl;
392 cout << "\t The transmitted codeword has all bits set to zero. The argument 'ldpc'" << endl;
393 cout << "\t determines how the LDPC code is constructed: either using a group structure," << endl;
394 cout << "\t or randomly, or a fixed small code with (N,K,k,j) = (4,4,3,3)." << endl << endl;
395
396 cout << "For all types except type=='LDPC', the factors have to be specified as well." << endl << endl;
397
398 cout << "EXPGAUSS factors (the default) are created by drawing all log-factor entries" << endl;
399 cout << "independently from a Gaussian with mean 0 and standard deviation <beta>." << endl << endl;
400
401 cout << "In case of pairwise interactions, one can also choose POTTS factors, for which" << endl;
402 cout << "the log-factors are simply delta functions multiplied by the strength <beta>." << endl << endl;
403
404 cout << "For pairwise interactions and binary variables, one can also use ISING factors." << endl;
405 cout << "Here variables x1...xN are assumed to be +1/-1--valued, and unary interactions" << endl;
406 cout << "are of the form exp(th*xi) with th drawn from a Gaussian distribution with mean" << endl;
407 cout << "<mean_th> and standard deviation <sigma_th>, and pairwise interactions are of the" << endl;
408 cout << "form exp(w*xi*xj) with w drawn from a Gaussian distribution with mean <mean_w>" << endl;
409 cout << "and standard deviation <sigma_w>." << endl;
410 return 1;
411 }
412
413 // Set default number of states
414 if( !vm.count("states") )
415 states = 2;
416
417 // Set default factor type
418 if( !vm.count("factors") )
419 ft = FactorType::EXPGAUSS;
420 // Check validness of factor type
421 if( ft == FactorType::POTTS )
422 if( type == HOI_TYPE )
423 throw "For factors=='POTTS', interactions should be pairwise (type!='HOI')";
424 if( ft == FactorType::ISING )
425 if( ((states != 2) || (type == HOI_TYPE)) )
426 throw "For factors=='ISING', variables should be binary (states==2) and interactions should be pairwise (type!='HOI')";
427
428 // Read random seed
429 if( !vm.count("seed") ) {
430 ifstream infile;
431 bool success;
432 infile.open( "/dev/urandom" );
433 success = infile.is_open();
434 if( success ) {
435 infile.read( (char *)&seed, sizeof(size_t) / sizeof(char) );
436 success = infile.good();
437 infile.close();
438 }
439 if( !success )
440 throw "Please specify random number seed.";
441 }
442 rnd_seed( seed );
443
444 // Set default periodicity
445 if( !vm.count("periodic") )
446 periodic = false;
447
448 // Store some options in a PropertySet object
449 PropertySet options;
450 if( vm.count("mean_th") )
451 options.Set("mean_th", mean_th);
452 if( vm.count("sigma_th") )
453 options.Set("sigma_th", sigma_th);
454 if( vm.count("mean_w") )
455 options.Set("mean_w", mean_w);
456 if( vm.count("sigma_w") )
457 options.Set("sigma_w", sigma_w);
458 if( vm.count("beta") )
459 options.Set("beta", beta);
460
461 // Output some comments
462 cout << "# Factor graph made by " << argv[0] << endl;
463 cout << "# type = " << type << endl;
464 cout << "# states = " << states << endl;
465
466 // The factor graph to be constructed
467 FactorGraph fg;
468
469 #define NEED_ARG(name, desc) do { if(!vm.count(name)) throw "Please specify " desc " with --" name; } while(0);
470
471 if( type == FULL_TYPE || type == DREG_TYPE || type == LOOP_TYPE || type == TREE_TYPE || type == GRID_TYPE || type == GRID3D_TYPE ) {
472 // Pairwise interactions
473
474 // Check command line options
475 if( type == GRID_TYPE ) {
476 NEED_ARG("n1", "width of grid");
477 NEED_ARG("n2", "height of grid");
478 N = n1 * n2;
479 } else if( type == GRID3D_TYPE ) {
480 NEED_ARG("n1", "width of grid");
481 NEED_ARG("n2", "height of grid");
482 NEED_ARG("n3", "depth of grid");
483 N = n1 * n2 * n3;
484 } else
485 NEED_ARG("N", "number of variables");
486
487 if( states > 2 || ft == FactorType::POTTS ) {
488 NEED_ARG("beta", "stddev of log-factor entries");
489 } else {
490 NEED_ARG("mean_w", "mean of pairwise interactions");
491 NEED_ARG("mean_th", "mean of unary interactions");
492 NEED_ARG("sigma_w", "stddev of pairwise interactions");
493 NEED_ARG("sigma_th", "stddev of unary interactions");
494 }
495
496 if( type == DREG_TYPE )
497 NEED_ARG("d", "connectivity (number of neighboring variables of each variable)");
498
499 // Build pairwise interaction graph
500 GraphAL G;
501 if( type == FULL_TYPE )
502 G = createGraphFull( N );
503 else if( type == DREG_TYPE )
504 G = createGraphRegular( N, d );
505 else if( type == LOOP_TYPE )
506 G = createGraphLoop( N );
507 else if( type == TREE_TYPE )
508 G = createGraphTree( N );
509 else if( type == GRID_TYPE )
510 G = createGraphGrid( n1, n2, periodic );
511 else if( type == GRID3D_TYPE )
512 G = createGraphGrid3D( n1, n2, n3, periodic );
513
514 // Construct factor graph from pairwise interaction graph
515 fg = createFG( G, ft, states, options );
516
517 // Output some additional comments
518 if( type == GRID_TYPE || type == GRID3D_TYPE ) {
519 cout << "# n1 = " << n1 << endl;
520 cout << "# n2 = " << n2 << endl;
521 if( type == GRID3D_TYPE )
522 cout << "# n3 = " << n3 << endl;
523 }
524 if( type == DREG_TYPE )
525 cout << "# d = " << d << endl;
526 cout << "# options = " << options << endl;
527 } else if( type == HOI_TYPE ) {
528 // Higher order interactions
529
530 // Check command line arguments
531 NEED_ARG("N", "number of variables");
532 NEED_ARG("K", "number of factors");
533 NEED_ARG("k", "number of variables per factor");
534 NEED_ARG("beta", "stddev of log-factor entries");
535
536 // Create higher-order interactions factor graph
537 do {
538 fg = createHOIFG( N, K, k, beta );
539 } while( !fg.isConnected() );
540
541 // Output some additional comments
542 cout << "# K = " << K << endl;
543 cout << "# k = " << k << endl;
544 cout << "# beta = " << beta << endl;
545 } else if( type == LDPC_TYPE ) {
546 // LDPC codes
547
548 // Check command line arguments
549 NEED_ARG("ldpc", "type of LDPC code");
550 NEED_ARG("noise", "bitflip probability for binary symmetric channel");
551
552 // Check more command line arguments (seperately for each LDPC type)
553 if( ldpc == LDPCType::RANDOM ) {
554 NEED_ARG("N", "number of variables");
555 NEED_ARG("K", "number of factors");
556 NEED_ARG("k", "number of variables per factor");
557 NEED_ARG("j", "number of parity checks per bit");
558 if( N * j != K * k )
559 throw "Parameters should satisfy N * j == K * k";
560 } else if( ldpc == LDPCType::GROUP ) {
561 NEED_ARG("prime", "prime number");
562 NEED_ARG("k", "number of variables per factor");
563 NEED_ARG("j", "number of parity checks per bit");
564
565 if( !isPrime(prime) )
566 throw "Parameter <prime> should be prime";
567 if( !((prime-1) % j == 0 ) )
568 throw "Parameters should satisfy (prime-1) % j == 0";
569 if( !((prime-1) % k == 0 ) )
570 throw "Parameters should satisfy (prime-1) % k == 0";
571
572 N = prime * k;
573 K = prime * j;
574 } else if( ldpc == LDPCType::SMALL ) {
575 N = 4;
576 K = 4;
577 j = 3;
578 k = 3;
579 }
580
581 // Output some additional comments
582 cout << "# N = " << N << endl;
583 cout << "# K = " << K << endl;
584 cout << "# j = " << j << endl;
585 cout << "# k = " << k << endl;
586 if( ldpc == LDPCType::GROUP )
587 cout << "# prime = " << prime << endl;
588 cout << "# noise = " << noise << endl;
589
590 // Construct likelihood and paritycheck factors
591 Real likelihood[4] = {1.0 - noise, noise, noise, 1.0 - noise};
592 Real *paritycheck = new Real[1 << k];
593 createParityCheck(paritycheck, k, 0.0);
594
595 // Create LDPC structure
596 BipartiteGraph ldpcG;
597 bool regular;
598 do {
599 if( ldpc == LDPCType::GROUP )
600 ldpcG = createGroupStructuredLDPCGraph( prime, j, k );
601 else if( ldpc == LDPCType::RANDOM )
602 ldpcG = createRandomBipartiteGraph( N, K, j, k );
603 else if( ldpc == LDPCType::SMALL )
604 ldpcG = createSmallLDPCGraph();
605
606 regular = true;
607 for( size_t i = 0; i < N; i++ )
608 if( ldpcG.nb1(i).size() != j )
609 regular = false;
610 for( size_t I = 0; I < K; I++ )
611 if( ldpcG.nb2(I).size() != k )
612 regular = false;
613 } while( !regular && !ldpcG.isConnected() );
614
615 // Convert to FactorGraph
616 vector<Factor> factors;
617 for( size_t I = 0; I < K; I++ ) {
618 VarSet vs;
619 for( size_t _i = 0; _i < k; _i++ ) {
620 size_t i = ldpcG.nb2(I)[_i];
621 vs |= Var( i, 2 );
622 }
623 factors.push_back( Factor( vs, paritycheck ) );
624 }
625 delete paritycheck;
626
627 // Generate noise vector
628 vector<char> noisebits(N,0);
629 size_t bitflips = 0;
630 for( size_t i = 0; i < N; i++ ) {
631 if( rnd_uniform() < noise ) {
632 noisebits[i] = 1;
633 bitflips++;
634 }
635 }
636 cout << "# bitflips = " << bitflips << endl;
637
638 // Simulate transmission of all-zero codeword
639 vector<char> input(N,0);
640 vector<char> output(N,0);
641 for( size_t i = 0; i < N; i++ )
642 output[i] = (input[i] + noisebits[i]) & 1;
643
644 // Add likelihoods
645 for( size_t i = 0; i < N; i++ )
646 factors.push_back( Factor(Var(i,2), likelihood + output[i]*2) );
647
648 // Construct Factor Graph
649 fg = FactorGraph( factors );
650 } else
651 throw "Invalid type";
652
653 // Output additional comments
654 cout << "# N = " << fg.nrVars() << endl;
655 cout << "# seed = " << seed << endl;
656
657 // Output factor graph
658 cout << fg;
659 } catch( const char *e ) {
660 /// Display error message
661 cerr << "Error: " << e << endl;
662 return 1;
663 }
664
665 return 0;
666 }