I am looking for some help on solve.ivp. Currently, the way solve.ivp is currently set up, is that it returns the differentials to be used in solve ivp. However, for the purposes of my research, I will need to return the intermediate variables that are used within the model code as well.
I have tried creating an empty list, and then appending the values created by the model to the list at certain time points, however, for some unknown reason this results in way too many time points being produced and initial values of the intermediate variables being ignored. (I.e, I'm starting body weight at 60gs but the data my method provided starts it way lower, and I only want data at time points at 1,2,3,4 and 5, but printing a len statement gives me 15 data values)
Any help in solving this issue, or if you have another completely different method to try and get the variables out model would be super helpful!
Also note this method provided here: Is there a way to store intermediate values when using solve_ivp? will not work for me as I need the whole model under one function.
Here is the code I am running, (dictionaries provided)
thisdict = {
"KS" : 0.01,
"AiPc" : 0.273,
"AePc" : 0.022,
"CiPc" : 0.18,
"CePc" : 0.023,
"CAA0" : 0.70,
"POB" : 0.55,
"BV" : 0.072,
"GENDER" : 1.0,
"SPECIES" : 1.0,
"AIPV" : 0.550,
"CIPV" : 0.50,
"AEPV" : 0.0525,
"AEPV" : 0.0525,
"CEPV" : 0.063,
}
Feeding = {
"FEEDAGE0": 0.0,
"FEEDAGE1" : 7.0,
"FEEDAGE2" : 35.0,
"FEEDAGE3" : 63.0,
"FEEDAGE4" : 106.0,
"FEEDDM0" : 0.880,
"FEEDDM1" : 0.874,
"FEEDDM2" : 0.89,
"FEEDDM3" : 0.90,
"FEEDDM4" : 0.90,
"FEEDCP0" : 0.210,
"FEEDCP1" : 0.203,
"FEEDCP2" : 0.187,
"FEEDCP3" : 0.180,
"FEEDCP4" : 0.170,
"FEEDCF0" : 0.0820,
"FEEDCF1" : 0.0860,
"FEEDCF2" : 0.0940,
"FEEDCF3" : 0.0873,
"FEEDCF4" : 0.090,
"FEEDME_KJG0" : 13.00,
"FEEDME_KJG1" : 13.30,
"FEEDME_KJG2" : 13.60,
"FEEDME_KJG3" : 13.70,
"FEEDME_KJG4" : 13.70,
"FEED_SIDAA0" : 0.90,
"FEED_SIDAA1" : 0.90,
"FEED_SIDAA2" : 0.90,
"FEED_SIDAA3" : 0.90,
"FEED_SIDAA4" : 0.90,
"FEED_SIDCF0" : 0.85,
"FEED_SIDCF1" : 0.85,
"FEED_SIDCF2" : 0.85,
"FEED_SIDCF3" : 0.85,
"FEED_SIDCF4" : 0.85,
"LYSDIET0" : 1.26,
"LYSDIET1" : 1.14,
"LYSDIET2" : 1.10,
"LYSDIET3" : 0.95,
"LYSDIET4" : 0.95,
"METDIET0" : 0.516,
"METDIET1" : 0.490,
"METDIET2" : 0.460,
"METDIET3" : 0.30,
"METDIET4" : 0.30,
"DROPVALUE" : 0.97,
"MVCO2" : 1.87,
"MVO2" : 1.35,
"NH2AANH2" : 1.4,
"AIEAA" : 6.5,
"ECOST_DEAMINATION" : 2.2,
"AY_AA" : 2.5,
"AY_FA" : 9.0,
"AY_GL" : 2.0,
"GL_AA" : 0.83,
"GL_FA" : 2.67,
"GL_TG" : 0.5,
"FA_TG" : 3.0,
"FA_AY" : 0.102,
"CO2_GLAY" : 2.0,
"CO2_AYCO2" : 2.0,
"O2_GLAY" : 2.0,
"O2_AYCO2" : 2.0,
"ATP_AAAY" : 14.0,
"ATP_FAAY" : 36.0,
"ATP_GlAY" : 14.0,
"ATP_AYFA" : 1.13,
"ATP_LCFA" : 2.26,
"ATP_AYCO2" : 12.0,
"ATP_LA" : 10.0,
"ATP_PA" : 5.0,
"ATP_PC" : 1.0,
"ATP_FDFA" : 2.15,
"ATP_AUUN" : 14.36,
"EATP_kjmol" : 79.5,
"EGL_KJG" : 15.7,
"EL_KJG" : 39.6,
"EP_KJG" : 23.8,
"ELD_KJG" : 56.0,
"EPD_KJG" : 50.0,
"KF" : 0.96,
"KP" : 0.65,
"KOX" : 0.02,
"MMAA" : 125.0,
"MMCO2" : 48.0,
"MMFA" : 282.0,
"MMGL" : 180.0,
"MMN" : 14.0,
"MMNH2" : 16.0,
"MMO2" : 32.0,
"MMTG" : 884.0,
}
IBODYWT_G = 60.0
#This is the inital empty lists
BodyWt_g_S = []
t_S = []
def Broiler_model(t, pools,):
#globals are being used here in attempt to sneak intermediate variables out of the model
global ADG_gd
global BodyWt_g
global Pc
global Pv
global dAAdt
global i
global BodyWt_g_S
global t_S
global PredFI
global FaaPc_AA
# this if else statement provides a counter to try and get unintegrated values out of the model at certain points
if t == 0:
i = 0
BodyWt_g = 60.0
else:
if np.any(t == times2):
i =+ 1
Pred_FI = 80.0
FaaPc_AA = 0.001
#this is where I tried setting body weight to 60 at t==0 but it doesn't work
if t == 0.0:
BodyWt_g = 60.0
Pc, AA, Pv, = pools
#dont worry about this, this is just a feeding program
if ( (t >= Feeding.get("FEEDAGE0")) & (t < Feeding.get("FEEDAGE1")) ):
feedCP = Feeding.get("FEEDCP0")
feedCF = Feeding.get("FEEDCF0")
feedME_KJg = Feeding.get("FEEDKJG0")
feedDM = Feeding.get("FEEDDM0")
SIDAA = Feeding.get("FEED_SIDAA0") # standardized ileal dig of AA, %
SIDCF = Feeding.get("FEED_SIDCF0") # standardized ileal dig of CF, %
LYSDIET = Feeding.get("LYSDIET0")
METDIET = Feeding.get("METDIET0")
elif ( (t >= Feeding.get("FEEDAGE1")) & (t < Feeding.get("FEEDAGE2")) ):
feedCP = Feeding.get("FEEDCP1")
feedCF = Feeding.get("FEEDCF1")
feedME_KJg = Feeding.get("FEEDME_KJG1")
feedDM = Feeding.get("FEEDDM1")
SIDAA = Feeding.get("FEED_SIDAA1") # standardized ileal dig of AA, %
SIDCF = Feeding.get("FEED_SIDCF1") # standardized ileal dig of CF, %
LYSDIET = Feeding.get("LYSDIET1")
METDIET = Feeding.get("METDIET1")
elif ( (t >= Feeding.get("FEEDAGE2")) & (t < Feeding.get("FEEDAGE3")) ):
feedCP = Feeding.get("FEEDCP2")
feedCF = Feeding.get("FEEDCF2")
feedME_KJg = Feeding.get("FEEDME_KJG2")
feedDM = Feeding.get("FEEDDM2")
SIDAA = Feeding.get("FEED_SIDAA2") # standardized ileal dig of AA, %
SIDCF = Feeding.get("FEED_SIDCF2") # standardized ileal dig of CF, %
LYSDIET = Feeding.get("LYSDIET2")
METDIET = Feeding.get("METDIET2")
elif ( (t >= Feeding.get("FEEDAGE3")) & (t < Feeding.get("FEEDAGE4")) ):
feedCP = Feeding.get("FEEDCP3")
feedCF = Feeding.get("FEEDCF3")
feedME_KJg = Feeding.get("FEEDME_KJG3")
feedDM = Feeding.get("FEEDDM3")
SIDAA = Feeding.get("FEED_SIDAA3") # standardized ileal dig of AA, %
SIDCF = Feeding.get("FEED_SIDCF3") # standardized ileal dig of CF, %
LYSDIET = Feeding.get("LYSDIET3")
METDIET = Feeding.get("METDIET3")
elif ( (t >= Feeding.get("FEEDAGE4")) ):
feedCP = Feeding.get("FEEDCP4")
feedCF = Feeding.get("FEEDCF4")
feedME_KJg = Feeding.get("FEEDME_KJG4")
feedDM = Feeding.get("FEEDDM4")
SIDAA = Feeding.get("FEED_SIDAA4") # standardized ileal dig of AA, %
SIDCF = Feeding.get("FEED_SIDCF4") # standardized ileal dig of CF, %
LYSDIET = Feeding.get("LYSDIET4")
METDIET = Feeding.get("METDIET4")
#some constants that need to calcuelated
DMI = Pred_FI * feedDM
plasma_vol = thisdict.get("BV") * thisdict.get("POB") * BodyWt_g/1000.0 # blood volume in L (JE)
cAA = AA /plasma_vol # conc of amino acids (JE)
#FSRPC Variables
n1 = 0.12713
CPM1 = 1.99642
Tcaa1 = 1.0 + ((CAA0/cAA) ** n1)
#FSRPC Variables
n2 = 1.0
CPM2 = 1.19385
Tcaa2 = 1.0 + ((CAA0/cAA) ** n2)
FSRPc = thisdict.get("KS") + thisdict.get("AiPc") * (math.exp(-1.0 * thisdict.get("AePc") * t)) * (CPM1/Tcaa1) #new/cut out of other eqn
#fractional synthesis rate for carcass protein, /d (JE)#
FDRPc = thisdict.get("KS") + (thisdict.get("CiPc") * math.exp(-1.0* thisdict.get("CePc") * t) * (CPM2/Tcaa2))
#fractional degradation rate for carcass protein, /d (JE) #new/cut out of orig eqn
FaaPc_AA1 = FSRPc * Pc
#Outputs
#----------------------------------------------
FPcAA_AA = FDRPc * Pc
dPcdt = FaaPc_AA1 - FPcAA_AA
n3 = 0.124593
CPM3 = 0.1999848
Tcaa3 = 1.0 + (CAA0/cAA) ** n3
n4 = 0.124593
CPM4 = 1.9999848
Tcaa4 = 1.0 + (CAA0/cAA) ** n4
# Auxillary
#----------------------------------------------
FSRPv = thisdict.get("KS") + (thisdict.get("AIPV") * (math.exp(-1.0 * thisdict.get("AEPV")) * t)) * (CPM3/Tcaa3)
# fractional synthesis rate for visera protein, /d (JE added)
FDRPv = thisdict.get("KS") + (thisdict.get("CIPV") * (math.exp(-1.0 * thisdict.get("CEPV")) * t)) * (CPM4/Tcaa4)
# fractional degradation rate for viscera protein, d (JE added)
# Inputs
#----------------------------------------------
FAAPvPv1 = FSRPv * Pv # (JE)
# Outputs
#----------------------------------------------
FPvAA_AA = FDRPv * Pv
# Differential equation
#----------------------------------------------
dPvdt = FAAPvPv1 - FPvAA_AA #(limAAfactor routine removed)
#AA POOL
# inputs
#----------------------------------------------------
Fdiet_AA = feedCP * Pred_FI * SIDAA
FPcAA_AA = FPcAA_AA
# Outputs
#----------------------------------------------------
FAAPc_AA = FaaPc_AA
FIE_AA = Feeding.get("AIEAA")*(DMI/1000.0)
dAAdt = FPcAA_AA + Fdiet_AA + FAAPvPv1 - FAAPc_AA - FIE_AA - FPvAA_AA
Carcass = Pc + Pv
BodyWt_g = Carcass
#Here is where I append the variables that increased as counter goes up
if np.any(t == times2):
BodyWt_g_S.append(BodyWt_g)
t_S.append(t)
return [dPcdt, dAAdt, dPvdt,]
#Initial Values
CAA0 = 0.70
ieBW = IBODYWT_G * 0.87
Pc0 = 0.752 * ieBW * 0.263
AA0 = 0.002
Pv0 = 0.185*ieBW*0.263
Y1 = [Pc0]
Y2 = [AA0]
Y3 = [Pv0]
#These contain varibles the model setup, that allows me to control timestep I THINK
# times2 can also be found here
start_t = 0
end_t = 5
times2 = np.linspace(start_t,end_t, num = end_t + 1)
time_step = 1e-3 #Time step you choose
times = np.linspace(start_t,end_t, int(end_t/time_step)+1)
tsymb = sp.symbols('t')
ysymb = sp.symbols('y')
for i in range(start_t,end_t,1):
for j in range(0,int(1/time_step),1):
tstart = i + j*time_step
tend = i + j*time_step + time_step
Y0 = [Y1[-1],Y2[-1],Y3[-1],]
answer = solve_ivp(lambda tsymb,ysymb : Broiler_model(tsymb,ysymb,), (tstart,tend,), y0=Y0, method = 'RK45', t_eval=(np.linspace(tstart,tend,1000)),)
Y1.append(answer.y[0][-1])
Y2.append(answer.y[1][-1])
Y3.append(answer.y[2][-1])
