1 /* This file is part of libDAI - http://www.libdai.org/
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.
7 * Copyright (C) 2006-2009 Joris Mooij [joris dot mooij at libdai dot org]
8 * Copyright (C) 2006-2007 Radboud University Nijmegen, The Netherlands
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>
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
;
33 // w should be upper triangular or lower triangular
34 void WTh2FG( const matrix
&w
, const vector
<Real
> &th
, FactorGraph
&fg
) {
36 vector
<Factor
> factors
;
39 DAI_ASSERT( (w
.size1() == N
) && (w
.size2() == N
) );
42 for( size_t i
= 0; i
< N
; i
++ )
43 vars
.push_back(Var(i
,2));
45 factors
.reserve( w
.nnz() + N
);
46 // walk through the sparse array structure
47 // this is similar to matlab sparse arrays
48 // index2 gives the column index (similar to ir in matlab)
49 // index1 gives the starting indices for each row (similar to jc in matlab)
51 for( size_t pos
= 0; pos
< w
.nnz(); pos
++ ) {
52 while( pos
== w
.index1_data()[i
+1] )
54 size_t j
= w
.index2_data()[pos
];
55 Real w_ij
= w
.value_data()[pos
];
56 factors
.push_back( BinaryFactor( vars
[i
], vars
[j
], w_ij
) );
58 for( size_t i
= 0; i
< N
; i
++ )
59 factors
.push_back( BinaryFactor( vars
[i
], th
[i
] ) );
61 fg
= FactorGraph( factors
.begin(), factors
.end(), vars
.begin(), vars
.end(), factors
.size(), vars
.size() );
65 void MakeHOIFG( size_t N
, size_t M
, size_t k
, Real sigma
, FactorGraph
&fg
) {
67 vector
<Factor
> factors
;
70 for( size_t i
= 0; i
< N
; i
++ )
71 vars
.push_back(Var(i
,2));
73 for( size_t I
= 0; I
< M
; I
++ ) {
75 while( vars
.size() < k
) {
77 size_t newind
= (size_t)(N
* rnd_uniform());
78 Var newvar
= Var(newind
, 2);
79 if( !vars
.contains( newvar
) ) {
85 factors
.push_back( RandomFactor( vars
, sigma
) );
88 fg
= FactorGraph( factors
.begin(), factors
.end(), vars
.begin(), vars
.end(), factors
.size(), vars
.size() );
92 void MakeFullFG( size_t N
, Real mean_w
, Real mean_th
, Real sigma_w
, Real sigma_th
, FactorGraph
&fg
) {
93 matrix
w(N
,N
,N
*(N
-1)/2);
94 vector
<Real
> th(N
,0.0);
96 for( size_t i
= 0; i
< N
; i
++ ) {
97 for( size_t j
= i
+1; j
< N
; j
++ )
98 w(i
,j
) = rnd_stdnormal() * sigma_w
+ mean_w
;
99 th
[i
] = rnd_stdnormal() * sigma_th
+ mean_th
;
106 void Make3DPotts( size_t n1
, size_t n2
, size_t n3
, size_t states
, Real beta
, FactorGraph
&fg
) {
108 vector
<Factor
> factors
;
110 for( size_t i1
= 0; i1
< n1
; i1
++ )
111 for( size_t i2
= 0; i2
< n2
; i2
++ )
112 for( size_t i3
= 0; i3
< n3
; i3
++ ) {
113 vars
.push_back( Var( i1
*n2
*n3
+ i2
*n3
+ i3
, states
) );
115 factors
.push_back( PottsFactor( vars
.back(), vars
[ (i1
-1)*n2
*n3
+ i2
*n3
+ i3
], beta
) );
117 factors
.push_back( PottsFactor( vars
.back(), vars
[ i1
*n2
*n3
+ (i2
-1)*n3
+ i3
], beta
) );
119 factors
.push_back( PottsFactor( vars
.back(), vars
[ i1
*n2
*n3
+ i2
*n3
+ (i3
-1) ], beta
) );
122 fg
= FactorGraph( factors
.begin(), factors
.end(), vars
.begin(), vars
.end(), factors
.size(), vars
.size() );
126 void MakeGridFG( long periodic
, size_t n
, Real mean_w
, Real mean_th
, Real sigma_w
, Real sigma_th
, FactorGraph
&fg
) {
130 vector
<Real
> th(N
,0.0);
132 for( size_t i
= 0; i
< n
; i
++ )
133 for( size_t j
= 0; j
< n
; j
++ ) {
134 if( i
+1 < n
|| periodic
)
135 w(i
*n
+j
, ((i
+1)%n
)*n
+j
) = rnd_stdnormal() * sigma_w
+ mean_w
;
136 if( j
+1 < n
|| periodic
)
137 w(i
*n
+j
, i
*n
+((j
+1)%n
)) = rnd_stdnormal() * sigma_w
+ mean_w
;
138 th
[i
*n
+j
] = rnd_stdnormal() * sigma_th
+ mean_th
;
145 void MakeGridNonbinaryFG( bool periodic
, size_t n
, size_t states
, Real beta
, FactorGraph
&fg
) {
149 vector
<Factor
> factors
;
152 for( size_t i
= 0; i
< N
; i
++ )
153 vars
.push_back(Var(i
, states
));
155 factors
.reserve( 2 * N
);
156 for( size_t i
= 0; i
< n
; i
++ ) {
157 for( size_t j
= 0; j
< n
; j
++ ) {
158 if( i
+1 < n
|| periodic
)
159 factors
.push_back( RandomFactor( VarSet( vars
[i
*n
+j
], vars
[((i
+1)%n
)*n
+j
] ), beta
) );
160 if( j
+1 < n
|| periodic
)
161 factors
.push_back( RandomFactor( VarSet( vars
[i
*n
+j
], vars
[i
*n
+((j
+1)%n
)] ), beta
) );
165 fg
= FactorGraph( factors
.begin(), factors
.end(), vars
.begin(), vars
.end(), factors
.size(), vars
.size() );
169 void MakeLoopFG( size_t N
, Real mean_w
, Real mean_th
, Real sigma_w
, Real sigma_th
, FactorGraph
&fg
) {
171 vector
<Real
> th(N
,0.0);
173 for( size_t i
= 0; i
< N
; i
++ ) {
174 w(i
, (i
+1)%N
) = rnd_stdnormal() * sigma_w
+ mean_w
;
175 th
[i
] = rnd_stdnormal() * sigma_th
+ mean_th
;
182 void MakeLoopNonbinaryFG( size_t N
, size_t states
, Real beta
, FactorGraph
&fg
) {
184 vector
<Factor
> factors
;
187 for( size_t i
= 0; i
< N
; i
++ )
188 vars
.push_back(Var(i
, states
));
190 factors
.reserve( N
);
191 for( size_t i
= 0; i
< N
; i
++ ) {
192 factors
.push_back( RandomFactor( VarSet( vars
[i
], vars
[(i
+1)%N
] ), beta
) );
195 fg
= FactorGraph( factors
.begin(), factors
.end(), vars
.begin(), vars
.end(), factors
.size(), vars
.size() );
199 void MakeTreeFG( size_t N
, Real mean_w
, Real mean_th
, Real sigma_w
, Real sigma_th
, FactorGraph
&fg
) {
201 vector
<Real
> th(N
,0.0);
203 for( size_t i
= 0; i
< N
; i
++ ) {
204 th
[i
] = rnd_stdnormal() * sigma_th
+ mean_th
;
206 size_t j
= rnd_int( 0, i
-1 );
207 w(i
,j
) = rnd_stdnormal() * sigma_w
+ mean_w
;
215 void MakeDRegFG( size_t N
, size_t d
, Real mean_w
, Real mean_th
, Real sigma_w
, Real sigma_th
, FactorGraph
&fg
) {
216 matrix
w(N
,N
,(d
*N
)/2);
217 vector
<Real
> th(N
,0.0);
219 Graph g
= RandomDRegularGraph( N
, d
);
220 foreach( const UEdge
&e
, g
)
221 w(e
.n1
, e
.n2
) = rnd_stdnormal() * sigma_w
+ mean_w
;
223 for( size_t i
= 0; i
< N
; i
++ )
224 th
[i
] = rnd_stdnormal() * sigma_th
+ mean_th
;
230 // N = number of variables
231 // n = size of variable neighborhoods
232 // K = number of factors
233 // k = size of factor neighborhoods
234 // asserts: N * n == K * k
235 BipartiteGraph
CreateRandomBipartiteGraph( size_t N
, size_t K
, size_t n
, size_t k
) {
238 DAI_ASSERT( N
* n
== K
* k
);
240 // build lists of degree-repeated vertex numbers
241 std::vector
<size_t> stubs1(N
*n
,0);
242 for( size_t i
= 0; i
< N
; i
++ )
243 for( size_t t
= 0; t
< n
; t
++ )
246 // build lists of degree-repeated vertex numbers
247 std::vector
<size_t> stubs2(K
*k
,0);
248 for( size_t I
= 0; I
< K
; I
++ )
249 for( size_t t
= 0; t
< k
; t
++ )
253 random_shuffle( stubs1
.begin(), stubs1
.end() );
254 random_shuffle( stubs2
.begin(), stubs2
.end() );
257 vector
<BipartiteGraph::Edge
> edges
;
258 edges
.reserve( N
*n
);
259 for( size_t e
= 0; e
< N
*n
; e
++ )
260 edges
.push_back( BipartiteGraph::Edge(stubs1
[e
], stubs2
[e
]) );
262 // finish construction
263 G
.construct( N
, K
, edges
.begin(), edges
.end() );
269 // Returns x**n % p, assuming p is prime
270 int powmod (int x
, int n
, int p
) {
272 for( int m
= 0; m
< n
; m
++ )
278 // Returns order of x in GF(p) with p prime
279 size_t order (int x
, int p
) {
281 DAI_ASSERT( x
!= 0 );
285 prod
= (prod
* x
) % p
;
287 } while( prod
!= 1 );
292 // Returns whether n is a prime number
293 bool isPrime (size_t n
) {
295 for( size_t k
= 2; (k
< n
) && result
; k
++ )
302 // Make a regular LDPC graph with N=6, j=2, K=4, k=3
303 BipartiteGraph
CreateSmallLDPCGraph() {
305 size_t N
=4, j
=3, K
=4; // k=3;
307 typedef BipartiteGraph::Edge Edge
;
309 edges
.reserve( N
*j
);
310 edges
.push_back( Edge(0,0) ); edges
.push_back( Edge(1,0) ); edges
.push_back( Edge(2,0) );
311 edges
.push_back( Edge(0,1) ); edges
.push_back( Edge(1,1) ); edges
.push_back( Edge(3,1) );
312 edges
.push_back( Edge(0,2) ); edges
.push_back( Edge(2,2) ); edges
.push_back( Edge(3,2) );
313 edges
.push_back( Edge(1,3) ); edges
.push_back( Edge(2,3) ); edges
.push_back( Edge(3,3) );
315 // finish construction
316 G
.construct( N
, K
, edges
.begin(), edges
.end() );
322 // Use construction described in "A Class of Group-Structured LDPC Codes"
323 // by R. M. Tanner, D. Sridhara and T. Fuja
324 // Proceedings of ICSTA, 2001
326 // Example parameters: (p,j,k) = (31,3,5)
327 // j and k must be divisors of p-1
328 BipartiteGraph
CreateGroupStructuredLDPCGraph( size_t p
, size_t j
, size_t k
) {
336 for( a
= 2; a
< p
; a
++ )
337 if( order(a
,p
) == k
)
339 DAI_ASSERT( a
!= p
);
340 for( b
= 2; b
< p
; b
++ )
341 if( order(b
,p
) == j
)
343 DAI_ASSERT( b
!= p
);
344 // cout << "# order(a=" << a << ") = " << order(a,p) << endl;
345 // cout << "# order(b=" << b << ") = " << order(b,p) << endl;
347 DAI_ASSERT( N
* n
== K
* k
);
349 typedef BipartiteGraph::Edge Edge
;
351 edges
.reserve( N
* n
);
353 for( size_t s
= 0; s
< j
; s
++ )
354 for( size_t t
= 0; t
< k
; t
++ ) {
355 size_t P
= (powmod(b
,s
,p
) * powmod(a
,t
,p
)) % p
;
356 for( size_t m
= 0; m
< p
; m
++ )
357 edges
.push_back( Edge(t
*p
+ m
, s
*p
+ ((m
+ P
) % p
)) );
360 // finish construction
361 G
.construct( N
, K
, edges
.begin(), edges
.end() );
367 // Make parity check table
368 void MakeParityCheck( Real
*result
, size_t n
, Real eps
) {
370 for( size_t i
= 0; i
< N
; i
++ ) {
372 for( size_t t
= 0; t
< n
; t
++ )
378 result
[i
] = 1.0 - eps
;
384 const char *FULL_TYPE
= "full";
385 const char *GRID_TYPE
= "grid";
386 const char *GRID_TORUS_TYPE
= "grid_torus";
387 const char *DREG_TYPE
= "dreg";
388 const char *HOI_TYPE
= "hoi";
389 const char *LDPC_RANDOM_TYPE
= "ldpc_random";
390 const char *LDPC_GROUP_TYPE
= "ldpc_group";
391 const char *LDPC_SMALL_TYPE
= "ldpc_small";
392 const char *LOOP_TYPE
= "loop";
393 const char *TREE_TYPE
= "tree";
394 const char *POTTS3D_TYPE
= "potts3d";
397 int main( int argc
, char *argv
[] ) {
399 size_t N
, K
, k
, d
, j
, n1
, n2
, n3
;
402 Real beta
, sigma_w
, sigma_th
, noise
, mean_w
, mean_th
;
406 // Declare the supported options.
407 po::options_description
desc("Allowed options");
409 ("help", "produce help message")
410 ("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'")
411 ("seed", po::value
<size_t>(&seed
), "random number seed (tries to read from /dev/urandom if not specified)")
412 ("N", po::value
<size_t>(&N
), "number of variables (not for type=='ldpc_small')")
413 ("n1", po::value
<size_t>(&n1
), "width of 3D grid (only for type=='potts3d')")
414 ("n2", po::value
<size_t>(&n2
), "height of 3D grid (only for type=='potts3d')")
415 ("n3", po::value
<size_t>(&n3
), "length of 3D grid (only for type=='potts3d')")
416 ("K", po::value
<size_t>(&K
), "number of factors\n\t(only for type=='hoi' and 'type=='ldpc_{random,group}')")
417 ("k", po::value
<size_t>(&k
), "number of variables per factor\n\t(only for type=='hoi' and type=='ldpc_{random,group}')")
418 ("d", po::value
<size_t>(&d
), "variable connectivity\n\t(only for type=='dreg')")
419 ("j", po::value
<size_t>(&j
), "number of parity checks per bit\n\t(only for type=='ldpc_{random,group}')")
420 ("prime", po::value
<size_t>(&prime
), "prime number for construction of LDPC code\n\t(only for type=='ldpc_group')")
421 ("beta", po::value
<Real
>(&beta
), "stddev of log-factor entries\n\t(only for type=='hoi', 'potts3d', 'grid' if states>2)")
422 ("mean_w", po::value
<Real
>(&mean_w
), "mean of pairwise interactions w_{ij}\n\t(not for type=='hoi', 'ldpc_*', 'potts3d')")
423 ("mean_th", po::value
<Real
>(&mean_th
), "mean of singleton interactions th_i\n\t(not for type=='hoi', 'ldpc_*', 'potts3d')")
424 ("sigma_w", po::value
<Real
>(&sigma_w
), "stddev of pairwise interactions w_{ij}\n\t(not for type=='hoi', 'ldpc_*', 'potts3d')")
425 ("sigma_th", po::value
<Real
>(&sigma_th
), "stddev of singleton interactions th_i\n\t(not for type=='hoi', 'ldpc_*', 'potts3d'")
426 ("noise", po::value
<Real
>(&noise
), "bitflip probability for binary symmetric channel (only for type=='ldpc')")
427 ("states", po::value
<size_t>(&states
), "number of states of each variable (should be 2 for all but type=='grid', 'grid_torus', 'loop', 'potts3d')")
430 po::variables_map vm
;
431 po::store(po::parse_command_line(argc
, argv
, desc
), vm
);
434 if( vm
.count("help") || !vm
.count("type") ) {
435 if( vm
.count("type") ) {
436 if( type
== FULL_TYPE
) {
437 cout
<< "Creates fully connected pairwise graphical model of <N> binary variables;" << endl
;
438 } else if( type
== GRID_TYPE
) {
439 cout
<< "Creates (non-periodic) 2D Ising grid of (approx.) <N> variables (which need not be binary);" << endl
;
440 } else if( type
== GRID_TORUS_TYPE
) {
441 cout
<< "Creates periodic 2D Ising grid of (approx.) <N> variables (which need not be binary);" << endl
;
442 } else if( type
== DREG_TYPE
) {
443 cout
<< "Creates random d-regular graph of <N> binary variables with uniform degree <d>" << endl
;
444 cout
<< "(where <d><N> should be even);" << endl
;
445 } else if( type
== LOOP_TYPE
) {
446 cout
<< "Creates a pairwise graphical model consisting of a single loop of" << endl
;
447 cout
<< "<N> variables (which need not be binary);" << endl
;
448 } else if( type
== TREE_TYPE
) {
449 cout
<< "Creates a pairwise, connected graphical model without cycles (i.e., a tree)" << endl
;
450 cout
<< "of <N> binary variables;" << endl
;
451 } else if( type
== HOI_TYPE
) {
452 cout
<< "Creates a random factor graph of <N> binary variables and" << endl
;
453 cout
<< "<K> factors, each factor being an interaction of <k> variables." << endl
;
454 cout
<< "The entries of the factors are exponentials of i.i.d. Gaussian" << endl
;
455 cout
<< "variables with mean 0 and standard deviation <beta>." << endl
;
456 } else if( type
== LDPC_RANDOM_TYPE
) {
457 cout
<< "Simulates LDPC decoding problem, using a LDPC code of <N> bits and <K> parity" << endl
;
458 cout
<< "checks, with <k> bits per check and <j> checks per bit, transmitted on a binary" << endl
;
459 cout
<< "symmetric channel with probability <noise> of flipping a bit. The transmitted" << endl
;
460 cout
<< "codeword has all bits set to zero. The LDPC code is randomly generated." << endl
;
461 } else if( type
== LDPC_GROUP_TYPE
) {
462 cout
<< "Simulates LDPC decoding problem, using a LDPC code of <N> bits and <K> parity" << endl
;
463 cout
<< "checks, with <k> bits per check and <j> checks per bit, transmitted on a binary" << endl
;
464 cout
<< "symmetric channel with probability <noise> of flipping a bit. The transmitted" << endl
;
465 cout
<< "codeword has all bits set to zero. The LDPC code is constructed (using group" << endl
;
466 cout
<< "theory) using a parameter <prime>; <j> and <k> should both be divisors of <prime>-1." << endl
;
467 } else if( type
== LDPC_SMALL_TYPE
) {
468 cout
<< "Simulates LDPC decoding problem, using a LDPC code of 4 bits and 4 parity" << endl
;
469 cout
<< "checks, with 3 bits per check and 3 checks per bit, transmitted on a binary" << endl
;
470 cout
<< "symmetric channel with probability <noise> of flipping a bit. The transmitted" << endl
;
471 cout
<< "codeword has all bits set to zero. The LDPC code is fixed." << endl
;
472 } else if( type
== POTTS3D_TYPE
) {
473 cout
<< "Builds 3D Potts model of size <n1>x<n2>x<n3> with nearest-neighbour Potts" << endl
;
474 cout
<< "interactions with <states> states and inverse temperature <beta>." << endl
;
476 cerr
<< "Unknown type (should be one of 'full', 'grid', 'grid_torus', 'dreg', 'loop', 'tree', 'hoi', 'ldpc_random', 'ldpc_group', 'ldpc_small', 'potts3d')" << endl
;
478 if( type
== FULL_TYPE
|| type
== GRID_TYPE
|| type
== GRID_TORUS_TYPE
|| type
== DREG_TYPE
|| type
== LOOP_TYPE
|| type
== TREE_TYPE
) {
479 if( type
== GRID_TYPE
|| type
== GRID_TORUS_TYPE
|| type
== LOOP_TYPE
) {
480 cout
<< "if <states> > 2: factor entries are exponents of Gaussians with mean 0 and standard deviation beta; otherwise," << endl
;
482 cout
<< "singleton interactions are Gaussian with mean <mean_th> and standard" << endl
;
483 cout
<< "deviation <sigma_th>; pairwise interactions are Gaussian with mean" << endl
;
484 cout
<< "<mean_w> and standard deviation <sigma_w>." << endl
;
487 cout
<< endl
<< desc
<< endl
;
491 if( !vm
.count("states") )
494 if( !vm
.count("seed") ) {
497 infile
.open( "/dev/urandom" );
498 success
= infile
.is_open();
500 infile
.read( (char *)&seed
, sizeof(size_t) / sizeof(char) );
501 success
= infile
.good();
505 throw "Please specify random number seed.";
511 cout
<< "# Factor graph made by " << argv
[0] << endl
;
512 cout
<< "# type = " << type
<< endl
;
514 if( type
== FULL_TYPE
) {
515 if( !vm
.count("N") || !vm
.count("mean_w") || !vm
.count("mean_th") || !vm
.count("sigma_w") || !vm
.count("sigma_th") )
516 throw "Please specify all required arguments";
517 MakeFullFG( N
, mean_w
, mean_th
, sigma_w
, sigma_th
, fg
);
519 cout
<< "# N = " << N
<< endl
;
520 cout
<< "# mean_w = " << mean_w
<< endl
;
521 cout
<< "# mean_th = " << mean_th
<< endl
;
522 cout
<< "# sigma_w = " << sigma_w
<< endl
;
523 cout
<< "# sigma_th = " << sigma_th
<< endl
;
524 } else if( type
== GRID_TYPE
|| type
== GRID_TORUS_TYPE
) {
525 #define NEED_ARG(name, desc) do { if(!vm.count(name)) throw "Please specify " desc " with --" name; } while(0);
527 NEED_ARG("N", "number of nodes");
528 NEED_ARG("beta", "stddev of log-factor entries");
530 NEED_ARG("N", "number of nodes");
531 NEED_ARG("mean_w", "mean of pairwise interactions");
532 NEED_ARG("mean_th", "mean of singleton interactions");
533 NEED_ARG("sigma_w", "stddev of pairwise interactions");
534 NEED_ARG("sigma_th", "stddev of singleton interactions");
537 size_t n
= (size_t)sqrt((long double)N
);
540 bool periodic
= false;
541 if( type
== GRID_TYPE
)
547 MakeGridNonbinaryFG( periodic
, n
, states
, beta
, fg
);
549 MakeGridFG( periodic
, n
, mean_w
, mean_th
, sigma_w
, sigma_th
, fg
);
551 cout
<< "# n = " << n
<< endl
;
552 cout
<< "# N = " << N
<< endl
;
555 cout
<< "# beta = " << beta
<< 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
;
562 } else if( type
== DREG_TYPE
) {
563 if( !vm
.count("N") || !vm
.count("mean_w") || !vm
.count("mean_th") || !vm
.count("sigma_w") || !vm
.count("sigma_th") || !vm
.count("d") )
564 throw "Please specify all required arguments";
566 MakeDRegFG( N
, d
, mean_w
, mean_th
, sigma_w
, sigma_th
, fg
);
568 cout
<< "# N = " << N
<< endl
;
569 cout
<< "# d = " << d
<< endl
;
570 cout
<< "# mean_w = " << mean_w
<< endl
;
571 cout
<< "# mean_th = " << mean_th
<< endl
;
572 cout
<< "# sigma_w = " << sigma_w
<< endl
;
573 cout
<< "# sigma_th = " << sigma_th
<< endl
;
574 } else if( type
== LOOP_TYPE
) {
576 if( !vm
.count("N") || !vm
.count("beta") )
577 throw "Please specify all required arguments";
579 if( !vm
.count("N") || !vm
.count("mean_w") || !vm
.count("mean_th") || !vm
.count("sigma_w") || !vm
.count("sigma_th") )
580 throw "Please specify all required arguments";
583 MakeLoopNonbinaryFG( N
, states
, beta
, fg
);
585 MakeLoopFG( N
, mean_w
, mean_th
, sigma_w
, sigma_th
, fg
);
587 cout
<< "# N = " << N
<< endl
;
590 cout
<< "# beta = " << beta
<< endl
;
592 cout
<< "# mean_w = " << mean_w
<< endl
;
593 cout
<< "# mean_th = " << mean_th
<< endl
;
594 cout
<< "# sigma_w = " << sigma_w
<< endl
;
595 cout
<< "# sigma_th = " << sigma_th
<< endl
;
597 } else if( type
== TREE_TYPE
) {
598 if( !vm
.count("N") || !vm
.count("mean_w") || !vm
.count("mean_th") || !vm
.count("sigma_w") || !vm
.count("sigma_th") )
599 throw "Please specify all required arguments";
600 MakeTreeFG( N
, mean_w
, mean_th
, sigma_w
, sigma_th
, fg
);
602 cout
<< "# N = " << N
<< endl
;
603 cout
<< "# mean_w = " << mean_w
<< endl
;
604 cout
<< "# mean_th = " << mean_th
<< endl
;
605 cout
<< "# sigma_w = " << sigma_w
<< endl
;
606 cout
<< "# sigma_th = " << sigma_th
<< endl
;
607 } else if( type
== HOI_TYPE
) {
608 if( !vm
.count("N") || !vm
.count("K") || !vm
.count("k") || !vm
.count("beta") )
609 throw "Please specify all required arguments";
611 MakeHOIFG( N
, K
, k
, beta
, fg
);
612 } while( !fg
.isConnected() );
614 cout
<< "# N = " << N
<< endl
;
615 cout
<< "# K = " << K
<< endl
;
616 cout
<< "# k = " << k
<< endl
;
617 cout
<< "# beta = " << beta
<< endl
;
618 } else if( type
== LDPC_RANDOM_TYPE
|| type
== LDPC_GROUP_TYPE
|| type
== LDPC_SMALL_TYPE
) {
619 if( !vm
.count("noise") )
620 throw "Please specify all required arguments";
622 if( type
== LDPC_RANDOM_TYPE
) {
623 if( !vm
.count("N") || !vm
.count("K") || !vm
.count("j") || !vm
.count("k") )
624 throw "Please specify all required arguments";
627 throw "Parameters should satisfy N * j == K * k";
628 } else if( type
== LDPC_GROUP_TYPE
) {
629 if( !vm
.count("prime") || !vm
.count("j") || !vm
.count("k") )
630 throw "Please specify all required arguments";
632 if( !isPrime(prime
) )
633 throw "Parameter <prime> should be prime";
634 if( !((prime
-1) % j
== 0 ) )
635 throw "Parameters should satisfy (prime-1) % j == 0";
636 if( !((prime
-1) % k
== 0 ) )
637 throw "Parameters should satisfy (prime-1) % k == 0";
641 } else if( type
== LDPC_SMALL_TYPE
) {
648 cout
<< "# N = " << N
<< endl
;
649 cout
<< "# K = " << K
<< endl
;
650 cout
<< "# j = " << j
<< endl
;
651 cout
<< "# k = " << k
<< endl
;
652 if( type
== LDPC_GROUP_TYPE
)
653 cout
<< "# prime = " << prime
<< endl
;
654 cout
<< "# noise = " << noise
<< endl
;
656 // p = 31, j = 3, k = 5
657 // p = 37, j = 3, k = 4
658 // p = 7 , j = 2, k = 3
659 // p = 29, j = 2, k = 4
661 // Construct likelihood and paritycheck factors
662 Real likelihood
[4] = {1.0 - noise
, noise
, noise
, 1.0 - noise
};
663 Real
*paritycheck
= new Real
[1 << k
];
664 MakeParityCheck(paritycheck
, k
, 0.0);
666 // Create LDPC structure
667 BipartiteGraph ldpcG
;
670 if( type
== LDPC_GROUP_TYPE
)
671 ldpcG
= CreateGroupStructuredLDPCGraph( prime
, j
, k
);
672 else if( type
== LDPC_RANDOM_TYPE
)
673 ldpcG
= CreateRandomBipartiteGraph( N
, K
, j
, k
);
674 else if( type
== LDPC_SMALL_TYPE
)
675 ldpcG
= CreateSmallLDPCGraph();
678 for( size_t i
= 0; i
< N
; i
++ )
679 if( ldpcG
.nb1(i
).size() != j
)
681 for( size_t I
= 0; I
< K
; I
++ )
682 if( ldpcG
.nb2(I
).size() != k
)
684 } while( !regular
&& !ldpcG
.isConnected() );
686 // Convert to FactorGraph
687 vector
<Factor
> factors
;
688 for( size_t I
= 0; I
< K
; I
++ ) {
690 for( size_t _i
= 0; _i
< k
; _i
++ ) {
691 size_t i
= ldpcG
.nb2(I
)[_i
];
694 factors
.push_back( Factor( vs
, paritycheck
) );
698 // Generate noise vector
699 vector
<char> noisebits(N
,0);
701 for( size_t i
= 0; i
< N
; i
++ ) {
702 if( rnd_uniform() < noise
) {
707 cout
<< "# bitflips = " << bitflips
<< endl
;
709 // Simulate transmission of all-zero codeword
710 vector
<char> input(N
,0);
711 vector
<char> output(N
,0);
712 for( size_t i
= 0; i
< N
; i
++ )
713 output
[i
] = (input
[i
] + noisebits
[i
]) & 1;
716 for( size_t i
= 0; i
< N
; i
++ )
717 factors
.push_back( Factor(Var(i
,2), likelihood
+ output
[i
]*2) );
719 // Construct Factor Graph
720 fg
= FactorGraph( factors
);
721 } else if( type
== POTTS3D_TYPE
) {
722 if( !vm
.count("n1") || !vm
.count("n2") || !vm
.count("n3") || !vm
.count("beta") || !vm
.count("states") )
723 throw "Please specify all required arguments";
724 Make3DPotts( n1
, n2
, n3
, states
, beta
, fg
);
726 cout
<< "# N = " << n1
*n2
*n3
<< endl
;
727 cout
<< "# n1 = " << n1
<< endl
;
728 cout
<< "# n2 = " << n2
<< endl
;
729 cout
<< "# n3 = " << n3
<< endl
;
730 cout
<< "# beta = " << beta
<< endl
;
731 cout
<< "# states = " << states
<< endl
;
733 throw "Invalid type";
736 cout
<< "# seed = " << seed
<< endl
;
738 } catch( const char *e
) {
739 cerr
<< "Error: " << e
<< endl
;