From 1ecb743bba4d857d879e0e54e95fde3141f687da Mon Sep 17 00:00:00 2001 From: Marco Tazzari Date: Wed, 14 Nov 2018 17:04:04 +0000 Subject: [PATCH 1/4] [utils] Make utils available as galario.utils --- python/CMakeLists.txt | 4 ++++ python/__init__.py.in | 1 + 2 files changed, 5 insertions(+) diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index 00f22cf..ca472e5 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -35,6 +35,10 @@ configure_file( "${CMAKE_CURRENT_SOURCE_DIR}/__init__.py.in" "${PYGALARIO_DIR}/__init__.py" ) +configure_file( + "${CMAKE_CURRENT_SOURCE_DIR}/utils.py" + "${PYGALARIO_DIR}/utils.py" + ) include(wrap_lib.cmake) wrap_lib() diff --git a/python/__init__.py.in b/python/__init__.py.in index 2852969..ede03c3 100644 --- a/python/__init__.py.in +++ b/python/__init__.py.in @@ -27,5 +27,6 @@ if HAVE_CUDA: from . import single from . import double +from . import utils from .double import arcsec, au, cgs_to_Jy, pc, deg \ No newline at end of file From 9fed4ac55c12f625dde25f379f17c7d923d3d4c1 Mon Sep 17 00:00:00 2001 From: Marco Tazzari Date: Wed, 14 Nov 2018 17:05:04 +0000 Subject: [PATCH 2/4] [docs] Add the utils functions relevant to the public to the docs --- docs/py-api.rst | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/docs/py-api.rst b/docs/py-api.rst index 54974fa..dfd79c5 100644 --- a/docs/py-api.rst +++ b/docs/py-api.rst @@ -67,6 +67,17 @@ Other useful functions .. autofunction:: galario.double.reduce_chi2 + +Utilities +--------- +The `galario.utils` module contains a number of general purpose functions as well as Python version of the compiled functions present in `galario.single` and `galario.double`. + +.. autofunction:: galario.utils.sweep_ref +.. autofunction:: galario.utils.apply_rotation +.. autofunction:: galario.utils.apply_phase_array +.. autofunction:: galario.utils.unique_part +.. autofunction:: galario.utils.assert_allclose + .. _galario_exceptions: Exceptions From bb7b8de1b7a799a1803ab01f692447dc9ca38398 Mon Sep 17 00:00:00 2001 From: Marco Tazzari Date: Wed, 14 Nov 2018 17:07:01 +0000 Subject: [PATCH 3/4] [utils,tests,speed_benchmark] Rename udat->u, vdat->v, fint->vis For consistency with libcommon.pyx, where only u, v, vis are used --- python/speed_benchmark.py | 14 +++---- python/test_galario.py | 80 +++++++++++++++++++-------------------- python/utils.py | 70 ++++++++++++++++++++-------------- 3 files changed, 88 insertions(+), 76 deletions(-) diff --git a/python/speed_benchmark.py b/python/speed_benchmark.py index 0ae486c..bab1e63 100644 --- a/python/speed_benchmark.py +++ b/python/speed_benchmark.py @@ -87,16 +87,16 @@ def setup_chi2Image(nxy, nsamples): PA = 80. maxuv_generator = 3e3 - udat, vdat = create_sampling_points(nsamples, maxuv_generator, dtype='float64') + u, v = create_sampling_points(nsamples, maxuv_generator, dtype='float64') x, _, w = generate_random_vis(nsamples, options.dtype) - _, _, maxuv = matrix_size(udat, vdat) + _, _, maxuv = matrix_size(u, v) dxy = 1 / maxuv # create model image (it happens to have 0 imaginary part) image_ref = create_reference_image(size=nxy, kernel='gaussian', dtype=options.dtype) - return image_ref, dxy, udat, vdat, x.real.copy(), x.imag.copy(), w, dRA, dDec, PA + return image_ref, dxy, u, v, x.real.copy(), x.imag.copy(), w, dRA, dDec, PA def setup_chi2Profile(nxy, nsamples): @@ -113,19 +113,19 @@ def setup_chi2Profile(nxy, nsamples): # generate the samples maxuv_generator = 3e3 - udat, vdat = create_sampling_points(nsamples, maxuv_generator, dtype=options.dtype) + u, v = create_sampling_points(nsamples, maxuv_generator, dtype=options.dtype) x, _, w = generate_random_vis(nsamples, options.dtype) - _, _, maxuv = matrix_size(udat, vdat) + _, _, maxuv = matrix_size(u, v) maxuv /= wle_m dxy = 1 / maxuv # compute the matrix size and maxuv - # nxy, dxy = g_double.get_image_size(udat/wle_m, vdat/wle_m) + # nxy, dxy = g_double.get_image_size(u/wle_m, v/wle_m) # compute radial profile intensity = radial_profile(Rmin, dR, nrad, profile_mode, dtype=options.dtype, gauss_width=150.*arcsec) - return intensity, Rmin, dR, nxy, dxy, udat/wle_m, vdat/wle_m, x.real.copy(), x.imag.copy(), w, dRA, dDec, inc, PA + return intensity, Rmin, dR, nxy, dxy, u/wle_m, v/wle_m, x.real.copy(), x.imag.copy(), w, dRA, dDec, inc, PA def do_timing(options, input_data, gpu=False, tpb=0, omp_num_threads=0): diff --git a/python/test_galario.py b/python/test_galario.py index 88d23e9..afc3f2b 100644 --- a/python/test_galario.py +++ b/python/test_galario.py @@ -128,10 +128,10 @@ def test_R2C_vs_C2C(nsamples, real_type, rtol, atol, acc_lib, pars): # generate the samples maxuv_generator = 3.e3 - udat, vdat = create_sampling_points(nsamples, maxuv_generator, dtype=real_type) + u, v = create_sampling_points(nsamples, maxuv_generator, dtype=real_type) # compute the matrix nxy and maxuv - _, minuv, maxuv = matrix_size(udat, vdat) + _, minuv, maxuv = matrix_size(u, v) du = maxuv/nxy # create model image (it happens to have 0 imaginary part) @@ -142,8 +142,8 @@ def test_R2C_vs_C2C(nsamples, real_type, rtol, atol, acc_lib, pars): PA *= deg dRA *= arcsec dDec *= arcsec - dRArot, dDecrot, urot, vrot = apply_rotation(PA, dRA, dDec, udat, vdat) - dRArot_g, dDecrot_g, urot_g, vrot_g = acc_lib.uv_rotate(PA, dRA, dDec, udat, vdat) + dRArot, dDecrot, urot, vrot = apply_rotation(PA, dRA, dDec, u, v) + dRArot_g, dDecrot_g, urot_g, vrot_g = acc_lib.uv_rotate(PA, dRA, dDec, u, v) np.testing.assert_allclose(dRArot, dRArot_g) np.testing.assert_allclose(dDecrot, dDecrot_g) @@ -162,7 +162,7 @@ def test_R2C_vs_C2C(nsamples, real_type, rtol, atol, acc_lib, pars): # CPU/GPU version (galario) dxy = 1./nxy/du - vis_galario = acc_lib.sampleImage(ref_real, dxy, udat, vdat, dRA=dRA, dDec=dDec, PA=PA) + vis_galario = acc_lib.sampleImage(ref_real, dxy, u, v, dRA=dRA, dDec=dDec, PA=PA) # check python c2c vs galario assert_allclose(vis_galario.real, vis_c2c_shifted.real, rtol, atol) @@ -183,16 +183,16 @@ def test_interpolate(size, real_type, complex_type, rtol, atol, acc_lib): maxuv = 1000. reference_image = create_reference_image(size=size, dtype=real_type) - udat, vdat = create_sampling_points(nsamples, maxuv/2.2) + u, v = create_sampling_points(nsamples, maxuv/2.2) # this factor has to be > than 2 because the matrix cover between -maxuv/2 to +maxuv/2, # therefore the sampling points have to be contained inside. - udat = udat.astype(real_type) - vdat = vdat.astype(real_type) + u = u.astype(real_type) + v = v.astype(real_type) # no rotation du = maxuv/size - uroti, vroti = uv_idx_r2c(udat, vdat, du, size/2.) + uroti, vroti = uv_idx_r2c(u, v, du, size/2.) uroti = uroti.astype(real_type) vroti = vroti.astype(real_type) @@ -202,7 +202,7 @@ def test_interpolate(size, real_type, complex_type, rtol, atol, acc_lib): ReInt = int_bilin_MT(ft.real, uroti, vroti) ImInt = int_bilin_MT(ft.imag, uroti, vroti) AmpInt = int_bilin_MT(np.abs(ft), uroti, vroti) - uneg = udat < 0. + uneg = u < 0. ImInt[uneg] *= -1. PhaseInt = np.angle(ReInt + 1j*ImInt) @@ -210,8 +210,8 @@ def test_interpolate(size, real_type, complex_type, rtol, atol, acc_lib): ImInt = AmpInt * np.sin(PhaseInt) complexInt = acc_lib.interpolate(ft, du, - udat.astype(real_type), - vdat.astype(real_type)) + u.astype(real_type), + v.astype(real_type)) assert_allclose(ReInt, complexInt.real, rtol, atol) assert_allclose(ImInt, complexInt.imag, rtol, atol) @@ -311,16 +311,16 @@ def test_apply_phase_vis(real_type, complex_type, rtol, atol, acc_lib, pars): # generate the samples nsamples = 10000 maxuv_generator = 3.e3 - udat, vdat = create_sampling_points(nsamples, maxuv_generator, dtype=real_type) + u, v = create_sampling_points(nsamples, maxuv_generator, dtype=real_type) # generate mock visibility values fint = np.zeros(nsamples, dtype=complex_type) fint.real = np.random.random(nsamples) * 10. fint.imag = np.random.random(nsamples) * 30. - fint_numpy = apply_phase_array(udat, vdat, fint.copy(), dRA, dDec) + fint_numpy = apply_phase_array(u, v, fint.copy(), dRA, dDec) - fint_shifted = acc_lib.apply_phase_vis(dRA, dDec, udat, vdat, fint) + fint_shifted = acc_lib.apply_phase_vis(dRA, dDec, u, v, fint) assert_allclose(fint_numpy.real, fint_shifted.real, rtol, atol) assert_allclose(fint_numpy.imag, fint_shifted.imag, rtol, atol) @@ -371,7 +371,7 @@ def model4(R): # u, v points maxuv_generator = 3e3 - udat, vdat = create_sampling_points(nsamples, maxuv_generator, + u, v = create_sampling_points(nsamples, maxuv_generator, dtype=real_type) nxy, dxy = 4096, 6.42956326721e-08 @@ -407,16 +407,16 @@ def model4(R): image_asym2[np.where(image_asym2 < 1e-10)] = 0. # Compute visibilities of ORIGINAL image with CURRENT GALARIO algorithm (only: origin='upper') - vis_C_upper_image_upper = acc_lib.sampleImage(image_asym, dxy, udat, vdat, dRA=0.5, dDec=-3., PA=10.) + vis_C_upper_image_upper = acc_lib.sampleImage(image_asym, dxy, u, v, dRA=0.5, dDec=-3., PA=10.) # Compute visibilities of ORIGINAL image with NEW algorithm, origin='upper' - vis_py_upper_image_upper = py_sampleImage(image_asym, dxy, udat, vdat, dRA=0.5, dDec=-3., PA=10., origin='upper') + vis_py_upper_image_upper = py_sampleImage(image_asym, dxy, u, v, dRA=0.5, dDec=-3., PA=10., origin='upper') # Compute visibilities of LOWER ORIGIN image with NEW algorithm, origin='lower' - vis_py_lower_image_lower = py_sampleImage(image_asym2, dxy, udat, vdat, dRA=0.5, dDec=-3., PA=10., origin='lower') + vis_py_lower_image_lower = py_sampleImage(image_asym2, dxy, u, v, dRA=0.5, dDec=-3., PA=10., origin='lower') # Compute with C implementation - vis_C_lower_image_lower = acc_lib.sampleImage(image_asym2, dxy, udat, vdat, dRA=0.5, dDec=-3., PA=10., origin='lower') + vis_C_lower_image_lower = acc_lib.sampleImage(image_asym2, dxy, u, v, dRA=0.5, dDec=-3., PA=10., origin='lower') # check that they produce all the same visibilities assert_allclose(vis_py_upper_image_upper, vis_C_lower_image_lower, atol=0., rtol=rtol) @@ -450,9 +450,9 @@ def test_all(nsamples, real_type, rtol, atol, acc_lib, pars): # generate the samples maxuv_generator = 3.e3 - udat, vdat = create_sampling_points(nsamples, maxuv_generator, dtype=real_type) + u, v = create_sampling_points(nsamples, maxuv_generator, dtype=real_type) - _, minuv, maxuv = matrix_size(udat, vdat) + _, minuv, maxuv = matrix_size(u, v) dxy = 1. / maxuv # pixel size (rad) # create intensity profile and model image @@ -466,15 +466,15 @@ def test_all(nsamples, real_type, rtol, atol, acc_lib, pars): reference_image = sweep_ref(intensity, Rmin, dR, nxy, nxy, dxy, inc, dtype_image=real_type) # test sampleImage - vis_py_sampleImage = py_sampleImage(reference_image, dxy, udat, vdat, PA=PA, dRA=dRA, dDec=dDec) - vis_g_sampleImage = acc_lib.sampleImage(reference_image, dxy, udat, vdat, PA=PA, dRA=dRA, dDec=dDec) + vis_py_sampleImage = py_sampleImage(reference_image, dxy, u, v, PA=PA, dRA=dRA, dDec=dDec) + vis_g_sampleImage = acc_lib.sampleImage(reference_image, dxy, u, v, PA=PA, dRA=dRA, dDec=dDec) assert_allclose(vis_py_sampleImage.real, vis_g_sampleImage.real, rtol=rtol, atol=atol) assert_allclose(vis_py_sampleImage.imag, vis_g_sampleImage.imag, rtol=rtol, atol=np.abs(np.mean(vis_g_sampleImage.real))*rtol) # test sampleProfile - vis_py_sampleProfile = py_sampleProfile(intensity.copy(), Rmin, dR, nxy, dxy, udat, vdat, inc=inc, dRA=dRA, dDec=dDec, PA=PA) - vis_g_sampleProfile = acc_lib.sampleProfile(intensity, Rmin, dR, nxy, dxy, udat, vdat, inc=inc, dRA=dRA, dDec=dDec, PA=PA) + vis_py_sampleProfile = py_sampleProfile(intensity.copy(), Rmin, dR, nxy, dxy, u, v, inc=inc, dRA=dRA, dDec=dDec, PA=PA) + vis_g_sampleProfile = acc_lib.sampleProfile(intensity, Rmin, dR, nxy, dxy, u, v, inc=inc, dRA=dRA, dDec=dDec, PA=PA) # check galario vs python implementation assert_allclose(vis_g_sampleProfile.real, vis_py_sampleProfile.real, rtol=rtol, atol=atol) @@ -487,12 +487,12 @@ def test_all(nsamples, real_type, rtol, atol, acc_lib, pars): # test chi2Image x, _, w = generate_random_vis(nsamples, real_type) - chi2_pychi2Image = py_chi2Image(reference_image, dxy, udat, vdat, x.real.copy(), x.imag.copy(), w, dRA=dRA, dDec=dDec) - chi2_g_chi2Image = acc_lib.chi2Image(reference_image, dxy, udat, vdat, x.real.copy(), x.imag.copy(), w, dRA=dRA, dDec=dDec) + chi2_pychi2Image = py_chi2Image(reference_image, dxy, u, v, x.real.copy(), x.imag.copy(), w, dRA=dRA, dDec=dDec) + chi2_g_chi2Image = acc_lib.chi2Image(reference_image, dxy, u, v, x.real.copy(), x.imag.copy(), w, dRA=dRA, dDec=dDec) # test chi2Profile - chi2_pychi2Profile = py_chi2Profile(intensity, Rmin, dR, nxy, dxy, udat, vdat, x.real.copy(), x.imag.copy(), w, inc=inc, dRA=dRA, dDec=dDec) - chi2_g_chi2Profile = acc_lib.chi2Profile(intensity, Rmin, dR, nxy, dxy, udat, vdat, x.real.copy(), x.imag.copy(), w, inc=inc, dRA=dRA, dDec=dDec) + chi2_pychi2Profile = py_chi2Profile(intensity, Rmin, dR, nxy, dxy, u, v, x.real.copy(), x.imag.copy(), w, inc=inc, dRA=dRA, dDec=dDec) + chi2_g_chi2Profile = acc_lib.chi2Profile(intensity, Rmin, dR, nxy, dxy, u, v, x.real.copy(), x.imag.copy(), w, inc=inc, dRA=dRA, dDec=dDec) # check galario vs python implementation assert_allclose(chi2_pychi2Profile, chi2_g_chi2Profile, rtol=rtol, atol=atol) @@ -516,10 +516,10 @@ def test_loss(nsamples, real_type, complex_type, rtol, atol, acc_lib, pars): # generate the samples maxuv_generator = 3.e3 - udat, vdat = create_sampling_points(nsamples, maxuv_generator, dtype=real_type) + u, v = create_sampling_points(nsamples, maxuv_generator, dtype=real_type) # compute the matrix size and maxuv - size, minuv, maxuv = matrix_size(udat, vdat) + size, minuv, maxuv = matrix_size(u, v) # create model complex image (it happens to have 0 imaginary part) reference_image = create_reference_image(size=size, dtype=real_type) @@ -558,7 +558,7 @@ def test_loss(nsamples, real_type, complex_type, rtol, atol, acc_lib, pars): # phase ### du = maxuv/size - uroti, vroti = uv_idx(udat, vdat, du, size/2.) + uroti, vroti = uv_idx(u, v, du, size/2.) ReInt = int_bilin_MT(py_shift_cmplx.real, uroti, vroti).astype(real_type) ImInt = int_bilin_MT(py_shift_cmplx.imag, uroti, vroti).astype(real_type) AmpInt = int_bilin_MT(np.abs(py_shift_cmplx), uroti, vroti).astype(real_type) @@ -566,8 +566,8 @@ def test_loss(nsamples, real_type, complex_type, rtol, atol, acc_lib, pars): fint = AmpInt * (np.cos(PhaseInt) + 1j*np.sin(PhaseInt)) fint_acc = fint.copy() - fint_shifted = apply_phase_array(udat, vdat, fint, dRA, dDec) - fint_acc_shifted = acc_lib.apply_phase_vis(dRA, dDec, udat, vdat, fint_acc) + fint_shifted = apply_phase_array(u, v, fint, dRA, dDec) + fint_acc_shifted = acc_lib.apply_phase_vis(dRA, dDec, u, v, fint_acc) # lose some absolute precision here --> not anymore. Really? check by decreasing rtol, atol @@ -580,12 +580,12 @@ def test_loss(nsamples, real_type, complex_type, rtol, atol, acc_lib, pars): ### # interpolation ### - uroti, vroti = uv_idx_r2c(udat, vdat, du, size/2.) + uroti, vroti = uv_idx_r2c(u, v, du, size/2.) ReInt = int_bilin_MT(py_shift_cmplx.real, uroti, vroti).astype(real_type) ImInt = int_bilin_MT(py_shift_cmplx.imag, uroti, vroti).astype(real_type) AmpInt = int_bilin_MT(np.abs(py_shift_cmplx), uroti, vroti).astype(real_type) - uneg = udat < 0. + uneg = u < 0. ImInt[uneg] *= -1. PhaseInt = np.angle(ReInt + 1j*ImInt) @@ -594,8 +594,8 @@ def test_loss(nsamples, real_type, complex_type, rtol, atol, acc_lib, pars): complexInt = acc_lib.interpolate(py_shift_cmplx.astype(complex_type, order='C'), du, - udat.astype(real_type), - vdat.astype(real_type)) + u.astype(real_type), + v.astype(real_type)) assert_allclose(ReInt, complexInt.real, rtol, atol) assert_allclose(ImInt, complexInt.imag, rtol, atol) @@ -604,7 +604,7 @@ def test_loss(nsamples, real_type, complex_type, rtol, atol, acc_lib, pars): # now all steps in one function # -> MT removed this because there is already a test for sample and here it is not clear what is the reference. ### - # sampled = acc_lib.sampleImage(ref_real, dRA, dDec, du, udat, vdat) + # sampled = acc_lib.sampleImage(ref_real, dRA, dDec, du, u, v) # # # a lot of precision lost. Why? --> not anymore # # rtol = 1 diff --git a/python/utils.py b/python/utils.py index 229085d..f0ff28e 100644 --- a/python/utils.py +++ b/python/utils.py @@ -35,7 +35,7 @@ "unique_part", "assert_allclose", "apply_rotation"] -def py_sampleImage(reference_image, dxy, udat, vdat, dRA=0., dDec=0., PA=0., origin='upper'): +def py_sampleImage(reference_image, dxy, u, v, dRA=0., dDec=0., PA=0., origin='upper'): """ Python implementation of sampleImage. @@ -61,8 +61,8 @@ def py_sampleImage(reference_image, dxy, udat, vdat, dRA=0., dDec=0., PA=0., ori cos_PA = np.cos(PA) sin_PA = np.sin(PA) - urot = udat * cos_PA - vdat * sin_PA - vrot = udat * sin_PA + vdat * cos_PA + urot = u * cos_PA - v * sin_PA + vrot = u * sin_PA + v * cos_PA dRArot = dRA * cos_PA - dDec * sin_PA dDecrot = dRA * sin_PA + dDec * cos_PA @@ -99,7 +99,7 @@ def py_sampleImage(reference_image, dxy, udat, vdat, dRA=0., dDec=0., PA=0., ori return vis -def py_sampleProfile(intensity, Rmin, dR, nxy, dxy, udat, vdat, dRA=0., dDec=0., PA=0, inc=0.): +def py_sampleProfile(intensity, Rmin, dR, nxy, dxy, u, v, dRA=0., dDec=0., PA=0, inc=0.): """ Python implementation of sampleProfile. @@ -129,29 +129,29 @@ def py_sampleProfile(intensity, Rmin, dR, nxy, dxy, udat, vdat, dRA=0., dDec=0., intensmap[nrow//2, ncol//2] = central_pixel(intensity, Rmin, dR, dxy) - vis = py_sampleImage(intensmap, dxy, udat, vdat, PA=PA, dRA=dRA, dDec=dDec) + vis = py_sampleImage(intensmap, dxy, u, v, PA=PA, dRA=dRA, dDec=dDec) return vis -def py_chi2Image(reference_image, dxy, udat, vdat, vis_obs_re, vis_obs_im, weights, dRA=0., dDec=0., PA=0.): +def py_chi2Image(reference_image, dxy, u, v, vis_obs_re, vis_obs_im, weights, dRA=0., dDec=0., PA=0.): """ Python implementation of chi2Image. """ - vis = py_sampleImage(reference_image, dxy, udat, vdat, PA=PA, dRA=dRA, dDec=dDec) + vis = py_sampleImage(reference_image, dxy, u, v, PA=PA, dRA=dRA, dDec=dDec) chi2 = np.sum(((vis.real - vis_obs_re)**2. + (vis.imag - vis_obs_im)**2.)*weights) return chi2 -def py_chi2Profile(intensity, Rmin, dR, nxy, dxy, udat, vdat, vis_obs_re, vis_obs_im, weights, dRA=0., dDec=0., PA=0, inc=0.): +def py_chi2Profile(intensity, Rmin, dR, nxy, dxy, u, v, vis_obs_re, vis_obs_im, weights, dRA=0., dDec=0., PA=0, inc=0.): """ Python implementation of chi2Profile. """ - vis = py_sampleProfile(intensity, Rmin, dR, nxy, dxy, udat, vdat, inc=inc, PA=PA, dRA=dRA, dDec=dDec) + vis = py_sampleProfile(intensity, Rmin, dR, nxy, dxy, u, v, inc=inc, PA=PA, dRA=dRA, dDec=dDec) chi2 = np.sum(((vis.real - vis_obs_re)**2. + (vis.imag - vis_obs_im)**2.)*weights) @@ -176,6 +176,8 @@ def central_pixel(I, Rmin, dR, dxy): """ Compute brightness in the central pixel as the average flux in the pixel. + TODO: add docs + """ # with quadrature method: tends to over-estimate it # area = np.pi*((dxy/2.)**2-Rmin**2) @@ -364,25 +366,25 @@ def create_sampling_points(nsamples, maxuv=1., dtype='float64'): return u.astype(dtype), v.astype(dtype) -def uv_idx(udat, vdat, du, half_size): +def uv_idx(u, v, du, half_size): """ For C2C transform. uv coordinates to pixel coordinates in range [0, npixels]. Assume image is square, same boundary in u and v direction. """ - return half_size + udat/du, half_size + vdat/du + return half_size + u/du, half_size + v/du -def uv_idx_r2c(udat, vdat, du, half_size): +def uv_idx_r2c(u, v, du, half_size): """ For R2C transform. uv coordinates to pixel coordinates in range [0, npixels]. Assume image is square, same boundary in u and v direction. """ - indu = np.abs(udat) / du - indv = half_size + vdat / du - uneg = udat < 0. - indv[uneg] = half_size - vdat[uneg] / du + indu = np.abs(u) / du + indv = half_size + v / du + uneg = u < 0. + indv[uneg] = half_size - v[uneg] / du return indu, indv @@ -407,12 +409,12 @@ def int_bilin_MT(f, x, y): return fint -def matrix_size(udat, vdat, **kwargs): +def matrix_size(u, v, **kwargs): maxuv_factor = kwargs.get('maxuv_factor', 4.8) minuv_factor = kwargs.get('minuv_factor', 4.) - uvdist = np.sqrt(udat**2 + vdat**2) + uvdist = np.sqrt(u**2 + v**2) maxuv = max(uvdist)*maxuv_factor minuv = min(uvdist)/minuv_factor @@ -424,7 +426,7 @@ def matrix_size(udat, vdat, **kwargs): return Nuv, minuv, maxuv -def apply_phase_array(u, v, fint, x0, y0): +def apply_phase_array(u, v, vis, x0, y0): """ Performs a translation in the real space by applying a phase shift in the Fourier space. This function applies the shift to data points sampling the Fourier transform of an image. @@ -433,7 +435,7 @@ def apply_phase_array(u, v, fint, x0, y0): ---------- u, v: 1D float array Coordinates of points in the Fourier space. units: observing wavelength - fint: 1D float array, complex + vis: 1D float array, complex Fourier Transform sampled in the (u, v) points. Re, Im, u, v must have the same length. x0, y0: floats, rad @@ -441,7 +443,7 @@ def apply_phase_array(u, v, fint, x0, y0): Returns ------- - fint_shifted: 1D float array, complex + vis_shifted: 1D float array, complex Phase-shifted of the Fourier Transform sampled in the (u, v) points. """ @@ -452,9 +454,9 @@ def apply_phase_array(u, v, fint, x0, y0): theta = u*x0 + v*y0 # apply the phase change - fint_shifted = fint * (np.cos(theta) + 1j*np.sin(theta)) + vis_shifted = vis * (np.cos(theta) + 1j * np.sin(theta)) - return fint_shifted + return vis_shifted def generate_random_vis(nsamples, dtype): @@ -469,15 +471,19 @@ def generate_random_vis(nsamples, dtype): return x, y, w -def apply_rotation(PA, dRA, dDec, udat, vdat): - """ Rotates the RA, Dec offsets and the udat and vdat coordinates by Position Angle PA """ +def apply_rotation(PA, dRA, dDec, u, v): + """ + Rotate the RA, Dec offsets and the u and v coordinates by Position Angle PA + + TODO: add docs + """ # PA: rad cos_PA = np.cos(PA) sin_PA = np.sin(PA) - urot = udat * cos_PA - vdat * sin_PA - vrot = udat * sin_PA + vdat * cos_PA + urot = u * cos_PA - v * sin_PA + vrot = u * sin_PA + v * cos_PA dRArot = dRA * cos_PA - dDec * sin_PA dDecrot = dRA * sin_PA + dDec * cos_PA @@ -486,12 +492,18 @@ def apply_rotation(PA, dRA, dDec, udat, vdat): def unique_part(array): - """Extract the unique part of a real-to-complex Fourier transform""" + """ + Extract the unique part of a real-to-complex Fourier transform + + """ return array[:, 0:int(array.shape[1]/2)+1] def assert_allclose(x, y, rtol=1e-10, atol=1e-8): - """Drop in replacement for `numpy.testing.assert_allclose` that shows the nonmatching elements""" + """ + Drop in replacement for `numpy.testing.assert_allclose` that shows the nonmatching elements + + """ if np.isscalar(x) and np.isscalar(y) == 1: return np.testing.assert_allclose(x, y, rtol=rtol, atol=atol) From 783f875f070bd689f6409073ea869276eb4a8431 Mon Sep 17 00:00:00 2001 From: Marco Tazzari Date: Wed, 14 Nov 2018 17:27:19 +0000 Subject: [PATCH 4/4] [pyx] Add alternative names for useful constants. Some people (including me) are more confortable with more explicit names that make it apparent the conversion, instead of relying on memory to remember that all units are in radians, or cm. --- python/__init__.py.in | 2 +- python/libcommon.pyx | 16 ++++++++++------ 2 files changed, 11 insertions(+), 7 deletions(-) diff --git a/python/__init__.py.in b/python/__init__.py.in index ede03c3..c2c6832 100644 --- a/python/__init__.py.in +++ b/python/__init__.py.in @@ -29,4 +29,4 @@ from . import single from . import double from . import utils -from .double import arcsec, au, cgs_to_Jy, pc, deg \ No newline at end of file +from .double import arcsec, au, cgs_to_Jy, pc, deg, arcsec_to_rad, deg_to_rad, au_to_cm, pc_to_cm \ No newline at end of file diff --git a/python/libcommon.pyx b/python/libcommon.pyx index 910e3ad..096eff9 100644 --- a/python/libcommon.pyx +++ b/python/libcommon.pyx @@ -30,6 +30,7 @@ include "galario_config.pxi" cimport galario_defs as cpp __all__ = ['arcsec', 'deg', 'cgs_to_Jy', 'pc', 'au', + 'arcsec_to_rad', 'deg_to_rad', 'cgs_to_Jy', 'pc_to_cm', 'au_to_cm', '_init', '_cleanup', 'set_v_origin', 'ngpus', 'use_gpu', 'threads', 'check_obs', 'check_image_size', 'get_image_size', @@ -40,12 +41,15 @@ __all__ = ['arcsec', 'deg', 'cgs_to_Jy', 'pc', 'au', # CONSTANTS -arcsec = 4.84813681109536e-06 # radians -deg = 0.017453292519943295 # radians -cgs_to_Jy = 1e23 # 1 Jy = 1.0e-23 erg/(s cm^2 Hz) -pc = 3.0856775815e18 # cm (IAU 2015 Resolution B2) -au = 1.49597870700e13 # cm (IAU 2012 Resolution B1) - +arcsec = 4.84813681109536e-06 # arcsecond to radians (1 arcsec = pi/3600/180 rad) +deg = 0.017453292519943295 # degree to radians (1 deg = pi/180 rad) +cgs_to_Jy = 1e23 # from cgs units to Jansky (1 Jy = 1.0e-23 erg/(s cm^2 Hz)) +pc = 3.0856775815e18 # parsec to cm (IAU 2015 Resolution B2) +au = 1.49597870700e13 # astronomical unit to cm (IAU 2012 Resolution B1) +arcsec_to_rad = arcsec # arcsecond to radians [alternative name for arcsec] +deg_to_rad = deg # degree to radians [alternative name for deg] +pc_to_cm = pc # parsec to cm [alternative name for pc] +au_to_cm = au # astronomical unit to cm [alternative name for au] cdef class ArrayWrapper: """Wrap an array allocated in C that has to be deleted by `free`.