diff --git a/CHANGELOG.md b/CHANGELOG.md index 091e31505..a9fffc395 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,8 +1,9 @@ # Changelog -## [Unreleased] +## [1.2.0] - 2020-7-20 ### Added +- Diffuse models in mock catalogs, via the analytic_diffuse module. - quantity_shared_bcast function to allow objects derived from astropy.units.Quantity to use shared memory broadcasting. - SkyModelData class to replace recarray conversion in pyradiosky. - quiet keyword for run_uvsim, to suppress stdout printing. @@ -12,14 +13,22 @@ - Benchmarking tools. ### Changed +- Switch to a new definition for the Counter class, without threading. +- Cleaned up unit tests - Use `at_frequencies` method to enable support for all pyradiosky spectral types. - Only do coherency calculation when the time changes - Only do beam eval when time, freq, or beam type changes. - The definition of the Airy beam now uses the exact value of c, not 3e8. ### Fixed +- Keep a UVBeam with more than two frequencies for tests, so the (default) cubic interpolation works. +- Use pytest hooks to ensure the profiler tests run last. +- Switched to using a remote-memory-access based counter class. Appears to have fixed bug in the Counter test. - Ensure that source positions, coherency matrices, and Jones matrices are updated at the right times. -- Error early if the task list is too long for gather +- Error early if the task list is too long for gather. + +### Deprecated +- Support for pyradiosky versions <0.1.0. ## [1.1.2] - 2020-2-14 diff --git a/ci/tests.yaml b/ci/tests.yaml index cf6a287f2..83f0a3d7e 100644 --- a/ci/tests.yaml +++ b/ci/tests.yaml @@ -19,5 +19,6 @@ dependencies: - scipy - setuptools_scm - pip: - - pyradiosky + - pyradiosky>=0.1.0 - lunarsky + - git+git://github.com/aelanman/analytic_diffuse diff --git a/docs/parameter_files.rst b/docs/parameter_files.rst index ba4d57a5b..29c059c66 100644 --- a/docs/parameter_files.rst +++ b/docs/parameter_files.rst @@ -149,8 +149,9 @@ Telescope Configuration This yaml file provides the telescope name, location in latitude/longitude/altitude in degrees/degrees/meters (respectively), and the `beam dictionary`. In this case, beam_id == 0 is the UVBeam file hera.uvbeam, beam_id == 1 - is an Airy disk with diameter 14 m, beam_id == 2 is a Gaussian beam with - sigma 0.03, and beam_id == 3 is another Airy beam. When specifying a shape + is an Airy disk with diameter 16 m, beam_id == 2 is a Gaussian beam with + sigma 0.03, beam_id == 3 is another Airy beam with diameter 12 m, 4 is a Gaussian + with diameter 14 m, and 5 is a Gaussian with with diameter 12 m. When specifying a shape parameter for a specific beam_id, the beam type needs to be specified using the type keyword (rather than on the same line as the beam_id) and then the shape keyword can be specified in the next line at the same indent level. diff --git a/environment.yml b/environment.yml index bc3701939..786fbea37 100644 --- a/environment.yml +++ b/environment.yml @@ -22,5 +22,6 @@ dependencies: - sphinx - pip: - pre-commit - - pyradiosky + - pyradiosky>=0.1.0 - lunarsky + - git+git://github.com/aelanman/analytic_diffuse diff --git a/pyuvsim/antenna.py b/pyuvsim/antenna.py index b76253cd1..e700c7396 100644 --- a/pyuvsim/antenna.py +++ b/pyuvsim/antenna.py @@ -103,6 +103,7 @@ def get_beam_jones(self, array, source_alt_az, frequency, reuse_spline=True, # interp_data has shape: # (Naxes_vec, Nspws, Nfeeds, 1 (freq), Ncomponents (source positions)) jones_matrix = np.zeros((2, 2, Ncomponents), dtype=np.complex) + # first axis is feed, second axis is theta, phi (opposite order of beam!) jones_matrix[0, 0] = interp_data[1, 0, 0, 0, :] jones_matrix[1, 1] = interp_data[0, 0, 1, 0, :] diff --git a/pyuvsim/data/28m_triangle_10time_10chan.uvfits b/pyuvsim/data/28m_triangle_10time_10chan.uvfits index 16856b170..cfb473d31 100644 --- a/pyuvsim/data/28m_triangle_10time_10chan.uvfits +++ b/pyuvsim/data/28m_triangle_10time_10chan.uvfits @@ -1,3 +1,4 @@ -SIMPLE = T / conforms to FITS standard BITPIX = -32 / array data type NAXIS = 7 / number of array dimensions NAXIS1 = 0 NAXIS2 = 3 NAXIS3 = 4 NAXIS4 = 10 NAXIS5 = 1 NAXIS6 = 1 NAXIS7 = 1 EXTEND = T GROUPS = T / has groups PCOUNT = 9 / number of parameters GCOUNT = 100 / number of groups PTYPE1 = 'UU ' PTYPE2 = 'VV ' PTYPE3 = 'WW ' PTYPE4 = 'DATE ' PTYPE5 = 'BASELINE' PTYPE6 = 'ANTENNA1' PTYPE7 = 'ANTENNA2' PTYPE8 = 'SUBARRAY' PTYPE9 = 'INTTIM ' PSCAL1 = 1.0 PZERO1 = 0.0 PSCAL2 = 1.0 PZERO2 = 0.0 PSCAL3 = 1.0 PZERO3 = 0.0 PSCAL4 = 1.0 PZERO4 = 2457458.25 PSCAL5 = 1.0 PZERO5 = 0.0 PSCAL6 = 1.0 PZERO6 = 0.0 PSCAL7 = 1.0 PZERO7 = 0.0 PSCAL8 = 1.0 PZERO8 = 0.0 PSCAL9 = 1.0 PZERO9 = 0.0 DATE-OBS= '2016-03-10T16:10:24.524' CTYPE2 = 'COMPLEX ' CRVAL2 = 1.0 CRPIX2 = 1.0 CDELT2 = 1.0 CTYPE3 = 'STOKES ' CRVAL3 = -8 CRPIX3 = 1.0 CDELT3 = 1 CTYPE4 = 'FREQ ' CRVAL4 = 100.5 CRPIX4 = 1.0 CDELT4 = 1.0 CTYPE5 = 'IF ' CRVAL5 = 1.0 CRPIX5 = 1.0 CDELT5 = 1.0 CTYPE6 = 'RA ' CRVAL6 = 72.63972830227426 CTYPE7 = 'DEC ' CRVAL7 = -30.74139849838783 BUNIT = 'Jy ' BSCALE = 1.0 BZERO = 0.0 OBJECT = 'None ' TELESCOP= '28m_triangle_10time_10chan.yaml' LAT = -30.72152777777791 LON = 21.42830555555556 ALT = 1073.000000009313 INSTRUME= 'HERA ' EPOCH = 2000.0 PHSFRAME= 'icrs ' OBS_PARA= 'param_10time_10chan_0.yaml' ANTENNA_= 'triangle_bl_layout.csv' HISTORY Test file made from param_10time_10chan_0.yaml Read/written with pyuvdaHISTORY ta version: 1.3.5. Git origin: https://github.com/HERA-Team/pyuvdata.giHISTORY t. Git hash: 62d1c7b11abfb158337fa1569f20e969e48ffec2. Git branch: masHISTORY ter. Git description: v1.3-318-g62d1c7b. END €€€½›ÜöC€€?€?€?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€3ÈŽô/4Ëë©@m½›ÜöD@@?€?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€2éù53ëK)Ñô%½›ÜöD@@@@?€?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€2êø¯.S¯³Ó½›ÜöD€ @€?€?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€€€½›ÜöD€@@?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€³Ž¦3ä* h½›ÜöD@€@@@?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€³ÐÍ®ÿƒ'³Òä½›ÜöD€@@€@?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€€€½›ÜöD@À@@@@?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€.ÿf.³ÐųÒö½›ÜöD€`@€@@?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€€€½›ÜöD€€@€@€?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€€€½›š6C€€?€?€?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€3ÈŽü/ -‹i®ŽCS½›š6D@@?€?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€2ê4“3å­¦«½›š6D@@@@?€?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€2ê”p.!ÿ³Ýb½›š6D€ @€?€?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€€€½›š6D€@@?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€³Žâ3ŸÀ.HI¸½›š6D@€@@@?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€³éÛ®Ä^‚³¹î½›š6D€@@€@?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€€€½›š6D@À@@@@?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€.?⛳Ðé³Ó½›š6D€`@€@@?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€€€½›š6D€€@€@€?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€€€½›WvC€€?€?€?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€3È.Á*¯ Ì:½›WvD@@?€?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€2êpO3Þö®&p7½›WvD@@@@?€?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€2ê06-àù}³ç¾½›WvD€ @€?€?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€€€½›WvD€@@?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€³òè3®ª.ÈRÀ½›WvD@€@@@?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€³ŽÑ®ˆðγ á½›WvD€@@€@?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€€€½›WvD@À@@@@?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€­þ˜£³Ðç³Òõ½›WvD€`@€@@?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€€€½›WvD€€@€@€?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€€€½›¶C€€?€?€?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€3ÈŽÖ.Z¯T‹¹½›¶D@@?€?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€2ê«f3ØÕ®y¿½›¶D@@@@?€?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€2éÌ&-}²î³ñü½›¶D€ @€?€?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€€€½›¶D€@@?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€³ä'3½—/QŽ›¶D@€@@@?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€³Žò® ¬³‡Ë½›¶D€@@€@?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€€€½›¶D@À@@@@?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€®ß±l³Ðæ³ÒÔ½›¶D€`@€@@?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€€€½›¶D€€@€@€?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€€€½šÑöC€€?€?€?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€3ÈŽÕ-H‚¯š:½šÑöD@@?€?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€2êç"3Ò±®¦8½šÑöD@@@@?€?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€2ég®,^ ܳüX½šÑöD€ @€?€?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€€€½šÑöD€@@?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€³Õ3Ì/He8½šÑöD@€@@@?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€³Ž4î­:I³n£½šÑöD€@@€@?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€€€½šÑöD@À@@@@?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€¯?m¿³Ðé³Ò½½šÑöD€`@€@@?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€€€½šÑöD€€@€@€?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€€€½š6C€€?€?€?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€3ÈŽ¼­ìh:¯±ü½š6D@@?€?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€2ë"3Ì‚®Ï៽š6D@@@@?€?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€2é—­ )—³Žš½š6D€ @€?€?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€€€½š6D€@@?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€³Æ3ÛT/z>C½š6D@€@@@?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€³ŽMÖ-¥Ón³U–½š6D€@@€@?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€€€½š6D@À@@@@?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€¯‡¹°³Ðé³Òœ½š6D€`@€@@?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€€€½š6D€€@€@€?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€€€½šLwC€€?€?€?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€3ÈŽy®6¯Ôf½šLwD@@?€?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€2ë]Ý3ÆM®ùܽšLwD@@@@?€?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€2蟭¨¦Ñ³Ž×½šLwD€ @€?€?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€€€½šLwD€@@?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€³·03ê0/–K½šLwD@€@@@?€A0€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€€?€³ŽfÃ.J(³= 2 and version[2] < 1): - with pytest.raises(TypeError, match='pyuvdata version >=2.0.1'): - antenna.get_beam_jones(array, altaz, 150e6) - else: - array.beam_list.spline_interp_opts = None - antenna.get_beam_jones(array, altaz, 150e6) + array.beam_list.spline_interp_opts = None + antenna.get_beam_jones(array, altaz, 150e6) + + +def test_jones_set_interp(cst_beam, hera_loc): + # check setting the interpolation method + + array_location = hera_loc + + beam = cst_beam.copy() + beam.freq_interp_kind = None + beam.interpolation_function = 'az_za_simple' + beam_list = pyuvsim.BeamList([beam]) + antenna1 = pyuvsim.Antenna('ant1', 1, np.array([0, 10, 0]), 0) + array = pyuvsim.Telescope('telescope_name', array_location, beam_list) + source_altaz = np.array([[0.0], [np.pi / 4.]]) + freq = 123e6 * units.Hz + + with pytest.raises(ValueError, match='freq_interp_kind must be set'): + antenna1.get_beam_jones(array, source_altaz, freq) + + jones = antenna1.get_beam_jones(array, source_altaz, freq, freq_interp_kind='cubic') + assert beam.freq_interp_kind == 'cubic' + jones0 = antenna1.get_beam_jones(array, source_altaz, freq) + jones1 = antenna1.get_beam_jones(array, source_altaz, freq, freq_interp_kind='linear') + assert beam.freq_interp_kind == 'linear' + jones2 = antenna1.get_beam_jones(array, source_altaz, freq) + + assert (np.all(jones2 == jones0) + and np.all(jones1 == jones) + and np.all(jones1 == jones0)) diff --git a/pyuvsim/tests/test_mpi.py b/pyuvsim/tests/test_mpi.py new file mode 100644 index 000000000..1eb2be743 --- /dev/null +++ b/pyuvsim/tests/test_mpi.py @@ -0,0 +1,198 @@ +# -*- mode: python; coding: utf-8 -* +# Copyright (c) 2018 Radio Astronomy Software Group +# Licensed under the 3-clause BSD License + +import numpy as np +import resource +import time +import pytest +import sys + +pytest.importorskip('mpi4py') # noqa +import mpi4py +mpi4py.rc.initialize = False # noqa +from mpi4py import MPI +from astropy import units +from astropy.coordinates import Latitude, Longitude + +import pyuvsim +from pyuvsim import mpi +from pyuvsim.astropy_interface import Time +import pyradiosky + + +@pytest.fixture(scope='module', autouse=True) +def start_mpi(): + mpi.start_mpi() + + +@pytest.fixture(scope='module') +def fake_tasks(): + sky = pyradiosky.SkyModel( + 'src', Longitude('10d'), Latitude('5d'), np.array([1, 0, 0, 0]) * units.Jy, + 'spectral_index', reference_frequency=np.array([100e6]) * units.Hz, + spectral_index=np.array([-0.74]) + ) + n_tasks = 30 + t0 = Time.now() + freq = 50 * units.Hz + objs = [pyuvsim.UVTask(sky, t0, freq, None, None, freq_i=ii) for ii in range(n_tasks)] + for ti, task in enumerate(objs): + task.visibility_vector = np.random.uniform(0, 3) + np.random.uniform(0, 3) * 1j + task.uvdata_index = (ti, 0, 0) + task.sources = 1 # Replace with something easier to compare later. + objs[ti] = task + + return objs + + +def test_mpi_version(): + assert MPI.VERSION == 3 + + +def test_mpi_funcs(): + assert mpi.get_rank() == MPI.COMM_WORLD.rank + assert mpi.get_Npus() == MPI.COMM_WORLD.size + assert isinstance(mpi.get_comm(), MPI.Intracomm) + assert isinstance(mpi.get_comm(), MPI.Intracomm) + assert isinstance(mpi.get_node_comm(), MPI.Intracomm) + + +def test_shared_mem(): + N = 200 + shape = (20, 10) + A = np.arange(N, dtype=float).reshape(shape) + + sA = mpi.shared_mem_bcast(A) + + # Equivalent to original + assert np.all(sA == A) + + # Not the original object: + assert hex(id(sA)) != hex(id(A)) + + # Shared array should be read-only + pytest.raises(ValueError, sA.itemset, 0, 3.0) + + +def test_mem_usage(): + # Check that the mpi-enabled memory check is consistent + # with a local memory check. + + # Also check that making a variable of a given size + # increases memory usage by the expected amount. + + scale = 1.0 + if 'linux' in sys.platform: + scale = 2**10 + + memory_usage_GiB = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss * scale / 2**30 + assert np.isclose(memory_usage_GiB, mpi.get_max_node_rss()) + incsize = 50 * 2**20 # 50 MiB + arr = bytearray(incsize) + time.sleep(1) + change = mpi.get_max_node_rss() - memory_usage_GiB + del arr + assert np.isclose(change, incsize / 2**30, atol=5e-2) + + +def test_mpi_counter(): + count = mpi.Counter() + N = 20 + for i in range(N): + count.next() + assert count.current_value() == N * mpi.world_comm.size + count.free() + + +@pytest.mark.parametrize('MAX_BYTES', [mpi.INT_MAX, 100]) +def test_big_gather(MAX_BYTES, fake_tasks): + + objs = fake_tasks + n_tasks = len(objs) + result, split_info = mpi.big_gather( + mpi.world_comm, objs, root=0, return_split_info=True, MAX_BYTES=MAX_BYTES + ) + if mpi.rank == 0: + assert all(objs[ii].freq_i == ii for ii in range(n_tasks)) + assert all(result[0][ii].freq == objs[ii].freq for ii in range(n_tasks)) + assert all(result[0][ii].time == objs[ii].time for ii in range(n_tasks)) + + # Compare with normal gather: + result2 = mpi.world_comm.gather(objs, root=0) + + if mpi.rank == 0: + assert result2 == result + + assert split_info['MAX_BYTES'] == MAX_BYTES + if MAX_BYTES < 200: + assert len(split_info['ranges']) > 1 + + +@pytest.mark.parametrize('MAX_BYTES', [mpi.INT_MAX, 100]) +def test_big_bcast(MAX_BYTES, fake_tasks): + + objs = fake_tasks + n_tasks = len(objs) + + result, split_info = mpi.big_bcast( + mpi.world_comm, objs, root=0, return_split_info=True, MAX_BYTES=MAX_BYTES + ) + + if mpi.rank == 0: + assert all(result[ii].freq_i == ii for ii in range(n_tasks)) + assert all(result[ii].freq == objs[ii].freq for ii in range(n_tasks)) + assert all(result[ii].time == objs[ii].time for ii in range(n_tasks)) + + # Compare with normal gather: + result2 = mpi.world_comm.bcast(objs, root=0) + + if mpi.rank == 0: + assert result2 == result + assert split_info['MAX_BYTES'] == MAX_BYTES + if MAX_BYTES < 200: + assert len(split_info['ranges']) > 1 + + +def test_big_bcast_gather_loop(fake_tasks): + + objs = fake_tasks + + broadcast = mpi.big_bcast(mpi.world_comm, objs, root=0, MAX_BYTES=35) + gathered = mpi.big_gather(mpi.world_comm, broadcast, root=0, MAX_BYTES=27) + + if mpi.rank == 0: + assert broadcast == gathered[0] + + +def test_sharedmem_bcast_with_quantities(): + # Use mpi.quantity_shared_bcast and check returned objects. + + lats = Latitude(np.linspace(-np.pi / 2, np.pi / 2, 10), 'rad') + freqs = np.linspace(100e6, 130e6, 15) * units.Hz + lat_return = mpi.quantity_shared_bcast(lats) + freq_return = mpi.quantity_shared_bcast(freqs) + + if mpi.rank == 0: + assert np.all(lat_return == lats) + assert np.all(freq_return == freqs) + assert np.all(freq_return.to("MHz") == freqs.to("MHz")) + + +def test_skymodeldata_share(): + # Test the SkyModelData share method. + sky = pyradiosky.SkyModel( + 'src', Longitude('10d'), Latitude('5d'), np.array([1, 0, 0, 0]) * units.Jy, + 'spectral_index', reference_frequency=np.array([100e6]) * units.Hz, + spectral_index=np.array([-0.74]) + ) + + smd = pyuvsim.simsetup.SkyModelData() + if mpi.rank == 0: + smd = pyuvsim.simsetup.SkyModelData(sky) + + smd.share() # Shares among processes. + + sky2 = smd.get_skymodel() + + assert sky2 == sky diff --git a/pyuvsim/tests/test_mpi_uvsim.py b/pyuvsim/tests/test_mpi_uvsim.py deleted file mode 100644 index 7ef86dcc0..000000000 --- a/pyuvsim/tests/test_mpi_uvsim.py +++ /dev/null @@ -1,311 +0,0 @@ -# -*- mode: python; coding: utf-8 -* -# Copyright (c) 2018 Radio Astronomy Software Group -# Licensed under the 3-clause BSD License - -import os - -import numpy as np -import yaml -import resource -import time -import pytest -import sys - -pytest.importorskip('mpi4py') # noqa -import mpi4py -mpi4py.rc.initialize = False # noqa -from mpi4py import MPI -from astropy import units -from astropy.coordinates import Latitude, Longitude - -from pyuvdata import UVData -from pyuvdata.data import DATA_PATH - -from pyuvsim.data import DATA_PATH as SIM_DATA_PATH -import pyuvsim -from pyuvsim import mpi -from pyuvsim.astropy_interface import Time -import pyuvsim.tests as simtest -import pyradiosky - - -cst_files = ['HERA_NicCST_150MHz.txt', 'HERA_NicCST_123MHz.txt'] -beam_files = [os.path.join(DATA_PATH, 'NicCSTbeams', f) for f in cst_files] -hera_miriad_file = os.path.join(DATA_PATH, 'hera_testfile') -EW_uvfits_file = os.path.join(SIM_DATA_PATH, '28mEWbl_1time_1chan.uvfits') -param_filenames = [ - os.path.join(SIM_DATA_PATH, 'test_config', 'param_10time_10chan_{}.yaml'.format(x)) - for x in range(4) -] # Five different test configs -singlesource_vot = os.path.join(SIM_DATA_PATH, 'single_source.vot') -singlesource_txt = os.path.join(SIM_DATA_PATH, 'single_source.txt') - - -@pytest.fixture -def fake_tasks(): - sky = pyradiosky.SkyModel( - 'src', Longitude('10d'), Latitude('5d'), [1, 0, 0, 0], 'spectral_index', - reference_frequency=np.array([100e6]) * units.Hz, spectral_index=np.array([-0.74]) - ) - n_tasks = 30 - t0 = Time.now() - freq = 50 * units.Hz - objs = [pyuvsim.UVTask(sky, t0, freq, None, None, freq_i=ii) for ii in range(n_tasks)] - for ti, task in enumerate(objs): - task.visibility_vector = np.random.uniform(0, 3) + np.random.uniform(0, 3) * 1j - task.uvdata_index = (ti, 0, 0) - task.sources = 1 # Replace with something easier to compare later. - objs[ti] = task - - return objs - - -def test_mpi_version(): - assert MPI.VERSION == 3 - - -def test_run_uvsim(): - hera_uv = UVData() - hera_uv.read_uvfits(EW_uvfits_file) - beam = simtest.make_cst_beams(freqs=[100e6, 123e6]) - beamfile = os.path.join(simtest.TESTDATA_PATH, "temp.uvbeam") - beam.write_beamfits(beamfile) - beam_list = pyuvsim.BeamList([beamfile]) - mock_keywords = {"Nsrcs": 3} - - with pytest.raises(TypeError, match='input_uv must be UVData object'): - pyuvsim.run_uvdata_uvsim('not_uvdata', beam_list, quiet=True) - - mpi.start_mpi() - catalog, mock_kwds = pyuvsim.simsetup.create_mock_catalog( - hera_uv.time_array[0], return_data=True, **mock_keywords - ) - uv_out = pyuvsim.run_uvdata_uvsim(hera_uv, beam_list, catalog=catalog) - rank = mpi.get_rank() - if rank == 0: - assert np.allclose(uv_out.data_array, hera_uv.data_array, atol=5e-3) - os.remove(beamfile) - - -@pytest.mark.filterwarnings("ignore:The frequency field is included in the recarray") -def test_run_paramfile_uvsim(): - # Test vot and txt catalogs for parameter simulation - - uv_ref = UVData() - uv_ref.read_uvfits(os.path.join(SIM_DATA_PATH, 'testfile_singlesource.uvfits')) - uv_ref.unphase_to_drift(use_ant_pos=True) - - param_filename = os.path.join(SIM_DATA_PATH, 'test_config', 'param_1time_1src_testcat.yaml') - with open(param_filename) as pfile: - params_dict = yaml.safe_load(pfile) - tempfilename = params_dict['filing']['outfile_name'] - - # This test obsparam file has "single_source.txt" as its catalog. - pyuvsim.uvsim.run_uvsim(param_filename) - - uv_new_txt = UVData() - with pytest.warns(UserWarning) as antdiam: - uv_new_txt.read_uvfits(tempfilename) - assert str(antdiam.pop().message).startswith('antenna_diameters is not set') - uv_new_txt.unphase_to_drift(use_ant_pos=True) - os.remove(tempfilename) - - param_filename = os.path.join(SIM_DATA_PATH, 'test_config', 'param_1time_1src_testvot.yaml') - pyuvsim.uvsim.run_uvsim(param_filename) - - uv_new_vot = UVData() - with pytest.warns(UserWarning) as antdiam: - uv_new_vot.read_uvfits(tempfilename) - assert str(antdiam.pop().message).startswith('antenna_diameters is not set') - uv_new_vot.unphase_to_drift(use_ant_pos=True) - os.remove(tempfilename) - uv_new_txt.history = uv_ref.history # History includes irrelevant info for comparison - uv_new_vot.history = uv_ref.history - uv_new_txt.object_name = uv_ref.object_name - uv_new_vot.object_name = uv_ref.object_name - assert uv_new_txt == uv_ref - assert uv_new_vot == uv_ref - - -@pytest.mark.filterwarnings("ignore:The frequency field is included in the recarray") -def test_run_paramdict_uvsim(): - # Running a simulation from parameter dictionary. - - params = pyuvsim.simsetup._config_str_to_dict( - os.path.join(SIM_DATA_PATH, 'test_config', 'param_1time_1src_testcat.yaml') - ) - - pyuvsim.run_uvsim(params, return_uv=True) - - -@pytest.mark.parametrize( - "spectral_type", - ["flat", "subband", "spectral_index"]) -def test_run_gleam_uvsim(spectral_type): - params = pyuvsim.simsetup._config_str_to_dict( - os.path.join(SIM_DATA_PATH, 'test_config', 'param_1time_1src_testgleam.yaml') - ) - params["sources"]["spectral_type"] = spectral_type - params["sources"].pop("min_flux") - params["sources"].pop("max_flux") - - pyuvsim.run_uvsim(params, return_uv=True) - - -def test_mpi_funcs(): - mpi.start_mpi() - assert mpi.get_rank() == 0 - assert mpi.get_Npus() == 1 - assert isinstance(mpi.get_comm(), MPI.Intracomm) - assert isinstance(mpi.get_comm(), MPI.Intracomm) - assert isinstance(mpi.get_node_comm(), MPI.Intracomm) - - -def test_shared_mem(): - mpi.start_mpi() - N = 200 - shape = (20, 10) - A = np.arange(N, dtype=float).reshape(shape) - - sA = mpi.shared_mem_bcast(A) - - # Equivalent to original - assert np.all(sA == A) - - # Not the original object: - assert hex(id(sA)) != hex(id(A)) - - # Shared array should be read-only - pytest.raises(ValueError, sA.itemset, 0, 3.0) - - -def test_mem_usage(): - # Check that the mpi-enabled memory check is consistent - # with a local memory check. - - # Also check that making a variable of a given size - # increases memory usage by the expected amount. - - mpi.start_mpi() - - scale = 1.0 - if 'linux' in sys.platform: - scale = 2**10 - - memory_usage_GiB = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss * scale / 2**30 - assert memory_usage_GiB == mpi.get_max_node_rss() - incsize = 50 * 2**20 # 50 MiB - arr = bytearray(incsize) - time.sleep(1) - change = mpi.get_max_node_rss() - memory_usage_GiB - assert np.isclose(change, incsize / 2**30, atol=5e-2) - del arr - - -@pytest.mark.skip -def test_mpi_counter(): - # This test should be run in parallel to check likely bug sources. - mpi.start_mpi() - - count = mpi.Counter() - N = 20 - for i in range(N): - count.next() - assert count.current_value() == N - count.free() - - -@pytest.mark.parametrize('MAX_BYTES', [mpi.INT_MAX, 100]) -def test_big_gather(MAX_BYTES, fake_tasks): - mpi.start_mpi() - - objs = fake_tasks - n_tasks = len(objs) - result, split_info = mpi.big_gather( - mpi.world_comm, objs, root=0, return_split_info=True, MAX_BYTES=MAX_BYTES - ) - assert all(objs[ii].freq_i == ii for ii in range(n_tasks)) - assert all(result[0][ii].freq == objs[ii].freq for ii in range(n_tasks)) - assert all(result[0][ii].time == objs[ii].time for ii in range(n_tasks)) - - # Compare with normal gather: - result2 = mpi.world_comm.gather(objs, root=0) - - assert result2 == result - - assert split_info['MAX_BYTES'] == MAX_BYTES - if MAX_BYTES < 200: - assert len(split_info['ranges']) > 1 - - -@pytest.mark.parametrize('MAX_BYTES', [mpi.INT_MAX, 100]) -def test_big_bcast(MAX_BYTES, fake_tasks): - mpi.start_mpi() - - objs = fake_tasks - n_tasks = len(objs) - - result, split_info = mpi.big_bcast( - mpi.world_comm, objs, root=0, return_split_info=True, MAX_BYTES=MAX_BYTES - ) - assert all(result[ii].freq_i == ii for ii in range(n_tasks)) - assert all(result[ii].freq == objs[ii].freq for ii in range(n_tasks)) - assert all(result[ii].time == objs[ii].time for ii in range(n_tasks)) - - # Compare with normal gather: - result2 = mpi.world_comm.bcast(objs, root=0) - - assert result2 == result - - assert split_info['MAX_BYTES'] == MAX_BYTES - if MAX_BYTES < 200: - assert len(split_info['ranges']) > 1 - - -def test_big_bcast_gather_loop(fake_tasks): - mpi.start_mpi() - - objs = fake_tasks - - broadcast = mpi.big_bcast(mpi.world_comm, objs, root=0, MAX_BYTES=35) - gathered = mpi.big_gather(mpi.world_comm, broadcast, root=0, MAX_BYTES=27) - - assert broadcast == gathered[0] - - -@pytest.mark.skipif('not pyuvsim.astropy_interface.hasmoon') -def test_sim_on_moon(): - from pyuvsim.astropy_interface import MoonLocation - param_filename = os.path.join(SIM_DATA_PATH, 'test_config', 'obsparam_tranquility_hex.yaml') - param_dict = pyuvsim.simsetup._config_str_to_dict(param_filename) - param_dict['select'] = {'redundant_threshold': 0.1} - uv_obj, beam_list, beam_dict = pyuvsim.initialize_uvdata_from_params(param_dict) - uv_obj.select(times=uv_obj.time_array[0]) - tranquility_base = MoonLocation.from_selenocentric(*uv_obj.telescope_location, 'meter') - - time = Time(uv_obj.time_array[0], format='jd', scale='utc') - sources, kwds = pyuvsim.create_mock_catalog( - time, array_location=tranquility_base, arrangement='zenith', Nsrcs=30, return_data=True - ) - print(uv_obj.extra_keywords['world']) - # Run simulation. - uv_out = pyuvsim.uvsim.run_uvdata_uvsim( - uv_obj, beam_list, beam_dict, catalog=sources, quiet=True - ) - assert np.allclose(uv_out.data_array[:, 0, :, 0], 0.5) - assert uv_out.extra_keywords['world'] == 'moon' - - -def test_sharedmem_bcast_with_quantities(): - # Use mpi.quantity_shared_bcast and check returned objects. - - mpi.start_mpi() - lats = Latitude(np.linspace(-np.pi / 2, np.pi / 2, 10), 'rad') - freqs = np.linspace(100e6, 130e6, 15) * units.Hz - lat_return = mpi.quantity_shared_bcast(lats) - freq_return = mpi.quantity_shared_bcast(freqs) - - assert np.all(lat_return == lats) - assert np.all(freq_return == freqs) - - assert np.all(freq_return.to("MHz") == freqs.to("MHz")) diff --git a/pyuvsim/tests/test_profiler.py b/pyuvsim/tests/test_profiler.py new file mode 100644 index 000000000..a3c294994 --- /dev/null +++ b/pyuvsim/tests/test_profiler.py @@ -0,0 +1,44 @@ +# -*- mode: python; coding: utf-8 -* +# Copyright (c) 2018 Radio Astronomy Software Group +# Licensed under the 3-clause BSD License + +from numpy import unique +import os +import atexit +import shutil +import pytest + +import pyuvsim +from pyuvsim.data import DATA_PATH as SIM_DATA_PATH + + +def profdata_dir_setup(tmpdir): + # The pytest temporary test directory doesn't work with the profiling outputs + # because they are only created after the exit functions are run, and pytest + # deletes its temporary directory before then. + # As a workaround, make a temporary directory for profiling data and register + # its deletion in the exit functions. + outpath = tmpdir.mkdir('temporary_profiling_data') + atexit.register(shutil.rmtree, outpath) + return outpath + + +def test_profiler(tmpdir): + line_profiler = pytest.importorskip('line_profiler') + outpath = profdata_dir_setup(tmpdir) + testprof_fname = str(outpath.join('time_profile.out')) + print(testprof_fname) + pyuvsim.profiling.set_profiler(outfile_prefix=testprof_fname, dump_raw=True) + with pytest.warns(UserWarning, match='Profiler already set'): + pyuvsim.profiling.set_profiler(outfile_prefix=testprof_fname[:-4], dump_raw=True) + param_filename = os.path.join(SIM_DATA_PATH, 'test_config', 'param_1time_1src_testcat.yaml') + pyuvsim.uvsim.run_uvsim(param_filename, return_uv=True) + time_profiler = pyuvsim.profiling.get_profiler() + time_profiler.disable_by_count() + assert isinstance(time_profiler, line_profiler.LineProfiler) + assert hasattr(time_profiler, 'rank') + assert hasattr(time_profiler, 'meta_file') + lstats = time_profiler.get_stats() + assert len(lstats.timings) != 0 + func_names = [k[2] for k in lstats.timings.keys()] + assert unique(func_names).tolist() == sorted(pyuvsim.profiling.default_profile_funcs) diff --git a/pyuvsim/tests/test_run.py b/pyuvsim/tests/test_run.py new file mode 100644 index 000000000..2502aeb41 --- /dev/null +++ b/pyuvsim/tests/test_run.py @@ -0,0 +1,233 @@ +# -*- mode: python; coding: utf-8 -* +# Copyright (c) 2020 Radio Astronomy Software Group +# Licensed under the 3-clause BSD License + +import numpy as np +import os +import pytest +import yaml +from astropy import units + +from pyuvdata import UVData +from pyradiosky.utils import jy_to_ksr, stokes_to_coherency + +import pyuvsim +from pyuvsim.astropy_interface import Time +from pyuvsim.data import DATA_PATH as SIM_DATA_PATH +from pyuvsim.analyticbeam import c_ms + + +pytest.importorskip('mpi4py') # noqa + + +@pytest.mark.filterwarnings("ignore:The frequency field is included in the recarray") +def test_run_paramfile_uvsim(): + # Test vot and txt catalogs for parameter simulation + # Compare to reference files. + + uv_ref = UVData() + uv_ref.read_uvfits(os.path.join(SIM_DATA_PATH, 'testfile_singlesource.uvfits')) + uv_ref.unphase_to_drift(use_ant_pos=True) + + param_filename = os.path.join(SIM_DATA_PATH, 'test_config', 'param_1time_1src_testcat.yaml') + with open(param_filename) as pfile: + params_dict = yaml.safe_load(pfile) + tempfilename = params_dict['filing']['outfile_name'] + + # This test obsparam file has "single_source.txt" as its catalog. + pyuvsim.uvsim.run_uvsim(param_filename) + + uv_new_txt = UVData() + with pytest.warns(UserWarning, match='antenna_diameters is not set'): + uv_new_txt.read_uvfits(tempfilename) + + uv_new_txt.unphase_to_drift(use_ant_pos=True) + os.remove(tempfilename) + + param_filename = os.path.join(SIM_DATA_PATH, 'test_config', 'param_1time_1src_testvot.yaml') + pyuvsim.uvsim.run_uvsim(param_filename) + + uv_new_vot = UVData() + with pytest.warns(UserWarning, match='antenna_diameters is not set'): + uv_new_vot.read_uvfits(tempfilename) + + uv_new_vot.unphase_to_drift(use_ant_pos=True) + os.remove(tempfilename) + uv_new_txt.history = uv_ref.history # History includes irrelevant info for comparison + uv_new_vot.history = uv_ref.history + uv_new_txt.object_name = uv_ref.object_name + uv_new_vot.object_name = uv_ref.object_name + assert uv_new_txt == uv_ref + assert uv_new_vot == uv_ref + + +@pytest.mark.parametrize('model', ['monopole', 'cosza', 'quaddome', 'monopole-nonflat']) +def test_analytic_diffuse(model, tmpdir): + # Generate the given model and simulate for a few baselines. + # Import from analytic_diffuse (consider moving to rasg_affiliates?) + pytest.importorskip('analytic_diffuse') + pytest.importorskip('astropy_healpix') + import analytic_diffuse + + modname = model + use_w = False + params = {} + if model == 'quaddome': + modname = 'polydome' + params['n'] = 2 + elif model == 'monopole-nonflat': + modname = 'monopole' + use_w = True + params['order'] = 30 # Expansion order for the non-flat monopole solution. + + # Making configuration files for this simulation. + template_path = os.path.join(SIM_DATA_PATH, 'test_config', 'obsparam_diffuse_sky.yaml') + obspar_path = str(tmpdir.join('obsparam_diffuse_sky.yaml')) + layout_path = str(tmpdir.join('threeant_layout.csv')) + herauniform_path = str(tmpdir.join('hera_uniform.yaml')) + + teleconfig = { + 'beam_paths': {0: 'uniform'}, + 'telescope_location': "(-30.72153, 21.42830, 1073.0)", + 'telescope_name': 'HERA' + } + if not use_w: + antpos_enu = np.array([[0, 0, 0], [0, 3, 0], [5, 0, 0]], dtype=float) + else: + antpos_enu = np.array([[0, 0, 0], [0, 3, 0], [0, 3, 5]], dtype=float) + + pyuvsim.simsetup._write_layout_csv( + layout_path, antpos_enu, np.arange(3).astype(str), np.arange(3) + ) + with open(herauniform_path, 'w') as ofile: + yaml.dump(teleconfig, ofile, default_flow_style=False) + + with open(template_path, 'r') as yfile: + obspar = yaml.safe_load(yfile) + obspar['telescope']['array_layout'] = layout_path + obspar['telescope']['telescope_config_name'] = herauniform_path + obspar['sources']['diffuse_model'] = modname + obspar['sources'].update(params) + obspar['filing']['outfile_name'] = 'diffuse_sim.uvh5' + obspar['filing']['output_format'] = 'uvh5' + obspar['filing']['outdir'] = str(tmpdir) + + with open(obspar_path, 'w') as ofile: + yaml.dump(obspar, ofile, default_flow_style=False) + + uv_out = pyuvsim.run_uvsim(obspar_path, return_uv=True) + # Convert from Jy to K sr + dat = uv_out.data_array[:, 0, 0, 0] * jy_to_ksr(uv_out.freq_array[0, 0]).value + # Evaluate the solution and compare to visibilities. + soln = analytic_diffuse.get_solution(modname) + uvw_lam = uv_out.uvw_array * uv_out.freq_array[0, 0] / c_ms + ana = soln(uvw_lam, **params) + assert np.allclose(ana / 2, dat, atol=1e-2) + + +@pytest.mark.filterwarnings("ignore:The frequency field is included in the recarray") +def test_run_paramdict_uvsim(): + # Running a simulation from parameter dictionary. + + params = pyuvsim.simsetup._config_str_to_dict( + os.path.join(SIM_DATA_PATH, 'test_config', 'param_1time_1src_testcat.yaml') + ) + + pyuvsim.run_uvsim(params, return_uv=True) + + +@pytest.mark.parametrize( + "spectral_type", + ["flat", "subband", "spectral_index"]) +def test_run_gleam_uvsim(spectral_type): + params = pyuvsim.simsetup._config_str_to_dict( + os.path.join(SIM_DATA_PATH, 'test_config', 'param_1time_1src_testgleam.yaml') + ) + params["sources"]["spectral_type"] = spectral_type + params["sources"].pop("min_flux") + params["sources"].pop("max_flux") + + pyuvsim.run_uvsim(params, return_uv=True) + + +@pytest.mark.parametrize( + "spectral_type", + ["subband", "spectral_index"]) +def test_zenith_spectral_sim(spectral_type, tmpdir): + # Make a power law source at zenith in three ways. + # Confirm that simulated visibilities match expectation. + + params = pyuvsim.simsetup._config_str_to_dict( + os.path.join(SIM_DATA_PATH, 'test_config', 'param_1time_1src_testcat.yaml') + ) + + alpha = -0.5 + ref_freq = 111e6 + Nfreqs = 20 + freqs = np.linspace(110e6, 115e6, Nfreqs) + freq_params = pyuvsim.simsetup.freq_array_to_params(freqs) + freqs = pyuvsim.simsetup.parse_frequency_params(freq_params)['freq_array'][0, :] + freqs *= units.Hz + spectrum = (freqs.value / ref_freq)**alpha + + source, kwds = pyuvsim.create_mock_catalog(Time.now(), arrangement='zenith', Nsrcs=1) + source.spectral_type = spectral_type + if spectral_type == 'spectral_index': + source.reference_frequency = np.array([ref_freq]) * units.Hz + source.spectral_index = np.array([alpha]) + else: + source.Nfreqs = Nfreqs + source.freq_array = freqs + source.stokes = np.repeat(source.stokes, Nfreqs, axis=1) + source.stokes[0, :, 0] *= spectrum + source.coherency_radec = stokes_to_coherency(source.stokes) + + catpath = str(tmpdir.join('spectral_test_catalog.txt')) + source.write_text_catalog(catpath) + params['sources'] = {"catalog" : catpath} + params['filing']['outdir'] = str(tmpdir) + params['freq'] = freq_params + params['time']['start_time'] = kwds['time'] + params['select'] = {'antenna_nums' : [1, 2]} + + uv_out = pyuvsim.run_uvsim(params, return_uv=True) + + for ii in range(uv_out.Nbls): + assert np.allclose(uv_out.data_array[ii, 0, :, 0], spectrum / 2) + + +def test_pol_error(): + # Check that running with a uvdata object without the proper polarizations will fail. + hera_uv = UVData() + + hera_uv.polarizations = ['xx'] + + with pytest.raises(ValueError, match='input_uv must have XX,YY,XY,YX polarization'): + pyuvsim.run_uvdata_uvsim(hera_uv, ['beamlist']) + + +def test_input_uv_error(): + with pytest.raises(TypeError, match="input_uv must be UVData object"): + pyuvsim.run_uvdata_uvsim(None, None) + + +@pytest.mark.skipif('not pyuvsim.astropy_interface.hasmoon') +def test_sim_on_moon(): + from pyuvsim.astropy_interface import MoonLocation + param_filename = os.path.join(SIM_DATA_PATH, 'test_config', 'obsparam_tranquility_hex.yaml') + param_dict = pyuvsim.simsetup._config_str_to_dict(param_filename) + param_dict['select'] = {'redundant_threshold': 0.1} + uv_obj, beam_list, beam_dict = pyuvsim.initialize_uvdata_from_params(param_dict) + uv_obj.select(times=uv_obj.time_array[0]) + tranquility_base = MoonLocation.from_selenocentric(*uv_obj.telescope_location, 'meter') + + time = Time(uv_obj.time_array[0], format='jd', scale='utc') + sources, kwds = pyuvsim.create_mock_catalog( + time, array_location=tranquility_base, arrangement='zenith', Nsrcs=30, return_data=True + ) + # Run simulation. + uv_out = pyuvsim.uvsim.run_uvdata_uvsim( + uv_obj, beam_list, beam_dict, catalog=sources, quiet=True + ) + assert np.allclose(uv_out.data_array[:, 0, :, 0], 0.5) + assert uv_out.extra_keywords['world'] == 'moon' diff --git a/pyuvsim/tests/test_simsetup.py b/pyuvsim/tests/test_simsetup.py index 48910cb55..2f8423c22 100644 --- a/pyuvsim/tests/test_simsetup.py +++ b/pyuvsim/tests/test_simsetup.py @@ -37,12 +37,10 @@ gleam_param_file = os.path.join(SIM_DATA_PATH, 'test_config', 'param_1time_1src_testgleam.yaml') -def test_mock_catalog_zenith_source(): +def test_mock_catalog_zenith_source(hera_loc): time = Time(2457458.65410, scale='utc', format='jd') - array_location = EarthLocation( - lat='-30d43m17.5s', lon='21d25m41.9s', height=1073. - ) + array_location = hera_loc source_coord = SkyCoord( alt=Angle(90 * units.deg), az=Angle(0 * units.deg), @@ -53,22 +51,20 @@ def test_mock_catalog_zenith_source(): ra = icrs_coord.ra dec = icrs_coord.dec - test_source = pyradiosky.SkyModel('src0', ra, dec, [1, 0, 0, 0], 'flat') + test_source = pyradiosky.SkyModel('src0', ra, dec, units.Quantity([1, 0, 0, 0], 'Jy'), 'flat') cat, mock_keywords = pyuvsim.create_mock_catalog(time, arrangement='zenith') assert cat == test_source -def test_mock_catalog_off_zenith_source(): +def test_mock_catalog_off_zenith_source(hera_loc): src_az = Angle('90.0d') src_alt = Angle('85.0d') time = Time(2457458.65410, scale='utc', format='jd') - array_location = EarthLocation( - lat='-30d43m17.5s', lon='21d25m41.9s', height=1073. - ) + array_location = hera_loc source_coord = SkyCoord( alt=src_alt, az=src_az, obstime=time, frame='altaz', location=array_location @@ -77,7 +73,7 @@ def test_mock_catalog_off_zenith_source(): ra = icrs_coord.ra dec = icrs_coord.dec - test_source = pyradiosky.SkyModel('src0', ra, dec, [1.0, 0, 0, 0], 'flat') + test_source = pyradiosky.SkyModel('src0', ra, dec, units.Quantity([1.0, 0, 0, 0], "Jy"), 'flat') cat, mock_keywords = pyuvsim.create_mock_catalog(time, arrangement='off-zenith', alt=src_alt.deg) @@ -85,14 +81,63 @@ def test_mock_catalog_off_zenith_source(): assert cat == test_source +def test_mock_diffuse_maps_errors(): + analytic_diffuse = pyuvsim.simsetup.analytic_diffuse + astropy_healpix = pyuvsim.simsetup.astropy_healpix + if (analytic_diffuse is not None) and (astropy_healpix is not None): + + # Error cases: + with pytest.raises(ValueError, match="Diffuse arrangement selected"): + pyuvsim.simsetup.create_mock_catalog(Time.now(), arrangement='diffuse') + + with pytest.warns(UserWarning, match="No nside chosen"): + pyuvsim.simsetup.create_mock_catalog( + Time.now(), arrangement='diffuse', diffuse_model='monopole' + ) + + else: + with pytest.raises(ValueError, match="analytic_diffuse and astropy_healpix"): + pyuvsim.simsetup.create_mock_catalog(Time.now(), arrangement='diffuse') + + +@pytest.mark.parametrize('model', [ + ['monopole', {}], + ['gauss', {'a': 0.05}], + ['polydome', {'n': 4}] + ]) +@pytest.mark.parametrize('location', ['earth', 'moon']) +def test_mock_diffuse_maps(model, hera_loc, apollo_loc, location): + analytic_diffuse = pytest.importorskip('analytic_diffuse') + pytest.importorskip('astropy_healpix') + if location == 'earth': + loc = hera_loc + else: + loc = apollo_loc + modname, modkwargs = model + map_nside = 128 + t0 = Time.now() + cat, kwds = pyuvsim.simsetup.create_mock_catalog( + t0, arrangement='diffuse', array_location=loc, + diffuse_model=modname, map_nside=map_nside, + diffuse_params=modkwargs + ) + + cat.update_positions(t0, loc) + + modfunc = analytic_diffuse.get_model(modname) + alt, az = cat.alt_az + za = np.pi / 2 - alt + + vals = modfunc(az, za, **modkwargs) + + assert cat.nside == map_nside + assert np.allclose(cat.stokes[0].to_value("K"), vals) + + def test_catalog_from_params(): # Pass in parameter dictionary as dict hera_uv = UVData() - with pytest.warns( - UserWarning, - match='Telescope 28m_triangle_10time_10chan.yaml is not in known_telescopes.' - ): - hera_uv.read_uvfits(triangle_uvfits_file) + hera_uv.read_uvfits(triangle_uvfits_file) source_dict = {} with pytest.raises(KeyError, match='No catalog defined.'): @@ -256,112 +301,114 @@ def test_vot_catalog_error(key_pop, message): pyuvsim.simsetup.initialize_catalog_from_params(param_dict, return_recarray=False)[0] -@pytest.mark.parametrize("config_num", [0, 2]) -def test_param_reader(config_num): +def test_param_reader(): pytest.importorskip('mpi4py') # Reading in various configuration files - param_filename = param_filenames[config_num] + param_filename = os.path.join(SIM_DATA_PATH, "test_config", "param_10time_10chan_0.yaml") hera_uv = UVData() - with pytest.warns(UserWarning) as warn: - hera_uv.read_uvfits(triangle_uvfits_file) - assert str(warn.pop().message).startswith('Telescope 28m_triangle_10time_10chan.yaml ' - 'is not in known_telescopes.') + hera_uv.read_uvfits(triangle_uvfits_file) + hera_uv.unphase_to_drift() hera_uv.telescope_name = 'HERA' - if config_num == 5: - hera_uv.select(bls=[(0, 1), (1, 2)]) time = Time(hera_uv.time_array[0], scale='utc', format='jd') sources, _ = pyuvsim.create_mock_catalog(time, arrangement='zenith', return_data=True) beam0 = UVBeam() beam0.read_beamfits(herabeam_default) + beam0.extra_keywords['beam_path'] = herabeam_default beam1 = pyuvsim.AnalyticBeam('uniform') beam2 = pyuvsim.AnalyticBeam('gaussian', sigma=0.02) beam3 = pyuvsim.AnalyticBeam('airy', diameter=14.6) beam_list = pyuvsim.BeamList([beam0, beam1, beam2, beam3]) + # To fill out other parameters in the UVBeam. + beam_list.set_str_mode() + beam_list.set_obj_mode() + beam_dict = {'ANT1': 0, 'ANT2': 1, 'ANT3': 2, 'ANT4': 3} - Ntasks = hera_uv.Nblts * hera_uv.Nfreqs - taskiter = pyuvsim.uvdata_to_task_iter(range(Ntasks), hera_uv, sources, - beam_list, beam_dict=beam_dict) - expected_uvtask_list = list(taskiter) - - # Check error conditions: - if config_num == 0: - params_bad = pyuvsim.simsetup._config_str_to_dict(param_filename) - bak_params_bad = copy.deepcopy(params_bad) - - # Missing config file info - params_bad['config_path'] = os.path.join( - SIM_DATA_PATH, 'nonexistent_directory', 'nonexistent_file' - ) - with pytest.raises(ValueError, match="nonexistent_directory is not a directory"): - pyuvsim.simsetup.initialize_uvdata_from_params(params_bad) - params_bad['config_path'] = os.path.join(SIM_DATA_PATH, "test_config") - params_bad['telescope']['array_layout'] = 'nonexistent_file' - with pytest.raises(ValueError, match="nonexistent_file from yaml does not exist"): - pyuvsim.simsetup.initialize_uvdata_from_params(params_bad) + # Error conditions: + params_bad = pyuvsim.simsetup._config_str_to_dict(param_filename) + bak_params_bad = copy.deepcopy(params_bad) - params_bad['telescope']['telescope_config_name'] = 'nonexistent_file' - with pytest.raises(ValueError, match="telescope_config_name file from yaml does not exist"): - pyuvsim.simsetup.initialize_uvdata_from_params(params_bad) + # Missing config file info + params_bad['config_path'] = os.path.join( + SIM_DATA_PATH, 'nonexistent_directory', 'nonexistent_file' + ) + with pytest.raises(ValueError, match="nonexistent_directory is not a directory"): + pyuvsim.simsetup.initialize_uvdata_from_params(params_bad) - # Missing beam keywords - params_bad = copy.deepcopy(bak_params_bad) + params_bad['config_path'] = os.path.join(SIM_DATA_PATH, "test_config") + params_bad['telescope']['array_layout'] = 'nonexistent_file' + with pytest.raises(ValueError, match="nonexistent_file from yaml does not exist"): + pyuvsim.simsetup.initialize_uvdata_from_params(params_bad) - params_bad['config_path'] = os.path.join(SIM_DATA_PATH, "test_config") + params_bad['telescope']['telescope_config_name'] = 'nonexistent_file' + with pytest.raises(ValueError, match="telescope_config_name file from yaml does not exist"): + pyuvsim.simsetup.initialize_uvdata_from_params(params_bad) - params_bad = copy.deepcopy(bak_params_bad) - params_bad['telescope']['telescope_config_name'] = os.path.join( - SIM_DATA_PATH, 'test_config', '28m_triangle_10time_10chan_gaussnoshape.yaml' - ) - with pytest.raises(KeyError, - match="Missing shape parameter for gaussian beam"): - pyuvsim.simsetup.initialize_uvdata_from_params(params_bad) + # Missing beam keywords + params_bad = copy.deepcopy(bak_params_bad) - params_bad['telescope']['telescope_config_name'] = os.path.join( - SIM_DATA_PATH, 'test_config', '28m_triangle_10time_10chan_nodiameter.yaml' - ) - with pytest.raises(KeyError, match="Missing diameter for airy beam."): - pyuvsim.simsetup.initialize_uvdata_from_params(params_bad) + params_bad['config_path'] = os.path.join(SIM_DATA_PATH, "test_config") - params_bad['telescope']['telescope_config_name'] = os.path.join( - SIM_DATA_PATH, 'test_config', '28m_triangle_10time_10chan_nofile.yaml' - ) - with pytest.raises(ValueError, match="Undefined beam model"): - pyuvsim.simsetup.initialize_uvdata_from_params(params_bad) + params_bad = copy.deepcopy(bak_params_bad) + params_bad['telescope']['telescope_config_name'] = os.path.join( + SIM_DATA_PATH, 'test_config', '28m_triangle_10time_10chan_gaussnoshape.yaml' + ) + with pytest.raises(KeyError, + match="Missing shape parameter for gaussian beam"): + pyuvsim.simsetup.initialize_uvdata_from_params(params_bad) + + params_bad['telescope']['telescope_config_name'] = os.path.join( + SIM_DATA_PATH, 'test_config', '28m_triangle_10time_10chan_nodiameter.yaml' + ) + with pytest.raises(KeyError, match="Missing diameter for airy beam."): + pyuvsim.simsetup.initialize_uvdata_from_params(params_bad) + + params_bad['telescope']['telescope_config_name'] = os.path.join( + SIM_DATA_PATH, 'test_config', '28m_triangle_10time_10chan_nofile.yaml' + ) + with pytest.raises(ValueError, match="Undefined beam model"): + pyuvsim.simsetup.initialize_uvdata_from_params(params_bad) # Check default configuration uv_obj, new_beam_list, new_beam_dict = pyuvsim.initialize_uvdata_from_params(param_filename) new_beam_list.set_obj_mode() + pyuvsim.simsetup._complete_uvdata(uv_obj, inplace=True) + with open(param_filename, 'r') as fhandle: param_dict = yaml.safe_load(fhandle) expected_ofilepath = pyuvsim.utils.write_uvdata( uv_obj, param_dict, return_filename=True, dryrun=True ) - ofilename = 'sim_results.uvfits' - if config_num == 1: - if os.path.isdir('tempdir'): - os.rmdir('tempdir') - ofilename = os.path.join('.', 'tempdir', ofilename) - else: - ofilename = os.path.join('.', ofilename) - assert ofilename == expected_ofilepath - - Ntasks = uv_obj.Nblts * uv_obj.Nfreqs - taskiter = pyuvsim.uvdata_to_task_iter( - range(Ntasks), hera_uv, sources, beam_list, beam_dict=beam_dict - ) - uvtask_list = list(taskiter) + assert './sim_results.uvfits' == expected_ofilepath + + # Spoof attributes that won't match. + uv_obj.antenna_names = uv_obj.antenna_names.tolist() + uv_obj.antenna_diameters = hera_uv.antenna_diameters + uv_obj.history = hera_uv.history + + uvfits_required_extra = [ + "_antenna_positions", + "_gst0", + "_rdate", + "_earth_omega", + "_dut1", + "_timesys", + ] + for attr in uvfits_required_extra: + param = getattr(uv_obj, attr) + if param.value is None: + param.value = param.spoof_val + setattr(uv_obj, attr, param) - # Tasks are not ordered in UVTask lists, so need to sort them. - uvtask_list = sorted(uvtask_list) - expected_uvtask_list = sorted(expected_uvtask_list) - assert uvtask_list == expected_uvtask_list + assert new_beam_dict == beam_dict + assert new_beam_list == beam_list + assert uv_obj == hera_uv def test_tele_parser(): @@ -654,7 +701,7 @@ def test_param_select_redundant(): @pytest.mark.parametrize('case', np.arange(6)) -def test_uvdata_keyword_init(case): +def test_uvdata_keyword_init(case, tmpdir): base_kwargs = { "antenna_layout_filepath": os.path.join(SIM_DATA_PATH, "test_config/triangle_bl_layout.csv"), @@ -728,12 +775,12 @@ def test_uvdata_keyword_init(case): new_kwargs['output_layout_filename'] = layout_fname new_kwargs['output_yaml_filename'] = obsparam_fname new_kwargs['array_layout'] = antpos_d - new_kwargs['path_out'] = simtest.TESTDATA_PATH + new_kwargs['path_out'] = str(tmpdir) new_kwargs['write_files'] = True uvd = pyuvsim.simsetup.initialize_uvdata_from_keywords(**new_kwargs) - layout_path = os.path.join(simtest.TESTDATA_PATH, layout_fname) - obsparam_path = os.path.join(simtest.TESTDATA_PATH, obsparam_fname) + layout_path = str(tmpdir.join(layout_fname)) + obsparam_path = str(tmpdir.join(obsparam_fname)) assert os.path.exists(layout_path) assert os.path.exists(obsparam_path) @@ -881,12 +928,12 @@ def test_mock_catalogs(): os.remove(fname) -def test_keyword_param_loop(): +def test_keyword_param_loop(tmpdir): # Check that yaml/csv files made by intialize_uvdata_from_keywords will work # on their own. layout_fname = 'temp_layout_kwdloop.csv' obsparam_fname = 'temp_obsparam_kwdloop.yaml' - path_out = simtest.TESTDATA_PATH + path_out = str(tmpdir) antpos_enu = np.ones(30).reshape((10, 3)) antnums = np.arange(10) antpos_d = dict(zip(antnums, antpos_enu)) @@ -907,7 +954,7 @@ def test_keyword_param_loop(): assert uv2 == uvd -def test_multi_analytic_beams(): +def test_multi_analytic_beams(tmpdir): # Test inline definitions of beam attributes. # eg. (in beam configuration file): # @@ -915,8 +962,8 @@ def test_multi_analytic_beams(): # 0 : airy, diameter=14 # 1 : airy, diameter=20 # 2 : gaussian, sigma=0.5 - par_fname = os.path.join(simtest.TESTDATA_PATH, 'test_teleconfig.yaml') - layout_fname = os.path.join(simtest.TESTDATA_PATH, 'test_layout_5ant.csv') + par_fname = str(tmpdir.join('test_teleconfig.yaml')) + layout_fname = str(tmpdir.join('test_layout_5ant.csv')) telescope_location = (-30.72152777777791, 21.428305555555557, 1073.0000000093132) telescope_name = 'SKA' @@ -945,7 +992,7 @@ def test_multi_analytic_beams(): param_dict = {'telescope_config_name': par_fname, 'array_layout': layout_fname} pdict, beam_list, beam_dict = pyuvsim.simsetup.parse_telescope_params( - param_dict, simtest.TESTDATA_PATH) + param_dict, str(tmpdir)) for i, nm in enumerate(names): bid = beam_ids[i] @@ -975,14 +1022,14 @@ def test_direct_fname(): os.remove("triangle_bl_layout.csv") -def test_beamlist_init(): +def test_beamlist_init_errors(): telescope_config_name = os.path.join(SIM_DATA_PATH, 'bl_lite_mixed.yaml') with open(telescope_config_name, 'r') as yf: telconfig = yaml.safe_load(yf) # The path for beam 0 is invalid, and it's not needed for this test. del telconfig['beam_paths'][0] - beam_ids = np.arange(1, 4) + beam_ids = np.arange(1, 6) bad_conf_0 = copy.deepcopy(telconfig) bad_conf_0['beam_paths'][1] = 1.35 @@ -1006,6 +1053,30 @@ def test_beamlist_init(): assert beam_list.spline_interp_opts is not None +def test_beamlist_init(): + telescope_config_name = os.path.join(SIM_DATA_PATH, 'bl_lite_mixed.yaml') + with open(telescope_config_name, 'r') as yf: + telconfig = yaml.safe_load(yf) + + telconfig['beam_paths'][0] = os.path.join(SIM_DATA_PATH, 'HERA_NicCST.uvbeam') + + beam_list = pyuvsim.simsetup._construct_beam_list(np.arange(6), telconfig) + beam_list.set_obj_mode() + + # How the beam attributes should turn out for this file: + assert isinstance(beam_list[0], UVBeam) + assert beam_list[1].type == 'airy' + assert beam_list[1].diameter == 16 + assert beam_list[2].type == 'gaussian' + assert beam_list[2].sigma == 0.03 + assert beam_list[3].type == 'airy' + assert beam_list[3].diameter == 12 + assert beam_list[4].type == 'gaussian' + assert beam_list[4].diameter == 14 + assert beam_list[5].type == 'gaussian' + assert beam_list[5].diameter == 12 + + def test_moon_lsts(): # Check that setting lsts for a Moon simulation works as expected. pytest.importorskip('lunarsky') @@ -1060,7 +1131,7 @@ def test_mock_catalog_moon(): assert mmock != emock -@pytest.fixture +@pytest.fixture(scope='module') def cat_with_some_pols(): # Mock catalog with a couple sources polarized. Nsrcs = 30 @@ -1082,7 +1153,7 @@ def cat_with_some_pols(): name=np.arange(Nsrcs).astype(str), ra=ra, dec=dec, - stokes=stokes, + stokes=stokes * units.Jy, spectral_type='full', freq_array=freqs ) @@ -1094,19 +1165,19 @@ def cat_with_some_pols(): def test_skymodeldata_with_quantity_stokes(unit, cat_with_some_pols): # Support for upcoming pyradiosky change setting SkyModel.stokes # to an astropy Quantity. - if unit == 'K': - pytest.importorskip('astropy_healpix') - path = os.path.join(SKY_DATA_PATH, 'healpix_disk.hdf5') - sky = pyradiosky.SkyModel() - sky.read_healpix_hdf5(path) - else: + if unit == 'Jy': sky = cat_with_some_pols - + else: + pytest.importorskip('analytic_diffuse') + sky, _ = pyuvsim.simsetup.create_mock_catalog( + Time.now(), arrangement='diffuse', diffuse_model='monopole', map_nside=16 + ) if not isinstance(sky.stokes, units.Quantity): sky.stokes *= units.Unit(unit) smd = pyuvsim.simsetup.SkyModelData(sky) assert np.all(sky.stokes.to_value(unit)[0] == smd.stokes_I) + assert smd.flux_unit == unit sky2 = smd.get_skymodel() if units.Quantity != pyradiosky.SkyModel()._stokes.expected_type: @@ -1187,3 +1258,19 @@ def test_skymodeldata_attr_bases(inds, cat_with_some_pols): assert smd_copy.ra.base is smd.ra.base assert smd_copy.dec.base is smd.dec.base assert smd_copy.stokes_I.base is smd.stokes_I.base + + +def test_set_lsts_errors(): + # Error cases on set_lsts function. + uv0 = UVData() + uv0.read_uvfits(longbl_uvfits_file) + uv0.lst_array = None + + uv0.extra_keywords['world'] = 'moon' + if not pyuvsim.astropy_interface.hasmoon: + with pytest.raises(ValueError, match="Cannot construct lsts for MoonLocation"): + pyuvsim.simsetup._set_lsts_on_uvdata(uv0) + + uv0.extra_keywords['world'] = 'tatooine' + with pytest.raises(ValueError, match="Invalid world tatooine."): + pyuvsim.simsetup._set_lsts_on_uvdata(uv0) diff --git a/pyuvsim/tests/test_telescope.py b/pyuvsim/tests/test_telescope.py index 36b2126e9..c4ae08251 100644 --- a/pyuvsim/tests/test_telescope.py +++ b/pyuvsim/tests/test_telescope.py @@ -37,6 +37,7 @@ def beam_objs(): return beams +@pytest.mark.filterwarnings('ignore:Achromatic gaussian') def test_convert_loop(beam_objs): beams = beam_objs beams[0].freq_interp_kind = 'linear' @@ -137,6 +138,19 @@ def test_string_mode(beam_objs): assert False +@pytest.mark.filterwarnings('ignore:Achromatic gaussian') +def test_comparison(beam_objs): + beamlist = pyuvsim.BeamList(beam_objs) + beamlist.set_str_mode() + + beamlist2 = pyuvsim.BeamList(beamlist._str_beam_list) + assert beamlist == beamlist2 + + beamlist.set_obj_mode() + beamlist2.set_obj_mode() + assert beamlist == beamlist2 + + def test_no_overwrite(beam_objs): # Ensure UVBeam keywords are not overwritten by BeamList.uvb_params # while in object mode. diff --git a/pyuvsim/tests/test_utils.py b/pyuvsim/tests/test_utils.py index eae9f8323..2e9dc939c 100644 --- a/pyuvsim/tests/test_utils.py +++ b/pyuvsim/tests/test_utils.py @@ -3,13 +3,11 @@ # Licensed under the 3-clause BSD License import os -import shutil import numpy as np import pytest from pyuvdata import UVData -import pyuvsim.tests as simtest from pyuvsim import utils as simutils from pyuvsim.data import DATA_PATH as SIM_DATA_PATH @@ -89,50 +87,32 @@ def test_altaz_za_az_errors(): simutils.zenithangle_azimuth_to_altaz(0, [0, np.pi / 2]) -def test_file_namer(): +@pytest.mark.parametrize('ext', ['.ext', '.uvfits', '.uvh5', '.yaml', '']) +def test_file_namer(tmpdir, ext): """ - File name incrementer utility + File name incrementer utility, with extensions. """ fnames = [] for i in range(111): - fname = os.path.join(simtest.TESTDATA_PATH, 'file_' + str(i)) + fname = str(tmpdir.join(f"file_{i}{ext}")) with open(fname, 'w') as f: f.write(' ') fnames.append(fname) existing_file = fnames[0] - new_filepath = simutils.check_file_exists_and_increment(existing_file) - for fn in fnames: - os.remove(fn) - assert new_filepath.endswith("_111") - - -def test_file_namer_extensions(): - """ - File name incrementer with specified extension - """ - os.mkdir(os.path.join(simtest.TESTDATA_PATH, 'tempfiles')) - fnames = [] - for i in range(111): - fname = os.path.join(simtest.TESTDATA_PATH, 'tempfiles', 'file_' + str(i) + '.ext') - with open(fname, 'w') as f: - f.write(' ') - fnames.append(fname) - existing_file = fnames[0] - new_filepath = simutils.check_file_exists_and_increment(existing_file, 'ext') - shutil.rmtree(os.path.join(simtest.TESTDATA_PATH, 'tempfiles')) - assert new_filepath.endswith("_111.ext") + if ext == '.ext': + new_filepath = simutils.check_file_exists_and_increment(existing_file, 'ext') + else: + new_filepath = simutils.check_file_exists_and_increment(existing_file) + assert new_filepath.endswith(f"_111{ext}") @pytest.mark.parametrize("save_format", [None, 'uvfits', 'miriad', 'uvh5']) -def test_write_uvdata(save_format): +def test_write_uvdata(save_format, tmpdir): """ Test function that defines filenames from parameter dict """ uv = UVData() - with pytest.warns(UserWarning) as telwarn: - uv.read_uvfits(triangle_uvfits_file) - assert str(telwarn.pop().message).startswith('Telescope 28m_triangle_10time_10chan.yaml ' - 'is not in known_telescopes.') + uv.read_uvfits(triangle_uvfits_file) - ofname = os.path.join(simtest.TESTDATA_PATH, 'test_file') + ofname = str(tmpdir.join('test_file')) filing_dict = {'outfile_name': ofname} expected_ofname = simutils.write_uvdata(uv, filing_dict, return_filename=True, @@ -141,39 +121,30 @@ def test_write_uvdata(save_format): if save_format == 'uvfits' or save_format is None: assert ofname + '.uvfits' == expected_ofname - os.remove(ofname + '.uvfits') elif save_format == 'uvh5': assert ofname + '.uvh5' == expected_ofname - os.remove(ofname + '.uvh5') else: assert ofname == expected_ofname - shutil.rmtree(ofname) -def test_write_error_with_no_format(): +def test_write_error_with_no_format(tmpdir): """Test write_uvdata will error if no format is given.""" uv = UVData() - with pytest.warns(UserWarning) as telwarn: - uv.read_uvfits(triangle_uvfits_file) - assert str(telwarn.pop().message).startswith('Telescope 28m_triangle_10time_10chan.yaml is ' - 'not in known_telescopes.') + uv.read_uvfits(triangle_uvfits_file) - ofname = os.path.join(simtest.TESTDATA_PATH, 'test_file') + ofname = str(tmpdir.join('test_file')) filing_dict = {'outfile_name': ofname} with pytest.raises(ValueError, match='Invalid output format. Options are'): simutils.write_uvdata(uv, filing_dict, return_filename=True, out_format='') -def test_file_format_in_filing_dict(): +def test_file_format_in_filing_dict(tmpdir): """Test file is written out when output_format is set in filing dict.""" uv = UVData() - with pytest.warns(UserWarning) as telwarn: - uv.read_uvfits(triangle_uvfits_file) - assert str(telwarn.pop().message).startswith('Telescope 28m_triangle_10time_10chan.yaml is ' - 'not in known_telescopes.') + uv.read_uvfits(triangle_uvfits_file) - ofname = os.path.join(simtest.TESTDATA_PATH, 'test_file') + ofname = str(tmpdir.join('test_file')) filing_dict = {'outfile_name': ofname} filing_dict['output_format'] = 'uvfits' expected_ofname = simutils.write_uvdata(uv, filing_dict, return_filename=True) diff --git a/pyuvsim/tests/test_uvsim.py b/pyuvsim/tests/test_uvsim.py index 3d827c23c..aac27c72e 100644 --- a/pyuvsim/tests/test_uvsim.py +++ b/pyuvsim/tests/test_uvsim.py @@ -2,32 +2,60 @@ # Copyright (c) 2018 Radio Astronomy Software Group # Licensed under the 3-clause BSD License +import numpy as np import copy import itertools import os +import warnings import astropy.constants as const -import numpy as np import pyuvdata.utils as uvutils from astropy import units from astropy.coordinates import Angle, SkyCoord, EarthLocation, Longitude, Latitude import pytest -from pyuvdata import UVData -from pyuvdata.data import DATA_PATH +from pyuvdata import UVData, UVBeam import pyradiosky import pyuvsim -import pyuvsim.tests as simtest import pyuvsim.utils as simutils from pyuvsim.data import DATA_PATH as SIM_DATA_PATH from pyuvsim.astropy_interface import Time -cst_files = ['HERA_NicCST_150MHz.txt', 'HERA_NicCST_123MHz.txt'] -beam_files = [os.path.join(DATA_PATH, 'NicCSTbeams', f) for f in cst_files] -hera_miriad_file = os.path.join(DATA_PATH, 'hera_testfile') EW_uvfits_file = os.path.join(SIM_DATA_PATH, '28mEWbl_1time_1chan.uvfits') EW_uvfits_10time10chan = os.path.join(SIM_DATA_PATH, '28mEWbl_10time_10chan.uvfits') longbl_uvfits_file = os.path.join(SIM_DATA_PATH, '5km_triangle_1time_1chan.uvfits') +herabeam_default = os.path.join(SIM_DATA_PATH, 'HERA_NicCST.uvbeam') + + +def multi_beams(): + beam0 = UVBeam() + beam0.read_beamfits(herabeam_default) + beam0.freq_interp_kind = 'cubic' + beam0.interpolation_function = 'az_za_simple' + beam1 = pyuvsim.AnalyticBeam('uniform') + with warnings.catch_warnings(): + warnings.simplefilter('ignore') + beam2 = pyuvsim.AnalyticBeam('gaussian', sigma=0.02) + beam3 = pyuvsim.AnalyticBeam('airy', diameter=14.6) + beam_list = pyuvsim.BeamList([beam0, beam1, beam2, beam3]) + + return beam_list + + +multi_beams = multi_beams() + + +@pytest.fixture(scope='module') +def triangle_pos(): + hera_uv = UVData() + hera_uv.read_uvfits(longbl_uvfits_file, + ant_str='cross') # consists of a right triangle of baselines with w term + hera_uv.unphase_to_drift(use_ant_pos=True) + + enu = hera_uv.get_ENU_antpos()[0] + uvw = hera_uv.uvw_array[:hera_uv.Nbls] + + return enu, uvw @pytest.fixture @@ -65,18 +93,18 @@ def uvobj_beams_srcs(): return uv_obj, beam_list, beam_dict, sources -def test_visibility_single_zenith_source(): +def test_visibility_single_zenith_source(cst_beam, hera_loc): """Test single zenith source.""" - beam0 = simtest.make_cst_beams() + beam0 = cst_beam beam1 = pyuvsim.AnalyticBeam('uniform') beam2 = pyuvsim.AnalyticBeam('gaussian', sigma=np.radians(10.0)) beam3 = pyuvsim.AnalyticBeam('airy', diameter=14.0) time = Time('2018-03-01 00:00:00', scale='utc') - array_location = EarthLocation(lat='-30d43m17.5s', lon='21d25m41.9s', - height=1073.) + array_location = hera_loc + time.location = array_location freq = (150e6 * units.Hz) @@ -100,11 +128,10 @@ def test_visibility_single_zenith_source(): assert np.allclose(visibility, np.array([.5, .5, 0, 0]), atol=5e-3) -def test_visibility_source_below_horizon(): +def test_visibility_source_below_horizon(cst_beam, hera_loc): time = Time('2018-03-01 00:00:00', scale='utc') - array_location = EarthLocation(lat='-30d43m17.5s', lon='21d25m41.9s', - height=1073.) + array_location = hera_loc time.location = array_location freq = (150e6 * units.Hz) @@ -118,7 +145,7 @@ def test_visibility_source_below_horizon(): baseline = pyuvsim.Baseline(antenna1, antenna2) - beam = simtest.make_cst_beams() + beam = cst_beam beam_list = pyuvsim.BeamList([beam]) array = pyuvsim.Telescope('telescope_name', array_location, beam_list) @@ -132,13 +159,12 @@ def test_visibility_source_below_horizon(): assert np.allclose(visibility, np.array([0, 0, 0, 0])) -def test_visibility_source_below_horizon_radec(): +def test_visibility_source_below_horizon_radec(cst_beam, hera_loc): # redo with RA/Dec defined source time = Time(2458098.27471265, format='jd') - array_location = EarthLocation( - lat='-30d43m17.5s', lon='21d25m41.9s', height=1073. - ) + array_location = hera_loc + time.location = array_location freq = (150e6 * units.Hz) @@ -146,14 +172,14 @@ def test_visibility_source_below_horizon_radec(): obstime=time, frame='icrs', location=array_location) source = pyradiosky.SkyModel('src_down', source_coord.ra, source_coord.dec, - np.array([1.0, 0, 0, 0]).reshape(4, 1), 'flat') + np.array([1.0, 0, 0, 0]).reshape(4, 1) * units.Jy, 'flat') antenna1 = pyuvsim.Antenna('ant1', 1, np.array([0, 0, 0]), 0) antenna2 = pyuvsim.Antenna('ant2', 2, np.array([107, 0, 0]), 0) baseline = pyuvsim.Baseline(antenna1, antenna2) - beam = simtest.make_cst_beams() + beam = cst_beam beam_list = pyuvsim.BeamList([beam]) array = pyuvsim.Telescope('telescope_name', array_location, beam_list) @@ -167,74 +193,29 @@ def test_visibility_source_below_horizon_radec(): assert np.allclose(visibility, np.array([0, 0, 0, 0])) -def test_visibility_single_zenith_source_uvdata(): - """Test single zenith source using test uvdata file.""" - hera_uv = UVData() - hera_uv.read_uvfits(EW_uvfits_file) - - time = Time(hera_uv.time_array[0], scale='utc', format='jd') - array_location = EarthLocation.from_geocentric( - *hera_uv.telescope_location, unit='m' - ) - freq = hera_uv.freq_array[0, 0] * units.Hz - - # get antennas positions into ENU - antpos = hera_uv.antenna_positions[0:2, :] + hera_uv.telescope_location - lat, lon, alt = hera_uv.telescope_location_lat_lon_alt - antpos = uvutils.ENU_from_ECEF(antpos, lat, lon, alt) - - antenna1 = pyuvsim.Antenna('ant1', 1, np.array(antpos[0, :]), 0) - antenna2 = pyuvsim.Antenna('ant2', 2, np.array(antpos[1, :]), 0) - - # setup the things that don't come from pyuvdata: - # make a source at zenith - time.location = array_location - source, _ = pyuvsim.create_mock_catalog(time, arrangement='zenith') - - beam = simtest.make_cst_beams(freqs=[100e6, 123e6]) - beam_list = pyuvsim.BeamList([beam]) - - baseline = pyuvsim.Baseline(antenna1, antenna2) - array = pyuvsim.Telescope('telescope_name', array_location, beam_list) - task = pyuvsim.UVTask(source, time, freq, baseline, array) - engine = pyuvsim.UVEngine(task) - - visibility = engine.make_visibility() - - assert np.allclose(visibility, np.array([.5, .5, 0, 0]), atol=5e-3) - - -def test_redundant_baselines(): +def test_redundant_baselines(cst_beam, hera_loc): """Check that two perfectly redundant baselines are truly redundant. """ - hera_uv = UVData() - hera_uv.read_uvfits(EW_uvfits_file, ant_str='cross') - hera_uv.unphase_to_drift() - + time = Time(2458098.27471265, format='jd') + array_location = hera_loc + time.location = array_location + freq = (150e6 * units.Hz) src_alt = Angle('85.0d') - time = Time(hera_uv.time_array[0], scale='utc', format='jd') - array_location = EarthLocation.from_geocentric( - *hera_uv.telescope_location, unit='m' - ) - freq = hera_uv.freq_array[0, 0] * units.Hz - - # get antennas positions into ENU - antpos, _ = hera_uv.get_ENU_antpos() + # Set up antenna positions in ENU: + antpos = np.array([[0, 0, 0], [28, 0, 0]], dtype=float) en_shift = [5., 5., 0] - antenna1 = pyuvsim.Antenna('ant1', 1, np.array(antpos[0, :]), 0) - antenna2 = pyuvsim.Antenna('ant2', 2, np.array(antpos[1, :]), 0) - antenna3 = pyuvsim.Antenna('ant3', 3, np.array(antpos[0, :]) + en_shift, 0) - antenna4 = pyuvsim.Antenna('ant4', 4, np.array(antpos[1, :]) + en_shift, 0) + antenna1 = pyuvsim.Antenna('ant1', 1, antpos[0, :], 0) + antenna2 = pyuvsim.Antenna('ant2', 2, antpos[1, :], 0) + antenna3 = pyuvsim.Antenna('ant3', 3, antpos[0, :] + en_shift, 0) + antenna4 = pyuvsim.Antenna('ant4', 4, antpos[1, :] + en_shift, 0) - # setup the things that don't come from pyuvdata: # make a source off zenith time.location = array_location source, _ = pyuvsim.create_mock_catalog(time, arrangement='off-zenith', alt=src_alt.deg) - beam = simtest.make_cst_beams(freqs=[100e6, 123e6]) - beam_list = pyuvsim.BeamList([beam]) + beam_list = pyuvsim.BeamList([cst_beam]) baseline1 = pyuvsim.Baseline(antenna1, antenna2) baseline2 = pyuvsim.Baseline(antenna3, antenna4) @@ -254,11 +235,14 @@ def test_redundant_baselines(): assert np.allclose(visibility1, visibility2) -def test_single_offzenith_source_uvfits(): - """Test single off-zenith source using test uvdata file.""" - hera_uv = UVData() - hera_uv.read_uvfits(EW_uvfits_file, ant_str='cross') - hera_uv.unphase_to_drift() +@pytest.mark.parametrize('beam', multi_beams) +def test_single_offzenith_source(beam, hera_loc): + """Test single off-zenith source.""" + + time = Time(2458098.27471265, format='jd') + array_location = hera_loc + time.location = array_location + freq = (123e6 * units.Hz) src_az = Angle('90.0d') src_alt = Angle('85.0d') @@ -268,14 +252,7 @@ def test_single_offzenith_source_uvfits(): src_m = np.cos(src_az.rad) * np.sin(src_za.rad) src_n = np.cos(src_za.rad) - time = Time(hera_uv.time_array[0], scale='utc', format='jd') - array_location = EarthLocation.from_geocentric( - *hera_uv.telescope_location, unit='m' - ) - freq = hera_uv.freq_array[0, 0] * units.Hz - - # get antennas positions into ENU - antpos, _ = hera_uv.get_ENU_antpos() + antpos = np.array([[0, 0, 0], [28, 0, 0]], dtype=float) antenna1 = pyuvsim.Antenna('ant1', 1, np.array(antpos[0, :]), 0) antenna2 = pyuvsim.Antenna('ant2', 2, np.array(antpos[1, :]), 0) @@ -296,7 +273,6 @@ def test_single_offzenith_source_uvfits(): assert np.isclose(src_lmn[1], src_m) assert np.isclose(src_lmn[2], src_n) - beam = simtest.make_cst_beams(freqs=[100e6, 123e6]) beam_list = pyuvsim.BeamList([beam]) baseline = pyuvsim.Baseline(antenna1, antenna2) @@ -308,7 +284,6 @@ def test_single_offzenith_source_uvfits(): visibility = engine.make_visibility() # analytically calculate visibility - beam.peak_normalize() beam.interpolation_function = 'az_za_simple' beam_za, beam_az = simutils.altaz_to_zenithangle_azimuth(src_alt.rad, src_az.rad) beam_za2, beam_az2 = simutils.altaz_to_zenithangle_azimuth(src_alt_az[0], src_alt_az[1]) @@ -329,7 +304,8 @@ def test_single_offzenith_source_uvfits(): beam_jones = antenna1.get_beam_jones(array, src_alt_az, freq) assert np.allclose(beam_jones, jones) - uvw_wavelength_array = hera_uv.uvw_array * units.m / const.c * freq.to('1/s') + uvw_array = units.Quantity([[28., 0., 0.]], unit='m') + uvw_wavelength_array = uvw_array / const.c * freq.to('1/s') # Remove source axis from jones matrix jones = jones.squeeze() vis_analytic = 0.5 * np.dot(jones, np.conj(jones).T) * np.exp( @@ -343,19 +319,16 @@ def test_single_offzenith_source_uvfits(): [vis_analytic[0, 0], vis_analytic[1, 1], vis_analytic[0, 1], vis_analytic[1, 0]] ) - assert np.allclose(baseline.uvw.to_value('m'), hera_uv.uvw_array[0:hera_uv.Nbls], atol=1e-4) - assert np.allclose(visibility, vis_analytic, atol=1e-4) + assert np.allclose(baseline.uvw.to_value('m'), uvw_array.value) + assert np.allclose(visibility, vis_analytic) -def test_offzenith_source_multibl_uvfits(): - """Test single off-zenith source using test uvdata file. - Calculate visibilities for a baseline triangle. - """ - hera_uv = UVData() - hera_uv.read_uvfits(longbl_uvfits_file, - ant_str='cross') # consists of a right triangle of baselines with w term - hera_uv.unphase_to_drift(use_ant_pos=True) +@pytest.mark.parametrize('beam', multi_beams) +def test_offzenith_source_multibl(beam, hera_loc, triangle_pos): + """Calculate visibilities for a baseline triangle of an off-zenith source.""" + enu_antpos, uvw_array = triangle_pos + time = Time(2458098.27471265, format='jd') src_az = Angle('90.0d') src_alt = Angle('85.0d') src_za = Angle('90.0d') - src_alt @@ -364,27 +337,18 @@ def test_offzenith_source_multibl_uvfits(): src_m = np.cos(src_az.rad) * np.sin(src_za.rad) src_n = np.cos(src_za.rad) - time = Time(hera_uv.time_array[0], scale='utc', format='jd') - array_location = EarthLocation.from_geocentric( - *hera_uv.telescope_location, unit='m' - ) - freq = hera_uv.freq_array[0, 0] * units.Hz + array_location = hera_loc + freq = (123e6 * units.Hz) - # get antennas positions into ENU - antpos, ants = hera_uv.get_ENU_antpos() - antenna1 = pyuvsim.Antenna('ant1', ants[0], np.array(antpos[0, :]), 0) - antenna2 = pyuvsim.Antenna('ant2', ants[1], np.array(antpos[1, :]), 0) - antenna3 = pyuvsim.Antenna('ant3', ants[2], np.array(antpos[2, :]), 0) + antenna1 = pyuvsim.Antenna('ant1', 0, np.array(enu_antpos[0, :]), 0) + antenna2 = pyuvsim.Antenna('ant2', 1, np.array(enu_antpos[1, :]), 0) + antenna3 = pyuvsim.Antenna('ant3', 2, np.array(enu_antpos[2, :]), 0) - # setup the things that don't come from pyuvdata: # make a source off zenith time.location = array_location - # create_mock_catalog uses azimuth of 90 source, _ = pyuvsim.create_mock_catalog(time, arrangement='off-zenith', alt=src_alt.deg) - beam = pyuvsim.AnalyticBeam('uniform') - - beam_list = [beam] + beam_list = pyuvsim.BeamList([beam]) baselines = [pyuvsim.Baseline(antenna2, antenna1), pyuvsim.Baseline(antenna3, antenna1), @@ -393,49 +357,45 @@ def test_offzenith_source_multibl_uvfits(): tasks = [pyuvsim.UVTask(source, time, freq, bl, array) for bl in baselines] visibilities = [] uvws = [] + engine = pyuvsim.UVEngine() for t in tasks: - engine = pyuvsim.UVEngine(t) + engine.set_task(t) visibilities.append(engine.make_visibility()) uvws.append(t.baseline.uvw) - uvws = np.array(uvws) - # analytically calculate visibilities + # analytically calculate visibilities beam.peak_normalize() beam.interpolation_function = 'az_za_simple' interpolated_beam, interp_basis_vector = beam.interp( - az_array=np.array([src_az.rad]), za_array=np.array([src_za.rad]), + az_array=np.array([0.0]), za_array=np.array([src_za.rad]), freq_array=np.array([freq.to_value('Hz')]) ) - jones = np.zeros((2, 2), dtype=np.complex64) - jones[0, 0] = interpolated_beam[1, 0, 0, 0, 0] - jones[1, 1] = interpolated_beam[0, 0, 1, 0, 0] - jones[1, 0] = interpolated_beam[1, 0, 1, 0, 0] - jones[0, 1] = interpolated_beam[0, 0, 0, 0, 0] + jones = np.zeros((2, 2, 1), dtype=np.complex64) + jones[0, 0] = interpolated_beam[1, 0, 0, 0, :] + jones[1, 1] = interpolated_beam[0, 0, 1, 0, :] + jones[1, 0] = interpolated_beam[1, 0, 1, 0, :] + jones[0, 1] = interpolated_beam[0, 0, 0, 0, :] - uvw_wavelength_array = hera_uv.uvw_array[0:hera_uv.Nbls] * units.m / const.c * freq.to('1/s') + uvw_wavelength_array = uvw_array * units.m / const.c * freq.to('1/s') visibilities_analytic = [] for u, v, w in uvw_wavelength_array: - vis = 0.5 * np.dot(jones, np.conj(jones).T) * np.exp( + vis = 0.5 * np.dot(jones[..., 0], np.conj(jones[..., 0]).T) * np.exp( 2j * np.pi * (u * src_l + v * src_m + w * src_n)) - visibilities_analytic.append(np.array([vis[0, 0], vis[1, 1], vis[1, 0], vis[0, 1]])) - - # the file used different phasing code than the test uses -- increase the tolerance - assert np.allclose(uvws, hera_uv.uvw_array[0:hera_uv.Nbls], atol=1e-4) + visibilities_analytic.append(np.array([vis[0, 0], vis[1, 1], vis[0, 1], vis[1, 0]])) - # the file used different phasing code than the test uses -- increase the tolerance - assert np.allclose(visibilities, visibilities_analytic, atol=1e-4) + assert np.allclose(uvws, uvw_array) + assert np.allclose(visibilities, visibilities_analytic) -def test_file_to_tasks(): +def test_file_to_tasks(cst_beam): hera_uv = UVData() hera_uv.read_uvfits(EW_uvfits_file) time = Time(hera_uv.time_array[0], scale='utc', format='jd') sources, _ = pyuvsim.create_mock_catalog(time, arrangement='zenith', Nsrcs=5, return_data=True) - beam = simtest.make_cst_beams() - beam_list = pyuvsim.BeamList([beam]) + beam_list = pyuvsim.BeamList([cst_beam]) Nblts = hera_uv.Nblts Nfreqs = hera_uv.Nfreqs @@ -514,8 +474,7 @@ def test_gather(): time = Time(hera_uv.time_array[0], scale='utc', format='jd') sources, _ = pyuvsim.create_mock_catalog(time, arrangement='zenith', return_data=True) - beam = simtest.make_cst_beams(freqs=[100e6, 123e6]) - beam_list = pyuvsim.BeamList([beam]) + beam_list = pyuvsim.BeamList([multi_beams[1]]) Nblts = hera_uv.Nblts Nfreqs = hera_uv.Nfreqs @@ -546,8 +505,7 @@ def test_local_task_gen(): time, arrangement='random', Nsrcs=5, return_data=True ) - beam = simtest.make_cst_beams(freqs=[100e6, 123e6]) - beam_list = pyuvsim.BeamList([beam]) + beam_list = pyuvsim.BeamList([multi_beams[1]]) Nblts = hera_uv.Nblts Nfreqs = hera_uv.Nfreqs @@ -585,19 +543,6 @@ def test_local_task_gen(): assert np.allclose(engine1.make_visibility(), engine0.make_visibility()) -def test_pol_error(): - # Check that running with a uvdata object without the proper polarizations will fail. - pytest.importorskip('mpi4py') - - hera_uv = UVData() - hera_uv.read_uvfits(EW_uvfits_file) - - hera_uv.select(polarizations=['xx']) - - with pytest.raises(ValueError, match='input_uv must have XX,YY,XY,YX polarization'): - pyuvsim.run_uvdata_uvsim(hera_uv, ['beamlist']) - - def test_task_coverage(): """ Check that the task ids generated in different scenarios @@ -720,35 +665,6 @@ def test_source_splitting(): pyuvsim.mpi.Npus_node = 1 -def test_get_beam_jones(): - # check setting the interpolation method - - array_location = EarthLocation(lat='-30d43m17.5s', lon='21d25m41.9s', height=1073.) - - beam = simtest.make_cst_beams(freqs=[100e6, 123e6]) - beam.freq_interp_kind = None - beam.interpolation_function = 'az_za_simple' - beam_list = pyuvsim.BeamList([beam]) - antenna1 = pyuvsim.Antenna('ant1', 1, np.array([0, 10, 0]), 0) - array = pyuvsim.Telescope('telescope_name', array_location, beam_list) - source_altaz = np.array([[0.0], [np.pi / 4.]]) - freq = 100e6 * units.Hz - - with pytest.raises(ValueError, match='freq_interp_kind must be set'): - antenna1.get_beam_jones(array, source_altaz, freq) - - jones = antenna1.get_beam_jones(array, source_altaz, freq, freq_interp_kind='cubic') - assert beam.freq_interp_kind == 'cubic' - jones0 = antenna1.get_beam_jones(array, source_altaz, freq) - jones1 = antenna1.get_beam_jones(array, source_altaz, freq, freq_interp_kind='linear') - assert beam.freq_interp_kind == 'linear' - jones2 = antenna1.get_beam_jones(array, source_altaz, freq) - - assert (np.all(jones2 == jones0) - and np.all(jones1 == jones) - and np.all(jones1 == jones0)) - - def test_quantity_reuse(uvobj_beams_srcs): # Check that the right quantities on the UVEngine are changed when # the time/frequency/antenna pair change. @@ -887,8 +803,8 @@ def test_fullfreq_check(uvobj_beams_srcs): freqs0 = np.linspace(100, 130, uv_obj.Nfreqs) * 1e6 * units.Hz freqs1 = uv_obj.freq_array[0, :] * units.Hz - stokes = np.zeros((4, uv_obj.Nfreqs, Nsrcs)) - stokes[0, :, :] = 1.0 + stokes = np.zeros((4, uv_obj.Nfreqs, Nsrcs)) * units.Jy + stokes[0, :, :] = 1.0 * units.Jy ra = Longitude(np.linspace(0, 2 * np.pi, Nsrcs), 'rad') dec = Latitude(np.linspace(-np.pi / 2, np.pi / 3, Nsrcs), 'rad') @@ -944,3 +860,15 @@ def test_moonloc_error(uvobj_beams_srcs): with pytest.raises(ValueError, match="If world keyword is set, it must "): next(pyuvsim.uvsim.uvdata_to_task_iter(range(5), uv_obj, sources, beam_list, beam_dict)) + + +def test_run_mpierr(): + params = pyuvsim.simsetup._config_str_to_dict( + os.path.join(SIM_DATA_PATH, 'test_config', 'param_1time_1src_testcat.yaml') + ) + if pyuvsim.mpi is None: + with pytest.raises(ImportError, match='You need mpi4py to use the uvsim module'): + pyuvsim.run_uvsim(params, return_uv=True) + + with pytest.raises(ImportError, match='You need mpi4py to use the uvsim module'): + pyuvsim.run_uvdata_uvsim(UVData(), ['beamlist']) diff --git a/pyuvsim/tests/test_z_profiler.py b/pyuvsim/tests/test_z_profiler.py deleted file mode 100644 index 81f10b5ad..000000000 --- a/pyuvsim/tests/test_z_profiler.py +++ /dev/null @@ -1,49 +0,0 @@ -# -*- mode: python; coding: utf-8 -* -# Copyright (c) 2018 Radio Astronomy Software Group -# Licensed under the 3-clause BSD License - -from numpy import unique -import os -import shutil -import atexit -import pytest - -import pyuvsim -from pyuvsim.data import DATA_PATH as SIM_DATA_PATH - -try: - from line_profiler import LineProfiler -except ImportError: - LineProfiler = False - - -def profdata_dir_setup(): - # The pytest temporary test directory doesn't work with the profiling outputs - # because they are only created after the exit functions are run, and pytest - # deletes its temporary directory before then. - # As a workaround, make a temporary directory for profiling data and register - # its deletion in the exit functions. - outpath = os.path.join(SIM_DATA_PATH, 'temporary_profiling_data') - os.mkdir(outpath) - atexit.register(shutil.rmtree, outpath) - return outpath - - -def test_profiler(): - if LineProfiler: - outpath = profdata_dir_setup() - testprof_fname = os.path.join(outpath, 'time_profile.out') - pyuvsim.profiling.set_profiler(outfile_prefix=testprof_fname, dump_raw=True) - with pytest.warns(UserWarning, match='Profiler already set'): - pyuvsim.profiling.set_profiler(outfile_prefix=testprof_fname[:-4], dump_raw=True) - param_filename = os.path.join(SIM_DATA_PATH, 'test_config', 'param_1time_1src_testcat.yaml') - pyuvsim.uvsim.run_uvsim(param_filename, return_uv=True) - time_profiler = pyuvsim.profiling.get_profiler() - time_profiler.disable_by_count() - assert isinstance(time_profiler, LineProfiler) - assert hasattr(time_profiler, 'rank') - assert hasattr(time_profiler, 'meta_file') - lstats = time_profiler.get_stats() - assert len(lstats.timings) != 0 - func_names = [k[2] for k in lstats.timings.keys()] - assert unique(func_names).tolist() == sorted(pyuvsim.profiling.default_profile_funcs) diff --git a/pyuvsim/utils.py b/pyuvsim/utils.py index 481e893c1..5e004176e 100644 --- a/pyuvsim/utils.py +++ b/pyuvsim/utils.py @@ -76,16 +76,21 @@ def altaz_to_zenithangle_azimuth(altitude, azimuth): """ Convert from astropy altaz convention to UVBeam az/za convention. - Args: - altitude: in radians - azimuth: in radians in astropy convention: East of North (N=0, E=pi/2) + Parameters + ---------- + altitude, azimuth: float or array of float + altitude above horizon + azimuth in radians in astropy convention: East of North (N=0, E=pi/2) - Returns: - zenith_angle in radians - azimuth in radians in uvbeam convention: North of East(East=0, North=pi/2) + Returns + ------- + zenith_angle: float or array of float + In radians + azimuth: float or array of float + In radians in uvbeam convention: North of East(East=0, North=pi/2) """ - input_alt = np.array(altitude) - input_az = np.array(azimuth) + input_alt = np.asarray(altitude) + input_az = np.asarray(azimuth) if input_alt.size != input_az.size: raise ValueError('number of altitude and azimuth values must match.') @@ -96,7 +101,7 @@ def altaz_to_zenithangle_azimuth(altitude, azimuth): wh_neg = np.where(new_azimuth < -1e-9) if wh_neg[0].size > 0: new_azimuth[wh_neg] = new_azimuth[wh_neg] + np.pi * 2 - else: + elif new_azimuth.size == 1: if new_azimuth < -1e-9: new_azimuth = new_azimuth + np.pi * 2 diff --git a/pyuvsim/uvsim.py b/pyuvsim/uvsim.py index ea8706db5..6c23bba1f 100644 --- a/pyuvsim/uvsim.py +++ b/pyuvsim/uvsim.py @@ -4,17 +4,13 @@ import numpy as np import yaml -import warnings from astropy.coordinates import EarthLocation import astropy.units as units from astropy.units import Quantity from astropy.constants import c as speed_of_light from pyuvdata import UVData -try: - from . import mpi -except ImportError: - mpi = None +from . import mpi from . import simsetup from . import utils as simutils from .antenna import Antenna @@ -138,21 +134,11 @@ def set_task(self, task): def apply_beam(self): """ Set apparent coherency from jones matrices and source coherency. """ - if not self.update_beams: - return - beam1_id, beam2_id = self.current_beam_pair sources = self.task.sources baseline = self.task.baseline - if not hasattr(sources, 'above_horizon'): - warnings.warn("SkyModel class lacks horizon cut on position and coherency calculations." - " This will slow evaluation considerably. Please update your pyradiosky" - " installation to the latest version." - , DeprecationWarning) - setattr(sources, "above_horizon", slice(None)) - if sources.alt_az is None: sources.update_positions(self.task.time, self.task.telescope.location) @@ -275,7 +261,6 @@ def uvdata_to_task_iter(task_ids, input_uv, catalog, beam_list, beam_dict, Nsky_ # Skymodel will now be passed in as a catalog array. if not isinstance(catalog, SkyModelData): raise TypeError("catalog must be a SkyModelData object.") - # Splitting the catalog for memory's sake. Nsrcs_total = catalog.Ncomponents if Nsky_parts > 1: @@ -320,14 +305,17 @@ def uvdata_to_task_iter(task_ids, input_uv, catalog, beam_list, beam_dict, Nsky_ time_array = Time(input_uv.time_array, scale='utc', format='jd', location=telescope.location) for src_i in src_iter: sky = catalog.get_skymodel(src_i) + if ( + sky.spectral_type == 'flat' + and sky.freq_array is None + and sky.reference_frequency is None + ): + sky.freq_array = freq_array[0] if sky.component_type == 'healpix' and hasattr(sky, 'healpix_to_point'): sky.healpix_to_point() - if sky.spectral_type != 'flat' and hasattr(sky, 'at_frequencies'): + if sky.spectral_type != 'flat': sky.at_frequencies(freq_array[0]) - elif sky.spectral_type == 'full': - if not np.allclose(sky.freq_array, freq_array): - raise ValueError("SkyModel spectral type is 'full', and " - "its frequencies do not match simulation frequencies.") + for task_index in task_ids: # Shape indicates slowest to fastest index. if not isinstance(task_index, tuple): diff --git a/setup.py b/setup.py index 9d55939fb..f2d576f0e 100644 --- a/setup.py +++ b/setup.py @@ -29,7 +29,7 @@ 'use_scm_version': {'local_scheme': branch_scheme}, 'include_package_data': True, 'install_requires': ['numpy>=1.15', 'scipy', 'astropy>=4.0', 'pyyaml', - 'pyradiosky', 'pyuvdata', 'setuptools_scm'], + 'pyradiosky>=0.1.0', 'pyuvdata', 'setuptools_scm'], 'test_requires': ['pytest'], 'classifiers': ['Development Status :: 5 - Production/Stable', 'Intended Audience :: Science/Research',