diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 1f16fc46..6d06d299 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -16,7 +16,7 @@ jobs: strategy: matrix: os: [ubuntu-latest, macos-latest] - python-version: [3.9, "3.10", "3.11", "3.12"] + python-version: ["3.10", "3.11", "3.12"] fail-fast: false steps: diff --git a/examples/UVWindow_tutorial.ipynb b/examples/UVWindow_tutorial.ipynb index d4391f7e..c3dc0ab4 100755 --- a/examples/UVWindow_tutorial.ipynb +++ b/examples/UVWindow_tutorial.ipynb @@ -119,13 +119,8 @@ "metadata": {}, "outputs": [], "source": [ - "from astropy.io import fits\n", "from astropy import units\n", - "from scipy.interpolate import interp2d, interp1d\n", - "from scipy import integrate\n", - "import sys, matplotlib, os, h5py\n", "import numpy as np\n", - "from itertools import product, combinations\n", "from astropy import units\n", "from multiprocessing import Pool\n", "import time, h5py\n", diff --git a/hera_pspec/pspecbeam.py b/hera_pspec/pspecbeam.py index 287e706b..38873555 100644 --- a/hera_pspec/pspecbeam.py +++ b/hera_pspec/pspecbeam.py @@ -79,7 +79,7 @@ def _compute_pspec_scalar(cosmo, beam_freqs, omega_ratio, pspec_freqs, X2Y = np.array([cosmo.X2Y(z, little_h=little_h) for z in redshifts]) if exact_norm: #Beam and spectral tapering are already taken into account in normalization. We only use averaged X2Y - scalar = integrate.trapz(X2Y, x=integration_freqs)/(np.abs(integration_freqs[-1]-integration_freqs[0])) + scalar = integrate.trapezoid(X2Y, x=integration_freqs)/(np.abs(integration_freqs[-1]-integration_freqs[0])) return scalar # Use linear interpolation to interpolate the frequency-dependent @@ -105,7 +105,7 @@ def _compute_pspec_scalar(cosmo, beam_freqs, omega_ratio, pspec_freqs, # Integrate to get scalar d_inv_scalar = dBpp_over_BpSq * dOpp_over_Op2 / X2Y - scalar = 1. / integrate.trapz(d_inv_scalar, x=integration_freqs) + scalar = 1. / integrate.trapezoid(d_inv_scalar, x=integration_freqs) return scalar diff --git a/hera_pspec/pspecdata.py b/hera_pspec/pspecdata.py index e74e765a..a60ee8af 100644 --- a/hera_pspec/pspecdata.py +++ b/hera_pspec/pspecdata.py @@ -198,9 +198,6 @@ def add(self, dsets, wgts, labels=None, dsets_std=None, cals=None, cal_flag=True for d, w, s in zip(dsets, wgts, dsets_std): if not isinstance(d, UVData): raise TypeError("Only UVData objects can be used as datasets.") - elif not d.future_array_shapes: - warnings.warn('Converting data to future_array_shapes...') - d.use_future_array_shapes() if not isinstance(w, UVData) and w is not None: raise TypeError("Only UVData objects (or None) can be used as " "weights.") @@ -322,7 +319,13 @@ def validate_datasets(self, verbose=True): # Check phase type phase_types = [] - for d in self.dsets: phase_types.append(d.phase_type) + for d in self.dsets: + if len(d.phase_center_catalog) > 1: + raise ValueError( + "phase_center_catalog should contain only one entry per dataset" + ) + + phase_types.append(next(iter(d.phase_center_catalog.values()))['cat_type']) if np.unique(phase_types).size > 1: raise ValueError("all datasets must have the same phase type " "(i.e. 'drift', 'phased', ...)\ncurrent phase " @@ -362,7 +365,7 @@ def check_key_in_dset(self, key, dset_ind): """ #FIXME: Fix this to enable label keys # get iterable - key = uvutils._get_iterable(key) + key = uvutils.tools._get_iterable(key) if isinstance(key, str): key = (key,) @@ -1721,7 +1724,7 @@ def get_unnormed_V(self, key1, key2, model='empirical', exact_norm=False, def get_analytic_covariance(self, key1, key2, M=None, exact_norm=False, pol=False, model='empirical', known_cov=None): - """ + r""" Calculates the auto-covariance matrix for both the real and imaginary parts of bandpowers (i.e., the q vectors and the p vectors). @@ -2242,7 +2245,7 @@ def get_Q_alt(self, mode: int, allow_fft=True, include_extension=False): ) def get_Q_alt_tensor(self, allow_fft=True, include_extension=False): - """ + r""" Gets the (Ndly, Nfreq, Nfreq) tensor of Q_alt matrices. Parameters @@ -2349,7 +2352,7 @@ def get_Q(self, mode): return Q_alt def p_hat(self, M, q): - """ + r""" Optimal estimate of bandpower p_alpha, defined as p_hat = M q_hat. Parameters @@ -2526,7 +2529,7 @@ def delays(self): def scalar(self, polpair, little_h=True, num_steps=2000, beam=None, taper_override='no_override', exact_norm=False): - """ + r""" Computes the scalar function to convert a power spectrum estimate in "telescope units" to cosmological units, using self.spw_range to set spectral window. @@ -2610,7 +2613,7 @@ def scalar(self, polpair, little_h=True, num_steps=2000, beam=None, def scalar_delay_adjustment(self, key1=None, key2=None, sampling=False, Gv=None, Hv=None): - """ + r""" Computes an adjustment factor for the pspec scalar. There are two reasons why this might be needed: @@ -3426,15 +3429,17 @@ def pspec(self, bls1, bls2, dsets, pols, n_dlys=None, # insert time and blpair info only once per blpair if i < 1 and j < 1: # insert time info + # inds is a slice if dset1 has rectangular blts. inds1 = dset1.antpair2ind(bl1, ordered=False) inds2 = dset2.antpair2ind(bl2, ordered=False) + time1.extend(dset1.time_array[inds1]) time2.extend(dset2.time_array[inds2]) lst1.extend(dset1.lst_array[inds1]) lst2.extend(dset2.lst_array[inds2]) # insert blpair info - blp_arr.extend(np.ones_like(inds1, int) \ + blp_arr.extend(np.ones_like(dset1.time_array[inds1], int) \ * uvputils._antnums_to_blpair(blp)) # insert into data and wgts integrations dictionaries @@ -3644,16 +3649,19 @@ def rephase_to_dset(self, dset_index=0, inplace=True): # iterate over dsets for i, dset in enumerate(dsets): + if len(dset.phase_center_catalog) > 1: + raise ValueError("Cannot deal with datasets with more than one phase center") + # don't rephase dataset we are using as our LST anchor if i == dset_index: # even though not phasing this dset, must set to match all other # dsets due to phasing-check validation - dset.phase_type = 'unknown' + dset.phase_center_catalog[0]['cat_type'] = 'unknown' continue # skip if dataset is not drift phased - if dset.phase_type != 'drift': - warnings.warn(f"Skipping dataset {i} b/c it isn't drift phased") + if dset.phase_center_catalog[0]['cat_type'] != 'unprojected': + warnings.warn(f"Skipping dataset {i} b/c it isn't unprojected") # convert UVData to DataContainers. Note this doesn't make # a copy of the data @@ -3685,7 +3693,7 @@ def rephase_to_dset(self, dset_index=0, inplace=True): # set phasing in UVData object to unknown b/c there isn't a single # consistent phasing for the entire data set. - dset.phase_type = 'unknown' + dset.phase_center_catalog[0]['cat_type'] = 'unknown' if inplace is False: return dsets @@ -4579,9 +4587,9 @@ def _load_cals(cnames, logf=None, verbose=True): # read data uvc = UVCal() if isinstance(cfile, str): - uvc.read_calfits(glob.glob(cfile), use_future_array_shapes=True) + uvc.read(glob.glob(cfile)) else: - uvc.read_calfits(cfile, use_future_array_shapes=True) + uvc.read(cfile) uvc.extra_keywords['filename'] = json.dumps(cfile) cals.append(uvc) diff --git a/hera_pspec/testing.py b/hera_pspec/testing.py index d26a988b..cd0059e5 100644 --- a/hera_pspec/testing.py +++ b/hera_pspec/testing.py @@ -447,8 +447,8 @@ def noise_sim( Nextend = int(Nextend) if Nextend > 0: assert ( - data.phase_type == "drift" - ), "data must be drift phased in order to extend along time axis" + data.phase_center_catalog[0]['cat_type'] == "unprojected" + ), "data must be unprojected in order to extend along time axis" data = copy.deepcopy(data) _data = copy.deepcopy(data) dt = np.median(np.diff(np.unique(_data.time_array))) diff --git a/hera_pspec/tests/Scalar_dev2.ipynb b/hera_pspec/tests/Scalar_dev2.ipynb index d627f721..f5da80bd 100644 --- a/hera_pspec/tests/Scalar_dev2.ipynb +++ b/hera_pspec/tests/Scalar_dev2.ipynb @@ -232,7 +232,7 @@ } ], "source": [ - "scalar = 1 / integrate.trapz(d_inv_scalar, pspec_freqs_Hz)\n", + "scalar = 1 / integrate.trapezoid(d_inv_scalar, pspec_freqs_Hz)\n", "print scalar#* (pspec_freqs_Hz[-1] - pspec_freqs_Hz[0])" ] }, diff --git a/hera_pspec/tests/test_pspecdata.py b/hera_pspec/tests/test_pspecdata.py index 86143082..9de35465 100644 --- a/hera_pspec/tests/test_pspecdata.py +++ b/hera_pspec/tests/test_pspecdata.py @@ -3,7 +3,7 @@ import numpy as np import pyuvdata as uv import os, copy, sys -from scipy.integrate import simps, trapz +from scipy.integrate import simpson, trapezoid from .. import pspecdata, pspecbeam, conversions, container, utils, testing, uvwindow from hera_pspec.data import DATA_PATH from pyuvdata import UVData, UVCal, utils as uvutils @@ -192,25 +192,25 @@ def test_init(self): pytest.raises(ValueError, ds.get_G, key, key) pytest.raises(ValueError, ds.get_H, key, key) - # Test conversion if not future_array_shapes - uvd_old = uv.UVData() - uvd_old.read_miriad(os.path.join(DATA_PATH, - "zen.2458042.17772.xx.HH.uvXA"), - use_future_array_shapes=False) - ds = pspecdata.PSpecData() - pytest.warns(UserWarning, ds.add, [uvd_old, uvd_old], [None, None]) def test_add_data(self): """ Test PSpecData add() """ uv = self.d[0] + # proper usage + ds1 = copy.deepcopy(self.ds) + ds1.add(dsets=[uv], wgts=None, labels=None) # test adding non list objects pytest.raises(TypeError, self.ds.add, 1, 1) # test adding non UVData objects pytest.raises(TypeError, self.ds.add, [1], None) pytest.raises(TypeError, self.ds.add, [uv], [1]) pytest.raises(TypeError, self.ds.add, [uv], None, dsets_std=[1]) + # test adding UVData object with old array shape + ds2 = copy.deepcopy(uv) + ds2.data_array = np.tile(uv.data_array, [1,1,1,1]) + pytest.raises(TypeError, self.ds.add, [uv], [1]) # test adding non UVCal for cals pytest.raises(TypeError, self.ds.add, [uv], None, cals=[1]) # test TypeError if dsets is dict but other inputs are not @@ -1077,9 +1077,9 @@ def test_parseval(self): ps_avg = np.fft.fftshift( np.mean(ps[0], axis=1) ) # Calculate integrals for Parseval's theorem - parseval_real = simps(gsq, x) - parseval_ft = dx**2. * simps(fsq, k) - parseval_phat = simps(ps_avg, tau) + parseval_real = simpson(gsq, x) + parseval_ft = dx**2. * simpson(fsq, k) + parseval_phat = simpson(ps_avg, tau) # Report on results for different ways of calculating Parseval integrals print "Parseval's theorem:" @@ -1186,6 +1186,14 @@ def test_validate_datasets(self): pytest.raises(ValueError, ds.validate_datasets) uvd2.phase_to_time(Time(2458042.5, format='jd')) ds.validate_datasets() + # phase_center_catalog should contain only one entry per dataset + uvd3 = copy.deepcopy(self.d[0]) + uvd3.phase_center_catalog[1] = uvd3.phase_center_catalog[0] + ds2 = pspecdata.PSpecData(dsets=[uvd, uvd3], wgts=[None, None]) + pytest.raises(ValueError, ds2.validate_datasets) + # phased data + uvd4 = copy.deepcopy(self.d[0]) + # test polarization ds.validate_pol((0,1), ('xx', 'xx')) @@ -1209,6 +1217,7 @@ def test_rephase_to_dset(self): uvp1 = ds.pspec(bls, bls, (0, 1), pols=('xx','xx'), verbose=False) # rephase and get pspec ds.rephase_to_dset(0) + ds2 = ds.rephase_to_dset(0, inplace=False) uvp2 = ds.pspec(bls, bls, (0, 1), pols=('xx','xx'), verbose=False) blp = (0, ((37,39),(37,39)), ('xx','xx')) assert np.isclose(np.abs(uvp2.get_data(blp)/uvp1.get_data(blp)), 1.0).min() @@ -1375,7 +1384,7 @@ def test_get_analytic_covariance(self): for i in range(2): new = copy.deepcopy(uvd) new.time_array += new.Ntimes * np.diff(np.unique(new.time_array))[0] - new.lst_array = uvutils.get_lst_for_time(new.time_array, *new.telescope_location_lat_lon_alt_degrees) + new.lst_array = uvutils.get_lst_for_time(new.time_array, telescope_loc=new.telescope.location) uvd += new # get redundant baselines @@ -1774,7 +1783,7 @@ def test_normalization(self): # taper window = windows.blackmanharris(len(freqs)) - NEB = Bp / trapz(window**2, x=freqs) + NEB = Bp / trapezoid(window**2, x=freqs) scalar = cosmo.X2Y(np.mean(cosmo.f2z(freqs))) * np.mean(OmegaP**2/OmegaPP) * Bp * NEB data1 = d1.get_data(bls1[0]) data2 = d2.get_data(bls2[0]) @@ -1812,13 +1821,13 @@ def test_broadcast_dset_flags(self): # test single integration being flagged within spw ds = pspecdata.PSpecData(dsets=[copy.deepcopy(uvd), copy.deepcopy(uvd)], wgts=[None, None]) - ds.dsets[0].flag_array[ds.dsets[0].antpair2ind(24, 25, ordered=False)[3], 600, 0] = True + ds.dsets[0].flag_array[ds.dsets[0].antpair2ind(24, 25, ordered=False), 600, 0][3] = True ds.broadcast_dset_flags(spw_ranges=[(400, 800)], time_thresh=0.25, unflag=False) assert ds.dsets[0].get_flags(24, 25)[3, 400:800].all() assert ds.dsets[0].get_flags(24, 25)[3, :].all() == False # test pspec run sets flagged integration to have zero weight - uvd.flag_array[uvd.antpair2ind(24, 25, ordered=False)[3], 400, :] = True + uvd.flag_array[uvd.antpair2ind(24, 25, ordered=False), 400, :][3] = True ds = pspecdata.PSpecData(dsets=[copy.deepcopy(uvd), copy.deepcopy(uvd)], wgts=[None, None]) ds.broadcast_dset_flags(spw_ranges=[(400, 450)], time_thresh=0.25) uvp = ds.pspec([(24, 25), (37, 38), (38, 39)], [(24, 25), (37, 38), (38, 39)], (0, 1), ('xx', 'xx'), @@ -1830,7 +1839,7 @@ def test_broadcast_dset_flags(self): # average spectra avg_uvp = uvp.average_spectra(blpair_groups=[sorted(np.unique(uvp.blpair_array))], time_avg=True, inplace=False) # repeat but change data in flagged portion - ds.dsets[0].data_array[uvd.antpair2ind(24, 25, ordered=False)[3], 400:450, :] *= 100 + ds.dsets[0].data_array[uvd.antpair2ind(24, 25, ordered=False), 400:450, :][3] *= 100 uvp2 = ds.pspec([(24, 25), (37, 38), (38, 39)], [(24, 25), (37, 38), (38, 39)], (0, 1), ('xx', 'xx'), spw_ranges=[(400, 450)], verbose=False) avg_uvp2 = uvp.average_spectra(blpair_groups=[sorted(np.unique(uvp.blpair_array))], time_avg=True, inplace=False) @@ -2004,8 +2013,8 @@ def test_pspec_run(): uvd.read_miriad(fnames[0], use_future_array_shapes=True) # interleave the data by hand, and add some flags in uvd.flag_array[:] = False - uvd.flag_array[uvd.antpair2ind(37, 38, ordered=False)[0], 10, 0] = True - uvd.flag_array[uvd.antpair2ind(37, 38, ordered=False)[:3], 15, 0] = True + uvd.flag_array[uvd.antpair2ind(37, 38, ordered=False), 10, 0][0] = True + uvd.flag_array[uvd.antpair2ind(37, 38, ordered=False), 15, 0][:3] = True uvd1 = uvd.select(times=np.unique(uvd.time_array)[::2], inplace=False) uvd2 = uvd.select(times=np.unique(uvd.time_array)[1::2], inplace=False) if os.path.exists("./out2.h5"): diff --git a/hera_pspec/tests/test_utils.py b/hera_pspec/tests/test_utils.py index 7140752c..82e1488a 100644 --- a/hera_pspec/tests/test_utils.py +++ b/hera_pspec/tests/test_utils.py @@ -64,19 +64,14 @@ def test_load_config(): assert cfg['pspec']['options']['foo'] == None -class Test_Utils(unittest.TestCase): +class Test_Utils: - def setUp(self): + def setup_class(self): # Load data into UVData object self.uvd = UVData() self.uvd.read_miriad(os.path.join(DATA_PATH, "zen.2458042.17772.xx.HH.uvXA"), use_future_array_shapes=True) - # without future array shapes - self.uvd2 = UVData() - self.uvd2.read_miriad(os.path.join(DATA_PATH, - "zen.2458042.17772.xx.HH.uvXA"), - use_future_array_shapes=False) # Load PSpecBeam object beamfile = os.path.join(DATA_PATH, 'HERA_NF_dipole_power.beamfits') self.beam = pspecbeam.PSpecBeamUV(beamfile) @@ -97,7 +92,7 @@ def test_spw_range_from_freqs(self): # Check that type errors and bounds errors are raised pytest.raises(AttributeError, utils.spw_range_from_freqs, np.arange(3), freq_range=(100e6, 110e6)) - for obj in [self.uvd, self.uvd2, self.uvp]: + for obj in [self.uvd, self.uvp]: pytest.raises(ValueError, utils.spw_range_from_freqs, obj, freq_range=(98e6, 110e6)) # lower bound pytest.raises(ValueError, utils.spw_range_from_freqs, obj, @@ -108,7 +103,7 @@ def test_spw_range_from_freqs(self): # Check that valid frequency ranges are returned # with and without future array shapes freq_list = [(100e6, 120e6), (120e6, 140e6), (140e6, 160e6)] - for obj in [self.uvd, self.uvd2]: + for obj in [self.uvd]: spw1 = utils.spw_range_from_freqs(obj, freq_range=(110e6, 130e6)) spw2 = utils.spw_range_from_freqs(obj, freq_range=freq_list) spw3 = utils.spw_range_from_freqs(obj, freq_range=(98e6, 120e6), @@ -136,7 +131,7 @@ def test_spw_range_from_redshifts(self): # Check that type errors and bounds errors are raised pytest.raises(AttributeError, utils.spw_range_from_redshifts, np.arange(3), z_range=(9.7, 12.1)) - for obj in [self.uvd, self.uvd2, self.uvp]: + for obj in [self.uvd, self.uvp]: pytest.raises(ValueError, utils.spw_range_from_redshifts, obj, z_range=(5., 8.)) # lower bound pytest.raises(ValueError, utils.spw_range_from_redshifts, obj, @@ -147,7 +142,7 @@ def test_spw_range_from_redshifts(self): # Check that valid frequency ranges are returned # with and without future array shapes z_list = [(6.5, 7.5), (7.5, 8.5), (8.5, 9.5)] - for obj in [self.uvd, self.uvd2]: + for obj in [self.uvd]: spw1 = utils.spw_range_from_redshifts(obj, z_range=(7., 8.)) spw2 = utils.spw_range_from_redshifts(obj, z_range=z_list) spw3 = utils.spw_range_from_redshifts(obj, z_range=(12., 14.), @@ -167,126 +162,123 @@ def test_spw_range_from_redshifts(self): def test_calc_blpair_reds(self): fname = os.path.join(DATA_PATH, 'zen.all.xx.LST.1.06964.uvA') - for future_array_shapes in [True, False]: - uvd = UVData() - uvd.read_miriad(fname, use_future_array_shapes=future_array_shapes) - - # basic execution - (bls1, bls2, blps, xants1, xants2, rgrps, lens, - angs) = utils.calc_blpair_reds(uvd, uvd, filter_blpairs=True, extra_info=True, - exclude_auto_bls=False, exclude_permutations=True) - assert len(bls1) == len(bls2) == 15 - assert blps == list(zip(bls1, bls2)) - assert xants1 == xants2 - assert len(xants1) == 42 - assert len(rgrps) == len(bls1) # assert rgrps matches bls1 shape - assert np.max(rgrps) == len(lens) - 1 # assert rgrps indexes lens / angs - - # test xant_flag_thresh - (bls1, bls2, blps, xants1, - xants2) = utils.calc_blpair_reds(uvd, uvd, filter_blpairs=True, exclude_auto_bls=True, exclude_permutations=True, - xant_flag_thresh=0.0) - assert len(bls1) == len(bls2) == 0 - - # test bl_len_range - (bls1, bls2, blps, xants1, - xants2) = utils.calc_blpair_reds(uvd, uvd, filter_blpairs=True, exclude_auto_bls=False, exclude_permutations=True, - bl_len_range=(0, 15.0)) - assert len(bls1) == len(bls2) == 12 - (bls1, bls2, blps, xants1, - xants2) = utils.calc_blpair_reds(uvd, uvd, filter_blpairs=True, exclude_auto_bls=True, exclude_permutations=True, - bl_len_range=(0, 15.0)) - assert len(bls1) == len(bls2) == 5 - assert np.all([bls1[i] != bls2[i] for i in range(len(blps))]) - - # test grouping - (bls1, bls2, blps, xants1, - xants2) = utils.calc_blpair_reds(uvd, uvd, filter_blpairs=True, exclude_auto_bls=False, exclude_permutations=True, - Nblps_per_group=2) - assert len(blps) == 10 - assert isinstance(blps[0], list) - assert blps[0] == [((24, 37), (25, 38)), ((24, 37), (24, 37))] - - # test baseline select on uvd - uvd2 = copy.deepcopy(uvd) - uvd2.select(bls=[(24, 25), (37, 38), (24, 39)]) - (bls1, bls2, blps, xants1, - xants2) = utils.calc_blpair_reds(uvd2, uvd2, filter_blpairs=True, exclude_auto_bls=True, exclude_permutations=True, - bl_len_range=(10.0, 20.0)) - assert blps == [((24, 25), (37, 38))] - - # test exclude_cross_bls - (bls1, bls2, blps, xants1, - xants2) = utils.calc_blpair_reds(uvd, uvd, filter_blpairs=True, exclude_cross_bls=True) - for bl1, bl2 in blps: - assert bl1 == bl2 - - # test exceptions - uvd2 = copy.deepcopy(uvd) - uvd2.antenna_positions[0] += 2 - pytest.raises(AssertionError, utils.calc_blpair_reds, uvd, uvd2) - pytest.raises(AssertionError, utils.calc_blpair_reds, uvd, uvd, exclude_auto_bls=True, exclude_cross_bls=True) + uvd = UVData() + uvd.read_miriad(fname, use_future_array_shapes=True) + + # basic execution + (bls1, bls2, blps, xants1, xants2, rgrps, lens, + angs) = utils.calc_blpair_reds(uvd, uvd, filter_blpairs=True, extra_info=True, + exclude_auto_bls=False, exclude_permutations=True) + assert len(bls1) == len(bls2) == 15 + assert blps == list(zip(bls1, bls2)) + assert xants1 == xants2 + assert len(xants1) == 42 + assert len(rgrps) == len(bls1) # assert rgrps matches bls1 shape + assert np.max(rgrps) == len(lens) - 1 # assert rgrps indexes lens / angs + + # test xant_flag_thresh + (bls1, bls2, blps, xants1, + xants2) = utils.calc_blpair_reds(uvd, uvd, filter_blpairs=True, exclude_auto_bls=True, exclude_permutations=True, + xant_flag_thresh=0.0) + assert len(bls1) == len(bls2) == 0 + + # test bl_len_range + (bls1, bls2, blps, xants1, + xants2) = utils.calc_blpair_reds(uvd, uvd, filter_blpairs=True, exclude_auto_bls=False, exclude_permutations=True, + bl_len_range=(0, 15.0)) + assert len(bls1) == len(bls2) == 12 + (bls1, bls2, blps, xants1, + xants2) = utils.calc_blpair_reds(uvd, uvd, filter_blpairs=True, exclude_auto_bls=True, exclude_permutations=True, + bl_len_range=(0, 15.0)) + assert len(bls1) == len(bls2) == 5 + assert np.all([bls1[i] != bls2[i] for i in range(len(blps))]) + + # test grouping + (bls1, bls2, blps, xants1, + xants2) = utils.calc_blpair_reds(uvd, uvd, filter_blpairs=True, exclude_auto_bls=False, exclude_permutations=True, + Nblps_per_group=2) + assert len(blps) == 10 + assert isinstance(blps[0], list) + assert blps[0] == [((24, 37), (25, 38)), ((24, 37), (24, 37))] + + # test baseline select on uvd + uvd2 = copy.deepcopy(uvd) + uvd2.select(bls=[(24, 25), (37, 38), (24, 39)]) + (bls1, bls2, blps, xants1, + xants2) = utils.calc_blpair_reds(uvd2, uvd2, filter_blpairs=True, exclude_auto_bls=True, exclude_permutations=True, + bl_len_range=(10.0, 20.0)) + assert blps == [((24, 25), (37, 38))] + + # test exclude_cross_bls + (bls1, bls2, blps, xants1, + xants2) = utils.calc_blpair_reds(uvd, uvd, filter_blpairs=True, exclude_cross_bls=True) + for bl1, bl2 in blps: + assert bl1 == bl2 + + # test exceptions + uvd2 = copy.deepcopy(uvd) + uvd2.antenna_positions[0] += 2 + pytest.raises(AssertionError, utils.calc_blpair_reds, uvd, uvd2) + pytest.raises(AssertionError, utils.calc_blpair_reds, uvd, uvd, exclude_auto_bls=True, exclude_cross_bls=True) def test_calc_blpair_reds_autos_only(self): # test include_crosscorrs selection option being set to false. fname = os.path.join(DATA_PATH, 'zen.all.xx.LST.1.06964.uvA') - for future_array_shapes in [True, False]: - uvd = UVData() - uvd.read_miriad(fname, use_future_array_shapes=future_array_shapes) - # basic execution - (bls1, bls2, blps, xants1, xants2, rgrps, lens, - angs) = utils.calc_blpair_reds(uvd, uvd, filter_blpairs=True, extra_info=True, - exclude_auto_bls=False, exclude_permutations=True, include_crosscorrs=False, - include_autocorrs=True) - assert len(bls1) > 0 - for bl1, bl2 in zip(bls1, bls2): - assert bl1[0] == bl1[1] - assert bl2[0] == bl2[1] + uvd = UVData() + uvd.read_miriad(fname, use_future_array_shapes=True) + # basic execution + (bls1, bls2, blps, xants1, xants2, rgrps, lens, + angs) = utils.calc_blpair_reds(uvd, uvd, filter_blpairs=True, extra_info=True, + exclude_auto_bls=False, exclude_permutations=True, include_crosscorrs=False, + include_autocorrs=True) + assert len(bls1) > 0 + for bl1, bl2 in zip(bls1, bls2): + assert bl1[0] == bl1[1] + assert bl2[0] == bl2[1] def test_get_delays(self): utils.get_delays(np.linspace(100., 200., 50)*1e6) def test_get_reds(self): fname = os.path.join(DATA_PATH, 'zen.all.xx.LST.1.06964.uvA') - for future_array_shapes in [True, False]: - uvd = UVData() - uvd.read_miriad(fname, read_data=False, use_future_array_shapes=future_array_shapes) - antpos, ants = uvd.get_ENU_antpos() - antpos_d = dict(list(zip(ants, antpos))) - - # test basic execution - xants = [0, 1, 2] - r, l, a = utils.get_reds(fname, xants=xants) - assert np.all([np.all([bl[0] not in xants and bl[1] not in xants for bl in _r]) for _r in r]) - assert len(r) == len(a) == len(l) - assert len(r) == 104 - - r2, l2, a2 = utils.get_reds(uvd, xants=xants) - _ = [np.testing.assert_array_equal(_r1, _r2) for _r1, _r2 in zip(r, r2)] - - r2, l2, a2 = utils.get_reds(antpos_d, xants=xants) - _ = [np.testing.assert_array_equal(_r1, _r2) for _r1, _r2 in zip(r, r2)] - - # restrict - bl_len_range = (14, 16) - bl_deg_range = (55, 65) - r, l, a = utils.get_reds(uvd, bl_len_range=bl_len_range, bl_deg_range=bl_deg_range) - assert (np.all([_l > bl_len_range[0] and _l < bl_len_range[1] for _l in l])) - assert (np.all([_a > bl_deg_range[0] and _a < bl_deg_range[1] for _a in a])) - - # min EW cut - r, l, a = utils.get_reds(uvd, bl_len_range=(14, 16), min_EW_cut=14) - assert len(l) == len(a) == 1 - assert np.isclose(a[0] % 180, 0, atol=1) - - # autos - r, l, a = utils.get_reds(fname, xants=xants, add_autos=True) - np.testing.assert_almost_equal(l[0], 0) - np.testing.assert_almost_equal(a[0], 0) - assert len(r) == 105 - - # Check errors when wrong types input - pytest.raises(TypeError, utils.get_reds, [1., 2.]) + uvd = UVData() + uvd.read_miriad(fname, read_data=False, use_future_array_shapes=True) + antpos, ants = uvd.get_ENU_antpos() + antpos_d = dict(list(zip(ants, antpos))) + + # test basic execution + xants = [0, 1, 2] + r, l, a = utils.get_reds(fname, xants=xants) + assert np.all([np.all([bl[0] not in xants and bl[1] not in xants for bl in _r]) for _r in r]) + assert len(r) == len(a) == len(l) + assert len(r) == 104 + + r2, l2, a2 = utils.get_reds(uvd, xants=xants) + _ = [np.testing.assert_array_equal(_r1, _r2) for _r1, _r2 in zip(r, r2)] + + r2, l2, a2 = utils.get_reds(antpos_d, xants=xants) + _ = [np.testing.assert_array_equal(_r1, _r2) for _r1, _r2 in zip(r, r2)] + + # restrict + bl_len_range = (14, 16) + bl_deg_range = (55, 65) + r, l, a = utils.get_reds(uvd, bl_len_range=bl_len_range, bl_deg_range=bl_deg_range) + assert (np.all([_l > bl_len_range[0] and _l < bl_len_range[1] for _l in l])) + assert (np.all([_a > bl_deg_range[0] and _a < bl_deg_range[1] for _a in a])) + + # min EW cut + r, l, a = utils.get_reds(uvd, bl_len_range=(14, 16), min_EW_cut=14) + assert len(l) == len(a) == 1 + assert np.isclose(a[0] % 180, 0, atol=1) + + # autos + r, l, a = utils.get_reds(fname, xants=xants, add_autos=True) + np.testing.assert_almost_equal(l[0], 0) + np.testing.assert_almost_equal(a[0], 0) + assert len(r) == 105 + + # Check errors when wrong types input + pytest.raises(TypeError, utils.get_reds, [1., 2.]) def test_get_reds_autos_only(self): fname = os.path.join(DATA_PATH, 'zen.all.xx.LST.1.06964.uvA') @@ -333,22 +325,21 @@ def test_uvd_to_Tsys(self): # PROPER USAGE # check different ways to call method work and are equivalent # with or without future array shapes - for obj in [self.uvd, self.uvd2]: - tsys_estimate = utils.uvd_to_Tsys(obj, self.beam) - tsys_estimate2 = utils.uvd_to_Tsys( - obj, - os.path.join(DATA_PATH, 'HERA_NF_dipole_power.beamfits') - ) - assert np.allclose(tsys_estimate.data_array, tsys_estimate2.data_array) - uvp2, _ = testing.build_vanilla_uvpspec(beam=self.beam) - tsys_estimate3 = utils.uvd_to_Tsys(obj, uvp2) - assert np.allclose(tsys_estimate.data_array, tsys_estimate3.data_array) - - # CHECK ERROR CALLS - # uvp called for beam has no beam information - pytest.raises(ValueError, utils.uvd_to_Tsys, obj, self.uvp) - # beam has wrong format - pytest.raises(ValueError, utils.uvd_to_Tsys, obj, 12.) + tsys_estimate = utils.uvd_to_Tsys(self.uvd, self.beam) + tsys_estimate2 = utils.uvd_to_Tsys( + self.uvd, + os.path.join(DATA_PATH, 'HERA_NF_dipole_power.beamfits') + ) + assert np.allclose(tsys_estimate.data_array, tsys_estimate2.data_array) + uvp2, _ = testing.build_vanilla_uvpspec(beam=self.beam) + tsys_estimate3 = utils.uvd_to_Tsys(self.uvd, uvp2) + assert np.allclose(tsys_estimate.data_array, tsys_estimate3.data_array) + + # CHECK ERROR CALLS + # uvp called for beam has no beam information + pytest.raises(ValueError, utils.uvd_to_Tsys, self.uvd, self.uvp) + # beam has wrong format + pytest.raises(ValueError, utils.uvd_to_Tsys, self.uvd, 12.) diff --git a/hera_pspec/tests/test_uvpspec.py b/hera_pspec/tests/test_uvpspec.py index ce418a73..109c88df 100644 --- a/hera_pspec/tests/test_uvpspec.py +++ b/hera_pspec/tests/test_uvpspec.py @@ -606,13 +606,11 @@ def test_get_exact_window_functions(self): ## tests - # give Gaussian beam as input - widths = -0.0343 * self.uvp_wf.freq_array/1e6 + 11.30 - gaussian_beam = uvwindow.FTBeam.gaussian(freq_array=self.uvp_wf.freq_array, - widths=widths, - pol='xx') - uvp.get_exact_window_functions(ftbeam=gaussian_beam, - spw_array=0, inplace=True, verbose=True) + # output exact_window functions but does not make them attributes + kperp_bins, kpara_bins, wf_array = uvp.get_exact_window_functions(ftbeam=self.ft_file, + inplace=False) + # check if result is the same with and without inplace + assert np.allclose(wf_array[0], uvp.window_function_array[0]) # obtain exact window functions for one spw only uvp.get_exact_window_functions(ftbeam=self.ft_file, spw_array=0, inplace=True, verbose=True) @@ -621,11 +619,13 @@ def test_get_exact_window_functions(self): ftbeam=self.ft_file, spw_array=2, inplace=True) - # output exact_window functions but does not make them attributes - kperp_bins, kpara_bins, wf_array = uvp.get_exact_window_functions(ftbeam=self.ft_file, - inplace=False) - # check if result is the same with and without inplace - assert np.all(wf_array[0]==uvp.window_function_array[0]) + # give Gaussian beam as input + widths = -0.0343 * self.uvp_wf.freq_array/1e6 + 11.30 + gaussian_beam = uvwindow.FTBeam.gaussian(freq_array=self.uvp_wf.freq_array, + widths=widths, + pol='xx') + uvp.get_exact_window_functions(ftbeam=gaussian_beam, + spw_array=0, inplace=True, verbose=True) def test_fold_spectra(self): uvp = copy.deepcopy(self.uvp) diff --git a/hera_pspec/tests/test_uvwindow.py b/hera_pspec/tests/test_uvwindow.py index 6d5f7240..84b98f88 100644 --- a/hera_pspec/tests/test_uvwindow.py +++ b/hera_pspec/tests/test_uvwindow.py @@ -58,7 +58,7 @@ def test_init(self): verbose=self.verbose, x_orientation=self.x_orientation) assert self.pol == test.pol - assert np.all(self.data == test.ft_beam) + assert np.allclose(self.data, test.ft_beam) # raise assertion error if data is not dim 3 pytest.raises(AssertionError, uvwindow.FTBeam, @@ -111,11 +111,11 @@ def test_from_file(self): # tests related to spw_range test1 = uvwindow.FTBeam.from_file(ftfile=self.ft_file, spw_range=self.spw_range) - assert np.all(test1.freq_array == self.freq_array) + assert np.allclose(test1.freq_array, self.freq_array) test2 = uvwindow.FTBeam.from_file(ftfile=self.ft_file, spw_range=None) - assert np.all(test2.freq_array == self.bandwidth) + assert np.allclose(test2.freq_array, self.bandwidth) pytest.raises(AssertionError, uvwindow.FTBeam.from_file, spw_range=(13), ftfile=self.ft_file) @@ -452,13 +452,13 @@ def test_get_cylindrical_wf(self): kpara_bins=None, return_bins='unweighted') # check normalisation - assert np.all(np.isclose(np.sum(cyl_wf, axis=(1, 2)), 1., atol=1e-3)) + assert np.allclose(np.sum(cyl_wf, axis=(1, 2)), 1., atol=1e-3) assert kperp.size == cyl_wf.shape[1] assert kpara.size == cyl_wf.shape[2] assert test.Nfreqs == cyl_wf.shape[0] # test the bins are recovered by get_kperp_bins and get_kpara_bins - assert np.all(kperp == test.get_kperp_bins(bl_len).value) - assert np.all(kpara == test.get_kpara_bins(test.freq_array).value) + assert np.allclose(kperp, test.get_kperp_bins(bl_len).value) + assert np.allclose(kpara, test.get_kpara_bins(test.freq_array).value) # test different key words @@ -467,8 +467,8 @@ def test_get_cylindrical_wf(self): kperp_bins=kperp*test.kunits, kpara_bins=None, return_bins='unweighted') - assert np.all(cyl_wf2 == cyl_wf) - assert np.all(kperp2 == kperp) # unweighted option to return_bins + assert np.allclose(cyl_wf2, cyl_wf) + assert np.allclose(kperp2, kperp) # unweighted option to return_bins # ValueError raised if kperp_bins not linearly spaced kperp_log = np.logspace(-2, 0, 100) pytest.raises(ValueError, test.get_cylindrical_wf, bl_len, @@ -480,8 +480,8 @@ def test_get_cylindrical_wf(self): kperp_bins=None, kpara_bins=kpara*test.kunits, return_bins='unweighted') - assert np.all(cyl_wf3 == cyl_wf) - assert np.all(kpara == kpara3) + assert np.allclose(cyl_wf3, cyl_wf) + assert np.allclose(kpara, kpara3) # ValueError raised if kpara_bins not linearly spaced kpara_log = np.logspace(-1, 1, 100) pytest.raises(ValueError, test.get_cylindrical_wf, bl_len, diff --git a/hera_pspec/uvpspec.py b/hera_pspec/uvpspec.py index acf6def8..b1313c5e 100644 --- a/hera_pspec/uvpspec.py +++ b/hera_pspec/uvpspec.py @@ -1291,9 +1291,13 @@ def get_ENU_bl_vecs(self): blvecs : array_like Array of baseline vectors. """ + latitude, longitude, altitude = uvutils.LatLonAlt_from_XYZ(self.telescope_location[None]) return uvutils.ENU_from_ECEF( - (self.bl_vecs + self.telescope_location), \ - *uvutils.LatLonAlt_from_XYZ(self.telescope_location[None])) + self.bl_vecs + self.telescope_location, + latitude=latitude, + longitude=longitude, + altitude=altitude + ) def read_from_group(self, grp, just_meta=False, spws=None, bls=None, @@ -1510,7 +1514,7 @@ def write_to_group(self, group, run_check=True): # Do Unicode/string conversions, as HDF5 struggles with them if k == 'labels': - this_attr = [np.string_(lbl) for lbl in this_attr] + this_attr = [np.bytes_(lbl) for lbl in this_attr] # Store attribute in group group.attrs[k] = this_attr diff --git a/hera_pspec/uvwindow.py b/hera_pspec/uvwindow.py index faf2ce4c..364ccaf5 100644 --- a/hera_pspec/uvwindow.py +++ b/hera_pspec/uvwindow.py @@ -7,7 +7,7 @@ import os import copy import time -from scipy.interpolate import interp2d +from scipy.interpolate import RegularGridInterpolator from pyuvdata import UVBeam, UVData from astropy import units from pathlib import Path @@ -569,9 +569,10 @@ def _interpolate_ft_beam(self, bl_len, ft_beam): # kperp values over frequency array for bl_len k = self._kperp4bl_freq(self.freq_array[i], bl_len, ngrid=ngrid) # interpolate the FT beam over values onto regular kgrid - A_real = interp2d(k, k, ft_beam[i, :, :], bounds_error=False, - fill_value=0.) - interp_ft_beam[:, :, i] = A_real(kgrid, kgrid) + A_real = RegularGridInterpolator( + (k,k), ft_beam[i], bounds_error=False, fill_value=0.0 + ) + interp_ft_beam[:, :, i] = A_real(np.array([kgrid, kgrid]).T) return interp_ft_beam, kperp_norm