/home/mark/model/software/ScrumPy/ScrumPy/Kinetic/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

"""

"""
M.G.Poolman 2001
Kinetic.Model - a user level interface to kinetic models
"""


import os,sys,  random, exceptions

from BasicInfo import Location, HeaderPath
from Data import DataSets
from Util import File,Seq,Set
#from ScrumPy import gui

import ScrumPy
import LowLevelKinMod, CGen

gui=None  # must be set by importing module
#TODO - remove gui from here, report errors etc back to calling functions
"""
Things relevant to user (treat as ReadOnly !!)
"""

# these 4 define valid values for the "code" argument to GetVec and GetNameVec
Param = LowLevelKinMod.Param
Conc = LowLevelKinMod.Conc
Vel = LowLevelKinMod.Vel
dMet = LowLevelKinMod.dMet



# a string representing time, avoids use of literals in code
TimeString = "Time"


LowLevel = LowLevelKinMod.LowLevelKinMod

ErrMsgTitle = "ScrumPy  error"


class Monitor(DataSets.DataSet):

    def __init__(self, model, Params=[], UpdateIfNotOK="Ignore"):

        if Params==[]:
            Params = model.GetVecNames(Param)

        self.model=model
        self.VarNames = ["Time"] + model.sm.cnames + model.sm.rnames
        self.Params = Params[:]
        self.UpdateIfNotOK=UpdateIfNotOK

        ItemNames = self.VarNames+Params

        DataSets.DataSet.__init__(self, ItemNames=ItemNames, SrcDic=model)


    def Update(self):

        if self.model.IsOK() or self.UpdateIfNotOK=="Yes":
            DataSets.DataSet.Update(self)
        elif self.UpdateIfNotOK=="Zero":
            vars = [0]*len(self.VarNames)
            pars = self.model.GetVals(self.Params)
            self.NewRow(vars+pars)
        elif self.UpdateIfNotOK!="Ignore":
            raise exceptions.ValueError, "Unknown value "+self.UpdateIfNotOK+" for UpdateIfNotOK"

        #else do nothing











def DefaultSimMon(Model, DataSet=None):
    pass

#TODO this needs moving into Util for both LP and Kinetic (and anything else thatmightr want it
from LP import Array
class DoubleArray(Array.BaseArray):
    _del = LowLevelKinMod.delete_doubleArray
    _get = LowLevelKinMod.doubleArray_getitem
    _set = LowLevelKinMod.doubleArray_setitem
    _new = LowLevelKinMod.new_doubleArray
    _cast = float
    _offset=0

    def _NoDel(self):
        pass

    def __init__(self,src=None, length=None, NoDel=False):

        if src == length == None:
            raise exceptions.ValueError, "src  and length can't both be None"

        if src == None:
            Array.BaseArray.__init__(self, length)

        elif type(src) == list:
             Array.BaseArray.__init__(self, src)

        else:
            self._array =  src
            self.size = length

        if NoDel:
            self._del = self._NoDel




class Model (LowLevel):
    """ instances of this class are the ScrumPy modellers main tool.
        House keeping methods are at the begining, user methods at the end
    """
    def __init__(self):


        LowLevel.__init__(self)
        self.Sim = LowLevel.Simulate   # a convenient handle  on the low level simulate when we overload it


        self.SimSSRetry = 4
        self.SimSSPoints = 128
        self.SimSSTincr = 0.125
        self.SimSSGrowth = 1.5
        self.NSteps = 1 # no of sim steps on last simulation

        self.SSTol = 1e-5

        self.BuildKinetic()



    def BuildKinetic(self): #

        ModelDesc = self.ModelDesc
        SrcFName = self.FileName+".c"
        self.MakeC(open(SrcFName,"w"))

        SoFName = self.FileName+".so"
        BuildSo  ="g++ -fPIC -shared -I " + HeaderPath + " -o " + SoFName + " " + SrcFName
        syserr = os.system(BuildSo)
        if syserr:
            gui.ErrorMsg("Failed to compile "+ self.ModelName +"\nsee Python terminal and stderr for clues")
            print "!!! ", BuildSo, " generated error ", str(syserr), "!!!\n"

        else:
            if not File.IsAbsPath(SoFName):                 # because dlib/ld.so need absolute path names
                SoFName = os.curdir + os.sep + SoFName
            LowLevel.Load(self, self.ModelName,SoFName)

            # we've generated, compiled and loaded the model, now to initialise it

            self.ZerVec(Conc)       # it's OK to init concentrations to zero by default
            self.NaNVec(Param)    # but not parameters, this ensures the user *will* find out

            ConsDic = ModelDesc.ConsDic
            nCons = len(ConsDic)
            if nCons >0:
                ParamNames = self.GetVecNames(Param)[:-nCons] # we will generate conservation params later
            else:
                ParamNames = self.GetVecNames(Param)

            IniDic = ModelDesc.InitDic
            IniNames = IniDic.keys() # initialise values as supplied in the Ini dictionary
            for Name in IniNames:
                self[Name] = IniDic[Name]

            self.InitCSums()  # sort out values for conserved totals

            Seq.RemoveFrom(ParamNames, IniNames)  # see which params are not initialised
            if len(ParamNames) > 0:
                gui.WarnMsg("Parameters not initialised" + str(ParamNames))

            MetNames = self.GetVecNames(Conc)
            Seq.RemoveFrom(MetNames, IniNames)  # see which mets are not initialised
            if len(MetNames) > 0:
                gui.WarnMsg("Metabolites not initialised" + str(MetNames))


            self.Eval()     # initialise velocities and dependent mets
            self.ResetInteg() # ensure integrator has correct start values
            self.SetRescue()


    def ConsDic(self):

        cm = self.ConsMoieties(self)
        for c in cm.cnames:
            cm.DelRow(c)
        return cm.ToColDict()






    def MakeC(self, out):
            md = self.ModelDesc
            sm = self.sm
            md.ConsDic = self.ConsDic()
            cons = md.ConsDic.keys()
            md.conpard = conpard = {}
            conparl = []
            for con in cons:
                self.sm.MakeLastRow(con)
                self.smexterns.MakeLastRow(con)
                conp = "CSUM_"+con
                conpard[con] = conp
                conparl.append(conp)


            md.Params = self.ParamNames = Set.Complement(md.IdList, sm.rnames+sm.cnames) + conparl
            #CGen.InitMD(md)
            out.write(CGen.Preamble)
            out.write("static double "+self.ModelName+"Time ; \n")
            out.write(CGen.VecDecls(md))
            out.write(CGen.Enums(md))
            out.write(CGen.NameArrays(md))
            out.write(CGen.DiffEqns(md))
            out.write(CGen.RateEqns(md))
            out.write(CGen.ConsEqns(md))
            out.write(CGen.Struct(md))



    def InitCSums(self):

        ConsDic = self.ModelDesc.ConsDic
        conpard = self.ModelDesc.conpard

        for k in ConsDic.keys():
            sum = self[k]
            tuplist = ConsDic[k]
            for tup in tuplist:
                sum += float(self[tup[0]]) * float(tup[1])

            self[conpard[k]] = sum




    def Eval(self, *args): # just ignore args
        LowLevel.Eval(self)


    def __del__(self):
        self.Destroy()

    def Destroy(self):
        LowLevel.Destroy(self)





    def __getitem__(self, name):
        return self.GetVal(name)

    def __setitem__(self, name, val):
        self.SetVal(name, val)

    def ReportBadPoints(self, Source, Points):
       """ do something gui at a later date """


       print "BadPoints from ", Source, ":\n", Points



       ### User methods start here  ################


    def GetVals(self, names):
        """ pre: names is a list of strings
           post: returns a list of the values of names """

        rv = []
        for n in names:
            rv.append(self[n])

        return rv

    def SetVals(self, names, vals):
        """ pre: names is a list of strings, vals is a list of float,
                 len(names) == len (vals)
           post: GetVals(names) == vals"""

        for i in range(len(names)):
            self[names[i]] = vals[i]


    def SetTime(self, time):
        """pre: time <=0.0
          post: self["Time"] == time
          NB: this ensures internal consistancy in the C library when time is changed.
              m["Time"] = t is best avoided."""

        self["Time"] = time
        self.ResetInteg()






    def GetVec(self,code):
        """ return a list of values of the vector associated with code """
        #return self.GetVals(self.GetVecNames(code))
        return DoubleArray(LowLevel.GetVec(self, code), LowLevel.SizeVec(self, code),NoDel=True)

    def SetVec(self, code, Vec):
        """ pre: len(Vec) == self.SizeVec(code), type(Vec[0..len(Vec)-1] == float
           post: Vec == self.GetVec(code) """

        return self.SetVals(self.GetVecNames(code),Vec)

    def GetVecNames(self, code):
        """ return a list of names of the values of the vector associated with code """
        vec = {
            Conc: self.sm.rnames,
            Vel: self.sm.cnames,
            Param:self.ParamNames}[code]
        return vec[:]


    def GetAllNames(self):
        """ return a list (of strings) of all identifiers in the model """
        return(
            self.GetVecNames(Vel) +
            self.GetVecNames(Conc) +
            self.GetVecNames(Param) +
            [TimeString]
        )

    def keys(self):
        return self.GetAllNames()
    """
    because models anyway have dictionary-like characteristics we supply
    a keys function for completeness. This is useful for the monitor
    classes which can then be updated from any dicionary-like instance
    """


    def Randomise(self, lorel=0.5, hirel=0.5, *args, **kwargs):

        for name,  val in self.ModelDesc.InitDic.items():
            self[name] = value * ( 1+random.uniform(-lorel, hirel))

        self.InitCSums()
        self.Eval()
        self.ResetInteg()


    def PingConc(self, pert=1e-6):
        """ perturb all concentrations by multiplying by a random number [1-pert..1+pert]"""

        concs = self.GetVec(Conc)
        for c in range(len(concs)):
            concs[c] *= random.uniform(1-pert,1+pert)




    def FindSS(self,  Retry=100):
        """ Use a combination of simulation and Newton to (attempt to) obtain a steady-state"""

        curtime =self["Time"]

        self.Eval()

        while self.GetSSErr() > self.SSTol *2 and Retry >0:
            self.Sim(self,1)
            Retry-=1
        self.NewtSS()

        self.SetTime(curtime) # restore time to starting pos.



    def PureSimSS(self, MaxSteps=100):
        """ Attempt to get to steady-state using only time course simulation """

        print "!! Deprecated PureSimSS() !!"
        # not a lot of point given current FindSS()
        self.SimToSS(MaxSteps, self.SSTol)


    def ScaledSensits(self, pname, Vars, Scaled=True,  SS=True, Pert=1e-3):

        if SS:
            def Eval():
                self.PingConc()
                self.FindSS()
                #if not self.IsOK():
                #    print "!!"
        else:
            Eval = self.CalcVels

        pval = self[pname]
        delta_p = pval*Pert
        TwodP = 2 * delta_p

        if Scaled:
            Eval()
            midvars = self.GetVals(Vars)   # mid of 3 point differential

        self[pname] += delta_p
        Eval()
        hivars = self.GetVals(Vars)         # high point of 3 point diff

        self[pname] -= TwodP
        Eval()
        lovars = self.GetVals(Vars)         # low point of 3 point diff

        self[pname] = pval

        dydx = map(lambda yh, yl:(yh-yl)/(TwodP),  hivars, lovars)
        if Scaled:
            dydx = map(lambda d, y:d*pval/y, dydx, midvars)

        return dydx





    def NewDataSet(self,Params=[],UpdateIfNotOK="Ignore"):
        """ UpdateIfNotOK == "Ignore" => only update if self.IsOK()
                             "Zero"   => update vars with [0,0..] but update Params if not self.IsOK()
                             "Yes"    => Always update"""
        #TODO: specify whether static or dynamic here also

        ds = Monitor(self,[],UpdateIfNotOK)
        ds.MakeFirstCol("Time")
        return ds



    def SetFromDataSet(self, ds, n):
        """ pre: n is an integer 0 <= n < len(ds)
           post: values of ds.keys() in self updated from the nth record of ds """

        for k in ds.cnames:
            self[k] = ds[k][n]



    def GetStateSaver(self):
        """ returns a StateSaver instance """


        class StateSaver(DataSets.DataSet):
            raise "Kinetic.Model.Model.StateSaver", "needs re-writing - DataSet API changed"

            def __init__(self, model):
                self.names = [TimeString]
                for comp in [Conc, Vel, Param]:
                    self.names = self.names+ model.GetVecNames(comp)

                self.model = model
                DataSets.DataSet.__init__(self, FromFile=0, NameList=self.names)

            #def Update(self):
            #    self.model.UpdateDataSet(self)

            def SaveFile(self, fname=None):
                DataSets.DataSet.SaveFile(self, fname=fname, InOrder=self.names)

        return  StateSaver(self)


    def GetState(self):
        rv = {}
        rv["Time"] = self["Time"]
        for vec in [ScrumPy.Conc, ScrumPy.Vel, ScrumPy.Param]:
            rv[vec] = self.GetVec(vec)
        return rv



    def SetFromState(self, state):
        """ pre: state = m.GetState(), self = m """

        self["Time"] = state["Time"]
        for vec in [ScrumPy.Conc, ScrumPy.Vel, ScrumPy.Param]:
            self.SetVec(vec, state[vec])