From 533847fed1603c6f40c115dbd38588438c6aa476 Mon Sep 17 00:00:00 2001 From: Simon Perkins Date: Wed, 5 Jul 2017 17:25:11 +0200 Subject: [PATCH 1/8] Compute parallactic angle in astropy Based on Bruce Merry's katsdpimager code. --- montblanc/util/parallactic_angles.py | 104 ++++++++++++++++++++++++++- 1 file changed, 103 insertions(+), 1 deletion(-) diff --git a/montblanc/util/parallactic_angles.py b/montblanc/util/parallactic_angles.py index e17f23e72..ca53615ed 100644 --- a/montblanc/util/parallactic_angles.py +++ b/montblanc/util/parallactic_angles.py @@ -80,4 +80,106 @@ def parallactic_angles(times, antenna_positions, field_centre): pm.posangle(fc_rad, zenith).get_value("rad") for rp in reference_positions ] - for t in times]) \ No newline at end of file + for t in times]) + + +from astropy.coordinates import EarthLocation, SkyCoord, AltAz, CIRS, Angle +from astropy.time import Time +from astropy import units + +def _parallactic_angle_astropy(times, antenna_positions, field_centre): + """ + Computes parallactic angles per timestep for the given + reference antenna position and field centre. + + Arguments: + times: ndarray + Array of unique times with shape (ntime,), + obtained from TIME column of MS table + antenna_positions: ndarray of shape (na, 3) + Antenna positions, obtained from POSITION + column of MS ANTENNA sub-table + field_centre : ndarray of shape (2,) + Field centre, should be obtained from MS PHASE_DIR + + Returns: + An array of parallactic angles per time-step + + """ + from astropy.coordinates import EarthLocation, SkyCoord, AltAz, CIRS + from astropy.time import Time + from astropy import units + + fc = field_centre + ap = antenna_positions + + FT = times.dtype + ntime = times.shape[0] + na = antenna_positions.shape[0] + + # Convert from MJD second to MJD + times = Time(times / 86400.00, format='mjd', scale='utc') + + ap = EarthLocation.from_geocentric(ap[:,0], ap[:,1], ap[:,2], unit='m') + fc = SkyCoord(ra=fc[0], dec=fc[1], unit=units.rad, frame='fk5') + pole = SkyCoord(ra=0, dec=90, unit=units.deg, frame='fk5') + + cirs_frame = CIRS(obstime=times) + pole_cirs = pole.transform_to(cirs_frame) + fc_cirs = fc.transform_to(cirs_frame) + + pa = units.Quantity(np.empty((ntime, na), dtype=FT), + dtype=FT, unit=units.rad, copy=False) + + for i, ant_pos in enumerate(ap): + altaz_frame = AltAz(location=ant_pos, obstime=times) + pole_altaz = pole_cirs.transform_to(altaz_frame) + fc_altaz = fc_cirs.transform_to(altaz_frame) + pa[:,i] = fc_altaz.position_angle(pole_altaz) + + return pa +# print pa + +if __name__ == "__main__": + # 5s second types from 12h00 midday on 1st Feb + times = np.arange('2017-02-01T12:00', '2017-02-01T16:00', dtype='datetime64[5s]') + ftimes = times.astype('datetime64[s]').astype(np.float64) + + # Westerbork antenna positions + antenna_positions = np.array([ + [ 3828763.10544699, 442449.10566454, 5064923.00777 ], + [ 3828746.54957258, 442592.13950824, 5064923.00792 ], + [ 3828729.99081359, 442735.17696417, 5064923.00829 ], + [ 3828713.43109885, 442878.2118934 , 5064923.00436 ], + [ 3828696.86994428, 443021.24917264, 5064923.00397 ], + [ 3828680.31391933, 443164.28596862, 5064923.00035 ], + [ 3828663.75159173, 443307.32138056, 5064923.00204 ], + [ 3828647.19342757, 443450.35604638, 5064923.0023 ], + [ 3828630.63486201, 443593.39226634, 5064922.99755 ], + [ 3828614.07606798, 443736.42941621, 5064923. ], + [ 3828609.94224429, 443772.19450029, 5064922.99868 ], + [ 3828601.66208572, 443843.71178407, 5064922.99963 ], + [ 3828460.92418735, 445059.52053929, 5064922.99071 ], + [ 3828452.64716351, 445131.03744105, 5064922.98793 ]], + dtype=np.float64) + + phase_centre = np.array([ 0. , 1.04719755], dtype=np.float64) + + import time + + t = time.time() + pa_astro = _parallactic_angle_astropy(ftimes, antenna_positions, phase_centre) + print "pa_astropy done in ", time.time() - t + + t = time.time() + pa_casa = parallactic_angles(ftimes, antenna_positions, phase_centre) + print "pa_casa done in ", time.time() - t + + pa_astro = Angle(pa_astro, unit=units.deg).wrap_at(180*units.deg) + pa_casa = Angle(pa_casa*units.rad, unit=units.deg).wrap_at(180*units.deg) + + for a, c in zip(pa_astro.flat, pa_casa.flat): + print a, c, (a-c).wrap_at(180*units.deg) + + + From e0ffa7a11adf2f4351b65108f5247c761482e6bd Mon Sep 17 00:00:00 2001 From: Simon Perkins Date: Wed, 5 Jul 2017 17:25:56 +0200 Subject: [PATCH 2/8] Use AZELGEO frame to match astropy code. But this breaks agreements with MeqTrees... --- montblanc/util/parallactic_angles.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/montblanc/util/parallactic_angles.py b/montblanc/util/parallactic_angles.py index ca53615ed..caafff46f 100644 --- a/montblanc/util/parallactic_angles.py +++ b/montblanc/util/parallactic_angles.py @@ -54,7 +54,7 @@ def parallactic_angles(times, antenna_positions, field_centre): try: # Create direction measure for the zenith - zenith = pm.direction('AZEL','0deg','90deg') + zenith = pm.direction('AZELGEO','0deg','90deg') except AttributeError as e: if pm is None: raise ImportError("python-casacore import failed") From 524a26e70834671a09369f2a56c20521bfe74ce6 Mon Sep 17 00:00:00 2001 From: Simon Perkins Date: Wed, 5 Jul 2017 17:43:45 +0200 Subject: [PATCH 3/8] Simplify astropy parallactic angle code. - Broadcast time and antenna positions against each other. --- montblanc/util/parallactic_angles.py | 29 ++++++++-------------------- 1 file changed, 8 insertions(+), 21 deletions(-) diff --git a/montblanc/util/parallactic_angles.py b/montblanc/util/parallactic_angles.py index caafff46f..51110704d 100644 --- a/montblanc/util/parallactic_angles.py +++ b/montblanc/util/parallactic_angles.py @@ -85,9 +85,10 @@ def parallactic_angles(times, antenna_positions, field_centre): from astropy.coordinates import EarthLocation, SkyCoord, AltAz, CIRS, Angle from astropy.time import Time + from astropy import units -def _parallactic_angle_astropy(times, antenna_positions, field_centre): +def _parallactic_angle_astropy(times, ap, fc): """ Computes parallactic angles per timestep for the given reference antenna position and field centre. @@ -96,10 +97,10 @@ def _parallactic_angle_astropy(times, antenna_positions, field_centre): times: ndarray Array of unique times with shape (ntime,), obtained from TIME column of MS table - antenna_positions: ndarray of shape (na, 3) + ap: ndarray of shape (na, 3) Antenna positions, obtained from POSITION column of MS ANTENNA sub-table - field_centre : ndarray of shape (2,) + fc : ndarray of shape (2,) Field centre, should be obtained from MS PHASE_DIR Returns: @@ -110,13 +111,6 @@ def _parallactic_angle_astropy(times, antenna_positions, field_centre): from astropy.time import Time from astropy import units - fc = field_centre - ap = antenna_positions - - FT = times.dtype - ntime = times.shape[0] - na = antenna_positions.shape[0] - # Convert from MJD second to MJD times = Time(times / 86400.00, format='mjd', scale='utc') @@ -128,17 +122,10 @@ def _parallactic_angle_astropy(times, antenna_positions, field_centre): pole_cirs = pole.transform_to(cirs_frame) fc_cirs = fc.transform_to(cirs_frame) - pa = units.Quantity(np.empty((ntime, na), dtype=FT), - dtype=FT, unit=units.rad, copy=False) - - for i, ant_pos in enumerate(ap): - altaz_frame = AltAz(location=ant_pos, obstime=times) - pole_altaz = pole_cirs.transform_to(altaz_frame) - fc_altaz = fc_cirs.transform_to(altaz_frame) - pa[:,i] = fc_altaz.position_angle(pole_altaz) - - return pa -# print pa + altaz_frame = AltAz(location=ap[None,:], obstime=times[:,None]) + pole_altaz = pole_cirs[:,None].transform_to(altaz_frame) + fc_altaz = fc_cirs[:,None].transform_to(altaz_frame) + return fc_altaz.position_angle(pole_altaz) if __name__ == "__main__": # 5s second types from 12h00 midday on 1st Feb From e0698ec816a4be40c2e85fab73b16d3dd583ffda Mon Sep 17 00:00:00 2001 From: Simon Perkins Date: Thu, 6 Jul 2017 10:45:31 +0200 Subject: [PATCH 4/8] Try not going through the CIRS frame Apparently this happens anyway when going to AltAz. --- montblanc/util/parallactic_angles.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/montblanc/util/parallactic_angles.py b/montblanc/util/parallactic_angles.py index 51110704d..25bad6b0a 100644 --- a/montblanc/util/parallactic_angles.py +++ b/montblanc/util/parallactic_angles.py @@ -118,13 +118,9 @@ def _parallactic_angle_astropy(times, ap, fc): fc = SkyCoord(ra=fc[0], dec=fc[1], unit=units.rad, frame='fk5') pole = SkyCoord(ra=0, dec=90, unit=units.deg, frame='fk5') - cirs_frame = CIRS(obstime=times) - pole_cirs = pole.transform_to(cirs_frame) - fc_cirs = fc.transform_to(cirs_frame) - altaz_frame = AltAz(location=ap[None,:], obstime=times[:,None]) - pole_altaz = pole_cirs[:,None].transform_to(altaz_frame) - fc_altaz = fc_cirs[:,None].transform_to(altaz_frame) + pole_altaz = pole.transform_to(altaz_frame) + fc_altaz = fc.transform_to(altaz_frame) return fc_altaz.position_angle(pole_altaz) if __name__ == "__main__": From a589c10b38980bde789bf84e4675162e7473b689 Mon Sep 17 00:00:00 2001 From: Simon Perkins Date: Thu, 6 Jul 2017 10:46:29 +0200 Subject: [PATCH 5/8] Revert "Try not going through the CIRS frame" This reverts commit 0c0849ca57f59aba933f10042567ea865db5c9ed. This makes things much slower, almost 10x. This is likely because the scalar pole and field centre are broadcast against all antenna positions and times. --- montblanc/util/parallactic_angles.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/montblanc/util/parallactic_angles.py b/montblanc/util/parallactic_angles.py index 25bad6b0a..51110704d 100644 --- a/montblanc/util/parallactic_angles.py +++ b/montblanc/util/parallactic_angles.py @@ -118,9 +118,13 @@ def _parallactic_angle_astropy(times, ap, fc): fc = SkyCoord(ra=fc[0], dec=fc[1], unit=units.rad, frame='fk5') pole = SkyCoord(ra=0, dec=90, unit=units.deg, frame='fk5') + cirs_frame = CIRS(obstime=times) + pole_cirs = pole.transform_to(cirs_frame) + fc_cirs = fc.transform_to(cirs_frame) + altaz_frame = AltAz(location=ap[None,:], obstime=times[:,None]) - pole_altaz = pole.transform_to(altaz_frame) - fc_altaz = fc.transform_to(altaz_frame) + pole_altaz = pole_cirs[:,None].transform_to(altaz_frame) + fc_altaz = fc_cirs[:,None].transform_to(altaz_frame) return fc_altaz.position_angle(pole_altaz) if __name__ == "__main__": From 09049ad211cbe488f9c58005bf57f9e906856fc2 Mon Sep 17 00:00:00 2001 From: Simon Perkins Date: Thu, 14 Sep 2017 11:04:01 +0200 Subject: [PATCH 6/8] Simplify include_pkg_dirs --- setup.py | 20 +++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/setup.py b/setup.py index 82d92135b..2541c5e29 100644 --- a/setup.py +++ b/setup.py @@ -21,6 +21,7 @@ import json import logging import os +from os.path import join as pjoin import sys #============== @@ -30,7 +31,7 @@ from install.install_log import log mb_path = 'montblanc' -mb_inc_path = os.path.join(mb_path, 'include') +mb_inc_path = pjoin(mb_path, 'include') #=================== # Detect readthedocs @@ -142,10 +143,23 @@ def include_pkg_dirs(): # OK, so everything starts with 'montblanc/' # Take everything after that ('include...') and # append a '/*.*' to it - pkg_dirs.append(os.path.join(root[l:], d, '*.*')) + pkg_dirs.append(pjoin(root[l:], d, '*.*')) return pkg_dirs + +def include_pkg_dirs(): + mb_inc_path = pjoin("montblanc", "include") + l = len("montblanc") + len(os.sep) + + exclude = set(['docs', '.git', '.svn']) + + return [pjoin(path, d, '*.*')[l:] for path, dirs, files + in os.walk(mb_inc_path, topdown=True) + for d in dirs if d not in exclude] + + + install_requires = [ 'attrdict >= 2.0.0', 'attrs >= 16.3.0', @@ -196,7 +210,7 @@ def include_pkg_dirs(): from install.versioning import maintain_version setup(name='montblanc', - version=maintain_version(os.path.join('montblanc', 'version.py')), + version=maintain_version(pjoin('montblanc', 'version.py')), description='GPU-accelerated RIME implementations.', long_description=readme(), url='http://github.com/ska-sa/montblanc', From 7166c5a3838e38d45a34150ba4d030097296fd4b Mon Sep 17 00:00:00 2001 From: Simon Perkins Date: Thu, 14 Sep 2017 11:04:17 +0200 Subject: [PATCH 7/8] Depend on pybind11 and cppimport --- setup.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/setup.py b/setup.py index 2541c5e29..e1598f7bc 100644 --- a/setup.py +++ b/setup.py @@ -181,10 +181,11 @@ def include_pkg_dirs(): else: # Add binary/C extension type packages install_requires += [ - 'astropy >= 1.3.0', + 'astropy >= 2.0.2', + 'cppimport >= 17.7.24', 'cerberus >= 1.1', 'numpy >= 1.11.3', - 'numexpr >= 2.6.1', + 'pybind11 >= 2.2.0', 'python-casacore >= 2.1.2', 'ruamel.yaml >= 0.15.22', "{} >= 1.3.0".format(tensorflow_package), From b5a04211b05b00eb05d055e84b6847a6ab8d7cad Mon Sep 17 00:00:00 2001 From: Simon Perkins Date: Thu, 14 Sep 2017 11:05:58 +0200 Subject: [PATCH 8/8] C++ pybind11 implementation of parallactic angles --- montblanc/util/casa_parangles.cpp | 93 ++++++++++++++++++++++++++++ montblanc/util/parallactic_angles.py | 25 ++++++++ setup.py | 3 +- 3 files changed, 119 insertions(+), 2 deletions(-) create mode 100644 montblanc/util/casa_parangles.cpp diff --git a/montblanc/util/casa_parangles.cpp b/montblanc/util/casa_parangles.cpp new file mode 100644 index 000000000..48d6b50cf --- /dev/null +++ b/montblanc/util/casa_parangles.cpp @@ -0,0 +1,93 @@ +/* +<% +setup_pybind11(cfg) +cfg['compiler_args'] = ['-std=c++11', '-fvisibility=hidden'] +cfg['libraries'] = ['casa_casa', 'casa_measures'] +%> +*/ + +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace py = pybind11; + +constexpr unsigned int flags = py::array::c_style | py::array::forcecast; + +template +py::array_t parallactic_angles( + py::array_t times, + py::array_t antenna_positions, + py::array_t phase_centre) +{ + py::gil_scoped_release release; + + int na = antenna_positions.shape(0); + int ntimes = times.shape(0); + + // Result array + py::array_t angles({ntimes, na}); + + std::vector itrf_antenna_positions; + itrf_antenna_positions.reserve(na); + + // Compute antenna positions in ITRF + for(int ant=0; ant(*times.data(time), "s"), + casa::MEpoch::UTC)); + + // For each antenna + for(int ant=0; ant convert(mzenith, casa::MDirection::J2000); + *angles.mutable_data(time, ant) = phase_dir.positionAngle(convert().getValue()); + } + } + + return angles; +} + +PYBIND11_MODULE(casa_parangles, m) { + m.doc() = "auto-compiled c++ extension"; + + m.def("parallactic_angles", ¶llactic_angles, py::return_value_policy::move); + m.def("parallactic_angles", ¶llactic_angles, py::return_value_policy::move); +} \ No newline at end of file diff --git a/montblanc/util/parallactic_angles.py b/montblanc/util/parallactic_angles.py index 51110704d..40128ae20 100644 --- a/montblanc/util/parallactic_angles.py +++ b/montblanc/util/parallactic_angles.py @@ -18,10 +18,16 @@ # You should have received a copy of the GNU General Public License # along with this program; if not, see . +import cppimport import numpy as np +import pybind11 import montblanc +cppimport.set_quiet(False) +mod = cppimport.imp("montblanc.util.casa_parangles") +mod = mod.util.casa_parangles + try: import pyrap.measures @@ -154,6 +160,25 @@ def _parallactic_angle_astropy(times, ap, fc): import time + t = time.clock() + cangles = mod.parallactic_angles(ftimes, antenna_positions, phase_centre) + print 'cangles dones in %f' % (time.clock() - t) + + t = time.time() + pa_astro = _parallactic_angle_astropy(ftimes, antenna_positions, phase_centre) + print "pa_astropy done in ", time.time() - t + + t = time.clock() + pangles = parallactic_angles(ftimes, antenna_positions, phase_centre) + print 'pangles dones in %f' % (time.clock() - t) + + assert np.allclose(cangles, pangles) + + import sys + sys.exit(0) + + import time + t = time.time() pa_astro = _parallactic_angle_astropy(ftimes, antenna_positions, phase_centre) print "pa_astropy done in ", time.time() - t diff --git a/setup.py b/setup.py index e1598f7bc..b9aa96a34 100644 --- a/setup.py +++ b/setup.py @@ -232,6 +232,5 @@ def include_pkg_dirs(): license='GPL2', install_requires=install_requires, packages=find_packages(), - package_data={'montblanc': include_pkg_dirs()}, - include_package_data=True, + package_data={'montblanc': include_pkg_dirs() + [pjoin("util", "*.cpp")] }, zip_safe=False)