From 6edf688702c6f196d6ae884566e5cc4735f17739 Mon Sep 17 00:00:00 2001 From: Hugo Karas Date: Wed, 15 May 2024 08:25:10 +0200 Subject: [PATCH 1/6] Smoothing of 2D Data --- autodeer/Relaxation.py | 40 ++++++++++++++++++++++++++++++++++------ 1 file changed, 34 insertions(+), 6 deletions(-) diff --git a/autodeer/Relaxation.py b/autodeer/Relaxation.py index 9d10d078..94e5d44d 100644 --- a/autodeer/Relaxation.py +++ b/autodeer/Relaxation.py @@ -6,6 +6,7 @@ from scipy.optimize import curve_fit from autodeer.sequences import Sequence import deerlab as dl +from scipy.interpolate import SmoothBivariateSpline # =========================================================================== @@ -365,8 +366,22 @@ def __init__(self, dataset, sequence: Sequence = None) -> None: dataset.epr.correctphasefull self.data = dataset.data self.dataset = dataset + + + def _smooth(self,*args,**kwargs): + train_x = self.axis[0].data + train_y = self.axis[1].data + train_z = self.data.real + train_x, train_y = np.meshgrid(train_x, train_y) + train_x = train_x.flatten() + train_y = train_y.flatten() + train_z = train_z.flatten() - def plot2D(self, contour=True, norm = 'Normal', axs=None, fig=None): + self.spline = SmoothBivariateSpline(train_x, train_y, train_z,*args,**kwargs) + self.data_smooth = self.spline(train_x, train_y,grid=False).reshape(self.data.shape) + return self.data_smooth + + def plot2D(self, contour=True,smooth=False, norm = 'Normal', axs=None, fig=None): """ Create a 2D plot of the 2D relaxation data. @@ -382,8 +397,13 @@ def plot2D(self, contour=True, norm = 'Normal', axs=None, fig=None): The figure to plot to, by default None """ + if smooth is True: + if not hasattr(self,'data_smooth'): + self._smooth() + data = self.data_smooth + else: + data = self.data.real - data = self.data.real if norm == 'Normal': data = data/np.max(data) elif norm == 'tau2': @@ -397,6 +417,7 @@ def plot2D(self, contour=True, norm = 'Normal', axs=None, fig=None): cmap = cm.get_cmap('Purples',lut=None) cmap_contour = cm.get_cmap('Greys_r',lut=None) + axs.pcolormesh(self.axis[0],self.axis[1],data,cmap=cmap) if contour is True: axs.contour(self.axis[0],self.axis[1],data, cmap=cmap_contour) @@ -426,7 +447,9 @@ def plot1D(self,axs=None,fig=None): elif axs is None: axs = fig.subplots(1,1) - data = self.data.real + if not hasattr(self,'data_smooth'): + self._smooth() + data = self.data_smooth data /= np.max(data) optimal_4p = np.argmax(data,axis=1) @@ -474,7 +497,10 @@ def find_optimal(self,type:str,SNR_target, target_time: float, target_step, aver averages = dataset.nAvgs * dataset.shots * dataset.nPcyc target_shrt = dataset.reptime * 1e-6 - data = np.abs(self.data) + if not hasattr(self,'data_smooth'): + self._smooth() + data = self.data_smooth + # data = np.abs(self.data) data /= np.max(data) calc_data = data @@ -509,9 +535,11 @@ def find_optimal(self,type:str,SNR_target, target_time: float, target_step, aver def optimal_tau1(self,tau2=None,): - + if not hasattr(self,'data_smooth'): + self._smooth() + data = self.data_smooth tau2_idx = np.argmin(np.abs(self.axis[1].data - tau2)) - self.tau1 = self.axis[0].data[np.argmax(self.data[:,tau2_idx])] + self.tau1 = self.axis[0].data[np.argmax(data[:,tau2_idx])] return self.tau1 From 21bd2f7448c12f91e717dd8252f234b92d383366 Mon Sep 17 00:00:00 2001 From: Hugo Karas Date: Wed, 15 May 2024 08:25:24 +0200 Subject: [PATCH 2/6] Minor Fixes --- autodeer/classes.py | 2 +- autodeer/dataset.py | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/autodeer/classes.py b/autodeer/classes.py index 67d186c0..e9ddbf45 100644 --- a/autodeer/classes.py +++ b/autodeer/classes.py @@ -114,7 +114,7 @@ def terminate_at(self, criterion, test_interval=2, keep_running=True, verbosity= # nAvgs = data.num_scans.value nAvgs = data.attrs['nAvgs'] except AttributeError: - print("WARNING: Dataset missing number of averages(nAvgs)!") + self.log.warning("WARNING: Dataset missing number of averages(nAvgs)!") nAvgs = 1 finally: if nAvgs < 1: diff --git a/autodeer/dataset.py b/autodeer/dataset.py index fae1e8d9..f7374474 100644 --- a/autodeer/dataset.py +++ b/autodeer/dataset.py @@ -166,7 +166,6 @@ class EPRAccessor: def __init__(self, xarray_obj): self._obj = xarray_obj - @property def save(self, filename,type='netCDF'): if type == 'netCDF': From c2887bc761963c816ec12f6ddf6ce929c8a1e22e Mon Sep 17 00:00:00 2001 From: Hugo Karas Date: Wed, 15 May 2024 13:34:45 +0200 Subject: [PATCH 3/6] Refocused2D Smotthing now done with SVD --- autodeer/Relaxation.py | 34 ++++++++++++++++++++++------------ 1 file changed, 22 insertions(+), 12 deletions(-) diff --git a/autodeer/Relaxation.py b/autodeer/Relaxation.py index 94e5d44d..3860f57e 100644 --- a/autodeer/Relaxation.py +++ b/autodeer/Relaxation.py @@ -6,7 +6,7 @@ from scipy.optimize import curve_fit from autodeer.sequences import Sequence import deerlab as dl -from scipy.interpolate import SmoothBivariateSpline +from scipy.linalg import svd # =========================================================================== @@ -368,17 +368,27 @@ def __init__(self, dataset, sequence: Sequence = None) -> None: self.dataset = dataset - def _smooth(self,*args,**kwargs): - train_x = self.axis[0].data - train_y = self.axis[1].data - train_z = self.data.real - train_x, train_y = np.meshgrid(train_x, train_y) - train_x = train_x.flatten() - train_y = train_y.flatten() - train_z = train_z.flatten() - - self.spline = SmoothBivariateSpline(train_x, train_y, train_z,*args,**kwargs) - self.data_smooth = self.spline(train_x, train_y,grid=False).reshape(self.data.shape) + def _smooth(self,elements=3): + """ + Used SVD to smooth the 2D data. + + Parameters + ---------- + elements : int, optional + The number of elements to use in the smoothing, by default 3 + + Returns + ------- + np.ndarray + The smoothed data. + """ + + U,E,V = svd(self.data.real) + E_mod = E.copy() + E_mod[elements:] = 0 + mod_data = U @ np.diag(E_mod) @ V + self.data_smooth = mod_data + return self.data_smooth def plot2D(self, contour=True,smooth=False, norm = 'Normal', axs=None, fig=None): From 2621ebe717d8fd079fb227f633414fd440678f24 Mon Sep 17 00:00:00 2001 From: Hugo Karas Date: Tue, 28 May 2024 12:44:17 +0200 Subject: [PATCH 4/6] Smoothing to EDFS --- autodeer/FieldSweep.py | 29 ++++++++++++++++++++++++++++- 1 file changed, 28 insertions(+), 1 deletion(-) diff --git a/autodeer/FieldSweep.py b/autodeer/FieldSweep.py index e89903be..ce489773 100644 --- a/autodeer/FieldSweep.py +++ b/autodeer/FieldSweep.py @@ -9,6 +9,8 @@ import deerlab as dl from xarray import DataArray from autodeer.colors import primary_colors, ReIm_colors +from scipy.interpolate import UnivariateSpline + def create_Nmodel(mwFreq): """Create the field sweep model for a Nitroxide spin system. @@ -162,8 +164,23 @@ def calc_noise_level(self,SNR_target=30): if level < 0.2: level = 0.2 return level - + def smooth(self,): + """ + Generates a smoothed version of the data using a 1D smoothing spline. + + Returns + ------- + np.ndarray + The smoothed data. + """ + smooth_spl = UnivariateSpline(self.axis, self.data) + smooth_spl_freq = UnivariateSpline(self.fs_x, self.data) + self.smooth_data = smooth_spl(self.axis) + self.func = smooth_spl + self.func_freq = smooth_spl_freq + return self.smooth_data + def fit(self, spintype='N', **kwargs): if spintype != 'N': @@ -245,6 +262,16 @@ def plot(self, norm: bool = True, axis: str = "field", axs=None, fig=None) -> Fi elif axis.lower() == 'freq': axs.plot(self.fs_x, np.flip(data), label='fit',c=primary_colors[0]) + axs.legend() + + elif hasattr(self,"smooth_data"): + if axis.lower() == 'field': + data = self.smooth_data / np.max(np.abs(self.smooth_data)) + axs.plot(self.axis, data, label='smooth',c=primary_colors[0]) + elif axis.lower() == 'freq': + data = self.func_freq(self.fs_x) + data /= np.max(np.abs(data)) + axs.plot(self.fs_x, data, label='smooth',c=primary_colors[0]) axs.legend() return fig From b8992479c8b68a4539d7f0b5214271e35e4f0231 Mon Sep 17 00:00:00 2001 From: Hugo Karas Date: Tue, 28 May 2024 12:44:40 +0200 Subject: [PATCH 5/6] Fixes to main.ui --- autodeer/gui/main.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/autodeer/gui/main.py b/autodeer/gui/main.py index fbb00130..de19cce8 100644 --- a/autodeer/gui/main.py +++ b/autodeer/gui/main.py @@ -691,10 +691,15 @@ def create_relax_figure(self): def refresh_relax_figure(self): - - if 'relax2D' in self.current_results: + if isinstance(self.relax_ax, np.ndarray): self.relax_ax[0].cla() self.relax_ax[1].cla() + else: + self.relax_ax.cla() + + + if 'relax2D' in self.current_results: + fig, axs = plt.subplots(2,1,figsize=(12.5, 6.28),layout='constrained',height_ratios=[2,1]) relax_canvas = FigureCanvas(fig) self.relax_canvas.figure.clear() @@ -705,7 +710,6 @@ def refresh_relax_figure(self): self.current_results['relax2D'].plot2D(axs=self.relax_ax[0],fig=fig) self.current_results['relax2D'].plot1D(axs=self.relax_ax[1],fig=fig) else: - self.relax_ax.cla() fig = self.relax_canvas.figure axs = self.relax_ax relax1D_results = [] @@ -909,9 +913,9 @@ def update_func(x): remaining_time = self.MaxTime.value() - ((time.time() - self.starttime) / (60*60)) self.correction_factor = ad.calc_correction_factor(self.current_results['quickdeer'],self.aim_MNR,self.aim_time) + main_log.info(f"Correction factor {self.correction_factor:.3f}") max_tau = self.current_results['relax'].find_optimal(SNR_target=45/(mod_depth*self.correction_factor), target_time=remaining_time, target_step=dt/1e3) main_log.info(f"Max tau {max_tau:.2f} us") - max_tau = 10.3 tau = np.min([rec_tau/2,max_tau]) self.deer_settings['tau1'] = ad.round_step(tau,self.waveform_precision/1e3) self.deer_settings['tau2'] = ad.round_step(tau,self.waveform_precision/1e3) From 78568688c6d34da8ba0737a4182dc70d58322070 Mon Sep 17 00:00:00 2001 From: Hugo Karas Date: Tue, 28 May 2024 12:45:33 +0200 Subject: [PATCH 6/6] Fixes in pulses --- autodeer/dataset.py | 2 +- autodeer/pulses.py | 39 ++++++++++++++++++++++++++------------- 2 files changed, 27 insertions(+), 14 deletions(-) diff --git a/autodeer/dataset.py b/autodeer/dataset.py index f7374474..795c568d 100644 --- a/autodeer/dataset.py +++ b/autodeer/dataset.py @@ -169,7 +169,7 @@ def __init__(self, xarray_obj): def save(self, filename,type='netCDF'): if type == 'netCDF': - self._obj.to_netcdf(f"{filename}.epr") + self._obj.to_netcdf(f"{filename}.h5",engine='h5netcdf',invalid_netcdf=True) @property def correctphase(self): diff --git a/autodeer/pulses.py b/autodeer/pulses.py index 7e4206e7..8a9a80e3 100644 --- a/autodeer/pulses.py +++ b/autodeer/pulses.py @@ -283,7 +283,7 @@ def exciteprofile(self, freqs=None): if isinstance(self, ChirpPulse): sweeprate = np.abs(self.final_freq.value-self.init_freq.value) / self.tp.value elif isinstance(self, HSPulse): - sweeprate = self.beta.value * self.BW.value / 2*self.tp.value + sweeprate = self.beta.value * np.abs(self.final_freq.value-self.init_freq.value) / 2*self.tp.value amp_factor = np.sqrt(2*np.pi*Qcrit*sweeprate)/(2*np.pi); @@ -314,7 +314,7 @@ def exciteprofile(self, freqs=None): Density = UPulse @ Density0 @ UPulse.conjugate().T Mag[0, iOffset] = -2 * (Sx @ Density.T).sum().real Mag[1, iOffset] = -2 * (Sy @ Density.T).sum().real - Mag[2, iOffset] = -2 * (Sz @ Density.T).sum().real + Mag[2, iOffset] = -2 * (Sz * Density.T).sum().real return Mag[0,:], Mag[1,:], Mag[2,:] @@ -801,21 +801,34 @@ def func(self, ax): BW = self.BW.value tp = ax.max() - ax.min() tcent = tp / 2 - + ti = ax nx = ax.shape[0] - beta_exp1 = np.log(beta*0.5**(1-order1)) / np.log(beta) - beta_exp2 = np.log(beta*0.5**(1-order2)) / np.log(beta) + # beta_exp1 = np.log(beta*0.5**(1-order1)) / np.log(beta) + # beta_exp2 = np.log(beta*0.5**(1-order2)) / np.log(beta) + # cut = round(nx/2) + # AM = np.zeros(nx) + # AM[0:cut] = 1/np.cosh( + # beta**beta_exp1 * (ax[0:cut]/tp)**order1) + # AM[cut:-1] = 1/np.cosh( + # beta**beta_exp2 * (ax[cut:-1]/tp)**order2) + + # FM = BW * cumulative_trapezoid(AM**2,ax,initial=0) /\ + # np.trapz(AM**2,ax) + self.init_freq.value + sech = lambda x: 1/np.cosh(x) cut = round(nx/2) - AM = np.zeros(nx) - AM[0:cut] = 1/np.cosh( - beta**beta_exp1 * (ax[0:cut]/tp)**order1) - AM[cut:-1] = 1/np.cosh( - beta**beta_exp2 * (ax[cut:-1]/tp)**order2) + AM = np.zeros_like(ti) + AM[:cut] = sech(beta*0.5*(2*ti[:cut]/tp)**order1) + AM[cut:] = sech(beta*0.5*(2*ti[cut:]/tp)**order2) - FM = BW * cumulative_trapezoid(AM**2,ax,initial=0) /\ - np.trapz(AM**2,ax) + self.init_freq.value - return AM, FM + BWinf = (self.final_freq.value - self.init_freq.value) / np.tanh(beta/2) + + freq = (BWinf/2) * np.tanh((beta/tp*ti)) + # phase = 2*np.pi*(BWinf/2) * (tp/beta) * np.log(np.cosh((beta/tp)*ti)) + + # total_phase = phase * 2* np.pi * np.mean([self.init_freq.value,self.final_freq.value]) + + return AM, freq # =============================================================================