From 3ba2a4da8160f39a91d0e2303e5d8849c3efd18f Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Wed, 19 Feb 2025 17:26:23 +0100 Subject: [PATCH 01/68] worked on pixave sph projections --- yt/utilities/lib/pixelization_routines.pyx | 279 ++++++++++++++++++++- 1 file changed, 278 insertions(+), 1 deletion(-) diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx index a6a9f23d9c7..464d86aa9e7 100644 --- a/yt/utilities/lib/pixelization_routines.pyx +++ b/yt/utilities/lib/pixelization_routines.pyx @@ -1124,7 +1124,11 @@ cdef class SPHKernelInterpolationTable: @cython.boundscheck(False) @cython.wraparound(False) @cython.cdivision(True) -def pixelize_sph_kernel_projection( +# SPH projection: integrate the SPH kernel over a pencil beam +# through the center of each pixel in the output grid. +# The pencil beam is perpendicular to the grid plane, along the +# line-of-sight direction +def pixelize_sph_kernel_projection_pencilbeam( np.float64_t[:, :] buff, np.uint8_t[:, :] mask, any_float[:] posx, @@ -1322,6 +1326,279 @@ def pixelize_sph_kernel_projection( return mask + +# SPH projection: integrate the SPH kernel over the rectangular prism +# defined by each pixel in the output grid, and the depth of the +# projection. +# Note that we do assume here that along the projection direction, the +# SPH kernel is entirely contained in the projection volume. +# to avoid adding too much mass as a consequence, we cut off particle +# inclusion in the projection entirely if a particle center is outside +# the center +/- depth range +# implementation is based on Borrow & Kelly (2021): arXiv:2106.05281, +# https://ui.adsabs.harvard.edu/abs/2021arXiv210605281B/abstract +# except we use the table-interpolation-based line-of-sight kernel +# integration to use the full 3D kernel shape, +# and at least for the first go, omit many of the speed optimisations +# additionally, this paper technically descibes projecting masses +# (i.e., resolution-element integrated quantities), +# while this function is for densities +# default nsample_min is from swiftsimio +@cython.initializedcheck(False) +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) +def pixelize_sph_kernel_projection_pixelave( + np.float64_t[:, :] buff, + np.uint8_t[:, :] mask, + any_float[:] posx, + any_float[:] posy, + any_float[:] posz, + any_float[:] hsml, + any_float[:] pmass, + any_float[:] pdens, + any_float[:] quantity_to_smooth, + bounds, + kernel_name="cubic", + weight_field=None, + _check_period = (1, 1, 1), + period=None, + np.int_t nsample_min = 32): + + cdef np.intp_t xsize, ysize + cdef np.float64_t x_min, x_max, y_min, y_max, z_min, z_max, prefactor_j + cdef np.int64_t xi, yi, x0, x1, y0, y1, xxi, yyi, nx, ny, nsubx, nsuby + cdef np.float64_t q_ij2, posx_diff, posy_diff, ih_j2, insubx, insuby + cdef np.float64_t x, y, dx, dy, idx, idy, ipixA, h_j2, px, py, pz + cdef np.float64_t period_x = 0, period_y = 0, period_z = 0 + cdef int i, j, ii, jj, kk + cdef np.float64_t[:] _weight_field + cdef int * xiter + cdef int * yiter + cdef int * ziter + cdef np.float64_t * xiterv + cdef np.float64_t * yiterv + cdef np.float64_t * ziterv + cdef np.int8_t[3] check_period + + if weight_field is not None: + _weight_field = weight_field + + if period is not None: + period_x = period[0] + period_y = period[1] + period_z = period[2] + for i in range(3): + check_period[i] = np.int8(_check_period[i]) + # we find the x and y range over which we have pixels and we find how many + # pixels we have in each dimension + xsize, ysize = buff.shape[0], buff.shape[1] + x_min = bounds[0] + x_max = bounds[1] + y_min = bounds[2] + y_max = bounds[3] + z_min = bounds[4] + z_max = bounds[5] + + dx = (x_max - x_min) / xsize + dy = (y_max - y_min) / ysize + + idx = 1.0/dx + idy = 1.0/dy + ipixA = idx * idy + + if kernel_name not in kernel_tables: + kernel_tables[kernel_name] = SPHKernelInterpolationTable(kernel_name) + cdef SPHKernelInterpolationTable itab = kernel_tables[kernel_name] + with nogil, parallel(): + # loop through every particle + # NOTE: this loop can be quite time consuming. However it is easily + # parallelizable in multiple ways, such as: + # 1) use multiple cores to process individual particles (the outer loop) + # 2) use multiple cores to process individual pixels for a given particle + # (the inner loops) + # Depending on the ratio of particles' "sphere of influence" (a.k.a. the smoothing + # length) to the physical width of the pixels, different parallelization + # strategies may yield different speed-ups. Strategy #1 works better in the case + # of lots of itty bitty particles. Strategy #2 works well when we have a + # not-very-large-number of reasonably large-compared-to-pixels particles. We + # currently employ #1 as its workload is more even and consistent, even though it + # comes with a price of an additional, per thread memory for storing the + # intermediate results. + + local_buff = malloc(sizeof(np.float64_t) * xsize * ysize) + xiterv = malloc(sizeof(np.float64_t) * 2) + yiterv = malloc(sizeof(np.float64_t) * 2) + ziterv = malloc(sizeof(np.float64_t) * 2) + xiter = malloc(sizeof(int) * 2) + yiter = malloc(sizeof(int) * 2) + ziter = malloc(sizeof(int) * 2) + xiter[0] = yiter[0] = ziter[0] = 0 + xiterv[0] = yiterv[0] = ziterv[0] = 0.0 + for i in range(xsize * ysize): + local_buff[i] = 0.0 + + for j in prange(0, posx.shape[0], schedule="dynamic"): + if j % 100000 == 0: + with gil: + PyErr_CheckSignals() + + xiter[1] = yiter[1] = ziter[1] = 999 + + if check_period[0] == 1: + if posx[j] - hsml[j] < x_min: + xiter[1] = +1 + xiterv[1] = period_x + elif posx[j] + hsml[j] > x_max: + xiter[1] = -1 + xiterv[1] = -period_x + if check_period[1] == 1: + if posy[j] - hsml[j] < y_min: + yiter[1] = +1 + yiterv[1] = period_y + elif posy[j] + hsml[j] > y_max: + yiter[1] = -1 + yiterv[1] = -period_y + if check_period[2] == 1: + if posz[j] - hsml[j] < z_min: + ziter[1] = +1 + ziterv[1] = period_z + elif posz[j] + hsml[j] > z_max: + ziter[1] = -1 + ziterv[1] = -period_z + + # we set the smoothing length squared with lower limit of the pixel + # Nope! that causes weird grid resolution dependences and increases + # total values when resolution elements have hsml < grid spacing + h_j2 = hsml[j]*hsml[j] + ih_j2 = 1.0/h_j2 + + prefactor_j = pmass[j] / pdens[j] / hsml[j]**2 * quantity_to_smooth[j] + if weight_field is not None: + prefactor_j *= _weight_field[j] + + # Discussion point: do we want the hsml margin on the z direction? + # it's consistent with Ray and Region selections, I think, + # but does tend to 'tack on' stuff compared to the nominal depth + for kk in range(2): + # discard if z is outside bounds + if ziter[kk] == 999: continue + pz = posz[j] + ziterv[kk] + ## removed hsml 'margin' in the projection direction to avoid + ## double-counting particles near periodic edges + ## and adding extra 'depth' to projections + #if (pz + hsml[j] < z_min) or (pz - hsml[j] > z_max): continue + if (pz < z_min) or (pz > z_max): continue + + for ii in range(2): + if xiter[ii] == 999: continue + px = posx[j] + xiterv[ii] + if (px + hsml[j] < x_min) or (px - hsml[j] > x_max): + continue + for jj in range(2): + if yiter[jj] == 999: continue + py = posy[j] + yiterv[jj] + if (py + hsml[j] < y_min) or (py - hsml[j] > y_max): + continue + + # case: particle is small, + # overlaps with a small number of pixels (max. 4) + if hsml[j] < 0.5 * dx and hsml[j] < 0.5 * dy: + # lower right corner indices + x0 = ((px - hsml[j] - x_min)*idx) + y0 = ((py - hsml[j] - y_min)*idy) + # sph particle lies entirely in one pixel + if (px + hsml[j] < x_min + (x0 + 1) * dx and + py + hsml[j] < y_min + (y0 + 1) * dy): + # check that we're not on the wrong + # periodicity iteration + if (x0 > 0 and x0 < xsize and + y0 > 0 and y0 < ysize): + # surface density + # = density * av. length + # = density * volume / area + local_buff[x0 + y0 * xsize] += ( + pmass[j] / pdens[j] * ipixA * + quantity_to_smooth[j] + ) + mask[x0, y0] = 1 + else: + continue + + # sph particle lies in 2--4 pixels + # use a pre-calculated grid: TODO + else: pass + else: + # subsample + + # here we find the pixels which this + # particle contributes to + # subsampling does not depend on whether a + # particle overlaps a boundary + x0 = ((px - hsml[j] - x_min)*idx) + x1 = ((px + hsml[j] - x_min)*idx) + nx = x1 - x0 + x0 = iclip(x0-1, 0, xsize) + x1 = iclip(x1+1, 0, xsize) + + y0 = ((py - hsml[j] - y_min)*idy) + y1 = ((py + hsml[j] - y_min)*idy) + ny = y1 - y0 + y0 = iclip(y0-1, 0, ysize) + y1 = iclip(y1+1, 0, ysize) + + # using the same sampling rate in x and y directions + # should be fine for (typical) square grids, but + # less than optimal for large aspect ratios + nsubx = ( + ( nsample_min / + nx) + + 1) + insubx = 1. / nsubx + nsubx = ( + ( nsample_min / + ny) + + 1) + insuby = 1. / nsuby + + # found pixels we deposit on, + # loop through those pixels + for xi in range(x0, x1): + for yi in range(y0, y1): + for xsubi in range(nsubx): + x = x_min + \ + (xi + (xsubi + 0.5) * insubx) * dx + posx_diff = px - x + posx_diff = posx_diff * posx_diff + if posx_diff > h_j2: continue + for ysubi in range(nsuby): + y = y_min + \ + (yi + (ysubi + 0.5) * insuby)\ + * dy + posy_diff = py - y + posy_diff = posy_diff * posy_diff + if posy_diff > h_j2: continue + + q_ij2 = (posx_diff + posy_diff) \ + * ih_j2 + if q_ij2 >= 1: continue + local_buff[xi + yi*xsize] += \ + prefactor_j * insuby * insubx \ + * itab.interpolate(q_ij2) + mask[xi, yi] = 1 + + with gil: + for xxi in range(xsize): + for yyi in range(ysize): + buff[xxi, yyi] += local_buff[xxi + yyi*xsize] + free(local_buff) + free(xiterv) + free(yiterv) + free(xiter) + free(yiter) + + return mask + @cython.boundscheck(False) @cython.wraparound(False) def interpolate_sph_positions_gather(np.float64_t[:] buff, From 81100cb970c293130faefd992b0e03d86ba1bc63 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Wed, 19 Feb 2025 17:27:20 +0100 Subject: [PATCH 02/68] free all allocated memory at end of pencil beam SPH projection --- yt/utilities/lib/pixelization_routines.pyx | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx index 464d86aa9e7..c6cceb0b4b5 100644 --- a/yt/utilities/lib/pixelization_routines.pyx +++ b/yt/utilities/lib/pixelization_routines.pyx @@ -1321,8 +1321,10 @@ def pixelize_sph_kernel_projection_pencilbeam( free(local_buff) free(xiterv) free(yiterv) + free(ziterv) free(xiter) free(yiter) + free(ziter) return mask @@ -1594,8 +1596,10 @@ def pixelize_sph_kernel_projection_pixelave( free(local_buff) free(xiterv) free(yiterv) + free(ziterv) free(xiter) free(yiter) + free(ziter) return mask From 7886d87822ea2592b439240749e67b4fcc79a538 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Wed, 19 Feb 2025 18:59:16 +0100 Subject: [PATCH 03/68] first go at on-axis pixel-ave SPH projection backend --- yt/utilities/lib/pixelization_routines.pyx | 75 ++++++++++++++++++---- 1 file changed, 64 insertions(+), 11 deletions(-) diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx index c6cceb0b4b5..dc4d8b0cecc 100644 --- a/yt/utilities/lib/pixelization_routines.pyx +++ b/yt/utilities/lib/pixelization_routines.pyx @@ -1370,8 +1370,9 @@ def pixelize_sph_kernel_projection_pixelave( cdef np.intp_t xsize, ysize cdef np.float64_t x_min, x_max, y_min, y_max, z_min, z_max, prefactor_j cdef np.int64_t xi, yi, x0, x1, y0, y1, xxi, yyi, nx, ny, nsubx, nsuby - cdef np.float64_t q_ij2, posx_diff, posy_diff, ih_j2, insubx, insuby - cdef np.float64_t x, y, dx, dy, idx, idy, ipixA, h_j2, px, py, pz + cdef np.float64_t q_ij2, posx_diff, posy_diff, ih_j2 + cdef np.float_64_t insubx, insuby, insample_min, xnorm, ynorm + cdef np.float64_t x, y, dx, dy, idx, idy, ipixA, subA, h_j2, px, py, pz cdef np.float64_t period_x = 0, period_y = 0, period_z = 0 cdef int i, j, ii, jj, kk cdef np.float64_t[:] _weight_field @@ -1412,6 +1413,22 @@ def pixelize_sph_kernel_projection_pixelave( if kernel_name not in kernel_tables: kernel_tables[kernel_name] = SPHKernelInterpolationTable(kernel_name) cdef SPHKernelInterpolationTable itab = kernel_tables[kernel_name] + + # pre-calculate kernel 1D integral values to use for small particles + kernel1dint = malloc( + sizeof(np.float64_t) * nsample_min * nsample_min) + # get line-of-sight kernel integral values + insample_min = 1.0 / nsample_min + for i in range(nsample_min): + xnorm = -1.0 + 2.0 * (i + 0.5) * insample_min + xnorm = xnorm * xnorm + for j in range(nsample_min): + ynorm = -1.0 + 2.0 * (i + 0.5) * insample_min + q_ij2 = xnorm + ynorm * ynorm + if q_ij2 >= 1.: + kernel1dint[nsample_min * i + j] = 0.0 + else: + kernel1dint[nsample_min * i + j] = itab.interpolate(q_ij2) with nogil, parallel(): # loop through every particle # NOTE: this loop can be quite time consuming. However it is easily @@ -1427,6 +1444,7 @@ def pixelize_sph_kernel_projection_pixelave( # currently employ #1 as its workload is more even and consistent, even though it # comes with a price of an additional, per thread memory for storing the # intermediate results. + # !! These comments were written for local_buff = malloc(sizeof(np.float64_t) * xsize * ysize) xiterv = malloc(sizeof(np.float64_t) * 2) @@ -1506,7 +1524,7 @@ def pixelize_sph_kernel_projection_pixelave( # case: particle is small, # overlaps with a small number of pixels (max. 4) if hsml[j] < 0.5 * dx and hsml[j] < 0.5 * dy: - # lower right corner indices + # lower left corner indices x0 = ((px - hsml[j] - x_min)*idx) y0 = ((py - hsml[j] - y_min)*idy) # sph particle lies entirely in one pixel @@ -1529,7 +1547,41 @@ def pixelize_sph_kernel_projection_pixelave( # sph particle lies in 2--4 pixels # use a pre-calculated grid: TODO - else: pass + else: # use pre-calculated kernel values + # pre-calculated grid has size + # (nsample_min, nsample_min) + # bin edge values -1 -- 1 + # in impact parameter / hsml units + + # where does the overlapped pixel edge + # fall in the save kernel space array? + xnorm = (x_min + (x0 + 1) * dx - px) \ + / hsml[j] + ynorm = (y_min + (y0 + 1) * dy - py) \ + / hsml[j] + # round to the nearest kernel grid edge + # (C float to int cast behaves like floor()) + xi = ( + 0.5 * (xnorm + 1.) * nsample_min + 0.5) + yi = ( + 0.5 * (ynorm + 1.) * nsample_min + 0.5) + # weight for each line integral: + # area of subsampled array in length units + # divided by output grid pixel size + prefactor_j *= 4.0 * hsml[j] * hsml[j] \ + * insample_min * insample_min \ + * ipixA + + for xxi in range(0, xi): + for yyi in range(0, yi): + if kernel1dint[xxi, yyi] == 0.: + # no erroneous mask = 1 + continue + local_buff[x0 + y0*xsize] += \ + prefactor_j \ + * kernel1dint[xxi, yyi] + mask[xi, yi] = 1 + else: # subsample @@ -1540,14 +1592,14 @@ def pixelize_sph_kernel_projection_pixelave( x0 = ((px - hsml[j] - x_min)*idx) x1 = ((px + hsml[j] - x_min)*idx) nx = x1 - x0 - x0 = iclip(x0-1, 0, xsize) - x1 = iclip(x1+1, 0, xsize) + x0 = iclip(x0 - 1, 0, xsize) + x1 = iclip(x1 + 1, 0, xsize) y0 = ((py - hsml[j] - y_min)*idy) y1 = ((py + hsml[j] - y_min)*idy) ny = y1 - y0 - y0 = iclip(y0-1, 0, ysize) - y1 = iclip(y1+1, 0, ysize) + y0 = iclip(y0 - 1, 0, ysize) + y1 = iclip(y1 + 1, 0, ysize) # using the same sampling rate in x and y directions # should be fine for (typical) square grids, but @@ -1556,12 +1608,12 @@ def pixelize_sph_kernel_projection_pixelave( ( nsample_min / nx) + 1) - insubx = 1. / nsubx + insubx = 1.0 / nsubx nsubx = ( ( nsample_min / ny) + 1) - insuby = 1. / nsuby + insuby = 1.0 / nsuby # found pixels we deposit on, # loop through those pixels @@ -1583,7 +1635,7 @@ def pixelize_sph_kernel_projection_pixelave( q_ij2 = (posx_diff + posy_diff) \ * ih_j2 - if q_ij2 >= 1: continue + if q_ij2 >= 1.0: continue local_buff[xi + yi*xsize] += \ prefactor_j * insuby * insubx \ * itab.interpolate(q_ij2) @@ -1594,6 +1646,7 @@ def pixelize_sph_kernel_projection_pixelave( for yyi in range(ysize): buff[xxi, yyi] += local_buff[xxi + yyi*xsize] free(local_buff) + free(kernel1dint) free(xiterv) free(yiterv) free(ziterv) From 3d461e2da328a50c123aa695a65396ca88a20c6c Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Thu, 20 Feb 2025 12:51:20 +0100 Subject: [PATCH 04/68] edit SPH pix-av proj. backend: 2x2 kernel case --- yt/utilities/lib/pixelization_routines.pyx | 164 +++++++++++++-------- 1 file changed, 106 insertions(+), 58 deletions(-) diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx index dc4d8b0cecc..d7b9559f934 100644 --- a/yt/utilities/lib/pixelization_routines.pyx +++ b/yt/utilities/lib/pixelization_routines.pyx @@ -1370,11 +1370,11 @@ def pixelize_sph_kernel_projection_pixelave( cdef np.intp_t xsize, ysize cdef np.float64_t x_min, x_max, y_min, y_max, z_min, z_max, prefactor_j cdef np.int64_t xi, yi, x0, x1, y0, y1, xxi, yyi, nx, ny, nsubx, nsuby - cdef np.float64_t q_ij2, posx_diff, posy_diff, ih_j2 - cdef np.float_64_t insubx, insuby, insample_min, xnorm, ynorm - cdef np.float64_t x, y, dx, dy, idx, idy, ipixA, subA, h_j2, px, py, pz + cdef np.float64_t q_ij2, posx_diff, posy_diff, ih_j2, qts_j + cdef np.float_64_t insubx, insuby, insmin, xnorm, ynorm, kv + cdef np.float64_t x, y, dx, dy, idx, idy, ipixA, h_j2, px, py, pz cdef np.float64_t period_x = 0, period_y = 0, period_z = 0 - cdef int i, j, ii, jj, kk + cdef int i, j, ii, jj, kk, ic, nsmin cdef np.float64_t[:] _weight_field cdef int * xiter cdef int * yiter @@ -1409,26 +1409,61 @@ def pixelize_sph_kernel_projection_pixelave( idx = 1.0/dx idy = 1.0/dy ipixA = idx * idy + nsmin = nsample_min if kernel_name not in kernel_tables: kernel_tables[kernel_name] = SPHKernelInterpolationTable(kernel_name) cdef SPHKernelInterpolationTable itab = kernel_tables[kernel_name] - # pre-calculate kernel 1D integral values to use for small particles - kernel1dint = malloc( - sizeof(np.float64_t) * nsample_min * nsample_min) - # get line-of-sight kernel integral values - insample_min = 1.0 / nsample_min - for i in range(nsample_min): - xnorm = -1.0 + 2.0 * (i + 0.5) * insample_min + # pre-calculate kernel integral values to use for small particles + # overlapping max. 2x2 output grid pixels. dimensions: + # (part. position x, part. position y, kernel x, kernel y) + # kernel x, kernel y: shape (2, 2), values for the output grid + # part. x, y positions: position of the center pixel boundary of + # the 2x2 output values in x and y directions. + # units: (boundary x/y - particle center x/y) / particle hsml + # value range: [-1.0, 1.0] + # (The volume integral of this kernel should be one. We could + # normalize it explicitly, but skip that step for now, since we + # also don't do this for larger particles.) + kern2by2 = malloc( + sizeof(np.float64_t) * (nsmin + 1) * (nsmin + 1) * 4) + # step 1: just get a big grid of the 1D integral values + kern1dint = malloc( + sizeof(np.float64_t) * nsmin * nsmin) + insmin = 1.0 / nsmin + for i in range(nsmin): + xnorm = -1.0 + 2.0 * (i + 0.5) * insmin xnorm = xnorm * xnorm - for j in range(nsample_min): - ynorm = -1.0 + 2.0 * (i + 0.5) * insample_min + for j in range(nsmin): + ynorm = -1.0 + 2.0 * (i + 0.5) * insmin q_ij2 = xnorm + ynorm * ynorm if q_ij2 >= 1.: - kernel1dint[nsample_min * i + j] = 0.0 + kern1dint[nsample_min * i + j] = 0.0 else: - kernel1dint[nsample_min * i + j] = itab.interpolate(q_ij2) + kern1dint[nsample_min * i + j] = itab.interpolate(q_ij2) + # step 2: integrate the 1D values for each grid edge position + for i in range(nsample_min + 1): + for j in range(nsample_min + 1): + # slowest 2 indices + ic = (4 * (nsample_min + 1) * (nsample_min + 1) * i + + 4 * (nsample_min + 1) * j) + kern2by2[ic] = 0.0 + kern2by2[ic + 1] = 0.0 + kern2by2[ic + 2] = 0.0 + kern2by2[ic + 3] = 0.0 + for ii in range(0, i): + for jj in range(0, j): + kern2by2[ic] += kern1dint[nsample_min * ii + jj] + for jj in range(j, nsample_min): + kern2by2[ic + 1] += kern1dint[nsample_min * ii + jj] + for ii in range(i, nsample_min): + for jj in range(0, j): + kern2by2[ic + 2] += kern1dint[nsample_min * ii + jj] + for jj in range(j, nsample_min): + kern2by2[ic + 3] += kern1dint[nsample_min * ii + jj] + free(kern1dint) + with nogil, parallel(): # loop through every particle # NOTE: this loop can be quite time consuming. However it is easily @@ -1444,7 +1479,8 @@ def pixelize_sph_kernel_projection_pixelave( # currently employ #1 as its workload is more even and consistent, even though it # comes with a price of an additional, per thread memory for storing the # intermediate results. - # !! These comments were written for + # !!! These comments were written for the line of sight integral + # !!! counterpart to this function, and not checked for this version local_buff = malloc(sizeof(np.float64_t) * xsize * ysize) xiterv = malloc(sizeof(np.float64_t) * 2) @@ -1490,13 +1526,16 @@ def pixelize_sph_kernel_projection_pixelave( # we set the smoothing length squared with lower limit of the pixel # Nope! that causes weird grid resolution dependences and increases # total values when resolution elements have hsml < grid spacing - h_j2 = hsml[j]*hsml[j] - ih_j2 = 1.0/h_j2 - - prefactor_j = pmass[j] / pdens[j] / hsml[j]**2 * quantity_to_smooth[j] + h_j2 = hsml[j] * hsml[j] + ih_j2 = 1.0 / h_j2 + + # for small particles, we use the particle to be smoothed + # directly, for larger ones, we need the whole prefactor + qts_j = quantity_to_smooth[j] if weight_field is not None: - prefactor_j *= _weight_field[j] - + qts_j *= _weight_field[j] + prefactor_j = pmass[j] / pdens[j] * ih_j2 * qts_j + # Discussion point: do we want the hsml margin on the z direction? # it's consistent with Ray and Region selections, I think, # but does tend to 'tack on' stuff compared to the nominal depth @@ -1522,11 +1561,11 @@ def pixelize_sph_kernel_projection_pixelave( continue # case: particle is small, - # overlaps with a small number of pixels (max. 4) + # overlaps with a small number of pixels (max. 2x2) if hsml[j] < 0.5 * dx and hsml[j] < 0.5 * dy: # lower left corner indices - x0 = ((px - hsml[j] - x_min)*idx) - y0 = ((py - hsml[j] - y_min)*idy) + x0 = ((px - hsml[j] - x_min) * idx) + y0 = ((py - hsml[j] - y_min) * idy) # sph particle lies entirely in one pixel if (px + hsml[j] < x_min + (x0 + 1) * dx and py + hsml[j] < y_min + (y0 + 1) * dy): @@ -1538,23 +1577,24 @@ def pixelize_sph_kernel_projection_pixelave( # = density * av. length # = density * volume / area local_buff[x0 + y0 * xsize] += ( - pmass[j] / pdens[j] * ipixA * - quantity_to_smooth[j] + pmass[j] / pdens[j] * ipixA * qts_j ) mask[x0, y0] = 1 else: continue # sph particle lies in 2--4 pixels - # use a pre-calculated grid: TODO + # use a pre-calculated grid else: # use pre-calculated kernel values # pre-calculated grid has size - # (nsample_min, nsample_min) + # (nsmin + 1, nsmin + 1, 2, 2) # bin edge values -1 -- 1 - # in impact parameter / hsml units + # in units: + # center of 2x2 output pix grid + # offset from particle center / hsml # where does the overlapped pixel edge - # fall in the save kernel space array? + # fall in the saved kernel space array? xnorm = (x_min + (x0 + 1) * dx - px) \ / hsml[j] ynorm = (y_min + (y0 + 1) * dy - py) \ @@ -1562,48 +1602,54 @@ def pixelize_sph_kernel_projection_pixelave( # round to the nearest kernel grid edge # (C float to int cast behaves like floor()) xi = ( - 0.5 * (xnorm + 1.) * nsample_min + 0.5) + 0.5 * (xnorm + 1.) * nsmin + 0.5) yi = ( - 0.5 * (ynorm + 1.) * nsample_min + 0.5) + 0.5 * (ynorm + 1.) * nsmin + 0.5) # weight for each line integral: # area of subsampled array in length units # divided by output grid pixel size - prefactor_j *= 4.0 * hsml[j] * hsml[j] \ - * insample_min * insample_min \ - * ipixA - - for xxi in range(0, xi): - for yyi in range(0, yi): - if kernel1dint[xxi, yyi] == 0.: + prefactor_j *= (4.0 * h_j2 + * insmin * insmin * ipixA) + # coarse (part. pos rel. to grid) index + ci = (4 * nsmin * nsmin * xi + + 4 * nsmin * yi) + for xxi in range(2): + # catch on the next periodic loop + if x0 + xxi < 0 or x0 + xxi > xsize: + continue + for yyi in range(2): + # catch on the next periodic loop + if y0 + yyi < 0 or y0 + yyi > ysize: + continue + kv = kern2by2[ci + 2 * xxi + yyi] + if kv == 0.: # no erroneous mask = 1 continue - local_buff[x0 + y0*xsize] += \ - prefactor_j \ - * kernel1dint[xxi, yyi] + local_buff[x0 + xxi + + (y0 + yyi) * xsize] += ( + prefactor_j * kv + ) mask[xi, yi] = 1 - + else: # subsample # here we find the pixels which this # particle contributes to # subsampling does not depend on whether a - # particle overlaps a boundary + # particle overlaps a (periodic) edge x0 = ((px - hsml[j] - x_min)*idx) x1 = ((px + hsml[j] - x_min)*idx) - nx = x1 - x0 + nx = x1 - x0 + 2 x0 = iclip(x0 - 1, 0, xsize) x1 = iclip(x1 + 1, 0, xsize) y0 = ((py - hsml[j] - y_min)*idy) y1 = ((py + hsml[j] - y_min)*idy) - ny = y1 - y0 + ny = y1 - y0 + 2 y0 = iclip(y0 - 1, 0, ysize) y1 = iclip(y1 + 1, 0, ysize) - # using the same sampling rate in x and y directions - # should be fine for (typical) square grids, but - # less than optimal for large aspect ratios nsubx = ( ( nsample_min / nx) @@ -1626,19 +1672,21 @@ def pixelize_sph_kernel_projection_pixelave( posx_diff = posx_diff * posx_diff if posx_diff > h_j2: continue for ysubi in range(nsuby): - y = y_min + \ - (yi + (ysubi + 0.5) * insuby)\ - * dy + y = (y_min + + (yi + (ysubi + 0.5) * insuby) + * dy + ) posy_diff = py - y posy_diff = posy_diff * posy_diff if posy_diff > h_j2: continue - q_ij2 = (posx_diff + posy_diff) \ - * ih_j2 + q_ij2 = ((posx_diff + posy_diff) + * ih_j2) if q_ij2 >= 1.0: continue - local_buff[xi + yi*xsize] += \ - prefactor_j * insuby * insubx \ + local_buff[xi + yi*xsize] += ( + prefactor_j * insuby * insubx * itab.interpolate(q_ij2) + ) mask[xi, yi] = 1 with gil: @@ -1646,7 +1694,7 @@ def pixelize_sph_kernel_projection_pixelave( for yyi in range(ysize): buff[xxi, yyi] += local_buff[xxi + yyi*xsize] free(local_buff) - free(kernel1dint) + free(kern2by2) free(xiterv) free(yiterv) free(ziterv) From f09b2236973664f9899a71e68fbd63c1ae7ab60a Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 20 Feb 2025 11:54:53 +0000 Subject: [PATCH 05/68] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- yt/utilities/lib/pixelization_routines.pyx | 48 +++++++++++----------- 1 file changed, 24 insertions(+), 24 deletions(-) diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx index d7b9559f934..31395fd32f4 100644 --- a/yt/utilities/lib/pixelization_routines.pyx +++ b/yt/utilities/lib/pixelization_routines.pyx @@ -1330,7 +1330,7 @@ def pixelize_sph_kernel_projection_pencilbeam( # SPH projection: integrate the SPH kernel over the rectangular prism -# defined by each pixel in the output grid, and the depth of the +# defined by each pixel in the output grid, and the depth of the # projection. # Note that we do assume here that along the projection direction, the # SPH kernel is entirely contained in the projection volume. @@ -1344,7 +1344,7 @@ def pixelize_sph_kernel_projection_pencilbeam( # and at least for the first go, omit many of the speed optimisations # additionally, this paper technically descibes projecting masses # (i.e., resolution-element integrated quantities), -# while this function is for densities +# while this function is for densities # default nsample_min is from swiftsimio @cython.initializedcheck(False) @cython.boundscheck(False) @@ -1416,7 +1416,7 @@ def pixelize_sph_kernel_projection_pixelave( cdef SPHKernelInterpolationTable itab = kernel_tables[kernel_name] # pre-calculate kernel integral values to use for small particles - # overlapping max. 2x2 output grid pixels. dimensions: + # overlapping max. 2x2 output grid pixels. dimensions: # (part. position x, part. position y, kernel x, kernel y) # kernel x, kernel y: shape (2, 2), values for the output grid # part. x, y positions: position of the center pixel boundary of @@ -1424,7 +1424,7 @@ def pixelize_sph_kernel_projection_pixelave( # units: (boundary x/y - particle center x/y) / particle hsml # value range: [-1.0, 1.0] # (The volume integral of this kernel should be one. We could - # normalize it explicitly, but skip that step for now, since we + # normalize it explicitly, but skip that step for now, since we # also don't do this for larger particles.) kern2by2 = malloc( sizeof(np.float64_t) * (nsmin + 1) * (nsmin + 1) * 4) @@ -1480,7 +1480,7 @@ def pixelize_sph_kernel_projection_pixelave( # comes with a price of an additional, per thread memory for storing the # intermediate results. # !!! These comments were written for the line of sight integral - # !!! counterpart to this function, and not checked for this version + # !!! counterpart to this function, and not checked for this version local_buff = malloc(sizeof(np.float64_t) * xsize * ysize) xiterv = malloc(sizeof(np.float64_t) * 2) @@ -1528,14 +1528,14 @@ def pixelize_sph_kernel_projection_pixelave( # total values when resolution elements have hsml < grid spacing h_j2 = hsml[j] * hsml[j] ih_j2 = 1.0 / h_j2 - - # for small particles, we use the particle to be smoothed + + # for small particles, we use the particle to be smoothed # directly, for larger ones, we need the whole prefactor qts_j = quantity_to_smooth[j] if weight_field is not None: qts_j *= _weight_field[j] prefactor_j = pmass[j] / pdens[j] * ih_j2 * qts_j - + # Discussion point: do we want the hsml margin on the z direction? # it's consistent with Ray and Region selections, I think, # but does tend to 'tack on' stuff compared to the nominal depth @@ -1560,7 +1560,7 @@ def pixelize_sph_kernel_projection_pixelave( if (py + hsml[j] < y_min) or (py - hsml[j] > y_max): continue - # case: particle is small, + # case: particle is small, # overlaps with a small number of pixels (max. 2x2) if hsml[j] < 0.5 * dx and hsml[j] < 0.5 * dy: # lower left corner indices @@ -1569,7 +1569,7 @@ def pixelize_sph_kernel_projection_pixelave( # sph particle lies entirely in one pixel if (px + hsml[j] < x_min + (x0 + 1) * dx and py + hsml[j] < y_min + (y0 + 1) * dy): - # check that we're not on the wrong + # check that we're not on the wrong # periodicity iteration if (x0 > 0 and x0 < xsize and y0 > 0 and y0 < ysize): @@ -1579,18 +1579,18 @@ def pixelize_sph_kernel_projection_pixelave( local_buff[x0 + y0 * xsize] += ( pmass[j] / pdens[j] * ipixA * qts_j ) - mask[x0, y0] = 1 + mask[x0, y0] = 1 else: - continue + continue # sph particle lies in 2--4 pixels # use a pre-calculated grid else: # use pre-calculated kernel values - # pre-calculated grid has size + # pre-calculated grid has size # (nsmin + 1, nsmin + 1, 2, 2) # bin edge values -1 -- 1 # in units: - # center of 2x2 output pix grid + # center of 2x2 output pix grid # offset from particle center / hsml # where does the overlapped pixel edge @@ -1605,27 +1605,27 @@ def pixelize_sph_kernel_projection_pixelave( 0.5 * (xnorm + 1.) * nsmin + 0.5) yi = ( 0.5 * (ynorm + 1.) * nsmin + 0.5) - # weight for each line integral: + # weight for each line integral: # area of subsampled array in length units # divided by output grid pixel size prefactor_j *= (4.0 * h_j2 * insmin * insmin * ipixA) - # coarse (part. pos rel. to grid) index + # coarse (part. pos rel. to grid) index ci = (4 * nsmin * nsmin * xi - + 4 * nsmin * yi) + + 4 * nsmin * yi) for xxi in range(2): # catch on the next periodic loop if x0 + xxi < 0 or x0 + xxi > xsize: - continue + continue for yyi in range(2): # catch on the next periodic loop if y0 + yyi < 0 or y0 + yyi > ysize: - continue + continue kv = kern2by2[ci + 2 * xxi + yyi] - if kv == 0.: + if kv == 0.: # no erroneous mask = 1 continue - local_buff[x0 + xxi + local_buff[x0 + xxi + (y0 + yyi) * xsize] += ( prefactor_j * kv ) @@ -1634,7 +1634,7 @@ def pixelize_sph_kernel_projection_pixelave( else: # subsample - # here we find the pixels which this + # here we find the pixels which this # particle contributes to # subsampling does not depend on whether a # particle overlaps a (periodic) edge @@ -1649,7 +1649,7 @@ def pixelize_sph_kernel_projection_pixelave( ny = y1 - y0 + 2 y0 = iclip(y0 - 1, 0, ysize) y1 = iclip(y1 + 1, 0, ysize) - + nsubx = ( ( nsample_min / nx) @@ -1679,7 +1679,7 @@ def pixelize_sph_kernel_projection_pixelave( posy_diff = py - y posy_diff = posy_diff * posy_diff if posy_diff > h_j2: continue - + q_ij2 = ((posx_diff + posy_diff) * ih_j2) if q_ij2 >= 1.0: continue From 190c2ff99bd92e2db46987e65629f1d3195c53b9 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Fri, 21 Feb 2025 16:04:13 +0100 Subject: [PATCH 06/68] first go at connecting the backend to the rest of the code --- .../data_selection_objects.py | 16 +++++- .../coordinates/cartesian_coordinates.py | 40 +++++++++++--- yt/utilities/answer_testing/framework.py | 4 +- yt/utilities/lib/pixelization_routines.pyx | 42 +++++++++------ yt/visualization/fits_image.py | 53 +++++++++++++++++-- yt/visualization/fixed_resolution.py | 15 ++++++ .../volume_rendering/off_axis_projection.py | 17 ++++++ 7 files changed, 159 insertions(+), 28 deletions(-) diff --git a/yt/data_objects/selection_objects/data_selection_objects.py b/yt/data_objects/selection_objects/data_selection_objects.py index fc47b693dc0..bcc95306c7a 100644 --- a/yt/data_objects/selection_objects/data_selection_objects.py +++ b/yt/data_objects/selection_objects/data_selection_objects.py @@ -567,7 +567,9 @@ def _get_pw(self, fields, center, width, origin, plot_type): pw._setup_plots() return pw - def to_frb(self, width, resolution, center=None, height=None, periodic=False): + def to_frb(self, width, resolution, center=None, height=None, + periodic=False, + pixelmeaning: {"pixelave", "pencilbeam"} = "pixelave"): r"""This function returns a FixedResolutionBuffer generated from this object. @@ -595,7 +597,15 @@ def to_frb(self, width, resolution, center=None, height=None, periodic=False): periodic : bool Should the returned Fixed Resolution Buffer be periodic? (default: False). + pixelmeaning: + "pixelav": a pixel represents an average surface density or + surface-density-weighted average across a pixel. + "pencilbeam": a pixel represents a column density or + column-density-weighted average integrated over a pencil + beam through the pixel center. + + Only applies to projections of SPH datasets. Returns ------- frb : :class:`~yt.visualization.fixed_resolution.FixedResolutionBuffer` @@ -659,7 +669,9 @@ def to_frb(self, width, resolution, center=None, height=None, periodic=False): center[yax] - height * 0.5, center[yax] + height * 0.5, ) - frb = FixedResolutionBuffer(self, bounds, resolution, periodic=periodic) + frb = FixedResolutionBuffer(self, bounds, resolution, + periodic=periodic, + pixelmeaning=pixelmeaning) return frb diff --git a/yt/geometry/coordinates/cartesian_coordinates.py b/yt/geometry/coordinates/cartesian_coordinates.py index 4fed1977c4e..f6bcb431bf5 100644 --- a/yt/geometry/coordinates/cartesian_coordinates.py +++ b/yt/geometry/coordinates/cartesian_coordinates.py @@ -13,7 +13,8 @@ pixelize_element_mesh_line, pixelize_off_axis_cartesian, pixelize_sph_kernel_cutting, - pixelize_sph_kernel_projection, + pixelize_sph_kernel_projection_pencilbeam, + pixelize_sph_kernel_projection_pixelave, pixelize_sph_kernel_slice, ) from yt.utilities.math_utils import compute_stddev_image @@ -170,6 +171,7 @@ def pixelize( size, antialias=True, periodic=True, + pixelmeaning: {"pixelave", "pencilbeam"} = "pixelave", *, return_mask=False, ): @@ -177,6 +179,16 @@ def pixelize( Method for pixelizing datasets in preparation for two-dimensional image plots. Relies on several sampling routines written in cython + + pixelmeaning: + "pixelav": a pixel represents an average surface density or + surface-density-weighted average across a pixel. + + "pencilbeam": a pixel represents a column density or + column-density-weighted average integrated over a pencil + beam through the pixel center. + + Only applies to SPH datasets. """ index = data_source.ds.index if hasattr(index, "meshes") and not isinstance( @@ -236,11 +248,12 @@ def pixelize( elif self.axis_id.get(dimension, dimension) is not None: buff, mask = self._ortho_pixelize( - data_source, field, bounds, size, antialias, dimension, periodic + data_source, field, bounds, size, antialias, dimension, + periodic, pixelmeaning ) else: buff, mask = self._oblique_pixelize( - data_source, field, bounds, size, antialias + data_source, field, bounds, size, antialias, pixelmeaning ) if return_mask: @@ -314,7 +327,8 @@ def pixelize_line(self, field, start_point, end_point, npoints): return arc_length, plot_values def _ortho_pixelize( - self, data_source, field, bounds, size, antialias, dim, periodic + self, data_source, field, bounds, size, antialias, dim, periodic, + pixelmeaning ): from yt.data_objects.construction_data_containers import YTParticleProj from yt.data_objects.selection_objects.slices import YTSlice @@ -396,6 +410,19 @@ def _ortho_pixelize( kernel_name = "cubic" if isinstance(data_source, YTParticleProj): # projection + # pick out projection function. _pencilbeam does have + # one extra optional parameter (minimum number of + # subsampling pixels), but we just use the default value + # for now. + if pixelmeaning == "pencilbeam": + pixelize_sph_kernel_projection = \ + pixelize_sph_kernel_projection_pencilbeam + elif pixelmeaning == "pixelave": + pixelize_sph_kernel_projection = \ + pixelize_sph_kernel_projection_pixelave + else: + raise NotImplementedError( + f"No pixelmeaning option {pixelmeaning}") weight = data_source.weight_field moment = data_source.moment le, re = data_source.data_source.get_bbox() @@ -431,7 +458,7 @@ def _ortho_pixelize( if weight is None: for chunk in proj_reg.chunks([], "io"): data_source._initialize_projected_units([field], chunk) - pixelize_sph_kernel_projection( + pixelize_sph_kernel_projection_pencilbeam( buff, mask_uint8, chunk[ptype, px_name].to("code_length"), @@ -667,7 +694,8 @@ def _ortho_pixelize( assert mask is None or mask.dtype == bool return buff, mask - def _oblique_pixelize(self, data_source, field, bounds, size, antialias): + def _oblique_pixelize(self, data_source, field, bounds, size, antialias, + pixelmeaning): from yt.data_objects.selection_objects.slices import YTCuttingPlane from yt.frontends.sph.data_structures import ParticleDataset from yt.frontends.stream.data_structures import StreamParticlesDataset diff --git a/yt/utilities/answer_testing/framework.py b/yt/utilities/answer_testing/framework.py index 6abd1dc366d..cf6caf12656 100644 --- a/yt/utilities/answer_testing/framework.py +++ b/yt/utilities/answer_testing/framework.py @@ -628,8 +628,10 @@ def __init__(self, ds_fn, axis, field, weight_field=None, obj_type=None): self.obj_type = obj_type def _get_frb(self, obj): + # pixelmeaning="pencilbeam" for backwards compatible SPH tests proj = self.ds.proj( - self.field, self.axis, weight_field=self.weight_field, data_source=obj + self.field, self.axis, weight_field=self.weight_field, + data_source=obj, pixelmeaning="pencilbeam", ) frb = proj.to_frb((1.0, "unitary"), 256) return proj, frb diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx index 31395fd32f4..44e3e1ee749 100644 --- a/yt/utilities/lib/pixelization_routines.pyx +++ b/yt/utilities/lib/pixelization_routines.pyx @@ -2494,7 +2494,8 @@ def off_axis_projection_SPH(np.float64_t[:] px, north_vector, weight_field=None, depth=None, - kernel_name="cubic"): + kernel_name="cubic", + pixelmeaning="pixelave"): # periodic: periodicity of the data set: # Do nothing in event of a 0 normal vector if np.allclose(normal_vector, 0.): @@ -2521,21 +2522,30 @@ def off_axis_projection_SPH(np.float64_t[:] px, # approach implemented for this along-axis projection method # would fail here check_period = np.array([0, 0, 0], dtype="int") - pixelize_sph_kernel_projection(projection_array, - mask, - px_rotated, - py_rotated, - pz_rotated, - smoothing_lengths, - particle_masses, - particle_densities, - quantity_to_smooth, - [rot_bounds_x0, rot_bounds_x1, - rot_bounds_y0, rot_bounds_y1, - rot_bounds_z0, rot_bounds_z1], - weight_field=weight_field, - _check_period=check_period, - kernel_name=kernel_name) + if pixelmeaning == "pencilbeam": + pixelize_sph_kernel_projection_pencilbeam( + projection_array, mask, + px_rotated, py_rotated, pz_rotated, + smoothing_lengths, particle_masses, particle_densities, + quantity_to_smooth, + [rot_bounds_x0, rot_bounds_x1, + rot_bounds_y0, rot_bounds_y1, + rot_bounds_z0, rot_bounds_z1], + weight_field=weight_field, _check_period=check_period, + kernel_name=kernel_name + ) + elif pixelmeaning == "pixelave": + pixelize_sph_kernel_projection_pixelave( + projection_array, mask, + px_rotated, py_rotated, pz_rotated, + smoothing_lengths, particle_masses, particle_densities, + quantity_to_smooth, + [rot_bounds_x0, rot_bounds_x1, + rot_bounds_y0, rot_bounds_y1, + rot_bounds_z0, rot_bounds_z1], + weight_field=weight_field, _check_period=check_period, + kernel_name=kernel_name + ) # like slice pixelization, but for off-axis planes def pixelize_sph_kernel_cutting( diff --git a/yt/visualization/fits_image.py b/yt/visualization/fits_image.py index 59d17acebea..18afd1c8b7b 100644 --- a/yt/visualization/fits_image.py +++ b/yt/visualization/fits_image.py @@ -111,6 +111,15 @@ def __init__( The dataset associated with the image(s), typically used to transfer metadata to the header(s). Does not need to be specified if *data* has a dataset as an attribute. + pixelmeaning: + "pixelav": a pixel represents an average surface density or + surface-density-weighted average across a pixel. + + "pencilbeam": a pixel represents a column density or + column-density-weighted average integrated over a pencil + beam through the pixel center. + + Only applies to SPH datasets. Examples -------- @@ -870,7 +879,9 @@ def sanitize_fits_unit(unit): def construct_image( - ds, axis, data_source, center, image_res, width, length_unit, origin="domain" + ds, axis, data_source, center, image_res, width, length_unit, + origin="domain", + pixelmeaning: {"pixelave", "pencilbeam"} = "pixelave", ): if width is None: width = ds.domain_width[axis_wcs[axis]] @@ -910,9 +921,12 @@ def construct_image( crval = np.zeros(2) if hasattr(data_source, "to_frb"): if is_sequence(axis): - frb = data_source.to_frb(width[0], (nx, ny), height=width[1]) + frb = data_source.to_frb(width[0], (nx, ny), height=width[1], + pixelmeaning=pixelmeaning) else: - frb = data_source.to_frb(width[0], (nx, ny), center=center, height=width[1]) + frb = data_source.to_frb(width[0], (nx, ny), center=center, + height=width[1], + pixelmeaning=pixelmeaning) elif isinstance(data_source, ParticleDummyDataSource): if hasattr(data_source, "normal_vector"): # If we have a normal vector, this means @@ -1138,6 +1152,15 @@ class FITSProjection(FITSImageData): for a weighted projection, moment = 1 (the default) corresponds to a weighted average. moment = 2 corresponds to a weighted standard deviation. + pixelmeaning: + "pixelav": a pixel represents an average surface density or + surface-density-weighted average across a pixel. + + "pencilbeam": a pixel represents a column density or + column-density-weighted average integrated over a pencil + beam through the pixel center. + + Only applies to SPH datasets. """ def __init__( @@ -1153,6 +1176,7 @@ def __init__( *, origin="domain", moment=1, + pixelmeaning: {"pixelave", "pencilbeam" } = "pixelave", **kwargs, ): fields = list(iter_fields(fields)) @@ -1170,6 +1194,7 @@ def __init__( width, length_unit, origin=origin, + pixelmeaning=pixelmeaning, ) super().__init__(frb, fields=fields, length_unit=lunit, wcs=w) @@ -1378,6 +1403,15 @@ class FITSParticleOffAxisProjection(FITSImageData): center coordinates will be the same as the center of the image as defined by the *center* keyword argument. If "image", then the center coordinates will be set to (0,0). Default: "domain" + pixelmeaning: + "pixelav": a pixel represents an average surface density or + surface-density-weighted average across a pixel. + + "pencilbeam": a pixel represents a column density or + column-density-weighted average integrated over a pencil + beam through the pixel center. + + Only applies to SPH datasets. """ def __init__( @@ -1396,6 +1430,7 @@ def __init__( field_parameters=None, data_source=None, north_vector=None, + pixelmeaning: {"pixelave", "pencilbeam"} = "pixelave", ): fields = list(iter_fields(fields)) center, dcenter = ds.coordinates.sanitize_center(center, None) @@ -1427,6 +1462,7 @@ def __init__( width, length_unit, origin="image", + pixelmeaning=pixelmeaning, ) super().__init__(frb, fields=fields, length_unit=lunit, wcs=w) @@ -1603,6 +1639,15 @@ class FITSOffAxisProjection(FITSImageData): for a weighted projection, moment = 1 (the default) corresponds to a weighted average. moment = 2 corresponds to a weighted standard deviation. + pixelmeaning: + "pixelav": a pixel represents an average surface density or + surface-density-weighted average across a pixel. + + "pencilbeam": a pixel represents a column density or + column-density-weighted average integrated over a pencil + beam through the pixel center. + + Only applies to SPH datasets. """ def __init__( @@ -1619,6 +1664,7 @@ def __init__( depth=(1.0, "unitary"), method="integrate", length_unit=None, + *, moment=1, ): @@ -1652,6 +1698,7 @@ def __init__( method=method, weight=weight_field, depth=depth, + pixelmeaning=pixelmeaning, ).swapaxes(0, 1) if moment == 2: diff --git a/yt/visualization/fixed_resolution.py b/yt/visualization/fixed_resolution.py index e5c6178e935..5333d7aa0eb 100644 --- a/yt/visualization/fixed_resolution.py +++ b/yt/visualization/fixed_resolution.py @@ -67,6 +67,15 @@ class FixedResolutionBuffer: periodic : boolean This can be true or false, and governs whether the pixelization will span the domain boundaries. + pixelmeaning: + "pixelav": a pixel represents an average surface density or + surface-density-weighted average across a pixel. + + "pencilbeam": a pixel represents a column density or + column-density-weighted average integrated over a pencil + beam through the pixel center. + + Only applies to SPH datasets. filters : list of FixedResolutionBufferFilter objects (optional) @@ -115,6 +124,7 @@ def __init__( buff_size, antialias=True, periodic=False, + pixelmeaning: {"pencilbeam", "pixelave"} = "pixelave", *, filters: list["FixedResolutionBufferFilter"] | None = None, ): @@ -127,6 +137,7 @@ def __init__( self.mask: dict[str, MaskT] = {} self.axis = data_source.axis self.periodic = periodic + self.pixelmeaning = pixelmeaning self._data_valid = False # import type here to avoid import cycles @@ -184,6 +195,7 @@ def _generate_image_and_mask(self, item) -> None: self.buff_size, int(self.antialias), return_mask=True, + pixelmeaning="pixelmeaning", ) buff = self._apply_filters(buff) @@ -255,6 +267,7 @@ def _get_info(self, item): info["length_unit"] = self.data_source.ds.length_unit info["length_to_cm"] = info["length_unit"].in_cgs().to_ndarray() info["center"] = self.data_source.center + info["pixelmeaning"] = self.pixelmeaning try: info["coord"] = self.data_source.coord @@ -653,6 +666,7 @@ def _generate_image_and_mask(self, item) -> None: north_vector=dd.north_vector, depth=dd.depth, method=dd.method, + pixelmeaning=self.pixelmeaning, ) if self.data_source.moment == 2: @@ -684,6 +698,7 @@ def _sq_field(field, data, item: FieldKey): north_vector=dd.north_vector, depth=dd.depth, method=dd.method, + pixelmeaning=self.pixelmeaning, ) buff = compute_stddev_image(buff2, buff) diff --git a/yt/visualization/volume_rendering/off_axis_projection.py b/yt/visualization/volume_rendering/off_axis_projection.py index 5de55f2fbbc..31ec10a5084 100644 --- a/yt/visualization/volume_rendering/off_axis_projection.py +++ b/yt/visualization/volume_rendering/off_axis_projection.py @@ -33,6 +33,7 @@ def off_axis_projection( depth=None, num_threads=1, method="integrate", + pixelmeaning: {"pencilbeam", "pixelave"} = "pixelave", ): r"""Project through a dataset, off-axis, and return the image plane. @@ -102,6 +103,16 @@ def off_axis_projection( This should only be used for uniform resolution grid datasets, as other datasets may result in unphysical images. or camera movements. + pixelmeaning: + "pixelav": a pixel represents an average surface density or + surface-density-weighted average across a pixel. + + "pencilbeam": a pixel represents a column density or + column-density-weighted average integrated over a pencil beam + through the pixel center. + + Only applies to SPH datasets. + Returns ------- image : array @@ -132,6 +143,9 @@ def off_axis_projection( "Only interpolated=False methods are currently implemented " "for off-axis-projections" ) + if pixelmeaning not in {"pencilbeam", "pixelave"}: + raise ValueError( + f"No pixelmeaning option {pixelmeaning} exists") data_source = data_source_or_all(data_source) @@ -261,6 +275,7 @@ def off_axis_projection( north, depth=depth, kernel_name=kernel_name, + pixelmeaning=pixelmeaning, ) # Assure that the path length unit is in the default length units @@ -303,6 +318,7 @@ def off_axis_projection( weight_field=chunk[weight].in_units(wounits), depth=depth, kernel_name=kernel_name, + pixelmeaning=pixelmeaning, ) for chunk in data_source.chunks([], "io"): @@ -324,6 +340,7 @@ def off_axis_projection( north, depth=depth, kernel_name=kernel_name, + pixelmeaning=pixelmeaning, ) normalization_2d_utility(buf, weight_buff) From 6168a1662ed3df382038e2f9937ff5a45d746fe7 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Fri, 21 Feb 2025 16:36:44 +0100 Subject: [PATCH 07/68] continued work on integrating backend --- yt/visualization/plot_window.py | 55 +++++++++++++++++++++++++++++++++ 1 file changed, 55 insertions(+) diff --git a/yt/visualization/plot_window.py b/yt/visualization/plot_window.py index 9a1bbc9e8e2..8506139f850 100644 --- a/yt/visualization/plot_window.py +++ b/yt/visualization/plot_window.py @@ -182,6 +182,15 @@ class PlotWindow(ImagePlotContainer, abc.ABC): window_size : float The size of the window on the longest axis (in units of inches), including the margins but not the colorbar. + pixelmeaning: + "pixelav": a pixel represents an average surface density or + surface-density-weighted average across a pixel. + + "pencilbeam": a pixel represents a column density or + column-density-weighted average integrated over a pencil + beam through the pixel center. + + Only applies to SPH datasets. """ @@ -199,6 +208,7 @@ def __init__( fontsize=18, aspect=None, setup=False, + pixelmeaning: {"pixelave", "pencilbeam"} = "pixelave", *, geometry: Geometry = Geometry.CARTESIAN, ) -> None: @@ -216,6 +226,7 @@ def __init__( self._axes_unit_names = None self._transform = None self._projection = None + self.pixelmeaning = pixelmeaning self.aspect = aspect skip = list(FixedResolutionBuffer._exclude_fields) + data_source._key_fields @@ -336,6 +347,7 @@ def _recreate_frb(self): self.antialias, periodic=self._periodic, filters=old_filters, + pixelmeaning=self.pixelmeaning, ) # At this point the frb has the valid bounds, size, aliasing, etc. @@ -763,6 +775,28 @@ def set_antialias(self, aa): aa : boolean """ self.antialias = aa + + @invalidate_data + def set_pixelmeaning( + self, + pixelmeaning: {"pixelave", "pencilbeam"} = "pencilbeam", + ): + """ + Change the SPH surface density calculation approach + + parameters + ---------- + pixelmeaning: + "pixelav": a pixel represents an average surface density or + surface-density-weighted average across a pixel. + + "pencilbeam": a pixel represents a column density or + column-density-weighted average integrated over a pencil + beam through the pixel center. + + Only applies to projections of SPH datasets. + """ + self.pixelmeaning = pixelmeaning @invalidate_data def set_buff_size(self, size): @@ -2014,7 +2048,15 @@ class AxisAlignedProjectionPlot(ProjectionPlot, PWViewerMPL): for a weighted projection, moment = 1 (the default) corresponds to a weighted average. moment = 2 corresponds to a weighted standard deviation. + pixelmeaning: + "pixelav": a pixel represents an average surface density or + surface-density-weighted average across a pixel. + "pencilbeam": a pixel represents a column density or + column-density-weighted average integrated over a pencil + beam through the pixel center. + + Only applies to SPH datasets. Examples -------- @@ -2048,6 +2090,7 @@ def __init__( window_size=8.0, buff_size=(800, 800), aspect=None, + pixelmeaning: {"pixelave", "pencilbeam"} = "pixelave", *, moment=1, ): @@ -2112,6 +2155,7 @@ def __init__( aspect=aspect, buff_size=buff_size, geometry=ds.geometry, + pixelmeaning=pixelmeaning, ) if axes_unit is None: axes_unit = get_axes_unit(width, ds) @@ -2427,6 +2471,15 @@ class OffAxisProjectionPlot(ProjectionPlot, PWViewerMPL): Size of the buffer to use for the image, i.e. the number of resolution elements used. Effectively sets a resolution limit to the image if buff_size is smaller than the finest gridding. + pixelmeaning: + "pixelav": a pixel represents an average surface density or + surface-density-weighted average across a pixel. + + "pencilbeam": a pixel represents a column density or + column-density-weighted average integrated over a pencil + beam through the pixel center. + + Only applies to SPH datasets. """ _plot_type = "OffAxisProjection" @@ -2455,6 +2508,7 @@ def __init__( moment=1, data_source=None, buff_size=(800, 800), + pixelmeaning: {"pixelave", "pencilbeam"} = "pixelave", ): if ds.geometry not in self._supported_geometries: raise NotImplementedError( @@ -2555,6 +2609,7 @@ def __init__( oblique=True, fontsize=fontsize, buff_size=buff_size, + pixelmeaning=pixelmeaning, ) if axes_unit is None: axes_unit = get_axes_unit(width, ds) From d7ec60638b44d85f8d1cdda495fa0ece174db0d9 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Fri, 21 Feb 2025 17:09:42 +0100 Subject: [PATCH 08/68] more backend integration, including for tests --- yt/geometry/coordinates/tests/test_sph_pixelization.py | 5 +++++ .../coordinates/tests/test_sph_pixelization_pytestonly.py | 1 + yt/visualization/fits_image.py | 2 +- yt/visualization/tests/test_offaxisprojection_pytestonly.py | 1 + 4 files changed, 8 insertions(+), 1 deletion(-) diff --git a/yt/geometry/coordinates/tests/test_sph_pixelization.py b/yt/geometry/coordinates/tests/test_sph_pixelization.py index 1d1c82be5c8..e7a141518b3 100644 --- a/yt/geometry/coordinates/tests/test_sph_pixelization.py +++ b/yt/geometry/coordinates/tests/test_sph_pixelization.py @@ -66,6 +66,7 @@ def test_sph_projection_basic1(): height=(2.5, "cm"), center=np.array([1.5, 1.5, 1.5]), periodic=False, + pixelmeaning="pencilbeam", ) out = frb.get_image(("gas", "density")) @@ -96,6 +97,7 @@ def test_sph_projection_basic2(): height=(2.5, "cm"), center=np.array([1.375, 1.375, 1.5]), periodic=False, + pixelmeaning="pencilbeam", ) out = frb.get_image(("gas", "density")) @@ -143,6 +145,7 @@ def refmass(i: int, j: int, k: int) -> float: massgenerator=refmass, unitrho=unitrho, bbox=bbox, + pixelmeaning="pencilbeam", ) return ds @@ -163,6 +166,7 @@ def getdata_test_gridproj2(): height=(2.5, "cm"), center=np.array([1.5, 1.5, 1.5]), periodic=False, + pixelmeaning="pencilbeam", ) out = frb.get_image(("gas", "density")) outlist.append(out) @@ -205,6 +209,7 @@ def test_sph_gridproj_reseffect2(): height=(3.0 - 2.0 * margin, "cm"), center=np.array([1.5, 1.5, 1.5]), periodic=False, + pixelmeaning="pencilbeam", ) out = frb.get_image(("gas", "density")) imgs[rl] = out diff --git a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py index ef12dac8242..e9dc3ff0aea 100644 --- a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py +++ b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py @@ -112,6 +112,7 @@ def makemasses(i, j, k): buff_size=(5, 5), center=center, data_source=source, + pixelmeaning="pencilbeam", ) img = prj.frb.data["gas", "density"] if weighted: diff --git a/yt/visualization/fits_image.py b/yt/visualization/fits_image.py index 18afd1c8b7b..d41c6864a31 100644 --- a/yt/visualization/fits_image.py +++ b/yt/visualization/fits_image.py @@ -1664,7 +1664,7 @@ def __init__( depth=(1.0, "unitary"), method="integrate", length_unit=None, - + pixelmeaning: {"pixelave", "pencilbeam"} = "pixelave", *, moment=1, ): diff --git a/yt/visualization/tests/test_offaxisprojection_pytestonly.py b/yt/visualization/tests/test_offaxisprojection_pytestonly.py index 246e5954c84..ac3070ace50 100644 --- a/yt/visualization/tests/test_offaxisprojection_pytestonly.py +++ b/yt/visualization/tests/test_offaxisprojection_pytestonly.py @@ -128,6 +128,7 @@ def makemasses(i, j, k): data_source=source, north_vector=northvector, depth=depth, + pixelmeaning="pencilbeam", ) img = prj.frb.data["gas", "density"] if weighted: From 7f06360b856d992a5c5bb57964b3eac2e0bd2c41 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Mon, 24 Feb 2025 16:26:05 +0100 Subject: [PATCH 09/68] debug until compiles (macOS) --- yt/utilities/lib/pixelization_routines.pyx | 230 +++++++++++---------- 1 file changed, 125 insertions(+), 105 deletions(-) diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx index 44e3e1ee749..91429fb0785 100644 --- a/yt/utilities/lib/pixelization_routines.pyx +++ b/yt/utilities/lib/pixelization_routines.pyx @@ -1360,54 +1360,63 @@ def pixelize_sph_kernel_projection_pixelave( any_float[:] pmass, any_float[:] pdens, any_float[:] quantity_to_smooth, - bounds, + any_float[:] bounds, kernel_name="cubic", weight_field=None, _check_period = (1, 1, 1), period=None, - np.int_t nsample_min = 32): - - cdef np.intp_t xsize, ysize - cdef np.float64_t x_min, x_max, y_min, y_max, z_min, z_max, prefactor_j - cdef np.int64_t xi, yi, x0, x1, y0, y1, xxi, yyi, nx, ny, nsubx, nsuby - cdef np.float64_t q_ij2, posx_diff, posy_diff, ih_j2, qts_j - cdef np.float_64_t insubx, insuby, insmin, xnorm, ynorm, kv - cdef np.float64_t x, y, dx, dy, idx, idy, ipixA, h_j2, px, py, pz + nsample_min=32): + + # shared variables (or used before parallel section) + cdef np.intp_t xsize, ysize, si, sj, sii, sjj, sci + cdef np.float64_t x_min, x_max, y_min, y_max, z_min, z_max + cdef np.float64_t insmin, xnorm, ynorm, q2 + cdef np.float64_t dx, dy, idx, idy, ipixA cdef np.float64_t period_x = 0, period_y = 0, period_z = 0 - cdef int i, j, ii, jj, kk, ic, nsmin + cdef int nsmin cdef np.float64_t[:] _weight_field + cdef np.int8_t[3] check_period + cdef np.int8_t weighted = 0 + # (hopefully?) private variables in parallel section + cdef np.float64_t x, y, px, py, pz + cdef np.float64_t h_j2, ih_j2, qts_j, prefactor_j + cdef np.float64_t q_ij2, posx_diff, posy_diff + cdef np.float64_t insubx, insuby, xn, yn, kv + cdef np.int64_t xi, yi, x0, x1, y0, y1, xxi, yyi + cdef np.int64_t nx, ny, nsubx, nsuby + cdef int i, j, ii, jj, kk, ci, xsubi, ysubi cdef int * xiter cdef int * yiter cdef int * ziter cdef np.float64_t * xiterv cdef np.float64_t * yiterv cdef np.float64_t * ziterv - cdef np.int8_t[3] check_period if weight_field is not None: _weight_field = weight_field + weighted = 1 if period is not None: period_x = period[0] period_y = period[1] period_z = period[2] - for i in range(3): - check_period[i] = np.int8(_check_period[i]) + for si in range(3): + check_period[si] = _check_period[si] # we find the x and y range over which we have pixels and we find how many # pixels we have in each dimension - xsize, ysize = buff.shape[0], buff.shape[1] - x_min = bounds[0] - x_max = bounds[1] - y_min = bounds[2] - y_max = bounds[3] - z_min = bounds[4] - z_max = bounds[5] - - dx = (x_max - x_min) / xsize - dy = (y_max - y_min) / ysize - - idx = 1.0/dx - idy = 1.0/dy + xsize, ysize = buff.shape[0], buff.shape[1] + x_min = bounds[0] + x_max = bounds[1] + y_min = bounds[2] + y_max = bounds[3] + z_min = bounds[4] + z_max = bounds[5] + + dx = (x_max - x_min) / ( xsize) + dy = (y_max - y_min) / ( ysize) + + idx = 1.0 / dx + idy = 1.0 / dy ipixA = idx * idy nsmin = nsample_min @@ -1427,62 +1436,71 @@ def pixelize_sph_kernel_projection_pixelave( # normalize it explicitly, but skip that step for now, since we # also don't do this for larger particles.) kern2by2 = malloc( - sizeof(np.float64_t) * (nsmin + 1) * (nsmin + 1) * 4) + sizeof(np.float64_t) * (nsmin + 1) * (nsmin + 1) * 4) # step 1: just get a big grid of the 1D integral values - kern1dint = malloc( - sizeof(np.float64_t) * nsmin * nsmin) + kern1dint = malloc(sizeof(np.float64_t) * nsmin * nsmin) insmin = 1.0 / nsmin - for i in range(nsmin): - xnorm = -1.0 + 2.0 * (i + 0.5) * insmin + for si in range(nsmin): + xnorm = -1.0 + 2.0 * (si + 0.5) * insmin xnorm = xnorm * xnorm - for j in range(nsmin): - ynorm = -1.0 + 2.0 * (i + 0.5) * insmin - q_ij2 = xnorm + ynorm * ynorm - if q_ij2 >= 1.: - kern1dint[nsample_min * i + j] = 0.0 + for sj in range(nsmin): + ynorm = -1.0 + 2.0 * (sj + 0.5) * insmin + q2 = xnorm + ynorm * ynorm + if q2 >= 1.: + kern1dint[nsmin * si + sj] = 0.0 else: - kern1dint[nsample_min * i + j] = itab.interpolate(q_ij2) + kern1dint[nsmin * si + sj] = itab.interpolate(q2) # step 2: integrate the 1D values for each grid edge position - for i in range(nsample_min + 1): - for j in range(nsample_min + 1): + for si in range(nsmin + 1): + for sj in range(nsmin + 1): # slowest 2 indices - ic = (4 * (nsample_min + 1) * (nsample_min + 1) * i - + 4 * (nsample_min + 1) * j) - kern2by2[ic] = 0.0 - kern2by2[ic + 1] = 0.0 - kern2by2[ic + 2] = 0.0 - kern2by2[ic + 3] = 0.0 - for ii in range(0, i): - for jj in range(0, j): - kern2by2[ic] += kern1dint[nsample_min * ii + jj] - for jj in range(j, nsample_min): - kern2by2[ic + 1] += kern1dint[nsample_min * ii + jj] - for ii in range(i, nsample_min): - for jj in range(0, j): - kern2by2[ic + 2] += kern1dint[nsample_min * ii + jj] - for jj in range(j, nsample_min): - kern2by2[ic + 3] += kern1dint[nsample_min * ii + jj] + sci = (4 * (nsmin + 1) * (nsmin + 1) * si + + 4 * (nsmin + 1) * sj) + kern2by2[sci] = 0.0 + kern2by2[sci + 1] = 0.0 + kern2by2[sci + 2] = 0.0 + kern2by2[sci + 3] = 0.0 + for sii in range(0, si): + for sjj in range(0, sj): + kern2by2[sci] += kern1dint[nsmin * sii + sjj] + for sjj in range(sj, nsmin): + kern2by2[sci + 1] += kern1dint[nsmin * sii + sjj] + for sii in range(si, nsmin): + for sjj in range(0, sj): + kern2by2[sci + 2] += kern1dint[nsmin * sii + sjj] + for sjj in range(sj, nsmin): + kern2by2[sci + 3] += kern1dint[nsmin * sii + sjj] free(kern1dint) with nogil, parallel(): # loop through every particle - # NOTE: this loop can be quite time consuming. However it is easily - # parallelizable in multiple ways, such as: - # 1) use multiple cores to process individual particles (the outer loop) - # 2) use multiple cores to process individual pixels for a given particle - # (the inner loops) - # Depending on the ratio of particles' "sphere of influence" (a.k.a. the smoothing - # length) to the physical width of the pixels, different parallelization - # strategies may yield different speed-ups. Strategy #1 works better in the case - # of lots of itty bitty particles. Strategy #2 works well when we have a - # not-very-large-number of reasonably large-compared-to-pixels particles. We - # currently employ #1 as its workload is more even and consistent, even though it - # comes with a price of an additional, per thread memory for storing the - # intermediate results. - # !!! These comments were written for the line of sight integral - # !!! counterpart to this function, and not checked for this version - - local_buff = malloc(sizeof(np.float64_t) * xsize * ysize) + # NOTE: this loop can be quite time consuming. However it is + # easily parallelizable in multiple ways, such as: + # 1) use multiple cores to process individual particles (the + # outer loop) + # 2) use multiple cores to process individual pixels for a + # given particle (the inner loops) + # Depending on the ratio of particles' "sphere of influence" + # (a.k.a. the smoothing length) to the physical width of the + # pixels, different parallelization strategies may yield + # different speed-ups. Strategy #1 works better in the case of + # lots of itty bitty particles. Strategy #2 works well when we + # have a not-very-large-number of reasonably + # large-compared-to-pixels particles. We currently employ #1 as + # its workload is more even and consistent, even though it + # comes with a price of an additional, per thread memory for + # storing the intermediate results. + # !!! These comments were written for the line of sight + # !!! integral counterpart to this function, and not checked + # for this version + # (also, I don't know if this is possible in cython, but in + # C + OpenMP without pixel subsampling, atomic read/write + # operations to the output array weren't much slower than + # using the memory-intensive ouput array buffer for each + # thread) + + local_buff = malloc(sizeof(np.float64_t) + * xsize * ysize) xiterv = malloc(sizeof(np.float64_t) * 2) yiterv = malloc(sizeof(np.float64_t) * 2) ziterv = malloc(sizeof(np.float64_t) * 2) @@ -1529,11 +1547,13 @@ def pixelize_sph_kernel_projection_pixelave( h_j2 = hsml[j] * hsml[j] ih_j2 = 1.0 / h_j2 - # for small particles, we use the particle to be smoothed + # for small particles, we use the quantity to be smoothed # directly, for larger ones, we need the whole prefactor - qts_j = quantity_to_smooth[j] - if weight_field is not None: - qts_j *= _weight_field[j] + qts_j = quantity_to_smooth[j] + if weighted: + # cython compiler complains when I use " *= ": + # "Cannot read reduction variable in loop body" + qts_j = qts_j * _weight_field[j] prefactor_j = pmass[j] / pdens[j] * ih_j2 * qts_j # Discussion point: do we want the hsml margin on the z direction? @@ -1590,23 +1610,22 @@ def pixelize_sph_kernel_projection_pixelave( # (nsmin + 1, nsmin + 1, 2, 2) # bin edge values -1 -- 1 # in units: - # center of 2x2 output pix grid - # offset from particle center / hsml + # (offset of center of 2x2 output pix + # grid from particle center) / hsml # where does the overlapped pixel edge # fall in the saved kernel space array? - xnorm = (x_min + (x0 + 1) * dx - px) \ - / hsml[j] - ynorm = (y_min + (y0 + 1) * dy - py) \ - / hsml[j] + xn = (x_min + (x0 + 1) * dx - px) / hsml[j] + yn = (y_min + (y0 + 1) * dy - py) / hsml[j] # round to the nearest kernel grid edge - # (C float to int cast behaves like floor()) + # (C float to int cast behaves like + # floor()) xi = ( - 0.5 * (xnorm + 1.) * nsmin + 0.5) + 0.5 * (xn + 1.) * nsmin + 0.5) yi = ( - 0.5 * (ynorm + 1.) * nsmin + 0.5) + 0.5 * (yn + 1.) * nsmin + 0.5) # weight for each line integral: - # area of subsampled array in length units + # area of subsampling pixel in length units # divided by output grid pixel size prefactor_j *= (4.0 * h_j2 * insmin * insmin * ipixA) @@ -1614,16 +1633,18 @@ def pixelize_sph_kernel_projection_pixelave( ci = (4 * nsmin * nsmin * xi + 4 * nsmin * yi) for xxi in range(2): - # catch on the next periodic loop if x0 + xxi < 0 or x0 + xxi > xsize: + # catch on the next periodic + # loop continue for yyi in range(2): - # catch on the next periodic loop if y0 + yyi < 0 or y0 + yyi > ysize: + # catch on the next + # periodic loop continue kv = kern2by2[ci + 2 * xxi + yyi] if kv == 0.: - # no erroneous mask = 1 + # don't set mask to 1 continue local_buff[x0 + xxi + (y0 + yyi) * xsize] += ( @@ -1632,34 +1653,34 @@ def pixelize_sph_kernel_projection_pixelave( mask[xi, yi] = 1 else: - # subsample + # subsample the kernel in each pixel # here we find the pixels which this # particle contributes to # subsampling does not depend on whether a # particle overlaps a (periodic) edge - x0 = ((px - hsml[j] - x_min)*idx) - x1 = ((px + hsml[j] - x_min)*idx) + x0 = ((px - hsml[j] - x_min) * idx) + x1 = ((px + hsml[j] - x_min) * idx) nx = x1 - x0 + 2 x0 = iclip(x0 - 1, 0, xsize) x1 = iclip(x1 + 1, 0, xsize) - y0 = ((py - hsml[j] - y_min)*idy) - y1 = ((py + hsml[j] - y_min)*idy) + y0 = ((py - hsml[j] - y_min) * idy) + y1 = ((py + hsml[j] - y_min) * idy) ny = y1 - y0 + 2 y0 = iclip(y0 - 1, 0, ysize) y1 = iclip(y1 + 1, 0, ysize) - nsubx = ( - ( nsample_min / + nsubx = ( + ( nsmin / nx) + 1) - insubx = 1.0 / nsubx - nsubx = ( - ( nsample_min / + insubx = 1.0 / ( nsubx) + nsuby = ( + ( nsmin / ny) + 1) - insuby = 1.0 / nsuby + insuby = 1.0 / ( nsuby) # found pixels we deposit on, # loop through those pixels @@ -1674,8 +1695,7 @@ def pixelize_sph_kernel_projection_pixelave( for ysubi in range(nsuby): y = (y_min + (yi + (ysubi + 0.5) * insuby) - * dy - ) + * dy) posy_diff = py - y posy_diff = posy_diff * posy_diff if posy_diff > h_j2: continue @@ -1690,9 +1710,9 @@ def pixelize_sph_kernel_projection_pixelave( mask[xi, yi] = 1 with gil: - for xxi in range(xsize): - for yyi in range(ysize): - buff[xxi, yyi] += local_buff[xxi + yyi*xsize] + for sxi in range(xsize): + for syi in range(ysize): + buff[sxi, syi] += local_buff[sxi + syi * xsize] free(local_buff) free(kern2by2) free(xiterv) From 9a0b3600084b8adbbba3568c6a35c32b0db6954d Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Mon, 24 Feb 2025 16:31:15 +0100 Subject: [PATCH 10/68] fix docstring typo --- yt/visualization/fits_image.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yt/visualization/fits_image.py b/yt/visualization/fits_image.py index d41c6864a31..fdbde7f462b 100644 --- a/yt/visualization/fits_image.py +++ b/yt/visualization/fits_image.py @@ -112,7 +112,7 @@ def __init__( to transfer metadata to the header(s). Does not need to be specified if *data* has a dataset as an attribute. pixelmeaning: - "pixelav": a pixel represents an average surface density or + "pixelave": a pixel represents an average surface density or surface-density-weighted average across a pixel. "pencilbeam": a pixel represents a column density or From dc8153b19e44535babb1572f9535660ea6e0e484 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Mon, 24 Feb 2025 16:32:13 +0100 Subject: [PATCH 11/68] more docstring typos --- yt/visualization/fits_image.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/yt/visualization/fits_image.py b/yt/visualization/fits_image.py index fdbde7f462b..0b6677af087 100644 --- a/yt/visualization/fits_image.py +++ b/yt/visualization/fits_image.py @@ -1153,7 +1153,7 @@ class FITSProjection(FITSImageData): weighted average. moment = 2 corresponds to a weighted standard deviation. pixelmeaning: - "pixelav": a pixel represents an average surface density or + "e": a pixel represents an average surface density or surface-density-weighted average across a pixel. "pencilbeam": a pixel represents a column density or @@ -1404,7 +1404,7 @@ class FITSParticleOffAxisProjection(FITSImageData): defined by the *center* keyword argument. If "image", then the center coordinates will be set to (0,0). Default: "domain" pixelmeaning: - "pixelav": a pixel represents an average surface density or + "pixelave": a pixel represents an average surface density or surface-density-weighted average across a pixel. "pencilbeam": a pixel represents a column density or @@ -1640,7 +1640,7 @@ class FITSOffAxisProjection(FITSImageData): weighted average. moment = 2 corresponds to a weighted standard deviation. pixelmeaning: - "pixelav": a pixel represents an average surface density or + "pixelave": a pixel represents an average surface density or surface-density-weighted average across a pixel. "pencilbeam": a pixel represents a column density or From df809a7980e03ebda18f15b5e9e4ce4bb80bfe16 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Mon, 24 Feb 2025 16:35:45 +0100 Subject: [PATCH 12/68] more docstring typos --- yt/visualization/plot_window.py | 8 ++++---- yt/visualization/volume_rendering/off_axis_projection.py | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/yt/visualization/plot_window.py b/yt/visualization/plot_window.py index 8506139f850..f634459a475 100644 --- a/yt/visualization/plot_window.py +++ b/yt/visualization/plot_window.py @@ -183,7 +183,7 @@ class PlotWindow(ImagePlotContainer, abc.ABC): The size of the window on the longest axis (in units of inches), including the margins but not the colorbar. pixelmeaning: - "pixelav": a pixel represents an average surface density or + "pixelave": a pixel represents an average surface density or surface-density-weighted average across a pixel. "pencilbeam": a pixel represents a column density or @@ -787,7 +787,7 @@ def set_pixelmeaning( parameters ---------- pixelmeaning: - "pixelav": a pixel represents an average surface density or + "pixelave": a pixel represents an average surface density or surface-density-weighted average across a pixel. "pencilbeam": a pixel represents a column density or @@ -2049,7 +2049,7 @@ class AxisAlignedProjectionPlot(ProjectionPlot, PWViewerMPL): weighted average. moment = 2 corresponds to a weighted standard deviation. pixelmeaning: - "pixelav": a pixel represents an average surface density or + "pixelave": a pixel represents an average surface density or surface-density-weighted average across a pixel. "pencilbeam": a pixel represents a column density or @@ -2472,7 +2472,7 @@ class OffAxisProjectionPlot(ProjectionPlot, PWViewerMPL): used. Effectively sets a resolution limit to the image if buff_size is smaller than the finest gridding. pixelmeaning: - "pixelav": a pixel represents an average surface density or + "pixelave": a pixel represents an average surface density or surface-density-weighted average across a pixel. "pencilbeam": a pixel represents a column density or diff --git a/yt/visualization/volume_rendering/off_axis_projection.py b/yt/visualization/volume_rendering/off_axis_projection.py index 31ec10a5084..82593fa7ef0 100644 --- a/yt/visualization/volume_rendering/off_axis_projection.py +++ b/yt/visualization/volume_rendering/off_axis_projection.py @@ -104,7 +104,7 @@ def off_axis_projection( datasets may result in unphysical images. or camera movements. pixelmeaning: - "pixelav": a pixel represents an average surface density or + "pixelave": a pixel represents an average surface density or surface-density-weighted average across a pixel. "pencilbeam": a pixel represents a column density or From 0ae4e73eda49befd85616a9bddc7a087f83f4662 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Mon, 24 Feb 2025 16:54:47 +0100 Subject: [PATCH 13/68] fix type annotations? --- .../selection_objects/data_selection_objects.py | 3 ++- yt/geometry/coordinates/cartesian_coordinates.py | 6 ++++-- yt/visualization/fits_image.py | 11 ++++++----- yt/visualization/fixed_resolution.py | 6 +++--- yt/visualization/plot_window.py | 10 +++++----- .../volume_rendering/off_axis_projection.py | 5 +++-- 6 files changed, 23 insertions(+), 18 deletions(-) diff --git a/yt/data_objects/selection_objects/data_selection_objects.py b/yt/data_objects/selection_objects/data_selection_objects.py index bcc95306c7a..a9517c7d9d2 100644 --- a/yt/data_objects/selection_objects/data_selection_objects.py +++ b/yt/data_objects/selection_objects/data_selection_objects.py @@ -7,6 +7,7 @@ import numpy as np from more_itertools import always_iterable +from typing import Literal from unyt import unyt_array from unyt.exceptions import UnitConversionError, UnitParseError @@ -569,7 +570,7 @@ def _get_pw(self, fields, center, width, origin, plot_type): def to_frb(self, width, resolution, center=None, height=None, periodic=False, - pixelmeaning: {"pixelave", "pencilbeam"} = "pixelave"): + pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave"): r"""This function returns a FixedResolutionBuffer generated from this object. diff --git a/yt/geometry/coordinates/cartesian_coordinates.py b/yt/geometry/coordinates/cartesian_coordinates.py index f6bcb431bf5..58588d2fa41 100644 --- a/yt/geometry/coordinates/cartesian_coordinates.py +++ b/yt/geometry/coordinates/cartesian_coordinates.py @@ -1,3 +1,5 @@ +from typing import Literal + import numpy as np from yt.data_objects.index_subobjects.unstructured_mesh import SemiStructuredMesh @@ -171,7 +173,7 @@ def pixelize( size, antialias=True, periodic=True, - pixelmeaning: {"pixelave", "pencilbeam"} = "pixelave", + pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave", *, return_mask=False, ): @@ -181,7 +183,7 @@ def pixelize( routines written in cython pixelmeaning: - "pixelav": a pixel represents an average surface density or + "pixelave": a pixel represents an average surface density or surface-density-weighted average across a pixel. "pencilbeam": a pixel represents a column density or diff --git a/yt/visualization/fits_image.py b/yt/visualization/fits_image.py index 0b6677af087..73b11ba904b 100644 --- a/yt/visualization/fits_image.py +++ b/yt/visualization/fits_image.py @@ -5,6 +5,7 @@ import numpy as np from more_itertools import first, mark_ends +from typing import Literal from yt._typing import FieldKey from yt.data_objects.construction_data_containers import YTCoveringGrid @@ -881,7 +882,7 @@ def sanitize_fits_unit(unit): def construct_image( ds, axis, data_source, center, image_res, width, length_unit, origin="domain", - pixelmeaning: {"pixelave", "pencilbeam"} = "pixelave", + pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave", ): if width is None: width = ds.domain_width[axis_wcs[axis]] @@ -1153,7 +1154,7 @@ class FITSProjection(FITSImageData): weighted average. moment = 2 corresponds to a weighted standard deviation. pixelmeaning: - "e": a pixel represents an average surface density or + "pixelave": a pixel represents an average surface density or surface-density-weighted average across a pixel. "pencilbeam": a pixel represents a column density or @@ -1176,7 +1177,7 @@ def __init__( *, origin="domain", moment=1, - pixelmeaning: {"pixelave", "pencilbeam" } = "pixelave", + pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave", **kwargs, ): fields = list(iter_fields(fields)) @@ -1430,7 +1431,7 @@ def __init__( field_parameters=None, data_source=None, north_vector=None, - pixelmeaning: {"pixelave", "pencilbeam"} = "pixelave", + pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave", ): fields = list(iter_fields(fields)) center, dcenter = ds.coordinates.sanitize_center(center, None) @@ -1664,7 +1665,7 @@ def __init__( depth=(1.0, "unitary"), method="integrate", length_unit=None, - pixelmeaning: {"pixelave", "pencilbeam"} = "pixelave", + pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave", *, moment=1, ): diff --git a/yt/visualization/fixed_resolution.py b/yt/visualization/fixed_resolution.py index 5333d7aa0eb..5dda14b7b3e 100644 --- a/yt/visualization/fixed_resolution.py +++ b/yt/visualization/fixed_resolution.py @@ -1,7 +1,7 @@ import sys import weakref from functools import partial -from typing import TYPE_CHECKING +from typing import Literal, TYPE_CHECKING import numpy as np @@ -68,7 +68,7 @@ class FixedResolutionBuffer: This can be true or false, and governs whether the pixelization will span the domain boundaries. pixelmeaning: - "pixelav": a pixel represents an average surface density or + "pixelave": a pixel represents an average surface density or surface-density-weighted average across a pixel. "pencilbeam": a pixel represents a column density or @@ -124,7 +124,7 @@ def __init__( buff_size, antialias=True, periodic=False, - pixelmeaning: {"pencilbeam", "pixelave"} = "pixelave", + pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave", *, filters: list["FixedResolutionBufferFilter"] | None = None, ): diff --git a/yt/visualization/plot_window.py b/yt/visualization/plot_window.py index f634459a475..0ace6fb500f 100644 --- a/yt/visualization/plot_window.py +++ b/yt/visualization/plot_window.py @@ -2,7 +2,7 @@ import sys from collections import defaultdict from numbers import Number -from typing import TYPE_CHECKING, Union +from typing import Literal, TYPE_CHECKING, Union import matplotlib import numpy as np @@ -208,7 +208,7 @@ def __init__( fontsize=18, aspect=None, setup=False, - pixelmeaning: {"pixelave", "pencilbeam"} = "pixelave", + pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave", *, geometry: Geometry = Geometry.CARTESIAN, ) -> None: @@ -779,7 +779,7 @@ def set_antialias(self, aa): @invalidate_data def set_pixelmeaning( self, - pixelmeaning: {"pixelave", "pencilbeam"} = "pencilbeam", + pixelmeaning: Literal["pixelave", "pencilbeam"] = "pencilbeam", ): """ Change the SPH surface density calculation approach @@ -2090,7 +2090,7 @@ def __init__( window_size=8.0, buff_size=(800, 800), aspect=None, - pixelmeaning: {"pixelave", "pencilbeam"} = "pixelave", + pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave", *, moment=1, ): @@ -2508,7 +2508,7 @@ def __init__( moment=1, data_source=None, buff_size=(800, 800), - pixelmeaning: {"pixelave", "pencilbeam"} = "pixelave", + pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave", ): if ds.geometry not in self._supported_geometries: raise NotImplementedError( diff --git a/yt/visualization/volume_rendering/off_axis_projection.py b/yt/visualization/volume_rendering/off_axis_projection.py index 82593fa7ef0..62dedebabfb 100644 --- a/yt/visualization/volume_rendering/off_axis_projection.py +++ b/yt/visualization/volume_rendering/off_axis_projection.py @@ -1,4 +1,5 @@ import numpy as np +from typing import Literal from yt.data_objects.api import ImageArray from yt.funcs import is_sequence, mylog @@ -33,7 +34,7 @@ def off_axis_projection( depth=None, num_threads=1, method="integrate", - pixelmeaning: {"pencilbeam", "pixelave"} = "pixelave", + pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave", ): r"""Project through a dataset, off-axis, and return the image plane. @@ -143,7 +144,7 @@ def off_axis_projection( "Only interpolated=False methods are currently implemented " "for off-axis-projections" ) - if pixelmeaning not in {"pencilbeam", "pixelave"}: + if pixelmeaning not in["pencilbeam", "pixelave"]: raise ValueError( f"No pixelmeaning option {pixelmeaning} exists") From ce5fa63863f55087bcdab1ebcdfb8a4e068dc44f Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 24 Feb 2025 15:56:59 +0000 Subject: [PATCH 14/68] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../data_selection_objects.py | 22 +++++++----- .../coordinates/cartesian_coordinates.py | 35 +++++++++++-------- yt/utilities/answer_testing/framework.py | 7 ++-- yt/utilities/lib/pixelization_routines.pyx | 22 ++++++------ yt/visualization/fits_image.py | 33 +++++++++++------ yt/visualization/fixed_resolution.py | 4 +-- yt/visualization/plot_window.py | 18 +++++----- .../volume_rendering/off_axis_projection.py | 14 ++++---- 8 files changed, 91 insertions(+), 64 deletions(-) diff --git a/yt/data_objects/selection_objects/data_selection_objects.py b/yt/data_objects/selection_objects/data_selection_objects.py index a9517c7d9d2..e99ebed4e5c 100644 --- a/yt/data_objects/selection_objects/data_selection_objects.py +++ b/yt/data_objects/selection_objects/data_selection_objects.py @@ -4,10 +4,10 @@ import uuid from collections import defaultdict from contextlib import contextmanager +from typing import Literal import numpy as np from more_itertools import always_iterable -from typing import Literal from unyt import unyt_array from unyt.exceptions import UnitConversionError, UnitParseError @@ -568,9 +568,15 @@ def _get_pw(self, fields, center, width, origin, plot_type): pw._setup_plots() return pw - def to_frb(self, width, resolution, center=None, height=None, - periodic=False, - pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave"): + def to_frb( + self, + width, + resolution, + center=None, + height=None, + periodic=False, + pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave", + ): r"""This function returns a FixedResolutionBuffer generated from this object. @@ -602,7 +608,7 @@ def to_frb(self, width, resolution, center=None, height=None, "pixelav": a pixel represents an average surface density or surface-density-weighted average across a pixel. - "pencilbeam": a pixel represents a column density or + "pencilbeam": a pixel represents a column density or column-density-weighted average integrated over a pencil beam through the pixel center. @@ -670,9 +676,9 @@ def to_frb(self, width, resolution, center=None, height=None, center[yax] - height * 0.5, center[yax] + height * 0.5, ) - frb = FixedResolutionBuffer(self, bounds, resolution, - periodic=periodic, - pixelmeaning=pixelmeaning) + frb = FixedResolutionBuffer( + self, bounds, resolution, periodic=periodic, pixelmeaning=pixelmeaning + ) return frb diff --git a/yt/geometry/coordinates/cartesian_coordinates.py b/yt/geometry/coordinates/cartesian_coordinates.py index 58588d2fa41..1611078edd7 100644 --- a/yt/geometry/coordinates/cartesian_coordinates.py +++ b/yt/geometry/coordinates/cartesian_coordinates.py @@ -186,7 +186,7 @@ def pixelize( "pixelave": a pixel represents an average surface density or surface-density-weighted average across a pixel. - "pencilbeam": a pixel represents a column density or + "pencilbeam": a pixel represents a column density or column-density-weighted average integrated over a pencil beam through the pixel center. @@ -250,8 +250,14 @@ def pixelize( elif self.axis_id.get(dimension, dimension) is not None: buff, mask = self._ortho_pixelize( - data_source, field, bounds, size, antialias, dimension, - periodic, pixelmeaning + data_source, + field, + bounds, + size, + antialias, + dimension, + periodic, + pixelmeaning, ) else: buff, mask = self._oblique_pixelize( @@ -329,8 +335,7 @@ def pixelize_line(self, field, start_point, end_point, npoints): return arc_length, plot_values def _ortho_pixelize( - self, data_source, field, bounds, size, antialias, dim, periodic, - pixelmeaning + self, data_source, field, bounds, size, antialias, dim, periodic, pixelmeaning ): from yt.data_objects.construction_data_containers import YTParticleProj from yt.data_objects.selection_objects.slices import YTSlice @@ -412,19 +417,20 @@ def _ortho_pixelize( kernel_name = "cubic" if isinstance(data_source, YTParticleProj): # projection - # pick out projection function. _pencilbeam does have - # one extra optional parameter (minimum number of + # pick out projection function. _pencilbeam does have + # one extra optional parameter (minimum number of # subsampling pixels), but we just use the default value - # for now. + # for now. if pixelmeaning == "pencilbeam": - pixelize_sph_kernel_projection = \ + pixelize_sph_kernel_projection = ( pixelize_sph_kernel_projection_pencilbeam + ) elif pixelmeaning == "pixelave": - pixelize_sph_kernel_projection = \ + pixelize_sph_kernel_projection = ( pixelize_sph_kernel_projection_pixelave + ) else: - raise NotImplementedError( - f"No pixelmeaning option {pixelmeaning}") + raise NotImplementedError(f"No pixelmeaning option {pixelmeaning}") weight = data_source.weight_field moment = data_source.moment le, re = data_source.data_source.get_bbox() @@ -696,8 +702,9 @@ def _ortho_pixelize( assert mask is None or mask.dtype == bool return buff, mask - def _oblique_pixelize(self, data_source, field, bounds, size, antialias, - pixelmeaning): + def _oblique_pixelize( + self, data_source, field, bounds, size, antialias, pixelmeaning + ): from yt.data_objects.selection_objects.slices import YTCuttingPlane from yt.frontends.sph.data_structures import ParticleDataset from yt.frontends.stream.data_structures import StreamParticlesDataset diff --git a/yt/utilities/answer_testing/framework.py b/yt/utilities/answer_testing/framework.py index cf6caf12656..29977c3e21b 100644 --- a/yt/utilities/answer_testing/framework.py +++ b/yt/utilities/answer_testing/framework.py @@ -630,8 +630,11 @@ def __init__(self, ds_fn, axis, field, weight_field=None, obj_type=None): def _get_frb(self, obj): # pixelmeaning="pencilbeam" for backwards compatible SPH tests proj = self.ds.proj( - self.field, self.axis, weight_field=self.weight_field, - data_source=obj, pixelmeaning="pencilbeam", + self.field, + self.axis, + weight_field=self.weight_field, + data_source=obj, + pixelmeaning="pencilbeam", ) frb = proj.to_frb((1.0, "unitary"), 256) return proj, frb diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx index 91429fb0785..8e1a2e25beb 100644 --- a/yt/utilities/lib/pixelization_routines.pyx +++ b/yt/utilities/lib/pixelization_routines.pyx @@ -1366,7 +1366,7 @@ def pixelize_sph_kernel_projection_pixelave( _check_period = (1, 1, 1), period=None, nsample_min=32): - + # shared variables (or used before parallel section) cdef np.intp_t xsize, ysize, si, sj, sii, sjj, sci cdef np.float64_t x_min, x_max, y_min, y_max, z_min, z_max @@ -1474,32 +1474,32 @@ def pixelize_sph_kernel_projection_pixelave( with nogil, parallel(): # loop through every particle - # NOTE: this loop can be quite time consuming. However it is + # NOTE: this loop can be quite time consuming. However it is # easily parallelizable in multiple ways, such as: # 1) use multiple cores to process individual particles (the # outer loop) - # 2) use multiple cores to process individual pixels for a + # 2) use multiple cores to process individual pixels for a # given particle (the inner loops) - # Depending on the ratio of particles' "sphere of influence" + # Depending on the ratio of particles' "sphere of influence" # (a.k.a. the smoothing length) to the physical width of the - # pixels, different parallelization strategies may yield + # pixels, different parallelization strategies may yield # different speed-ups. Strategy #1 works better in the case of # lots of itty bitty particles. Strategy #2 works well when we - # have a not-very-large-number of reasonably + # have a not-very-large-number of reasonably # large-compared-to-pixels particles. We currently employ #1 as # its workload is more even and consistent, even though it - # comes with a price of an additional, per thread memory for + # comes with a price of an additional, per thread memory for # storing the intermediate results. - # !!! These comments were written for the line of sight + # !!! These comments were written for the line of sight # !!! integral counterpart to this function, and not checked # for this version # (also, I don't know if this is possible in cython, but in - # C + OpenMP without pixel subsampling, atomic read/write + # C + OpenMP without pixel subsampling, atomic read/write # operations to the output array weren't much slower than # using the memory-intensive ouput array buffer for each # thread) - local_buff = malloc(sizeof(np.float64_t) + local_buff = malloc(sizeof(np.float64_t) * xsize * ysize) xiterv = malloc(sizeof(np.float64_t) * 2) yiterv = malloc(sizeof(np.float64_t) * 2) @@ -1618,7 +1618,7 @@ def pixelize_sph_kernel_projection_pixelave( xn = (x_min + (x0 + 1) * dx - px) / hsml[j] yn = (y_min + (y0 + 1) * dy - py) / hsml[j] # round to the nearest kernel grid edge - # (C float to int cast behaves like + # (C float to int cast behaves like # floor()) xi = ( 0.5 * (xn + 1.) * nsmin + 0.5) diff --git a/yt/visualization/fits_image.py b/yt/visualization/fits_image.py index 73b11ba904b..ee1ff4771e5 100644 --- a/yt/visualization/fits_image.py +++ b/yt/visualization/fits_image.py @@ -2,10 +2,10 @@ import sys from functools import partial from numbers import Number as numeric_type +from typing import Literal import numpy as np from more_itertools import first, mark_ends -from typing import Literal from yt._typing import FieldKey from yt.data_objects.construction_data_containers import YTCoveringGrid @@ -116,7 +116,7 @@ def __init__( "pixelave": a pixel represents an average surface density or surface-density-weighted average across a pixel. - "pencilbeam": a pixel represents a column density or + "pencilbeam": a pixel represents a column density or column-density-weighted average integrated over a pencil beam through the pixel center. @@ -880,7 +880,13 @@ def sanitize_fits_unit(unit): def construct_image( - ds, axis, data_source, center, image_res, width, length_unit, + ds, + axis, + data_source, + center, + image_res, + width, + length_unit, origin="domain", pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave", ): @@ -922,12 +928,17 @@ def construct_image( crval = np.zeros(2) if hasattr(data_source, "to_frb"): if is_sequence(axis): - frb = data_source.to_frb(width[0], (nx, ny), height=width[1], - pixelmeaning=pixelmeaning) + frb = data_source.to_frb( + width[0], (nx, ny), height=width[1], pixelmeaning=pixelmeaning + ) else: - frb = data_source.to_frb(width[0], (nx, ny), center=center, - height=width[1], - pixelmeaning=pixelmeaning) + frb = data_source.to_frb( + width[0], + (nx, ny), + center=center, + height=width[1], + pixelmeaning=pixelmeaning, + ) elif isinstance(data_source, ParticleDummyDataSource): if hasattr(data_source, "normal_vector"): # If we have a normal vector, this means @@ -1157,7 +1168,7 @@ class FITSProjection(FITSImageData): "pixelave": a pixel represents an average surface density or surface-density-weighted average across a pixel. - "pencilbeam": a pixel represents a column density or + "pencilbeam": a pixel represents a column density or column-density-weighted average integrated over a pencil beam through the pixel center. @@ -1408,7 +1419,7 @@ class FITSParticleOffAxisProjection(FITSImageData): "pixelave": a pixel represents an average surface density or surface-density-weighted average across a pixel. - "pencilbeam": a pixel represents a column density or + "pencilbeam": a pixel represents a column density or column-density-weighted average integrated over a pencil beam through the pixel center. @@ -1644,7 +1655,7 @@ class FITSOffAxisProjection(FITSImageData): "pixelave": a pixel represents an average surface density or surface-density-weighted average across a pixel. - "pencilbeam": a pixel represents a column density or + "pencilbeam": a pixel represents a column density or column-density-weighted average integrated over a pencil beam through the pixel center. diff --git a/yt/visualization/fixed_resolution.py b/yt/visualization/fixed_resolution.py index 5dda14b7b3e..3e04ea2c482 100644 --- a/yt/visualization/fixed_resolution.py +++ b/yt/visualization/fixed_resolution.py @@ -1,7 +1,7 @@ import sys import weakref from functools import partial -from typing import Literal, TYPE_CHECKING +from typing import TYPE_CHECKING, Literal import numpy as np @@ -71,7 +71,7 @@ class FixedResolutionBuffer: "pixelave": a pixel represents an average surface density or surface-density-weighted average across a pixel. - "pencilbeam": a pixel represents a column density or + "pencilbeam": a pixel represents a column density or column-density-weighted average integrated over a pencil beam through the pixel center. diff --git a/yt/visualization/plot_window.py b/yt/visualization/plot_window.py index 0ace6fb500f..b22f4d4254c 100644 --- a/yt/visualization/plot_window.py +++ b/yt/visualization/plot_window.py @@ -2,7 +2,7 @@ import sys from collections import defaultdict from numbers import Number -from typing import Literal, TYPE_CHECKING, Union +from typing import TYPE_CHECKING, Literal, Union import matplotlib import numpy as np @@ -186,7 +186,7 @@ class PlotWindow(ImagePlotContainer, abc.ABC): "pixelave": a pixel represents an average surface density or surface-density-weighted average across a pixel. - "pencilbeam": a pixel represents a column density or + "pencilbeam": a pixel represents a column density or column-density-weighted average integrated over a pencil beam through the pixel center. @@ -775,12 +775,12 @@ def set_antialias(self, aa): aa : boolean """ self.antialias = aa - + @invalidate_data def set_pixelmeaning( - self, - pixelmeaning: Literal["pixelave", "pencilbeam"] = "pencilbeam", - ): + self, + pixelmeaning: Literal["pixelave", "pencilbeam"] = "pencilbeam", + ): """ Change the SPH surface density calculation approach @@ -790,7 +790,7 @@ def set_pixelmeaning( "pixelave": a pixel represents an average surface density or surface-density-weighted average across a pixel. - "pencilbeam": a pixel represents a column density or + "pencilbeam": a pixel represents a column density or column-density-weighted average integrated over a pencil beam through the pixel center. @@ -2052,7 +2052,7 @@ class AxisAlignedProjectionPlot(ProjectionPlot, PWViewerMPL): "pixelave": a pixel represents an average surface density or surface-density-weighted average across a pixel. - "pencilbeam": a pixel represents a column density or + "pencilbeam": a pixel represents a column density or column-density-weighted average integrated over a pencil beam through the pixel center. @@ -2475,7 +2475,7 @@ class OffAxisProjectionPlot(ProjectionPlot, PWViewerMPL): "pixelave": a pixel represents an average surface density or surface-density-weighted average across a pixel. - "pencilbeam": a pixel represents a column density or + "pencilbeam": a pixel represents a column density or column-density-weighted average integrated over a pencil beam through the pixel center. diff --git a/yt/visualization/volume_rendering/off_axis_projection.py b/yt/visualization/volume_rendering/off_axis_projection.py index 62dedebabfb..e760ba739eb 100644 --- a/yt/visualization/volume_rendering/off_axis_projection.py +++ b/yt/visualization/volume_rendering/off_axis_projection.py @@ -1,6 +1,7 @@ -import numpy as np from typing import Literal +import numpy as np + from yt.data_objects.api import ImageArray from yt.funcs import is_sequence, mylog from yt.geometry.oct_geometry_handler import OctreeIndex @@ -108,12 +109,12 @@ def off_axis_projection( "pixelave": a pixel represents an average surface density or surface-density-weighted average across a pixel. - "pencilbeam": a pixel represents a column density or + "pencilbeam": a pixel represents a column density or column-density-weighted average integrated over a pencil beam through the pixel center. - + Only applies to SPH datasets. - + Returns ------- image : array @@ -144,9 +145,8 @@ def off_axis_projection( "Only interpolated=False methods are currently implemented " "for off-axis-projections" ) - if pixelmeaning not in["pencilbeam", "pixelave"]: - raise ValueError( - f"No pixelmeaning option {pixelmeaning} exists") + if pixelmeaning not in ["pencilbeam", "pixelave"]: + raise ValueError(f"No pixelmeaning option {pixelmeaning} exists") data_source = data_source_or_all(data_source) From 3868858fcdfe29b0ce36e039101fb30e86a9b317 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Mon, 24 Feb 2025 18:42:22 +0100 Subject: [PATCH 15/68] first go on-axes pixelave proj test setup --- .../tests/test_sph_pixelization_pytestonly.py | 143 ++++++++++++++++++ yt/testing.py | 57 ++++++- 2 files changed, 199 insertions(+), 1 deletion(-) diff --git a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py index e9dc3ff0aea..8a8ef313be8 100644 --- a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py +++ b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py @@ -11,6 +11,7 @@ fake_random_sph_ds, fake_sph_flexible_grid_ds, integrate_kernel, + pixelintegrate_kernel, ) @@ -160,6 +161,148 @@ def makemasses(i, j, k): # print(img.v) assert_rel_equal(expected_out, img.v, 5) +@pytest.mark.parametrize("weighted", [True, False]) +@pytest.mark.parametrize("periodic", [True, False]) +@pytest.mark.parametrize("depth", [None, (1.0, "cm"), (5.0, "cm")]) +@pytest.mark.parametrize("shiftcenter", [False, True]) +@pytest.mark.parametrize("axis", [0, 1, 2]) +@pytest.mark.parametrize("hsmlfac", [0.2, 1.0, 2.0]) +def test_sph_projection_pixelave_alongaxes( + axis: int, + shiftcenter: bool, + periodic: bool, + depth: float, + weighted: bool, + hsmlfac: float, +) -> None: + if shiftcenter: + center = unyt.unyt_array(np.array((0.625, 0.525, 0.425)), "cm") + else: + center = unyt.unyt_array(np.array((1.5, 1.5, 1.5)), "cm") + bbox = unyt.unyt_array(np.array([[0.0, 3.0], [0.0, 3.0], [0.0, 3.0]]), "cm") + hsml_factor = hsmlfac + unitrho = 1.5 + + if axis == 2: + ax0 = 0 + ax1 = 1 + elif axis == 0: + ax0 = 1 + ax1 = 2 + elif axis == 1: + ax0 = 2 + ax1 = 0 + buff_size = (5, 5) + # test correct centering, particle selection + def makemasses(i, j, k): + if i == j == k == 1: + return 2.0 + else: + return 1.0 + + # result shouldn't depend explicitly on the center if we re-center + # the data, unless we get cut-offs in the non-periodic case + ds = fake_sph_flexible_grid_ds( + hsml_factor=hsml_factor, + nperside=3, + periodic=periodic, + offsets=np.full(3, 0.5), + massgenerator=makemasses, + unitrho=unitrho, + bbox=bbox.v, + recenter=center.v, + ) + if depth is None: + source = ds.all_data() + else: + depth = unyt.unyt_quantity(*depth) + le = np.array(ds.domain_left_edge) + re = np.array(ds.domain_right_edge) + le[axis] = center[axis] - 0.5 * depth + re[axis] = center[axis] + 0.5 * depth + cen = 0.5 * (le + re) + reg = YTRegion(center=cen, left_edge=le, right_edge=re, ds=ds) + source = reg + + # we don't actually want a plot, it's just a straightforward, + # common way to get an frb / image array + if weighted: + toweight_field = ("gas", "density") + else: + toweight_field = None + prj = yt.ProjectionPlot( + ds, + axis, + ("gas", "density"), + width=(2.5, "cm"), + weight_field=toweight_field, + buff_size=buff_size, + center=center, + data_source=source, + pixelmeaning="pixelave", + ) + img = prj.frb.data[("gas", "density")] + + projwidth = unyt.unyt_array(np.array((2.5,)), "cm") + pixedges_x = np.linspace(center[ax0] - 0.5 * projwidth, + center[ax0] + 0.5 * projwidth, + buff_size[0] + 1) + pixedges_y = np.linspace(center[ax1] - 0.5 * projwidth, + center[ax1] + 0.5 * projwidth, + buff_size[1] + 1) + _depth = depth if depth is not None else np.inf + zlims = (center[axis] - 0.5 * _depth, + center[axis] + 0.5 * _depth) + ad = ds.all_data() + pcenxy = np.array( + [ + (ad[("gas", "xyz"[ax0])]).to("cm"), + (ad[("gas", "xyz"[ax1])]).to("cm"), + ] + ).T + projz = (ad[("gas", "xyz"[axis])]).to("cm") + hsml = (ad[("gas", "smoothing_length")]).to("cm") + mass = (ad[("gas", "mass")]).to("g") + density = (ad[("gas", "density")]).to("g * cm**-3") + baseline = np.zeros(buff_size) + if periodic: + periodxy = (bbox[ax0][1] - bbox[ax0][0], + bbox[ax1][1] - bbox[ax1][0]) + periodz = bbox[axis][1] - bbox[axis][0] + else: + periodxy = (None, None) + for i in range(buff_size[0]): + for j in range(buff_size[1]): + pixelxy = np.array([pixedges_x[i], pixedges_x[i + 1], + pixedges_y[j], pixedges_y[j + 1]]) + weightsum = 0. + weightedsum = 0. + for p in range(len(pcenxy.shape[0])): + # check if the particle is excluded along z + if projz[p] < zlims[0] or projz[p] > zlims[1]: + if periodz is None: + continue + else: + dz = projz - center[axis].to("cm") + dz += 0.5 * periodz + dz %= periodz + dz -= periodz + if np.abs(dz) > 0.5 * depth: + continue + + kernint = pixelintegrate_kernel( + cubicspline_python, pixelxy, pcenxy[p], hsml[p], + nsample=100, periodxy=periodxy, + ) + weightsum += kernint * mass[p] + if weighted: + weightedsum += kernint * mass[p] * density[p] + if weighted: + baseline[i, j] += weightedsum / weightsum + else: + baseline[i, j] += weightsum + assert_rel_equal(baseline, img.v, 5) + @pytest.mark.parametrize("periodic", [True, False]) @pytest.mark.parametrize("shiftcenter", [False, True]) diff --git a/yt/testing.py b/yt/testing.py index 5a1c82f1e2f..bb54de774e0 100644 --- a/yt/testing.py +++ b/yt/testing.py @@ -141,7 +141,7 @@ def integrate_kernel( Returns: -------- the integral of the SPH kernel function. - units: 1 / units of b and hsml + units: 1 / (units of b and hsml)**2 """ pre = 1.0 / hsml**2 x = b / hsml @@ -154,6 +154,61 @@ def integrate_kernel( integral = np.sum(spv * dx, axis=0) return np.atleast_1d(pre * integral) +def pixelintegrate_kernel( + kernelfunc: Callable[[float], float], + pixelxy: np.array[float], + pcenxy: np.array[float], + hsml: float, + nsample: int = 100, + periodxy: tuple[float | None, float | None] = (True, True), +) -> float: + """ + integrates a kernel function over a rectangular prism. Assumes + The kernel support is entirely within the prism in the + line-of-sight direction. Does not handle periodicity. + + Parameters: + ----------- + kernelfunc: + the kernel function to integrate + pixelxy: + corners of the pixel in the output grid; + (xmin, xmax, ymin, ymax); [length units] + pcenxy: + sph particle center in the output grid plane + (x, y) [same units as pixelxy] + hsml: + smoothing length [same units as pixelxy] + nsample: + number of subsamples used within the pixel + in each grid direction + + Returns: + -------- + the integral of the SPH kernel function. + units: (1 / input length units)**2 + """ + xedges = np.linspace(pixelxy[0], pixelxy[1], nsample + 1) + xcens = 0.5 * (xedges[:-1] + xedges[1:]) + dx = np.diff(xedges) + yedges = np.linspace(pixelxy[2], pixelxy[3], nsample + 1) + ycens = 0.5 * (yedges[:-1] + yedges[1:]) + dy = np.diff(yedges) + + xoff = xcens - pcenxy[0] + if periodxy[0] is not None: + xoff += 0.5 * periodxy[0] + xoff %= periodxy[0] + xoff -= 0.5 * periodxy[0] + yoff = ycens - pcenxy[1] + if periodxy[1] is not None: + yoff += 0.5 * periodxy[1] + yoff %= periodxy[1] + yoff -= 0.5 * periodxy[1] + bpars = np.sqrt(xoff[:, np.newaxis]**2 + yoff[np.newaxis, :]**2) + lineints = integrate_kernel(kernelfunc, bpars, hsml) + volint = np.sum(lineints * dx[:, np.newaxis] * dy[np.newaxis, :]) + return volint / (pixelxy[1] - pixelxy[0]) / (pixelxy[3] - pixelxy[2]) _zeroperiods = np.array([0.0, 0.0, 0.0]) From 593b0ed237be88102527c6f8d115b89b172aab49 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Tue, 25 Feb 2025 16:58:11 +0100 Subject: [PATCH 16/68] function annotation fixes? --- yt/testing.py | 28 +++++++++++++++++++--------- 1 file changed, 19 insertions(+), 9 deletions(-) diff --git a/yt/testing.py b/yt/testing.py index bb54de774e0..366e4b61a0c 100644 --- a/yt/testing.py +++ b/yt/testing.py @@ -156,19 +156,22 @@ def integrate_kernel( def pixelintegrate_kernel( kernelfunc: Callable[[float], float], - pixelxy: np.array[float], - pcenxy: np.array[float], + pixelxy: np.ndarray[float], + pcenxy: np.ndarray[float], hsml: float, nsample: int = 100, periodxy: tuple[float | None, float | None] = (True, True), ) -> float: """ - integrates a kernel function over a rectangular prism. Assumes - The kernel support is entirely within the prism in the - line-of-sight direction. Does not handle periodicity. + integrates a kernel function over a rectangular prism. + + Returns (kernel_volume_fraction) / pixel area, + where the kernel_volume_fraction is the volume integral of the + kernel over the rectangular prism defined by the pixel, divided + by the volume integral of the kernel over its entire support. - Parameters: - ----------- + Parameters + ---------- kernelfunc: the kernel function to integrate pixelxy: @@ -183,10 +186,17 @@ def pixelintegrate_kernel( number of subsamples used within the pixel in each grid direction - Returns: + Returns -------- - the integral of the SPH kernel function. + integral of the SPH kernel function / pixel area units: (1 / input length units)**2 + + Note + ---- + The calculation assumes the kernel support is entirely + within the rectangular prism defined by the pixel in the + line-of-sight direction, as the backend functions for SPH + projections do. """ xedges = np.linspace(pixelxy[0], pixelxy[1], nsample + 1) xcens = 0.5 * (xedges[:-1] + xedges[1:]) From f8f52df65ff063ce5fd59fd6e6a747bf06dc95c5 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Tue, 25 Feb 2025 16:58:29 +0100 Subject: [PATCH 17/68] raise error if invalid pixelmeaning arg --- yt/utilities/lib/pixelization_routines.pyx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx index 8e1a2e25beb..31fb9339d4c 100644 --- a/yt/utilities/lib/pixelization_routines.pyx +++ b/yt/utilities/lib/pixelization_routines.pyx @@ -2566,6 +2566,8 @@ def off_axis_projection_SPH(np.float64_t[:] px, weight_field=weight_field, _check_period=check_period, kernel_name=kernel_name ) + else: + raise ValueError(f"no pixelmeaning option '{pixelmeaning}'") # like slice pixelization, but for off-axis planes def pixelize_sph_kernel_cutting( From cc2510a7e04159bf25894f712031d416796cc5e0 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 25 Feb 2025 16:19:57 +0000 Subject: [PATCH 18/68] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../tests/test_sph_pixelization_pytestonly.py | 41 +++++++++++-------- yt/testing.py | 14 ++++--- 2 files changed, 31 insertions(+), 24 deletions(-) diff --git a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py index 8a8ef313be8..732438b7341 100644 --- a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py +++ b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py @@ -161,6 +161,7 @@ def makemasses(i, j, k): # print(img.v) assert_rel_equal(expected_out, img.v, 5) + @pytest.mark.parametrize("weighted", [True, False]) @pytest.mark.parametrize("periodic", [True, False]) @pytest.mark.parametrize("depth", [None, (1.0, "cm"), (5.0, "cm")]) @@ -193,6 +194,7 @@ def test_sph_projection_pixelave_alongaxes( ax0 = 2 ax1 = 0 buff_size = (5, 5) + # test correct centering, particle selection def makemasses(i, j, k): if i == j == k == 1: @@ -244,15 +246,14 @@ def makemasses(i, j, k): img = prj.frb.data[("gas", "density")] projwidth = unyt.unyt_array(np.array((2.5,)), "cm") - pixedges_x = np.linspace(center[ax0] - 0.5 * projwidth, - center[ax0] + 0.5 * projwidth, - buff_size[0] + 1) - pixedges_y = np.linspace(center[ax1] - 0.5 * projwidth, - center[ax1] + 0.5 * projwidth, - buff_size[1] + 1) + pixedges_x = np.linspace( + center[ax0] - 0.5 * projwidth, center[ax0] + 0.5 * projwidth, buff_size[0] + 1 + ) + pixedges_y = np.linspace( + center[ax1] - 0.5 * projwidth, center[ax1] + 0.5 * projwidth, buff_size[1] + 1 + ) _depth = depth if depth is not None else np.inf - zlims = (center[axis] - 0.5 * _depth, - center[axis] + 0.5 * _depth) + zlims = (center[axis] - 0.5 * _depth, center[axis] + 0.5 * _depth) ad = ds.all_data() pcenxy = np.array( [ @@ -266,17 +267,17 @@ def makemasses(i, j, k): density = (ad[("gas", "density")]).to("g * cm**-3") baseline = np.zeros(buff_size) if periodic: - periodxy = (bbox[ax0][1] - bbox[ax0][0], - bbox[ax1][1] - bbox[ax1][0]) + periodxy = (bbox[ax0][1] - bbox[ax0][0], bbox[ax1][1] - bbox[ax1][0]) periodz = bbox[axis][1] - bbox[axis][0] else: periodxy = (None, None) for i in range(buff_size[0]): for j in range(buff_size[1]): - pixelxy = np.array([pixedges_x[i], pixedges_x[i + 1], - pixedges_y[j], pixedges_y[j + 1]]) - weightsum = 0. - weightedsum = 0. + pixelxy = np.array( + [pixedges_x[i], pixedges_x[i + 1], pixedges_y[j], pixedges_y[j + 1]] + ) + weightsum = 0.0 + weightedsum = 0.0 for p in range(len(pcenxy.shape[0])): # check if the particle is excluded along z if projz[p] < zlims[0] or projz[p] > zlims[1]: @@ -289,10 +290,14 @@ def makemasses(i, j, k): dz -= periodz if np.abs(dz) > 0.5 * depth: continue - + kernint = pixelintegrate_kernel( - cubicspline_python, pixelxy, pcenxy[p], hsml[p], - nsample=100, periodxy=periodxy, + cubicspline_python, + pixelxy, + pcenxy[p], + hsml[p], + nsample=100, + periodxy=periodxy, ) weightsum += kernint * mass[p] if weighted: @@ -302,7 +307,7 @@ def makemasses(i, j, k): else: baseline[i, j] += weightsum assert_rel_equal(baseline, img.v, 5) - + @pytest.mark.parametrize("periodic", [True, False]) @pytest.mark.parametrize("shiftcenter", [False, True]) diff --git a/yt/testing.py b/yt/testing.py index 366e4b61a0c..5772f9e1818 100644 --- a/yt/testing.py +++ b/yt/testing.py @@ -154,6 +154,7 @@ def integrate_kernel( integral = np.sum(spv * dx, axis=0) return np.atleast_1d(pre * integral) + def pixelintegrate_kernel( kernelfunc: Callable[[float], float], pixelxy: np.ndarray[float], @@ -163,8 +164,8 @@ def pixelintegrate_kernel( periodxy: tuple[float | None, float | None] = (True, True), ) -> float: """ - integrates a kernel function over a rectangular prism. - + integrates a kernel function over a rectangular prism. + Returns (kernel_volume_fraction) / pixel area, where the kernel_volume_fraction is the volume integral of the kernel over the rectangular prism defined by the pixel, divided @@ -174,11 +175,11 @@ def pixelintegrate_kernel( ---------- kernelfunc: the kernel function to integrate - pixelxy: + pixelxy: corners of the pixel in the output grid; (xmin, xmax, ymin, ymax); [length units] pcenxy: - sph particle center in the output grid plane + sph particle center in the output grid plane (x, y) [same units as pixelxy] hsml: smoothing length [same units as pixelxy] @@ -193,7 +194,7 @@ def pixelintegrate_kernel( Note ---- - The calculation assumes the kernel support is entirely + The calculation assumes the kernel support is entirely within the rectangular prism defined by the pixel in the line-of-sight direction, as the backend functions for SPH projections do. @@ -215,11 +216,12 @@ def pixelintegrate_kernel( yoff += 0.5 * periodxy[1] yoff %= periodxy[1] yoff -= 0.5 * periodxy[1] - bpars = np.sqrt(xoff[:, np.newaxis]**2 + yoff[np.newaxis, :]**2) + bpars = np.sqrt(xoff[:, np.newaxis] ** 2 + yoff[np.newaxis, :] ** 2) lineints = integrate_kernel(kernelfunc, bpars, hsml) volint = np.sum(lineints * dx[:, np.newaxis] * dy[np.newaxis, :]) return volint / (pixelxy[1] - pixelxy[0]) / (pixelxy[3] - pixelxy[2]) + _zeroperiods = np.array([0.0, 0.0, 0.0]) From 9bd50cec0fac6b2f5f971bcfd4bc8361d25c654d Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Tue, 25 Feb 2025 18:01:56 +0100 Subject: [PATCH 19/68] bugfix: pass on pixelmeaning value, not string "pixelmeaning" --- yt/visualization/fixed_resolution.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yt/visualization/fixed_resolution.py b/yt/visualization/fixed_resolution.py index 3e04ea2c482..316fa411f96 100644 --- a/yt/visualization/fixed_resolution.py +++ b/yt/visualization/fixed_resolution.py @@ -195,7 +195,7 @@ def _generate_image_and_mask(self, item) -> None: self.buff_size, int(self.antialias), return_mask=True, - pixelmeaning="pixelmeaning", + pixelmeaning=self.pixelmeaning, ) buff = self._apply_filters(buff) From 3e9d5bca2c4c2cd43bf10a32fa923050a420f766 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Tue, 25 Feb 2025 18:05:53 +0100 Subject: [PATCH 20/68] remove wrongfully added pixelmeaning arg --- yt/geometry/coordinates/tests/test_sph_pixelization.py | 1 - 1 file changed, 1 deletion(-) diff --git a/yt/geometry/coordinates/tests/test_sph_pixelization.py b/yt/geometry/coordinates/tests/test_sph_pixelization.py index e7a141518b3..95bf5bc7c5f 100644 --- a/yt/geometry/coordinates/tests/test_sph_pixelization.py +++ b/yt/geometry/coordinates/tests/test_sph_pixelization.py @@ -145,7 +145,6 @@ def refmass(i: int, j: int, k: int) -> float: massgenerator=refmass, unitrho=unitrho, bbox=bbox, - pixelmeaning="pencilbeam", ) return ds From adfc60fadc01acafac81ac991e14ee5947c30a39 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Tue, 25 Feb 2025 18:24:44 +0100 Subject: [PATCH 21/68] first try: test for both pixelmeaning opts. --- .../tests/test_off_axis_SPH.py | 57 +++++++++++++------ 1 file changed, 41 insertions(+), 16 deletions(-) diff --git a/yt/visualization/volume_rendering/tests/test_off_axis_SPH.py b/yt/visualization/volume_rendering/tests/test_off_axis_SPH.py index 57f6f85d522..3ed798f1f38 100644 --- a/yt/visualization/volume_rendering/tests/test_off_axis_SPH.py +++ b/yt/visualization/volume_rendering/tests/test_off_axis_SPH.py @@ -1,20 +1,25 @@ import numpy as np from numpy.testing import assert_almost_equal +import pytest from yt.testing import fake_sph_orientation_ds, requires_module -from yt.utilities.lib.pixelization_routines import pixelize_sph_kernel_projection +from yt.utilities.lib.pixelization_routines import ( + pixelize_sph_kernel_projection_pencilbeam, + pixelize_sph_kernel_projection_pixelave +) from yt.utilities.on_demand_imports import _scipy from yt.visualization.volume_rendering import off_axis_projection as OffAP spatial = _scipy.spatial ndimage = _scipy.ndimage - -def test_no_rotation(): +@pytest.mark.parametrize("pixelmeaning", ["pixelave", "pencilbeam"]) +def test_no_rotation(pixelmeaning): """Determines if a projection processed through off_axis_projection with no rotation will give the same image buffer if processed directly through pixelize_sph_kernel_projection + (done for both pixelmeaning options) """ normal_vector = [0.0, 0.0, 1.0] resolution = (64, 64) @@ -36,16 +41,23 @@ def test_no_rotation(): buf2 = np.zeros(resolution) mask = np.ones_like(buf2, dtype="uint8") buf1 = OffAP.off_axis_projection( - ds, center, normal_vector, width, resolution, ("gas", "density") + ds, center, normal_vector, width, resolution, ("gas", "density"), + pixelmeaning=pixelmeaning, ) - pixelize_sph_kernel_projection( + if pixelmeaning == "pixelave": + onaxisfunc = pixelize_sph_kernel_projection_pixelave + elif pixelmeaning == "pencilbeam": + onaxisfunc = pixelize_sph_kernel_projection_pencilbeam + + onaxisfunc( buf2, mask, px, py, pz, hsml, mass, density, quantity_to_smooth, bounds ) assert_almost_equal(buf1.ndarray_view(), buf2) @requires_module("scipy") -def test_basic_rotation_1(): +@pytest.mark.parametrize("pixelmeaning", ["pixelave", "pencilbeam"]) +def test_basic_rotation_1(pixelmeaning): """All particles on Z-axis should now be on the negative Y-Axis fake_sph_orientation has three z-axis particles, so there should be three y-axis particles after rotation @@ -76,12 +88,14 @@ def test_basic_rotation_1(): resolution, ("gas", "density"), north_vector=north_vector, + pixelmeaning=pixelmeaning, ) find_compare_maxima(expected_maxima, buf1, resolution, width) @requires_module("scipy") -def test_basic_rotation_2(): +@pytest.mark.parametrize("pixelmeaning", ["pixelave", "pencilbeam"]) +def test_basic_rotation_2(pixelmeaning): """Rotation of x-axis onto z-axis. All particles on z-axis should now be on the negative x-Axis fake_sph_orientation has three z-axis particles, so there should be three x-axis particles after rotation @@ -115,12 +129,14 @@ def test_basic_rotation_2(): resolution, ("gas", "density"), north_vector=north_vector, + pixelmeaning=pixelmeaning, ) find_compare_maxima(expected_maxima, buf1, resolution, width) @requires_module("scipy") -def test_basic_rotation_3(): +@pytest.mark.parametrize("pixelmeaning", ["pixelave", "pencilbeam"]) +def test_basic_rotation_3(pixelmeaning): """Rotation of z-axis onto negative z-axis. All fake particles on z-axis should now be of the negative z-Axis. fake_sph_orientation has three z-axis particles, @@ -145,13 +161,15 @@ def test_basic_rotation_3(): center = (left_edge + right_edge) / 2 width = right_edge - left_edge buf1 = OffAP.off_axis_projection( - ds, center, normal_vector, width, resolution, ("gas", "density") + ds, center, normal_vector, width, resolution, ("gas", "density"), + pixelmeaning=pixelmeaning, ) find_compare_maxima(expected_maxima, buf1, resolution, width) @requires_module("scipy") -def test_basic_rotation_4(): +@pytest.mark.parametrize("pixelmeaning", ["pixelave", "pencilbeam"]) +def test_basic_rotation_4(pixelmeaning): """Rotation of x-axis to z-axis and original z-axis to y-axis with the use of the north_vector. All fake particles on z-axis should now be on the y-Axis. All fake particles on the x-axis should now be on the z-axis, and @@ -185,12 +203,14 @@ def test_basic_rotation_4(): resolution, ("gas", "density"), north_vector=north_vector, + pixelmeaning=pixelmeaning, ) find_compare_maxima(expected_maxima, buf1, resolution, width) @requires_module("scipy") -def test_center_1(): +@pytest.mark.parametrize("pixelmeaning", ["pixelave", "pencilbeam"]) +def test_center_1(pixelmeaning): """Change the center to [0, 3, 0] Every point will be shifted by 3 in the y-domain With this, we should not be able to see any of the y-axis particles @@ -214,13 +234,15 @@ def test_center_1(): center = [0.0, 3.0, 0.0] width = right_edge - left_edge buf1 = OffAP.off_axis_projection( - ds, center, normal_vector, width, resolution, ("gas", "density") + ds, center, normal_vector, width, resolution, ("gas", "density"), + pixelmeaning=pixelmeaning, ) find_compare_maxima(expected_maxima, buf1, resolution, width) @requires_module("scipy") -def test_center_2(): +@pytest.mark.parametrize("pixelmeaning", ["pixelave", "pencilbeam"]) +def test_center_2(pixelmeaning): """Change the center to [0, -1, 0] Every point will be shifted by 1 in the y-domain With this, we should not be able to see any of the y-axis particles @@ -241,13 +263,15 @@ def test_center_2(): center = [0.0, -1.0, 0.0] width = right_edge - left_edge buf1 = OffAP.off_axis_projection( - ds, center, normal_vector, width, resolution, ("gas", "density") + ds, center, normal_vector, width, resolution, ("gas", "density"), + pixelmeaning=pixelmeaning, ) find_compare_maxima(expected_maxima, buf1, resolution, width) @requires_module("scipy") -def test_center_3(): +@pytest.mark.parametrize("pixelmeaning", ["pixelave", "pencilbeam"]) +def test_center_3(pixelmeaning): """Change the center to the left edge, or [0, -8, 0] Every point will be shifted by 8 in the y-domain With this, we should not be able to see anything ! @@ -265,7 +289,8 @@ def test_center_3(): (right_edge[2] - left_edge[2]), ] buf1 = OffAP.off_axis_projection( - ds, center, normal_vector, width, resolution, ("gas", "density") + ds, center, normal_vector, width, resolution, ("gas", "density"), + pixelmeaning=pixelmeaning, ) find_compare_maxima(expected_maxima, buf1, resolution, width) From cc39df4135a192fc8c7ba15a809664b8b44fa254 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Tue, 25 Feb 2025 18:35:46 +0100 Subject: [PATCH 22/68] add pixelmeaning arg to avoid errors in pixelize calls --- yt/geometry/coordinates/coordinate_handler.py | 3 +++ yt/geometry/coordinates/cylindrical_coordinates.py | 8 ++++++++ yt/geometry/coordinates/geographic_coordinates.py | 8 ++++++++ yt/geometry/coordinates/spherical_coordinates.py | 8 ++++++++ 4 files changed, 27 insertions(+) diff --git a/yt/geometry/coordinates/coordinate_handler.py b/yt/geometry/coordinates/coordinate_handler.py index 159c4d0e028..28dcb2af484 100644 --- a/yt/geometry/coordinates/coordinate_handler.py +++ b/yt/geometry/coordinates/coordinate_handler.py @@ -156,6 +156,7 @@ def pixelize( size, antialias=True, periodic=True, + pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave", *, return_mask: Literal[False], ) -> "np.ndarray[Any, np.dtype[np.float64]]": ... @@ -170,6 +171,7 @@ def pixelize( size, antialias=True, periodic=True, + pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave", *, return_mask: Literal[True], ) -> tuple[ @@ -186,6 +188,7 @@ def pixelize( size, antialias=True, periodic=True, + pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave", *, return_mask=False, ): diff --git a/yt/geometry/coordinates/cylindrical_coordinates.py b/yt/geometry/coordinates/cylindrical_coordinates.py index 37c95211e5f..9c7385e50dd 100644 --- a/yt/geometry/coordinates/cylindrical_coordinates.py +++ b/yt/geometry/coordinates/cylindrical_coordinates.py @@ -1,4 +1,5 @@ from functools import cached_property +from typing import Literal import numpy as np @@ -87,9 +88,16 @@ def pixelize( size, antialias=True, periodic=False, + pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave", *, return_mask=False, ): + """ + Parameters + ---------- + pixelmeaning: ignored, argument meant for cartesian SPH data + + """ # Note that above, we set periodic by default to be *false*. This is # because our pixelizers, at present, do not handle periodicity # correctly, and if you change the "width" of a cylindrical plot, it diff --git a/yt/geometry/coordinates/geographic_coordinates.py b/yt/geometry/coordinates/geographic_coordinates.py index 036a28ac757..03b31a20e6b 100644 --- a/yt/geometry/coordinates/geographic_coordinates.py +++ b/yt/geometry/coordinates/geographic_coordinates.py @@ -1,4 +1,5 @@ import numpy as np +from typing import Literal import unyt from yt.utilities.lib.pixelization_routines import pixelize_cartesian, pixelize_cylinder @@ -228,9 +229,16 @@ def pixelize( size, antialias=True, periodic=True, + pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave", *, return_mask=False, ): + """ + Parameters + ---------- + pixelmeaning: ignored, argument meant for cartesian SPH data + + """ if self.axis_name[dimension] in ("latitude", "longitude"): buff, mask = self._cyl_pixelize( data_source, field, bounds, size, antialias, dimension diff --git a/yt/geometry/coordinates/spherical_coordinates.py b/yt/geometry/coordinates/spherical_coordinates.py index 165287393a2..3c15b7ac77f 100644 --- a/yt/geometry/coordinates/spherical_coordinates.py +++ b/yt/geometry/coordinates/spherical_coordinates.py @@ -1,4 +1,5 @@ from functools import cached_property +from typing import Literal import numpy as np @@ -90,9 +91,16 @@ def pixelize( size, antialias=True, periodic=True, + pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave", *, return_mask=False, ): + """ + Parameters + ---------- + pixelmeaning: ignored, argument meant for cartesian SPH data + + """ self.period name = self.axis_name[dimension] if name == "r": From fc626686340ad9111ef76556577668b68560c995 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Tue, 25 Feb 2025 18:45:01 +0100 Subject: [PATCH 23/68] add pixelmeaning arg to to_frb to avoid errors --- yt/data_objects/selection_objects/slices.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/yt/data_objects/selection_objects/slices.py b/yt/data_objects/selection_objects/slices.py index fd510737f30..77c7807e919 100644 --- a/yt/data_objects/selection_objects/slices.py +++ b/yt/data_objects/selection_objects/slices.py @@ -1,4 +1,5 @@ import numpy as np +from typing import Literal from yt.data_objects.selection_objects.data_selection_objects import ( YTSelectionContainer, @@ -313,7 +314,8 @@ def to_pw(self, fields=None, center="center", width=None, axes_unit=None): pw._setup_plots() return pw - def to_frb(self, width, resolution, height=None, periodic=False): + def to_frb(self, width, resolution, height=None, periodic=False, + pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave"): r"""This function returns a FixedResolutionBuffer generated from this object. @@ -338,6 +340,7 @@ def to_frb(self, width, resolution, height=None, periodic=False): periodic : boolean This can be true or false, and governs whether the pixelization will span the domain boundaries. + pixelmeaning: ignored, arg. meant for SPH projections Returns ------- From 6308cf4e40a5b7f2854686450f2ff20741e582d4 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Tue, 25 Feb 2025 18:52:47 +0100 Subject: [PATCH 24/68] add pixelmeaning arg. for ParticleImageBuffer to avoid errors --- yt/visualization/fixed_resolution.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/yt/visualization/fixed_resolution.py b/yt/visualization/fixed_resolution.py index 316fa411f96..5a81deb9cd6 100644 --- a/yt/visualization/fixed_resolution.py +++ b/yt/visualization/fixed_resolution.py @@ -717,6 +717,10 @@ class ParticleImageBuffer(FixedResolutionBuffer): that supports particle plots. It splats points onto an image buffer. + parameters + ---------- + pixelmeaning: ignored; arg. is meant for SPH projection plots + """ def __init__( @@ -726,6 +730,7 @@ def __init__( buff_size, antialias=True, periodic=False, + pixelmeaning="pixelave", *, filters=None, ): From 4833b910bf99bc8906f4724aa115b21077fb944a Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 25 Feb 2025 18:01:09 +0000 Subject: [PATCH 25/68] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- yt/data_objects/selection_objects/slices.py | 13 ++++-- .../coordinates/cylindrical_coordinates.py | 2 +- .../coordinates/geographic_coordinates.py | 5 ++- .../coordinates/spherical_coordinates.py | 2 +- .../tests/test_off_axis_SPH.py | 44 ++++++++++++++----- 5 files changed, 49 insertions(+), 17 deletions(-) diff --git a/yt/data_objects/selection_objects/slices.py b/yt/data_objects/selection_objects/slices.py index 77c7807e919..ed690517241 100644 --- a/yt/data_objects/selection_objects/slices.py +++ b/yt/data_objects/selection_objects/slices.py @@ -1,6 +1,7 @@ -import numpy as np from typing import Literal +import numpy as np + from yt.data_objects.selection_objects.data_selection_objects import ( YTSelectionContainer, YTSelectionContainer2D, @@ -314,8 +315,14 @@ def to_pw(self, fields=None, center="center", width=None, axes_unit=None): pw._setup_plots() return pw - def to_frb(self, width, resolution, height=None, periodic=False, - pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave"): + def to_frb( + self, + width, + resolution, + height=None, + periodic=False, + pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave", + ): r"""This function returns a FixedResolutionBuffer generated from this object. diff --git a/yt/geometry/coordinates/cylindrical_coordinates.py b/yt/geometry/coordinates/cylindrical_coordinates.py index 9c7385e50dd..9ec03d998e7 100644 --- a/yt/geometry/coordinates/cylindrical_coordinates.py +++ b/yt/geometry/coordinates/cylindrical_coordinates.py @@ -96,7 +96,7 @@ def pixelize( Parameters ---------- pixelmeaning: ignored, argument meant for cartesian SPH data - + """ # Note that above, we set periodic by default to be *false*. This is # because our pixelizers, at present, do not handle periodicity diff --git a/yt/geometry/coordinates/geographic_coordinates.py b/yt/geometry/coordinates/geographic_coordinates.py index 03b31a20e6b..befe5347427 100644 --- a/yt/geometry/coordinates/geographic_coordinates.py +++ b/yt/geometry/coordinates/geographic_coordinates.py @@ -1,5 +1,6 @@ -import numpy as np from typing import Literal + +import numpy as np import unyt from yt.utilities.lib.pixelization_routines import pixelize_cartesian, pixelize_cylinder @@ -237,7 +238,7 @@ def pixelize( Parameters ---------- pixelmeaning: ignored, argument meant for cartesian SPH data - + """ if self.axis_name[dimension] in ("latitude", "longitude"): buff, mask = self._cyl_pixelize( diff --git a/yt/geometry/coordinates/spherical_coordinates.py b/yt/geometry/coordinates/spherical_coordinates.py index 3c15b7ac77f..1a52c304f0e 100644 --- a/yt/geometry/coordinates/spherical_coordinates.py +++ b/yt/geometry/coordinates/spherical_coordinates.py @@ -99,7 +99,7 @@ def pixelize( Parameters ---------- pixelmeaning: ignored, argument meant for cartesian SPH data - + """ self.period name = self.axis_name[dimension] diff --git a/yt/visualization/volume_rendering/tests/test_off_axis_SPH.py b/yt/visualization/volume_rendering/tests/test_off_axis_SPH.py index 3ed798f1f38..7f2d359a764 100644 --- a/yt/visualization/volume_rendering/tests/test_off_axis_SPH.py +++ b/yt/visualization/volume_rendering/tests/test_off_axis_SPH.py @@ -1,11 +1,11 @@ import numpy as np -from numpy.testing import assert_almost_equal import pytest +from numpy.testing import assert_almost_equal from yt.testing import fake_sph_orientation_ds, requires_module from yt.utilities.lib.pixelization_routines import ( pixelize_sph_kernel_projection_pencilbeam, - pixelize_sph_kernel_projection_pixelave + pixelize_sph_kernel_projection_pixelave, ) from yt.utilities.on_demand_imports import _scipy from yt.visualization.volume_rendering import off_axis_projection as OffAP @@ -13,6 +13,7 @@ spatial = _scipy.spatial ndimage = _scipy.ndimage + @pytest.mark.parametrize("pixelmeaning", ["pixelave", "pencilbeam"]) def test_no_rotation(pixelmeaning): """Determines if a projection processed through @@ -41,7 +42,12 @@ def test_no_rotation(pixelmeaning): buf2 = np.zeros(resolution) mask = np.ones_like(buf2, dtype="uint8") buf1 = OffAP.off_axis_projection( - ds, center, normal_vector, width, resolution, ("gas", "density"), + ds, + center, + normal_vector, + width, + resolution, + ("gas", "density"), pixelmeaning=pixelmeaning, ) if pixelmeaning == "pixelave": @@ -49,9 +55,7 @@ def test_no_rotation(pixelmeaning): elif pixelmeaning == "pencilbeam": onaxisfunc = pixelize_sph_kernel_projection_pencilbeam - onaxisfunc( - buf2, mask, px, py, pz, hsml, mass, density, quantity_to_smooth, bounds - ) + onaxisfunc(buf2, mask, px, py, pz, hsml, mass, density, quantity_to_smooth, bounds) assert_almost_equal(buf1.ndarray_view(), buf2) @@ -161,7 +165,12 @@ def test_basic_rotation_3(pixelmeaning): center = (left_edge + right_edge) / 2 width = right_edge - left_edge buf1 = OffAP.off_axis_projection( - ds, center, normal_vector, width, resolution, ("gas", "density"), + ds, + center, + normal_vector, + width, + resolution, + ("gas", "density"), pixelmeaning=pixelmeaning, ) find_compare_maxima(expected_maxima, buf1, resolution, width) @@ -234,7 +243,12 @@ def test_center_1(pixelmeaning): center = [0.0, 3.0, 0.0] width = right_edge - left_edge buf1 = OffAP.off_axis_projection( - ds, center, normal_vector, width, resolution, ("gas", "density"), + ds, + center, + normal_vector, + width, + resolution, + ("gas", "density"), pixelmeaning=pixelmeaning, ) find_compare_maxima(expected_maxima, buf1, resolution, width) @@ -263,7 +277,12 @@ def test_center_2(pixelmeaning): center = [0.0, -1.0, 0.0] width = right_edge - left_edge buf1 = OffAP.off_axis_projection( - ds, center, normal_vector, width, resolution, ("gas", "density"), + ds, + center, + normal_vector, + width, + resolution, + ("gas", "density"), pixelmeaning=pixelmeaning, ) find_compare_maxima(expected_maxima, buf1, resolution, width) @@ -289,7 +308,12 @@ def test_center_3(pixelmeaning): (right_edge[2] - left_edge[2]), ] buf1 = OffAP.off_axis_projection( - ds, center, normal_vector, width, resolution, ("gas", "density"), + ds, + center, + normal_vector, + width, + resolution, + ("gas", "density"), pixelmeaning=pixelmeaning, ) find_compare_maxima(expected_maxima, buf1, resolution, width) From 30b8f11e0d555117074c6b305c8997b0cec8303b Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Tue, 25 Feb 2025 19:16:18 +0100 Subject: [PATCH 26/68] fix function annotation error --- yt/testing.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/yt/testing.py b/yt/testing.py index 5772f9e1818..d088f0710d5 100644 --- a/yt/testing.py +++ b/yt/testing.py @@ -157,8 +157,8 @@ def integrate_kernel( def pixelintegrate_kernel( kernelfunc: Callable[[float], float], - pixelxy: np.ndarray[float], - pcenxy: np.ndarray[float], + pixelxy: np.ndarray, + pcenxy: np.ndarray, hsml: float, nsample: int = 100, periodxy: tuple[float | None, float | None] = (True, True), From 42bf76664bc6a53f6dd2619e53851298091856a2 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Tue, 25 Feb 2025 19:45:21 +0100 Subject: [PATCH 27/68] cython function wants bounds to be an array --- yt/visualization/volume_rendering/off_axis_projection.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yt/visualization/volume_rendering/off_axis_projection.py b/yt/visualization/volume_rendering/off_axis_projection.py index e760ba739eb..e7cc7551142 100644 --- a/yt/visualization/volume_rendering/off_axis_projection.py +++ b/yt/visualization/volume_rendering/off_axis_projection.py @@ -246,7 +246,7 @@ def off_axis_projection( re = data_source.ds.domain_right_edge.to("code_length").d x_min, y_min, z_min = le x_max, y_max, z_max = re - bounds = [x_min, x_max, y_min, y_max, z_min, z_max] + bounds = np.array([x_min, x_max, y_min, y_max, z_min, z_max]) # only need (rotated) x/y widths _width = (width.to("code_length").d)[:2] finfo = data_source.ds.field_info[item] From 2a3e0e1dc9f5876ae0abdb3d1f07cf39efc405ed Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Tue, 25 Feb 2025 20:09:10 +0100 Subject: [PATCH 28/68] pytest parametrize used -> ignore with nose --- nose_ignores.txt | 1 + tests/tests.yaml | 1 + 2 files changed, 2 insertions(+) diff --git a/nose_ignores.txt b/nose_ignores.txt index 20013886fba..075699d81e7 100644 --- a/nose_ignores.txt +++ b/nose_ignores.txt @@ -48,5 +48,6 @@ --ignore-file=test_disks\.py --ignore-file=test_offaxisprojection_pytestonly\.py --ignore-file=test_sph_pixelization_pytestonly\.py +--ignore-file=test_off_axis_SPH\.py --ignore-file=test_time_series\.py --ignore-file=test_cf_radial_pytest\.py diff --git a/tests/tests.yaml b/tests/tests.yaml index ae5515bca98..0e92d3f7408 100644 --- a/tests/tests.yaml +++ b/tests/tests.yaml @@ -231,5 +231,6 @@ other_tests: - "--ignore-file=test_offaxisprojection_pytestonly\\.py" - "--ignore-file=test_sph_pixelization_pytestonly\\.py" - "--ignore-file=test_cf_radial_pytest\\.py" + - "--ignore-file=test_off_axis_SPH\\.py" cookbook: - 'doc/source/cookbook/tests/test_cookbook.py' From df40e05ae72f68711ceb2c718015df962ba39104 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Tue, 25 Feb 2025 20:16:41 +0100 Subject: [PATCH 29/68] just be less strict on pixelave bounds input type --- 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 31fb9339d4c..b8dd73344f0 100644 --- a/yt/utilities/lib/pixelization_routines.pyx +++ b/yt/utilities/lib/pixelization_routines.pyx @@ -1360,7 +1360,7 @@ def pixelize_sph_kernel_projection_pixelave( any_float[:] pmass, any_float[:] pdens, any_float[:] quantity_to_smooth, - any_float[:] bounds, + bounds, kernel_name="cubic", weight_field=None, _check_period = (1, 1, 1), From 105018152a2738b7c74f107afe1c422863bcac20 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Wed, 26 Feb 2025 16:14:09 +0100 Subject: [PATCH 30/68] free shared kern2by2 *after* parallel section --- yt/utilities/lib/pixelization_routines.pyx | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx index b8dd73344f0..0f05e1ece22 100644 --- a/yt/utilities/lib/pixelization_routines.pyx +++ b/yt/utilities/lib/pixelization_routines.pyx @@ -1713,15 +1713,17 @@ def pixelize_sph_kernel_projection_pixelave( for sxi in range(xsize): for syi in range(ysize): buff[sxi, syi] += local_buff[sxi + syi * xsize] + # free memory in (hopefully) private variables assigned in + # each thread free(local_buff) - free(kern2by2) free(xiterv) free(yiterv) free(ziterv) free(xiter) free(yiter) free(ziter) - + # free memory in shared variable + free(kern2by2) return mask @cython.boundscheck(False) From 1642c3aaf930e1a2eea29a390487ca2542fb14cf Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Wed, 26 Feb 2025 18:10:22 +0100 Subject: [PATCH 31/68] debugged until test runs to assert --- .../tests/test_sph_pixelization_pytestonly.py | 31 ++++++++++--------- yt/testing.py | 17 ++++++---- 2 files changed, 28 insertions(+), 20 deletions(-) diff --git a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py index 732438b7341..bec58ff257c 100644 --- a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py +++ b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py @@ -161,12 +161,11 @@ def makemasses(i, j, k): # print(img.v) assert_rel_equal(expected_out, img.v, 5) - -@pytest.mark.parametrize("weighted", [True, False]) +@pytest.mark.parametrize("axis", [0, 1, 2]) +@pytest.mark.parametrize("shiftcenter", [False, True]) @pytest.mark.parametrize("periodic", [True, False]) @pytest.mark.parametrize("depth", [None, (1.0, "cm"), (5.0, "cm")]) -@pytest.mark.parametrize("shiftcenter", [False, True]) -@pytest.mark.parametrize("axis", [0, 1, 2]) +@pytest.mark.parametrize("weighted", [True, False]) @pytest.mark.parametrize("hsmlfac", [0.2, 1.0, 2.0]) def test_sph_projection_pixelave_alongaxes( axis: int, @@ -252,7 +251,8 @@ def makemasses(i, j, k): pixedges_y = np.linspace( center[ax1] - 0.5 * projwidth, center[ax1] + 0.5 * projwidth, buff_size[1] + 1 ) - _depth = depth if depth is not None else np.inf + _depth = (depth if depth is not None + else unyt.unyt_array(np.array((np.inf,)), "cm")) zlims = (center[axis] - 0.5 * _depth, center[axis] + 0.5 * _depth) ad = ds.all_data() pcenxy = np.array( @@ -267,18 +267,20 @@ def makemasses(i, j, k): density = (ad[("gas", "density")]).to("g * cm**-3") baseline = np.zeros(buff_size) if periodic: - periodxy = (bbox[ax0][1] - bbox[ax0][0], bbox[ax1][1] - bbox[ax1][0]) + periodxy = ((bbox[ax0][1] - bbox[ax0][0]).to("cm").v, + (bbox[ax1][1] - bbox[ax1][0]).to("cm").v) periodz = bbox[axis][1] - bbox[axis][0] else: periodxy = (None, None) for i in range(buff_size[0]): for j in range(buff_size[1]): - pixelxy = np.array( - [pixedges_x[i], pixedges_x[i + 1], pixedges_y[j], pixedges_y[j + 1]] - ) + pixelxy = unyt.unyt_array(np.array( + [pixedges_x[i], pixedges_x[i + 1], + pixedges_y[j], pixedges_y[j + 1]] + ), "cm") weightsum = 0.0 weightedsum = 0.0 - for p in range(len(pcenxy.shape[0])): + for p in range(pcenxy.shape[0]): # check if the particle is excluded along z if projz[p] < zlims[0] or projz[p] > zlims[1]: if periodz is None: @@ -288,15 +290,15 @@ def makemasses(i, j, k): dz += 0.5 * periodz dz %= periodz dz -= periodz - if np.abs(dz) > 0.5 * depth: + if np.abs(dz) > 0.5 * _depth: continue kernint = pixelintegrate_kernel( cubicspline_python, - pixelxy, + pixelxy.to("cm").v, pcenxy[p], - hsml[p], - nsample=100, + (hsml[p].to("cm")).v, + nsample=50, periodxy=periodxy, ) weightsum += kernint * mass[p] @@ -306,6 +308,7 @@ def makemasses(i, j, k): baseline[i, j] += weightedsum / weightsum else: baseline[i, j] += weightsum + baseline[np.isnan(baseline)] = 0. assert_rel_equal(baseline, img.v, 5) diff --git a/yt/testing.py b/yt/testing.py index d088f0710d5..b13cbcec744 100644 --- a/yt/testing.py +++ b/yt/testing.py @@ -124,6 +124,7 @@ def integrate_kernel( ], b: float | npt.NDArray[np.floating], hsml: float | npt.NDArray[np.floating], + nsample: int = 500, ) -> npt.NDArray[np.floating]: """ integrates a kernel function over a line passing entirely @@ -137,6 +138,8 @@ def integrate_kernel( impact parameter hsml: smoothing length [same units as impact parameter] + nsample: + number of points to sample for the integral Returns: -------- @@ -145,9 +148,10 @@ def integrate_kernel( """ pre = 1.0 / hsml**2 x = b / hsml + x[x >= 1.] = 1. # kernel is zero at 1. and > 1. xmax = np.sqrt(1.0 - x**2) xmin = -1.0 * xmax - xe = np.linspace(xmin, xmax, 500) # shape: 500, x.shape + xe = np.linspace(xmin, xmax, nsample) # shape: 500, x.shape xc = 0.5 * (xe[:-1, ...] + xe[1:, ...]) dx = np.diff(xe, axis=0) spv = kernelfunc(np.sqrt(xc**2 + x**2)) @@ -160,7 +164,7 @@ def pixelintegrate_kernel( pixelxy: np.ndarray, pcenxy: np.ndarray, hsml: float, - nsample: int = 100, + nsample: int = 50, periodxy: tuple[float | None, float | None] = (True, True), ) -> float: """ @@ -201,10 +205,10 @@ def pixelintegrate_kernel( """ xedges = np.linspace(pixelxy[0], pixelxy[1], nsample + 1) xcens = 0.5 * (xedges[:-1] + xedges[1:]) - dx = np.diff(xedges) + dx = np.diff(xedges, axis=0) yedges = np.linspace(pixelxy[2], pixelxy[3], nsample + 1) ycens = 0.5 * (yedges[:-1] + yedges[1:]) - dy = np.diff(yedges) + dy = np.diff(yedges, axis=0) xoff = xcens - pcenxy[0] if periodxy[0] is not None: @@ -216,8 +220,9 @@ def pixelintegrate_kernel( yoff += 0.5 * periodxy[1] yoff %= periodxy[1] yoff -= 0.5 * periodxy[1] - bpars = np.sqrt(xoff[:, np.newaxis] ** 2 + yoff[np.newaxis, :] ** 2) - lineints = integrate_kernel(kernelfunc, bpars, hsml) + bpars = np.sqrt(xoff[:, np.newaxis, ...] ** 2 + + yoff[np.newaxis, :, ...] ** 2) + lineints = integrate_kernel(kernelfunc, bpars, hsml, nsample=nsample) volint = np.sum(lineints * dx[:, np.newaxis] * dy[np.newaxis, :]) return volint / (pixelxy[1] - pixelxy[0]) / (pixelxy[3] - pixelxy[2]) From 429874a4365713644ed3b9c222d896122841688f Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 26 Feb 2025 17:17:06 +0000 Subject: [PATCH 32/68] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../tests/test_sph_pixelization_pytestonly.py | 22 +++++++++++-------- yt/testing.py | 7 +++--- 2 files changed, 16 insertions(+), 13 deletions(-) diff --git a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py index bec58ff257c..5a92f8b93b3 100644 --- a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py +++ b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py @@ -161,6 +161,7 @@ def makemasses(i, j, k): # print(img.v) assert_rel_equal(expected_out, img.v, 5) + @pytest.mark.parametrize("axis", [0, 1, 2]) @pytest.mark.parametrize("shiftcenter", [False, True]) @pytest.mark.parametrize("periodic", [True, False]) @@ -251,8 +252,7 @@ def makemasses(i, j, k): pixedges_y = np.linspace( center[ax1] - 0.5 * projwidth, center[ax1] + 0.5 * projwidth, buff_size[1] + 1 ) - _depth = (depth if depth is not None - else unyt.unyt_array(np.array((np.inf,)), "cm")) + _depth = depth if depth is not None else unyt.unyt_array(np.array((np.inf,)), "cm") zlims = (center[axis] - 0.5 * _depth, center[axis] + 0.5 * _depth) ad = ds.all_data() pcenxy = np.array( @@ -267,17 +267,21 @@ def makemasses(i, j, k): density = (ad[("gas", "density")]).to("g * cm**-3") baseline = np.zeros(buff_size) if periodic: - periodxy = ((bbox[ax0][1] - bbox[ax0][0]).to("cm").v, - (bbox[ax1][1] - bbox[ax1][0]).to("cm").v) + periodxy = ( + (bbox[ax0][1] - bbox[ax0][0]).to("cm").v, + (bbox[ax1][1] - bbox[ax1][0]).to("cm").v, + ) periodz = bbox[axis][1] - bbox[axis][0] else: periodxy = (None, None) for i in range(buff_size[0]): for j in range(buff_size[1]): - pixelxy = unyt.unyt_array(np.array( - [pixedges_x[i], pixedges_x[i + 1], - pixedges_y[j], pixedges_y[j + 1]] - ), "cm") + pixelxy = unyt.unyt_array( + np.array( + [pixedges_x[i], pixedges_x[i + 1], pixedges_y[j], pixedges_y[j + 1]] + ), + "cm", + ) weightsum = 0.0 weightedsum = 0.0 for p in range(pcenxy.shape[0]): @@ -308,7 +312,7 @@ def makemasses(i, j, k): baseline[i, j] += weightedsum / weightsum else: baseline[i, j] += weightsum - baseline[np.isnan(baseline)] = 0. + baseline[np.isnan(baseline)] = 0.0 assert_rel_equal(baseline, img.v, 5) diff --git a/yt/testing.py b/yt/testing.py index b13cbcec744..0cff1fb646b 100644 --- a/yt/testing.py +++ b/yt/testing.py @@ -124,7 +124,7 @@ def integrate_kernel( ], b: float | npt.NDArray[np.floating], hsml: float | npt.NDArray[np.floating], - nsample: int = 500, + nsample: int = 500, ) -> npt.NDArray[np.floating]: """ integrates a kernel function over a line passing entirely @@ -148,7 +148,7 @@ def integrate_kernel( """ pre = 1.0 / hsml**2 x = b / hsml - x[x >= 1.] = 1. # kernel is zero at 1. and > 1. + x[x >= 1.0] = 1.0 # kernel is zero at 1. and > 1. xmax = np.sqrt(1.0 - x**2) xmin = -1.0 * xmax xe = np.linspace(xmin, xmax, nsample) # shape: 500, x.shape @@ -220,8 +220,7 @@ def pixelintegrate_kernel( yoff += 0.5 * periodxy[1] yoff %= periodxy[1] yoff -= 0.5 * periodxy[1] - bpars = np.sqrt(xoff[:, np.newaxis, ...] ** 2 - + yoff[np.newaxis, :, ...] ** 2) + bpars = np.sqrt(xoff[:, np.newaxis, ...] ** 2 + yoff[np.newaxis, :, ...] ** 2) lineints = integrate_kernel(kernelfunc, bpars, hsml, nsample=nsample) volint = np.sum(lineints * dx[:, np.newaxis] * dy[np.newaxis, :]) return volint / (pixelxy[1] - pixelxy[0]) / (pixelxy[3] - pixelxy[2]) From 01371d1c73192fe0d89166390f292b1efdae47cc Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Fri, 28 Feb 2025 11:20:37 +0100 Subject: [PATCH 33/68] fix errors in pencilbeam tests (types in testing funcs) --- yt/testing.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/yt/testing.py b/yt/testing.py index 0cff1fb646b..afbbf6cf531 100644 --- a/yt/testing.py +++ b/yt/testing.py @@ -147,7 +147,7 @@ def integrate_kernel( units: 1 / (units of b and hsml)**2 """ pre = 1.0 / hsml**2 - x = b / hsml + x = np.asarray(b / hsml) x[x >= 1.0] = 1.0 # kernel is zero at 1. and > 1. xmax = np.sqrt(1.0 - x**2) xmin = -1.0 * xmax @@ -160,7 +160,7 @@ def integrate_kernel( def pixelintegrate_kernel( - kernelfunc: Callable[[float], float], + kernelfunc: Callable[[float | np.ndarray], np.ndarray], pixelxy: np.ndarray, pcenxy: np.ndarray, hsml: float, From ad6f3774c30521577719cf881340c578296ba990 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Fri, 28 Feb 2025 11:22:20 +0100 Subject: [PATCH 34/68] add debugging print lines --- .../tests/test_sph_pixelization_pytestonly.py | 16 +++++++-- yt/utilities/lib/pixelization_routines.pyx | 34 +++++++++++++++---- 2 files changed, 40 insertions(+), 10 deletions(-) diff --git a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py index 5a92f8b93b3..17ebf7a6bc0 100644 --- a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py +++ b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py @@ -156,9 +156,12 @@ def makemasses(i, j, k): if (not periodic) and shiftcenter: expected_out[:1, :] = 0.0 expected_out[:, :1] = 0.0 - # print(axis, shiftcenter, depth, periodic, weighted) - # print(expected_out) - # print(img.v) + #print(f"axis: {axis}, shiftcenter: {shiftcenter}, " + # f"depth: {depth}, periodic: {periodic}, weighted: {weighted}") + #print("expected:") + #print(expected_out) + #print("got:") + #print(img.v) assert_rel_equal(expected_out, img.v, 5) @@ -313,6 +316,13 @@ def makemasses(i, j, k): else: baseline[i, j] += weightsum baseline[np.isnan(baseline)] = 0.0 + print(f"axis: {axis}, shiftcenter: {shiftcenter}, " + f"depth: {depth}, periodic: {periodic}, weighted: {weighted}, " + f"hsmlfac: {hsmlfac}") + print("expected:") + print(baseline) + print("got:") + print(img.v) assert_rel_equal(baseline, img.v, 5) diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx index 0f05e1ece22..83869594fb7 100644 --- a/yt/utilities/lib/pixelization_routines.pyx +++ b/yt/utilities/lib/pixelization_routines.pyx @@ -1385,12 +1385,14 @@ def pixelize_sph_kernel_projection_pixelave( cdef np.int64_t xi, yi, x0, x1, y0, y1, xxi, yyi cdef np.int64_t nx, ny, nsubx, nsuby cdef int i, j, ii, jj, kk, ci, xsubi, ysubi - cdef int * xiter - cdef int * yiter - cdef int * ziter - cdef np.float64_t * xiterv - cdef np.float64_t * yiterv - cdef np.float64_t * ziterv + #cdef int * xiter + #cdef int * yiter + #cdef int * ziter + #cdef np.float64_t * xiterv + #cdef np.float64_t * yiterv + #cdef np.float64_t * ziterv + + print("Line 1395 pixelization_routines.pyx : function call, cdef block OK") if weight_field is not None: _weight_field = weight_field @@ -1423,6 +1425,7 @@ def pixelize_sph_kernel_projection_pixelave( if kernel_name not in kernel_tables: kernel_tables[kernel_name] = SPHKernelInterpolationTable(kernel_name) cdef SPHKernelInterpolationTable itab = kernel_tables[kernel_name] + print("Line 1428 pixelization_routines.pyx : basic setup seems to work") # pre-calculate kernel integral values to use for small particles # overlapping max. 2x2 output grid pixels. dimensions: @@ -1440,6 +1443,7 @@ def pixelize_sph_kernel_projection_pixelave( # step 1: just get a big grid of the 1D integral values kern1dint = malloc(sizeof(np.float64_t) * nsmin * nsmin) insmin = 1.0 / nsmin + print("Line 1446 pixelization_routines.pyx : kern2by2, kern1dint allocated") for si in range(nsmin): xnorm = -1.0 + 2.0 * (si + 0.5) * insmin xnorm = xnorm * xnorm @@ -1450,6 +1454,7 @@ def pixelize_sph_kernel_projection_pixelave( kern1dint[nsmin * si + sj] = 0.0 else: kern1dint[nsmin * si + sj] = itab.interpolate(q2) + print("Line 1446 pixelization_routines.pyx : kern1dint calculated") # step 2: integrate the 1D values for each grid edge position for si in range(nsmin + 1): for sj in range(nsmin + 1): @@ -1470,7 +1475,9 @@ def pixelize_sph_kernel_projection_pixelave( kern2by2[sci + 2] += kern1dint[nsmin * sii + sjj] for sjj in range(sj, nsmin): kern2by2[sci + 3] += kern1dint[nsmin * sii + sjj] + print("Line 1478 pixelization_routines.pyx : kern2by2 calculated") free(kern1dint) + print("Line 1480 pixelization_routines.pyx : kern1dint freed") with nogil, parallel(): # loop through every particle @@ -1498,7 +1505,8 @@ def pixelize_sph_kernel_projection_pixelave( # operations to the output array weren't much slower than # using the memory-intensive ouput array buffer for each # thread) - + with gil: + print("Line 1509 pixelization_routines.pyx : parallel started") local_buff = malloc(sizeof(np.float64_t) * xsize * ysize) xiterv = malloc(sizeof(np.float64_t) * 2) @@ -1511,11 +1519,16 @@ def pixelize_sph_kernel_projection_pixelave( xiterv[0] = yiterv[0] = ziterv[0] = 0.0 for i in range(xsize * ysize): local_buff[i] = 0.0 + with gil: + print("Line 1522 pixelization_routines.pyx : thread vars set up") for j in prange(0, posx.shape[0], schedule="dynamic"): if j % 100000 == 0: with gil: PyErr_CheckSignals() + with gil: + print("Line 1530 pixelization_routines.pyx :" + f"loop over particles started; at index {j}") xiter[1] = yiter[1] = ziter[1] = 999 @@ -1710,11 +1723,15 @@ def pixelize_sph_kernel_projection_pixelave( mask[xi, yi] = 1 with gil: + print("Line 1726 pixelization_routines.pyx :" + f"loop over particles done in this thread") for sxi in range(xsize): for syi in range(ysize): buff[sxi, syi] += local_buff[sxi + syi * xsize] # free memory in (hopefully) private variables assigned in # each thread + print("Line 1733 pixelization_routines.pyx :" + f"thread buffer added to output grid") free(local_buff) free(xiterv) free(yiterv) @@ -1722,6 +1739,9 @@ def pixelize_sph_kernel_projection_pixelave( free(xiter) free(yiter) free(ziter) + with gil: + print("Line 1743 pixelization_routines.pyx :" + f"thread variables freed") # free memory in shared variable free(kern2by2) return mask From b0201fa78114d843af9fa77885d4b15eef126602 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Fri, 28 Feb 2025 11:24:36 +0100 Subject: [PATCH 35/68] shorter func name --- .../coordinates/tests/test_sph_pixelization_pytestonly.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py index 17ebf7a6bc0..62566fe6ec5 100644 --- a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py +++ b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py @@ -171,7 +171,7 @@ def makemasses(i, j, k): @pytest.mark.parametrize("depth", [None, (1.0, "cm"), (5.0, "cm")]) @pytest.mark.parametrize("weighted", [True, False]) @pytest.mark.parametrize("hsmlfac", [0.2, 1.0, 2.0]) -def test_sph_projection_pixelave_alongaxes( +def test_sph_proj_pixelave_alongaxes( axis: int, shiftcenter: bool, periodic: bool, From ba9e7071a6a3588b2934491066e8df1601011c7e Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Fri, 28 Feb 2025 11:38:38 +0100 Subject: [PATCH 36/68] use pixel-meaning-dependent projection function --- yt/geometry/coordinates/cartesian_coordinates.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yt/geometry/coordinates/cartesian_coordinates.py b/yt/geometry/coordinates/cartesian_coordinates.py index 1611078edd7..2b750d87f5c 100644 --- a/yt/geometry/coordinates/cartesian_coordinates.py +++ b/yt/geometry/coordinates/cartesian_coordinates.py @@ -466,7 +466,7 @@ def _ortho_pixelize( if weight is None: for chunk in proj_reg.chunks([], "io"): data_source._initialize_projected_units([field], chunk) - pixelize_sph_kernel_projection_pencilbeam( + pixelize_sph_kernel_projection( buff, mask_uint8, chunk[ptype, px_name].to("code_length"), From f1ed43c47e72fda398c6bc0a97ccbdc20088cbaf Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Fri, 28 Feb 2025 13:23:58 +0100 Subject: [PATCH 37/68] lots of debug msgs, some actually debugging --- .../tests/test_sph_pixelization_pytestonly.py | 28 +++--- yt/utilities/lib/pixelization_routines.pyx | 89 +++++++++++++------ 2 files changed, 79 insertions(+), 38 deletions(-) diff --git a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py index 62566fe6ec5..7614c430954 100644 --- a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py +++ b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py @@ -165,12 +165,12 @@ def makemasses(i, j, k): assert_rel_equal(expected_out, img.v, 5) -@pytest.mark.parametrize("axis", [0, 1, 2]) -@pytest.mark.parametrize("shiftcenter", [False, True]) -@pytest.mark.parametrize("periodic", [True, False]) -@pytest.mark.parametrize("depth", [None, (1.0, "cm"), (5.0, "cm")]) -@pytest.mark.parametrize("weighted", [True, False]) -@pytest.mark.parametrize("hsmlfac", [0.2, 1.0, 2.0]) +@pytest.mark.parametrize("axis", [0]) # , 1, 2 +@pytest.mark.parametrize("shiftcenter", [False]) #, True +@pytest.mark.parametrize("periodic", [True]) # , False +@pytest.mark.parametrize("depth", [None]) # , (1.0, "cm"), (5.0, "cm") +@pytest.mark.parametrize("weighted", [False]) # , True +@pytest.mark.parametrize("hsmlfac", [0.2]) # , 1.0, 2.0 def test_sph_proj_pixelave_alongaxes( axis: int, shiftcenter: bool, @@ -180,7 +180,7 @@ def test_sph_proj_pixelave_alongaxes( hsmlfac: float, ) -> None: if shiftcenter: - center = unyt.unyt_array(np.array((0.625, 0.525, 0.425)), "cm") + center = unyt.unyt_array(np.array((0.6, 0.5, 0.4)), "cm") else: center = unyt.unyt_array(np.array((1.5, 1.5, 1.5)), "cm") bbox = unyt.unyt_array(np.array([[0.0, 3.0], [0.0, 3.0], [0.0, 3.0]]), "cm") @@ -200,7 +200,7 @@ def test_sph_proj_pixelave_alongaxes( # test correct centering, particle selection def makemasses(i, j, k): - if i == j == k == 1: + if i == 0 and j == 1 and k == 2: return 2.0 else: return 1.0 @@ -247,6 +247,8 @@ def makemasses(i, j, k): pixelmeaning="pixelave", ) img = prj.frb.data[("gas", "density")] + print("FRB info:") + print(prj.frb._get_info(("gas", "density"))) projwidth = unyt.unyt_array(np.array((2.5,)), "cm") pixedges_x = np.linspace( @@ -308,13 +310,15 @@ def makemasses(i, j, k): nsample=50, periodxy=periodxy, ) - weightsum += kernint * mass[p] + weightsum += kernint * mass[p].to("g").v if weighted: - weightedsum += kernint * mass[p] * density[p] + weightedsum += (kernint * mass[p].to("g").v + * density[p].to("g * cm**-3").v) if weighted: - baseline[i, j] += weightedsum / weightsum + baseline[i, j] += weightedsum[()] / weightsum[()] else: - baseline[i, j] += weightsum + baseline[i, j] += weightsum[()] + baseline = baseline.T # projections use (y, x) image axes order baseline[np.isnan(baseline)] = 0.0 print(f"axis: {axis}, shiftcenter: {shiftcenter}, " f"depth: {depth}, periodic: {periodic}, weighted: {weighted}, " diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx index 83869594fb7..2692028e900 100644 --- a/yt/utilities/lib/pixelization_routines.pyx +++ b/yt/utilities/lib/pixelization_routines.pyx @@ -1366,6 +1366,12 @@ def pixelize_sph_kernel_projection_pixelave( _check_period = (1, 1, 1), period=None, nsample_min=32): + # axes x, y, and z here need not be the corresponding axes in + # the simulation. Projections are always along the "z" axis here, + # bounds are in this function's x, y, z order. The calling python + # function should ensure these line up with the line of sight + # (here z) and projection plane (x, y) coordinate orders in the + # (possible rotated) simulation coordinates. # shared variables (or used before parallel section) cdef np.intp_t xsize, ysize, si, sj, sii, sjj, sci @@ -1385,15 +1391,14 @@ def pixelize_sph_kernel_projection_pixelave( cdef np.int64_t xi, yi, x0, x1, y0, y1, xxi, yyi cdef np.int64_t nx, ny, nsubx, nsuby cdef int i, j, ii, jj, kk, ci, xsubi, ysubi - #cdef int * xiter - #cdef int * yiter - #cdef int * ziter - #cdef np.float64_t * xiterv - #cdef np.float64_t * yiterv - #cdef np.float64_t * ziterv - + cdef int * xiter + cdef int * yiter + cdef int * ziter + cdef np.float64_t * xiterv + cdef np.float64_t * yiterv + cdef np.float64_t * ziterv print("Line 1395 pixelization_routines.pyx : function call, cdef block OK") - + print(f"check_period: {check_period}") if weight_field is not None: _weight_field = weight_field weighted = 1 @@ -1402,8 +1407,13 @@ def pixelize_sph_kernel_projection_pixelave( period_x = period[0] period_y = period[1] period_z = period[2] + print(f"period_x: {period_x}, period_y: {period_y}, period_z: {period_z}") + else: + print("period input was None") for si in range(3): check_period[si] = _check_period[si] + print(f"check_period: ({check_period[0]}, {check_period[1]}, {check_period[2]})") + # we find the x and y range over which we have pixels and we find how many # pixels we have in each dimension xsize, ysize = buff.shape[0], buff.shape[1] @@ -1443,7 +1453,7 @@ def pixelize_sph_kernel_projection_pixelave( # step 1: just get a big grid of the 1D integral values kern1dint = malloc(sizeof(np.float64_t) * nsmin * nsmin) insmin = 1.0 / nsmin - print("Line 1446 pixelization_routines.pyx : kern2by2, kern1dint allocated") + print("pixelization_routines.pyx : kern2by2, kern1dint allocated") for si in range(nsmin): xnorm = -1.0 + 2.0 * (si + 0.5) * insmin xnorm = xnorm * xnorm @@ -1454,7 +1464,7 @@ def pixelize_sph_kernel_projection_pixelave( kern1dint[nsmin * si + sj] = 0.0 else: kern1dint[nsmin * si + sj] = itab.interpolate(q2) - print("Line 1446 pixelization_routines.pyx : kern1dint calculated") + print("pixelization_routines.pyx : kern1dint calculated") # step 2: integrate the 1D values for each grid edge position for si in range(nsmin + 1): for sj in range(nsmin + 1): @@ -1475,9 +1485,9 @@ def pixelize_sph_kernel_projection_pixelave( kern2by2[sci + 2] += kern1dint[nsmin * sii + sjj] for sjj in range(sj, nsmin): kern2by2[sci + 3] += kern1dint[nsmin * sii + sjj] - print("Line 1478 pixelization_routines.pyx : kern2by2 calculated") + print("pixelization_routines.pyx : kern2by2 calculated") free(kern1dint) - print("Line 1480 pixelization_routines.pyx : kern1dint freed") + print("pixelization_routines.pyx : kern1dint freed") with nogil, parallel(): # loop through every particle @@ -1506,7 +1516,7 @@ def pixelize_sph_kernel_projection_pixelave( # using the memory-intensive ouput array buffer for each # thread) with gil: - print("Line 1509 pixelization_routines.pyx : parallel started") + print("pixelization_routines.pyx : parallel started") local_buff = malloc(sizeof(np.float64_t) * xsize * ysize) xiterv = malloc(sizeof(np.float64_t) * 2) @@ -1520,14 +1530,14 @@ def pixelize_sph_kernel_projection_pixelave( for i in range(xsize * ysize): local_buff[i] = 0.0 with gil: - print("Line 1522 pixelization_routines.pyx : thread vars set up") + print("pixelization_routines.pyx : thread vars set up") for j in prange(0, posx.shape[0], schedule="dynamic"): if j % 100000 == 0: with gil: PyErr_CheckSignals() with gil: - print("Line 1530 pixelization_routines.pyx :" + print("pixelization_routines.pyx :" f"loop over particles started; at index {j}") xiter[1] = yiter[1] = ziter[1] = 999 @@ -1597,18 +1607,28 @@ def pixelize_sph_kernel_projection_pixelave( # overlaps with a small number of pixels (max. 2x2) if hsml[j] < 0.5 * dx and hsml[j] < 0.5 * dy: # lower left corner indices - x0 = ((px - hsml[j] - x_min) * idx) - y0 = ((py - hsml[j] - y_min) * idy) + # floor: need to round down negative vals + # (int cast -> neg. vals. closer to zero) + x0 = math.floor( + (px - hsml[j] - x_min) * idx) + y0 = math.floor( + (py - hsml[j] - y_min) * idy) # sph particle lies entirely in one pixel if (px + hsml[j] < x_min + (x0 + 1) * dx and py + hsml[j] < y_min + (y0 + 1) * dy): # check that we're not on the wrong # periodicity iteration - if (x0 > 0 and x0 < xsize and - y0 > 0 and y0 < ysize): + with gil: + print(f"part {j}: has 1-cell overlap" + f" ({x0}, {y0})") + if (x0 >= 0 and x0 < xsize and + y0 >= 0 and y0 < ysize): # surface density # = density * av. length # = density * volume / area + with gil: + print(f"part {j}:" + " 1-cell overlap add") local_buff[x0 + y0 * xsize] += ( pmass[j] / pdens[j] * ipixA * qts_j ) @@ -1632,10 +1652,11 @@ def pixelize_sph_kernel_projection_pixelave( yn = (y_min + (y0 + 1) * dy - py) / hsml[j] # round to the nearest kernel grid edge # (C float to int cast behaves like - # floor()) - xi = ( + # floor() for input > 0, we need to round + # down negative numbers for periodic cases) + xi = math.floor( 0.5 * (xn + 1.) * nsmin + 0.5) - yi = ( + yi = math.floor( 0.5 * (yn + 1.) * nsmin + 0.5) # weight for each line integral: # area of subsampling pixel in length units @@ -1645,6 +1666,9 @@ def pixelize_sph_kernel_projection_pixelave( # coarse (part. pos rel. to grid) index ci = (4 * nsmin * nsmin * xi + 4 * nsmin * yi) + with gil: + print(f"part {j}: has 2x2-cell overlap;" + f"lower-left corner ({x0}, {y0})") for xxi in range(2): if x0 + xxi < 0 or x0 + xxi > xsize: # catch on the next periodic @@ -1659,6 +1683,9 @@ def pixelize_sph_kernel_projection_pixelave( if kv == 0.: # don't set mask to 1 continue + with gil: + print(f"part {j}: 2x2-cell" + "overlap add") local_buff[x0 + xxi + (y0 + yyi) * xsize] += ( prefactor_j * kv @@ -1672,18 +1699,26 @@ def pixelize_sph_kernel_projection_pixelave( # particle contributes to # subsampling does not depend on whether a # particle overlaps a (periodic) edge - x0 = ((px - hsml[j] - x_min) * idx) - x1 = ((px + hsml[j] - x_min) * idx) + x0 = math.floor( + (px - hsml[j] - x_min) * idx) + x1 = math.floor( + (px + hsml[j] - x_min) * idx) nx = x1 - x0 + 2 x0 = iclip(x0 - 1, 0, xsize) x1 = iclip(x1 + 1, 0, xsize) - y0 = ((py - hsml[j] - y_min) * idy) - y1 = ((py + hsml[j] - y_min) * idy) + y0 = math.floor( + (py - hsml[j] - y_min) * idy) + y1 = math.floor( + (py + hsml[j] - y_min) * idy) ny = y1 - y0 + 2 y0 = iclip(y0 - 1, 0, ysize) y1 = iclip(y1 + 1, 0, ysize) + with gil: + print(f"part {j}: has multi-cell overlap:" + f" {x0}--{x1}, {y0}--{y1}") + nsubx = ( ( nsmin / nx) @@ -1716,6 +1751,8 @@ def pixelize_sph_kernel_projection_pixelave( q_ij2 = ((posx_diff + posy_diff) * ih_j2) if q_ij2 >= 1.0: continue + with gil: + print(f"part {j}: multi-cell overlap add") local_buff[xi + yi*xsize] += ( prefactor_j * insuby * insubx * itab.interpolate(q_ij2) From ff0aa45da91ef5d6fb8586715941bc2927c4bb1b Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Fri, 28 Feb 2025 14:20:20 +0100 Subject: [PATCH 38/68] attempt to fix indexing in output nd kern2by2 --- yt/utilities/lib/pixelization_routines.pyx | 38 +++++++++++----------- 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx index 2692028e900..b8d91859dbd 100644 --- a/yt/utilities/lib/pixelization_routines.pyx +++ b/yt/utilities/lib/pixelization_routines.pyx @@ -1311,7 +1311,8 @@ def pixelize_sph_kernel_projection_pencilbeam( # see equation 32 of the SPLASH paper # now we just use the kernel projection - local_buff[xi + yi*xsize] += prefactor_j * itab.interpolate(q_ij2) + local_buff[xi + yi*xsize] += ( + prefactor_j * itab.interpolate(q_ij2)) mask[xi, yi] = 1 with gil: @@ -1461,30 +1462,30 @@ def pixelize_sph_kernel_projection_pixelave( ynorm = -1.0 + 2.0 * (sj + 0.5) * insmin q2 = xnorm + ynorm * ynorm if q2 >= 1.: - kern1dint[nsmin * si + sj] = 0.0 + kern1dint[si + nsmin * sj] = 0.0 else: - kern1dint[nsmin * si + sj] = itab.interpolate(q2) + kern1dint[si + nsmin * sj] = itab.interpolate(q2) print("pixelization_routines.pyx : kern1dint calculated") # step 2: integrate the 1D values for each grid edge position for si in range(nsmin + 1): for sj in range(nsmin + 1): # slowest 2 indices - sci = (4 * (nsmin + 1) * (nsmin + 1) * si - + 4 * (nsmin + 1) * sj) + sci = (4 * (nsmin + 1) * si + + 4 * (nsmin + 1) * (nsmin + 1) * sj) kern2by2[sci] = 0.0 kern2by2[sci + 1] = 0.0 kern2by2[sci + 2] = 0.0 kern2by2[sci + 3] = 0.0 for sii in range(0, si): for sjj in range(0, sj): - kern2by2[sci] += kern1dint[nsmin * sii + sjj] + kern2by2[sci] += kern1dint[sii + nsmin * sjj] for sjj in range(sj, nsmin): - kern2by2[sci + 1] += kern1dint[nsmin * sii + sjj] + kern2by2[sci + 1] += kern1dint[sii + nsmin * sjj] for sii in range(si, nsmin): for sjj in range(0, sj): - kern2by2[sci + 2] += kern1dint[nsmin * sii + sjj] + kern2by2[sci + 2] += kern1dint[sii + nsmin * sjj] for sjj in range(sj, nsmin): - kern2by2[sci + 3] += kern1dint[nsmin * sii + sjj] + kern2by2[sci + 3] += kern1dint[sii + nsmin * sjj] print("pixelization_routines.pyx : kern2by2 calculated") free(kern1dint) print("pixelization_routines.pyx : kern1dint freed") @@ -1616,11 +1617,11 @@ def pixelize_sph_kernel_projection_pixelave( # sph particle lies entirely in one pixel if (px + hsml[j] < x_min + (x0 + 1) * dx and py + hsml[j] < y_min + (y0 + 1) * dy): - # check that we're not on the wrong - # periodicity iteration with gil: print(f"part {j}: has 1-cell overlap" f" ({x0}, {y0})") + # check that we're not on the wrong + # periodicity iteration if (x0 >= 0 and x0 < xsize and y0 >= 0 and y0 < ysize): # surface density @@ -1629,7 +1630,7 @@ def pixelize_sph_kernel_projection_pixelave( with gil: print(f"part {j}:" " 1-cell overlap add") - local_buff[x0 + y0 * xsize] += ( + local_buff[x0 * ysize + y0] += ( pmass[j] / pdens[j] * ipixA * qts_j ) mask[x0, y0] = 1 @@ -1664,18 +1665,17 @@ def pixelize_sph_kernel_projection_pixelave( prefactor_j *= (4.0 * h_j2 * insmin * insmin * ipixA) # coarse (part. pos rel. to grid) index - ci = (4 * nsmin * nsmin * xi - + 4 * nsmin * yi) + ci = (4 * xi * (nsmin + 1) + 4 * yi) with gil: print(f"part {j}: has 2x2-cell overlap;" f"lower-left corner ({x0}, {y0})") for xxi in range(2): - if x0 + xxi < 0 or x0 + xxi > xsize: + if x0 + xxi < 0 or x0 + xxi >= xsize: # catch on the next periodic # loop continue for yyi in range(2): - if y0 + yyi < 0 or y0 + yyi > ysize: + if y0 + yyi < 0 or y0 + yyi >= ysize: # catch on the next # periodic loop continue @@ -1686,8 +1686,8 @@ def pixelize_sph_kernel_projection_pixelave( with gil: print(f"part {j}: 2x2-cell" "overlap add") - local_buff[x0 + xxi - + (y0 + yyi) * xsize] += ( + local_buff[(x0 + xxi) * ysize + + (y0 + yyi)] += ( prefactor_j * kv ) mask[xi, yi] = 1 @@ -1753,7 +1753,7 @@ def pixelize_sph_kernel_projection_pixelave( if q_ij2 >= 1.0: continue with gil: print(f"part {j}: multi-cell overlap add") - local_buff[xi + yi*xsize] += ( + local_buff[xi * ysize + yi] += ( prefactor_j * insuby * insubx * itab.interpolate(q_ij2) ) From 57a84764963705538d39d65286bbc7d5f504f02c Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Fri, 28 Feb 2025 14:20:45 +0100 Subject: [PATCH 39/68] matched due to indexing error? --- .../coordinates/tests/test_sph_pixelization_pytestonly.py | 1 - 1 file changed, 1 deletion(-) diff --git a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py index 7614c430954..ab3a5de15e6 100644 --- a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py +++ b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py @@ -318,7 +318,6 @@ def makemasses(i, j, k): baseline[i, j] += weightedsum[()] / weightsum[()] else: baseline[i, j] += weightsum[()] - baseline = baseline.T # projections use (y, x) image axes order baseline[np.isnan(baseline)] = 0.0 print(f"axis: {axis}, shiftcenter: {shiftcenter}, " f"depth: {depth}, periodic: {periodic}, weighted: {weighted}, " From 63708bd7c805bbe06b0f55fef1ddc94bbb0bd446 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Fri, 28 Feb 2025 14:47:42 +0100 Subject: [PATCH 40/68] print statement edits --- yt/utilities/lib/pixelization_routines.pyx | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx index b8d91859dbd..6b3afc7910a 100644 --- a/yt/utilities/lib/pixelization_routines.pyx +++ b/yt/utilities/lib/pixelization_routines.pyx @@ -1398,7 +1398,7 @@ def pixelize_sph_kernel_projection_pixelave( cdef np.float64_t * xiterv cdef np.float64_t * yiterv cdef np.float64_t * ziterv - print("Line 1395 pixelization_routines.pyx : function call, cdef block OK") + print("pixelization_routines.pyx : function call, cdef block OK") print(f"check_period: {check_period}") if weight_field is not None: _weight_field = weight_field @@ -1408,12 +1408,14 @@ def pixelize_sph_kernel_projection_pixelave( period_x = period[0] period_y = period[1] period_z = period[2] - print(f"period_x: {period_x}, period_y: {period_y}, period_z: {period_z}") + print(f"period_x: {period_x}, period_y: {period_y}," + f" period_z: {period_z}") else: print("period input was None") for si in range(3): check_period[si] = _check_period[si] - print(f"check_period: ({check_period[0]}, {check_period[1]}, {check_period[2]})") + print(f"check_period: ({check_period[0]}, {check_period[1]}," + f" {check_period[2]})") # we find the x and y range over which we have pixels and we find how many # pixels we have in each dimension @@ -1436,7 +1438,7 @@ def pixelize_sph_kernel_projection_pixelave( if kernel_name not in kernel_tables: kernel_tables[kernel_name] = SPHKernelInterpolationTable(kernel_name) cdef SPHKernelInterpolationTable itab = kernel_tables[kernel_name] - print("Line 1428 pixelization_routines.pyx : basic setup seems to work") + print("pixelization_routines.pyx : basic setup seems to work") # pre-calculate kernel integral values to use for small particles # overlapping max. 2x2 output grid pixels. dimensions: @@ -1761,14 +1763,14 @@ def pixelize_sph_kernel_projection_pixelave( with gil: print("Line 1726 pixelization_routines.pyx :" - f"loop over particles done in this thread") + "loop over particles done in this thread") for sxi in range(xsize): for syi in range(ysize): buff[sxi, syi] += local_buff[sxi + syi * xsize] # free memory in (hopefully) private variables assigned in # each thread - print("Line 1733 pixelization_routines.pyx :" - f"thread buffer added to output grid") + print("pixelization_routines.pyx :" + "thread buffer added to output grid") free(local_buff) free(xiterv) free(yiterv) @@ -1777,10 +1779,11 @@ def pixelize_sph_kernel_projection_pixelave( free(yiter) free(ziter) with gil: - print("Line 1743 pixelization_routines.pyx :" - f"thread variables freed") + print("pixelization_routines.pyx : thread variables freed") # free memory in shared variable + print("pixelization_routines.pyx : parallel part is over") free(kern2by2) + print("pixelization_routines.pyx : kern2by2 freed") return mask @cython.boundscheck(False) From cc0adef8f6b67d2c5bde79fc5b5eb84283459ff8 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Fri, 28 Feb 2025 15:41:49 +0100 Subject: [PATCH 41/68] a bit more test debugging --- .../tests/test_sph_pixelization_pytestonly.py | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py index ab3a5de15e6..fdbf2c1104b 100644 --- a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py +++ b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py @@ -165,12 +165,12 @@ def makemasses(i, j, k): assert_rel_equal(expected_out, img.v, 5) -@pytest.mark.parametrize("axis", [0]) # , 1, 2 -@pytest.mark.parametrize("shiftcenter", [False]) #, True -@pytest.mark.parametrize("periodic", [True]) # , False -@pytest.mark.parametrize("depth", [None]) # , (1.0, "cm"), (5.0, "cm") -@pytest.mark.parametrize("weighted", [False]) # , True -@pytest.mark.parametrize("hsmlfac", [0.2]) # , 1.0, 2.0 +@pytest.mark.parametrize("axis", [0, 1, 2]) +@pytest.mark.parametrize("shiftcenter", [False, True]) +@pytest.mark.parametrize("periodic", [True, False]) +@pytest.mark.parametrize("depth", [None, (1.0, "cm"), (5.0, "cm")]) +@pytest.mark.parametrize("weighted", [False, True ]) +@pytest.mark.parametrize("hsmlfac", [0.2, 0.5, 1.0]) def test_sph_proj_pixelave_alongaxes( axis: int, shiftcenter: bool, @@ -246,6 +246,7 @@ def makemasses(i, j, k): data_source=source, pixelmeaning="pixelave", ) + print("ProjectionPlot returned") img = prj.frb.data[("gas", "density")] print("FRB info:") print(prj.frb._get_info(("gas", "density"))) @@ -295,7 +296,7 @@ def makemasses(i, j, k): if periodz is None: continue else: - dz = projz - center[axis].to("cm") + dz = projz[p] - center[axis].to("cm") dz += 0.5 * periodz dz %= periodz dz -= periodz From 467c50974898389d180fc013757aa17660adf2cc Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Fri, 28 Feb 2025 16:17:09 +0100 Subject: [PATCH 42/68] style edit --- yt/geometry/coordinates/cartesian_coordinates.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/yt/geometry/coordinates/cartesian_coordinates.py b/yt/geometry/coordinates/cartesian_coordinates.py index 2b750d87f5c..ce82a34bdb3 100644 --- a/yt/geometry/coordinates/cartesian_coordinates.py +++ b/yt/geometry/coordinates/cartesian_coordinates.py @@ -430,7 +430,8 @@ def _ortho_pixelize( pixelize_sph_kernel_projection_pixelave ) else: - raise NotImplementedError(f"No pixelmeaning option {pixelmeaning}") + raise NotImplementedError( + f"No pixelmeaning option {pixelmeaning}") weight = data_source.weight_field moment = data_source.moment le, re = data_source.data_source.get_bbox() From 9635ae115dfa709f87313200b94f1041deff8534 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Fri, 28 Feb 2025 16:28:09 +0100 Subject: [PATCH 43/68] output grid axis order is (y, x) --- yt/utilities/lib/pixelization_routines.pyx | 43 +++++++++++----------- 1 file changed, 22 insertions(+), 21 deletions(-) diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx index 6b3afc7910a..0649efc2017 100644 --- a/yt/utilities/lib/pixelization_routines.pyx +++ b/yt/utilities/lib/pixelization_routines.pyx @@ -1392,12 +1392,12 @@ def pixelize_sph_kernel_projection_pixelave( cdef np.int64_t xi, yi, x0, x1, y0, y1, xxi, yyi cdef np.int64_t nx, ny, nsubx, nsuby cdef int i, j, ii, jj, kk, ci, xsubi, ysubi - cdef int * xiter - cdef int * yiter - cdef int * ziter - cdef np.float64_t * xiterv - cdef np.float64_t * yiterv - cdef np.float64_t * ziterv + #cdef int * xiter + #cdef int * yiter + #cdef int * ziter + #cdef np.float64_t * xiterv + #cdef np.float64_t * yiterv + #cdef np.float64_t * ziterv print("pixelization_routines.pyx : function call, cdef block OK") print(f"check_period: {check_period}") if weight_field is not None: @@ -1414,8 +1414,8 @@ def pixelize_sph_kernel_projection_pixelave( print("period input was None") for si in range(3): check_period[si] = _check_period[si] - print(f"check_period: ({check_period[0]}, {check_period[1]}," - f" {check_period[2]})") + print(f"check_period: ({check_period[0]}, {check_period[1]}," + f" {check_period[2]})") # we find the x and y range over which we have pixels and we find how many # pixels we have in each dimension @@ -1540,7 +1540,7 @@ def pixelize_sph_kernel_projection_pixelave( with gil: PyErr_CheckSignals() with gil: - print("pixelization_routines.pyx :" + print("pixelization_routines.pyx : " f"loop over particles started; at index {j}") xiter[1] = yiter[1] = ziter[1] = 999 @@ -1632,7 +1632,7 @@ def pixelize_sph_kernel_projection_pixelave( with gil: print(f"part {j}:" " 1-cell overlap add") - local_buff[x0 * ysize + y0] += ( + local_buff[x0 + y0 * xsize] += ( pmass[j] / pdens[j] * ipixA * qts_j ) mask[x0, y0] = 1 @@ -1670,7 +1670,7 @@ def pixelize_sph_kernel_projection_pixelave( ci = (4 * xi * (nsmin + 1) + 4 * yi) with gil: print(f"part {j}: has 2x2-cell overlap;" - f"lower-left corner ({x0}, {y0})") + f" lower-left corner ({x0}, {y0})") for xxi in range(2): if x0 + xxi < 0 or x0 + xxi >= xsize: # catch on the next periodic @@ -1687,9 +1687,9 @@ def pixelize_sph_kernel_projection_pixelave( continue with gil: print(f"part {j}: 2x2-cell" - "overlap add") - local_buff[(x0 + xxi) * ysize - + (y0 + yyi)] += ( + " overlap add") + local_buff[(x0 + xxi) + + (y0 + yyi) * xsize] += ( prefactor_j * kv ) mask[xi, yi] = 1 @@ -1754,8 +1754,9 @@ def pixelize_sph_kernel_projection_pixelave( * ih_j2) if q_ij2 >= 1.0: continue with gil: - print(f"part {j}: multi-cell overlap add") - local_buff[xi * ysize + yi] += ( + print(f"part {j}: multi-cell" + " overlap add") + local_buff[xi + yi * xsize] += ( prefactor_j * insuby * insubx * itab.interpolate(q_ij2) ) @@ -1766,11 +1767,11 @@ def pixelize_sph_kernel_projection_pixelave( "loop over particles done in this thread") for sxi in range(xsize): for syi in range(ysize): - buff[sxi, syi] += local_buff[sxi + syi * xsize] - # free memory in (hopefully) private variables assigned in - # each thread + buff[syi, sxi] += local_buff[sxi + syi * xsize] print("pixelization_routines.pyx :" "thread buffer added to output grid") + # free memory in (hopefully) private variables assigned in + # each thread free(local_buff) free(xiterv) free(yiterv) @@ -1778,8 +1779,8 @@ def pixelize_sph_kernel_projection_pixelave( free(xiter) free(yiter) free(ziter) - with gil: - print("pixelization_routines.pyx : thread variables freed") + #with gil: + print("pixelization_routines.pyx : thread variables freed") # free memory in shared variable print("pixelization_routines.pyx : parallel part is over") free(kern2by2) From 80c5bc6bd4d27d606644bc3542ebbe8d00838691 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Fri, 28 Feb 2025 16:36:59 +0100 Subject: [PATCH 44/68] fix domain-edge error non-periodic test --- .../coordinates/tests/test_sph_pixelization_pytestonly.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py index fdbf2c1104b..009389312fa 100644 --- a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py +++ b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py @@ -225,6 +225,9 @@ def makemasses(i, j, k): re = np.array(ds.domain_right_edge) le[axis] = center[axis] - 0.5 * depth re[axis] = center[axis] + 0.5 * depth + if not periodic: + le[axis] = max(le[axis], ds.domain_left_edge[axis]) + re[axis] = min(re[axis], ds.domain_right_edge[axis]) cen = 0.5 * (le + re) reg = YTRegion(center=cen, left_edge=le, right_edge=re, ds=ds) source = reg From 2f3f613d5e7767ba7e572ab6c208b16a6f89b316 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Fri, 28 Feb 2025 17:16:39 +0100 Subject: [PATCH 45/68] fix bug in test --- .../coordinates/tests/test_sph_pixelization_pytestonly.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py index 009389312fa..022ff66d6b8 100644 --- a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py +++ b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py @@ -169,7 +169,7 @@ def makemasses(i, j, k): @pytest.mark.parametrize("shiftcenter", [False, True]) @pytest.mark.parametrize("periodic", [True, False]) @pytest.mark.parametrize("depth", [None, (1.0, "cm"), (5.0, "cm")]) -@pytest.mark.parametrize("weighted", [False, True ]) +@pytest.mark.parametrize("weighted", [False, True]) @pytest.mark.parametrize("hsmlfac", [0.2, 0.5, 1.0]) def test_sph_proj_pixelave_alongaxes( axis: int, @@ -283,6 +283,7 @@ def makemasses(i, j, k): periodz = bbox[axis][1] - bbox[axis][0] else: periodxy = (None, None) + periodz = None for i in range(buff_size[0]): for j in range(buff_size[1]): pixelxy = unyt.unyt_array( From 6c4631684d535d06ed5e4f44f52d2fea78e848a1 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Fri, 28 Feb 2025 17:29:01 +0100 Subject: [PATCH 46/68] pixelave projection gets segfaults even if I free nothing --- 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 0649efc2017..53d91c9b4ae 100644 --- a/yt/utilities/lib/pixelization_routines.pyx +++ b/yt/utilities/lib/pixelization_routines.pyx @@ -1770,7 +1770,7 @@ def pixelize_sph_kernel_projection_pixelave( buff[syi, sxi] += local_buff[sxi + syi * xsize] print("pixelization_routines.pyx :" "thread buffer added to output grid") - # free memory in (hopefully) private variables assigned in + # free memory in (hopefully?) private variables assigned in # each thread free(local_buff) free(xiterv) From 76a7e69a3f1af8d50db1cd5036f30f70f72c333b Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Fri, 28 Feb 2025 18:11:08 +0100 Subject: [PATCH 47/68] more index revision --- yt/utilities/lib/pixelization_routines.pyx | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx index 53d91c9b4ae..2fb380e6e53 100644 --- a/yt/utilities/lib/pixelization_routines.pyx +++ b/yt/utilities/lib/pixelization_routines.pyx @@ -1464,30 +1464,29 @@ def pixelize_sph_kernel_projection_pixelave( ynorm = -1.0 + 2.0 * (sj + 0.5) * insmin q2 = xnorm + ynorm * ynorm if q2 >= 1.: - kern1dint[si + nsmin * sj] = 0.0 + kern1dint[si * nsmin + sj] = 0.0 else: - kern1dint[si + nsmin * sj] = itab.interpolate(q2) + kern1dint[si * nsmin + sj] = itab.interpolate(q2) print("pixelization_routines.pyx : kern1dint calculated") # step 2: integrate the 1D values for each grid edge position for si in range(nsmin + 1): for sj in range(nsmin + 1): # slowest 2 indices - sci = (4 * (nsmin + 1) * si - + 4 * (nsmin + 1) * (nsmin + 1) * sj) + sci = 4 * (nsmin + 1) * si + 4 * sj kern2by2[sci] = 0.0 kern2by2[sci + 1] = 0.0 kern2by2[sci + 2] = 0.0 kern2by2[sci + 3] = 0.0 for sii in range(0, si): for sjj in range(0, sj): - kern2by2[sci] += kern1dint[sii + nsmin * sjj] + kern2by2[sci] += kern1dint[sii * nsmin + sjj] for sjj in range(sj, nsmin): - kern2by2[sci + 1] += kern1dint[sii + nsmin * sjj] + kern2by2[sci + 1] += kern1dint[sii * nsmin + sjj] for sii in range(si, nsmin): for sjj in range(0, sj): - kern2by2[sci + 2] += kern1dint[sii + nsmin * sjj] + kern2by2[sci + 2] += kern1dint[sii * nsmin + sjj] for sjj in range(sj, nsmin): - kern2by2[sci + 3] += kern1dint[sii + nsmin * sjj] + kern2by2[sci + 3] += kern1dint[sii * nsmin + sjj] print("pixelization_routines.pyx : kern2by2 calculated") free(kern1dint) print("pixelization_routines.pyx : kern1dint freed") @@ -1692,7 +1691,7 @@ def pixelize_sph_kernel_projection_pixelave( + (y0 + yyi) * xsize] += ( prefactor_j * kv ) - mask[xi, yi] = 1 + mask[x0 + xxi, y0 + yyi] = 1 else: # subsample the kernel in each pixel From 85ade08a5587f8235ad2503bf0721264c14c1038 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Fri, 28 Feb 2025 18:26:15 +0100 Subject: [PATCH 48/68] tests fail but at least the segfaults are gone --- yt/utilities/lib/pixelization_routines.pyx | 47 ---------------------- 1 file changed, 47 deletions(-) diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx index 2fb380e6e53..12f7bc5dbbe 100644 --- a/yt/utilities/lib/pixelization_routines.pyx +++ b/yt/utilities/lib/pixelization_routines.pyx @@ -1398,8 +1398,6 @@ def pixelize_sph_kernel_projection_pixelave( #cdef np.float64_t * xiterv #cdef np.float64_t * yiterv #cdef np.float64_t * ziterv - print("pixelization_routines.pyx : function call, cdef block OK") - print(f"check_period: {check_period}") if weight_field is not None: _weight_field = weight_field weighted = 1 @@ -1408,14 +1406,8 @@ def pixelize_sph_kernel_projection_pixelave( period_x = period[0] period_y = period[1] period_z = period[2] - print(f"period_x: {period_x}, period_y: {period_y}," - f" period_z: {period_z}") - else: - print("period input was None") for si in range(3): check_period[si] = _check_period[si] - print(f"check_period: ({check_period[0]}, {check_period[1]}," - f" {check_period[2]})") # we find the x and y range over which we have pixels and we find how many # pixels we have in each dimension @@ -1438,7 +1430,6 @@ def pixelize_sph_kernel_projection_pixelave( if kernel_name not in kernel_tables: kernel_tables[kernel_name] = SPHKernelInterpolationTable(kernel_name) cdef SPHKernelInterpolationTable itab = kernel_tables[kernel_name] - print("pixelization_routines.pyx : basic setup seems to work") # pre-calculate kernel integral values to use for small particles # overlapping max. 2x2 output grid pixels. dimensions: @@ -1456,7 +1447,6 @@ def pixelize_sph_kernel_projection_pixelave( # step 1: just get a big grid of the 1D integral values kern1dint = malloc(sizeof(np.float64_t) * nsmin * nsmin) insmin = 1.0 / nsmin - print("pixelization_routines.pyx : kern2by2, kern1dint allocated") for si in range(nsmin): xnorm = -1.0 + 2.0 * (si + 0.5) * insmin xnorm = xnorm * xnorm @@ -1467,7 +1457,6 @@ def pixelize_sph_kernel_projection_pixelave( kern1dint[si * nsmin + sj] = 0.0 else: kern1dint[si * nsmin + sj] = itab.interpolate(q2) - print("pixelization_routines.pyx : kern1dint calculated") # step 2: integrate the 1D values for each grid edge position for si in range(nsmin + 1): for sj in range(nsmin + 1): @@ -1487,9 +1476,7 @@ def pixelize_sph_kernel_projection_pixelave( kern2by2[sci + 2] += kern1dint[sii * nsmin + sjj] for sjj in range(sj, nsmin): kern2by2[sci + 3] += kern1dint[sii * nsmin + sjj] - print("pixelization_routines.pyx : kern2by2 calculated") free(kern1dint) - print("pixelization_routines.pyx : kern1dint freed") with nogil, parallel(): # loop through every particle @@ -1517,8 +1504,6 @@ def pixelize_sph_kernel_projection_pixelave( # operations to the output array weren't much slower than # using the memory-intensive ouput array buffer for each # thread) - with gil: - print("pixelization_routines.pyx : parallel started") local_buff = malloc(sizeof(np.float64_t) * xsize * ysize) xiterv = malloc(sizeof(np.float64_t) * 2) @@ -1531,16 +1516,11 @@ def pixelize_sph_kernel_projection_pixelave( xiterv[0] = yiterv[0] = ziterv[0] = 0.0 for i in range(xsize * ysize): local_buff[i] = 0.0 - with gil: - print("pixelization_routines.pyx : thread vars set up") for j in prange(0, posx.shape[0], schedule="dynamic"): if j % 100000 == 0: with gil: PyErr_CheckSignals() - with gil: - print("pixelization_routines.pyx : " - f"loop over particles started; at index {j}") xiter[1] = yiter[1] = ziter[1] = 999 @@ -1618,9 +1598,6 @@ def pixelize_sph_kernel_projection_pixelave( # sph particle lies entirely in one pixel if (px + hsml[j] < x_min + (x0 + 1) * dx and py + hsml[j] < y_min + (y0 + 1) * dy): - with gil: - print(f"part {j}: has 1-cell overlap" - f" ({x0}, {y0})") # check that we're not on the wrong # periodicity iteration if (x0 >= 0 and x0 < xsize and @@ -1628,9 +1605,6 @@ def pixelize_sph_kernel_projection_pixelave( # surface density # = density * av. length # = density * volume / area - with gil: - print(f"part {j}:" - " 1-cell overlap add") local_buff[x0 + y0 * xsize] += ( pmass[j] / pdens[j] * ipixA * qts_j ) @@ -1667,9 +1641,6 @@ def pixelize_sph_kernel_projection_pixelave( * insmin * insmin * ipixA) # coarse (part. pos rel. to grid) index ci = (4 * xi * (nsmin + 1) + 4 * yi) - with gil: - print(f"part {j}: has 2x2-cell overlap;" - f" lower-left corner ({x0}, {y0})") for xxi in range(2): if x0 + xxi < 0 or x0 + xxi >= xsize: # catch on the next periodic @@ -1684,9 +1655,6 @@ def pixelize_sph_kernel_projection_pixelave( if kv == 0.: # don't set mask to 1 continue - with gil: - print(f"part {j}: 2x2-cell" - " overlap add") local_buff[(x0 + xxi) + (y0 + yyi) * xsize] += ( prefactor_j * kv @@ -1716,10 +1684,6 @@ def pixelize_sph_kernel_projection_pixelave( y0 = iclip(y0 - 1, 0, ysize) y1 = iclip(y1 + 1, 0, ysize) - with gil: - print(f"part {j}: has multi-cell overlap:" - f" {x0}--{x1}, {y0}--{y1}") - nsubx = ( ( nsmin / nx) @@ -1752,9 +1716,6 @@ def pixelize_sph_kernel_projection_pixelave( q_ij2 = ((posx_diff + posy_diff) * ih_j2) if q_ij2 >= 1.0: continue - with gil: - print(f"part {j}: multi-cell" - " overlap add") local_buff[xi + yi * xsize] += ( prefactor_j * insuby * insubx * itab.interpolate(q_ij2) @@ -1762,13 +1723,9 @@ def pixelize_sph_kernel_projection_pixelave( mask[xi, yi] = 1 with gil: - print("Line 1726 pixelization_routines.pyx :" - "loop over particles done in this thread") for sxi in range(xsize): for syi in range(ysize): buff[syi, sxi] += local_buff[sxi + syi * xsize] - print("pixelization_routines.pyx :" - "thread buffer added to output grid") # free memory in (hopefully?) private variables assigned in # each thread free(local_buff) @@ -1778,12 +1735,8 @@ def pixelize_sph_kernel_projection_pixelave( free(xiter) free(yiter) free(ziter) - #with gil: - print("pixelization_routines.pyx : thread variables freed") # free memory in shared variable - print("pixelization_routines.pyx : parallel part is over") free(kern2by2) - print("pixelization_routines.pyx : kern2by2 freed") return mask @cython.boundscheck(False) From d375ae8132fa81b486dc80af6389045235365254 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 28 Feb 2025 13:23:25 +0000 Subject: [PATCH 49/68] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../tests/test_sph_pixelization_pytestonly.py | 23 +++++++++++-------- yt/utilities/lib/pixelization_routines.pyx | 2 +- 2 files changed, 14 insertions(+), 11 deletions(-) diff --git a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py index 022ff66d6b8..95228237d6c 100644 --- a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py +++ b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py @@ -156,12 +156,12 @@ def makemasses(i, j, k): if (not periodic) and shiftcenter: expected_out[:1, :] = 0.0 expected_out[:, :1] = 0.0 - #print(f"axis: {axis}, shiftcenter: {shiftcenter}, " + # print(f"axis: {axis}, shiftcenter: {shiftcenter}, " # f"depth: {depth}, periodic: {periodic}, weighted: {weighted}") - #print("expected:") - #print(expected_out) - #print("got:") - #print(img.v) + # print("expected:") + # print(expected_out) + # print("got:") + # print(img.v) assert_rel_equal(expected_out, img.v, 5) @@ -317,16 +317,19 @@ def makemasses(i, j, k): ) weightsum += kernint * mass[p].to("g").v if weighted: - weightedsum += (kernint * mass[p].to("g").v - * density[p].to("g * cm**-3").v) + weightedsum += ( + kernint * mass[p].to("g").v * density[p].to("g * cm**-3").v + ) if weighted: baseline[i, j] += weightedsum[()] / weightsum[()] else: baseline[i, j] += weightsum[()] baseline[np.isnan(baseline)] = 0.0 - print(f"axis: {axis}, shiftcenter: {shiftcenter}, " - f"depth: {depth}, periodic: {periodic}, weighted: {weighted}, " - f"hsmlfac: {hsmlfac}") + print( + f"axis: {axis}, shiftcenter: {shiftcenter}, " + f"depth: {depth}, periodic: {periodic}, weighted: {weighted}, " + f"hsmlfac: {hsmlfac}" + ) print("expected:") print(baseline) print("got:") diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx index 12f7bc5dbbe..5aee1219fdd 100644 --- a/yt/utilities/lib/pixelization_routines.pyx +++ b/yt/utilities/lib/pixelization_routines.pyx @@ -1370,7 +1370,7 @@ def pixelize_sph_kernel_projection_pixelave( # axes x, y, and z here need not be the corresponding axes in # the simulation. Projections are always along the "z" axis here, # bounds are in this function's x, y, z order. The calling python - # function should ensure these line up with the line of sight + # function should ensure these line up with the line of sight # (here z) and projection plane (x, y) coordinate orders in the # (possible rotated) simulation coordinates. From ed579779abda6fab9603df380c56bda939656c90 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 28 Feb 2025 17:47:57 +0000 Subject: [PATCH 50/68] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- yt/geometry/coordinates/cartesian_coordinates.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/yt/geometry/coordinates/cartesian_coordinates.py b/yt/geometry/coordinates/cartesian_coordinates.py index ce82a34bdb3..2b750d87f5c 100644 --- a/yt/geometry/coordinates/cartesian_coordinates.py +++ b/yt/geometry/coordinates/cartesian_coordinates.py @@ -430,8 +430,7 @@ def _ortho_pixelize( pixelize_sph_kernel_projection_pixelave ) else: - raise NotImplementedError( - f"No pixelmeaning option {pixelmeaning}") + raise NotImplementedError(f"No pixelmeaning option {pixelmeaning}") weight = data_source.weight_field moment = data_source.moment le, re = data_source.data_source.get_bbox() From 529ad6be31577e8cf88bd1ef90c0d392a598fdd9 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Wed, 5 Mar 2025 15:15:13 +0100 Subject: [PATCH 51/68] transposed later in coordinate handler -> normal x, y axis order in backend --- yt/utilities/lib/pixelization_routines.pyx | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx index 5aee1219fdd..6e6a4baf01a 100644 --- a/yt/utilities/lib/pixelization_routines.pyx +++ b/yt/utilities/lib/pixelization_routines.pyx @@ -1605,7 +1605,7 @@ def pixelize_sph_kernel_projection_pixelave( # surface density # = density * av. length # = density * volume / area - local_buff[x0 + y0 * xsize] += ( + local_buff[x0 * ysize + y0] += ( pmass[j] / pdens[j] * ipixA * qts_j ) mask[x0, y0] = 1 @@ -1655,8 +1655,8 @@ def pixelize_sph_kernel_projection_pixelave( if kv == 0.: # don't set mask to 1 continue - local_buff[(x0 + xxi) - + (y0 + yyi) * xsize] += ( + local_buff[(x0 + xxi) * ysize + + (y0 + yyi)] += ( prefactor_j * kv ) mask[x0 + xxi, y0 + yyi] = 1 @@ -1716,7 +1716,7 @@ def pixelize_sph_kernel_projection_pixelave( q_ij2 = ((posx_diff + posy_diff) * ih_j2) if q_ij2 >= 1.0: continue - local_buff[xi + yi * xsize] += ( + local_buff[xi * ysize + yi] += ( prefactor_j * insuby * insubx * itab.interpolate(q_ij2) ) @@ -1725,7 +1725,7 @@ def pixelize_sph_kernel_projection_pixelave( with gil: for sxi in range(xsize): for syi in range(ysize): - buff[syi, sxi] += local_buff[sxi + syi * xsize] + buff[sxi, syi] += local_buff[sxi * ysize + syi] # free memory in (hopefully?) private variables assigned in # each thread free(local_buff) From f23b7e20e6fad9fbcd50572204ae93225dd64d69 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Wed, 5 Mar 2025 15:23:53 +0100 Subject: [PATCH 52/68] update test: match projection plot (y, x) axis order --- .../coordinates/tests/test_sph_pixelization_pytestonly.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py index 95228237d6c..70c3b1604fa 100644 --- a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py +++ b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py @@ -334,7 +334,7 @@ def makemasses(i, j, k): print(baseline) print("got:") print(img.v) - assert_rel_equal(baseline, img.v, 5) + assert_rel_equal(baseline.T, img.v, 5) @pytest.mark.parametrize("periodic", [True, False]) From 4d534397e83299fa7ebda0131afb1bebd39caf9e Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Wed, 5 Mar 2025 15:56:08 +0100 Subject: [PATCH 53/68] fixed *some* problems in 2x2 kernel stamps (indexing) --- yt/utilities/lib/pixelization_routines.pyx | 3 +++ 1 file changed, 3 insertions(+) diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx index 6e6a4baf01a..9383ed069e9 100644 --- a/yt/utilities/lib/pixelization_routines.pyx +++ b/yt/utilities/lib/pixelization_routines.pyx @@ -1634,6 +1634,9 @@ def pixelize_sph_kernel_projection_pixelave( 0.5 * (xn + 1.) * nsmin + 0.5) yi = math.floor( 0.5 * (yn + 1.) * nsmin + 0.5) + # in case the pixel edge was beyond hsml + xi = iclip(xi, 0, nsmin) + yi = iclip(yi, 0, nsmin) # weight for each line integral: # area of subsampling pixel in length units # divided by output grid pixel size From d7d0c5519392fa532ffd975d90f2d6b35d970325 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Wed, 5 Mar 2025 16:37:44 +0100 Subject: [PATCH 54/68] bugfix: for 2x2 stamps, was adjusting prefector_j in each periodicity loop, not just once --- yt/utilities/lib/pixelization_routines.pyx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx index 9383ed069e9..1a7e7f0c281 100644 --- a/yt/utilities/lib/pixelization_routines.pyx +++ b/yt/utilities/lib/pixelization_routines.pyx @@ -1386,7 +1386,7 @@ def pixelize_sph_kernel_projection_pixelave( cdef np.int8_t weighted = 0 # (hopefully?) private variables in parallel section cdef np.float64_t x, y, px, py, pz - cdef np.float64_t h_j2, ih_j2, qts_j, prefactor_j + cdef np.float64_t h_j2, ih_j2, qts_j, prefactor_j, kernnorm_j cdef np.float64_t q_ij2, posx_diff, posy_diff cdef np.float64_t insubx, insuby, xn, yn, kv cdef np.int64_t xi, yi, x0, x1, y0, y1, xxi, yyi @@ -1640,8 +1640,8 @@ def pixelize_sph_kernel_projection_pixelave( # weight for each line integral: # area of subsampling pixel in length units # divided by output grid pixel size - prefactor_j *= (4.0 * h_j2 - * insmin * insmin * ipixA) + kernnorm_j = (4.0 * h_j2 + * insmin * insmin * ipixA) # coarse (part. pos rel. to grid) index ci = (4 * xi * (nsmin + 1) + 4 * yi) for xxi in range(2): @@ -1660,7 +1660,7 @@ def pixelize_sph_kernel_projection_pixelave( continue local_buff[(x0 + xxi) * ysize + (y0 + yyi)] += ( - prefactor_j * kv + prefactor_j * kernnorm_j * kv ) mask[x0 + xxi, y0 + yyi] = 1 From a3725fc3a697840def2e6e74fe71b545029d475a Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Wed, 5 Mar 2025 16:55:38 +0100 Subject: [PATCH 55/68] fewer test cases for faster iteration --- .../tests/test_sph_pixelization_pytestonly.py | 23 ++++++++++++------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py index 70c3b1604fa..315b74d6d75 100644 --- a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py +++ b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py @@ -164,13 +164,14 @@ def makemasses(i, j, k): # print(img.v) assert_rel_equal(expected_out, img.v, 5) - -@pytest.mark.parametrize("axis", [0, 1, 2]) -@pytest.mark.parametrize("shiftcenter", [False, True]) +# I'm not sure if this warrants a "slow" marker, but the brute-force +# python/numpy kernel integration does take a while to run. +@pytest.mark.parametrize("axis", [0]) #, 1, 2 +@pytest.mark.parametrize("shiftcenter", [True, False]) @pytest.mark.parametrize("periodic", [True, False]) -@pytest.mark.parametrize("depth", [None, (1.0, "cm"), (5.0, "cm")]) -@pytest.mark.parametrize("weighted", [False, True]) -@pytest.mark.parametrize("hsmlfac", [0.2, 0.5, 1.0]) +@pytest.mark.parametrize("depth", [None, (1.0, "cm")]) #, (5.0, "cm") +@pytest.mark.parametrize("weighted", [False]) #, True +@pytest.mark.parametrize("hsmlfac", [0.2, 1.0]) #0.5, def test_sph_proj_pixelave_alongaxes( axis: int, shiftcenter: bool, @@ -288,7 +289,8 @@ def makemasses(i, j, k): for j in range(buff_size[1]): pixelxy = unyt.unyt_array( np.array( - [pixedges_x[i], pixedges_x[i + 1], pixedges_y[j], pixedges_y[j + 1]] + [pixedges_x[i], pixedges_x[i + 1], + pixedges_y[j], pixedges_y[j + 1]] ), "cm", ) @@ -318,7 +320,8 @@ def makemasses(i, j, k): weightsum += kernint * mass[p].to("g").v if weighted: weightedsum += ( - kernint * mass[p].to("g").v * density[p].to("g * cm**-3").v + kernint * mass[p].to("g").v + * density[p].to("g * cm**-3").v ) if weighted: baseline[i, j] += weightedsum[()] / weightsum[()] @@ -334,6 +337,10 @@ def makemasses(i, j, k): print(baseline) print("got:") print(img.v) + # in a few single-particle projection tests, mass conservation + # seems to be good to about 4 decimal places; the brute-force + # kernel integration for the baseline image may also introduce + # errors assert_rel_equal(baseline.T, img.v, 5) From 0502ba46e37667a9e2a9d2ed31d198dd902de12a Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Wed, 5 Mar 2025 17:01:38 +0100 Subject: [PATCH 56/68] try higher kernel sampling for baseline? --- .../coordinates/tests/test_sph_pixelization_pytestonly.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py index 315b74d6d75..82c28b37941 100644 --- a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py +++ b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py @@ -314,7 +314,7 @@ def makemasses(i, j, k): pixelxy.to("cm").v, pcenxy[p], (hsml[p].to("cm")).v, - nsample=50, + nsample=100, periodxy=periodxy, ) weightsum += kernint * mass[p].to("g").v From dfbe1ace27e32c94b0c0d858f9b2ec68575b3d67 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Mon, 24 Mar 2025 17:16:24 +0100 Subject: [PATCH 57/68] more mucking around with tests --- .../coordinates/tests/test_sph_pixelization_pytestonly.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py index 82c28b37941..088956975f4 100644 --- a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py +++ b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py @@ -168,8 +168,8 @@ def makemasses(i, j, k): # python/numpy kernel integration does take a while to run. @pytest.mark.parametrize("axis", [0]) #, 1, 2 @pytest.mark.parametrize("shiftcenter", [True, False]) -@pytest.mark.parametrize("periodic", [True, False]) -@pytest.mark.parametrize("depth", [None, (1.0, "cm")]) #, (5.0, "cm") +@pytest.mark.parametrize("periodic", [False]) # True, +@pytest.mark.parametrize("depth", [None]) #, (1.0, "cm"), (5.0, "cm") @pytest.mark.parametrize("weighted", [False]) #, True @pytest.mark.parametrize("hsmlfac", [0.2, 1.0]) #0.5, def test_sph_proj_pixelave_alongaxes( @@ -197,7 +197,7 @@ def test_sph_proj_pixelave_alongaxes( elif axis == 1: ax0 = 2 ax1 = 0 - buff_size = (5, 5) + buff_size = (4, 4) # test correct centering, particle selection def makemasses(i, j, k): From c51a4b2ed2c6fdba87221956789a1dce95304dba Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Fri, 28 Mar 2025 19:54:33 +0100 Subject: [PATCH 58/68] speed up tests, on-axis test pass with lowered precision reqs. --- .../tests/test_sph_pixelization_pytestonly.py | 194 ++++++++---------- 1 file changed, 87 insertions(+), 107 deletions(-) diff --git a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py index 088956975f4..475198e8c38 100644 --- a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py +++ b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py @@ -1,6 +1,7 @@ import numpy as np import pytest import unyt +from numpy.typing import NDArray import yt from yt.data_objects.selection_objects.region import YTRegion @@ -164,22 +165,31 @@ def makemasses(i, j, k): # print(img.v) assert_rel_equal(expected_out, img.v, 5) -# I'm not sure if this warrants a "slow" marker, but the brute-force -# python/numpy kernel integration does take a while to run. -@pytest.mark.parametrize("axis", [0]) #, 1, 2 +def resreduce(arr: NDArray, outshape: tuple[int, int]) -> NDArray: + ''' + gets an array of size outshape by averaging the pixel values in arr + ''' + insize = arr.shape + tempi0 = insize[0] // outshape[0] + tempi1 = insize[1] // outshape[1] + temp = arr.reshape(outshape[0], tempi0, outshape[1], tempi1) + out = np.sum(temp, axis=(1, 3)) / (tempi0 * tempi1) + return out + +@pytest.mark.parametrize("axis", [0, 1, 2]) @pytest.mark.parametrize("shiftcenter", [True, False]) -@pytest.mark.parametrize("periodic", [False]) # True, -@pytest.mark.parametrize("depth", [None]) #, (1.0, "cm"), (5.0, "cm") -@pytest.mark.parametrize("weighted", [False]) #, True -@pytest.mark.parametrize("hsmlfac", [0.2, 1.0]) #0.5, +@pytest.mark.parametrize("periodic", [True, False]) +@pytest.mark.parametrize("depth", [None, (1.0, "cm"), (5.0, "cm")]) +@pytest.mark.parametrize("hsmlfac", [0.2, 0.5, 1.0]) def test_sph_proj_pixelave_alongaxes( axis: int, shiftcenter: bool, periodic: bool, depth: float, - weighted: bool, hsmlfac: float, ) -> None: + # weighted and unweighted tested in one round: + # need weight map to downsample the baseline weighted map if shiftcenter: center = unyt.unyt_array(np.array((0.6, 0.5, 0.4)), "cm") else: @@ -187,17 +197,8 @@ def test_sph_proj_pixelave_alongaxes( bbox = unyt.unyt_array(np.array([[0.0, 3.0], [0.0, 3.0], [0.0, 3.0]]), "cm") hsml_factor = hsmlfac unitrho = 1.5 - - if axis == 2: - ax0 = 0 - ax1 = 1 - elif axis == 0: - ax0 = 1 - ax1 = 2 - elif axis == 1: - ax0 = 2 - ax1 = 0 buff_size = (4, 4) + subsample_size = (4 * 256, 4 * 256) # test correct centering, particle selection def makemasses(i, j, k): @@ -235,113 +236,92 @@ def makemasses(i, j, k): # we don't actually want a plot, it's just a straightforward, # common way to get an frb / image array - if weighted: - toweight_field = ("gas", "density") - else: - toweight_field = None - prj = yt.ProjectionPlot( + testprj = yt.ProjectionPlot( ds, axis, ("gas", "density"), width=(2.5, "cm"), - weight_field=toweight_field, + weight_field=None, buff_size=buff_size, center=center, data_source=source, pixelmeaning="pixelave", ) - print("ProjectionPlot returned") - img = prj.frb.data[("gas", "density")] - print("FRB info:") - print(prj.frb._get_info(("gas", "density"))) - - projwidth = unyt.unyt_array(np.array((2.5,)), "cm") - pixedges_x = np.linspace( - center[ax0] - 0.5 * projwidth, center[ax0] + 0.5 * projwidth, buff_size[0] + 1 + testprj_wtd = yt.ProjectionPlot( + ds, + axis, + ("gas", "density"), + width=(2.5, "cm"), + weight_field=("gas", "density"), + buff_size=buff_size, + center=center, + data_source=source, + pixelmeaning="pixelave", ) - pixedges_y = np.linspace( - center[ax1] - 0.5 * projwidth, center[ax1] + 0.5 * projwidth, buff_size[1] + 1 + testimg = testprj.frb.data[("gas", "density")] + testimg_wtd = testprj_wtd.frb.data[("gas", "density")] + + baseprj = yt.ProjectionPlot( + ds, + axis, + ("gas", "density"), + width=(2.5, "cm"), + weight_field=None, + buff_size=subsample_size, + center=center, + data_source=source, + pixelmeaning="pencilbeam", ) - _depth = depth if depth is not None else unyt.unyt_array(np.array((np.inf,)), "cm") - zlims = (center[axis] - 0.5 * _depth, center[axis] + 0.5 * _depth) - ad = ds.all_data() - pcenxy = np.array( - [ - (ad[("gas", "xyz"[ax0])]).to("cm"), - (ad[("gas", "xyz"[ax1])]).to("cm"), - ] - ).T - projz = (ad[("gas", "xyz"[axis])]).to("cm") - hsml = (ad[("gas", "smoothing_length")]).to("cm") - mass = (ad[("gas", "mass")]).to("g") - density = (ad[("gas", "density")]).to("g * cm**-3") - baseline = np.zeros(buff_size) - if periodic: - periodxy = ( - (bbox[ax0][1] - bbox[ax0][0]).to("cm").v, - (bbox[ax1][1] - bbox[ax1][0]).to("cm").v, - ) - periodz = bbox[axis][1] - bbox[axis][0] - else: - periodxy = (None, None) - periodz = None - for i in range(buff_size[0]): - for j in range(buff_size[1]): - pixelxy = unyt.unyt_array( - np.array( - [pixedges_x[i], pixedges_x[i + 1], - pixedges_y[j], pixedges_y[j + 1]] - ), - "cm", - ) - weightsum = 0.0 - weightedsum = 0.0 - for p in range(pcenxy.shape[0]): - # check if the particle is excluded along z - if projz[p] < zlims[0] or projz[p] > zlims[1]: - if periodz is None: - continue - else: - dz = projz[p] - center[axis].to("cm") - dz += 0.5 * periodz - dz %= periodz - dz -= periodz - if np.abs(dz) > 0.5 * _depth: - continue - - kernint = pixelintegrate_kernel( - cubicspline_python, - pixelxy.to("cm").v, - pcenxy[p], - (hsml[p].to("cm")).v, - nsample=100, - periodxy=periodxy, - ) - weightsum += kernint * mass[p].to("g").v - if weighted: - weightedsum += ( - kernint * mass[p].to("g").v - * density[p].to("g * cm**-3").v - ) - if weighted: - baseline[i, j] += weightedsum[()] / weightsum[()] - else: - baseline[i, j] += weightsum[()] - baseline[np.isnan(baseline)] = 0.0 + baseprj = yt.ProjectionPlot( + ds, + axis, + ("gas", "density"), + width=(2.5, "cm"), + weight_field=None, + buff_size=subsample_size, + center=center, + data_source=source, + pixelmeaning="pencilbeam", + ) + baseprj_wtd = yt.ProjectionPlot( + ds, + axis, + ("gas", "density"), + width=(2.5, "cm"), + weight_field=("gas", "density"), + buff_size=subsample_size, + center=center, + data_source=source, + pixelmeaning="pencilbeam", + ) + _baseimg = baseprj.frb.data[("gas", "density")] + baseimg = resreduce(_baseimg, buff_size) + _baseimg_wtd = baseprj_wtd.frb.data[("gas", "density")] + baseimg_wtd = resreduce(_baseimg * _baseimg_wtd, buff_size) / baseimg + baseimg_wtd[baseimg == 0.] = 0. # handle 0./0. NaNs + print( f"axis: {axis}, shiftcenter: {shiftcenter}, " - f"depth: {depth}, periodic: {periodic}, weighted: {weighted}, " + f"depth: {depth}, periodic: {periodic}, " f"hsmlfac: {hsmlfac}" ) print("expected:") - print(baseline) + print(baseimg.v) print("got:") - print(img.v) + print(testimg.v) + print("expected (weighted):") + print(baseimg_wtd.v) + print("got (weighted):") + print(testimg_wtd.v) # in a few single-particle projection tests, mass conservation - # seems to be good to about 4 decimal places; the brute-force - # kernel integration for the baseline image may also introduce - # errors - assert_rel_equal(baseline.T, img.v, 5) + # seems to be good to about 4 decimal places; agreement with + # subsampling with multiple particles is not quite at 4. + # pixel values seem to converge more slowly though, so we test + # mass conservation explicitly, and pixel agreement at low + # precision + assert_rel_equal(baseimg.v, testimg.v, 1) + assert_rel_equal(np.sum(baseimg.v), np.sum(testimg.v), 2) + assert_rel_equal(baseimg_wtd.v, testimg_wtd.v, 1) @pytest.mark.parametrize("periodic", [True, False]) From 6d7eda467dd634797ec941c5111dd20bc5745a15 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Fri, 28 Mar 2025 21:07:44 +0100 Subject: [PATCH 59/68] fix div. by 0 error test failure? --- .../coordinates/tests/test_sph_pixelization_pytestonly.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py index 475198e8c38..40bc20c9bf6 100644 --- a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py +++ b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py @@ -297,7 +297,9 @@ def makemasses(i, j, k): _baseimg = baseprj.frb.data[("gas", "density")] baseimg = resreduce(_baseimg, buff_size) _baseimg_wtd = baseprj_wtd.frb.data[("gas", "density")] - baseimg_wtd = resreduce(_baseimg * _baseimg_wtd, buff_size) / baseimg + _divimg = np.copy(baseimg.v) + _divimg[_divimg == 0.] = -1. # avoid div. by 0 test failures + baseimg_wtd = resreduce(_baseimg * _baseimg_wtd.v, buff_size) / _divimg baseimg_wtd[baseimg == 0.] = 0. # handle 0./0. NaNs print( From dd0d9f50317c83076a45d5cc7f13cf1798c2836b Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Tue, 1 Apr 2025 15:37:16 +0200 Subject: [PATCH 60/68] remove function not used in updated test --- yt/testing.py | 67 --------------------------------------------------- 1 file changed, 67 deletions(-) diff --git a/yt/testing.py b/yt/testing.py index afbbf6cf531..82c71bc29e1 100644 --- a/yt/testing.py +++ b/yt/testing.py @@ -159,73 +159,6 @@ def integrate_kernel( return np.atleast_1d(pre * integral) -def pixelintegrate_kernel( - kernelfunc: Callable[[float | np.ndarray], np.ndarray], - pixelxy: np.ndarray, - pcenxy: np.ndarray, - hsml: float, - nsample: int = 50, - periodxy: tuple[float | None, float | None] = (True, True), -) -> float: - """ - integrates a kernel function over a rectangular prism. - - Returns (kernel_volume_fraction) / pixel area, - where the kernel_volume_fraction is the volume integral of the - kernel over the rectangular prism defined by the pixel, divided - by the volume integral of the kernel over its entire support. - - Parameters - ---------- - kernelfunc: - the kernel function to integrate - pixelxy: - corners of the pixel in the output grid; - (xmin, xmax, ymin, ymax); [length units] - pcenxy: - sph particle center in the output grid plane - (x, y) [same units as pixelxy] - hsml: - smoothing length [same units as pixelxy] - nsample: - number of subsamples used within the pixel - in each grid direction - - Returns - -------- - integral of the SPH kernel function / pixel area - units: (1 / input length units)**2 - - Note - ---- - The calculation assumes the kernel support is entirely - within the rectangular prism defined by the pixel in the - line-of-sight direction, as the backend functions for SPH - projections do. - """ - xedges = np.linspace(pixelxy[0], pixelxy[1], nsample + 1) - xcens = 0.5 * (xedges[:-1] + xedges[1:]) - dx = np.diff(xedges, axis=0) - yedges = np.linspace(pixelxy[2], pixelxy[3], nsample + 1) - ycens = 0.5 * (yedges[:-1] + yedges[1:]) - dy = np.diff(yedges, axis=0) - - xoff = xcens - pcenxy[0] - if periodxy[0] is not None: - xoff += 0.5 * periodxy[0] - xoff %= periodxy[0] - xoff -= 0.5 * periodxy[0] - yoff = ycens - pcenxy[1] - if periodxy[1] is not None: - yoff += 0.5 * periodxy[1] - yoff %= periodxy[1] - yoff -= 0.5 * periodxy[1] - bpars = np.sqrt(xoff[:, np.newaxis, ...] ** 2 + yoff[np.newaxis, :, ...] ** 2) - lineints = integrate_kernel(kernelfunc, bpars, hsml, nsample=nsample) - volint = np.sum(lineints * dx[:, np.newaxis] * dy[np.newaxis, :]) - return volint / (pixelxy[1] - pixelxy[0]) / (pixelxy[3] - pixelxy[2]) - - _zeroperiods = np.array([0.0, 0.0, 0.0]) From 2007da43c0fa65b957fc7060a89f559eaf5e8027 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Tue, 1 Apr 2025 15:57:29 +0200 Subject: [PATCH 61/68] added function to reduce image res by averaging --- yt/testing.py | 38 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) diff --git a/yt/testing.py b/yt/testing.py index 82c71bc29e1..5050d263c2d 100644 --- a/yt/testing.py +++ b/yt/testing.py @@ -18,6 +18,7 @@ import numpy.typing as npt from more_itertools import always_iterable from numpy.random import RandomState +from numpy.typing import NDArray from unyt.exceptions import UnitOperationError from yt._maintenance.deprecation import issue_deprecation_warning @@ -159,6 +160,43 @@ def integrate_kernel( return np.atleast_1d(pre * integral) +def resreduce_image( + arr: NDArray[np.float32], + outshape: tuple[int, int] +) -> NDArray[np.float32]: + ''' + gets an image array of size outshape by averaging the pixel values + in arr over contiguous blocks in each array dimension. + + Parameters: + ----------- + arr: 2d-array + the input array for averaging, size (Nx, Ny) + outshape: + the shape of the output array + + Returns: + -------- + an array of shape `outshape` + + Notes: + ------ + values in `arr` are averaged over contiguous blocks of size + (Nsx, Nsy), where Nsx = Nx // outshape[0], + and Nsy = Ny // outshape[1]. + ''' + insize = arr.shape + if insize[0] % outshape[0] != 0 or insize[1] % outshape[1] != 0: + raise ValueError("desired output array shape must use integer" + "divisors of the input array shape. " + f"desired {outshape} doesn't divide " + f"input array shape {arr.shape()}") + tempi0 = insize[0] // outshape[0] + tempi1 = insize[1] // outshape[1] + temp = arr.reshape(outshape[0], tempi0, outshape[1], tempi1) + out = np.sum(temp, axis=(1, 3)) / (tempi0 * tempi1) + return out + _zeroperiods = np.array([0.0, 0.0, 0.0]) From a8f5a9217227d68c77ad0b6664a6202114a7920e Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Tue, 1 Apr 2025 16:28:31 +0200 Subject: [PATCH 62/68] added off-axis projection test for pixave; passes --- .../test_offaxisprojection_pytestonly.py | 139 ++++++++++++++++++ 1 file changed, 139 insertions(+) diff --git a/yt/visualization/tests/test_offaxisprojection_pytestonly.py b/yt/visualization/tests/test_offaxisprojection_pytestonly.py index ac3070ace50..a0a566de631 100644 --- a/yt/visualization/tests/test_offaxisprojection_pytestonly.py +++ b/yt/visualization/tests/test_offaxisprojection_pytestonly.py @@ -3,6 +3,7 @@ import unyt from numpy.testing import assert_allclose +from yt.data_objects.selection_objects.region import YTRegion from yt.testing import ( assert_rel_equal, cubicspline_python, @@ -11,6 +12,7 @@ fake_sph_grid_ds, integrate_kernel, requires_module_pytest, + resreduce_image, ) from yt.visualization.api import OffAxisProjectionPlot, ProjectionPlot @@ -172,6 +174,143 @@ def makemasses(i, j, k): assert_rel_equal(expected_out, img.v, 4) +@pytest.mark.parametrize("axis", [(0., 1e-4, 1.), (0.2, 0., -0.8)]) +@pytest.mark.parametrize("north_vector", [None, (1.0e-4, 1.0, 0.0)]) +@pytest.mark.parametrize("shiftcenter", [True, False]) +@pytest.mark.parametrize("periodic", [True, False]) +@pytest.mark.parametrize("depth", [None, (1.0, "cm"), (5.0, "cm")]) +@pytest.mark.parametrize("hsmlfac", [0.2, 0.5, 1.0]) +def test_sph_proj_pixelave_offaxes( + axis: tuple[float, float, float], + north_vector: tuple[float, float, float] | None, + shiftcenter: bool, + periodic: bool, + depth: float, + hsmlfac: float, +) -> None: + # weighted and unweighted tested in one round: + # need weight map to downsample the baseline weighted map + if shiftcenter: + center = unyt.unyt_array(np.array((0.6, 0.5, 0.4)), "cm") + else: + center = unyt.unyt_array(np.array((1.5, 1.5, 1.5)), "cm") + bbox = unyt.unyt_array(np.array([[0.0, 3.0], [0.0, 3.0], [0.0, 3.0]]), + "cm") + hsml_factor = hsmlfac + unitrho = 1.5 + buff_size = (4, 4) + subsample_size = (4 * 256, 4 * 256) + + # test correct centering, particle selection + def makemasses(i, j, k): + if i == 0 and j == 1 and k == 2: + return 2.0 + else: + return 1.0 + + # result shouldn't depend explicitly on the center if we re-center + # the data, unless we get cut-offs in the non-periodic case + ds = fake_sph_flexible_grid_ds( + hsml_factor=hsml_factor, + nperside=3, + periodic=periodic, + offsets=np.full(3, 0.5), + massgenerator=makemasses, + unitrho=unitrho, + bbox=bbox.v, + recenter=center.v, + ) + source = ds.all_data() + + # we don't actually want a plot, it's just a straightforward, + # common way to get an frb / image array + testprj = ProjectionPlot( + ds, + axis, + ("gas", "density"), + width=(2.5, "cm"), + weight_field=None, + north_vector=north_vector, + buff_size=buff_size, + center=center, + data_source=source, + pixelmeaning="pixelave", + depth=depth, + ) + testprj_wtd = ProjectionPlot( + ds, + axis, + ("gas", "density"), + width=(2.5, "cm"), + weight_field=("gas", "density"), + north_vector=north_vector, + buff_size=buff_size, + center=center, + data_source=source, + pixelmeaning="pixelave", + depth=depth, + ) + testimg = testprj.frb.data[("gas", "density")] + testimg_wtd = testprj_wtd.frb.data[("gas", "density")] + + baseprj = ProjectionPlot( + ds, + axis, + ("gas", "density"), + width=(2.5, "cm"), + weight_field=None, + buff_size=subsample_size, + north_vector=north_vector, + center=center, + data_source=source, + pixelmeaning="pencilbeam", + depth=depth, + ) + baseprj_wtd = ProjectionPlot( + ds, + axis, + ("gas", "density"), + width=(2.5, "cm"), + weight_field=("gas", "density"), + buff_size=subsample_size, + north_vector=north_vector, + center=center, + data_source=source, + pixelmeaning="pencilbeam", + depth=depth, + ) + _baseimg = baseprj.frb.data[("gas", "density")] + baseimg = resreduce_image(_baseimg, buff_size) + _baseimg_wtd = baseprj_wtd.frb.data[("gas", "density")] + _divimg = np.copy(baseimg.v) + _divimg[_divimg == 0.] = -1. # avoid div. by 0 test failures + baseimg_wtd = (resreduce_image(_baseimg * _baseimg_wtd.v, buff_size) + / _divimg) + + print( + f"axis: {axis}, shiftcenter: {shiftcenter}, " + f"depth: {depth}, periodic: {periodic}, " + f"hsmlfac: {hsmlfac}" + ) + print("expected:") + print(baseimg.v) + print("got:") + print(testimg.v) + print("expected (weighted):") + print(baseimg_wtd.v) + print("got (weighted):") + print(testimg_wtd.v) + # in a few single-particle projection tests, mass conservation + # seems to be good to about 4 decimal places; agreement with + # subsampling with multiple particles is not quite at 4. + # pixel values seem to converge more slowly though, so we test + # mass conservation explicitly, and pixel agreement at low + # precision + assert_rel_equal(baseimg.v, testimg.v, 1) + assert_rel_equal(np.sum(baseimg.v), np.sum(testimg.v), 2) + assert_rel_equal(baseimg_wtd.v, testimg_wtd.v, 1) + + _diag_dist = np.sqrt(3.0) # diagonal distance of a cube with length 1. # each case is depth, center, expected integrated distance _cases_to_test = [ From 395527dde0bab4d88de4e9cc7730d3f3b6bad885 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Tue, 1 Apr 2025 16:28:56 +0200 Subject: [PATCH 63/68] moved resreduce to testing.py, renamed --- .../tests/test_sph_pixelization_pytestonly.py | 30 ++++--------------- 1 file changed, 5 insertions(+), 25 deletions(-) diff --git a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py index 40bc20c9bf6..decf325101e 100644 --- a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py +++ b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py @@ -12,7 +12,7 @@ fake_random_sph_ds, fake_sph_flexible_grid_ds, integrate_kernel, - pixelintegrate_kernel, + resreduce_image, ) @@ -165,16 +165,7 @@ def makemasses(i, j, k): # print(img.v) assert_rel_equal(expected_out, img.v, 5) -def resreduce(arr: NDArray, outshape: tuple[int, int]) -> NDArray: - ''' - gets an array of size outshape by averaging the pixel values in arr - ''' - insize = arr.shape - tempi0 = insize[0] // outshape[0] - tempi1 = insize[1] // outshape[1] - temp = arr.reshape(outshape[0], tempi0, outshape[1], tempi1) - out = np.sum(temp, axis=(1, 3)) / (tempi0 * tempi1) - return out + @pytest.mark.parametrize("axis", [0, 1, 2]) @pytest.mark.parametrize("shiftcenter", [True, False]) @@ -261,17 +252,6 @@ def makemasses(i, j, k): testimg = testprj.frb.data[("gas", "density")] testimg_wtd = testprj_wtd.frb.data[("gas", "density")] - baseprj = yt.ProjectionPlot( - ds, - axis, - ("gas", "density"), - width=(2.5, "cm"), - weight_field=None, - buff_size=subsample_size, - center=center, - data_source=source, - pixelmeaning="pencilbeam", - ) baseprj = yt.ProjectionPlot( ds, axis, @@ -295,12 +275,12 @@ def makemasses(i, j, k): pixelmeaning="pencilbeam", ) _baseimg = baseprj.frb.data[("gas", "density")] - baseimg = resreduce(_baseimg, buff_size) + baseimg = resreduce_image(_baseimg, buff_size) _baseimg_wtd = baseprj_wtd.frb.data[("gas", "density")] _divimg = np.copy(baseimg.v) _divimg[_divimg == 0.] = -1. # avoid div. by 0 test failures - baseimg_wtd = resreduce(_baseimg * _baseimg_wtd.v, buff_size) / _divimg - baseimg_wtd[baseimg == 0.] = 0. # handle 0./0. NaNs + baseimg_wtd = (resreduce_image(_baseimg * _baseimg_wtd.v, buff_size) + / _divimg) print( f"axis: {axis}, shiftcenter: {shiftcenter}, " From d73e07e62f6470f632c7f42718735129d4071d3a Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Tue, 1 Apr 2025 16:43:43 +0200 Subject: [PATCH 64/68] add test for pixave mass conservation --- .../tests/test_sph_pixelization_pytestonly.py | 73 +++++++++++++++++++ 1 file changed, 73 insertions(+) diff --git a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py index decf325101e..79482e43804 100644 --- a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py +++ b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py @@ -5,6 +5,7 @@ import yt from yt.data_objects.selection_objects.region import YTRegion +from yt.loaders import load_particles from yt.testing import ( assert_rel_equal, cubicspline_python, @@ -306,6 +307,78 @@ def makemasses(i, j, k): assert_rel_equal(baseimg_wtd.v, testimg_wtd.v, 1) +def test_massconservation_pixave(): + bbox = np.array([[0., 3.], [0., 3.], [0., 3.]]) + pz = 1.5 + periodic = True + periods = 3. + nrandom = 50 + + # random centers, three pixel size + hsml sets + # test three cases in the backend routine + rng = np.random.default_rng() + pxs = rng.uniform(0., 3., nrandom) + pys = rng.uniform(0., 3., nrandom) + # particles < pixels: overlap 1, 2x1, or 2x2 pixels + hsmls1 = rng.uniform(0.01, 0.05, nrandom) + outsize1 = (7, 7) + # particles ~ pixels: typical 2x2 overlaps + hsmls2 = rng.uniform(0.1, 0.25, nrandom) + outsize2 = (7, 7) + # particles > pixels: + hsmls3 = rng.uniform(0.07, 0.8, nrandom) + outsize3 = (50, 50) + + for i, (hsmls, outsize) in enumerate([(hsmls1, outsize1), + (hsmls2, outsize2), + (hsmls3, outsize3), + ]): + + masses_rel = [] + for i in range(nrandom): + data = { + "particle_position_x": (np.array([pxs[i]]), "cm"), + "particle_position_y": (np.array([pys[i]]), "cm"), + "particle_position_z": (np.array([pz]), "cm"), + "particle_mass": (np.array([1.]), "g"), + "particle_velocity_x": (np.zeros(1), "cm/s"), + "particle_velocity_y": (np.zeros(1), "cm/s"), + "particle_velocity_z": (np.zeros(1), "cm/s"), + "smoothing_length": (np.ones(1) * hsmls[i], "cm"), + "density": (np.array([1.5]), "g/cm**3"), + } + ds = load_particles( + data=data, + bbox=bbox, + periodicity=(periodic,) * 3, + length_unit=1.0, + mass_unit=1.0, + time_unit=1.0, + velocity_unit=1.0, + ) + ds.kernel_name = "cubic" + + proj = ds.proj(("gas", "density"), 2) + frb = proj.to_frb( + width=(3., "cm"), + resolution=outsize, + height=(3., "cm"), + center=np.array([1.5, 1.5, 1.5]), + periodic=True, + pixelmeaning="pixelave", + ) + out = frb.get_image(("gas", "density")) + mass_in = 1. + mass_out = np.sum(out.v) * periods**2 / np.prod(outsize) + masses_rel.append(mass_out / mass_in - 1.) + masses_rel = np.array(masses_rel) + minrel = np.min(masses_rel) + maxrel = np.max(masses_rel) + print(f"Mass conservation test pixel/hsml set {i}:" + " mass conservation deviations fractions" + f" {minrel:.2e} -- {maxrel:.2e} ") + assert np.all(np.abs(masses_rel) < 1e-4) + @pytest.mark.parametrize("periodic", [True, False]) @pytest.mark.parametrize("shiftcenter", [False, True]) @pytest.mark.parametrize("zoff", [0.0, 0.1, 0.5, 1.0]) From b5360a0b0008928e4346d87ccbbb3c408b64de89 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 1 Apr 2025 14:52:14 +0000 Subject: [PATCH 65/68] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../tests/test_sph_pixelization_pytestonly.py | 56 ++++++++++--------- yt/testing.py | 18 +++--- .../test_offaxisprojection_pytestonly.py | 17 +++--- 3 files changed, 46 insertions(+), 45 deletions(-) diff --git a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py index 79482e43804..27e9953c621 100644 --- a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py +++ b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py @@ -1,7 +1,6 @@ import numpy as np import pytest import unyt -from numpy.typing import NDArray import yt from yt.data_objects.selection_objects.region import YTRegion @@ -167,12 +166,11 @@ def makemasses(i, j, k): assert_rel_equal(expected_out, img.v, 5) - @pytest.mark.parametrize("axis", [0, 1, 2]) @pytest.mark.parametrize("shiftcenter", [True, False]) @pytest.mark.parametrize("periodic", [True, False]) @pytest.mark.parametrize("depth", [None, (1.0, "cm"), (5.0, "cm")]) -@pytest.mark.parametrize("hsmlfac", [0.2, 0.5, 1.0]) +@pytest.mark.parametrize("hsmlfac", [0.2, 0.5, 1.0]) def test_sph_proj_pixelave_alongaxes( axis: int, shiftcenter: bool, @@ -180,7 +178,7 @@ def test_sph_proj_pixelave_alongaxes( depth: float, hsmlfac: float, ) -> None: - # weighted and unweighted tested in one round: + # weighted and unweighted tested in one round: # need weight map to downsample the baseline weighted map if shiftcenter: center = unyt.unyt_array(np.array((0.6, 0.5, 0.4)), "cm") @@ -279,9 +277,8 @@ def makemasses(i, j, k): baseimg = resreduce_image(_baseimg, buff_size) _baseimg_wtd = baseprj_wtd.frb.data[("gas", "density")] _divimg = np.copy(baseimg.v) - _divimg[_divimg == 0.] = -1. # avoid div. by 0 test failures - baseimg_wtd = (resreduce_image(_baseimg * _baseimg_wtd.v, buff_size) - / _divimg) + _divimg[_divimg == 0.0] = -1.0 # avoid div. by 0 test failures + baseimg_wtd = resreduce_image(_baseimg * _baseimg_wtd.v, buff_size) / _divimg print( f"axis: {axis}, shiftcenter: {shiftcenter}, " @@ -301,46 +298,48 @@ def makemasses(i, j, k): # subsampling with multiple particles is not quite at 4. # pixel values seem to converge more slowly though, so we test # mass conservation explicitly, and pixel agreement at low - # precision + # precision assert_rel_equal(baseimg.v, testimg.v, 1) assert_rel_equal(np.sum(baseimg.v), np.sum(testimg.v), 2) assert_rel_equal(baseimg_wtd.v, testimg_wtd.v, 1) def test_massconservation_pixave(): - bbox = np.array([[0., 3.], [0., 3.], [0., 3.]]) + bbox = np.array([[0.0, 3.0], [0.0, 3.0], [0.0, 3.0]]) pz = 1.5 periodic = True - periods = 3. + periods = 3.0 nrandom = 50 - + # random centers, three pixel size + hsml sets # test three cases in the backend routine rng = np.random.default_rng() - pxs = rng.uniform(0., 3., nrandom) - pys = rng.uniform(0., 3., nrandom) + pxs = rng.uniform(0.0, 3.0, nrandom) + pys = rng.uniform(0.0, 3.0, nrandom) # particles < pixels: overlap 1, 2x1, or 2x2 pixels hsmls1 = rng.uniform(0.01, 0.05, nrandom) outsize1 = (7, 7) # particles ~ pixels: typical 2x2 overlaps hsmls2 = rng.uniform(0.1, 0.25, nrandom) outsize2 = (7, 7) - # particles > pixels: + # particles > pixels: hsmls3 = rng.uniform(0.07, 0.8, nrandom) outsize3 = (50, 50) - - for i, (hsmls, outsize) in enumerate([(hsmls1, outsize1), - (hsmls2, outsize2), - (hsmls3, outsize3), - ]): + for i, (hsmls, outsize) in enumerate( + [ + (hsmls1, outsize1), + (hsmls2, outsize2), + (hsmls3, outsize3), + ] + ): masses_rel = [] for i in range(nrandom): data = { "particle_position_x": (np.array([pxs[i]]), "cm"), "particle_position_y": (np.array([pys[i]]), "cm"), "particle_position_z": (np.array([pz]), "cm"), - "particle_mass": (np.array([1.]), "g"), + "particle_mass": (np.array([1.0]), "g"), "particle_velocity_x": (np.zeros(1), "cm/s"), "particle_velocity_y": (np.zeros(1), "cm/s"), "particle_velocity_z": (np.zeros(1), "cm/s"), @@ -360,25 +359,28 @@ def test_massconservation_pixave(): proj = ds.proj(("gas", "density"), 2) frb = proj.to_frb( - width=(3., "cm"), + width=(3.0, "cm"), resolution=outsize, - height=(3., "cm"), + height=(3.0, "cm"), center=np.array([1.5, 1.5, 1.5]), periodic=True, pixelmeaning="pixelave", ) out = frb.get_image(("gas", "density")) - mass_in = 1. + mass_in = 1.0 mass_out = np.sum(out.v) * periods**2 / np.prod(outsize) - masses_rel.append(mass_out / mass_in - 1.) + masses_rel.append(mass_out / mass_in - 1.0) masses_rel = np.array(masses_rel) minrel = np.min(masses_rel) maxrel = np.max(masses_rel) - print(f"Mass conservation test pixel/hsml set {i}:" - " mass conservation deviations fractions" - f" {minrel:.2e} -- {maxrel:.2e} ") + print( + f"Mass conservation test pixel/hsml set {i}:" + " mass conservation deviations fractions" + f" {minrel:.2e} -- {maxrel:.2e} " + ) assert np.all(np.abs(masses_rel) < 1e-4) + @pytest.mark.parametrize("periodic", [True, False]) @pytest.mark.parametrize("shiftcenter", [False, True]) @pytest.mark.parametrize("zoff", [0.0, 0.1, 0.5, 1.0]) diff --git a/yt/testing.py b/yt/testing.py index 5050d263c2d..813a8e381ba 100644 --- a/yt/testing.py +++ b/yt/testing.py @@ -161,10 +161,9 @@ def integrate_kernel( def resreduce_image( - arr: NDArray[np.float32], - outshape: tuple[int, int] + arr: NDArray[np.float32], outshape: tuple[int, int] ) -> NDArray[np.float32]: - ''' + """ gets an image array of size outshape by averaging the pixel values in arr over contiguous blocks in each array dimension. @@ -184,19 +183,22 @@ def resreduce_image( values in `arr` are averaged over contiguous blocks of size (Nsx, Nsy), where Nsx = Nx // outshape[0], and Nsy = Ny // outshape[1]. - ''' + """ insize = arr.shape if insize[0] % outshape[0] != 0 or insize[1] % outshape[1] != 0: - raise ValueError("desired output array shape must use integer" - "divisors of the input array shape. " - f"desired {outshape} doesn't divide " - f"input array shape {arr.shape()}") + raise ValueError( + "desired output array shape must use integer" + "divisors of the input array shape. " + f"desired {outshape} doesn't divide " + f"input array shape {arr.shape()}" + ) tempi0 = insize[0] // outshape[0] tempi1 = insize[1] // outshape[1] temp = arr.reshape(outshape[0], tempi0, outshape[1], tempi1) out = np.sum(temp, axis=(1, 3)) / (tempi0 * tempi1) return out + _zeroperiods = np.array([0.0, 0.0, 0.0]) diff --git a/yt/visualization/tests/test_offaxisprojection_pytestonly.py b/yt/visualization/tests/test_offaxisprojection_pytestonly.py index a0a566de631..8e773375da6 100644 --- a/yt/visualization/tests/test_offaxisprojection_pytestonly.py +++ b/yt/visualization/tests/test_offaxisprojection_pytestonly.py @@ -3,7 +3,6 @@ import unyt from numpy.testing import assert_allclose -from yt.data_objects.selection_objects.region import YTRegion from yt.testing import ( assert_rel_equal, cubicspline_python, @@ -174,12 +173,12 @@ def makemasses(i, j, k): assert_rel_equal(expected_out, img.v, 4) -@pytest.mark.parametrize("axis", [(0., 1e-4, 1.), (0.2, 0., -0.8)]) +@pytest.mark.parametrize("axis", [(0.0, 1e-4, 1.0), (0.2, 0.0, -0.8)]) @pytest.mark.parametrize("north_vector", [None, (1.0e-4, 1.0, 0.0)]) @pytest.mark.parametrize("shiftcenter", [True, False]) @pytest.mark.parametrize("periodic", [True, False]) @pytest.mark.parametrize("depth", [None, (1.0, "cm"), (5.0, "cm")]) -@pytest.mark.parametrize("hsmlfac", [0.2, 0.5, 1.0]) +@pytest.mark.parametrize("hsmlfac", [0.2, 0.5, 1.0]) def test_sph_proj_pixelave_offaxes( axis: tuple[float, float, float], north_vector: tuple[float, float, float] | None, @@ -188,14 +187,13 @@ def test_sph_proj_pixelave_offaxes( depth: float, hsmlfac: float, ) -> None: - # weighted and unweighted tested in one round: + # weighted and unweighted tested in one round: # need weight map to downsample the baseline weighted map if shiftcenter: center = unyt.unyt_array(np.array((0.6, 0.5, 0.4)), "cm") else: center = unyt.unyt_array(np.array((1.5, 1.5, 1.5)), "cm") - bbox = unyt.unyt_array(np.array([[0.0, 3.0], [0.0, 3.0], [0.0, 3.0]]), - "cm") + bbox = unyt.unyt_array(np.array([[0.0, 3.0], [0.0, 3.0], [0.0, 3.0]]), "cm") hsml_factor = hsmlfac unitrho = 1.5 buff_size = (4, 4) @@ -283,9 +281,8 @@ def makemasses(i, j, k): baseimg = resreduce_image(_baseimg, buff_size) _baseimg_wtd = baseprj_wtd.frb.data[("gas", "density")] _divimg = np.copy(baseimg.v) - _divimg[_divimg == 0.] = -1. # avoid div. by 0 test failures - baseimg_wtd = (resreduce_image(_baseimg * _baseimg_wtd.v, buff_size) - / _divimg) + _divimg[_divimg == 0.0] = -1.0 # avoid div. by 0 test failures + baseimg_wtd = resreduce_image(_baseimg * _baseimg_wtd.v, buff_size) / _divimg print( f"axis: {axis}, shiftcenter: {shiftcenter}, " @@ -305,7 +302,7 @@ def makemasses(i, j, k): # subsampling with multiple particles is not quite at 4. # pixel values seem to converge more slowly though, so we test # mass conservation explicitly, and pixel agreement at low - # precision + # precision assert_rel_equal(baseimg.v, testimg.v, 1) assert_rel_equal(np.sum(baseimg.v), np.sum(testimg.v), 2) assert_rel_equal(baseimg_wtd.v, testimg_wtd.v, 1) From 7ac336f7a02bb41505acd44abcbb6eb47127bd6a Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Mon, 7 Apr 2025 16:54:31 +0200 Subject: [PATCH 66/68] fix tuple annotation? --- yt/testing.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/yt/testing.py b/yt/testing.py index 813a8e381ba..3cd1d1635d8 100644 --- a/yt/testing.py +++ b/yt/testing.py @@ -10,7 +10,7 @@ from functools import wraps from importlib.util import find_spec from shutil import which -from typing import TYPE_CHECKING, TypeVar +from typing import Tuple, TYPE_CHECKING, TypeVar from unittest import SkipTest import matplotlib @@ -161,7 +161,7 @@ def integrate_kernel( def resreduce_image( - arr: NDArray[np.float32], outshape: tuple[int, int] + arr: NDArray[np.float32], outshape: Tuple[int, int] ) -> NDArray[np.float32]: """ gets an image array of size outshape by averaging the pixel values From e4c908bf3a0f414dd20e010c9db84454fe5be997 Mon Sep 17 00:00:00 2001 From: nastasha-w <35771073+nastasha-w@users.noreply.github.com> Date: Mon, 7 Apr 2025 16:57:52 +0200 Subject: [PATCH 67/68] update docstring --- yt/data_objects/selection_objects/slices.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/yt/data_objects/selection_objects/slices.py b/yt/data_objects/selection_objects/slices.py index ed690517241..bdfe9dbe105 100644 --- a/yt/data_objects/selection_objects/slices.py +++ b/yt/data_objects/selection_objects/slices.py @@ -347,7 +347,8 @@ def to_frb( periodic : boolean This can be true or false, and governs whether the pixelization will span the domain boundaries. - pixelmeaning: ignored, arg. meant for SPH projections + pixelmeaning: ignored + argument meant for SPH projections Returns ------- From 3befdb79651c2bc3a47aac388cd035e8d15c39cb Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 7 Apr 2025 15:07:25 +0000 Subject: [PATCH 68/68] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- yt/testing.py | 6 +++--- yt/utilities/lib/pixelization_routines.pyx | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/yt/testing.py b/yt/testing.py index 3cd1d1635d8..f76e08935be 100644 --- a/yt/testing.py +++ b/yt/testing.py @@ -10,7 +10,7 @@ from functools import wraps from importlib.util import find_spec from shutil import which -from typing import Tuple, TYPE_CHECKING, TypeVar +from typing import TYPE_CHECKING, TypeVar from unittest import SkipTest import matplotlib @@ -125,7 +125,7 @@ def integrate_kernel( ], b: float | npt.NDArray[np.floating], hsml: float | npt.NDArray[np.floating], - nsample: int = 500, + nsample: int = 500, ) -> npt.NDArray[np.floating]: """ integrates a kernel function over a line passing entirely @@ -161,7 +161,7 @@ def integrate_kernel( def resreduce_image( - arr: NDArray[np.float32], outshape: Tuple[int, int] + arr: NDArray[np.float32], outshape: tuple[int, int] ) -> NDArray[np.float32]: """ gets an image array of size outshape by averaging the pixel values diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx index 1a7e7f0c281..f9957dcbff5 100644 --- a/yt/utilities/lib/pixelization_routines.pyx +++ b/yt/utilities/lib/pixelization_routines.pyx @@ -1408,7 +1408,7 @@ def pixelize_sph_kernel_projection_pixelave( period_z = period[2] for si in range(3): check_period[si] = _check_period[si] - + # we find the x and y range over which we have pixels and we find how many # pixels we have in each dimension xsize, ysize = buff.shape[0], buff.shape[1]