Merge pull request #3 from JuliaApproximation/dl/3d

3D support

Author: dlfivefifty

Project.toml

@@ -43,7 +43,7 @@ FastTransforms = "≥ 0.4.2"
 FillArrays = "≥ 0.5.0"
 InfiniteArrays = "0.1"
 IntervalSets = "≥ 0.3.1"
-LazyArrays = "0.8.0"
+LazyArrays = "0.8, 0.9"
 LowRankApprox = "0.2"
 SpecialFunctions = "≥ 0.7.0"
 StaticArrays = "≥ 0.8.3"

src/Domains/ProductDomain.jl

@@ -19,15 +19,15 @@ function pushappendpts!(ret, xx, pts)
         push!(ret,Vec(xx...))
     else
         for x in pts[1]
-            pushappendpts!(ret,(xx...,x),pts[2:end])
+            pushappendpts!(ret,(xx...,x...),pts[2:end])
         end
     end
     ret
 end
 
 function checkpoints(d::ProductDomain)
-    pts=map(checkpoints,d.domains)
-    ret=Vector{Vec{length(d.domains),float(mapreduce(eltype,promote_type,d.domains))}}(undef, 0)
+    pts = checkpoints.(d.domains)
+    ret=Vector{Vec{sum(dimension.(d.domains)),float(promote_type(eltype.(eltype.(d.domains))...))}}(undef, 0)
 
     pushappendpts!(ret,(),pts)
     ret

src/LinearAlgebra/helper.jl

@@ -716,18 +716,6 @@ conv(y::SVector{1}, x::AbstractVector) = y[1]*x
 conv(x::AbstractFill, y::SVector{1}) = x*y[1]
 conv(y::SVector{1}, x::AbstractFill) = y[1]*x
 conv(x::AbstractFill, y::AbstractFill) = DSP.conv(x, y)
-function conv(x::AbstractFill, y::AbstractVector)
-    isinf(length(x)) || return DSP.conv(x,y)
-    @assert length(y) == 1
-    x*y[1]
-end
-function conv(y::AbstractVector, x::AbstractFill)
-    isinf(length(x)) || return DSP.conv(y,x)
-    @assert length(y) == 1
-    y[1]*x
-end
-
-
 
 
 ## BlockInterlacer

src/Multivariate/TensorSpace.jl

@@ -246,7 +246,7 @@ TensorSpace(A::ProductDomain) = TensorSpace(tuple(map(Space,A.domains)...))
 ⊗(A::Space,B::TensorSpace) = TensorSpace(A,B.spaces...)
 ⊗(A::Space,B::Space) = TensorSpace(A,B)
 
-domain(f::TensorSpace) = mapreduce(domain,×,f.spaces)
+domain(f::TensorSpace) = ×(domain.(f.spaces)...)
 Space(sp::ProductDomain) = TensorSpace(sp)
 
 setdomain(sp::TensorSpace, d::ProductDomain) = TensorSpace(setdomain.(factors(sp), factors(d)))
@@ -368,9 +368,15 @@ function plan_transform(sp::TensorSpace, ::Type{T}, n::Integer) where {T}
     TransformPlan(sp,((plan_transform(sp.spaces[1],T,N),N),
                     (plan_transform(sp.spaces[2],T,M),M)),
                 Val{false})
-end   
+end  
+
+function plan_transform!(sp::TensorSpace, ::Type{T}, n::Integer) where {T}
+    P = plan_transform(sp, T, n)
+    TransformPlan(sp, P.plan, Val{true})
+end
 
 plan_transform(sp::TensorSpace, v::AbstractVector) = plan_transform(sp,eltype(v),length(v))
+plan_transform!(sp::TensorSpace, v::AbstractVector) = plan_transform!(sp,eltype(v),length(v))
 
 function plan_itransform(sp::TensorSpace, v::AbstractVector{T}) where {T}
     N,M = size(totensor(sp, v)) # wasteful
@@ -380,13 +386,13 @@ function plan_itransform(sp::TensorSpace, v::AbstractVector{T}) where {T}
 end   
 
 
-function *(T::TransformPlan{<:Any,<:TensorSpace,true},v::AbstractVector) 
+function *(T::TransformPlan{TT,<:TensorSpace,true},v::AbstractVector) where TT # need where TT
     N,M = T.plan[1][2],T.plan[2][2]
     V=reshape(v,N,M)
     fromtensor(T.space,T*V)
 end
 
-*(T::ITransformPlan{<:Any,<:TensorSpace,true},v::AbstractVector)  =
+*(T::ITransformPlan{TT,<:TensorSpace,true},v::AbstractVector) where TT  =
     vec(T*totensor(T.space,v))
 
 

src/Operators/Operator.jl

@@ -461,6 +461,7 @@ macro wrapperstructure(Wrap)
              $ret
 
              $func(D::$Wrap,k::Integer) = $func(D.op,k)
+             $func(A::$Wrap,i::ApproxFunBase.Infinity) = $func(D.op,k)
          end
      end
 

src/Space.jl

@@ -98,6 +98,7 @@ isambiguous(sp::Space) = isambiguous(rangetype(sp))
 
 #TODO: should it default to canonicalspace?
 points(d::Space,n) = points(domain(d),n)
+points(d::Space) = points(d, dimension(d))
 
 
 

src/constructors.jl

@@ -38,10 +38,10 @@ end
 
 # default_Fun is the default constructor, based on evaluation and transforms
 # last argument is whether to splat or not
-default_Fun(::Type{T},f,d::Space{ReComp},pts::AbstractVector,::Type{Val{true}}) where {T,ReComp} =
+default_Fun(::Type{T},f,d::Space{ReComp},pts::AbstractArray,::Type{Val{true}}) where {T,ReComp} =
     Fun(d,transform(d,T[f(x...) for x in pts]))
 
-default_Fun(::Type{T},f,d::Space{ReComp},pts::AbstractVector,::Type{Val{false}}) where {T,ReComp} =
+default_Fun(::Type{T},f,d::Space{ReComp},pts::AbstractArray,::Type{Val{false}}) where {T,ReComp} =
     Fun(d,transform(d,broadcast!(f, similar(pts, T), pts)))
 
 

test/SpacesTest.jl

@@ -1,5 +1,5 @@
 using ApproxFunBase, Test
-import ApproxFunBase: PointSpace, HeavisideSpace, PiecewiseSegment
+import ApproxFunBase: PointSpace, HeavisideSpace, PiecewiseSegment, dimension, Vec, checkpoints
 
 @testset "Spaces" begin
     @testset "PointSpace" begin
@@ -34,6 +34,10 @@ import ApproxFunBase: PointSpace, HeavisideSpace, PiecewiseSegment
 
         S = HeavisideSpace([-1.0,0.0,1.0])
         @test Derivative(S) === Derivative(S,1)
+
+        a = HeavisideSpace(0:0.25:1)
+        @test dimension(a) == 4
+        @test @inferred(points(a)) == 0.125:0.25:0.875
     end
 
     @testset "DiracDelta integration and differentiation" begin
@@ -57,4 +61,22 @@ import ApproxFunBase: PointSpace, HeavisideSpace, PiecewiseSegment
         @test h(2) == 0.3+1im
         @test h(3) == 3.3+1im
     end
+
+    @testset "Multivariate" begin
+        a = HeavisideSpace(0:0.25:1)
+        @test @inferred(dimension(a^2)) == dimension(a)^2
+        @test @inferred(domain(a^2)) == domain(a)^2
+        @test @inferred(points(a^2)) == vec(Vec.(points(a), points(a)'))
+        @test  @inferred(checkpoints(a^2)) == vec(Vec.(checkpoints(a)', checkpoints(a)))
+        
+        aa2 = TensorSpace(a , a^2)
+        @test dimension(aa2) == dimension(a)^3
+        @test @inferred(domain(aa2)) == domain(a)^3
+        @test @inferred(points(aa2)) == vec(Vec.(points(a), points(a)', reshape(points(a), 1,1,4)))
+        @test  @inferred(checkpoints(aa2)) == vec(Vec.(reshape(checkpoints(a), 1,1,length(checkpoints(a))), checkpoints(a)', checkpoints(a)))
+
+        @test dimension(a^3) == dimension(a)^3
+        @test @inferred(domain(a^3)) == domain(a)^3
+        @test_broken @inferred(points(a^3)) == vec(Vec.(points(a), points(a)', reshape(points(a), 1,1,4)))
+    end
 end