Update main from 2D_I_GUI

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

autodeer/Relaxation.py

@@ -6,6 +6,7 @@
 from scipy.optimize import curve_fit
 from autodeer.sequences import Sequence
 import deerlab as dl
+from scipy.linalg import svd
 # ===========================================================================
 
 
@@ -365,8 +366,32 @@ def __init__(self, dataset, sequence: Sequence = None) -> None:
         dataset = dataset.epr.correctphasefull
         self.data = dataset.data
         self.dataset = dataset
+
+
+    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, norm = 'Normal', axs=None, fig=None):
+    def plot2D(self, contour=True,smooth=False, norm = 'Normal', axs=None, fig=None):
         """
         Create a 2D plot of the 2D relaxation data.
 
@@ -382,8 +407,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 +427,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 +457,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 +507,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 +545,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            
 
 

autodeer/classes.py

@@ -113,8 +113,9 @@ def terminate_at(self, criterion, test_interval=2, keep_running=True, verbosity=
             try:
                 # nAvgs = data.num_scans.value
                 nAvgs = data.attrs['nAvgs']
+
             except AttributeError or KeyError:
-                print("WARNING: Dataset missing number of averages(nAvgs)!")
+                self.log.warning("WARNING: Dataset missing number of averages(nAvgs)!")
                 nAvgs = 1
             finally:
                 if nAvgs < 1:

autodeer/dataset.py

@@ -166,11 +166,10 @@ class EPRAccessor:
     def __init__(self, xarray_obj):
         self._obj = xarray_obj
 
-    @property
     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):

autodeer/gui/main.py

@@ -913,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)

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
 
 # =============================================================================