/home/mark/model/software/ScrumPy/ScrumPy/Structural/Model.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,  exceptions


from Util import Seq, DynMatrix, Set, Types,Tree
from Structural import StoMat, EnzSubset, ReacTree
from Parser import lexer
import ElModes

BadElMsg="""
!! Element types of stoichiometry matrix !!
!!      *MUST* be arbitrary rational     !!
!!              Giving up.               !!
"""

def SMFromDesc(Desc):

    ElType =  Desc.OptDic["ElType"]
    StoDic = Desc.StoDic
    MetList =Desc.MetList

    sm = StoMat.StoMat(
        cnames=StoDic.keys(),
        rnames=MetList[:],
        Conv=ElType
    )

    for reac in sm.cnames:
        sto = StoDic[reac]
        sm.SetSto(reac, sto)

    return sm


class Model:

    def Destroy(self):
        pass

    def BuildStructural(self):
        """pre: ModelDesc is correctly filled in instance of ScrumPy.Parser.Parser.ModelDescriptor """

        ModelDesc = self.ModelDesc

        self.smexterns = SMFromDesc(ModelDesc)
        ExtList = ModelDesc.ExtList[:]
        if ModelDesc.OptDic["AutoExtern"]:
            ExtList = Set.Union(ExtList,self.smexterns.Unbals())

        self.sm = self.smexterns.Copy()
        for e in ExtList:
            try:
                self.sm.DelRow(e)
            except:
                sys.stderr.write("! Warning: "+e+" declared external, but not used !\n")

        self.sm.Externs =  ExtList[:]
        self.smexterns.Externs = ExtList[:]
        if ModelDesc.OptDic["DeQuote"]:
            self.sm.DeQuote()
            self.smexterns.DeQuote()

        self.sm.Irrevs = self.sm.GetIrrevs()
        self.smexterns.Irrevs = self.sm.GetIrrevs()

        self.RemoveUnused()





    def DelReactions(self,Reactions):
        """ pre: Reactions is a list of strings
          post: reactions removed from self """
        for r in Reactions:
            if not r in self.sm.cnames:
                sys.stderr.write("Warning: tried to delete non-existent reaction, " + r + "\n")
            else:
                self.sm.DelReac(r)
                self.smexterns.DelReac(r)


        self.sm.ZapZeroes()
        self.smexterns.ZapZeroes()
        self.smexterns.Externs = Set.Intersect(self.smexterns.Externs, self.smexterns.rnames)
        self.sm.Externs = self.smexterns.Externs[:]


    def AddReactions(self, StoDicDic):
        print "ScrumPy.Structural.Model.AddReactions due for removal"
        for reac in StoDicDic.keys():
           sto = StoDicDic[reac].copy()
           self.smexterns.NewReaction(reac, sto)
           for met in sto.keys():
               if met in self.sm.Externs:
                   del(sto[met])
           self.sm.NewReaction(reac, sto)




    def RemoveUnused(self):

        Unused = Set.Complement(self.smexterns.cnames, self.sm.cnames)

        if len(Unused) >0:
            print "The following reactions affect no internal metabolites, removing them:\n",Unused
        for u in Unused:
            if u in self.smexterns.cnames:
                self.smexterns.DelCol(u)
            if u in self.sm.cnames:
                self.sm.DelCol(u)



    def IsIrrev(self, name):
        """pre: name is a reaction name
          post: self.IsIrrev(name) => name is an irreversible reaction """

        return self.sm.IsIrrev(name)


    def NetSto(self,vels):
        "return a matrix of the net stoichs represented by the rate information in vels"

        rv = DynMatrix.matrix()
        sm = self.smexterns.ToNumPyMtx(conv=int)
        v = vels.ToNumPyMtx(conv=int)
        res = sm * v
        rv.FromNumPyMtx(res)
        rv.rnames = self.smexterns.rnames
        rv.cnames = vels.cnames
        return rv


    def Externals(self):
        return Set.Complement(self.smexterns.rnames, self.sm.rnames)

    def Transporters(self,NamesOnly=False):
        seen = {}
        rv = []
        exs = self.Externals()
        for e in exs:
            for tx in self.smexterns.InvolvedWith(e):
                if not seen.has_key(tx):
                    seen[tx] = 1
                    rv.append([tx, Set.Intersect(exs,self.smexterns.Reactants(tx))])

        if NamesOnly:
            return seen.keys()
            return rv


    def FindIsoforms(self):
        stodic = {}
        for reac in self.smexterns.cnames:
            strsto = str(self.smexterns.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 DelIsoforms(self):
        isos = self.FindIsoforms()
        rv = {}
        for iso in isos:
            rv[iso[0]] = iso[1:]
            for reac in iso[1:]:
                self.sm.DelReac(reac)
                self.smexterns.DelReac(reac)
        return rv


    def FindOrphans(self):
        raise DeprecationWarning,  "use model.OrphanMets() instead"


    def OrphanMets(self):
        return self.smexterns.OrphanMets()


    def Core(self):
        """ pre: True
           post: A copy of the external stoichiometry matrix with all internal orphans
                 and associated reactions removed. """

        if self.smexterns.IsDeQuoted():
            self.smexterns.ReQuote()
            smx = self.Core()
            self.smexterns.DeQuote()
            smx.DeQuote()
        else:
            smx = self.smexterns.Copy()
            orphans  = Set.Complement(smx.OrphanMets(), smx.Externs)
            while len(orphans) > 0:
                print "len(o)", len(orphans)
                for met in orphans:
                    if not met in smx.Externs:
                        for rn in smx.InvolvedWith(met):
                            smx.DelReac(rn)
                        smx.DelRow(met)

                smx.ZapZeroes()
                orphans  = Set.Complement(smx.OrphanMets(), smx.Externs)

        return smx


    def DeadReactions(self):
        rv  = []
        k = self.sm.NullSpace()
        for r in k.rnames:
            if Seq.AllEq(k[r],0):
                rv.append(r)
        return rv

    def DisconnectsK(self, thresh=1e-6,exts=False):
        """ pre: True
           post: return list of lists of disconnected subsystems by null-space analysis of self.sm
                    exts => use the external stoichiometry matrix
                    first item is possible empty list of dead reactions
        """
        if exts:
            k = self.smexterns.OrthNullSpace()
        else:
            k = self.sm.OrthNullSpace()
        deads = []
        for reac in k.rnames[:]:
            if Seq.AllMatch(k[reac], lambda x:abs(x)<thresh):
                deads.append(reac)
                k.DelRow(reac)
        return k.NonZeroR([deads],thresh)


    def MaxCycles(self):
        return self.DisconnectsK(exts=1)[1:]





    def EnzSubsets(self):
        """ pre: True
            post: return EnzSubset.EssDic(self)
                     see EnzSubset module for details
        """

        return EnzSubset.EssDic(self)



    def ElModes(self):

        ssm, sm, smx = self.EnzSubsets().CondensedSMs()
        if len(ssm.cnames) <2: # 0 or 1 ssubsets means 0 or 1 ems
            elmodes = ssm
            stos = smx
        else:
            ems = ElModes.GetElModes(sm)
            #ems.Transpose()
            elmodes = ssm.Mul(ems)

        return ElModes.ModesDB(elmodes,self)



    def ConsMoieties(self,raise_err=True):

        smt = self.sm.Copy(tx = 1)
        ksmt = smt.NullSpace()
        deps = []
        for c in ksmt.cnames:
            for r in ksmt.rnames:
                if ksmt[r,c] != 0 and Seq.NumOfMatches(ksmt[r], lambda x: x!=0) ==1:
                    deps.append(r)
                    break

        if len(deps) != len(ksmt.cnames):
           print "! missed deps should have", len(smtk.cnames), "found", len(deps)
           if raise_er:
               raise exceptions.ValueError, str(len(smtk.cnames)-len(deps)) + "missing dependencies"

        for d in deps:
            ksmt.MakeFirstRow(d)
            c = deps.index(d)
            val = ksmt[0,c]
            ksmt.cnames[c] = d
            ksmt.DivCol(c,k=val)


        return ksmt



##
##    def ConsMoieties(self):
##        self.sm.Transpose()
##        rv  = self.sm.NullSpace()
##        if rv != None:
##            rv.RowReorder(self.sm.cnames)
##            rv.Transpose()
##            for r in range(len(rv)):
##                rv.rnames[r] = "CONS_"+str(r)
##                if  Seq.NumOfMatches(rv[r], lambda x: x<0) > Seq.NumOfMatches(rv[r], lambda x: x>0):# more -ve than +ve elements
##                    rv.MulRow(r, k = -1)
##
##            self.sm.Transpose()
##
##        return rv
##
    def ReacTree(self):

        def diff(a,b):
            return 1- abs(Seq.CosTheta(a,b))

        k = self.sm.OrthNullSpace()
        rdm = k.RowDiffMtx(diff)
        t= rdm.ToNJTree()
        t.rdm = rdm
        return t