QuasiFill (#63)

* Start QuasiFill

* Fix bugs

* add fill algebra and broadcasting

* more tests

* more tests pass

* tests pass

* Update test_quasifill.jl

* fix stack overflow

* ones, zeros, one, zero, fill

* StaticArray indexing

* Support hcat with quasivec and quasimat

* Update test_quasifill.jl

* union cardinality

* Update quasifill.jl

Author: dlfivefifty

Project.toml

@@ -1,24 +1,29 @@
 name = "QuasiArrays"
 uuid = "c4ea9172-b204-11e9-377d-29865faadc5c"
 authors = ["Sheehan Olver <solver@mac.com>"]
-version = "0.4.10"
+version = "0.5.0"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
 DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
+FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
 ArrayLayouts = "0.5, 0.6"
 DomainSets = "0.4"
+FillArrays = "0.11"
 LazyArrays = "0.20, 0.21"
+StaticArrays = "1"
 julia = "1.5"
 
 [extras]
 Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Base64", "Test", "IntervalSets"]
+test = ["Base64", "IntervalSets", "Random", "Test"]

src/QuasiArrays.jl

@@ -1,5 +1,5 @@
 module QuasiArrays
-using Base, LinearAlgebra, LazyArrays, ArrayLayouts, DomainSets
+using Base, LinearAlgebra, LazyArrays, ArrayLayouts, DomainSets, FillArrays, StaticArrays
 import Base: getindex, size, axes, axes1, length, ==, isequal, iterate, CartesianIndices, LinearIndices,
                 Indices, IndexStyle, getindex, setindex!, parent, vec, convert, similar, copy, copyto!, zero,
                 map, eachindex, eltype, first, last, firstindex, lastindex, in, reshape, permutedims, all,
@@ -8,7 +8,7 @@ import Base: @_inline_meta, DimOrInd, OneTo, @_propagate_inbounds_meta, @_noinli
                 DimsInteger, error_if_canonical_getindex, @propagate_inbounds, _return_type,
                 _maybetail, tail, _getindex, _maybe_reshape, index_ndims, _unsafe_getindex,
                 index_shape, to_shape, unsafe_length, @nloops, @ncall, unalias, _unaliascopy,
-                to_index, to_indices, _to_subscript_indices, _splatmap, dataids, 
+                to_index, to_indices, _to_subscript_indices, _splatmap, dataids,
                 compute_stride1, compute_offset1, fill_to_length
 import Base: Slice, IdentityUnitRange, ScalarIndex, RangeIndex, view, viewindexing, mightalias, ensure_indexable, index_dimsum,
                 check_parent_index_match, reindex, _isdisjoint, unsafe_indices, _unsafe_ind2sub,
@@ -24,19 +24,22 @@ import Base: exp, log, sqrt,
           acos, asin, atan, acsc, asec, acot,
           acosh, asinh, atanh, acsch, asech, acoth
 import Base: Array, Matrix, Vector
-import Base: union, intersect
+import Base: union, intersect, sort, sort!
+import Base: conj, real, imag
+import Base: sum, cumsum, diff
+import Base: ones, zeros, one, zero, fill
 
 import Base.Broadcast: materialize, materialize!, BroadcastStyle, AbstractArrayStyle, Style, broadcasted, Broadcasted, Unknown,
-                        newindex, broadcastable, preprocess, _eachindex, _broadcast_getindex,
-                        DefaultArrayStyle, axistype, throwdm, instantiate, combine_eltypes, eltypes                   
+                        newindex, broadcastable, preprocess, _eachindex, _broadcast_getindex, broadcast_shape,
+                        DefaultArrayStyle, axistype, throwdm, instantiate, combine_eltypes, eltypes
 
 import LinearAlgebra: transpose, adjoint, checkeltype_adjoint, checkeltype_transpose, Diagonal,
                         AbstractTriangular, pinv, inv, promote_leaf_eltypes, power_by_squaring,
-                        integerpow, schurpow, tr, factorize, copy_oftype
+                        integerpow, schurpow, tr, factorize, copy_oftype, rank
 
-import ArrayLayouts: indextype, concretize
+import ArrayLayouts: indextype, concretize, fillzeros, OnesLayout, AbstractFillLayout, FillLayout, ZerosLayout
 import LazyArrays: MemoryLayout, UnknownLayout, Mul, ApplyLayout, BroadcastLayout,
-                    InvOrPInv, ApplyStyle, AbstractLazyLayout, LazyLayout, 
+                    InvOrPInv, ApplyStyle, AbstractLazyLayout, LazyLayout,
                     MulStyle, MulAddStyle, LazyArrayApplyStyle, combine_mul_styles, DefaultArrayApplyStyle,
                     Applied, flatten, _flatten, arguments, _mat_mul_arguments, _vec_mul_arguments, _mul_arguments,
                     rowsupport, colsupport, tuple_type_memorylayouts, applylayout, broadcastlayout,
@@ -50,7 +53,7 @@ export AbstractQuasiArray, AbstractQuasiMatrix, AbstractQuasiVector, materialize
        QuasiArray, QuasiMatrix, QuasiVector, QuasiDiagonal, Inclusion,
        QuasiAdjoint, QuasiTranspose, ApplyQuasiArray, ApplyQuasiMatrix, ApplyQuasiVector,
        BroadcastQuasiArray, BroadcastQuasiMatrix, BroadcastQuasiVector, indextype,
-       QuasiKron, quasikron, UnionVcat, SubQuasiArray
+       QuasiKron, quasikron, UnionVcat, SubQuasiArray, QuasiFill, QuasiZeros, QuasiOnes, QuasiEye
 
 import Base.Broadcast: broadcast_preserving_zero_d
 
@@ -70,6 +73,7 @@ const AbstractQuasiOrArray{T} = Union{AbstractArray{T},AbstractQuasiArray{T}}
 
 
 cardinality(d) = length(d)
+cardinality(d::UnionDomain) = +(cardinality.(d.domains)...)
 
 size(A::AbstractQuasiArray) = map(cardinality, axes(A))
 axes(A::AbstractQuasiArray) = error("Override axes for $(typeof(A))")
@@ -96,6 +100,7 @@ include("matmul.jl")
 include("inv.jl")
 include("quasiadjtrans.jl")
 include("quasidiagonal.jl")
+include("quasifill.jl")
 include("dense.jl")
 
 include("quasikron.jl")

src/abstractquasiarray.jl

@@ -353,8 +353,6 @@ function copymutable(a::AbstractQuasiArray)
     copyto!(similar(a), a)
 end
 
-zero(x::AbstractQuasiArray{T}) where {T} = fill!(similar(x), zero(T))
-
 ## iteration support for arrays by iterating over `eachindex` in the array ##
 # Allows fast iteration by default for both IndexLinear and IndexCartesian arrays
 

src/indices.jl

@@ -111,8 +111,13 @@ to_quasi_index(::Type{IND}, I::AbstractQuasiArray{<:AbstractArray}) where IND<:A
 to_quasi_index(A, IND, i) = to_quasi_index(IND,i)
 
 to_indices(A::AbstractQuasiArray, inds, ::Tuple{}) = ()
+# if inds is empty then we are indexing past indices, so array-like
+to_indices(A::AbstractQuasiArray, ::Tuple{}, I::Tuple{Any,Vararg{Any}}) = 
+    (@_inline_meta; (to_index(A, I[1]), to_indices(A, (), tail(I))...))
 to_indices(A::AbstractQuasiArray, inds, I::Tuple{Any,Vararg{Any}}) =
     (@_inline_meta; (to_quasi_index(A, eltype(inds[1]), I[1]), to_indices(A, _maybetail(inds), tail(I))...))
+@inline to_indices(A::AbstractQuasiArray, inds, I::Tuple{CartesianIndex, Vararg{Any}}) =
+    to_indices(A, inds, (I[1].I..., tail(I)...))
 @inline to_indices(A::AbstractQuasiArray, inds, I::Tuple{Colon, Vararg{Any}}) =
     (uncolon(inds, I), to_indices(A, _maybetail(inds), tail(I))...)
 
@@ -177,6 +182,7 @@ size(S::Inclusion) = (cardinality(S.domain),)
 length(S::Inclusion) = cardinality(S.domain)
 unsafe_length(S::Inclusion) = cardinality(S.domain)
 cardinality(S::Inclusion) = cardinality(S.domain)
+measure(x) = cardinality(x) # TODO: Inclusion(0:0.5:1) should have 
 getindex(S::Inclusion{T}, i::T) where T = (@_inline_meta; @boundscheck checkbounds(S, i); convert(T,i))
 getindex(S::Inclusion{T}, i::AbstractArray{T}) where T = (@_inline_meta; @boundscheck checkbounds(S, i); convert(AbstractArray{T},i))
 getindex(S::Inclusion, i::Inclusion) = (@_inline_meta; @boundscheck checkbounds(S, i); copy(S))
@@ -197,6 +203,7 @@ checkindex(::Type{Bool}, inds::Inclusion{T}, i::T) where T = i ∈ inds
 checkindex(::Type{Bool}, inds::Inclusion, i) = i ⊆ inds
 checkindex(::Type{Bool}, inds::Inclusion, ::Colon) = true
 checkindex(::Type{Bool}, inds::Inclusion, ::Inclusion) = true
+
 function __checkindex(::Type{Bool}, inds::Inclusion, I::AbstractArray)
     @_inline_meta
     b = true
@@ -213,9 +220,6 @@ checkindex(::Type{Bool}, inds::Inclusion{T}, I::AbstractArray{T}) where T<:Abstr
 checkindex(::Type{Bool}, inds::Inclusion{T}, I::AbstractArray{<:AbstractArray}) where T<:AbstractArray = 
     __checkindex(Bool, inds, convert(AbstractArray{T}, I))
 
-@propagate_inbounds _affine_checkindex(inds, r) = isempty(r) | (checkindex(Bool, inds, Base.to_indices(inds, (first(r),))...) & checkindex(Bool, inds,  Base.to_indices(inds, (last(r),))...))
-@propagate_inbounds checkindex(::Type{Bool}, inds::Inclusion, r::AbstractRange) = _affine_checkindex(inds, r)
-
 
 checkindex(::Type{Bool}, indx::Inclusion, I::AbstractVector{Bool}) = indx == axes1(I)
 checkindex(::Type{Bool}, indx::Inclusion, I::AbstractArray{Bool}) = false

src/lazyquasiarrays.jl

@@ -79,6 +79,7 @@ MemoryLayout(M::Type{ApplyQuasiArray{T,N,F,Args}}) where {T,N,F,Args} =
 copy(A::Applied{LazyQuasiArrayApplyStyle}) = ApplyQuasiArray(A)
 copy(A::Applied{<:AbstractQuasiArrayApplyStyle}) = QuasiArray(A)
 QuasiArray(A::Applied) = QuasiArray(ApplyQuasiArray(A))
+QuasiArray(A::Applied{<:Any,typeof(inv)}) = inv(QuasiArray(A.args[1]))
 
 
 ####

src/matmul.jl

@@ -89,9 +89,10 @@ end
 
 struct QuasiArrayLayout <: MemoryLayout end
 MemoryLayout(::Type{<:AbstractQuasiArray}) = QuasiArrayLayout()
-
+Array(M::Mul{<:Any,<:Any,<:Any,<:AbstractQuasiMatrix}) = eltype(M)[M[k,j] for k in axes(M)[1], j in axes(M)[2]]
+Array(M::Mul{<:Any,<:Any,<:Any,<:AbstractQuasiVector}) = eltype(M)[M[k] for k in axes(M)[1]]
 _quasi_mul(M, _) = QuasiArray(M)
-_quasi_mul(M, ::NTuple{N,OneTo{Int}}) where N = Array(QuasiArray(M))
+_quasi_mul(M, ::NTuple{N,OneTo{Int}}) where N = Array(M)
 copy(M::Mul{QuasiArrayLayout,QuasiArrayLayout}) = _quasi_mul(M, axes(M))
 copy(M::Mul{QuasiArrayLayout}) = _quasi_mul(M, axes(M))
 copy(M::Mul{<:Any,QuasiArrayLayout}) = _quasi_mul(M, axes(M))
@@ -100,6 +101,9 @@ copy(M::Mul{ApplyLayout{typeof(\)},QuasiArrayLayout}) = ApplyQuasiArray(M)
 copy(M::Mul{QuasiArrayLayout,<:AbstractLazyLayout}) = ApplyQuasiArray(M)
 
 
+LazyArrays._vec_mul_view(a::AbstractQuasiVector, kr, ::Colon) = view(a, kr)
+
+
 ###
 # Scalar special case, simplifies x * A and A * x
 # Find an AbstractArray to latch on to by commuting

src/multidimensional.jl

@@ -9,6 +9,7 @@ module QuasiIteratorsMD
     using Base: IndexLinear, IndexCartesian, AbstractCartesianIndex, fill_to_length, tail
     using Base.Iterators: Reverse
     import QuasiArrays: AbstractQuasiArray, AbstractQuasiVector, domain, AbstractQuasiOrVector
+    import StaticArrays: StaticArray
 
     export QuasiCartesianIndex, QuasiCartesianIndices
 
@@ -423,6 +424,10 @@ to_indices(A::AbstractQuasiArray, I::Tuple{Any}) = (@_inline_meta; to_indices(A,
     _, indstail = IteratorsMD.split(inds, Val(N))
     (to_index(A, I[1]), to_indices(A, indstail, tail(I))...)
 end
+Base.to_indices(A::AbstractQuasiArray, I::Tuple{Vararg{Union{Integer, CartesianIndex, StaticArray{<:Tuple,Int}}}}) =
+    to_indices(A, axes(A), I)
+
+@inline index_shape(A::AbstractQuasiArray, rest...) = (axes(A)..., index_shape(rest...)...)
 
 ### From abstractarray.jl: Internal multidimensional indexing definitions ###
 

src/quasiarray.jl

@@ -8,6 +8,9 @@ struct QuasiArray{T,N,AXES<:NTuple{N,AbstractVector}} <: AbstractQuasiArray{T,N}
     end
 end
 
+QuasiArray{T,N,AXES}(par::AbstractArray{<:Any,N}, axes::AXES) where {T,N,AXES<:NTuple{N,AbstractVector}} =
+    QuasiArray{T,N,AXES}(convert(AbstractArray{T,N}, par)::AbstractArray{T,N}, axes)
+
 const QuasiMatrix{T,AXES<:NTuple{2,AbstractVector}} = QuasiArray{T,2,AXES}
 const QuasiVector{T,AXES<:Tuple{AbstractVector}} = QuasiArray{T,1,AXES}
 
@@ -28,6 +31,39 @@ QuasiMatrix(par::AbstractMatrix{T}, axes::NTuple{2,AbstractVector}) where T =
 QuasiVector(par::AbstractVector{T}, axes::Tuple{AbstractVector}) where T = 
     QuasiArray{T,1,typeof(axes)}(par, axes)
 
+QuasiArray{T}(par::AbstractArray{<:Any,N}, axes::NTuple{N,AbstractVector}) where {T,N} = 
+    QuasiArray{T,N,typeof(axes)}(par, axes)
+QuasiMatrix{T}(par::AbstractMatrix, axes::NTuple{2,AbstractVector}) where T = 
+    QuasiMatrix{T,typeof(axes)}(par, axes)
+QuasiVector{T}(par::AbstractVector, axes::Tuple{AbstractVector}) where T = 
+    QuasiVector{T,typeof(axes)}(par, axes)
+
+QuasiArray(par::AbstractArray{T,N}, axes::Vararg{AbstractVector,N}) where {T,N} = 
+    QuasiArray(par, axes)
+QuasiMatrix(par::AbstractMatrix{T}, axes::Vararg{AbstractVector,2}) where T = 
+    QuasiMatrix(par, axes)
+QuasiVector(par::AbstractVector{T}, axes::AbstractVector) where T = 
+    QuasiVector(par, (axes,))
+
+QuasiArray{T}(par::AbstractArray{<:Any,N}, axes::Vararg{AbstractVector,N}) where {T,N} = 
+    QuasiArray{T}(par, axes)
+QuasiMatrix{T}(par::AbstractMatrix, axes::Vararg{AbstractVector,2}) where T = 
+    QuasiMatrix{T}(par, axes)
+QuasiVector{T}(par::AbstractVector, axes::AbstractVector) where T = 
+    QuasiVector{T}(par, (axes,))
+
+
+
+QuasiMatrix(λ::UniformScaling, ax::NTuple{2,AbstractVector}) = QuasiMatrix(Matrix(λ,map(length,ax)...), ax)
+QuasiMatrix{T}(λ::UniformScaling, ax::NTuple{2,AbstractVector}) where T = QuasiMatrix{T}(Matrix(λ,map(length,ax)...), ax)
+QuasiMatrix(λ::UniformScaling, ax1::AbstractVector, ax2::AbstractVector) = QuasiMatrix(λ, (ax1, ax2))
+QuasiMatrix{T}(λ::UniformScaling, ax1::AbstractVector, ax2::AbstractVector) where T = QuasiMatrix{T}(λ, (ax1, ax2))
+
+
+QuasiArray(par::AbstractArray) = QuasiArray(par, axes(par))
+QuasiMatrix(par::AbstractArray) = QuasiMatrix(par, axes(par))
+QuasiVector(par::AbstractArray) = QuasiVector(par, axes(par))
+
 QuasiArray{T,N,AXES}(par::AbstractArray{T,N}, axes::NTuple{N,AbstractQuasiOrVector}) where {T,N,AXES} =
     QuasiArray{T,N,AXES}(par, map(domain,axes))
 QuasiArray{T,N}(par::AbstractArray{T,N}, axes::NTuple{N,AbstractQuasiOrVector}) where {T,N} =
@@ -45,9 +81,9 @@ QuasiVector(par::AbstractVector{T}, axes::AbstractArray) where {T} =
     QuasiVector(par, (axes,))
 
 
-QuasiArray(a::AbstractQuasiArray) = QuasiArray(Array(a), axes(a))
-QuasiArray{T}(a::AbstractQuasiArray) where T = QuasiArray(Array{T}(a), axes(a))
-QuasiArray{T,N}(a::AbstractQuasiArray{<:Any,N}) where {T,N} = QuasiArray(Array{T}(a), axes(a))
+QuasiArray(a::AbstractQuasiArray) = QuasiArray(a[map(collect,axes(a))...], axes(a))
+QuasiArray{T}(a::AbstractQuasiArray) where T = QuasiArray(convert(AbstractArray{T},a[map(collect,axes(a))...]), axes(a))
+QuasiArray{T,N}(a::AbstractQuasiArray{<:Any,N}) where {T,N} = QuasiArray(convert(AbstractArray{T,N},a[map(collect,axes(a))...]), axes(a))
 QuasiArray{T,N,AXES}(a::AbstractQuasiArray{<:Any,N}) where {T,N,AXES} = QuasiArray{T,N,AXES}(Array{T}(a), axes(a))
 QuasiMatrix(a::AbstractQuasiMatrix) = QuasiArray(a)
 QuasiVector(a::AbstractQuasiVector) = QuasiArray(a)
@@ -61,14 +97,17 @@ convert(::Type{AbstractQuasiArray{T}}, a::QuasiArray) where T = convert(QuasiArr
 convert(::Type{AbstractQuasiArray{T,N}}, a::QuasiArray{T,N}) where {T,N} = a
 convert(::Type{AbstractQuasiArray{T}}, a::QuasiArray{T}) where T = a
 
+AbstractQuasiArray{T,N}(a::QuasiArray) where {T,N} = QuasiArray{T,N}(a)
+AbstractQuasiArray{T}(a::QuasiArray) where T = QuasiArray{T}(a)
+
 
 _inclusion(d::Slice) = d
 _inclusion(d::OneTo) = d
 _inclusion(d::IdentityUnitRange{<:Integer}) = d
-_inclusion(d::AbstractUnitRange{<:Integer}) = Slice(d)
+_inclusion(d::AbstractUnitRange{<:Integer}) = IdentityUnitRange(d)
 _inclusion(d) = Inclusion(d)
 
-axes(A::QuasiArray) = _inclusion.(A.axes)
+axes(A::QuasiArray) = map(_inclusion,A.axes)
 parent(A::QuasiArray) = A.parent
 
 @propagate_inbounds @inline function _getindex(::Type{IND}, A::QuasiArray, I::IND) where IND
@@ -89,3 +128,9 @@ end
 
 ^(A::QuasiMatrix, p::Number) = _quasimatrix_pow(A, p)
 ^(A::QuasiMatrix, p::Integer) = _quasimatrix_pow(A, p)
+
+
+==(A::QuasiArray{T,N,NTuple{N,OneTo{Int}}}, B::AbstractArray{V,N}) where {T,V,N} =
+    axes(A) == axes(B) && A.parent == B
+==(B::AbstractArray{V,N}, A::QuasiArray{T,N,NTuple{N,OneTo{Int}}}) where {T,V,N} =
+    A == B

src/quasiconcat.jl

@@ -7,16 +7,24 @@ ApplyStyle(::typeof(hcat), ::Type{<:AbstractQuasiArray}...) = LazyQuasiArrayAppl
 
 axes(f::ApplyQuasiMatrix{<:Any,typeof(hcat)}) = (axes(f.args[1],1), oneto(sum(size.(f.args,2))))
 
+_vecormat_getindex(A::AbstractQuasiMatrix,k,ξ) = A[k,ξ]
+_vecormat_getindex(A::AbstractQuasiVector,k,ξ) = A[k]
+@inline LazyArrays._view_hcat(a::AbstractQuasiVector, kr, jr::Colon) = view(a, kr)
+# TODO: generalize
+@inline LazyArrays._view_hcat(a::AbstractQuasiVector, kr::Number, jr) = a[kr]
+
 function getindex(f::ApplyQuasiMatrix{T,typeof(hcat)}, k::Number, j::Number) where T
     ξ = j
     for A in f.args
         n = size(A,2)
-        ξ ≤ n && return T(A[k,ξ])::T
+        ξ ≤ n && return T(_vecormat_getindex(A,k,ξ))::T
         ξ -= n
     end
     throw(BoundsError(f, (k,j)))
 end
 
+\(A::AbstractQuasiArray, H::ApplyQuasiMatrix{<:Any,typeof(hcat)}) = hcat((Ref(A) .\ H.args)...)
+
 """
     UnionVcat
 

src/quasidiagonal.jl

@@ -35,7 +35,7 @@ julia> QuasiDiagonal(A)
 """
 QuasiDiagonal(A::AbstractQuasiMatrix) = QuasiDiagonal(diag(A))
 
-Diagonal(A::AbstractQuasiArray) = QuasiDiagonal(A)
+@deprecate Diagonal(A::AbstractQuasiArray) QuasiDiagonal(A)
 
 
 """
@@ -63,8 +63,11 @@ QuasiDiagonal{T}(D::QuasiDiagonal{T}) where {T} = D
 QuasiDiagonal{T}(D::QuasiDiagonal) where {T} = QuasiDiagonal{T}(D.diag)
 
 AbstractQuasiMatrix{T}(D::QuasiDiagonal) where {T} = QuasiDiagonal{T}(D)
-Matrix(D::QuasiDiagonal) = diagm(0 => D.diag)
-Array(D::QuasiDiagonal) = Matrix(D)
+AbstractQuasiArray{T}(D::QuasiDiagonal) where T = AbstractQuasiMatrix{T}(D)
+
+convert(::Type{QuasiDiagonal}, A::QuasiDiagonal) = A
+convert(::Type{QuasiDiagonal{T}}, A::QuasiDiagonal{T}) where T = A
+convert(::Type{QuasiDiagonal{T}}, A::QuasiDiagonal) where T = QuasiDiagonal(convert(AbstractQuasiVector{T}, A.diag))
 
 # For D<:QuasiDiagonal, similar(D[, neweltype]) should yield a QuasiDiagonal matrix.
 # On the other hand, similar(D, [neweltype,] shape...) should yield a sparse matrix.
@@ -196,3 +199,7 @@ for f in (:exp, :log, :sqrt,
     @eval $f(D::QuasiDiagonal) = QuasiDiagonal($f.(D.diag))
 end
 
+
+permutedims(D::QuasiDiagonal) = D
+inv(D::QuasiDiagonal) = QuasiDiagonal(inv.(D.diag))
+copy(D::QuasiDiagonal) = QuasiDiagonal(copy(D.diag))