genericmatrix.py

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

                        is not invertable.
        """

        workingCopy = self.Copy()
        result = self.MakeSimilarMatrix(self.Size(),'i')
        workingCopy.LowerGaussianElim(result)
        workingCopy.UpperInverse(result)
        return result

    def Determinant(self):
        """
        Function:       Determinant
        Description:    Returns the determinant of the matrix or raises
                        a ValueError if the matrix is not square.
        """
        if (self.rows != self.cols):
            raise ValueError, 'matrix not square'
        workingCopy = self.Copy()
        result = self.MakeSimilarMatrix(self.Size(),'i')
        workingCopy.LowerGaussianElim(result)
        det = self.identityElement
        for i in range(self.rows):
            det = det * workingCopy.data[i][i]
        return det

    def LUP(self):
        """
        Function:       (l,u,p) = self.LUP()
        Purpose:        Compute the LUP decomposition of self.
        Description:    This function returns three matrices
                        l, u, and p such that p * self = l * u
                        where l, u, and p have the following properties:

                        l is lower triangular with ones on the diagonal
                        u is upper triangular 
                        p is a permutation matrix.

                        The idea behind the algorithm is to first
                        do Gaussian elimination to obtain an upper
                        triangular matrix u and lower triangular matrix
                        r such that r * self = u, then by inverting r to
                        get l = r ^-1 we obtain self = r^-1 * u = l * u.
                        Note tha since r is lower triangular its
                        inverse must also be lower triangular.

                        Where does the p come in?  Well, with some
                        matrices our technique doesn't work due to
                        zeros appearing on the diagonal of r.  So we
                        apply some permutations to the orginal to
                        prevent this.
                        
        """
        upper = self.Copy()
        resultInv = self.MakeSimilarMatrix(self.Size(),'i')
        perm = self.MakeSimilarMatrix((self.rows,self.rows),'i')

        (rowIndex,colIndex) = (0,0)
        lastRow = self.rows - 1
        lastCol = self.cols - 1
        while( rowIndex < lastRow and colIndex < lastCol ):
            leader = upper.FindRowLeader(rowIndex,colIndex)
            if (leader < 0):
                colIndex = colIndex+1
                continue            
            if (leader != rowIndex):
                upper.SwapRows(leader,rowIndex)
                resultInv.SwapRows(leader,rowIndex)
                perm.SwapRows(leader,rowIndex)
            (rowIndex,colIndex) = (
                upper.PartialLowerGaussElim(rowIndex,colIndex,resultInv))

        lower = self.MakeSimilarMatrix((self.rows,self.rows),'i')
        resultInv.LowerGaussianElim(lower)
        resultInv.UpperInverse(lower)
        # possible optimization: due perm*lower explicitly without
        # relying on the * operator.
        return (perm*lower, upper, perm)
    
    def Solve(self,b):
        """
        Solve(self,b):

        b:   A list.
        
        Returns the values of x such that Ax = b.

        This is done using the LUP decomposition by 
        noting that Ax = b implies PAx = Pb implies LUx = Pb.
        First we solve for Ly = Pb and then we solve Ux = y.
        The following is an example of how to use Solve:

>>> # Floating point example
>>> import genericmatrix
>>> A = genericmatrix.GenericMatrix(size=(2,5),str=lambda x: '%.4f' % x)
>>> A.SetRow(0,[0.0, 0.0, 0.160, 0.550, 0.280])
>>> A.SetRow(1,[0.0, 0.0, 0.745, 0.610, 0.190])
>>> A
<matrix
 0.0000 0.0000 0.1600 0.5500 0.2800
 0.0000 0.0000 0.7450 0.6100 0.1900>
>>> b = [0.975, 0.350]
>>> x = A.Solve(b)
>>> z = A.LeftMulColumnVec(x)
>>> diff = reduce(lambda xx,yy: xx+yy,map(lambda aa,bb: abs(aa-bb),b,z))
>>> diff > 1e-6
0
>>> # Boolean example
>>> XOR = lambda x,y: x^y
>>> AND = lambda x,y: x&y
>>> DIV = lambda x,y: x  
>>> m=GenericMatrix(size=(3,6),zeroElement=0,identityElement=1,add=XOR,mul=AND,sub=XOR,div=DIV)
>>> m.SetRow(0,[1,0,0,1,0,1])
>>> m.SetRow(1,[0,1,1,0,1,0])
>>> m.SetRow(2,[0,1,0,1,1,0])
>>> b = [0, 1, 1]
>>> x = m.Solve(b)
>>> z = m.LeftMulColumnVec(x)
>>> z
[0, 1, 1]

        """
        assert self.cols >= self.rows
        
        (L,U,P) = self.LUP()
        Pb = P.LeftMulColumnVec(b)
        y = [0]*len(Pb)
        for row in range(L.rows):
            y[row] = Pb[row]
            for i in range(row+1,L.rows):
                Pb[i] = L.sub(Pb[i],L.mul(L[i,row],Pb[row]))
        x = [0]*self.cols
        curRow = self.rows-1

        for curRow in range(len(y)-1,-1,-1):
            col = U.FindColLeader(curRow,0)
            assert col > -1
            x[col] = U.div(y[curRow],U[curRow,col])
            y[curRow] = x[col]
            for i in range(0,curRow):
                y[i] = U.sub(y[i],U.mul(U[i,col],y[curRow]))
        return x


def DotProduct(mul,add,x,y):
    """
    Function:    DotProduct(mul,add,x,y)
    Description: Return the dot product of lists x and y using mul and
                 add as the multiplication and addition operations.
    """
    assert len(x) == len(y), 'sizes do not match'
    return reduce(add,map(mul,x,y))

class GenericMatrixTester:
    def DoTests(self,numTests,sizeList):
        """
        Function:       DoTests(numTests,sizeList)

        Description:    For each test, run numTests tests for square
                        matrices with the sizes in sizeList.
        """

        for size in sizeList:
            self.RandomInverseTest(size,numTests)
            self.RandomLUPTest(size,numTests)
            self.RandomSolveTest(size,numTests)
            self.RandomDetTest(size,numTests)


    def MakeRandom(self,s):
        import random 
        r = GenericMatrix(size=s,fillMode=lambda x,y: random.random(),
                          equalsZero = lambda x: abs(x) < 1e-6)
        return r

    def MatAbs(self,m):
        r = -1
        (N,M) = m.Size()
        for i in range(0,N):
            for j in range(0,M):
                if (abs(m[i,j]) > r):
                    r = abs(m[i,j])
        return r

    def RandomInverseTest(self,s,n):
        ident = GenericMatrix(size=(s,s),fillMode='i')
        for i in range(n):
            m = self.MakeRandom((s,s))
            assert self.MatAbs(ident - m * m.Inverse()) < 1e-6, (
                'offender = ' + `m`)

    def RandomLUPTest(self,s,n):
        ident = GenericMatrix(size=(s,s),fillMode='i')
        for i in range(n):
            m = self.MakeRandom((s,s))
            (l,u,p) = m.LUP()
            assert self.MatAbs(p*m - l*u) < 1e-6, 'offender = ' + `m`

    def RandomSolveTest(self,s,n):
        import random
        if (s <= 1):
            return
        extraEquations=3
        
        for i in range(n):
            m = self.MakeRandom((s,s+extraEquations))
            for j in range(extraEquations):
                colToKill = random.randrange(s+extraEquations)
                for r in range(m.rows):
                    m[r,colToKill] = 0.0
            b = map(lambda x: random.random(), range(s))
            x = m.Solve(b)
            z = m.LeftMulColumnVec(x)
            diff = reduce(lambda xx,yy:xx+yy, map(lambda aa,bb:abs(aa-bb),b,z))
            assert diff < 1e-6, ('offenders: m = ' + `m` + '\nx = ' + `x`
                                 + '\nb = ' + `b` + '\ndiff = ' + `diff`)

    def RandomDetTest(self,s,n):
        for i in range(n):
            m1 = self.MakeRandom((s,s))
            m2 = self.MakeRandom((s,s))
            prod = m1 * m2
            assert (abs(m1.Determinant() * m2.Determinant()
                        - prod.Determinant() )
                    < 1e-6), 'offenders = ' + `m1` + `m2`


license_doc = """
  This code was originally written by Emin Martinian ([e-mail]).
  You may copy, modify, redistribute in source or binary form as long
  as credit is given to the original author.  Specifically, please
  include some kind of comment or docstring saying that Emin Martinian
  was one of the original authors.  Also, if you publish anything based
  on this work, it would be nice to cite the original author and any
  other contributers.

  There is NO WARRANTY for this software just as there is no warranty
  for GNU software (although this is not GNU software).  Specifically
  we adopt the same policy towards warranties as the GNU project:

  BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW.  EXCEPT WHEN
OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES
PROVIDE THE PROGRAM 'AS IS' WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED
OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.  THE ENTIRE RISK AS
TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU.  SHOULD THE
PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,
REPAIR OR CORRECTION.

  IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,
INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING
OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED
TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY
YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER
PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
POSSIBILITY OF SUCH DAMAGES.
"""

testing_doc = """
The GenericMatrixTester class contains some simple
testing functions such as RandomInverseTest, RandomLUPTest,
RandomSolveTest, and RandomDetTest which generate random floating
point values and test the appropriate routines.  The simplest way to
run these tests is via

>>> import genericmatrix
>>> t = genericmatrix.GenericMatrixTester()
>>> t.DoTests(100,[1,2,3,4,5,10])

# runs 100 tests each for sizes 1-5, 10
# note this may take a few minutes

If any problems occur, assertion errors are raised.  Otherwise
nothing is returned.  Note that you can also use the doctest
package to test all the python examples in the documentation
by typing 'python ffield.py' or 'python -v ffield.py' at the
command line.
"""


# The following code is used to make the doctest package
# check examples in docstrings when you enter

__test__ = {
    'testing_doc' : testing_doc
}

def _test():
    import doctest, genericmatrix
    return doctest.testmod(genericmatrix)

if __name__ == "__main__":
    _test()
    print 'Tests Passed.'