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

Delete pkOptimizer.py

parent d37fa77b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
pk_optimizer/pkOptimizer.py deleted 100644 → 0
+ 0
− 216
View file @ d37fa77b
#!/usr/bin/env python
# coding: utf-8
# In[1]:
from scipy.stats import gamma
from scipy.integrate import odeint
from scipy.optimize import fmin as fmin
import os
import csv
import re
import math as math
import numpy as np
import matplotlib.pyplot as plt
class pkOptimizer:
"""The pkOptimizer object is an optimizer for parameters in pk models."""
def __init__ (self, wd):#, Flow = 1/60, Vp = 0.05, Visf = 0.15, PS = 1/60):
"""Initializes the model with initial guess parameter values for flow, Vp, Visf, and PS.
Parameters
----------
wd : string
wd is the working directory path
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.
time : double[]
list of all timepoints
aorta : double[]
concentration of tracer in aorta (input function)
myo : double[]
concentration of tracer in myocardial tissue (Cisf)
opt : float[]
opt is the optimized parameters to fit curve
"""
self.wd = wd
#self.Flow = Flow
#self.Vp = Vp
#self.Visf = Visf
#self.PS = PS
self.time = []
self.aorta = []
self.myo = []
self.opt = []
def getData(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 (Cisf)
Returns
-------
time : double[]
list of all timepoints
aorta : double[]
concentration of tracer in aorta (input function)
myo : double[]
concentration of tracer in myocardial tissue (Cisf)
"""
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, t = np.arange(0,25), ymax = 10, tmax = 10, alpha = 2, delay = 5):
"""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 = t
ymax = ymax
tmax = tmax
alpha = alpha
delay = delay
y = np.zeros(np.size(t)); #preallocate output
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
#gamma_var_error
def MSE(self, param = [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(param) < 1:
ymax = 10; tmax = 10; alpha = 2,; delay = 5
elif len(param) < 2:
ymax = param[0]; tmax = 10; alpha = 2,; delay = 5
elif len(param) < 3:
ymax = param[0]; tmax = param[1]; alpha = 2,; delay = 5
elif len(param) < 4:
ymax = param[0]; tmax = param[1]; alpha = param[2],; delay = 5
else:
# Mean squared error (MSE) between data and gamma variate with given parameters
ymax = param[0]
tmax = param[1]
alpha = param[2]
delay = param[3]
MSE = 0
if tmax <= 0 or ymax <= 10 or delay < 0 or alpha<0 or alpha>1000 or tmax>1000:
MSE = 1000000; #just return a big number
return MSE
else:
model_vals = self.gamma_var(self.time, ymax, tmax, alpha, delay);
for i in range(len(self.myo)):
MSE = (self.aorta[i]-model_vals[i])**2 + MSE
MSE = MSE / len(self.myo)
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
self.opt = fmin(self.MSE, initGuesses, maxiter = 1000)
return self.opt.round(2)
def plotOpt(self):
"""Plots the original data to the fitted curve."""
plt.plot(self.time, self.aorta, 'bo', label='data')
#plt.plot(t, y, 'b-', label='data')
y = self.gamma_var(np.arange(0,25,0.01), self.opt[0], self.opt[1], self.opt[2], self.opt[3])
plt.plot(np.arange(0,25,0.01), y, 'r-')
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