/home/mark/model/software/ScrumPy/ScrumPy/Structural/StoMat.py

"""

ScrumPy -- Metabolic Modelling with Python

Copyright Mark Poolman 1995 - 2002

 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

"""

import sys, types

from Util import DynMatrix, Seq,Set,Sci,Format, Types
# NB "types" is the standard module, "Types" contains some local type defs

from Parser.lexer import t_Rever, t_Irrev, t_BackIrrev
from Parser.Parser import Sto
# TODO: Parser.Sto and ModelDesc need a module of their own

# Exceptions
NoSuchMet = "No such metabolite"
NoSuchreac = "No such reaction"
DupReac = "Duplicate reaction"

class StoMat(DynMatrix.matrix):

    def __init__(self, r=0, c=0, **kwargs):

        DynMatrix.matrix.__init__(self,r,c, **kwargs)
        self.RevProps = {}
        self.Externs = []
        #self.Irrevs = self.GetIrrevs()


        # TODO: StoMat needs a proper copy constructer


    def DelRow(self, r):

        if type(r) == types.IntType:
            r = self.rnames[r]

        if r in self.Externs:
            self.Externs.remove(r)

        DynMatrix.matrix.DelRow(self,r)


    def DelReac(self, c):

        if type(c) == types.IntType:
            c = self.cnames[c]

        del self.RevProps[c]
        DynMatrix.matrix.DelCol(self,c)


    def DeQuote(self):

        DynMatrix.matrix.DeQuote(self)

        revpropdq = {}   # reversibility property dictionary dequited
        for k in self.RevProps:
            revpropdq[self.QuoteD[k]] = self.RevProps[k]
        self.RevProps = revpropdq

        xnamesd = map(Format.StripQuote,self.Externs) # have to externs by hand
        for n in range(len(xnamesd)):                 # as internal sms won't otherwise
            self.QuoteD[xnamesd[n]] = self.Externs[n] # have Externs in QuoteD
            self.QuoteD[self.Externs[n]] = xnamesd[n]
        self.Externs = xnamesd


    def ReQuote(self):

        revpropdq = {}
        for k in self.RevProps:
            revpropdq[self.QuoteD[k]] = self.RevProps[k]
        self.RevProps = revpropdq
        self.Externs = map(lambda x:self.QuoteD[x], self.Externs)
        DynMatrix.matrix.ReQuote(self)



    def MakeExtern(self, name):

        self.DelRow(name)
        self.Externs.append(name)

    def GetIrrevs(self):
       return filter(lambda reac: self.RevProps[reac] != t_Rever,  self.cnames)

    def MakeIrrev(self,name):

        self.RevProps[name] = t_Irrev


    def MakeRevers(self,name):

        self.RevProps[name] = t_Rever


    def InvolvedWith(self, name, crit=None):
        if crit==None:
            def crit(x):
                return x!=0


        if name in  self.cnames:
            list = self.GetCol(name)
            namelist = self.rnames
        elif name in self.rnames:
            list = self.GetRow(name)
            namelist = self.cnames
        else:
            print "name ", name, "not found"
            list = []

        rv = {}
        for idx in range(len(list)):

            if crit( list[idx]):
                rv[namelist[idx]] = self.Conv(list[idx])

        return (rv)


    def Products(self, Reaction):
        return self.InvolvedWith(Reaction, lambda x: x>0).keys()

    def Substrates(self, Reaction):
        return self.InvolvedWith(Reaction, lambda x: x<0).keys()

    def Reactants(self, Reaction):
        """ List of reactant names involved with Reaction """
        return self.InvolvedWith(Reaction).keys()

    def Transporters(self):
        rv = {}
        for met in self.Externs:
            rv.update(self.InvolvedWith(met))
        return rv.keys()


    def OrphanMets(self):

        rv = []

        for met in self.rnames:
            if len(self.InvolvedWith(met))==1:
                rv.append(met)

        return rv


    def Delete(self, name):
        if name in self.cnames:
            self.DelReac(name)
            for name in self.rnames:                 # ugly, find nicer way
                if Seq.AllEq(self.GetRow(name), 0):
                    self.DelRow(name)
        elif name in self.rnames:
            self.DelRow(name)
            for name in self.cnames:
                if Seq.AllEq(self.GetCol(name), 0):
                    self.DelReac(name)
        else:
            print "can't find ", name


    def NewReaction(self, reac, StoDic,RevProp=t_Rever):

        if reac in self.cnames:
            raise DupReac

        self.NewCol(name=reac)
        for met in StoDic.keys():
            if not met in self.rnames:
                self.NewRow(name=met)
            self[met,reac]=StoDic[met]

        self.RevProps[reac] = RevProp


    def SetSto(self, reac,  Sto,  reset = False):
        """ pre: reac in self.cnames, Sto = Parser.Sto(...)"""

        if reset:
            self.MulCol(reac, k=0)

        for met in Sto.lhs:
            coeff = met[0]
            name = met[1]
            self[name,reac] = coeff
            self[name,reac] *= -1 # cos coeff is a string until it gets into self

        for met in Sto.rhs:
            coeff = met[0]
            name = met[1]
            self[name,reac] = coeff

        self.RevProps[reac] = Sto.Revers



    def FirstReaction(self, Name, lhs, rhs, Revers):

        print "Deprecated FirstReaction() ?"
        self.cnames.append(Name)
        for tup in lhs:
            if tup[0] in self.rnames:
                 raise "Duplicate Metabolite"
            self.rnames.append(tup[0])
            self.rows.append([self.Conv(-tup[1])])

        for tup in rhs:
            if tup[0] in self.rnames:
                 raise "DuplicateMetabolite"
            self.rnames.append(tup[0])
            self.rows.append([self.Conv(tup[1])])

        if not Revers:
                self.MakeIrrev(Name)


    def AddReaction(self, Name, lhs, rhs, Revers=1):

        sys.stderr.write("Deprecated: AddReaction - Use NewReaction\n")
        raise "Deprecated method", "AddReaction"


    def Copy(self, *args,**kwargs):

        rv =  DynMatrix.matrix.Copy(self,*args, **kwargs)
        rv.RevProps = dict(self.RevProps)

        rv.Externs = self.Externs[:]

        return rv


    def ReacToStr(self,reac):

        def d2l(d):
            rv = []
            for k in d.keys():
                rv.append(" ".join((str(d[k]),k)))
            return rv

        substs = self.InvolvedWith(reac, lambda x: x<0)
        for k in substs.keys():
            substs[k] *= -1
        substs = d2l(substs)


        prods = d2l(self.InvolvedWith(reac, lambda x: x>0))

        sep = self.RevProps[reac].join((" "," "))


        rv = reac+":\n\t" + " + ".join(substs) + sep + " + ".join(prods) + "\n\t~\n"

        if len(substs) ==0 or len(prods)==0:
            print reac, "missing substrate or product !"
            rv = "#"+ "\n#".join(rv.split("\n"))+"\n"

        return rv

    def IntegiseReacs(self):
        for reac in self.cnames:
            col = self.GetCol(reac)
            if Seq.AllEq(col, 0):
                self.DelCol(reac)
            else:
                #Seq.Normalise(col,KeepSign=False)
                self.SetCol(reac, Seq.Integise(col))


    def NiceOrder(self):

        txs = self.Transporters()
        interns = Set.Complement(self.cnames,txs)
        txs.sort()
        interns.sort()
        self.ColReorder(txs+interns)


    def ElTypeAsScrumPy(self):
       
        return Types.TypeToDirective[type(self.Conv(1))]




    def ToScrumPy(self, file=None, FirstLine="Structural()\n",Externs=None, NiceOrder=False):

        if self.IsDeQuoted():
            FirstLine += "DeQuote()\n"
            self.ReQuote()
            rv = self.ToScrumPy(file, FirstLine,Externs, NiceOrder)
            self.DeQuote()
            return rv

        if Externs==None:
            Externs=self.Externs[:]

        for ex in Externs[:]:
            if ex.upper()[:2] == "X_":
                Externs.remove(ex)

        self2 = self.Copy()

        self2.ZapZeroes()

        if self.Conv != Types.ArbRat:
            FirstLine += "ElType("+self.ElTypeAsScrumPy() + ")\n"
        if NiceOrder:
            self2.NiceOrder()

        Header = FirstLine
        Body = "\n"
        if len(Externs) >0:
            Header += "External("+ ", ".join(Externs)+")\n\n"
        for reac in self.cnames:
            if len(self2.InvolvedWith(reac))==0:
                Header += "# Empty Reaction : "+ reac + "\n"
            else:
                Body += self2.ReacToStr(reac)

        rv = Header + Body
        if file != None:
            if not hasattr(file,"write"):
                file = open(file,"w")
            file.write(rv)
        else:
            return rv




    def IsIrrev(self,name):
        return  self.RevProps[name] != t_Rever

    def Connectedness(self, met):
        return Seq.NumOfMatches(self[met], lambda x: x != 0)

    def Unbals(self):
        " return list of mets that do not have >0 producer and >0 consumer "

        rv = []

        for m in self.rnames:
            row = self[m]
            if not (Seq.HasMatches(row, lambda x: x>0) and Seq.HasMatches(row, lambda x: x<0)):
                rv.append(m)
        return rv


    def SubSys(self, reacs):

        if self.IsDeQuoted():
            self.ReQuote()
            reacsq  = []
            for reac in reacs:
                reacq = '"'+reac+'"'
                if reacq in self.cnames:
                    reacsq.append(reacq)
                else:
                    reacsq.append(reac)

            rv = self.SubSys(reacsq)
            rv.DeQuote()
            self.DeQuote()

        else:

            rv = StoMat(rnames=self.rnames, cnames=reacs,Conv=self.Conv)

            for reac in reacs:
                stod = self.InvolvedWith(reac)
                for met in stod:
                    rv[met,reac] = stod[met]

            rv.ZapZeroes()
            rv.Externs = Set.Intersect(rv.rnames,self.Externs)
            for reac in rv.cnames:
                rv.RevProps[reac] = self.RevProps[reac]

            # TODO - rnamesq etc. no longer needed ? check and remove

            if self.IsDeQuoted():
                for atr in "rnamesq",  "cnamesq", "Externsq", "Irrevsq":
                    setattr(rv, atr, map(self.FindQuoted, getattr(self, atr[:-1])))


        return rv




    def FindIsoforms(self):
        stodic = {}
        for reac in self.cnames:
            strsto = str(self.InvolvedWith(reac))
            if stodic.has_key(strsto):
                stodic[strsto].append(reac)
            else:
                stodic[strsto] = [reac]
        rv = []
        for reac in stodic.values():
            if len(reac)>1:
                rv.append(reac)

        return rv




    def OrthSubSysts(self,tol=1e-6):
        """ a list of lists of reaction names of the orthogonal subsystems in self,
        The first, possibly empty, list contains the dead reactions """

        k = self.OrthNullSpace()
        dead = []
        for reac in k.rnames[:]:
            if max(map(abs, k[reac])) <tol:
                dead.append(reac)
                k.DelRow(reac)

        rv = [dead]
                
        while len(k.rnames) >0:
           ref = k.rnames[0]
           curss = [ref]

           for reac in k.rnames[1:]:
               if not Sci.AreOrthogonal(k[ref],k[reac],tol):
                   curss.append(reac)
                   k.DelRow(reac)

           k.DelRow(ref)
           rv.append(curss)
        return rv


    def ReactionNJTree(self):
        k = self.OrthNullSpace()
        t = k.RowDiffMtx(Sci.Theta).ToNJTree()
        t.Consolidate(thresh= 1e-8)
        return t

    def AdjMtx(self):
        """ The integer adjacency matrix of self """

        rv = DynMatrix.matrix(Conv = int, rnames = self.rnames, cnames = self.cnames)
        for met in self.rnames:
            reacs = self.InvolvedWith(met).keys()
            for r in reacs:
                for neighbour in self.InvolvedWith(r).keys():
                    rv[r,neighbour] = 1
        return rv