let's see

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)

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):

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')]