From f06761ffbc53343961acb75673aff9110b4db48b Mon Sep 17 00:00:00 2001 From: chrishavlin Date: Fri, 4 Oct 2024 11:27:21 -0500 Subject: [PATCH 1/7] release gil in pixelize_cylinder --- yt/utilities/lib/pixelization_routines.pyx | 118 +++++++++++---------- 1 file changed, 60 insertions(+), 58 deletions(-) diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx index a207473e9ea..2702a89005d 100644 --- a/yt/utilities/lib/pixelization_routines.pyx +++ b/yt/utilities/lib/pixelization_routines.pyx @@ -566,6 +566,7 @@ def pixelize_cylinder(np.float64_t[:,:] buff, cdef np.float64_t r_inc, theta_inc cdef np.float64_t costheta, sintheta cdef int i, i1, pi, pj + cdef np.float64_t npPI = np.pi cdef int imin, imax imin = np.asarray(radius).argmin() @@ -611,65 +612,66 @@ def pixelize_cylinder(np.float64_t[:,:] buff, rbounds[0] = 0.0 r_inc = 0.5 * fmin(dx, dy) - for i in range(radius.shape[0]): - r0 = radius[i] - theta0 = theta[i] - dr_i = dradius[i] - dtheta_i = dtheta[i] - # Skip out early if we're offsides, for zoomed in plots - if r0 + dr_i < rbounds[0] or r0 - dr_i > rbounds[1]: - continue - theta_i = theta0 - dtheta_i - theta_inc = r_inc / (r0 + dr_i) - - while theta_i < theta0 + dtheta_i: - r_i = r0 - dr_i - costheta = math.cos(theta_i) - sintheta = math.sin(theta_i) - while r_i < r0 + dr_i: - if rmax <= r_i: + with nogil: + for i in range(radius.shape[0]): + r0 = radius[i] + theta0 = theta[i] + dr_i = dradius[i] + dtheta_i = dtheta[i] + # Skip out early if we're offsides, for zoomed in plots + if r0 + dr_i < rbounds[0] or r0 - dr_i > rbounds[1]: + continue + theta_i = theta0 - dtheta_i + theta_inc = r_inc / (r0 + dr_i) + + while theta_i < theta0 + dtheta_i: + r_i = r0 - dr_i + costheta = math.cos(theta_i) + sintheta = math.sin(theta_i) + while r_i < r0 + dr_i: + if rmax <= r_i: + r_i += r_inc + continue + y = r_i * costheta + x = r_i * sintheta + pi = ((x - x0)/dx) + pj = ((y - y0)/dy) + if pi >= 0 and pi < buff.shape[0] and \ + pj >= 0 and pj < buff.shape[1]: + # we got a pixel that intersects the grid cell + # now check that this pixel doesn't go beyond the data domain + xp = x0 + pi*dx + yp = y0 + pj*dy + corners[0] = xp*xp + yp*yp + corners[1] = xp*xp + (yp+dy)**2 + corners[2] = (xp+dx)**2 + yp*yp + corners[3] = (xp+dx)**2 + (yp+dy)**2 + prbounds[0] = prbounds[1] = corners[3] + for i1 in range(3): + prbounds[0] = fmin(prbounds[0], corners[i1]) + prbounds[1] = fmax(prbounds[1], corners[i1]) + prbounds[0] = math.sqrt(prbounds[0]) + prbounds[1] = math.sqrt(prbounds[1]) + + corners[0] = math.atan2(xp, yp) + corners[1] = math.atan2(xp, yp+dy) + corners[2] = math.atan2(xp+dx, yp) + corners[3] = math.atan2(xp+dx, yp+dy) + ptbounds[0] = ptbounds[1] = corners[3] + for i1 in range(3): + ptbounds[0] = fmin(ptbounds[0], corners[i1]) + ptbounds[1] = fmax(ptbounds[1], corners[i1]) + + # shift to a [0, PI] interval + ptbounds[0] = ptbounds[0] % (2*npPI) + ptbounds[1] = ptbounds[1] % (2*npPI) + + if prbounds[0] >= rmin and prbounds[1] <= rmax and \ + ptbounds[0] >= tmin and ptbounds[1] <= tmax: + buff[pi, pj] = field[i] + mask[pi, pj] = 1 r_i += r_inc - continue - y = r_i * costheta - x = r_i * sintheta - pi = ((x - x0)/dx) - pj = ((y - y0)/dy) - if pi >= 0 and pi < buff.shape[0] and \ - pj >= 0 and pj < buff.shape[1]: - # we got a pixel that intersects the grid cell - # now check that this pixel doesn't go beyond the data domain - xp = x0 + pi*dx - yp = y0 + pj*dy - corners[0] = xp*xp + yp*yp - corners[1] = xp*xp + (yp+dy)**2 - corners[2] = (xp+dx)**2 + yp*yp - corners[3] = (xp+dx)**2 + (yp+dy)**2 - prbounds[0] = prbounds[1] = corners[3] - for i1 in range(3): - prbounds[0] = fmin(prbounds[0], corners[i1]) - prbounds[1] = fmax(prbounds[1], corners[i1]) - prbounds[0] = math.sqrt(prbounds[0]) - prbounds[1] = math.sqrt(prbounds[1]) - - corners[0] = math.atan2(xp, yp) - corners[1] = math.atan2(xp, yp+dy) - corners[2] = math.atan2(xp+dx, yp) - corners[3] = math.atan2(xp+dx, yp+dy) - ptbounds[0] = ptbounds[1] = corners[3] - for i1 in range(3): - ptbounds[0] = fmin(ptbounds[0], corners[i1]) - ptbounds[1] = fmax(ptbounds[1], corners[i1]) - - # shift to a [0, PI] interval - ptbounds[0] = ptbounds[0] % (2*np.pi) - ptbounds[1] = ptbounds[1] % (2*np.pi) - - if prbounds[0] >= rmin and prbounds[1] <= rmax and \ - ptbounds[0] >= tmin and ptbounds[1] <= tmax: - buff[pi, pj] = field[i] - mask[pi, pj] = 1 - r_i += r_inc - theta_i += theta_inc + theta_i += theta_inc if return_mask: return mask_arr.astype("bool") From 52698375a9175adf4f1c9f1fa3bcf44249ac789d Mon Sep 17 00:00:00 2001 From: Chris Havlin Date: Thu, 10 Oct 2024 10:16:07 -0500 Subject: [PATCH 2/7] pre-compute twoPI Co-authored-by: Corentin Cadiou --- yt/utilities/lib/pixelization_routines.pyx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx index 2702a89005d..2d5bf2d6fce 100644 --- a/yt/utilities/lib/pixelization_routines.pyx +++ b/yt/utilities/lib/pixelization_routines.pyx @@ -663,8 +663,8 @@ def pixelize_cylinder(np.float64_t[:,:] buff, ptbounds[1] = fmax(ptbounds[1], corners[i1]) # shift to a [0, PI] interval - ptbounds[0] = ptbounds[0] % (2*npPI) - ptbounds[1] = ptbounds[1] % (2*npPI) + ptbounds[0] = ptbounds[0] % twoPI + ptbounds[1] = ptbounds[1] % twoPI if prbounds[0] >= rmin and prbounds[1] <= rmax and \ ptbounds[0] >= tmin and ptbounds[1] <= tmax: From ad92f037f7c66f997078760b9919bb40ae402df7 Mon Sep 17 00:00:00 2001 From: Chris Havlin Date: Thu, 10 Oct 2024 10:17:18 -0500 Subject: [PATCH 3/7] precompute twoPI Co-authored-by: Corentin Cadiou --- yt/utilities/lib/pixelization_routines.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx index 2d5bf2d6fce..53ea75144f3 100644 --- a/yt/utilities/lib/pixelization_routines.pyx +++ b/yt/utilities/lib/pixelization_routines.pyx @@ -566,7 +566,7 @@ def pixelize_cylinder(np.float64_t[:,:] buff, cdef np.float64_t r_inc, theta_inc cdef np.float64_t costheta, sintheta cdef int i, i1, pi, pj - cdef np.float64_t npPI = np.pi + cdef np.float64_t twoPI = 2 * np.pi cdef int imin, imax imin = np.asarray(radius).argmin() From 3e1ef2c51387437c563f6365aa3a66bd99b2ee1d Mon Sep 17 00:00:00 2001 From: chrishavlin Date: Thu, 10 Oct 2024 14:55:14 -0500 Subject: [PATCH 4/7] handle theta<0 in shift to 0,2pi --- yt/utilities/lib/pixelization_routines.pyx | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx index 53ea75144f3..194ce7af6d4 100644 --- a/yt/utilities/lib/pixelization_routines.pyx +++ b/yt/utilities/lib/pixelization_routines.pyx @@ -81,6 +81,8 @@ cdef extern from "pixelization_constants.hpp": int WEDGE_NF np.uint8_t wedge_face_defs[MAX_NUM_FACES][2][2] +cdef extern from "numpy/npy_math.h": + double NPY_PI @cython.cdivision(True) @cython.boundscheck(False) @@ -566,7 +568,7 @@ def pixelize_cylinder(np.float64_t[:,:] buff, cdef np.float64_t r_inc, theta_inc cdef np.float64_t costheta, sintheta cdef int i, i1, pi, pj - cdef np.float64_t twoPI = 2 * np.pi + cdef np.float64_t twoPI = 2 * NPY_PI cdef int imin, imax imin = np.asarray(radius).argmin() @@ -578,7 +580,6 @@ def pixelize_cylinder(np.float64_t[:,:] buff, imax = np.asarray(theta).argmax() tmin = theta[imin] - dtheta[imin] tmax = theta[imax] + dtheta[imax] - cdef np.ndarray[np.uint8_t, ndim=2] mask_arr = np.zeros_like(buff, dtype="uint8") cdef np.uint8_t[:, :] mask = mask_arr @@ -664,7 +665,11 @@ def pixelize_cylinder(np.float64_t[:,:] buff, # shift to a [0, PI] interval ptbounds[0] = ptbounds[0] % twoPI + if ptbounds[0] < 0: + ptbounds[0] += NPY_PI ptbounds[1] = ptbounds[1] % twoPI + if ptbounds[1] < 0: + ptbounds[1] += NPY_PI if prbounds[0] >= rmin and prbounds[1] <= rmax and \ ptbounds[0] >= tmin and ptbounds[1] <= tmax: From 16c873759f6ebd864ae43a04713da214f7fa52cd Mon Sep 17 00:00:00 2001 From: chrishavlin Date: Mon, 21 Oct 2024 15:43:33 -0500 Subject: [PATCH 5/7] adjust the 0,2pi modulo calculation --- yt/utilities/lib/pixelization_routines.pyx | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx index 194ce7af6d4..da941beb986 100644 --- a/yt/utilities/lib/pixelization_routines.pyx +++ b/yt/utilities/lib/pixelization_routines.pyx @@ -663,13 +663,13 @@ def pixelize_cylinder(np.float64_t[:,:] buff, ptbounds[0] = fmin(ptbounds[0], corners[i1]) ptbounds[1] = fmax(ptbounds[1], corners[i1]) - # shift to a [0, PI] interval - ptbounds[0] = ptbounds[0] % twoPI - if ptbounds[0] < 0: - ptbounds[0] += NPY_PI - ptbounds[1] = ptbounds[1] % twoPI - if ptbounds[1] < 0: - ptbounds[1] += NPY_PI + # shift to a [0, 2*PI] interval + # note: in cython, % is an alias to fmod and the sign + # of the returned value matches the sign of the first + # argument, so need to offset by 2pi to ensure we + # always get a positive result in [0,2pi] + ptbounds[0] = math.fmod(ptbounds[0]+twoPI, twoPI) + ptbounds[1] = math.fmod(ptbounds[1]+twoPI, twoPI) if prbounds[0] >= rmin and prbounds[1] <= rmax and \ ptbounds[0] >= tmin and ptbounds[1] <= tmax: From 6c3fce579458173138e522f165712cffc589035e Mon Sep 17 00:00:00 2001 From: Chris Havlin Date: Tue, 22 Oct 2024 11:29:44 -0500 Subject: [PATCH 6/7] update note on fmod --- yt/utilities/lib/pixelization_routines.pyx | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx index da941beb986..67ac42ced8f 100644 --- a/yt/utilities/lib/pixelization_routines.pyx +++ b/yt/utilities/lib/pixelization_routines.pyx @@ -664,10 +664,9 @@ def pixelize_cylinder(np.float64_t[:,:] buff, ptbounds[1] = fmax(ptbounds[1], corners[i1]) # shift to a [0, 2*PI] interval - # note: in cython, % is an alias to fmod and the sign - # of the returned value matches the sign of the first - # argument, so need to offset by 2pi to ensure we - # always get a positive result in [0,2pi] + # note: with fmod, the sign of the returned value + # matches the sign of the first argument, so need + # to offset by 2pi to ensure a positive result in [0, 2pi] ptbounds[0] = math.fmod(ptbounds[0]+twoPI, twoPI) ptbounds[1] = math.fmod(ptbounds[1]+twoPI, twoPI) From 93b75d7be643151edec27f06997bb38c91e8e755 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 22 Oct 2024 16:47:20 +0000 Subject: [PATCH 7/7] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- yt/utilities/lib/pixelization_routines.pyx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx index 67ac42ced8f..a6a9f23d9c7 100644 --- a/yt/utilities/lib/pixelization_routines.pyx +++ b/yt/utilities/lib/pixelization_routines.pyx @@ -664,8 +664,8 @@ def pixelize_cylinder(np.float64_t[:,:] buff, ptbounds[1] = fmax(ptbounds[1], corners[i1]) # shift to a [0, 2*PI] interval - # note: with fmod, the sign of the returned value - # matches the sign of the first argument, so need + # note: with fmod, the sign of the returned value + # matches the sign of the first argument, so need # to offset by 2pi to ensure a positive result in [0, 2pi] ptbounds[0] = math.fmod(ptbounds[0]+twoPI, twoPI) ptbounds[1] = math.fmod(ptbounds[1]+twoPI, twoPI)