formlp.py

著名的ldpc编解码的资料及源代码

"""
This file implements the ECCLPFormer class and QuantLPFormer to
decode a low density parity check (LDPC) error correcting code or
quantize via the dual of an LDPC code using a linear programming
relaxation.  The basic idea is that you instantiate an ECCLPFormer or
QuantLPFormer and give it the code parameters.  Then you tell it to
form an LP and either solve it or print it out.

Before using this file you need to download and install the freely
available GNU Linear Programming Kit (GLPK).  If you install the
executable glpsol from GLPK in a weird place that is not in the path
seen by python, set the default for the LPSolver variable in the
__init__ method of the LPFormer base class to point to the glpsol
executable.

Some examples of how to use the ECCLPFormer and QuantLPFormer are shown
below.  While this file can stand alone (provided you have GLPK), it
is most useful when used with the other pycodes routines developed
by Emin Martinian and available from http://freshmeat.net/projects/pycodes.

If you have these other files available you can have python test the
examples provided below by using the command 'python FormLP.py'.
These tests can range from less than a minute to a 3GHz Linux system
to around 5 minutes on an older G3 Mac running OS 10.2.

######################################################################
#                                                                    #
# Start error correction examples                                    #
#                                                                    #
######################################################################


# The following is a simple example of how to use the ECCLPFormer class
# for a parity check matrix representing the code shown below.
#
# y0  y1     y2  y3
#  =  =      =   =
#  \  |\    /|  /
#   \ | \  / | /
#    \|  \/  |/
#     +   +  +
#
# The only two codewords for this code are 0000 and 1111.
# First we generat the ECCLPFormer class, then we use it solve for
# the optimal y0,y1,y2,y3 given the received data [1,0,1,1].
# This received data corresponds to sending 1111 and getting an
# error on the second bit.  The LP decoder correcterly decoes
# to the answer y0,y1,y2,y3 = 1,1,1,1.
#
>>> from FormLP import *
>>> r = ECCLPFormer(4,1,[[0,1],[1,2],[2,3]])
>>> r.FormLP([1,0,1,1])
>>> (v,s,o) = r.SolveLP()
>>> print v
[1.0, 1.0, 1.0, 1.0]

# Next we do LP decoding for a medium size Gallager code assuming
# that the all zeros codeword was transmitted.  Feldman, Karger,
# and Wainwright argue that analyzing things assuming the all-0
# codeword was sent over a binary symmetric channel is valid provided
# the LP satisfies certain conditions (see their 2003 CISS paper
# for more details).  IMPORTANT: the all-0 assumption works for
# analyzing things sent over a BSC but *NOT* over an erasure channel.
# The following test takes about a minute to run on a Mac G3.
>>> N = 1200
>>> K = 600
>>> numErrors = 90 # error rate of 7.5%
>>> from FormLP import *
>>> from CodeMaker import *
>>> from random import *
>>> regL = make_H_gallager(N,3,6)
>>> origSource = [0]*N
>>> recSource = list(origSource)
>>> i = 0
>>> while (i < numErrors):
... 	index = randrange(N)
... 	if (0 == recSource[index]):
... 	    recSource[randrange(N)] = 1
... 	    i = i+1
... 
>>> r = ECCLPFormer(N,K,regL)
>>> r.FormLP(recSource)
>>> (v,s,o) = r.SolveLP()
>>> errors = map(lambda x,y: int(x) != int(y), origSource,v)
>>> print 'num errors = ',  errors.count(1)
num errors =  0

######################################################################
#                                                                    #
# Start quantization examples                                        #
#                                                                    #
######################################################################


# The following is a simple example of how to use the QuantLPFormer class
# for a generator matrix representing the code shown below.
#
# y0  y1     y2  y3
#  +  +      +   +
#  \  |\    /|  /
#   \ | \  / | /
#    \|  \/  |/
#    x0  x1  x2
#
# First we generat the QuantLPFormer class, then we use it solve for
# the optimal x0,x1,x2 when y0,y1,y2,y3=[1,0,0,1].  The answer turns 
# out to be x0,x1,x2 = 1,1,1.
#
>>> from FormLP import *
>>> r = QuantLPFormer(4,3,[[0],[0,1],[1,2],[2]])
>>> r.FormLP([1,0,0,1])
>>> (v,s,o) = r.SolveLP()
>>> print v
[1.0, 1.0, 1.0]

# In the following example we take the dual of a (7,4) Hamming code using the
# built in function TransposeCodeMatrix and then use that as the generator
# matrix for quantization.  In this example we quantize the
# sequence 0,*,*,*,*,*,1 where the *'s represent don't cares which can
# be reconstructed to either 0 or 1.
# 
>>> r = QuantLPFormer(7,3,TransposeCodeMatrix(7,4,[[0,1,2,4],[1,2,3,5],[2,3,4,6]]))
>>> r.FormLP([0,.5,.5,.5,.5,.5,1])
>>> (v,s,o) = r.SolveLP()
>>> print v
[0.0, 0.0, 1.0]


#
# Iteratively solve an quantization LP for a medium size code.  
# This does not seem to work all that well, but none of the other LP
# relaxations does much better at quantization either.
#
>>> N = 300
>>> K = 150
>>> numErase = 180
>>> numIter = 1000
>>> from FormLP import *
>>> from CodeMaker import *
>>> from random import *
>>> regL = make_H_gallager(N,3,6)
>>> source = map(lambda x: round(random()),range(N))
>>> for i in range(numErase):
... 	source[randrange(N)] = 0.5
... 
>>> r = QuantLPFormer(N,K,TransposeCodeMatrix(N,K,regL))
>>> r.FormLP(source)
>>> (v,s,o) = r.IterSolveLP(numIter,verbose=0)
>>> from encoders import *
>>> recon = EncodeFromLinkArray(map(lambda x: int(x),v),N,regL)
>>> diffs = map(lambda x,y: x != 0.5 and x != y, source,recon)
>>> print 'num flips = ',  diffs.count(1)
num flips =  25
>>> if (25 != diffs.count(1)):
... 	print 'failure may be due to diffs w/glpsol on different platforms'
>>>

# The following example illustrates what can go wrong with the LP
# relaxation in doing quantization.  Choosing v = [1.0,1.0,0.0] would
# reconstruct the source perfectly in all the unerased positions
# (places where the source is not 0.5).  But the LP relaxation produces
# the vector [1.0/3.0, 1.0/3.0, 1.0/3.0].  First, this 'solution' is
# not even binary, and second even rounding the bits would not give the
# right answer.

>>> from FormLP import *
>>> from CodeMaker import *
>>> from random import *
>>> hammingCode = [[0,1,2,4],[1,2,3,5],[2,3,4,6]]
>>> source = [0.5, 0.0, 0.0, 0.5, 1.0, 0.5, 0.5]
>>> r = QuantLPFormer(7,3,TransposeCodeMatrix(7,4,hammingCode))
>>> r.FormErasureLP(source)
>>> (v,s,o) = r.SolveLP()
>>> from encoders import *
>>> recon = EncodeFromLinkArray(map(lambda x: int(x),v),7,hammingCode)
>>> diffs = map(lambda x,y: x != 0.5 and x != y, source,recon)
>>> print 'num flips = ', diffs.count(1)
num flips =  1

"""

import string, tempfile, commands, re

def q(a,b,c):
    if a:
        return b
    else:
        return c

def CountOnes(numBitsToCheck,z):
    mask = 1
    count = 0
    for i in range(numBitsToCheck):
        if (mask & z):
            count = count + 1
        mask = mask << 1
    return count 

def OddParity(numBitsToCheck,z):
    return CountOnes(numBitsToCheck,z) % 2

def EvenParity(numBitsToCheck,z):
    return not OddParity(numBitsToCheck,z)

def TransposeCodeMatrix(N,K,codeMatrix):
    """
    Take a parity check code matrix with N columns and N-K rows and
    transposes it.

    The input codeMatrix is a list of N-K lists representing the rows
    of the parity check matrix.  Specifically, codeMatrix[i][j]=t
    means that the jth one in the ith row appears in column t.

    This function returns a list of N lists representing the transpose.
    So result[i][j]=t means that the jth one in the ith column appears
    in row t.

    This function is useful for taking a code produced by the file
    CodeMaker.py and putting it into a form suitable for the
    QuantLPFormer class.
    """
    assert N-K == len(codeMatrix)

    result = range(N)
    for i in range(N):    
        result[i] = []
    for row in range(N-K):
        for col in codeMatrix[row]:
            result[col].append(row)
    return result

class LPFormer:
    """
    The LPFormer class is a base class for the ECCLPFormer and
    QuantLPFormer classes.  You should use one of those and *not*
    use LPFormer class directly.
    """
    
    def __init__(self,N,K,codeMatrix,
                 LPSolver = 'glpsol',
                 LPSolveOpts = ' --nomip '):
        """
        N:          Number of outputs.
        K:          Number of inputs.
        codeMatrix: A representation of the generator matrix for the code
                    as a list of lists.  Specifically, codeMatrix[i] should
                    be a list of the variables which are summed modulo 2
                    to produce the ith output.
        LPSolver:   Path to the executable for the GNU Linear Programming Kit.
        LPSolveOpts:Options to GLPK.
        """
        self.N = N
        self.K = K
        self.codeMatrix = codeMatrix
        self.LPSolver = LPSolver
        self.LPSolveOpts = LPSolveOpts

        self.Reset()

    def Reset(self):
        self.output = []
        self.channelData = []


    def SolveLP(self):
        """
        Solve the linear program developed by calling the FormLP method
        (i.e., make sure you call FormLP before calling this).

        The return value is a 3-tuple (v,s,o) where v represents the
        (possibly fractional) values for the K variables, s represents
        the status returned by the LPSolver, and o is the text output
        of the LPSolver.
        """
        
        resultFile = tempfile.mktemp()
        modelFile = tempfile.mktemp()
        fd = open(modelFile,'w')
        fd.write(self.GetLP())
        fd.close()
        (s,o) = commands.getstatusoutput(self.LPSolver + self.LPSolveOpts +