From 3dbbc7c5bb1a87089a1d9c29e4481b4559cd446e Mon Sep 17 00:00:00 2001 From: Joris Mooij Date: Thu, 18 Sep 2008 10:40:52 +0200 Subject: [PATCH] Now compiles also with Visual Studio 2008 under Windows (still buggy!) - Merged utils/createfg from SVN head - Added Makefile.win for Visual Studio and fixed some bugs; HAK doesn't work yet, because of problems with NANs. --- Makefile.win | 242 ++++++++++++ include/dai/factorgraph.h | 1 - include/dai/lc.h | 2 +- include/dai/mr.h | 4 +- include/dai/util.h | 24 +- include/dai/var.h | 2 +- src/factorgraph.cpp | 3 +- src/util.cpp | 17 +- src/x2x.cpp | 2 +- tests/test.cpp | 5 +- utils/createfg.cpp | 805 +++++++++++++++++++++++++++++--------- 11 files changed, 911 insertions(+), 196 deletions(-) create mode 100755 Makefile.win diff --git a/Makefile.win b/Makefile.win new file mode 100755 index 0000000..7a12f2e --- /dev/null +++ b/Makefile.win @@ -0,0 +1,242 @@ +# Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl] +# Radboud University Nijmegen, The Netherlands +# +# This file is part of libDAI. +# +# libDAI is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 2 of the License, or +# (at your option) any later version. +# +# libDAI is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with libDAI; if not, write to the Free Software +# Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + + +# Enable/disable various approximate inference methods +WITH_BP=true +WITH_MF=true +WITH_HAK=true +WITH_LC=true +WITH_TREEEP=true +WITH_JTREE=true +WITH_MR=true +# Build with debug info? +DEBUG=true +# Build matlab interface? +#WITH_MATLAB=true +# New/old matlab version? +NEW_MATLAB=true +# Windows or linux (default)? +WINDOWS=true + +# Directories +INC=include/dai +SRC=src +LIB=lib + +# For some reason, we have to set the library path, although it is in the environment +LIBS=/LIBPATH:"C:\Program Files\Microsoft Visual Studio 9.0\VC\ATLMFC\LIB" /LIBPATH:"C:\Program Files\Microsoft Visual Studio 9.0\VC\LIB" /LIBPATH:"C:\Program Files\Microsoft SDKs\Windows\v6.0A\lib" +# We use the BOOST Program Options library +BOOSTLIBS=/LIBPATH:C:\boost_1_36_0\stage\lib + +# Compile using GNU C++ Compiler +CC=cl +OBJEXT=obj +LIBEXT=lib + +# Flags for the C++ compiler +CCFLAGS=/I./include /IC:\boost_1_36_0 /EHsc /Ox +!IFDEF DEBUG +CCFLAGS=$(CCFLAGS) /Zi /DDAI_DEBUG +!ELSE +CCFLAGS=$(CCFLAGS) +!ENDIF + +!IFDEF WINDOWS +CCFLAGS=$(CCFLAGS) /DWINDOWS +!ENDIF + +!IFDEF WITH_BP +CCFLAGS=$(CCFLAGS) /DWITH_BP +OBJECTS=$(OBJECTS) bp.obj +!ENDIF +!IFDEF WITH_MF +CCFLAGS=$(CCFLAGS) /DWITH_MF +OBJECTS=$(OBJECTS) mf.obj +!ENDIF +!IFDEF WITH_HAK +CCFLAGS=$(CCFLAGS) /DWITH_HAK +OBJECTS=$(OBJECTS) hak.obj +!ENDIF +!IFDEF WITH_LC +CCFLAGS=$(CCFLAGS) /DWITH_LC +OBJECTS=$(OBJECTS) lc.obj +!ENDIF +!IFDEF WITH_TREEEP +CCFLAGS=$(CCFLAGS) /DWITH_TREEEP +OBJECTS=$(OBJECTS) treeep.obj +!ENDIF +!IFDEF WITH_JTREE +CCFLAGS=$(CCFLAGS) /DWITH_JTREE +OBJECTS=$(OBJECTS) jtree.obj +!ENDIF +!IFDEF WITH_MR +CCFLAGS=$(CCFLAGS) /DWITH_MR +OBJECTS=$(OBJECTS) mr.obj +!ENDIF + +!IFDEF WITH_MATLAB +# Replace the following by the directory where Matlab has been installed +MATLABDIR = /opt/matlab/bin +# Replace the following with the extension of compiled MEX files on this platform, e.g. .mexglx for x86 +MEXEXT = .mexglx +MEX = $(MATLABDIR)/mex +MEXFLAGS = -I. +!IFDEF DEBUG +MEXFLAGS = $(MEXFLAGS) -g /DDAI_DEBUG +!ENDIF +!IFDEF NEW_MATLAB +MEXFLAGS = $(MEXFLAGS) -largeArrayDims +!ELSE +MEXFLAGS = $(MEXFLAGS) /DSMALLMEM +!ENDIF +!ENDIF + +HEADERS = $(INC)/bipgraph.h $(INC)/diffs.h $(INC)/index.h $(INC)/var.h $(INC)/factor.h $(INC)/varset.h $(INC)/prob.h $(INC)/daialg.h $(INC)/properties.h $(INC)/alldai.h $(INC)/enum.h $(INC)/x2x.h + +TARGETS = tests utils $(LIB)/libdai.lib example.exe +# testregression disabled, it uses sh and awk +# doc disabled, it uses doxygen +!IFDEF WITH_MATLAB +TARGETS = $(TARGETS) matlabs +!ENDIF +all : $(TARGETS) + +matlabs : matlab/dai.$(MEXEXT) matlab/dai_readfg.$(MEXEXT) matlab/dai_writefg.$(MEXEXT) matlab/dai_removeshortloops.$(MEXEXT) matlab/dai_potstrength.$(MEXEXT) + +$(LIB)/libdai.lib : bipgraph.obj daialg.obj alldai.obj clustergraph.obj factorgraph.obj properties.obj regiongraph.obj util.obj weightedgraph.obj x2x.obj $(OBJECTS) + lib /out:$(LIB)/libdai.lib bipgraph.obj daialg.obj alldai.obj clustergraph.obj factorgraph.obj properties.obj regiongraph.obj util.obj weightedgraph.obj x2x.obj $(OBJECTS) + +tests : tests/test.exe + +utils : utils/createfg.exe utils/fg2dot.exe utils/remove_short_loops.exe utils/fginfo.exe + +testregression : tests/test + echo Testing...this can take a while... + cd tests; time ./testregression; cd .. + +doc : $(INC)/*.h $(SRC)/*.cpp doxygen.conf + doxygen doxygen.conf + +clean : + del *.obj example.exe matlab\*.$(MEXEXT) matlab\*.obj tests/test.exe utils\fg2dot.exe utils\createfg.exe utils\remove_short_loops.exe utils\fginfo.exe $(LIB)\libdai.lib + + +bipgraph.obj : $(SRC)/bipgraph.cpp $(HEADERS) + $(CC) $(CCFLAGS) -c $(SRC)/bipgraph.cpp + +daialg.obj : $(SRC)/daialg.cpp $(HEADERS) + $(CC) $(CCFLAGS) -c $(SRC)/daialg.cpp + +bp.obj : $(SRC)/bp.cpp $(INC)/bp.h $(HEADERS) + $(CC) $(CCFLAGS) -c $(SRC)/bp.cpp + +lc.obj : $(SRC)/lc.cpp $(INC)/lc.h $(HEADERS) + $(CC) $(CCFLAGS) -c $(SRC)/lc.cpp + +mf.obj : $(SRC)/mf.cpp $(INC)/mf.h $(HEADERS) + $(CC) $(CCFLAGS) -c $(SRC)/mf.cpp + +factorgraph.obj : $(SRC)/factorgraph.cpp $(INC)/factorgraph.h $(HEADERS) + $(CC) $(CCFLAGS) -c $(SRC)/factorgraph.cpp + +util.obj : $(SRC)/util.cpp $(INC)/util.h $(HEADERS) + $(CC) $(CCFLAGS) -c $(SRC)/util.cpp + +regiongraph.obj : $(SRC)/regiongraph.cpp $(INC)/regiongraph.h $(HEADERS) + $(CC) $(CCFLAGS) -c $(SRC)/regiongraph.cpp + +hak.obj : $(SRC)/hak.cpp $(INC)/hak.h $(HEADERS) $(INC)/regiongraph.h + $(CC) $(CCFLAGS) -c $(SRC)/hak.cpp + +clustergraph.obj : $(SRC)/clustergraph.cpp $(INC)/clustergraph.h $(HEADERS) + $(CC) $(CCFLAGS) -c $(SRC)/clustergraph.cpp + +jtree.obj : $(SRC)/jtree.cpp $(INC)/jtree.h $(HEADERS) $(INC)/weightedgraph.h $(INC)/clustergraph.h $(INC)/regiongraph.h + $(CC) $(CCFLAGS) -c $(SRC)/jtree.cpp + +treeep.obj : $(SRC)/treeep.cpp $(INC)/treeep.h $(HEADERS) $(INC)/weightedgraph.h $(INC)/clustergraph.h $(INC)/regiongraph.h $(INC)/jtree.h + $(CC) $(CCFLAGS) -c $(SRC)/treeep.cpp + +weightedgraph.obj : $(SRC)/weightedgraph.cpp $(INC)/weightedgraph.h $(HEADERS) + $(CC) $(CCFLAGS) -c $(SRC)/weightedgraph.cpp + +mr.obj : $(SRC)/mr.cpp $(INC)/mr.h $(HEADERS) + $(CC) $(CCFLAGS) -c $(SRC)/mr.cpp + +properties.obj : $(SRC)/properties.cpp $(HEADERS) + $(CC) $(CCFLAGS) -c $(SRC)/properties.cpp + +alldai.obj : $(SRC)/alldai.cpp $(HEADERS) + $(CC) $(CCFLAGS) -c $(SRC)/alldai.cpp + +x2x.obj : $(SRC)/x2x.cpp $(HEADERS) + $(CC) $(CCFLAGS) -c $(SRC)/x2x.cpp + + +# EXAMPLE +########## + +example.exe : example.cpp $(HEADERS) $(LIB)/libdai.lib + $(CC) $(CCFLAGS) /Feexample example.cpp /link $(LIB)/libdai.lib $(LIBS) + + +# TESTS +######## + +tests/test.exe : tests/test.cpp $(HEADERS) $(LIB)/libdai.lib + $(CC) $(CCFLAGS) /Fetests/test tests/test.cpp /link $(LIB)/libdai.lib $(LIBS) $(BOOSTLIBS) + + +# MATLAB INTERFACE +################### + +matlab/dai.$(MEXEXT) : matlab/dai.cpp $(HEADERS) matlab/matlab.obj $(LIB)/libdai.lib + $(MEX) $(MEXFLAGS) -o matlab/dai matlab/dai.cpp matlab/matlab.obj $(LIB)/libdai.lib + +matlab/dai_readfg.$(MEXEXT) : matlab/dai_readfg.cpp $(HEADERS) factorgraph.obj matlab/matlab.obj + $(MEX) $(MEXFLAGS) -o matlab/dai_readfg matlab/dai_readfg.cpp factorgraph.obj matlab/matlab.obj + +matlab/dai_writefg.$(MEXEXT) : matlab/dai_writefg.cpp $(HEADERS) factorgraph.obj matlab/matlab.obj + $(MEX) $(MEXFLAGS) -o matlab/dai_writefg matlab/dai_writefg.cpp factorgraph.obj matlab/matlab.obj + +matlab/dai_removeshortloops.$(MEXEXT) : matlab/dai_removeshortloops.cpp $(HEADERS) factorgraph.obj matlab/matlab.obj + $(MEX) $(MEXFLAGS) -o matlab/dai_removeshortloops matlab/dai_removeshortloops.cpp factorgraph.obj matlab/matlab.obj + +matlab/dai_potstrength.$(MEXEXT) : matlab/dai_potstrength.cpp $(HEADERS) matlab/matlab.obj + $(MEX) $(MEXFLAGS) -o matlab/dai_potstrength matlab/dai_potstrength.cpp matlab/matlab.obj + +matlab/matlab.obj : matlab/matlab.cpp matlab/matlab.h $(HEADERS) + $(MEX) $(MEXFLAGS) -outdir matlab -c matlab/matlab.cpp + + +# UTILS +######## + +utils/createfg.exe : utils/createfg.cpp $(HEADERS) factorgraph.obj weightedgraph.obj util.obj bipgraph.obj + $(CC) $(CCFLAGS) /Feutils/createfg utils/createfg.cpp factorgraph.obj weightedgraph.obj util.obj bipgraph.obj /link $(LIBS) $(BOOSTLIBS) + +utils/fg2dot.exe : utils/fg2dot.cpp $(HEADERS) factorgraph.obj + $(CC) $(CCFLAGS) /Feutils/fg2dot utils/fg2dot.cpp factorgraph.obj /link $(LIBS) + +utils/remove_short_loops.exe : utils/remove_short_loops.cpp $(HEADERS) factorgraph.obj + $(CC) $(CCFLAGS) /Feutils/remove_short_loops utils/remove_short_loops.cpp factorgraph.obj /link $(LIBS) + +utils/fginfo.exe : utils/fginfo.cpp $(HEADERS) factorgraph.obj bipgraph.obj + $(CC) $(CCFLAGS) /Feutils/fginfo utils/fginfo.cpp factorgraph.obj bipgraph.obj /link $(LIBS) diff --git a/include/dai/factorgraph.h b/include/dai/factorgraph.h index 77f5f8a..a502ae8 100644 --- a/include/dai/factorgraph.h +++ b/include/dai/factorgraph.h @@ -25,7 +25,6 @@ #include #include -#include #include #include diff --git a/include/dai/lc.h b/include/dai/lc.h index 231b043..7ded83e 100644 --- a/include/dai/lc.h +++ b/include/dai/lc.h @@ -81,7 +81,7 @@ class LC : public DAIAlgFG { std::string identify() const; Factor belief (const Var &n) const { return( _beliefs[findVar(n)] ); } - Factor belief (const VarSet &/*ns*/) const { assert( 0 == 1 ); } + Factor belief (const VarSet &/*ns*/) const { assert( 0 == 1 ); return Factor(); } std::vector beliefs() const { return _beliefs; } Complex logZ() const { return NAN; } void CalcBelief (size_t i); diff --git a/include/dai/mr.h b/include/dai/mr.h index 153b83a..eb0726d 100644 --- a/include/dai/mr.h +++ b/include/dai/mr.h @@ -71,7 +71,7 @@ class MR : public DAIAlgFG { void solveM(); double run(); Factor belief( const Var &n ) const; - Factor belief( const VarSet &/*ns*/ ) const { assert( 0 == 1 ); } + Factor belief( const VarSet &/*ns*/ ) const { assert( 0 == 1 ); return Factor(); } std::vector beliefs() const; Complex logZ() const { return NAN; } void init() { assert( checkProperties() ); } @@ -89,7 +89,7 @@ class MR : public DAIAlgFG { void sum_subs(size_t j, sub_nb A, double *sum_even, double *sum_odd); double sign(double a) { return (a >= 0) ? 1.0 : -1.0; } - MR* clone() const { assert( 0 == 1 ); } + MR* clone() const { return new MR(*this); } bool checkProperties(); diff --git a/include/dai/util.h b/include/dai/util.h index d83c204..4c093e0 100644 --- a/include/dai/util.h +++ b/include/dai/util.h @@ -30,24 +30,34 @@ #include #include + #ifdef WINDOWS - #include // for _isnan (?) - #include // for atanh - #include // for log1p -// #include + #include // for NAN + #include +#else + #include #endif #define foreach BOOST_FOREACH +#ifdef WINDOWS + bool isnan( double x ); + double atanh( double x ); + double log1p( double x ); +#endif + + namespace dai { #ifdef WINDOWS - bool isnan( double x ) inline { return _isnan( x ); } - using boost::math:atanh; - using boost::math::log1p; + template + class hash_map : public std::map {}; +#else + template + class hash_map : public std::tr1::unordered_map {}; #endif diff --git a/include/dai/var.h b/include/dai/var.h index e9deda1..7deaa95 100644 --- a/include/dai/var.h +++ b/include/dai/var.h @@ -31,7 +31,7 @@ namespace dai { /// Represents a discrete variable. -/* It contains the label of the variable, an integer-valued +/** It contains the label of the variable, an integer-valued * unique ID for that variable, and the number of possible * values ("states") of the variable. */ diff --git a/src/factorgraph.cpp b/src/factorgraph.cpp index b887316..38a6b5a 100644 --- a/src/factorgraph.cpp +++ b/src/factorgraph.cpp @@ -27,7 +27,6 @@ #include #include #include -#include #include #include @@ -62,7 +61,7 @@ FactorGraph::FactorGraph( const std::vector &P ) : G(), _undoProbs(), _n /// Part of constructors (creates edges, neighbours and adjacency matrix) void FactorGraph::createGraph( size_t nrEdges ) { // create a mapping for indices - std::tr1::unordered_map hashmap; + hash_map hashmap; for( size_t i = 0; i < vars.size(); i++ ) hashmap[var(i).label()] = i; diff --git a/src/util.cpp b/src/util.cpp index 18596ee..728f62f 100644 --- a/src/util.cpp +++ b/src/util.cpp @@ -27,6 +27,9 @@ #ifdef WINDOWS #include + #include // for atanh + #include // for log1p + #include // for _isnan #else // Assume POSIX compliant system. We need the following for querying the CPU time for this process #include @@ -34,6 +37,19 @@ #endif +#ifdef WINDOWS +bool isnan( double x ) { + return _isnan( x ); +} +double atanh( double x ) { + return boost::math::atanh( x ); +} +double log1p( double x ) { + return boost::math::log1p( x ); +} +#endif + + namespace dai { @@ -50,7 +66,6 @@ double toc() { #endif } - // This is a typedef for a random number generator. // Try boost::mt19937 or boost::ecuyer1988 instead of boost::minstd_rand typedef boost::minstd_rand _rnd_gen_type; diff --git a/src/x2x.cpp b/src/x2x.cpp index 51d026e..49b5352 100644 --- a/src/x2x.cpp +++ b/src/x2x.cpp @@ -24,7 +24,7 @@ CHANGES: #include #include - +#include namespace x2x { diff --git a/tests/test.cpp b/tests/test.cpp index 007e7f3..2ac19e4 100644 --- a/tests/test.cpp +++ b/tests/test.cpp @@ -25,6 +25,7 @@ #include #include #include +#include #include #include #include @@ -148,10 +149,10 @@ pair parseMethod( const string &_s, const map if( pos == string::npos ) throw "Malformed method"; size_t n = 0; - for( ; n < sizeof(DAINames) / sizeof(string); n++ ) + for( ; strlen( DAINames[n] ) != 0; n++ ) if( name == DAINames[n] ) break; - if( n == sizeof(DAINames) / sizeof(string) ) + if( strlen( DAINames[n] ) == 0 ) throw "Unknown inference algorithm"; stringstream ss; diff --git a/utils/createfg.cpp b/utils/createfg.cpp index cc5fe80..4feb946 100644 --- a/utils/createfg.cpp +++ b/utils/createfg.cpp @@ -20,8 +20,13 @@ #include +#include +#include #include +#include #include +#include +#include #include #include #include @@ -30,12 +35,53 @@ using namespace std; using namespace dai; namespace po = boost::program_options; +typedef boost::numeric::ublas::compressed_matrix matrix; +typedef matrix::value_array_type::const_iterator matrix_vcit; +typedef matrix::index_array_type::const_iterator matrix_icit; + + +Factor BinaryFactor( const Var &n, double field ) { + assert( n.states() == 2 ); + double buf[2]; + buf[0] = exp(-field); + buf[1] = exp(field); + return Factor(n, &buf[0]); +} + + +Factor BinaryFactor( const Var &n1, const Var &n2, double coupling ) { + assert( n1.states() == 2 ); + assert( n2.states() == 2 ); + assert( n1 != n2 ); + double buf[4]; + buf[0] = (buf[3] = exp(coupling)); + buf[1] = (buf[2] = exp(-coupling)); + return Factor(n1 | n2, &buf[0]); +} + + +Factor RandomFactor( const VarSet &ns, double beta ) { + Factor fac( ns ); + for( size_t t = 0; t < fac.states(); t++ ) + fac[t] = exp(rnd_stdnormal() * beta); + return fac; +} + + +Factor PottsFactor( const Var &n1, const Var &n2, double beta ) { + Factor fac( n1 | n2, 1.0 ); + assert( n1.states() == n2.states() ); + for( size_t s = 0; s < n1.states(); s++ ) + fac[ s * (n1.states() + 1) ] = exp(beta); + return fac; +} void MakeHOIFG( size_t N, size_t M, size_t k, double sigma, FactorGraph &fg ) { vector vars; vector factors; + vars.reserve(N); for( size_t i = 0; i < N; i++ ) vars.push_back(Var(i,2)); @@ -45,286 +91,689 @@ void MakeHOIFG( size_t N, size_t M, size_t k, double sigma, FactorGraph &fg ) { do { size_t newind = (size_t)(N * rnd_uniform()); Var newvar = Var(newind, 2); - if( !(vars.contains( newvar )) ) { + if( !vars.contains( newvar ) ) { vars |= newvar; break; } } while( 1 ); } - Factor newfac(vars); - for( size_t t = 0; t < newfac.states(); t++ ) - newfac[t] = exp(rnd_stdnormal() * sigma); - factors.push_back(newfac); + factors.push_back( RandomFactor( vars, sigma ) ); } - fg = FactorGraph(factors); -}; + fg = FactorGraph( factors.begin(), factors.end(), vars.begin(), vars.end(), factors.size(), vars.size() ); +} -void MakeFullFG( size_t N, double sigma_w, double sigma_th, string type, FactorGraph &fg ) { - vector vars; +// w should be upper triangular or lower triangular +void WTh2FG( const matrix &w, const vector &th, FactorGraph &fg ) { + vector vars; vector factors; - double w[N][N]; - double th[N]; - double buf[4]; + size_t N = th.size(); + assert( (w.size1() == N) && (w.size2() == N) ); + vars.reserve(N); for( size_t i = 0; i < N; i++ ) vars.push_back(Var(i,2)); - - for( size_t i = 0; i < N; i++ ) - for( size_t j = 0; j < N; j++ ) - w[i][j] = 0.0; + factors.reserve( w.nnz() + N ); + // walk through the sparse array structure + // this is similar to matlab sparse arrays + // index2 gives the column index (similar to ir in matlab) + // index1 gives the starting indices for each row (similar to jc in matlab) + size_t i = 0; + for( size_t pos = 0; pos < w.nnz(); pos++ ) { + while( pos == w.index1_data()[i+1] ) + i++; + size_t j = w.index2_data()[pos]; + double w_ij = w.value_data()[pos]; + factors.push_back( BinaryFactor( vars[i], vars[j], w_ij ) ); + } for( size_t i = 0; i < N; i++ ) - for( size_t j = i+1; j < N; j++ ) { - w[i][j] = rnd_stdnormal() * sigma_w; - if( type == "fe" ) - w[i][j] = fabs(w[i][j]); - else if( type == "af" ) - w[i][j] = -fabs(w[i][j]); - w[j][i] = w[i][j]; - buf[0] = (buf[3] = exp(w[i][j])); - buf[1] = (buf[2] = exp(-w[i][j])); - factors.push_back(Factor(VarSet(vars[i],vars[j]),buf)); - } + factors.push_back( BinaryFactor( vars[i], th[i] ) ); + + fg = FactorGraph( factors.begin(), factors.end(), vars.begin(), vars.end(), factors.size(), vars.size() ); +} + + +void MakeFullFG( size_t N, double mean_w, double mean_th, double sigma_w, double sigma_th, FactorGraph &fg ) { + matrix w(N,N,N*(N-1)/2); + vector th(N,0.0); for( size_t i = 0; i < N; i++ ) { - th[i] = rnd_stdnormal() * sigma_th; - buf[0] = exp(th[i]); - buf[1] = exp(-th[i]); - factors.push_back(Factor(vars[i],buf)); + for( size_t j = i+1; j < N; j++ ) + w(i,j) = rnd_stdnormal() * sigma_w + mean_w; + th[i] = rnd_stdnormal() * sigma_th + mean_th; } - fg = FactorGraph(factors); -}; + WTh2FG( w, th, fg ); +} -void MakeGridFG( long periodic, long n, double sigma_w, double sigma_th, string type, FactorGraph &fg ) { +void Make3DPotts( size_t n1, size_t n2, size_t n3, size_t states, double beta, FactorGraph &fg ) { vector vars; vector factors; + + for( size_t i1 = 0; i1 < n1; i1++ ) + for( size_t i2 = 0; i2 < n2; i2++ ) + for( size_t i3 = 0; i3 < n3; i3++ ) { + vars.push_back( Var( i1*n2*n3 + i2*n3 + i3, states ) ); + if( i1 ) + factors.push_back( PottsFactor( vars.back(), vars[ (i1-1)*n2*n3 + i2*n3 + i3 ], beta ) ); + if( i2 ) + factors.push_back( PottsFactor( vars.back(), vars[ i1*n2*n3 + (i2-1)*n3 + i3 ], beta ) ); + if( i3 ) + factors.push_back( PottsFactor( vars.back(), vars[ i1*n2*n3 + i2*n3 + (i3-1) ], beta ) ); + } - long N = n*n; + fg = FactorGraph( factors.begin(), factors.end(), vars.begin(), vars.end(), factors.size(), vars.size() ); +} - double w[N][N]; - double th[N]; - double buf[4]; - for( long i = 0; i < N; i++ ) - vars.push_back(Var(i,2)); - - for( long i = 0; i < N; i++ ) - for( long j = 0; j < N; j++ ) - w[i][j] = 0.0; +void MakeGridFG( long periodic, size_t n, double mean_w, double mean_th, double sigma_w, double sigma_th, FactorGraph &fg ) { + size_t N = n*n; - for( long i = 0; i < n; i++ ) - for( long j = 0; j < n; j++ ) { + matrix w(N,N,2*N); + vector th(N,0.0); + + for( size_t i = 0; i < n; i++ ) + for( size_t j = 0; j < n; j++ ) { if( i+1 < n || periodic ) - w[i*n+j][((i+1)%n)*n+j] = 1.0; - if( i > 0 || periodic ) - w[i*n+j][((i+n-1)%n)*n+j] = 1.0; + w(i*n+j, ((i+1)%n)*n+j) = rnd_stdnormal() * sigma_w + mean_w; if( j+1 < n || periodic ) - w[i*n+j][i*n+((j+1)%n)] = 1.0; - if( j > 0 || periodic ) - w[i*n+j][i*n+((j+n-1)%n)] = 1.0; + w(i*n+j, i*n+((j+1)%n)) = rnd_stdnormal() * sigma_w + mean_w; + th[i*n+j] = rnd_stdnormal() * sigma_th + mean_th; } - for( long i = 0; i < N; i++ ) - for( long j = i+1; j < N; j++ ) - if( w[i][j] ) { - w[i][j] = rnd_stdnormal() * sigma_w; - if( type == "fe" ) - w[i][j] = fabs(w[i][j]); - else if( type == "af" ) - w[i][j] = -fabs(w[i][j]); - w[j][i] = w[i][j]; - buf[0] = (buf[3] = exp(w[i][j])); - buf[1] = (buf[2] = exp(-w[i][j])); - factors.push_back(Factor(VarSet(vars[i],vars[j]),buf)); - } + WTh2FG( w, th, fg ); +} + + +void MakeGridNonbinaryFG( bool periodic, size_t n, size_t states, double beta, FactorGraph &fg ) { + size_t N = n*n; + + vector vars; + vector factors; + + vars.reserve(N); + for( size_t i = 0; i < N; i++ ) + vars.push_back(Var(i, states)); + + factors.reserve( 2 * N ); + for( size_t i = 0; i < n; i++ ) { + for( size_t j = 0; j < n; j++ ) { + if( i+1 < n || periodic ) + factors.push_back( RandomFactor( vars[i*n+j] | vars[((i+1)%n)*n+j], beta ) ); + if( j+1 < n || periodic ) + factors.push_back( RandomFactor( vars[i*n+j] | vars[i*n+((j+1)%n)], beta ) ); + } + } + + fg = FactorGraph( factors.begin(), factors.end(), vars.begin(), vars.end(), factors.size(), vars.size() ); +} - for( long i = 0; i < N; i++ ) { - th[i] = rnd_stdnormal() * sigma_th; - buf[0] = exp(th[i]); - buf[1] = exp(-th[i]); - factors.push_back(Factor(vars[i],buf)); + +void MakeLoopFG( size_t N, double mean_w, double mean_th, double sigma_w, double sigma_th, FactorGraph &fg ) { + matrix w(N,N,N); + vector th(N,0.0); + + for( size_t i = 0; i < N; i++ ) { + w(i, (i+1)%N) = rnd_stdnormal() * sigma_w + mean_w; + th[i] = rnd_stdnormal() * sigma_th + mean_th; } - fg = FactorGraph(factors); -}; + WTh2FG( w, th, fg ); +} -void MakeDRegFG( size_t N, size_t d, double sigma_w, double sigma_th, string type, FactorGraph &fg ) { - vector vars; +void MakeLoopNonbinaryFG( size_t N, size_t states, double beta, FactorGraph &fg ) { + vector vars; vector factors; - double w[N][N]; - double th[N]; - double buf[4]; - - for( size_t i = 0; i < N; i++ ) - vars.push_back(Var(i,2)); - + vars.reserve(N); for( size_t i = 0; i < N; i++ ) - for( size_t j = 0; j < N; j++ ) - w[i][j] = 0.0; + vars.push_back(Var(i, states)); - UEdgeVec g = RandomDRegularGraph( N, d ); - for( size_t i = 0; i < g.size(); i++ ) { - w[g[i].n1][g[i].n2] = 1.0; - w[g[i].n2][g[i].n1] = 1.0; + factors.reserve( N ); + for( size_t i = 0; i < N; i++ ) { + factors.push_back( RandomFactor( vars[i] | vars[(i+1)%N], beta ) ); } - - for( size_t i = 0; i < N; i++ ) - for( size_t j = i+1; j < N; j++ ) - if( w[i][j] ) { - w[i][j] = rnd_stdnormal() * sigma_w; - if( type == "fe" ) - w[i][j] = fabs(w[i][j]); - else if( type == "af" ) - w[i][j] = -fabs(w[i][j]); - w[j][i] = w[i][j]; - buf[0] = (buf[3] = exp(w[i][j])); - buf[1] = (buf[2] = exp(-w[i][j])); - factors.push_back(Factor(VarSet(vars[i],vars[j]),buf)); - } + + fg = FactorGraph( factors.begin(), factors.end(), vars.begin(), vars.end(), factors.size(), vars.size() ); +} + + +void MakeTreeFG( size_t N, double mean_w, double mean_th, double sigma_w, double sigma_th, FactorGraph &fg ) { + matrix w(N,N,N-1); + vector th(N,0.0); for( size_t i = 0; i < N; i++ ) { - th[i] = rnd_stdnormal() * sigma_th; - buf[0] = exp(th[i]); - buf[1] = exp(-th[i]); - factors.push_back(Factor(vars[i],buf)); + th[i] = rnd_stdnormal() * sigma_th + mean_th; + if( i > 0 ) { + size_t j = rnd_int( 0, i-1 ); + w(i,j) = rnd_stdnormal() * sigma_w + mean_w; + } } - fg = FactorGraph(factors); -}; + WTh2FG( w, th, fg ); +} + + +void MakeDRegFG( size_t N, size_t d, double mean_w, double mean_th, double sigma_w, double sigma_th, FactorGraph &fg ) { + matrix w(N,N,(d*N)/2); + vector th(N,0.0); + + UEdgeVec g = RandomDRegularGraph( N, d ); + for( size_t i = 0; i < g.size(); i++ ) + w(g[i].n1, g[i].n2) = rnd_stdnormal() * sigma_w + mean_w; + + for( size_t i = 0; i < N; i++ ) + th[i] = rnd_stdnormal() * sigma_th + mean_th; + + WTh2FG( w, th, fg ); +} + +// N = number of variables +// n = size of variable neighborhoods +// K = number of factors +// k = size of factor neighborhoods +// asserts: N * n == K * k +BipartiteGraph CreateRandomBipartiteGraph( size_t N, size_t K, size_t n, size_t k ) { + BipartiteGraph G; -const char *HOITYPE = "hoi"; -const char *FULLTYPE = "full"; -const char *GRIDTYPE = "grid"; -const char *DREGTYPE = "dreg"; + assert( N * n == K * k ); + // build lists of degree-repeated vertex numbers + std::vector stubs1(N*n,0); + for( size_t i = 0; i < N; i++ ) + for( size_t t = 0; t < n; t++ ) + stubs1[i*n + t] = i; + + // build lists of degree-repeated vertex numbers + std::vector stubs2(K*k,0); + for( size_t I = 0; I < K; I++ ) + for( size_t t = 0; t < k; t++ ) + stubs2[I*k + t] = I; + + // shuffle lists + random_shuffle( stubs1.begin(), stubs1.end() ); + random_shuffle( stubs2.begin(), stubs2.end() ); -// Old usages: -// create_full_fg -// create_grid_fg -// create_dreg_fg + // add edges + vector edges; + edges.reserve( N*n ); + for( size_t e = 0; e < N*n; e++ ) + edges.push_back( BipartiteGraph::Edge(stubs1[e], stubs2[e]) ); + + // finish construction + G.create( N, K, edges.begin(), edges.end() ); + + return G; +} + + +// Returns x**n % p, assuming p is prime +int powmod (int x, int n, int p) { + int y = 1; + for( int m = 0; m < n; m++ ) + y = (x * y) % p; + return y; +} + + +// Returns order of x in GF(p) with p prime +size_t order (int x, int p) { + x = x % p; + assert( x != 0 ); + size_t n = 0; + size_t prod = 1; + do { + prod = (prod * x) % p; + n++; + } while( prod != 1 ); + return n; +} + + +// Returns whether n is a prime number +bool isPrime (size_t n) { + bool result = true; + for( size_t k = 2; (k < n) && result; k++ ) + if( n % k == 0 ) + result = false; + return result; +} + + +// Make a regular LDPC graph with N=6, j=2, K=4, k=3 +BipartiteGraph CreateSmallLDPCGraph() { + BipartiteGraph G; + size_t N=4, j=3, K=4; // k=3; + + typedef BipartiteGraph::Edge Edge; + vector edges; + edges.reserve( N*j ); + edges.push_back( Edge(0,0) ); edges.push_back( Edge(1,0) ); edges.push_back( Edge(2,0) ); + edges.push_back( Edge(0,1) ); edges.push_back( Edge(1,1) ); edges.push_back( Edge(3,1) ); + edges.push_back( Edge(0,2) ); edges.push_back( Edge(2,2) ); edges.push_back( Edge(3,2) ); + edges.push_back( Edge(1,3) ); edges.push_back( Edge(2,3) ); edges.push_back( Edge(3,3) ); + + // finish construction + G.create( N, K, edges.begin(), edges.end() ); + + return G; +} + + +// Use construction described in "A Class of Group-Structured LDPC Codes" +// by R. M. Tanner, D. Sridhara and T. Fuja +// Proceedings of ICSTA, 2001 +// +// Example parameters: (p,j,k) = (31,3,5) +// j and k must be divisors of p-1 +BipartiteGraph CreateGroupStructuredLDPCGraph( size_t p, size_t j, size_t k ) { + BipartiteGraph G; + + size_t n = j; + size_t N = p * k; + size_t K = p * j; + + size_t a, b; + for( a = 2; a < p; a++ ) + if( order(a,p) == k ) + break; + assert( a != p ); + for( b = 2; b < p; b++ ) + if( order(b,p) == j ) + break; + assert( b != p ); +// cout << "# order(a=" << a << ") = " << order(a,p) << endl; +// cout << "# order(b=" << b << ") = " << order(b,p) << endl; + + assert( N * n == K * k ); + + typedef BipartiteGraph::Edge Edge; + vector edges; + edges.reserve( N * n ); + + for( size_t s = 0; s < j; s++ ) + for( size_t t = 0; t < k; t++ ) { + size_t P = (powmod(b,s,p) * powmod(a,t,p)) % p; + for( size_t m = 0; m < p; m++ ) + edges.push_back( Edge(t*p + m, s*p + ((m + P) % p)) ); + } + + // finish construction + G.create( N, K, edges.begin(), edges.end() ); + + return G; +} + + +// Make parity check table +void MakeParityCheck( double *result, size_t n, double eps ) { + size_t N = 1 << n; + for( size_t i = 0; i < N; i++ ) { + size_t c = 0; + for( size_t t = 0; t < n; t++ ) + if( i & (1 << t) ) + c ^= 1; + if( c ) + result[i] = eps; + else + result[i] = 1.0 - eps; + } + return; +} + + +const char *FULL_TYPE = "full"; +const char *GRID_TYPE = "grid"; +const char *GRID_TORUS_TYPE = "grid_torus"; +const char *DREG_TYPE = "dreg"; +const char *HOI_TYPE = "hoi"; +const char *LDPC_RANDOM_TYPE = "ldpc_random"; +const char *LDPC_GROUP_TYPE = "ldpc_group"; +const char *LDPC_SMALL_TYPE = "ldpc_small"; +const char *LOOP_TYPE = "loop"; +const char *TREE_TYPE = "tree"; +const char *POTTS3D_TYPE = "potts3d"; int main( int argc, char *argv[] ) { try { - size_t N, M, k, d; - size_t periodic; + size_t N, K, k, d, j, n1, n2, n3; + size_t prime; size_t seed; - double beta, sigma_w, sigma_th; - string type, subtype; + double beta, sigma_w, sigma_th, noise, mean_w, mean_th; + string type; + size_t states = 2; // Declare the supported options. po::options_description desc("Allowed options"); desc.add_options() ("help", "produce help message") - ("type", po::value(&type), "factor graph type:\n\t'full', 'grid', 'dreg' or 'hoi'") - ("seed", po::value(&seed), "random number seed") - ("subtype", po::value(&subtype), "interactions type:\n\t'sg', 'fe' or 'af'\n\t(ignored for type=='hoi')") - ("N", po::value(&N), "number of (binary) variables") - ("M", po::value(&M), "number of factors\n\t(only for type=='hoi')") - ("k", po::value(&k), "connectivity of the factors\n\t(only for type=='hoi')") + ("type", po::value(&type), "factor graph type:\n\t'full', 'grid', 'grid_torus', 'dreg', 'loop', 'tree', 'hoi', 'ldpc_random', 'ldpc_group', 'ldpc_small', 'potts3d'") + ("seed", po::value(&seed), "random number seed (tries to read from /dev/urandom if not specified)") + ("N", po::value(&N), "number of variables (not for type=='ldpc_small')") + ("n1", po::value(&n1), "width of 3D grid (only for type=='potts3d')") + ("n2", po::value(&n2), "height of 3D grid (only for type=='potts3d')") + ("n3", po::value(&n3), "length of 3D grid (only for type=='potts3d')") + ("K", po::value(&K), "number of factors\n\t(only for type=='hoi' and 'type=='ldpc_{random,group}')") + ("k", po::value(&k), "number of variables per factor\n\t(only for type=='hoi' and type=='ldpc_{random,group}')") ("d", po::value(&d), "variable connectivity\n\t(only for type=='dreg')") - ("beta", po::value(&beta), "stddev of log-factor entries\n\t(only for type=='hoi')") - ("sigma_w", po::value(&sigma_w), "stddev of pairwise interactions w_{ij}\n\t(ignored for type=='hoi')") - ("sigma_th", po::value(&sigma_th), "stddev of singleton interactions th_i\n\t(ignored for type=='hoi')") - ("periodic", po::value(&periodic), "0/1 corresponding to nonperiodic/periodic grid\n\t(only for type=='grid')") + ("j", po::value(&j), "number of parity checks per bit\n\t(only for type=='ldpc_{random,group}')") + ("prime", po::value(&prime), "prime number for construction of LDPC code\n\t(only for type=='ldpc_group')") + ("beta", po::value(&beta), "stddev of log-factor entries\n\t(only for type=='hoi', 'potts3d', 'grid' if states>2)") + ("mean_w", po::value(&mean_w), "mean of pairwise interactions w_{ij}\n\t(not for type=='hoi', 'ldpc_*', 'potts3d')") + ("mean_th", po::value(&mean_th), "mean of singleton interactions th_i\n\t(not for type=='hoi', 'ldpc_*', 'potts3d')") + ("sigma_w", po::value(&sigma_w), "stddev of pairwise interactions w_{ij}\n\t(not for type=='hoi', 'ldpc_*', 'potts3d')") + ("sigma_th", po::value(&sigma_th), "stddev of singleton interactions th_i\n\t(not for type=='hoi', 'ldpc_*', 'potts3d'") + ("noise", po::value(&noise), "bitflip probability for binary symmetric channel (only for type=='ldpc')") + ("states", po::value(&states), "number of states of each variable (should be 2 for all but type=='grid', 'grid_torus', 'loop', 'potts3d')") ; - po::variables_map vm; + po::variables_map vm; po::store(po::parse_command_line(argc, argv, desc), vm); - po::notify(vm); + po::notify(vm); if( vm.count("help") || !vm.count("type") ) { if( vm.count("type") ) { - if( type == HOITYPE ) { + if( type == FULL_TYPE ) { + cout << "Creates fully connected pairwise graphical model of binary variables;" << endl; + } else if( type == GRID_TYPE ) { + cout << "Creates (non-periodic) 2D Ising grid of (approx.) variables (which need not be binary);" << endl; + } else if( type == GRID_TORUS_TYPE ) { + cout << "Creates periodic 2D Ising grid of (approx.) variables (which need not be binary);" << endl; + } else if( type == DREG_TYPE ) { + cout << "Creates random d-regular graph of binary variables with uniform degree " << endl; + cout << "(where should be even);" << endl; + } else if( type == LOOP_TYPE ) { + cout << "Creates a pairwise graphical model consisting of a single loop of" << endl; + cout << " variables (which need not be binary);" << endl; + } else if( type == TREE_TYPE ) { + cout << "Creates a pairwise, connected graphical model without cycles (i.e., a tree)" << endl; + cout << "of binary variables;" << endl; + } else if( type == HOI_TYPE ) { cout << "Creates a random factor graph of binary variables and" << endl; - cout << " factors, each factor being an interaction of variables." << endl; + cout << " factors, each factor being an interaction of variables." << endl; cout << "The entries of the factors are exponentials of i.i.d. Gaussian" << endl; cout << "variables with mean 0 and standard deviation ." << endl; - } else if( type == FULLTYPE ) { - cout << "Creates fully connected pairwise graphical model of variables;" << endl; - } else if( type == GRIDTYPE ) { - cout << "Creates 2D Ising grid (periodic if !=0) of (approx.) variables;" << endl; - } else if( type == DREGTYPE ) { - cout << "Creates random d-regular graph of nodes with uniform degree " << endl; - cout << "(where should be even);" << endl; + } else if( type == LDPC_RANDOM_TYPE ) { + cout << "Simulates LDPC decoding problem, using a LDPC code of bits and parity" << endl; + cout << "checks, with bits per check and checks per bit, transmitted on a binary" << endl; + cout << "symmetric channel with probability of flipping a bit. The transmitted" << endl; + cout << "codeword has all bits set to zero. The LDPC code is randomly generated." << endl; + } else if( type == LDPC_GROUP_TYPE ) { + cout << "Simulates LDPC decoding problem, using a LDPC code of bits and parity" << endl; + cout << "checks, with bits per check and checks per bit, transmitted on a binary" << endl; + cout << "symmetric channel with probability of flipping a bit. The transmitted" << endl; + cout << "codeword has all bits set to zero. The LDPC code is constructed (using group" << endl; + cout << "theory) using a parameter ; and should both be divisors of -1." << endl; + } else if( type == LDPC_SMALL_TYPE ) { + cout << "Simulates LDPC decoding problem, using a LDPC code of 4 bits and 4 parity" << endl; + cout << "checks, with 3 bits per check and 3 checks per bit, transmitted on a binary" << endl; + cout << "symmetric channel with probability of flipping a bit. The transmitted" << endl; + cout << "codeword has all bits set to zero. The LDPC code is fixed." << endl; + } else if( type == POTTS3D_TYPE ) { + cout << "Builds 3D Potts model of size xx with nearest-neighbour Potts" << endl; + cout << "interactions with states and inverse temperature ." << endl; } else - cerr << "Unknown type (should be one of 'full', 'grid', 'dreg' or 'hoi')" << endl; + cerr << "Unknown type (should be one of 'full', 'grid', 'grid_torus', 'dreg', 'loop', 'tree', 'hoi', 'ldpc_random', 'ldpc_group', 'ldpc_small', 'potts3d')" << endl; - if( type == FULLTYPE || type == GRIDTYPE || type == DREGTYPE ) { - cout << "singleton interactions are Gaussian with mean 0 and standard" << endl; - cout << "deviation ; pairwise interactions are Gaussian with mean 0" << endl; - cout << "and standard deviation if =='sg', absolute value" << endl; - cout << "is taken if =='fe' and a minus sign is added if =='af'." << endl; - } + if( type == FULL_TYPE || type == GRID_TYPE || type == GRID_TORUS_TYPE || type == DREG_TYPE || type == LOOP_TYPE || type == TREE_TYPE ) { + if( type == GRID_TYPE || type == GRID_TORUS_TYPE || type == LOOP_TYPE ) { + cout << "if > 2: factor entries are exponents of Gaussians with mean 0 and standard deviation beta; otherwise," << endl; + } + cout << "singleton interactions are Gaussian with mean and standard" << endl; + cout << "deviation ; pairwise interactions are Gaussian with mean" << endl; + cout << " and standard deviation ." << endl; + } } cout << endl << desc << endl; return 1; } - if( !vm.count("seed") ) - throw "Please specify random number seed."; + if( !vm.count("states") ) + states = 2; + + if( !vm.count("seed") ) { + ifstream infile; + bool success; + infile.open( "/dev/urandom" ); + success = infile.is_open(); + if( success ) { + infile.read( (char *)&seed, sizeof(size_t) / sizeof(char) ); + success = infile.good(); + infile.close(); + } + if( !success ) + throw "Please specify random number seed."; + } rnd_seed( seed ); -// srand( gsl_rng_default_seed ); FactorGraph fg; cout << "# Factor graph made by " << argv[0] << endl; cout << "# type = " << type << endl; - if( type == HOITYPE ) { - if( !vm.count("N") || !vm.count("M") || !vm.count("k") || !vm.count("beta") ) - throw "Please specify all required arguments"; - do { - MakeHOIFG( N, M, k, beta, fg ); - } while( !fg.G.isConnected() ); - - cout << "# N = " << N << endl; - cout << "# M = " << M << endl; - cout << "# k = " << k << endl; - cout << "# beta = " << beta << endl; - } else if( type == FULLTYPE ) { - if( !vm.count("N") || !vm.count("sigma_w") || !vm.count("sigma_th") || !vm.count("subtype") ) + if( type == FULL_TYPE ) { + if( !vm.count("N") || !vm.count("mean_w") || !vm.count("mean_th") || !vm.count("sigma_w") || !vm.count("sigma_th") ) throw "Please specify all required arguments"; - MakeFullFG( N, sigma_w, sigma_th, subtype, fg ); + MakeFullFG( N, mean_w, mean_th, sigma_w, sigma_th, fg ); cout << "# N = " << N << endl; + cout << "# mean_w = " << mean_w << endl; + cout << "# mean_th = " << mean_th << endl; cout << "# sigma_w = " << sigma_w << endl; cout << "# sigma_th = " << sigma_th << endl; - cout << "# subtype = " << subtype << endl; - } else if( type == GRIDTYPE ) { - if( !vm.count("N") || !vm.count("sigma_w") || !vm.count("sigma_th") || !vm.count("subtype") || !vm.count("periodic") ) - throw "Please specify all required arguments"; + } else if( type == GRID_TYPE || type == GRID_TORUS_TYPE ) { + if( states > 2 ) { + if( !vm.count("N") || !vm.count("beta") ) + throw "Please specify all required arguments"; + } else { + if( !vm.count("N") || !vm.count("mean_w") || !vm.count("mean_th") || !vm.count("sigma_w") || !vm.count("sigma_th") ) + throw "Please specify all required arguments"; + } size_t n = (size_t)sqrt((long double)N); N = n * n; - MakeGridFG( periodic, n, sigma_w, sigma_th, subtype, fg ); - - cout << "# periodic = " << periodic << endl; + bool periodic = false; + if( type == GRID_TYPE ) + periodic = false; + else + periodic = true; + + if( states > 2 ) + MakeGridNonbinaryFG( periodic, n, states, beta, fg ); + else + MakeGridFG( periodic, n, mean_w, mean_th, sigma_w, sigma_th, fg ); + cout << "# n = " << n << endl; cout << "# N = " << N << endl; - cout << "# sigma_w = " << sigma_w << endl; - cout << "# sigma_th = " << sigma_th << endl; - cout << "# subtype = " << subtype << endl; - } else if( type == DREGTYPE ) { - if( !vm.count("N") || !vm.count("sigma_w") || !vm.count("sigma_th") || !vm.count("subtype") || !vm.count("d") ) + + if( states > 2 ) + cout << "# beta = " << beta << endl; + else { + cout << "# mean_w = " << mean_w << endl; + cout << "# mean_th = " << mean_th << endl; + cout << "# sigma_w = " << sigma_w << endl; + cout << "# sigma_th = " << sigma_th << endl; + } + } else if( type == DREG_TYPE ) { + if( !vm.count("N") || !vm.count("mean_w") || !vm.count("mean_th") || !vm.count("sigma_w") || !vm.count("sigma_th") || !vm.count("d") ) throw "Please specify all required arguments"; - MakeDRegFG( N, d, sigma_w, sigma_th, subtype, fg ); + MakeDRegFG( N, d, mean_w, mean_th, sigma_w, sigma_th, fg ); cout << "# N = " << N << endl; cout << "# d = " << d << endl; + cout << "# mean_w = " << mean_w << endl; + cout << "# mean_th = " << mean_th << endl; + cout << "# sigma_w = " << sigma_w << endl; + cout << "# sigma_th = " << sigma_th << endl; + } else if( type == LOOP_TYPE ) { + if( states > 2 ) { + if( !vm.count("N") || !vm.count("beta") ) + throw "Please specify all required arguments"; + } else { + if( !vm.count("N") || !vm.count("mean_w") || !vm.count("mean_th") || !vm.count("sigma_w") || !vm.count("sigma_th") ) + throw "Please specify all required arguments"; + } + if( states > 2 ) + MakeLoopNonbinaryFG( N, states, beta, fg ); + else + MakeLoopFG( N, mean_w, mean_th, sigma_w, sigma_th, fg ); + + cout << "# N = " << N << endl; + + if( states > 2 ) + cout << "# beta = " << beta << endl; + else { + cout << "# mean_w = " << mean_w << endl; + cout << "# mean_th = " << mean_th << endl; + cout << "# sigma_w = " << sigma_w << endl; + cout << "# sigma_th = " << sigma_th << endl; + } + } else if( type == TREE_TYPE ) { + if( !vm.count("N") || !vm.count("mean_w") || !vm.count("mean_th") || !vm.count("sigma_w") || !vm.count("sigma_th") ) + throw "Please specify all required arguments"; + MakeTreeFG( N, mean_w, mean_th, sigma_w, sigma_th, fg ); + + cout << "# N = " << N << endl; + cout << "# mean_w = " << mean_w << endl; + cout << "# mean_th = " << mean_th << endl; cout << "# sigma_w = " << sigma_w << endl; cout << "# sigma_th = " << sigma_th << endl; - cout << "# subtype = " << subtype << endl; + } else if( type == HOI_TYPE ) { + if( !vm.count("N") || !vm.count("K") || !vm.count("k") || !vm.count("beta") ) + throw "Please specify all required arguments"; + do { + MakeHOIFG( N, K, k, beta, fg ); + } while( !fg.G.isConnected() ); + + cout << "# N = " << N << endl; + cout << "# K = " << K << endl; + cout << "# k = " << k << endl; + cout << "# beta = " << beta << endl; + } else if( type == LDPC_RANDOM_TYPE || type == LDPC_GROUP_TYPE || type == LDPC_SMALL_TYPE ) { + if( !vm.count("noise") ) + throw "Please specify all required arguments"; + + if( type == LDPC_RANDOM_TYPE ) { + if( !vm.count("N") || !vm.count("K") || !vm.count("j") || !vm.count("k") ) + throw "Please specify all required arguments"; + + if( N * j != K * k ) + throw "Parameters should satisfy N * j == K * k"; + } else if( type == LDPC_GROUP_TYPE ) { + if( !vm.count("prime") || !vm.count("j") || !vm.count("k") ) + throw "Please specify all required arguments"; + + if( !isPrime(prime) ) + throw "Parameter should be prime"; + if( !((prime-1) % j == 0 ) ) + throw "Parameters should satisfy (prime-1) % j == 0"; + if( !((prime-1) % k == 0 ) ) + throw "Parameters should satisfy (prime-1) % k == 0"; + + N = prime * k; + K = prime * j; + } else if( type == LDPC_SMALL_TYPE ) { + N = 4; + K = 4; + j = 3; + k = 3; + } + + cout << "# N = " << N << endl; + cout << "# K = " << K << endl; + cout << "# j = " << j << endl; + cout << "# k = " << k << endl; + if( type == LDPC_GROUP_TYPE ) + cout << "# prime = " << prime << endl; + cout << "# noise = " << noise << endl; + + // p = 31, j = 3, k = 5 + // p = 37, j = 3, k = 4 + // p = 7 , j = 2, k = 3 + // p = 29, j = 2, k = 4 + + // Construct likelihood and paritycheck factors + double likelihood[4] = {1.0 - noise, noise, noise, 1.0 - noise}; + double *paritycheck = new double[1 << k]; + MakeParityCheck(paritycheck, k, 0.0); + + // Create LDPC structure + BipartiteGraph ldpcG; + bool regular; + do { + if( type == LDPC_GROUP_TYPE ) + ldpcG = CreateGroupStructuredLDPCGraph( prime, j, k ); + else if( type == LDPC_RANDOM_TYPE ) + ldpcG = CreateRandomBipartiteGraph( N, K, j, k ); + else if( type == LDPC_SMALL_TYPE ) + ldpcG = CreateSmallLDPCGraph(); + + regular = true; + for( size_t i = 0; i < N; i++ ) + if( ldpcG.nb1(i).size() != j ) + regular = false; + for( size_t I = 0; I < K; I++ ) + if( ldpcG.nb2(I).size() != k ) + regular = false; + } while( !regular && !ldpcG.isConnected() ); + + // Convert to FactorGraph + vector factors; + for( size_t I = 0; I < K; I++ ) { + VarSet vs; + for( size_t _i = 0; _i < k; _i++ ) { + size_t i = ldpcG.nb2(I)[_i]; + vs |= Var( i, 2 ); + } + factors.push_back( Factor( vs, paritycheck ) ); + } + delete paritycheck; + + // Generate noise vector + vector noisebits(N,0); + size_t bitflips = 0; + for( size_t i = 0; i < N; i++ ) { + if( rnd_uniform() < noise ) { + noisebits[i] = 1; + bitflips++; + } + } + cout << "# bitflips = " << bitflips << endl; + + // Simulate transmission of all-zero codeword + vector input(N,0); + vector output(N,0); + for( size_t i = 0; i < N; i++ ) + output[i] = (input[i] + noisebits[i]) & 1; + + // Add likelihoods + for( size_t i = 0; i < N; i++ ) + factors.push_back( Factor(Var(i,2), likelihood + output[i]*2) ); + + // Construct Factor Graph + fg = FactorGraph( factors ); + } else if( type == POTTS3D_TYPE ) { + if( !vm.count("n1") || !vm.count("n2") || !vm.count("n3") || !vm.count("beta") || !vm.count("states") ) + throw "Please specify all required arguments"; + Make3DPotts( n1, n2, n3, states, beta, fg ); + + cout << "# N = " << n1*n2*n3 << endl; + cout << "# n1 = " << n1 << endl; + cout << "# n2 = " << n2 << endl; + cout << "# n3 = " << n3 << endl; + cout << "# beta = " << beta << endl; + cout << "# states = " << states << endl; + } else { + throw "Invalid type"; } cout << "# seed = " << seed << endl; -- 2.20.1