Define values for ArraySpace (#688)

`values(f)` errored for an array-valued `Fun`, which broke `iszero`
since it started consulting `values` in #635. Three separate causes:

* `itransform` was only defined for `VectorSpace`, so a `MatrixSpace`
  fell through to the generic `Space` method and hit
  `checkcanonicalspace`, which errors as an `ArraySpace` of canonical
  spaces is its own canonical space. Widen the method to `ArraySpace`
  and reshape the values to the shape of the space.

* The component coefficients may be views, so `pad!` may not resize
  them. Use `pad` instead.

* `_values` asserted a `Vector{float(T)}` return for any space over a
  numeric domain, which includes array-valued ones. Restrict the
  assertion to spaces with a scalar range type.

Also fix `points(f::Fun{<:ArraySpace})`, which passed the total number
of interlaced coefficients where a per-component count is expected, so
that `values(f) == f.(points(f))`.

Fixes #647

Co-authored-by: Claude Opus 5 (1M context) <noreply@anthropic.com>

Author: jishnub

src/Fun.jl

@@ -426,7 +426,8 @@ true
 """
 values(f::Fun,dat...) = _values(f.space, f.coefficients, dat...)
 _values(sp, v, dat...) = itransform(sp, v, dat...)
-_values(sp::UnivariateSpace, v::Vector{T}, dat...) where {T<:Number} =
+# the return type may be asserted only for scalar-valued spaces
+_values(sp::Space{<:Domain{<:Number},<:Number}, v::Vector{T}, dat...) where {T<:Number} =
     itransform(sp, v, dat...)::Vector{float(T)}
 
 """

src/Spaces/ArraySpace.jl

@@ -74,6 +74,11 @@ setdomain(A::ArraySpace,d::Domain) = ArraySpace(map(sp->setdomain(sp,d),A.spaces
 
 #TODO: rework for different spaces
 points(d::ArraySpace,n) = points(d.spaces[1],n)
+# points(::ArraySpace, n) forwards n to a component space, so n must be the number of
+# coefficients per component. The generic points(f::Fun) would pass ncoefficients(f),
+# which counts the interlaced coefficients of all the components together.
+# The count below matches the padding in itransform, so that values(f) == f.(points(f)).
+points(f::Fun{<:ArraySpace}) = points(space(f), maximum(ncoefficients, vec(f), init=0))
 
 
 transform(AS::ArraySpace{SS,1},vals::AbstractVector{Vector{V}}) where {SS,V} =
@@ -98,10 +103,12 @@ transform(AS::VectorSpace{SS},vals::AbstractVector{AV}) where {SS,AV<:AbstractVe
 transform(AS::VectorSpace{SS},vals::AbstractVector{SVector{V,n}}) where {SS,n,V} =
     transform(AS,map(Vector,vals))
 
-function itransform(AS::VectorSpace,cfs::AbstractVector)
+function itransform(AS::ArraySpace,cfs::AbstractVector)
     vf = vec(Fun(AS, cfs))
     n = maximum(ncoefficients, vf)
-    vcat.(values.(pad!.(vf, n))...)
+    vals = vcat.(values.(pad.(vf, n))...)
+    # the values of an array-valued Fun are arrays of the same shape
+    [reshape(v, size(AS)) for v in vals]
 end
 
 

test/SpacesTest.jl

@@ -395,5 +395,30 @@ using Test
         A = ApproxFunBase.ArraySpace(empty!([PointSpace(1:3)]))
         @test length(A) == 0
         @test ApproxFunBase.dimension(A) == 0
+
+        @testset "values" begin
+            # the two spaces exercise different paths: an ArraySpace of
+            # HeavisideSpaces lies over a numeric domain and has a trivial
+            # interlacer, so that the component coefficients are views
+            @testset for S in (PointSpace(1:3), ApproxFunBase.HeavisideSpace([-1.0, -0.5, 0.0, 1.0]))
+                @testset for sz in ((2,), (2,2))
+                    A = ApproxFunBase.ArraySpace(S, sz...)
+                    # 3 coefficients per component, so that each component
+                    # may be evaluated on the full grid
+                    n = 3length(A)
+                    f = Fun(A, Float64.(1:n))
+                    v = values(f)
+                    @test all(x -> size(x) == sz, v)
+                    @test v == f.(points(f))
+                    @test !iszero(f)
+                    @test iszero(Fun(A, zeros(n)))
+                end
+            end
+
+            # static array space
+            A = ApproxFunBase.ArraySpace(PointSpace(1:3), Val((2,2)))
+            f = Fun(A, Float64.(1:3length(A)))
+            @test values(f) == f.(points(f))
+        end
     end
 end