added fitting function variation

Author: Jawad Chowdhury

BGlib/be/analysis/utils/be_loop.py

@@ -18,6 +18,9 @@
 from scipy.special import erf, erfinv
 import warnings
 from scipy.optimize import curve_fit
+from sklearn.metrics import r2_score
+from sklearn.gaussian_process import GaussianProcessRegressor
+from sklearn.gaussian_process.kernels import Matern
 
 # switching32 = np.dtype([('V+', np.float32),
 #                         ('V-', np.float32),
@@ -1012,7 +1015,6 @@ def intersection(L1, L2):
         except Exception as e:
             print('Error: ', e)                  
     ##    
-    from sklearn.metrics import r2_score
     best_guess = init_guess_coef_vec.copy()
     try:
         best_guess, _ = curve_fit(loop_fit_function, vdc, pr_vec, p0=init_guess_coef_vec, maxfev=5000)
@@ -1039,6 +1041,228 @@ def intersection(L1, L2):
     ##
     return best_guess        
 
+
+def generate_deepGP_guess(vdc, pr_vec, show_plots=False):
+    """
+    Given a single unfolded loop and centroid, return the best-fit parameter guess.
+    We start with an initial estimate based on the loop centroid and intersection points.
+    Then, we refine it by running the fitting program multiple times with randomized perturbations
+    around the initial guess and keeping the parameters that yield the lowest fitting error.
+
+    Parameters
+    -----------
+    vdc : 1D numpy array
+        DC offsets
+    pr_vec : 1D numpy array
+        Piezoresponse or unfolded loop
+    show_plots : Boolean (Optional. Default = False)
+        Whether or not the plot the convex hull, centroid, intersection points
+
+    Returns
+    -----------------
+    init_guess_coef_vec : 1D Numpy array
+        Fit guess coefficient vector
+    """
+
+    points = np.transpose(np.array([np.squeeze(vdc), pr_vec]))  # [points,axis]
+
+    geom_centroid, geom_area = calculate_loop_centroid(points[:, 0], points[:, 1])
+
+    hull = ConvexHull(points)
+
+    """
+    Now we need to find the intersection points on the N,S,E,W
+    the simplex of the complex hull is essentially a set of line equations.
+    We need to find the two lines (top and bottom) or (left and right) that
+    interect with the vertical / horizontal lines passing through the geometric centroid
+    """
+
+    def find_intersection(A, B, C, D):
+        """
+        Finds the coordinates where two line segments intersect
+
+        Parameters
+        ------------
+        A, B, C, D : Tuple or 1D list or 1D numpy array
+            (x,y) coordinates of the points that define the two line segments AB and CD
+
+        Returns
+        ----------
+        obj : None or tuple
+            None if not intersecting. (x,y) coordinates of intersection
+        """
+
+        def ccw(A, B, C):
+            """Credit - StackOverflow"""
+            return (C[1] - A[1]) * (B[0] - A[0]) > (B[1] - A[1]) * (C[0] - A[0])
+
+        def line(p1, p2):
+            """Credit - StackOverflow"""
+            A = (p1[1] - p2[1])
+            B = (p2[0] - p1[0])
+            C = (p1[0] * p2[1] - p2[0] * p1[1])
+            return A, B, -C
+
+        def intersection(L1, L2):
+            """
+            Finds the intersection of two lines (NOT line segments).
+            Credit - StackOverflow
+            """
+            D = L1[0] * L2[1] - L1[1] * L2[0]
+            Dx = L1[2] * L2[1] - L1[1] * L2[2]
+            Dy = L1[0] * L2[2] - L1[2] * L2[0]
+            if D != 0:
+                x = Dx / D
+                y = Dy / D
+                return x, y
+            else:
+                return None
+
+        if ((ccw(A, C, D) is not ccw(B, C, D)) and (ccw(A, B, C) is not ccw(A, B, D))) is False:
+            return None
+        else:
+            return intersection(line(A, B), line(C, D))
+
+    # start and end coordinates of each line segment defining the convex hull
+    outline_1 = np.zeros((hull.simplices.shape[0], 2), dtype=float)
+    outline_2 = np.zeros((hull.simplices.shape[0], 2), dtype=float)
+    for index, pair in enumerate(hull.simplices):
+        outline_1[index, :] = points[pair[0]]
+        outline_2[index, :] = points[pair[1]]
+
+    """Find the coordinates of the points where the vertical line through the
+    centroid intersects with the convex hull"""
+    y_intersections = []
+    for pair in range(outline_1.shape[0]):
+        x_pt = find_intersection(outline_1[pair], outline_2[pair],
+                                 [geom_centroid[0], hull.min_bound[1]],
+                                 [geom_centroid[0], hull.max_bound[1]])
+        if x_pt is not None:
+            y_intersections.append(x_pt)
+
+    '''
+    Find the coordinates of the points where the horizontal line through the
+    centroid intersects with the convex hull
+    '''
+    x_intersections = []
+    for pair in range(outline_1.shape[0]):
+        x_pt = find_intersection(outline_1[pair], outline_2[pair],
+                                 [hull.min_bound[0], geom_centroid[1]],
+                                 [hull.max_bound[0], geom_centroid[1]])
+        if x_pt is not None:
+            x_intersections.append(x_pt)
+
+    '''
+    Default values if not intersections can be found.
+    '''
+    if len(y_intersections) < 2:
+        min_y_intercept = min(pr_vec)
+        max_y_intercept = max(pr_vec)
+    else:
+        min_y_intercept = min(y_intersections[0][1], y_intersections[1][1])
+        max_y_intercept = max(y_intersections[0][1], y_intersections[1][1])
+
+    if len(x_intersections) < 2:
+        min_x_intercept = min(vdc) / 2.0
+        max_x_intercept = max(vdc) / 2.0
+    else:
+        min_x_intercept = min(x_intersections[0][0], x_intersections[1][0])
+        max_x_intercept = max(x_intersections[0][0], x_intersections[1][0])
+
+    # Only the first four parameters use the information from the intercepts
+    # a3, a4 are swapped in Stephen's figure. That was causing the branches to swap during fitting
+    # the a3, a4 are fixed now below:
+    init_guess_coef_vec = np.zeros(shape=9)
+    init_guess_coef_vec[0] = min_y_intercept
+    init_guess_coef_vec[1] = max_y_intercept - min_y_intercept
+    init_guess_coef_vec[2] = min_x_intercept
+    init_guess_coef_vec[3] = max_x_intercept
+    init_guess_coef_vec[4] = 0
+    init_guess_coef_vec[5] = 2  # 0.5
+    init_guess_coef_vec[6] = 2  # 0.2
+    init_guess_coef_vec[7] = 2  # 1.0
+    init_guess_coef_vec[8] = 2  # 0.2
+
+    if show_plots:
+        try:
+            fig, ax = plt.subplots()
+            ax.plot(points[:, 0], points[:, 1], 'o')
+            ax.plot(geom_centroid[0], geom_centroid[1], 'r*')
+            ax.plot([geom_centroid[0], geom_centroid[0]], [hull.max_bound[1], hull.min_bound[1]], 'g')
+            ax.plot([hull.min_bound[0], hull.max_bound[0]], [geom_centroid[1], geom_centroid[1]], 'g')
+            for simplex in hull.simplices:
+                ax.plot(points[simplex, 0], points[simplex, 1], 'k')
+            ax.plot(x_intersections[0][0], x_intersections[0][1], 'r*')
+            ax.plot(x_intersections[1][0], x_intersections[1][1], 'r*')
+            ax.plot(y_intersections[0][0], y_intersections[0][1], 'r*')
+            ax.plot(y_intersections[1][0], y_intersections[1][1], 'r*')
+            ax.plot(vdc, loop_fit_function(vdc, *init_guess_coef_vec))
+        except Exception as e:
+            print('Error: ', e)                  
+    ##
+    X = []
+    y = []    
+    #
+    best_guess = init_guess_coef_vec.copy()
+    best_err = np.inf
+    best_r2 = -np.inf
+    def objective(params):
+        """Returns error for a given parameter vector."""
+        try:
+            guess, _ = curve_fit(loop_fit_function, vdc, pr_vec, p0=params, maxfev=5000)
+            pred = loop_fit_function(vdc, *guess)
+            err = np.sum((pr_vec - pred) ** 2)
+            r2 = r2_score(pr_vec, pred)
+            return err, r2, guess
+        except RuntimeError:
+            return np.inf, -np.inf, params
+    #
+    for _ in range(10):
+        p0_rand = init_guess_coef_vec * (1 + 0.3 * np.random.randn(len(init_guess_coef_vec)))
+        err, r2, guess = objective(p0_rand)
+        if np.isfinite(err):        
+            X.append(p0_rand)
+            y.append(err)
+        if err < best_err:
+            best_guess, best_err, best_r2 = guess, err, r2
+        if best_r2 >= 0.95:
+            break
+    #       
+    kernel = Matern(length_scale=1.0, nu=2.5)
+    gp = GaussianProcessRegressor(kernel=kernel, alpha=1e-6, normalize_y=True)
+    #
+    for trial_no in range(50):
+        if best_r2 >= 0.95:
+            break
+        X_arr = np.array(X)
+        y_arr = np.array(y)
+        y_arr = np.nan_to_num(y_arr, nan=1e6, posinf=1e6, neginf=1e6)
+        if len(y_arr) == 0:
+            X_arr = np.array([init_guess_coef_vec])
+            y_arr = np.array([1e6])        
+            
+        gp.fit(X_arr, y_arr)
+        scale = np.maximum(np.abs(best_guess), 1e-3)
+        candidates = best_guess * (1 + 0.3 * np.random.randn(100, len(best_guess)) / scale)
+
+        mu, sigma = gp.predict(candidates, return_std=True)
+        
+        # Choose the most promising candidate (Expected Improvement heuristic)
+        acquisition = mu - 1.0 * sigma  # exploit low mean, explore high uncertainty
+        next_p = candidates[np.argmin(acquisition)]
+    
+        err, r2, guess = objective(next_p)
+        if np.isfinite(err):
+            X.append(next_p)
+            y.append(err)
+
+        if err < best_err:
+            best_guess, best_err, best_r2 = guess, err, r2    
+    ##
+    return best_guess      
+
+        
+
 ###############################################################################
 
 

BGlib/misc/bg_gui_basic.py

@@ -30,7 +30,7 @@
 from BGlib.be.analysis.utils.sidpy_sho_fitter import SHOestimateGuess, SHOestimateGuess, SHO_fit_flattened
 from PyQt5.QtGui import QTextCursor
 from PyQt5.QtCore import QObject, pyqtSignal
-from BGlib.be.analysis.utils.be_loop import projectLoop, loop_fit_function, generate_guess, generate_shallow_guess, generate_deep_guess, calc_switching_coef_vec
+from BGlib.be.analysis.utils.be_loop import projectLoop, loop_fit_function, generate_guess, generate_shallow_guess, generate_deep_guess, generate_deepGP_guess, calc_switching_coef_vec
 
 
 class EmittingStream(QObject):
@@ -434,7 +434,8 @@ def on_do_loop_guess(self):
             d_guess_fn = {
                 'Basic': generate_guess,
                 'Shallow': generate_shallow_guess,
-                'Deep': generate_deep_guess
+                'Deep': generate_deep_guess,
+                'DeepGP': generate_deepGP_guess
             }
             ##
             self.do_loop_fit_button.setEnabled(False)