From 9bcbcea57a16237ea67299d03774167cfce30826 Mon Sep 17 00:00:00 2001 From: msilvafe Date: Tue, 30 Jan 2024 17:58:22 -0500 Subject: [PATCH 01/14] Initial commit for bias_wave lots to go... --- sodetlib/operations/bias_wave.py | 326 +++++++++++++++++++++++++++++++ 1 file changed, 326 insertions(+) create mode 100644 sodetlib/operations/bias_wave.py diff --git a/sodetlib/operations/bias_wave.py b/sodetlib/operations/bias_wave.py new file mode 100644 index 00000000..0d37cbe1 --- /dev/null +++ b/sodetlib/operations/bias_wave.py @@ -0,0 +1,326 @@ +import time +import os +import traceback +import numpy as np +import sodetlib as sdl +import matplotlib.pyplot as plt +import scipy.optimize +from sodetlib.operations import bias_steps + +np.seterr(all='ignore') + +def play_bias_wave(S, bias_group, freqs_wave, amp_wave, duration, + dc_bias=None): + """ + UPDATE THE DOCSTRING + """ + start_times, stop_times = [], [] + + if dc_bias is None: + dc_bias = S.get_tes_bias_bipolar_array()[bias_group] + + for freq in freqs_wave: + S.log(f"BL sine wave with bg={bias_group}, freq={freq}") + S.play_sine_tes(bias_group, amp_wave, freq, dc_amp=dc_bias) + start_times.append(time.time()) + time.sleep(duration) + stop_times.append(time.time()) + S.set_rtm_arb_waveform_enable(0) + S.set_tes_bias_bipolar(bias_group, dc_bias) + return start_times, stop_times + +class BiasWaveAnalysis: + """ + UPDATE THE DOCSTRING + """ + def __init__(self, S=None, cfg=None, bgs=None, run_kwargs=None): + self._S = S + self._cfg = cfg + + self.bgs = bgs + self.am = None + # WHAT DO I REPLACE THIS WITH?? + # self.edge_idxs = None + + # DOES BEN USE THIS? + # self.transition_range = None + + if S is not None: + self.meta = sdl.get_metadata(S, cfg) + self.stream_id = cfg.stream_id + + if run_kwargs is None: + run_kwargs = {} + self.run_kwargs = run_kwargs + self.high_current_mode = run_kwargs.get("high_current_mode", True) + + def save(self, path=None): + data = {} + saved_fields = [ + # Run data and metadata + 'bands', 'channels', 'sid', 'meta', 'run_kwargs', 'start', 'stop', + 'high_current_mode', 'start_times', 'stop_times', + # Bgmap data + 'bgmap', 'polarity', + # Step data and fits + 'resp_times', 'mean_resp', 'step_resp', + # Step fit data + 'step_fit_tmin', 'step_fit_popts', 'step_fit_pcovs', + 'tau_eff', + # Det param data + 'transition_range', 'Ibias', 'Vbias', 'dIbias', 'dVbias', 'dItes', + 'R0', 'I0', 'Pj', 'Si', + # From IV's + 'R_n_IV', 'Rfrac', + ] + + for f in saved_fields: + if not hasattr(self, f): + print(f"WARNING: field {f} does not exist... " + "defaulting to None") + data[f] = getattr(self, f, None) + + if path is not None: + np.save(path, data, allow_pickle=True) + self.filepath = path + else: + self.filepath = sdl.validate_and_save( + 'bias_wave_analysis.npy', data, S=self._S, cfg=self._cfg, + make_path=True + ) + + @classmethod + def load(cls, filepath): + self = cls() + data = np.load(filepath, allow_pickle=True).item() + for k, v in data.items(): + setattr(self, k, v) + self.filepath = filepath + return self + + def run_analysis(self, assignment_thresh=0.3, arc=None, base_dir='/data/so/timestreams', + step_window=0.03, fit_tmin=1.5e-3, transition=None, R0_thresh=30e-3, + save=False, bg_map_file=None): + """ + UPDATE THE DOCSTRING + """ + self._load_am(arc=arc, base_dir=base_dir) + self._split_frequencies() + if bg_map_file is not None: + self.bgmap, self.polarity = sdl.load_bgmap( + self.bands, self.channels, bg_map_file) + else: + self.bgmap, self.polarity = sdl.load_bgmap( + self.bands, self.channels, self.meta['bgmap_file']) + + # GOT TO HERE + self._get_wave_response() + self._compute_dc_params(R0_thresh=R0_thresh) + + # Load R_n from IV + self.R_n_IV = np.full(self.nchans, np.nan) + # Rfrac determined from R0 and R_n + self.Rfrac = np.full(self.nchans, np.nan) + if self.meta['iv_file'] is not None: + if os.path.exists(self.meta['iv_file']): + iva = iv.IVAnalysis.load(self.meta['iv_file']) + chmap = sdl.map_band_chans( + self.bands, self.channels, iva.bands, iva.channels + ) + self.R_n_IV = iva.R_n[chmap] + self.R_n_IV[chmap == -1] = np.nan + self.Rfrac = self.R0 / self.R_n_IV + + if create_bg_map and save_bg_map and self._S is not None: + # Write bgmap after compute_dc_params because bg-assignment + # will be un-set if resistance estimation is too high. + ts = str(int(time.time())) + data = { + 'bands': self.bands, + 'channels': self.channels, + 'sid': self.sid, + 'meta': self.meta, + 'bgmap': self.bgmap, + 'polarity': self.polarity, + } + path = os.path.join('/data/smurf_data/bias_group_maps', + ts[:5], + self.meta['stream_id'], + f'{ts}_bg_map.npy') + os.makedirs(os.path.dirname(path), exist_ok=True) + sdl.validate_and_save(path, data, S=self._S, cfg=self._cfg, + register=True, make_path=False) + self._cfg.dev.update_experiment({'bgmap_file': path}, + update_file=True) + + + self._fit_tau_effs(tmin=fit_tmin) + + if save: + self.save() + + def _load_am(self, arc=None, base_dir='/data/so/timestreams', fix_timestamps=True): + """ + Attempts to load the axis manager from the sid or return one that's + already loaded. Also sets the `abs_chans` array. + """ + if self.am is None: + if arc: + self.am = arc.load_data(self.start, self.stop, stream_id=self.meta['stream_id']) + else: + self.am = sdl.load_session(self.meta['stream_id'], self.sid, + base_dir=base_dir) + + # Fix up timestamp jitter from timestamping in software + if fix_timestamps: + fsamp, t0 = np.polyfit(self.am.primary['FrameCounter'], + self.am.timestamps, 1) + self.am.timestamps = t0 + self.am.primary['FrameCounter']*fsamp + if "det_info" in self.am: + self.bands = self.am.det_info.smurf.band + self.channels = self.am.det_info.smurf.channel + else: + self.bands = self.am.ch_info.band + self.channels = self.am.ch_info.channel + self.abs_chans = self.bands*512 + self.channels + self.nbgs = len(self.am.biases) + self.nchans = len(self.am.signal) + return self.am + + def _split_frequencies(self, am=None): + """ + UPDATE THE DOCSTRING + """ + if am is None: + am = self.am + + self.start_idxs = np.full(np.shape(self.start_times), np.nan) + self.stop_idxs = np.full(np.shape(self.stop_times), np.nan) + + for bg, bias in enumerate(am.biases): + if np.isnan(self.start_times[bg,:]): + continue + for i, [start, stop] in enumerate(zip(self.start_times[bg,:], self.stop_times[bg,:])): + self.start_idxs[bg, i] = np.argmin(np.abs(am.timestamps - start)) + self.stop_idxs[bg, i] = np.argmin(np.abs(am.timestamps - stop)) + + return self.start_idxs, self.stop_idxs + + def _get_wave_response(self, step_window=0.03, pts_before_step=20, + restrict_to_bg_sweep=False, am=None): + """ + UPDATE THE DOCSTRING + """ + if am is None: + am = self.am + + nchans = len(am.signal) + nbgs = len(am.biases) + n_freqs = np.shape(self.start_idxs)[-1] + npts = np.nanmin(self.stop_idxs-self.start_idxs) + + sigs = np.full((nchans, n_freqs, npts), np.nan) + ts = np.full((nbgs, npts), np.nan) + + A_per_rad = self.meta['pA_per_phi0'] / (2*np.pi) * 1e-12 + for bg in np.unique(self.bgmap): + if bg == -1: + continue + rcs = np.where(self.bgmap == bg)[0] + for i, si in enumerate(self.start_idxs[bg]): + s = slice(si, si + npts) + if np.isnan(ts[bg]).all(): + ts[bg, :] = am.timestamps[s] - am.timestamps[si] + sigs[rcs, i, :] = am.signal[rcs, s] * A_per_rad + + # NEED TO ADD BACK IN POLARITY STUFF HERE DERIVED FROM BGMAP + self.resp_times = ts + self.step_resp = sigs.T + self.mean_resp = np.nanmean(sigs, axis=1) + + return ts, sigs + +@sdl.set_action() +def take_bias_waves(S, cfg, bgs=None, amp_wave=0.05, freqs_wave=[23.0], + duration=10, high_current_mode=True, hcm_wait_time=3, + run_analysis=True, analysis_kwargs=None, dacs='pos', + channel_mask=None, g3_tag=None, stream_subtype='bias_waves', + enable_compression=False, plot_rfrac=True, show_plots=False): + """ + UPDATE DOCSTRING + """ + if bgs is None: + bgs = cfg.dev.exp['active_bgs'] + bgs = np.atleast_1d(bgs) + + # Dumb way to get all run kwargs, but we probably want to save these in + # data object + run_kwargs = { + 'bgs': bgs, 'amp_wave': amp_wave, + 'freqs_wave': freqs_wave, 'duration': duration, + 'high_current_mode': high_current_mode, + 'hcm_wait_time': hcm_wait_time, 'run_analysis': run_analysis, + 'analysis_kwargs': analysis_kwargs, 'channel_mask': channel_mask, + } + + initial_ds_factor = S.get_downsample_factor() + initial_filter_disable = S.get_filter_disable() + initial_dc_biases = S.get_tes_bias_bipolar_array() + + try: + dc_biases = initial_dc_biases + init_current_mode = sdl.get_current_mode_array(S) + if high_current_mode: + dc_biases = dc_biases / S.high_low_current_ratio + amp_wave /= S.high_low_current_ratio + sdl.set_current_mode(S, bgs, 1) + S.log(f"Waiting {hcm_wait_time} sec after switching to hcm") + time.sleep(hcm_wait_time) + + bwa = BiasWaveAnalysis(S, cfg, bgs, run_kwargs=run_kwargs) + + bwa.sid = sdl.stream_g3_on( + S, tag=g3_tag, channel_mask=channel_mask, downsample_factor=1, + filter_disable=True, subtype=stream_subtype, enable_compression=enable_compression + ) + + bwa.start_times = np.full((12, len(freqs_wave)), np.nan) + bwa.stop_times = np.full((12, len(freqs_wave)), np.nan) + + bwa.start = time.time() + + for bg in bgs: + bwa.start_times[bg,:], bwa.stop_times[bg,:] = play_bias_wave(S, cfg, bg, + freqs_wave, + amp_wave, duration) + + bwa.stop = time.time() + + finally: + sdl.stream_g3_off(S) + + # Restores current mode to initial values + sdl.set_current_mode(S, np.where(init_current_mode == 0)[0], 0) + sdl.set_current_mode(S, np.where(init_current_mode == 1)[0], 1) + + S.set_downsample_factor(initial_ds_factor) + S.set_filter_disable(initial_filter_disable) + + if run_analysis: + S.log("Running bias wave analysis") + try: + if analysis_kwargs is None: + analysis_kwargs = {} + bwa.run_analysis(save=True, **analysis_kwargs) + if plot_rfrac: + fig, _ = bias_steps.plot_Rfrac(bwa) + path = sdl.make_filename(S, 'bw_rfrac_summary.png', plot=True) + fig.savefig(path) + S.pub.register_file(path, 'bw_rfrac_summary', plot=True, format='png') + if not show_plots: + plt.close(fig) + except Exception: + print(f"Bias wave analysis failed with exception:") + print(traceback.format_exc()) + + return bwa \ No newline at end of file From 6329c67e0c71d50a82e8b83ced1ef177f8436a0e Mon Sep 17 00:00:00 2001 From: msilvafe Date: Wed, 31 Jan 2024 19:19:51 -0500 Subject: [PATCH 02/14] Finish up functions required to take and analyze dc parameters. --- sodetlib/operations/bias_wave.py | 306 ++++++++++++++++++++++++++----- 1 file changed, 256 insertions(+), 50 deletions(-) diff --git a/sodetlib/operations/bias_wave.py b/sodetlib/operations/bias_wave.py index 0d37cbe1..e2e96b06 100644 --- a/sodetlib/operations/bias_wave.py +++ b/sodetlib/operations/bias_wave.py @@ -5,6 +5,7 @@ import sodetlib as sdl import matplotlib.pyplot as plt import scipy.optimize +from scipy.signal import welch from sodetlib.operations import bias_steps np.seterr(all='ignore') @@ -12,7 +13,27 @@ def play_bias_wave(S, bias_group, freqs_wave, amp_wave, duration, dc_bias=None): """ - UPDATE THE DOCSTRING + Play a sine wave on the bias group. + + Args + ---- + bias_group : int + The bias group + freq_wave : float array + List of frequencies to play. Unit = Hz. + amp_wave : float + Amplitude of sine wave to use. Unit = Volts. + duration : float + Duration each sine wave is played + dc_bias : float, optional + Offset voltage of sine wave. Unit = Volts. + + Returns + ------- + start_times : float list + unix timestamp for the beginning of each sine wave. + stop_times : float list + unix timestamp for the end of each sine wave. """ start_times, stop_times = [], [] @@ -29,9 +50,42 @@ def play_bias_wave(S, bias_group, freqs_wave, amp_wave, duration, S.set_tes_bias_bipolar(bias_group, dc_bias) return start_times, stop_times +def get_amplitudes(f_c, x, fs = 4000, N = 12000, window = 'hann'): + """ + Function for calculating the amplitude of a sine wave. + + Args + ---- + f_c : float + Target frequency to calculate sine wave amplitude at. + x : np.array + Data to analyze + fs : int + Sample rate. Unit = samples/second. + N : int + Number of samples to calculate FFT on. + window : str + Window function to use. See scipy.signal.get_window for a list of + valid window functions to use. Default is ``hann``. + + Returns + ------- + a_peak : float + Amplitudes of sine waves. Shape is len(x) + + COULD ADD HERE OPTION TO GET AMPLITUDE FROM LINEAR REGRESSION INSTEAD OF FFT/Periodogram. + """ + f, p = welch(x, fs = fs, nperseg = N, + scaling = 'spectrum',return_onesided=True, + window = window) + a = np.sqrt(p) + a_rms = a[:, f == f_c] + a_peak = np.sqrt(2)*a_rms + return a_peak + class BiasWaveAnalysis: """ - UPDATE THE DOCSTRING + UPDATE THE DOCSTRING...Maybe we can ask Ben to do this one. """ def __init__(self, S=None, cfg=None, bgs=None, run_kwargs=None): self._S = S @@ -39,11 +93,8 @@ def __init__(self, S=None, cfg=None, bgs=None, run_kwargs=None): self.bgs = bgs self.am = None - # WHAT DO I REPLACE THIS WITH?? - # self.edge_idxs = None - - # DOES BEN USE THIS? - # self.transition_range = None + + self.transition_range = None if S is not None: self.meta = sdl.get_metadata(S, cfg) @@ -63,10 +114,10 @@ def save(self, path=None): # Bgmap data 'bgmap', 'polarity', # Step data and fits - 'resp_times', 'mean_resp', 'step_resp', - # Step fit data - 'step_fit_tmin', 'step_fit_popts', 'step_fit_pcovs', - 'tau_eff', + 'resp_times', 'mean_resp', 'wave_resp', + # Add in tau fit stuff here. + # 'step_fit_tmin', 'step_fit_popts', 'step_fit_pcovs', + # 'tau_eff', # Det param data 'transition_range', 'Ibias', 'Vbias', 'dIbias', 'dVbias', 'dItes', 'R0', 'I0', 'Pj', 'Si', @@ -98,11 +149,31 @@ def load(cls, filepath): self.filepath = filepath return self - def run_analysis(self, assignment_thresh=0.3, arc=None, base_dir='/data/so/timestreams', - step_window=0.03, fit_tmin=1.5e-3, transition=None, R0_thresh=30e-3, - save=False, bg_map_file=None): + def run_analysis(self, arc=None, base_dir='/data/so/timestreams', + R0_thresh=30e-3, save=False, bg_map_file=None): """ - UPDATE THE DOCSTRING + Analyzes data taken with take_bias_waves. + + Parameters + ---------- + arc (optional, G3tSmurf): + G3tSmurf archive. If specified, will attempt to load + axis-manager using archive instead of sid. + base_dir (optiional, str): + Base directory where timestreams are stored. Defaults to + /data/so/timestreams. + R0_thresh (float): + Any channel with resistance greater than R0_thresh will be + unassigned from its bias group under the assumption that it's + crosstalk + save (bool): + If true will save the analysis to a npy file. + bg_map_file (optional, path): + If create_bg_map is false and this file is not None, use this file + to load the bg_map. + + CURRENTLY ONLY INCLUDES CALCULATION OF DC PARAMETERS. + NO TIME CONSTANT FITS. """ self._load_am(arc=arc, base_dir=base_dir) self._split_frequencies() @@ -131,30 +202,8 @@ def run_analysis(self, assignment_thresh=0.3, arc=None, base_dir='/data/so/times self.R_n_IV[chmap == -1] = np.nan self.Rfrac = self.R0 / self.R_n_IV - if create_bg_map and save_bg_map and self._S is not None: - # Write bgmap after compute_dc_params because bg-assignment - # will be un-set if resistance estimation is too high. - ts = str(int(time.time())) - data = { - 'bands': self.bands, - 'channels': self.channels, - 'sid': self.sid, - 'meta': self.meta, - 'bgmap': self.bgmap, - 'polarity': self.polarity, - } - path = os.path.join('/data/smurf_data/bias_group_maps', - ts[:5], - self.meta['stream_id'], - f'{ts}_bg_map.npy') - os.makedirs(os.path.dirname(path), exist_ok=True) - sdl.validate_and_save(path, data, S=self._S, cfg=self._cfg, - register=True, make_path=False) - self._cfg.dev.update_experiment({'bgmap_file': path}, - update_file=True) - - - self._fit_tau_effs(tmin=fit_tmin) + # Can add in function for fitting time constant from multifrequency data here. + # self._fit_tau_effs(tmin=fit_tmin) if save: self.save() @@ -189,7 +238,16 @@ def _load_am(self, arc=None, base_dir='/data/so/timestreams', fix_timestamps=Tru def _split_frequencies(self, am=None): """ - UPDATE THE DOCSTRING + Gets indices for each sine wave frequency on each bias group. + + Returns + ------- + start_idxs : float array + Array of indices for the start of each frequency sine wave. + Shape (n_bias_groups, n_frequencies). + stop_idxs : float array + Array of indices for the end of each frequency sine wave. + Shape (n_bias_groups, n_frequencies). """ if am is None: am = self.am @@ -206,9 +264,19 @@ def _split_frequencies(self, am=None): return self.start_idxs, self.stop_idxs - def _get_wave_response(self, step_window=0.03, pts_before_step=20, - restrict_to_bg_sweep=False, am=None): + def _get_wave_response(self, am=None): """ + Splits up full axis manager into each sine wave. + Here we enforce the number of points in each sine wave to be equal. + + Returns + ------- + ts : float array + Array of timestamps for each sine wave period + Shape (n_bias_groups, n_frequencies, n_pts_in_sine_wave). + sigs : float array + Array of sine response data. + Shape (n_detectors, n_frequencies, n_pts_in_sine_wave). UPDATE THE DOCSTRING """ if am is None: @@ -220,7 +288,7 @@ def _get_wave_response(self, step_window=0.03, pts_before_step=20, npts = np.nanmin(self.stop_idxs-self.start_idxs) sigs = np.full((nchans, n_freqs, npts), np.nan) - ts = np.full((nbgs, npts), np.nan) + ts = np.full((nbgs, n_freqs, npts), np.nan) A_per_rad = self.meta['pA_per_phi0'] / (2*np.pi) * 1e-12 for bg in np.unique(self.bgmap): @@ -230,24 +298,161 @@ def _get_wave_response(self, step_window=0.03, pts_before_step=20, for i, si in enumerate(self.start_idxs[bg]): s = slice(si, si + npts) if np.isnan(ts[bg]).all(): - ts[bg, :] = am.timestamps[s] - am.timestamps[si] + ts[bg, i, :] = am.timestamps[s] - am.timestamps[si] sigs[rcs, i, :] = am.signal[rcs, s] * A_per_rad - # NEED TO ADD BACK IN POLARITY STUFF HERE DERIVED FROM BGMAP self.resp_times = ts - self.step_resp = sigs.T - self.mean_resp = np.nanmean(sigs, axis=1) + self.wave_resp = (sigs.T * self.polarity).T + self.mean_resp = (np.nanmean(sigs, axis=1).T * self.polarity).T return ts, sigs + + def _compute_dc_params(self, R0_thresh=30e-3): + """ + Calculates Ibias, dIbias, and dItes from axis manager, and then + runs the DC param calc to estimate R0, I0, Pj, etc. + If multiple frequency are taken, the DC Params are only calculated + off of the minimum frequency in the array of frequencies. + + Args: + transition: (tuple) + Range of voltage bias values (in low-cur units) where the + "in-transition" resistance calculation should be used. If True, + or False, will use in-transition or normal calc for all + channels. Will default to ``cfg.dev.exp['transition_range']`` or + (1, 8) if that does not exist or if self._cfg is not set. + R0_thresh: (float) + Any channel with resistance greater than R0_thresh will be + unassigned from its bias group under the assumption that it's + crosstalk + + Saves: + Ibias: + Array of shape (nbgs) containing the DC bias current for each + bias group + Vbias: + Array of shape (nbgs) containing the DC bias voltage for each + bias group + dIbias: + Array of shape (nbgs) containing the wave current amplitude for + each bias group + dVbias: + Array of shape (nbgs) containing the wave voltage amplitude for + each bias group + dItes: + Array of shape (nchans) containing the wave current amplitude for + each detector. + """ + nbgs = len(self.am.biases) + nchans = len(am.signal) + npts = np.nanmin(self.stop_idxs-self.start_idxs) + + Ibias = np.full(nbgs, np.nan) + dIbias = np.full(nbgs, 0.0, dtype=float) + dItes = np.full(nchans, np.nan) + + # Compute Ibias and dIbias + bias_line_resistance = self.meta['bias_line_resistance'] + high_low_current_ratio = self.meta['high_low_current_ratio'] + rtm_bit_to_volt = self.meta['rtm_bit_to_volt'] + amp_per_bit = 2 * rtm_bit_to_volt / bias_line_resistance + if self.high_current_mode: + amp_per_bit *= high_low_current_ratio + + for bg in range(nbgs): + if len(self.start_idxs[bg]) == 0: + continue + rcs = np.where(self.bgmap == bg)[0] + s = slice(self.start_idxs[bg][0], self.start_idxs[bg][0] + npts) + Ibias[bg] = np.nanmean(self.am.biases[bg, s]) * amp_per_bit + dIbias[bg] = get_amplitudes(self.run_kwargs['freqs_wave'][0], + self.am.biases[bg, s]) * amp_per_bit + dItes[rcs] = get_amplitudes(self.run_kwargs['freqs_wave'][0], + self.wave_resp[rcs,0,:]) + + self.Ibias = Ibias + self.Vbias = Ibias * bias_line_resistance + self.dIbias = dIbias + self.dVbias = dIbias * bias_line_resistance + self.dItes = dItes + + R0, I0, Pj = bias_steps._compute_R0_I0_Pj() + + Si = -1./(I0 * (R0 - self.meta['R_sh'])) + + # If resistance is too high, most likely crosstalk so just reset + # bg mapping and det params + if R0_thresh is not None: + m = np.abs(R0) > R0_thresh + self.bgmap[m] = -1 + for arr in [R0, I0, Pj, Si]: + arr[m] = np.nan + + self.R0 = R0 + self.I0 = I0 + self.Pj = Pj + self.Si = Si + + return R0, I0, Pj, Si @sdl.set_action() def take_bias_waves(S, cfg, bgs=None, amp_wave=0.05, freqs_wave=[23.0], duration=10, high_current_mode=True, hcm_wait_time=3, - run_analysis=True, analysis_kwargs=None, dacs='pos', - channel_mask=None, g3_tag=None, stream_subtype='bias_waves', + run_analysis=True, analysis_kwargs=None, channel_mask=None, + g3_tag=None, stream_subtype='bias_waves', enable_compression=False, plot_rfrac=True, show_plots=False): """ - UPDATE DOCSTRING + Takes bias wave data at the current DC voltage. Assumes bias lines + are already in low-current mode (if they are in high-current this will + not run correction). This function runs bias waves and returns a + BiasWaveAnalysis object, which can be used to easily view and re-analyze + data. + + Parameters: + S (SmurfControl): + Pysmurf control instance + cfg (DetConfig): + Detconfig instance + bgs ( int, list, optional): + Bias groups to run steps on, defaulting to all 12. Note that the + bias-group mapping generated by the bias step analysis will be + restricted to the bgs set here so if you only run with a small + subset of bias groups, the map might not be correct. + amp_wave (float): + Bias wave amplitude voltage in Low-current-mode units. + This will be divided by the high-low-ratio before running the steps + in high-current mode. + freqs_wave (float): + List of frequencies to take bias wave data at. + duration (float): + Duration in seconds of bias wave frequency + high_current_mode (bool): + If true, switches to high-current-mode. If False, leaves in LCM + which runs through the bias-line filter, so make sure you + extend the step duration to be like >2 sec or something + hcm_wait_time (float): + Time to wait after switching to high-current-mode. + run_analysis (bool): + If True, will attempt to run the analysis to calculate DC params + and tau_eff. If this fails, the analysis object will + still be returned but will not contain all analysis results. + analysis_kwargs (dict, optional): + Keyword arguments to be passed to the BiasWaveAnalysis run_analysis + function. + channel_mask : np.ndarray, optional + Mask containing absolute smurf-channels to write to disk + g3_tag: string, optional + Tag to attach to g3 stream. + stream_subtype : optional, string + Stream subtype for this operation. This will default to 'bias_waves'. + enable_compression: bool, optional + If True, will tell the smurf-streamer to compress G3Frames. Defaults + to False because this dominates frame-processing time for high + data-rate streams. + plot_rfrac : bool + Create rfrac plot, publish it, and save it. Default is True. + show_plots : bool + Show plot in addition to saving when running interactively. Default is False. """ if bgs is None: bgs = cfg.dev.exp['active_bgs'] @@ -255,6 +460,7 @@ def take_bias_waves(S, cfg, bgs=None, amp_wave=0.05, freqs_wave=[23.0], # Dumb way to get all run kwargs, but we probably want to save these in # data object + freqs_wave = np.sort(np.asarray(freqs_wave)) # Enforces lowest frequency first. run_kwargs = { 'bgs': bgs, 'amp_wave': amp_wave, 'freqs_wave': freqs_wave, 'duration': duration, From c0814d4c83aae2eb1f155f2ac95001fc8110bae0 Mon Sep 17 00:00:00 2001 From: msilvafe Date: Thu, 1 Feb 2024 14:31:21 -0500 Subject: [PATCH 03/14] Fix issue in iv fn. --- sodetlib/operations/iv.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sodetlib/operations/iv.py b/sodetlib/operations/iv.py index a8e70c73..04e349f0 100644 --- a/sodetlib/operations/iv.py +++ b/sodetlib/operations/iv.py @@ -645,7 +645,7 @@ def overbias_and_sweep(bgs, cool_voltage=None): stop_times[bgs, i] = time.time() if ivcfg.high_current_mode: - biases /= S.high_low_current_ratio + ivcfg.biases /= S.high_low_current_ratio try: sid = sdl.stream_g3_on(S, tag=ivcfg.g3_tag, subtype='iv') if ivcfg.run_serially: From 65ca1d367f57f6434030a4ab45e44cd2a386a6cd Mon Sep 17 00:00:00 2001 From: cryo Date: Wed, 21 Feb 2024 16:36:19 -0500 Subject: [PATCH 04/14] Debuggging from testing on satp3 --- sodetlib/operations/bias_wave.py | 73 +++++++++++++++++++++++--------- 1 file changed, 53 insertions(+), 20 deletions(-) diff --git a/sodetlib/operations/bias_wave.py b/sodetlib/operations/bias_wave.py index e2e96b06..fa1e8f1f 100644 --- a/sodetlib/operations/bias_wave.py +++ b/sodetlib/operations/bias_wave.py @@ -6,11 +6,11 @@ import matplotlib.pyplot as plt import scipy.optimize from scipy.signal import welch -from sodetlib.operations import bias_steps +from sodetlib.operations import bias_steps, iv np.seterr(all='ignore') -def play_bias_wave(S, bias_group, freqs_wave, amp_wave, duration, +def play_bias_wave(S, cfg, bias_group, freqs_wave, amp_wave, duration, dc_bias=None): """ Play a sine wave on the bias group. @@ -42,7 +42,9 @@ def play_bias_wave(S, bias_group, freqs_wave, amp_wave, duration, for freq in freqs_wave: S.log(f"BL sine wave with bg={bias_group}, freq={freq}") - S.play_sine_tes(bias_group, amp_wave, freq, dc_amp=dc_bias) + S.play_sine_tes(bias_group=bias_group, + tone_amp=amp_wave, + tone_freq=freq, dc_amp=dc_bias) start_times.append(time.time()) time.sleep(duration) stop_times.append(time.time()) @@ -75,11 +77,15 @@ def get_amplitudes(f_c, x, fs = 4000, N = 12000, window = 'hann'): COULD ADD HERE OPTION TO GET AMPLITUDE FROM LINEAR REGRESSION INSTEAD OF FFT/Periodogram. """ + x = np.atleast_2d(x) f, p = welch(x, fs = fs, nperseg = N, scaling = 'spectrum',return_onesided=True, window = window) a = np.sqrt(p) - a_rms = a[:, f == f_c] + # Ran into problem with f == f_c when nperseg and len(x) are not right. + # Ben to add a function to check this or enforce correct choice somehow. + idx = np.argmin(np.abs(f-f_c)) + a_rms = a[:, idx] a_peak = np.sqrt(2)*a_rms return a_peak @@ -255,8 +261,8 @@ def _split_frequencies(self, am=None): self.start_idxs = np.full(np.shape(self.start_times), np.nan) self.stop_idxs = np.full(np.shape(self.stop_times), np.nan) - for bg, bias in enumerate(am.biases): - if np.isnan(self.start_times[bg,:]): + for bg, bias in enumerate(am.biases[:12]): + if np.all(np.isnan(self.start_times[bg,:])): continue for i, [start, stop] in enumerate(zip(self.start_times[bg,:], self.stop_times[bg,:])): self.start_idxs[bg, i] = np.argmin(np.abs(am.timestamps - start)) @@ -283,9 +289,9 @@ def _get_wave_response(self, am=None): am = self.am nchans = len(am.signal) - nbgs = len(am.biases) + nbgs = 12 n_freqs = np.shape(self.start_idxs)[-1] - npts = np.nanmin(self.stop_idxs-self.start_idxs) + npts = np.nanmin(self.stop_idxs-self.start_idxs).astype(int) sigs = np.full((nchans, n_freqs, npts), np.nan) ts = np.full((nbgs, n_freqs, npts), np.nan) @@ -295,11 +301,10 @@ def _get_wave_response(self, am=None): if bg == -1: continue rcs = np.where(self.bgmap == bg)[0] - for i, si in enumerate(self.start_idxs[bg]): - s = slice(si, si + npts) - if np.isnan(ts[bg]).all(): - ts[bg, i, :] = am.timestamps[s] - am.timestamps[si] - sigs[rcs, i, :] = am.signal[rcs, s] * A_per_rad + for i, si in enumerate(self.start_idxs[bg].astype(int)): + if np.isnan(ts[bg,i]).all(): + ts[bg, i, :] = am.timestamps[si:si+npts] - am.timestamps[si] + sigs[rcs, i, :] = am.signal[rcs, si:si+npts] * A_per_rad self.resp_times = ts self.wave_resp = (sigs.T * self.polarity).T @@ -343,8 +348,8 @@ def _compute_dc_params(self, R0_thresh=30e-3): Array of shape (nchans) containing the wave current amplitude for each detector. """ - nbgs = len(self.am.biases) - nchans = len(am.signal) + nbgs = 12 + nchans = len(self.am.signal) npts = np.nanmin(self.stop_idxs-self.start_idxs) Ibias = np.full(nbgs, np.nan) @@ -363,10 +368,11 @@ def _compute_dc_params(self, R0_thresh=30e-3): if len(self.start_idxs[bg]) == 0: continue rcs = np.where(self.bgmap == bg)[0] - s = slice(self.start_idxs[bg][0], self.start_idxs[bg][0] + npts) + s = slice(int(self.start_idxs[bg][0]), + int(self.start_idxs[bg][0] + npts)) Ibias[bg] = np.nanmean(self.am.biases[bg, s]) * amp_per_bit dIbias[bg] = get_amplitudes(self.run_kwargs['freqs_wave'][0], - self.am.biases[bg, s]) * amp_per_bit + self.am.biases[bg, s])[0] * amp_per_bit dItes[rcs] = get_amplitudes(self.run_kwargs['freqs_wave'][0], self.wave_resp[rcs,0,:]) @@ -376,7 +382,7 @@ def _compute_dc_params(self, R0_thresh=30e-3): self.dVbias = dIbias * bias_line_resistance self.dItes = dItes - R0, I0, Pj = bias_steps._compute_R0_I0_Pj() + R0, I0, Pj = self._compute_R0_I0_Pj() Si = -1./(I0 * (R0 - self.meta['R_sh'])) @@ -395,6 +401,33 @@ def _compute_dc_params(self, R0_thresh=30e-3): return R0, I0, Pj, Si + def _compute_R0_I0_Pj(self): + """ + Computes the DC params R0 I0 and Pj + + Carbon copy from BiasStepAnalysis + """ + Ib = self.Ibias[self.bgmap] + dIb = self.dIbias[self.bgmap] + dItes = self.dItes + Ib[self.bgmap == -1] = np.nan + dIb[self.bgmap == -1] = np.nan + dIrat = dItes / dIb + + R_sh = self.meta["R_sh"] + + I0 = np.zeros_like(dIrat) + I0_nontransition = Ib * dIrat + I0_transition = Ib * dIrat / (2 * dIrat - 1) + I0[dIrat>0] = I0_nontransition[dIrat>0] + I0[dIrat<0] = I0_transition[dIrat<0] + + Pj = I0 * R_sh * (Ib - I0) + R0 = Pj / I0**2 + R0[I0 == 0] = 0 + + return R0, I0, Pj + @sdl.set_action() def take_bias_waves(S, cfg, bgs=None, amp_wave=0.05, freqs_wave=[23.0], duration=10, high_current_mode=True, hcm_wait_time=3, @@ -460,7 +493,7 @@ def take_bias_waves(S, cfg, bgs=None, amp_wave=0.05, freqs_wave=[23.0], # Dumb way to get all run kwargs, but we probably want to save these in # data object - freqs_wave = np.sort(np.asarray(freqs_wave)) # Enforces lowest frequency first. + freqs_wave = np.sort(np.atleast_1d(freqs_wave)) # Enforces lowest frequency first. run_kwargs = { 'bgs': bgs, 'amp_wave': amp_wave, 'freqs_wave': freqs_wave, 'duration': duration, @@ -529,4 +562,4 @@ def take_bias_waves(S, cfg, bgs=None, amp_wave=0.05, freqs_wave=[23.0], print(f"Bias wave analysis failed with exception:") print(traceback.format_exc()) - return bwa \ No newline at end of file + return bwa From ddf97d8f96042bcbef50338536e500d774dabe45 Mon Sep 17 00:00:00 2001 From: bkelle <51238483+bkelle@users.noreply.github.com> Date: Thu, 22 Feb 2024 13:57:42 -0500 Subject: [PATCH 05/14] Update bias_wave.py Fixed a few of the requested changes from Jack: Deleted unused __init__ removed some unused arg docstrings changed start_idxs and stop_idxs assignment to be ints --- sodetlib/operations/bias_wave.py | 36 +++++++------------------------- 1 file changed, 7 insertions(+), 29 deletions(-) diff --git a/sodetlib/operations/bias_wave.py b/sodetlib/operations/bias_wave.py index fa1e8f1f..28a03568 100644 --- a/sodetlib/operations/bias_wave.py +++ b/sodetlib/operations/bias_wave.py @@ -61,7 +61,7 @@ def get_amplitudes(f_c, x, fs = 4000, N = 12000, window = 'hann'): f_c : float Target frequency to calculate sine wave amplitude at. x : np.array - Data to analyze + Data to analyze. Either can be shape nsamp or ndet x nsamp. fs : int Sample rate. Unit = samples/second. N : int @@ -93,23 +93,6 @@ class BiasWaveAnalysis: """ UPDATE THE DOCSTRING...Maybe we can ask Ben to do this one. """ - def __init__(self, S=None, cfg=None, bgs=None, run_kwargs=None): - self._S = S - self._cfg = cfg - - self.bgs = bgs - self.am = None - - self.transition_range = None - - if S is not None: - self.meta = sdl.get_metadata(S, cfg) - self.stream_id = cfg.stream_id - - if run_kwargs is None: - run_kwargs = {} - self.run_kwargs = run_kwargs - self.high_current_mode = run_kwargs.get("high_current_mode", True) def save(self, path=None): data = {} @@ -121,7 +104,7 @@ def save(self, path=None): 'bgmap', 'polarity', # Step data and fits 'resp_times', 'mean_resp', 'wave_resp', - # Add in tau fit stuff here. + # Add in tau fit stuff here. The below commented out params are anticipated to be reported for tau analysis. # 'step_fit_tmin', 'step_fit_popts', 'step_fit_pcovs', # 'tau_eff', # Det param data @@ -258,8 +241,8 @@ def _split_frequencies(self, am=None): if am is None: am = self.am - self.start_idxs = np.full(np.shape(self.start_times), np.nan) - self.stop_idxs = np.full(np.shape(self.stop_times), np.nan) + self.start_idxs = np.full(np.shape(self.start_times), -1) + self.stop_idxs = np.full(np.shape(self.stop_times), -1) for bg, bias in enumerate(am.biases[:12]): if np.all(np.isnan(self.start_times[bg,:])): @@ -291,7 +274,7 @@ def _get_wave_response(self, am=None): nchans = len(am.signal) nbgs = 12 n_freqs = np.shape(self.start_idxs)[-1] - npts = np.nanmin(self.stop_idxs-self.start_idxs).astype(int) + npts = np.nanmin(self.stop_idxs-self.start_idxs) sigs = np.full((nchans, n_freqs, npts), np.nan) ts = np.full((nbgs, n_freqs, npts), np.nan) @@ -301,7 +284,7 @@ def _get_wave_response(self, am=None): if bg == -1: continue rcs = np.where(self.bgmap == bg)[0] - for i, si in enumerate(self.start_idxs[bg].astype(int)): + for i, si in enumerate(self.start_idxs[bg]): if np.isnan(ts[bg,i]).all(): ts[bg, i, :] = am.timestamps[si:si+npts] - am.timestamps[si] sigs[rcs, i, :] = am.signal[rcs, si:si+npts] * A_per_rad @@ -320,12 +303,6 @@ def _compute_dc_params(self, R0_thresh=30e-3): off of the minimum frequency in the array of frequencies. Args: - transition: (tuple) - Range of voltage bias values (in low-cur units) where the - "in-transition" resistance calculation should be used. If True, - or False, will use in-transition or normal calc for all - channels. Will default to ``cfg.dev.exp['transition_range']`` or - (1, 8) if that does not exist or if self._cfg is not set. R0_thresh: (float) Any channel with resistance greater than R0_thresh will be unassigned from its bias group under the assumption that it's @@ -371,6 +348,7 @@ def _compute_dc_params(self, R0_thresh=30e-3): s = slice(int(self.start_idxs[bg][0]), int(self.start_idxs[bg][0] + npts)) Ibias[bg] = np.nanmean(self.am.biases[bg, s]) * amp_per_bit + dIbias[bg] = get_amplitudes(self.run_kwargs['freqs_wave'][0], self.am.biases[bg, s])[0] * amp_per_bit dItes[rcs] = get_amplitudes(self.run_kwargs['freqs_wave'][0], From 3b4fc4c1a93c3a561de27dde3be8fa8986afb0dc Mon Sep 17 00:00:00 2001 From: bkelle <51238483+bkelle@users.noreply.github.com> Date: Thu, 22 Feb 2024 14:07:15 -0500 Subject: [PATCH 06/14] Update 2 bias_wave.py whoops misunderstood Jack's comment on the __init__ modification. Fixing it to align with what he meant here! --- sodetlib/operations/bias_wave.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/sodetlib/operations/bias_wave.py b/sodetlib/operations/bias_wave.py index 28a03568..0f0204f1 100644 --- a/sodetlib/operations/bias_wave.py +++ b/sodetlib/operations/bias_wave.py @@ -93,7 +93,22 @@ class BiasWaveAnalysis: """ UPDATE THE DOCSTRING...Maybe we can ask Ben to do this one. """ + def __init__(self, S=None, cfg=None, bgs=None, run_kwargs=None): + self._S = S + self._cfg = cfg + self.bgs = bgs + self.am = None + + if S is not None: + self.meta = sdl.get_metadata(S, cfg) + self.stream_id = cfg.stream_id + + if run_kwargs is None: + run_kwargs = {} + self.run_kwargs = run_kwargs + self.high_current_mode = run_kwargs.get("high_current_mode", True) + def save(self, path=None): data = {} saved_fields = [ From e162554c612fe92e1df1509bd4e0b1d832ded37a Mon Sep 17 00:00:00 2001 From: bkelle <51238483+bkelle@users.noreply.github.com> Date: Fri, 23 Feb 2024 13:37:19 -0500 Subject: [PATCH 07/14] Update bias_wave.py adding in sign convention for phase on dIrat calculation. Negative sign = in transition. This is a response to Jack's earlier comments! --- sodetlib/operations/bias_wave.py | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/sodetlib/operations/bias_wave.py b/sodetlib/operations/bias_wave.py index 0f0204f1..4ec42846 100644 --- a/sodetlib/operations/bias_wave.py +++ b/sodetlib/operations/bias_wave.py @@ -403,9 +403,22 @@ def _compute_R0_I0_Pj(self): Ib = self.Ibias[self.bgmap] dIb = self.dIbias[self.bgmap] dItes = self.dItes + + #obtain sign of phase information to determine dIrat sign + vects = np.atleast_2d(dIb) + I = np.linalg.inv(np.tensordot(vects, vects, (1, 1))) + coeffs = np.zeros(self.am.dets.count) + + for di in range(self.am.dets.count): + c = np.matmul(np.atleast_2d(dItes[di]), vects.T) + c = np.dot(I, c.T).T + coeffs[di] = c[0] + coeffs = np.sign(coeffs) + Ib[self.bgmap == -1] = np.nan dIb[self.bgmap == -1] = np.nan - dIrat = dItes / dIb + + dIrat = coeffs * (dItes / dIb) R_sh = self.meta["R_sh"] From 1d46da88d63a084e77ca0dd421e70a0c8a748071 Mon Sep 17 00:00:00 2001 From: bkelle <51238483+bkelle@users.noreply.github.com> Date: Fri, 23 Feb 2024 13:41:51 -0500 Subject: [PATCH 08/14] Update bias_wave.py phase sign calculation slight modification to previous commit in the phase sign calculation code. --- sodetlib/operations/bias_wave.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sodetlib/operations/bias_wave.py b/sodetlib/operations/bias_wave.py index 4ec42846..b4bb6430 100644 --- a/sodetlib/operations/bias_wave.py +++ b/sodetlib/operations/bias_wave.py @@ -407,9 +407,9 @@ def _compute_R0_I0_Pj(self): #obtain sign of phase information to determine dIrat sign vects = np.atleast_2d(dIb) I = np.linalg.inv(np.tensordot(vects, vects, (1, 1))) - coeffs = np.zeros(self.am.dets.count) + coeffs = np.zeros(len(self.am.signal)) - for di in range(self.am.dets.count): + for di in range(len(self.am.signal)): c = np.matmul(np.atleast_2d(dItes[di]), vects.T) c = np.dot(I, c.T).T coeffs[di] = c[0] From 1950532cde6f213416a70c64abe3f419597449b8 Mon Sep 17 00:00:00 2001 From: bkelle <51238483+bkelle@users.noreply.github.com> Date: Mon, 26 Feb 2024 11:16:48 -0500 Subject: [PATCH 09/14] write chunked bias signal to bwa object included code to write chunked bias signal to bwa object. Also removed some phase calc code that probably wont work, in advance of testing said code separately. --- sodetlib/operations/bias_wave.py | 22 +++++++--------------- 1 file changed, 7 insertions(+), 15 deletions(-) diff --git a/sodetlib/operations/bias_wave.py b/sodetlib/operations/bias_wave.py index b4bb6430..63fd4e69 100644 --- a/sodetlib/operations/bias_wave.py +++ b/sodetlib/operations/bias_wave.py @@ -117,8 +117,8 @@ def save(self, path=None): 'high_current_mode', 'start_times', 'stop_times', # Bgmap data 'bgmap', 'polarity', - # Step data and fits - 'resp_times', 'mean_resp', 'wave_resp', + # Step data and fits, including chunked bias data + 'resp_times', 'mean_resp', 'wave_resp', 'wave_biases', # Add in tau fit stuff here. The below commented out params are anticipated to be reported for tau analysis. # 'step_fit_tmin', 'step_fit_popts', 'step_fit_pcovs', # 'tau_eff', @@ -292,6 +292,7 @@ def _get_wave_response(self, am=None): npts = np.nanmin(self.stop_idxs-self.start_idxs) sigs = np.full((nchans, n_freqs, npts), np.nan) + biases = np.full((nchans, n_freqs, npts), np.nan) ts = np.full((nbgs, n_freqs, npts), np.nan) A_per_rad = self.meta['pA_per_phi0'] / (2*np.pi) * 1e-12 @@ -303,12 +304,14 @@ def _get_wave_response(self, am=None): if np.isnan(ts[bg,i]).all(): ts[bg, i, :] = am.timestamps[si:si+npts] - am.timestamps[si] sigs[rcs, i, :] = am.signal[rcs, si:si+npts] * A_per_rad + biases[rcs, i, :] = am.biases[rcs, si:si+npts] self.resp_times = ts self.wave_resp = (sigs.T * self.polarity).T self.mean_resp = (np.nanmean(sigs, axis=1).T * self.polarity).T + self.wave_biases = biases - return ts, sigs + return ts, sigs, biases def _compute_dc_params(self, R0_thresh=30e-3): """ @@ -403,22 +406,11 @@ def _compute_R0_I0_Pj(self): Ib = self.Ibias[self.bgmap] dIb = self.dIbias[self.bgmap] dItes = self.dItes - - #obtain sign of phase information to determine dIrat sign - vects = np.atleast_2d(dIb) - I = np.linalg.inv(np.tensordot(vects, vects, (1, 1))) - coeffs = np.zeros(len(self.am.signal)) - - for di in range(len(self.am.signal)): - c = np.matmul(np.atleast_2d(dItes[di]), vects.T) - c = np.dot(I, c.T).T - coeffs[di] = c[0] - coeffs = np.sign(coeffs) Ib[self.bgmap == -1] = np.nan dIb[self.bgmap == -1] = np.nan - dIrat = coeffs * (dItes / dIb) + dIrat = dItes / dIb R_sh = self.meta["R_sh"] From 3bdb5ff1fb4ea0444a6b8aa52a9cdb82e52ed057 Mon Sep 17 00:00:00 2001 From: bkelle <51238483+bkelle@users.noreply.github.com> Date: Mon, 26 Feb 2024 11:51:23 -0500 Subject: [PATCH 10/14] fix writing of chunked bias signal to BiasWaveAnalysis object --- sodetlib/operations/bias_wave.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sodetlib/operations/bias_wave.py b/sodetlib/operations/bias_wave.py index 63fd4e69..c3bc07a9 100644 --- a/sodetlib/operations/bias_wave.py +++ b/sodetlib/operations/bias_wave.py @@ -292,7 +292,7 @@ def _get_wave_response(self, am=None): npts = np.nanmin(self.stop_idxs-self.start_idxs) sigs = np.full((nchans, n_freqs, npts), np.nan) - biases = np.full((nchans, n_freqs, npts), np.nan) + biases = np.full((nbgs, n_freqs, npts), np.nan) ts = np.full((nbgs, n_freqs, npts), np.nan) A_per_rad = self.meta['pA_per_phi0'] / (2*np.pi) * 1e-12 @@ -304,7 +304,7 @@ def _get_wave_response(self, am=None): if np.isnan(ts[bg,i]).all(): ts[bg, i, :] = am.timestamps[si:si+npts] - am.timestamps[si] sigs[rcs, i, :] = am.signal[rcs, si:si+npts] * A_per_rad - biases[rcs, i, :] = am.biases[rcs, si:si+npts] + biases[bg, i, :] = am.biases[bg, si:si+npts] self.resp_times = ts self.wave_resp = (sigs.T * self.polarity).T From 45c4a346e71382355a69132e8b9b9096b11a3fd9 Mon Sep 17 00:00:00 2001 From: bkelle <51238483+bkelle@users.noreply.github.com> Date: Tue, 27 Feb 2024 23:00:28 -0500 Subject: [PATCH 11/14] Update bias_wave.py with deprojection method for amplitude/phase calc I have tested ~I think~ all parts of these updates in a notebook. Using new get amplitude method of deprojection for recovering amplitude and phase data. I didn't yet remove the old get_amplitudes function, it would be trivia to delete it if that's better practice. --- sodetlib/operations/bias_wave.py | 63 ++++++++++++++++++++++++++++---- 1 file changed, 56 insertions(+), 7 deletions(-) diff --git a/sodetlib/operations/bias_wave.py b/sodetlib/operations/bias_wave.py index c3bc07a9..7c3d0e58 100644 --- a/sodetlib/operations/bias_wave.py +++ b/sodetlib/operations/bias_wave.py @@ -74,8 +74,6 @@ def get_amplitudes(f_c, x, fs = 4000, N = 12000, window = 'hann'): ------- a_peak : float Amplitudes of sine waves. Shape is len(x) - - COULD ADD HERE OPTION TO GET AMPLITUDE FROM LINEAR REGRESSION INSTEAD OF FFT/Periodogram. """ x = np.atleast_2d(x) f, p = welch(x, fs = fs, nperseg = N, @@ -89,6 +87,41 @@ def get_amplitudes(f_c, x, fs = 4000, N = 12000, window = 'hann'): a_peak = np.sqrt(2)*a_rms return a_peak + def get_amplitudes_deprojection(f_c, x, ts): + """ + Function for calculating the amplitude and phase of a wave by deprojecting sine and cosine components of desired frequency. + + Args + ---- + f_c : float + Target frequency to calculate wave amplitude at. + x : np.array + Data to analyze. Either can be shape nsamp or ndet x nsamp. + ts : np.array + Timestamps of data to analyze. Shape is nsamp. + + Returns + ------- + a_peak : float + Amplitudes of sine wave of the desired frequency. Shape is len(x) + phase : float + Phase of sine wave of the desired frequency. Shape is len(x) + """ + + vects = np.zeros((2, len(ts)), dtype='float32') + vects[0, :] = np.sin(2*np.pi*f_c*ts) + vects[1, :] = np.cos(2*np.pi*f_c*ts) + I = np.linalg.inv(np.tensordot(vects, vects, (1, 1))) + coeffs = np.matmul(x, vects.T) + coeffs = np.dot(I, coeffs.T).T + coeffs = np.atleast_2d(coeffs) + + a_peak = np.sqrt(coeffs[:,0]**2 + coeffs[:,1]**2) + phase = np.arctan2(coeffs[:,0], coeffs[:,1]) + + + return a_peak, phase + class BiasWaveAnalysis: """ UPDATE THE DOCSTRING...Maybe we can ask Ben to do this one. @@ -351,6 +384,9 @@ def _compute_dc_params(self, R0_thresh=30e-3): dIbias = np.full(nbgs, 0.0, dtype=float) dItes = np.full(nchans, np.nan) + dIbias_phase = np.full(nbgs, 0.0, dtype=float) + dItes_phase = np.full(nchans, np.nan) + # Compute Ibias and dIbias bias_line_resistance = self.meta['bias_line_resistance'] high_low_current_ratio = self.meta['high_low_current_ratio'] @@ -366,17 +402,21 @@ def _compute_dc_params(self, R0_thresh=30e-3): s = slice(int(self.start_idxs[bg][0]), int(self.start_idxs[bg][0] + npts)) Ibias[bg] = np.nanmean(self.am.biases[bg, s]) * amp_per_bit + + dIbias[bg], dIbias_phase[bg] = get_amplitudes_deprojection(self.run_kwargs['freqs_wave'][0], + self.am.biases[bg, s], self.resp_times[bg, 0, :]) + dIbias[bg] = dIbias[bg] * amp_per_bit - dIbias[bg] = get_amplitudes(self.run_kwargs['freqs_wave'][0], - self.am.biases[bg, s])[0] * amp_per_bit - dItes[rcs] = get_amplitudes(self.run_kwargs['freqs_wave'][0], - self.wave_resp[rcs,0,:]) + dItes[rcs], dItes_phase[rcs] = get_amplitudes_deprojection(self.run_kwargs['freqs_wave'][0], + self.wave_resp[rcs,0,:], self.resp_times[bg, 0, :]) self.Ibias = Ibias self.Vbias = Ibias * bias_line_resistance self.dIbias = dIbias + self.dIbias_phase = dIbias_phase self.dVbias = dIbias * bias_line_resistance self.dItes = dItes + self.dItes_phase = dItes_phase R0, I0, Pj = self._compute_R0_I0_Pj() @@ -406,11 +446,20 @@ def _compute_R0_I0_Pj(self): Ib = self.Ibias[self.bgmap] dIb = self.dIbias[self.bgmap] dItes = self.dItes + + dIb_phase = self.dIbias_phase[self.bgmap] + dItes_phase = self.dItes_phase Ib[self.bgmap == -1] = np.nan dIb[self.bgmap == -1] = np.nan + + #sign of phase response, is there a better way to do it? + rel_phase = dItes_phase - dIb_phase + rel_phase = np.isclose(abs(rel_phase), np.pi, rtol = 1e-1) + rel_phase_sign = np.full(rel_phase.shape, 1.0) + rel_phase_sign[np.where(rel_phase == True)] = -1.0 - dIrat = dItes / dIb + dIrat = rel_phase_sign * (dItes / dIb) R_sh = self.meta["R_sh"] From 78f659c2bda646dcae6582f1add93fcb256fa71d Mon Sep 17 00:00:00 2001 From: bkelle <51238483+bkelle@users.noreply.github.com> Date: Wed, 28 Feb 2024 10:51:33 -0500 Subject: [PATCH 12/14] couple quick indentation fixes --- sodetlib/operations/bias_wave.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/sodetlib/operations/bias_wave.py b/sodetlib/operations/bias_wave.py index 7c3d0e58..9ff6d3c1 100644 --- a/sodetlib/operations/bias_wave.py +++ b/sodetlib/operations/bias_wave.py @@ -87,7 +87,7 @@ def get_amplitudes(f_c, x, fs = 4000, N = 12000, window = 'hann'): a_peak = np.sqrt(2)*a_rms return a_peak - def get_amplitudes_deprojection(f_c, x, ts): +def get_amplitudes_deprojection(f_c, x, ts): """ Function for calculating the amplitude and phase of a wave by deprojecting sine and cosine components of desired frequency. @@ -107,7 +107,7 @@ def get_amplitudes_deprojection(f_c, x, ts): phase : float Phase of sine wave of the desired frequency. Shape is len(x) """ - + vects = np.zeros((2, len(ts)), dtype='float32') vects[0, :] = np.sin(2*np.pi*f_c*ts) vects[1, :] = np.cos(2*np.pi*f_c*ts) @@ -115,13 +115,12 @@ def get_amplitudes_deprojection(f_c, x, ts): coeffs = np.matmul(x, vects.T) coeffs = np.dot(I, coeffs.T).T coeffs = np.atleast_2d(coeffs) - + a_peak = np.sqrt(coeffs[:,0]**2 + coeffs[:,1]**2) phase = np.arctan2(coeffs[:,0], coeffs[:,1]) - return a_peak, phase - + class BiasWaveAnalysis: """ UPDATE THE DOCSTRING...Maybe we can ask Ben to do this one. From a8002a437b4e623dd9bfb053433fd9af6010b53d Mon Sep 17 00:00:00 2001 From: bkelle <51238483+bkelle@users.noreply.github.com> Date: Thu, 29 Feb 2024 07:31:44 -0500 Subject: [PATCH 13/14] added docstring for BiasWaveAnalysis --- sodetlib/operations/bias_wave.py | 81 +++++++++++++++++++++++++++++++- 1 file changed, 79 insertions(+), 2 deletions(-) diff --git a/sodetlib/operations/bias_wave.py b/sodetlib/operations/bias_wave.py index 9ff6d3c1..38eb0d85 100644 --- a/sodetlib/operations/bias_wave.py +++ b/sodetlib/operations/bias_wave.py @@ -123,7 +123,84 @@ def get_amplitudes_deprojection(f_c, x, ts): class BiasWaveAnalysis: """ - UPDATE THE DOCSTRING...Maybe we can ask Ben to do this one. + Container to manage analysis of bias waves taken with the take_bias_waves + function. The main function is ``run_analysis`` and will do a series of + analysis procedures to calculate DC detector parameters from a single bias + wave frequency. The method will be extended to calculate tau_eff from multiple + frequencies. Procedure: + + - Loads an axis manager with all the data + - Finds locations of bias wave frequencies for each bias group + - Gets detector response to each bias wave + - Computes DC params R0, I0, Pj, Si from step responses + - Later: fit multiple TES amplitudes with one-pole filter for tau_eff measurement. + + Most analysis inputs and products will be saved to a npy file so they can + be loaded and re-analyzed easily on another computer like simons1. + + To load data from an saved step file, you can run:: + + bwa = BiasWaveAnalysis.load() + + Attributes: + bands (array (int) of shape (ndets)): + Smurf bands for each resonator. + channels (array (int) of shape (ndets)): + Smurf channel for each resonator. + meta (dict) : + Dictionary of smurf metadata. + sid (int) : + Session-id of streaming session + run_kwargs : + input kwargs + start, stop (float): + start and stop time of all steps + high_current_mode (Bool) : + If high-current-mode was used + start_times, stop_times (array (float) of shape (nbgs, nfreqs)) : + Arrays of start and stop times for the start of each frequency sine wave. + bgmap (array (int) of shape (nchans)): + Map from readout channel to assigned bias group. -1 means not + assigned (that the assignment threshold was not met for any of the + 12 bgs) + polarity (array (int) of shape (ndets)): + The sign of each channel. + resp_times (array (float) shape (nbgs, nfreqs, npts)): + Shared timestamps for each of the wave responses in and + with respect to the bg-wave location, with the wave + occuring at t=0. + mean_resp (array (float) shape (nchans, nfreqs, npts)): + Wave response averaged accross all bias steps for a given channel + in Amps. + wave_resp (array (float) shape (nchans, nfreqs, npts)): + Each individual wave response for a given channel in amps. + wave_biases (array (float) shape (nbgs, nfreqs, npts)): + Each individual bias wave for a given channel. + Ibias (array (float) shape (nbgs)): + DC bias current of each bias group (amps). DC Params are + only calculated off of the minimum frequency in the array of frequencies. + Vbias: + DC bias voltage of each bias group (volts in low-current mode) + dIbias (array (float) shape (nbgs)): + Wave current amplitude for each bias group (amps) + dVbias (array (float) shape (nbgs)): + Wave voltage amplitude for each bias group (volts in low-current mode) + dItes (array (float) shape (nchans)): + Array of tes wave response height for each channel (amps) + R0 (array (float) shape (nchans)): + Computed TES resistances for each channel (ohms). + I0 (array (float) shape (nchans)): + Computed TES currents for each channel (amps) + Pj (array (float) shape (nchans)): + Bias power computed for each channel + Si (array (float) shape (nchans)): + Responsivity computed for each channel + R_n_IV (array (float) shape (nchans)): + Array of normal resistances for each channel pulled from IV + in the device cfg. + Rfrac (array (float) shape (nchans)): + Rfrac of each channel, determined from R0 and the channel's normal + resistance. """ def __init__(self, S=None, cfg=None, bgs=None, run_kwargs=None): self._S = S @@ -155,7 +232,7 @@ def save(self, path=None): # 'step_fit_tmin', 'step_fit_popts', 'step_fit_pcovs', # 'tau_eff', # Det param data - 'transition_range', 'Ibias', 'Vbias', 'dIbias', 'dVbias', 'dItes', + 'Ibias', 'Vbias', 'dIbias', 'dVbias', 'dItes', 'R0', 'I0', 'Pj', 'Si', # From IV's 'R_n_IV', 'Rfrac', From 7bf6df683abecf40f5b58c936c5efa8019c368f7 Mon Sep 17 00:00:00 2001 From: bkelle <51238483+bkelle@users.noreply.github.com> Date: Thu, 29 Feb 2024 07:36:19 -0500 Subject: [PATCH 14/14] adjust tolerance for determination of transition/non-transition phase --- sodetlib/operations/bias_wave.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sodetlib/operations/bias_wave.py b/sodetlib/operations/bias_wave.py index 38eb0d85..f80e1b27 100644 --- a/sodetlib/operations/bias_wave.py +++ b/sodetlib/operations/bias_wave.py @@ -531,7 +531,7 @@ def _compute_R0_I0_Pj(self): #sign of phase response, is there a better way to do it? rel_phase = dItes_phase - dIb_phase - rel_phase = np.isclose(abs(rel_phase), np.pi, rtol = 1e-1) + rel_phase = np.isclose(abs(rel_phase), np.pi, atol=3e-1, rtol = 1e-1) #arbitrary cutoff? rel_phase_sign = np.full(rel_phase.shape, 1.0) rel_phase_sign[np.where(rel_phase == True)] = -1.0