first FEM solve working

Author: dlfivefifty

src/ContinuumArrays.jl

@@ -8,7 +8,7 @@ import BandedMatrices: AbstractBandedLayout
 
 include("QuasiArrays/QuasiArrays.jl")
 using .QuasiArrays
-import .QuasiArrays: _length, checkindex, Adjoint, Transpose, slice, QSlice
+import .QuasiArrays: cardinality, checkindex, Adjoint, Transpose, slice, QSlice, SubQuasiArray
 
 export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative
 
@@ -20,7 +20,7 @@ struct AlephInfinity{N} <: Integer end
 const ℵ₁ = AlephInfinity{1}()
 
 
-_length(::AbstractInterval) = ℵ₁
+cardinality(::AbstractInterval) = ℵ₁
 *(ℵ::AlephInfinity) = ℵ
 
 

src/QuasiArrays/QuasiArrays.jl

@@ -2,13 +2,14 @@ module QuasiArrays
 using Base, LinearAlgebra, LazyArrays
 import Base: getindex, size, axes, length, ==, isequal, iterate, CartesianIndices, LinearIndices,
                 Indices, IndexStyle, getindex, setindex!, parent, vec, convert, similar, zero,
-                map, eachindex, eltype, first, last
+                map, eachindex, eltype, first, last, firstindex, lastindex
 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, Slice, unalias
 import Base: ViewIndex, Slice, ScalarIndex, RangeIndex, view, viewindexing, ensure_indexable, index_dimsum,
-                check_parent_index_match, reindex, _isdisjoint
+                check_parent_index_match, reindex, _isdisjoint, unsafe_indices,
+                parentindices
 import Base: *, /, \, +, -, inv
 
 import Base.Broadcast: materialize
@@ -25,15 +26,15 @@ AbstractQuasiMatrix{T} = AbstractQuasiArray{T,2}
 AbstractQuasiVecOrMat{T} = Union{AbstractQuasiVector{T}, AbstractQuasiMatrix{T}}
 
 
-_length(d) = length(d)
+cardinality(d) = length(d)
 
-size(A::AbstractQuasiArray) = _length.(axes(A))
+size(A::AbstractQuasiArray) = cardinality.(axes(A))
 axes(A::AbstractQuasiArray) = error("Override axes for $(typeof(A))")
 
 include("indices.jl")
 include("abstractquasiarray.jl")
 include("multidimensional.jl")
-include("subarray.jl")
+include("subquasiarray.jl")
 include("matmul.jl")
 
 include("adjtrans.jl")

src/QuasiArrays/indices.jl

@@ -157,6 +157,7 @@ struct QSlice{T,AX} <: AbstractQuasiVector{T}
 end
 QSlice(axis) = QSlice{eltype(axis),typeof(axis)}(axis)
 QSlice(S::QSlice) = S
+==(A::QSlice, B::QSlice) = A.axis == B.axis
 axes(S::QSlice) = (S,)
 unsafe_indices(S::QSlice) = (S,)
 axes1(S::QSlice) = S
@@ -169,7 +170,10 @@ last(S::QSlice) = last(S.axis)
 size(S::QSlice) = (length(S.axis),)
 length(S::QSlice) = length(S.axis)
 unsafe_length(S::QSlice) = unsafe_length(S.axis)
+cardinality(S::QSlice) = cardinality(S.axis)
 getindex(S::QSlice, i::Real) = (@_inline_meta; @boundscheck checkbounds(S, i); i)
 getindex(S::QSlice, i::AbstractVector{<:Real}) = (@_inline_meta; @boundscheck checkbounds(S, i); i)
 show(io::IO, r::QSlice) = print(io, "QSlice(", r.axis, ")")
 iterate(S::QSlice, s...) = iterate(S.axis, s...)
+
+checkindex(::Type{Bool}, inds::QSlice, i::Real) = i ∈ inds.axis

src/QuasiArrays/matmul.jl

@@ -19,6 +19,16 @@ function getindex(M::QuasiMatMulVec, k::Real)
     ret
 end
 
+function getindex(M::QuasiMatMulVec, k::AbstractArray)
+    A,B = M.factors
+    ret = zeros(eltype(M),length(k))
+    @inbounds for j = 1:size(A,2)
+        ret .+= Mul(view(A,k,j) , B[j])
+    end
+    ret
+end
+
+
 QuasiMatMulMat{styleA, styleB, T, V} = QuasiArrayMulArray{styleA, styleB, 2, 2, T, V}
 QuasiMatMulQuasiMat{styleA, styleB, T, V} = QuasiArrayMulQuasiArray{styleA, styleB, 2, 2, T, V}
 
@@ -28,5 +38,5 @@ QuasiMatMulQuasiMat{styleA, styleB, T, V} = QuasiArrayMulQuasiArray{styleA, styl
 *(A::AbstractArray, B::AbstractQuasiArray) = materialize(Mul(A,B))
 inv(A::AbstractQuasiArray) = materialize(Inv(A))
 
-*(A::AbstractQuasiArray, B::Mul) = Mul(A, B.factors...)
-*(A::Mul, B::AbstractQuasiArray) = Mul(A.factors..., B)
+*(A::AbstractQuasiArray, B::Mul) = materialize(Mul(A, B.factors...))
+*(A::Mul, B::AbstractQuasiArray) = materialize(Mul(A.factors..., B))

src/QuasiArrays/multidimensional.jl

@@ -7,6 +7,11 @@ end
 @inline to_indices(A, inds, I::Tuple{Colon, Vararg{Any}}) =
     (uncolon(inds, I), to_indices(A, _maybetail(inds), tail(I))...)
 
+@inline index_dimsum(::AbstractQuasiArray{Bool}, I...) = (true, index_dimsum(I...)...)
+@inline function index_dimsum(::AbstractQuasiArray{<:Any,N}, I...) where N
+    (ntuple(x->true, Val(N))..., index_dimsum(I...)...)
+end
+
 slice(d::AbstractVector) = Slice(d)
 slice(d) = QSlice(d)
 

src/QuasiArrays/subarray.jl renamed to src/QuasiArrays/subquasiarray.jl

@@ -33,7 +33,7 @@ check_parent_index_match(parent::AbstractQuasiArray{T,N}, ::NTuple{N, Bool}) whe
 viewindexing(I::Tuple{AbstractQuasiArray, Vararg{Any}}) = IndexCartesian()
 
 # Simple utilities
-size(V::SubQuasiArray) = (@_inline_meta; map(n->unsafe_length(n), axes(V)))
+size(V::SubQuasiArray) = (@_inline_meta; map(n->cardinality(n), axes(V)))
 
 similar(V::SubQuasiArray, T::Type, dims::Dims) = similar(V.parent, T, dims)
 
@@ -84,6 +84,22 @@ end
 unsafe_view(V::SubQuasiArray, I::Vararg{ViewIndex,N}) where {N} =
     (@_inline_meta; _maybe_reindex(V, I))
 
+_maybe_reindex(V, I) = (@_inline_meta; _maybe_reindex(V, I, I))
+# _maybe_reindex(V, I, ::Tuple{AbstractArray{<:AbstractCartesianIndex}, Vararg{Any}}) =
+#     (@_inline_meta; SubArray(V, I))
+# # But allow arrays of CartesianIndex{1}; they behave just like arrays of Ints
+# _maybe_reindex(V, I, A::Tuple{AbstractArray{<:AbstractCartesianIndex{1}}, Vararg{Any}}) =
+#     (@_inline_meta; _maybe_reindex(V, I, tail(A)))
+_maybe_reindex(V, I, A::Tuple{Any, Vararg{Any}}) = (@_inline_meta; _maybe_reindex(V, I, tail(A)))
+
+_subarray(A::AbstractArray, idxs) = SubArray(A, idxs)
+_subarray(A::AbstractQuasiArray, idxs) = SubQuasiArray(A, idxs)
+
+function _maybe_reindex(V, I, ::Tuple{})
+    @_inline_meta
+    @inbounds idxs = to_indices(V.parent, reindex(V, V.indices, I))
+    _subarray(V.parent, idxs)
+end
 ## Re-indexing is the heart of a view, transforming A[i, j][x, y] to A[i[x], j[y]]
 #
 # Recursively look through the heads of the parent- and sub-indices, considering
@@ -115,42 +131,42 @@ end
 
 # In general, we simply re-index the parent indices by the provided ones
 SlowSubQuasiArray{T,N,P,I} = SubQuasiArray{T,N,P,I,false}
-function getindex(V::SlowSubQuasiArray{T,N}, I::Vararg{Int,N}) where {T,N}
+function getindex(V::SlowSubQuasiArray{T,N}, I::Vararg{Real,N}) where {T,N}
     @_inline_meta
     @boundscheck checkbounds(V, I...)
     @inbounds r = V.parent[reindex(V, V.indices, I)...]
     r
 end
 
 FastSubQuasiArray{T,N,P,I} = SubQuasiArray{T,N,P,I,true}
-function getindex(V::FastSubQuasiArray, i::Int)
+function getindex(V::FastSubQuasiArray, i::Real)
     @_inline_meta
     @boundscheck checkbounds(V, i)
     @inbounds r = V.parent[V.offset1 + V.stride1*i]
     r
 end
 # We can avoid a multiplication if the first parent index is a Colon or AbstractUnitRange
 FastContiguousSubQuasiArray{T,N,P,I<:Tuple{Union{Slice, AbstractUnitRange}, Vararg{Any}}} = SubQuasiArray{T,N,P,I,true}
-function getindex(V::FastContiguousSubQuasiArray, i::Int)
+function getindex(V::FastContiguousSubQuasiArray, i::Real)
     @_inline_meta
     @boundscheck checkbounds(V, i)
     @inbounds r = V.parent[V.offset1 + i]
     r
 end
 
-function setindex!(V::SlowSubQuasiArray{T,N}, x, I::Vararg{Int,N}) where {T,N}
+function setindex!(V::SlowSubQuasiArray{T,N}, x, I::Vararg{Real,N}) where {T,N}
     @_inline_meta
     @boundscheck checkbounds(V, I...)
     @inbounds V.parent[reindex(V, V.indices, I)...] = x
     V
 end
-function setindex!(V::FastSubQuasiArray, x, i::Int)
+function setindex!(V::FastSubQuasiArray, x, i::Real)
     @_inline_meta
     @boundscheck checkbounds(V, i)
     @inbounds V.parent[V.offset1 + V.stride1*i] = x
     V
 end
-function setindex!(V::FastContiguousSubQuasiArray, x, i::Int)
+function setindex!(V::FastContiguousSubQuasiArray, x, i::Real)
     @_inline_meta
     @boundscheck checkbounds(V, i)
     @inbounds V.parent[V.offset1 + i] = x
@@ -239,7 +255,7 @@ _pointer(V::SubQuasiArray, i::Int) = pointer(V, Base._ind2sub(axes(V), i))
 axes(S::SubQuasiArray) = (@_inline_meta; _indices_sub(S, S.indices...))
 _indices_sub(S::SubQuasiArray) = ()
 _indices_sub(S::SubQuasiArray, ::Real, I...) = (@_inline_meta; _indices_sub(S, I...))
-function _indices_sub(S::SubQuasiArray, i1::AbstractQuasiArray, I...)
+function _indices_sub(S::SubQuasiArray, i1::Union{AbstractQuasiArray,AbstractArray}, I...)
     @_inline_meta
     (unsafe_indices(i1)..., _indices_sub(S, I...)...)
 end

src/bases/splines.jl

@@ -86,8 +86,8 @@ function copyto!(dest::Mul2{<:Any,<:Any,<:HeavisideSpline},
     bandwidths(A) == (0,1) || throw(ArgumentError("Not implemented"))
 
     d = diff(x)
-    A[band(0)] .= (-).(d)
-    A[band(1)] .= d
+    A[band(0)] .= inv.((-).(d))
+    A[band(1)] .= inv.(d)
 
     dest
 end

src/operators.jl

@@ -20,8 +20,8 @@ end
 
 
 function materialize(M::Mul2{<:Any,<:Any,<:QuasiArrays.Adjoint{<:Any,<:DiracDelta},<:AbstractQuasiVector})
-    A, B = M.A, M.B
-    axes(A,2) == axes(A,1) || throw(DimensionMismatch())
+    A, B = M.factors
+    axes(A,2) == axes(B,1) || throw(DimensionMismatch())
     B[parent(A).x]
 end
 
@@ -31,7 +31,7 @@ function materialize(M::Mul2{<:Any,<:Any,<:QuasiArrays.Adjoint{<:Any,<:DiracDelt
     B[parent(A).x,:]
 end
 
-struct Derivative{T,A} <: AbstractQuasiVector{T}
+struct Derivative{T,A} <: AbstractQuasiMatrix{T}
     axis::A
 end
 
@@ -40,3 +40,11 @@ Derivative(axis) = Derivative{Float64}(axis)
 
 axes(D::Derivative) = (D.axis, D.axis)
 ==(a::Derivative, b::Derivative) = a.axis == b.axis
+
+
+function materialize(M::Mul2{<:Any,<:Any,<:Derivative,<:SubQuasiArray})
+    A, B = M.factors
+    axes(A,2) == axes(B,1) || throw(DimensionMismatch())
+    P = parent(B)
+    (Derivative(axes(P,1))*P)[parentindices(B)...]
+end

test/runtests.jl

@@ -1,5 +1,6 @@
 using ContinuumArrays, LazyArrays, IntervalSets, FillArrays, LinearAlgebra, BandedMatrices, Test
     import ContinuumArrays: ℵ₁, materialize
+    import ContinuumArrays.QuasiArrays: SubQuasiArray
 
 @testset "DiracDelta" begin
     δ = DiracDelta(-1..3)
@@ -31,7 +32,6 @@ end
 
     @test_throws BoundsError H[[0.1,2.1], 1]
 
-
     f = H*[1,2]
     @test axes(f) == (1.0..3.0,)
     @test f[1.1] ≈ 1
@@ -89,18 +89,107 @@ end
     DL = D*L
     @test M.factors == tuple(D', (D*L).factors...)
 
-    @test materialize(Mul(L', D', D, L)) == materialize(L'D'*D*L) ==
+    @test materialize(Mul(L', D', D, L)) == (L'D'*D*L) ==
         [1.0 -1 0; -1.0 2.0 -1.0; 0.0 -1.0 1.0]
 
     @test materialize(Mul(L', D', D, L)) isa BandedMatrix
-    @test materialize(L'D'*D*L) isa BandedMatrix
+    @test (L'D'*D*L) isa BandedMatrix
 
     @test bandwidths(materialize(L'D'*D*L)) == (1,1)
 end
 
 @testset "Views" begin
     L = LinearSpline(0:2)
     @test view(L,0.1,1)[1] == L[0.1,1]
+
+    L = LinearSpline(0:2)
+    B1 = L[:,1]
+    @test B1 isa SubQuasiArray{Float64,1}
+    @test size(B1) == (ℵ₁,)
+    @test B1[0.1] == L[0.1,1]
+    @test_throws BoundsError B1[2.2]
+
+    B = L[:,1:2]
+    @test B isa SubQuasiArray{Float64,2}
+    @test B[0.1,:] == L[0.1,1:2]
+
+    B = L[:,2:end-1]
+    @test B[0.1,:] == [0.1]
+end
+
+A = randn(4,4)
+@which lastindex(A,2)
+
+@time begin
+L = LinearSpline(range(0,stop=1,length=10))[:,2:end-1]
+D = Derivative(axes(L,1))
+
+D*L
+
+Derivative(0..1)*parent(L)
+
+
+
+
+M = Mul(D,L)
+A, B = M.factors
+axes(A,2) == axes(B,1) || throw(DimensionMismatch())
+P = parent(B)
+(Derivative(axes(P,1))*P)[parentindices(P)...]
+
+@which axes(D,2)
+
+axes(D,2)
+D*L
+A = -(L'D'D*L)
+f = L*exp.(L.points)
+
+A \ (L'f)
+
+u = A[2:end-1,2:end-1] \ (L'f)[2:end-1]
+
+cond(A[2:end-1,2:end-1])
+
+A[2:end-1,2:end-1] *u - (L'f)[2:end-1]
+
+v = L*[0; u; 0]
+
+
+v[0.2]
+using Plots
+plot(0:0.01:1,getindex.(Ref(v),0:0.01:1))
+    plot!(u1)
+ui
+
+using ApproxFun
+
+x = Fun(0..1)
+    u1 = [Dirichlet(Chebyshev(0..1)); ApproxFun.Derivative()^2] \ [[0,0], exp(x)]
+
+plot(u1)
+
+u_ex = L*u1.(L.points)
+
+xx = 0:0.01:1;
+    plot(xx,getindex.(Ref(D*u_ex),xx))
+
+getindex.(Ref(D*u_ex),xx)
+plot!(u1')
+
+f[0.1]-exp(0.1)
+
+(D*v)'*(D*v)
+
+
+
+L'f
+
+
+
+
+
+(L'f)
+
+(L'f)
+
 end
-L = LinearSpline(0:2)
-L[:,1]