"""
PySCeS - Python Simulator for Cellular Systems (http://pysces.sourceforge.net)
Copyright (C) 2004-2024 B.G. Olivier, J.M. Rohwer, J.-H.S Hofmeyr all rights reserved,
Brett G. Olivier (bgoli@users.sourceforge.net)
Triple-J Group for Molecular Cell Physiology
Stellenbosch University, South Africa.
Permission to use, modify, and distribute this software is given under the
terms of the PySceS (BSD style) license. See LICENSE.txt that came with
this distribution for specifics.
NO WARRANTY IS EXPRESSED OR IMPLIED. USE AT YOUR OWN RISK.
Brett G. Olivier
"""
from __future__ import division, print_function
from __future__ import absolute_import
from __future__ import unicode_literals
from pysces.version import __version__
__doc__ = '''
PyscesLink
----------
Interfaces to external software and API's, has replaced the PySCeS contrib classes.
'''
try:
input = raw_input # Py2 compatibility
except NameError:
pass
# for METATOOLlink
import os, re, io
# for SBWWebLink
try:
from urllib.request import (
ProxyHandler,
build_opener,
HTTPHandler,
install_opener,
urlopen,
) # Py 3
except ImportError:
from urllib2 import (
ProxyHandler,
build_opener,
HTTPHandler,
install_opener,
urlopen,
) # Py 2
try:
from urllib.parse import urlencode # Py 3
except ImportError:
from urllib import urlencode # Py 2
[docs]class SBWlink(object):
"""Generic access for local SBW services using SBWPython """
sbw = None
psbw = None
sbwModuleProxy = None
moduleDict = None
modules = None
def __init__(self):
try:
import SBW as SBW
import SBW.psbw as psbw
## reload(SBW)
## reload(psbw)
self.sbw = SBW
self.psbw = SBW.psbw
self.sbwModuleProxy = SBW.sbwModuleProxy
self.moduleDict = SBW.sbwModuleProxy.moduleDict
self.modules = []
for m in self.moduleDict:
if self.moduleDict[m].pythonName not in ['python']:
self.SBW_exposeAll(self.moduleDict[m])
self.modules.append(self.moduleDict[m].pythonName)
setattr(self, self.moduleDict[m].pythonName, self.moduleDict[m])
print('\nSBWlink established.')
except Exception as ex:
print(ex)
print('\nSBWlink not established.')
[docs] def SBW_exposeAll(self, module):
for s in module.services:
s = getattr(module, s)
for m in s.methods:
getattr(s, m)
[docs] def SBW_getActiveModules(self):
idlist = []
namelst = []
for id in self.psbw.getModuleIdList():
idlist.append(id)
namelst.append(self.psbw.getModuleName(id))
for id in list(self.moduleDict.keys()):
if id not in idlist:
self.moduleDict.pop(id)
for name in range(len(self.modules) - 1, -1, -1):
if self.modules[name] not in namelst:
delattr(self, self.modules[name])
self.modules.pop(name)
for name in namelst:
if name not in self.modules:
self.SBW_loadModule(name)
return namelst
[docs] def SBW_loadModule(self, module_name):
ans = 'Y'
if module_name[-3:] == 'GUI':
ans = input(
"Warning! This may hang the console\nPress 'Y' to continue: "
)
if ans == 'Y':
module_id = self.psbw.SBWGetModuleInstance(module_name)
assert module_id != None, '\nUnknow module, %s' % module_name
module = self.sbwModuleProxy.ModuleProxy(module_id)
self.SBW_exposeAll(module)
if module_id not in self.moduleDict:
print(
'<PySCeS_SBW> Adding '
+ module.pythonName
+ ' to ModuleProxy (id='
+ str(module_id)
+ ')'
)
self.moduleDict.update({module_id: module})
if module.pythonName not in self.modules:
print('<PySCeS_SBW> Adding ' + module.pythonName + ' to SBWlink')
self.modules.append(module.pythonName)
setattr(self, module.pythonName, module)
else:
print('\nModule %s not loaded' % module_name)
[docs]class SBWLayoutWebLink(object):
"""Enables access to DrawNetwork and SBMLLayout web services at www.sys-bio.org"""
sbwhost = '128.208.17.26'
sbml = None
sbmllayout = None
svg = None
DEBUGMODE = False
DEBUGLEVEL = 1
DRAWNETWORKLOADED = False
LAYOUTMODULELOADED = False
[docs] def setProxy(self, **kwargs):
"""Set as many proxy settings as you need. You may supply a user name without
a password in which case you will be prompted to enter one (once) when required
(NO guarantees, implied or otherwise, on password security AT ALL). Arguments can be:
user = 'daUser',
pwd = 'daPassword',
host = 'proxy.paranoid.net',
port = 3128
"""
proxy_info = {}
for k in list(kwargs.keys()):
proxy_info.update({k: kwargs[k]})
if 'user' in proxy_info and 'pwd' not in proxy_info:
proxy_info.update({'pwd': getpass.getpass()})
proxy_support = ProxyHandler(
{"http": "http://%(user)s:%(pwd)s@%(host)s:%(port)d" % proxy_info}
)
opener = build_opener(proxy_support, HTTPHandler)
install_opener(opener)
del proxy_info, proxy_support
[docs] def loadSBMLFileFromDisk(self, File, Dir=None):
if Dir != None:
path = os.path.join(Dir, File)
else:
path = File
if os.path.exists(path):
self.sbmllayout = None
self.svg = None
self.DRAWNETWORKLOADED = False
self.LAYOUTMODULELOADED = False
sbmlF = open(path, 'r')
self.sbml = sbmlF.read()
sbmlF.close()
return True
else:
print("%s is an invalid path" % path)
return False
[docs] def loadSBMLFromString(self, str):
self.sbmllayout = None
self.svg = None
self.DRAWNETWORKLOADED = False
self.LAYOUTMODULELOADED = False
self.sbml = str
return True
[docs] def urlGET(self, host, urlpath):
url = 'http://%s%s' % (host, urlpath)
con = urlopen(url)
resp = con.read()
if self.DEBUGMODE:
print(con.headers)
if self.DEBUGMODE and self.DEBUGLEVEL == 2:
print(resp)
con.close()
return resp
[docs] def urlPOST(self, host, urlpath, data):
assert type(data) == dict, '\nData must be a dictionary'
url = 'http://%s%s' % (host, urlpath)
con = urlopen(url, urlencode(data))
resp = con.read()
if self.DEBUGMODE:
print(con.headers)
if self.DEBUGMODE and self.DEBUGLEVEL == 2:
print(resp)
con.close()
return resp
[docs] def getVersion(self):
print('Inspector.getVersion()')
ver = self.urlGET(self.sbwhost, '/generate/Inspector.asmx/getVersion')
ver = ver.replace('<?xml version="1.0" encoding="utf-8"?>', '')
ver = ver.replace('<string xmlns="http://www.sys-bio.org/">', '')
ver = ver.replace('</string>', '')
return ver
[docs] def drawNetworkLoadSBML(self):
print('DrawNetwork.loadSBML()')
assert self.sbml != None, '\nNo SBML file loaded'
data = {'var0': self.sbml}
self.DRAWNETWORKLOADED = True
return self.urlPOST(self.sbwhost, '/generate/DrawNetwork.asmx/loadSBML', data)
[docs] def drawNetworkGetSBMLwithLayout(self):
print('DrawNetwork.getSBML()')
assert self.DRAWNETWORKLOADED, '\nSBML not loaded into DrawNetwork module'
sbml = self.urlGET(self.sbwhost, '/generate/DrawNetwork.asmx/getSBML')
sbml = sbml.replace('>', '>')
sbml = sbml.replace('<', '<')
sbml = sbml.replace(
'''<string xmlns="http://www.sys-bio.org/"><?xml version="1.0" encoding="utf-8"?>''',
'',
)
sbml = sbml.replace('</string>', '')
self.sbmllayout = sbml
return True
[docs] def layoutModuleLoadSBML(self):
print('SBMLLayoutModule.loadSBML()')
assert self.sbmllayout != None, '\nNo SBML Layout loaded'
data = {'var0': self.sbmllayout}
self.LAYOUTMODULELOADED = True
return self.urlPOST(
self.sbwhost, '/generate/SBMLLayoutModule.asmx/loadSBML', data
)
[docs] def layoutModuleGetSVG(self):
assert self.LAYOUTMODULELOADED, '\nSBML not loaded into SBMLLayout module'
svg = self.urlGET(self.sbwhost, '/generate/SBMLLayoutModule.asmx/getSVG')
svg = svg.replace('>', '>')
svg = svg.replace('<', '<')
svg = svg.replace('''<string xmlns="http://www.sys-bio.org/">''', '')
svg = svg.replace('''<?xml version="1.0" encoding="utf-8"?>''', '')
svg = svg.replace('</string>', '')
self.svg = svg
return True
[docs] def getSBML(self):
return self.sbml
[docs] def getSBMLlayout(self):
return self.sbmllayout
[docs] def getSVG(self):
return self.svg
'''
_HAVE_STOMPY = False
_STOMPY_LOAD_ERROR = ''
try:
## import stompy
import stochpy as stompy
_HAVE_STOMPY = True
except Exception, ex:
_STOMPY_LOAD_ERROR = '%s' % ex
_HAVE_STOMPY = False
class StomPyInterface(object):
"""
StomPy interface to PySCeS this may move to pysces.link in the future
"""
SSA = None
SSA_REACTIONS = None
SSA_SPECIES = None
stompy = None
_MOD2PSC = None
TMP_FNAME = None
TMP_PATH = None
MODEL_PATH = None
OUTPUT_PATH = None
STP_IS_TIME_SIM = False
STP_METHOD = 'Direct'
STP_TIMEEND = 1
STP_TRAJ = 1
STP_INTERACTIVE = True
STP_TRACK_PROPENSITIES = True
STP_WAITING_TIMES = True
STP_STEPS = 10
STP_INITIAL_SPECIES = True
STP_KEEP_PSC_FILES = False
def __init__(self, model_path, output_path):
"""
An interface class to the StomPy stochastic simulator
- *model_path* the default PySCeS model directory
- *output_path* the default PySCeS output directory
"""
self.stompy = stompy
self.OUTPUT_PATH = output_path
self.MODEL_PATH = model_path
self.TMP_PATH = os.path.join(model_path, 'orca')
self._MOD2PSC = interface.writeMod2PSC
def setProperty(self, **kwargs):
"""
Sets a StomPy simulation parameter
- *method* [default='Direct'] select simulation algorithm
- *trajectories* [default=1]
- *interactive* [default=True]
- *track_propensities* [default=True]
- *steps* [default=10]
"""
## print kwargs
if kwargs.has_key('method'):
self.STP_METHOD = kwargs['method']
## print '%s = %s' % ('method', kwargs['method'])
if kwargs.has_key('trajectories'):
self.STP_TRAJ = kwargs['trajectories']
self.STP_TRAJ = 1 # TODO I need to look into this
## print 'Currently only single trajectories are supported via the PySCeS interface'
## print '%s = %s' % ('trajectories', self.STP_TRAJ)
if kwargs.has_key('interactive'):
self.STP_INTERACTIVE = kwargs['interactive']
## print '%s = %s' % ('interactive', kwargs['interactive'])
if kwargs.has_key('track_propensities'):
self.STP_TRACK_PROPENSITIES = kwargs['track_propensities']
## print '%s = %s' % ('track_propensities', kwargs['track_propensities'])
if kwargs.has_key('steps'):
self.STP_STEPS = kwargs['steps']
## print '%s = %s' % ('steps', kwargs['steps'])
if kwargs.has_key('species_initial'):
self.STP_INITIAL_SPECIES = kwargs['initial_species']
## print '%s = %s' % ('initial_species', kwargs['initial_species'])
if kwargs.has_key('keep_psc_files'):
self.STP_KEEP_PSC_FILES = kwargs['keep_psc_files']
## print '%s = %s' % ('keep_psc_files', kwargs['keep_psc_files'])
def initModelFromMod(self, pscmod, iValues=False):
"""
Initialise a StomPy SSA instance from a PySCeS model.
- *pscmod* an initialised PySCeS model
- *iValues* [default=False] use initial values (not current)
"""
self.TMP_FNAME = str(time.time()).replace('.','')+'.psc'
if self.STP_INITIAL_SPECIES:
for s in pscmod.species:
setattr(pscmod, s, pscmod.__sDict__[s]['initial'])
self._MOD2PSC(pscmod, self.TMP_FNAME, self.TMP_PATH, iValues=iValues)
self.SSA = self.stompy.SSA(Method=self.STP_METHOD, File=self.TMP_FNAME, dir=self.TMP_PATH, IsInteractive=self.STP_INTERACTIVE)
self.SSA.Trajectories(self.STP_TRAJ)
self.SSA_REACTIONS = self.SSA.SSA.rate_names
self.SSA_SPECIES = self.SSA.SSA.species
if self.STP_TRACK_PROPENSITIES:
self.SSA.TrackPropensities()
try:
print os.path.join(self.TMP_PATH, self.TMP_FNAME)
if not self.STP_KEEP_PSC_FILES and self.TMP_PATH != None and self.TMP_FNAME != None:
os.remove(os.path.join(self.TMP_PATH, self.TMP_FNAME))
except:
print 'Could not delete intermediatery StomPy PSC file: %s' % os.path.join(self.TMP_PATH, self.TMP_FNAME)
self.TMP_FNAME = None
print 'StomPy model ... initialised.'
def runTimeSimulation(self, pscmod, endtime=None, method='Direct', iValues=False):
"""
Run a stochastic simulation
- *pscmod* and instanitiated PySCeS model
- *endtime* [default=1] the end time **Note: this could take a long time i.e. generate ahuge amount of steps**
- *method* [default='Direct'] select the simulation method, one of:
- *Direct*
- *FirstReactionMethod*
- *NextReactionMethod*
- *TauLeaping*
- *iValues* [default=False] use initial values (not current)
"""
if method not in ['Direct','FirstReactionMethod','NextReactionMethod','TauLeaping']:
print 'Method: %s does not exist using - Direct' % method
self.STP_METHOD = 'Direct'
else:
self.STP_METHOD = method
if endtime != None:
self.STP_TIMEEND = endtime
self.initModelFromMod(pscmod, iValues=iValues)
## self.SSA.Timesteps(self.STP_STEPS)
self.SSA.Endtime(self.STP_TIMEEND)
self.SSA.Run()
self.STP_IS_TIME_SIM = True
## self.SSA.PlotTimeSim()
print 'StomPy time simulation ... done.'
# TODO STOCHPY
## if self.SSA.SSA.output[0][-1] == '':
## self.SSA.SSA.output[0][-1] = 0.0
## sim_dat = numpy.array(self.SSA.SSA.output, 'd')
## pscmod.data_stochsim = IntegrationStochasticDataObj()
## pscmod.data_stochsim.setTime(sim_dat[:,0])
## pscmod.data_stochsim.setSpecies(sim_dat[:,1:-1], self.SSA_SPECIES)
pscmod.data_stochsim = self.SSA.data_stochsim
if self.STP_WAITING_TIMES:
wtimes, wt_lbls = self.getWaitingtimesData(reactions=None,lbls=True)
pscmod.data_stochsim.setWaitingtimes(wtimes, wt_lbls)
if self.STP_TRACK_PROPENSITIES:
pscmod.data_stochsim.setPropensities(self.SSA.SSA.propensities_output)
pscmod.data_stochsim.TYPE_INFO = 'Stochastic'
def runStepSimulation(self, pscmod, steps=None, method='Direct', iValues=False):
"""
Run a stochastic simulation
- *pscmod* and instanitiated PySCeS model
- *steps* [default=10] the number of steps to simulate
- *method* [default='Direct'] select the simulation method, one of:
- *Direct*
- *FirstReactionMethod*
- *NextReactionMethod*
- *TauLeaping*
- *iValues* [default=False] use initial values (not current)
"""
if method not in ['Direct','FirstReactionMethod','NextReactionMethod','TauLeaping']:
print 'Method: %s does not exist using - Direct' % method
self.STP_METHOD = 'Direct'
else:
self.STP_METHOD = method
if steps != None:
self.STP_STEPS = steps
self.initModelFromMod(pscmod, iValues=iValues)
self.SSA.Timesteps(self.STP_STEPS)
## self.SSA.Endtime(self.STP_TIMEEND)
self.SSA.Run()
self.STP_IS_TIME_SIM = False
print 'StomPy step simulation ... done.'
## print self.SSA.SSA.output[0]
## print self.SSA.SSA.output[1]
## print self.SSA.SSA.output[-1]
## header_line = self.SSA.SSA.output.pop(0)
## if self.SSA.SSA.output[0][-1] == '':
## self.SSA.SSA.output[0][-1] = 0.0
## sim_dat = numpy.array(self.SSA.SSA.output, 'd')
## pscmod.data_stochsim = IntegrationStochasticDataObj()
## pscmod.data_stochsim.setTime(sim_dat[:,0])
## pscmod.data_stochsim.setSpecies(sim_dat[:,1:-1], self.SSA_SPECIES)
pscmod.data_stochsim = self.SSA.data_stochsim
if self.STP_WAITING_TIMES:
wtimes, wt_lbls = self.getWaitingtimesData(reactions=None,lbls=True)
pscmod.data_stochsim.setWaitingtimes(wtimes, wt_lbls)
if self.STP_TRACK_PROPENSITIES:
pscmod.data_stochsim.setPropensities(self.SSA.SSA.propensities_output)
pscmod.data_stochsim.TYPE_INFO = 'Stochastic'
def getWaitingtimesData(self,reactions=None,lbls=False):
"""
Plots the waiting times for each reaction in the model. Makes use of ObtainWaitingtimes to derive the waiting times out of the SSA output.
Input:
- *reactions* [default=0] a list of reactions to plot defualts to all reactions
- *traj* [default=0] trajectory to plot (defaults to first one)
- *lbls* [default=False] if True return (data_array, column_labels) otherwise just data_array
This method is derived from StomPy 0.9 (http://stompy.sf.net) Analysis.py
"""
if self.SSA.IsTauLeaping:
print 'INFO: waiting times not available when method is Tau Leaping'
if not lbls:
return None
else:
return None, None
self.SSA.GetWaitingtimes()
if reactions == None:
reactions = self.SSA_REACTIONS
vect = []
vect_lbls = []
for r in reactions:
if r in self.SSA_REACTIONS:
vect.append(self.SSA_REACTIONS.index(r)+1)
vect_lbls.append('wt'+str(r))
else:
print "INFO: '%s' is not a valid reaction name" % r
OUTPUT = []
## for t in range(len(self.SSA.data_stochsim.waiting_times)):
T_OUTPUT = []
for i in vect:
if self.SSA.data_stochsim.waiting_times.has_key(i):
waiting_time = self.SSA.data_stochsim.waiting_times[i]
if len(waiting_time) > 1: # At least 2 waiting times are necessary per reaction
T_OUTPUT.append(self.stompy.modules.Analysis.LogBin(waiting_time, 1.5)) # Create logarithmic bins
else:
T_OUTPUT.append(None)
else:
T_OUTPUT.append(None)
OUTPUT.append(T_OUTPUT)
if not lbls:
return OUTPUT
else:
return OUTPUT, vect_lbls
class IntegrationStochasticDataObj(object):
"""
This class is specifically designed to store the
results of a stochastic time simulation
It has methods for setting the Time, Labels, Species and Propensity data and
getting Time, Species and Rate (including time) arrays. However, of more use:
- getOutput(*args) feed this method species/rate labels and it will return
an array of [time, sp1, r1, ....]
- getDataAtTime(time) the data generated at time point "time".
- getDataInTimeInterval(time, bounds=None) more intelligent version of the above
returns an array of all data points where: time-bounds <= time <= time+bounds
"""
time = None
waiting_times = None
species = None
propensities = None
xdata = None
time_label = 'Time'
waiting_time_labels = None
species_labels = None
propensities_labels = None
xdata_labels = None
HAS_SPECIES = False
HAS_WAITING_TIMES = False
HAS_PROPENSITIES = False
HAS_TIME = False
HAS_XDATA = False
IS_VALID = True
TYPE_INFO = 'Stochastic'
def setLabels(self, species):
"""
Set the species
- *species* a list of species labels
"""
self.species_labels = species
def setTime(self, time, lbl=None):
"""
Set the time vector
- *time* a 1d array of time points
- *lbl* [default=None] is "Time" set as required
"""
self.time = time.reshape(len(time), 1)
self.HAS_TIME = True
if lbl != None:
self.time_label = lbl
def setSpecies(self, species, lbls=None):
"""
Set the species array
- *species* an array of species vs time data
- *lbls* [default=None] a list of species labels
"""
self.species = species
self.HAS_SPECIES = True
if lbls != None:
self.species_labels = lbls
def setWaitingtimes(self, waiting_times, lbls=None):
"""
Set the `waiting_times` this data structure is not an array but a nested list of: waiting time log bins per reaction per trajectory::
waiting_times = [traj_1, ..., traj_n]
traj_1 = [wt_J1, ..., wt_Jn] # in order of SSA_REACTIONS
wt_J1 = (xval, yval, nbin)
xval =[x_1, ..., x_n]
yval =[y_1, ..., y_n]
nbin = n
- *waiting_times* a list of waiting times
- *lbls* [default=None] a list of matching reaction names
"""
self.waiting_times = waiting_times
self.HAS_WAITING_TIMES = True
if lbls != None:
self.waiting_time_labels = lbls
def setPropensities(self, propensities, lbls=None):
"""
Sets an array of propensities.
- *propensities* a list of propensities
- *lbls* [default=None] a list of matching reaction names
"""
if lbls == None:
LB = copy.copy(propensities[0])
lbls = LB[1:]
lbls = ['p'+str(r) for r in lbls]
P_ARR = numpy.zeros((len(propensities), len(propensities[0])-1), 'd')
P_ARR[-1,:] = numpy.NaN
for r in range(1, P_ARR.shape[0]):
P_ARR[r, :] = propensities[r][1:]
self.propensities = P_ARR
self.HAS_PROPENSITIES = True
if lbls != None:
self.propensities_labels = lbls
## print self.propensities_labels
## print self.propensities
def setXData(self, xdata, lbls=None):
"""
Sets an array of extra simulation data
- *xdata* an array of xdata vs time
- *lbls* [default=None] a list of xdata labels
"""
self.xdata = xdata
self.HAS_XDATA = True
if lbls != None:
self.xdata_labels = lbls
def getTime(self, lbls=False):
"""
Return the time vector
- *lbls* [default=False] return only the time array or optionally both the time array and time label
"""
output = None
if self.HAS_TIME:
output = self.time.reshape(len(self.time),)
if not lbls:
return output
else:
return output, [self.time_label]
def getSpecies(self, lbls=False):
"""
Return an array fo time+species
- *lbls* [default=False] return only the time+species array or optionally both the data array and a list of column label
"""
output = None
if self.HAS_SPECIES:
output = numpy.hstack((self.time, self.species))
labels = [self.time_label]+self.species_labels
else:
output = self.time
labels = [self.time_label]
if not lbls:
return output
else:
return output, labels
def getWaitingtimes(self, lbls=False, traj=[]):
"""
Return waiting times, time+waiting_time array
- *lbls* [default=False] return only the time+waiting_time array or optionally both the data array and a list of column label
- *traj* [default=[0]] return the firs or trajectories defined in this list
"""
output = None
labels = None
if self.HAS_WAITING_TIMES:
output = []
if len(traj) == 0:
traj = range(len(self.waiting_times))
## if len(traj) == 1:
## output = self.waiting_times[0]
## else:
for t in traj:
output.append(self.waiting_times[t])
labels = self.waiting_time_labels
else:
output = []
labels = []
if not lbls:
return output
else:
return output, labels
def getPropensities(self, lbls=False):
"""
Return time+propensity array
- *lbls* [default=False] return only the time+propensity array or optionally both the data array and a list of column label
"""
#assert self.propensities != None, "\nNo propensities"
output = None
if self.HAS_PROPENSITIES:
print self.time.shape
print self.propensities.shape
output = numpy.hstack((self.time, self.propensities))
labels = [self.time_label]+self.propensities_labels
else:
output = self.time
labels = [self.time_label]
if not lbls:
return output
else:
return output, labels
def getXData(self, lbls=False):
"""
Return time+xdata array
- *lbls* [default=False] return only the time+xdata array or optionally both the data array and a list of column label
"""
output = None
if self.HAS_XDATA:
output = numpy.hstack((self.time, self.xdata))
labels = [self.time_label]+self.xdata_labels
else:
output = self.time
labels = [self.time_label]
if not lbls:
return output
else:
return output, labels
def getDataAtTime(self, time):
"""
Return all data generated at "time"
- *time* the required exact time point
"""
#TODO add rate rule data
t = None
sp = None
ra = None
ru = None
xd = None
temp_t = self.time.reshape(len(self.time),)
for tt in range(len(temp_t)):
if temp_t[tt] == time:
t = tt
if self.HAS_SPECIES:
sp = self.species.take([tt], axis=0)
if self.HAS_PROPENSITIES:
ru = self.propensities.take([tt], axis=0)
if self.HAS_XDATA:
xd = self.xdata.take([tt], axis=0)
break
output = None
if t is not None:
output = numpy.array([[temp_t[t]]])
if sp is not None:
output = numpy.hstack((output,sp))
if ra is not None:
output = numpy.hstack((output,ra))
if ru is not None:
output = numpy.hstack((output,ru))
if xd is not None:
output = numpy.hstack((output,xd))
return output
def getDataInTimeInterval(self, time, bounds=None):
"""
Returns an array of all data in interval: time-bounds <= time <= time+bounds
where bound defaults to stepsize
- *time* the interval midpoint
- *bounds* [default=None] interval halfspan defaults to stepsize
"""
temp_t = self.time.reshape(len(self.time),)
if bounds == None:
bounds = temp_t[1] - temp_t[0]
c1 = (temp_t >= time-bounds)
c2 = (temp_t <= time+bounds)
print 'Searching (%s:%s:%s)' % (time-bounds, time, time+bounds)
t = []
sp = None
ra = None
for tt in range(len(c1)):
if c1[tt] and c2[tt]:
t.append(tt)
output = None
if len(t) > 0:
output = self.time.take(t)
output = output.reshape(len(output),1)
if self.HAS_SPECIES and self.HAS_TIME:
output = numpy.hstack((output, self.species.take(t, axis=0)))
if self.HAS_PROPENSITIES:
output = numpy.hstack((output, self.propensities.take(t, axis=0)))
if self.HAS_XDATA:
output = numpy.hstack((output, self.xdata.take(t, axis=0)))
return output
def getAllSimData(self,lbls=False):
"""
Return an array of time + all available simulation data
- *lbls* [default=False] return only the data array or (data array, list of labels)
"""
labels = [self.time_label]
if self.HAS_SPECIES and self.HAS_TIME:
output = numpy.hstack((self.time, self.species))
labels += self.species_labels
if self.HAS_PROPENSITIES:
output = numpy.hstack((output, self.propensities))
labels += self.propensities_labels
if self.HAS_XDATA:
output = numpy.hstack((output, self.xdata))
labels += self.xdata_labels
if not lbls:
return output
else:
return output, labels
def getSimData(self, *args, **kwargs):
"""
Feed this method species/xdata labels and it will return an array of [time, sp1, ....]
- 'speces_l', 'xdatal' ...
- *lbls* [default=False] return only the data array or (data array, list of labels)
"""
output = self.time
if kwargs.has_key('lbls'):
lbls = kwargs['lbls']
else:
lbls = False
lout = [self.time_label]
for roc in args:
if self.HAS_SPECIES and roc in self.species_labels:
lout.append(roc)
output = numpy.hstack((output, self.species.take([self.species_labels.index(roc)], axis=-1)))
if self.HAS_PROPENSITIES and roc in self.propensities_labels:
lout.append(roc)
output = numpy.hstack((output, self.propensities.take([self.propensities_labels.index(roc)], axis=-1)))
if self.HAS_XDATA and roc in self.xdata_labels:
lout.append(roc)
output = numpy.hstack((output, self.xdata.take([self.xdata_labels.index(roc)], axis=-1)))
if not lbls:
return output
else:
return output, lout
class PysMod:
#STOMPY INSERT START
def StochSimPlot(self, plot='species', filename=None, title=None, log=None, format='points'):
"""
Plot the Stochastic simulation results, uses the new UPI pysces.plt interface:
- *plot* [default='species'] output to plot, can be one of:
- 'all' species and propensities
- 'species' species
- 'waiting_times' waiting_times
- 'propensities' propensities
- `['S1', 'R1', ]` a list of model attributes ('species')
- *filename* [default=None] if defined file is exported to filename
- *title* [default=None] the plot title
- *log* [default=None] use log axis for axis 'x', 'y', 'xy'
- *format* [default='lines'] use UPI or backend specific keys
"""
data = None
labels = None
allowedplots = ['all', 'species', 'propensities','waiting_times']
## allowedplots = ['all', 'species', 'waiting_times']
if type(plot) != list and plot not in allowedplots:
raise RuntimeError, '\nPlot must be one of %s not \"%s\"' % (str(allowedplots), plot)
if plot == 'all':
## data, labels = self.data_stochsim.getSpecies(lbls=True)
data, labels = self.data_stochsim.getAllSimData(lbls=True)
elif plot == 'species':
data, labels = self.data_stochsim.getSpecies(lbls=True)
elif plot == 'propensities':
data, labels = self.data_stochsim.getPropensities(lbls=True)
## data, labels = self.data_stochsim.getRates(lbls=True)
elif plot == 'waiting_times':
dataLst, labels = self.data_stochsim.getWaitingtimes(lbls=True)
format='points'
## data, labels = self.data_stochsim.getRates(lbls=True)
else:
plot = [at for at in plot if at in self.__species__+[self.data_stochsim.time_label]+self.data_stochsim.propensities_labels]
kwargs = {'lbls' : True}
print plot
if len(plot) > 0:
data, labels = self.data_stochsim.getSimData(*plot, **kwargs)
del allowedplots
xu = 'Time (%(multiplier)s x %(kind)s x 10**%(scale)s)**%(exponent)s' % self.__uDict__['time']
if plot == 'waiting_times':
xu = 'Inter-arrival time (%s)' % xu
xrng_start = 0.1
xrng_end = 0.1
yrng_start = 0.1
yrng_end = 0.1
for wt in range(len(dataLst)):
for d in range(len(dataLst[wt])):
D = dataLst[wt][d]
if plt.__USE_MATPLOTLIB__ and d > 0:
plt.m.hold(True)
if D != None and len(D[0]) > 0 and len(D[1]) > 0:
data = numpy.vstack([D[0], D[1]]).transpose()
if min(D[0]) < xrng_start and min(D[0]) > 0.0:
xrng_start = min(D[0])
if max(D[0]) > xrng_end:
xrng_end = max(D[0])
if min(D[1]) < yrng_start and min(D[1]) > 0.0:
yrng_start = min(D[1])
if max(D[1]) > yrng_end:
yrng_end = max(D[1])
plt.plotLines(data, 0, [1], titles=['Time']+[labels[d]], formats=[format])
plt.setRange('x', xrng_start*0.8, xrng_end*1.2)
plt.setRange('y', yrng_start*0.8, yrng_end*1.2)
if plt.__USE_MATPLOTLIB__:
plt.m.hold(False)
else:
plt.plotLines(data, 0, range(1, data.shape[1]), titles=labels, formats=[format])
# set the x-axis range so that it is original range + 0.2*sim_end
# this is a sceintifcally dtermned amount of space that is needed for the title at the
# end of the line :-) - brett 20040209
RngTime = self.data_stochsim.getTime()
end = RngTime[-1] + 0.2*RngTime[-1]
plt.setRange('x', RngTime[0], end)
del RngTime
# For now StochPy results are plotted as Amounts directly from StochPy
M = 'Amount'
## if self.__KeyWords__['Output_In_Conc']:
## M = 'Concentration'
## else:
## M = 'Amount (%(multiplier)s x %(kind)s x 10**%(scale)s)**%(exponent)s' % self.__uDict__['substance']
if plot == 'all':
yl = 'Amount, propensities'
elif plot == 'propensities':
yl = 'Propensities'
elif plot == 'waiting_times':
yl = 'Frequency'
if log == None:
log = 'xy'
elif plot == 'species':
yl = '%s' % M
else:
yl = 'User defined'
plt.setAxisLabel('x', xu)
plt.setAxisLabel('y', yl)
if log != None:
plt.setLogScale(log)
if title == None:
plt.setGraphTitle('PySCeS/StochPy simulation (' + self.ModelFile + ') ' + time.strftime("%a, %d %b %Y %H:%M:%S"))
else:
plt.setGraphTitle(title)
plt.replot()
if filename != None:
plt.export(filename, directory=self.ModelOutput, type='png')
def doStochSim(self,end=10,mode='steps',method='Direct',trajectories=1):
"""
doStochSim(end=10, mode='steps', method='Direct')
Run a stochastic simulation for until `end` is reached. This can be either steps or end time (which could be a *HUGE* number of steps).
Arguments:
- *end* [default=10] simulation end (steps or time)
- *mode* [default='steps'] simulation mode, can be one of:
- *steps* total number of steps to simulate
- *time* simulate until time is reached
- *method* [default='Direct'] stochastic algorithm, can be one of:
- Direct
- FirstReactionMethod
- NextReactionMethod
- TauLeaping
"""
if method not in ['Direct', 'FirstReactionMethod','NextReactionMethod','TauLeaping']:
print 'Method "%s" not recognised using: "Direct"' % method
method = 'Direct'
if mode not in ['steps','time']:
print 'Mode "%s" not recognised using: "steps"' % mode
mode = 'steps'
stompy_track_propensities = True
stompy_keep_psc_files = False
self.__STOMPY__.setProperty(method=method, trajectories=trajectories, interactive=True, track_propensities=stompy_track_propensities, keep_psc_files=stompy_keep_psc_files)
if mode == 'time':
self.__STOMPY__.runTimeSimulation(self, endtime=end, method=method)
else:
self.__STOMPY__.runStepSimulation(self, steps=end, method=method)
def doStochSimPlot(self, end=10.0, mode='steps', method='Direct', plot='species', fmt='points', log=None, filename=None):
"""
doStochSimPlot(end=10.0, mode='steps', method='Direct', plot='species', fmt='points', log=None, filename=None)
Run a stochastic simulation for until `end` is reached and plot the results. This can be either steps or end time (which could be a *HUGE* number of steps).
Arguments:
- *end* [default=10] simulation end (steps or time)
- *mode* [default='steps'] simulation mode, can be one of:
- *steps* total number of 'steps' to simulate
- *time* simulate until 'time' is reached
- *method* [default='Direct'] stochastic algorithm, can be one of:
- Direct
- FirstReactionMethod
- NextReactionMethod
- TauLeaping
- *plot* [default='species'] output to plot, can be one of:
- 'all' species and propensities
- 'species' species
- 'waiting_times' waiting_times
- 'propensities' propensities
- `['S1', 'R1', ]` a list of model attributes ('species')
- *filename* [default=None] if defined file is exported to filename
- *title* [default=None] the plot title
- *log* [default=None] use log axis for axis 'x', 'y', 'xy'
- *fmt* [default='lines'] use UPI or backend specific keys
"""
self.doStochSim(end=end, mode=mode, method=method,trajectories=1)
self.StochSimPlot(plot='species', filename=filename, log=log, format=fmt)
#STOMPY INSERT START
if not _HAVE_STOMPY:
def nofunc(self, *args, **kwargs):
print '\nStochastic simulation not available, please download/install *StomPy* from: http://stompy.sf.net\n'
PysMod.doStochSim = nofunc
PysMod.doStochSimPlot = nofunc
PysMod.StochSimPlot = nofunc
'''