Apply a pre-commit config/format all/document LuaMap/cubic interpolation

[davschneller]

Also format all of the code.
And add the LuaMap to the docs.
Add cubic interpolation.

[Thomas-Ulrich]

Ok, maybe there is still room for improvement ^^
left cubic, right linear:

[davschneller]

Hmm, looks like it. :) Smoothes a lot seemingly. For clarity — is the linear version computed with the same easi version? (i.e. this one) Just to exclude that there's something more fundamentally wrong.

I'll look into the code. ...

[davschneller]

One more question; is the picture supposed to show something close to the initial setup, or after the simulation is done? If it's the latter, then how does the initial projection look like by chance? (otherwise, the left part somehow looks like the bottom left part of the right image, magnified — maybe that'll help)

[Thomas-Ulrich]

I've verified the cubic interpolation code using some vibe-coded scripts (Gemini 3.6 flash):

generate_grids.py

import os
import netCDF4 as nc
import numpy as np

# Ensure target folder exists
output_dir = "ASAGI_files"
os.makedirs(output_dir, exist_ok=True)


def build_1d_linear():
    fn = os.path.join(output_dir, "test_1d_grid.nc")
    with nc.Dataset(fn, "w", format="NETCDF3_64BIT_OFFSET") as ds:
        nx = 11
        ds.createDimension("x", nx)
        x_var = ds.createVariable("x", "f4", ("x",))
        data_var = ds.createVariable("data", "f4", ("x",))

        x_vals = np.linspace(0.0, 10.0, nx)
        x_var[:] = x_vals
        data_var[:] = 2.0 * x_vals
    print(f"Created {fn}")


def build_1d_cubic():
    fn = os.path.join(output_dir, "test_cubic_grid.nc")
    with nc.Dataset(fn, "w", format="NETCDF3_64BIT_OFFSET") as ds:
        nx = 11
        ds.createDimension("x", nx)
        x_var = ds.createVariable("x", "f4", ("x",))
        data_var = ds.createVariable("data", "f4", ("x",))

        x_vals = np.linspace(0.0, 10.0, nx)
        x_var[:] = x_vals
        data_var[:] = x_vals**3
    print(f"Created {fn}")


def build_2d_cubic():
    fn = os.path.join(output_dir, "test_2d_cubic_grid.nc")
    with nc.Dataset(fn, "w", format="NETCDF3_64BIT_OFFSET") as ds:
        nx, ny = 11, 11
        ds.createDimension("x", nx)
        ds.createDimension("y", ny)
        x_var = ds.createVariable("x", "f4", ("x",))
        y_var = ds.createVariable("y", "f4", ("y",))
        data_var = ds.createVariable("data", "f4", ("y", "x"))

        x_vals = np.linspace(0.0, 10.0, nx)
        y_vals = np.linspace(0.0, 10.0, ny)
        x_var[:] = x_vals
        y_var[:] = y_vals

        X, Y = np.meshgrid(x_vals, y_vals)
        data_var[:, :] = X**3 + Y**3
    print(f"Created {fn}")


if __name__ == "__main__":
    build_1d_linear()
    build_1d_cubic()
    build_2d_cubic()

test_asagi_all.yaml

!Switch
[val_1d_lin]: !AffineMap
  matrix:
    x: [1.0, 0.0, 0.0]
  translation:
    x: 0.0
  components: !Any
    - !ASAGI
        file: ./ASAGI_files/test_1d_grid.nc
        parameters: [val_1d_lin]
        var: data
        interpolation: linear
    - !ConstantMap
        map:
          val_1d_lin: -999.0

[val_1d_cub]: !AffineMap
  matrix:
    x: [1.0, 0.0, 0.0]
  translation:
    x: 0.0
  components: !Any
    - !ASAGI
        file: ./ASAGI_files/test_cubic_grid.nc
        parameters: [val_1d_cub]
        var: data
        interpolation: cubic
    - !ConstantMap
        map:
          val_1d_cub: -999.0

[val_2d_cub]: !AffineMap
  matrix:
    x: [1.0, 0.0, 0.0]
    y: [0.0, 1.0, 0.0]
  translation:
    x: 0.0
    y: 0.0
  components: !Any
    - !ASAGI
        file: ./ASAGI_files/test_2d_cubic_grid.nc
        parameters: [val_2d_cub]
        var: data
        interpolation: cubic
    - !ConstantMap
        map:
          val_2d_cub: -999.0

and

verify_all.py

import easi
import numpy as np

# Test positions [x, y, z]
points = np.array(
    [
        [2.5, 0.0, 0.0],  # 1D Linear (interior) & 1D Cubic (interior)
        [0.5, 0.0, 0.0],  # 1D Cubic boundary (Linear fallback)
        [2.5, 3.5, 0.0],  # 2D Bicubic interior
        [0.5, 9.5, 0.0],  # 2D Boundary (Bilinear fallback)
    ],
    dtype=np.float64,
)
tags = np.array([1, 1, 1, 1])

yaml_file = "test_asagi_all.yaml"

# 1. Test 1D Linear (x = 2.5)
res_lin = easi.evaluate_model(points[:1], tags[:1], ["val_1d_lin"], yaml_file)
val_lin = res_lin["val_1d_lin"][0]
exp_lin = 2.0 * 2.5
assert np.isclose(val_lin, exp_lin), f"1D Linear failed: {val_lin} != {exp_lin}"
print(f"[PASS] 1D Linear (x=2.5): Evaluated={val_lin:.4f} | Expected={exp_lin:.4f}")

# 2. Test 1D Cubic Interior (x = 2.5) & Boundary Fallback (x = 0.5)
res_cub = easi.evaluate_model(points[:2], tags[:2], ["val_1d_cub"], yaml_file)

# 2a. Interior
val_cub_int = res_cub["val_1d_cub"][0]
exp_cub_int = 2.5**3  # 15.625
assert np.isclose(
    val_cub_int, exp_cub_int, atol=1e-6
), f"1D Cubic interior failed: {val_cub_int} != {exp_cub_int}"
print(
    f"[PASS] 1D Cubic Interior (x=2.5): Evaluated={val_cub_int:.4f} | Expected={exp_cub_int:.4f}"
)

# 2b. Boundary Fallback
val_cub_bnd = res_cub["val_1d_cub"][1]
exp_cub_bnd = 0.5  # Linear fallback between (0,0) and (1,1)
assert np.isclose(
    val_cub_bnd, exp_cub_bnd, atol=1e-6
), f"1D Cubic boundary failed: {val_cub_bnd} != {exp_cub_bnd}"
print(
    f"[PASS] 1D Cubic Boundary (x=0.5): Evaluated={val_cub_bnd:.4f} | Expected={exp_cub_bnd:.4f}"
)

# 3. Test 2D Cubic (Interior & Bilinear boundary fallback)
res_2d = easi.evaluate_model(points[2:], tags[2:], ["val_2d_cub"], yaml_file)

# 3a. 2D Interior
val_2d_int = res_2d["val_2d_cub"][0]
exp_2d_int = 2.5**3 + 3.5**3  # 58.5
assert np.isclose(val_2d_int, exp_2d_int, atol=1e-6)
print(
    f"[PASS] 2D Bicubic Interior (2.5, 3.5): Evaluated={val_2d_int:.4f} | Expected={exp_2d_int:.4f}"
)

# 3b. 2D Boundary Fallback
val_2d_bnd = res_2d["val_2d_cub"][1]
# 2D Boundary Fallback (z = x^3 + y^3 at x=0.5, y=9.5 -> 865.0)
# - (0.5, 9.5) lies in outer boundary cells for both axes (x in [0,1], y in [9,10]).
# - ASAGI falls back to Bilinear Interpolation across the 4 surrounding nodes:
#   (0,9)=729, (1,9)=730, (0,10)=1000, (1,10)=1001.
# - Step 1: Linear interp along x=0.5 at y=9 -> 729.5; at y=10 -> 1000.5.
# - Step 2: Linear interp along y=9.5 -> (729.5 + 1000.5) / 2 = 865.0.
exp_2d_bnd = 865.0
assert np.isclose(val_2d_bnd, exp_2d_bnd, atol=1e-6)
print(
    f"[PASS] 2D Bilinear Boundary (0.5, 9.5): Evaluated={val_2d_bnd:.4f} | Expected={exp_2d_bnd:.4f}"
)

print("\nAll 1D and 2D linear/cubic ASAGI tests passed successfully!")

and everything works fine, also near the boundaries.

(venv-stack24.5-3.10) di73yeq4@login01:/hppfs/work/pn49ha/di73yeq4/test_asagi> python verify_all.py 
[PASS] 1D Linear (x=2.5): Evaluated=5.0000 | Expected=5.0000
[PASS] 1D Cubic Interior (x=2.5): Evaluated=15.6250 | Expected=15.6250
[PASS] 1D Cubic Boundary (x=0.5): Evaluated=0.5000 | Expected=0.5000
[PASS] 2D Bicubic Interior (2.5, 3.5): Evaluated=58.5000 | Expected=58.5000
[PASS] 2D Bilinear Boundary (0.5, 9.5): Evaluated=865.0000 | Expected=865.0000

We could add these tests to the CI, but this would require compiling a docker image with easi in the CI.

[davschneller]

Great to see it work.

We could add these tests to the CI, but this would require compiling a docker image with easi in the CI.

We already have tests with easi+ASAGI in here — so adding these shouldn't be too hard as well. In fact, the AI (Claude) has already added a cubic test to the current easi CI setup with this PR. Here's the check file: https://github.com/SeisSol/easi/blob/a8d58b3ed6efb0611f9ab70d7a3c7b19279535ca/tests/101_asagi_nearest.cpp .

[Thomas-Ulrich]

Yes, I was looking through the changes and also discovered these tests.
Well, we have everything triple-checked then :) (and also via the Python bindings).

[Thomas-Ulrich]

Ha yes, I thought it was unusual to apply the formatting to the "external" folder.
(but fine for me if intentional).

[Thomas-Ulrich]

LGTM
(easi 1.7.0 or 1.6.3?)

[Thomas-Ulrich]

by the way, you could mention x["y"] can also be written x.y?

[davschneller]

Yes, I'll add that one.

[davschneller]

LGTM (easi 1.7.0 or 1.6.3?)

We might even want to go for 1.7.0; the new interpolation type should be enough as new feature IMO.