Initial version of ContinuumArrays, with own implementation of AbstractAxisArray

Author: dlfivefifty

REQUIRE

@@ -0,0 +1,2 @@
+julia 0.7
+IntervalSets

examples/BSplines.jl

@@ -0,0 +1,70 @@
+using ContinuumArrays, LazyArrays, IntervalSets
+import ContinuumArrays: AbstractAxisMatrix, ℵ₀
+import Base: axes, getindex, convert
+
+struct Spline{order,T} <: AbstractAxisMatrix{T}
+    points::Vector{T}
+end
+
+const LinearSpline{T} = Spline{1,T}
+const HeavisideSpline{T} = Spline{0,T}
+
+Spline{o}(pts::AbstractVector{T}) where {o,T} = Spline{o,float(T)}(pts)
+
+axes(B::Spline{o}) where o = (first(B.points)..last(B.points), Base.OneTo(length(B.points)-o-1))
+
+function getindex(B::LinearSpline{T}, x::Real, k::Int) where T
+    x ∈ axes(B,1) && 1 ≤ k ≤ size(B,2)|| throw(BoundsError())
+
+    p = B.points
+    n = length(p)
+
+    k > 1 && x ≤ p[k-1] && return zero(T)
+    k < n && x ≥ p[k+1] && return zero(T)
+    x == p[k] && return one(T)
+    x < p[k] && return (x-p[k-1])/(p[k]-p[k-1])
+    return (x-p[k+1])/(p[k]-p[k+1]) # x ≥ p[k]
+end
+
+function getindex(B::HeavisideSpline{T}, x::Real, k::Int) where T
+    x ∈ axes(B,1) && 1 ≤ k ≤ size(B,2)|| throw(BoundsError())
+
+    p = B.points
+    n = length(p)
+
+    x < p[k] && return zero(T)
+    k < n && x > p[k+1] && return zero(T)
+    return one(T)
+end
+
+# getindex(B::LinearSpline, ::Colon, k::Int) = Mul(B, [Zeros{Int}(k-1); 1; Zeros{Int}(size(B,2)-k)])
+
+# function convert(::Type{SymTridiagonal}, AB::Mul{<:Any,<:Any,<:Any,<:ContinuumArrays.Adjoint{<:Any,<:LinearSpline},<:LinearSpline})
+#     Ac,B = AB.A, AB.B
+#     A = parent(Ac)
+#     @assert A.points == B.points
+#     x = A.points
+#     SymTridiagonal(x, x/2) # TODO fix
+# end
+#
+# materialize(M::Mul{<:Any,<:Any,<:Any,<:ContinuumArrays.Adjoint{<:Any,<:LinearSpline},<:LinearSpline}) =
+#     convert(SymTridiagonal, M)
+
+## tests
+
+B = HeavisideSpline([1,2,3])
+@test size(B) == (ℵ₀, 2)
+
+@test_throws BoundsError B[0.1, 1]
+@test B[1.1,1] === 1.0
+@test B[2.1,1] === 0.0
+@test B[1.1,2] === 0.0
+@test B[2.1,2] === 1.0
+@test_throws BoundsError B[2.1,3]
+@test_throws BoundsError B[3.1,2]
+
+@test all(B[[1.1,2.1], 1] .=== [1.0,0.0])
+@test all(B[1.1,1:2] .=== [1.0,0.0])
+@test all(B[[1.1,2.1], 1:2] .=== [1.0 0.0; 0.0 1.0])
+
+@test_throws BoundsError B[[0.1,2.1], 1]

examples/BSplinesSphericalHarmonics.jl

@@ -1,14 +1,17 @@
 SH = SphericalHarmonics()
-Spl = BSplines(0.0:0.01:1)
-S = SH ⊗ Spl
+CS = CubicSplines(0.0:0.01:1)
+QS = QuadraticSplines(0.0:0.01:1)
+LS = LinearSplines(0.0:0.01:1)
+S = SH ⊗ CS
 
 λ = @. -(1:∞)*((1:∞) +1)
 L = BlockDiagonal(Diagonal.(Fill.(λ,1:2:∞))) # create a Block matrix with growing block sizes
 Δ_s = Mul(SH, L, inv(SH)) # the Laplace–Boltrami operator
-D_r = Mul(Spl, SymTridiagonal(...), inv(Spl))
-R = Mul(Spl, Diagonal(Spl.points), inv(Spl))
+D_r = Mul(QS, SymTridiagonal(...), inv(CS))
+R = Mul(QS, Banded(...), inv(QS)) # Diagonal(Spl.points)
+D̃_r = Mul(LS, SymTridiagonal(...), inv(QS))
 
-Δ = I ⊗ (D_r*R^2*D_r) + Δ_s ⊗ I # creates Laplacian as a Mul(S, ..., inv(S)) operator
+Δ = I ⊗ (D̃_r*R^2*D_r) + Δ_s ⊗ I # creates Laplacian as a Mul(S, ..., inv(S)) operator
 
 # finite_dimensional case
 N = 100
@@ -31,4 +34,8 @@ MPI_f = Mul(S_N, BlockVector{Float64,MPIVector{Float64}}(undef, blocklengths(bac
 MPI_f .= f # automatically generates data remotely, via MPI version of FastTransforms
 
 MPI_v = similar(MPI_f)
+
+C = Mul(SH ⊗ QS, I ⊗ (QS' * LS), inv(SH ⊗ LS))
+MPI_C = ...
+
 MPI_v .= Mul(MPI_Δ_N, MPI_f) # applies the Laplacian and fills in missing information.

src/ContinuumArrays.jl

@@ -0,0 +1,30 @@
+module ContinuumArrays
+using Base, LinearAlgebra, IntervalSets
+import Base: getindex, size, axes, length, ==, isequal, iterate, CartesianIndices, LinearIndices,
+                Indices, IndexStyle, getindex, setindex!, parent, vec, convert, similar, zero,
+                map, eachindex
+import Base: @_inline_meta, DimOrInd, OneTo, @_propagate_inbounds_meta, @_noinline_meta,
+                DimsInteger, error_if_canonical_getindex, @propagate_inbounds, _return_type, _default_type,
+                _maybetail, tail, _getindex, _maybe_reshape, index_ndims, _unsafe_getindex,
+                index_shape, to_shape, unsafe_length, @nloops, @ncall
+
+
+import LinearAlgebra: transpose, adjoint
+
+abstract type AbstractAxisArray{T,N} end
+AbstractAxisVector{T} = AbstractAxisArray{T,1}
+AbstractAxisMatrix{T} = AbstractAxisArray{T,2}
+AbstractAxisVecOrMat{T} = Union{AbstractAxisVector{T}, AbstractAxisMatrix{T}}
+
+struct ℵ₀ <: Number end
+_length(::AbstractInterval) = ℵ₀
+_length(d) = length(d)
+
+size(A::AbstractAxisArray) = _length.(axes(A))
+axes(A::AbstractAxisArray) = error("Override axes for $(typeof(A))")
+
+include("indices.jl")
+include("abstractaxisarray.jl")
+include("adjtrans.jl")
+include("multidimensional.jl")
+end