Accept vectors in PiecewiseSpace

src/ApproxFunBase.jl

@@ -162,6 +162,28 @@ function _IteratorSize(::Type{T}) where {T<:Tuple}
     any(x -> x isa Base.IsInfinite, s) ? Base.IsInfinite() : Base.HasLength()
 end
 
+_IteratorSize(::Type{<:AbstractVector{<:AbstractFill{<:Any,<:Any,
+        <:Tuple{Vararg{InfRanges}}}}}) = Base.IsInfinite()
+_IteratorSize(::Type{<:AbstractVector{<:InfRanges}}) = Base.IsInfinite()
+_IteratorSize(::Type{<:AbstractVector{<:Union{Vector{<:Integer},
+                                    SVector{<:Any,<:Integer}}}}) = Base.HasLength()
+_IteratorSize(::Type{<:AbstractVector}) = Base.SizeUnknown()
+
+tuple_to_SVector(t::Tuple) = SVector(t)
+tuple_to_SVector(x) = x
+_vcat_toabsvec(args...) = mapreduce(tuple_to_SVector, vcat, args)
+
+# Concatenation that preserves how the pieces are stored. Each argument is either a
+# container of spaces or a single space. Tuples stay tuples, so the number of pieces and
+# the type of each one remain in the type; as soon as a runtime-length vector is involved
+# the result is a vector.
+_ascontainer(x::Tuple) = x
+_ascontainer(x::AbstractVector) = x
+_ascontainer(x) = (x,)
+_vcat_containers(a::Tuple, b::Tuple) = (a..., b...)
+_vcat_containers(a, b) = _vcat_toabsvec(a, b)
+_vcat_preservecontainer(a, b) = _vcat_containers(_ascontainer(a), _ascontainer(b))
+
 include("LinearAlgebra/LinearAlgebra.jl")
 include("Fun.jl")
 include("Domains/Domains.jl")

src/LinearAlgebra/helper.jl

@@ -673,15 +673,13 @@ conv(x::AbstractFill, y::AbstractFill) = DSP.conv(x, y)
 # TODO: cache sums
 
 
-struct BlockInterlacer{DMS<:Tuple{Vararg{AbstractVector{Int}}}}
+struct BlockInterlacer{DMS}
     blocks::DMS
 end
 
 
 const TrivialInterlacer{d,Ax} = BlockInterlacer{<:NTuple{d,Ones{Int,1,Tuple{Ax}}}}
 
-BlockInterlacer(v::AbstractVector) = BlockInterlacer(Tuple(v))
-
 eltype(::Type{<:BlockInterlacer}) = Tuple{Int,Int}
 
 dimensions(b::BlockInterlacer) = map(sum,b.blocks)
@@ -708,10 +706,16 @@ function done(it::BlockInterlacer,st)
     return true
 end
 
-iterate(it::BlockInterlacer) =
+iterate(it::BlockInterlacer{<:Tuple}) =
     iterate(it, (1,1,ntuple(_ -> tuple(), length(it.blocks)),
             ntuple(zero,length(it.blocks))))
 
+iterate(it::BlockInterlacer{<:AbstractVector}) =
+    iterate(it, (1,1, [tuple() for _ in 1:length(it.blocks)],
+            [zero(i) for i in 1:length(it.blocks)]))
+
+_setindex(coll, v, ind) = (coll[ind] = v; coll)
+_setindex(coll::Tuple, v, ind) = Base.setindex(coll, v, ind)
 function iterate(it::BlockInterlacer, (N,k,blkst,lngs))
     done(it, (N,k,blkst,lngs)) && return nothing
 
@@ -731,7 +735,7 @@ function iterate(it::BlockInterlacer, (N,k,blkst,lngs))
         return iterate(it,(N+1,1,blkst,lngs))
     end
 
-    lngs = Base.setindex(lngs, lngs[N]+1, N)
+    lngs = _setindex(lngs, lngs[N]+1, N)
     return (N,lngs[N]),(N,k+1,blkst,lngs)
 end
 

src/Operators/banded/PermutationOperator.jl

@@ -35,8 +35,9 @@ end
 
 # the permutation that rearranges a to be b
 multiplyperm(b,a) = Int[a[bk] for bk in b]
-perm(a::Vector,b::Vector) = multiplyperm(invperm(sortperm(b)),sortperm(a))
-perm(a::Tuple,b::Tuple) = perm(collect(a),collect(b))
+# a and b may be stored either as tuples or as vectors, and the two may be mixed
+perm(a::Union{Tuple,AbstractVector}, b::Union{Tuple,AbstractVector}) =
+    multiplyperm(invperm(sortperm(convert_vector(b))), sortperm(convert_vector(a)))
 
 
 struct NegateEven{T,DS,RS} <: Operator{T}

src/Operators/general/InterlaceOperator.jl

@@ -229,11 +229,21 @@ InterlaceOperator(ops::AbstractArray, ds=NoSpace, rs=ds) =
     InterlaceOperator(Array{Operator{promote_eltypeof(ops)}, ndims(ops)}(ops), ds, rs)
 
 
+# `map` over a `Diagonal` runs the function on the structural zeros as well, which
+# for operator eltypes either errors or silently drops the diagonal structure,
+# depending on how the diagonal is stored. Map over the diagonal itself instead.
+_mapops(f, ops) = map(f, ops)
+_mapops(f, ops::Diagonal) = Diagonal(map(f, parent(ops)))
+
+# a dense array of operators, with the structural zeros of a `Diagonal` filled in
+_denseops(ops) = ops
+_denseops(ops::Diagonal) = [ops[k,j] for k in axes(ops,1), j in axes(ops,2)]
+
 function convert(::Type{Operator{T}},S::InterlaceOperator) where T
     if T == eltype(S)
         S
     else
-        ops = map(x -> convert(Operator{T},x), S.ops)
+        ops = _mapops(x -> convert(Operator{T},x), S.ops)
         InterlaceOperator(ops,domainspace(S),rangespace(S),
                             S.domaininterlacer,S.rangeinterlacer,S.bandwidths,
                             S.blockbandwidths, S.israggedbelow)
@@ -428,7 +438,7 @@ function blockbanded_interlace_convert!(S,ret)
                 KR_size = Block.(Int(first(KR)):min(Int(last(KR)),blocksize(op,1)))
                 JR_size = Block.(Int(first(JR)):min(Int(last(JR)),blocksize(op,2)))
                 BlockBandedMatrix(view(op, KR_size, JR_size))
-            end, parent(S).ops)
+            end, _denseops(parent(S).ops))
 
     for J=blockaxes(ret,2),K=blockcolrange(ret,J)
         Bs=view(ret,K,J)

src/Spaces/ProductSpaceOperators.jl

@@ -187,6 +187,14 @@ end
 
 
 
+# combine the pieces of two direct sums one by one, keeping tuples as tuples
+_mapspaces(f,S1::Tuple,S2::Tuple) = map(f,S1,S2)
+_mapspaces(f,S1,S2) = [f(S1[k],S2[k]) for k=1:length(S1)]
+
+# The pieces of the two spaces may be stored in different containers, and a tuple never
+# compares equal to a vector however its elements compare, so compare the pieces.
+_samespaces(S1,S2) = length(S1) == length(S2) && all(((a,b),) -> a == b, zip(S1,S2))
+
 for (OPrule,OP) in ((:conversion_rule,:conversion_type),(:maxspace_rule,:maxspace),
                         (:union_rule,:union))
     for TYP in (:SumSpace,:PiecewiseSpace)
@@ -198,11 +206,11 @@ for (OPrule,OP) in ((:conversion_rule,:conversion_type),(:maxspace_rule,:maxspac
                 NoSpace()
             elseif canonicalspace(S1sp) == canonicalspace(S2sp)  # this sorts S1 and S2
                 S1sp ≤ S2sp ? S1sp : S2sp  # choose smallest space by sorting
-            elseif cs1 == cs2
+            elseif _samespaces(cs1,cs2)
                 # we can just map down
                 # $TYP(map($OP,S1.spaces,S2.spaces))
                 # this is commented out due to Issue #13261
-                newspaces = [$OP(S1[k],S2[k]) for k=1:length(S1)]
+                newspaces = _mapspaces($OP,S1,S2)
                 if any(b->b==NoSpace(),newspaces)
                     NoSpace()
                 else
@@ -211,7 +219,8 @@ for (OPrule,OP) in ((:conversion_rule,:conversion_type),(:maxspace_rule,:maxspac
             elseif sort(collect(cs1)) == sort(collect(cs2))
                 # sort S1
                 p=perm(cs1,cs2)
-                $OP($TYP(S1[p]),S2sp)
+                # a permutation that reorders nothing would recurse forever
+                p == 1:length(S1) ? NoSpace() : $OP($TYP(S1[p]),S2sp)
             elseif length(S1) == length(S2) == 2  &&
                     $OP(S1[1],S2[1]) != NoSpace() &&
                     $OP(S1[2],S2[2]) != NoSpace()

src/Spaces/SumSpace.jl

@@ -50,49 +50,67 @@ end
 SumSpace(sp::Tuple) = SumSpace{typeof(sp),domaintype(first(sp)),
                                 mapreduce(rangetype,promote_type,sp)}(sp)
 
+SumSpace(A::SumSpace,B::SumSpace) = SumSpace(tuple(A.spaces...,B.spaces...))
+
+SumSpace(A::Space,B::SumSpace) = SumSpace(tuple(A,B.spaces...))
+SumSpace(A::SumSpace,B::Space) = SumSpace(tuple(A.spaces...,B))
+SumSpace(A::Space...) = SumSpace(A)
+SumSpace(sp::AbstractArray) = SumSpace(tuple(sp...))
+
+canonicalspace(A::SumSpace) = SumSpace(sort(collect(A.spaces)))
+
+# TODO: Fix this Hack
+SumSpace(A::ConstantSpace{AnyDomain}, B::ConstantSpace{AnyDomain}) = error("Should not happen")
+SumSpace(A::SumSpace, B::ConstantSpace{AnyDomain}) = SumSpace(A, setdomain(B, domain(A)))
+SumSpace(B::ConstantSpace{AnyDomain}, A::SumSpace) = SumSpace(setdomain(B, domain(A)), A)
+SumSpace(A::Space, B::ConstantSpace{AnyDomain}) = SumSpace(A, setdomain(B, domain(A)))
+SumSpace(B::ConstantSpace{AnyDomain}, A::Space) = SumSpace(setdomain(B, domain(A)), A)
 
 struct PiecewiseSpace{SV,D<:UnionDomain,R} <: DirectSumSpace{SV,D,R}
     spaces::SV
     PiecewiseSpace{SV,D,R}(dom::AnyDomain) where {SV,D,R} =
         new{SV,D,R}(tuple(map(typ->typ(dom),SV.parameters)...))
     PiecewiseSpace{SV,D,R}(dom::UnionDomain) where {SV,D,R} =
         new{SV,D,R}(tuple(map((typ,dom)->typ(dom),SV.parameters,dom.domains)...))
-    PiecewiseSpace{SV,D,R}(sp::Tuple) where {SV,D,R} =
+    PiecewiseSpace{SV,D,R}(sp::SV) where {SV,D,R} =
         new{SV,D,R}(sp)
 end
 
 function _PiecewiseSpace(sp)
     PiecewiseSpace{typeof(sp),typeof(UnionDomain(map(domain,sp))),
                    mapreduce(rangetype,promote_type,sp)}(sp)
 end
-function PiecewiseSpace(spin::Tuple)
-    sp=tuple(union(spin)...)  # remove duplicates
-    _PiecewiseSpace(sp)
+# `union` removes duplicates, but it always returns a `Vector`, which would drop the
+# number of pieces and their individual types from the type of a tuple of spaces.
+# Duplicates are rare here -- `union(::Space, ::Space)` already short-circuits equal
+# spaces before a `PiecewiseSpace` is built -- so keep the tuple when nothing was
+# removed, and only fall back to the vector that `union` returned when it was.
+_uniquespaces(sp::AbstractVector) = union(sp)
+function _uniquespaces(sp::Tuple)
+    u = union(sp)
+    length(u) == length(sp) ? sp : u
 end
 
-PiecewiseSpace(spin::Set) = PiecewiseSpace(collect(spin))
+function PiecewiseSpace(spacesin::Union{Tuple{Vararg{Space}}, AbstractVector{<:Space}})
+    _PiecewiseSpace(_uniquespaces(spacesin))
+end
 
+PiecewiseSpace(spacesin::Set) = PiecewiseSpace(collect(spacesin))
 
+PiecewiseSpace(A::PiecewiseSpace, B::PiecewiseSpace) =
+    PiecewiseSpace(_vcat_preservecontainer(A.spaces, B.spaces))
 
-for TYP in (:SumSpace,:PiecewiseSpace)
-    @eval begin
-        $TYP(A::$TYP,B::$TYP) = $TYP(tuple(A.spaces...,B.spaces...))
-
-        $TYP(A::Space,B::$TYP) = $TYP(tuple(A,B.spaces...))
-        $TYP(A::$TYP,B::Space) = $TYP(tuple(A.spaces...,B))
-        $TYP(A::Space...) = $TYP(A)
-        $TYP(sp::AbstractArray) = $TYP(tuple(sp...))
+PiecewiseSpace(A::Space,B::PiecewiseSpace) = PiecewiseSpace(_vcat_preservecontainer(A, B.spaces))
+PiecewiseSpace(A::PiecewiseSpace,B::Space) = PiecewiseSpace(_vcat_preservecontainer(A.spaces, B))
+PiecewiseSpace(A::Space...) = PiecewiseSpace(A)
 
-        canonicalspace(A::$TYP) = $TYP(sort(collect(A.spaces)))
-    end
-end
+# `sort` of a tuple is only defined if every piece has the same type, so fall back to
+# permuting by a sorted vector for a mix of spaces. Either way a tuple stays a tuple.
+_sortspaces(sp::AbstractVector) = sort(sp)
+_sortspaces(sp::Tuple{T,Vararg{T}}) where {T} = sort(sp)
+_sortspaces(sp::Tuple) = sp[sortperm(convert_vector(sp))]
 
-# TODO: Fix this Hack
-SumSpace(A::ConstantSpace{AnyDomain}, B::ConstantSpace{AnyDomain}) = error("Should not happen")
-SumSpace(A::SumSpace, B::ConstantSpace{AnyDomain}) = SumSpace(A, setdomain(B, domain(A)))
-SumSpace(B::ConstantSpace{AnyDomain}, A::SumSpace) = SumSpace(setdomain(B, domain(A)), A)
-SumSpace(A::Space, B::ConstantSpace{AnyDomain}) = SumSpace(A, setdomain(B, domain(A)))
-SumSpace(B::ConstantSpace{AnyDomain}, A::Space) = SumSpace(setdomain(B, domain(A)), A)
+canonicalspace(A::PiecewiseSpace) = PiecewiseSpace(_sortspaces(A.spaces))
 
 pieces(sp::PiecewiseSpace) = sp.spaces
 piece(s::Space,k) = pieces(s)[k]
@@ -112,8 +130,16 @@ setdomain(A::PiecewiseSpace,d::UnionDomain) =
 
 
 
-function spacescompatible(A::S,B::S) where S<:DirectSumSpace
-    if ncomponents(A) != ncomponents(B)
+# Two direct sums are compatible if they are the same kind of sum and their pieces
+# match up. The pieces may be stored as a tuple in one and as a vector in the other,
+# so this must not require the two spaces to have the same type.
+_directsumkind(::SumSpace) = SumSpace
+_directsumkind(::PiecewiseSpace) = PiecewiseSpace
+
+function spacescompatible(A::DirectSumSpace,B::DirectSumSpace)
+    if _directsumkind(A) !== _directsumkind(B)
+        false
+    elseif ncomponents(A) != ncomponents(B)
         false
     else
         ret=true
@@ -140,7 +166,18 @@ end
 Base.promote_rule(::Type{Fun{SumSpace{SV,D,R}}},::Type{T}) where {SV,D,R,T<:Number} =
     promote_rule(VFun{SumSpace{SV,D,R},Float64},T)
 
-function Base.promote_rule(::Type{Fun{PiecewiseSpace{SV,D,R},V,VV}},::Type{T}) where {SV,D,R,V,VV,T<:Number}
+# pieces stored in a vector: there is no per-piece type to promote, only the element type
+function Base.promote_rule(::Type{Fun{PiecewiseSpace{SV,D,R},V,VV}},
+        ::Type{T}) where {SV<:AbstractVector,D,R,V,VV,T<:Number}
+    newf = promote_type(VFun{eltype(SV),V},T)
+    if newf == Fun
+        Fun
+    else
+        VFun{PiecewiseSpace{Vector{newf.parameters[1]},D,R},promote_type(V,T)}
+    end
+end
+
+function Base.promote_rule(::Type{Fun{PiecewiseSpace{SV,D,R},V,VV}},::Type{T}) where {SV<:Tuple,D,R,V,VV,T<:Number}
     # if any doesn't support promoting, just leave unpromoted
 
     newfsp=map(s->promote_type(VFun{s,V},T),SV.parameters)

test/SpacesTest.jl

@@ -1,5 +1,5 @@
 using ApproxFunBase
-using ApproxFunBase: PointSpace, HeavisideSpace, PiecewiseSegment, dimension, SVector, checkpoints, AnyDomain
+using ApproxFunBase: PointSpace, HeavisideSpace, PiecewiseSegment, dimension, SVector, checkpoints, AnyDomain, SumSpace
 using BandedMatrices: rowrange, colrange, BandedMatrix
 using DomainSets: Point
 using LinearAlgebra
@@ -445,4 +445,44 @@ using Test
             @test values(f) == f.(points(f))
         end
     end
+
+    @testset "PiecewiseSpace" begin
+        ps = PiecewiseSpace([PointSpace(1:2), PointSpace(3:4)])
+        @test PiecewiseSpace(ps, ps) == ps
+        @test ApproxFunBase.canonicalspace(ps) == ps
+        d = domain(ps)
+        @test all(x -> x in d, 1:4)
+        @test all(x -> !(x in d), 5:6)
+        d2 = domain(PiecewiseSpace(ps, PointSpace(5:6)))
+        @test all(x -> x in d2, 5:6)
+        d2 = domain(PiecewiseSpace(PointSpace(5:6), ps))
+        @test all(x -> x in d2, 5:6)
+
+        ps2 = PiecewiseSpace([PointSpace(3:4), PointSpace(1:2)])
+        @test ApproxFunBase.canonicalspace(ps2) == ps
+
+        # a Set has no order, so canonicalize before comparing
+        @test ApproxFunBase.canonicalspace(
+            PiecewiseSpace(Set([PointSpace(1:2), PointSpace(3:4)]))) == ps
+    end
+
+    @testset "SumSpace" begin
+        a, b, c = PointSpace(1:2), PointSpace(3:4), PointSpace(5:6)
+        ss = SumSpace(a, b)
+        @test components(ss) == (a, b)
+
+        @test components(SumSpace(ss, c)) == (a, b, c)
+        @test components(SumSpace(c, ss)) == (c, a, b)
+        @test components(SumSpace(ss, SumSpace(c, a))) == (a, b, c, a)
+        @test SumSpace([a, b]) == ss
+
+        @test ApproxFunBase.canonicalspace(SumSpace(b, a)) == ApproxFunBase.canonicalspace(ss)
+
+        # a ConstantSpace{AnyDomain} argument inherits the domain of the other space
+        @test domain(components(SumSpace(ss, ConstantSpace()))[end]) == domain(ss)
+        @test domain(components(SumSpace(a,  ConstantSpace()))[end]) == domain(a)
+        @test domain(components(SumSpace(ConstantSpace(), ss))[1])   == domain(ss)
+        @test domain(components(SumSpace(ConstantSpace(), a))[1])    == domain(a)
+        @test_throws ErrorException SumSpace(ConstantSpace(), ConstantSpace())
+    end
 end

test/runtests.jl

@@ -14,6 +14,7 @@ using LinearAlgebra
 using LowRankMatrices
 using Random
 using SpecialFunctions
+using StaticArrays
 using Test
 
 @testset "Project quality" begin
@@ -86,6 +87,31 @@ end
             @test first(C, 10) == C[1:10] == B[1:10] == first(B, 10)
             @test C[2:10][1:2:end] == B[2:10][1:2:end]
         end
+
+        a = [1,2]
+        sa = SVector(1,2)
+        f = Fill(2,∞)
+        @test Base.IteratorSize(ApproxFunBase.BlockInterlacer((f, f))) == Base.IsInfinite()
+        @test Base.IteratorSize(ApproxFunBase.BlockInterlacer([f, f])) == Base.IsInfinite()
+        @test Base.IteratorSize(ApproxFunBase.BlockInterlacer((1:∞, 1:∞))) == Base.IsInfinite()
+        @test Base.IteratorSize(ApproxFunBase.BlockInterlacer([1:∞, 1:∞])) == Base.IsInfinite()
+        @test Base.IteratorSize(ApproxFunBase.BlockInterlacer([sa, sa])) == Base.HasLength()
+        @test Base.IteratorSize(ApproxFunBase.BlockInterlacer([a, a])) == Base.HasLength()
+        @test Base.IteratorSize(ApproxFunBase.BlockInterlacer((a, a))) == Base.HasLength()
+
+        b1 = ApproxFunBase.BlockInterlacer([Fill(2,2), 1:2])
+        b2 = ApproxFunBase.BlockInterlacer((Fill(2,2), 1:2))
+        @test collect(b1) == collect(b2)
+
+        # a bare Space is one component, so interlacer must wrap its
+        # blocklengths in a 1-element collection
+        ps = ApproxFunBase.PointSpace(1:3)
+        @test ApproxFunBase.interlacer(ps).blocks == (blocklengths(ps),)
+        cs = ApproxFunBase.ContinuousSpace(ApproxFunBase.PiecewiseSegment([1.0,2.0,3.0]))
+        itc = ApproxFunBase.interlacer(cs)
+        @test itc.blocks == (blocklengths(cs),)
+        @test Base.IteratorSize(itc) == Base.IsInfinite()
+        @test collect(Iterators.take(itc, 4)) == [(1,k) for k in 1:4]
     end
 
     @testset "issue #94" begin