I've created two programs to run a stochastic simulation on a system of chemical reactions. In the program I've got a function that's meant to update the elements of an array with the derivative of the changing molecule numbers popul_num and the stochastic rate constant of each reaction stoch_rate In the first program the function looks as follows:
popul_num = np.array([1.0E9, 0, 0])
stoch_rate = np.array([1.0, 0.002, 0.5, 0.04])
def update_array(popul_num, stoch_rate):
"""Specific to this model
will need to change if different model
implements equaiton 24 of the Gillespie paper"""
# calcualte in seperate varaible then pass it into the array
s_derviative = stoch_rate[1]*(2*popul_num[0] -1)/2
b = np.array([[1.0, 0.0, 0.0], [s_derviative, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.4, 0.0]])
return b
This function returns b which is an array of shape(4, 3)
In the next program I've added more reactions and more reactants and the function is as follows:
popul_num = np.array([1.0E5, 3.0E5, 0.0, 1.0E5, 0.0, 0.0, 0.0, 0.0])
stoch_rate = np.array([0.015, 0.00016, 0.5, 0.002, 0.002, 0.8])
def update_array(popul_num, stoch_rate):
"""Specific to this model
will need to change if different model
implements equaiton 24 of the Gillespie paper"""
s_derivative = stoch_rate[0]*popul_num[1]*((popul_num[1] - 1)/2) # derivative with respect to S is a function of X
x_derivative = stoch_rate[0]*popul_num[0]*((2*popul_num[1] - 1)/2) # derivative with respect to X is a function of S
r_derivative = stoch_rate[1]*((popul_num[3]*(popul_num[3]))/2)
r2_derivative = stoch_rate[2]*popul_num[4] # derivative with respect to R is a function of Y type = numpy.float64
y_derivative = stoch_rate[3]*popul_num[3] # derivative with respect to Y is a function of R type = numpy.float64
x2_derivative = stoch_rate[4]*((popul_num[1] - 1)/2)
b = np.array([[x_derivative, s_derivative, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0,
r_derivative, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, r2_derivative, y_derivative, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, stoch_rate[3], 0.0, 0.0, 0.0, 0.0], [x2_derivative, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, stoch_rate[5], 0.0]])
b.reshape((6,8))
print("Shape b:\n", b.shape)
return b
Only this returns an array of shape(6,) and I need it to be a 2D array of shape(6, 8) I've tried using the reshape() method but this results in the following error:
ValueError: cannot reshape array of size 6 into shape (6,8)
Which is thrown on the line where I call the reshape() command
I don't understand whats different about the second function meaning it doesn't return a 2D array?
Cheers
reshapecall in. I don't really understand why I'd have to call it in the first place and why thearrayisn't 2D automatically.