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' diff --git a/yt/data_objects/selection_objects/data_selection_objects.py b/yt/data_objects/selection_objects/data_selection_objects.py index fc47b693dc0..e99ebed4e5c 100644 --- a/yt/data_objects/selection_objects/data_selection_objects.py +++ b/yt/data_objects/selection_objects/data_selection_objects.py @@ -4,6 +4,7 @@ import uuid from collections import defaultdict from contextlib import contextmanager +from typing import Literal import numpy as np from more_itertools import always_iterable @@ -567,7 +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): + 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. @@ -595,7 +604,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 +676,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/data_objects/selection_objects/slices.py b/yt/data_objects/selection_objects/slices.py index fd510737f30..bdfe9dbe105 100644 --- a/yt/data_objects/selection_objects/slices.py +++ b/yt/data_objects/selection_objects/slices.py @@ -1,3 +1,5 @@ +from typing import Literal + import numpy as np from yt.data_objects.selection_objects.data_selection_objects import ( @@ -313,7 +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): + 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 +347,8 @@ 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 + argument meant for SPH projections Returns ------- diff --git a/yt/geometry/coordinates/cartesian_coordinates.py b/yt/geometry/coordinates/cartesian_coordinates.py index 4fed1977c4e..2b750d87f5c 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 @@ -13,7 +15,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 +173,7 @@ def pixelize( size, antialias=True, periodic=True, + pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave", *, return_mask=False, ): @@ -177,6 +181,16 @@ def pixelize( Method for pixelizing datasets in preparation for two-dimensional image plots. Relies on several sampling routines written in cython + + 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 + 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 +250,18 @@ 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 +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 + 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 +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 + # 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() @@ -667,7 +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): + 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/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..9ec03d998e7 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..befe5347427 100644 --- a/yt/geometry/coordinates/geographic_coordinates.py +++ b/yt/geometry/coordinates/geographic_coordinates.py @@ -1,3 +1,5 @@ +from typing import Literal + import numpy as np import unyt @@ -228,9 +230,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..1a52c304f0e 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": diff --git a/yt/geometry/coordinates/tests/test_sph_pixelization.py b/yt/geometry/coordinates/tests/test_sph_pixelization.py index 1d1c82be5c8..95bf5bc7c5f 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")) @@ -163,6 +165,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 +208,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..27e9953c621 100644 --- a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py +++ b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py @@ -4,6 +4,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, @@ -11,6 +12,7 @@ fake_random_sph_ds, fake_sph_flexible_grid_ds, integrate_kernel, + resreduce_image, ) @@ -112,6 +114,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: @@ -154,12 +157,230 @@ 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(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) +@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]) +def test_sph_proj_pixelave_alongaxes( + axis: int, + 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, + ) + 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 + 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 + + # we don't actually want a plot, it's just a straightforward, + # common way to get an frb / image array + testprj = yt.ProjectionPlot( + ds, + axis, + ("gas", "density"), + width=(2.5, "cm"), + weight_field=None, + buff_size=buff_size, + center=center, + data_source=source, + pixelmeaning="pixelave", + ) + 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", + ) + 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_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_image(_baseimg, buff_size) + _baseimg_wtd = baseprj_wtd.frb.data[("gas", "density")] + _divimg = np.copy(baseimg.v) + _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}, " + 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) + + +def test_massconservation_pixave(): + bbox = np.array([[0.0, 3.0], [0.0, 3.0], [0.0, 3.0]]) + pz = 1.5 + periodic = True + 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.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: + 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.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"), + "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.0, "cm"), + resolution=outsize, + 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.0 + mass_out = np.sum(out.v) * periods**2 / np.prod(outsize) + 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} " + ) + 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 5a1c82f1e2f..f76e08935be 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 @@ -124,6 +125,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,17 +139,20 @@ def integrate_kernel( impact parameter hsml: smoothing length [same units as impact parameter] + nsample: + number of points to sample for the integral 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 + 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 - 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)) @@ -155,6 +160,45 @@ 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]) diff --git a/yt/utilities/answer_testing/framework.py b/yt/utilities/answer_testing/framework.py index 6abd1dc366d..29977c3e21b 100644 --- a/yt/utilities/answer_testing/framework.py +++ b/yt/utilities/answer_testing/framework.py @@ -628,8 +628,13 @@ 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 a6a9f23d9c7..f9957dcbff5 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, @@ -1307,7 +1311,8 @@ def pixelize_sph_kernel_projection( # 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: @@ -1317,11 +1322,426 @@ def pixelize_sph_kernel_projection( free(local_buff) free(xiterv) free(yiterv) + free(ziterv) free(xiter) free(yiter) + free(ziter) 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, + 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 + 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 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, 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 + 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 + 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 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 + 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 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 si in range(nsmin): + xnorm = -1.0 + 2.0 * (si + 0.5) * insmin + xnorm = xnorm * xnorm + for sj in range(nsmin): + ynorm = -1.0 + 2.0 * (sj + 0.5) * insmin + q2 = xnorm + ynorm * ynorm + if q2 >= 1.: + kern1dint[si * nsmin + sj] = 0.0 + else: + kern1dint[si * nsmin + sj] = itab.interpolate(q2) + # 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 * 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] + for sjj in range(sj, nsmin): + 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] + for sjj in range(sj, nsmin): + kern2by2[sci + 3] += kern1dint[sii * nsmin + 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 + # (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) + 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 + + # 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 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? + # 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. 2x2) + if hsml[j] < 0.5 * dx and hsml[j] < 0.5 * dy: + # lower left corner indices + # 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): + # surface density + # = density * av. length + # = density * volume / area + local_buff[x0 * ysize + y0] += ( + 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 + else: # use pre-calculated kernel values + # pre-calculated grid has size + # (nsmin + 1, nsmin + 1, 2, 2) + # bin edge values -1 -- 1 + # in units: + # (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? + 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() for input > 0, we need to round + # down negative numbers for periodic cases) + xi = math.floor( + 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 + 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): + 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: + # catch on the next + # periodic loop + continue + kv = kern2by2[ci + 2 * xxi + yyi] + if kv == 0.: + # don't set mask to 1 + continue + local_buff[(x0 + xxi) * ysize + + (y0 + yyi)] += ( + prefactor_j * kernnorm_j * kv + ) + mask[x0 + xxi, y0 + yyi] = 1 + + else: + # 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 = 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 = 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) + + nsubx = ( + ( nsmin / + nx) + + 1) + insubx = 1.0 / ( nsubx) + nsuby = ( + ( nsmin / + ny) + + 1) + insuby = 1.0 / ( 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.0: continue + local_buff[xi * ysize + yi] += ( + prefactor_j * insuby * insubx + * itab.interpolate(q_ij2) + ) + mask[xi, yi] = 1 + + with gil: + for sxi in range(xsize): + for syi in range(ysize): + buff[sxi, syi] += local_buff[sxi * ysize + syi] + # free memory in (hopefully?) private variables assigned in + # each thread + free(local_buff) + free(xiterv) + free(yiterv) + free(ziterv) + free(xiter) + free(yiter) + free(ziter) + # free memory in shared variable + free(kern2by2) + return mask + @cython.boundscheck(False) @cython.wraparound(False) def interpolate_sph_positions_gather(np.float64_t[:] buff, @@ -2112,7 +2532,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.): @@ -2139,21 +2560,32 @@ 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 + ) + else: + raise ValueError(f"no pixelmeaning option '{pixelmeaning}'") # 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..ee1ff4771e5 100644 --- a/yt/visualization/fits_image.py +++ b/yt/visualization/fits_image.py @@ -2,6 +2,7 @@ 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 @@ -111,6 +112,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: + "pixelave": 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 +880,15 @@ 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: Literal["pixelave", "pencilbeam"] = "pixelave", ): if width is None: width = ds.domain_width[axis_wcs[axis]] @@ -910,9 +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]) + 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 +1164,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: + "pixelave": 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 +1188,7 @@ def __init__( *, origin="domain", moment=1, + pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave", **kwargs, ): fields = list(iter_fields(fields)) @@ -1170,6 +1206,7 @@ def __init__( width, length_unit, origin=origin, + pixelmeaning=pixelmeaning, ) super().__init__(frb, fields=fields, length_unit=lunit, wcs=w) @@ -1378,6 +1415,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: + "pixelave": 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 +1442,7 @@ def __init__( field_parameters=None, data_source=None, north_vector=None, + pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave", ): fields = list(iter_fields(fields)) center, dcenter = ds.coordinates.sanitize_center(center, None) @@ -1427,6 +1474,7 @@ def __init__( width, length_unit, origin="image", + pixelmeaning=pixelmeaning, ) super().__init__(frb, fields=fields, length_unit=lunit, wcs=w) @@ -1603,6 +1651,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: + "pixelave": 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 +1676,7 @@ def __init__( depth=(1.0, "unitary"), method="integrate", length_unit=None, + pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave", *, moment=1, ): @@ -1652,6 +1710,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..5a81deb9cd6 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 TYPE_CHECKING, Literal import numpy as np @@ -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: + "pixelave": 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: Literal["pixelave", "pencilbeam"] = "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=self.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) @@ -702,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__( @@ -711,6 +730,7 @@ def __init__( buff_size, antialias=True, periodic=False, + pixelmeaning="pixelave", *, filters=None, ): diff --git a/yt/visualization/plot_window.py b/yt/visualization/plot_window.py index 9a1bbc9e8e2..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 TYPE_CHECKING, Union +from typing import TYPE_CHECKING, Literal, Union import matplotlib import numpy as np @@ -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: + "pixelave": 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: Literal["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. @@ -764,6 +776,28 @@ def set_antialias(self, aa): """ self.antialias = aa + @invalidate_data + def set_pixelmeaning( + self, + pixelmeaning: Literal["pixelave", "pencilbeam"] = "pencilbeam", + ): + """ + Change the SPH surface density calculation approach + + parameters + ---------- + 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 + 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): """Sets a new buffer size for the fixed resolution buffer @@ -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: + "pixelave": 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: Literal["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: + "pixelave": 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: Literal["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) diff --git a/yt/visualization/tests/test_offaxisprojection_pytestonly.py b/yt/visualization/tests/test_offaxisprojection_pytestonly.py index 246e5954c84..8e773375da6 100644 --- a/yt/visualization/tests/test_offaxisprojection_pytestonly.py +++ b/yt/visualization/tests/test_offaxisprojection_pytestonly.py @@ -11,6 +11,7 @@ fake_sph_grid_ds, integrate_kernel, requires_module_pytest, + resreduce_image, ) from yt.visualization.api import OffAxisProjectionPlot, ProjectionPlot @@ -128,6 +129,7 @@ def makemasses(i, j, k): data_source=source, north_vector=northvector, depth=depth, + pixelmeaning="pencilbeam", ) img = prj.frb.data["gas", "density"] if weighted: @@ -171,6 +173,141 @@ def makemasses(i, j, k): assert_rel_equal(expected_out, img.v, 4) +@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]) +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.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}, " + 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 = [ diff --git a/yt/visualization/volume_rendering/off_axis_projection.py b/yt/visualization/volume_rendering/off_axis_projection.py index 5de55f2fbbc..e7cc7551142 100644 --- a/yt/visualization/volume_rendering/off_axis_projection.py +++ b/yt/visualization/volume_rendering/off_axis_projection.py @@ -1,3 +1,5 @@ +from typing import Literal + import numpy as np from yt.data_objects.api import ImageArray @@ -33,6 +35,7 @@ def off_axis_projection( depth=None, num_threads=1, method="integrate", + pixelmeaning: Literal["pixelave", "pencilbeam"] = "pixelave", ): r"""Project through a dataset, off-axis, and return the image plane. @@ -102,6 +105,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: + "pixelave": 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 +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") data_source = data_source_or_all(data_source) @@ -231,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] @@ -261,6 +276,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 +319,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 +341,7 @@ def off_axis_projection( north, depth=depth, kernel_name=kernel_name, + pixelmeaning=pixelmeaning, ) normalization_2d_utility(buf, weight_buff) 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..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,8 +1,12 @@ import numpy as np +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 +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 @@ -10,11 +14,13 @@ 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 +42,26 @@ 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") - ) - pixelize_sph_kernel_projection( - buf2, mask, px, py, pz, hsml, mass, density, quantity_to_smooth, bounds + ds, + center, + normal_vector, + width, + resolution, + ("gas", "density"), + pixelmeaning=pixelmeaning, ) + 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 +92,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 +133,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 +165,20 @@ 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 +212,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 +243,20 @@ 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 +277,20 @@ 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 +308,13 @@ 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)