diff --git a/autodeer/dataset.py b/autodeer/dataset.py index 914c578..1f8f2c3 100644 --- a/autodeer/dataset.py +++ b/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 diff --git a/autodeer/gui/autoDEER_worker.py b/autodeer/gui/autoDEER_worker.py index 83c8dc7..869286c 100644 --- a/autodeer/gui/autoDEER_worker.py +++ b/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) diff --git a/autodeer/gui/main.py b/autodeer/gui/main.py index b02612c..c3d385d 100644 --- a/autodeer/gui/main.py +++ b/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): diff --git a/autodeer/hardware/ETH_awg.py b/autodeer/hardware/ETH_awg.py index 2f8bad6..6866266 100644 --- a/autodeer/hardware/ETH_awg.py +++ b/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): diff --git a/test/test_dataset.py b/test/test_dataset.py index a9f80fb..3653524 100644 --- a/test/test_dataset.py +++ b/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 \ No newline at end of file + 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] +