Exact interpolation is not correct at low grid resolutions

[JBorrow]

I was attempting to use sarracen to re-create some figures from https://arxiv.org/abs/2106.05281 as part of my JOSS review and I am running into some issues with the 'exact' backend.

The following image shows a gradient in density that is purposefully misaligned with the background grid. Your base setup does excellently on this test, producing less than 0.1% error throughout. But the 'exact=True' setup causes some issues:

gradients_2d.pdf

It looks even worse when looking at the convergence (though it does converge super fast). I suspect there may be a bug?

test_convergence.pdf

For what it's worth, I could never get the method to work when implementing in swiftismio, hence why we went the route of subsampling.

I am also interested in how you manage to keep the error low when the grid is poorly resolved; is the numerical method documented somewhere?

Code (apologies for the poor quality; note subsampled <-> exact, original <-> fast):

"""
Sets up a gradient in particles, and visualises it with
different grid spacings.
"""

import numpy as np
import matplotlib.pyplot as plt

import sarracen
from numba import njit

from matplotlib.colors import LogNorm, Normalize
from swiftsimio.visualisation.projection_backends.fast import kernel_gamma
from swiftsimio.visualisation.projection_backends.fast import scatter as scatter_fast
from swiftsimio.visualisation.projection_backends.subsampled import (
    scatter as scatter_sub,
)
from swiftsimio.visualisation.projection_backends.histogram import (
    scatter as scatter_hist,
)
from swiftascmaps import lover

image_width = 7.3  # inches

# Internal units, from 0 -> 1 on edges
particle_spacing_lower = 1e-3
particle_spacing_upper = 10 ** (-2.5)
particle_spacing_gradient = particle_spacing_upper - particle_spacing_lower

@njit()
def get_x_y_h():
    dx = particle_spacing_lower
    x = 0.0
    y = 0.0
    particle_x = []
    particle_y = []
    smoothings = []

    while x < np.sqrt(2) + 0.3:
        y = 0.0

        while y < np.sqrt(2) + 0.3:
            particle_x.append(x)
            particle_y.append(y)

            y += dx

            smoothings.append(dx)  # * kernel_gamma)

        x += dx
        dx = particle_spacing_gradient * x + particle_spacing_lower

    particle_x = np.array(particle_x)  # - 0.5 * np.sqrt(2)
    particle_y = np.array(particle_y)  # - 0.5 * np.sqrt(2)
    smoothings = np.array(smoothings)

    #    return particle_x, particle_y, smoothings

    particle_x -= 0.5 * np.sqrt(2)
    particle_y -= 0.5 * np.sqrt(2)

    rot_by = np.pi / 4
    rot_x = particle_x * np.cos(rot_by) - particle_y * np.sin(rot_by) + 0.5 * np.sqrt(2)
    rot_y = particle_x * np.sin(rot_by) + particle_y * np.cos(rot_by) + 0.5 * np.sqrt(2)

    # Yes, I do hate geometry.
    offset_left_corner = 1 / (2 + np.sqrt(2)) + 0.05

    return rot_x - offset_left_corner, rot_y - offset_left_corner, smoothings

particle_x, particle_y, smoothings = get_x_y_h()

fast = {}
sub = {}
hist = {}

resolutions = [64, 128, 256]

for pixels in resolutions:
    # Build sarracen data frame
    data = sarracen.SarracenDataFrame(
        data=dict(
            x=particle_x,
            y=particle_y,
            #z=np.zeros_like(particle_x),
            m=np.ones_like(particle_x),
            h=smoothings,
            ndim=2,
        )
    )
    data.params = dict(hfact = 1.2, mass=1)
    data.calc_density()

    fast[pixels] = data.sph_interpolate("rho", x_pixels=pixels, y_pixels=pixels, xlim=[0, 1], ylim=[0,1])[
        1:-1, 1:-1
    ]

    sub[pixels] = data.sph_interpolate(
        "rho", x_pixels=pixels, y_pixels=pixels, exact=True, xlim=[0,1], ylim=[0, 1]
    )[1:-1, 1:-1]

    hist[pixels] = scatter_hist(
        x=particle_x,
        y=particle_y,
        m=np.ones_like(particle_x),
        h=smoothings,
        res=pixels + 2,
    )[1:-1, 1:-1]

# min = np.min([np.min(v) for v in sub.values()])
# max = np.max([np.max(v) for v in sub.values()])
min = np.min(sub[resolutions[-1]][np.nonzero(sub[resolutions[-1]])])
max = np.max(sub[resolutions[-1]])

fig, axes = plt.subplots(
    3,
    3,
    sharex=True,
    sharey=True,
    squeeze=True,
    figsize=(image_width, image_width + 0.57),
    constrained_layout=True
)

for row, style in zip(axes, [hist, fast, sub]):
    for ax, dpi in zip(row, resolutions):
        ax.tick_params(
            axis="both",
            which="both",
            **{k: False for k in ["bottom", "top", "left", "right"]},
            **{f"label{k}": False for k in ["bottom", "top", "left", "right"]},
        )

        img = ax.imshow(
            np.log10(style[dpi]),
            origin="lower",
            norm=Normalize(vmin=np.log10(min), vmax=np.log10(max), clip=True),
            cmap=lover,
            extent=[0, 1, 0, 1],
            interpolation="nearest",
        )

for ax, res in zip(axes[0], resolutions):
    ax.text(0.5, 0.975, f"${res}^2$", color="black", ha="center", va="top")

for ax, method in zip(axes[:, 0], ["Histogram", "Base", "Exact"]):
    ax.text(0.025, 0.5, method, color="black", rotation=90, ha="left", va="center")

cb = fig.colorbar(
    img,
    label="Projected Density",
    orientation="horizontal",
    ticks=[-1],
    ax=axes.ravel().tolist(),
    pad=0.005,
)
cb.set_ticks([])
fig.savefig("gradients_2d.pdf", dpi=300)

[JBorrow]

https://github.com/openjournals/joss-reviews/issues/5263

[ttricco]

@JBorrow This is awesome. I'm really intrigued by your analysis. I'm going to dig deeper into this, as I have an idea of what may be causing the low resolution results you've found.

The numerical method is documented insofar as it is in Petkova, Laibe & Bonnell (2018). That gives the general form of the equations for polyhedra (Voronoi cells). We are using this method applied to regular cubes (cubic grid cells). Petkova and Price derived the specific form of the equations for that case and implemented into Splash. To my knowledge, there is no documentation on that set of equations specifically.

[JBorrow]

So something that always concerned me about the Petkova paper is that they show that there are still errors when using their exact framework (though they are smaller). At some level, given that it is a static problem, there should be no need to converge (we know the exact answer already).

Also - I was actually referring to the numerical method for your 'fast' case. Does this do the typical splash thing of setting a minimal smoothing length for all gas particles based upon the grid spacing?

[ttricco]

The normal case is described in more detail in the API: https://sarracen.readthedocs.io/en/latest/api/sarracen.SarracenDataFrame.render.html#sarracen.SarracenDataFrame.render

The normal case does not set a minimum smoothing length. Something to consider -- as is your blitting approach from the paper you linked (thanks, was not aware of this).

[JBorrow]

So in Splash, as far as I am aware, h is set as max(h, pixel_size / 2): https://ui.adsabs.harvard.edu/abs/2007PASA...24..159P/abstract (eqn 13). I am somewhat surprised then that you are able to get such accurate interpolations, that's very cool.

One thing that I have noticed about both your implementations is that if a kernel has no overlap with any pixels, then it doesn't get included. See this pathalogical case:

pathalogical_cases.pdf

Code:

"""
Plots of a couple of pathalogical cases.
"""

import matplotlib.pyplot as plt
import numpy as np

#plt.style.use("../style.mplstyle")

from matplotlib.patches import Circle
from matplotlib.colors import LogNorm

from swiftsimio.visualisation.projection_backends import backends
from swiftsimio.visualisation.projection import kernel_gamma

from swiftascmaps import nineteen_eighty_nine

import sarracen

# Sarracen uses cubic spline
kernel_gamma = 1.778002

data_1 = sarracen.SarracenDataFrame(
    data=dict(
        x=np.array([0.525]),
        y=np.array([0.525]),
        #z=np.zeros_like(particle_x),
        m=np.array([1.0]),
        h=np.array([0.03]),
        ndim=2,
    )
)
data_1.params = dict(hfact = 1.0, mass=1)
data_1.calc_density()

data_2 = sarracen.SarracenDataFrame(
    data=dict(
        x=np.array([0.6]),
        y=np.array([0.6]),
        #z=np.zeros_like(particle_x),
        m=np.array([1.0]),
        h=np.array([0.2]),
        ndim=2,
    )
)
data_2.params = dict(hfact = 1.0, mass=1)
data_2.calc_density()

def divide_grid_res_by_two(grid):
    return 0.25 * (
        (grid[0::2, 0::2] + grid[1::2, 0::2]) + (grid[0::2, 1::2] + grid[1::2, 1::2])
    )

converged_answer_one = data_1.sph_interpolate("rho", x_pixels=1024, y_pixels=1024, xlim=[0, 1], ylim=[0,1])
converged_answer_two = data_2.sph_interpolate("rho", x_pixels=1024, y_pixels=1024, xlim=[0, 1], ylim=[0,1])

plt.imsave("test.png", plt.get_cmap()(plt.Normalize()(converged_answer_two)))

while converged_answer_one.size > 16:
    converged_answer_one = divide_grid_res_by_two(converged_answer_one)

while converged_answer_two.size > 16:
    converged_answer_two = divide_grid_res_by_two(converged_answer_two)

fig, ax = plt.subplots(
    3, 2, sharex=True, sharey=True, squeeze=True, figsize=(3.65, 3.65*1.5 + 0.43), constrained_layout=True
)

ax[0][0].text(0.025, 0.5, "Base", ha="left", va="center", rotation=90, transform=ax[0][0].transAxes)
ax[1][0].text(0.025, 0.5, "Exact", ha="left", va="center", rotation=90, transform=ax[1][0].transAxes)
ax[2][0].text(0.025, 0.5, "True", ha="left", va="center", rotation=90, transform=ax[2][0].transAxes)

imaging_args = dict(
    origin="lower",
    norm=LogNorm(vmin=1e-3, vmax=1e1),
    extent=[0, 1, 0, 1],
    cmap=nineteen_eighty_nine,
)
text_args = dict(x=0.025, y=0.025, ha="left", va="bottom")

def chi_square(o, e):
    return float(((o[e != 0] - e[e != 0]) ** 2 / e[e != 0] ** 2).sum())

results = [
    [
        data_1.sph_interpolate("rho", x_pixels=4, y_pixels=4, xlim=[0, 1], ylim=[0,1]),
        data_1.sph_interpolate("rho", x_pixels=4, y_pixels=4, xlim=[0, 1], ylim=[0,1], exact=True),
        converged_answer_one,
    ],
    [
        data_2.sph_interpolate("rho", x_pixels=4, y_pixels=4, xlim=[0, 1], ylim=[0,1]),
        data_2.sph_interpolate("rho", x_pixels=4, y_pixels=4, xlim=[0, 1], ylim=[0,1], exact=True),
        converged_answer_two,
    ],
]

expected = [
    [converged_answer_one, converged_answer_one, converged_answer_one],
    [converged_answer_two, converged_answer_two, converged_answer_two],
]

for row, row_res, row_exp in zip(ax.T, results, expected):
    for axis, res, exp in zip(row, row_res, row_exp):
        mappable = axis.imshow(res.T, **imaging_args)
        axis.text(**text_args, s=f"$\\chi^2=$\n${chi_square(res, exp):.3f}$")

cb = fig.colorbar(
    mappable,
    label="Projected Density",
    orientation="horizontal",
    ticks=[],
    ax=ax.ravel().tolist(),
    pad=0.005,
)

for data, row in zip([data_1, data_2, data_1, data_2], ax.T):
    for axis in row:
        c = Circle(
            (data.x, data.y),
            data.h * kernel_gamma,
            color="k",
            fill=False,
            linestyle="dashed",
        )
        axis.add_artist(c)
        axis.tick_params(
            axis="both",
            which="both",
            bottom=False,
            top=False,
            left=False,
            right=False,
            labelbottom=False,
            labelleft=False,
        )

cb.set_ticks([])
cb.ax.tick_params(size=0, which="both")

fig.savefig("pathalogical_cases.pdf")

[AndrewHarris709]

My apologies for the long delay! The underlying issue was related to what you observed above, any pixels that did not have the center overlapped by a particular particle had that particular "contribution" excluded. This didn't matter for our fast interpolation method, but completely broke our exact interpolation method at low resolutions. #54 should resolve this.

With #54 included, the two codes you posted now produce the following results:

[JBorrow]

Why does the bottom figure still have missing contributions for the top left case?

[JBorrow]

In that case, there is no overlap with any of the kernel with any of the cell centers. But there is still overlap between the pixel and the cells (as shown perfectly by the exact case!)

[ttricco]

It depends upon what you consider a pixel. For our fast case, the "pixel" is just considered to be the point at the centre of the cell. Only particles whose smoothing lengths overlap with that point will contribute.

Depending upon your pixel resolution, some particles can "fall between the cracks", but I'm not convinced this is a problem per se. Splash implements an option to set the minimum smoothing length based on the grid size so that all particles are included. This can help represent structure below the grid resolution. Though artificially inflating particles lowers the particle resolution, and diverges from the actual resolution of the calculation, so there is a tradeoff.

We've specifically chosen not to use this minimum smoothing length approach in our 'fast' case so that our default interpolation most closely aligns with the standard SPH interpolation. We think this is what most people are familiar with and would naively expect. Adding this minimum smoothing length as an option that the user can select is something we will likely add though.

[JBorrow]

The problem with this technique is that in adaptive applications (e.g. galaxy formation, stellar formation), a massive fraction of the mass can 'fall between the cracks'. If you take all of your particle masses, and sum them up, the result should be the same as the whole grid density grid area (modulo any grid interpolation errors). At present, in your implementation, you are potentially missing out on a huge fraction of the mass (you are actually neglecting almost all the mass where h <~ 0.5 grid_spacing). It also introduces a term where the choice of grid resolution directly corresponds to a (massive, increasing with decreasing grid resolution) potential error in the density estimation.

For instance, consider the case where you have an adaptive simulation with a structure like a galaxy. Let's say the galaxy has a typical diameter of 50 kpc, disk thickness 1kpc, and we'd like to make an image with a typical 512^2 resolution (i.e. each pixel represents a 100 pc^2 region). I'll also use a typical ISM density of 1e3 nh/cm^3.

If I simulate this galaxy with a particle mass resolution of 10^5 Msun, I will have ~2000 particles in this pixel. A pretty big number of them won't overlap with the center of the pixel. If I do the same, with a mass resolution of 10^3 Msun, I now have ~200000 particles in this pixel! Now my measured density with your method will be tiny.

That's exactly what I've done here in this image (albeit with a far lower resolution galaxy). The left panel removes the explicit case we have in our 'fast' backend (we treat cases with zero overlap as PIC, i.e. we add their entire mass onto the cell their center lives within). Without the explicit correction, notice how much less dense the center is!

The worst thing about this error is that it is introduced into the area where you actually have the most information about the simulation...

error_image.pdf

[JBorrow]

Another way to think about this is taking a uniform density of particles.

Your method as it stands will, as the resolution of the grid decreases relative to the background particles, effectively act as a random sampling filter. What this corresponds to is the image below:

test_images.pdf

When it should look like this:

test_images_swio.pdf

As histograms:

resolution_test_uniform.pdf

v.s. expected result (notice as resolution decreases, and as we average over more and more particles, the distribution of pixel densities gets thinner, towards the mean density in the box (vertical blue line).

resolution_test_uniform_swio.pdf

[ttricco]

Thanks for being so thorough on this point, Josh. You've definitely convinced me that we need to move this min smoothing length up in our priority list.

What I'm going to do is close this current issue, as we have a PR merging in to address the original issue highlighted with exact interpolation, and then I'll open a new issue focused on this issue with our fast interpolation.