Merge pull request #379 from HERA-Team/allow_combine_interleaved_pspec

Unique time-keys  in `UVPSpec` objects based on `time_1_array` and `time_2_array` instead of the current `time_avg_array`

Author: philbull

hera_pspec/grouping.py

@@ -540,10 +540,11 @@ def average_spectra(uvp_in, blpair_groups=None, time_avg=False,
                         for bl in bl_arr])
 
     # Assign arrays and metadata to UVPSpec object
-    uvp.Ntimes = len(np.unique(time_avg_arr))
-    uvp.Nblpairts = len(time_avg_arr)
+    uvp.Ntimes = len(np.unique(np.hstack([time_1, time_2])))
+    uvp.Nbltpairs = len(time_avg_arr)
     uvp.Nblpairs = len(np.unique(blpair_arr))
     uvp.Nbls = len(bl_arr)
+    uvp.Ntpairs = len(set((t1, t2) for t1, t2 in zip(time_1, time_2)))
 
     # Baselines
     uvp.bl_array = bl_arr
@@ -733,7 +734,7 @@ def spherical_average(uvp_in, kbins, bin_widths, blpair_groups=None, time_avg=Fa
         E = np.zeros((uvp.Ntimes, Ndlyblps, Ndlys, uvp.Npols), dtype=np.float64)
 
 
-        # get kperps for this spw: shape (Nblpairts,)
+        # get kperps for this spw: shape (Nbltpairs,)
         kperps = uvp.get_kperps(spw, little_h=True)
 
         # get kparas for this spw: shape (Ndlys,)
@@ -914,7 +915,7 @@ def spherical_average(uvp_in, kbins, bin_widths, blpair_groups=None, time_avg=Fa
     blp = uvp.blpair_array[0]
     blp_inds = uvp.blpair_to_indices(blp)
     uvp.blpair_array = uvp.blpair_array[blp_inds]
-    uvp.Nblpairts = uvp.Ntimes
+    uvp.Nbltpairs = uvp.Ntpairs
     uvp.Nblpairs = 1
     bl_array = np.unique([uvp.antnums_to_bl(an) for an in uvp.blpair_to_antnums(blp)])
     uvp.bl_vecs = np.asarray([uvp.bl_vecs[np.argmin(uvp.bl_array - bl)] for bl in bl_array])

hera_pspec/plot.py

@@ -384,6 +384,7 @@ def delay_waterfall(uvp, blpairs, spw, pol, component='abs-real',
     """
     import matplotlib
     import matplotlib.pyplot as plt
+    import matplotlib.axes
 
     # assert component
     assert component in ['real', 'abs', 'imag', 'abs-real', 'abs-imag'], "Can't parse specified component {}".format(component)
@@ -1076,4 +1077,4 @@ def _round_sigfig(x, up=True):
     if up:
         return np.ceil(10**sigfigs * x) / 10**sigfigs
     else:
-        return np.floor(10**sigfigs * x) / 10**sigfigs
+        return np.floor(10**sigfigs * x) / 10**sigfigs

hera_pspec/pspecdata.py

@@ -243,6 +243,7 @@ def add(self, dsets, wgts, labels=None, dsets_std=None, cals=None, cal_flag=True
                     break
 
         # Store no. frequencies and no. times
+        # We are still only supporting dsets with same number of times
         self.Nfreqs = self.dsets[0].Nfreqs
         self.Ntimes = self.dsets[0].Ntimes
 
@@ -326,8 +327,8 @@ def validate_datasets(self, verbose=True):
 
         # Check phase centers if phase type is phased
         if 'phased' in set(phase_types):
-            phase_ra = [d.phase_center_ra_degrees for d in self.dsets]
-            phase_dec = [d.phase_center_dec_degrees for d in self.dsets]
+            phase_ra = [d.phase_center_app_ra_degrees for d in self.dsets]
+            phase_dec = [d.phase_center_app_dec_degrees for d in self.dsets]
             max_diff_ra = np.max( [np.diff(d)
                                    for d in itertools.combinations(phase_ra, 2)])
             max_diff_dec = np.max([np.diff(d)
@@ -3119,6 +3120,7 @@ def pspec(self, bls1, bls2, dsets, pols, n_dlys=None,
         pols = _pols
 
         # initialize empty lists
+        # We should really pre-allocate arrays instead.
         data_array = odict()
         wgt_array = odict()
         integration_array = odict()
@@ -3484,8 +3486,15 @@ def pspec(self, bls1, bls2, dsets, pols, n_dlys=None,
                                      % (2*np.pi)
         uvp.blpair_array = np.array(blp_arr)
         uvp.Nblpairs = len(np.unique(blp_arr))
-        uvp.Ntimes = len(np.unique(time1))
-        uvp.Nblpairts = len(time1)
+        # Ntimes in a uvpspec object now means the total number of times.
+        # In all possible datasets that could be produced by pspecdata this
+        # is equal to Ntpairs or 2 x Ntpairs if we interleave
+        # but we have given upspec the capability to have this
+        # not necessarily be the same.
+        # Ntimes could still be the number of "average times" which might be useful for noise purposes. 
+        uvp.Ntimes = len(np.unique(np.hstack([uvp.time_1_array, uvp.time_2_array])))
+        uvp.Nbltpairs = len(set([(blp, t1, t2) for blp, t1, t2 in zip(uvp.blpair_array, uvp.time_1_array, uvp.time_2_array)]))
+        uvp.Ntpairs = len(set([(t1, t2) for t1, t2 in zip(uvp.time_1_array, uvp.time_2_array)]))
         bls_arr = sorted(set(bls_arr))
         uvp.bl_array = np.array([uvp.antnums_to_bl(bl) for bl in bls_arr])
         antpos = dict(zip(dset1.antenna_numbers, dset1.antenna_positions))
@@ -3516,9 +3525,9 @@ def pspec(self, bls1, bls2, dsets, pols, n_dlys=None,
         label1 = self.labels[self.dset_idx(dsets[0])]
         label2 = self.labels[self.dset_idx(dsets[1])]
         uvp.labels = sorted(set([label1, label2]))
-        uvp.label_1_array = np.ones((uvp.Nspws, uvp.Nblpairts, uvp.Npols), int) \
+        uvp.label_1_array = np.ones((uvp.Nspws, uvp.Nbltpairs, uvp.Npols), int) \
                             * uvp.labels.index(label1)
-        uvp.label_2_array = np.ones((uvp.Nspws, uvp.Nblpairts, uvp.Npols), int) \
+        uvp.label_2_array = np.ones((uvp.Nspws, uvp.Nbltpairs, uvp.Npols), int) \
                             * uvp.labels.index(label2)
         uvp.labels = np.array(uvp.labels, str)
         uvp.r_params = uvputils.compress_r_params(r_params)

hera_pspec/testing.py

@@ -61,6 +61,8 @@ def build_vanilla_uvpspec(beam=None):
     time_array = np.repeat(time_array, Nblpairs)
     time_1_array = time_array
     time_2_array = time_array
+    Ntpairs = len(set([(t1, t2) for t1, t2 in zip(time_1_array, time_2_array)]))
+    Nbltpairs = len(set([(blp, t1, t2) for blp, t1, t2 in zip(blpair_array, time_1_array, time_2_array)]))
     lst_array = np.repeat(lst_array, Nblpairs)
     lst_1_array = lst_array
     lst_2_array = lst_array
@@ -153,7 +155,8 @@ def build_vanilla_uvpspec(beam=None):
         "Nspwfreqs",
         "Nspws",
         "Nblpairs",
-        "Nblpairts",
+        "Nbltpairs",
+        "Ntpairs",
         "Npols",
         "Ndlys",
         "Nbls",

hera_pspec/tests/test_grouping.py

@@ -130,7 +130,7 @@ def test_average_spectra(self):
             uvp.set_stats("noise", key, error)
 
         # Add the simple error bar (all are set to be one) to stat_array
-        errs = np.ones((uvp.Ntimes, uvp.Ndlys))
+        errs = np.ones((uvp.Ntpairs, uvp.Ndlys))
         for key in keys:
             uvp.set_stats("simple", key, errs)
         blpair_groups = [blpairs]
@@ -158,7 +158,7 @@ def test_average_spectra(self):
         # is 1/sqrt{N} times the error bar on one single sample.
         averaged_stat = uvp_avg_simple_wgts.stats_array["simple"][0][0, 0, 0]
         initial_stat = uvp.stats_array["simple"][0][0, 0, 0] \
-                     / np.sqrt(uvp.Ntimes) / np.sqrt(len(blpairs))
+                     / np.sqrt(uvp.Ntpairs) / np.sqrt(len(blpairs))
         assert np.all(np.isclose(initial_stat, averaged_stat))
 
         # For non-uniform weights, we test the error bar on the average power
@@ -173,7 +173,7 @@ def test_average_spectra(self):
         # it matches initial over sqrt(Nblpairs - 1)
         uvp_inf_var = copy.deepcopy(uvp)
         initial_stat = uvp.get_stats('simple', (0, blpairs[0], 'xx'))
-        inf_var_stat = np.ones((uvp_inf_var.Ntimes, uvp_inf_var.Ndlys)) * np.inf
+        inf_var_stat = np.ones((uvp_inf_var.Ntpairs, uvp_inf_var.Ndlys)) * np.inf
         uvp_inf_var.set_stats('simple', (0, blpairs[1], 'xx'), inf_var_stat)
         uvp_inf_var_avg = uvp_inf_var.average_spectra(blpair_groups=blpair_groups,
                                                       error_weights='simple',
@@ -185,7 +185,7 @@ def test_average_spectra(self):
         # and check that averaged stat for that time is inf (not zero)
         uvp_inf_var = copy.deepcopy(uvp)
         initial_stat = uvp.get_stats('simple', (0, blpairs[0], 'xx'))
-        inf_var_stat = np.ones((uvp_inf_var.Ntimes, uvp_inf_var.Ndlys))
+        inf_var_stat = np.ones((uvp_inf_var.Ntpairs, uvp_inf_var.Ndlys))
         inf_var_stat[0] = np.inf
         for blp in blpairs:
             uvp_inf_var.set_stats('simple', (0, blp, 'xx'), inf_var_stat)
@@ -220,8 +220,8 @@ def test_average_spectra(self):
                                                 normalize_weights=True,
                                                 inplace=False,
                                                 add_to_history='')
-        assert uvp_time_avg.Nblpairts == uvp_time_avg.Nblpairs
-        assert uvp_time_avg.window_function_array[0].shape[0] == uvp_time_avg.Nblpairts
+        assert uvp_time_avg.Nbltpairs == uvp_time_avg.Nblpairs
+        assert uvp_time_avg.window_function_array[0].shape[0] == uvp_time_avg.Nbltpairs
         blpair_groups, blpair_lens, _ = uvp.get_red_blpairs()
 
         # redundant average
@@ -234,7 +234,7 @@ def test_average_spectra(self):
                                                 normalize_weights=True,
                                                 inplace=False,
                                                 add_to_history='')
-        assert uvp_red_avg.Nblpairts == uvp_red_avg.Ntimes
+        assert uvp_red_avg.Nbltpairs == uvp_red_avg.Ntpairs
 
         # both + error_weights
         keys = uvp.get_all_keys()
@@ -292,8 +292,8 @@ def test_bootstrap_average_blpairs(self):
                                                         blpair_groups,
                                                         time_avg=True)
         assert uvp1[0].Nblpairs == 1
-        assert uvp1[0].Ntimes == self.uvp.Ntimes
-        assert uvp2[0].Ntimes == 1
+        assert uvp1[0].Ntpairs == self.uvp.Ntpairs
+        assert uvp2[0].Ntpairs == 1
 
         # Total of weights assigned should equal total no. of blpairs
         assert np.sum(wgts) == np.array(blpair_groups).size
@@ -308,7 +308,7 @@ def test_bootstrap_average_blpairs(self):
         _blpairs = list(np.unique(self.uvp.blpair_array)[:3])
         uvp3 = self.uvp.select(spws=0, inplace=False, blpairs=_blpairs)
 
-        Nt = uvp3.Ntimes
+        Nt = uvp3.Ntpairs
         uvp3.data_array[0][Nt:2*Nt] = uvp3.data_array[0][:Nt]
         uvp3.data_array[0][2*Nt:] = uvp3.data_array[0][:Nt]
         uvp3.integration_array[0][Nt:2*Nt] = uvp3.integration_array[0][:Nt]
@@ -438,7 +438,7 @@ def test_bootstrap_run():
 
     # assert average only has one time and 3 blpairs
     uvp_avg = psc.get_pspec("grp1", "uvp_avg")
-    assert uvp_avg.Ntimes == 1
+    assert uvp_avg.Ntpairs == 1
     assert uvp_avg.Nblpairs == 3
 
     # check avg file history
@@ -507,11 +507,11 @@ def test_spherical_average():
 
     # insert cov_array and stats_array
     uvp.cov_model = 'empirical'
-    uvp.cov_array_real = {s: np.repeat(np.repeat(np.eye(uvp.Ndlys, dtype=np.float64)[None, : , :, None], uvp.Nblpairts, 0), uvp.Npols, -1)
+    uvp.cov_array_real = {s: np.repeat(np.repeat(np.eye(uvp.Ndlys, dtype=np.float64)[None, : , :, None], uvp.Nbltpairs, 0), uvp.Npols, -1)
                         for s in range(uvp.Nspws)}
-    uvp.cov_array_imag = {s: np.repeat(np.repeat(np.eye(uvp.Ndlys, dtype=np.float64)[None, : , :, None], uvp.Nblpairts, 0), uvp.Npols, -1)
+    uvp.cov_array_imag = {s: np.repeat(np.repeat(np.eye(uvp.Ndlys, dtype=np.float64)[None, : , :, None], uvp.Nbltpairs, 0), uvp.Npols, -1)
                         for s in range(uvp.Nspws)}
-    uvp.stats_array = {'err': {s: np.ones((uvp.Nblpairts, uvp.Ndlys, uvp.Npols), dtype=np.complex128)
+    uvp.stats_array = {'err': {s: np.ones((uvp.Nbltpairs, uvp.Ndlys, uvp.Npols), dtype=np.complex128)
                                   for s in range(uvp.Nspws)}}
 
     # try a spherical average
@@ -539,8 +539,8 @@ def test_spherical_average():
         assert np.isclose(np.sqrt(sph.cov_array_real[spw])[:, range(Nk), range(Nk)], 1/np.sqrt(A[spw].sum(axis=1))).all()
         assert np.isclose(sph.stats_array['err'][spw], 1/np.sqrt(A[spw].sum(axis=1))).all()
 
-    # bug check: time_avg_array was not down-selected to new Nblpairts
-    assert sph.time_avg_array.size == sph.Nblpairts
+    # bug check: time_avg_array was not down-selected to new Nbltpairs
+    assert sph.time_avg_array.size == sph.Nbltpairs
 
     # bug check: cov_array_imag was not updated
     assert sph.cov_array_real[0].shape == sph.cov_array_imag[0].shape
@@ -552,7 +552,7 @@ def test_spherical_average():
 
     # try time average
     sph = grouping.spherical_average(uvp, kbins, bin_widths, time_avg=True)
-    assert sph.Ntimes == 1
+    assert sph.Ntpairs == 1
 
     # try weighting by stats_array
     sph = grouping.spherical_average(uvp, kbins, bin_widths, error_weights='err')

hera_pspec/tests/test_noise.py

@@ -178,7 +178,7 @@ def test_analytic_noise():
     frac_ratio = (uvp.stats_array["P_SN"][0] - uvp.stats_array["P_N"][0]) / uvp.stats_array["P_N"][0]
     dlys = uvp.get_dlys(0) * 1e9
     select = np.abs(dlys) > 3000
-    assert np.abs(frac_ratio[:, select].mean()) < 1 / np.sqrt(uvp.Nblpairts)
+    assert np.abs(frac_ratio[:, select].mean()) < 1 / np.sqrt(uvp.Nbltpairs)
 
     # test single time
     uvp.select(times=uvp.time_avg_array[:1], inplace=True)
@@ -187,7 +187,7 @@ def test_analytic_noise():
     frac_ratio = (uvp.stats_array["P_SN"][0] - uvp.stats_array["P_N"][0]) / uvp.stats_array["P_N"][0]
     dlys = uvp.get_dlys(0) * 1e9
     select = np.abs(dlys) > 3000
-    assert np.abs(frac_ratio[:, select].mean()) < 1 / np.sqrt(uvp.Nblpairts)
+    assert np.abs(frac_ratio[:, select].mean()) < 1 / np.sqrt(uvp.Nbltpairs)
 
     # clean up
     os.remove(os.path.join(DATA_PATH, "test_uvd.uvh5"))

hera_pspec/tests/test_plot.py

@@ -108,7 +108,7 @@ def test_plot_average(self):
         # Average over baseline-pairs but keep the time bins intact
         f2 = plot.delay_spectrum(self.uvp, [blps,], spw=0, pol=('xx','xx'),
                                   average_blpairs=True, average_times=False)
-        elements = [(matplotlib.lines.Line2D, self.uvp.Ntimes),]
+        elements = [(matplotlib.lines.Line2D, self.uvp.Ntpairs),]
         assert axes_contains(f2.axes[0], elements)
         plt.close(f2)
 
@@ -215,7 +215,9 @@ def test_delay_spectrum_misc(self):
         for i in range(3):
             # build-up a large uvpspec object
             _uvp = copy.deepcopy(uvp)
-            _uvp.time_avg_array += (i+1)**2
+            _uvp.time_avg_array += 0 # don't change avg time to make sure this is ok.
+            _uvp.time_1_array += (i+1)**2 
+            _uvp.time_2_array += (i+1)**2
             uvp = uvp + _uvp
         pytest.raises(ValueError, plot.delay_spectrum, uvp, uvp.get_blpairs(), 0, 'xx')
 

hera_pspec/tests/test_uvpspec.py

@@ -53,9 +53,9 @@ def _add_optionals(self, uvp):
         uvp.stats_array = odict({stat: odict()})
         for spw in uvp.spw_array:
             ndlys = uvp.get_spw_ranges(spw)[0][-1]
-            uvp.cov_array_real[spw] = np.empty((uvp.Nblpairts, ndlys, ndlys, uvp.Npols), np.float64)
-            uvp.cov_array_imag[spw] = np.empty((uvp.Nblpairts, ndlys, ndlys, uvp.Npols), np.float64)
-            uvp.stats_array[stat][spw] = np.empty((uvp.Nblpairts, ndlys, uvp.Npols), np.complex128)
+            uvp.cov_array_real[spw] = np.empty((uvp.Nbltpairs, ndlys, ndlys, uvp.Npols), np.float64)
+            uvp.cov_array_imag[spw] = np.empty((uvp.Nbltpairs, ndlys, ndlys, uvp.Npols), np.float64)
+            uvp.stats_array[stat][spw] = np.empty((uvp.Nbltpairs, ndlys, uvp.Npols), np.complex128)
         return uvp
 
     def test_param(self):
@@ -591,7 +591,7 @@ def test_get_exact_window_functions(self):
         uvp.get_exact_window_functions(ftbeam=self.ft_file,
                                        inplace=True)
         assert uvp.exact_windows
-        assert uvp.window_function_array[0].shape[0] == uvp.Nblpairts
+        assert uvp.window_function_array[0].shape[0] == uvp.Nbltpairs
         # if not exact window function, array dim is 4
         assert uvp.window_function_array[0].ndim == 5
 
@@ -1028,5 +1028,13 @@ def test_backwards_compatibility_read():
     # test read in of a static test file dated 8/2019
     uvp = uvpspec.UVPSpec()
     uvp.read_hdf5(os.path.join(DATA_PATH, 'test_uvp.h5'))
+    for dattr in uvp._meta_deprecated:
+        with pytest.raises(AttributeError) as excinfo:
+            raise AttributeError("'UVPSpec' object has no attribute")
+        assert "'UVPSpec' object has no attribute" in str(excinfo.value)
+    for dattr in uvp._meta_dsets_deprecated:
+        with pytest.raises(AttributeError) as excinfo:
+            raise AttributeError("'UVPSpec' object has no attribute")
+        assert "'UVPSpec' object has no attribute" in str(excinfo.value)
     # assert check does not fail
     uvp.check()

hera_pspec/tests/test_uvpspec_utils.py

@@ -237,7 +237,7 @@ def test_conj_blpair():
     pytest.raises(ValueError, uvputils._conj_blpair, 102101103104, which='foo')
 
 
-def test_fast_is_in():
+def test_is_in():
     blps = [ 102101103104, 102101103104, 102101103104, 102101103104,
              101102104103, 101102104103, 101102104103, 101102104103,
              102101104103, 102101104103, 102101104103, 102101104103 ]
@@ -246,7 +246,7 @@ def test_fast_is_in():
               0.1, 0.15, 0.3, 0.3, ]
     src_blpts = np.array(list(zip(blps, times)))
 
-    assert uvputils._fast_is_in(src_blpts, [(101102104103, 0.2)])[0]
+    assert uvputils._is_in(src_blpts, [(101102104103, 0.2)])[0]
 
 
 def test_fast_lookup_blpairts():