Update StateDiagramV2.py : less import, add domain detection and sigmoid Fit

StateDiagramV2.py

 import matplotlib.pyplot as plt
-import matplotlib.gridspec as gridspec
-
 
 import numpy as np
 
-from contourpy import contour_generator
 import itertools as itt
 import functools as ft
-from types import MethodType
 from collections import namedtuple,OrderedDict
 
 import os
-import sys
 
-from scipy.optimize import curve_fit,least_squares
-from scipy.interpolate import RegularGridInterpolator
-from sklearn.cluster import OPTICS, Birch
+from scipy.optimize import curve_fit
 from sklearn.neighbors import LocalOutlierFactor
 from scipy.ndimage import gaussian_filter
 
-from PIL import Image
 import cv2
 
 import tifffile as TIFF
 import json
 from tabulate import tabulate
 
-import time
 from tqdm import tqdm
 
 scriptFolder = os.path.dirname(os.path.abspath(__file__))
@@ -296,7 +287,7 @@ imageDictIn = normImages
 
 #%% display images
 
-# keyFilterFunc = lambda x:x.Initial_State=="P+" and x.QWP==45.0 and x.Number_pulse==1000 # Filter images to display based on their propreties 
+# keyFilterFunc = lambda x:x.Initial_State=="P+" and x.QWP==45.0 and x.Number_pulse==100 # Filter images to display based on their propreties 
 # sortKey = lambda x:x.Power # Sort displayed images based on their propreties 
 # plotDictListIm(normImages,keyFilter=keyFilterFunc,sortKey=sortKey)
 
@@ -467,36 +458,36 @@ circles = mapDictList(getCircles,imageDictIn)
 
 #%% plot effective radius
 
-# radiusSelect = OUTERCERCLE # select circle radius to plot can be INNERCERCLE or OUTERCERCLE
+radiusSelect = OUTERCERCLE # select circle radius to plot can be INNERCERCLE or OUTERCERCLE
 
-# groups = ["QWP","Number_pulse"] # group data based on the propreties
-# # groups = ["Compressor","QWP","Number_pulse"] # group data based on the propreties
+groups = ["QWP","Number_pulse"] # group data based on the propreties
+# groups = ["Compressor","QWP","Number_pulse"] # group data based on the propreties
 
-# subFilterFunc = lambda x:True # Filter data to display based on group parameter list
+subFilterFunc = lambda x:True # Filter data to display based on group parameter list
 
-# keyfunc = lambda x:[getattr(x, key) for key in groups] # get keys of the grouped data
-# sortedKeys = list(sorted(circles.keys(),key=keyfunc)) # sort dict keys !! necessary for itt.groupby
+keyfunc = lambda x:[getattr(x, key) for key in groups] # get keys of the grouped data
+sortedKeys = list(sorted(circles.keys(),key=keyfunc)) # sort dict keys !! necessary for itt.groupby
 
 
-# for k, g in itt.groupby(sortedKeys, keyfunc): # group data with commun keyfunc output (more info : https://docs.python.org/3/library/itertools.html#itertools.groupby)
-#     if subFilterFunc(k): # apply filter
-#         figax = plt.subplots(figsize=(10,7)) # create new plot
-#         g = list(g) # convert from grouper type to list (if not the generator will empty itself after the first passage)
-#         keyfunc2 = lambda x:x.Initial_State # key for the sub grouping 
-#         sortedKeys2 = list(sorted(g,key=keyfunc2)) # sort sub group key !! necessary for itt.groupby
-#         for k2, g2 in itt.groupby(sortedKeys2, keyfunc2): # sub-group grouped data with commun keyfunc2 output the data will be on the same plot but with different label
-#             g2 = list(g2) # convert from grouper type to list (if not the generator will empty itself after the first passage)
-#             # plot the sub grouped data point 
-#             plotDictListValue(mapDictList(lambda prop,x:x[radiusSelect], circles),  # create a sub dictionary of circle only selecting one of the circles
-#                               Xaxis="Power",                                        # name of the proprety to use as the X axis
-#                               keyFilter=lambda x:x in g2,                           # only select point of the sub group
-#                               figax=figax,                                          # figure and axis to use
-#                               logScale=True,                                        # logaritmic Y scaling
-#                               label=k2,                                             # use the sub-grouping key as label
-#                               title='\n'.join([str(a)+" = "+str(b) for a,b in zip(groups,k)]) # format title using the outer group keys
-#                               )
+for k, g in itt.groupby(sortedKeys, keyfunc): # group data with commun keyfunc output (more info : https://docs.python.org/3/library/itertools.html#itertools.groupby)
+    if subFilterFunc(k): # apply filter
+        figax = plt.subplots(figsize=(10,7)) # create new plot
+        g = list(g) # convert from grouper type to list (if not the generator will empty itself after the first passage)
+        keyfunc2 = lambda x:x.Initial_State # key for the sub grouping 
+        sortedKeys2 = list(sorted(g,key=keyfunc2)) # sort sub group key !! necessary for itt.groupby
+        for k2, g2 in itt.groupby(sortedKeys2, keyfunc2): # sub-group grouped data with commun keyfunc2 output the data will be on the same plot but with different label
+            g2 = list(g2) # convert from grouper type to list (if not the generator will empty itself after the first passage)
+            # plot the sub grouped data point 
+            plotDictListValue(mapDictList(lambda prop,x:x[radiusSelect], circles),  # create a sub dictionary of circle only selecting one of the circles
+                              Xaxis="Power",                                        # name of the proprety to use as the X axis
+                              keyFilter=lambda x:x in g2,                           # only select point of the sub group
+                              figax=figax,                                          # figure and axis to use
+                              logScale=True,                                        # logaritmic Y scaling
+                              label=k2,                                             # use the sub-grouping key as label
+                              title='\n'.join([str(a)+" = "+str(b) for a,b in zip(groups,k)]) # format title using the outer group keys
+                              )
 
-#         plt.show()
+        plt.show()
 
 #%% export radius
 
@@ -595,7 +586,7 @@ for k, g in itt.groupby(sortedKeys, keyfunc): # group data with commun keyfunc o
         plt.title(k) # use group key as title
         plt.show()
         
-        fitParam[k] = np.array([popt[0],popt[1]/np.pi/popt[0]**2]) # convert to sigma and fluence threshold and store the fit parameters into the fitParam dictonary
+        fitParam[k] = np.array([np.abs(popt[0]),popt[1]/np.pi/popt[0]**2]) # convert to sigma and fluence threshold and store the fit parameters into the fitParam dictonary
     except RuntimeError as err: # if fit fails
         print(err) # print error
         
@@ -613,3 +604,74 @@ paramName = ",".join(parameters) + ",Sigma,Fth"
 units = ",".join([next(filter(lambda x:x['Decription']==propDef[param][0],communParam))['units'] for param in parameters]) +\
     ',m,J/m2'
 np.savetxt(os.path.join(saveDir,'GaussFit.csv'),data,header=paramName+'\n'+units,fmt='%s',delimiter=",") # export to .csv as str array and with the constructed header
+
+
+#%% get sigmoid threshold
+
+fitTanhFunc = lambda x,offset,w:0.5*(1+np.tanh((x-offset)/w)) # function to fit
+
+domainDetectionCount = {}
+sigmoidFit = {}
+
+radiusSelect = OUTERCERCLE # select circle radius to plot can be INNERCERCLE or OUTERCERCLE
+
+groups = [param for param in propDef.keys() if param != "Power"] # group data based on the propreties except Power
+subParams = namedtuple('subParams', groups) # create a subparameter nametuple class to store selected propreties
+
+keyfunc = lambda x:subParams(**{key:getattr(x, key) for key in groups}) # get keys of the grouped data
+sortedKeys = list(sorted(circles.keys(),key=keyfunc)) # sort dict keys !! necessary for itt.groupby
+
+for k, g in itt.groupby(sortedKeys, keyfunc): # group data with commun keyfunc output (more info : https://docs.python.org/3/library/itertools.html#itertools.groupby)
+    g = list(sorted(g,key=lambda x:x.Power)) # convert and sort from grouper type to list (if not the generator will empty itself after the first passage)
+    for k2 in g: # for each key 
+        domainDetectionCount[k2] = np.sum(list(map(lambda x:x[radiusSelect]>0,circles[k2]))) # count the number of non zero selected radius
+    data = np.array([[k2.Power/1E3/repRate,domainDetectionCount[k2]/len(circles[k2])] for k2 in g]) # calculate data to fit Energy pulse [J] ; number of count
+    try:
+        popt, pcov = curve_fit(fitTanhFunc, data[:,0], data[:,1],p0=(np.min(data[:,0]),1),maxfev=1000000) # fit with the previously defined function fitTanhFunc
+        perr = np.sqrt(np.diag(pcov)) # calculate fit standard deviation
+        plt.plot(data[:,0], data[:,1],".") # plot to be fitted data
+        newEp = np.linspace(min(data[:,0]),max(data[:,0]),500) # X for ploting of the the fited function calculated from the min and max of the data
+        newN = [fitTanhFunc(p,*popt) for p in newEp] # Y for ploting of the the fited function calculaded at every X 
+        plt.plot(newEp,newN) # plot fited function
+        plt.title(k) # use group key as title
+        plt.show()
+        
+        sigmoidFit[k] = [*popt,*perr] # store fit parameters and standard deviation
+        
+    except RuntimeError as err: # if fit fails
+        print(err,k) # print error
+        plt.plot(data[:,0], data[:,1],".") # plot data without fit
+        plt.title(k) # use group key as title
+        plt.show()
+
+#%% export domain detection count
+
+sortKey = lambda x:[getattr(x[0], p) for p in x[0]._fields if p != "Power"] + [x[0].Power] # function to sort parameters with power as last sort element
+dataMapper = lambda x : list(x[0]) + [x[0].Power/1E3/repRate , x[1]] # function to extract parameters and data from the domainDetectionCount dictonary
+
+# map domainDetectionCount dictonary to export data
+data = list(map(dataMapper,sorted(domainDetectionCount.items(),key=sortKey)))
+
+# format the columns name
+paramName = ",".join(next(iter(domainDetectionCount.keys()))._fields) + ",Pulse Energy,Count"
+
+# get the units of the selected column parameters and append the fixed ouput unit for count
+units = ",".join([next(filter(lambda x:x['Decription']==propDef[param][0],communParam))['units'] for param in next(iter(domainDetectionCount.keys()))._fields]) + ',J,'
+
+np.savetxt(os.path.join(saveDir,'DetectionCount.csv'),data,header=paramName+'\n'+units,fmt='%s',delimiter=",") # export to .csv as str array and with the constructed header
+
+#%% export sigmoid fit parameters
+
+parameters = [param for param in propDef.keys() if param != "Power"] # list of propreties to use as column to discribe the output
+
+# retrive propreties to fill the columns and append the fit parameters
+data = np.array(sorted([[getattr(k,param) for param in parameters] + v for k,v in sorted(sigmoidFit.items(),key=lambda x:x[0])]))
+
+# format the columns name
+paramName = ",".join(parameters) +',EpOffset,EpOffsetStd,Width,WidthStd'
+
+# get the units of the selected column parameters and append the fixed ouput unit for power offset and width
+units = ",".join([next(filter(lambda x:x['Decription']==propDef[param][0],communParam))['units'] for param in parameters]) +\
+    ',J'*4
+    
+np.savetxt(os.path.join(saveDir,'sigmoidFit.csv'),data,header=paramName+'\n'+units,fmt='%s',delimiter=",") # export to .csv as str array and with the constructed header