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

import sys
from ScrumPy.Util.Tree import Tree
from ScrumPy.Util import Set,Sci,Seq
import StoMat


def ListToName(l):
    l.sort()
    return "_".join(l)

class Smholder:
    def __init__(self, sm, xsm):
        self.sm = sm
        self.smexterns = xsm
        

class PartialMode:
    def __init__(self, sm,Balances = [],Uses = [],Irrev=False):
        
        names = sm.cnames[:]
        names.sort()
        self.Name = "_".join(names)
        self.sm = sm
        self.Balances = Balances
        self.Irrev = Irrev
        self.Uses = Uses
        

    def DefinedSto(self):
        return len(self.sm.cnames)==1

    def __str__(self):
        return ", ".join(
            (
                self.Name,
                str(self.Substrates()),
                str(self.Products())
                )
        )

    def __repr__(self):
        return self.Name

    def Substrates(self):
        rv = []
        for reac in self.sm.cnames:
            rv = Set.Union(rv,self.sm.InvolvedWith(reac,lambda x: x<0).keys())
        return rv


    def Products(self):
        rv = []
        for reac in self.sm.cnames:
            rv = Set.Union(rv,self.sm.InvolvedWith(reac,lambda x: x>0).keys())
        return rv

            
    def Reactants(self):
        return self.sm.rnames[:]


    def CombinedSM(self, Other, UseMets=[]):

        cols = self.sm.cnames + Other.sm.cnames

        rv = StoMat.StoMat(Conv=self.sm.Conv,rnames=UseMets,cnames=cols)
        for srcsm in self.sm,Other.sm:
            rows = Set.Intersect(srcsm.rnames, UseMets)
            for row in rows:
                for col in srcsm.cnames:       
                    rv[row,col]=srcsm[row,col]
        return rv

    def SubSysModes(self, Other):

        sys.stderr.write( "Subsys modes\n")

        sReactants= self.Reactants()
        oReactants =Other.Reactants()
        Interns = Set.Intersect(sReactants, oReactants)
        if len(Interns)==0:
            return []
        AllMets = Set.Union(sReactants, oReactants)
        Intsm = self.CombinedSM(Other, Interns)
        Extsm = self.CombinedSM(Other, AllMets)

        Subtree = ReacTree(Smholder(Intsm,Extsm))
        Subtree.ElModes()

        print "Subsys modes end"
        return Subtree.ems.values() + Subtree.ParModes.values()

        

    def CombinedWith(self, Other):
        if Set.DoesIntersect(self.Balances, Other.Reactants()) or Set.DoesIntersect(self.Reactants(), Other.Balances):
            return []

        

        sReactants = self.Reactants()
        oReactants = Other.Reactants()
        Interns = Set.Intersect(sReactants, oReactants)
        Externs = Set.Difference(sReactants,oReactants)
        Balances = Set.Union(self.Balances, Other.Balances)
        LeafReacs = Set.Union(self.Uses,Other.Uses)
        
        print "Combining", self.Name, Other.Name, "sharing", Interns
        UnDefSto = True

        if Interns != []:
            #if not (self.DefinedSto() and Other.DefinedSto()):
            #    return self.SubSysModes(Other)
            
            Intsm = self.CombinedSM(Other, Interns)  
            print "intsm\n", Intsm,"\n"
            k = Intsm.NullSpace()
            print "k\n", k,"\n"
            if k != None:
                kdim = len(k.cnames)
                UnDefSto = False
                Extsm = self.CombinedSM(Other, Externs)
                if Externs !=[]:
                    stoich = Sci.mul(Extsm,k).tolist()                                                               # stoich of the new mode
                    Extsm.cnames = k.cnames
                    Extsm.rows = stoich
                    print "extsm orig\n", Extsm
                    name = ListToName(self.sm.cnames + Other.sm.cnames)
                    if kdim == 1:
                         Extsm.cnames = [name]
                    else:
                        Extsm.cnames = []
                        for cn in k.NamesNonZeroesInCols():
                            Extsm.cnames.append(ListToName(cn))
                Balances = Set.Union(Balances, Intsm.rnames)
                
        if UnDefSto:
            Extsm = self.CombinedSM(Other,Externs+Interns)
        Extsm.ZapZeroes()
        if UnDefSto: sys.stderr.write("!! ")
        print "final sto", Extsm    
        return [PartialMode(Extsm,Balances,LeafReacs)]


            

               
class ReacTree(Tree):
    def __init__(self, Model=None,*args,**kwargs):
        Tree.__init__(self,*args,**kwargs)
        self.Model = Model
        if Model != None:
            self.FromNJStr(self.Model.sm.ReactionNJTree(),Constructor=ReacTree)
            self.SetModelData()


    def __str__(self):
        if self.IsLeaf():
            return str(self.UserData)
        return Tree.__str__(self)    
          

    def SetModelData(self,bid=0):
        if self.Parent !=None:
            self.Model = self.Parent.Model
            self.bid = bid  # binary id
        else:
            self.bid=1

        self.strbid = str(bid)
            
        chbid=0
        for ch in self.Children:
            ch.SetModelData(self.bid*2+chbid)
            chbid+=1
            


    def HasLeaf(self, name):
        return Tree.HasLeaf(self,name,lambda leaf,name: leaf.UserData==name)



    def Reactions(self):
        """ pre: True
          post: list of reactions in this tree """

        try:
            return self.__reactions
        except:
            if self.IsLeaf():
                self.__reactions = [self.UserData] # UserData set by base class FromNJStr
            else:
                self.__reactions = []
                for l in  self.AllLeaves():
                    self.__reactions.extend(l.Reactions())
        return self.__reactions


    def Metabolites(self):

        try:
            return self.__metabolites
        except:
            self.__metabolites = []
            for reac in self.Reactions():
                mets = self.Model.smexterns.InvolvedWith(reac).keys()
                self.__metabolites= Set.Union(self.__metabolites,mets)
        return self.__metabolites



    def ExtMetabolites(self,Globals=True):
        """pre: True
         post: list of internal metabolites used by self external to the model or used by non-descendent nodes """

        mets = self.Metabolites() # metabolites we use
        if self.IsLeaf():   # if we are a leaf node, all our mets are external
            return mets
        root = self.Root()
        rv = []
        for met in mets:
            subtree =root.Uses([met])    # minimal subtree that uses met
            if not (subtree == self or  subtree.IsSub(self)):
                rv.append(met)                  # met used outside of self
            elif Globals and met in self.Model.sm.Externs:
                rv.append(met)
        return rv



    def Contains(self,Reactions=[]):
        """ pre: self.HasLeaves(Reactions), len(Reactions)>1
           post: minimal subtree containing Reactions"""


        if self[0].HasLeaf(Reactions[0]):
            ch = self[0]
        elif self[1].HasLeaf(Reactions[0]):
            ch = self[1]
        else:
            return self

        for r in Reactions[1:]:
            if not ch.HasLeaf(r):
                return self
        return ch.Contains(Reactions)


    def Uses(self,Metabolites=[]):
        """pre: Metabolite is used by at least one reaction under self
          post: minimal subtree of reactions using Metabolite """

        reacs = []
        for met in Metabolites:
            reacs.extend(self.Model.smexterns.InvolvedWith(met).keys())
        return self.Root().Contains(reacs)



       
    def GetDoneParents(self):
        """ pre: True
          post: list of nodes for whom all children have a stoichiometry matrix """


        if self.DoneModes:
            return []
        if self[0].DoneModes and self[1].DoneModes:
            return [self]

        
        return self[0].GetDoneParents() + self[1].GetDoneParents()



        
    def PromoteCombinations(self):
 
        pms0 = self[0].ParModes
        pms1 = self[1].ParModes
        pms = self.ParModes
        ems = self.Root().ems
        Externs = self.Model.sm.Externs
                                                                                  
        for reac0 in pms0.keys():                                                   # now find and combine usable pairs of reactions
            pm0 = pms0[reac0]
            for reac1 in pms1.keys():
                modes = []
                pm1 = pms1[reac1]
                if not pm0.DefinedSto() or not pm1.DefinedSto():
                    if len(Set.Union(pm0.Uses,pm1.Uses)) < len(self.Model.sm.cnames):
                        modes = pm0.SubSysModes(pm1)
                else:
                    modes = pm1.CombinedWith(pm0)
                        
                for pmnew in modes:
                    if Set.IsSubset(pmnew.Reactants(), Externs):
                        dest = ems
                        sys.stderr.write("Elementary " + pmnew.Name +"\n")
                    else:
                        dest = pms
                    dest[pmnew.Name] = pmnew
                
                
                                            
    def PromoteChildren(self):

        ParModes = self.ParModes
        ExtMets = self.ExtMetabolites()
        print  ExtMets
        
        for ch in self[0],self[1]:                           # promote children that can produce/consume metabolites external to this node
            chpms = ch.ParModes
            for reac in chpms.keys():
                chpm = chpms[reac]
                mets = chpm.Reactants()
 
                if Set.DoesIntersect(mets, ExtMets):                # involved with at least one metabolite external to ourself
                    ParModes[reac] = ch.ParModes[reac]     # record it as one of our partial modes
   


    def CalcModes(self):

        print 
        if self.IsRoot(): print "Root" 
     
        self.ParModes = {}                                               # store partial modes in terms of leaf reactions
        self.ElType = self.Model.sm.Conv
       
        self.PromoteChildren()
        self.PromoteCombinations()
  
        self.DoneModes=True

        
                            

    def ElModes(self):
        """ pre: self.IsRoot()
          post: return pair of matrices sto, mo - stoichiometry matrix and reaction matrix of elmodes """
  
        self.ems = {}
        sm = self.Model.smexterns

        def SetFlags(self):
            if self.IsLeaf():
                self.DoneModes=True
            else:
                self.DoneModes=False
                SetFlags(self[0])
                SetFlags(self[1])

        SetFlags(self)
            
        for l in self.AllLeaves():
            reac = str(l)
            l.ParModes={reac:
                        PartialMode(
                            sm.SubSys([reac]),
                            Irrev = sm.IsIrrev(reac),
                            Uses= [reac]
                        )
                    }
                           
        while not self.DoneModes:
            for lp in self.GetDoneParents():
               lp.CalcModes()