<jupyter_start><jupyter_text>Use py-earth to smoothen data and peakutils to find peaks<jupyter_code>matx_filename = 'threegaussian_broadpeak.txt'
datanm, datatime, dataz_matx = loaddata(matx_filename)
## get the peak position dataframe of true data set
true_df, smooth_df = earth_peak_matrix(datanm, dataz_matx, 0.1, 0, 10)
from syntheticdata_threegaussian_noshift_broadpeak import gaussian, monoexp, data_matrix
def make_data_matrix(nm_array, time_coeff_array, spectrum):
    """
    generates a two-way data matrix based on a known 
    spectrum and the spectrum's decay
    
    Args:
        nm_array: wavelength array
        time_coeff_array: an array that describes the decay
        spectrum: an array of the spectrum
    
    Returns:
        data_matrix: a matrix that contains the spectrum at each time
    """
    
    data_matrix = np.empty((np.shape(nm_array)[0], np.shape(time_coeff_array)[0]))
    for i in range(np.shape(time_coeff_array)[0]):
        data_matrix[:, i] = time_coeff_array[i] * spectrum
    
    return data_matrix
#create time array  
time = np.arange(0, 5000, 1)

#create wavelength array
nm = np.arange(850, 1600, 1)

a1 = 1
#center and FWHM of the gaussian 
x0_1 = 950
sigma_1 = 30

#life-time of the gaussian
tau1 = 10


#create a second gaussian
a2 = 0.3
x0_2 = 1300
sigma_2 = 100
tau2 = 5000


#create gaussian #3
a3 = 0.2
x0_3 = 1100
sigma_3 = 600
tau3 = 5000


#generate a gaussian model
species_1 = gaussian(nm, a1, x0_1, sigma_1)

#generate an array of time-coefficients 
#describing a mono-exponential decay with a given lifetime 
time_coeff_1 = monoexp(time, tau1)

#generate a data matrix that contains a gaussian model at each
#time and decays mono-exponentially
data_matrix_1 = make_data_matrix(nm, time_coeff_1, species_1)


#generate a second data matrix that contains a gaussian model
#at each time and decays mono-exponentially
species_2 = gaussian(nm, a2, x0_2, sigma_2)
time_coeff_2 = monoexp(time, tau2)
data_matrix_2 = make_data_matrix(nm, time_coeff_2, species_2)


#generate the third data matrix that contains a gaussian model
#at each time and decays mono-exponentially
species_3 = gaussian(nm, a3, x0_3, sigma_3)
time_coeff_3 = monoexp(time, tau3)
data_matrix_3 = make_data_matrix(nm, time_coeff_3, species_3)



#generate a two-gaussian mixture model by adding 
#the two gaussians above
data_matrix_sum = data_matrix_1 + data_matrix_2 + data_matrix_3
def plot_timeslice(nm, time, select_times, data):
    select_times_idx = [find_nearest(time, x) for x in select_times]
    num_times = len(select_times_idx)
    data_timeslice = np.empty((len(nm), num_times))
    
    plt.figure()
    plt.xlabel('Wavelength (nm)')
    plt.ylabel('Intensity (a.u.)')
    for i in range(num_times):
        data_timeslice_i = data[:, select_times_idx[i]]
        data_timeslice[:, i] = data_timeslice_i
        plt.plot(nm, data_timeslice_i, label = select_times[i])
    plt.legend()
    return data_timeslice

def find_nearest(array,value):
    idx = (np.abs(array-value)).argmin()
    return idx

def get_timeslice(nm, time, select_time, data):
    select_times_idx = find_nearest(time, select_time)
    data_timeslice = data[:, select_times_idx]
    data_timeslice_norm = data_timeslice / np.max(data_timeslice)
    
    return data_timeslice_norm<jupyter_output><empty_output><jupyter_text>Plot time-slices<jupyter_code>#select time points to plot
times = [1, 10, 100, 1000, 3000, 5000]

#plot 1d time-slice 
timeslice_matx = plot_timeslice(datanm, datatime, times, data_matrix_sum)<jupyter_output><empty_output><jupyter_text>Plot normalized individual spectrum<jupyter_code>time = 1000
timeslice_matx_1 = get_timeslice(datanm, datatime, time, data_matrix_1)
timeslice_matx_2 = get_timeslice(datanm, datatime, time, data_matrix_2)
timeslice_matx_3 = get_timeslice(datanm, datatime, time, data_matrix_3)
timeslice_matx_sum = get_timeslice(datanm, datatime, time, data_matrix_sum)


plt.figure()
plt.plot(datanm, timeslice_matx_1)
plt.plot(datanm, timeslice_matx_2)
plt.plot(datanm, timeslice_matx_3)
plt.plot(datanm, timeslice_matx_sum)
plt.legend(('feature 1', 'feature 2', 'feature 3', 'sum'))<jupyter_output><empty_output><jupyter_text>Peak-finding doesn't work in this scenario<jupyter_code>true_df
smooth_df<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup