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

"""
Classes and functions dealing with enzyme subsets
Main user interface is class EssDic
functions before this are private(ish) helpers for EssDic
functions following are public
"""




# import pickle - No ! pickle **WAY** too slow for large models
import sys
from Util import DynMatrix,Set,Seq,Types
from Structural import StoMat


def EssList(m,k=None):
    """ pre: m = ScrumPy.Model()
                 k = None || k = m.sm.NullSpace()
         post: list subsets, each subset is a list of (name, value) tuples
                  first item is (possibly empty) subset of dead reactions
    """

    if k == None:
        k = m.sm.NullSpace()
    else:
        k = k.Copy()

    k.ZapZeroes()
    One = k.Conv(1)
    deads = Set.Complement(m.sm.cnames, k.rnames)
    ess = [map(lambda x:(x,One), deads)]

    while len(k) !=0:
        ss = [(k.rnames[0],One)]
        row = k[0]
        k.DelRow(0)
        for r in k.rnames[:]:
            ap = Seq.AreParallel(k[r],row)
            if ap[0]:
                ss.append((r,ap[1]))
                k.DelRow(r)
        ess.append(ss)

    return ess


def EssMtx(m,essl = None):
    """ pre: m = ScrumPy.Model()
                 essl == None || EssAsList(m)
       post:  matrix of enz subsets """

    essl = EssList(m)[1:]
    rle = range(len(essl))
    ssnames = map(lambda x: "Ess_"+str(x), rle)
    rnames  = m.sm.cnames[:]
    rv = StoMat.StoMat(rnames = rnames, cnames = ssnames)
    for ssi in rle:                  # for each subset index
        ss = essl[ssi]              # the subset
        for reac,val in ss:
            rv[reac,ssi] = val
    return rv



class EssDic(dict):
    """ Dictionary subclass mapping enzyme subsets to reactions and vice-versa
          pretty much anything useful with subsets should be done via this class
    """

    def __init__(self,m,essl=None):
        """ pre: m = ScrumPy.Model
                     essl = None || essl = EssList(m)
                     (only use essl for development/debugging purposes)
             post: OK. will write messages on sys.stderr if inconsistent irreversibility
                      specifications found, but will not act on this
        """

        if essl == None:
            essl = EssList(m)

        self.Model = m
        self.ReacRevProps = dict(m.sm.RevProps)
        self.RevProps = {}
        self.Rev ={}
        self.Inconsistants = [] # list of subsets with inconsitent irreversible reactions

        self._AddSS("DeadReacs",essl[0])
        n = 1
        for ss in essl[1:]:
            if len(ss)==1:
                name = ss[0][0]
            else:
                name = "Ess_"+str(n)
                n+= 1
            self._AddSS(name,ss)


    def _AddSS(self,key,ss):
        """ Private """

        val = {}
        rev = True    # Reversibility of subset
        for s in ss:
            reac = s[0]
            val[reac] = s[1]
            self.Rev[reac] = key # reverse lookup: reac->ss
            if self.ReacRevProps[reac] != StoMat.t_Rever:
                rev = False

        self[key]= val
        if rev:
            self.RevProps[key] = StoMat.t_Rever
        else:
            self.RevProps[key] = StoMat.t_Irrev
            self.RevCheck(key)




    def RevCheck(self,  key):

        def WrongWay(ss):
            for reac in ss:
                coeff = ss[reac]
                revcrit = self.ReacRevProps[reac]
                if (coeff < 0 and revcrit == StoMat.t_Irrev) or (
                    coeff > 0 and revcrit == StoMat.t_BackIrrev):
                       return True
            return False


        ss = self[key]

        if WrongWay(ss):
            for reac in ss.keys():
                ss[reac] *= -1

        if WrongWay(ss):
            sys.stderr.write("!\n! Inconsistent irreversible reactions in subset "+ key + "\n")
            sys.stderr.write("! " + str(ss.keys()) + "\n!\n")
            self.Inconsistants.append(key)


    def RevCheckOld(self, key):
        """ pre: self.has_key(key) (privateish)
            post: checks for irreversible reactions in self[key],
                     ensures correct sign  of coeffs, and marks as irrev  if irreversible,
                     reports inconsistent irreversbility in self[key]
        """

        subset = self[key]
        irrevs = []
        negate = False

        for reac in subset.keys():
            if reac in self.Model.sm.Irrevs:
                irrevs.append(reac)
                if  subset[reac] <0:
                    negate = True

        if negate:
            for reac in subset.keys():
                subset[reac] *= -1

            error = False
            for reac in  irrevs:
                if subset[reac] < 0:
                    error = True

            if error:
                sys.stderr.write("!\n! Inconsistent irreversible reactions in subset "+ key + "\n")
                sys.stderr.write("! " + str(subset.keys()) + "\n!\n")
                self.Inconsistants.append(key)
                del self.Irrevs[key]





    def ReacToSS(self,reac):
        """ pre: reac present in the model
            post: name of subset of which reac is a member
        """
        return self.Rev[reac]



    def SameSS(self,reac):
        """ pre: self.ReacToSS(reac)
            post: return the names of reactions in the same subset (including reac)
        """
        return self[self.ReacToSS(reac)].keys()




    def ToMtx(self,ExclDeadIncs=False):
        """ pre: True
            post: subsets in matrix form, ExclDeadIncs => don't include dead or inconsistant subsets in the matrix
        """

        rnames  = self.Model.sm.cnames[:]
        ssnames = self.keys()
        if ExclDeadIncs:
            for i in self.Inconsistants+["DeadReacs"]:
                ssnames.remove(i)

        rv = StoMat.StoMat(rnames = rnames, cnames =ssnames)
        for ssn in ssnames:
            ss = self[ssn]
            for reac in ss.keys():
                rv[reac,ssn] = self[ssn][reac]
        rv.RevProps = dict(self.ReacRevProps)

        return rv


    def CondensedSMs(self):
        """ pre: True
           post: (self.ToMtx(), Condensed Internal SM, Condensed External SM) """


        ssm = self.ToMtx(True)
        m = self.Model
        sm = m.sm.Mul(ssm)
        sm.ZapZeroes()

        sm.RevProps = dict(self.RevProps)
        smx = m.smexterns.Mul(ssm)
        smx.ZapZeroes()
        smx.RevProps = dict(self.RevProps)

        return ssm, sm, smx



    def SaveCondensedModel(self,fname=None):
        """ pre: fname != None => fname can be opened for writing
                     else:  current working directory is writable
            post: current model is in condensed form saved as a ScrumPy file
                     with a cPickled self.
                     as fname_Condensed and fname_Condensed_esd.pic
                     if fname was None model.fname will be used
        """

        m = self.Model

        mfname = fname or m.FileName + "_Condensed"
        out = open(mfname,"w")
        print "Saving condensed model", mfname

        smx = m.smexterns.Mul(self.ToMtx(True))


        smx.RevProps = self.RevProps()
        for c in range(len(smx.cnames)):
            cname = smx.cnames[c]
            if len(self[cname])>1:
                reac = "_".join([cname]+self[cname].keys())
                smx.RevProps[reac] = self.RevProps[cname]
                smx.cnames[c] = reac

        smx.ToScrumPy(out,Externs=self.Model.Externals())



    def ToList(self):

        keys = self.keys()
        keys.remove("DeadReacs")
        rv = [self["DeadReacs"].items()]
        for k in keys:
            rv.append(map(lambda x:x[0],self[k].items()))
        return rv

    def SaveFile(self, FileName=None):
        """ pre: FileName != None => open(FileName,"w") will succeed
           post: self saved as Pickle to FileName if given
                    or to self.Model.FileName +" _esd.pic" if FileName != None"""

        print "Don't use 'SaveFile()'"
        return

        if FileName == None:
            FileName = self.Model.FileName + "_esd.pic"

        model = self.Model
        self.Model = self.Model.FileName

        for v in self.values():
            for k in v.keys():
                v[k] = str(v[k])

        out = open(FileName,"w")
        pickle.dump(self,out)

        self.Model = model

        for v in self.values():
            for k in v.keys():
                v[k] = Types.ArbRat(v[k])



    def IsIrrev(self,ss):
        """ pre: True
           post: self.IsIrrev(ss) => ss is irreversible
        """
        return  self.RevProps[ss] != StoMat.t_Rever


    def Irrev_ss(self):
        """ pre: True
            post: list of irreversible subsets
        """
        return map(lambda reac:self.RevProps[reac],  filter(self.IsIrrev,  self.keys()))



def EsdFromFile(FileName, LoadModel=False):
    """ pre: Esd.SaveFile(FileName)
        post: loads and returns previously saved esd
                LoadModel => model is loaded too
    """

    inp = open(FileName)
    rv = pickle.load(inp)


    for v in rv.values():
        for k in v.keys():
            v[k] = Types.ArbRat(v[k])


    if LoadModel:
        rv.Model = ScrumPy.Model(rv.Model)
    return rv