Feature Request: Support for Exclusive Periodic Boundaries in ND Interpolators

[henry2004y]

Hi again!

FastInterpolations.jl is quickly evolving. I'm going to seek help on my common use cases of interpolating a field from numerical simulation outputs. As the boundary conditions are only supported for part of the methods (e.g. not for linear_interp), I am using WrapExtrap() currently in my workflow.

The grid setup in the following demo is a common pattern in Finite Volume simulations and periodic signal processing, where the grid points represent unique cell centers (or unique samples) and the total grid span is one short of the full period. In Interpolations.jl, this is distinguished by the gridstyle (OnGrid(), OnCell()) Interpolations.jl

Grid Configuration

• Coordinates: x = [0.5, 1.5, 2.5] ($N=3$ points).
• Values: y = [1.0, 2.0, 3.0].
• Spacing: $\Delta x = 1.0$.
• Desired Period: $L = N \times \Delta x = 3.0$.

This setup describes a domain of length 3.0 (e.g., $[0.0, 3.0]$) divided into 3 cells of width $1.0$. The data points are located at the cell centers. In a periodic system of this type:

• The point $x = 3.0$ is equivalent to $x = 0.0$.
• The point $x = 3.5$ is equivalent to $x = 0.5$ (the first cell center).
• Expected Value at $x=3.5$: 1.0.

The Issue with Standard Wrapping

By default, FastInterpolations.jl's Extrap(:wrap) uses the distance between the first and last points as the period ($L_{\text{default}} = x_{\text{end}} - x_{\text{start}} = 2.0$).

• Standard Result: $x=3.5 \pmod{2.0}$ maps to $x=1.5$.
• Returned Value: 2.0 (Incorrect for this grid type).

Expected native support

While cubic_interp supports this via PeriodicBC(endpoint=:exclusive, period=3.0), the generic and ND methods (linear_interp, cardinal_interp, etc.) do not currently provide a unified keyword-based way to specify the Exclusive Period $L$. Is it possible to have a similar support with Extrap(:warp) without the need of manual padding? Or if you have a better idea, let me know!

using FastInterpolations

println("--- Exclusive Periodic Grid Demo ---")

# 1. Setup an "Exclusive" Grid (e.g., cell centers)
# Grid: [0.5, 1.5, 2.5], Values: [1.0, 2.0, 3.0]
# Expected Period L = 3.0 (not 2.0! Next point at 3.5 should wrap to 0.5)
x = [0.5, 1.5, 2.5]
y = [1.0, 2.0, 3.0]
dx = 1.0
L_expected = 3.0

println("Original Grid: ", x)
println("Original Values: ", y)
println("Expected Period (N*dx): ", L_expected)

# 2. Standard FastInterpolations wrapping (Inclusive)
# This uses period = last(x) - first(x) = 2.0
itp_standard = linear_interp(x, y; extrap = Extrap(:wrap))
val_standard = itp_standard(3.5)

println("\nStandard Wrap (Inclusive):")
println("  Query at x=3.5 maps to: ", 0.5 + mod(3.5 - 0.5, 2.5 - 0.5))
println("  Value at x=3.5: ", val_standard)

# 3. Exclusive Wrap via Padding
# We add a virtual point at x=3.5 with value y=1.0
x_padded = [0.5, 1.5, 2.5, 3.5]
y_padded = [1.0, 2.0, 3.0, 1.0]

itp_exclusive = linear_interp(x_padded, y_padded; extrap = Extrap(:wrap))
val_exclusive = itp_exclusive(3.5)

println("\nExclusive Wrap (via Padding):")
println("  Query at x=3.5 maps to: ", 0.5 + mod(3.5 - 0.5, 3.5 - 0.5))
println("  Value at x=3.5: ", val_exclusive, " (Success!)")

# 4. Same demo for Cardinal (Order 3) Interpolation
# This is even more important for splines to ensure C2 continuity across the boundary
itp_cardinal = cardinal_interp(x_padded, y_padded; extrap = Extrap(:wrap))
val_cardinal = itp_cardinal(3.5)

println("\nCardinal Exclusive (Order 3):")
println("  Value at x=3.5: ", val_cardinal)

# 5. Native FastInterpolations PeriodicBC (Order 3)
# Note: For non-uniform grids (including Vector), we need to specify the period.
itp_periodicBC = cubic_interp(x, y; bc = PeriodicBC(endpoint = :exclusive, period = 3.0))
val_periodicBC = itp_periodicBC(3.5)

println("\nNative PeriodicBC (endpoint=:exclusive, period=3.0):")
println("  Value at x=3.5: ", val_periodicBC)

Original Grid: [0.5, 1.5, 2.5]
Original Values: [1.0, 2.0, 3.0]
Expected Period (N*dx): 3.0

Standard Wrap (Inclusive):
  Query at x=3.5 maps to: 1.5
  Value at x=3.5: 2.0

Exclusive Wrap (via Padding):
  Query at x=3.5 maps to: 0.5
  Value at x=3.5: 1.0 (Success!)

Cardinal Exclusive (Order 3):
  Value at x=3.5: 1.0

Native PeriodicBC (endpoint=:exclusive, period=3.0):
  Value at x=3.5: 1.0

[mgyoo86]

Hi again!
Thanks for the detailed write-up and the concrete example — really appreciate it!

I'm still thinking about what the best API design would be— want to get it right without being too invasive.
The direction I'm currently leaning toward is extending PeriodicBC support to more methods (linear, constant, PCHIP/Cardinal/Akima) rather than adding parameters to WrapExtrap.

The reason is that the role of WrapExtrap is purely simple and geometric — it maps queries back into the domain via mod, with no assumptions about the data (values don't even need to be continuous).
On the other hand, PeriodicBC imposes an actual constraint to the system, making the interpolation continuous across the periodic boundary. Currently only cubic_interp supports it, by solving a different system (cyclic tridiagonal) to enforce C2-continuity. Other methods haven't needed this kind of enforcement until now.
(See 📖 PeriodicBC docs for more details)

So the plan is to extend PeriodicBC to each method where it makes sense (linear, constant, and the Hermite family). This would be especially meaningful for the cubic Hermite family (PCHIP, Cardinal, Akima). Right now, the cubic Hermite family's boundary slopes are computed with reduced stencils (one-sided FD, extrapolated secants), so they're not truly C1-continuous at the boundary. With PeriodicBC, the slope computation would use periodic stencils that wrap around, giving genuine C1-continuity.

I'm currently working on another feature release first, and will revisit this afterwards.
Will update here as progress is made :)

[mgyoo86]

Hello!
This is just a general progress update, not a release yet.

The master branch now supports the exclusive PeriodicBC for constant_interp and linear_interp (some advanced APIs are not covered yet).

If you have time, please try it on master and let me know if this works for you.

Here is an MWE:

using FastInterpolations

x = [0.5, 1.5, 2.5]
y = [1.0, 2.0, 3.0]
dx = 1.0
L_expected = 3.0

# 1. Interpolant path
# This internally extends x & y to x=[0.5, 1.5, 2.5, 3.5] and y=[1.0, 2.0, 3.0, 1.0])
# Create the interpolation object with exclusive PeriodicBC
itp = linear_interp(x, y; bc = PeriodicBC(endpoint=:exclusive, period=L_expected))
itp(3.5) # Should wrap to 0.5 and return 1.0

# 2. One-shot path
# Zero-memory copy or overhead
# Evaluation at 3.5 with the same exclusive PeriodicBC
linear_interp(x, y, 3.5; bc = PeriodicBC(endpoint=:exclusive, period=L_expected)) # Should also return 1.0

Note that the interpolant path simply extends the grid and data, since it allocates new memory and copies from the original data anyway.

On the other hand, the one-shot path does not allocate new memory — it uses the original grid and data in a periodic manner (which required major refactoring).

I'll extend the exclusive PeriodicBC to other methods, along with some performance improvements.
This is taking some time, but it's worth implementing.

I'll update you again when everything is ready!

[henry2004y]

Thanks for the update! A rework like this would take months without AI, but gladly we are in a new era now.

[mgyoo86]

A rework like this would take months without AI, but gladly we are in a new era now.

This is so true! I can’t imagine myself manually reading/typing hundreds of lines just for a refactor 😂

Still designing clean, extensible API layers, while addressing performance issues, is not easy, but AI is definetly very helpful!

[mgyoo86]

@henry2004y
A new version v0.4.10 is finally released,
and it now supports the PeriodicBC for all methods except for the quadratic one.
📖 PeriodicBC docs

It took much longer than expected due to massive refactorings 😅
Please give it a try, and let me know if you find any issues.
Thank you!