CP extending data with merge initial commit

autodeer/dataset.py

@@ -285,9 +285,57 @@ def sequence(self):
         return sequence
 
 
+    def merge(self,other):
+        """
+        Merge two datasets into one dataset.
+
+        Handles the following cases:
+        1. Both datasets have the same parameters but different axes and are 1D
+        """
+
+        # Check if the datasets have the same parameters
+
+        dataarray1: xr.DataArray = self._obj
+        dataarray2: xr.DataArray = other
+
+        # check both axes are 1D
+        if len(dataarray1.dims) != 1 or len(dataarray2.dims) != 1:
+            raise ValueError("Both datasets must be 1D")
+
+        keys_check = [
+            'B','LO','reptime','shots','nAvgs','nPcyc','pcyc_name',
+        ]
+
+        for key in keys_check:
+            if key in dataarray1.attrs and key in dataarray2.attrs:
+                if dataarray1.attrs[key] != dataarray2.attrs[key]:
+                    raise ValueError(f"Datasets have different values for {key}, cannot merge")
+            elif key in dataarray1.attrs:
+                print(f"Parameter {key} not found in dataset 2")
+            elif key in dataarray2.attrs:
+                print(f"Parameter {key} not found in dataset 1")
+
+
+        new_data = np.concatenate((dataarray1.data,dataarray2.data),axis=0)
 
+        new_coords = {}
+        for key, coord in dataarray1.coords.items():
+            new_coords[key] = (coord.dims,np.concatenate((coord,dataarray2.coords[key]),axis=0))
+
+        # Sort based on the first coord
+        first_coord = new_coords[list(new_coords.keys())[0]][1]
+        sort_dir = first_coord[-1] - first_coord[0] # check if ascending or descending
+        sort_idx = np.argsort(first_coord)
+        if sort_dir < 0:
+            sort_idx = np.flip(sort_idx)
+
+        new_data = new_data[sort_idx]
+        for key in new_coords.keys():
+            new_coords[key] = (new_coords[key][0],new_coords[key][1][sort_idx])
+
+        new_dataset = xr.DataArray(new_data, dims=dataarray1.dims, coords=new_coords, attrs=dataarray1.attrs)
 
-
+        return new_dataset
 
 
 

autodeer/gui/autoDEER_worker.py

@@ -148,7 +148,7 @@ def run_fsweep(self):
             pi2_pulse=p90, pi_pulse=p180,
         )
 
-        self.interface.launch(fsweep,savename=self.savename("EDFS_Q"),IFgain=1)
+        self.interface.launch(fsweep,savename=self.savename("EDFS_Q"),)
         self.signals.status.emit('Field-sweep running')
         self.interface.terminate_at(SNRCriteria(150),test_interval=self.test_interval)
         while self.interface.isrunning():
@@ -170,7 +170,7 @@ def run_respro(self):
             pi2_pulse=p90, pi_pulse=p180,fwidth=0.15
         )
 
-        self.interface.launch(RPseq,savename=self.savename("ResPro"),IFgain=1)
+        self.interface.launch(RPseq,savename=self.savename("ResPro"),)
         self.signals.status.emit('Resonator Profile running')
 
         self.interface.terminate_at(SNRCriteria(5),test_interval=self.test_interval)
@@ -187,19 +187,20 @@ def run_respro(self):
 
 
 
-    def run_CP_relax(self,dt=10):
+    def run_CP_relax(self,dt=10,tmin=0.7,averages=30,autoStop=True,autoIFGain=True,*kwargs):
         '''
         Initialise the runner function for relaxation. 
         '''
         self.signals.status.emit('Running Carr-Purcell Experiment')
+        print(f"Running Carr-Purcell Experiment with tmin: {tmin} us and dt: {dt} us")
         LO = self.LO
         gyro = self.gyro
         reptime = self.reptime
         shots = int(100*self.noise_mode)
         shots = np.min([shots,25])
         relax = DEERSequence(
-            B=LO/gyro, LO=LO,reptime=reptime,averages=30,shots=shots,
-            tau1=0.3, tau2=0.3, tau3=0.2, dt=15,
+            B=LO/gyro, LO=LO,reptime=reptime,averages=averages,shots=shots,
+            tau1=tmin/2, tau2=tmin/2, tau3=0.2, dt=15,
             exc_pulse=self.pulses['exc_pulse'], ref_pulse=self.pulses['ref_pulse'],
             pump_pulse=self.pulses['pump_pulse'], det_event=self.pulses['det_event']
             )
@@ -209,18 +210,22 @@ def run_CP_relax(self,dt=10):
         else:
             relax.select_pcyc("DC")
             relax.shots.value *= 8 
-        relax._estimate_time();
+        relax._estimate_time()
 
         # relax.pulses[1].scale.value = 0
         # relax.pulses[3].scale.value = 0
-        self.interface.launch(relax,savename=self.savename("CP"),IFgain=1)
-        self.interface.terminate_at(SNRCriteria(50),test_interval=self.test_interval)
+        if not autoIFGain:
+            autoIFGain = self.interface.IFgain
+
+        self.interface.launch(relax,savename=self.savename("CP"),IFgain=autoIFGain)
+        if autoStop:
+            self.interface.terminate_at(SNRCriteria(50),test_interval=self.test_interval)
         while self.interface.isrunning():
             time.sleep(self.updaterate)
         self.signals.relax_result.emit(self.interface.acquire_dataset())
         self.signals.status.emit('Carr-Purcell experiment complete')
 
-    def run_T2_relax(self,dt=10):
+    def run_T2_relax(self,dt=10,tmin=0.5,averages=30,autoStop=True,autoIFGain=True,*kwargs):
         self.signals.status.emit('Running T2 experiment')
         LO = self.LO
         gyro = self.gyro
@@ -229,12 +234,16 @@ def run_T2_relax(self,dt=10):
         shots = np.min([shots,25])
 
         seq = T2RelaxationSequence(
-            B=LO/gyro, LO=LO,reptime=reptime,averages=30,shots=shots,
-            step=dt,dim=500,pi2_pulse=self.pulses['exc_pulse'],
+            B=LO/gyro, LO=LO,reptime=reptime,averages=averages,shots=shots,
+            start=tmin*1e3,step=dt,dim=500,pi2_pulse=self.pulses['exc_pulse'],
             pi_pulse=self.pulses['ref_pulse'], det_event=self.pulses['det_event'])
 
-        self.interface.launch(seq,savename=self.savename("T2_Q"),IFgain=1)
-        self.interface.terminate_at(SNRCriteria(50),test_interval=self.test_interval)
+        if not autoIFGain:
+            autoIFGain = self.interface.IFgain
+
+        self.interface.launch(seq,savename=self.savename("T2_Q"),IFgain=autoIFGain)
+        if autoStop:
+            self.interface.terminate_at(SNRCriteria(50),test_interval=self.test_interval)
         while self.interface.isrunning():
             time.sleep(self.updaterate)
         self.signals.T2_result.emit(self.interface.acquire_dataset())
@@ -255,7 +264,7 @@ def run_2D_relax(self):
             pi_pulse=self.pulses['ref_pulse'], det_event=self.pulses['det_event'])
 
 
-        self.interface.launch(seq,savename=self.savename("2D_DEC"),IFgain=2)
+        self.interface.launch(seq,savename=self.savename("2D_DEC"),)
         self.interface.terminate_at(SNRCriteria(15),test_interval=self.test_interval,verbosity=2,)
         while self.interface.isrunning():
             time.sleep(self.updaterate)
@@ -361,7 +370,7 @@ def run_deer(self,end_criteria,signal, dt=16,shot=50,averages=1000,):
 
         deer._estimate_time();
 
-        self.interface.launch(deer,savename=self.savename(savename_type,savename_suffix),IFgain=2)
+        self.interface.launch(deer,savename=self.savename(savename_type,savename_suffix),)
         time.sleep(30) # Always wait for the experiment to properly start
         with threadpool_limits(limits=self.cores, user_api='blas'):
             self.interface.terminate_at(end_criteria,verbosity=2,test_interval=self.test_interval) # Change criteria backagain
@@ -378,7 +387,7 @@ def run_reptime_opt(self):
         n_shots = int(np.min([int(50*self.noise_mode),50]))
         scan = ReptimeScan(B=LO/self.gyro, LO=LO,reptime=reptime_guess, reptime_max=20e3, averages=10, shots=n_shots,
                            pi2_pulse=p90, pi_pulse=p180)
-        self.interface.launch(scan,savename=f"{self.samplename}_reptimescan",IFgain=1)
+        self.interface.launch(scan,savename=f"{self.samplename}_reptimescan",)
         self.interface.terminate_at(SNRCriteria(45),verbosity=2,test_interval=self.test_interval)
         while self.interface.isrunning():
             time.sleep(self.updaterate)

autodeer/gui/main.py

@@ -724,6 +724,10 @@ def update_relax(self, dataset=None):
         if dataset is None:
             dataset = self.current_data['relax']
         else:
+            if 'relax' in self.current_data:
+                # attempt to merge datasets
+                self.current_data['relax'].epr.merge(dataset)
+
             self.current_data['relax'] = dataset
 
         worker = Worker(relax_process, dataset)
@@ -992,43 +996,47 @@ def check_T2(self, fitresult):
         # Check if the T2 measurment is too short. 
 
         test_result = fitresult.check_decay()
+        test_dt = fitresult.axis[1].values - fitresult.axis[0].values
+        test_dt *= 1e3
 
         if test_result == 0:
             if self.waitCondition is not None: # Wake up the runner thread
                 self.waitCondition.wakeAll()
                 return None
         elif test_result == -1: # The trace needs to be longer
-            test_dt = fitresult.axis[1] - fitresult.axis[0]
-            test_dt *= 1e3
             new_dt = ad.round_step(test_dt*2,self.waveform_precision)
         elif test_result == 1: # The trace needs to be shorter
-            test_dt = fitresult.axis[1] - fitresult.axis[0]
-            test_dt *= 1e3
             new_dt = ad.round_step(test_dt/2,self.waveform_precision)
+
+        new_tmin = fitresult.axis[-1].values
+        new_tmin += new_dt*1e-3
+        nAvgs = fitresult.dataset.attrs['nAvgs']
 
         if self.worker is not None:
-            self.worker.run_T2_relax(dt=new_dt)
+            self.worker.run_T2_relax(dt=new_dt,tmin=new_tmin,averages=nAvgs,autoStop=False)
 
     def check_CP(self, fitresult):
         # Check if the CP measurment is too short. 
 
         test_result = fitresult.check_decay()
 
+        test_dt = fitresult.axis[1].values - fitresult.axis[0].values
+        test_dt *= 1e3
         if test_result == 0:
             if self.waitCondition is not None: # Wake up the runner thread
                 self.waitCondition.wakeAll()
                 return None
         elif test_result == -1: # The trace needs to be longer
-            test_dt = fitresult.axis[1] - fitresult.axis[0]
-            test_dt *= 1e3
             new_dt = ad.round_step(test_dt*2,self.waveform_precision)
         elif test_result == 1: # The trace needs to be shorter
-            test_dt = fitresult.axis[1] - fitresult.axis[0]
-            test_dt *= 1e3
             new_dt = ad.round_step(test_dt/2,self.waveform_precision)
+
+        new_tmin = fitresult.axis[-1].values
+        new_tmin += new_dt*1e-3
+        nAvgs = fitresult.dataset.attrs['nAvgs']
 
         if self.worker is not None:
-            self.worker.run_CP_relax(dt=new_dt)
+            self.worker.run_CP_relax(dt=new_dt,tmin=new_tmin*2,averages=nAvgs,autoStop=False)
 
 
     def refresh_T2(self, fitresult):

autodeer/hardware/ETH_awg.py

@@ -173,47 +173,50 @@ def acquire_dataset_from_matlab(self, verbosity=0,**kwargs):
                 return data
         raise RuntimeError("Data was unable to be retrieved")
 
-    def launch(self, sequence: Sequence , savename: str, *args,**kwargs):
+    def launch(self, sequence: Sequence , savename: str, IFgain=None, *args,**kwargs):
 
-        test_IF = True
-        while test_IF:
-            self.terminate()
-            print(self.IFgain)
-            self.launch_withIFGain(sequence,savename,self.IFgain)
-            scan_time = sequence._estimate_time() / sequence.averages.value
-            check_1stScan = True
-            while check_1stScan:
-                dataset = self.acquire_dataset()
-                time.sleep(np.min([scan_time//10,2]))
-
-                dig_level = dataset.attrs['diglevel'] / (2**11 *sequence.shots.value* sequence.pcyc_dets.shape[0])
-                pos_levels = dig_level * self.IFgain_options / self.IFgain_options[self.IFgain]
-                pos_levels[pos_levels > 0.85] = 0
-                if dig_level == 0:
-                    continue
-                if (pos_levels[self.IFgain] > 0.85) or  (pos_levels[self.IFgain] < 0.03):
-                    best_IFgain = np.argmax(pos_levels)
-                else:
-                    best_IFgain = self.IFgain
-
-                if np.all(pos_levels==0):
-                    log.critical('Saturation detected with IF gain 0. Please check the power levels.')
-                    raise RuntimeError('Saturation detected with IF gain 0. Please check the power levels.')
-                elif (best_IFgain < self.IFgain) and (dataset.nAvgs == 0):
-                    new_IFgain = self.IFgain - 1
-                    log.warning(f"IF gain changed from {self.IFgain} to {new_IFgain}")
-                    self.IFgain = new_IFgain
-                    check_1stScan = False
-                elif (best_IFgain != self.IFgain) and (dataset.nAvgs >= 1):
-                    new_IFgain = self.IFgain +1
-                    log.warning(f"IF gain changed from {self.IFgain} to {new_IFgain}")
-                    self.IFgain = new_IFgain
-                    check_1stScan = False
-                elif dataset.nAvgs >=1:
-                    log.debug(f"IF gain {self.IFgain} is optimal")
-                    check_1stScan = False
-                    test_IF = False
-
+        if IFgain is None:
+            test_IF = True
+            while test_IF:
+                self.terminate()
+                print(self.IFgain)
+                self.launch_withIFGain(sequence,savename,self.IFgain)
+                scan_time = sequence._estimate_time() / sequence.averages.value
+                check_1stScan = True
+                while check_1stScan:
+                    dataset = self.acquire_dataset()
+                    time.sleep(np.min([scan_time//10,2]))
+
+                    dig_level = dataset.attrs['diglevel'] / (2**11 *sequence.shots.value* sequence.pcyc_dets.shape[0])
+                    pos_levels = dig_level * self.IFgain_options / self.IFgain_options[self.IFgain]
+                    pos_levels[pos_levels > 0.85] = 0
+                    if dig_level == 0:
+                        continue
+                    if (pos_levels[self.IFgain] > 0.85) or  (pos_levels[self.IFgain] < 0.03):
+                        best_IFgain = np.argmax(pos_levels)
+                    else:
+                        best_IFgain = self.IFgain
+
+                    if np.all(pos_levels==0):
+                        log.critical('Saturation detected with IF gain 0. Please check the power levels.')
+                        raise RuntimeError('Saturation detected with IF gain 0. Please check the power levels.')
+                    elif (best_IFgain < self.IFgain) and (dataset.nAvgs == 0):
+                        new_IFgain = self.IFgain - 1
+                        log.warning(f"IF gain changed from {self.IFgain} to {new_IFgain}")
+                        self.IFgain = new_IFgain
+                        check_1stScan = False
+                    elif (best_IFgain != self.IFgain) and (dataset.nAvgs >= 1):
+                        new_IFgain = self.IFgain +1
+                        log.warning(f"IF gain changed from {self.IFgain} to {new_IFgain}")
+                        self.IFgain = new_IFgain
+                        check_1stScan = False
+                    elif dataset.nAvgs >=1:
+                        log.debug(f"IF gain {self.IFgain} is optimal")
+                        check_1stScan = False
+                        test_IF = False
+        elif IFgain is not None:
+
+            self.launch_withIFGain(sequence,savename,IFgain)
 
 
     def launch_withIFGain(self, sequence , savename: str, IFgain: int = 0):

test/test_dataset.py

@@ -1,7 +1,7 @@
 from autodeer.classes import *
 from autodeer.dataset import *
 from autodeer.sequences import *
-from autodeer.hardware.dummy import _simulate_field_sweep
+from autodeer.hardware.dummy import _simulate_field_sweep, _simulate_T2
 import numpy as np
 from scipy.signal import hilbert
 import pytest
@@ -86,4 +86,22 @@ def test_calc_snr(dataset_from_sequence):
 def test_MeasurementTime(dataset_from_sequence):
     dataset = dataset_from_sequence
     dset = dataset.epr.correctphase
-    assert dset.epr.MeasurementTime == 0.6
+    assert dset.epr.MeasurementTime == 0.6
+
+def test_merge_dataset(dataset_from_sequence):
+    rng = np.random.default_rng(seed=0)
+    def build_dataset(tmin):
+        seq = T2RelaxationSequence(
+            B=12200, LO=34.4, reptime=3e3, averages=1, shots=100, start=tmin)
+        t, data = _simulate_T2(seq,0.2)
+        data += rng.normal(0,0.04,data.shape) + 1j*rng.normal(0,0.04,data.shape)
+        extra_params = {"nAvgs":1}
+        dset = create_dataset_from_sequence(data,seq, extra_params)
+        return dset
+
+    dset1 = build_dataset(500)
+    dset2 = build_dataset(dset1.tau[-1].values*1e3+10)
+
+    dset = dset1.epr.merge(dset2)
+    assert dset.data.shape[0] == dset1.data.shape[0] + dset2.data.shape[0]
+