import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
from scipy.integrate import trapz
from scipy import optimize
import warnings

def one_peak_gaussian_fit(x,y):
    y, denominator = __norm_gaussian(x,y)
    popt_gauss = __one_peak_gaussian_fit(x, y)
    fit_data = __one_peak_gaussian(x, *popt_gauss)
    return fit_data*denominator, popt_gauss

def __one_peak_gaussian_fit(x_ori,y_ori):
    """Do not use pico & Nano scale
        fit_data = __one_peak_gaussian(x_axis,popt_gauss[0],popt_gauss[1],popt_gauss[2])
        fit_data = __one_peak_gaussian(x_axis,*popt_gauss)
    """
    x = x_ori.copy()
    y = y_ori.copy() 
    x, y = __interpolate(x, y)
    x,y = __clip_negative(x, y)

    amp = y.max()
    cen = x[np.where(y == amp)][0]
    sigma = np.std(y)
    warnings.filterwarnings("error", category=UserWarning)
    try:
        popt_gauss, _ = optimize.curve_fit(__one_peak_gaussian, x, y, p0=[amp, cen, sigma])
    except UserWarning as e:
        popt_gauss = (amp,cen,sigma)
    except RuntimeError:
        popt_gauss = (amp,cen,sigma)
    return popt_gauss

def __norm_gaussian(x,y):
    # integral = 1
    y = np.abs(y)
    eps = 7/3. - 4/3. -1 
    denominator = np.abs(trapz(x,y)) + eps
    return y / denominator, denominator

def __one_peak_gaussian(x, amp,cen,sigma):
    eps = 7/.3 - 4/.3 - 1
    return amp *( 1 / (eps + sigma * (np.sqrt(2 * np.pi)))) * (np.exp((-1.0 / 2.0) * (((x-cen) / (sigma + eps))**2)))

def __clip_negative(x,y):
    clip_index = np.where(y > 0)
    return x[clip_index], y[clip_index]

def __interpolate(x,y, num:int=10000):
    fun = interp1d(x,y,kind='cubic')
    x_new = np.linspace(x.min(), x.max(),num=num,endpoint=True)
    y_new = fun(x_new)
    return x_new, y_new

if __name__ == "__main__":
    # Create synthetic data
    x = np.linspace(-200, 200, 500)
    A_true = 1.0
    mu_true = 0.0
    sigma_true = 2.0
    y = __one_peak_gaussian(x, A_true, mu_true, sigma_true)
    noise = 0.001 * np.random.normal(size=x.size)
    y_noisy = y + noise
    
    fit_data, popt = one_peak_gaussian_fit(x, y_noisy)
    
    plt.figure(figsize=(8,6))
    plt.plot(x, y_noisy, 'b.', label='Noisy Data')
    plt.plot(x, fit_data, 'r-', label='Fitted Gaussian')
    plt.title('One Peak Gaussian Fit')
    plt.legend()
    plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
from scipy.integrate import trapz
from scipy import optimize
import warnings

def one_peak_gaussian_fit(x,y):
    y, denominator = __norm_gaussian(x,y)
    popt_gauss = __one_peak_gaussian_fit(x, y)
    fit_data = __one_peak_gaussian(x, *popt_gauss)
    return fit_data*denominator, popt_gauss

def __one_peak_gaussian_fit(x_ori,y_ori):
    """Do not use pico & Nano scale
        fit_data = __one_peak_gaussian(x_axis,popt_gauss[0],popt_gauss[1],popt_gauss[2])
        fit_data = __one_peak_gaussian(x_axis,*popt_gauss)
    """
    x = x_ori.copy()
    y = y_ori.copy() 
    x, y = __interpolate(x, y)
    x,y = __clip_negative(x, y)

    amp = y.max()
    cen = x[np.where(y == amp)][0]
    sigma = np.std(y)
    warnings.filterwarnings("error", category=UserWarning)
    try:
        popt_gauss, _ = optimize.curve_fit(__one_peak_gaussian, x, y, p0=[amp, cen, sigma])
    except UserWarning as e:
        popt_gauss = (amp,cen,sigma)
    except RuntimeError:
        popt_gauss = (amp,cen,sigma)
    return popt_gauss

def __norm_gaussian(x,y):
    # integral = 1
    y = np.abs(y)
    eps = 7/3. - 4/3. -1 
    denominator = np.abs(trapz(x,y)) + eps
    return y / denominator, denominator

def __one_peak_gaussian(x, amp,cen,sigma):
    eps = 7/.3 - 4/.3 - 1
    return amp *( 1 / (eps + sigma * (np.sqrt(2 * np.pi)))) * (np.exp((-1.0 / 2.0) * (((x-cen) / (sigma + eps))**2)))

def __clip_negative(x,y):
    clip_index = np.where(y > 0)
    return x[clip_index], y[clip_index]

def __interpolate(x,y, num:int=10000):
    fun = interp1d(x,y,kind='cubic')
    x_new = np.linspace(x.min(), x.max(),num=num,endpoint=True)
    y_new = fun(x_new)
    return x_new, y_new

if __name__ == "__main__":
    # Create synthetic data
    x = np.linspace(-200, 200, 500)
    A_true = 1.0
    mu_true = 0.0
    sigma_true = 2.0
    y = __one_peak_gaussian(x, A_true, mu_true, sigma_true)
    noise = 0.001 * np.random.normal(size=x.size)
    y_noisy = y + noise
    
    fit_data, popt = one_peak_gaussian_fit(x, y_noisy)
    
    plt.figure(figsize=(8,6))
    plt.plot(x, y_noisy, 'b.', label='Noisy Data')
    plt.plot(x, fit_data, 'r-', label='Fitted Gaussian')
    plt.title('One Peak Gaussian Fit')
    plt.legend()
    plt.show()