Using the new formulation of the Cluster class

jupyter_notebooks/gIV_sIV_Nat_Comm_2017.ipynb

@@ -582,16 +582,10 @@
 #fig.tight_layout()
 """fig.savefig(os.path.join(figure_folder, 'capacitance_map.png'), format='png', dpi=300)
 fig.savefig(os.path.join(figure_folder, 'capacitance_map.pdf'), format='pdf',bbox_inches = 'tight', pad_inches = 2.0)"""
 ```
 
-%% Cell type:code id: tags:
-
-``` python
-cap_vec
-```
-
 %% Cell type:markdown id: tags:
 
 ## Extracting Switching Current
 * The forward switching current is calculated as the difference between the forward and reverse currents for positive biases. 
 * The reverse switching current is calculated as the difference between the reverse and the forward currents for negative biases
@@ -952,17 +946,17 @@
 P is half the full charge (-Ps to +Ps = 2Ps)
 
 %% Cell type:code id: tags:
 
 ``` python
-cap_diameter = 500E-9 #Diameter of cap
+cap_diameter = 500E-9  # Diameter of cap in m
 
 cap_area = np.pi*(cap_diameter/2)**2
 P_forward = 0.5*(Q_integrated[:,:,0]/(cap_area))  # Polarization in C/m^2, 
-P_forward = P_forward*1E6*1E-4 #Polarization in uC (*1E6) /cm^2 (*1E-4)
+P_forward = P_forward*1E6*1E-4  # Polarization in uC (*1E6) /cm^2 (*1E-4)
 P_reverse = 0.5*(Q_integrated[:,:,1]/(cap_area)) #Polarization in C/m^2
-P_reverse = P_reverse*1E6*1E-4 #Polarization in uC (*1E6) /cm^2 (*1E-4)
+P_reverse = P_reverse*1E6*1E-4  # Polarization in uC (*1E6) /cm^2 (*1E-4)
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Now visualize
@@ -1094,49 +1088,49 @@
 Then fit to a Gaussian, with the constraint that the area equals that of the switched fraction
 
 %% Cell type:code id: tags:
 
 ``` python
-#Define some functions
+# Define some functions
 def gauss1(x,*p1):
     A1,x1,sigma1= p1
     return (A1*np.exp(-((x-x1)**2)/(2*sigma1*sigma1)))
 
 # here you include the penalization factor
 def residuals(p,x,y):
-    #Calcualte the area given the whole xvector
+    # Calcualte the area given the whole xvector
     fitted_gaussian = gauss1(xvec, *p) #Gaussian, for whole x
     integral = np.trapz(fitted_gaussian,xvec) #Area, for whole x
     penalization = abs(integral - area_of_Rayleigh)*1E4 #c.f. with area of Rayleigh
     return y - gauss1(x, *p) - penalization
 
 coef_disorder_fit = np.zeros(shape=(P_switching_fwd.shape[0],3))
 middle_ind_mat= np.zeros(shape=(P_switching_fwd.shape[0],1))
 
-#Fitting guesses
+# Fitting guesses
 p0 = [1, 4E-4, 1E-4] 
 lb = [0,0,0]
 ub=[1,1,1]
 bounds = (lb,ub)
 
 xvec = np.linspace(0,tstep*time_ratio_fwd, len(xvec_Pr_fwd[:])) #xvec
 
 for pix_ind, yvec in enumerate(P_switching_fwd):
 #If it's a switching event
     if pix_ind in switching_points[0]:
-        #Calculate the area
-        #del area_of_Rayleigh
+        # Calculate the area
+        # del area_of_Rayleigh
         global area_of_Rayleigh 
         area_of_Rayleigh = np.trapz(yvec, xvec)
         cum_area = scipy.integrate.cumtrapz(yvec,xvec)
 
-        #Take the point where the area is half of the total as the cutoff point
+        # Take the point where the area is half of the total as the cutoff point
         middle_index = np.argsort((cum_area - 0.5*area_of_Rayleigh)**2)[0]
-        #Smaller x,yvecs
+        # Smaller x,yvecs
         new_x, new_y = xvec[middle_index:],yvec[middle_index:]
         
-        #Apply the area constraint
+        # Apply the area constraint
         if len(new_x)>4:
             popt2, pcov2 = scipy.optimize.leastsq(func=residuals, x0=p0, args=(new_x,new_y))   
             y_fit2 = gauss1(xvec, *popt2)
             coef_disorder_fit[pix_ind,:] = popt2
             middle_ind_mat[pix_ind,:] = middle_index
@@ -1149,11 +1143,11 @@
 At a specific location
 
 %% Cell type:code id: tags:
 
 ``` python
-#row, col = 115, 51
+# row, col = 115, 51
 row, col = 110, 134
 pix_ind = row * num_cols + col
 
 xvec = np.linspace(0,tstep*time_ratio_fwd, len(xvec_Pr_fwd[:]))
 yvec = P_switching_fwd[pix_ind,:] #yvec
@@ -1336,19 +1330,20 @@
 
 %% Cell type:code id: tags:
 
 ``` python
 num_clusts = 6
+estimator = KMeans(n_clusters=num_clusts)
+
+proc = px.Cluster(h5_siv_raw, estimator)
 
-h5_siv_kmeans_grp = px.hdf_utils.check_for_old(h5_siv_raw, 'Cluster', 
-                                               {'n_clusters':num_clusts, 'n_jobs':1})
-if h5_siv_kmeans_grp is None:
+if proc.duplicate_h5_groups is None:
     print('Performing K-means now')
-    clusterer = px.processing.Cluster(h5_siv_raw, 'KMeans', n_clusters=num_clusts)
-    h5_siv_kmeans_grp = clusterer.do_cluster()
+    h5_siv_kmeans_grp = proc.compute()
 else:
     print('Taking previous results')
+    h5_siv_kmeans_grp = proc.duplicate_h5_groups[-1]
 
 # fig_km, axes_km = px.plot_utils.plot_cluster_h5_group(h5_siv_kmeans_grp, siv_y_label)
 
 fig, axes = plt.subplots(ncols=2, figsize=(11.5,5))
 _ = px.plot_utils.plot_map(axes[0], np.reshape(h5_siv_kmeans_grp['Labels'][()],(siv_num_rows, siv_num_cols)),