import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

# Set a seed for the random number generator so we get the same random numbers each time
np.random.seed(42)

# Create fake x-data
x = np.arange(10)
# Create fake y-data
a = 4.5
b = 0.5
c = 20
y = a * np.exp(-b * x) + c 
y = y + np.random.normal(loc=0, scale=0.3, size=len(x))  # Add noise

def plot_sc(x,y):
    fig, ax = plt.subplots()
    ax.scatter(x, y, c="green")
    ax.set_xlabel("X")
    ax.set_ylabel("y")
    #ax.legend()
    plt.show()

plot_sc(x,y)

# Set a seed for the random number generator so we get the same random numbers each time
np.random.seed(6354)

# Create fake x-data
x1 = np.arange(10)
# Create fake y-data
a = 4.5
b = 0.5
c = 0  #set c=0
y1 = a * np.exp(b * x) 
y1 = y1 + np.random.normal(loc=0, scale=10, size=len(x))  # Add noise

plot_sc(x1,y1)

# Fit a polynomial of degree 1 (a linear function) to the data
p = np.polyfit(x1, np.log(y1), 1)

# Convert the polynomial back into an exponential
a = np.exp(p[1])
b = p[0]
x1_fit = np.linspace(np.min(x1), np.max(x1), 100)
y1_fit = a*np.exp(b * x1_fit)

ax = plt.axes()
ax.scatter(x1, y1, label='Raw data')
ax.plot(x1_fit, y1_fit, 'k', label='Fitted curve')
ax.set_title('Using polyfit() to fit an exponential function')
ax.set_ylabel('y')
ax.set_ylim(0, 500)
ax.set_xlabel('X')
ax.legend()

# Fit a weighted polynomial of degree 1 (a linear function) to the data
p = np.polyfit(x1, np.log(y1), 1, w=np.sqrt(y1))

# Convert the polynomial back into an exponential
a = np.exp(p[1])
b = p[0]
x1_fit_wt = np.linspace(np.min(x1), np.max(x1), 100)
y1_fit_wt = a * np.exp(b * x1_fit_wt)

# Plot
ax = plt.axes()
ax.scatter(x1, y1, label='Raw data')
ax.plot(x1_fit, y1_fit, 'k', label='Fitted curve, unweighted')
ax.plot(x1_fit_wt, y1_fit_wt, 'k--', label='Fitted curve, weighted')
ax.set_title('Using polyfit() to fit an exponential function')
ax.set_ylabel('y')
ax.set_ylim(0, 500)
ax.set_xlabel('X')
ax.legend()

# Set a seed for the random number generator so we get the same random numbers each time
np.random.seed(6354)

# Create fake x-data
x2 = np.arange(10)
# Create fake y-data
a = 4.5
b = 0.5
c = 50  
y2 = a * np.exp(b * x) + c
y2 = y2 + np.random.normal(loc=0, scale=10, size=len(x))  # Add noise

plot_sc(x2,y2)

# Have an initial guess as to what the values of the parameters are
a0 = 5
b0 = 0.6
c0 = 40

from scipy.optimize import curve_fit

# Fit the function a * np.exp(b * t) + c to x and y
popt, pcov = curve_fit(lambda t, a, b, c: a * np.exp(b * t) + c, x2, y2, p0=(a0,b0,c0), maxfev=5000)

# Create the fitted curve
x2_fit = np.linspace(np.min(x2), np.max(x2), 100)
y2_fit = popt[0] * np.exp(popt[1] * x2_fit) + popt[2]

ax = plt.axes()
ax.scatter(x2, y2, label='Raw data')
ax.plot(x2_fit, y2_fit, 'k', label='Fitted curve')
ax.set_title('Using curve_fit() to fit an exponential function')
ax.set_ylabel('y')
ax.set_ylim(0, 500)
ax.set_xlabel('X')
ax.legend()

print(popt)

# Set a seed for the random number generator so we get the same random numbers each time
np.random.seed(64)

# Create fake x-data
x3 = np.arange(10)
# Create fake y-data
a = 4.5
b = - 0.5
c = 10
y3 = a * np.exp(b * x) + c
y3 = y3 + np.random.normal(loc=0, scale=0.1, size=len(x3))  # Add noise

plot_sc(x3,y3)

# Have an initial guess as to what the values of the parameters are
a0 = 5
b0 = - 0.6
c0 = 12

# Fit the function a * np.exp(b * t) + c to x and y
popt, pcov = curve_fit(lambda t, a, b, c: a * np.exp(b * t) + c, x3, y3, p0=(a0,b0,c0), maxfev=5000)

# Create the fitted curve
x3_fit = np.linspace(np.min(x3), np.max(x3), 100)
y3_fit = popt[0] * np.exp(popt[1] * x3_fit) + popt[2]

ax = plt.axes()
ax.scatter(x3, y3, label='Raw data')
ax.plot(x3_fit, y3_fit, 'k', label='Fitted curve')
ax.set_title('Fitting an negative exponential function')
ax.set_ylabel('y')
ax.set_ylim(9.5, 15.5)
ax.set_xlabel('X')
ax.legend()

popt

def func(x,a,b,c):
    return a * np.exp(b * x) + c

from lmfit import Model

model = Model(func)

model.independent_vars

model.param_names

params = model.make_params(a=5.0, b=-0.6, c=12)
result = model.fit(y3, params, x=x3)

print(result.fit_report())

y4 = np.array([146., 94., 79.,62., 41.])
x4 = np.array([375, 470, 528, 625,880])

def fun(x, alpha1, alpha2, A, B):
    return A*x**(-alpha1)

uncertainty = abs(0.16 + np.random.normal(size=x.size, scale=0.05))
uncertainty 

from lmfit import minimize, Parameters

def residual(params, x, data):
    alpha1 = params['alpha1']
    A = params['A']

    
    model = A*x**(-alpha1)
    
    return (data-model)

params = Parameters()
params.add('alpha1', value=1.0, min=0.8, max=10.0)

params.add('A', value=100.0, min=100.0, max=1.0e10)

out = minimize(residual, params, args=(x4, y4))

out.params

x4_fit = np.linspace(np.min(x4), np.max(x4), 100)
y4_fit = 15423.092*x4_fit**(-1.40000000) + 1.0000e+29*x4_fit**(-13.0000000)

ax = plt.axes()
ax.scatter(x4, y4, label='Raw data')
ax.plot(x4_fit, y4_fit, 'k', label='Fitted curve')
ax.set_title('Fitting an power curve')
ax.set_ylabel('y')
ax.set_xlabel('X')
ax.legend()


y4_B = 1.0000e+10*x4_fit**(-10.0000000)

ax = plt.axes()
ax.scatter(x4, y4, label='Raw data')
ax.plot(x4_fit, y4_fit, 'k', label='Fitted curve')
ax.plot(x4_fit, y4_B, 'k--', label='Fit 2nd Component')
ax.set_title('Fitting an power curve')
ax.set_ylabel('y')
ax.set_xlabel('X')
ax.legend()

y4 = np.array([146., 94., 79.,62., 41.])
x4 = np.array([375., 470., 528., 625.,880.])

from scipy.optimize import differential_evolution

def func(parameters, *data):

    #we have 3 parameters which will be passed as parameters and
    #"experimental" x,y which will be passed as data

    a1,a2,A1,A2 = parameters
    x,y = data

    result = 0

    for i in range(len(x)):
        result += (A1*x[i]**(-a1) + A2*x[i]**(-a2) - y[i])**2

    return result**0.5

bounds = [(0.8, 1.4), (1.5,15.0), (1000.0, 1.0e10), (1000.0, 1.0e30)]


#packing "experimental" data into args
args = (x4,y4)

result = differential_evolution(func, bounds, args=args)

result.x

def func1(parameters, *data):

    #we have 3 parameters which will be passed as parameters and
    #"experimental" x,y which will be passed as data

    a1,A1= parameters
    x,y = data

    result = 0

    for i in range(len(x)):
        result += (A1*x[i]**(-a1) -y4[i])**2

    return result

bounds = [(0.5, 10), (100.0, 1.0e10)]


#packing "experimental" data into args
args = (x4,y4)

result = differential_evolution(func1, bounds, args=args)

result.x



alpha1 = [optimizer(bounds)[0] for i in range(20)]
alpha2 = [optimizer(bounds)[1] for i in range(20)]

min(alpha2), max(alpha2)

x4_fit = np.linspace(np.min(x4), np.max(x4), 100)

alpha1 = result.x[0]
alpha2 = result.x[1]
A = result.x[2]
B = result.x[3]

y4_A = A*x4_fit**(-alpha1)
y4_B = B*x4_fit**(-alpha2)

ax = plt.axes()
# ax.scatter(x4, y4, label='Raw data')
ax.plot(x4_fit, y4_A, 'k', label='Fit curve 1st Component')
ax.plot(x4_fit, y4_B, 'k--', label='Fit curve 2nd Component')
ax.set_title('Fitting an power curve')
ax.set_ylabel('y')
ax.set_xlabel('X')
ax.legend()

y4_fit = result.x[1]*x4_fit**(-result.x[0])

ax = plt.axes()
ax.scatter(x4, y4, label='Raw data')
ax.plot(x4_fit, y4_fit, 'k', label='Fit curve 1st Component')
#ax.plot(x4_fit, y4_B, 'k--', label='Fit curve 2nd Component')
ax.set_title('Fitting an power curve')
ax.set_ylabel('y')
ax.set_xlabel('X')
ax.legend()