diff --git a/BGlib/be/analysis/utils/sidpy_sho_fitter.py b/BGlib/be/analysis/utils/sidpy_sho_fitter.py new file mode 100644 index 0000000..6893992 --- /dev/null +++ b/BGlib/be/analysis/utils/sidpy_sho_fitter.py @@ -0,0 +1,191 @@ +#Sidpy fitting +''' +This file will contain a class that is used for SHO fitting for sidpy datasets +for now I am just dumping the code from the notebook, but essentially this should be it's own class that is instantiated with the BEPS dataset +THen it should be able to do sho fitting and maybe loop fitting. + +''' + +# import numpy as np +import time +import h5py + +import pyNSID +import matplotlib.pyplot as plt +import numba + +import sidpy +#Let's open up a sample dataset and see... + +import SciFiReaders as sr +from scipy.optimize import curve_fit +import numpy as np + + + +def SHO_fit_flattened(wvec,*p): + Amp, w_0, Q, phi=p[0],p[1],p[2],p[3] + func = Amp * np.exp(1.j * phi) * w_0 ** 2 / (wvec ** 2 - 1j * wvec * w_0 / Q - w_0 ** 2) + return np.hstack([np.real(func),np.imag(func)]) + +def SHO_fit_abs(wvec,*p): + Amp, w_0, Q, phi=p[0],p[1],p[2],p[3] + func = Amp * np.exp(1.j * phi) * w_0 ** 2 / (wvec ** 2 - 1j * wvec * w_0 / Q - w_0 ** 2) + return np.abs(func) + +def my_guess_fn(freq_vec,ydata): + ydata = np.array(ydata) + amp_guess = np.abs(ydata)[np.argmax(np.abs(ydata))] + Q_guess = 50 + max_min_ratio = np.max(abs(ydata)) / np.min(abs(ydata)) + phi_guess = np.angle(ydata)[np.argmax(np.abs(ydata))] + w_guess = freq_vec[np.argmax(np.abs(ydata))] + + #Let's just run some Q values to find the closest one + Q_values = [5,10,20,50,100,200,500] + err_vals = [] + for q_val in Q_values: + p_test = [amp_guess/q_val, w_guess, q_val, phi_guess] + func_out = SHO_fit_flattened(freq_vec,*p_test) + complex_output = func_out[:len(func_out)//2] + 1j*func_out[(len(func_out)//2):] + amp_output = np.abs(complex_output) + err = np.mean((amp_output - np.abs(ydata))**2) + err_vals.append(err) + Q_guess = Q_values[np.argmin(err_vals)] + p0 = [amp_guess/Q_guess, w_guess, Q_guess, phi_guess] + return p0 + + +#Complex Gaussian Guess function +from numpy import exp, abs, sqrt, sum, real, imag, arctan2, append + +def SHOestimateGuess(w_vec, resp_vec, num_points=5): + """ + Generates good initial guesses for fitting + + Parameters + ------------ + w_vec : 1D numpy array or list + Vector of BE frequencies + resp_vec : 1D complex numpy array or list + BE response vector as a function of frequency + num_points : (Optional) unsigned int + Quality factor of the SHO peak + + Returns + --------- + retval : tuple + SHO fit parameters arranged as amplitude, frequency, quality factor, phase + """ + + ii = np.argsort(abs(resp_vec))[::-1] + + a_mat = np.array([]) + e_vec = np.array([]) + + for c1 in range(num_points): + for c2 in range(c1 + 1, num_points): + w1 = w_vec[ii[c1]] + w2 = w_vec[ii[c2]] + X1 = real(resp_vec[ii[c1]]) + X2 = real(resp_vec[ii[c2]]) + Y1 = imag(resp_vec[ii[c1]]) + Y2 = imag(resp_vec[ii[c2]]) + + denom = (w1 * (X1 ** 2 - X1 * X2 + Y1 * (Y1 - Y2)) + w2 * (-X1 * X2 + X2 ** 2 - Y1 * Y2 + Y2 ** 2)) + if denom > 0: + a = ((w1 ** 2 - w2 ** 2) * (w1 * X2 * (X1 ** 2 + Y1 ** 2) - w2 * X1 * (X2 ** 2 + Y2 ** 2))) / denom + b = ((w1 ** 2 - w2 ** 2) * (w1 * Y2 * (X1 ** 2 + Y1 ** 2) - w2 * Y1 * (X2 ** 2 + Y2 ** 2))) / denom + c = ((w1 ** 2 - w2 ** 2) * (X2 * Y1 - X1 * Y2)) / denom + d = (w1 ** 3 * (X1 ** 2 + Y1 ** 2) - + w1 ** 2 * w2 * (X1 * X2 + Y1 * Y2) - + w1 * w2 ** 2 * (X1 * X2 + Y1 * Y2) + + w2 ** 3 * (X2 ** 2 + Y2 ** 2)) / denom + + if d > 0: + a_mat = append(a_mat, [a, b, c, d]) + + A_fit = abs(a + 1j * b) / d + w0_fit = sqrt(d) + Q_fit = -sqrt(d) / c + phi_fit = arctan2(-b, -a) + + H_fit = A_fit * w0_fit ** 2 * exp(1j * phi_fit) / ( + w_vec ** 2 - 1j * w_vec * w0_fit / Q_fit - w0_fit ** 2) + + e_vec = append(e_vec, + sum((real(H_fit) - real(resp_vec)) ** 2) + + sum((imag(H_fit) - imag(resp_vec)) ** 2)) + if a_mat.size > 0: + a_mat = a_mat.reshape(-1, 4) + + weight_vec = (1 / e_vec) ** 4 + w_sum = sum(weight_vec) + + a_w = sum(weight_vec * a_mat[:, 0]) / w_sum + b_w = sum(weight_vec * a_mat[:, 1]) / w_sum + c_w = sum(weight_vec * a_mat[:, 2]) / w_sum + d_w = sum(weight_vec * a_mat[:, 3]) / w_sum + + A_fit = abs(a_w + 1j * b_w) / d_w + w0_fit = sqrt(d_w) + Q_fit = -sqrt(d_w) / c_w + phi_fit = np.arctan2(-b_w, -a_w) + + H_fit = A_fit * w0_fit ** 2 * exp(1j * phi_fit) / (w_vec ** 2 - 1j * w_vec * w0_fit / Q_fit - w0_fit ** 2) + + if np.std(abs(resp_vec)) / np.std(abs(resp_vec - H_fit)) < 1.2 or w0_fit < np.min(w_vec) or w0_fit > np.max( + w_vec): + p0 = SHOfastGuess(w_vec, resp_vec) + else: + p0 = np.array([A_fit, w0_fit, Q_fit, phi_fit]) + else: + p0 = SHOfastGuess(w_vec, resp_vec) + + return p0 + +def SHOfastGuess(w_vec, resp_vec, qual_factor=200): + """ + Default SHO guess from the maximum value of the response + + Parameters + ------------ + w_vec : 1D numpy array or list + Vector of BE frequencies + resp_vec : 1D complex numpy array or list + BE response vector as a function of frequency + qual_factor : float + Quality factor of the SHO peak + + Returns + ------- + retval : 1D numpy array + SHO fit parameters arranged as [amplitude, frequency, quality factor, phase] + """ + amp_vec = abs(resp_vec) + i_max = int(len(resp_vec) / 2) + return np.array([np.mean(amp_vec) / qual_factor, w_vec[i_max], qual_factor, np.angle(resp_vec[i_max])]) + + +#Now let's fit them all with sidpy +#Let's try sidpy fitter +#Instantiate the SidFitter class + +p0 = SHOestimateGuess(freq_vec, ydata) + +lb = [1E-6, freq_vec.min(), 50, -2*np.pi] +ub = [1E-3, freq_vec.max(), 500, 2*np.pi] + + +fitter = sidpy.proc.fitter.SidFitter(beps_small, SHO_fit_flattened,num_workers=1, + guess_fn = SHOestimateGuess,ind_dims=[0,1,3,4], + threads=1, return_cov=False, return_fit=False, return_std=False, + km_guess=True,num_fit_parms = 4, n_clus = 5) + +n_workers = 8 + #for n_workers in [2,4,8]: + +fitter = sidpy.proc.fitter.SidFitter(beps_small, SHO_fit_flattened,num_workers=n_workers, + guess_fn = SHOestimateGuess,ind_dims=[0,1,3,4], + threads=1, return_cov=False, return_fit=False, return_std=False, + km_guess=True,num_fit_parms = 4, n_clus = 4) diff --git a/BGlib/be/translators/labview_h5_patcher.py b/BGlib/be/translators/labview_h5_patcher.py index 6ccdf75..64664e9 100644 --- a/BGlib/be/translators/labview_h5_patcher.py +++ b/BGlib/be/translators/labview_h5_patcher.py @@ -25,6 +25,46 @@ if sys.version_info.major == 3: unicode = str +def fix_BE_problem(h5_main): + """ + This function will fix dodgy BEline files that are output from some versions of V3 + This bug causes incorrect sizes of the spectroscopic dimensions. + WARNING: This will irretrievably change the h5 file! + Input: + - h5_main: main dataset + Output: + None. The file is modified in-place. + """ + spec_inds = h5_main.parent['Spectroscopic_Indices'] + spec_vals = h5_main.parent['Spectroscopic_Values'] + pos_inds = h5_main.parent['Position_Indices'] + pos_vals = h5_main.parent['Position_Values'] + + freq_vec = np.unique(spec_vals) + + spec_inds_new = np.zeros((1, len(freq_vec))) + spec_inds_new[0,:] = np.arange(len(freq_vec)) + spec_vals_new = np.zeros_like(spec_inds_new) + spec_vals_new[0,:] = freq_vec + + del h5_main.parent['Spectroscopic_Indices'] + del h5_main.parent['Spectroscopic_Values'] + h5_file = h5_main.file + + h5_file.flush() + + h5_spec_inds = h5_main.parent.create_dataset(name='Spectroscopic_Indices', shape = (1, len(freq_vec)), data = spec_inds_new) + h5_spec_vals = h5_main.parent.create_dataset(name='Spectroscopic_Values', shape = (1, len(freq_vec)), data = spec_vals_new) + + h5_spec_inds.attrs['labels'] = ['Frequency'] + h5_spec_inds.attrs['units'] = ['Hz'] + + h5_spec_vals.attrs['labels'] = ['Frequency'] + h5_spec_vals.attrs['units'] = ['Hz'] + + h5_file.flush() + + return class LabViewH5Patcher(Translator): """ @@ -34,6 +74,11 @@ class LabViewH5Patcher(Translator): """ def __init__(self): + """ + Patches the hdf5 files from the LabView V3 data aquisition software to meet the + standards of the Pycroscopy data format. + Output: Sidpy.Translator object + """ super(LabViewH5Patcher, self).__init__() def _parse_file_path(self, input_path): diff --git a/BGlib/be/viz/be_viz_utils.py b/BGlib/be/viz/be_viz_utils.py index 8a8ceb1..8dcb035 100644 --- a/BGlib/be/viz/be_viz_utils.py +++ b/BGlib/be/viz/be_viz_utils.py @@ -130,7 +130,11 @@ def __plot_loops_maps(ac_vec, resp_mat, grp_name, win_title, spec_var_title, mea num_rows = int(np.floor((np.sqrt(h5_main.shape[0])))) num_cols = int(np.reshape(h5_main, [num_rows, -1, h5_main.shape[1]]).shape[1]) else: - num_rows, num_cols = h5_main.pos_dim_sizes + if len(h5_main.pos_dim_sizes)==1: + num_rows=h5_main.pos_dim_sizes[0] + num_cols = 1 + else: + num_rows, num_cols = h5_main.pos_dim_sizes try: h5_spec_vals = h5_file[get_attr(h5_main, 'Spectroscopic_Values')]