This isn’t an ODE question, per see. It’s more of a referencing issue, I think (but I might be wrong). I’ve copied the code below. I can’t seem to sweep parameters for an ODE function that I would like to run. Any advice/insight on how to fix this problem would be greatly appreciated! Thank you very much in advance!
CONTEXTUAL CODE:
import numpy as np
from scipy.integrate import odeint # odeint allows to run ODEs
import matplotlib.pyplot as plt
def model(x, t):
# Definitions. The 'x' input is of the form [B0, G0, S0, F0]. So they need to be allocated to the appropriate variable here.
B = x[0]
G = x[1]
S = x[2]
F = x[3]
# ODEs.
dBdt = l - d*G*(B/(B+K0)) - aB*B
dGdt = pG1*(B**n/((B**n)+(K1**n))) + pG2*(S**n/((S**n)+(K2**n))) - aG*G
dSdt = pS*((F**n)/((F**n)+(K4**n))) - aS*S
dFdt = pF*((K3**n)/((K3**n)+(B**n))) - aF*F
return [dBdt, dGdt, dSdt, dFdt]
# Parameters for 'model'
# This is a list of the parameters that are used in the 'model' function above.
pG1 = 0.25
pG2 = 0.25
pF = 0.25
pS = 0.25
aB = 0.25
aG = 0.25
aF = 0.25
aS = 0.25
K0 = 0.4
K1 = 0.5
K2 = 0.3
K3 = 0.45
K4 = 0.35
n = 3
l = 0.25
d = 1.5
n1 = 3
n2 = 3
n3 = 3
n4 = 3
# Initial conditions for the ODE
model_x0 = [1,0,0,0] # this will be entered as an input to the 'model' function
#Defining the timeline of the model
model_t = np.linspace (0, 50, 200)
def sweep(param, p_low, p_high, values):
B = np.array([])
parameter_values = np.linspace(p_low, p_high, values)
for parameter_value in parameter_values:
param = parameter_value # **THIS IS THE KEY SECTION, I THINK. 'param' isn't referencing the variable that is being given in the argument of the call**
model_result = odeint(model, model_x0, model_t)
temp= np.array(model_result[:,0])
B = np.append(B, temp, axis=0)
return tuple(B)
When I test the sweep with two values for ‘pG1’ (they should give different outputs):
test = sweep(pG1, 0, 0.8, 2) test1 = test[:200] test2 = test[200:] test1==test2
This outputs True. And it shouldn’t.
Out[10]: True
Advertisement
Answer
Based on my understanding you basically want to sweep the variable in this scenario pG1 through your ode. THe primary mistake is that the ODE is not accepting the values. odeint allows for odeint(model,model_init,t, args=(a,b,c)) according to the docs. Seeing as you are initializing the parameters globally, it doesnt actually work since the variable isn’t changed in the initial ode. I am not an expert at it but I got a working version with some changes to your code. Pretty sure there is a more elegant way of doing this which I hope someone can contribute.
import numpy as np
from scipy.integrate import odeint # odeint allows to run ODEs
import matplotlib.pyplot as plt
def model(x, t,param_value): #we modify the ode model slightly to allow the passing of parameters.
pG1 = param_value['pG1'] #since its a dict we can access the parameter like this. I used a dict because if you want to sweep for all the other parameters its much easier like this I think.
# Definitions. The 'x' input is of the form [B0, G0, S0, F0]. So they need to be allocated to the appropriate variable here.
B = x[0]
G = x[1]
S = x[2]
F = x[3]
# ODEs.
dBdt = l - d*G*(B/(B+K0)) - aB*B
dGdt = pG1*(B**n/((B**n)+(K1**n))) + pG2*(S**n/((S**n)+(K2**n))) - aG*G
dSdt = pS*((F**n)/((F**n)+(K4**n))) - aS*S
dFdt = pF*((K3**n)/((K3**n)+(B**n))) - aF*F
#print(pG1) # debugged here to see if the value actually changed
return [dBdt, dGdt, dSdt, dFdt]
# Parameters for 'model'
# This is a list of the parameters that are used in the 'model' function above.
pG1 = 0.25 #Defining the parameters like this defines it once. It does not change in your original code. model takes the values here but there is no call to change it.
pG2 = 0.25
pF = 0.25
pS = 0.25
aB = 0.25
aG = 0.25
aF = 0.25
aS = 0.25
K0 = 0.4
K1 = 0.5
K2 = 0.3
K3 = 0.45
K4 = 0.35
n = 3
l = 0.25
d = 1.5
n1 = 3
n2 = 3
n3 = 3
n4 = 3
param = {'pG1':pG1,'pG2':pG2,'pF':pF,'pS':pS,'aB':aB,'aG':aG,'aF':aF,'aS':aS,'K0':K0,'K1':K1,'K2':K2,'K3':K3,'K4':K4,'n':n,'l':l,'d':d,'n1':n1,'n2':n2,'n3':n3,'n4':n4} # Here we put all your parameters defined into a nice dict.
# Initial conditions for the ODE
model_x0 = [1,0,0,0] # this will be entered as an input to the 'model' function
#Defining the timeline of the model
model_t = np.linspace (0, 50, 200)
def sweep(p_name, p_low, p_high, values, param): #note i changed the input name. So we pass the variable name we want to sweep here.
B = np.array([])
parameter_values = np.linspace(p_low, p_high, values)
for parameter_value in parameter_values:
param[p_name] = parameter_value # Here we use the linspace values to 'sweep' the parameter value that we are manipulating
model_result = odeint(model, model_x0, model_t,args = (param,)) #here we pass the 'new' parameters into the actual ode.
temp= np.array(model_result[:,0])
B = np.append(B, temp, axis=0)
return tuple(B)
test = sweep('pG1', 0, 0.8, 2,param) #so we pass the 'default' parameters
test1 = test[:200]
test2 = test[200:]
print(test1==test2)
>>False
This is definitely a hacky way to do it but it works. I am sure someone with more experience using odeint and the numpy/scipy package can give you an easier/cleaner way to do this. You can extend this for all the parameters if you so wish. I would not recommend globals however.