# 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()