# Created by Nick Fidalgo, Chris Gatesman, Phillip Less, Brett Peters, Samantha Sandwick, and Alex Shetlerimport numpy as np # Data-fitting algorithm import matplotlib.pyplot as plt import math import scipy.optimize import pandas as pd from matplotlib import pyplot import ntpath as nt # Equations # dA = -(beta) A # dB = -(beta) B # dD = (beta) A # dE = (sigma) D + (eta) B -(phi) E #Currently unused def eulers_old(beta, sigma, eta, phi, t, xi): A = np.zeros(len(t)) B = np.zeros(len(t)) D = np.zeros(len(t)) E = np.zeros(len(t)) A[0] = 9.1 B[0] = 174000000 D[0] = 0 E[0] = 84000 for i in range(0, len(t) - 1): A[i + 1] = math.exp(480.407) * math.exp(-67.63 * i) B[i + 1] = - beta * A[i] + B[i] D[i + 1] = beta * A[i] + D[i] E[i + 1] = sigma * D[i] + eta * B[i] - phi * E[i] + E[i] return A, B, D, E # Solved equation E(t), currently unused def equation_E(sigma, beta, A0, alpha, phi, eta, t, constantE): return sigma - eta * ((beta * A0 * math.exp(-alpha * t)) / (alpha(-alpha + phi))) + constantE # Solved equation D(t), currently unused def equation_D(beta, A0, alpha, t, constantD): return (beta * A0 * math.exp(-alpha * t)) / alpha + constantD # Solved equation B(t), currently unused def equation_B(beta, A0, alpha, t, constantB): return (-beta * A0 * math.exp(-alpha * t)) / alpha + constantB # Solved equation A(t), currently USED as the function to fit A(t) def equation_A(t, A0, alpha): return A0 * np.exp(-alpha * t) """ Eulers with open parameters for fitting. Non-finished """ def eulers_E(t, sigma,v): #List of unknowns: mu, k50, kmax, sigma, w #A = 46511.17195379639 * math.exp(- 0.23492878746836654 * t) A = np.zeros(len(t)) B = np.zeros(len(t)) D = np.zeros(len(t)) E = np.zeros(len(t)) F = np.zeros(len(t)) w = 1 B0 = 174 E0= 45 B[0] = 174 E[0]=45 F[0]=1 F0 = 1 #A_tilde = A[i]/A0 ####CHECK THIS kmax = .00202 k50 = 55.33 xi =1.8741 * 10**(-4) gamma = 483 phi = .35 mu = 15.75 A0 = 46511.17195379639 kmax_hat = kmax*A0**w k50_hat = sigma*k50**w gamma_hat= gamma*(B0/E0) xi_gamma = xi*gamma*(B0/E0)*(1/sigma) xi_gamma_B0 = xi*gamma*(B0/F0)*(1/sigma) phi_hat = phi/sigma D_hat = (D0/B0) A_hat = sigma*A_tilde**w*A0**w xi_hat = xi/sigma mu_hat = mu/sigma mu_F0 = (mu*F0)/(sigma*E0) A0_hat = sigma*A0**w for i in range(0,len(t)-1): A[i] = ( A0 * math.exp(- 0.23492878746836654 * t[i])) B[i+1] = (-kmax_hat*(A[i]/A0)**w)/(k50_hat+A0_hat*(A[i]/A0)^w)+ B[i] D[i+1] = ((kmax*A0**w) * A[i]**w)/(sigma*k50**w+A[i]**w)*B[i] - sigma*D[i] + D[i] F[i+1] = sigma*gamma*(B0/E0)*D[i] + xi*(1/sigma) * gamma*(B0/E0)*B[i]-mu*(1/sigma)*F[i] + F[i] E[i+1] = mu*(1/sigma)*F[i] - phi*(1/sigma) * E[i] + E[i] return E """ Main function for fitting the equations file_name: The name of the file BEFORE .csv iteration: Currently unused but can be activated to fit several functions of several data sets at once would need slight integration to use again. Default value is 0. fig: The figure is declared outside the function so it can be used in iterations if needed ax: Same thing as fig, without a loop doing multiple iterations this is simply declared and input directly before function use """ def interpolation(file_name, iteration, fig, ax): path = r"C:\Users\AlexS\PycharmProjects\Capstone\data" + "\\" + file_name + r".csv" df = pd.read_csv(path, header=None) equation = equation_A # Distinguish from the CSV the 'x' and 'y' values from the columns t = df.iloc[:, 0] y = df.iloc[:, 1] # Bulk of the optimization popt, cov = scipy.optimize.curve_fit(equation, t, y) newA0, newAlpha = popt # Calculating r^2 residuals = y - equation(t, newA0, newAlpha) ss_res = np.sum(residuals ** 2) ss_tot = np.sum((y - np.mean(y)) ** 2) r_squared = 1 - (ss_res / ss_tot) # Begin Text Display and Graph pyplot.scatter(t, y, c="r") equation_plotting = np.linspace(-.5, max(t), num=50) pyplot.plot(equation_plotting, equation(equation_plotting, newA0, newAlpha)) # Semilog Plot Code # pyplot.semilogy(equation_plotting,equation(equation_plotting,newA0, newAlpha)) pyplot.title("Serum T-DM1 Drug Profile", fontsize=20) print("") print(newA0, "* exp(-", newAlpha, "* t)") # Square root of covariance is standard deviation print("A0 =", newA0, " +/- ", np.sqrt(cov[0, 0])) print("alpha = ", newAlpha, " +/- ", np.sqrt(cov[1, 1])) print("r^2 = ", r_squared) pyplot.legend(["Original Data", "Fitted Curve"], fontsize=18) pyplot.ylabel("Mean Serum T-DM1 concentration (ng\\mL)", fontsize=18) pyplot.xlabel("Time (days)", fontsize=18) print("") return 0 def eulers_fitting(file_name): path = r"C:\Users\AlexS\PycharmProjects\Capstone\data" + "\\" + file_name + r".csv" df = pd.read_csv(path, header=None) # Distinguish from the CSV the 'x' and 'y' values from the columns t = df.iloc[:, 0] y = df.iloc[:, 1] # Bulk of the optimization popt, cov = scipy.optimize.curve_fit(eulers_E, t, y, bounds=(0.1,6)) newsigma, newv = popt # Begin Text Display and Graph pyplot.scatter(t, y, c="r") equation_plotting = np.linspace(0.01, max(t), num=50) pyplot.plot(equation_plotting, eulers_E(equation_plotting, newsigma, newv)) # Semilog Plot Code # pyplot.semilogy(equation_plotting,equation(equation_plotting,newA0, newAlpha)) pyplot.title("E(t) function", fontsize=20) print("") print("Sigma = ", newsigma) print("v= ", newv) pyplot.legend(["Original Data", "Fitted Curve"], fontsize=18) pyplot.ylabel("ALT Concentration", fontsize=18) pyplot.xlabel("Time (days)", fontsize=18) print("") return 0 #Currently unused def interpolation_loop(): A = ["3", "6", "1_2", "2_4", "3_6", "4_8"] fig, ax = pyplot.subplots(len(A)) counter = 0 for i in A: interpolation(i, counter, fig, ax) counter += 1 pyplot.show() #Currently unused def plotting(A, B, D, E, t): xaxis = range(0, len(t) - 1) fig, ax = plt.subplots(4) ax[0].plot(A, label="A(t)-Drug Concentration") ax[0].legend() ax[1].plot(B, label="B(t)- Healthy Liver Cells") ax[1].legend() ax[2].plot(D, label="D(t) - Damaged Liver Cells") ax[2].legend() ax[3].plot(E, label="E(t) - ALT Concentration") ax[3].legend() plt.show() #Currently unused def run_simulation(): # Parameters beta = 500 sigma = 500 eta = 500 phi = 50 dt = 1 T = 5 t = np.linspace(0, T, int(T / dt) + 1) xi = 50 Equations = eulers(beta, sigma, eta, phi, t, xi) print(Equations[0]) print(Equations[1]) print(Equations[2]) print(Equations[3]) plotting(Equations[0], Equations[1], Equations[2], Equations[3], t) if __name__ == '__main__': fig, ax = pyplot.subplots(1) interpolation("Final", 0, fig, ax) #pyplot.show() #eulers_fitting("ALT") pyplot.show()