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

import StoMat, Decompose
from clib.ElModes import ElModes
from Util import   Set, Seq, Sorter,  Types, DynMatrix
from Data import Histo



neg = lambda x: x<0
pos = lambda x: x>0

def deb(msg):
    sys.__stderr__.write("ElModes.py" + msg + "\n")

class ModesDB:      # a data base type thing, to allow searching of elementary modes etc.

    def __init__(self, tabl,srcmodel, sto=None):

        self.model = srcmodel
        if tabl.__class__ == Tableau:
            self.mo = tabl.AsStoMat()
        else:
            self.mo = tabl

        if sto != None:
            self.sto = sto
        else:
            self.sto =  srcmodel.smexterns.Mul(self.mo)

            self.sto.RevProps = self.mo.RevProps = RevProps = {}
            for e in self.mo.cnames:
                reacs = self.mo.InvolvedWith(e).keys()
                RevProps[e] = StoMat.t_Rever
                for r in reacs:
                    if not srcmodel.sm.RevProps[r] == StoMat.t_Rever:
                        RevProps[e] = StoMat.t_Irrev
                        break






    def _retmosto(self):

        mo  = StoMat.StoMat(len(self.mo))
        mo.rnames = self.mo.rnames[:]
        return mo


    def __len__(self):
        if len(self.mo) >0:
            return len(self.mo[0])
        return 0

    def __add__(self, other):
        mo = self.mo.Copy()
        sto = self.sto.Copy()

        mo.AugRow(other.mo)
        sto.AugRow(other.sto)

        return ModesDB(mo,sto,self.model,False)

    def NormByName(self,names):
        """ Normalise all modes and stos such that they consume 1 mol of the first
            species in  names that they utilise """

        for mo in self.mo.rnames:
            for name in names:
                coeff = self.sto[mo,name]
                if coeff <0:
                    self.mo.DivRow(mo, k=-coeff)
                    self.sto.DivRow(mo, k=-coeff)
                    break





    def Copy(self):
        return ModesDB(self.mo.Copy(), self.sto.Copy(), self.model,False)

    def AddMode(self,StoDic,MoDic,Name=None):

        if Name == None:
            Name = "ElMode_"+str(len(self))

        for dic,mat in (StoDic, self.sto),(MoDic,self.mo):
            mat.NewRow(name=Name)
            for k in dic.keys():
                if not k in mat.cnames:
                    mat.NewCol(name=k)
                mat[Name,k] = dic[k]




    def DelMode(self,mode):
        self.mo.DelRow(mode)
        self.sto.DelRow(mode)

    def Stos(self):
        "return a string repesentation of stoichiometries"
        return "".join(map(self.sto.ReacToStr, self.sto.cnames))

    #def Modes(self):
    #    "return a string representation of the elementary modes"
    #    return self.mo.ToStr()

    def Modes(self):
        retlist = []
        for e in self.mo.cnames:
            molist = [e]
            rd = self.mo.InvolvedWith(e)
            for r in rd:
                molist.append(" ".join((str(rd[r]),r)))
            retlist.append(", ".join(molist))
        return "\n\n".join(retlist)


    def LenHist(self,nbins=20):
        """ histogram of mode lengths. """

        lens  = map(lambda em: self.LenOf(em), self.mo.cnames)
        return Histo.Histo(lens,nbins=nbins)



    def _rowsearch(self, name, crit, mtx, neg=0):

        if neg:
            return self._rowsearch(name, lambda x,fun=crit: not fun(x), mtx)

        matches = mtx.InvolvedWith(name, crit).keys()
        mo = self.mo.SubSys(matches)
        sto = self.sto.SubSys(matches)

        return ModesDB(mo, self.model, sto)


    def Consumes(self, name, neg=0):
        return self._rowsearch(name, lambda x: x<0, self.sto, neg)

    def Produces(self,name, neg=0):
        return self._rowsearch(name, lambda x: x>0, self.sto, neg)

    def Used(self, name, neg=0):
        return self._rowsearch(name, lambda x: x!=0, self.sto, neg)

    def Unused(self, name, neg=0):
        return self._rowsearch(name, lambda x: x==0, self.sto, neg)

    def PosFlux(self,name, neg=0):
        return self._rowsearch(name, lambda x: x>0, self.mo, neg)

    def NegFlux(self, name, neg=0):
        return self._rowsearch(name, lambda x: x<0, self.mo, neg)

    def HasFlux(self, name, neg=0):
        return self._rowsearch(name, lambda x: x!=0, self.mo, neg)

    def NoFlux(self, name, neg=0):
        return self._rowsearch(name, lambda x: x==0, self.mo, neg)

    def PartSto(self, names):
        rv = DynMatrix.matrix()
        for name in names:
            rv.NewCol(self.sto.GetCol(name), name)
        rv.rnames = self.sto.rnames

    def ReacsOf(self, m):
        return self.mo.InvolvedWith(m)


    def ConsumesAnyOf(self, names):

        rvd = {}
        for name in names:
            rvd.update( self.sto.InvolvedWith(name, neg))

        ems = rvd.keys()
        mo = self.mo.SubSys(ems)
        sto = self.sto.SubSys(ems)

        return ModesDB(mo,  self.model, sto)

    def ConsumesOnly(self, names):

        rv = self.ConsumesAnyOf(names)
        for em in rv.sto.cnames[:]:
            if len(rv.sto.InvolvedWith(em, neg))!=1:
                rv.mo.DelCol(em)
                rv.sto.DelCol(em)

        return rv



    def ProducesAnyOf(self, names):
        rvd = {}
        for name in names:
            rvd.update( self.sto.InvolvedWith(name, pos))

        ems = rvd.keys()
        mo = self.mo.SubSys(ems)
        sto = self.sto.SubSys(ems)

        return ModesDB(mo,  self.model, sto)

    def ProducesOnly(self, names):

        rv = self.ProducesAnyOf(names)
        for em in rv.sto.cnames[:]:
            if len(rv.sto.InvolvedWith(em, pos))!=1:
                rv.mo.DelCol(em)
                rv.sto.DelCol(em)

        return rv


    def Integise(self):
        """ convert rational -> integer representation"""

        self.sto.IntegiseC()
        self.mo.IntegiseC()

        #Integise = lambda s: map(int,  Seq.Integise(s))
        #if self.sto.Conv == Types.ArbRat:
        #    self.sto.rows = map(Integise,  self.sto.rows)
        #    self.sto.Conv = int
        #    self.mo.rows = map(Integise,  self.mo.rows)
        #    self.mo.Conv = int
        #else:
        #    print >>sys.__stderr__,  "Can only integise if matrix elements are rational - ignoring."



    def Futile(self):
        """Strictly internal cycles - ie no net stoichiometry"""

        ems = []
        for e in self.sto.cnames:
            if Seq.AllEq(self.sto.GetCol(e),0):
                ems.append(e)
        mo = self.mo.SubSys(ems)
        sto = self.sto.SubSys(ems)
        return ModesDB(mo,  self.model,  sto)

    def GetDupStos(self):
        return self.sto.DupRows()


    def DelDupStos(self):
        """ remove elmodes with identical stoichs """

        for dups in self.GetDupStos():
            for d in dups[1:]:
                self.sto.DelRow(d)
                self.mo.DelRow(d)

    def MergeDupStos(self):
        """ merge elmodes with duplicate stoichs to form non-elementary modes """

        dups = self.GetDupStos()
        for dup in dups:
            for d in dup[1:]:
                self.sto.DelRow(d)
            self.mo.AddAndRemoveRows(dup)




    def Decompose(self,DataSet,**opts):

        return Decompose.Decompose(self, DataSet, **opts)


    def StoOf(self, ElMo):
        """ returns a string repn of the net Sto of ElMo"""

        return self.sto.ReacToStr(ElMo)



    def LenOf(self, mo):
        """ The number of reactions in mode, mo """

        return  len(filter(lambda x: x!=0,self.mo.GetCol(mo)))


    def Shortest(self):
        rv = []
        best = len(self.mo[0])
        for mode in self.mo.rnames:
            length = self.LenOf(mode)
            if length < best:
                rv = [mode]
                best = length
            elif length == best:
                rv.append(mode)
        return rv



    def Irrevs(self):
        """ pre: true:
           post: e.Irrves() == a posibly empty list of irreverible mode names
        """

        return filter(lambda em:self.IsIrrev(em),  self.mo.cnames)


    def IsIrrev(self, em):

        return self.mo.RevProps[em] == "->"









class Tableau(StoMat.StoMat): #DynMatrix.matrix):

    def __init__(self, sm):
        """pre: sm = StoMat.StoMat(...) """

        #StoMat.StoMat.__init__(self, Conv=Types.ArbRat, FromMtx=sm)
        ## FromMtx will over-ride Conv, hence we do it as a separate step. The C module
        ## *MUST* have ArbRats as ElType.
        StoMat.StoMat.__init__(self, FromMtx=sm)
        if self.Conv != Types.ArbRat:
            self.ChangeType(Types.ArbRat)
        self.RevProps.update(sm.RevProps)
        self.Transpose()
        # start with a transposed copy of the stoichiometry matrix

        id = DynMatrix.GetIdent(len(sm.cnames))           # get ident matrix of same size
        id.cnames = sm.cnames[:]                          # ident col names will be reaction names
        self.AugCol(id)                                   # augment transposed stomat with ident

        self.nMets = len(sm.rnames)                       # Make data available in convenient form
        self.RowLen = len(self.rows[0])                   # for the C library
        self.IrrevIdxs = []                               # indices of irreversible reactions in self
        for name in sm.GetIrrevs():                          # from those in the original sm
            self.IrrevIdxs.append(self.rnames.index(name)+self.nMets)
            self.IrrevLen = len(self.IrrevIdxs)
        self.Results = []                                   # where the C lib will deposit the results




    def GetElModes(self):

        ElModes.ElModes(self)  # ElModes is the low-level C entity

        self.cnames = []
        self.rows=[]
        self.RevProps = {}

        print "" , # very wierd - Idle ? - next line sometimes excepts without *some* write to stdout/stderr
        for i in self.rnames:
            self.rows.append([])

        for r in range(len(self.Results)):
            cname = "ElMo_" + str(r)
            self.NewCol(map(
                self.Conv,
                self.Results[r][self.nMets:]
              ),
              cname
            )






    def AsStoMat(self):
        rv =StoMat.StoMat(FromMtx=self)
        rv.RevProps=dict(self.RevProps)
        return rv


def GetElModes(sm):
    t = Tableau(sm)
    t.GetElModes()
    return t


def GetModesDB(m):

    t = GetElModes(m.sm)
    return ModesDB(t, m)




def FancyGetElModes(sm):


    #def ColSort(Col, mtx):
    #    return Seq.Angle(mtx.GetCol(0), mtx.GetCol(Col))

    def ModSort(row, mtx):
        return Seq.Mod(mtx[row])
    sort = Sorter.Sorter()
    for r in sm.rnames:
        sort.append(r)
    sorted = sort.Sort(ModSort,mtx=sm)
    print sorted
    sm.RowReorder(sorted)
    t = Tableau(sm)

    t.GetElModes()

    return t