Skip to content

Commit

Permalink
Merge pull request #173 from HERA-Team/combine_wgts_dimensions
Browse files Browse the repository at this point in the history 
Fix dimensionality of wgt_array in combine_uvpspec
  • Loading branch information
@nkern
nkern committed Sep 14, 2018
2 parents 4093bb9 + ee35efc commit b129d47
Show file tree
Hide file tree
Showing 4 changed files with 93 additions and 47 deletions.
4 changes: 2 additions & 2 deletions hera_pspec/container.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -406,10 +406,10 @@ def combine_psc_spectra(psc, groups=None, dset_split_str='_x_', ext_split_str='_
for uvp in to_merge:
if uvp != spc:
del psc.data[grp][uvp]
except:
except Exception as exc:
# merge failed, so continue
if verbose:
print "uvp merge failed for spectra {}/{}".format(grp, spc)
print "uvp merge failed for spectra {}/{}, exception: {}".format(grp, spc, exc)


def get_combine_psc_spectra_argparser():
Expand Down
84 changes: 53 additions & 31 deletions hera_pspec/testing.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,8 @@

def build_vanilla_uvpspec(beam=None):
"""
Build an example vanilla UVPSpec object from scratch, with all necessary metadata.
Build an example vanilla UVPSpec object from scratch, with all necessary
metadata.
Parameters
----------
Expand Down Expand Up @@ -54,7 +55,8 @@ def build_vanilla_uvpspec(beam=None):
spw_freq_array = np.tile(np.arange(Nspws), Nfreqs)
spw_dly_array = np.tile(np.arange(Nspws), Ndlys)
spw_array = np.arange(Nspws)
freq_array = np.repeat(np.linspace(100e6, 105e6, Nfreqs, endpoint=False), Nspws)
freq_array = np.repeat(np.linspace(100e6, 105e6, Nfreqs, endpoint=False),
Nspws)
dly_array = np.repeat(utils.get_delays(freq_array, n_dlys=Ndlys), Nspws)
pol_array = np.array([-5])
Npols = len(pol_array)
Expand Down Expand Up @@ -84,26 +86,31 @@ def build_vanilla_uvpspec(beam=None):
store_cov=True
cosmo = conversions.Cosmo_Conversions()

data_array, wgt_array, integration_array, nsample_array, cov_array = {}, {}, {}, {}, {}
data_array, wgt_array = {}, {}
integration_array, nsample_array, cov_array = {}, {}, {}
for s in spw_array:
data_array[s] = np.ones((Nblpairts, Ndlys, Npols), dtype=np.complex) \
* blpair_array[:, None, None] / 1e9
wgt_array[s] = np.ones((Nblpairts, Ndlys, 2, Npols), dtype=np.float)
wgt_array[s] = np.ones((Nblpairts, Nfreqs, 2, Npols), dtype=np.float)
# NB: The wgt_array has dimensions Nfreqs rather than Ndlys; it has the
# dimensions of the input visibilities, not the output delay spectra
integration_array[s] = np.ones((Nblpairts, Npols), dtype=np.float)
nsample_array[s] = np.ones((Nblpairts, Npols), dtype=np.float)
cov_array[s] =np.moveaxis(np.array([[np.identity(Ndlys,dtype=np.complex)\
for m in range(Nblpairts)] for n in range(Npols)]),0,-1)
for m in range(Nblpairts)]
for n in range(Npols)]), 0, -1)


params = ['Ntimes', 'Nfreqs', 'Nspws', 'Nspwdlys', 'Nspwfreqs', 'Nspws', 'Nblpairs', 'Nblpairts',
'Npols', 'Ndlys', 'Nbls', 'blpair_array', 'time_1_array',
'time_2_array', 'lst_1_array', 'lst_2_array', 'spw_array',
params = ['Ntimes', 'Nfreqs', 'Nspws', 'Nspwdlys', 'Nspwfreqs', 'Nspws',
'Nblpairs', 'Nblpairts', 'Npols', 'Ndlys', 'Nbls',
'blpair_array', 'time_1_array', 'time_2_array',
'lst_1_array', 'lst_2_array', 'spw_array',
'dly_array', 'freq_array', 'pol_array', 'data_array', 'wgt_array',
'integration_array', 'bl_array', 'bl_vecs', 'telescope_location',
'vis_units', 'channel_width', 'weighting', 'history', 'taper', 'norm',
'git_hash', 'nsample_array', 'time_avg_array', 'lst_avg_array',
'cosmo', 'scalar_array', 'labels', 'norm_units', 'labels', 'label_1_array',
'label_2_array','store_cov','cov_array', 'spw_dly_array', 'spw_freq_array']
'vis_units', 'channel_width', 'weighting', 'history', 'taper',
'norm', 'git_hash', 'nsample_array', 'time_avg_array',
'lst_avg_array', 'cosmo', 'scalar_array', 'labels', 'norm_units',
'labels', 'label_1_array', 'label_2_array', 'store_cov',
'cov_array', 'spw_dly_array', 'spw_freq_array']

if beam is not None:
params += ['OmegaP', 'OmegaPP', 'beam_freqs']
Expand All @@ -117,7 +124,9 @@ def build_vanilla_uvpspec(beam=None):
return uvp, cosmo


def uvpspec_from_data(data, bl_grps, data_std=None, spw_ranges=None, beam=None, taper='none', cosmo=None, verbose=False):
def uvpspec_from_data(data, bl_grps, data_std=None, spw_ranges=None,
beam=None, taper='none', cosmo=None, n_dlys=None,
verbose=False):
"""
Build an example UVPSpec object from a visibility file and PSpecData.
Expand All @@ -127,26 +136,34 @@ def uvpspec_from_data(data, bl_grps, data_std=None, spw_ranges=None, beam=None,
This can be a UVData object or a string filepath to a miriad file.
bl_grps : list
This is a list of baseline groups (e.g. redundant groups) to form blpairs from.
This is a list of baseline groups (e.g. redundant groups) to form
blpairs from.
Ex: [[(24, 25), (37, 38), ...], [(24, 26), (37, 39), ...], ... ]
data_std: UVData object or str or None
Can be UVData object or a string filepath to a miriad file.
data_std: UVData object or str, optional
Can be UVData object or a string filepath to a miriad file.
Default: None.
spw_ranges : list
List of spectral window tuples. See PSpecData.pspec docstring for details.
spw_ranges : list, optional
List of spectral window tuples. See PSpecData.pspec docstring for
details. Default: None.
beam : PSpecBeamBase subclass or str
beam : PSpecBeamBase subclass or str, optional
This can be a subclass of PSpecBeamBase of a string filepath to a
UVBeam healpix map.
UVBeam healpix map. Default: None.
taper : string
Optional tapering applied to the data before OQE.
taper : str, optional
Optional tapering applied to the data before OQE. Default: 'none'.
cosmo : Cosmo_Conversions object
Cosmology object.
n_dlys : int, optional
Number of delay bins to use. Default: None (uses as many delay bins as
frequency channels).
verbose : bool
if True, report feedback to standard output
verbose : bool, optional
if True, report feedback to standard output. Default: False.
Returns
-------
Expand Down Expand Up @@ -181,20 +198,25 @@ def uvpspec_from_data(data, bl_grps, data_std=None, spw_ranges=None, beam=None,
beam.cosmo = cosmo

# instantiate pspecdata
ds = pspecdata.PSpecData(dsets=[uvd, uvd], dsets_std=[uvd_std, uvd_std], wgts=[None, None], labels=['d1', 'd2'], beam=beam)
ds = pspecdata.PSpecData(dsets=[uvd, uvd], dsets_std=[uvd_std, uvd_std],
wgts=[None, None], labels=['d1', 'd2'], beam=beam)

# get blpair groups
assert isinstance(bl_grps, list), "bl_grps must be a list"
if not isinstance(bl_grps[0], list): bl_grps = [bl_grps]
assert np.all([isinstance(blgrp, list) for blgrp in bl_grps]), "bl_grps must be fed as a list of lists"
assert np.all([isinstance(blgrp[0], tuple) for blgrp in bl_grps]), "bl_grps must be fed as a list of lists of tuples"
assert np.all([isinstance(blgrp, list) for blgrp in bl_grps]), \
"bl_grps must be fed as a list of lists"
assert np.all([isinstance(blgrp[0], tuple) for blgrp in bl_grps]), \
"bl_grps must be fed as a list of lists of tuples"
bls1, bls2 = [], []
for blgrp in bl_grps:
_bls1, _bls2, _ = utils.construct_blpairs(blgrp, exclude_auto_bls=True, exclude_permutations=True)
_bls1, _bls2, _ = utils.construct_blpairs(blgrp, exclude_auto_bls=True,
exclude_permutations=True)
bls1.extend(_bls1)
bls2.extend(_bls2)

# run pspec
uvp = ds.pspec(bls1, bls2, (0, 1), (pol, pol), input_data_weight='identity', spw_ranges=spw_ranges,
taper=taper, verbose=verbose,store_cov=store_cov)
uvp = ds.pspec(bls1, bls2, (0, 1), (pol, pol), input_data_weight='identity',
spw_ranges=spw_ranges, taper=taper, verbose=verbose,
store_cov=store_cov, n_dlys=n_dlys)
return uvp
46 changes: 34 additions & 12 deletions hera_pspec/tests/test_uvpspec.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ def test_get_funcs(self):
nt.assert_almost_equal(d[0,0], (101.1021011020000001+0j))
# get_wgts
w = self.uvp.get_wgts((0, ((1, 2), (1, 2)), 'xx'))
nt.assert_equal(w.shape, (10, 30, 2))
nt.assert_equal(w.shape, (10, 50, 2)) # should have Nfreq dim, not Ndlys
nt.assert_true(w.dtype == np.float)
nt.assert_equal(w[0,0,0], 1.0)
# get_integrations
Expand Down Expand Up @@ -504,16 +504,29 @@ def test_combine_uvpspec(self):
uvp3.pol_array[0] = -7
out = uvp1 + uvp2 + uvp3
nt.assert_equal(out.Npols, 3)

# Test whether n_dlys != Nfreqs works
uvp4 = testing.uvpspec_from_data(uvd, bls, beam=beam,
spw_ranges=[(20, 30), (60, 90)],
n_dlys=[5, 15])
uvp4b = copy.deepcopy(uvp4)
uvp4b.pol_array[0] = -6
out = uvpspec.combine_uvpspec([uvp4, uvp4b], verbose=False)


def test_combine_uvpspec_std(self):
# setup uvp build
uvd = UVData()
uvd_std = UVData()
uvd.read_miriad(os.path.join(DATA_PATH, 'zen.even.xx.LST.1.28828.uvOCRSA'))
uvd_std.read_miriad(os.path.join(DATA_PATH,'zen.even.xx.LST.1.28828.uvOCRSA'))
beam = pspecbeam.PSpecBeamUV(os.path.join(DATA_PATH, "HERA_NF_dipole_power.beamfits"))
uvd_std.read_miriad(
os.path.join(DATA_PATH,'zen.even.xx.LST.1.28828.uvOCRSA'))
beam = pspecbeam.PSpecBeamUV(
os.path.join(DATA_PATH, "HERA_NF_dipole_power.beamfits"))
bls = [(37, 38), (38, 39), (52, 53)]
uvp1 = testing.uvpspec_from_data(uvd, bls, data_std=uvd_std, spw_ranges=[(20, 24), (64, 68)], beam=beam)
uvp1 = testing.uvpspec_from_data(uvd, bls, data_std=uvd_std,
spw_ranges=[(20, 24), (64, 68)],
beam=beam)
# test failure due to overlapping data
uvp2 = copy.deepcopy(uvp1)
nt.assert_raises(AssertionError, uvpspec.combine_uvpspec, [uvp1, uvp2])
Expand All @@ -527,24 +540,32 @@ def test_combine_uvpspec_std(self):
nt.assert_raises(AssertionError, uvpspec.combine_uvpspec, [uvp1, uvp2])

# test partial data overlap failure
uvp2 = testing.uvpspec_from_data(uvd, [(37, 38), (38, 39), (53, 54)], data_std=uvd_std, spw_ranges=[(20, 24), (64, 68)], beam=beam)
uvp2 = testing.uvpspec_from_data(uvd, [(37, 38), (38, 39), (53, 54)],
data_std=uvd_std,
spw_ranges=[(20, 24), (64, 68)],
beam=beam)
nt.assert_raises(AssertionError, uvpspec.combine_uvpspec, [uvp1, uvp2])
uvp2 = testing.uvpspec_from_data(uvd, bls, spw_ranges=[(20, 24), (64, 68)], data_std=uvd_std, beam=beam)
uvp2 = testing.uvpspec_from_data(uvd, bls,
spw_ranges=[(20, 24), (64, 68)],
data_std=uvd_std, beam=beam)
nt.assert_raises(AssertionError, uvpspec.combine_uvpspec, [uvp1, uvp2])
uvp2 = copy.deepcopy(uvp1)
uvp2.pol_array[0] = -6
uvp2 = uvpspec.combine_uvpspec([uvp1, uvp2], verbose=False)
nt.assert_raises(AssertionError, uvpspec.combine_uvpspec, [uvp1, uvp2])

# test concat across spw
uvp2 = testing.uvpspec_from_data(uvd, bls, spw_ranges=[(85, 91)], data_std=uvd_std, beam=beam)
uvp2 = testing.uvpspec_from_data(uvd, bls, spw_ranges=[(85, 91)],
data_std=uvd_std, beam=beam)
out = uvpspec.combine_uvpspec([uvp1, uvp2], verbose=False)
nt.assert_equal(out.Nspws, 3)
nt.assert_equal(out.Nfreqs, 14)
nt.assert_equal(out.Nspwdlys, 14)

# test concat across blpairts
uvp2 = testing.uvpspec_from_data(uvd, [(53, 54), (67, 68)], spw_ranges=[(20, 24), (64, 68)], data_std=uvd_std, beam=beam)
uvp2 = testing.uvpspec_from_data(uvd, [(53, 54), (67, 68)],
spw_ranges=[(20, 24), (64, 68)],
data_std=uvd_std, beam=beam)
out = uvpspec.combine_uvpspec([uvp1, uvp2], verbose=False)
nt.assert_equal(out.Nblpairs, 4)
nt.assert_equal(out.Nbls, 5)
Expand All @@ -558,17 +579,17 @@ def test_combine_uvpspec_std(self):
nt.assert_raises(AssertionError, uvpspec.combine_uvpspec, [uvp1, uvp2])

# test feed as strings
uvp1 = testing.uvpspec_from_data(uvd, bls, spw_ranges=[(20, 30)], data_std=uvd_std, beam=beam)
uvp1 = testing.uvpspec_from_data(uvd, bls, spw_ranges=[(20, 30)],
data_std=uvd_std, beam=beam)
uvp2 = copy.deepcopy(uvp1)
uvp2.pol_array[0] = -6
uvp1.write_hdf5('uvp1.hdf5', overwrite=True)
uvp2.write_hdf5('uvp2.hdf5', overwrite=True)
out = uvpspec.combine_uvpspec(['uvp1.hdf5', 'uvp2.hdf5'], verbose=False)
nt.assert_true(out.Npols, 2)
for ff in ['uvp1.hdf5', 'uvp2.hdf5']:
if os.path.exists(ff):
os.remove(ff)

if os.path.exists(ff): os.remove(ff)

# test UVPSpec __add__
uvp3 = copy.deepcopy(uvp1)
uvp3.pol_array[0] = -7
Expand All @@ -591,3 +612,4 @@ def test_conj_blpair():
blpair = uvputils._conj_blpair(101102103104, which='both')
nt.assert_equal(blpair, 102101104103)
nt.assert_raises(ValueError, uvputils._conj_blpair, 102101103104, which='foo')

6 changes: 4 additions & 2 deletions hera_pspec/uvpspec.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1820,7 +1820,9 @@ def combine_uvpspec(uvps, verbose=True):
# Initialize new arrays
u.data_array[i] = np.empty((Nblpairts, spw[3], Npols), np.complex128)
u.integration_array[i] = np.empty((Nblpairts, Npols), np.float64)
u.wgt_array[i] = np.empty((Nblpairts, spw[3], 2, Npols), np.float64)
u.wgt_array[i] = np.empty((Nblpairts, spw[2], 2, Npols), np.float64)
# spw[2] == Nfreqs (wgt_array is not resampled if Ndlys != Nfreqs,
# so needs to keep this shape)
u.nsample_array[i] = np.empty((Nblpairts, Npols), np.float64)
if store_cov:
u.cov_array[i] = np.empty((Nblpairts, spw[3], spw[3], Npols),
Expand Down Expand Up @@ -2074,7 +2076,7 @@ def get_uvp_overlap(uvps, just_meta=True, verbose=True):
unique_spws : list
List of unique spectral window tuples (spw_freq_start, spw_freq_end,
spw_Nfreqs) across all input uvps
spw_Nfreqs, spw_Ndlys) across all input uvps
unique_blpts : list
List of unique baseline-pair-time tuples (blpair_integer,
Expand Down

0 comments on commit b129d47

Please sign in to comment.