Module pk_two_comp
The pk2Comp object is a two compartment PK model that outputs graphs of concentration of tracer over time.
Expand source code
#!/usr/bin/env python
# coding: utf-8
# In[1]:
"""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
class pk_two_comp:
"""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
----------
wd : path
Absolute path name to current directory. Defaults to ./Data
filename : String
Name of the data file you'd like to access
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)
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.
ymax : int
Magnitude of Gamma-var peak.
tmax : double
Time of which highest peak in Gamma-var appears
alpha : double
Scale factor
delay : double
Delay to start Gamma-var curve.
"""
# Subject Data
if os.path.basename(os.path.normpath(pathlib.Path().absolute())) != 'Data':
self.wd = pathlib.Path('Data').absolute()
else:
self.wd = wd
if not isinstance(filename, str):
raise ValueError('Filename must be a string')
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.deriv_sol = np.array([])
self.fit_myo = np.array([])
def get_data(self, filename):
"""Imports data from all .csv files in directory.
Parameters
----------
filename : string
Name of the data file you'd like to access
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)
"""
self.time = []
self.aorta = []
self.myo = []
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(r'\d+[.]+\d+|\d+').findall(data[i+1][0])[0]))
self.aorta.append(
float(re.compile(r'\d+[.]+\d+|\d+').findall(data[i+1][1])[0]))
self.myo.append(
float(re.compile(r'\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)
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 input_mse(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 input_func_fit(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.input_mse, 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 output_mse(self, guess):
"""Calculates Mean squared error (MSE) between data and
output derivs 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 output_func_fit(self, initGuesses):
"""Uses fmin algorithm (Nelder-Mead simplex algorithm) to minimize
loss function (MSE) of output function.
Parameters
----------
initGuesses : ndarray[]
Array of initial guesses containing:
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
-------
opt : double[]
optimized parameters
"""
# Mean squared error (MSE) between data and gamma variate with given parameters
opt1 = fmin(self.output_mse, initGuesses, maxiter=10000)
self.flow = opt1[0].round(4)
self.visf = opt1[1].round(4)
self.baseline = opt1[2].round(4)
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.input_func_fit([250, 7, 4, 0])
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.output_func_fit([.011418, .62, self.myo[0]])
self.deriv_sol = solve_ivp(self.derivs, [0, 30],
[0, 0], t_eval=self.time)
self.fit_myo = self.deriv_sol.y[0] + self.deriv_sol.y[1]
plt.plot(self.time, self.fit_myo + self.baseline, 'm-')
Classes
class pk_two_comp (wd=WindowsPath('C:/Users/Ethan/OneDrive - Michigan State University/MSU/Classwork/Computational Modeling/Models/Data'), filename='CTPERF005_stress.csv')
-
The pk2Comp object is a two compartment PK model that outputs graphs of concentration of tracer over time.
Initializes the model with default parameter values for flow, Vp, Visf, and PS. Parameters
wd
:path
- Absolute path name to current directory. Defaults to ./Data
filename
:String
- Name of the data file you'd like to access
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)
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.
ymax
:int
- Magnitude of Gamma-var peak.
tmax
:double
- Time of which highest peak in Gamma-var appears
alpha
:double
- Scale factor
delay
:double
- Delay to start Gamma-var curve.
Expand source code
class pk_two_comp: """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 ---------- wd : path Absolute path name to current directory. Defaults to ./Data filename : String Name of the data file you'd like to access 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) 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. ymax : int Magnitude of Gamma-var peak. tmax : double Time of which highest peak in Gamma-var appears alpha : double Scale factor delay : double Delay to start Gamma-var curve. """ # Subject Data if os.path.basename(os.path.normpath(pathlib.Path().absolute())) != 'Data': self.wd = pathlib.Path('Data').absolute() else: self.wd = wd if not isinstance(filename, str): raise ValueError('Filename must be a string') 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.deriv_sol = np.array([]) self.fit_myo = np.array([]) def get_data(self, filename): """Imports data from all .csv files in directory. Parameters ---------- filename : string Name of the data file you'd like to access 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) """ self.time = [] self.aorta = [] self.myo = [] 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(r'\d+[.]+\d+|\d+').findall(data[i+1][0])[0])) self.aorta.append( float(re.compile(r'\d+[.]+\d+|\d+').findall(data[i+1][1])[0])) self.myo.append( float(re.compile(r'\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) 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 input_mse(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 input_func_fit(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.input_mse, 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 output_mse(self, guess): """Calculates Mean squared error (MSE) between data and output derivs 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 output_func_fit(self, initGuesses): """Uses fmin algorithm (Nelder-Mead simplex algorithm) to minimize loss function (MSE) of output function. Parameters ---------- initGuesses : ndarray[] Array of initial guesses containing: 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 ------- opt : double[] optimized parameters """ # Mean squared error (MSE) between data and gamma variate with given parameters opt1 = fmin(self.output_mse, initGuesses, maxiter=10000) self.flow = opt1[0].round(4) self.visf = opt1[1].round(4) self.baseline = opt1[2].round(4) 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.input_func_fit([250, 7, 4, 0]) 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.output_func_fit([.011418, .62, self.myo[0]]) self.deriv_sol = solve_ivp(self.derivs, [0, 30], [0, 0], t_eval=self.time) self.fit_myo = self.deriv_sol.y[0] + self.deriv_sol.y[1] plt.plot(self.time, self.fit_myo + self.baseline, 'm-')
Methods
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.
Expand source code
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 gamma_var(self, time=array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24]), 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.
Expand source code
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) else: y = round((ymax*tmax**(-alpha)*math.exp(alpha))*(t-delay) ** alpha*math.exp(-alpha*(t-delay)/tmax), 3) return y
def get_data(self, filename)
-
Imports data from all .csv files in directory. Parameters
filename
:string
- Name of the data file you'd like to access
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)
Expand source code
def get_data(self, filename): """Imports data from all .csv files in directory. Parameters ---------- filename : string Name of the data file you'd like to access 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) """ self.time = [] self.aorta = [] self.myo = [] 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(r'\d+[.]+\d+|\d+').findall(data[i+1][0])[0])) self.aorta.append( float(re.compile(r'\d+[.]+\d+|\d+').findall(data[i+1][1])[0])) self.myo.append( float(re.compile(r'\d+[.]+\d+|\d+').findall(data[i+1][2])[0])) return self.time, self.aorta, self.myo
def input_func_fit(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
Expand source code
def input_func_fit(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.input_mse, initGuesses, maxiter=1000) self.ymax = opt[0] self.tmax = opt[1] self.alpha = opt[2] self.delay = opt[3] return opt.round(2)
def input_mse(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
Expand source code
def input_mse(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 main(self)
-
Expand source code
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.input_func_fit([250, 7, 4, 0]) 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.output_func_fit([.011418, .62, self.myo[0]]) self.deriv_sol = solve_ivp(self.derivs, [0, 30], [0, 0], t_eval=self.time) self.fit_myo = self.deriv_sol.y[0] + self.deriv_sol.y[1] plt.plot(self.time, self.fit_myo + self.baseline, 'm-')
def output_func_fit(self, initGuesses)
-
Uses fmin algorithm (Nelder-Mead simplex algorithm) to minimize loss function (MSE) of output function. Parameters
initGuesses
:ndarray
[]-
Array of initial guesses containing: 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
opt
:double
[]- optimized parameters
Expand source code
def output_func_fit(self, initGuesses): """Uses fmin algorithm (Nelder-Mead simplex algorithm) to minimize loss function (MSE) of output function. Parameters ---------- initGuesses : ndarray[] Array of initial guesses containing: 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 ------- opt : double[] optimized parameters """ # Mean squared error (MSE) between data and gamma variate with given parameters opt1 = fmin(self.output_mse, initGuesses, maxiter=10000) self.flow = opt1[0].round(4) self.visf = opt1[1].round(4) self.baseline = opt1[2].round(4) return opt1.round(4)
def output_mse(self, guess)
-
Calculates Mean squared error (MSE) between data and output derivs 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
Expand source code
def output_mse(self, guess): """Calculates Mean squared error (MSE) between data and output derivs 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