Skip to content
Snippets Groups Projects
Commit e14a6fc0 authored by LeGuen Yann's avatar LeGuen Yann
Browse files

BaseCode

parent aa9259ff
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 781 additions and 0 deletions
StateDiagram.py 0 → 100644
+ 781
− 0
View file @ e14a6fc0
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
# from mayavi import mlab
import numpy as np
from contourpy import contour_generator
import itertools as itt
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 PIL import Image
import cv2
import tifffile as TIFF
import json
from tabulate import tabulate
scriptFolder = os.path.dirname(os.path.abspath(__file__))
#%% base variable
mainDir = scriptFolder
dirList = [
r"\24-11-15\1kHz\QWP45deg",
r"\24-11-16\1kHz\QWP0deg",
r"\24-11-16\1kHz\QWP45deg",
r"\24-11-16\1kHz\QWP315deg",
r"\24-11-17\1kHz\QWP0deg"
]
saveDir = mainDir + "/output/"
if not os.path.exists(saveDir):
os.makedirs(saveDir)
pxSizeX = 3.287E-6 # m/px bin by 4
pxSizeY = 1.156E-6 # m/px bin by 4
repRate = 1e3 #Hz
LOWSTATE = "P-"
HIGHSTATE = "P+"
TAG_INDEX = 699
#%% def funct
def getImageTagsDict(FilePath):
Tags = []
with TIFF.TiffFile(FilePath) as tif:
for page in tif.pages:
Tags.append({tag['Decription']:{k:tag[k] for k in ['Short','units','value']} for tag in json.loads(page.tags[TAG_INDEX].value)})
return Tags
def getImageTagsList(FilePath):
Tags = []
with TIFF.TiffFile(FilePath) as tif:
for page in tif.pages:
Tags.append(json.loads(page.tags[TAG_INDEX].value))
return Tags
def mapDictList(f,d):
return {k: [f(k,a) for a in v] for k, v in d.items()}
def cropIm(im,ROI):
return im[int(ROI[1]):int(ROI[1]+ROI[3]), int(ROI[0]):int(ROI[0]+ROI[2])]
def plotIm(im,fig,ax,title=None,aspectRatio=None,colorBar=None):
imx = ax.imshow(im,aspect=aspectRatio,cmap="bwr",vmin=0,vmax=1)
ax.set_title(title)
if colorBar==True:
fig.colorbar(imx)
def plotDictListIm(d,keyFilter=None,aspectRatio=None):
keyList = filter(keyFilter,d.keys())
for prop in keyList:
for im in d[prop]:
fig, ax = plt.subplots()
plotIm(im,fig,ax,title=prop,aspectRatio=aspectRatio,colorBar=True)
plt.show()
def plotDictListValue(d,figax=None,keyFilter=None,Xaxis=None,logScale=False,label=None,title=None):
keyList = list(filter(keyFilter,d.keys()))
if Xaxis!=None:
keyList = sorted(keyList,key=lambda x:getattr(x, Xaxis))
keyLookUp = []
X,Y = [],[]
for i,prop in enumerate(keyList):
for j,val in enumerate(d[prop]):
if Xaxis!=None:
X.append(getattr(prop, Xaxis))
else:
X.append(i)
Y.append(val)
keyLookUp.append((prop,j))
if figax==None:
figax = plt.subplots()
figax[1].plot(X,Y,'o',markersize=2,picker=5,label=label)
if logScale:
figax[1].set_yscale('log')
figax[1].set_title(title)
if label!=None:
figax[1].legend(bbox_to_anchor=(1.05, 1.0), loc='upper left')
figax[0].tight_layout()
return figax[0], figax[1], keyLookUp
def plotContourList(cs,fig,ax,title=None,color=None):
if color==None:
color="k"
for c in cs:
ax.plot(c[:,0],c[:,1],c=color)
ax.set_title(title)
def plotDictListContour(d,keyFilter=None,sortKey=None):
keyList = list(filter(keyFilter,d.keys()))
colorCycle = plt.rcParams['axes.prop_cycle']()
color = lambda x:next(colorCycle)["color"]
if sortKey!=None:
keyList = list(sorted(keyList,key=sortKey))
cmap = plt.cm.get_cmap('viridis')
cmin,cmax=min(map(sortKey,keyList)),max(map(sortKey,keyList))
color = lambda x:cmap((sortKey(x)-cmin)/(cmax-cmin))
fig, ax = plt.subplots()
for prop in keyList:
plotContourList(d[prop],fig,ax,color=color(prop))
ax.invert_yaxis()
ax.axis('scaled')
return fig, ax, keyList
def area(vs):
a = 0
x0,y0 = vs[0]
for [x1,y1] in vs[1:]:
dx = x1-x0
dy = y1-y0
a += 0.5*(y0*dx - x0*dy)
x0 = x1
y0 = y1
return abs(a)
#%% get data set propreties
communParam = None
for subDir in dirList:
imageDir = mainDir + subDir
imageFilenameList = list(filter(lambda x:x.endswith('.tiff'),os.listdir(imageDir)))
*_, communParam = itt.accumulate([getImageTagsList(os.path.join(imageDir,filename))[0] for filename in imageFilenameList],
lambda x,y:[a if a in y else {**a, "value":"\x1b[35m\x1b[1mvariable\x1b[0m"} for a in x],
initial=communParam)
print(tabulate([OrderedDict([(k,p[k]) for k in ['Decription','Short','value','units']]) for p in communParam],
colalign=('right',),
tablefmt="rst",
headers="keys"))
#%% open image and background
propDef = {
"Power" : ('Pump power',float),
"QWP" : ('Pump QWP angle',float),
"Compressor" : ('Compressor',int),
"Initial_State" : ('Initial state',str),
"Number_pulse" : ('Number of pulse',int),
# "Wavelength" : ('Pump wavelength',int)
}
Params = namedtuple('Params', list(propDef.keys()))
keyfunc = lambda x:[Params(**{k:v[1](getImageTagsDict(os.path.join(imageDir,x))[0][v[0]]["value"]) for k,v in propDef.items()}),"Bkg" in x]
openfunc = lambda x:np.array(Image.open(os.path.join(imageDir,x)))
rawImages = {}
bkgImage = {}
for subDir in dirList:
imageDir = mainDir + subDir
imageFilenameList = list(sorted(filter(lambda x:x.endswith('.tiff'),os.listdir(imageDir)),key=keyfunc))
for k, g in itt.groupby(imageFilenameList, keyfunc):
if k[1]:
bkgImage[k[0]] = list(map(openfunc,g))
else:
rawImages[k[0]] = list(map(openfunc,g))
#%% display images
# keyFilterFunc = lambda x:x.Initial_State=="P+" and x.Power==201.0
# plotDictListIm(rawImages,keyFilter=keyFilterFunc)
#%% fix background
bkgSumFunc = lambda prop,im : np.mean(im)
bkgSum = mapDictList(bkgSumFunc,bkgImage)
#%%
data = bkgSum
Xaxis = "Power"
groups = [
# "QWP",
"Initial_State",
# "Compressor",
"Number_pulse"
]
subFilterFunc = lambda x:True
keyfunc = lambda x:[getattr(x, key) for key in groups]
sortedKeys = list(sorted(data.keys(),key=keyfunc))
figax = plt.subplots(figsize=(10,7))
for k, g in itt.groupby(sortedKeys, keyfunc):
if subFilterFunc(k):
g = list(g)
plotDictListValue(data,
Xaxis=Xaxis,
keyFilter=lambda x:x in g and subFilterFunc(x),
figax=figax,
logScale=False,
label='\n'.join([str(a)+" = "+str(b) for a,b in zip(groups,k)]))
plt.show()
#%% norm images
normFunc = lambda prop,im : cv2.resize((im.astype('float32') - np.average(bkgImage[prop._replace(Initial_State=LOWSTATE)],axis=0).astype('float32'))\
/(np.average(bkgImage[prop._replace(Initial_State=HIGHSTATE)],axis=0).astype('float32') - np.average(bkgImage[prop._replace(Initial_State=LOWSTATE)],axis=0).astype('float32')),
dsize=(0,0),fx=pxSizeX/pxSizeY,fy=1)
normImages = mapDictList(normFunc,rawImages)
#%% norm images
# normFunc = lambda prop,im : (prop.Initial_State==HIGHSTATE) + (im.astype('float32') - np.average(bkgImage[prop],axis=0).astype('float32'))\
# /(np.average(bkgImage[prop._replace(Initial_State=HIGHSTATE)],axis=0).astype('float32') - np.average(bkgImage[prop._replace(Initial_State=LOWSTATE)],axis=0).astype('float32'))
# normImages = mapDictList(normFunc,rawImages)
#%% display images
keyFilterFunc = lambda x:x.Initial_State=="P+" and x.Compressor==170000 and x.QWP==315.00 and x.Number_pulse==100
plotDictListIm(normImages,keyFilter=keyFilterFunc)
#%% Special plot
keys = [Params(Power=130.0, QWP=315.0, Compressor=210000, Initial_State='P-', Number_pulse=10000),
Params(Power=130.0, QWP=315.0, Compressor=210000, Initial_State='P+', Number_pulse=10000),
Params(Power=130.0, QWP=0.0, Compressor=210000, Initial_State='P-', Number_pulse=10000),
Params(Power=130.0, QWP=0.0, Compressor=210000, Initial_State='P+', Number_pulse=10000),
Params(Power=130.0, QWP=45.0, Compressor=210000, Initial_State='P-', Number_pulse=10000),
Params(Power=130.0, QWP=45.0, Compressor=210000, Initial_State='P+', Number_pulse=10000)]
fig,axs = plt.subplots(2,3,figsize=(12,8))
for i,prop in enumerate(keys):
center = np.mean(np.where(np.logical_or(((prop.Initial_State == 'P-') and (normImages[prop][0]>=0.4)),
((prop.Initial_State == 'P+') and (normImages[prop][0]<=0.6)))),axis=1)
if np.any(np.isnan(center)):
center = np.array(normImages[prop][0].shape)/2
crop = np.s_[int(center[0]-150):int(center[0]+150),int(center[1]-150):int(center[1]+150)]
im = axs[i%2][i//2].imshow(normImages[prop][0][crop],cmap="bwr",vmin=0,vmax=1,interpolation='nearest', aspect='auto')
axs[i%2][i//2].tick_params(left = False,
right = False,
labelleft = False,
labelbottom = False,
bottom = False)
axs[0][0].set_ylabel('Init - (P-)', fontsize=20)
axs[1][0].set_ylabel('Init + (P+)', fontsize=20)
axs[0][0].set_title('σ- (QWP=315deg)', fontsize=20)
axs[0][1].set_title('L (QWP=0deg)', fontsize=20)
axs[0][2].set_title('σ+ (QWP=45deg)', fontsize=20)
cb_ax = fig.add_axes([1.03, 0.05, 0.06, 0.87])
cbar = fig.colorbar(im, cax=cb_ax)
cbar.ax.tick_params(labelsize=20)
plt.tight_layout()
plt.show()
#%% crop set
# cv2.namedWindow("output", cv2.WINDOW_NORMAL)
# ROIBox = cv2.selectROI("output",cv2.cvtColor(normImages[sorted(normImages.keys())[-1]][1],cv2.COLOR_GRAY2RGB))
# cv2.destroyAllWindows()
# cropImages = mapDictList(lambda prop,im:cropIm(im,ROIBox),normImages)
#%% display images
# keyFilterFunc = lambda x:x.QWP==45.0
# plotDictListIm(cropImages,keyFilter=keyFilterFunc)
#%% bin images
# del(binImages)
# threshold = 0.3
# binFunc = lambda prop,im : cv2.threshold(im,threshold,1,cv2.THRESH_BINARY)[1] if prop.Initial_State == LOWSTATE else\
# cv2.threshold(im,1-threshold,1,cv2.THRESH_BINARY_INV)[1]
# binImages = mapDictList(binFunc,normImages)
#%% sum value images
threshold = 0.5
imageSumFunc = lambda prop,im : np.sum(cv2.threshold(im,threshold,1,cv2.THRESH_BINARY)[1] if prop.Initial_State == LOWSTATE else\
cv2.threshold(im,1-threshold,1,cv2.THRESH_BINARY_INV)[1])
sumImages = mapDictList(imageSumFunc,normImages)
#%% sum plot
# subFilterFunc = lambda x:(x==[45.0, 'P-'] or x==[315.0, 'P+']) #and x.Power>52.5
# subFilterFunc = lambda x:x[2]==150000 #and x.Power>52.5
subFilterFunc = lambda x:True
X = np.linspace(min([k.Power for k in sumImages.keys()]),max([k.Power for k in sumImages.keys()]),1000)
# EmptySumThreshold = lambda x:max(1000*np.log((x-130)/10),10)
# OutSumThreshold = lambda x:max(10000*np.log((x-110)/10),10)
EmptySumThreshold = lambda x:10
groups = ["QWP","Initial_State","Compressor"]
keyfunc = lambda x:[getattr(x, key) for key in groups]
sortedKeys = list(sorted(sumImages.keys(),key=keyfunc))
figax = plt.subplots(figsize=(10,7))
for k, g in itt.groupby(sortedKeys, keyfunc):
if subFilterFunc(k):
g = list(g)
plotDictListValue(sumImages,
Xaxis="Power",
keyFilter=lambda x:x in g and subFilterFunc(x),
figax=figax,
logScale=True,
label='\n'.join([str(a)+" = "+str(b) for a,b in zip(groups,k)]))
# figax[1].plot(X,[EmptySumThreshold(x) for x in X],c="red")
# figax[1].plot(X,[OutSumThreshold(x) for x in X],c="red")
figax[1].axhline(y=EmptySumThreshold(k),c="red")
plt.show()
#%% interactive plot sum
# def on_close(event):
# fig.canvas.stop_event_loop()
# def on_pick(event,keylist,linesList,labelList,legend,d):
# line = event.artist
# if line in legend.get_lines():
# ind = legend.get_lines().index(line)
# origline = linesList[ind]
# visible = not origline.get_visible()
# origline.set_visible(visible)
# line.set_alpha(1.0 if visible else 0.2)
# fig.canvas.draw()
# else:
# ind2 = event.ind[0]
# ind1 = labelList.index(event.artist.get_label())
# k = keylist[ind1][ind2]
# # print(k)
# newfig, newax = plt.subplots()
# im = d[k[0]][k[1]]
# plotIm(im,newfig,newax,title=k)
# newfig.show()
# Xaxis = "Power"
# klList = []
# lineList = []
# labelList = []
# %matplotlib auto
# fig, ax = plt.subplots()
# for QWPpos in set(map(lambda x:x.QWP,sumImages.keys())):
# keyFilterFunc = lambda x:x.QWP==QWPpos
# fig, ax, kl = plotDictListValue(sumImages,
# keyFilter=keyFilterFunc,
# Xaxis="Power",
# figax=(fig, ax),
# logScale=True,
# label=QWPpos)
# klList.append(kl)
# labelList.append(str(QWPpos))
# # ax.axhline(y=EmptySumThreshold,c="red")
# ax.plot(X,[EmptySumThreshold(x) for x in X],c="red")
# ax.plot(X,[OutSumThreshold(x) for x in X],c="red")
# leg = ax.legend(bbox_to_anchor=(1.05, 1.0), loc='upper left')
# for legline in leg.get_lines():
# legline.set_picker(7)
# on_pick_kl = lambda x: on_pick(x, klList,ax.get_lines(),labelList,leg,normImages)
# cid = fig.canvas.mpl_connect('pick_event', on_pick_kl)
# fig.canvas.mpl_connect('close_event', on_close)
# fig.canvas.start_event_loop(timeout = -1)
# fig.canvas.mpl_disconnect(cid)
# %matplotlib inline
#%% filter out empty images
filteredImages = {}
for prop in sumImages.keys():
filterfunc = lambda x:np.sum(cv2.threshold(x,threshold,1,cv2.THRESH_BINARY)[1] if prop.Initial_State == LOWSTATE else\
cv2.threshold(x,1-threshold,1,cv2.THRESH_BINARY_INV)[1]) > EmptySumThreshold(prop.Power)
listIm = list(filter(filterfunc,normImages[prop]))
if len(listIm) > 0 :
filteredImages[prop] = listIm
#%% Separate Set
switchPair = [{"QWP":45.0,"Initial_State":"P-"},
{"QWP":315.0,"Initial_State":"P+"}]
noSwitchPair = [{"QWP":315.0,"Initial_State":"P-"},
{"QWP":45.0,"Initial_State":"P+"}]
SwitchImages = dict(filter(lambda x:any([subParam.items() <= x[0]._asdict().items() for subParam in switchPair]),
filteredImages.items()))
noSwitchImages = dict(filter(lambda x:any([subParam.items() <= x[0]._asdict().items() for subParam in noSwitchPair]),
filteredImages.items()))
linPolImages = dict(filter(lambda x:{"QWP":0.0}.items() <= x[0]._asdict().items() ,
filteredImages.items()))
#%%
contourFunc = lambda prop,im : sorted(contour_generator(z=im).lines(0.5), key=len, reverse=True)[0]
cont = mapDictList(contourFunc,SwitchImages)
#%%
def inscribeCircleFunc(prop,im):
binContour = cv2.fillPoly(np.zeros(im.shape, dtype=np.uint8),pts=[sorted(contour_generator(z=im).lines(0.5), key=len, reverse=True)[0].astype("int32")],color=255)
_, radius, _, center = cv2.minMaxLoc(cv2.distanceTransform(binContour,cv2.DIST_L2, cv2.DIST_MASK_PRECISE))
return radius
# cv2.drawContours(np.zeros(im.shape, dtype=np.uint8), sorted(contour_generator(z=im).lines(0.5), key=len, reverse=True),
# 0, 255, -1)
circleRadius = mapDictList(inscribeCircleFunc,SwitchImages)
#%% plot effective radius
subFilterFunc = lambda x:True
groups = ["Compressor","QWP","Initial_State","Number_pulse"]
keyfunc = lambda x:[getattr(x, key) for key in groups]
sortedKeys = list(sorted(circleRadius.keys(),key=keyfunc))
figax = plt.subplots(figsize=(10,7))
for k, g in itt.groupby(sortedKeys, keyfunc):
if subFilterFunc(k):
g = list(g)
plotDictListValue(circleRadius,
Xaxis="Power",
keyFilter=lambda x:x in g,
figax=figax,
logScale=True,
# label='\n'.join([str(a)+" = "+str(b) for a,b in zip(groups,k)])
)
plt.show()
#%% display contours
# keyFilterFunc = lambda x:200<=x.Power and x.Initial_State=="P-"
# plotDictListContour(cont,sortKey=lambda x:x.Power)
# plt.show()
#%% get effective radius
contourRadiusFunc = lambda prop,c:np.sqrt(area(c)*pxSizeY**2/np.pi)
contRadius = mapDictList(contourRadiusFunc,cont)
#%% plot effective radius
subFilterFunc = lambda x:True
groups = ["Compressor","QWP","Initial_State","Number_pulse"]
keyfunc = lambda x:[getattr(x, key) for key in groups]
sortedKeys = list(sorted(contRadius.keys(),key=keyfunc))
figax = plt.subplots(figsize=(10,7))
for k, g in itt.groupby(sortedKeys, keyfunc):
if subFilterFunc(k):
g = list(g)
plotDictListValue(contRadius,
Xaxis="Power",
keyFilter=lambda x:x in g,
figax=figax,
logScale=True,
# label='\n'.join([str(a)+" = "+str(b) for a,b in zip(groups,k)])
)
plt.show()
#%% export effective radius
data = np.array(sorted([[k.Power/repRate,k.Initial_State,r] for k,v in contRadius.items() for r in v]))
header = "Ep,InitialState,r"
np.savetxt(os.path.join(saveDir,'effectiveRadius'+dirList[0][1:]+'.csv'),data,header=header,fmt='%s',delimiter=",")
#%% def gauss radius
Rfunc = lambda Ep,s,Fth:s*np.sqrt(np.log(Ep/(np.pi*s**2*Fth)))
#%% fit radius
fitParam = {}
groups = ["Compressor","QWP","Initial_State","Number_pulse"]
subParams = namedtuple('subParams', groups)
keyfunc = lambda x:subParams(**{key:getattr(x, key) for key in groups})
sortedKeys = list(sorted(circleRadius.keys(),key=keyfunc))
for k, g in itt.groupby(sortedKeys, keyfunc):
try:
EpRList = np.array([[k.Power*1E-3/repRate,pxSizeY*r] for k in g for r in circleRadius[k]])
popt, pcov = curve_fit(Rfunc, EpRList[:,0], EpRList[:,1],p0=(1e-5,100),maxfev=1000000)
plt.plot(EpRList[:,0], EpRList[:,1],".")
newEp = np.linspace(min(EpRList[:,0]),max(EpRList[:,0]),500)
newR = [Rfunc(ep,*popt) for ep in newEp]
plt.plot(newEp,newR)
plt.title(k)
plt.show()
fitParam[k] = popt
except RuntimeError as err:
print(err)
plt.plot(EpRList[:,0], EpRList[:,1],".")
plt.title(k)
plt.show()
#%% def gauss fit
OneOverEpFunc = lambda r,s,Eth:(1/Eth)*np.exp(-(r/s)**2)
#%% fit radius
fitParam = {}
groups = ["Compressor","QWP","Initial_State","Number_pulse"]
subParams = namedtuple('subParams', groups)
keyfunc = lambda x:subParams(**{key:getattr(x, key) for key in groups})
sortedKeys = list(sorted(circleRadius.keys(),key=keyfunc))
for k, g in itt.groupby(sortedKeys, keyfunc):
try:
oneOverEpRList = np.array([[repRate/k.Power/1E-3,pxSizeY*r] for k in g for r in circleRadius[k]])
symDataList = np.vstack((np.flip(oneOverEpRList*[1,-1],axis=0),oneOverEpRList))
popt, pcov = curve_fit(OneOverEpFunc, symDataList[:,1], symDataList[:,0],p0=(1e-4,1E-6),maxfev=10000)
plt.plot(symDataList[:,1], symDataList[:,0],".")
newR = np.linspace(min(symDataList[:,1]),max(symDataList[:,1]),100)
newOneOverEp = [OneOverEpFunc(r,*popt) for r in newR]
plt.plot(newR,newOneOverEp)
plt.title(k)
plt.show()
fitParam[k] = popt
except RuntimeError as err:
print(err)
#%% fit radius
fitParam = {}
groups = ["Compressor","QWP","Initial_State","Number_pulse"]
subParams = namedtuple('subParams', groups)
keyfunc = lambda x:subParams(**{key:getattr(x, key) for key in groups})
sortedKeys = list(sorted(contRadius.keys(),key=keyfunc))
for k, g in itt.groupby(sortedKeys, keyfunc):
try:
EpRList = np.array([[k.Power*1E-3/repRate,r] for k in g for r in contRadius[k]])
popt, pcov = curve_fit(Rfunc, EpRList[:,0], EpRList[:,1],p0=(1e-4,50),maxfev=10000)
plt.plot(EpRList[:,0], EpRList[:,1],".")
newEp = np.linspace(min(EpRList[:,0]),max(EpRList[:,0]),100)
newR = [Rfunc(ep,*popt) for ep in newEp]
plt.plot(newEp,newR)
plt.title(k)
plt.show()
fitParam[k] = popt
except RuntimeError as err:
print(err)
#%% def gauss radius
Rfunc = lambda Ep,Fth:(1.18E-4)*np.sqrt(np.log(Ep/(np.pi*(1.18E-4)**2*Fth)))
#%% fit radius
fitParam = {}
groups = ["Compressor","QWP","Initial_State","Number_pulse"]
subParams = namedtuple('subParams', groups)
keyfunc = lambda x:subParams(**{key:getattr(x, key) for key in groups})
sortedKeys = list(sorted(contRadius.keys(),key=keyfunc))
for k, g in itt.groupby(sortedKeys, keyfunc):
try:
EpRList = np.array([[k.Power*1E-3/repRate,r] for k in g for r in contRadius[k]])
popt, pcov = curve_fit(Rfunc, EpRList[:,0], EpRList[:,1],p0=(50),maxfev=10000)
plt.plot(EpRList[:,0], EpRList[:,1],".")
newEp = np.linspace(min(EpRList[:,0]),max(EpRList[:,0]),100)
newR = [Rfunc(ep,*popt) for ep in newEp]
plt.plot(newEp,newR)
plt.title(k)
plt.show()
fitParam[k] = popt
except RuntimeError as err:
print(err)
#%% plot and export fit param
# data = fitParam
# print(f"Sigma\t=\t{fitParam[0]:.3e}\tm\nFth\t\t=\t{fitParam[1]:.3e}\tJ/m2")
# header = "sigma,Fth"
# np.savetxt(os.path.join(saveDir,'switchingThreshold'+dirList[0][1:]+'.csv'),data,header=header,fmt='%s',delimiter=",")
#%% plot and export fit param
data = np.array([list(k) + [f'{v[0]:.3E}',f'{v[1]:.3E}'] for k,v in fitParam.items()])
header = ','.join(next(iter(fitParam.keys()))._fields) + ',Sigma,Fth'
filePath = os.path.join(saveDir,'SwitchThreshold.csv')
with open(filePath, 'a+') as f:
f.write("\n")
np.savetxt(f,data,header=header,fmt='%s',delimiter=",")
#%%
thresholdMD = 0.4
validPointThreshold = 50
def EnclosingCircleFunc(prop,im):
contours = contour_generator(z=im).lines(thresholdMD if prop.Initial_State == LOWSTATE else 1-thresholdMD)
if len(contours) > 0 :
contoursPoint = np.vstack(contours).astype("int32")
validContourPoint = [p for p in contoursPoint if np.linalg.norm(p - np.mean(contoursPoint,axis=0)) <= validPointThreshold]
(x_axis,y_axis),radius = cv2.minEnclosingCircle(np.array(validContourPoint))
return radius
circleRadiusNoSwitch = mapDictList(EnclosingCircleFunc,noSwitchImages)
circleRadiusLinear = mapDictList(EnclosingCircleFunc,linPolImages)
#%% plot effective radius
subFilterFunc = lambda x:True
groups = ["Compressor","QWP","Initial_State","Number_pulse"]
keyfunc = lambda x:[getattr(x, key) for key in groups]
sortedKeys = list(sorted(circleRadiusNoSwitch.keys(),key=keyfunc))
figax = plt.subplots(figsize=(10,7))
for k, g in itt.groupby(sortedKeys, keyfunc):
if subFilterFunc(k):
# figax = plt.subplots(figsize=(10,7))
g = list(g)
plotDictListValue(circleRadiusNoSwitch,
Xaxis="Power",
keyFilter=lambda x:x in g,
figax=figax,
logScale=True,
# label='\n'.join([str(a)+" = "+str(b) for a,b in zip(groups,k)])
)
figax[1].axhline(y=validPointThreshold,c="red")
plt.show()
#%% plot effective radius
subFilterFunc = lambda x:True
groups = ["Compressor","QWP","Initial_State","Number_pulse"]
keyfunc = lambda x:[getattr(x, key) for key in groups]
sortedKeys = list(sorted(circleRadiusLinear.keys(),key=keyfunc))
figax = plt.subplots(figsize=(10,7))
for k, g in itt.groupby(sortedKeys, keyfunc):
if subFilterFunc(k):
# figax = plt.subplots(figsize=(10,7))
g = list(g)
plotDictListValue(circleRadiusLinear,
Xaxis="Power",
keyFilter=lambda x:x in g,
figax=figax,
logScale=True,
# label='\n'.join([str(a)+" = "+str(b) for a,b in zip(groups,k)])
)
figax[1].axhline(y=validPointThreshold,c="red")
plt.show()
#%% def gauss radius
Rfunc = lambda Ep,s,Fth:s*np.sqrt(np.log(Ep/(np.pi*s**2*Fth)))
#%% fit radius
fitParamLinear = {}
groups = ["Compressor","QWP","Initial_State","Number_pulse"]
subParams = namedtuple('subParams', groups)
keyfunc = lambda x:subParams(**{key:getattr(x, key) for key in groups})
sortedKeys = list(sorted(circleRadiusLinear.keys(),key=keyfunc))
for k, g in itt.groupby(sortedKeys, keyfunc):
try:
EpRList = np.array([[key.Power*1E-3/repRate,pxSizeY*r] for key in g for r in circleRadiusLinear[key]])
popt, pcov = curve_fit(Rfunc, EpRList[:,0], EpRList[:,1],p0=(1e-4,100),maxfev=1000000)
plt.plot(EpRList[:,0], EpRList[:,1],".")
newEp = np.linspace(min(EpRList[:,0]),max(EpRList[:,0]),100)
newR = [Rfunc(ep,*popt) for ep in newEp]
plt.plot(newEp,newR)
plt.title(k)
plt.show()
fitParamLinear[k] = popt
except RuntimeError as err:
print(err)
#%% plot and export fit param
data = np.array([list(k) + [f'{v[0]:.3E}',f'{v[1]:.3E}'] for k,v in fitParamLinear.items()])
header = ','.join(next(iter(fitParamLinear.keys()))._fields) + ',Sigma,Fmd'
filePath = os.path.join(saveDir,'LinearMDThreshold.csv')
with open(filePath, 'a+') as f:
f.write("\n")
np.savetxt(f,data,header=header,fmt='%s',delimiter=",")
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment