feat(holography) : add quality_flag = 4 for CHIME noise source off (#23)

* feat(holograph) : add quality_flag = 4 for CHIME noise source off
(most commonly because 26 m in non-CHIME observing state)

* style: blacken

Co-authored-by: tristpinsm <[email protected]>

ch_util/andata.py

@@ -2549,7 +2549,7 @@ def _unwrap_fpga_counts(data):
 
     # Calculate the length of an FPGA count and the time it takes to wrap
     seconds_per_count = 2.0 * nfreq / (samp_freq_MHz * 1e6)
-    wrap_time = 2 ** 32.0 * seconds_per_count
+    wrap_time = 2**32.0 * seconds_per_count
 
     # Estimate the FPGA initial zero time from the timestamp in the acquisition
     # name, if the acq name is not there, or of the correct format just silently return
@@ -2567,7 +2567,7 @@ def _unwrap_fpga_counts(data):
     num_wraps = np.round((last_wrap - acq_start) / wrap_time).astype(np.uint64)
 
     # Correct the FPGA counts by adding on the counts lost by wrapping
-    fpga_corrected = time_map["fpga_count"] + num_wraps * 2 ** 32
+    fpga_corrected = time_map["fpga_count"] + num_wraps * 2**32
 
     # Create an array to represent the new time dataset, and fill in the corrected values
     _time_dtype = [("fpga_count", np.uint64), ("ctime", np.float64)]
@@ -2604,7 +2604,7 @@ def _timestamp_from_fpga_cpu(cpu_s, cpu_us, fpga_counts):
             diff_cpu = timestamp_cpu[sl] - mean_cpu
             diff_fpga = fpga_counts[sl] - mean_fpga
             slope_num += np.sum(diff_cpu * diff_fpga)
-            slope_den += np.sum(diff_fpga ** 2)
+            slope_den += np.sum(diff_fpga**2)
         slope = slope_num / slope_den
         # Calculate offset in each section.
         for ii in range(len(edge_inds) - 1):
@@ -3522,7 +3522,7 @@ def _insert_gains(data, input_sel):
                 g_real, g_imag = g_data[1:-1:2], g_data[2:-1:2]
                 g_exp = g_data[-1]
 
-                g_full = (g_real + 1.0j * g_imag) * 2 ** g_exp
+                g_full = (g_real + 1.0j * g_imag) * 2**g_exp
 
                 # Select frequencies that are loaded from the file
                 g_sel = g_full[fsel]

ch_util/cal_utils.py

@@ -652,7 +652,7 @@ def _fit(
         # Iterate to obtain model estimate for amplitude
         for kk in range(niter):
 
-            wk = w0 * model_amp ** 2
+            wk = w0 * model_amp**2
 
             if window is not None:
 
@@ -681,7 +681,7 @@ def _fit(
         if np.isnan(center):
             raise RuntimeError("No peak found.")
 
-        wf = w0 * model_amp ** 2
+        wf = w0 * model_amp**2
         if window is not None:
             wf *= (np.abs(ha - center) <= window).astype(np.float)
 
@@ -1028,7 +1028,7 @@ def fit_jac(x, *param):
             jac[:, 0] = tools.invert_no_zero(peak_amplitude) * model
             jac[:, 1] = 8.0 * np.log(2.0) * dx * tools.invert_no_zero(fwhm) ** 2 * model
             jac[:, 2] = (
-                8.0 * np.log(2.0) * dx ** 2 * tools.invert_no_zero(fwhm) ** 3 * model
+                8.0 * np.log(2.0) * dx**2 * tools.invert_no_zero(fwhm) ** 3 * model
             )
             jac[:, 3:] = (
                 self._vander(x, self.poly_deg_phi) * dmodel_dphase[:, np.newaxis]
@@ -1093,7 +1093,7 @@ def _jacobian_amp(self, ha, elementwise=False):
         jac = np.zeros(shp, dtype=ha.dtype)
         jac[slc + (0,)] = tools.invert_no_zero(peak_amplitude)
         jac[slc + (1,)] = 8.0 * np.log(2.0) * dha * tools.invert_no_zero(fwhm) ** 2
-        jac[slc + (2,)] = 8.0 * np.log(2.0) * dha ** 2 * tools.invert_no_zero(fwhm) ** 3
+        jac[slc + (2,)] = 8.0 * np.log(2.0) * dha**2 * tools.invert_no_zero(fwhm) ** 3
 
         return jac
 
@@ -1880,7 +1880,7 @@ def fit_histogram(
 
     # Define the fitting function
     def gauss(x, peak, mu, sigma):
-        return peak * np.exp(-((x - mu) ** 2) / (2.0 * sigma ** 2))
+        return peak * np.exp(-((x - mu) ** 2) / (2.0 * sigma**2))
 
     # Perform the fit
     par, var_par = curve_fit(
@@ -2103,7 +2103,7 @@ def interpolate_gain(freq, gain, weight, flag=None, length_scale=30.0):
             interp_gain[test, ii] = (ypred[:, 0] + ytrain_mu[:, 0]) + 1.0j * (
                 ypred[:, 1] + ytrain_mu[:, 1]
             )
-            interp_weight[test, ii] = tools.invert_no_zero(err_ypred ** 2)
+            interp_weight[test, ii] = tools.invert_no_zero(err_ypred**2)
 
         else:
             # No valid data

ch_util/chan_monitor.py

@@ -381,7 +381,7 @@ def gaus(x, A, mu, sig2):
             ]
         )
         # 1 deg => dt = 1*(pi/180)*(24*3600/2*pi) = 240s
-        sigsqr0 = 1.0 / (4.0 * np.pi ** 2 * (240.0 * np.cos(dec)) ** 2)
+        sigsqr0 = 1.0 / (4.0 * np.pi**2 * (240.0 * np.cos(dec)) ** 2)
         p0 = np.array(
             [
                 [[A0[ii, jj], mu0[ii, jj], sigsqr0] for jj in range(self.Npr)]

ch_util/finder.py

@@ -1361,7 +1361,7 @@ def print_results_summary(self):
                     )
                 size_q = size_q.where(cond & info_cond)
                 try:
-                    s = float(size_q.scalar()) / 1024 ** 2  # MB.
+                    s = float(size_q.scalar()) / 1024**2  # MB.
                 except TypeError:
                     s = 0
                 info = (interval_number, acq.name, offset, length, n_files, s)

ch_util/fluxcat.py

@@ -204,7 +204,7 @@ def fit(self, freq, flux, eflux, flag=None):
             if not isinstance(self.stats, dict):
                 self.stats = {}
             self.stats["ndof"] = len(x) - self.nparam
-            self.stats["chisq"] = np.sum(resid ** 2 / np.diag(C))
+            self.stats["chisq"] = np.sum(resid**2 / np.diag(C))
 
         # Return results
         return self.param, self.param_cov, self.stats
@@ -216,7 +216,7 @@ def _get_x(self, freq):
     @staticmethod
     def _vandermonde(x, nparam):
 
-        return np.vstack(tuple([x ** rank for rank in range(nparam)])).T
+        return np.vstack(tuple([x**rank for rank in range(nparam)])).T
 
     @staticmethod
     def _fit_func(x, *param):
@@ -236,7 +236,7 @@ def _deriv_fit_func(x, *param):
             )
         )
 
-        dfdp = np.array([z * x ** rank for rank in range(len(param))])
+        dfdp = np.array([z * x**rank for rank in range(len(param))])
         dfdp[0] /= param[0]
 
         return dfdp
@@ -584,7 +584,7 @@ def plot(self, legend=True, catalog=True, residuals=False):
 
         fplot = np.logspace(*xrng, num=nplot)
 
-        xrng = [10.0 ** xx for xx in xrng]
+        xrng = [10.0**xx for xx in xrng]
 
         if residuals:
             flux_hat = self.predict_flux(self.freq)

ch_util/holography.py

@@ -29,6 +29,7 @@
 QUALITY_GOOD = 0
 QUALITY_OFFSOURCE = 1
 QUALITY_BADGATING = 2
+QUALITY_NOISEOFF = 4
 ONSOURCE_DIST_TO_FLAG = 0.1
 
 # Tables in the for tracking Holography observations
@@ -74,6 +75,7 @@ class HolographyObservation(base_model):
     )  # maximum of 64 fields. If we need more, use BigBitField
     off_source = quality_flag.flag(QUALITY_OFFSOURCE)
     bad_gating = quality_flag.flag(QUALITY_BADGATING)
+    noise_off = quality_flag.flag(QUALITY_NOISEOFF)
 
     notes = pw.TextField(null=True)
 
@@ -230,7 +232,13 @@ def parse_post_report(cls, post_report_file):
 
     @classmethod
     def create_from_ant_logs(
-        cls, logs, verbose=False, onsource_dist=0.1, notes=None, **kwargs
+        cls,
+        logs,
+        verbose=False,
+        onsource_dist=0.1,
+        notes=None,
+        quality_flag=0,
+        **kwargs,
     ):
         """
         Read John Galt Telescope log files and create an entry in the
@@ -300,6 +308,7 @@ def create_from_ant_logs(
                             "Questionable on source. Mean, STD(offset) : "
                             "{:.3f}, {:.3f}. {}".format(meanoffset, stdoffset, noteout)
                         )
+                    obs["quality_flag"] |= quality_flag
                     if verbose:
                         print(
                             "Times in .ANT log    : {} {}".format(

ch_util/ni_utils.py

@@ -879,7 +879,7 @@ def ni_als(R, g0, Gamma, Upsilon, maxsteps, abs_tol, rel_tol, weighted_als=True)
             LA.pinv(np.dot(C, np.dot(G.conj().T, W_pow2GC)).conj() * W_pow2),
             np.dot(
                 ktrprod(W_pow2GC, W_pow2).conj().T,
-                (R - N).reshape(Nchannels ** 2, order="F"),
+                (R - N).reshape(Nchannels**2, order="F"),
             ),
         )
 

ch_util/plot.py

@@ -328,7 +328,7 @@ def _coerce_data_shape(
     elif part_sel == "imag":
         data = data.imag
     elif part_sel == "mag":
-        data = (data.real ** 2 + data.imag ** 2) ** (0.5)
+        data = (data.real**2 + data.imag**2) ** (0.5)
     elif part_sel == "phase":
         data = np.arctan(data.imag / data.real)
     elif part_sel == "complex":

ch_util/timing.py

@@ -234,7 +234,7 @@ def _interpret_and_read(cls, acq_files, start, stop, datasets, out_group):
         asc = amp_ref * _amplitude_scaling(freq[:, np.newaxis, np.newaxis])
 
         alpha = np.sum(weight * asc * damp, axis=0) * tools.invert_no_zero(
-            np.sum(weight * asc ** 2, axis=0)
+            np.sum(weight * asc**2, axis=0)
         )
 
         for tt, obj in enumerate(objs):
@@ -788,8 +788,8 @@ def get_tau(self, timestamp, ignore_amp=False, interp="linear", extrap_limit=Non
             # Extract the weights
             weight_corrected = _weight_propagation_addition(
                 self.weight_tau[:],
-                self.weight_alpha[:] / self.amp_to_delay ** 2,
-                self.weight_alpha[iref, np.newaxis, :] / self.amp_to_delay ** 2,
+                self.weight_alpha[:] / self.amp_to_delay**2,
+                self.weight_alpha[iref, np.newaxis, :] / self.amp_to_delay**2,
             )
 
             # Interpolate to the requested times
@@ -1170,7 +1170,7 @@ def get_gain(self, freq, inputs, timestamp, **kwargs):
                 # The delay for each input is a linear combination of the
                 # delay from the noise source inputs
                 tau = np.matmul(C, tau)
-                vartau = np.matmul(C ** 2, vartau)
+                vartau = np.matmul(C**2, vartau)
 
             else:
                 # Find the reference for the requested inputs
@@ -1182,7 +1182,7 @@ def get_gain(self, freq, inputs, timestamp, **kwargs):
 
                 tau = np.matmul(C, tau) - sumC * tau[iref, :]
 
-                vartau = np.matmul(C ** 2, vartau) + sumC ** 2 * vartau[iref, :]
+                vartau = np.matmul(C**2, vartau) + sumC**2 * vartau[iref, :]
 
             # Check if we need to correct the delay using the noise source amplitude
             if self.has_amplitude and self.has_coeff_alpha:
@@ -1199,7 +1199,7 @@ def get_gain(self, freq, inputs, timestamp, **kwargs):
                 # amplitude from the noise source inputs
                 tau += np.matmul(Calpha, alpha)
 
-                vartau += np.matmul(Calpha ** 2, varalpha)
+                vartau += np.matmul(Calpha**2, varalpha)
 
             # Scale by 2 pi nu to convert to gain
             gain = np.exp(
@@ -2179,12 +2179,12 @@ def construct_delay_template(
 
         # Construct alpha template
         alpha = np.sum(weight_amp * asc * damp, axis=0) * tools.invert_no_zero(
-            np.sum(weight_amp * asc ** 2, axis=0)
+            np.sum(weight_amp * asc**2, axis=0)
         )
 
         # Construct delay template
         tau = np.sum(weight_phi * omega * dphi, axis=0) * tools.invert_no_zero(
-            np.sum(weight_phi * omega ** 2, axis=0)
+            np.sum(weight_phi * omega**2, axis=0)
         )
 
         # Calculate amplitude residuals
@@ -2234,8 +2234,8 @@ def construct_delay_template(
     num_freq = np.sum(weight_amp > 0.0, axis=0, dtype=np.int)
 
     # Calculate the uncertainties on the fit parameters
-    weight_tau = np.sum(weight_phi * omega ** 2, axis=0)
-    weight_alpha = np.sum(weight_amp * asc ** 2, axis=0)
+    weight_tau = np.sum(weight_phi * omega**2, axis=0)
+    weight_alpha = np.sum(weight_amp * asc**2, axis=0)
 
     # Calculate the average delay over this period using non-linear
     # least squares that is insensitive to phase wrapping
@@ -2501,7 +2501,7 @@ def model_poly_phase(freq, *param):
 
     model_phase = np.zeros_like(freq)
     for pp, par in enumerate(param):
-        model_phase += par * x ** pp
+        model_phase += par * x**pp
 
     model_phase = model_phase % (2.0 * np.pi)
     model_phase -= 2.0 * np.pi * (model_phase > np.pi)
@@ -2547,7 +2547,7 @@ def _func_poly_phase(freq, *param):
 
     model_phase = np.zeros(x.size, dtype=x.dtype)
     for pp, par in enumerate(param):
-        model_phase += par * x ** pp
+        model_phase += par * x**pp
 
     return np.concatenate((np.cos(model_phase), np.sin(model_phase)))
 
@@ -2644,7 +2644,7 @@ def _interpolation_linear(x, y, var, xeval):
     a2 = adx1 * norm
 
     yeval = a1 * y[ind1] + a2 * y[ind2]
-    weval = tools.invert_no_zero(a1 ** 2 * var[ind1] + a2 ** 2 * var[ind2])
+    weval = tools.invert_no_zero(a1**2 * var[ind1] + a2**2 * var[ind2])
 
     return yeval, weval
 

scripts/update_psrcat.py

@@ -100,7 +100,7 @@
                 if np.ma.is_masked(spec[m + "_ERR"])
                 else 1e-3 * spec[m + "_ERR"]
             )
-            entry.add_measurement(float(m[1:]), 1e-3 * spec[m], err, True, u"ATNF")
+            entry.add_measurement(float(m[1:]), 1e-3 * spec[m], err, True, "ATNF")
 
     entry.fit_model()