From a2c508fa5d9a70b7f64e5a2f590f1e4810fe53c8 Mon Sep 17 00:00:00 2001 From: Marco Cammarata Date: Fri, 25 Nov 2016 11:03:07 +0100 Subject: [PATCH] various changes: some naming, added function for averaging and changed doShots to accept list of calibs --- xanes_analyzeRun.py | 157 +++++++++++++++++++++++++++++--------------- 1 file changed, 103 insertions(+), 54 deletions(-) diff --git a/xanes_analyzeRun.py b/xanes_analyzeRun.py index 2b2b608..2a12552 100644 --- a/xanes_analyzeRun.py +++ b/xanes_analyzeRun.py @@ -101,6 +101,54 @@ def showShots(im1,im2): a[0].plot(p1) a[1].plot(p2) +def ratioOfAverage(p1,p2,threshold=0.03): + """ + p1 and p2 are the energy spectrum. if 2D the first index has to be the shot number + calculate median ratio taking into account only regions where p1 and p2 are > 5% of the max """ + # check if they are 2D + if p1.ndim == 1: + p1 = p1[np.newaxis,:] + p2 = p2[np.newaxis,:] + # w1 and w2 are the weights + w1 = p1.copy(); w2 = p2.copy() + if threshold is not None: + # weights will be set to zero if intensity is smaller than 5% of max + # for each shots, get maximum + m1 = np.nanmax(p1,axis=1); m2 = np.nanmax(p2,axis=1) + # find where each spectrum is smaller than threshold*max_for_that_shot; they will be masked out + idx1 = p1 < (m1[:,np.newaxis]*threshold) + idx2 = p2 < (m2[:,np.newaxis]*threshold) + w1[idx1]=0 + w2[idx2]=0 + # using masked array because some pixel will have zero shots contributing + av1 = np.ma.average(p1,axis=0,weights=w1) + av1[av1.mask] = np.nan + av2 = np.ma.average(p2,axis=0,weights=w2) + av2[av2.mask] = np.nan + return av2/av1 + +def medianRatio(p1,p2,threshold=0.03): + """ + p1 and p2 are the energy spectrum. if 2D the first index has to be the shot number + calculate median ratio taking into account only regions where p1 and p2 are > 5% of the max """ + # check if they are 2D + if p1.ndim == 1: + p1 = p1[np.newaxis,:] + p2 = p2[np.newaxis,:] + p1 = np.ma.asarray( p1.copy() ) + p2 = np.ma.asarray( p2.copy() ) + if threshold is not None: + m1 = np.nanmax(p1,axis=1); m2 = np.nanmax(p2,axis=1) + # find where each spectrum is smaller than threshold*max_for_that_shot; they will be masked out + idx1 = p1 < (m1[:,np.newaxis]*threshold) + idx2 = p2 < (m2[:,np.newaxis]*threshold) + idx = idx1 & idx2 + p1.mask = idx + p2.mask = idx + ratio = p2/p1 + return np.ma.average(ratio,axis=0,weights=p1) + + class AnalyzeRun(object): def __init__(self,run,initAlign="auto",swapx=g_swapx,swapy=g_swapy): @@ -125,7 +173,8 @@ class AnalyzeRun(object): d = self.data self.spec1 = d.spec1 ; # spec1 is the one that is moved - self.spec2 = d.spec2 ; + self.spec2 = d.spec2 ; + self.E = alignment.defaultE try: self.loadTransform(initAlign) @@ -136,10 +185,15 @@ class AnalyzeRun(object): else: self.initAlign = initAlign - def getShot(self,shot=0,calib=None,bkgSub="line",roi=g_roi_height): + def getShots(self,shots=0,calib=None,bkgSub="line",roi=g_roi_height): + if shots == "all": + if calib != None: + shots = slice(0,self.nShotsPerCalib[calib]) + else: + shots = slice(0,self.data.spec1.nShots) # read data - im1 = self.spec1.getShots(shot,calib=calib) - im2 = self.spec2.getShots(shot,calib=calib) + im1 = self.spec1.getShots(shots,calib=calib) + im2 = self.spec2.getShots(shots,calib=calib) # subtractBkg bkg im1 = alignment.subtractBkg(im1,bkg_type=bkgSub) im2 = alignment.subtractBkg(im2,bkg_type=bkgSub) @@ -174,31 +228,40 @@ class AnalyzeRun(object): gui.save(fname) - def analyzeScan(self,initpars=None,nShotsPerCalib=20,nC=None,doFit=False,fitEveryCalib=False,nSaveImg=5): - """ this is a comment """ + def analyzeScan(self,initpars=None,shots=slice(0,30),calibs="all",nImagesToFit=0,nSaveImg=5): + """ nImagesToFit: number of images to Fit per calibcycle, (int or "all") """ if initpars is None: initpars= self.initAlign - if nC is None: nC = self.nCalib - out = [] - for i in range(nC): - if nShotsPerCalib == 'all': - shots = slice(self.nShotsPerCalib[i]) + if calibs == "all": calibs=list(range(self.nCalib)) + if isinstance(calibs,slice): calibs=list(range(self.nCalib))[calibs] + nC = len(calibs) + for ic,calib in enumerate(calibs): + if nImagesToFit == "all": + nToFit = self.nShotsPerCalib[calib] else: - shots = slice(nShotsPerCalib) - s1,s2 = self.getShot(shots,calib=i) - if fitEveryCalib is not False: - ret,bestTransf = alignment.doShots(s1[:fitEveryCalib],s2[:fitEveryCalib],doFit=True,\ + nToFit = nImagesToFit + #print("Memory available 1",x3py.toolsOS.memAvailable()) + s1,s2 = self.getShots(shots,calib=calib) + #print("Memory available 2",x3py.toolsOS.memAvailable()) + if nImagesToFit > 0: + ret,bestTransf = alignment.doShots(s1[:nToFit],s2[:nToFit],doFit=True,\ initpars=initpars,nSaveImg=nSaveImg,returnBestTransform=True); initpars = bestTransf; self.initAlign=bestTransf - ret = alignment.doShots(s1,s2,initpars=initpars,doFit=doFit,nSaveImg=nSaveImg) - self.results[i] = ret - print("Calib cycle %d -> %.3f (best FOM: %.2f)" % (i,self.scanpos[i],np.min(ret.fom))) - out.append(ret) - return out + if nToFit < s1.shape[0]: + ret2 = alignment.doShots(s1[nToFit:],s2[nToFit:],initpars=initpars,doFit=False,nSaveImg=0) + if ret is None: + ret = ret2 + else: + ret = alignment.unravel_results( (ret,ret2) ) + #print("Memory available 3",x3py.toolsOS.memAvailable()) + self.results[calib] = ret + #print("Memory available 4",x3py.toolsOS.memAvailable()) + print("Calib cycle %d/%d -> %.3f (best FOM: %.2f)" % (ic,nC,self.scanpos[ic],np.nanmin(ret.fom))) + return [self.results[c] for c in calibs] def doShot(self,shot=0,calib=None,initpars=None,im1=None,im2=None,doFit=True,show=False,showInit=False,save=False,savePlot="auto"): if initpars is None: initpars= self.initAlign if (im1 is None) or (im2 is None): - im1,im2 = self.getShot(shot,calib=calib); im1=im1[0]; im2=im2[0] + im1,im2 = self.getShots(shot,calib=calib); im1=im1[0]; im2=im2[0] r = alignment.doShot(im1,im2,initpars,doFit=doFit,show=showInit) im1 = r.im1 im2 = r.im2 @@ -219,7 +282,7 @@ class AnalyzeRun(object): """ if initpars is None: initpars= self.initAlign if shots == "all": shots = slice(self.nShotsPerCalib[calib]) - s1,s2 = self.getShot(shots,calib=calib) + s1,s2 = self.getShots(shots,calib=calib) ret,transformForBestFit = alignment.doShots(s1,s2,initpars=initpars,doFit=doFit,\ returnBestTransform=True,nSaveImg=nSaveImg) if doFit: self.initAlign = transformForBestFit @@ -234,48 +297,34 @@ class AnalyzeRun(object): if len(self.results) == 0: print("self.results are empty, returning without saving") if not os.path.isdir(g_folder_out): os.makedirs(g_folder_out) if fname == "auto": - fname = g_folder_out+"/run%04d_analysis.h5" % self.run + fname = g_folder_out+"/run%04d_analysis" % self.run if os.path.exists(fname) and not overwrite: print("File %s exists, **NOT** saving, use overwrite=True is you want ..."%fname) return - if os.path.exists(fname) and overwrite: os.unlink(fname) - print("Saving results to %s"%fname) - h = h5py.File(fname) - h["roi1"] = (self.roi1.start,self.roi1.stop) - h["roi2"] = (self.roi2.start,self.roi2.stop) - h["scanmot0"] = self.data.scan.scanmotor0 - h["scanpos0"] = self.data.scan.scanmotor0_values - if hasattr(self.data.scan,"scanmotor1"): - h["scanmot1"] = self.data.scan.scanmotor1 - h["scanpos1"] = self.data.scan.scanmotor1_values - #h["transform"] = self.initAlign - for (c,v) in self.results.items(): - cname = "calib%04d/" % c if isinstance(c,int) else "calib%s/" % c - for p,vv in mc.objToDict(v).items(): - # cannot save in hfd5 certain python objects - if p == "fit_result" or p.find("final_transform")==0: - continue - if isinstance(vv,dict): - for pname,parray in vv.items(): - name = cname + p + "/" + pname - h[name] = parray - else: - h[cname + p] = vv - h.close() + def load(self,fname="auto"): + if fname == "auto": fname = g_folder_out+"/run%04d_analysis.npz" % self.run + temp = np.load(fname) + self.results = temp["results"].item() + temp.close() - def saveTransform(self,fname="auto",transform=None): + def _auto_transform_name(self,run=None,calib=None): + if run is None: run = self.run + fname = g_folder_init+"/run%04d_transform" % run + if calib is not None: + fname = fname + "_c%03d" % calib + return fname + ".npy" + + def saveTransform(self,fname="auto",calib=None,transform=None): if transform is None: transform = self.initAlign - if fname == "auto": - fname = g_folder_init+"/run%04d_transform.npy" % self.run + if fname == "auto": fname = self._auto_transform_name(calib=calib) print("Saving roi and transformation parameter to %s"%fname) alignment.saveAlignment(fname,self.initAlign,self.roi1,self.roi2) - def loadTransform(self,fname="auto"): + def loadTransform(self,fname="auto", calib=None): if isinstance(fname,dict): raise FileNotFoundError - if fname == "auto": - fname = g_folder_init+"/run%04d_transform.npy" % self.run - elif isinstance(fname,int): + if fname == "auto": fname = self._auto_transform_name(calib=calib) + if isinstance(fname,int): fname = g_folder_init+"/run%04d_transform.npy" % fname if not os.path.exists(fname): print("Asked to read %s, but it does not exist"%fname) temp = np.load(fname).item()