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 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..c2c6832 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 .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`. 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)