WIP: add pixel-averaging option for SPH projections

[nastasha-w]

PR Summary

This PR will add an option (to be set as default) for SPH projections. The current implementation integrates the specified field over a line of sight through the center of each pixel in the grid the simulation is projected onto. This PR will add the option to average the specified field over the pixel area instead. The output units will be the same: input field units per unit pixel area (non-weighted), or input field units (weighted). This should close issue #4993.

The implementation is based on Borrow & Kelly (2021).

This PR should add this projection option for axis-aligned and off-axis projections.

PR Checklist

• New features are documented, with docstrings and narrative docs
• Adds tests: unit tests for the same cases as the old/current SPH projections.
• Adds test: mass conservation for different smoothing length and pixel size combinations.
• Qualitative test: convergence between SPH projection methods as pixel size approaches 0.

[nastasha-w]

pre-commit.ci autofix

[nastasha-w]

The math behind the pixel-averaging SPH projection implementation

I thought this might be useful to get some early feedback, and as a first draft for the documentation.

In SPH simulations, a field $q(\vec{x})$, which depends on position $\vec{x}$, is sampled by the SPH particles $i$ such that

$$q(\vec{x}) \approx \sum_{i} \frac{m_i}{\rho_i} q_i W\left(\vec{x} - \vec{x_i}, h_i\right),$$

where $m_i$, $\rho_i$, $\vec{x_i}$, and $h_i$ are the mass, density, position, and smoothing length of particle $i$, respectively. $W$ is the kernel function used in a given simulation. Its value is $>0$ only for arguments $<1$. We can write this kernel function as

$$W\left(\left| \vec{x} - \vec{x_i}\right|, h_i \right) = \frac{1}{{h_i}^3} U\left(\frac{\left| \vec{x} - \vec{x_i}\right|}{h_i}\right),$$

where U is a function with a dimensionless output (and input). The volume integral $\int_V U(x) dx = 1$ (see e.g., #4781). Note the input unit change between $W$ and $U$. Since both sides should integrate to the same quantity,

$$W(x, y, z, h_i) \, dx \, dy \, dz = U\left(\sqrt{{x'}^2 + {y'}^2 + {z'}^2}\right) dx' dy' dz',$$

where the cartesian coordinates $(x, y, z) = h_i (x', y', z')$, with the particle centered on $(x, y, z) = (0, 0, 0)$.

For the pixel-averaging projections, we want to integrate the field $q$ over the rectangular prism defined by each pixel and the projection depth, then divide by the pixel area $A_{\text{pix}}$. We can split this into an integral along the line of sight direction ($\int dl$) and across the pixel area ($\int dA$). As with the pixel center line-of-sight integrals, we assume the SPH particle lies entirely within the volume in the line-of-sight direction. (However, we restrict the particles to project to those whose centers lie strictly within projection depth region, i.e., without a smoothing length margin, meaning we do not generally bias $\Sigma_q$ high. Projection regions that overlap only at a common edge along the line of sight will only share a particle if it lies exactly on the projection region boundary.) This assumption is generally reasonable, since a projection with a depth comparable to typical SPH particle sizes would measure something more akin to a Slice than the line-of-sight aggregate quantities (or weighted averages) that Projections are meant to capture. However, if the goal is to capture features on the finest resolved scales of a simulation, the pixel size will be comparable to the size of the smaller SPH particles in the volume. For a typical astrophysical simulation analyzed with yt, density contrasts are large, meaning that the larger SPH particles in the volume will intersect many pixels in the projection grid.

The assumption that integrals in the line-of-sight direction are limited by the SPH particle size, rather than by the projection domain, allows us to speed up the integral calculations. We get the line-of-sight integrals by interpolating tabulated values for different impact parameters, since $\int_{-\infty}^{\infty} U(\sqrt{b^2 + {l'}^2}) dl'$ depends only on $b$, the impact parameter of the line of sight relative to the particle center, divided by the smoothing length (see e.g., #4781).

For the projected quantity $\Sigma_q$ in one pixel of the output grid, we then get:

$$\displaylines{\Sigma_q = \frac{1}{A_{\text{pix}}} \int dA \int dl \, q(\vec{x})\\\ \approx \frac{1}{A_{\text{pix}}} \int dA \int dl \, \sum_{i} \frac{m_i}{\rho_i} q_i W\left(\frac{\left| \vec{x} - \vec{x_i}\right|}{h_i}\right)\\\ = \frac{1}{A_{\text{pix}}} \sum_{i} \frac{m_i}{\rho_i {h_i}^3} q_i \int dA \int dl' \, U\left(\frac{\left| \vec{x} - \vec{x_i}\right|}{h_i}\right) \text{ where } l' = \frac{l}{h_i}\\\ = \frac{1}{A_{\text{pix}}} \sum_{i} \frac{m_i}{\rho_i {h_i}^2} q_i \int dA \int dl' \, U\left(\sqrt{{b_i}^2 + {l'}^2}\right) \text{ where } b_i \text{ is the impact parameter of the sightline in units of } h_i. }$$

We can get the 1-D line-of-sight integrals $\int dl' U(\sqrt{{b_i}^2 + {l'}^2})$ by linearly interpolating a table of the integral values as a function of $b_i$, ranging from 0 to 1.

For the integrals $\int dA$, we take a few different approaches, depending on the relative size of the SPH particle and the output grid pixels. These strategies are taken from Borrow & Kelly (2021).

If a particle $i$ falls entirely within one pixel,

we take $dA = dxdy = {h_i}^2 dx'dy'$, where $x', y'$ are the smoothing-length-normalized positions in the output grid
axis directions. The area we are integrating over is limited by the extent of the kernel, not that of the pixel. We get

$$\displaylines{\frac{1}{A_{\text{pix}}} \frac{m_i}{\rho_i {h_i}^2} q_i \int dA \int dl' \, U\left(\sqrt{{b_i}^2 + {l'}^2}\right) \\\ = \frac{1}{A_{\text{pix}}} \frac{m_i}{\rho_i {h_i}^2} q_i {h_i}^2 \int dx' dy' \int dl' \, U\left(\sqrt{{x'}^2 + {y'}^2 + {l'}^2}\right)\\\ = \frac{m_i}{\rho_i A_{\text{pix}}} q_i, \text{since } \int dV \, U(\vec{x'}) = 1.}$$

In other words, we just multiply the quantity $q_i$ by the particle volume and divide by the pixel area. Qualitatively, $\frac{m_i}{\rho_i} q_i$ is the volume integral of the field $q$ over particle $i$, e.g., the mass of $i$ if $q$ is the density.
The contribution of a small particle to each pixel is therefore also what we'd use to make a surface density map by histogramming a particle-integrated property, e.g., total mass per unit pixel area. This agreement forms a good sanity check to this approach.

If a particle $i$ falls within multiple pixels,

the area we are integrating over is generally limited by both the extent of the kernel and that of the pixel. In this case, we approximate the area integral with a finite sum. For this, we define 'subpixels', labelled by index $j$, and each having an area $A_j$. $A_{\text{pix}}$ is still the area of each pixel in the output grid. We evaluate the line integral for a line of sight through the center of each subpixel. For this, we get

$$\displaylines{\frac{1}{A_{\text{pix}}} \frac{m_i}{\rho_i {h_i}^2} q_i \int dA \int dl' \, U\left(\sqrt{{b_i}^2 + {l'}^2}\right) \\\ \approx \frac{1}{A_{\text{pix}}} \frac{m_i}{\rho_i {h_i}^2} q_i \sum_{\text{subpixels } j} A_{j} \int dl' \, U\left(\sqrt{{b_{i, j}}^2 + {l'}^2}\right) \\\ = \frac{m_i}{\rho_i {h_i}^2} q_i \sum_{\text{subpixels } j} \frac{A_{j}}{A_{\text{pix}}} \int dl' \, U\left(\sqrt{{b_{i, j}}^2 + {l'}^2}\right) \\\ = \frac{m_i}{\rho_i {h_i}^2} q_i \sum_{\text{subpixels } j} \frac{1}{n_{\text{sub}, x} n_{\text{sub}, y}} \int dl' \, U\left(\sqrt{{b_{i, j}}^2 + {l'}^2}\right),}$$

where $n_{\text{sub}, x/y}$ is the number of subpixels per output grid pixel in the output grid axis directions $x$ and $y$.
Following Borrow & Kelly (2021), we set a minimum number of subpixels to sample the kernel with. For a given particle $i$, each output grid pixel is subsampled with the same number of pixels, such that the total meets the minimum requirement.

Qualitatively, $\frac{m_i}{\rho_i} q_i$ is the integral contribution of particle $i$ to the field $q$. The rest of the equation should then encapsulate the overlapped particle volume fraction, divided by the pixel area. We can see the subpixel sum here as an (approximate) area-weighted average over the pixel of the dimensionless particle path length. Using $dA = h_i dA'$ and denoting the average with angle brackets, we get

$$\displaylines{ \frac{1}{{h_j}^2} \sum_{\text{subpixels } j} \frac{1}{n_{\text{sub}, x} n_{\text{sub}, y}} \int dl' \, U\left(\sqrt{{b_{i, j}}^2 + {l'}^2}\right) \\ \approx \frac{1}{{h_j}^2} \left\langle \int dl' \, U\left(\sqrt{{b_{i, j}}^2 + {l'}^2}\right) \right\rangle_{\text{area}} \\\ = \frac{1}{A_{\text{pix}}} \int dA' \int dl' \, U\left(\sqrt{{b_{i, j}}^2 + {(l / h_i)}^2}\right),}$$

which works out since the dimensionless kernel integral over the pixel is indeed the fraction of particle $i$'s volume contained within that pixel (or rather, the rectangular prism it defines).

For small particles, which overlap at most 2x2 output grid pixels,

we use a time-saving optimization from Borrow & Kelly (2021). If we define $A = {h_i}^2 A'$, we can write the equation for multiple-pixel-overlaps case as

$$\displaylines{\frac{1}{A_{\text{pix}}} \frac{m_i}{\rho_i {h_i}^2} q_i \int dA \int dl' \, U\left(\sqrt{{b_i}^2 + {l'}^2}\right) \\\ = \frac{m_i}{\rho_i {h_i}^2} q_i \int dA \frac{1}{A_{\text{pix}}} \int dl' \, U\left(\sqrt{{b_{i}}^2 + {l'}^2}\right) \\\ = \frac{m_i}{\rho_i {h_i}^2} q_i \int dA' \frac{{h_j}^2}{A_{\text{pix}}} \int dl' \, U\left(\sqrt{{b_{i}}^2 + {l'}^2}\right). \\\ = \frac{m_i}{\rho_i A_{\text{pix}}} q_i \int dA' \int dl' \, U\left(\sqrt{{b_{i}}^2 + {l'}^2}\right).}$$

For a given particle $i$, the dimensionless integral $\int dA' \int dl' , U\left(\sqrt{{b_{i}}^2 + {l'}^2}\right)$ is the fraction of the particle volume enclosed with the output grid pixel. If the number of pixels a particle overlaps is small, as in $2\times2$ at most, we can tabulate the values in each of the four pixels as a function of the $h_i$-normalized offset (in both output grid axis directions) between the $2\times2$ grid center and the particle position.

We tabulate these values for normalized offsets in the range [-1, 1] in both axis directions. We approximate the integral $\int dA'$ by subsampling the 1-D kernel integrals with the minimum number of subpixels in each dimension. This leaves us with a (minimum + 1, minimum + 1) grid of normalized offsets. For these small particles, we calculate their contribution to each of the 4 pixels they might overlap by using the tabulated volume fractions for the nearest normalized offset position in each dimension. This is similar to the 1-pixel-overlap case, except now we use the kernel integral over each pixel to divide the integral contribution of the SPH particle $\frac{m_i}{\rho_i} q_i$ over the output grid pixels.

(This could maybe be improved by interpolating between tabulated values instead? I'm not sure how much the exact offset differences matter compared to the finite sum approximation though.)

[nastasha-w]

pre-commit.ci autofix

[chummels]

I don't yet have time to review this, but I just wanted to comment. Thank you for this incredibly important work on getting the ProjectionPlot functionality upgraded. I think this will have a positive impact on the work of hundreds to thousands of people. Looking forward to discussing further soon!

[nastasha-w]

pre-commit.ci autofix

[nastasha-w]

pre-commit.ci autofix

[nastasha-w]

pre-commit.ci autofix

[nastasha-w]

pre-commit.ci autofix

[nastasha-w]

pre-commit.ci autofix

[nastasha-w]

I was hitting a bit of a wall with testing: the baseline and yt image projections looked similar, but there were differences larger than I'd typically be comfortable with (~25%). Part of the issue was likely that these were comparisons to brute-forces calculations, integrating each kernel over each output grid pixel by subsampling the pixel by a factor of 50 in each direction, and using 50 samples along the line-of-sight direction as well. This is pretty slow though, and in some cases a worse sampling rate than in the yt method new to this PR: pixel-averaging.

I switched to testing against the pencil-beam in a smaller test; this approach gets me better subsampling faster, by simply projecting onto a higher-resolution grid than the testing target, then averaging the values in blocks of those pixels. Based on this testing, I think getting the pixel values right to high precision with this method is maybe just kind of hard...

For this first test, I downsampled pixel values from grids of increasing size, then calculated the total projected mass by adding up the pixel values and multiplying by the area per pixel. The x axis shows the subsampling factor along each grid axis. These projections are the ones shown in the next plot: a simple test dataset with particles of different sizes, and different placements relative to the grid.

The main conclusion here is that the pixel averaging method does a pretty good job of conserving mass: results are accurate to about 1e-4 fractionally even at the lowest resolution.

However, the values of individual pixels do not seem to converge as rapidly. The image below is messy, as I tend to make them, but convergence does require many images. The top two rows show the pixel averaging method: first, projections of the same dataset, subsampled by the factors listed above, which increase from left to right. Black dashed lines indicate the extent of the seven SPH particles in the test dataset; this dataset was periodic in the x-axis direction, but not in the y-axis direction. The second row shows the difference between each subsampled image and the image at the right, with the highest subsampling factor. The labels at the top of the different images show the maximum absolute value of the relative difference (pixel value) / (pixel value, max. subsampling) - 1. The lowest-resolution images here show ~10%-level differences with the ones subsampled at higher resolution, which is not as good as I had hoped.

For comparison, the next two rows show the same convergence test, but for the pencil-beam projections. Those definitely don't do as well, with large differences at the lowest subsampling resolution. However, even at a 64x64 subsampling rate, there are still ~1% differences here. The smallest smoothing length in the test data is 0.2, so this subsampling rate corresponds to 0.4 * 64 / 1 = 26 pixels per dimension across the kernel. The minimum in the pixel-averaging method's backend is 32 pixels per dimension.

The bottom row compares the pencil-beam and pixel-averaging projections at different subsampling resolutions. The images for 128x128 subsampling pixels are identical, because the subsampling pixels here are small enough that the smallest particles in the test dataset are resolved at the minimum level. This also holds for higher resolutions I tried.

tldr:

• I want to switch my baseline image in the tests to one that using the existing pencil-beam projections, since those run faster than my initial tests, and should do the exact same thing.
• Getting converged pixel values is hard. Using a larger subsampling resolution would probably help, but of course, this comes with a speed trade-off. It also makes testing harder, in the sense that setting wider tolerances might also let some genuine problems slip through the cracks.

@JBorrow @neutrinoceros @chummels , any advice on what sort level of convergence we want for these pixel values? As a point of reference, I think the volume integrals are limited to a 1e-4 -- 1e-5 in 1 level accuracy by the kernel line integral approximation used in both projection methods.

[nastasha-w]

pre-commit.ci autofix

[nastasha-w]

I've taken a look at convergence of the methods with pixel size; here I showed this for a test dataset, but I did also take a look at some actual simulation data. With grids up to 800x800 pixels, it seems like cosmological volumes and zoom simulations projected over their whole parent boxes, don't converge well (yet). That's not unexpected given that the pixels are still larger than many resolution elements here, but it is worth keeping in mind.

This gizmo Milky-Way disk dataset shows some convergence, even though it's no quite there yet at 800x800 pixels.

The Arepo cluster projections show pretty good convergence, I think:

[nastasha-w]

I still need to write the narrative docs to go with this PR, but I would appreciate if people could start taking a look at the implementation part of this PR. I've run the tests and done the further tests that I planned to for this, but the pixel values don't agree as well with those obtained from very high pixel subsampling rates as I had hoped. @JBorrow and @AshKelly, I would particularly appreciate feedback from you on whether this is a bug that needs fixing, or whether these values just converge slowly.

I'm hoping I can resolve the type check failure by rebasing after #5137 is done and merged. [grammar edit]

[chrishavlin]

One quick comment not related to the details of the implementation:

This PR will add an option (to be set as default)

I could be wrong, but this feels like a pretty big change to drop on folks by default in the next feature release (4.5). IMO it'd be a smoother transition to not enable this by default initially: i.e., I'd prefer that next release 4.5 introduce and advertise the functionality and then the subsequent release (4.6 or 5.0?) switch it to default.

[nastasha-w]

Good point! It would indeed be better to set the default to "pencilbeam" for now.

[JBorrow]

I've taken a brief look at the code and I can't immediately see if there is something wrong.

However, if you are using the subsampling strategy that we demonstrated in Borrow & Kelly (2021), there should be no 'convergence' - you should get the same answer irrespective of your pixel size. I believe there must hence be a bug in the code. You can see in our Fig 9 that for a similar particle setup to you we are able to constrain reconstruction errors to approximately less than 0.01% in all cases.

You could test this yourself against our reference implementation included in swiftsimio. You can find information on the lower level API here: https://swiftsimio.readthedocs.io/en/latest/visualisation/projection.html#lower-level-api

To use that, you would do something like:

# Rescale to 0, 1... then
from swiftsimio.visualisation.projection_backends.subsampled import scatter

parameters = {
  "x": x,
  "y": y,
  "m": masses,
  "h": hsml
}

output = {
  x: scatter(**parameters, res=x) for x in (100, 200, 400, 800)
}

Also consider comparing against subsampled_extreme.

For what it's worth, I believe we use a different strategy than you to determine how to subsample pixels. We define a minimum number of kernel evaluations (one-sided), this is by default 16. So if our kernel overlaps less than 32x32 pixels, we oversample those pixels. If it overlaps 16x16 < x < 32x32, each pixel is evaluated 4 (2x2) times. For 8x8 < 16x16 we do each pixel 16 (4x4) times, and so forth: https://swiftsimio.readthedocs.io/en/latest/_modules/swiftsimio/visualisation/projection_backends/subsampled.html#scatter

[nastasha-w]

@JBorrow Thanks! I was doing 1 sample per pixel if 32x32 pixels were overlapped, I'll see if a higher minimum number of evaluations works better. If not, I'll have to see if there's a deeper problem.

[nastasha-w]

Note: I rebased this branch to main after Chris Havelin's type annotation fixes in the SPH projection tests were merged

[nastasha-w]

pre-commit.ci autofix