Add VectorInterface package extension

Implements the VectorInterface.jl API (scalartype, zerovector, scale,
add families and inner) for Fun, so Fun can be used directly as a
vector type with packages like KrylovKit without type piracy. In-place
variants operate on the coefficient vectors of Funs in compatible
spaces; the !! variants fall back to out-of-place versions on space
mismatch or scalar type promotion.

Assisted-by: Claude Fable 5 <noreply@anthropic.com>

Author: simeonschaub

Project.toml

@@ -20,9 +20,11 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [weakdeps]
 DualNumbers = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74"
+VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8"
 
 [extensions]
 ApproxFunDualNumbersExt = "DualNumbers"
+ApproxFunVectorInterfaceExt = "VectorInterface"
 
 [compat]
 AbstractFFTs = "1.0"
@@ -48,6 +50,7 @@ Reexport = "1.0"
 SpecialFunctions = "1.1, 2"
 StaticArrays = "1"
 Test = "1"
+VectorInterface = "0.4, 0.5, 0.6"
 julia = "1.10"
 
 [extras]
@@ -61,6 +64,7 @@ IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
 ParallelTestRunner = "d3525ed8-44d0-4b2c-a655-542cee43accc"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8"
 
 [targets]
-test = ["Aqua", "BandedMatrices", "BlockBandedMatrices", "Documenter", "DualNumbers", "FFTW", "ParallelTestRunner", "Random", "Test"]
+test = ["Aqua", "BandedMatrices", "BlockBandedMatrices", "Documenter", "DualNumbers", "FFTW", "ParallelTestRunner", "Random", "Test", "VectorInterface"]

ext/ApproxFunVectorInterfaceExt.jl

@@ -0,0 +1,86 @@
+module ApproxFunVectorInterfaceExt
+
+using ApproxFun: Fun, space, coefficients, innerproduct
+using ApproxFun.ApproxFunBase: spacescompatible
+using VectorInterface: VectorInterface, One, scalartype, zerovector, zerovector!,
+    scale, scale!, scale!!, add, add!
+
+# Funs in compatible spaces share a coefficient basis, so the vector operations
+# can act directly on the (possibly differently truncated) coefficient vectors.
+# The in-place (`!`) variants require compatible spaces and mutable coefficient
+# vectors; the `!!` variants fall back to the out-of-place versions otherwise.
+
+VectorInterface.scalartype(::Type{<:Fun{<:Any,T}}) where {T} = scalartype(T)
+
+function _check_compatible(y::Fun, x::Fun)
+    spacescompatible(y, x) || throw(ArgumentError(
+        "in-place operations require Funs in compatible spaces, got $(space(y)) and $(space(x))"))
+    return nothing
+end
+
+_resizable(c) = c isa Vector
+
+function VectorInterface.zerovector(f::Fun, ::Type{S}) where {S<:Number}
+    return Fun(space(f), zerovector(coefficients(f), S))
+end
+VectorInterface.zerovector!(f::Fun) = (zerovector!(coefficients(f)); f)
+function VectorInterface.zerovector!!(f::Fun)
+    c = coefficients(f)
+    c′ = VectorInterface.zerovector!!(c)
+    return c′ === c ? f : Fun(space(f), c′)
+end
+
+VectorInterface.scale(f::Fun, α::Number) = Fun(space(f), scale(coefficients(f), α))
+VectorInterface.scale!(f::Fun, α::Number) = (scale!(coefficients(f), α); f)
+function VectorInterface.scale!!(f::Fun, α::Number)
+    c = coefficients(f)
+    c′ = scale!!(c, α)
+    return c′ === c ? f : Fun(space(f), c′)
+end
+
+function VectorInterface.scale!(y::Fun, x::Fun, α::Number)
+    _check_compatible(y, x)
+    cy, cx = coefficients(y), coefficients(x)
+    length(cy) == length(cx) || resize!(cy, length(cx))
+    cy .= scale.(cx, α)
+    return y
+end
+function VectorInterface.scale!!(y::Fun, x::Fun, α::Number)
+    if spacescompatible(y, x) &&
+            (length(coefficients(y)) == length(coefficients(x)) || _resizable(coefficients(y))) &&
+            VectorInterface.promote_scale(x, α) <: scalartype(y)
+        return scale!(y, x, α)
+    end
+    return scale(x, α)
+end
+
+VectorInterface.add(y::Fun, x::Fun, α::Number, β::Number) = scale(y, β) + scale(x, α)
+
+function VectorInterface.add!(y::Fun, x::Fun, α::Number, β::Number)
+    _check_compatible(y, x)
+    cy, cx = coefficients(y), coefficients(x)
+    ny, nx = length(cy), length(cx)
+    if ny < nx
+        resize!(cy, nx)
+        for i in (ny + 1):nx
+            @inbounds cy[i] = scale(cx[i], α)
+        end
+        add!(view(cy, 1:ny), view(cx, 1:ny), α, β)
+    else
+        add!(view(cy, 1:nx), cx, α, β)
+        β === One() || scale!(view(cy, (nx + 1):ny), β)
+    end
+    return y
+end
+function VectorInterface.add!!(y::Fun, x::Fun, α::Number, β::Number)
+    if spacescompatible(y, x) &&
+            (length(coefficients(y)) >= length(coefficients(x)) || _resizable(coefficients(y))) &&
+            VectorInterface.promote_add(y, x, α, β) <: scalartype(y)
+        return add!(y, x, α, β)
+    end
+    return add(y, x, α, β)
+end
+
+VectorInterface.inner(f::Fun, g::Fun) = innerproduct(f, g)
+
+end

test/VectorInterfaceTest.jl

@@ -0,0 +1,151 @@
+module VectorInterfaceTest
+
+using ApproxFun
+using Test
+using VectorInterface
+using VectorInterface: One, Zero
+
+include(joinpath(@__DIR__, "testutils.jl"))
+
+@verbose @testset "VectorInterface" begin
+    f = Fun(x -> exp(x) * cos(x))
+    g = Fun(x -> sin(3x))
+    pts = range(-1, 1; length=10)
+
+    @testset "scalartype" begin
+        @test scalartype(f) == Float64
+        @test scalartype(typeof(f)) == Float64
+        @test scalartype(Fun(x -> cis(x), Fourier())) == ComplexF64
+        @test scalartype(Fun(Chebyshev(), Float32[0, 1])) == Float32
+    end
+
+    @testset "zerovector" begin
+        z = zerovector(f)
+        @test z isa Fun
+        @test space(z) == space(f)
+        @test iszero(coefficients(z))
+        @test scalartype(z) == Float64
+
+        z = zerovector(f, ComplexF64)
+        @test iszero(coefficients(z))
+        @test scalartype(z) == ComplexF64
+
+        h = copy(f)
+        @test zerovector!(h) === h
+        @test iszero(coefficients(h))
+
+        h = copy(f)
+        @test zerovector!!(h) === h
+        @test iszero(coefficients(h))
+    end
+
+    @testset "scale" begin
+        h = scale(f, 2.5)
+        @test h(0.3) ≈ 2.5 * f(0.3)
+        @test coefficients(f) == coefficients(Fun(x -> exp(x) * cos(x)))
+
+        h = scale(f, 2im)
+        @test scalartype(h) == ComplexF64
+        @test h(0.3) ≈ 2im * f(0.3)
+
+        h = copy(f)
+        @test scale!(h, 3.0) === h
+        @test h(0.3) ≈ 3 * f(0.3)
+        @test scale!(h, One()) === h
+        @test h(0.3) ≈ 3 * f(0.3)
+
+        h = copy(f)
+        @test scale!!(h, 3.0) === h
+        @test h(0.3) ≈ 3 * f(0.3)
+
+        h = copy(f)
+        h2 = scale!!(h, im)
+        @test h2 !== h
+        @test scalartype(h2) == ComplexF64
+        @test h2(0.3) ≈ im * f(0.3)
+
+        h = copy(f)
+        @test scale!(h, g, 2.0) === h
+        @test ncoefficients(h) == ncoefficients(g)
+        @test h(0.3) ≈ 2 * g(0.3)
+
+        h = copy(f)
+        @test scale!!(h, g, 2.0) === h
+        @test h(0.3) ≈ 2 * g(0.3)
+
+        h = copy(f)
+        h2 = scale!!(h, g, im)
+        @test h2 !== h
+        @test h2(0.3) ≈ im * g(0.3)
+
+        @test_throws ArgumentError scale!(copy(f), Fun(sin, 0..1), 2.0)
+    end
+
+    @testset "add" begin
+        for (α, β) in ((2.0, 3.0), (2.0, One()), (One(), One()), (0.5, Zero()))
+            h = add(copy(f), g, α, β)
+            @test all(x -> isapprox(h(x), β * f(x) + α * g(x); atol=1e-12), pts)
+
+            h = copy(f)
+            @test add!(h, g, α, β) === h
+            @test all(x -> isapprox(h(x), β * f(x) + α * g(x); atol=1e-12), pts)
+
+            h = copy(g)
+            @test add!(h, f, α, β) === h
+            @test all(x -> isapprox(h(x), β * g(x) + α * f(x); atol=1e-12), pts)
+
+            h = copy(f)
+            @test add!!(h, g, α, β) === h
+            @test all(x -> isapprox(h(x), β * f(x) + α * g(x); atol=1e-12), pts)
+        end
+
+        h = add(copy(f), g)
+        @test all(x -> isapprox(h(x), f(x) + g(x); atol=1e-12), pts)
+        h = add(copy(f), g, 2.0)
+        @test all(x -> isapprox(h(x), f(x) + 2g(x); atol=1e-12), pts)
+        h = add!(copy(f), g)
+        @test all(x -> isapprox(h(x), f(x) + g(x); atol=1e-12), pts)
+
+        h = copy(f)
+        h2 = add!!(h, g, im, 1.0)
+        @test h2 !== h
+        @test scalartype(h2) == ComplexF64
+        @test all(x -> isapprox(h2(x), f(x) + im * g(x); atol=1e-12), pts)
+
+        f2 = Fun(x -> x^2, 0..1)
+        g2 = Fun(x -> cospi(2x), Fourier(0..1))
+        @test_throws ArgumentError add!(copy(f2), g2, 1.0, 1.0)
+        h = add!!(copy(f2), g2, 2.0, 3.0)
+        @test all(x -> isapprox(h(x), 3 * f2(x) + 2 * g2(x); atol=1e-10), range(0.1, 0.9; length=5))
+        h = add(copy(f2), g2, 2.0, 3.0)
+        @test all(x -> isapprox(h(x), 3 * f2(x) + 2 * g2(x); atol=1e-10), range(0.1, 0.9; length=5))
+    end
+
+    @testset "inner and norm" begin
+        @test inner(f, g) ≈ sum(f * g)
+        @test norm(f) ≈ sqrt(inner(f, f))
+
+        u = Fun(x -> cis(x), Fourier())
+        v = Fun(x -> sin(x) * cis(2x), Fourier())
+        @test inner(u, v) ≈ sum(conj(u) * v)
+        @test inner(v, u) ≈ conj(inner(u, v))
+        @test norm(u) ≈ sqrt(real(inner(u, u)))
+    end
+
+    @testset "Gram-Schmidt orthogonalization" begin
+        basis = [Fun(x -> x^k) for k in 0:3]
+        ortho = typeof(f)[]
+        for b in basis
+            v = copy(b)
+            for q in ortho
+                v = add!!(v, q, -inner(q, v) / inner(q, q), One())
+            end
+            push!(ortho, scale!!(v, 1 / norm(v)))
+        end
+        for (i, p) in enumerate(ortho), (j, q) in enumerate(ortho)
+            @test inner(p, q) ≈ (i == j ? 1.0 : 0.0) atol=1e-10
+        end
+    end
+end
+
+end # module