Skip to content

(wip) Release utils [do not merge] #161

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
mtazzari wants to merge 4 commits into
base: master
Choose a base branch
from
Open

(wip) Release utils [do not merge] #161

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
1ecb743
[utils] Make utils available as galario.utils
mtazzari Nov 14, 2018
9fed4ac
[docs] Add the utils functions relevant to the public to the docs
mtazzari Nov 14, 2018
bb7b8de
[utils,tests,speed_benchmark] Rename udat->u, vdat->v, fint->vis
mtazzari Nov 14, 2018
783f875
[pyx] Add alternative names for useful constants.
mtazzari Nov 14, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 11 additions & 0 deletions docs/py-api.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
4 changes: 4 additions & 0 deletions python/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand Down
3 changes: 2 additions & 1 deletion python/__init__.py.in
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
from .double import arcsec, au, cgs_to_Jy, pc, deg, arcsec_to_rad, deg_to_rad, au_to_cm, pc_to_cm
16 changes: 10 additions & 6 deletions python/libcommon.pyx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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',
Expand All @@ -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`.
Expand Down
14 changes: 7 additions & 7 deletions python/speed_benchmark.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand All @@ -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):
Expand Down
80 changes: 40 additions & 40 deletions python/test_galario.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -202,16 +202,16 @@ 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)

ReInt = AmpInt * np.cos(PhaseInt)
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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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)
Expand Down Expand Up @@ -558,16 +558,16 @@ 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)
PhaseInt = np.angle(ReInt + 1j*ImInt)

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
Expand All @@ -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)

Expand All @@ -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)
Expand All @@ -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
Expand Down
Loading