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