How to include numbers we need in a list which is generated by some stochastic algorithm

I need to implement a stochastic algorithm that provides as output the times and the states at the corresponding time points of a dynamic system. We include randomness in defining the time points by retrieving a random number from the uniform distribution. What I want to do, is to find the state at time points 0,1,2,…,24. Given the randomness of the algorithm, the time points 1, 2, 3,…,24 are not necessarily hit. We my include rounding at two decimal places, but even with rounding I can not find/insert all of these time points. The question is, how to change the code so as to be able to include in the list of the time points the numbers 1, 2,…, 24 while preserving the stochasticity of the algorithm ? Thanks for any suggestion.

import numpy as np
import random
import math as m

np.random.seed(seed = 5)

# Stoichiometric matrix
S = np.array([(-1, 0), (1, -1)])

# Reaction parameters
ke = 0.3; ka = 0.5
k = [ke, ka]

# Initial state vector at time t0

X1 = [200]; X2 = [0]

# We will update it for each time.
X = [X1, X2]

# Initial time is t0 = 0, which we will update.
t = [0]
# End time
tfinal = 24

# The propensity vector R concerning the last/updated value of time
def ReactionRates(k, X1, X2):
    R = np.zeros((2,1))
    R[0] = k[1] * X1[-1]
    R[1] = k[0] * X2[-1] 
    return R   

# We implement the Gillespie (SSA) algorithm

while True:

    # Reaction propensities/rates
    R = ReactionRates(k,X1,X2)  
    propensities =  R
    propensities_sum = sum(R)[0]
    if propensities_sum == 0:
        break
    # we include randomness 
    u1 = np.random.uniform(0,1) 

    delta_t = (1/propensities_sum) * m.log(1/u1)
    if t[-1] + delta_t > tfinal:
        break
    t.append(t[-1] + delta_t)

    b = [0,R[0], R[1]]
    u2 = np.random.uniform(0,1)

# Choose j
    lambda_u2 = propensities_sum * u2
    for j in range(len(b)):
        if sum(b[0:j-1+1]) < lambda_u2 <= sum(b[1:j+1]):
        break # out of for j
    # make j zero based
    j -= 1
# We update the state vector
    X1.append(X1[-1] + S.T[j][0])
    X2.append(X2[-1] + S.T[j][1])

# round t values
t = [round(tt,2) for tt in t]
print("The time steps:", t)
print("The second component of the state vector:", X2)

Answer

After playing with your model, I conclude, that interpolation works fine.

Basically, just append the following lines to your code:

ts = np.arange(tfinal+1)
xs = np.interp(ts, t, X2)

and if you have matplotlib installed, you can visualize using

import matplotlib.pyplot as plt

plt.plot(t, X2)
plt.plot(ts, xs)
plt.show()

Advertisement

Answer