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 }