/home/mark/model/software/ScrumPy/ScrumPy/Util/DynMatrix.py

"""

ScrumPy -- Metabolic Modelling with Python

Copyright Mark Poolman 1995 - 2002 Copy

 This file is part of ScrumPy.

    ScrumPy 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.

    ScrumPy 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 ScrumPy; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

"""
"""
M. Poolman 10/01/01
Dynamic matrix class
"""




import  types,math,sys,exceptions
Pi = math.pi
PiOver2 = Pi/2


import  Sci,Seq,Sorter,Types,Format, LaTeX,  Tree

try:
    from ScrumPy import wxGUI as GUI
    HaveGUI = True
except:
    HaveGUI = False

DefaultConv = Types.ArbRat
DynMatrixError = "Dynamic Matrix Error"
Tol = 1e-9 # tolerance for floating point error checks etc.

def DepWarn(Str):
    sys.stderr.write( "!! Deprecated " + Str+" !!" )


def GetIdent(n, Conv=DefaultConv, mul=1):
    rv = matrix(nrow=n,ncol=n,Conv=Conv)
    rv.MakeIdent(mul)
    return rv

def FromVec(vec, Conv=None):

    if not Conv:
        vtype = type(vec[0]).__name__
        try:
            Conv = {
                'float' : float,
                'int' : int,
                'long int' : long
                }[vtype]
        except:
            Conv = DefaultConv

    size = len(vec)
    rv = matrix(size,size,Conv)
    for n in range(size):
        rv[n,n] = Conv(vec[n])

    return rv






class matrix:

    def __init__(self, nrow=0, ncol=0,rnames=[],cnames=[], Conv=DefaultConv, FromMtx=None):

        self.Conv = Conv
        self.rnames = []
        self.cnames = []
        self.rows = []


        if FromMtx != None:
            try:
                self.FromDynMatrix(FromMtx)
                return
            except:
                try:
                    self.FromNumPyMtx(FromMtx)
                except:
                    self.FromList(FromMtx)

            if len(rnames) != 0 and len(rnames) != len(self):
                raise exceptions.IndexError, "len(rnames) != len(FromMtx)"
            else:
                self.rnames = rnames[:]

            if len(cnames) != 0 and len(cnames) != len(self.rows[0]):
                raise exceptions.IndexError, "len(cnames) != len(FromMtx[0])"
            else:
                self.cnames = cnames[:]

            return

        if len(rnames) != 0:
            self.rnames = rnames[:]
        else:
            for r in range(nrow):
                self.rnames.append("row"+str(r))

        if len(cnames) != 0:
            self.cnames = cnames[:]
        else:
            for c in range(ncol):
                self.cnames.append("col"+str(c))

        self.MakeZero()



    def _getrc(self,rc):
        if type(rc) == types.TupleType:
            r = rc[0]
            c = rc[1]
        else:
            r = rc
            c = None

        return r,c


    def __len__(self):
        return len(self.rows)


    def __getitem__(self,rc):
        r,c = self._getrc(rc)

        try:
            row =  self.rows[r]
        except:
            try:
                row =  self.rows[self.rnames.index(r)]
            except:
                raise IndexError, r

        if  c==None:
            return row

        try:
            return row[c]
        except:
            try:
                return row[self.cnames.index(c)]
            except:
                raise IndexError, c


    def __setitem__(self,rc, val):

        r,c = self._getrc(rc)

        if type(r) == types.StringType:
            try:
                r = self.rnames.index(r)
            except:
                raise IndexError, r

        if  c==None:
            self.rows[r] = map(self.Conv, val)
        else:
            if type(c) == types.StringType:
                try:
                    c = self.cnames.index(c)
                except:
                    raise IndexError, c

            self.rows[r][c] = self.Conv(val)


    def __getslice__(self,lo,hi):
        rv = matrix(Conv=self.Conv)
        rv.rows = self.rows[lo:hi]
        rv.rnames = self.rnames[lo:hi]
        rv.cnames = self.cnames[:]
        return rv


    def __str__(self):
        retstr = "    "
        for cname in self.cnames:
            retstr += str(cname) + " "
        retstr += "\n"
        for r in range(len(self.rnames)):
            retstr += "\n " + str(self.rnames[r])
            for el in self[r]:
                retstr +=" "+str(el)+ " "

        return retstr


    def __repr__(self):
        try:
            return str(self)
        except:
            return "can't generate string rep'n for this matrix"

    def Show(self):
        if HaveGUI:
            GUI.Mtx.SimpleView(self).Show()
        else:
            print self


    def index(self,row):

        for i in range(len(self.rows)):
            if Seq.AreEqual(self.rows[i], row):
                return i

        raise "ValueError", "matrix.index(row): row not in matrix"



    def DeQuote(self):
        """ replace double quoted strings in row and col names with unquoted
            (preserves mappings in self.QuoteD) """

        self.QuoteD = {}
        rnamesd = map(Format.StripQuote,self.rnames)
        cnamesd = map(Format.StripQuote,self.cnames)

        for n in range(len(rnamesd)):
            self.QuoteD[rnamesd[n]] = self.rnames[n]
            self.QuoteD[self.rnames[n]] = rnamesd[n]

        for n in range(len(cnamesd)):
            self.QuoteD[cnamesd[n]] = self.cnames[n]
            self.QuoteD[self.cnames[n]] = cnamesd[n]

        self.rnames = rnamesd
        self.cnames = cnamesd



    def IsDeQuoted(self):
        """ bool: self is DeQuote()ed """

        return hasattr(self, "QuoteD")


    def ReQuote(self):
        """ undo DeQuote """

        if  self.IsDeQuoted():

            for n in range(len(self.rnames)):
                name = self.rnames[n]
                if self.QuoteD.has_key(name):  # because might have added something sunce DeQuote
                    self.rnames[n] = self.QuoteD[name]

            for n in range(len(self.cnames)):
                name = self.cnames[n]
                if self.QuoteD.has_key(name):
                    self.cnames[n] = self.QuoteD[name]

            del self.QuoteD

    def ResetRowNames(self):
        """ pre: True
           post: self.rnames==["Row_0".."Row_n"] """

        self.rnames = map(lambda n: "Row_"+str(n),range(len(self)))


    def ChangeType(self, Conv):
        self.Conv = Conv
        for r in self.rnames:
            self[r] = self.MakeSeqSelfType(self[r])

    def MakeSeqSelfType(self, seq):

        return map(self.Conv, seq)


    def MakeZero(self):
        Conv = self.Conv
        self.rows = []
        for n in range(len(self.rnames)):
            self.rows.append(map(self.Conv, [0]*len(self.cnames)))


    def Clear(self):
        "erase all row and col data"

        self.cnames, self.rnames, self.rows = [],[],[]



    def Copy(self, Conv=None,tx=False):#,_Super=None):
        """ pre: True
           post: Copy of self, with Element type Conv (self.Conv if Conv==None)
                    tx => Copy is transposed """

        if Conv==None:
            Conv=self.Conv


        rv = self.__class__(Conv=Conv)
        rv.cnames = self.cnames[:]
        rv.rnames = self.rnames[:]
        if Conv==self.Conv:
            for row in self.rows:
                rv.rows.append(row[:])
        else:
            for row in self.rows:
                rv.rows.append(map(Conv, row))

        if tx:
            rv.Transpose()

        return rv

    def ReadFile(self, FileName):
        """ pre: (FileName != None) -> (File conforms to FileSpec)
           post: contents of self are those of FileName  """

        f = open(FileName, "r")
        self.Clear()
        Titles = f.readline()
        for t in Titles.split():
            self.NewCol(name=t)

        for row in f.readlines():
            row = row.split()
            try:
                self.NewRow(row[1:], row[0])
            except:
                print "problem with row", row, "ignoring it"



    def WriteFile(self,File,WithHash=False,WithRNames=True,Delim="\t"):
        """ pre: (ColOrder == None) || (complete and exclusive list of column headings of self)
           post: FileName contains FileSpec representation of self
                 (ColOrder != None) => Columns in ColOrder
                 ELSE Exception """

        if type(File) == types.StringType:
            f = open(File, "w")
        else:
            f = File

        first = {True:["#"],
                 False:[""]
                 }[WithHash]

        print >>f, Delim.join(first+self.cnames)

        if WithRNames:
            for n in range(len(self)):
                print >>f, Delim.join([self.rnames[n]]+map(str,self.rows[n]))
        else:
            for n in len(self):
                print >>f, Delim.join(map(str,self.rows[n]))



    def ToFile(self, fname):
        f = open(fname,"w")
        f.write(str(self))

    def ToList(self):
        """ return a list of lists of tuples. Each list corresponds to a  row,
             each tuple is (name,value) of the column and value of the non zero elements within that row """

        rv = []
        for row in self:
            curlist = []
            curidxs = Seq.IdxListOfMatches(row, lambda x: x!=0)
            for idx in curidxs:
                curlist.append((self.cnames[idx], row[idx]))
            rv.append(curlist)

        return rv



    def ToDict(self):
        """ return dictionary in which keys are row names and values lists of tuples (colnames, elval) for elval != 0 """

        l = self.ToList()
        rv = {}
        for idx in range(len(self)):
            rv[self.rnames[idx]] = l[idx]
        return rv


    def UpTri(self):
        """ return upper triangular elements of self as a single list """

        rv = []
        for i in range(len(self)):
            rv.extend(self[i][i+1:])
        return rv


    def FromDicDic(self,d):
        """ pre: d is a dictionary of rows,
        in which each value is a dictionary mapping col names to values
            post: self'[r,c] == d[r,c] for non-zero values, all other values in self are zero """

        self.rnames = []
        self.rows = []
        self.cnames = []

        for kr in d.keys():              # kr - key of row
            self.NewRow(name=kr)
            rd = d[kr]                      # rd - row dictionary
            for kc in rd.keys():        # kc - key of col
                if not kc in self.cnames:
                    self.NewCol(name=kc)

                self[kr,kc] = d[kr][kc]




    def ToColList(self):
        """ as ToList, but returns list of columns """

        rv = []
        for cn in self.cnames:
            col = self.GetCol(cn)
            curidxs = Seq.IdxListOfMatches(col, lambda x: x!=0)
            curlist = []
            for idx in curidxs:
                curlist.append((self.rnames[idx],col[idx]))
            rv.append(curlist)
        return rv


    def ToColDict(self):
        """ return dictionary in which keys are col names and values lists of tuples (rownames, elval) for elval != 0 """

        l = self.ToColList()
        rv = {}
        for idx in range(len(self.cnames)):
            rv[self.cnames[idx]] = l[idx]
        return rv





    def ToSciMtx(self,Conv=None):

        sys.stderr.write("ToSciMtx: why do we need this ?\n")

        if not Conv:
            return Sci.matrix(self.rows)

        return Sci.matrix(self.Copy(Conv).rows)


    def FromNumPyMtx(self, mtx, Conv=None):
        sys.stderr.write("FromNumPyMtx: use from SciMtx ?\n")
        return self.FromSciMtx(mtx,  Conv)

    def FromSciMtx(self, mtx, Conv=None):

        sys.stderr.write("FromSciMtx: why do we need this ?\n")

        nrow,ncol = mtx.shape
        if Conv!=None:
            self.Conv = Conv
        self.cnames = map(lambda x: "col"+str(x),range(ncol))
        self.rnames = map(lambda x: "row"+str(x),range(nrow))

        self.rows = mtx.tolist() # ! no - do this in the Sci module.






    def ToOctaveMtx(self, name, Conv = float):
        """ pre: True
           post: returns a string that when printed assigns the octave matrix, name, to current contents
                    using Conv as a type converter
        """
        rv = name + "= ["
        for row in self.rows:
            rv += "\n  ["
            for el in  row[:-1]:
                rv += str(Conv(el)) + ", "
            rv += str(Conv(row[-1])) + "]"
        rv += "\n]\n"

        return rv


    def ToCMtx(self, name, Conv = float):
        """ pre: True
           post: returns a string that when printed declares a C matrix, name, of current contents
                    type will need editing for non-default Conv, which as ever, is the type converter
        """

        c_nrows = name + "_nRows"
        c_ncols = name + "_nCols"
        strlenrow  = str(len(self.rows))
        strlencol = str(len(self.rows[0]))
        rv = "  const int "

        rv += c_nrows + " = " + strlenrow + ", " +   c_ncols + " = " + strlencol + ";"

        rv += "\n  double " +  name + "[" + c_nrows + "]" + "[" + c_ncols + "] = {\n"
        for row in self.rows[:-1]:
            rv += "\t{"
            for el in row[:-1]:
                rv += str(Conv(el)) + ", "
            rv += str(Conv(row[-1])) + "},\n"

        row = self.rows[-1]
        rv += "\t{"
        for el in row[:-1]:
            rv += str(Conv(el)) + ", "
        rv += str(Conv(row[-1])) + "} \n  } ;"

        return rv


    def ToStr(self):
        """ return  an easy-to-interpret string repn of self, reflecting information content
             rather than matrix notation -NOT- the same as __str__ or __repr__"""

        rv = ""
        list = self.ToList()
        for row in list:
            rowstr = "    "
            for el in row:
                elstr = str(el[1]) + " " + el[0]
                if el[1] < 0:
                    rowstr = elstr + "  " + rowstr
                else:
                    rowstr += "  " + elstr
            rv += rowstr + "\n"
        return rv


    def ToTexTab(self, RowLabels=True):
        "return a string representation in Latex tabular form"

        if RowLabels:
            colheads = ["Row"] + self.cnames
        else:
            colheads = self.cnames[:]

        ncols = len(colheads)
        rv = LaTeX.Tabular(ncols)
        rv.AddRow(colheads)

        if self.Conv == float:
            TeXConv = Float2Tex
        else:
            TeXConv = str

        if RowLabels:
            for r in self.rnames:
                rv.AddRow([r]+map(TeXConv,self[r]))
        else:
            for r in self.rows:
                rv.AddRow(map(TeXConv,r))

        return str(rv)



    def FromDynMatrix(self,other):

        self.Conv = other.Conv
        self.cnames = other.cnames[:]
        self.rnames = other.rnames[:]
        self.rows = []
        for r in other.rows:
            self.rows.append(r[:])




    def InitRows(self, rnames):
        """ pre: rnames is a list of strings
           post: self is an empty (r x 0) matrix, using rnames as row names """
        self.rnames = rnames[:]
        self.rows = []



    def Transpose(self):
        NewRows = []
        for c in range(len(self.cnames)):
           NewRows.append(self.GetCol(c))
        self.rows = NewRows

        tmp = self.cnames
        self.cnames = self.rnames
        self.rnames = tmp


    def RowReorder(self, NewOrder):
        tmp = matrix(Conv=self.Conv,ncol = len(self.cnames))

        for n in NewOrder:
            if type(n)==int:
                name = self.rnames[n]
            else:
                name = n
            tmp.NewRow(self[n], name)

        self.rows = tmp.rows
        self.rnames = tmp.rnames


    def ColReorder(self, NewOrder):
        self.Transpose()
        self.RowReorder(NewOrder)
        self.Transpose()

    def RandRowOrder(self):
        self.RowReorder(Seq.RandOrder(self.rnames))

    def RandColOrder(self):
        self.ColReorder(Seq.RandOrder(self.cnames))

    def RandOrder(self):
        self.RandRowOrder()
        self.RandColOrder()


    def SortBy(self,cname,Descend=False):

        origc = self.GetCol(cname)
        sortc = sorted(origc,reverse=Descend)
        mark = min(origc)-1  # ie a number known not to be in the list

        neworder=[]
        for x in sortc:
            idx = origc.index(x)
            neworder.append(idx)
            origc[idx] = mark

        self.RowReorder(neworder)




    def RowNamesFromIdxs(self,idxs=[]):
        rv = []
        for i in idxs:
            rv.append(self.rnames[i])
        return rv


    def ColNamesFromIdxs(self,idxs=[]):
        rv = []
        for i in idxs:
            rv.append(self.cnames[i])
        return rv


    def RowNZIdxs(self,r,Abs=True):
        "list of indices to nonzero elements of r, sorted by ascending (Abs => absolute) element value "

        return Seq.NZIdxs(self[r],Abs)


    def ColNZIdxs(self,c,Abs=True):
        "list of indices to nonzero elements of c, sorted by ascending (Abs => absolute) element value "

        return Seq.NZIdxs(self.GetCol(c),Abs)



    def Sort(self, fun, **fundat):
        """ sort rows in ascending order of fun(row,mtx,fundat) where:
            fun returns a numeric type
            row is an index into mtx
            mtx is self
            fundat is optional named args"""

        s = Sorter.Sorter()
        for r in self.rnames:
            s.append(r)
        self.RowReorder(s.Sort(fun, mtx=self, **fundat))




    def MakeIdent(self, mul):

        r = 0
        for row in self.rows:
            row = map(lambda x: 0, row)
            row[r] = mul
            self.rows[r] = self.MakeSeqSelfType(row)
            r = r + 1

    def SetEl(self, row, col, val):

        DepWarn("DynMatrix.matrix.SetEl()")

        if type(row) != types.IntType:
            row = self.rnames.index(row)
        if type(col) != types.IntType:
            col = self.cnames.index(col)

        self.rows[row][col] = self.Conv(val)


    def NewEl(self, rname, cname, val):
        """ pre: true
           post: self[rname,cname] == val
                 (rows and cols created if needed)
        """

        if not rname in self.rnames:
            self.NewRow(name=rname)
        if not cname in self.cnames:
            self.NewCol(name=cname)
        self[rname,cname] = self.Conv(val)



    def NewCol(self,  col=None, name=None):

        if name == None:
            name = "col" + str(len(self.cnames)+1)

        self.cnames.append(name)

        if col == None:
            col = [0] *len(self.rnames)

        if len(self.rnames)==0:
            for i in range(len(col)):
                self.rows.append([self.Conv(col[i])])
                self.rnames.append("row_" + str(i))

        elif len(col) != len(self.rnames):
            raise exceptions.IndexError, "wrong length for new col"

        else:
            for i in range(len(self.rows)):
                self[i].append(self.Conv(col[i]))




    def DelCol(self, c=0):

        if type(c) == types.StringType:
            c = self.cnames.index(c)

        del self.cnames[c]
        for row in self.rows:
            del row[c]


    def GetCol(self, c):
        if type(c) == types.StringType:
            c = self.cnames.index(c)
        col = []
        for row in self.rows:
            col.append(row[c])

        return col


    def SetCol(self, c, vals):
        if type(c) == types.StringType:
            c = self.cnames.index(c)
        vals = self.MakeSeqSelfType(vals)
        r = 0
        for row in self.rows:
            row[c] = vals[r]
            r = r + 1


    def AddCol(self, c, vals=None, k = None):
        #DepWarn("AddCol")
        if vals !=None:
            self.SetCol(c, map(lambda x,y: x+y, self.GetCol(c), vals))
        else:
            self.SetCol(c, map(lambda x: x+k, self.GetCol(c)))

    def SubCol(self, c, vals=None, k=None):
        #DepWarn("SubCol")
        if vals != None:
            self.SetCol(c, map(lambda x,y: x-y, self.GetCol(c), vals))
        else:
            self.SetCol(c, map(lambda x: x-k, self.GetCol(c)))

    def MulCol(self, c, vals=None, k=None):
        #DepWarn("MulCol")
        if vals!=None:
            self.SetCol(c, map(lambda x,y: x*y, self.GetCol(c), vals))
        else:
            self.SetCol(c, map(lambda x: x*k, self.GetCol(c)))

    def DivCol(self, c, vals=None, k=None):
        #DepWarn("DivCol")
        if vals!=None:
            self.SetCol(c, map(lambda x,y: x/y, self.GetCol(c), vals))
        else:
            self.SetCol(c, map(lambda x: x/k, self.GetCol(c)))


    def FunCol(self, c, Fun):
        self.SetCol(c, map(Fun, self.GetCol(c)))

    def MaxInCol(self, c):
        return max(self.GetCol(c))

    def MinInCol(self, c):
        return min(self.GetCol(c))

    def PosOfValsInCol(self, c, val):
        return Seq.PosOfValsInSeq(self.GetCol(c), val)

    def NZeroesInCol(self,c):
        return len(self.PosOfValsInCol(c,0))


    def NonZeroesInCols(self):
        """ returns a list vector, one element per col,
        each element being the number of non-zero elements in the coresponding col """

        rv = []
        for c in range(len(self.cnames)):
            rv.append(self.NZeroesInCol(c))

        return rv

    def NamesNonZeroesInCols(self):
        """ return a list (one el per col) lists of row names that have no zero element in corresponding col """

        rows = self.rows
        rnames = self.rnames
        rv = []

        for c in range(len(self.cnames)):
            res = []
            for r in range(len(self)):
                if rows[r][c] != 0:
                    res.append(rnames[r])
            rv.append(res)
        return rv


    def MaxInRange(self,names):
        """ pre: self.rnames == self.cnames, Set.IsSubset(names, self.rnames)
           post: greatest value of self[names[x],names[y]] """

        rv = 0.0
        for n1 in names:
            for n2 in names[names.index(n1):]:
                if self[n1,n2] > rv:
                    rv = self[n1,n2]

        return rv




    def SwapCol(self, c1, c2):

        if type(c1).__name__ == 'string':
            c1 = self.cnames.index(c1)

        if type(c2).__name__ == 'string':
            c2 = self.cnames.index(c2)

        tmp = self.GetCol(c1)
        self.SetCol(c1, self.GetCol(c2))
        self.SetCol(c2, tmp)

        tmp = self.cnames[c1]
        self.cnames[c1] = self.cnames[c2]
        self.cnames[c2] = tmp


    def AugCol(self, m):
        for c in range(len(m.cnames)):
            self.NewCol(m.GetCol(c),m.cnames[c])

    def MakeLastCol(self, c):
        if type(c) == types.StringType:
            cname = c
            c = self.cnames.index(c)
        else:
            cname = self.cnames[c]

        col = self.GetCol(c)
        self.DelCol(c)
        self.NewCol(col, cname)


    def MakeFirstCol(self, c):
        if type(c) == types.StringType:
            cname = c
            c = self.cnames.index(c)
        else:
            cname = self.cnames[c]

        col = self.GetCol(c)
        self.DelCol(c)
        self.cnames.insert(0,cname)
        for idx in range(len(col)):
            self.rows[idx].insert(0,col[idx])





    def NewRow(self,  r=None, name=None):

        if name == None:
            name = "row" + str(len(self.rows)+1)

        self.rnames.append(name)

        if r == None:
            r = [0] *len(self.cnames)
        elif len(r) != len(self.cnames):
            raise exceptions.IndexError, "wrong length for new row"

        self.rows.append(self.MakeSeqSelfType(r))


    def DelRow(self, r):
        if type(r) == types.StringType:
            r = self.rnames.index(r)

        del self.rows[r]
        del self.rnames[r]


    def GetRow(self, r):
        if type(r) == types.StringType:
            r = self.rnames.index(r)
        return(self.rows[r][:])


    def SetRow(self, r, vals):
        if type(r) == types.StringType:
            r = self.rnames.index(r)
        self.rows[r] = self.MakeSeqSelfType(vals)


    def AddRow(self, r, vals=None, k=None):
        if vals and k:
            self[r] = map(lambda x,y: x+y, self[r], map(lambda x,mul=k: x*mul,vals))
        elif vals:
            self[r] = map(lambda x,y: x+y, self[r], vals)
        else:
            self[r] = map(lambda x,add=k: x+add,vals)


    def InsertRow(self, pos, Val=0, Name=None):

        r = map(lambda x,v=Val: v, range(0,len(self.cnames)))
        self.rows.insert(pos, self.MakeSeqSelfType(r))
        if Name:
            self.rnames.insert(pos,Name)


    def SubRow(self, r, vals):
        DepWarn("SubRow")
        self.rows[r] = map(lambda x,y: x-y, self.rows[r], vals)


    def MulRow(self, r, vals=None, k=None):
        if vals and k:
            self[r] = map(lambda x,y: x*y, self[r], map(lambda x,mul=k: x*mul,vals))
        elif vals:
            self[r] = map(lambda x,y: x*y, self[r], vals)
        else:
            self[r] = map(lambda x,mul=k: x*mul,self[r])


    def DivRow(self, r, vals=None, k=None):
        if type(r) != types.IntType:
            r = self.rnames.index(r)
        if vals and k:
            self.rows[r] = map(lambda x,y: x/y, self.rows[r], map(lambda x,mul=k: x*mul,vals))
        elif vals:
            self.rows[r] = map(lambda x,y: x/y, self.rows[r], vals)
        else:
            self.rows[r] = map(lambda x,den=k: x/den,self.rows[r])


    def AddAndRemoveRows(self, rows):
        """ pre: rows is subset of self.rnames
           post: self'[rows[0]] """

        r0 = rows[0]
        for r in rows[1:]:
            self.AddRow(r0,self[r])
        for r in rows[1:]:
            self.DelRow(r)

    def FunRow(self, r, fun):
        self.rows[r] = map(fun, self.rows[r])

    def MaxInRow(self, r):
        return max(self[r])

    def MinInRow(self, c):
        return min(self[r])

    def PosOfValsInRow(self, r, val):
        return Seq.PosOfValsInSeq(self.rows[r], val)



    def FunAllRows(self, fun):
        for r in range(len(self.rows)):
            self.rows[r] = fun(self.rows[r])

    def DupRows(self):
        """ pre: True
           post: self.DupRows()  a (possibly empty) list of lists of names of rows with identical elements
        """

        rv = []
        rows = self.rows
        rnames = self.rnames
        ridxs = range(len(self.rows))
        while len(ridxs)>0:
            r = ridxs[0]
            duplist = [r]
            for r2 in ridxs[1:]:
                if Seq.AreEqual(rows[r], rows[r2]):
                    duplist.append(r2)

            for d in duplist:
                del ridxs[ridxs.index(d)]
            if len(duplist)> 1:
                rv.append(map(lambda x,rn = self.rnames: rn[x],duplist))

        return rv



    def DelDupRows(self):
        """pre: True
          post: len(self.DupRows())==0
        """

        for dups in self.DupRows():
            for d in dups[1:]:  # leave the first one !!
                self.DelRow(d)


    def MaxCountRow(self, val):
        """ return the row name with the greatest number of elements of value val
            in the event of a tie the first matching row name is returned
            if no rows contain val the first row name is returned """

        rv = self.rnames[0]
        best = 0
        for r in self.rnames:
            hits = self[r].count(val)
            if hits > best:
                best = hits
                rv = r
        return rv



    def RowsWithSingCols(self):
        "return list of idxs of rows with at least one non-zero element unique to the column"
        rv = []
        for c in self.cnames:
            col = self.GetCol(c)
            idx = Seq.IdxOfUnique(col, lambda x: x!=0)
            if idx != None and not idx in rv:
                rv.append(idx)
        return rv


    def SwapRow(self, r1, r2):

        if type(r1).__name__ == 'string':
            r1 = self.rnames.index(r1)

        if type(r2).__name__ == 'string':
            r2 = self.rnames.index(r2)

        tmp = self.rows[r1]
        self.rows[r1] = self.rows[r2]
        self.rows[r2] = tmp

        tmp = self.rnames[r1]
        self.rnames[r1] = self.rnames[r2]
        self.rnames[r2] = tmp


    def MakeFirstRow(self, r):

        row = self[r]
        if type(r)==types.StringType:
            name = r
        else:
            name = self.rnames[r]

        self.DelRow(r)
        self.rows.insert(0,row)
        self.rnames.insert(0,name)


    def MakeLastRow(self, r):

        row = self[r]
        if type(r)==types.StringType:
            name = r
        else:
            name = self.rnames[r]

        self.DelRow(r)
        self.rows.append(row)
        self.rnames.append(name)


    def AugRow(self, m):
        for r in range(len(m.rnames)):
            self.NewRow(m.rows[r], m.rnames[r] )


    def NonZeroR(self, seed=[], thresh=1e-6):
        """pre: mtx has no empty rows, and numrical elements
          post: return a seed + list of list of rownames for which r(a,b)> thresh
                  where r is sample correlaation coeff,
                  ** self is empty  -DESTRUCTIVE ! ** """


        lenself = len(self)

        if lenself==0:
            return seed
        elif lenself==1:
            return seed+self.rnames[:]

        FirstRow = self.rnames[0]
        rv = [FirstRow]
        for ThisRow in self.rnames[1:]:
            if abs(Seq.CosTheta(self[FirstRow],self[ThisRow])) > thresh:
                rv.append(ThisRow)
                self.DelRow(ThisRow)
        self.DelRow(FirstRow)
        return self.NonZeroR(seed+[rv], thresh)



    def RowDiffMtxFull(self, fun= None, Conv=float):

        rv = matrix(rnames=self.rnames, cnames=self.rnames, Conv=Conv)
        lenself = len(self)
        for r1 in range(lenself):
            for r2 in range(lenself):
                rv[r1,r2] = fun(self[r1],self[r2])

        return rv




    def RowDiffMtx(self, fun=None,Conv=float):
        """ pre:  fun(self[r1],self[r2]) for 0 <= r1,r2  < len(self) returns type Conv
                     fun(r1,r2) == fun(r2,r1)
           post: returns M such that M[r1,r2] = fun(self[r1],self[r2])
            """

        if fun == None:
            fun = lambda x,y: abs(Sci.CosTheta(x,y))

        rv = matrix(rnames=self.rnames, cnames=self.rnames, Conv=Conv)
        lenself = len(self)
        for r1 in range(lenself):
            for r2 in range(r1, lenself):
                rv.rows[r1][r2] = rv.rows[r2][r1] = fun(self[r1],self[r2])

        return rv


    def RowDiffMtx2(self,fun):
        """ as above, but fun takes self, r1, and r2 as args, and no type conversion is done """

        rv = matrix(rnames=self.rnames, cnames=self.rnames, Conv=self.Conv)
        lenself = len(self)
        for r1 in range(lenself):
            for r2 in range(r1+1, lenself):
                rv[r1,r2] = rv[r2,r1] =  fun(self,r1,r2)

        return rv


    def ToNJTree(self):
        """pre:  self is square difference matrix
          post:  A string representing the NJ tree, based on row differences, as calculated by fun """

        def FindNearestRows(mtx):

            rv = mtx.rnames[0] ,mtx.cnames[1],mtx[0,1]
            if len(mtx) == 2:
                return rv

            for r in range(len(mtx)):
                for c in range(r):
                    el =  mtx[r,c]
                    if el == 0:
                       return mtx.rnames[r],mtx.cnames[c],0
                    else:
                        if el < rv[2]:
                            rv = mtx.rnames[r],mtx.cnames[c],mtx[r,c]
            return rv

        nodes = {}
        mtx = self.Copy()
        for r in mtx.rnames:
            nodes[r] = t = Tree.Tree()
            t.add_feature("name", r)

        while len(mtx) >1:
            print ".",
            m0,m1,dist = FindNearestRows(mtx)
            ch0 = nodes[m0]
            ch0.add_feature("dist", dist)
            ch1 = nodes[m1]
            ch1.add_feature("dist", dist)

            p = Tree.Tree()
            p.add_child(ch0)
            p.add_child(ch1)

            pname = m0+m1
            nodes[pname] = p

            mtx.NewRow(name=pname)
            mtx.NewCol(name=pname)


            for i in range(len(mtx)):
                mtx[pname,i] = mtx[i,pname] = (mtx[m0,i]+mtx[m1,i])/2

            for m in m0, m1:
                mtx.DelRow(m)
                mtx.DelCol(m)
                del nodes[m]

        return nodes.values()[0]





    def ToNJTreeOld(self):
        """pre:  self is square difference matrix
          post:  A string representing the NJ tree, based on row differences, as calculated by fun """

        def FindNearestRows(mtx):

            rv = 0,1,mtx[0,1]
            if len(mtx) == 2:
                return rv

            for r in range(len(mtx)):
                for c in range(r):
                    el =  mtx[r,c]
                    if el == 0:
                       return r,c,0
                    else:
                        if el < rv[2]:
                            rv = r,c,mtx[r,c]
            return rv

        while len(self) >1:
            m0,m1,v = FindNearestRows(self)
            strv = ":"+str(v)
            self.rnames[m0] = "(" + self.rnames[m0] + strv + "," + self.rnames[m1] + strv+")"
            for i in range(len(self)):
                self[m0,i] = self[i,m0] = (self[m0,i]+self[m1,i])/2

            self.DelRow(m1)
            self.DelCol(m1)

        return self.rnames[0] + ";"


    def IntegiseR(self):
        """pre: Seq.Integise(r) for all r in self.rows
         post: all denominators = 1, ratios unchanged

        !!! Will Fail **SILENTLY** if self[r,c].GetNum() > max(int) !!!

            """

        self.Conv=int
        for r in self.rnames:
            self[r] = map(int,Seq.Integise(self[r]))

    def IntegiseC(self):

        self.Conv=int
        for c in self.cnames:
            self.SetCol(c, map(int,Seq.Integise(self.GetCol(c))))


    def ZapZeroes(self):

        for r in range(len(self)-1, -1, -1):
            if Seq.AllEq(self[r],0):
                self.DelRow(r)

    def Mul(self, other):
        try:
            rows = (Sci.matrix(self.rows) *Sci.matrix(other.rows)).tolist()
            rv = self.__class__(rnames=self.rnames, cnames=other.cnames,Conv=self.Conv)
            rv.rows  = rows
            return rv
        except: # problem in Sci.mul, fall back to slower method
            print "Matrix.Mul() ?",
            selfr = len(self)
            otherr = len(other)
            selfc = len(self.cnames)
            otherc = len(other.cnames)

            if  selfc != otherr:
                raise exceptions.ArithmeticError, "inconsistant matrices for multiplication"

            rv = self.__class__(selfr,otherc,Conv=self.Conv)
            other.Transpose()
            for i in range(selfr):
                for j in range(otherc):
                    rv.rows[i][j] = sum(map(lambda x,y:x*y,self.rows[i],other.rows[j]))
            other.Transpose()
            rv.cnames = other.cnames[:]
            rv.rnames = self.rnames[:]
            return rv


    def OrthNullSpace(self):
        "orthogonal nullspace as matrix of floats"

        try:
            N = Sci.matrix(self.rows,float)
            K = Sci.null(N)
            Z = N*K
            Fail = Z.max() > Sci.TOL or Z.min() < -Sci.TOL
            if Fail:
                print "OrthNullSpace out of tol"
        except:
            Fail = True   # possbility of unknown exceptions from low level math libs
            print "Sci.null() fails !"

        if not Fail:
            rv = matrix(FromMtx=K,Conv=float)
            rv.rnames = self.cnames[:]
            rv.cnames =map(lambda x: "c_"+str(x), range(len(rv[0])))
        else:
            print "falling back to (slow) GaussJ + GramS\n"
            if not self.Conv == Types.ArbRat:
                rv = self.Copy(Conv=Types.ArbRat).NullSpace()
            else:
                rv = self.NullSpace()
            rv.Orthogonalise()

        return rv


    def Orthogonalise(self):
        """ Use Gram-Schmidt orthogonalisation to make self orthognal"""

        self.Transpose() # transposing lets us use more efficient row (aot col) operations
        self.ChangeType(float)

        for j in range(len(self)):
            vj = self[j][:]
            for i in range(j):
                ui = self[i]
                vjui = -Seq.DotProd(vj,ui)
                self.AddRow(j, map(lambda x: x*vjui, ui))
            NormJ = math.sqrt(sum(map(lambda x:x*x,self[j])))
            self.DivRow(j,k=NormJ)
        self.Transpose()


    def PseudoInverse(self,  AsSci=False):

        if self.Conv != float:
            rows = self.Copy(float).rows
        else:
            rows = self.rows

        rvrows = Sci.pinv(rows)
        if AsSci:
            rv = rvrows
        else:
            nrows, ncols = len(rvrows), len(rvrows[0])
            rv = self.__class__(nrows, ncols, Conv=float)
            rv.rows = rvrows.tolist
        return rv


    def SetPivot(self,Row):

        GotOne = 0
        vals=self.GetCol(Row)

        vals = map(abs, vals)[Row:]
        PivVal = max(vals)
        if PivVal != 0:
            GotOne = 1
            PivIdx = Row+vals.index(PivVal)
            if PivIdx != Row:
                self.SwapRow(PivIdx,Row)

        return GotOne


    def ElimCol(self, Row):
        "pre: self.GetEl(Row,Row)==1"

        RowVals = self.rows[Row]
        for r in range(len(self.rows)):
            if r != Row and self.rows[r][Row] != 0:
                mul = -self.rows[r][Row]
                self.AddRow(r, RowVals, mul)


    def RedRowEch(self):
        """ put self in reduced row echelon form,
        uses col swaps in pivoting """

        CurRow = 0
        done = 0
        if self.Conv == float:
            Flatten = Seq.Flatten
        else:
            Flatten = lambda x,y: 0

        while not done:
            if Seq.AllEq(self.rows[CurRow], 0):
                self.DelRow(CurRow)
            else:
                if self.SetPivot(CurRow):
                    self.DivRow(CurRow,vals=None, k=self.rows[CurRow][CurRow])
                    self.ElimCol(CurRow)
                    Flatten(self[CurRow], 1e-6)
                    CurRow += 1
                else:
                    self.MakeLastCol(CurRow)

            done = CurRow == len(self.rows)


    def NullSpace(self):

        rv = self.Copy(Conv=DefaultConv)
        rv.RedRowEch()
        nsrnames = rv.cnames[:]   # because we do col xchs in RedRowEch names get reordered and rv is othrewise mangled
        slice = len(rv.rows)
        if len(rv.cnames) != slice:
            for r in range(len(rv.rows)):
                rv.rows[r] = rv.rows[r][slice:]
            rv.cnames = map(lambda x: "c_"+str(x), range(len(rv[0])))
            id = GetIdent(len(rv.cnames),Conv=self.Conv, mul = -1)
            rv.AugRow(id)
            rv.rnames = nsrnames
            rv.RowReorder(self.cnames) # so now  (self * rv) == 0
        else:
            # special case, invertible matrix, zero dimensional NullSpace
            rv.cnames = []
            rv.rnames = self.cnames[:]
            rv.rows = map(lambda x:[], rv.rnames)
        return rv