From 7997d23da15246ac358c27c5183d3e31995f5fbe Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Tue, 2 Nov 2021 22:11:25 -0500 Subject: [PATCH 01/23] Added sampleUnstructuredImage function. Currently uses interpolation from scipy. --- python/libcommon.pyx | 95 +++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 94 insertions(+), 1 deletion(-) diff --git a/python/libcommon.pyx b/python/libcommon.pyx index 5420888..27bdab9 100644 --- a/python/libcommon.pyx +++ b/python/libcommon.pyx @@ -25,6 +25,8 @@ from cpython cimport PyObject, Py_INCREF import numpy as np np.import_array() +from scipy.interpolate import LinearNDInterpolator + include "galario_config.pxi" cimport galario_defs as cpp @@ -33,7 +35,7 @@ __all__ = ['arcsec', 'deg', 'cgs_to_Jy', 'pc', 'au', '_init', '_cleanup', 'set_v_origin', 'ngpus', 'use_gpu', 'threads', 'check_obs', 'check_image_size', 'get_image_size', - 'sampleImage', 'sampleProfile', 'chi2Image', 'chi2Profile', + 'sampleImage', 'sampleUnstructuredImage', 'sampleProfile', 'chi2Image', 'chi2Profile', 'get_coords_meshgrid', 'sweep', 'uv_rotate', 'interpolate', 'apply_phase_vis', 'reduce_chi2', '_fft2d', '_fftshift', '_fftshift_axis0'] @@ -443,6 +445,97 @@ def sampleImage(dreal[:,::1] image, dxy, dreal[::1] u, dreal[::1] v, return vis +def sampleUnstructuredImage(dreal[::1] x, dreal[::1] y, dreal[:,::1] image, + nxy, dxy, dreal[::1] u, dreal[::1] v, + dRA=0., dDec=0., PA=0., check=False, origin='upper'): + """ + Compute the synthetic visibilities of a model image at the specified (u, v) locations. + + The 2D surface brightness in `image` is Fourier transformed and sampled in the + (u, v) locations given in the `u` and `v` arrays. + + Typical call signature:: + + vis = sampleImage(image, dxy, u, v, dRA=0, dDec=0, PA=0, check=False, origin='upper') + + Parameters + ---------- + x : 1D array_like, float + List of x coordinates at which intensities are known. + **units**: rad + y : 1D array_like, float + List of y coordinates at which intensities are known. + **units**: rad + image : 1D array_like, float + Array containing the surface brightness of the model. + Assume the x-axis (R.A.) increases from right (West) to left (East) + and the y-axis (Dec.) increases from bottom (South) to top (North). + `nxy` must be even. + **units**: Jy/st + nxy : int + Number of pixels to use for the interpolated gridded image. + dxy : float + Size of the image cell in the interpolated image, assumed equal in both x and y direction. + **units**: rad + u : array_like, float + u coordinate of the visibility points where the FT has to be sampled. + **units**: wavelength + v : array_like, float + v coordinate of the visibility points where the FT has to be sampled. + The length of v must be equal to the length of u. + **units**: wavelength + dRA : float, optional + R.A. offset w.r.t. the phase center by which the image is translated. + If dRA > 0 translate the image towards the left (East). Default is 0. + **units**: rad + dDec : float, optional + Dec. offset w.r.t. the phase center by which the image is translated. + If dDec > 0 translate the image towards the top (North). Default is 0. + **units**: rad + PA : float, optional + Position Angle, defined East of North. Default is 0. + **units**: rad + check : bool, optional + If True, check whether `image` and `dxy` satisfy Nyquist criterion for + computing the synthetic visibilities in the (u, v) locations provided. + Additionally check that the (u, v) points fall in the image to avoid + segmentation violations. Default is False since the check might take + time. For executions where speed is important, set to False. + origin : ['upper' | 'lower'], optional + Set the [0,0] pixel index of the matrix in the upper left or lower left corner of the axes. + It follows the same convention as in matplotlib `matshow` and `imshow` commands. + Declination axis and the matrix y axis are parallel for `origin='lower'`, anti-parallel for `origin='upper'`. + The central pixel corresponding to the (RA, Dec) = (0, 0) is always [Nxy/2, Nxy/2]. + For more details see the Technical Requirements page in the online docs. + + Returns + ------- + vis : array_like, complex + Synthetic visibilities sampled in the (u, v) locations given in `u` and `v`. + **units**: Jy + + """ + #nxy = image.shape[0] + + interp = LinearNDInterpolator(list(zip(x, y)), image, fill_value=0) + + _, _, grid_x, grid_y, _ = get_coords_meshgrid(nxy, nxy, dxy, origin=origin) + cdef dreal[:,::1] new_image = interp(grid_x, grid_y) * dxy**2 + + duv = 1 / (dxy*nxy) + + if check: + check_image_size(u, v, nxy, dxy, duv) + + vis = np.zeros(len(u), dtype=complex_dtype) + v_origin = set_v_origin(origin) + cpp._sample_image(nxy, nxy, &new_image[0,0], v_origin, dRA, dDec, duv, PA, len(u), &u[0], &v[0], np.PyArray_DATA(vis)) + + return vis + + + + def sampleProfile(dreal[::1] intensity, Rmin, dR, nxy, dxy, dreal[::1] u, dreal[::1] v, dRA=0., dDec=0., PA=0., inc=0., check=False): """ From e21e2491b074d6a12ef92b4920b883da4459ad8f Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Wed, 3 Nov 2021 06:48:17 -0500 Subject: [PATCH 02/23] The image array in sampleUnstructuredImage is 1D. --- python/libcommon.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/libcommon.pyx b/python/libcommon.pyx index 27bdab9..85ba1f9 100644 --- a/python/libcommon.pyx +++ b/python/libcommon.pyx @@ -445,7 +445,7 @@ def sampleImage(dreal[:,::1] image, dxy, dreal[::1] u, dreal[::1] v, return vis -def sampleUnstructuredImage(dreal[::1] x, dreal[::1] y, dreal[:,::1] image, +def sampleUnstructuredImage(dreal[::1] x, dreal[::1] y, dreal[::1] image, nxy, dxy, dreal[::1] u, dreal[::1] v, dRA=0., dDec=0., PA=0., check=False, origin='upper'): """ From cbc2ef3ef53e97b1fce8e5b8265e0e5fab31b466 Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Wed, 3 Nov 2021 14:33:04 -0500 Subject: [PATCH 03/23] Added binning when pixels are super-sampled. --- python/libcommon.pyx | 45 +++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 42 insertions(+), 3 deletions(-) diff --git a/python/libcommon.pyx b/python/libcommon.pyx index 85ba1f9..37e5550 100644 --- a/python/libcommon.pyx +++ b/python/libcommon.pyx @@ -26,6 +26,7 @@ import numpy as np np.import_array() from scipy.interpolate import LinearNDInterpolator +from scipy.spatial import Voronoi, ConvexHull include "galario_config.pxi" @@ -446,7 +447,7 @@ def sampleImage(dreal[:,::1] image, dxy, dreal[::1] u, dreal[::1] v, def sampleUnstructuredImage(dreal[::1] x, dreal[::1] y, dreal[::1] image, - nxy, dxy, dreal[::1] u, dreal[::1] v, + int nxy, dreal dxy, dreal[::1] u, dreal[::1] v, dRA=0., dDec=0., PA=0., check=False, origin='upper'): """ Compute the synthetic visibilities of a model image at the specified (u, v) locations. @@ -515,13 +516,51 @@ def sampleUnstructuredImage(dreal[::1] x, dreal[::1] y, dreal[::1] image, **units**: Jy """ - #nxy = image.shape[0] + # Use scipy to inerpolate onto a regular grid. interp = LinearNDInterpolator(list(zip(x, y)), image, fill_value=0) - _, _, grid_x, grid_y, _ = get_coords_meshgrid(nxy, nxy, dxy, origin=origin) + cdef dreal[:,::1] grid_x, grid_y + grid_x_1D, grid_y_1D, grid_x, grid_y, _ = get_coords_meshgrid(nxy, nxy, \ + dxy, origin=origin) cdef dreal[:,::1] new_image = interp(grid_x, grid_y) * dxy**2 + # In pixels where we oversample, average instead of interpolate in case the + # intensity is varying quickly over the pixel. And use the volume of the + # associated Voronoi cell to weight each point being averaged. + cdef int[::1] i, j + cdef int[:,::1] npoints = np.zeros((nxy, nxy), dtype=np.dtype('i')) + cdef dreal[:,::1] binned_image = np.zeros((nxy, nxy)) + cdef dreal[:,::1] binned_weights = np.zeros((nxy, nxy)) + cdef int k, l, m + cdef int nx = x.shape[0] + + i = ((x - grid_x_1D.min()) / dxy + 0.5).astype(np.dtype('i')) + j = ((y - grid_y_1D.min()) / dxy + 0.5).astype(np.dtype('i')) + + vor = Voronoi(list(zip(x, y))) + cdef dreal[::1] vol = np.zeros(vor.npoints)+1 + for k, reg_num in enumerate(vor.point_region): + indices = vor.regions[reg_num] + if -1 in indices: + vol[k] = np.inf + else: + vol[k] = ConvexHull(vor.vertices[indices]).volume + + with nogil: + for k in range(nx): + if j[k] >= 0 and j[k] < nxy and i[k] >= 0 and i[k] < nxy: + npoints[j[k],i[k]] += 1 + binned_image[j[k],i[k]] += image[k]*vol[k] + binned_weights[j[k],i[k]] += vol[k] + + for l in range(nxy): + for m in range(nxy): + if npoints[l,m] > 1: + new_image[l,m] = binned_image[l,m] / binned_weights[l,m] * \ + dxy**2 + + # Now pick back up with what is typically done for regular grids. duv = 1 / (dxy*nxy) if check: From 043e8cbdea364b7b26838da1a918248c413ab043 Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Wed, 3 Nov 2021 20:50:13 -0500 Subject: [PATCH 04/23] Add the option to recycle weights to cut time on spectral line transforms. --- python/libcommon.pyx | 25 +++++++++++++++---------- 1 file changed, 15 insertions(+), 10 deletions(-) diff --git a/python/libcommon.pyx b/python/libcommon.pyx index 37e5550..04d6401 100644 --- a/python/libcommon.pyx +++ b/python/libcommon.pyx @@ -448,7 +448,8 @@ def sampleImage(dreal[:,::1] image, dxy, dreal[::1] u, dreal[::1] v, def sampleUnstructuredImage(dreal[::1] x, dreal[::1] y, dreal[::1] image, int nxy, dreal dxy, dreal[::1] u, dreal[::1] v, - dRA=0., dDec=0., PA=0., check=False, origin='upper'): + dRA=0., dDec=0., PA=0., check=False, origin='upper', \ + dreal[::1] vol=None, return_weights=False): """ Compute the synthetic visibilities of a model image at the specified (u, v) locations. @@ -538,14 +539,15 @@ def sampleUnstructuredImage(dreal[::1] x, dreal[::1] y, dreal[::1] image, i = ((x - grid_x_1D.min()) / dxy + 0.5).astype(np.dtype('i')) j = ((y - grid_y_1D.min()) / dxy + 0.5).astype(np.dtype('i')) - vor = Voronoi(list(zip(x, y))) - cdef dreal[::1] vol = np.zeros(vor.npoints)+1 - for k, reg_num in enumerate(vor.point_region): - indices = vor.regions[reg_num] - if -1 in indices: - vol[k] = np.inf - else: - vol[k] = ConvexHull(vor.vertices[indices]).volume + if vol is None: + vor = Voronoi(list(zip(x, y))) + vol = np.zeros(vor.npoints) + for k, reg_num in enumerate(vor.point_region): + indices = vor.regions[reg_num] + if -1 in indices: + vol[k] = np.inf + else: + vol[k] = ConvexHull(vor.vertices[indices]).volume with nogil: for k in range(nx): @@ -570,7 +572,10 @@ def sampleUnstructuredImage(dreal[::1] x, dreal[::1] y, dreal[::1] image, v_origin = set_v_origin(origin) cpp._sample_image(nxy, nxy, &new_image[0,0], v_origin, dRA, dDec, duv, PA, len(u), &u[0], &v[0], np.PyArray_DATA(vis)) - return vis + if return_weights: + return vis, vol + else: + return vis From d34629828e8c5f01f5cf7e064eb5e7b376a0679b Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Thu, 4 Nov 2021 16:38:21 -0500 Subject: [PATCH 05/23] Fixed some upper/lower difference issues. --- python/libcommon.pyx | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/python/libcommon.pyx b/python/libcommon.pyx index 04d6401..142732a 100644 --- a/python/libcommon.pyx +++ b/python/libcommon.pyx @@ -518,6 +518,10 @@ def sampleUnstructuredImage(dreal[::1] x, dreal[::1] y, dreal[::1] image, """ + if origin == "upper": + for n in range(y.size): + y[n] *= -1 + # Use scipy to inerpolate onto a regular grid. interp = LinearNDInterpolator(list(zip(x, y)), image, fill_value=0) @@ -536,8 +540,11 @@ def sampleUnstructuredImage(dreal[::1] x, dreal[::1] y, dreal[::1] image, cdef int k, l, m cdef int nx = x.shape[0] - i = ((x - grid_x_1D.min()) / dxy + 0.5).astype(np.dtype('i')) - j = ((y - grid_y_1D.min()) / dxy + 0.5).astype(np.dtype('i')) + i = ((x - grid_x_1D.max()) / -dxy + 0.5).astype(np.dtype('i')) + if origin == "upper": + j = ((y - grid_y_1D.max()) / -dxy + 0.5).astype(np.dtype('i')) + elif origin == "lower": + j = ((y - grid_y_1D.min()) / dxy + 0.5).astype(np.dtype('i')) if vol is None: vor = Voronoi(list(zip(x, y))) @@ -562,6 +569,10 @@ def sampleUnstructuredImage(dreal[::1] x, dreal[::1] y, dreal[::1] image, new_image[l,m] = binned_image[l,m] / binned_weights[l,m] * \ dxy**2 + if origin == "upper": + for n in range(y.size): + y[n] *= -1 + # Now pick back up with what is typically done for regular grids. duv = 1 / (dxy*nxy) From ce3598f2a4749e26d141540676c54efedac8af4d Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Thu, 11 Nov 2021 16:53:08 -0600 Subject: [PATCH 06/23] First go at C++ version of sample_unstructured_image. Slow because it uses brute force to find which triangle each point is in. --- .gitmodules | 3 + delaunator-cpp | 1 + python/galario_defs.pxd | 1 + python/libcommon.pyx | 91 ++++++++++++++++++++++++++- src/CMakeLists.txt | 2 +- src/galario.cpp | 135 ++++++++++++++++++++++++++++++++++++++++ src/galario.h | 2 + src/galario_py.h | 2 + 8 files changed, 235 insertions(+), 2 deletions(-) create mode 100644 .gitmodules create mode 160000 delaunator-cpp diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..b970088 --- /dev/null +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule "delaunator-cpp"] + path = delaunator-cpp + url = https://github.com/abellgithub/delaunator-cpp.git diff --git a/delaunator-cpp b/delaunator-cpp new file mode 160000 index 0000000..6f28799 --- /dev/null +++ b/delaunator-cpp @@ -0,0 +1 @@ +Subproject commit 6f2879967bc96a9bcdbacf418e560e9f2e170ace diff --git a/python/galario_defs.pxd b/python/galario_defs.pxd index fd06835..305d7e3 100644 --- a/python/galario_defs.pxd +++ b/python/galario_defs.pxd @@ -23,6 +23,7 @@ cdef extern from "galario_py.h" namespace "galario": # Main user functions void _sample_profile(int nr, void* intensity, dreal Rmin, dreal dR, dreal dxy, int nxy, dreal inc, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, void* u, void* v, void* vis) except + void _sample_image(int nx, int ny, void* image, dreal v_origin, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, void* u, void* v, void* vis) except + + void _sample_unstructured_image(void* x, void* y, int nx, int ny, dreal dxy, int ni, void* image, dreal v_origin, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, void* u, void* v, void* vis) except + dreal _chi2_profile(int nr, void* intensity, dreal Rmin, dreal dR, dreal dxy, int nxy, dreal inc, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, void* u, void* v, void* vis_obs_re, void* vis_obs_im, void* vis_obs_w) except + dreal _chi2_image(int nx, int ny, void* image, dreal v_origin, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, void* u, void* v, void* vis_obs_re, void* vis_obs_im, void* vis_obs_w) except + void _sweep(int nr, void* intensity, dreal Rmin, dreal dR, int nxy, dreal dxy, dreal inc, void* image) except + diff --git a/python/libcommon.pyx b/python/libcommon.pyx index 142732a..13a3633 100644 --- a/python/libcommon.pyx +++ b/python/libcommon.pyx @@ -36,7 +36,7 @@ __all__ = ['arcsec', 'deg', 'cgs_to_Jy', 'pc', 'au', '_init', '_cleanup', 'set_v_origin', 'ngpus', 'use_gpu', 'threads', 'check_obs', 'check_image_size', 'get_image_size', - 'sampleImage', 'sampleUnstructuredImage', 'sampleProfile', 'chi2Image', 'chi2Profile', + 'sampleImage', 'sampleUnstructuredImage', 'sampleUnstructuredImageCPP', 'sampleProfile', 'chi2Image', 'chi2Profile', 'get_coords_meshgrid', 'sweep', 'uv_rotate', 'interpolate', 'apply_phase_vis', 'reduce_chi2', '_fft2d', '_fftshift', '_fftshift_axis0'] @@ -583,10 +583,99 @@ def sampleUnstructuredImage(dreal[::1] x, dreal[::1] y, dreal[::1] image, v_origin = set_v_origin(origin) cpp._sample_image(nxy, nxy, &new_image[0,0], v_origin, dRA, dDec, duv, PA, len(u), &u[0], &v[0], np.PyArray_DATA(vis)) + """ if return_weights: return vis, vol else: return vis + """ + return new_image + + +def sampleUnstructuredImageCPP(dreal[::1] x, dreal[::1] y, dreal[::1] image, + int nxy, dreal dxy, dreal[::1] u, dreal[::1] v, + dRA=0., dDec=0., PA=0., check=False, origin='upper', \ + dreal[::1] vol=None, return_weights=False): + """ + Compute the synthetic visibilities of a model image at the specified (u, v) locations. + + The 2D surface brightness in `image` is Fourier transformed and sampled in the + (u, v) locations given in the `u` and `v` arrays. + + Typical call signature:: + + vis = sampleImage(image, dxy, u, v, dRA=0, dDec=0, PA=0, check=False, origin='upper') + + Parameters + ---------- + x : 1D array_like, float + List of x coordinates at which intensities are known. + **units**: rad + y : 1D array_like, float + List of y coordinates at which intensities are known. + **units**: rad + image : 1D array_like, float + Array containing the surface brightness of the model. + Assume the x-axis (R.A.) increases from right (West) to left (East) + and the y-axis (Dec.) increases from bottom (South) to top (North). + `nxy` must be even. + **units**: Jy/st + nxy : int + Number of pixels to use for the interpolated gridded image. + dxy : float + Size of the image cell in the interpolated image, assumed equal in both x and y direction. + **units**: rad + u : array_like, float + u coordinate of the visibility points where the FT has to be sampled. + **units**: wavelength + v : array_like, float + v coordinate of the visibility points where the FT has to be sampled. + The length of v must be equal to the length of u. + **units**: wavelength + dRA : float, optional + R.A. offset w.r.t. the phase center by which the image is translated. + If dRA > 0 translate the image towards the left (East). Default is 0. + **units**: rad + dDec : float, optional + Dec. offset w.r.t. the phase center by which the image is translated. + If dDec > 0 translate the image towards the top (North). Default is 0. + **units**: rad + PA : float, optional + Position Angle, defined East of North. Default is 0. + **units**: rad + check : bool, optional + If True, check whether `image` and `dxy` satisfy Nyquist criterion for + computing the synthetic visibilities in the (u, v) locations provided. + Additionally check that the (u, v) points fall in the image to avoid + segmentation violations. Default is False since the check might take + time. For executions where speed is important, set to False. + origin : ['upper' | 'lower'], optional + Set the [0,0] pixel index of the matrix in the upper left or lower left corner of the axes. + It follows the same convention as in matplotlib `matshow` and `imshow` commands. + Declination axis and the matrix y axis are parallel for `origin='lower'`, anti-parallel for `origin='upper'`. + The central pixel corresponding to the (RA, Dec) = (0, 0) is always [Nxy/2, Nxy/2]. + For more details see the Technical Requirements page in the online docs. + + Returns + ------- + vis : array_like, complex + Synthetic visibilities sampled in the (u, v) locations given in `u` and `v`. + **units**: Jy + + """ + + # Now pick back up with what is typically done for regular grids. + duv = 1 / (dxy*nxy) + + if check: + check_image_size(u, v, nxy, dxy, duv) + + vis = np.zeros(len(u), dtype=complex_dtype) + v_origin = set_v_origin(origin) + cpp._sample_unstructured_image(&x[0], &y[0], nxy, nxy, dxy, len(x), &image[0], v_origin, dRA, dDec, duv, PA, len(u), &u[0], &v[0], np.PyArray_DATA(vis)) + + return vis + diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 5c16505..953a52b 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -71,7 +71,7 @@ OPTION(GALARIO_TIMING "Output timing of selected functions. For testing only. De foreach(t IN ITEMS galario_single galario) target_link_libraries(${t} ${FFTW3_LIBRARIES}) - target_include_directories(${t} PUBLIC ${FFTW3_INCLUDE_DIRS} ${FFTW3_INCLUDE_DIR_PARALLEL}) + target_include_directories(${t} PUBLIC ${FFTW3_INCLUDE_DIRS} ${FFTW3_INCLUDE_DIR_PARALLEL} ../delaunator-cpp/include) if(GALARIO_TIMING) target_compile_definitions(${t} PRIVATE GALARIO_TIMING=1) endif() diff --git a/src/galario.cpp b/src/galario.cpp index 3d31e88..0c282cb 100644 --- a/src/galario.cpp +++ b/src/galario.cpp @@ -19,6 +19,7 @@ #include "galario.h" #include "galario_py.h" +#include // full function makes code hard to read #define tpb galario::threads() @@ -1351,6 +1352,96 @@ void sample_h(int nx, int ny, dcomplex* data, const dreal v_origin, dreal dRA, d #endif +namespace galario { +/** + * Interpolate from an unstructured image onto a regular grid. + */ +dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x, const dreal* y, const dreal* realdata, dreal v_origin) { + // Set up the Delauney triangulation. + + std::vector coords; + + for (int i=0; i < ni; i++) { + coords.push_back(x[i]); + coords.push_back(y[i]); + } + + delaunator::Delaunator d(coords); + + // Create an image including the appropriate coordinates. + auto gx = static_cast(malloc(sizeof(dreal)*nx)); + auto gy = static_cast(malloc(sizeof(dreal)*nx)); + auto image = static_cast(malloc(sizeof(dreal)*nx*ny)); + + for (int i = 0; i < nx; i++) + gx[i] = (0.5 - i * 1./nx) * nx * dxy; + for (int i = 0; i < ny; i++) + gy[i] = (0.5 - i * 1./ny) * ny * dxy * v_origin; + + // Now loop through the pixels in the image pixels, find the triangle each point is in, and interpolate. + for (int i = 0; i < nx; i++) { + for (int j = 0; j < ny; j++) { + bool found_triangle = false; + // First, find the triangle that this point is in. + for (int k = 0; k < d.triangles.size(); k+=3) { + int ia = d.triangles[k]; + double ax = x[ia]; + double ay = y[ia]; + + int ib = d.triangles[k+1]; + double bx = x[ib]; + double by = y[ib]; + + int ic = d.triangles[k+2]; + double cx = x[ic]; + double cy = y[ic]; + + double vbx = bx - ax; + double vby = by - ay; + double vcx = cx - ax; + double vcy = cy - ay; + + double det_vv2 = gx[i]*vcy - gy[j]*vcx; + double det_v0v2 = ax*vcy - ay*vcx; + double det_v1v2 = vbx*vcy - vby*vcx; + double det_vv1 = gx[i]*vby - gy[j]*vbx; + double det_v0v1 = ax*vby - ay*vbx; + + double a = (det_vv2 - det_v0v2) / det_v1v2; + double b = -(det_vv1 - det_v0v1) / det_v1v2; + + // We've found the right triangle, now interpolate. + if ((a > 0) & (b > 0) & (a + b < 1)) { + double wa = ((by - cy)*(gx[i] - cx) + (cx - bx)*(gy[j] - cy)) / + ((by - cy)*(ax - cx) + (cx - bx)*(ay - cy)); + double wb = ((cy - ay)*(gx[i] - cx) + (ax - cx)*(gy[j] - cy)) / + ((by - cy)*(ax - cx) + (cx - bx)*(ay - cy)); + double wc = 1 - wa - wb; + + image[i * nx + j] = (wa*realdata[ia] + wb*realdata[ib] + wc*realdata[ic])*dxy*dxy; + + found_triangle = true; + break; + } + } + + // If no triangle was found, the point is outside the area with data so set to 0. + if (not found_triangle) + image[i * nx + j] = 0.; + } + } + + // Now copy to an image. + auto buffer = copy_input(nx, ny, image); + + return buffer; +} + +void* _interpolate_to_image(int nx, int ny, int ni, dreal dxy, void* x, void *y, void* realdata, dreal v_origin) { + return interpolate_to_image(nx, ny, ni, dxy, static_cast(x), static_cast(y), static_cast(realdata), v_origin); +} +} + namespace galario { /** @@ -1396,6 +1487,50 @@ void _sample_image(int nx, int ny, void* data, dreal v_origin, dreal dRA, dreal sample_image(nx, ny, static_cast(data), v_origin, dRA, dDec, duv, PA, nd, static_cast(u), static_cast(v), static_cast(vis_int)); } +/** + * return result in `vis_int` + */ +void sample_unstructured_image(const dreal* realx, const dreal* realy, int nx, int ny, dreal dxy, int ni, const dreal* realdata, dreal v_origin, dreal dRA, dreal dDec, dreal duv, + const dreal PA, int nd, const dreal* u, const dreal* v, dcomplex* vis_int) { + CPUTimer t_start; + + // Initialization for uv_idx and interpolate + CHECK_INPUT(nx); + +/*#ifdef __CUDACC__ + GPUTimer t_total; + CudaMemory vis_int_d(nd); + + auto data_d = copy_input_d(nx, ny, realdata); + + // do the actual computation + sample_d(nx, ny, data_d.ptr, v_origin, dRA, dDec, nd, duv, PA, u, v, vis_int_d.ptr); + + // retrieve interpolated values + CCheck(cudaDeviceSynchronize()); + + GPUTimer t; + vis_int_d.Retrieve(vis_int); + t.Elapsed("sample_image::vis_int_ D->H"); + + t_total.Elapsed("sample_image_tot"); +#else*/ + CPUTimer t; + + auto data = interpolate_to_image(nx, ny, ni, dxy, realx, realy, realdata, v_origin); t.Elapsed("sample_image::copy_input"); + + sample_h(nx, ny, data, v_origin, dRA, dDec, nd, duv, PA, u, v, vis_int); + + t = CPUTimer(); galario_free(data); t.Elapsed("sample_image::free_data"); +//#endif + t_start.Elapsed("sample_image_tot"); +} + +void _sample_unstructured_image(void* x, void* y, int nx, int ny, dreal dxy, int ni, void* data, dreal v_origin, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, void* u, void* v, void* vis_int) { + sample_unstructured_image(static_cast(x), static_cast(y), nx, ny, dxy, ni, static_cast(data), v_origin, dRA, dDec, duv, PA, nd, static_cast(u), static_cast(v), static_cast(vis_int)); +} + + /** * return result in `vis_int` diff --git a/src/galario.h b/src/galario.h index 9a3c4ab..f9ce326 100644 --- a/src/galario.h +++ b/src/galario.h @@ -27,6 +27,7 @@ namespace galario { void sample_profile(int nr, const dreal* intensity, dreal Rmin, dreal dR, dreal dxy, int nxy, dreal inc, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, const dreal *u, const dreal *v, dcomplex *vis_int); void sample_image(int nx, int ny, const dreal* image, const dreal v_origin, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, const dreal* u, const dreal* v, dcomplex* vis_int); +void sample_unstructured_image(const dreal* x, const dreal *y, int nx, int ny, dreal dxy, int ni, const dreal* image, const dreal v_origin, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, const dreal* u, const dreal* v, dcomplex* vis_int); dreal chi2_profile(int nr, const dreal* intensity, dreal Rmin, dreal dR, dreal dxy, int nxy, dreal inc, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, const dreal *u, const dreal *v, const dreal *vis_obs_re, const dreal *vis_obs_im, const dreal *weights); @@ -36,6 +37,7 @@ void uv_rotate(dreal PA, dreal dRA, dreal dDec, dreal* dRArot, dreal* dDecrot, i /* Interface for the experts */ dcomplex* copy_input(int nx, int ny, const dreal* image); +dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x, const dreal* y, const dreal* data, dreal v_origin); void galario_free(void*); void fft2d(int nx, int ny, dcomplex* image); void fftshift(int nx, int ny, dcomplex* image); diff --git a/src/galario_py.h b/src/galario_py.h index 1321fa8..cf96db4 100644 --- a/src/galario_py.h +++ b/src/galario_py.h @@ -30,6 +30,7 @@ namespace galario { void _sample_profile(int nr, void *intensity, dreal Rmin, dreal dR, dreal dxy, int nxy, dreal inc, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, void *u, void *v, void *vis_int); void _sample_image(int nx, int ny, void* data, dreal v_origin, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, void* u, void* v, void* vis_int); +void _sample_unstructured_image(void* x, void* y, int nx, int ny, dreal dxy, int ni, void* data, dreal v_origin, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, void* u, void* v, void* vis_int); dreal _chi2_profile(int nr, void *intensity, dreal Rmin, dreal dR, dreal dxy, int nxy, dreal inc, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, void *u, void *v, void *vis_obs_re, void *vis_obs_im, void *weights); dreal _chi2_image(int nx, int ny, void* data, dreal v_origin, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, void* u, void* v, void* vis_obs_re, void* vis_obs_im, void* weights); @@ -38,6 +39,7 @@ void _uv_rotate(dreal PA, dreal dRA, dreal dDec, void* dRArot, void* dDecrot, in /* Interface for the experts */ void* _copy_input(int nx, int ny, void* realdata); +void* interpolate_to_image(int nx, int ny, int ni, dreal dxy, void* x, void* y, void* data, dreal v_origin); void _fft2d(int nx, int ny, void* data); void _fftshift(int nx, int ny, void* data); void _fftshift_axis0(int nx, int ncol, void* data); From 84cd1854daa9d8a6020ae39d5f4caeb371404438 Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Fri, 12 Nov 2021 17:14:35 -0600 Subject: [PATCH 07/23] Added a directed walk to speed up triangle finding Also further fixes like passing a pointer to the Delauney triangulation and adding some timing tests. --- python/libcommon.pyx | 3 - src/galario.cpp | 228 ++++++++++++++++++++++++++++++++++++------- src/timer.h | 25 +++++ 3 files changed, 218 insertions(+), 38 deletions(-) create mode 100644 src/timer.h diff --git a/python/libcommon.pyx b/python/libcommon.pyx index 13a3633..c4449be 100644 --- a/python/libcommon.pyx +++ b/python/libcommon.pyx @@ -583,13 +583,10 @@ def sampleUnstructuredImage(dreal[::1] x, dreal[::1] y, dreal[::1] image, v_origin = set_v_origin(origin) cpp._sample_image(nxy, nxy, &new_image[0,0], v_origin, dRA, dDec, duv, PA, len(u), &u[0], &v[0], np.PyArray_DATA(vis)) - """ if return_weights: return vis, vol else: return vis - """ - return new_image def sampleUnstructuredImageCPP(dreal[::1] x, dreal[::1] y, dreal[::1] image, diff --git a/src/galario.cpp b/src/galario.cpp index 0c282cb..6b8c861 100644 --- a/src/galario.cpp +++ b/src/galario.cpp @@ -20,6 +20,7 @@ #include "galario.h" #include "galario_py.h" #include +#include "timer.h" // full function makes code hard to read #define tpb galario::threads() @@ -1352,6 +1353,130 @@ void sample_h(int nx, int ny, dcomplex* data, const dreal v_origin, dreal dRA, d #endif +/** + * Find the index of the triangle that a point is in using brute force. + */ +int find_triangle_bruteforce(delaunator::Delaunator *d, const dreal *x, const dreal *y, dreal gx, dreal gy) { + + bool found_triangle = false; + int which_triangle = -1; + + // Loop through all the triangles to brute force-find which one a point is in. + for (int k = 0; k < d->triangles.size(); k+=3) { + int ia = d->triangles[k]; + double ax = x[ia]; + double ay = y[ia]; + + int ib = d->triangles[k+1]; + double bx = x[ib]; + double by = y[ib]; + + int ic = d->triangles[k+2]; + double cx = x[ic]; + double cy = y[ic]; + + double vbx = bx - ax; + double vby = by - ay; + double vcx = cx - ax; + double vcy = cy - ay; + + double det_vv2 = gx*vcy - gy*vcx; + double det_v0v2 = ax*vcy - ay*vcx; + double det_v1v2 = vbx*vcy - vby*vcx; + double det_vv1 = gx*vby - gy*vbx; + double det_v0v1 = ax*vby - ay*vbx; + + double a = (det_vv2 - det_v0v2) / det_v1v2; + double b = -(det_vv1 - det_v0v1) / det_v1v2; + + // We've found the right triangle, now interpolate. + if ((a > 0) & (b > 0) & (a + b < 1)) { + which_triangle = k; + found_triangle = true; + break; + } + } + + return which_triangle; +} + +/** + * Find which triangle a point is in using a directed walk. + */ +int find_triangle_directedwalk(delaunator::Delaunator *d, const dreal *x, const dreal *y, dreal gx, dreal gy, int start, int* last_good, double *time) { + int which_triangle = -2; + int count = 0; + dreal eps = 1.0e-3; + bool found_triangle = false; + //TCREATE(moo); TCLEAR(moo); + //TSTART(moo); + while (count < d->triangles.size() / (3*4)) { + int ia = d->triangles[start]; + double ax = x[ia]; + double ay = y[ia]; + int ib = d->triangles[start+1]; + double bx = x[ib]; + double by = y[ib]; + int ic = d->triangles[start+2]; + double cx = x[ic]; + double cy = y[ic]; + + double wa = ((by - cy)*(gx - cx) + (cx - bx)*(gy - cy)) / + ((by - cy)*(ax - cx) + (cx - bx)*(ay - cy)); + double wb = ((cy - ay)*(gx - cx) + (ax - cx)*(gy - cy)) / + ((by - cy)*(ax - cx) + (cx - bx)*(ay - cy)); + double wc = 1 - wa - wb; + + if (wa < -eps) { + start = d->halfedges[start+1]; + } else if (wb < -eps) { + start = d->halfedges[start+2]; + } else if (wc < -eps) { + start = d->halfedges[start+0]; + } else { + which_triangle = start; + found_triangle = true; + } + + if (start >= 0) { + start = start - start%3; + *last_good = start; + } + else + which_triangle = start; + + if ((found_triangle) or (which_triangle == -1)) + break; + + count++; + } + // TSTOP(moo); + //if (count > 0) printf("gx = %f, gy = %f, count = %d \n", gx, gy, count); + //*time += TGIVE(moo); + + return which_triangle; +} + +/** + * First try to find the triangle index using a directed walk, and if that fails switch to brute force. + */ +int find_triangle(delaunator::Delaunator *d, const dreal *x, const dreal *y, dreal gx, dreal gy, int start, int* last_good, double* time) { +#ifdef GALARIO_TIMING + TCREATE(boo); TCLEAR(boo); TSTART(boo); +#endif + int which_triangle = find_triangle_directedwalk(d, x, y, gx, gy, start, last_good, time); + if (which_triangle == -2) { + printf("Switching to brute force \n"); + which_triangle = find_triangle_bruteforce(d, x, y, gx, gy); + } +#ifdef GALARIO_TIMING + TSTOP(boo); *time += TGIVE(boo); +#endif + + return which_triangle; +} + + namespace galario { /** * Interpolate from an unstructured image onto a regular grid. @@ -1361,12 +1486,26 @@ dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x std::vector coords; + dreal xmin = std::numeric_limits::max(); dreal xmax = -std::numeric_limits::max(); + dreal ymin = std::numeric_limits::max(); dreal ymax = -std::numeric_limits::max(); for (int i=0; i < ni; i++) { coords.push_back(x[i]); coords.push_back(y[i]); + + if (x[i] > xmax) xmax = x[i]; + if (x[i] < xmin) xmin = x[i]; + if (y[i] > ymax) ymax = y[i]; + if (y[i] < ymin) ymin = y[i]; } +#ifdef GALARIO_TIMING + TCREATE(moo); TCLEAR(moo); TSTART(moo); +#endif delaunator::Delaunator d(coords); +#ifdef GALARIO_TIMING + TSTOP(moo); + printf("Time to triangulate %f \n", TGIVE(moo)); +#endif // Create an image including the appropriate coordinates. auto gx = static_cast(malloc(sizeof(dreal)*nx)); @@ -1378,58 +1517,70 @@ dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x for (int i = 0; i < ny; i++) gy[i] = (0.5 - i * 1./ny) * ny * dxy * v_origin; + int which_triangle = 0; + int last_triangle = 0; + int col_start_triangle = -1; + double time = 0.; +#ifdef GALARIO_TIMING + TCLEAR(moo); +#endif // Now loop through the pixels in the image pixels, find the triangle each point is in, and interpolate. for (int i = 0; i < nx; i++) { + if ((i > 0) and (col_start_triangle > -1)) { + which_triangle = col_start_triangle; + last_triangle = col_start_triangle; + col_start_triangle = -1; + } for (int j = 0; j < ny; j++) { - bool found_triangle = false; - // First, find the triangle that this point is in. - for (int k = 0; k < d.triangles.size(); k+=3) { - int ia = d.triangles[k]; + // Check whether the triangle is out of the triangulation. + if ((gx[i] > xmin) and (gx[i] < xmax) and (gy[j] > ymin) and (gy[j] < ymax)) { + // Find which triangle this grid point is in. +#ifdef GALARIO_TIMING + TSTART(moo); +#endif + which_triangle = find_triangle(&d, x, y, gx[i], gy[j], which_triangle, &last_triangle, &time); +#ifdef GALARIO_TIMING + TSTOP(moo); +#endif + } + else + which_triangle = -1; + + // We've found the right triangle, now interpolate. + if (which_triangle > -1) { + int ia = d.triangles[which_triangle]; double ax = x[ia]; double ay = y[ia]; - - int ib = d.triangles[k+1]; + int ib = d.triangles[which_triangle+1]; double bx = x[ib]; double by = y[ib]; - - int ic = d.triangles[k+2]; + int ic = d.triangles[which_triangle+2]; double cx = x[ic]; double cy = y[ic]; - double vbx = bx - ax; - double vby = by - ay; - double vcx = cx - ax; - double vcy = cy - ay; - - double det_vv2 = gx[i]*vcy - gy[j]*vcx; - double det_v0v2 = ax*vcy - ay*vcx; - double det_v1v2 = vbx*vcy - vby*vcx; - double det_vv1 = gx[i]*vby - gy[j]*vbx; - double det_v0v1 = ax*vby - ay*vbx; - - double a = (det_vv2 - det_v0v2) / det_v1v2; - double b = -(det_vv1 - det_v0v1) / det_v1v2; - - // We've found the right triangle, now interpolate. - if ((a > 0) & (b > 0) & (a + b < 1)) { - double wa = ((by - cy)*(gx[i] - cx) + (cx - bx)*(gy[j] - cy)) / - ((by - cy)*(ax - cx) + (cx - bx)*(ay - cy)); - double wb = ((cy - ay)*(gx[i] - cx) + (ax - cx)*(gy[j] - cy)) / - ((by - cy)*(ax - cx) + (cx - bx)*(ay - cy)); - double wc = 1 - wa - wb; + double wa = ((by - cy)*(gx[i] - cx) + (cx - bx)*(gy[j] - cy)) / + ((by - cy)*(ax - cx) + (cx - bx)*(ay - cy)); + double wb = ((cy - ay)*(gx[i] - cx) + (ax - cx)*(gy[j] - cy)) / + ((by - cy)*(ax - cx) + (cx - bx)*(ay - cy)); + double wc = 1 - wa - wb; - image[i * nx + j] = (wa*realdata[ia] + wb*realdata[ib] + wc*realdata[ic])*dxy*dxy; + image[i * nx + j] = (wa*realdata[ia] + wb*realdata[ib] + wc*realdata[ic])*dxy*dxy; - found_triangle = true; - break; + if (col_start_triangle == -1) { + col_start_triangle = last_triangle; } } - // If no triangle was found, the point is outside the area with data so set to 0. - if (not found_triangle) + else { image[i * nx + j] = 0.; + which_triangle = last_triangle; + } } } +#ifdef GALARIO_TIMING + printf("Time to calculate barycentric coords %f \n", time); + printf("Time to find triangles %f \n", TGIVE(moo)); +#endif // Now copy to an image. auto buffer = copy_input(nx, ny, image); @@ -1517,7 +1668,14 @@ void sample_unstructured_image(const dreal* realx, const dreal* realy, int nx, i #else*/ CPUTimer t; - auto data = interpolate_to_image(nx, ny, ni, dxy, realx, realy, realdata, v_origin); t.Elapsed("sample_image::copy_input"); +#ifdef GALARIO_TIMING + TCREATE(moo); TCLEAR(moo); TSTART(moo); +#endif + auto data = interpolate_to_image(nx, ny, ni, dxy, realx, realy, realdata, v_origin); t.Elapsed("sample_image::interpolate_to_grid"); +#ifdef GALARIO_TIMING + TSTOP(moo); + printf("Time to interpolate: %f \n", TGIVE(moo)); +#endif sample_h(nx, ny, data, v_origin, dRA, dDec, nd, duv, PA, u, v, vis_int); diff --git a/src/timer.h b/src/timer.h new file mode 100644 index 0000000..033e8fe --- /dev/null +++ b/src/timer.h @@ -0,0 +1,25 @@ +#include +#include + +struct timeval tuse; + +#define CPU_TIME gettimeofday( &tuse, (struct timezone *)0 ); + +#define TCREATE(x) \ + double __timerseconds##x=0; double __timerstartseconds##x=0; \ + double __timerusec##x=0; double __timerstartusec##x=0; + +#define TCLEAR(x) {__timerseconds##x = 0; __timerusec##x = 0; } + +#define TSTART(x) { CPU_TIME; \ + __timerstartseconds##x = tuse.tv_sec; \ + __timerstartusec##x = tuse.tv_usec; } + +#define TSTOP(x) { CPU_TIME; \ + __timerseconds##x += (tuse.tv_sec - __timerstartseconds##x); \ + __timerusec##x += (tuse.tv_usec - __timerstartusec##x); } + +#define TTIME(str,x) \ + printf("%s %6.6f seconds \n", str, __timerseconds##x+__timerusec##x*1.0e-6); + +#define TGIVE(x) (__timerseconds##x+__timerusec##x*1.0e-6) From 82bc077e10e3c61c92451e3322fbf99524e9563f Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Fri, 12 Nov 2021 18:32:25 -0600 Subject: [PATCH 08/23] When there are multiple triangles in a binned cell, average. --- src/galario.cpp | 82 +++++++++++++++++++++++++++++++++++++++---------- 1 file changed, 65 insertions(+), 17 deletions(-) diff --git a/src/galario.cpp b/src/galario.cpp index 6b8c861..6800da4 100644 --- a/src/galario.cpp +++ b/src/galario.cpp @@ -1507,6 +1507,50 @@ dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x printf("Time to triangulate %f \n", TGIVE(moo)); #endif + // For each triangle, calculate the centroid and which grid cell it falls in. + auto tx = static_cast(malloc(sizeof(dreal)*d.triangles.size()/3)); + auto ty = static_cast(malloc(sizeof(dreal)*d.triangles.size()/3)); + auto tf = static_cast(malloc(sizeof(dreal)*d.triangles.size()/3)); + auto ta = static_cast(malloc(sizeof(dreal)*d.triangles.size()/3)); + + auto itx = static_cast(malloc(sizeof(int)*d.triangles.size()/3)); + auto ity = static_cast(malloc(sizeof(int)*d.triangles.size()/3)); + + dreal gx_max = 0.5*nx*dxy; + dreal gy_max = 0.5*ny*dxy*v_origin; + + for (int i = 0; i < d.triangles.size()/3; i++) { + tx[i] = (x[d.triangles[3*i]] + x[d.triangles[3*i+1]] + x[d.triangles[3*i+2]]) / 3.; + ty[i] = (y[d.triangles[3*i]] + y[d.triangles[3*i+1]] + y[d.triangles[3*i+2]]) / 3.; + tf[i] = (realdata[d.triangles[3*i]] + realdata[d.triangles[3*i+1]] + realdata[d.triangles[3*i+2]]) / 3.; + ta[i] = std::fabs((y[d.triangles[3*i+1]] - y[d.triangles[3*i]]) * (x[d.triangles[3*i+2]] - x[d.triangles[3*i+1]]) - + (x[d.triangles[3*i+1]] - x[d.triangles[3*i]]) * (y[d.triangles[3*i+2]] - y[d.triangles[3*i+1]])); + + itx[i] = trunc((tx[i] - gx_max) / (-dxy) + 0.5); + ity[i] = trunc((ty[i] - gy_max) / (-dxy*v_origin) + 0.5); + } + + // Create an image that bins the triangles weighted by their area. + auto binned_image = static_cast(malloc(sizeof(dreal)*nx*ny)); + auto binned_weights = static_cast(malloc(sizeof(dreal)*nx*ny)); + auto npoints = static_cast(malloc(sizeof(int)*nx*ny)); + + for (int i = 0; i < nx; i++) { + for (int j = 0; j < ny; j++) { + binned_image[i,j] = 0; + binned_weights[i,j] = 0; + npoints[i,j] = 0; + } + } + + for (int i = 0; i < d.triangles.size()/3; i++) { + if ((itx[i] >= 0) and (itx[i] < nx) and (ity[i] >= 0) and (ity[i] < ny)) { + npoints[itx[i] * nx + ity[i]] += 1; + binned_image[itx[i] * nx + ity[i]] += tf[i] * ta[i]; + binned_weights[itx[i] * nx + ity[i]] += ta[i]; + } + } + // Create an image including the appropriate coordinates. auto gx = static_cast(malloc(sizeof(dreal)*nx)); auto gy = static_cast(malloc(sizeof(dreal)*nx)); @@ -1548,23 +1592,27 @@ dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x // We've found the right triangle, now interpolate. if (which_triangle > -1) { - int ia = d.triangles[which_triangle]; - double ax = x[ia]; - double ay = y[ia]; - int ib = d.triangles[which_triangle+1]; - double bx = x[ib]; - double by = y[ib]; - int ic = d.triangles[which_triangle+2]; - double cx = x[ic]; - double cy = y[ic]; - - double wa = ((by - cy)*(gx[i] - cx) + (cx - bx)*(gy[j] - cy)) / - ((by - cy)*(ax - cx) + (cx - bx)*(ay - cy)); - double wb = ((cy - ay)*(gx[i] - cx) + (ax - cx)*(gy[j] - cy)) / - ((by - cy)*(ax - cx) + (cx - bx)*(ay - cy)); - double wc = 1 - wa - wb; - - image[i * nx + j] = (wa*realdata[ia] + wb*realdata[ib] + wc*realdata[ic])*dxy*dxy; + if (npoints[i * nx + j] > 1) { + image[i * nx + j] = binned_image[i * nx + j] / binned_weights[i * nx + j] * dxy * dxy; + } else { + int ia = d.triangles[which_triangle]; + double ax = x[ia]; + double ay = y[ia]; + int ib = d.triangles[which_triangle+1]; + double bx = x[ib]; + double by = y[ib]; + int ic = d.triangles[which_triangle+2]; + double cx = x[ic]; + double cy = y[ic]; + + double wa = ((by - cy)*(gx[i] - cx) + (cx - bx)*(gy[j] - cy)) / + ((by - cy)*(ax - cx) + (cx - bx)*(ay - cy)); + double wb = ((cy - ay)*(gx[i] - cx) + (ax - cx)*(gy[j] - cy)) / + ((by - cy)*(ax - cx) + (cx - bx)*(ay - cy)); + double wc = 1 - wa - wb; + + image[i * nx + j] = (wa*realdata[ia] + wb*realdata[ib] + wc*realdata[ic])*dxy*dxy; + } if (col_start_triangle == -1) { col_start_triangle = last_triangle; From 058b59aa95fcd91a15ff84e25d28bb18e6026b13 Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Fri, 12 Nov 2021 19:02:48 -0600 Subject: [PATCH 09/23] Get the orientation of sampleUnstructuredCPP correct --- src/galario.cpp | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/src/galario.cpp b/src/galario.cpp index 6800da4..d210a21 100644 --- a/src/galario.cpp +++ b/src/galario.cpp @@ -1481,7 +1481,12 @@ namespace galario { /** * Interpolate from an unstructured image onto a regular grid. */ -dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x, const dreal* y, const dreal* realdata, dreal v_origin) { +dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x, const dreal* realy, const dreal* realdata, dreal v_origin) { + // Flip y to get the orientation correct. + auto y = static_cast(malloc(sizeof(dreal)*ni)); + for (int i = 0; i < ni; i++) + y[i] = (-1*v_origin)*realy[i]; + // Set up the Delauney triangulation. std::vector coords; @@ -1569,20 +1574,20 @@ dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x TCLEAR(moo); #endif // Now loop through the pixels in the image pixels, find the triangle each point is in, and interpolate. - for (int i = 0; i < nx; i++) { + for (int i = 0; i < ny; i++) { if ((i > 0) and (col_start_triangle > -1)) { which_triangle = col_start_triangle; last_triangle = col_start_triangle; col_start_triangle = -1; } - for (int j = 0; j < ny; j++) { + for (int j = 0; j < nx; j++) { // Check whether the triangle is out of the triangulation. - if ((gx[i] > xmin) and (gx[i] < xmax) and (gy[j] > ymin) and (gy[j] < ymax)) { + if ((gx[j] > xmin) and (gx[j] < xmax) and (gy[i] > ymin) and (gy[i] < ymax)) { // Find which triangle this grid point is in. #ifdef GALARIO_TIMING TSTART(moo); #endif - which_triangle = find_triangle(&d, x, y, gx[i], gy[j], which_triangle, &last_triangle, &time); + which_triangle = find_triangle(&d, x, y, gx[j], gy[i], which_triangle, &last_triangle, &time); #ifdef GALARIO_TIMING TSTOP(moo); #endif @@ -1605,9 +1610,9 @@ dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x double cx = x[ic]; double cy = y[ic]; - double wa = ((by - cy)*(gx[i] - cx) + (cx - bx)*(gy[j] - cy)) / + double wa = ((by - cy)*(gx[j] - cx) + (cx - bx)*(gy[i] - cy)) / ((by - cy)*(ax - cx) + (cx - bx)*(ay - cy)); - double wb = ((cy - ay)*(gx[i] - cx) + (ax - cx)*(gy[j] - cy)) / + double wb = ((cy - ay)*(gx[j] - cx) + (ax - cx)*(gy[i] - cy)) / ((by - cy)*(ax - cx) + (cx - bx)*(ay - cy)); double wc = 1 - wa - wb; From 34aa7790d4295240e3c0ba3343227e9b4269d338 Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Sat, 13 Nov 2021 00:36:25 -0600 Subject: [PATCH 10/23] Fixes to get binned images working. Also added further timing tests and make sure to clean up. --- src/galario.cpp | 89 ++++++++++++++++++++++++++++++++++++------------- 1 file changed, 65 insertions(+), 24 deletions(-) diff --git a/src/galario.cpp b/src/galario.cpp index d210a21..d973444 100644 --- a/src/galario.cpp +++ b/src/galario.cpp @@ -1461,17 +1461,17 @@ int find_triangle_directedwalk(delaunator::Delaunator *d, const dreal *x, const * First try to find the triangle index using a directed walk, and if that fails switch to brute force. */ int find_triangle(delaunator::Delaunator *d, const dreal *x, const dreal *y, dreal gx, dreal gy, int start, int* last_good, double* time) { -#ifdef GALARIO_TIMING - TCREATE(boo); TCLEAR(boo); TSTART(boo); -#endif +//#ifdef GALARIO_TIMING +// TCREATE(boo); TCLEAR(boo); TSTART(boo); +//#endif int which_triangle = find_triangle_directedwalk(d, x, y, gx, gy, start, last_good, time); if (which_triangle == -2) { printf("Switching to brute force \n"); which_triangle = find_triangle_bruteforce(d, x, y, gx, gy); } -#ifdef GALARIO_TIMING - TSTOP(boo); *time += TGIVE(boo); -#endif +//#ifdef GALARIO_TIMING +// TSTOP(boo); *time += TGIVE(boo); +//#endif return which_triangle; } @@ -1489,6 +1489,9 @@ dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x // Set up the Delauney triangulation. +#ifdef GALARIO_TIMING + TCREATE(moo); TCLEAR(moo); TSTART(moo); +#endif std::vector coords; dreal xmin = std::numeric_limits::max(); dreal xmax = -std::numeric_limits::max(); @@ -1503,16 +1506,17 @@ dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x if (y[i] < ymin) ymin = y[i]; } -#ifdef GALARIO_TIMING - TCREATE(moo); TCLEAR(moo); TSTART(moo); -#endif delaunator::Delaunator d(coords); #ifdef GALARIO_TIMING TSTOP(moo); - printf("Time to triangulate %f \n", TGIVE(moo)); + printf(" Time to triangulate %f \n", TGIVE(moo)); #endif // For each triangle, calculate the centroid and which grid cell it falls in. + +#ifdef GALARIO_TIMING + TCLEAR(moo); TSTART(moo); +#endif auto tx = static_cast(malloc(sizeof(dreal)*d.triangles.size()/3)); auto ty = static_cast(malloc(sizeof(dreal)*d.triangles.size()/3)); auto tf = static_cast(malloc(sizeof(dreal)*d.triangles.size()/3)); @@ -1540,21 +1544,25 @@ dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x auto binned_weights = static_cast(malloc(sizeof(dreal)*nx*ny)); auto npoints = static_cast(malloc(sizeof(int)*nx*ny)); - for (int i = 0; i < nx; i++) { - for (int j = 0; j < ny; j++) { - binned_image[i,j] = 0; - binned_weights[i,j] = 0; - npoints[i,j] = 0; + for (int i = 0; i < ny; i++) { + for (int j = 0; j < nx; j++) { + binned_image[i*nx+j] = 0; + binned_weights[i*nx+j] = 0; + npoints[i*nx+j] = 0; } } for (int i = 0; i < d.triangles.size()/3; i++) { if ((itx[i] >= 0) and (itx[i] < nx) and (ity[i] >= 0) and (ity[i] < ny)) { - npoints[itx[i] * nx + ity[i]] += 1; - binned_image[itx[i] * nx + ity[i]] += tf[i] * ta[i]; - binned_weights[itx[i] * nx + ity[i]] += ta[i]; + npoints[ity[i] * nx + itx[i]] += 1; + binned_image[ity[i] * nx + itx[i]] += tf[i] * ta[i]; + binned_weights[ity[i] * nx + itx[i]] += ta[i]; } } +#ifdef GALARIO_TIMING + TSTOP(moo); + printf(" Time to create binned image: %f \n", TGIVE(moo)); +#endif // Create an image including the appropriate coordinates. auto gx = static_cast(malloc(sizeof(dreal)*nx)); @@ -1571,7 +1579,8 @@ dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x int col_start_triangle = -1; double time = 0.; #ifdef GALARIO_TIMING - TCLEAR(moo); + TCLEAR(moo); TSTART(moo); + TCREATE(boo); TCLEAR(boo); TSTART(boo); #endif // Now loop through the pixels in the image pixels, find the triangle each point is in, and interpolate. for (int i = 0; i < ny; i++) { @@ -1585,11 +1594,11 @@ dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x if ((gx[j] > xmin) and (gx[j] < xmax) and (gy[i] > ymin) and (gy[i] < ymax)) { // Find which triangle this grid point is in. #ifdef GALARIO_TIMING - TSTART(moo); + TSTART(boo); #endif which_triangle = find_triangle(&d, x, y, gx[j], gy[i], which_triangle, &last_triangle, &time); #ifdef GALARIO_TIMING - TSTOP(moo); + TSTOP(boo); #endif } else @@ -1631,12 +1640,37 @@ dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x } } #ifdef GALARIO_TIMING - printf("Time to calculate barycentric coords %f \n", time); - printf("Time to find triangles %f \n", TGIVE(moo)); + TSTOP(moo); + TSTOP(boo); + //printf("Time to calculate barycentric coords %f \n", time); + printf(" Time to interpolate onto the grid. %f \n", TGIVE(moo)); + printf(" Time to find triangles %f \n", TGIVE(boo)); #endif // Now copy to an image. +#ifdef GALARIO_TIMING + TCLEAR(moo); TSTART(moo); +#endif auto buffer = copy_input(nx, ny, image); +#ifdef GALARIO_TIMING + TSTOP(moo); + //printf("Time to calculate barycentric coords %f \n", time); + printf(" Time to copy to complex. %f \n", TGIVE(moo)); +#endif + + // Clean up. + +#ifdef GALARIO_TIMING + TCLEAR(moo); TSTART(moo); +#endif + galario_free(y); galario_free(tx); galario_free(ty); galario_free(tf); galario_free(ta); galario_free(itx); galario_free(ity); + galario_free(binned_image); galario_free(binned_weights); galario_free(npoints); + galario_free(gx); galario_free(gy); galario_free(image); +#ifdef GALARIO_TIMING + TSTOP(moo); + //printf("Time to calculate barycentric coords %f \n", time); + printf(" Time to clean up. %f \n", TGIVE(moo)); +#endif return buffer; } @@ -1727,10 +1761,17 @@ void sample_unstructured_image(const dreal* realx, const dreal* realy, int nx, i auto data = interpolate_to_image(nx, ny, ni, dxy, realx, realy, realdata, v_origin); t.Elapsed("sample_image::interpolate_to_grid"); #ifdef GALARIO_TIMING TSTOP(moo); - printf("Time to interpolate: %f \n", TGIVE(moo)); + printf("Total time to interpolate: %f \n", TGIVE(moo)); #endif +#ifdef GALARIO_TIMING + TCLEAR(moo); TSTART(moo); +#endif sample_h(nx, ny, data, v_origin, dRA, dDec, nd, duv, PA, u, v, vis_int); +#ifdef GALARIO_TIMING + TSTOP(moo); + printf("Total time to FFT and sample on (u,v): %f \n", TGIVE(moo)); +#endif t = CPUTimer(); galario_free(data); t.Elapsed("sample_image::free_data"); //#endif From 0a3275466e1fbc48245700de67c030876a43053a Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Sat, 13 Nov 2021 00:45:00 -0600 Subject: [PATCH 11/23] Clean up the Cython wrapper sampleUnstructuredImageCPP --- python/libcommon.pyx | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/python/libcommon.pyx b/python/libcommon.pyx index c4449be..589bcf7 100644 --- a/python/libcommon.pyx +++ b/python/libcommon.pyx @@ -591,8 +591,7 @@ def sampleUnstructuredImage(dreal[::1] x, dreal[::1] y, dreal[::1] image, def sampleUnstructuredImageCPP(dreal[::1] x, dreal[::1] y, dreal[::1] image, int nxy, dreal dxy, dreal[::1] u, dreal[::1] v, - dRA=0., dDec=0., PA=0., check=False, origin='upper', \ - dreal[::1] vol=None, return_weights=False): + dRA=0., dDec=0., PA=0., check=False, origin='upper'): """ Compute the synthetic visibilities of a model image at the specified (u, v) locations. From d9921b96426fe7577f3582db90a2abea13622416 Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Sat, 13 Nov 2021 09:00:45 -0600 Subject: [PATCH 12/23] Added OpenMP parallelization to interpolate_to_image --- src/galario.cpp | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/src/galario.cpp b/src/galario.cpp index d973444..3e53a00 100644 --- a/src/galario.cpp +++ b/src/galario.cpp @@ -1484,6 +1484,7 @@ namespace galario { dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x, const dreal* realy, const dreal* realdata, dreal v_origin) { // Flip y to get the orientation correct. auto y = static_cast(malloc(sizeof(dreal)*ni)); + #pragma omp parallel for for (int i = 0; i < ni; i++) y[i] = (-1*v_origin)*realy[i]; @@ -1528,6 +1529,7 @@ dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x dreal gx_max = 0.5*nx*dxy; dreal gy_max = 0.5*ny*dxy*v_origin; + #pragma omp parallel for for (int i = 0; i < d.triangles.size()/3; i++) { tx[i] = (x[d.triangles[3*i]] + x[d.triangles[3*i+1]] + x[d.triangles[3*i+2]]) / 3.; ty[i] = (y[d.triangles[3*i]] + y[d.triangles[3*i+1]] + y[d.triangles[3*i+2]]) / 3.; @@ -1544,6 +1546,7 @@ dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x auto binned_weights = static_cast(malloc(sizeof(dreal)*nx*ny)); auto npoints = static_cast(malloc(sizeof(int)*nx*ny)); + #pragma omp parallel for collapse(2) for (int i = 0; i < ny; i++) { for (int j = 0; j < nx; j++) { binned_image[i*nx+j] = 0; @@ -1553,6 +1556,8 @@ dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x } for (int i = 0; i < d.triangles.size()/3; i++) { + // Note: cant do this in parallel because two threads could access same + // grid cell at the same time. Locking made this very slow. if ((itx[i] >= 0) and (itx[i] < nx) and (ity[i] >= 0) and (ity[i] < ny)) { npoints[ity[i] * nx + itx[i]] += 1; binned_image[ity[i] * nx + itx[i]] += tf[i] * ta[i]; @@ -1569,20 +1574,25 @@ dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x auto gy = static_cast(malloc(sizeof(dreal)*nx)); auto image = static_cast(malloc(sizeof(dreal)*nx*ny)); + #pragma omp parallel for for (int i = 0; i < nx; i++) gx[i] = (0.5 - i * 1./nx) * nx * dxy; + #pragma omp parallel for for (int i = 0; i < ny; i++) gy[i] = (0.5 - i * 1./ny) * ny * dxy * v_origin; - int which_triangle = 0; - int last_triangle = 0; - int col_start_triangle = -1; - double time = 0.; #ifdef GALARIO_TIMING TCLEAR(moo); TSTART(moo); TCREATE(boo); TCLEAR(boo); TSTART(boo); #endif + #pragma omp parallel + { + int which_triangle = 0; + int last_triangle = 0; + int col_start_triangle = -1; + double time = 0.; // Now loop through the pixels in the image pixels, find the triangle each point is in, and interpolate. + #pragma omp for schedule(static) for (int i = 0; i < ny; i++) { if ((i > 0) and (col_start_triangle > -1)) { which_triangle = col_start_triangle; @@ -1639,6 +1649,7 @@ dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x } } } + } #ifdef GALARIO_TIMING TSTOP(moo); TSTOP(boo); From 0e70bd6b1fd6e4a3ae602556a6b3881b6107936c Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Sat, 13 Nov 2021 10:23:59 -0600 Subject: [PATCH 13/23] Use unordered map for binned image to save time and memory --- src/galario.cpp | 41 +++++++++++++++++++++++------------------ 1 file changed, 23 insertions(+), 18 deletions(-) diff --git a/src/galario.cpp b/src/galario.cpp index 3e53a00..22cc0b2 100644 --- a/src/galario.cpp +++ b/src/galario.cpp @@ -21,6 +21,7 @@ #include "galario_py.h" #include #include "timer.h" +#include // full function makes code hard to read #define tpb galario::threads() @@ -1541,29 +1542,34 @@ dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x ity[i] = trunc((ty[i] - gy_max) / (-dxy*v_origin) + 0.5); } - // Create an image that bins the triangles weighted by their area. - auto binned_image = static_cast(malloc(sizeof(dreal)*nx*ny)); - auto binned_weights = static_cast(malloc(sizeof(dreal)*nx*ny)); - auto npoints = static_cast(malloc(sizeof(int)*nx*ny)); - - #pragma omp parallel for collapse(2) - for (int i = 0; i < ny; i++) { - for (int j = 0; j < nx; j++) { - binned_image[i*nx+j] = 0; - binned_weights[i*nx+j] = 0; - npoints[i*nx+j] = 0; - } - } + std::unordered_map binned_image; + std::unordered_map binned_weights; + std::unordered_map npoints; for (int i = 0; i < d.triangles.size()/3; i++) { // Note: cant do this in parallel because two threads could access same // grid cell at the same time. Locking made this very slow. if ((itx[i] >= 0) and (itx[i] < nx) and (ity[i] >= 0) and (ity[i] < ny)) { - npoints[ity[i] * nx + itx[i]] += 1; - binned_image[ity[i] * nx + itx[i]] += tf[i] * ta[i]; - binned_weights[ity[i] * nx + itx[i]] += ta[i]; + if (npoints.find(ity[i] * nx + itx[i]) == npoints.end()) { + npoints[ity[i] * nx + itx[i]] = 1; + binned_image[ity[i] * nx + itx[i]] = tf[i] * ta[i]; + binned_weights[ity[i] * nx + itx[i]] = ta[i]; + } else { + npoints[ity[i] * nx + itx[i]] += 1; + binned_image[ity[i] * nx + itx[i]] += tf[i] * ta[i]; + binned_weights[ity[i] * nx + itx[i]] += ta[i]; + } } } + + std::unordered_map::iterator it = npoints.begin(); + while (it != npoints.end()) { + // Erase any places where npoints = 1 + if (it->first <= 1) + it = npoints.erase(it); + else + it++; + } #ifdef GALARIO_TIMING TSTOP(moo); printf(" Time to create binned image: %f \n", TGIVE(moo)); @@ -1616,7 +1622,7 @@ dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x // We've found the right triangle, now interpolate. if (which_triangle > -1) { - if (npoints[i * nx + j] > 1) { + if (npoints.find(i * nx + j) != npoints.end()) { image[i * nx + j] = binned_image[i * nx + j] / binned_weights[i * nx + j] * dxy * dxy; } else { int ia = d.triangles[which_triangle]; @@ -1675,7 +1681,6 @@ dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x TCLEAR(moo); TSTART(moo); #endif galario_free(y); galario_free(tx); galario_free(ty); galario_free(tf); galario_free(ta); galario_free(itx); galario_free(ity); - galario_free(binned_image); galario_free(binned_weights); galario_free(npoints); galario_free(gx); galario_free(gy); galario_free(image); #ifdef GALARIO_TIMING TSTOP(moo); From 13e46d7c19fb95feeff271c5a4bc293139ff3ad4 Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Sat, 18 Dec 2021 18:23:23 -0600 Subject: [PATCH 14/23] Add timing for the python version of unstructured transform. --- python/libcommon.pyx | 23 +++++++++++++++++++---- 1 file changed, 19 insertions(+), 4 deletions(-) diff --git a/python/libcommon.pyx b/python/libcommon.pyx index 589bcf7..21009d5 100644 --- a/python/libcommon.pyx +++ b/python/libcommon.pyx @@ -22,6 +22,7 @@ from cpython cimport PyObject, Py_INCREF # Numpy must be initialized. When using numpy from C or Cython you must # _always_ do that, or you will have segfaults +import time import numpy as np np.import_array() @@ -523,12 +524,17 @@ def sampleUnstructuredImage(dreal[::1] x, dreal[::1] y, dreal[::1] image, y[n] *= -1 # Use scipy to inerpolate onto a regular grid. + t1 = time. time() interp = LinearNDInterpolator(list(zip(x, y)), image, fill_value=0) + t2 = time. time() + print(" Time to triangulate "+str(t2-t1)) cdef dreal[:,::1] grid_x, grid_y grid_x_1D, grid_y_1D, grid_x, grid_y, _ = get_coords_meshgrid(nxy, nxy, \ dxy, origin=origin) cdef dreal[:,::1] new_image = interp(grid_x, grid_y) * dxy**2 + t3 = time.time() + print(" Time to do scipy interpolation "+str(t3-t2)) # In pixels where we oversample, average instead of interpolate in case the # intensity is varying quickly over the pixel. And use the volume of the @@ -539,6 +545,8 @@ def sampleUnstructuredImage(dreal[::1] x, dreal[::1] y, dreal[::1] image, cdef dreal[:,::1] binned_weights = np.zeros((nxy, nxy)) cdef int k, l, m cdef int nx = x.shape[0] + t4 = time.time() + print(" Time to create binned images "+str(t4 - t3)) i = ((x - grid_x_1D.max()) / -dxy + 0.5).astype(np.dtype('i')) if origin == "upper": @@ -569,6 +577,10 @@ def sampleUnstructuredImage(dreal[::1] x, dreal[::1] y, dreal[::1] image, new_image[l,m] = binned_image[l,m] / binned_weights[l,m] * \ dxy**2 + t5 = time.time() + print(" Time to incorporate binned image "+str(t5 - t4)) + print("Total time to interpolate "+str(t5 - t1)) + if origin == "upper": for n in range(y.size): y[n] *= -1 @@ -582,11 +594,14 @@ def sampleUnstructuredImage(dreal[::1] x, dreal[::1] y, dreal[::1] image, vis = np.zeros(len(u), dtype=complex_dtype) v_origin = set_v_origin(origin) cpp._sample_image(nxy, nxy, &new_image[0,0], v_origin, dRA, dDec, duv, PA, len(u), &u[0], &v[0], np.PyArray_DATA(vis)) + t6 = time.time() + print("Total time to FFT and sample on (u,v): "+str(t6 - t5)) - if return_weights: - return vis, vol - else: - return vis + #if return_weights: + # return vis, vol + #else: + # return vis + return vis def sampleUnstructuredImageCPP(dreal[::1] x, dreal[::1] y, dreal[::1] image, From 9c2c4524e8447dd1b0be925e629725227d9b1ef5 Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Sat, 18 Dec 2021 18:27:05 -0600 Subject: [PATCH 15/23] Add a cassert import since delaunator seems to need it... --- src/galario.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/galario.cpp b/src/galario.cpp index 22cc0b2..7b317ec 100644 --- a/src/galario.cpp +++ b/src/galario.cpp @@ -19,6 +19,7 @@ #include "galario.h" #include "galario_py.h" +#include #include #include "timer.h" #include From 45ba94da89ee46656bb36e82d204be9a3de61eee Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Wed, 9 Mar 2022 09:23:35 -0600 Subject: [PATCH 16/23] Put the bulk of the work from interpolate_to_image into three functions and one helper: - triangulate_h - bin_triangles_h - interpolate_or_bin_to_image_h - interpolate_on_triangle_h --- src/galario.cpp | 268 +++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 264 insertions(+), 4 deletions(-) diff --git a/src/galario.cpp b/src/galario.cpp index 7b317ec..5f367a9 100644 --- a/src/galario.cpp +++ b/src/galario.cpp @@ -1479,6 +1479,240 @@ int find_triangle(delaunator::Delaunator *d, const dreal *x, const dreal *y, dre } +/** + * Run the Delauney triangulation. + */ +delaunator::Delaunator triangulate_h(int ni, const dreal* x, const dreal* y, dreal v_origin) { + // Set up the Delauney triangulation. + +#ifdef GALARIO_TIMING + TCREATE(moo); TCLEAR(moo); TSTART(moo); +#endif + std::vector coords; + + dreal xmin = std::numeric_limits::max(); dreal xmax = -std::numeric_limits::max(); + dreal ymin = std::numeric_limits::max(); dreal ymax = -std::numeric_limits::max(); + for (int i=0; i < ni; i++) { + coords.push_back(x[i]); + coords.push_back(y[i]); + + if (x[i] > xmax) xmax = x[i]; + if (x[i] < xmin) xmin = x[i]; + if (y[i] > ymax) ymax = y[i]; + if (y[i] < ymin) ymin = y[i]; + } + + delaunator::Delaunator d(coords); +#ifdef GALARIO_TIMING + TSTOP(moo); + printf(" Time to triangulate %f \n", TGIVE(moo)); +#endif + + return d; +} + +/** + * For each triangle, calculate the centroid and which grid cell it falls in. + */ +void bin_triangles_h(int nx, int ny, dreal dxy, const dreal *x, const dreal *y, const dreal *realdata, delaunator::Delaunator &d, std::unordered_map &binned_image, + std::unordered_map &binned_weights, std::unordered_map &npoints, dreal v_origin) { +#ifdef GALARIO_TIMING + TCREATE(moo); TCLEAR(moo); TSTART(moo); +#endif + auto tx = static_cast(malloc(sizeof(dreal)*d.triangles.size()/3)); + auto ty = static_cast(malloc(sizeof(dreal)*d.triangles.size()/3)); + auto tf = static_cast(malloc(sizeof(dreal)*d.triangles.size()/3)); + auto ta = static_cast(malloc(sizeof(dreal)*d.triangles.size()/3)); + + auto itx = static_cast(malloc(sizeof(int)*d.triangles.size()/3)); + auto ity = static_cast(malloc(sizeof(int)*d.triangles.size()/3)); + + dreal gx_max = 0.5*nx*dxy; + dreal gy_max = 0.5*ny*dxy*v_origin; + + #pragma omp parallel for + for (int i = 0; i < d.triangles.size()/3; i++) { + tx[i] = (x[d.triangles[3*i]] + x[d.triangles[3*i+1]] + x[d.triangles[3*i+2]]) / 3.; + ty[i] = (y[d.triangles[3*i]] + y[d.triangles[3*i+1]] + y[d.triangles[3*i+2]]) / 3.; + tf[i] = (realdata[d.triangles[3*i]] + realdata[d.triangles[3*i+1]] + realdata[d.triangles[3*i+2]]) / 3.; + ta[i] = std::fabs((y[d.triangles[3*i+1]] - y[d.triangles[3*i]]) * (x[d.triangles[3*i+2]] - x[d.triangles[3*i+1]]) - + (x[d.triangles[3*i+1]] - x[d.triangles[3*i]]) * (y[d.triangles[3*i+2]] - y[d.triangles[3*i+1]])); + + itx[i] = trunc((tx[i] - gx_max) / (-dxy) + 0.5); + ity[i] = trunc((ty[i] - gy_max) / (-dxy*v_origin) + 0.5); + } + + for (int i = 0; i < d.triangles.size()/3; i++) { + // Note: cant do this in parallel because two threads could access same + // grid cell at the same time. Locking made this very slow. + if ((itx[i] >= 0) and (itx[i] < nx) and (ity[i] >= 0) and (ity[i] < ny)) { + if (npoints.find(ity[i] * nx + itx[i]) == npoints.end()) { + npoints[ity[i] * nx + itx[i]] = 1; + binned_image[ity[i] * nx + itx[i]] = tf[i] * ta[i]; + binned_weights[ity[i] * nx + itx[i]] = ta[i]; + } else { + npoints[ity[i] * nx + itx[i]] += 1; + binned_image[ity[i] * nx + itx[i]] += tf[i] * ta[i]; + binned_weights[ity[i] * nx + itx[i]] += ta[i]; + } + } + } + + std::unordered_map::iterator it = npoints.begin(); + while (it != npoints.end()) { + // Erase any places where npoints = 1 + if (it->first <= 1) + it = npoints.erase(it); + else + it++; + } + + free(tx); free(ty); free(tf); free(ta); free(itx); free(ity); + +#ifdef GALARIO_TIMING + TSTOP(moo); + printf(" Time to create binned image: %f \n", TGIVE(moo)); +#endif +} + +/** + * Do the interpolation onto a single point in a single triangle. + */ +double interpolate_on_triangle_h(delaunator::Delaunator &d, int which_triangle, const dreal *x, const dreal *y, const dreal *realdata, dreal gx, dreal gy) { + int ia = d.triangles[which_triangle]; + double ax = x[ia]; + double ay = y[ia]; + int ib = d.triangles[which_triangle+1]; + double bx = x[ib]; + double by = y[ib]; + int ic = d.triangles[which_triangle+2]; + double cx = x[ic]; + double cy = y[ic]; + + double wa = ((by - cy)*(gx - cx) + (cx - bx)*(gy - cy)) / + ((by - cy)*(ax - cx) + (cx - bx)*(ay - cy)); + double wb = ((cy - ay)*(gx - cx) + (ax - cx)*(gy - cy)) / + ((by - cy)*(ax - cx) + (cx - bx)*(ay - cy)); + double wc = 1 - wa - wb; + + return wa*realdata[ia] + wb*realdata[ib] + wc*realdata[ic]; +} + +/** + * Interpolate when the triangles are bigger than the grid cells, and use the binned image when triangles are smaller. + */ +dreal* interpolate_or_bin_to_image_h(int nx, int ny, int ni, dreal dxy, const dreal* x, const dreal* y, const dreal* realdata, dreal v_origin, + delaunator::Delaunator &d, std::unordered_map &binned_image, std::unordered_map &binned_weights, + std::unordered_map &npoints) { + + // Get the max and min x and y values from the triangulation. + dreal xmin = std::numeric_limits::max(); dreal xmax = -std::numeric_limits::max(); + dreal ymin = std::numeric_limits::max(); dreal ymax = -std::numeric_limits::max(); + for (int i=0; i < ni; i++) { + if (x[i] > xmax) xmax = x[i]; + if (x[i] < xmin) xmin = x[i]; + if (y[i] > ymax) ymax = y[i]; + if (y[i] < ymin) ymin = y[i]; + } + + // Create an image including the appropriate coordinates. + auto gx = static_cast(malloc(sizeof(dreal)*nx)); + auto gy = static_cast(malloc(sizeof(dreal)*nx)); + auto image = static_cast(malloc(sizeof(dreal)*nx*ny)); + + #pragma omp parallel for + for (int i = 0; i < nx; i++) + gx[i] = (0.5 - i * 1./nx) * nx * dxy; + #pragma omp parallel for + for (int i = 0; i < ny; i++) + gy[i] = (0.5 - i * 1./ny) * ny * dxy * v_origin; + +#ifdef GALARIO_TIMING + TCREATE(moo); TCLEAR(moo); TSTART(moo); + TCREATE(boo); TCLEAR(boo); TSTART(boo); +#endif + #pragma omp parallel + { + int which_triangle = 0; + int last_triangle = 0; + int col_start_triangle = -1; + double time = 0.; + + // Now loop through the pixels in the image pixels, find the triangle each point is in, and interpolate. + #pragma omp for schedule(static) + for (int i = 0; i < ny; i++) { + if ((i > 0) and (col_start_triangle > -1)) { + which_triangle = col_start_triangle; + last_triangle = col_start_triangle; + col_start_triangle = -1; + } + for (int j = 0; j < nx; j++) { + // Check whether the triangle is out of the triangulation. + if ((gx[j] > xmin) and (gx[j] < xmax) and (gy[i] > ymin) and (gy[i] < ymax)) { + // Find which triangle this grid point is in. +#ifdef GALARIO_TIMING + TSTART(boo); +#endif + which_triangle = find_triangle(&d, x, y, gx[j], gy[i], which_triangle, &last_triangle, &time); +#ifdef GALARIO_TIMING + TSTOP(boo); +#endif + } + else + which_triangle = -1; + + // We've found the right triangle, now interpolate. + if (which_triangle > -1) { + if (npoints.find(i * nx + j) != npoints.end()) { + image[i * nx + j] = binned_image[i * nx + j] / binned_weights[i * nx + j] * dxy * dxy; + } else { + /*int ia = d.triangles[which_triangle]; + double ax = x[ia]; + double ay = y[ia]; + int ib = d.triangles[which_triangle+1]; + double bx = x[ib]; + double by = y[ib]; + int ic = d.triangles[which_triangle+2]; + double cx = x[ic]; + double cy = y[ic]; + + double wa = ((by - cy)*(gx[j] - cx) + (cx - bx)*(gy[i] - cy)) / + ((by - cy)*(ax - cx) + (cx - bx)*(ay - cy)); + double wb = ((cy - ay)*(gx[j] - cx) + (ax - cx)*(gy[i] - cy)) / + ((by - cy)*(ax - cx) + (cx - bx)*(ay - cy)); + double wc = 1 - wa - wb; + + image[i * nx + j] = (wa*realdata[ia] + wb*realdata[ib] + wc*realdata[ic])*dxy*dxy;*/ + image[i * nx + j] = interpolate_on_triangle_h(d, which_triangle, x, y, realdata, gx[i], gy[i])*dxy*dxy; + } + + if (col_start_triangle == -1) { + col_start_triangle = last_triangle; + } + } + // If no triangle was found, the point is outside the area with data so set to 0. + else { + image[i * nx + j] = 0.; + which_triangle = last_triangle; + } + } + } + } + + // Clean up + free(gx); free(gy); + +#ifdef GALARIO_TIMING + TSTOP(moo); + TSTOP(boo); + //printf("Time to calculate barycentric coords %f \n", time); + printf(" Time to interpolate onto the grid. %f \n", TGIVE(moo)); + printf(" Time to find triangles %f \n", TGIVE(boo)); +#endif + + return image; +} + namespace galario { /** * Interpolate from an unstructured image onto a regular grid. @@ -1492,6 +1726,7 @@ dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x // Set up the Delauney triangulation. + /* #ifdef GALARIO_TIMING TCREATE(moo); TCLEAR(moo); TSTART(moo); #endif @@ -1514,9 +1749,13 @@ dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x TSTOP(moo); printf(" Time to triangulate %f \n", TGIVE(moo)); #endif + */ + + delaunator::Delaunator d = triangulate_h(ni, x, y, v_origin); // For each triangle, calculate the centroid and which grid cell it falls in. + /* #ifdef GALARIO_TIMING TCLEAR(moo); TSTART(moo); #endif @@ -1541,12 +1780,15 @@ dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x itx[i] = trunc((tx[i] - gx_max) / (-dxy) + 0.5); ity[i] = trunc((ty[i] - gy_max) / (-dxy*v_origin) + 0.5); - } + }*/ std::unordered_map binned_image; std::unordered_map binned_weights; std::unordered_map npoints; + bin_triangles_h(nx, ny, dxy, x, y, realdata, d, binned_image, binned_weights, npoints, v_origin); + + /* for (int i = 0; i < d.triangles.size()/3; i++) { // Note: cant do this in parallel because two threads could access same // grid cell at the same time. Locking made this very slow. @@ -1575,12 +1817,28 @@ dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x TSTOP(moo); printf(" Time to create binned image: %f \n", TGIVE(moo)); #endif + */ + + /* + // Get the max and min x and y values from the triangulation. + dreal xmin = std::numeric_limits::max(); dreal xmax = -std::numeric_limits::max(); + dreal ymin = std::numeric_limits::max(); dreal ymax = -std::numeric_limits::max(); + for (int i=0; i < ni; i++) { + if (x[i] > xmax) xmax = x[i]; + if (x[i] < xmin) xmin = x[i]; + if (y[i] > ymax) ymax = y[i]; + if (y[i] < ymin) ymin = y[i]; + } // Create an image including the appropriate coordinates. auto gx = static_cast(malloc(sizeof(dreal)*nx)); auto gy = static_cast(malloc(sizeof(dreal)*nx)); auto image = static_cast(malloc(sizeof(dreal)*nx*ny)); + */ + + auto image = interpolate_or_bin_to_image_h(nx, ny, ni, dxy, x, y, realdata, v_origin, d, binned_image, binned_weights, npoints); + /* #pragma omp parallel for for (int i = 0; i < nx; i++) gx[i] = (0.5 - i * 1./nx) * nx * dxy; @@ -1598,6 +1856,7 @@ dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x int last_triangle = 0; int col_start_triangle = -1; double time = 0.; + // Now loop through the pixels in the image pixels, find the triangle each point is in, and interpolate. #pragma omp for schedule(static) for (int i = 0; i < ny; i++) { @@ -1664,10 +1923,11 @@ dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x printf(" Time to interpolate onto the grid. %f \n", TGIVE(moo)); printf(" Time to find triangles %f \n", TGIVE(boo)); #endif + */ // Now copy to an image. #ifdef GALARIO_TIMING - TCLEAR(moo); TSTART(moo); + TCREATE(moo); TCLEAR(moo); TSTART(moo); #endif auto buffer = copy_input(nx, ny, image); #ifdef GALARIO_TIMING @@ -1681,8 +1941,8 @@ dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x #ifdef GALARIO_TIMING TCLEAR(moo); TSTART(moo); #endif - galario_free(y); galario_free(tx); galario_free(ty); galario_free(tf); galario_free(ta); galario_free(itx); galario_free(ity); - galario_free(gx); galario_free(gy); galario_free(image); + galario_free(y); //galario_free(tx); galario_free(ty); galario_free(tf); galario_free(ta); galario_free(itx); galario_free(ity); + /*galario_free(gx); galario_free(gy);*/ galario_free(image); #ifdef GALARIO_TIMING TSTOP(moo); //printf("Time to calculate barycentric coords %f \n", time); From 9e823fc5cbcbfddfb50562b3c49b1c839e317664 Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Wed, 9 Mar 2022 13:31:37 -0600 Subject: [PATCH 17/23] Renamed interpolate_to_image => unstructured_to_grid_h Also cleaned up some commented out code and some small rearranging. --- src/galario.cpp | 258 ++++------------------------------------------- src/galario.h | 1 - src/galario_py.h | 1 - 3 files changed, 18 insertions(+), 242 deletions(-) diff --git a/src/galario.cpp b/src/galario.cpp index 5f367a9..dcf6aca 100644 --- a/src/galario.cpp +++ b/src/galario.cpp @@ -1617,7 +1617,7 @@ dreal* interpolate_or_bin_to_image_h(int nx, int ny, int ni, dreal dxy, const dr // Create an image including the appropriate coordinates. auto gx = static_cast(malloc(sizeof(dreal)*nx)); - auto gy = static_cast(malloc(sizeof(dreal)*nx)); + auto gy = static_cast(malloc(sizeof(dreal)*ny)); auto image = static_cast(malloc(sizeof(dreal)*nx*ny)); #pragma omp parallel for @@ -1666,23 +1666,6 @@ dreal* interpolate_or_bin_to_image_h(int nx, int ny, int ni, dreal dxy, const dr if (npoints.find(i * nx + j) != npoints.end()) { image[i * nx + j] = binned_image[i * nx + j] / binned_weights[i * nx + j] * dxy * dxy; } else { - /*int ia = d.triangles[which_triangle]; - double ax = x[ia]; - double ay = y[ia]; - int ib = d.triangles[which_triangle+1]; - double bx = x[ib]; - double by = y[ib]; - int ic = d.triangles[which_triangle+2]; - double cx = x[ic]; - double cy = y[ic]; - - double wa = ((by - cy)*(gx[j] - cx) + (cx - bx)*(gy[i] - cy)) / - ((by - cy)*(ax - cx) + (cx - bx)*(ay - cy)); - double wb = ((cy - ay)*(gx[j] - cx) + (ax - cx)*(gy[i] - cy)) / - ((by - cy)*(ax - cx) + (cx - bx)*(ay - cy)); - double wc = 1 - wa - wb; - - image[i * nx + j] = (wa*realdata[ia] + wb*realdata[ib] + wc*realdata[ic])*dxy*dxy;*/ image[i * nx + j] = interpolate_on_triangle_h(d, which_triangle, x, y, realdata, gx[i], gy[i])*dxy*dxy; } @@ -1713,11 +1696,10 @@ dreal* interpolate_or_bin_to_image_h(int nx, int ny, int ni, dreal dxy, const dr return image; } -namespace galario { /** * Interpolate from an unstructured image onto a regular grid. */ -dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x, const dreal* realy, const dreal* realdata, dreal v_origin) { +dreal* unstructured_to_grid_h(int nx, int ny, int ni, dreal dxy, const dreal* x, const dreal* realy, const dreal* realdata, dreal v_origin) { // Flip y to get the orientation correct. auto y = static_cast(malloc(sizeof(dreal)*ni)); #pragma omp parallel for @@ -1725,236 +1707,22 @@ dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x y[i] = (-1*v_origin)*realy[i]; // Set up the Delauney triangulation. - - /* -#ifdef GALARIO_TIMING - TCREATE(moo); TCLEAR(moo); TSTART(moo); -#endif - std::vector coords; - - dreal xmin = std::numeric_limits::max(); dreal xmax = -std::numeric_limits::max(); - dreal ymin = std::numeric_limits::max(); dreal ymax = -std::numeric_limits::max(); - for (int i=0; i < ni; i++) { - coords.push_back(x[i]); - coords.push_back(y[i]); - - if (x[i] > xmax) xmax = x[i]; - if (x[i] < xmin) xmin = x[i]; - if (y[i] > ymax) ymax = y[i]; - if (y[i] < ymin) ymin = y[i]; - } - - delaunator::Delaunator d(coords); -#ifdef GALARIO_TIMING - TSTOP(moo); - printf(" Time to triangulate %f \n", TGIVE(moo)); -#endif - */ - delaunator::Delaunator d = triangulate_h(ni, x, y, v_origin); // For each triangle, calculate the centroid and which grid cell it falls in. - - /* -#ifdef GALARIO_TIMING - TCLEAR(moo); TSTART(moo); -#endif - auto tx = static_cast(malloc(sizeof(dreal)*d.triangles.size()/3)); - auto ty = static_cast(malloc(sizeof(dreal)*d.triangles.size()/3)); - auto tf = static_cast(malloc(sizeof(dreal)*d.triangles.size()/3)); - auto ta = static_cast(malloc(sizeof(dreal)*d.triangles.size()/3)); - - auto itx = static_cast(malloc(sizeof(int)*d.triangles.size()/3)); - auto ity = static_cast(malloc(sizeof(int)*d.triangles.size()/3)); - - dreal gx_max = 0.5*nx*dxy; - dreal gy_max = 0.5*ny*dxy*v_origin; - - #pragma omp parallel for - for (int i = 0; i < d.triangles.size()/3; i++) { - tx[i] = (x[d.triangles[3*i]] + x[d.triangles[3*i+1]] + x[d.triangles[3*i+2]]) / 3.; - ty[i] = (y[d.triangles[3*i]] + y[d.triangles[3*i+1]] + y[d.triangles[3*i+2]]) / 3.; - tf[i] = (realdata[d.triangles[3*i]] + realdata[d.triangles[3*i+1]] + realdata[d.triangles[3*i+2]]) / 3.; - ta[i] = std::fabs((y[d.triangles[3*i+1]] - y[d.triangles[3*i]]) * (x[d.triangles[3*i+2]] - x[d.triangles[3*i+1]]) - - (x[d.triangles[3*i+1]] - x[d.triangles[3*i]]) * (y[d.triangles[3*i+2]] - y[d.triangles[3*i+1]])); - - itx[i] = trunc((tx[i] - gx_max) / (-dxy) + 0.5); - ity[i] = trunc((ty[i] - gy_max) / (-dxy*v_origin) + 0.5); - }*/ - std::unordered_map binned_image; std::unordered_map binned_weights; std::unordered_map npoints; bin_triangles_h(nx, ny, dxy, x, y, realdata, d, binned_image, binned_weights, npoints, v_origin); - /* - for (int i = 0; i < d.triangles.size()/3; i++) { - // Note: cant do this in parallel because two threads could access same - // grid cell at the same time. Locking made this very slow. - if ((itx[i] >= 0) and (itx[i] < nx) and (ity[i] >= 0) and (ity[i] < ny)) { - if (npoints.find(ity[i] * nx + itx[i]) == npoints.end()) { - npoints[ity[i] * nx + itx[i]] = 1; - binned_image[ity[i] * nx + itx[i]] = tf[i] * ta[i]; - binned_weights[ity[i] * nx + itx[i]] = ta[i]; - } else { - npoints[ity[i] * nx + itx[i]] += 1; - binned_image[ity[i] * nx + itx[i]] += tf[i] * ta[i]; - binned_weights[ity[i] * nx + itx[i]] += ta[i]; - } - } - } - - std::unordered_map::iterator it = npoints.begin(); - while (it != npoints.end()) { - // Erase any places where npoints = 1 - if (it->first <= 1) - it = npoints.erase(it); - else - it++; - } -#ifdef GALARIO_TIMING - TSTOP(moo); - printf(" Time to create binned image: %f \n", TGIVE(moo)); -#endif - */ - - /* - // Get the max and min x and y values from the triangulation. - dreal xmin = std::numeric_limits::max(); dreal xmax = -std::numeric_limits::max(); - dreal ymin = std::numeric_limits::max(); dreal ymax = -std::numeric_limits::max(); - for (int i=0; i < ni; i++) { - if (x[i] > xmax) xmax = x[i]; - if (x[i] < xmin) xmin = x[i]; - if (y[i] > ymax) ymax = y[i]; - if (y[i] < ymin) ymin = y[i]; - } - - // Create an image including the appropriate coordinates. - auto gx = static_cast(malloc(sizeof(dreal)*nx)); - auto gy = static_cast(malloc(sizeof(dreal)*nx)); - auto image = static_cast(malloc(sizeof(dreal)*nx*ny)); - */ - + // Interpolate or bin, as appropriate to get to an image. auto image = interpolate_or_bin_to_image_h(nx, ny, ni, dxy, x, y, realdata, v_origin, d, binned_image, binned_weights, npoints); - /* - #pragma omp parallel for - for (int i = 0; i < nx; i++) - gx[i] = (0.5 - i * 1./nx) * nx * dxy; - #pragma omp parallel for - for (int i = 0; i < ny; i++) - gy[i] = (0.5 - i * 1./ny) * ny * dxy * v_origin; - -#ifdef GALARIO_TIMING - TCLEAR(moo); TSTART(moo); - TCREATE(boo); TCLEAR(boo); TSTART(boo); -#endif - #pragma omp parallel - { - int which_triangle = 0; - int last_triangle = 0; - int col_start_triangle = -1; - double time = 0.; - - // Now loop through the pixels in the image pixels, find the triangle each point is in, and interpolate. - #pragma omp for schedule(static) - for (int i = 0; i < ny; i++) { - if ((i > 0) and (col_start_triangle > -1)) { - which_triangle = col_start_triangle; - last_triangle = col_start_triangle; - col_start_triangle = -1; - } - for (int j = 0; j < nx; j++) { - // Check whether the triangle is out of the triangulation. - if ((gx[j] > xmin) and (gx[j] < xmax) and (gy[i] > ymin) and (gy[i] < ymax)) { - // Find which triangle this grid point is in. -#ifdef GALARIO_TIMING - TSTART(boo); -#endif - which_triangle = find_triangle(&d, x, y, gx[j], gy[i], which_triangle, &last_triangle, &time); -#ifdef GALARIO_TIMING - TSTOP(boo); -#endif - } - else - which_triangle = -1; - - // We've found the right triangle, now interpolate. - if (which_triangle > -1) { - if (npoints.find(i * nx + j) != npoints.end()) { - image[i * nx + j] = binned_image[i * nx + j] / binned_weights[i * nx + j] * dxy * dxy; - } else { - int ia = d.triangles[which_triangle]; - double ax = x[ia]; - double ay = y[ia]; - int ib = d.triangles[which_triangle+1]; - double bx = x[ib]; - double by = y[ib]; - int ic = d.triangles[which_triangle+2]; - double cx = x[ic]; - double cy = y[ic]; - - double wa = ((by - cy)*(gx[j] - cx) + (cx - bx)*(gy[i] - cy)) / - ((by - cy)*(ax - cx) + (cx - bx)*(ay - cy)); - double wb = ((cy - ay)*(gx[j] - cx) + (ax - cx)*(gy[i] - cy)) / - ((by - cy)*(ax - cx) + (cx - bx)*(ay - cy)); - double wc = 1 - wa - wb; - - image[i * nx + j] = (wa*realdata[ia] + wb*realdata[ib] + wc*realdata[ic])*dxy*dxy; - } - - if (col_start_triangle == -1) { - col_start_triangle = last_triangle; - } - } - // If no triangle was found, the point is outside the area with data so set to 0. - else { - image[i * nx + j] = 0.; - which_triangle = last_triangle; - } - } - } - } -#ifdef GALARIO_TIMING - TSTOP(moo); - TSTOP(boo); - //printf("Time to calculate barycentric coords %f \n", time); - printf(" Time to interpolate onto the grid. %f \n", TGIVE(moo)); - printf(" Time to find triangles %f \n", TGIVE(boo)); -#endif - */ - - // Now copy to an image. -#ifdef GALARIO_TIMING - TCREATE(moo); TCLEAR(moo); TSTART(moo); -#endif - auto buffer = copy_input(nx, ny, image); -#ifdef GALARIO_TIMING - TSTOP(moo); - //printf("Time to calculate barycentric coords %f \n", time); - printf(" Time to copy to complex. %f \n", TGIVE(moo)); -#endif - // Clean up. + free(y); -#ifdef GALARIO_TIMING - TCLEAR(moo); TSTART(moo); -#endif - galario_free(y); //galario_free(tx); galario_free(ty); galario_free(tf); galario_free(ta); galario_free(itx); galario_free(ity); - /*galario_free(gx); galario_free(gy);*/ galario_free(image); -#ifdef GALARIO_TIMING - TSTOP(moo); - //printf("Time to calculate barycentric coords %f \n", time); - printf(" Time to clean up. %f \n", TGIVE(moo)); -#endif - - return buffer; -} - -void* _interpolate_to_image(int nx, int ny, int ni, dreal dxy, void* x, void *y, void* realdata, dreal v_origin) { - return interpolate_to_image(nx, ny, ni, dxy, static_cast(x), static_cast(y), static_cast(realdata), v_origin); -} + return image; } @@ -2035,7 +1803,7 @@ void sample_unstructured_image(const dreal* realx, const dreal* realy, int nx, i #ifdef GALARIO_TIMING TCREATE(moo); TCLEAR(moo); TSTART(moo); #endif - auto data = interpolate_to_image(nx, ny, ni, dxy, realx, realy, realdata, v_origin); t.Elapsed("sample_image::interpolate_to_grid"); + auto data = unstructured_to_grid_h(nx, ny, ni, dxy, realx, realy, realdata, v_origin); t.Elapsed("sample_image::interpolate_to_grid"); #ifdef GALARIO_TIMING TSTOP(moo); printf("Total time to interpolate: %f \n", TGIVE(moo)); @@ -2044,13 +1812,23 @@ void sample_unstructured_image(const dreal* realx, const dreal* realy, int nx, i #ifdef GALARIO_TIMING TCLEAR(moo); TSTART(moo); #endif - sample_h(nx, ny, data, v_origin, dRA, dDec, nd, duv, PA, u, v, vis_int); + auto image = copy_input(nx, ny, data); +#ifdef GALARIO_TIMING + TSTOP(moo); + //printf("Time to calculate barycentric coords %f \n", time); + printf("Time to copy to complex. %f \n", TGIVE(moo)); +#endif + +#ifdef GALARIO_TIMING + TCLEAR(moo); TSTART(moo); +#endif + sample_h(nx, ny, image, v_origin, dRA, dDec, nd, duv, PA, u, v, vis_int); #ifdef GALARIO_TIMING TSTOP(moo); printf("Total time to FFT and sample on (u,v): %f \n", TGIVE(moo)); #endif - t = CPUTimer(); galario_free(data); t.Elapsed("sample_image::free_data"); + t = CPUTimer(); galario_free(data); galario_free(image); t.Elapsed("sample_image::free_data"); //#endif t_start.Elapsed("sample_image_tot"); } diff --git a/src/galario.h b/src/galario.h index f9ce326..68a00e5 100644 --- a/src/galario.h +++ b/src/galario.h @@ -37,7 +37,6 @@ void uv_rotate(dreal PA, dreal dRA, dreal dDec, dreal* dRArot, dreal* dDecrot, i /* Interface for the experts */ dcomplex* copy_input(int nx, int ny, const dreal* image); -dcomplex* interpolate_to_image(int nx, int ny, int ni, dreal dxy, const dreal* x, const dreal* y, const dreal* data, dreal v_origin); void galario_free(void*); void fft2d(int nx, int ny, dcomplex* image); void fftshift(int nx, int ny, dcomplex* image); diff --git a/src/galario_py.h b/src/galario_py.h index cf96db4..39a360c 100644 --- a/src/galario_py.h +++ b/src/galario_py.h @@ -39,7 +39,6 @@ void _uv_rotate(dreal PA, dreal dRA, dreal dDec, void* dRArot, void* dDecrot, in /* Interface for the experts */ void* _copy_input(int nx, int ny, void* realdata); -void* interpolate_to_image(int nx, int ny, int ni, dreal dxy, void* x, void* y, void* data, dreal v_origin); void _fft2d(int nx, int ny, void* data); void _fftshift(int nx, int ny, void* data); void _fftshift_axis0(int nx, int ncol, void* data); From 8629a61af5ed23065e489be86827a283ce27a018 Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Wed, 9 Mar 2022 14:42:01 -0600 Subject: [PATCH 18/23] Moved back to the built in galario timer. --- src/galario.cpp | 90 ++++++------------------------------------------- src/timer.h | 25 -------------- 2 files changed, 10 insertions(+), 105 deletions(-) delete mode 100644 src/timer.h diff --git a/src/galario.cpp b/src/galario.cpp index dcf6aca..a82568f 100644 --- a/src/galario.cpp +++ b/src/galario.cpp @@ -21,7 +21,6 @@ #include "galario_py.h" #include #include -#include "timer.h" #include // full function makes code hard to read @@ -1410,8 +1409,6 @@ int find_triangle_directedwalk(delaunator::Delaunator *d, const dreal *x, const int count = 0; dreal eps = 1.0e-3; bool found_triangle = false; - //TCREATE(moo); TCLEAR(moo); - //TSTART(moo); while (count < d->triangles.size() / (3*4)) { int ia = d->triangles[start]; double ax = x[ia]; @@ -1452,9 +1449,6 @@ int find_triangle_directedwalk(delaunator::Delaunator *d, const dreal *x, const count++; } - // TSTOP(moo); - //if (count > 0) printf("gx = %f, gy = %f, count = %d \n", gx, gy, count); - //*time += TGIVE(moo); return which_triangle; } @@ -1463,17 +1457,11 @@ int find_triangle_directedwalk(delaunator::Delaunator *d, const dreal *x, const * First try to find the triangle index using a directed walk, and if that fails switch to brute force. */ int find_triangle(delaunator::Delaunator *d, const dreal *x, const dreal *y, dreal gx, dreal gy, int start, int* last_good, double* time) { -//#ifdef GALARIO_TIMING -// TCREATE(boo); TCLEAR(boo); TSTART(boo); -//#endif int which_triangle = find_triangle_directedwalk(d, x, y, gx, gy, start, last_good, time); if (which_triangle == -2) { printf("Switching to brute force \n"); which_triangle = find_triangle_bruteforce(d, x, y, gx, gy); } -//#ifdef GALARIO_TIMING -// TSTOP(boo); *time += TGIVE(boo); -//#endif return which_triangle; } @@ -1485,9 +1473,6 @@ int find_triangle(delaunator::Delaunator *d, const dreal *x, const dreal *y, dre delaunator::Delaunator triangulate_h(int ni, const dreal* x, const dreal* y, dreal v_origin) { // Set up the Delauney triangulation. -#ifdef GALARIO_TIMING - TCREATE(moo); TCLEAR(moo); TSTART(moo); -#endif std::vector coords; dreal xmin = std::numeric_limits::max(); dreal xmax = -std::numeric_limits::max(); @@ -1503,10 +1488,6 @@ delaunator::Delaunator triangulate_h(int ni, const dreal* x, const dreal* y, dre } delaunator::Delaunator d(coords); -#ifdef GALARIO_TIMING - TSTOP(moo); - printf(" Time to triangulate %f \n", TGIVE(moo)); -#endif return d; } @@ -1516,9 +1497,6 @@ delaunator::Delaunator triangulate_h(int ni, const dreal* x, const dreal* y, dre */ void bin_triangles_h(int nx, int ny, dreal dxy, const dreal *x, const dreal *y, const dreal *realdata, delaunator::Delaunator &d, std::unordered_map &binned_image, std::unordered_map &binned_weights, std::unordered_map &npoints, dreal v_origin) { -#ifdef GALARIO_TIMING - TCREATE(moo); TCLEAR(moo); TSTART(moo); -#endif auto tx = static_cast(malloc(sizeof(dreal)*d.triangles.size()/3)); auto ty = static_cast(malloc(sizeof(dreal)*d.triangles.size()/3)); auto tf = static_cast(malloc(sizeof(dreal)*d.triangles.size()/3)); @@ -1568,11 +1546,6 @@ void bin_triangles_h(int nx, int ny, dreal dxy, const dreal *x, const dreal *y, } free(tx); free(ty); free(tf); free(ta); free(itx); free(ity); - -#ifdef GALARIO_TIMING - TSTOP(moo); - printf(" Time to create binned image: %f \n", TGIVE(moo)); -#endif } /** @@ -1627,10 +1600,6 @@ dreal* interpolate_or_bin_to_image_h(int nx, int ny, int ni, dreal dxy, const dr for (int i = 0; i < ny; i++) gy[i] = (0.5 - i * 1./ny) * ny * dxy * v_origin; -#ifdef GALARIO_TIMING - TCREATE(moo); TCLEAR(moo); TSTART(moo); - TCREATE(boo); TCLEAR(boo); TSTART(boo); -#endif #pragma omp parallel { int which_triangle = 0; @@ -1648,30 +1617,21 @@ dreal* interpolate_or_bin_to_image_h(int nx, int ny, int ni, dreal dxy, const dr } for (int j = 0; j < nx; j++) { // Check whether the triangle is out of the triangulation. - if ((gx[j] > xmin) and (gx[j] < xmax) and (gy[i] > ymin) and (gy[i] < ymax)) { + if ((gx[j] > xmin) and (gx[j] < xmax) and (gy[i] > ymin) and (gy[i] < ymax)) // Find which triangle this grid point is in. -#ifdef GALARIO_TIMING - TSTART(boo); -#endif which_triangle = find_triangle(&d, x, y, gx[j], gy[i], which_triangle, &last_triangle, &time); -#ifdef GALARIO_TIMING - TSTOP(boo); -#endif - } else which_triangle = -1; // We've found the right triangle, now interpolate. if (which_triangle > -1) { - if (npoints.find(i * nx + j) != npoints.end()) { + if (npoints.find(i * nx + j) != npoints.end()) image[i * nx + j] = binned_image[i * nx + j] / binned_weights[i * nx + j] * dxy * dxy; - } else { + else image[i * nx + j] = interpolate_on_triangle_h(d, which_triangle, x, y, realdata, gx[i], gy[i])*dxy*dxy; - } - if (col_start_triangle == -1) { + if (col_start_triangle == -1) col_start_triangle = last_triangle; - } } // If no triangle was found, the point is outside the area with data so set to 0. else { @@ -1685,14 +1645,6 @@ dreal* interpolate_or_bin_to_image_h(int nx, int ny, int ni, dreal dxy, const dr // Clean up free(gx); free(gy); -#ifdef GALARIO_TIMING - TSTOP(moo); - TSTOP(boo); - //printf("Time to calculate barycentric coords %f \n", time); - printf(" Time to interpolate onto the grid. %f \n", TGIVE(moo)); - printf(" Time to find triangles %f \n", TGIVE(boo)); -#endif - return image; } @@ -1707,17 +1659,17 @@ dreal* unstructured_to_grid_h(int nx, int ny, int ni, dreal dxy, const dreal* x, y[i] = (-1*v_origin)*realy[i]; // Set up the Delauney triangulation. - delaunator::Delaunator d = triangulate_h(ni, x, y, v_origin); + OPENMPTIME(delaunator::Delaunator d = triangulate_h(ni, x, y, v_origin), "unstructured_to_grid::triangulation"); // For each triangle, calculate the centroid and which grid cell it falls in. std::unordered_map binned_image; std::unordered_map binned_weights; std::unordered_map npoints; - bin_triangles_h(nx, ny, dxy, x, y, realdata, d, binned_image, binned_weights, npoints, v_origin); + OPENMPTIME(bin_triangles_h(nx, ny, dxy, x, y, realdata, d, binned_image, binned_weights, npoints, v_origin), "unstructured_to_grid::bin_trixels"); // Interpolate or bin, as appropriate to get to an image. - auto image = interpolate_or_bin_to_image_h(nx, ny, ni, dxy, x, y, realdata, v_origin, d, binned_image, binned_weights, npoints); + OPENMPTIME(auto image = interpolate_or_bin_to_image_h(nx, ny, ni, dxy, x, y, realdata, v_origin, d, binned_image, binned_weights, npoints), "unstructured_to_grid::generate_gridded_image"); // Clean up. free(y); @@ -1798,35 +1750,13 @@ void sample_unstructured_image(const dreal* realx, const dreal* realy, int nx, i t_total.Elapsed("sample_image_tot"); #else*/ - CPUTimer t; -#ifdef GALARIO_TIMING - TCREATE(moo); TCLEAR(moo); TSTART(moo); -#endif - auto data = unstructured_to_grid_h(nx, ny, ni, dxy, realx, realy, realdata, v_origin); t.Elapsed("sample_image::interpolate_to_grid"); -#ifdef GALARIO_TIMING - TSTOP(moo); - printf("Total time to interpolate: %f \n", TGIVE(moo)); -#endif + auto data = unstructured_to_grid_h(nx, ny, ni, dxy, realx, realy, realdata, v_origin); -#ifdef GALARIO_TIMING - TCLEAR(moo); TSTART(moo); -#endif - auto image = copy_input(nx, ny, data); -#ifdef GALARIO_TIMING - TSTOP(moo); - //printf("Time to calculate barycentric coords %f \n", time); - printf("Time to copy to complex. %f \n", TGIVE(moo)); -#endif + CPUTimer t; + auto image = copy_input(nx, ny, data); t.Elapsed("sample_image::copy_input"); -#ifdef GALARIO_TIMING - TCLEAR(moo); TSTART(moo); -#endif sample_h(nx, ny, image, v_origin, dRA, dDec, nd, duv, PA, u, v, vis_int); -#ifdef GALARIO_TIMING - TSTOP(moo); - printf("Total time to FFT and sample on (u,v): %f \n", TGIVE(moo)); -#endif t = CPUTimer(); galario_free(data); galario_free(image); t.Elapsed("sample_image::free_data"); //#endif diff --git a/src/timer.h b/src/timer.h deleted file mode 100644 index 033e8fe..0000000 --- a/src/timer.h +++ /dev/null @@ -1,25 +0,0 @@ -#include -#include - -struct timeval tuse; - -#define CPU_TIME gettimeofday( &tuse, (struct timezone *)0 ); - -#define TCREATE(x) \ - double __timerseconds##x=0; double __timerstartseconds##x=0; \ - double __timerusec##x=0; double __timerstartusec##x=0; - -#define TCLEAR(x) {__timerseconds##x = 0; __timerusec##x = 0; } - -#define TSTART(x) { CPU_TIME; \ - __timerstartseconds##x = tuse.tv_sec; \ - __timerstartusec##x = tuse.tv_usec; } - -#define TSTOP(x) { CPU_TIME; \ - __timerseconds##x += (tuse.tv_sec - __timerstartseconds##x); \ - __timerusec##x += (tuse.tv_usec - __timerstartusec##x); } - -#define TTIME(str,x) \ - printf("%s %6.6f seconds \n", str, __timerseconds##x+__timerusec##x*1.0e-6); - -#define TGIVE(x) (__timerseconds##x+__timerusec##x*1.0e-6) From 21c95e061c230096e160f916a7a0a22692e6cdb9 Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Wed, 9 Mar 2022 14:46:36 -0600 Subject: [PATCH 19/23] Moved functions for finding triangles to have names with _h --- src/galario.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/galario.cpp b/src/galario.cpp index a82568f..25ac9fd 100644 --- a/src/galario.cpp +++ b/src/galario.cpp @@ -1357,7 +1357,7 @@ void sample_h(int nx, int ny, dcomplex* data, const dreal v_origin, dreal dRA, d /** * Find the index of the triangle that a point is in using brute force. */ -int find_triangle_bruteforce(delaunator::Delaunator *d, const dreal *x, const dreal *y, dreal gx, dreal gy) { +int find_triangle_bruteforce_h(delaunator::Delaunator *d, const dreal *x, const dreal *y, dreal gx, dreal gy) { bool found_triangle = false; int which_triangle = -1; @@ -1404,7 +1404,7 @@ int find_triangle_bruteforce(delaunator::Delaunator *d, const dreal *x, const dr /** * Find which triangle a point is in using a directed walk. */ -int find_triangle_directedwalk(delaunator::Delaunator *d, const dreal *x, const dreal *y, dreal gx, dreal gy, int start, int* last_good, double *time) { +int find_triangle_directedwalk_h(delaunator::Delaunator *d, const dreal *x, const dreal *y, dreal gx, dreal gy, int start, int* last_good, double *time) { int which_triangle = -2; int count = 0; dreal eps = 1.0e-3; @@ -1456,11 +1456,11 @@ int find_triangle_directedwalk(delaunator::Delaunator *d, const dreal *x, const /** * First try to find the triangle index using a directed walk, and if that fails switch to brute force. */ -int find_triangle(delaunator::Delaunator *d, const dreal *x, const dreal *y, dreal gx, dreal gy, int start, int* last_good, double* time) { - int which_triangle = find_triangle_directedwalk(d, x, y, gx, gy, start, last_good, time); +int find_triangle_h(delaunator::Delaunator *d, const dreal *x, const dreal *y, dreal gx, dreal gy, int start, int* last_good, double* time) { + int which_triangle = find_triangle_directedwalk_h(d, x, y, gx, gy, start, last_good, time); if (which_triangle == -2) { printf("Switching to brute force \n"); - which_triangle = find_triangle_bruteforce(d, x, y, gx, gy); + which_triangle = find_triangle_bruteforce_h(d, x, y, gx, gy); } return which_triangle; @@ -1619,7 +1619,7 @@ dreal* interpolate_or_bin_to_image_h(int nx, int ny, int ni, dreal dxy, const dr // Check whether the triangle is out of the triangulation. if ((gx[j] > xmin) and (gx[j] < xmax) and (gy[i] > ymin) and (gy[i] < ymax)) // Find which triangle this grid point is in. - which_triangle = find_triangle(&d, x, y, gx[j], gy[i], which_triangle, &last_triangle, &time); + which_triangle = find_triangle_h(&d, x, y, gx[j], gy[i], which_triangle, &last_triangle, &time); else which_triangle = -1; From 4de775f3caa93a9c017cd01cf1b15ece0379de4b Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Wed, 9 Mar 2022 15:12:18 -0600 Subject: [PATCH 20/23] Removed old Cython version of sampleUnstructured --- python/libcommon.pyx | 159 +------------------------------------------ 1 file changed, 1 insertion(+), 158 deletions(-) diff --git a/python/libcommon.pyx b/python/libcommon.pyx index 21009d5..8fa6ae8 100644 --- a/python/libcommon.pyx +++ b/python/libcommon.pyx @@ -37,7 +37,7 @@ __all__ = ['arcsec', 'deg', 'cgs_to_Jy', 'pc', 'au', '_init', '_cleanup', 'set_v_origin', 'ngpus', 'use_gpu', 'threads', 'check_obs', 'check_image_size', 'get_image_size', - 'sampleImage', 'sampleUnstructuredImage', 'sampleUnstructuredImageCPP', 'sampleProfile', 'chi2Image', 'chi2Profile', + 'sampleImage', 'sampleUnstructuredImage', 'sampleProfile', 'chi2Image', 'chi2Profile', 'get_coords_meshgrid', 'sweep', 'uv_rotate', 'interpolate', 'apply_phase_vis', 'reduce_chi2', '_fft2d', '_fftshift', '_fftshift_axis0'] @@ -448,163 +448,6 @@ def sampleImage(dreal[:,::1] image, dxy, dreal[::1] u, dreal[::1] v, def sampleUnstructuredImage(dreal[::1] x, dreal[::1] y, dreal[::1] image, - int nxy, dreal dxy, dreal[::1] u, dreal[::1] v, - dRA=0., dDec=0., PA=0., check=False, origin='upper', \ - dreal[::1] vol=None, return_weights=False): - """ - Compute the synthetic visibilities of a model image at the specified (u, v) locations. - - The 2D surface brightness in `image` is Fourier transformed and sampled in the - (u, v) locations given in the `u` and `v` arrays. - - Typical call signature:: - - vis = sampleImage(image, dxy, u, v, dRA=0, dDec=0, PA=0, check=False, origin='upper') - - Parameters - ---------- - x : 1D array_like, float - List of x coordinates at which intensities are known. - **units**: rad - y : 1D array_like, float - List of y coordinates at which intensities are known. - **units**: rad - image : 1D array_like, float - Array containing the surface brightness of the model. - Assume the x-axis (R.A.) increases from right (West) to left (East) - and the y-axis (Dec.) increases from bottom (South) to top (North). - `nxy` must be even. - **units**: Jy/st - nxy : int - Number of pixels to use for the interpolated gridded image. - dxy : float - Size of the image cell in the interpolated image, assumed equal in both x and y direction. - **units**: rad - u : array_like, float - u coordinate of the visibility points where the FT has to be sampled. - **units**: wavelength - v : array_like, float - v coordinate of the visibility points where the FT has to be sampled. - The length of v must be equal to the length of u. - **units**: wavelength - dRA : float, optional - R.A. offset w.r.t. the phase center by which the image is translated. - If dRA > 0 translate the image towards the left (East). Default is 0. - **units**: rad - dDec : float, optional - Dec. offset w.r.t. the phase center by which the image is translated. - If dDec > 0 translate the image towards the top (North). Default is 0. - **units**: rad - PA : float, optional - Position Angle, defined East of North. Default is 0. - **units**: rad - check : bool, optional - If True, check whether `image` and `dxy` satisfy Nyquist criterion for - computing the synthetic visibilities in the (u, v) locations provided. - Additionally check that the (u, v) points fall in the image to avoid - segmentation violations. Default is False since the check might take - time. For executions where speed is important, set to False. - origin : ['upper' | 'lower'], optional - Set the [0,0] pixel index of the matrix in the upper left or lower left corner of the axes. - It follows the same convention as in matplotlib `matshow` and `imshow` commands. - Declination axis and the matrix y axis are parallel for `origin='lower'`, anti-parallel for `origin='upper'`. - The central pixel corresponding to the (RA, Dec) = (0, 0) is always [Nxy/2, Nxy/2]. - For more details see the Technical Requirements page in the online docs. - - Returns - ------- - vis : array_like, complex - Synthetic visibilities sampled in the (u, v) locations given in `u` and `v`. - **units**: Jy - - """ - - if origin == "upper": - for n in range(y.size): - y[n] *= -1 - - # Use scipy to inerpolate onto a regular grid. - t1 = time. time() - interp = LinearNDInterpolator(list(zip(x, y)), image, fill_value=0) - t2 = time. time() - print(" Time to triangulate "+str(t2-t1)) - - cdef dreal[:,::1] grid_x, grid_y - grid_x_1D, grid_y_1D, grid_x, grid_y, _ = get_coords_meshgrid(nxy, nxy, \ - dxy, origin=origin) - cdef dreal[:,::1] new_image = interp(grid_x, grid_y) * dxy**2 - t3 = time.time() - print(" Time to do scipy interpolation "+str(t3-t2)) - - # In pixels where we oversample, average instead of interpolate in case the - # intensity is varying quickly over the pixel. And use the volume of the - # associated Voronoi cell to weight each point being averaged. - cdef int[::1] i, j - cdef int[:,::1] npoints = np.zeros((nxy, nxy), dtype=np.dtype('i')) - cdef dreal[:,::1] binned_image = np.zeros((nxy, nxy)) - cdef dreal[:,::1] binned_weights = np.zeros((nxy, nxy)) - cdef int k, l, m - cdef int nx = x.shape[0] - t4 = time.time() - print(" Time to create binned images "+str(t4 - t3)) - - i = ((x - grid_x_1D.max()) / -dxy + 0.5).astype(np.dtype('i')) - if origin == "upper": - j = ((y - grid_y_1D.max()) / -dxy + 0.5).astype(np.dtype('i')) - elif origin == "lower": - j = ((y - grid_y_1D.min()) / dxy + 0.5).astype(np.dtype('i')) - - if vol is None: - vor = Voronoi(list(zip(x, y))) - vol = np.zeros(vor.npoints) - for k, reg_num in enumerate(vor.point_region): - indices = vor.regions[reg_num] - if -1 in indices: - vol[k] = np.inf - else: - vol[k] = ConvexHull(vor.vertices[indices]).volume - - with nogil: - for k in range(nx): - if j[k] >= 0 and j[k] < nxy and i[k] >= 0 and i[k] < nxy: - npoints[j[k],i[k]] += 1 - binned_image[j[k],i[k]] += image[k]*vol[k] - binned_weights[j[k],i[k]] += vol[k] - - for l in range(nxy): - for m in range(nxy): - if npoints[l,m] > 1: - new_image[l,m] = binned_image[l,m] / binned_weights[l,m] * \ - dxy**2 - - t5 = time.time() - print(" Time to incorporate binned image "+str(t5 - t4)) - print("Total time to interpolate "+str(t5 - t1)) - - if origin == "upper": - for n in range(y.size): - y[n] *= -1 - - # Now pick back up with what is typically done for regular grids. - duv = 1 / (dxy*nxy) - - if check: - check_image_size(u, v, nxy, dxy, duv) - - vis = np.zeros(len(u), dtype=complex_dtype) - v_origin = set_v_origin(origin) - cpp._sample_image(nxy, nxy, &new_image[0,0], v_origin, dRA, dDec, duv, PA, len(u), &u[0], &v[0], np.PyArray_DATA(vis)) - t6 = time.time() - print("Total time to FFT and sample on (u,v): "+str(t6 - t5)) - - #if return_weights: - # return vis, vol - #else: - # return vis - return vis - - -def sampleUnstructuredImageCPP(dreal[::1] x, dreal[::1] y, dreal[::1] image, int nxy, dreal dxy, dreal[::1] u, dreal[::1] v, dRA=0., dDec=0., PA=0., check=False, origin='upper'): """ From c329898cefe8958c00aaefca6e2e9f8dbf41fe2c Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Wed, 9 Mar 2022 19:20:43 -0600 Subject: [PATCH 21/23] Fixed an orientation issue. --- src/galario.cpp | 13 ++----------- 1 file changed, 2 insertions(+), 11 deletions(-) diff --git a/src/galario.cpp b/src/galario.cpp index 25ac9fd..4880cb3 100644 --- a/src/galario.cpp +++ b/src/galario.cpp @@ -1628,7 +1628,7 @@ dreal* interpolate_or_bin_to_image_h(int nx, int ny, int ni, dreal dxy, const dr if (npoints.find(i * nx + j) != npoints.end()) image[i * nx + j] = binned_image[i * nx + j] / binned_weights[i * nx + j] * dxy * dxy; else - image[i * nx + j] = interpolate_on_triangle_h(d, which_triangle, x, y, realdata, gx[i], gy[i])*dxy*dxy; + image[i * nx + j] = interpolate_on_triangle_h(d, which_triangle, x, y, realdata, gx[j], gy[i])*dxy*dxy; if (col_start_triangle == -1) col_start_triangle = last_triangle; @@ -1651,13 +1651,7 @@ dreal* interpolate_or_bin_to_image_h(int nx, int ny, int ni, dreal dxy, const dr /** * Interpolate from an unstructured image onto a regular grid. */ -dreal* unstructured_to_grid_h(int nx, int ny, int ni, dreal dxy, const dreal* x, const dreal* realy, const dreal* realdata, dreal v_origin) { - // Flip y to get the orientation correct. - auto y = static_cast(malloc(sizeof(dreal)*ni)); - #pragma omp parallel for - for (int i = 0; i < ni; i++) - y[i] = (-1*v_origin)*realy[i]; - +dreal* unstructured_to_grid_h(int nx, int ny, int ni, dreal dxy, const dreal* x, const dreal* y, const dreal* realdata, dreal v_origin) { // Set up the Delauney triangulation. OPENMPTIME(delaunator::Delaunator d = triangulate_h(ni, x, y, v_origin), "unstructured_to_grid::triangulation"); @@ -1671,9 +1665,6 @@ dreal* unstructured_to_grid_h(int nx, int ny, int ni, dreal dxy, const dreal* x, // Interpolate or bin, as appropriate to get to an image. OPENMPTIME(auto image = interpolate_or_bin_to_image_h(nx, ny, ni, dxy, x, y, realdata, v_origin, d, binned_image, binned_weights, npoints), "unstructured_to_grid::generate_gridded_image"); - // Clean up. - free(y); - return image; } From b5978431759fd1d873d1f9bb497124558ed1e91c Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Wed, 9 Mar 2022 19:22:17 -0600 Subject: [PATCH 22/23] Added some test scripts. --- test/test_centering.py | 79 ++++++++++++++++++++++++++++++ test/test_circle.py | 77 +++++++++++++++++++++++++++++ test/test_orientation.py | 101 +++++++++++++++++++++++++++++++++++++++ 3 files changed, 257 insertions(+) create mode 100755 test/test_centering.py create mode 100755 test/test_circle.py create mode 100755 test/test_orientation.py diff --git a/test/test_centering.py b/test/test_centering.py new file mode 100755 index 0000000..bf751e4 --- /dev/null +++ b/test/test_centering.py @@ -0,0 +1,79 @@ +#!/usr/bin/env python3 + +from galario import double +import matplotlib.pyplot as plt +import matplotlib.tri as tri +import pyDOE +import numpy + +# Make an image. + +grid = pyDOE.lhs(2, samples=400) + +r = grid[:,0] * 2 +phi = grid[:,1] * 2*numpy.pi + +grid = pyDOE.lhs(2, samples=2000) + +r = numpy.hstack((r, grid[:,0]*0.04 + 0.98)) +phi = numpy.hstack((phi, grid[:,1]*2*numpy.pi)) + +x = r * numpy.cos(phi) +y = r * numpy.sin(phi) + +flux = numpy.where(r < 1., 1., 0.) + +# Do the Fourier transform with TrIFT + +u, v = numpy.meshgrid(numpy.linspace(-3.,3.,100),numpy.linspace(-3.,3.,100)) + +u = u.reshape((u.size,)) +v = v.reshape((v.size,)) + +vis = double.sampleUnstructuredImage(x, y, flux, 4096, 0.02, u, v, 0.5, 0.25) + +# Now shift the image manually. + +x += 0.5 +y += 0.25 + +vvis = double.sampleUnstructuredImage(x, y, flux, 4096, 0.02, u, v, 0., 0.) + +# Plot the image to make sure we did the correct shifting. + +triang = tri.Triangulation(x, y) + +fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(4,4)) + +ax.tripcolor(triang, flux, "ko-") +ax.triplot(triang, "k.-", linewidth=0.1, markersize=0.1) + +ax.set_aspect("equal") + +ax.set_xlim(1.6,-1.6) +ax.set_ylim(-1.6,1.6) + +ax.set_xlabel("x", fontsize=14) +ax.set_ylabel("y", fontsize=14) + +ax.tick_params(labelsize=14) + +plt.show() + +# Finally, plot the visibilities. + +fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(8,8)) + +ax[0,0].scatter(u, v, c=vis.real, marker=".") +ax[0,1].scatter(u, v, c=vis.imag, marker=".") + +ax[1,0].scatter(u, v, c=vvis.real, marker=".") +ax[1,1].scatter(u, v, c=vvis.imag, marker=".") + +for a in ax.flatten(): + a.set_xlabel("u", fontsize=14) + a.set_ylabel("v", fontsize=14) + + a.tick_params(labelsize=14) + +plt.show() diff --git a/test/test_circle.py b/test/test_circle.py new file mode 100755 index 0000000..880c57c --- /dev/null +++ b/test/test_circle.py @@ -0,0 +1,77 @@ +#!/usr/bin/env python3 + +import matplotlib.pyplot as plt +import matplotlib.tri as tri +import scipy.special +import pyDOE +import numpy +import time + +from galario import double + +# Make an image. + +grid = pyDOE.lhs(2, samples=100) + +r = grid[:,0] * 2 +phi = grid[:,1] * 2*numpy.pi + +grid = pyDOE.lhs(2, samples=2000) + +r = numpy.hstack((r, grid[:,0]*0.04 + 0.98)) +phi = numpy.hstack((phi, grid[:,1]*2*numpy.pi)) + +x = r * numpy.cos(phi) +y = r * numpy.sin(phi) + +flux = numpy.where(r < 1., 1., 0.) + +# Plot the image. + +triang = tri.Triangulation(x, y) + +plt.tripcolor(triang, flux, "ko-") +plt.triplot(triang, "k.-", linewidth=0.1, markersize=0.1) + +plt.axes().set_aspect("equal") + +plt.xlim(-1.1,1.1) +plt.ylim(-1.1,1.1) + +plt.xlabel("x", fontsize=14) +plt.ylabel("y", fontsize=14) + +plt.axes().tick_params(labelsize=14) + +plt.show() + +# Do the Fourier transform with TrIFT + +u = numpy.linspace(0.001,10.,1000) +v = numpy.repeat(0., 1000) + +t1 = time.time() +vis = double.sampleUnstructuredImage(x, y, flux, 4096, 0.02, u, v, 0.25, 0.25) +t2 = time.time() +print(t2 - t1) + +# Calculate the analytic result. + +vis_analytic = scipy.special.jv(1, 2*numpy.pi*u) / u * numpy.exp(2*numpy.pi*\ + 1j*(0.25*u + 0.25*v)) + +# Finally, plot the visibilities. + +plt.plot(u, vis.real, "k.-", label="Unstructured Fourier Transform") +plt.plot(u, vis_analytic.real, "r-", label="Analytic Solution") + +plt.xlabel("u", fontsize=14) +plt.ylabel("Real Component", fontsize=14) + +plt.legend(fontsize=14) + +plt.axes().tick_params(labelsize=14) + +plt.subplots_adjust(left=0.17, right=0.95, top=0.99) + +plt.show() diff --git a/test/test_orientation.py b/test/test_orientation.py new file mode 100755 index 0000000..731a028 --- /dev/null +++ b/test/test_orientation.py @@ -0,0 +1,101 @@ +#!/usr/bin/env python3 + +from galario import double +import matplotlib.pyplot as plt +import matplotlib.tri as tri +import pyDOE +import numpy + +# Make an image. + +grid = pyDOE.lhs(2, samples=400) + +r = grid[:,0] * 2 +phi = grid[:,1] * 2*numpy.pi + +grid = pyDOE.lhs(2, samples=2000) + +r = numpy.hstack((r, grid[:,0]*0.04 + 0.98)) +phi = numpy.hstack((phi, grid[:,1]*2*numpy.pi)) + +x = r * numpy.cos(phi) * numpy.cos(numpy.pi/3) +y = r * numpy.sin(phi) + +flux = numpy.where(r < 1., y - y.min(), 0.) + +pa = numpy.pi/4 + +xp = x * numpy.cos(-pa) - y * numpy.sin(-pa) +yp = x * numpy.sin(-pa) + y * numpy.cos(-pa) + +x = xp +y = yp + +# Also make a traditional image to compare with. + +xx, yy = numpy.meshgrid(numpy.linspace(15.,-15.,1024, endpoint=False), \ + numpy.linspace(15.,-15.,1024, endpoint=False)) + +xp = xx * numpy.cos(pa) - yy * numpy.sin(pa) +yp = xx * numpy.sin(pa) + yy * numpy.cos(pa) + +rr = numpy.sqrt((xp/numpy.cos(numpy.pi/3))**2 + yp**2) + +fflux = numpy.where(rr < 1., yp - y.min(), 0.) + +# Plot the image. + +triang = tri.Triangulation(x, y) + +fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(8,4)) + +ax[0].tripcolor(triang, flux, "ko-") +ax[0].triplot(triang, "k.-", linewidth=0.1, markersize=0.1) + +ax[1].imshow(fflux, interpolation="nearest") + +for i in range(1): + ax[i].set_aspect("equal") + + ax[i].set_xlim(1.1,-1.1) + ax[i].set_ylim(-1.1,1.1) + + ax[i].set_xlabel("x", fontsize=14) + ax[i].set_ylabel("y", fontsize=14) + + ax[i].tick_params(labelsize=14) + +plt.show() + +# Do the Fourier transform with TrIFT + +u, v = numpy.meshgrid(numpy.linspace(-3.,3.,100),numpy.linspace(-3.,3.,100)) + +u = u.reshape((u.size,)) +v = v.reshape((v.size,)) + +vis = double.sampleUnstructuredImage(x, y, flux, 1024*4, 0.025/4, u, v, 0., 0.) + +# Do the Fourier transform with GALARIO. + +dxy = abs(xx[0,1] - xx[0,0]) + +vvis = double.sampleImage(fflux, dxy, u, v) + +# Finally, plot the visibilities. + +fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(8,8)) + +ax[0,0].scatter(u, v, c=vis.real/vis.real.max(), marker=".") +ax[0,1].scatter(u, v, c=vis.imag/vis.imag.max(), marker=".") + +ax[1,0].scatter(u, v, c=vvis.real/vvis.real.max(), marker=".") +ax[1,1].scatter(u, v, c=vvis.imag/vvis.imag.max(), marker=".") + +for a in ax.flatten(): + a.set_xlabel("u", fontsize=14) + a.set_ylabel("v", fontsize=14) + + a.tick_params(labelsize=14) + +plt.show() From b66965aec99617d443cd763fb742dbda4845bc55 Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Thu, 10 Mar 2022 10:20:57 -0600 Subject: [PATCH 23/23] Added a chi2_unstructured_image function. --- python/galario_defs.pxd | 1 + python/libcommon.pyx | 104 ++++++++++++++++++++++++++++++++++++++++ src/galario.cpp | 59 +++++++++++++++++++++++ src/galario.h | 1 + src/galario_py.h | 1 + 5 files changed, 166 insertions(+) diff --git a/python/galario_defs.pxd b/python/galario_defs.pxd index 305d7e3..c214325 100644 --- a/python/galario_defs.pxd +++ b/python/galario_defs.pxd @@ -26,6 +26,7 @@ cdef extern from "galario_py.h" namespace "galario": void _sample_unstructured_image(void* x, void* y, int nx, int ny, dreal dxy, int ni, void* image, dreal v_origin, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, void* u, void* v, void* vis) except + dreal _chi2_profile(int nr, void* intensity, dreal Rmin, dreal dR, dreal dxy, int nxy, dreal inc, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, void* u, void* v, void* vis_obs_re, void* vis_obs_im, void* vis_obs_w) except + dreal _chi2_image(int nx, int ny, void* image, dreal v_origin, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, void* u, void* v, void* vis_obs_re, void* vis_obs_im, void* vis_obs_w) except + + dreal _chi2_unstructured_image(void* x, void* y, int nx, int ny, dreal dxy, int ni, void* data, const dreal v_origin, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, void* u, void* v, void* vis_obs_re, void* vis_obs_im, void* weights) except + void _sweep(int nr, void* intensity, dreal Rmin, dreal dR, int nxy, dreal dxy, dreal inc, void* image) except + void _uv_rotate(dreal PA, dreal dRA, dreal dDec, void* dRArot, void* dDecrot, int nd, void* u, void* v, void* urot, void* vrot) except + diff --git a/python/libcommon.pyx b/python/libcommon.pyx index 8fa6ae8..d9a0f71 100644 --- a/python/libcommon.pyx +++ b/python/libcommon.pyx @@ -717,6 +717,110 @@ def chi2Image(dreal[:,::1] image, dxy, dreal[::1] u, dreal[::1] v, return cpp._chi2_image(image.shape[0], image.shape[1], &image[0,0], v_origin, dRA, dDec, duv, PA, len(u), &u[0], &v[0], &vis_obs_re[0], &vis_obs_im[0], &vis_obs_w[0]) +def chi2UnstructuredImage(dreal[::1] x, dreal[::1] y, dreal[::1] image, + int nxy, dreal dxy, dreal[::1] u, dreal[::1] v, + dreal[::1] vis_obs_re, dreal[::1] vis_obs_im, dreal[::1] vis_obs_w, + dRA=0., dDec=0., PA=0., check=False, origin='upper'): + """ + Compute the chi square of a model unstructured image given the observed visibilities. + + The chi square is computed from the observed and synthetic visibilities as: + + .. math:: + + \chi^2 = \sum_{j=1}^N w_j * [(Re V_{obs\ j}-Re V_{mod\ j})^2 + (Im V_{obs\ j}-Im V_{mod\ j})^2] + + where :math:`V_{mod}` are the synthetic visibilities, which are computed internally + as in :func:`.sampleUnstructuredImage`. + + Typical call signature:: + + chi2 = chi2UnstructuredImage(x, y, image, nxy, dxy, u, v, vis_obs_re, vis_obs_im, vis_obs_w, + dRA=0, dDec=0, PA=0, check=False, origin='upper') + + Parameters + ---------- + x : 1D array_like, float + List of x coordinates at which intensities are known. + **units**: rad + y : 1D array_like, float + List of y coordinates at which intensities are known. + **units**: rad + image : 1D array_like, float + Array containing the surface brightness of the model. + Assume the x-axis (R.A.) increases from right (West) to left (East) + and the y-axis (Dec.) increases from bottom (South) to top (North). + `nxy` must be even. + **units**: Jy/st + nxy : int + Number of pixels to use for the interpolated gridded image. + dxy : float + Size of the image cell in the interpolated image, assumed equal in both x and y direction. + **units**: rad + u : array_like, float + u coordinate of the visibility points where the FT has to be sampled. + **units**: wavelength + v : array_like, float + v coordinate of the visibility points where the FT has to be sampled. + The length of `v` must be equal to the length of `u`. + **units**: wavelength + vis_obs_re : array_like, float + Real part of the observed visibilities. + **units**: Jy + vis_obs_im: array_like, float + Imaginary part of the observed visibilities. + The length of `vis_obs_im` must be equal to the length of `vis_obs_re`. + **units**: Jy + vis_obs_w: array_like, float + Weight associated to the observed visibilities. + The length of `vis_obs_w` must be equal to the length of `vis_obs_re`. + **units**: + dRA : float, optional + R.A. offset w.r.t. the phase center by which the image is translated. + If dRA > 0 translate the image towards the left (East). Default is 0. + **units**: rad + dDec : float, optional + Dec. offset w.r.t. the phase center by which the image is translated. + If dDec > 0 translate the image towards the top (North). Default is 0. + **units**: rad + PA : float, optional + Position Angle, defined East of North. Default is 0. + **units**: rad + check : bool, optional + If True, check whether `image` and `dxy` satisfy Nyquist criterion for + computing the synthetic visibilities in the (u, v) locations provided. + Additionally check that the (u, v) points fall in the image to avoid + segmentation violations. Default is False since the check might take + time. For executions where speed is important, set to False. + origin : ['upper' | 'lower'], optional + Set the [0,0] pixel index of the matrix in the upper left or lower left corner of the axes. + It follows the same convention as in matplotlib `matshow` and `imshow` commands. + Declination axis and the matrix y axis are parallel for `origin='lower'`, anti-parallel for `origin='upper'`. + The central pixel corresponding to the (RA, Dec) = (0, 0) is always [Nxy/2, Nxy/2]. + For more details see the Technical Requirements page in the online docs. + + Returns + ------- + chi2: float + The chi square, not normalized. + + See also + -------- + :func:`.sampleImage` + + """ + check_obs(vis_obs_re, vis_obs_im, vis_obs_w, u=u, v=v) + + duv = 1 / (dxy*nxy) + + if check: + check_image_size(u, v, nxy, dxy, duv) + + v_origin = set_v_origin(origin) + + return cpp._chi2_unstructured_image(&x[0], &y[0], nxy, nxy, dxy, len(x), &image[0], v_origin, dRA, dDec, duv, PA, len(u), &u[0], &v[0], &vis_obs_re[0], &vis_obs_im[0], &vis_obs_w[0]) + + def chi2Profile(dreal[::1] intensity, Rmin, dR, nxy, dxy, dreal[::1] u, dreal[::1] v, dreal[::1] vis_obs_re, dreal[::1] vis_obs_im, dreal[::1] vis_obs_w, dRA=0., dDec=0., PA=0., inc=0., check=False): diff --git a/src/galario.cpp b/src/galario.cpp index 4880cb3..e35cc62 100644 --- a/src/galario.cpp +++ b/src/galario.cpp @@ -1975,6 +1975,65 @@ dreal _chi2_image(int nx, int ny, void* realdata, const dreal v_origin, dreal dR static_cast(weights)); } +dreal chi2_unstructured_image(const dreal* realx, const dreal* realy, int nx, int ny, dreal dxy, int ni, const dreal* realdata, const dreal v_origin, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, const dreal* u, const dreal* v, const dreal* vis_obs_re, const dreal* vis_obs_im, const dreal* weights) { + CPUTimer t_start; + + CHECK_INPUTXY(nx, ny); + dreal chi2 = 0; +#ifdef __CUDACC__ + GPUTimer t; + // ################################ + // ### ALLOCATION, INITIALIZATION ### + // ################################ + + /* async memory copy: + TODO copy memory asynchronously or create streams to define dependencies + use nonzero cudaStream_t + kernel<<< blocks, threads, bytes=0, stream =! 0>>>(); + + all cufft calls are asynchronous, can specify the stream explicitly (cf. doc) + same for cublas + draw dependcies on paper: first thing is to do fft while other data is transferred + + While the FFT etc. are calculated, we can copy over the weights and observed values. + */ + // reserve memory for the interpolated values + //CudaMemory vis_int_d(nd); + //t.Elapsed("chi2_image::malloc_vis_int"); + + // Initialization for comparison and chi square computation + /* allocate and copy observational data */ + /*CudaMemory vis_obs_re_d(nd, vis_obs_re); + CudaMemory vis_obs_im_d(nd, vis_obs_im); + CudaMemory weights_d(nd, weights); + t.Elapsed("chi2_image::copy_observations"); + + auto data_d = copy_input_d(nx, ny, realdata); + + sample_d(nx, ny, data_d.ptr, v_origin, dRA, dDec, nd, duv, PA, u, v, vis_int_d.ptr); + chi2 = reduce_chi2_d(nd, vis_obs_re_d.ptr, vis_obs_im_d.ptr, vis_int_d.ptr, weights_d.ptr);*/ +#else + CPUTimer t; + + auto vis_int = reinterpret_cast(FFTW(alloc_complex)(nd)); t.Elapsed("chi2_imag::fftw_alloc"); + sample_unstructured_image(realx, realy, nx, ny, dxy, ni, realdata, v_origin, dRA, dDec, duv, PA, nd, u, v, vis_int); + + chi2 = reduce_chi2(nd, vis_obs_re, vis_obs_im, vis_int, weights); + + t = CPUTimer(); galario_free(vis_int); t.Elapsed("chi2_imag::free_vis_int"); +#endif + t_start.Elapsed("chi2_image_tot"); + flush_timing(); + + return chi2; +} + +dreal _chi2_unstructured_image(void* realx, void* realy, int nx, int ny, dreal dxy, int ni, void* realdata, const dreal v_origin, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, void* u, void* v, void* vis_obs_re, void* vis_obs_im, void* weights) { + return chi2_unstructured_image(static_cast(realx), static_cast(realy), nx, ny, dxy, ni, static_cast(realdata), v_origin, dRA, dDec, duv, PA, nd, + static_cast(u), static_cast(v), static_cast(vis_obs_re), static_cast(vis_obs_im), + static_cast(weights)); +} + dreal chi2_profile(int nr, dreal *const intensity, dreal Rmin, dreal dR, dreal dxy, int nxy, dreal inc, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, const dreal *u, const dreal *v, const dreal *vis_obs_re, const dreal *vis_obs_im, const dreal *weights) { diff --git a/src/galario.h b/src/galario.h index 68a00e5..1787edc 100644 --- a/src/galario.h +++ b/src/galario.h @@ -32,6 +32,7 @@ dreal chi2_profile(int nr, const dreal* intensity, dreal Rmin, dreal dR, dreal d dreal dDec, dreal duv, dreal PA, int nd, const dreal *u, const dreal *v, const dreal *vis_obs_re, const dreal *vis_obs_im, const dreal *weights); dreal chi2_image(int nx, int ny, const dreal* image, const dreal v_origin, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, const dreal* u, const dreal* v, const dreal* vis_obs_re, const dreal* vis_obs_im, const dreal* weights); +dreal chi2_unstructured_image(const dreal* realx, const dreal* realy, int nx, int ny, dreal dxy, int ni, const dreal* realdata, const dreal v_origin, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, const dreal* u, const dreal* v, const dreal* vis_obs_re, const dreal* vis_obs_im, const dreal* weights); void sweep(int nr, const dreal* intensity, dreal Rmin, dreal dR, int nxy, dreal dxy, dreal inc, dcomplex *image); void uv_rotate(dreal PA, dreal dRA, dreal dDec, dreal* dRArot, dreal* dDecrot, int nd, const dreal* u, const dreal* v, dreal* urot, dreal* vrot); diff --git a/src/galario_py.h b/src/galario_py.h index 39a360c..fca8cca 100644 --- a/src/galario_py.h +++ b/src/galario_py.h @@ -34,6 +34,7 @@ void _sample_unstructured_image(void* x, void* y, int nx, int ny, dreal dxy, int dreal _chi2_profile(int nr, void *intensity, dreal Rmin, dreal dR, dreal dxy, int nxy, dreal inc, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, void *u, void *v, void *vis_obs_re, void *vis_obs_im, void *weights); dreal _chi2_image(int nx, int ny, void* data, dreal v_origin, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, void* u, void* v, void* vis_obs_re, void* vis_obs_im, void* weights); +dreal _chi2_unstructured_image(void* realx, void* realy, int nx, int ny, dreal dxy, int ni, void* realdata, const dreal v_origin, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, void* u, void* v, void* vis_obs_re, void* vis_obs_im, void* weights); void _sweep(int nr, void *intensity, dreal Rmin, dreal dR, int nxy, dreal dxy, dreal inc, void *image); void _uv_rotate(dreal PA, dreal dRA, dreal dDec, void* dRArot, void* dDecrot, int nd, void* u, void* v, void* urot, void* vrot);