Generalized @simplify to arbitrary amount of arguments

[jagot]

With this, I can also write the simplified version of the weak Laplacian:

@simplify function *(Ac::AdjointFEDVROrRestricted, Dc::QuasiAdjoint{<:Any,<:Derivative},
                     D::Derivative, B::FEDVROrRestricted)
    # Do something
end

I noticed however, in restricted bases, i.e. with SubQuasiArray{T,N,<:Basis}, I get a dense matrix back, instead of BlockSkylineMatrix, since the simplify framework calculate the derivative in the unrestricted basis and then indexes it afterwards. I think there are two ways to fix this:

1. Pass the restricted bases to Base.copy generated by @simplify and let the user worry about handling the restrictions correctly (which is what I've done before; hence the FEDVROrRestricted in the signature above)
2. Make sure slicing of a BlockSkylineMatrix returns a BlockSkylineMatrix and not a dense one. I guess that is only well-defined if one is "shaving" off from the outside, i.e. no non-contiguous slicing.

To illustrate what's going on:

julia> B = FEDVR(1.0:3, 4)
FEDVR{Float64} basis with 2 elements on 1.0..3.0

julia> B̃ = B[:,2:end-1]
FEDVR{Float64} basis with 2 elements on 1.0..3.0, restricted to basis functions 2..6 ⊂ 1..7

julia> D = Derivative(axes(B,1))
Derivative{Float64,IntervalSets.Interval{:closed,:closed,Float64}}(Inclusion(1.0..3.0))

julia> B'D'D*B
3×3-blocked 7×7 BlockBandedMatrices.BlockSkylineMatrix{Float64,Array{Float64,1},BlockBandedMatrices.BlockSkylineSizes{Tuple{BlockArrays.BlockedUnitRange{Array{Int64,1}},BlockArrays.BlockedUnitRange{Array{Int64,1}}},Array{Int64,1},Array{Int64,1},BandedMatrices.BandedMatrix{Int64,Array{Int64,2},Base.OneTo{Int64}},Array{Int64,1}}}:
 -52.0       26.1803    -3.81966  │    1.41421  │     ⋅          ⋅          ⋅
  26.1803   -20.0       10.0      │   -2.70091  │     ⋅          ⋅          ⋅
  -3.81966   10.0      -20.0      │   18.5123   │     ⋅          ⋅          ⋅
 ─────────────────────────────────┼─────────────┼─────────────────────────────────
   1.41421   -2.70091   18.5123   │  -52.0      │   18.5123    -2.70091    1.41421
 ─────────────────────────────────┼─────────────┼─────────────────────────────────
    ⋅          ⋅          ⋅       │   18.5123   │  -20.0       10.0       -3.81966
    ⋅          ⋅          ⋅       │   -2.70091  │   10.0      -20.0       26.1803
    ⋅          ⋅          ⋅       │    1.41421  │   -3.81966   26.1803   -52.0

julia> B̃'D'D*B̃
5×5 Array{Float64,2}:
 -20.0       10.0      -2.70091    0.0       0.0
  10.0      -20.0      18.5123     0.0       0.0
  -2.70091   18.5123  -52.0       18.5123   -2.70091
   0.0        0.0      18.5123   -20.0      10.0
   0.0        0.0      -2.70091   10.0     -20.0

[codecov]

Codecov Report

Merging #37 into master will decrease coverage by 1.06%.
The diff coverage is n/a.

@@            Coverage Diff             @@
##           master      #37      +/-   ##
==========================================
- Coverage   75.21%   74.14%   -1.07%     
==========================================
  Files           4        4              
  Lines         347      352       +5     
==========================================
  Hits          261      261              
- Misses         86       91       +5     

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 75d7dd9...e429d8d. Read the comment docs.

[jagot]

Actually, in light of #36, I think I'd prefer an alternative interface where you'd deal with the restrictions explicitly; if you're anyway only going to work on a subinterval, it's a waste of effort to calculate the full materialization and then crop the matrix.

[dlfivefifty]

I feel like the @simplify macro is a hack waiting for a better design....

[jagot]

Sure, I was thinking about writing a @materialize macro that would give the user some more control, but I did not get around to it yet.

Happy New Year!

[jagot]

I came up with this small @materialize DSL (which could go into LazyArrays.jl, since it does not depend on anything Quasi):

# For unparametrized destination types
generate_copyto!_signature(dest, dest_type::Symbol, Msig) =
    :(Base.copyto!($(dest)::$(dest_type), applied_obj::$(Msig)))

# For parametrized destination types
function generate_copyto!_signature(dest, dest_type::Expr, Msig)
    dest_type.head == :curly ||
        throw(ArgumentError("Invalid destination specification $(dest)::$(dest_type)"))
    :(Base.copyto!($(dest)::$(dest_type), applied_obj::$(Msig)) where {$(dest_type.args[2:end]...)})
end

function generate_copyto!(body, factor_names, Msig)
    body.head == :(->) ||
        throw(ArgumentError("Invalid copyto! specification"))
    body.args[1].head == :(::) ||
        throw(ArgumentError("Invalid destination specification $(body.args[1])"))
    (dest,dest_type) = body.args[1].args
    copyto!_signature = generate_copyto!_signature(dest, dest_type, Msig)
    f_body = quote
        axes($dest) == axes(applied_obj) || throw(DimensionMismatch("axes must be same"))
        $(factor_names) = applied_obj.args
        $(body.args[2].args...)
        $(dest)
    end
    Expr(:function, copyto!_signature, f_body)
end

# Generate three functions: similar, copyto!, and materialize for a
# LazyArrays.Applied object.
macro materialize(expr)
    expr.head == :function || expr.head == :(=) || error("Must start with a function")
    @assert expr.args[1].head == :call
    op = expr.args[1].args[1]
    factor_types = :(<:Tuple{})
    factor_names = :(())

    bodies = filter(e -> !(e isa LineNumberNode), expr.args[2].args)
    length(bodies) < 2 && throw(ArgumentError("At least two blocks required (one similar and at least one copyto!)"))

    # Generate Applied signature
    for arg in expr.args[1].args[2:end]
        arg isa Expr && arg.head == :(::) ||
            throw(ArgumentError("Invalid argument specification $(arg)"))
        arg_name, arg_typ = arg.args
        push!(factor_types.args[1].args, :(<:$(arg_typ)))
        push!(factor_names.args, arg_name)
    end
    Msig = :(LazyArrays.Applied{<:Any, typeof($op), $(factor_types)})

    sim_body = first(bodies)
    sim_body.head == :(->) ||
        throw(ArgumentError("Invalid similar specification"))
    T = first(sim_body.args)
    
    copytos! = map(body -> generate_copyto!(body, factor_names, Msig), bodies[2:end])

    f = quote
        function Base.similar(applied_obj::$Msig, ::Type{$T}=eltype(applied_obj)) where $T
            $(factor_names) = applied_obj.args
            $(first(bodies).args[2])
        end
        
        $(copytos!...)

        materialize(applied_obj::$Msig) = copyto!(similar(applied_obj, eltype(applied_obj)), applied_obj)
    end
    esc(f)
end

Perhaps you find a bit too "magical", but with it, I can do

# These type-aliases I would like to see in ContinuumArrays.jl
const RestrictionMatrix = BandedMatrix{<:Int, <:FillArrays.Ones}
const RestrictedQuasiArray{T,N,B<:Basis} = SubQuasiArray{T,N,B}
const AdjointRestrictedQuasiArray{T,N,B<:Basis} = QuasiAdjoint{T,<:RestrictedQuasiArray{T,N,B}}

const BasisOrRestricted{B<:Basis} = Union{B,<:RestrictedQuasiArray{<:Any,<:Any,<:B}}
const AdjointBasisOrRestricted{B<:Basis} = Union{<:QuasiAdjoint{<:Any,B},<:AdjointRestrictedQuasiArray{<:Any,<:Any,<:B}}

indices(B::AbstractFiniteDifferences) = axes(B)
indices(B::RestrictedQuasiArray{<:Any,2,<:AbstractFiniteDifferences}) = B.indices

@materialize function *(Ac::AdjointBasisOrRestricted{<:AbstractFiniteDifferences},
                        D::QuasiDiagonal,
                        B::BasisOrRestricted{<:AbstractFiniteDifferences})
    T -> begin # generates similar
        A = parent(Ac)
        parent(A) == parent(B) ||
            throw(ArgumentError("Cannot multiply functions on different grids"))
        Ai,Bi = last(indices(A)), last(indices(B))
        if Ai == Bi # "diagonal" restriction
            Diagonal(Vector{T}(undef, length(Ai)))
        else
            m,n = length(Ai),length(Bi)
            offset = Ai[1]-Bi[1]
            BandedMatrix{T}(undef, (m,n), (-offset,offset))
        end
    end
    dest::Diagonal{T} -> begin # generate copyto!(dest::Diagonal{T}, ...) where T
        dest.diag .= 1
    end
    dest::BandedMatrix{T} -> begin # generate copyto!(dest::BandedMatrix{T}, ...) where T
        dest.data .= 3
    end
end

Using this, we get

N = 10
ρ = 0.25

R = FiniteDifferences(N, ρ)
R̃ = R[:, 3:8]
r = axes(R, 1)
x = QuasiDiagonal(r)

julia> materialize(applied(*, R', x, R))
10×10 Diagonal{Float64,Array{Float64,1}}:
 1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0   ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0

julia> materialize(applied(*, R̃', x, R̃))
6×6 Diagonal{Float64,Array{Float64,1}}:
 1.0   ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅   1.0   ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅   1.0   ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅   1.0   ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅   1.0   ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅   1.0

julia> materialize(applied(*, R̃', x, R))
6×10 BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}:
  ⋅    ⋅   3.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅   3.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅   3.0   ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅   3.0   ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   3.0   ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   3.0   ⋅    ⋅

julia> materialize(applied(*, R', x, R̃))
10×6 BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}:
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
 3.0   ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅   3.0   ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅   3.0   ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅   3.0   ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅   3.0   ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅   3.0
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅

For some reason, it does not play well with infix or apply:

julia> R̃'x*R
ERROR: MethodError: no method matching Int64(::ContinuumArrays.AlephInfinity{1})
Closest candidates are:
  Int64(::T) where T<:Number at boot.jl:718
  Int64(::Union{Bool, Int32, Int64, UInt32, UInt64, UInt8, Int128, Int16, Int8, UInt128, UInt16}) at boot.jl:710
  Int64(::Ptr) at boot.jl:720
  ...
Stacktrace:
 [1] convert(::Type{Int64}, ::ContinuumArrays.AlephInfinity{1}) at ./number.jl:7
 [2] convert(::Type{Tuple{Vararg{Int64,N} where N}}, ::Tuple{ContinuumArrays.AlephInfinity{1}}) at ./essentials.jl:312 (repeats 2 times)
 [3] Array{Float64,N} where N(::UndefInitializer, ::Tuple{Int64,ContinuumArrays.AlephInfinity{1}}) at ./baseext.jl:24
 [4] similar(::Applied{LazyArrays.DefaultArrayApplyStyle,typeof(*),Tuple{QuasiAdjoint{Float64,QuasiArrays.SubQuasiArray{Float64,2,FiniteDifferences{Float64,Int64},Tuple{Inclusion{Float64,Interval{:closed,:closed,Float64}},UnitRange{Int64}},false}},QuasiDiagonal{Float64,Inclusion{Float64,Interval{:closed,:closed,Float64}}}}}, ::Type{Float64}, ::Tuple{Base.OneTo{Int64},Inclusion{Float64,Interval{:closed,:closed,Float64}}}) at /Users/jagot/.julia/packages/LazyArrays/14GOk/src/lazyapplying.jl:65
 [5] similar(::Applied{LazyArrays.DefaultArrayApplyStyle,typeof(*),Tuple{QuasiAdjoint{Float64,QuasiArrays.SubQuasiArray{Float64,2,FiniteDifferences{Float64,Int64},Tuple{Inclusion{Float64,Interval{:closed,:closed,Float64}},UnitRange{Int64}},false}},QuasiDiagonal{Float64,Inclusion{Float64,Interval{:closed,:closed,Float64}}}}}, ::Type{Float64}) at /Users/jagot/.julia/packages/LazyArrays/14GOk/src/lazyapplying.jl:66
 [6] similar(::Applied{LazyArrays.DefaultArrayApplyStyle,typeof(*),Tuple{QuasiAdjoint{Float64,QuasiArrays.SubQuasiArray{Float64,2,FiniteDifferences{Float64,Int64},Tuple{Inclusion{Float64,Interval{:closed,:closed,Float64}},UnitRange{Int64}},false}},QuasiDiagonal{Float64,Inclusion{Float64,Interval{:closed,:closed,Float64}}}}}) at /Users/jagot/.julia/packages/LazyArrays/14GOk/src/lazyapplying.jl:67
 [7] copy(::Applied{LazyArrays.DefaultArrayApplyStyle,typeof(*),Tuple{QuasiAdjoint{Float64,QuasiArrays.SubQuasiArray{Float64,2,FiniteDifferences{Float64,Int64},Tuple{Inclusion{Float64,Interval{:closed,:closed,Float64}},UnitRange{Int64}},false}},QuasiDiagonal{Float64,Inclusion{Float64,Interval{:closed,:closed,Float64}}}}}) at /Users/jagot/.julia/packages/LazyArrays/14GOk/src/linalg/mul.jl:101
 [8] materialize at /Users/jagot/.julia/packages/LazyArrays/14GOk/src/lazyapplying.jl:48 [inlined]
 [9] apply at /Users/jagot/.julia/packages/LazyArrays/14GOk/src/lazyapplying.jl:45 [inlined]
 [10] *(::QuasiAdjoint{Float64,QuasiArrays.SubQuasiArray{Float64,2,FiniteDifferences{Float64,Int64},Tuple{Inclusion{Float64,Interval{:closed,:closed,Float64}},UnitRange{Int64}},false}}, ::QuasiDiagonal{Float64,Inclusion{Float64,Interval{:closed,:closed,Float64}}}) at /Users/jagot/.julia/packages/QuasiArrays/d7Sua/src/matmul.jl:58
 [11] top-level scope at REPL[11]:1

julia> apply(*, R̃', x, R)
ApplyQuasiArray{Float64,2,typeof(*),Tuple{QuasiAdjoint{Float64,QuasiArrays.SubQuasiArray{Float64,2,FiniteDifferences{Float64,Int64},Tuple{Inclusion{Float64,Interval{:closed,:closed,Float64}},UnitRange{Int64}},false}},QuasiDiagonal{Float64,Inclusion{Float64,Interval{:closed,:closed,Float64}}},FiniteDifferences{Float64,Int64}}}(*, (QuasiAdjoint{Float64,QuasiArrays.SubQuasiArray{Float64,2,FiniteDifferences{Float64,Int64},Tuple{Inclusion{Float64,Interval{:closed,:closed,Float64}},UnitRange{Int64}},false}}(Finite differences basis {Float64} on 0.0..2.75 with 10 points spaced by Δx = 0.25, restricted to basis functions 3..8 ⊂ 1..10), QuasiDiagonal{Float64,Inclusion{Float64,Interval{:closed,:closed,Float64}}}(Inclusion(0.0..2.75)), Finite differences basis {Float64} on 0.0..2.75 with 10 points spaced by Δx = 0.25))

[dlfivefifty]

I think this is out of date

[jagot]

Probably, yes.