(refac): Relax periodic endpoint check + PeriodicBC(check=false) opt-out

docs/src/boundary-conditions/overview.md

@@ -109,7 +109,7 @@ PeriodicBC()     # S(x) = S(x + τ) with C² continuity
 | `ZeroCurvBC()` | S''=0 at both ends | Zero-curvature assumption |
 | `ZeroSlopeBC()` | S'=0 at both ends | Flat endpoints |
 | `BCPair(...)` | Custom at each end | Known derivatives |
-| `PeriodicBC()` | True periodicity (inclusive) | Cyclic data with `y[1] == y[end]` |
+| `PeriodicBC()` | True periodicity (inclusive) | Cyclic data with `y[1] ≈ y[end]` |
 | `PeriodicBC(endpoint=:exclusive)` | True periodicity (exclusive) | FFT grids, `[0, 2π)` data |
 | `CubicFit()` | 4-point polynomial fit | Exact for cubic polynomials |
 

docs/src/boundary-conditions/periodicbc.md

@@ -12,7 +12,7 @@ These two features sound similar but solve fundamentally different problems:
 |---|---|---|
 | **What it does** | Solves a cyclic tridiagonal system (Sherman-Morrison) so the spline is **C² continuous** at the period boundary | Maps out-of-domain queries back into `[x₁, xₙ]` via modular arithmetic |
 | **Smoothness** | ``S, S', S''`` all match at the wrap point | No smoothness guarantee — may have jumps in value, slope, or curvature |
-| **Data requirement** | `y[1] == y[end]` (inclusive) or `endpoint=:exclusive` | None |
+| **Data requirement** | `y[1] ≈ y[end]` (inclusive) or `endpoint=:exclusive` | None |
 | **Works with** | Cubic splines only | Any interpolation method |
 | **Use case** | Physically periodic signals (angles, phases, Fourier-sampled data) | Quick "repeat" behavior without physical periodicity |
 
@@ -27,16 +27,20 @@ These two features sound similar but solve fundamentally different problems:
 
 ### Inclusive Mode (Default)
 
-The grid includes the repeated start point at the end: `y[1] == y[end]` (exact equality required). This is the standard definition in most spline libraries.
+The grid includes the repeated start point at the end: `y[1] ≈ y[end]` (validated via `isapprox` with `atol = 8eps(T)` noise floor). This is the standard definition in most spline libraries.
 
 ```julia
 # Grid covers [0, 2π], with repeated endpoint
 x = range(0, 2π, 65)  # 65 points, last point is 2π
-y = cos.(x)           # y[1] == y[end] (cos(0) == cos(2π) == 1.0)
+y = sin.(x)           # sin(0) ≈ sin(2π) — passes isapprox check
 
 itp = cubic_interp(x, y; bc=PeriodicBC())
 ```
 
+!!! tip "Scaled data"
+    For scaled data like `1e6 .* sin.(x)` where noise exceeds the `8eps` floor,
+    either set `y[end] = y[1]` explicitly or use `PeriodicBC(check=false)` to skip validation.
+
 ### Exclusive Mode
 
 The grid contains only **unique** points. The theoretical "next" point would be the start of the next period. This is common in FFT-based applications or simulation grids.

docs/src/interpolation/cubic.md

@@ -85,7 +85,7 @@ cubic_interp(x, y, 1.0; bc=ZeroCurvBC())             # zero curvature at endpoin
 cubic_interp(x, y, 1.0; bc=ZeroSlopeBC())            # flat endpoints
 cubic_interp(x, y, 1.0; bc=BCPair(Deriv1(1), Deriv2(0)))  # custom
 
-# Periodic (closed curve) - requires y[1] == y[end]
+# Periodic (closed curve) - requires y[1] ≈ y[end]
 cubic_interp(x, y, 1.0; bc=PeriodicBC())
 
 # In-place evaluation (zero allocation)

docs/src/migration/to_v0.4.md

@@ -20,24 +20,24 @@ Series([y1, y2])            # vector of vectors
 Series(hcat(y1, y2))        # matrix (columns = series)
 ```
 
-## 2. `PeriodicBC()` Strict Endpoint Check
+## 2. `PeriodicBC()` Endpoint Check
 
-The default `:inclusive` mode now requires `y[1] == y[end]` (exact equality) instead of `isapprox`. This catches silent data errors from floating-point arithmetic:
+**v0.4.0–v0.4.5**: The default `:inclusive` mode required `y[1] == y[end]` (strict bitwise equality). This was overly strict for computed data — e.g., `sin(0) != sin(2π)` due to floating-point round-off.
+
+**v0.4.6+**: Relaxed to `isapprox` with `atol = 8eps(T)` noise floor. Typical computed periodic data now passes without manual fixup:
 
 ```julia
 t = range(0, 2π, 101)
-y = sin.(t)
-# sin(2π) ≈ -2.4e-16, NOT exactly 0.0
-
-cubic_interp(t, y; bc=PeriodicBC())
-# ERROR: y[1] != y[end] for inclusive PeriodicBC
+y = sin.(t)               # sin(0) ≈ sin(2π) — passes (diff ~1 eps)
+cubic_interp(t, y; bc=PeriodicBC())  # works (would ERROR in v0.4.0–v0.4.5)
 ```
 
-Fix by explicitly ensuring the endpoint matches:
+For scaled data where noise exceeds `8eps`, use `check=false` or set the endpoint explicitly:
 
 ```julia
-y[end] = y[1]   # force exact equality
-cubic_interp(t, y; bc=PeriodicBC())  # works
+y_scaled = 1e6 .* sin.(t)
+cubic_interp(t, y_scaled; bc=PeriodicBC(check=false))  # skip validation
+# or: y_scaled[end] = y_scaled[1]
 ```
 
 Alternatively, use `:exclusive` mode if your data does not include the repeated endpoint:

src/core/bc_types.jl

@@ -173,7 +173,7 @@ BCPair(t::Tuple{L, R}) where {L <: PointBC, R <: PointBC} =
 
 
 """
-    PeriodicBC{E, P} <: AbstractBC
+    PeriodicBC{E, P, C} <: AbstractBC
 
 Periodic boundary condition: S(x_0) = S(x_n), S'(x_0) = S'(x_n), S''(x_0) = S''(x_n)
 
@@ -182,6 +182,7 @@ Internally, periodic BC uses Sherman-Morrison solver with `PeriodicData{T}` for
 # Type Parameters
 - `E::Symbol`: `:inclusive` or `:exclusive` (compile-time endpoint convention)
 - `P`: `Nothing` (inclusive or auto-infer) or `<:AbstractFloat` (explicit period)
+- `C::Bool`: whether to validate `y[1] ≈ y[end]` at construction time (default `true`)
 
 # Endpoint Conventions
 - **Inclusive** (`endpoint=:inclusive`, default): `y[1] ≈ y[end]` required (standard convention)
@@ -200,22 +201,27 @@ itp = cubic_interp(x, y; bc=PeriodicBC(endpoint=:exclusive, period=2π))
 # Exclusive with Range grid (period auto-inferred)
 x = range(0, step=2π/64, length=64)
 itp = cubic_interp(x, sin.(x); bc=PeriodicBC(endpoint=:exclusive))
+
+# Skip endpoint check (e.g., scaled data where noise > 8eps)
+itp = cubic_interp(x, 1e6 .* sin.(x); bc=PeriodicBC(check=false))
 ```
 """
-struct PeriodicBC{E, P} <: AbstractBC
+struct PeriodicBC{E, P, C} <: AbstractBC
     period::P         # Nothing or AbstractFloat
-    function PeriodicBC{E, P}(period::P) where {E, P}
+    function PeriodicBC{E, P, C}(period::P) where {E, P, C}
         E isa Symbol || error("PeriodicBC type parameter E must be a Symbol")
         E in (:inclusive, :exclusive) || error("PeriodicBC type parameter E must be :inclusive or :exclusive")
-        return new{E, P}(period)
+        C isa Bool || error("PeriodicBC type parameter C must be a Bool")
+        return new{E, P, C}(period)
     end
 end
 
-# Accessor for endpoint (from type parameter, zero-cost)
+# Accessors (from type parameters, zero-cost)
 @inline endpoint(::PeriodicBC{E}) where {E} = E
+@inline periodic_check(::PeriodicBC{E, P, C}) where {E, P, C} = C
 
 # Keyword constructor with validation (also serves as zero-arg constructor via defaults)
-function PeriodicBC(; endpoint::Symbol = :inclusive, period::Union{Real, Nothing} = nothing)
+function PeriodicBC(; endpoint::Symbol = :inclusive, period::Union{Real, Nothing} = nothing, check::Bool = true)
     endpoint in (:inclusive, :exclusive) || throw(
         ArgumentError(
             "endpoint must be :inclusive or :exclusive, got :$endpoint"
@@ -227,14 +233,14 @@ function PeriodicBC(; endpoint::Symbol = :inclusive, period::Union{Real, Nothing
                 "period is not applicable for endpoint=:inclusive (y[1]≈y[end] convention)"
             )
         )
-        return PeriodicBC{:inclusive, Nothing}(nothing)
+        return PeriodicBC{:inclusive, Nothing, check}(nothing)
     else # :exclusive
         if period !== nothing
             p = float(period)
             p > 0 || throw(ArgumentError("period must be positive, got $period"))
-            return PeriodicBC{:exclusive, typeof(p)}(p)
+            return PeriodicBC{:exclusive, typeof(p), check}(p)
         else
-            return PeriodicBC{:exclusive, Nothing}(nothing)  # infer from Range at build time
+            return PeriodicBC{:exclusive, Nothing, check}(nothing)  # infer from Range at build time
         end
     end
 end

src/core/periodic.jl

@@ -51,18 +51,49 @@ end
 """
     _check_periodic_endpoints(y::AbstractVector)
 
-Validate that `y[1] == y[end]` for periodic boundary conditions (inclusive endpoint).
+Validate that `y[1] ≈ y[end]` for periodic boundary conditions (inclusive endpoint).
 Called once at construction time (zero runtime overhead).
 
-Uses strict `==` equality — no approximate comparison. This is universal for all
-value types (scalars, vectors, duck-typed custom types) without requiring `norm`,
-`isapprox`, or any tolerance parameters.
+Three-tier dispatch based on element type:
 
-If your data is computed (e.g., `sin.(range(0, 2π, n))`), set `y[end] = y[1]`
-explicitly to ensure exact periodicity.
+- **`AbstractFloat`**: `isapprox` with `atol = 8eps(T)` — handles both relative differences
+  (e.g., `cos(0) ≈ cos(2π)`) and near-zero noise floor (e.g., `sin(0)` vs `sin(2π)`).
+  The `8eps` constant is compile-time folded (zero overhead vs plain `isapprox`).
+- **`Complex{<:AbstractFloat}`**: same, using `eps(real(T))`.
+- **Other `_PromotableValue`** (Integer, Rational): `isapprox` with default tolerances.
+- **Duck types** (Dual, SVector, ...): strict `==` (isapprox semantics not guaranteed).
+
+!!! note "Scaled near-zero endpoints"
+    `atol = 8eps` covers direct evaluations (e.g., `sin.(x)`), but not scaled
+    variants (e.g., `1e6 .* sin.(x)` where noise ≈ 1e6·eps). For those cases,
+    set `y[end] = y[1]` explicitly, or use `PeriodicBC(check=false)` to skip
+    this validation.
 
 Throws `ArgumentError` if endpoints differ.
 """
+@inline function _check_periodic_endpoints(bc::PeriodicBC, y::AbstractVector)
+    periodic_check(bc) || return nothing
+    _check_periodic_endpoints(y)
+    return nothing
+end
+
+@inline function _check_periodic_endpoints(y::AbstractVector{T}) where {T <: AbstractFloat}
+    isapprox(first(y), last(y); atol = 8 * eps(T)) || _throw_periodic_endpoint_error(first(y), last(y))
+    return nothing
+end
+
+@inline function _check_periodic_endpoints(y::AbstractVector{Complex{T}}) where {T <: AbstractFloat}
+    isapprox(first(y), last(y); atol = 8 * eps(T)) || _throw_periodic_endpoint_error(first(y), last(y))
+    return nothing
+end
+
+# Integer, Rational: isapprox with default tolerances (effectively ==)
+@inline function _check_periodic_endpoints(y::AbstractVector{<:_PromotableValue})
+    isapprox(first(y), last(y)) || _throw_periodic_endpoint_error(first(y), last(y))
+    return nothing
+end
+
+# Duck-type fallback: strict == (isapprox not guaranteed for arbitrary types)
 @inline function _check_periodic_endpoints(y::AbstractVector)
     _extract_primal(first(y)) == _extract_primal(last(y)) || _throw_periodic_endpoint_error(first(y), last(y))
     return nothing
@@ -71,10 +102,11 @@ end
 @noinline function _throw_periodic_endpoint_error(y1, yn)
     throw(
         ArgumentError(
-            "PeriodicBC (inclusive endpoint) requires y[1] == y[end], " *
+            "PeriodicBC (inclusive endpoint) requires y[1] ≈ y[end], " *
                 "got y[1]=$y1, y[end]=$yn. " *
-                "Tip: set y[end] = y[1] to ensure exact periodicity, or use " *
-                "PeriodicBC(endpoint=:exclusive) if your data does not repeat the first point."
+                "Tip: set y[end] = y[1] explicitly, use " *
+                "PeriodicBC(endpoint=:exclusive) if your data does not repeat the first point, " *
+                "or PeriodicBC(check=false) to skip this validation."
         )
     )
 end
@@ -93,9 +125,10 @@ end
 @noinline function _throw_periodic_nd_error(d, v_first, v_last)
     throw(
         ArgumentError(
-            "Periodic BC on dim $d requires data[1,...] == data[end,...], " *
+            "Periodic BC on dim $d requires data[1,...] ≈ data[end,...], " *
                 "but found data[1,...]=$v_first, data[end,...]=$v_last. " *
-                "Tip: set the last slice equal to the first along dim $d."
+                "Tip: set the last slice equal to the first along dim $d, " *
+                "or use PeriodicBC(check=false) to skip this validation."
         )
     )
 end
@@ -170,8 +203,8 @@ Used so that `itp.bc` always carries the actual period for display/introspection
 Uses the inner constructor directly to bypass keyword-constructor validation
 (which rejects `period` for inclusive BCs).
 """
-@inline _with_resolved_period(::PeriodicBC{E}, period::T) where {E, T} =
-    PeriodicBC{E, T}(period)
+@inline _with_resolved_period(::PeriodicBC{E, <:Any, C}, period::T) where {E, T, C} =
+    PeriodicBC{E, T, C}(period)
 
 """
     _extend_exclusive(x, y, bc::PeriodicBC) -> (x_ext, y_ext)

src/cubic/cubic_interpolant.jl

@@ -290,7 +290,7 @@ so the pool memory can be safely reused after this function returns.
         search::AbstractSearchPolicy = AutoSearch()
     ) where {Tg <: AbstractFloat, Tv}
     x, y = _prepare_periodic(x, y, bc)
-    _check_periodic_endpoints(y)
+    _check_periodic_endpoints(bc, y)
     cache = _get_cubic_cache(x, PeriodicBC(), autocache)
     tmp_z = similar!(pool, y)
     _solve_system!(tmp_z, cache, y, cache.bc_config)

src/cubic/cubic_oneshot.jl

@@ -166,7 +166,7 @@ manages their lifetime. Follows the `_create_spacing_pooled(pool, ...)` pattern.
     end
 
     # ── Solve periodic tridiagonal system ──
-    _check_periodic_endpoints(y_p)
+    _check_periodic_endpoints(bc, y_p)
     cache = _get_cubic_cache(x_p, PeriodicBC(), autocache)
     z = acquire!(pool, Tv, length(y_p))
     _solve_system!(z, cache, y_p, cache.bc_config)

src/cubic/nd/cubic_nd_build.jl

@@ -251,11 +251,32 @@ the index expressions into the method body at specialization time.
     first_idx = [d == D ? 1 : idx_vars[d] for d in 1:N]
     last_idx = [d == D ? :n_D : idx_vars[d] for d in 1:N]
 
-    # Inner comparison body: strict == equality, no tolerance parameters
-    check = quote
-        v1 = @inbounds data[$(first_idx...)]
-        vn = @inbounds data[$(last_idx...)]
-        v1 == vn || _throw_periodic_nd_error($D, v1, vn)
+    # Inner comparison body: isapprox for AbstractFloat/Complex, strict == for others
+    if Tv <: AbstractFloat
+        check = quote
+            v1 = @inbounds data[$(first_idx...)]
+            vn = @inbounds data[$(last_idx...)]
+            isapprox(v1, vn; atol = 8 * eps($Tv)) || _throw_periodic_nd_error($D, v1, vn)
+        end
+    elseif Tv <: Complex{<:AbstractFloat}
+        RT = real(Tv)
+        check = quote
+            v1 = @inbounds data[$(first_idx...)]
+            vn = @inbounds data[$(last_idx...)]
+            isapprox(v1, vn; atol = 8 * eps($RT)) || _throw_periodic_nd_error($D, v1, vn)
+        end
+    elseif Tv <: _PromotableValue
+        check = quote
+            v1 = @inbounds data[$(first_idx...)]
+            vn = @inbounds data[$(last_idx...)]
+            isapprox(v1, vn) || _throw_periodic_nd_error($D, v1, vn)
+        end
+    else
+        check = quote
+            v1 = @inbounds data[$(first_idx...)]
+            vn = @inbounds data[$(last_idx...)]
+            v1 == vn || _throw_periodic_nd_error($D, v1, vn)
+        end
     end
 
     # Wrap in nested loops over all dims except D (outermost = N, innermost = 1)
@@ -374,7 +395,7 @@ end
     # in the data (it is added by _prepare_periodic_nd/_prepare_periodic_nd_pooled after
     # this validation).  Checking data[1] ≈ data[end] on unextended exclusive data would
     # produce false positives for perfectly valid periodic inputs.
-    if bcs[D] isa PeriodicBC{:inclusive}
+    if bcs[D] isa PeriodicBC{:inclusive} && periodic_check(bcs[D])
         _check_periodic_data_noalloc!(data, Val(D), Tg)
     end
     polyfit_deg = get_polyfit_degree(bcs[D])

src/tensor_product/tensor_product_interpolant.jl

@@ -13,7 +13,7 @@
 # the existing constructors already accept Union{Single, NTuple} for these kwargs.
 
 function _interp_nd_dispatch(
-        grids, data, methods::Tuple{<:CubicInterp, Vararg{<:CubicInterp}}, coeffs, extrap, search
+        grids, data, methods::Tuple{CubicInterp, Vararg{CubicInterp}}, coeffs, extrap, search
     )
     bcs = map(m -> m.bc, methods)
     return cubic_interp(grids, data; bc = bcs, extrap = extrap, search = search, coeffs = coeffs)
@@ -26,14 +26,14 @@ function _interp_nd_dispatch(
 end
 
 function _interp_nd_dispatch(
-        grids, data, methods::Tuple{<:QuadraticInterp, Vararg{<:QuadraticInterp}}, ::Any, extrap, search
+        grids, data, methods::Tuple{QuadraticInterp, Vararg{QuadraticInterp}}, ::Any, extrap, search
     )
     bcs = map(m -> m.bc, methods)
     return quadratic_interp(grids, data; bc = bcs, extrap = extrap, search = search)
 end
 
 function _interp_nd_dispatch(
-        grids, data, methods::Tuple{<:ConstantInterp, Vararg{<:ConstantInterp}}, ::Any, extrap, search
+        grids, data, methods::Tuple{ConstantInterp, Vararg{ConstantInterp}}, ::Any, extrap, search
     )
     sides = map(m -> m.side, methods)
     return constant_interp(grids, data; side = sides, extrap = extrap, search = search)

src/tensor_product/tensor_product_oneshot.jl

@@ -117,7 +117,7 @@ end
 
 function _interp_nd_oneshot_dispatch(
         grids, data, query,
-        methods::Tuple{<:CubicInterp, Vararg{<:CubicInterp}},
+        methods::Tuple{CubicInterp, Vararg{CubicInterp}},
         deriv, extrap, search, hints,
     )
     bcs = map(m -> m.bc, methods)
@@ -134,7 +134,7 @@ end
 
 function _interp_nd_oneshot_dispatch(
         grids, data, query,
-        methods::Tuple{<:QuadraticInterp, Vararg{<:QuadraticInterp}},
+        methods::Tuple{QuadraticInterp, Vararg{QuadraticInterp}},
         deriv, extrap, search, hints,
     )
     bcs = map(m -> m.bc, methods)
@@ -143,7 +143,7 @@ end
 
 function _interp_nd_oneshot_dispatch(
         grids, data, query,
-        methods::Tuple{<:ConstantInterp, Vararg{<:ConstantInterp}},
+        methods::Tuple{ConstantInterp, Vararg{ConstantInterp}},
         deriv, extrap, search, hints,
     )
     sides = map(m -> m.side, methods)
@@ -177,7 +177,7 @@ end
 
 function _interp_nd_oneshot_batch_dispatch!(
         output, grids, data, queries,
-        methods::Tuple{<:CubicInterp, Vararg{<:CubicInterp}},
+        methods::Tuple{CubicInterp, Vararg{CubicInterp}},
         deriv, extrap, search, hints,
     )
     bcs = map(m -> m.bc, methods)
@@ -194,7 +194,7 @@ end
 
 function _interp_nd_oneshot_batch_dispatch!(
         output, grids, data, queries,
-        methods::Tuple{<:QuadraticInterp, Vararg{<:QuadraticInterp}},
+        methods::Tuple{QuadraticInterp, Vararg{QuadraticInterp}},
         deriv, extrap, search, hints,
     )
     bcs = map(m -> m.bc, methods)
@@ -203,7 +203,7 @@ end
 
 function _interp_nd_oneshot_batch_dispatch!(
         output, grids, data, queries,
-        methods::Tuple{<:ConstantInterp, Vararg{<:ConstantInterp}},
+        methods::Tuple{ConstantInterp, Vararg{ConstantInterp}},
         deriv, extrap, search, hints,
     )
     sides = map(m -> m.side, methods)