Source code for pysces.PyscesLink

"""
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











'''
_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

'''