-
Notifications
You must be signed in to change notification settings - Fork 3
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
11 changed files
with
783 additions
and
204 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,168 @@ | ||
| #!/usr/bin/env python2 | ||
| import numpy as np | ||
| import copy, operator, itertools | ||
| from collections import OrderedDict as odict | ||
| from hera_pspec import uvpspec, pspecdata, conversions, pspecbeam | ||
| from pyuvdata import UVData | ||
|
|
||
|
|
||
| def build_vanilla_uvpspec(beam=None): | ||
| """ | ||
| Build an example vanilla UVPSpec object from scratch, with all necessary metadata. | ||
| Parameters | ||
| ---------- | ||
| beam : PSpecBeamBase subclass | ||
| Returns | ||
| ------- | ||
| uvp : UVPSpec object | ||
| """ | ||
| uvp = uvpspec.UVPSpec() | ||
|
|
||
| Ntimes = 10 | ||
| Nfreqs = 50 | ||
| Ndlys = Nfreqs | ||
| Nspws = 1 | ||
| Nspwdlys = Nspws * Nfreqs | ||
|
|
||
| # [((1, 2), (1, 2)), ((2, 3), (2, 3)), ((1, 3), (1, 3))] | ||
| blpairs = [1002001002, 2003002003, 1003001003] | ||
| bls = [1002, 2003, 1003] | ||
| Nbls = len(bls) | ||
| Nblpairs = len(blpairs) | ||
| Nblpairts = Nblpairs * Ntimes | ||
|
|
||
| blpair_array = np.tile(blpairs, Ntimes) | ||
| bl_array = np.array(bls) | ||
| bl_vecs = np.array([[ 5.33391548e+00, -1.35907816e+01, -7.91624188e-09], | ||
| [ -8.67982998e+00, 4.43554478e+00, -1.08695203e+01], | ||
| [ -3.34591450e+00, -9.15523687e+00, -1.08695203e+01]]) | ||
| time_array = np.repeat(np.linspace(2458042.1, 2458042.2, Ntimes), Nblpairs) | ||
| time_1_array = time_array | ||
| time_2_array = time_array | ||
| lst_array = np.repeat(np.ones(Ntimes, dtype=np.float), Nblpairs) | ||
| lst_1_array = lst_array | ||
| lst_2_array = lst_array | ||
| time_avg_array = time_array | ||
| lst_avg_array = lst_array | ||
| spws = np.arange(Nspws) | ||
| spw_array = np.tile(spws, Ndlys) | ||
| freq_array = np.repeat(np.linspace(100e6, 105e6, Nfreqs, endpoint=False), Nspws) | ||
| dly_array = np.fft.fftshift(np.repeat(np.fft.fftfreq(Nfreqs, np.median(np.diff(freq_array))), Nspws)) | ||
| pol_array = np.array([-5]) | ||
| Npols = len(pol_array) | ||
| vis_units = 'unknown' | ||
| norm_units = 'Hz str' | ||
| weighting = 'identity' | ||
| channel_width = np.median(np.diff(freq_array)) | ||
| history = 'example' | ||
| taper = "none" | ||
| norm = "I" | ||
| git_hash = "random" | ||
| scalar_array = np.ones((Nspws, Npols), np.float) | ||
| label1 = 'red' | ||
| label2 = 'blue' | ||
| labels = np.array([label1, label2]) | ||
| label_1_array = np.ones((Nspws, Nblpairts, Npols), np.int) * 0 | ||
| label_2_array = np.ones((Nspws, Nblpairts, Npols), np.int) * 1 | ||
| if beam is not None: | ||
| OmegaP, OmegaPP = beam.get_Omegas(beam.primary_beam.polarization_array[0]) | ||
| beam_freqs = beam.beam_freqs | ||
|
|
||
| # HERA coordinates in Karoo Desert, SA | ||
| telescope_location = np.array([5109325.85521063, | ||
| 2005235.09142983, | ||
| -3239928.42475397]) | ||
|
|
||
| cosmo = conversions.Cosmo_Conversions() | ||
|
|
||
| data_array, wgt_array, integration_array, nsample_array = {}, {}, {}, {} | ||
| for s in spws: | ||
| 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) | ||
| integration_array[s] = np.ones((Nblpairts, Npols), dtype=np.float) | ||
| nsample_array[s] = np.ones((Nblpairts, Npols), dtype=np.float) | ||
|
|
||
| params = ['Ntimes', 'Nfreqs', 'Nspws', 'Nspwdlys', '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'] | ||
|
|
||
| if beam is not None: | ||
| params += ['OmegaP', 'OmegaPP', 'beam_freqs'] | ||
|
|
||
| # Set all parameters | ||
| for p in params: | ||
| setattr(uvp, p, locals()[p]) | ||
|
|
||
| uvp.check() | ||
|
|
||
| return uvp, cosmo | ||
|
|
||
| def uvpspec_from_data(data, bls, spw_ranges=None, beam=None, taper='none', cosmo=None, verbose=False): | ||
| """ | ||
| Build an example UVPSpec object from a visibility file and PSpecData. | ||
| Parameters | ||
| ---------- | ||
| data : UVData object or str | ||
| This can be a UVData object or a string filepath to a miriad file. | ||
| bls : list | ||
| This is a list of at least 2 baseline tuples. | ||
| Ex: [(24, 25), (37, 38), ...] | ||
| spw_ranges : list | ||
| List of spectral window tuples. See PSpecData.pspec docstring for details. | ||
| beam : PSpecBeamBase subclass or str | ||
| This can be a subclass of PSpecBeamBase of a string filepath to a | ||
| UVBeam healpix map. | ||
| taper : string | ||
| Optional tapering applied to the data before OQE. | ||
| cosmo : Cosmo_Conversions object | ||
| verbose : bool | ||
| if True, report feedback to standard output | ||
| Returns | ||
| ------- | ||
| uvp : UVPSpec object | ||
| """ | ||
| # load data | ||
| if isinstance(data, str): | ||
| uvd = UVData() | ||
| uvd.read_miriad(data) | ||
| elif isinstance(data, UVData): | ||
| uvd = data | ||
|
|
||
| # get pol | ||
| pol = uvd.polarization_array[0] | ||
|
|
||
| # load beam | ||
| if isinstance(beam, str): | ||
| beam = pspecbeam.PSpecBeamUV(beam, cosmo=cosmo) | ||
| if beam is not None and cosmo is not None: | ||
| beam.cosmo = cosmo | ||
|
|
||
| # instantiate pspecdata | ||
| ds = pspecdata.PSpecData(dsets=[uvd, uvd], wgts=[None, None], labels=['d1', 'd2'], beam=beam) | ||
|
|
||
| # get red bls | ||
| bls1, bls2, _ = pspecdata.construct_blpairs(bls, exclude_auto_bls=True) | ||
|
|
||
| # run pspec | ||
| uvp = ds.pspec(bls1, bls2, (0, 1), (pol, pol), input_data_weight='identity', spw_ranges=spw_ranges, | ||
| taper=taper, verbose=verbose) | ||
|
|
||
| return uvp | ||
|
|
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,36 @@ | ||
| import nose.tools as nt | ||
| from hera_pspec.data import DATA_PATH | ||
| from hera_pspec import testing, uvpspec, conversions, pspecbeam | ||
| import os | ||
| from pyuvdata import UVData | ||
| import numpy as np | ||
|
|
||
|
|
||
| def test_build_vanilla_uvpspec(): | ||
| uvp, cosmo = testing.build_vanilla_uvpspec() | ||
| nt.assert_true(isinstance(uvp, uvpspec.UVPSpec)) | ||
| nt.assert_true(isinstance(cosmo, conversions.Cosmo_Conversions)) | ||
| nt.assert_equal(uvp.cosmo, cosmo) | ||
|
|
||
| beam = pspecbeam.PSpecBeamUV(os.path.join(DATA_PATH, 'NF_HERA_Beams.beamfits')) | ||
| uvp, cosmo = testing.build_vanilla_uvpspec(beam=beam) | ||
| beam_OP = beam.get_Omegas(uvp.pol_array[0])[0] | ||
| nt.assert_equal(beam_OP.tolist(), uvp.OmegaP.tolist()) | ||
|
|
||
| def test_uvpspec_from_data(): | ||
| fname = os.path.join(DATA_PATH, "zen.even.xx.LST.1.28828.uvOCRSA") | ||
| uvd = UVData() | ||
| uvd.read_miriad(fname) | ||
| beamfile = os.path.join(DATA_PATH, 'NF_HERA_Beams.beamfits') | ||
| beam = pspecbeam.PSpecBeamUV(beamfile) | ||
|
|
||
| uvp = testing.uvpspec_from_data(fname, [(37, 38), (38, 39), (52, 53), (53, 54)], beam=beam) | ||
| nt.assert_equal(uvp.Nfreqs, 150) | ||
| nt.assert_equal(np.unique(uvp.blpair_array).tolist(), [37038038039, 37038052053, 37038053054, 38039037038, | ||
| 38039052053, 38039053054, 52053037038, 52053038039, | ||
| 52053053054, 53054037038, 53054038039, 53054052053]) | ||
| uvp2 = testing.uvpspec_from_data(uvd, [(37, 38), (38, 39), (52, 53), (53, 54)], beam=beamfile) | ||
| uvp.history = '' | ||
| uvp2.history = '' | ||
| nt.assert_equal(uvp, uvp2) | ||
|
|
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.