Module pk_one_comp
The pk1Comp object is a one compartment PK model that outputs graphs of mass of tracer over time.
Expand source code
#!/usr/bin/env python
# coding: utf-8
# In[1]:
"""The pk1Comp object is a one compartment PK
model that outputs graphs of mass of tracer over time."""
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_one_comp:
"""The pk1Comp object is a one compartment PK
model that outputs graphs of mass 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 :
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
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(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)
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 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, init_guesses):
"""Uses fmin algorithm (Nelder-Mead simplex algorithm) to
minimize loss function (MSE) of input function.
Parameters
----------
init_guesses : 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, init_guesses, 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
-------
dc_dt : double[]
contains the derivative of concentrations with respect to time.
"""
# Define value of input function conc_in
conc_in = self.gamma_var(time, self.ymax, self.tmax, self.alpha, self.delay)
# Unpack the current values of the variables we wish to "update" from the curr_vals list
myo_vals = curr_vals
# Right-hand side of odes, which are used to computer the derivative
dc_dt = self.flow*(conc_in - myo_vals)/self.vol_plasma
return dc_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.vol_plasma = guess[1]
self.baseline = guess[2]
mse = 0
if self.flow <= 0 or self.flow >= 25 or self.vol_plasma > 100 or self.vol_plasma < 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], 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] - sol.y[0][i])**2 + mse
mse = mse / len(self.myo)
return mse
def output_func_fit(self, init_guesses):
"""Uses fmin algorithm (Nelder-Mead simplex algorithm) to minimize
loss function (MSE) of output function.
Parameters
----------
init_guesses : 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, init_guesses, maxiter=10000)
self.flow = opt1[0].round(4)
self.vol_plasma = opt1[1].round(4)
self.baseline = opt1[2].round(4)
return opt1.round(4)
def main(self):
"""Run this for demonstration of what pk_one_comp can do"""
# 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])
print(opt)
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([.03102, .462, self.myo[0]])
print(opt2)
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_one_comp (wd=WindowsPath('C:/Users/Ethan/OneDrive - Michigan State University/MSU/Classwork/Computational Modeling/Models/Data'), filename='CTPERF005_stress.csv')
-
The pk1Comp object is a one compartment PK model that outputs graphs of mass 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 : 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/(mLmin). 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/(gmin). 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_one_comp: """The pk1Comp object is a one compartment PK model that outputs graphs of mass 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 : 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 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(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) 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 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, init_guesses): """Uses fmin algorithm (Nelder-Mead simplex algorithm) to minimize loss function (MSE) of input function. Parameters ---------- init_guesses : 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, init_guesses, 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 ------- dc_dt : double[] contains the derivative of concentrations with respect to time. """ # Define value of input function conc_in conc_in = self.gamma_var(time, self.ymax, self.tmax, self.alpha, self.delay) # Unpack the current values of the variables we wish to "update" from the curr_vals list myo_vals = curr_vals # Right-hand side of odes, which are used to computer the derivative dc_dt = self.flow*(conc_in - myo_vals)/self.vol_plasma return dc_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.vol_plasma = guess[1] self.baseline = guess[2] mse = 0 if self.flow <= 0 or self.flow >= 25 or self.vol_plasma > 100 or self.vol_plasma < 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], 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] - sol.y[0][i])**2 + mse mse = mse / len(self.myo) return mse def output_func_fit(self, init_guesses): """Uses fmin algorithm (Nelder-Mead simplex algorithm) to minimize loss function (MSE) of output function. Parameters ---------- init_guesses : 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, init_guesses, maxiter=10000) self.flow = opt1[0].round(4) self.vol_plasma = opt1[1].round(4) self.baseline = opt1[2].round(4) return opt1.round(4) def main(self): """Run this for demonstration of what pk_one_comp can do""" # 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]) print(opt) 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([.03102, .462, self.myo[0]]) print(opt2) 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
dc_dt
:double
[]- contains the derivative of concentrations 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 ------- dc_dt : double[] contains the derivative of concentrations with respect to time. """ # Define value of input function conc_in conc_in = self.gamma_var(time, self.ymax, self.tmax, self.alpha, self.delay) # Unpack the current values of the variables we wish to "update" from the curr_vals list myo_vals = curr_vals # Right-hand side of odes, which are used to computer the derivative dc_dt = self.flow*(conc_in - myo_vals)/self.vol_plasma return dc_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) return y 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
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 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(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, init_guesses)
-
Uses fmin algorithm (Nelder-Mead simplex algorithm) to minimize loss function (MSE) of input function. Parameters
init_guesses
: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, init_guesses): """Uses fmin algorithm (Nelder-Mead simplex algorithm) to minimize loss function (MSE) of input function. Parameters ---------- init_guesses : 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, init_guesses, 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)
-
Run this for demonstration of what pk_one_comp can do
Expand source code
def main(self): """Run this for demonstration of what pk_one_comp can do""" # 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]) print(opt) 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([.03102, .462, self.myo[0]]) print(opt2) 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, init_guesses)
-
Uses fmin algorithm (Nelder-Mead simplex algorithm) to minimize loss function (MSE) of output function. Parameters
init_guesses
: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, init_guesses): """Uses fmin algorithm (Nelder-Mead simplex algorithm) to minimize loss function (MSE) of output function. Parameters ---------- init_guesses : 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, init_guesses, maxiter=10000) self.flow = opt1[0].round(4) self.vol_plasma = 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.vol_plasma = guess[1] self.baseline = guess[2] mse = 0 if self.flow <= 0 or self.flow >= 25 or self.vol_plasma > 100 or self.vol_plasma < 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], 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] - sol.y[0][i])**2 + mse mse = mse / len(self.myo) return mse