Skip to content
Snippets Groups Projects
Commit 67f4ef34 authored by Tu, Ethan's avatar Tu, Ethan
Browse files

Delete pk2Comp.py

parent 78df2f91
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 0 additions and 401 deletions
pk_optimizer/pk2Comp.py deleted 100644 → 0
+ 0
− 401
View file @ 78df2f91
"""The pk2Comp object is a two compartment PK model
that outputs graphs of concentration of tracer over time.
"""
#!/usr/bin/env python
# coding: utf-8
# In[1]:
import pathlib
import os
import csv
import re
import math
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from scipy.optimize import fmin
# %matplotlib inline
# np.set_printoptions(threshold=sys.maxsize)
class pk2Comp:
"""The pk2Comp object is a two compartment PK model
that outputs graphs of concentration of tracer over time.
"""
def __init__(self, wd=pathlib.Path('Data').absolute(), filename='CTPERF005_stress.csv'):
"""Initializes the model with default parameter values for flow, Vp, Visf, and PS.
Parameters
----------
time : double[]
list of all timepoints
aorta : double[]
concentration of tracer in aorta (input function)
myo : double[]
concentration of tracer in myocardial tissue (conc_isf)
Flow : double
Flow is the flow of plasma through the blood vessel in mL/(mL*min). Defaults to 1/60.
Vp : double
Vp is the volume of plasma in mL. Defaults to 0.05.
Visf : double
Visf is the volume of interstitial fluid in mL. Defaults to 0.15.
PS : double
PS is the permeability-surface area constant in mL/(g*min). Defaults to 1/60.
"""
# Subject Data
self.wd = wd
self.filename = filename
self.time = []
self.aorta = []
self.myo = []
# Declare Variables for initial conditions
# Compartment variables to be fitted
self.flow = 1/60
self.visf = 0.15
self.baseline = 60
# Other Compartmental Modelvariables
self.perm_surf = 0.35
self.vol_plasma = 0.10
# Solved ode
self.sol = []
# Gamma variables
self.ymax = 250
self.tmax = 6.5
self.alpha = 2.5
self.delay = 0
self.f2 = np.array([])
self.f3 = np.array([])
def Get_data(self, filename):
"""Imports data from all .csv files in directory.
Parameters
----------
wd : str
wd is the working directory path
Attributes
----------
time : double[]
list of all timepoints
aorta : double[]
concentration of tracer in aorta (input function)
myo : double[]
concentration of tracer in myocardial tissue (conc_isf)
Returns
-------
time : double[]
list of all timepoints
aorta : double[]
concentration of tracer in aorta (input function)
myo : double[]
concentration of tracer in myocardial tissue (conc_isf)
"""
os.chdir(self.wd)
# File not found error
if not os.path.isfile(filename):
raise ValueError(
"Input file does not exist: {0}. I'll quit now.".format(filename))
data = list(csv.reader(open(filename), delimiter='\t'))
for i in range(12):
self.time.append(
float(re.compile('\d+[.]+\d+|\d+').findall(data[i+1][0])[0]))
self.aorta.append(
float(re.compile('\d+[.]+\d+|\d+').findall(data[i+1][1])[0]))
self.myo.append(
float(re.compile('\d+[.]+\d+|\d+').findall(data[i+1][2])[0]))
return self.time, self.aorta, self.myo
# gamma_var distribution curve
def gamma_var(self, time=np.arange(0, 25), ymax=10, tmax=10, alpha=2, delay=0):
"""Creates a gamma variate probability density function with given alpha,
location, and scale values.
Parameters
----------
t : double[]
array of timepoints
ymax : double
peak y value of gamma distribution
tmax : double
location of 50th percentile of function
alpha : double
scale parameter
delay : double
time delay of which to start gamma distribution
Returns
-------
y : double[]
probability density function of your gamma variate.
"""
# Following Madsen 1992 simplified parameterization for gamma variate
t = time
self.ymax = ymax
self.tmax = tmax
self.alpha = alpha
self.delay = delay
y = np.zeros(np.size(t)) # preallocate output
# For odeint, checks if t input is array or float
if isinstance(t, (list, np.ndarray)):
for i in range(np.size(y)):
if t[i] < delay:
y[i] = 0
else:
y[i] = round((ymax*tmax**(-alpha)*math.exp(alpha))*(t[i]-delay)
** alpha*math.exp(-alpha*(t[i]-delay)/tmax), 3)
return y
else:
y = round((ymax*tmax**(-alpha)*math.exp(alpha))*(t-delay)
** alpha*math.exp(-alpha*(t-delay)/tmax), 3)
return y
# gamma_var_error
def inputMSE(self, guess=[10, 10, 2, 5]):
"""Calculates Mean squared error (MSE) between data and
gamma variate with given parameters values.
Parameters
----------
param : ndarray[]
time : double[]
array of timepoints
ymax : double
peak y value of gamma distribution
tmax : double
location of 50th percentile of function
alpha : double
scale parameter
delay : double
time delay of which to start gamma distribution
Returns
-------
MSE : double
Mean squared error
"""
if len(guess) < 1:
self.ymax = 10
self.tmax = 10
self.alpha = 2
self.delay = 5
elif len(guess) < 2:
self.ymax = guess[0]
self.tmax = 10
self.alpha = 2
self.delay = 5
elif len(guess) < 3:
self.ymax = guess[0]
self.tmax = guess[1]
self.alpha = 2
self.delay = 5
elif len(guess) < 4:
self.ymax = guess[0]
self.tmax = guess[1]
self.alpha = guess[2]
self.delay = 5
else:
# Mean squared error (MSE) between data and gamma variate with given parameters
self.ymax = guess[0]
self.tmax = guess[1]
self.alpha = guess[2]
self.delay = guess[3]
mse = 0
if self.tmax <= 0 or self.ymax <= 10 or self.delay < 0 or self.alpha < 0 \
or self.alpha > 1000 or self.tmax > 1000:
mse = 1000000 # just return a big number
else:
model_vals = self.gamma_var(
self.time, self.ymax, self.tmax, self.alpha, self.delay)
for i in range(len(self.aorta)):
mse = (self.aorta[i] - model_vals[i])**2 + mse
mse = mse / len(self.aorta)
return round(mse, 3)
def inputFuncFit(self, initGuesses):
"""Uses fmin algorithm (Nelder-Mead simplex algorithm) to
minimize loss function (MSE) of input function.
Parameters
----------
initGuesses : ndarray[]
Array of initial guesses containing:
time : double[]
array of timepoints
ymax : double
peak y value of gamma distribution
tmax : double
location of 50th percentile of function
alpha : double
scale parameter
delay : double
time delay of which to start gamma distribution
Returns
-------
opt : double[]
optimized parameters
"""
# Mean squared error (MSE) between data and gamma variate with given parameters
opt = fmin(self.inputMSE, initGuesses, maxiter=1000)
self.ymax = opt[0]
self.tmax = opt[1]
self.alpha = opt[2]
self.delay = opt[3]
return opt.round(2)
# Derivative function
def derivs(self, time, curr_vals):
"""Finds derivatives of ODEs.
Parameters
----------
curr_vals : double[]
curr_vals it he current values of the variables we wish to
"update" from the curr_vals list.
time : double[]
time is our time array from 0 to tmax with timestep dt.
Returns
-------
dconc_plasma_dt : double[]
contains the derivative of concentration in plasma with respect to time.
dconc_isf_dt : double[]
contains the derivative of concentration in interstitial fluid with respect to time.
"""
# Unpack the current values of the variables we wish to "update" from the curr_vals list
conc_plasma, conc_isf = curr_vals
# Define value of input function conc_in
conc_in = self.gamma_var(time, self.ymax, self.tmax,
self.alpha, self.delay)
# Right-hand side of odes, which are used to computer the derivative
dconc_plasma_dt = (self.flow/self.vol_plasma)*(conc_in - conc_plasma) + (self.perm_surf/self.vol_plasma)*(conc_isf - conc_plasma)
dconc_isf_dt = (self.perm_surf/self.visf)*(conc_plasma - conc_isf)
return dconc_plasma_dt, dconc_isf_dt
def outputMSE(self, guess):
"""Calculates Mean squared error (MSE) between data and
gamma variate with given parameters values.
Parameters
----------
guess : ndarray[]
Flow : double
Flow is the flow of plasma through the blood vessel in mL/(mL*min).
Defaults to 1/60.
Vp : double
Vp is the volume of plasma in mL. Defaults to 0.05.
Visf : double
Visf is the volume of interstitial fluid in mL. Defaults to 0.15.
PS : double
PS is the permeability-surface area constant in mL/(g*min). Defaults to 1/60.
Returns
-------
MSE : double
Mean squared error
"""
self.flow = guess[0]
self.visf = guess[1]
self.baseline = guess[2]
mse = 0
if self.flow <= 0 or self.flow >= 25 or self.visf > 100 or self.visf < 0 \
or self.baseline > 150 or self.baseline < 0:
mse = 100000 # just return a big number
else:
sol = solve_ivp(self.derivs, [0, 30], [0, 0], t_eval=self.time)
MBF = sol.y[0] + sol.y[1]
temp = np.asarray(self.myo) - self.baseline
for i in range(len(self.myo)):
mse = (temp[i] - MBF[i])**2 + mse
mse = mse / len(self.myo)
return mse
def outputFuncFit(self, initGuesses):
"""Uses fmin algorithm (Nelder-Mead simplex algorithm) to minimize
loss function (MSE) of input function.
Parameters
----------
initGuesses : ndarray[]
Array of initial guesses containing:
time : double[]
array of timepoints
ymax : double
peak y value of gamma distribution
tmax : double
location of 50th percentile of function
alpha : double
scale parameter
delay : double
time delay of which to start gamma distribution
Returns
-------
opt : double[]
optimized parameters
"""
# Mean squared error (MSE) between data and gamma variate with given parameters
opt1 = fmin(self.outputMSE, initGuesses, maxiter=10000)
self.flow = opt1[0]
self.visf = opt1[1]
self.baseline = opt1[2]
return opt1 # .round(4)
def main(self):
# Gets data from file
self.Get_data(self.filename)
# Plots original data
plt.plot(self.time, self.aorta, 'bo')
plt.plot(self.time, self.myo, 'ro')
# Fit gamma_var input function and plots it
opt = self.inputFuncFit([250, 7, 4, 0])
print(opt)
print(self.ymax, self.tmax, self.alpha, self.delay)
plt.plot(np.arange(0, 25, 0.01), self.gamma_var(np.arange(0, 25, 0.01),
opt[0], opt[1], opt[2], opt[3]), 'k-')
# Fit uptake function and plot it
opt2 = self.outputFuncFit([.011418, .62, self.myo[0]])
print('myo is ', self.myo[0])
print(opt2)
print(self.flow, self.visf, self.baseline)
print('time is ', self.time)
self.f2 = solve_ivp(self.derivs, [0, 30], [0, 0], t_eval=self.time)
self.f3 = self.f2.y[0] + self.f2.y[1]
plt.plot(self.time, self.f3 + self.baseline, 'm-')
# In[2]:
PK = pk2Comp()
PK.main()
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment