From 7413876d899aa5d967a14919f9cf5a85e4cebccf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20Hakkel?= Date: Thu, 11 Apr 2024 18:16:29 +0200 Subject: [PATCH] Storage type (#23) Co-authored-by: Tamas Hakkel --- Project.toml | 2 +- src/calculus/AdjointOperator.jl | 4 +- src/calculus/AffineAdd.jl | 18 +++--- src/calculus/Ax_mul_Bx.jl | 19 +++--- src/calculus/Ax_mul_Bxt.jl | 21 +++---- src/calculus/Axt_mul_Bx.jl | 21 +++---- src/calculus/BroadCast.jl | 19 +++--- src/calculus/Compose.jl | 17 +++--- src/calculus/DCAT.jl | 62 +++++++++---------- src/calculus/HCAT.jl | 86 +++++++++++++-------------- src/calculus/HadamardProd.jl | 17 +++--- src/calculus/Jacobian.jl | 38 ++++++------ src/calculus/Reshape.jl | 10 ++-- src/calculus/Scale.jl | 20 +++---- src/calculus/Sum.jl | 20 +++---- src/calculus/VCAT.jl | 84 +++++++++++++------------- src/linearoperators/Conv.jl | 36 +++++++---- src/linearoperators/DCT.jl | 12 ++-- src/linearoperators/DiagOp.jl | 4 +- src/linearoperators/Eye.jl | 2 +- src/linearoperators/Filt.jl | 4 +- src/linearoperators/FiniteDiff.jl | 62 ++++++++----------- src/linearoperators/GetIndex.jl | 14 ++--- src/linearoperators/IRDFT.jl | 10 ++-- src/linearoperators/LBFGS.jl | 14 +++-- src/linearoperators/LMatrixOp.jl | 18 +++--- src/linearoperators/MIMOFilt.jl | 10 ++-- src/linearoperators/MatrixOp.jl | 10 ++-- src/linearoperators/MyLinOp.jl | 4 +- src/linearoperators/RDFT.jl | 10 ++-- src/linearoperators/Variation.jl | 19 +++--- src/linearoperators/Xcorr.jl | 6 +- src/linearoperators/ZeroPad.jl | 24 ++++---- src/linearoperators/Zeros.jl | 12 ++-- src/nonlinearoperators/Atan.jl | 6 +- src/nonlinearoperators/Cos.jl | 6 +- src/nonlinearoperators/Exp.jl | 6 +- src/nonlinearoperators/Pow.jl | 6 +- src/nonlinearoperators/Sech.jl | 6 +- src/nonlinearoperators/Sigmoid.jl | 8 +-- src/nonlinearoperators/Sin.jl | 6 +- src/nonlinearoperators/SoftMax.jl | 14 +++-- src/nonlinearoperators/SoftPlus.jl | 6 +- src/nonlinearoperators/Tanh.jl | 6 +- src/properties.jl | 95 +++++++++++++++++++++++------- src/syntax.jl | 11 +--- 46 files changed, 468 insertions(+), 437 deletions(-) diff --git a/Project.toml b/Project.toml index 4bfb91f..4a86f9a 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "AbstractOperators" uuid = "d9c5613a-d543-52d8-9afd-8f241a8c3f1c" -version = "0.3" +version = "0.4" [deps] AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c" diff --git a/src/calculus/AdjointOperator.jl b/src/calculus/AdjointOperator.jl index 2c53237..c507e87 100644 --- a/src/calculus/AdjointOperator.jl +++ b/src/calculus/AdjointOperator.jl @@ -3,7 +3,7 @@ export AdjointOperator """ `AdjointOperator(A::AbstractOperator)` -Shorthand constructor: +Shorthand constructor: `'(A::AbstractOperator)` @@ -19,7 +19,7 @@ julia> [DFT(10); DCT(10)]' """ struct AdjointOperator{T <: AbstractOperator} <: AbstractOperator A::T - function AdjointOperator(A::T) where {T<:AbstractOperator} + function AdjointOperator(A::T) where {T<:AbstractOperator} is_linear(A) == false && error("Cannot transpose a nonlinear operator. You might use `jacobian`") new{T}(A) end diff --git a/src/calculus/AffineAdd.jl b/src/calculus/AffineAdd.jl index 8428559..b767858 100644 --- a/src/calculus/AffineAdd.jl +++ b/src/calculus/AffineAdd.jl @@ -3,7 +3,7 @@ export AffineAdd """ `AffineAdd(A::AbstractOperator, d, [sign = true])` -Affine addition to `AbstractOperator` with an array or scalar `d`. +Affine addition to `AbstractOperator` with an array or scalar `d`. Use `sign = false` to perform subtraction. @@ -26,17 +26,17 @@ true struct AffineAdd{L <: AbstractOperator, D <: Union{AbstractArray, Number}, S} <: AbstractOperator A::L d::D - function AffineAdd(A::L, d::D, sign::Bool = true) where {L, D <: AbstractArray} - if size(d) != size(A,1) + function AffineAdd(A::L, d::D, sign::Bool = true) where {L, D <: AbstractArray} + if size(d) != size(A,1) throw(DimensionMismatch("codomain size of $A not compatible with array `d` of size $(size(d))")) end - if eltype(d) != codomainType(A) + if eltype(d) != codomainType(A) error("cannot tilt opertor having codomain type $(codomainType(A)) with array of type $(eltype(d))") end new{L,D,sign}(A,d) end # scalar - function AffineAdd(A::L, d::D, sign::Bool = true) where {L, D <: Number} + function AffineAdd(A::L, d::D, sign::Bool = true) where {L, D <: Number} if typeof(d) <: Complex && codomainType(A) <: Real error("cannot tilt opertor having codomain type $(codomainType(A)) with array of type $(eltype(d))") end @@ -46,12 +46,12 @@ end # Mappings # array -function mul!(y::DD, T::AffineAdd{L, D, true}, x) where {L <: AbstractOperator, DD, D} +function mul!(y::DD, T::AffineAdd{L, D, true}, x) where {L <: AbstractOperator, DD, D} mul!(y,T.A,x) y .+= T.d end -function mul!(y::DD, T::AffineAdd{L, D, false}, x) where {L <: AbstractOperator, DD, D} +function mul!(y::DD, T::AffineAdd{L, D, false}, x) where {L <: AbstractOperator, DD, D} mul!(y,T.A,x) y .-= T.d end @@ -70,7 +70,7 @@ is_null(L::AffineAdd) = is_null(L.A) is_eye(L::AffineAdd) = is_diagonal(L.A) is_diagonal(L::AffineAdd) = is_diagonal(L.A) is_invertible(L::AffineAdd) = is_invertible(L.A) -is_AcA_diagonal(L::AffineAdd) = is_AcA_diagonal(L.A) +is_AcA_diagonal(L::AffineAdd) = is_AcA_diagonal(L.A) is_AAc_diagonal(L::AffineAdd) = is_AAc_diagonal(L.A) is_full_row_rank(L::AffineAdd) = is_full_row_rank(L.A) is_full_column_rank(L::AffineAdd) = is_full_column_rank(L.A) @@ -90,7 +90,7 @@ sign(T::AffineAdd{L,D, true}) where {L,D} = 1 function permute(T::AffineAdd{L,D,S}, p::AbstractVector{Int}) where {L,D,S} A = permute(T.A,p) - return AffineAdd(A,T.d,S) + return AffineAdd(A,T.d,S) end displacement(A::AffineAdd{L,D,true}) where {L,D} = A.d .+ displacement(A.A) diff --git a/src/calculus/Ax_mul_Bx.jl b/src/calculus/Ax_mul_Bx.jl index 46832f9..ea04d14 100644 --- a/src/calculus/Ax_mul_Bx.jl +++ b/src/calculus/Ax_mul_Bx.jl @@ -60,13 +60,10 @@ end # Constructors function Ax_mul_Bx(A::AbstractOperator,B::AbstractOperator) - s,t = size(A,1), codomainType(A) - bufA = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)...) - s,t = size(B,1), codomainType(B) - bufB = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)...) - bufC = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)...) - s,t = size(A,2), domainType(A) - bufD = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)...) + bufA = allocateInCodomain(A) + bufB = allocateInCodomain(B) + bufC = allocateInCodomain(B) + bufD = allocateInDomain(A) Ax_mul_Bx(A,B,bufA,bufB,bufC,bufD) end @@ -95,16 +92,16 @@ end size(P::Union{Ax_mul_Bx,Ax_mul_BxJac}) = ((size(P.A,1)[1],size(P.B,1)[2]),size(P.A,2)) -fun_name(L::Union{Ax_mul_Bx,Ax_mul_BxJac}) = fun_name(L.A)*"*"*fun_name(L.B) +fun_name(L::Union{Ax_mul_Bx,Ax_mul_BxJac}) = fun_name(L.A)*"*"*fun_name(L.B) domainType(L::Union{Ax_mul_Bx,Ax_mul_BxJac}) = domainType(L.A) codomainType(L::Union{Ax_mul_Bx,Ax_mul_BxJac}) = codomainType(L.A) # utils -function permute(P::Ax_mul_Bx{L1,L2,C,D}, +function permute(P::Ax_mul_Bx{L1,L2,C,D}, p::AbstractVector{Int}) where {L1,L2,C,D <:ArrayPartition} - Ax_mul_Bx(permute(P.A,p),permute(P.B,p),P.bufA,P.bufB,P.bufC,ArrayPartition(P.bufD.x[p]) ) + Ax_mul_Bx(permute(P.A,p),permute(P.B,p),P.bufA,P.bufB,P.bufC,ArrayPartition(P.bufD.x[p])) end -remove_displacement(P::Ax_mul_Bx) = +remove_displacement(P::Ax_mul_Bx) = Ax_mul_Bx(remove_displacement(P.A), remove_displacement(P.B), P.bufA, P.bufB, P.bufC, P.bufD) diff --git a/src/calculus/Ax_mul_Bxt.jl b/src/calculus/Ax_mul_Bxt.jl index 7e0f834..f2eaa0b 100644 --- a/src/calculus/Ax_mul_Bxt.jl +++ b/src/calculus/Ax_mul_Bxt.jl @@ -38,11 +38,11 @@ struct Ax_mul_Bxt{ bufD::D function Ax_mul_Bxt(A::L1, B::L2, bufA::C, bufB::C, bufC::C, bufD::D) where {L1,L2,C,D} if ndims(A,1) == 1 - if size(A) != size(B) + if size(A) != size(B) throw(DimensionMismatch("Cannot compose operators")) end elseif ndims(A,1) == 2 && ndims(B,1) == 2 && size(A,2) == size(B,2) - if size(A,1)[2] != size(B,1)[2] + if size(A,1)[2] != size(B,1)[2] throw(DimensionMismatch("Cannot compose operators")) end else @@ -68,13 +68,10 @@ end # Constructors function Ax_mul_Bxt(A::AbstractOperator,B::AbstractOperator) - s,t = size(A,1), codomainType(A) - bufA = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)...) - bufC = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)...) - s,t = size(B,1), codomainType(B) - bufB = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)...) - s,t = size(A,2), domainType(A) - bufD = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)...) + bufA = allocateInCodomain(A) + bufB = allocateInCodomain(B) + bufC = allocateInCodomain(A) + bufD = allocateInDomain(A) Ax_mul_Bxt(A,B,bufA,bufB,bufC,bufD) end @@ -103,16 +100,16 @@ end size(P::Union{Ax_mul_Bxt,Ax_mul_BxtJac}) = ((size(P.A,1)[1],size(P.B,1)[1]),size(P.A,2)) -fun_name(L::Union{Ax_mul_Bxt,Ax_mul_BxtJac}) = fun_name(L.A)*"*"*fun_name(L.B) +fun_name(L::Union{Ax_mul_Bxt,Ax_mul_BxtJac}) = fun_name(L.A)*"*"*fun_name(L.B) domainType(L::Union{Ax_mul_Bxt,Ax_mul_BxtJac}) = domainType(L.A) codomainType(L::Union{Ax_mul_Bxt,Ax_mul_BxtJac}) = codomainType(L.A) # utils -function permute(P::Ax_mul_Bxt{L1,L2,C,D}, +function permute(P::Ax_mul_Bxt{L1,L2,C,D}, p::AbstractVector{Int}) where {L1,L2,C,D <:ArrayPartition} Ax_mul_Bxt(permute(P.A,p),permute(P.B,p),P.bufA,P.bufB,P.bufC,ArrayPartition(P.bufD.x[p]) ) end -remove_displacement(P::Ax_mul_Bxt) = +remove_displacement(P::Ax_mul_Bxt) = Ax_mul_Bxt(remove_displacement(P.A), remove_displacement(P.B), P.bufA, P.bufB, P.bufC, P.bufD) diff --git a/src/calculus/Axt_mul_Bx.jl b/src/calculus/Axt_mul_Bx.jl index 2a951b3..dcbc846 100644 --- a/src/calculus/Axt_mul_Bx.jl +++ b/src/calculus/Axt_mul_Bx.jl @@ -38,11 +38,11 @@ struct Axt_mul_Bx{N, bufD::D function Axt_mul_Bx(A::L1, B::L2, bufA::C, bufB::C, bufC::C, bufD::D) where {L1,L2,C,D} if ndims(A,1) == 1 - if size(A) != size(B) + if size(A) != size(B) throw(DimensionMismatch("Cannot compose operators")) end elseif ndims(A,1) == 2 && ndims(B,1) == 2 && size(A,2) == size(B,2) - if size(A,1)[1] != size(B,1)[1] + if size(A,1)[1] != size(B,1)[1] throw(DimensionMismatch("Cannot compose operators")) end else @@ -69,13 +69,10 @@ end # Constructors function Axt_mul_Bx(A::AbstractOperator,B::AbstractOperator) - s,t = size(A,1), codomainType(A) - bufA = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)...) - bufC = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)...) - s,t = size(B,1), codomainType(B) - bufB = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)...) - s,t = size(A,2), domainType(A) - bufD = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)...) + bufA = allocateInCodomain(A) + bufB = allocateInCodomain(B) + bufC = allocateInCodomain(A) + bufD = allocateInDomain(A) Axt_mul_Bx(A,B,bufA,bufB,bufC,bufD) end @@ -122,16 +119,16 @@ end size(P::Union{Axt_mul_Bx{1},Axt_mul_BxJac{1}}) = ((1,),size(P.A,2)) size(P::Union{Axt_mul_Bx{2},Axt_mul_BxJac{2}}) = ((size(P.A,1)[2],size(P.B,1)[2]),size(P.A,2)) -fun_name(L::Union{Axt_mul_Bx,Axt_mul_BxJac}) = fun_name(L.A)*"*"*fun_name(L.B) +fun_name(L::Union{Axt_mul_Bx,Axt_mul_BxJac}) = fun_name(L.A)*"*"*fun_name(L.B) domainType(L::Union{Axt_mul_Bx,Axt_mul_BxJac}) = domainType(L.A) codomainType(L::Union{Axt_mul_Bx,Axt_mul_BxJac}) = codomainType(L.A) # utils -function permute(P::Axt_mul_Bx{N,L1,L2,C,D}, +function permute(P::Axt_mul_Bx{N,L1,L2,C,D}, p::AbstractVector{Int}) where {N,L1,L2,C,D <:ArrayPartition} Axt_mul_Bx(permute(P.A,p),permute(P.B,p),P.bufA,P.bufB,P.bufC,ArrayPartition(P.bufD.x[p]) ) end -remove_displacement(P::Axt_mul_Bx) = +remove_displacement(P::Axt_mul_Bx) = Axt_mul_Bx(remove_displacement(P.A), remove_displacement(P.B), P.bufA, P.bufB, P.bufC, P.bufD) diff --git a/src/calculus/BroadCast.jl b/src/calculus/BroadCast.jl index b1bac4f..be62828 100644 --- a/src/calculus/BroadCast.jl +++ b/src/calculus/BroadCast.jl @@ -21,9 +21,9 @@ julia> B*[1.;2.] ``` """ -struct BroadCast{N, - L <: AbstractOperator, - T <: AbstractArray, +struct BroadCast{N, + L <: AbstractOperator, + T <: AbstractArray, D <: AbstractArray, M, C <: NTuple{M,Colon}, @@ -36,14 +36,14 @@ struct BroadCast{N, cols::C idxs::I - function BroadCast(A::L,dim_out::NTuple{N,Int},bufC::T, bufD::D) where {N, - L<:AbstractOperator, + function BroadCast(A::L,dim_out::NTuple{N,Int},bufC::T, bufD::D) where {N, + L<:AbstractOperator, T<:AbstractArray, D<:AbstractArray } Base.Broadcast.check_broadcast_shape(dim_out,size(A,1)) if size(A,1) != (1,) - M = length(size(A,1)) + M = length(size(A,1)) cols = ([Colon() for i = 1:M]...,) idxs = CartesianIndices((dim_out[M+1:end]...,)) new{N,L,T,D,M,typeof(cols),typeof(idxs)}(A,dim_out,bufC,bufD,cols,idxs) @@ -52,14 +52,13 @@ struct BroadCast{N, idxs = CartesianIndices((1,)) new{N,L,T,D,M,NTuple{0,Colon},typeof(idxs)}(A,dim_out,bufC,bufD,(),idxs) end - end end # Constructors BroadCast(A::L, dim_out::NTuple{N,Int}) where {N,L<:AbstractOperator} = -BroadCast(A, dim_out, zeros(codomainType(A),size(A,1)), zeros(domainType(A),size(A,2)) ) +BroadCast(A, dim_out, allocateInCodomain(A), allocateInDomain(A)) # Mappings @@ -82,7 +81,7 @@ end function mul!(y::CC, A::AdjointOperator{BroadCast{N,L,T,D,0,C,I}}, b::DD) where {N,L,T,D,C,I,CC,DD} R = A.A fill!(y, 0.) - bii = zeros(eltype(b),1) + bii = allocateInCodomain(R.A) for bi in b bii[1] = bi mul!(R.bufD, R.A', bii) @@ -92,7 +91,7 @@ function mul!(y::CC, A::AdjointOperator{BroadCast{N,L,T,D,0,C,I}}, b::DD) where end #TODO make this more general -#length(dim_out) == size(A,1) e.g. a .= b; size(a) = (m,n) size(b) = (1,n) matrix out, column in +#length(dim_out) == size(A,1) e.g. a .= b; size(a) = (m,n) size(b) = (1,n) matrix out, column in function mul!(y::CC, A::AdjointOperator{BroadCast{2,L,T,D,2,C,I}}, b::DD) where {L,T,D,C,I,CC,DD} R = A.A fill!(y, 0.) diff --git a/src/calculus/Compose.jl b/src/calculus/Compose.jl index f05b113..f369989 100644 --- a/src/calculus/Compose.jl +++ b/src/calculus/Compose.jl @@ -3,9 +3,9 @@ export Compose """ `Compose(A::AbstractOperator,B::AbstractOperator)` -Shorthand constructor: +Shorthand constructor: -`A*B` +`A*B` Compose different `AbstractOperator`s. Notice that the domain and codomain of the operators `A` and `B` must match, i.e. `size(A,2) == size(B,1)` and `domainType(A) == codomainType(B)`. @@ -28,19 +28,22 @@ end function Compose(L1::AbstractOperator, L2::AbstractOperator) if size(L1,2) != size(L2,1) - throw(DimensionMismatch("cannot compose operators")) + throw(DimensionMismatch("cannot compose operators with different domain and codomain sizes")) end if domainType(L1) != codomainType(L2) - throw(DomainError()) + throw(DomainError((domainType(L1),codomainType(L2)), "cannot compose operators with different domain and codomain types")) end - Compose( L1, L2, Array{domainType(L1)}(undef,size(L2,1)) ) + if domainStorageType(L1) != codomainStorageType(L2) + throw(DomainError((domainStorageType(L1),codomainStorageType(L2)), "cannot compose operators with different input and output storage types")) + end + Compose(L1, L2, allocateInCodomain(L2)) end Compose(L1::AbstractOperator,L2::AbstractOperator,buf::AbstractArray) = -Compose( (L2,L1), (buf,)) +Compose((L2,L1), (buf,)) Compose(L1::Compose, L2::AbstractOperator,buf::AbstractArray) = -Compose( (L2,L1.A...), (buf,L1.buf...)) +Compose((L2,L1.A...), (buf,L1.buf...)) Compose(L1::AbstractOperator,L2::Compose, buf::AbstractArray) = Compose((L2.A...,L1), (L2.buf...,buf)) diff --git a/src/calculus/DCAT.jl b/src/calculus/DCAT.jl index 53432c3..3767d18 100644 --- a/src/calculus/DCAT.jl +++ b/src/calculus/DCAT.jl @@ -40,14 +40,14 @@ function DCAT(A::Vararg{AbstractOperator}) # build H.idxs idxD,idxC = (),() - for d in (1,2) + for d in (1,2) idxs = [] K = 0 for i in eachindex(ndoms.(A,d)) if ndoms(A[i],d) == 1 # flatten operator K += 1 push!(idxs,K) - else # stacked operator + else # stacked operator idxs = push!(idxs,(collect(K+1:K+ndoms(A[i],d))...,) ) for ii = 1:ndoms(A[i],d) K += 1 @@ -62,37 +62,37 @@ function DCAT(A::Vararg{AbstractOperator}) end # Mappings -@generated function mul!(yy::ArrayPartition, - H::DCAT{N,L,P1,P2}, - bb::ArrayPartition) where {N,L,P1,P2} +@generated function mul!(yy::ArrayPartition, + H::DCAT{N,L,P1,P2}, + bb::ArrayPartition) where {N,L,P1,P2} # extract stuff ex = :(y = yy.x; b = bb.x ) for i = 1:N - if fieldtype(P2,i) <: Int - # flatten operator - # build mul!(y[H.idxC[i]], H.A[i], b) + if fieldtype(P2,i) <: Int + # flatten operator + # build mul!(y[H.idxC[i]], H.A[i], b) yy = :(y[H.idxC[$i]]) else - # stacked operator + # stacked operator # build mul!(( y[H.idxC[i][1]], y[H.idxC[i][2]] ... ), H.A[i], b) yy = [ :(y[H.idxC[$i][$ii]]) for ii in eachindex(fieldnames(fieldtype(P2,i)))] yy = :( ArrayPartition( $(yy...) ) ) end - if fieldtype(P1,i) <: Int - # flatten operator - # build mul!(H.buf, H.A[i], b[H.idxD[i]]) + if fieldtype(P1,i) <: Int + # flatten operator + # build mul!(H.buf, H.A[i], b[H.idxD[i]]) bb = :(b[H.idxD[$i]]) else - # stacked operator + # stacked operator # build mul!(H.buf, H.A[i],( b[H.idxD[i][1]], b[H.idxD[i][2]] ... )) bb = [ :(b[H.idxD[$i][$ii]]) for ii in eachindex(fieldnames(fieldtype(P1,i))) ] bb = :( ArrayPartition( $(bb...) ) ) end - + ex = :($ex; mul!($yy,H.A[$i],$bb)) end @@ -101,37 +101,37 @@ end end -@generated function mul!(yy::ArrayPartition, +@generated function mul!(yy::ArrayPartition, A::AdjointOperator{DCAT{N,L,P1,P2}}, - bb::ArrayPartition) where {N,L,P1,P2} + bb::ArrayPartition) where {N,L,P1,P2} # extract stuff ex = :(H = A.A; y = yy.x; b = bb.x ) for i = 1:N - if fieldtype(P1,i) <: Int - # flatten operator - # build mul!(y[H.idxD[i]], H.A[i]', b) + if fieldtype(P1,i) <: Int + # flatten operator + # build mul!(y[H.idxD[i]], H.A[i]', b) yy = :(y[H.idxD[$i]]) else - # stacked operator + # stacked operator # build mul!(( y[H.idxD[i][1]], y[H.idxD[i][2]] ... ), H.A[i]', b) yy = [ :(y[H.idxD[$i][$ii]]) for ii in eachindex(fieldnames(fieldtype(P1,i)))] yy = :( ArrayPartition( $(yy...) )) end - if fieldtype(P2,i) <: Int - # flatten operator - # build mul!(H.buf, H.A[i]', b[H.idxC[i]]) + if fieldtype(P2,i) <: Int + # flatten operator + # build mul!(H.buf, H.A[i]', b[H.idxC[i]]) bb = :(b[H.idxC[$i]]) else - # stacked operator + # stacked operator # build mul!(H.buf, H.A[i]',( b[H.idxC[i][1]], b[H.idxC[i][2]] ... )) bb = [ :(b[H.idxC[$i][$ii]]) for ii in eachindex(fieldnames(fieldtype(P2,i)))] bb = :( ArrayPartition( $(bb...) ) ) end - + ex = :($ex; mul!($yy,H.A[$i]',$bb)) end @@ -141,13 +141,13 @@ end end # Properties -size(H::DCAT) = size(H,1),size(H,2) +size(H::DCAT) = size(H,1),size(H,2) -function size(H::DCAT, i::Int) +function size(H::DCAT, i::Int) sz = [] for s in size.(H.A,i) - eltype(s) <: Int ? push!(sz,s) : push!(sz,s...) + eltype(s) <: Int ? push!(sz,s) : push!(sz,s...) end p = vcat([[idx... ] for idx in (i == 1 ? H.idxC : H.idxD) ]...) invpermute!(sz,p) @@ -158,13 +158,13 @@ end fun_name(L::DCAT) = length(L.A) == 2 ? "["*fun_name(L.A[1])*",0;0,"*fun_name(L.A[2])*"]" : "DCAT" -function domainType(H::DCAT) +function domainType(H::DCAT) domain = vcat([typeof(d)<:Tuple ? [d...] : d for d in domainType.(H.A)]...) p = vcat([[idx... ] for idx in H.idxD]...) invpermute!(domain,p) return (domain...,) end -function codomainType(H::DCAT) +function codomainType(H::DCAT) codomain = vcat([typeof(d)<:Tuple ? [d...] : d for d in codomainType.(H.A)]...) p = vcat([[idx... ] for idx in H.idxC]...) invpermute!(codomain,p) @@ -185,7 +185,7 @@ is_full_column_rank(L::DCAT) = all(is_full_column_rank.(L.A)) function permute(H::DCAT{N,L,P1,P2}, p::AbstractVector{Int}) where {N,L,P1,P2} - unfolded = vcat([[idx... ] for idx in H.idxD]...) + unfolded = vcat([[idx... ] for idx in H.idxD]...) invpermute!(unfolded,p) new_part = () diff --git a/src/calculus/HCAT.jl b/src/calculus/HCAT.jl index 3a9ab88..1c1df51 100644 --- a/src/calculus/HCAT.jl +++ b/src/calculus/HCAT.jl @@ -3,11 +3,11 @@ export HCAT """ `HCAT(A::AbstractOperator...)` -Shorthand constructors: +Shorthand constructors: -`[A1 A2 ...]` +`[A1 A2 ...]` -`hcat(A...)` +`hcat(A...)` Horizontally concatenate `AbstractOperator`s. Notice that all the operators must share the same codomain dimensions and type, e.g. `size(A1,1) == size(A2,1)` and `codomainType(A1) == codomainType(A2)`. @@ -23,7 +23,7 @@ HCAT ℝ^10 ℝ^10 ℝ^10 -> ℝ^10 ``` -To evaluate `HCAT` operators multiply them with a `Tuple` of `AbstractArray` of the correct dimensions and type. +To evaluate `HCAT` operators multiply them with a `Tuple` of `AbstractArray` of the correct dimensions and type. ```julia julia> H*ArrayPartition(ones(10),ones(10)) @@ -35,20 +35,20 @@ julia> H*ArrayPartition(ones(10),ones(10)) ``` """ -struct HCAT{N, # number of AbstractOperator +struct HCAT{N, # number of AbstractOperator L <: NTuple{N,AbstractOperator}, P <: Tuple, C <: AbstractArray, } <: AbstractOperator A::L # tuple of AbstractOperators - idxs::P # indices - # H = HCAT(Eye(n),HCAT(Eye(n),Eye(n))) has H.idxs = (1,2,3) + idxs::P # indices + # H = HCAT(Eye(n),HCAT(Eye(n),Eye(n))) has H.idxs = (1,2,3) # `AbstractOperators` are flatten - # H = HCAT(Eye(n),Compose(MatrixOp(randn(n,n)),HCAT(Eye(n),Eye(n)))) + # H = HCAT(Eye(n),Compose(MatrixOp(randn(n,n)),HCAT(Eye(n),Eye(n)))) # has H.idxs = (1,(2,3)) # `AbstractOperators` are stack buf::C # buffer memory - function HCAT(A::L, idxs::P, buf::C) where {N, + function HCAT(A::L, idxs::P, buf::C) where {N, L <: NTuple{N,AbstractOperator}, P <: Tuple, C} @@ -66,7 +66,7 @@ function HCAT(A::Vararg{AbstractOperator}) if any((<:).(typeof.(A),HCAT)) #there are HCATs in A AA = () for a in A - if typeof(a) <: HCAT # flatten + if typeof(a) <: HCAT # flatten AA = (AA...,a.A...) else # stack AA = (AA...,a) @@ -74,18 +74,16 @@ function HCAT(A::Vararg{AbstractOperator}) end # use buffer from HCAT in A buf = A[findfirst( (<:).(typeof.(A),HCAT) ) ].buf - else + else AA = A - s = size(AA[1],1) - t = codomainType(AA[1]) # generate buffer - buf = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)) + buf = allocateInCodomain(AA[1]) end return HCAT(AA, buf) end -function HCAT(AA::NTuple{N,AbstractOperator}, buf::C) where {N,C} +function HCAT(AA::NTuple{N,AbstractOperator}, buf::C) where {N,C} if N == 1 return AA[1] else @@ -96,7 +94,7 @@ function HCAT(AA::NTuple{N,AbstractOperator}, buf::C) where {N,C} if ndoms(AA[i],2) == 1 # flatten operator K += 1 push!(idxs,K) - else # stacked operator + else # stacked operator idxs = push!(idxs,(collect(K+1:K+ndoms(AA[i],2))...,) ) for ii = 1:ndoms(AA[i],2) K += 1 @@ -113,12 +111,12 @@ HCAT(A::AbstractOperator) = A @generated function mul!(y::C, H::HCAT{N,L,P,C}, b::DD) where {N,L,P,C,DD <: ArrayPartition} ex = :() - if fieldtype(P,1) <: Int - # flatten operator - # build mul!(y, H.A[1], b.x[H.idxs[1]]) + if fieldtype(P,1) <: Int + # flatten operator + # build mul!(y, H.A[1], b.x[H.idxs[1]]) bb = :(b.x[H.idxs[1]]) else - # stacked operator + # stacked operator # build mul!(y, H.A[1],ArrayPartition( b.x[H.idxs[1][1]], b.x[H.idxs[1][2]] ... )) bb = [ :(b.x[H.idxs[1][$ii]]) for ii in eachindex(fieldnames(fieldtype(P,1)))] bb = :( ArrayPartition($(bb...)) ) @@ -126,12 +124,12 @@ HCAT(A::AbstractOperator) = A ex = :($ex; mul!(y,H.A[1],$bb)) # write on y for i = 2:N - if fieldtype(P,i) <: Int - # flatten operator - # build mul!(H.buf, H.A[i], b.x[H.idxs[i]]) + if fieldtype(P,i) <: Int + # flatten operator + # build mul!(H.buf, H.A[i], b.x[H.idxs[i]]) bb = :(b.x[H.idxs[$i]]) else - # stacked operator + # stacked operator # build mul!(H.buf, H.A[i],( b.x[H.idxs[i][1]], b.x[H.idxs[i][2]] ... )) bb = [ :( b.x[H.idxs[$i][$ii]] ) for ii in eachindex(fieldnames(fieldtype(P,i)))] bb = :( ArrayPartition( $(bb...) ) ) @@ -147,12 +145,12 @@ end @generated function mul!(y::DD, A::AdjointOperator{HCAT{N,L,P,C}}, b::C) where {N,L,P,C,DD <: ArrayPartition} ex = :(H = A.A) for i = 1:N - if fieldtype(P,i) <: Int - # flatten operator - # build mul!(y.x[H.idxs[i]], H.A[i]', b) + if fieldtype(P,i) <: Int + # flatten operator + # build mul!(y.x[H.idxs[i]], H.A[i]', b) yy = :(y.x[H.idxs[$i]]) else - # stacked operator + # stacked operator # build mul!(ArrayPartition( y[.xH.idxs[i][1]], y.x[H.idxs[i][2]] ... ), H.A[i]', b) yy = [ :(y.x[H.idxs[$i][$ii]]) for ii in eachindex(fieldnames(fieldtype(P,i)))] yy = :(ArrayPartition( $(yy...) ) ) @@ -167,12 +165,12 @@ end @generated function mul_skipZeros!(y::C, H::HCAT{N,L,P,C}, b::DD) where {N,L,P,C,DD <: ArrayPartition} ex = :() - if fieldtype(P,1) <: Int - # flatten operator - # build mul!(y, H.A[1], b.x[H.idxs[1]]) + if fieldtype(P,1) <: Int + # flatten operator + # build mul!(y, H.A[1], b.x[H.idxs[1]]) bb = :(b.x[H.idxs[1]]) else - # stacked operator + # stacked operator # build mul!(y, H.A[1],ArrayPartition( b.x[H.idxs[1][1]], b.x[H.idxs[1][2]] ... )) bb = [ :(b.x[H.idxs[1][$ii]]) for ii in eachindex(fieldnames(fieldtype(P,1)))] bb = :( ArrayPartition($(bb...)) ) @@ -181,12 +179,12 @@ end for i = 2:N if !(fieldtype(L,i) <: Zeros) - if fieldtype(P,i) <: Int - # flatten operator - # build mul!(H.buf, H.A[i], b.x[H.idxs[i]]) + if fieldtype(P,i) <: Int + # flatten operator + # build mul!(H.buf, H.A[i], b.x[H.idxs[i]]) bb = :(b.x[H.idxs[$i]]) else - # stacked operator + # stacked operator # build mul!(H.buf, H.A[i],( b.x[H.idxs[i][1]], b.x[H.idxs[i][2]] ... )) bb = [ :( b.x[H.idxs[$i][$ii]] ) for ii in eachindex(fieldnames(fieldtype(P,i)))] bb = :( ArrayPartition( $(bb...) ) ) @@ -204,12 +202,12 @@ end ex = :(H = A.A) for i = 1:N if !(fieldtype(L,i) <: Zeros) - if fieldtype(P,i) <: Int - # flatten operator - # build mul!(y.x[H.idxs[i]], H.A[i]', b) + if fieldtype(P,i) <: Int + # flatten operator + # build mul!(y.x[H.idxs[i]], H.A[i]', b) yy = :(y.x[H.idxs[$i]]) else - # stacked operator + # stacked operator # build mul!(ArrayPartition( y[.xH.idxs[i][1]], y.x[H.idxs[i][2]] ... ), H.A[i]', b) yy = [ :(y.x[H.idxs[$i][$ii]]) for ii in eachindex(fieldnames(fieldtype(P,i)))] yy = :(ArrayPartition( $(yy...) ) ) @@ -223,10 +221,10 @@ end # Properties -function size(H::HCAT) +function size(H::HCAT) size_in = [] for s in size.(H.A,2) - eltype(s) <: Int ? push!(size_in,s) : push!(size_in,s...) + eltype(s) <: Int ? push!(size_in,s) : push!(size_in,s...) end p = vcat([[idx... ] for idx in H.idxs]...) invpermute!(size_in,p) @@ -236,7 +234,7 @@ end fun_name(L::HCAT) = length(L.A) == 2 ? "["*fun_name(L.A[1])*","*fun_name(L.A[2])*"]" : "HCAT" -function domainType(H::HCAT) +function domainType(H::HCAT) domain = vcat([typeof(d)<:Tuple ? [d...] : d for d in domainType.(H.A)]...) p = vcat([[idx... ] for idx in H.idxs]...) invpermute!(domain,p) @@ -252,7 +250,7 @@ diag_AAc(L::HCAT) = (+).(diag_AAc.(L.A)...) # utils function permute(H::HCAT, p::AbstractVector{Int}) - unfolded = vcat([[idx... ] for idx in H.idxs]...) + unfolded = vcat([[idx... ] for idx in H.idxs]...) invpermute!(unfolded,p) new_part = () diff --git a/src/calculus/HadamardProd.jl b/src/calculus/HadamardProd.jl index b3e8ad1..862a978 100644 --- a/src/calculus/HadamardProd.jl +++ b/src/calculus/HadamardProd.jl @@ -37,7 +37,7 @@ struct HadamardProd{ bufB::C bufD::D function HadamardProd(A::L1, B::L2, bufA::C, bufB::C, bufD::D) where {L1,L2,C,D} - if size(A) != size(B) + if size(A) != size(B) throw(DimensionMismatch("Cannot compose operators")) end new{L1,L2,C,D}(A,B,bufA,bufB,bufD) @@ -59,12 +59,9 @@ end # Constructors function HadamardProd(A::AbstractOperator,B::AbstractOperator) - s,t = size(A,1), codomainType(A) - bufA = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)...) - s,t = size(B,1), codomainType(B) - bufB = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)...) - s,t = size(A,2), domainType(A) - bufD = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)...) + bufA = allocateInCodomain(A) + bufB = allocateInCodomain(B) + bufD = allocateInDomain(A) HadamardProd(A,B,bufA,bufB,bufD) end @@ -94,16 +91,16 @@ end size(P::Union{HadamardProd,HadamardProdJac}) = (size(P.A,1),size(P.A,2)) -fun_name(L::Union{HadamardProd,HadamardProdJac}) = fun_name(L.A)*".*"*fun_name(L.B) +fun_name(L::Union{HadamardProd,HadamardProdJac}) = fun_name(L.A)*".*"*fun_name(L.B) domainType(L::Union{HadamardProd,HadamardProdJac}) = domainType(L.A) codomainType(L::Union{HadamardProd,HadamardProdJac}) = codomainType(L.A) # utils -function permute(P::HadamardProd{L1,L2,C,D}, +function permute(P::HadamardProd{L1,L2,C,D}, p::AbstractVector{Int}) where {L1,L2,C,D <:ArrayPartition} HadamardProd(permute(P.A,p),permute(P.B,p),P.bufA,P.bufB,ArrayPartition(P.bufD.x[p]) ) end -remove_displacement(P::HadamardProd) = +remove_displacement(P::HadamardProd) = HadamardProd(remove_displacement(P.A), remove_displacement(P.B), P.bufA, P.bufB, P.bufD) diff --git a/src/calculus/Jacobian.jl b/src/calculus/Jacobian.jl index 9f599d9..b051a5f 100644 --- a/src/calculus/Jacobian.jl +++ b/src/calculus/Jacobian.jl @@ -3,7 +3,7 @@ export Jacobian """ `Jacobian(A::AbstractOperator,x)` -Shorthand constructor: +Shorthand constructor: `jacobian(A::AbstractOperator,x)` @@ -24,61 +24,61 @@ struct Jacobian{T <: NonLinearOperator, X<:AbstractArray} <: LinearOperator x::X end -#Jacobian of LinearOperator +#Jacobian of LinearOperator Jacobian(L::T,x::X) where {T<:LinearOperator,X<:AbstractArray} = L #Jacobian of Scale -Jacobian(S::T2,x::AbstractArray) where {T,L,T2<:Scale{T,L}} = Scale(S.coeff,Jacobian(S.A,x)) +Jacobian(S::T2,x::AbstractArray) where {T,L,T2<:Scale{T,L}} = Scale(S.coeff,Jacobian(S.A,x)) ##Jacobian of DCAT -function Jacobian(H::DCAT{N,L,P1,P2},b) where {N,L,P1,P2} +function Jacobian(H::DCAT{N,L,P1,P2},b) where {N,L,P1,P2} x = b.x A = () c = 0 for (k, idx) in enumerate(H.idxD) - if length(idx) == 1 - A = (A...,jacobian(H.A[k],x[idx])) + if length(idx) == 1 + A = (A...,jacobian(H.A[k],x[idx])) else xx = ([x[i] for i in idx]...,) - A = (A...,jacobian(H.A[k],xx)) + A = (A...,jacobian(H.A[k],xx)) end end DCAT(A,H.idxD,H.idxC) end #Jacobian of HCAT -function Jacobian(H::HCAT{N,L,P,C},b::ArrayPartition) where {N,L,P,C,D} +function Jacobian(H::HCAT{N,L,P,C},b::ArrayPartition) where {N,L,P,C} x = b.x A = () for (k, idx) in enumerate(H.idxs) - if length(idx) == 1 - A = (A...,jacobian(H.A[k],x[idx])) + if length(idx) == 1 + A = (A...,jacobian(H.A[k],x[idx])) else xx = ArrayPartition([x[i] for i in idx]...,) - A = (A...,jacobian(H.A[k],xx)) + A = (A...,jacobian(H.A[k],xx)) end end HCAT(A,H.idxs,H.buf) end #Jacobian of VCAT -function Jacobian(V::VCAT{N,L,P,C},x::AbstractArray) where {N,L,P,C,D} +function Jacobian(V::VCAT{N,L,P,C},x::AbstractArray) where {N,L,P,C} JJ = ([Jacobian(a,x) for a in V.A]...,) VCAT(JJ, V.idxs, V.buf) end -#Jacobian of Compose -function Jacobian(L::Compose, x::X) where {X<:AbstractArray} +#Jacobian of Compose +function Jacobian(L::Compose, x::X) where {X<:AbstractArray} Compose(Jacobian.(L.A,(x,L.buf...)),L.buf) end -function Jacobian(L::Compose, x::X) where {N,X<:NTuple{N,AbstractArray}} +function Jacobian(L::Compose, x::X) where {N,X<:NTuple{N,AbstractArray}} Compose(Jacobian.(L.A,(x,L.buf...)),L.buf) end #Jacobian of Reshape -Jacobian(R::Reshape{N,L},x::AbstractArray) where {N,L} = Reshape(Jacobian(R.A,x),R.dim_out) +Jacobian(R::Reshape{N,L},x::AbstractArray) where {N,L} = Reshape(Jacobian(R.A,x),R.dim_out) #Jacobian of Sum -Jacobian(S::Sum{K,C,D},x::D) where {K,C,D} = +Jacobian(S::Sum{K,C,D},x::D) where {K,C,D} = Sum(([Jacobian(a,x) for a in S.A]...,),S.bufC,S.bufD) #Jacobian of Transpose -Jacobian(T::Transpose{A}, x::AbstractArray) where {A <: AbstractOperator} = T +Jacobian(T::Transpose{A}, x::AbstractArray) where {A <: AbstractOperator} = T #Jacobian of BroadCast -Jacobian(B::A, x::AbstractArray) where {A <: BroadCast} = BroadCast(Jacobian(B.A,x),B.dim_out) +Jacobian(B::A, x::AbstractArray) where {A <: BroadCast} = BroadCast(Jacobian(B.A,x),B.dim_out) #Jacobian of AffineAdd Jacobian(B::A, x) where {A <: AffineAdd} = Jacobian(B.A,x) diff --git a/src/calculus/Reshape.jl b/src/calculus/Reshape.jl index 25d5ef2..928e9d8 100644 --- a/src/calculus/Reshape.jl +++ b/src/calculus/Reshape.jl @@ -4,9 +4,9 @@ export Reshape """ `Reshape(A::AbstractOperator, dim_out...)` -Shorthand constructor: +Shorthand constructor: -`reshape(A, idx...)` +`reshape(A, idx...)` Reshape the codomain dimensions of an `AbstractOperator`. @@ -61,11 +61,11 @@ is_linear( R::Reshape) = is_linear(R.A) is_null( R::Reshape) = is_null(R.A) is_eye( R::Reshape) = is_eye(R.A) is_diagonal( R::Reshape) = is_diagonal(R.A) -is_AcA_diagonal(R::Reshape) = is_AcA_diagonal(R.A) +is_AcA_diagonal(R::Reshape) = is_AcA_diagonal(R.A) is_AAc_diagonal(R::Reshape) = is_AAc_diagonal(R.A) is_orthogonal( R::Reshape) = is_orthogonal( R.A) is_invertible( R::Reshape) = is_invertible(R.A) -is_full_row_rank( R::Reshape) = is_full_row_rank( R.A) +is_full_row_rank( R::Reshape) = is_full_row_rank( R.A) is_full_column_rank( R::Reshape) = is_full_column_rank( R.A) fun_name(R::Reshape) = "¶"*fun_name(R.A) @@ -73,5 +73,5 @@ remove_displacement(R::Reshape) = Reshape(remove_displacement(R.A), R.dim_out) function permute(T::Reshape{N,L}, p::AbstractVector{Int}) where {N,L} A = AbstractOperators.permute(T.A,p) - return Reshape(A,T.dim_out) + return Reshape(A,T.dim_out) end diff --git a/src/calculus/Scale.jl b/src/calculus/Scale.jl index 4e1876b..36e2270 100644 --- a/src/calculus/Scale.jl +++ b/src/calculus/Scale.jl @@ -3,9 +3,9 @@ export Scale """ `Scale(α::Number,A::AbstractOperator)` -Shorthand constructor: +Shorthand constructor: -`*(α::Number,A::AbstractOperator)` +`*(α::Number,A::AbstractOperator)` Scale an `AbstractOperator` by a factor of `α`. @@ -16,7 +16,7 @@ julia> A = FiniteDiff((10,2)) julia> S = Scale(10,A) αδx ℝ^(10, 2) -> ℝ^(9, 2) -julia> 10*A #shorthand +julia> 10*A #shorthand αℱc ℝ^10 -> ℝ^10 ``` @@ -30,11 +30,11 @@ end # Constructors -function Scale(coeff::T, L::R) where {T <: RealOrComplex, R <: AbstractOperator} +function Scale(coeff::T, L::R) where {T <: RealOrComplex, R <: AbstractOperator} coeff_conj = conj(coeff) coeff, coeff_conj = promote(coeff, coeff_conj) cT = codomainType(L) - isCodomainReal = typeof(cT) <: Tuple ? all([t <: Real for t in cT]) : cT <: Real + isCodomainReal = typeof(cT) <: Tuple ? all([t <: Real for t in cT]) : cT <: Real if isCodomainReal && T <: Complex error("Cannot Scale AbstractOperator with real codomain with complex scalar. Use `DiagOp` instead.") end @@ -43,10 +43,10 @@ end # Special Constructors # scale of scale -Scale(coeff::T2, L::S) where {T1 <: RealOrComplex, - T2 <: RealOrComplex, - R <: AbstractOperator, - S <: Scale{T1, R}}= +Scale(coeff::T2, L::S) where {T1 <: RealOrComplex, + T2 <: RealOrComplex, + R <: AbstractOperator, + S <: Scale{T1, R}}= Scale(*(promote(coeff,L.coeff)...), L.A) # scale of DiagOp Scale(coeff::T,L::DiagOp) where {T<:RealOrComplex} = DiagOp(coeff*diag(L)) @@ -91,7 +91,7 @@ is_null(L::Scale) = is_null(L.A) is_eye(L::Scale) = is_diagonal(L.A) is_diagonal(L::Scale) = is_diagonal(L.A) is_invertible(L::Scale) = L.coeff == 0 ? false : is_invertible(L.A) -is_AcA_diagonal(L::Scale) = is_AcA_diagonal(L.A) +is_AcA_diagonal(L::Scale) = is_AcA_diagonal(L.A) is_AAc_diagonal(L::Scale) = is_AAc_diagonal(L.A) is_full_row_rank(L::Scale) = is_full_row_rank(L.A) is_full_column_rank(L::Scale) = is_full_column_rank(L.A) diff --git a/src/calculus/Sum.jl b/src/calculus/Sum.jl index d8db14d..9565fd0 100644 --- a/src/calculus/Sum.jl +++ b/src/calculus/Sum.jl @@ -22,14 +22,8 @@ end Sum(A::AbstractOperator) = A function Sum(A::Vararg{AbstractOperator}) - s = size(A[1],1) - t = codomainType(A[1]) - bufC = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)...) - - s = size(A[1],2) - t = domainType(A[1]) - bufD = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)...) - + bufC = allocateInCodomain(A[1]) + bufD = allocateInDomain(A[1]) return Sum(A, bufC, bufD) end @@ -85,17 +79,17 @@ fun_name(S::Sum) = length(S.A) == 2 ? fun_name(S.A[1])*"+"*fun_name(S.A[2]) : "Σ" -is_linear(L::Sum) = all(is_linear.(L.A)) -is_null(L::Sum) = all(is_null.(L.A)) -is_diagonal(L::Sum) = all(is_diagonal.(L.A)) +is_linear(L::Sum) = all(is_linear.(L.A)) +is_null(L::Sum) = all(is_null.(L.A)) +is_diagonal(L::Sum) = all(is_diagonal.(L.A)) is_full_row_rank(L::Sum) = any(is_full_row_rank.(L.A)) is_full_column_rank(L::Sum) = any(is_full_column_rank.(L.A)) diag(L::Sum) = (+).(diag.(L.A)...,) # utils -function permute(S::Sum, p::AbstractVector{Int}) - AA = ([permute(A,p) for A in S.A]...,) +function permute(S::Sum, p::AbstractVector{Int}) + AA = ([permute(A,p) for A in S.A]...,) return Sum(AA, S.bufC, ArrayPartition(S.bufD.x[p]...)) end diff --git a/src/calculus/VCAT.jl b/src/calculus/VCAT.jl index 0749d59..0a9c5bc 100644 --- a/src/calculus/VCAT.jl +++ b/src/calculus/VCAT.jl @@ -4,11 +4,11 @@ export VCAT """ `VCAT(A::AbstractOperator...)` -Shorthand constructors: +Shorthand constructors: -`[A1; A2 ...]` +`[A1; A2 ...]` -`vcat(A...)` +`vcat(A...)` Vertically concatenate `AbstractOperator`s. Notice that all the operators must share the same domain dimensions and type, e.g. `size(A1,2) == size(A2,2)` and `domainType(A1) == domainType(A2)`. @@ -32,20 +32,20 @@ julia> V*ones(3) ``` """ -struct VCAT{N, # number of AbstractOperator +struct VCAT{N, # number of AbstractOperator L <: NTuple{N,AbstractOperator}, P <: Tuple, C <: AbstractArray, } <: AbstractOperator A::L # tuple of AbstractOperators - idxs::P # indices - # H = VCAT(Eye(n),VCAT(Eye(n),Eye(n))) has H.idxs = (1,2,3) + idxs::P # indices + # H = VCAT(Eye(n),VCAT(Eye(n),Eye(n))) has H.idxs = (1,2,3) # `AbstractOperators` are flatten - # H = VCAT(Eye(n),Compose(MatrixOp(randn(n,n)),VCAT(Eye(n),Eye(n)))) + # H = VCAT(Eye(n),Compose(MatrixOp(randn(n,n)),VCAT(Eye(n),Eye(n)))) # has H.idxs = (1,(2,3)) # `AbstractOperators` are stack buf::C # buffer memory - function VCAT(A::L, idxs::P, buf::C) where {N, + function VCAT(A::L, idxs::P, buf::C) where {N, L <: NTuple{N,AbstractOperator}, P <: Tuple, C <: AbstractArray, @@ -64,7 +64,7 @@ function VCAT(A::Vararg{AbstractOperator}) if any((<:).(typeof.(A),VCAT)) #there are VCATs in A AA = () for a in A - if typeof(a) <: VCAT # flatten + if typeof(a) <: VCAT # flatten AA = (AA...,a.A...) else # stack AA = (AA...,a) @@ -72,18 +72,16 @@ function VCAT(A::Vararg{AbstractOperator}) end # use buffer from VCAT in A buf = A[findfirst( (<:).(typeof.(A),VCAT) ) ].buf - else + else AA = A - s = size(AA[1],2) - t = domainType(AA[1]) # generate buffer - buf = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)) + buf = allocateInDomain(AA[1]) end return VCAT(AA, buf) end -function VCAT(AA::NTuple{N,AbstractOperator}, buf::C) where {N,C} +function VCAT(AA::NTuple{N,AbstractOperator}, buf::C) where {N,C} if N == 1 return AA[1] else @@ -94,7 +92,7 @@ function VCAT(AA::NTuple{N,AbstractOperator}, buf::C) where {N,C} if ndoms(AA[i],1) == 1 # flatten operator K += 1 push!(idxs,K) - else # stacked operator + else # stacked operator idxs = push!(idxs,(collect(K+1:K+ndoms(AA[i],2))...,) ) for ii = 1:ndoms(AA[i],1) K += 1 @@ -111,12 +109,12 @@ VCAT(A::AbstractOperator) = A @generated function mul!(y::C, A::AdjointOperator{VCAT{N,L,P,C}}, b::DD) where {N,L,P,C,DD <: ArrayPartition} ex = :(H = A.A) - if fieldtype(P,1) <: Int - # flatten operator - # build mul!(y, H.A[1]', b.x[H.idxs[1]]) + if fieldtype(P,1) <: Int + # flatten operator + # build mul!(y, H.A[1]', b.x[H.idxs[1]]) bb = :(b.x[H.idxs[1]]) else - # stacked operator + # stacked operator # build mul!(y, H.A[1]',ArrayPartition( b.x[H.idxs[1][1]], b.x[H.idxs[1][2]] ... )) bb = [ :(b.x[H.idxs[1][$ii]]) for ii in eachindex(fieldnames(fieldtype(P,1)))] bb = :( ArrayPartition($(bb...)) ) @@ -124,12 +122,12 @@ VCAT(A::AbstractOperator) = A ex = :($ex; mul!(y,H.A[1]',$bb)) # write on y for i = 2:N - if fieldtype(P,i) <: Int - # flatten operator - # build mul!(H.buf, H.A[i]', b.x[H.idxs[i]]) + if fieldtype(P,i) <: Int + # flatten operator + # build mul!(H.buf, H.A[i]', b.x[H.idxs[i]]) bb = :(b.x[H.idxs[$i]]) else - # stacked operator + # stacked operator # build mul!(H.buf, H.A[i]',( b.x[H.idxs[i][1]], b.x[H.idxs[i][2]] ... )) bb = [ :( b.x[H.idxs[$i][$ii]] ) for ii in eachindex(fieldnames(fieldtype(P,i)))] bb = :( ArrayPartition( $(bb...) ) ) @@ -145,12 +143,12 @@ end @generated function mul!(y::DD, H::VCAT{N,L,P,C}, b::C) where {N,L,P,C,DD <: ArrayPartition} ex = :() for i = 1:N - if fieldtype(P,i) <: Int - # flatten operator - # build mul!(y.x[H.idxs[i]], H.A[i], b) + if fieldtype(P,i) <: Int + # flatten operator + # build mul!(y.x[H.idxs[i]], H.A[i], b) yy = :(y.x[H.idxs[$i]]) else - # stacked operator + # stacked operator # build mul!(ArrayPartition( y[.xH.idxs[i][1]], y.x[H.idxs[i][2]] ... ), H.A[i], b) yy = [ :(y.x[H.idxs[$i][$ii]]) for ii in eachindex(fieldnames(fieldtype(P,i)))] yy = :(ArrayPartition( $(yy...) ) ) @@ -165,12 +163,12 @@ end @generated function mul_skipZeros!(y::C, A::AdjointOperator{VCAT{N,L,P,C}}, b::DD) where {N,L,P,C,DD <: ArrayPartition} ex = :(H = A.A) - if fieldtype(P,1) <: Int - # flatten operator - # build mul!(y, H.A[1]', b.x[H.idxs[1]]) + if fieldtype(P,1) <: Int + # flatten operator + # build mul!(y, H.A[1]', b.x[H.idxs[1]]) bb = :(b.x[H.idxs[1]]) else - # stacked operator + # stacked operator # build mul!(y, H.A[1]',ArrayPartition( b.x[H.idxs[1][1]], b.x[H.idxs[1][2]] ... )) bb = [ :(b.x[H.idxs[1][$ii]]) for ii in eachindex(fieldnames(fieldtype(P,1)))] bb = :( ArrayPartition($(bb...)) ) @@ -179,12 +177,12 @@ end for i = 2:N if !(fieldtype(L,i) <: Zeros) - if fieldtype(P,i) <: Int - # flatten operator - # build mul!(H.buf, H.A[i]', b.x[H.idxs[i]]) + if fieldtype(P,i) <: Int + # flatten operator + # build mul!(H.buf, H.A[i]', b.x[H.idxs[i]]) bb = :(b.x[H.idxs[$i]]) else - # stacked operator + # stacked operator # build mul!(H.buf, H.A[i]',( b.x[H.idxs[i][1]], b.x[H.idxs[i][2]] ... )) bb = [ :( b.x[H.idxs[$i][$ii]] ) for ii in eachindex(fieldnames(fieldtype(P,i)))] bb = :( ArrayPartition( $(bb...) ) ) @@ -202,12 +200,12 @@ end ex = :() for i = 1:N if !(fieldtype(L,i) <: Zeros) - if fieldtype(P,i) <: Int - # flatten operator - # build mul!(y.x[H.idxs[i]], H.A[i], b) + if fieldtype(P,i) <: Int + # flatten operator + # build mul!(y.x[H.idxs[i]], H.A[i], b) yy = :(y.x[H.idxs[$i]]) else - # stacked operator + # stacked operator # build mul!(ArrayPartition( y[.xH.idxs[i][1]], y.x[H.idxs[i][2]] ... ), H.A[i], b) yy = [ :(y.x[H.idxs[$i][$ii]]) for ii in eachindex(fieldnames(fieldtype(P,i)))] yy = :(ArrayPartition( $(yy...) ) ) @@ -221,10 +219,10 @@ end # Properties -function size(H::VCAT) +function size(H::VCAT) size_out = [] for s in size.(H.A,1) - eltype(s) <: Int ? push!(size_out,s) : push!(size_out,s...) + eltype(s) <: Int ? push!(size_out,s) : push!(size_out,s...) end p = vcat([[idx... ] for idx in H.idxs]...) invpermute!(size_out,p) @@ -234,7 +232,7 @@ end fun_name(L::VCAT) = length(L.A) == 2 ? "["*fun_name(L.A[1])*";"*fun_name(L.A[2])*"]" : "VCAT" -function codomainType(H::VCAT) +function codomainType(H::VCAT) codomain = vcat([typeof(d)<:Tuple ? [d...] : d for d in codomainType.(H.A)]...) p = vcat([[idx... ] for idx in H.idxs]...) invpermute!(codomain,p) @@ -252,7 +250,7 @@ diag_AcA(L::VCAT) = (+).(diag_AcA.(L.A)...,) function permute(H::VCAT{N,L,P,C}, p::AbstractVector{Int}) where {N,L,P,C} - unfolded = vcat([[idx... ] for idx in H.idxs]...) + unfolded = vcat([[idx... ] for idx in H.idxs]...) invpermute!(unfolded,p) new_part = () diff --git a/src/linearoperators/Conv.jl b/src/linearoperators/Conv.jl index e5dfd1f..6dfc98c 100644 --- a/src/linearoperators/Conv.jl +++ b/src/linearoperators/Conv.jl @@ -5,12 +5,12 @@ export Conv `Conv(x::AbstractVector, h::AbstractVector)` -Creates a `LinearOperator` which, when multiplied with an array `x::AbstractVector`, returns the convolution between `x` and `h`. Uses `conv` and hence FFT algorithm. +Creates a `LinearOperator` which, when multiplied with an array `x::AbstractVector`, returns the convolution between `x` and `h`. Uses `conv` and hence FFT algorithm. """ struct Conv{T, H <: AbstractVector{T}, - Hc <: AbstractVector{Complex{T}}, + Hc <: AbstractVector, } <: LinearOperator dim_in::Tuple{Int} h::H @@ -23,15 +23,25 @@ end # Constructors +isTypeReal(::Type{T}) where {T} = T <: Real + ###standard constructor -function Conv(DomainType::Type, dim_in::Tuple{Int}, h::H) where {H<:AbstractVector} +function Conv(DomainType::Type, dim_in::Tuple{Int}, h::H) where {H<:AbstractVector} eltype(h) != DomainType && error("eltype(h) is $(eltype(h)), should be $(DomainType)") - buf = zeros(dim_in[1]+length(h)-1) - R = plan_rfft(buf) - buf_c1 = zeros(Complex{eltype(h)}, div(dim_in[1]+length(h)-1,2)+1) - buf_c2 = zeros(Complex{eltype(h)}, div(dim_in[1]+length(h)-1,2)+1) - I = plan_irfft(buf_c1,dim_in[1]+length(h)-1) + if isTypeReal(DomainType) + buf = zeros(DomainType,dim_in[1]+length(h)-1) + R = plan_rfft(buf) + buf_c1 = zeros(Complex{DomainType}, div(dim_in[1]+length(h)-1,2)+1) + buf_c2 = zeros(Complex{DomainType}, div(dim_in[1]+length(h)-1,2)+1) + I = plan_irfft(buf_c1,dim_in[1]+length(h)-1) + else + buf = zeros(DomainType,dim_in[1]+length(h)-1) + R = plan_fft(buf) + buf_c1 = zeros(DomainType, div(dim_in[1]+length(h)-1,2)+1) + buf_c2 = zeros(DomainType, div(dim_in[1]+length(h)-1,2)+1) + I = plan_ifft(buf_c1,dim_in[1]+length(h)-1) + end Conv{DomainType, H, typeof(buf_c1)}(dim_in,h,buf,buf_c1,buf_c2,R,I) end @@ -43,11 +53,11 @@ Conv(x::H, h::H) where {H} = Conv(eltype(x), size(x), h) function mul!(y::H, A::Conv{T,H}, b::H) where {T, H} #y .= conv(A.h,b) #naive implementation for i in eachindex(A.buf) - A.buf[i] = i <= length(A.h) ? A.h[i] : zero(T) + A.buf[i] = i <= length(A.h) ? A.h[i] : zero(T) end mul!(A.buf_c1, A.R, A.buf) for i in eachindex(A.buf) - A.buf[i] = i <= length(b) ? b[i] : zero(T) + A.buf[i] = i <= length(b) ? b[i] : zero(T) end mul!(A.buf_c2, A.R, A.buf) A.buf_c2 .*= A.buf_c1 @@ -59,11 +69,11 @@ function mul!(y::H, L::AdjointOperator{C}, b::H) where {T, H, C <: Conv{T,H}} #y .= xcorr(b,L.A.h)[size(L.A,1)[1]:end-length(L.A.h)+1] #naive implementation for i in eachindex(L.A.buf) ii = length(L.A.buf)-i+1 - L.A.buf[ii] = i <= length(L.A.h) ? L.A.h[i] : zero(T) + L.A.buf[ii] = i <= length(L.A.h) ? L.A.h[i] : zero(T) end mul!(L.A.buf_c1, L.A.R, L.A.buf) for i in eachindex(L.A.buf) - L.A.buf[i] = b[i] + L.A.buf[i] = b[i] end mul!(L.A.buf_c2, L.A.R, L.A.buf) L.A.buf_c2 .*= L.A.buf_c1 @@ -79,7 +89,7 @@ end domainType(L::Conv{T}) where {T} = T codomainType(L::Conv{T}) where {T} = T -#TODO find out a way to verify this, +#TODO find out a way to verify this, is_full_row_rank(L::Conv) = true is_full_column_rank(L::Conv) = true diff --git a/src/linearoperators/DCT.jl b/src/linearoperators/DCT.jl index 5b0a455..9b77412 100644 --- a/src/linearoperators/DCT.jl +++ b/src/linearoperators/DCT.jl @@ -9,14 +9,14 @@ abstract type CosineTransform{N,C,T1,T2} <: LinearOperator end `DCT(x::AbstractArray)` -Creates a `LinearOperator` which, when multiplied with an array `x::AbstractArray{N}`, returns the `N`-dimensional Inverse Discrete Cosine Transform of `x`. +Creates a `LinearOperator` which, when multiplied with an array `x::AbstractArray{N}`, returns the `N`-dimensional Inverse Discrete Cosine Transform of `x`. ```julia julia> DCT(Complex{Float64},(10,10)) -ℱc ℂ^(10, 10) -> ℂ^(10, 10) +ℱc ℂ^(10, 10) -> ℂ^(10, 10) julia> DCT(10,10) -ℱc ℝ^(10, 10) -> ℂ^(10, 10) +ℱc ℝ^(10, 10) -> ℂ^(10, 10) julia> A = DCT(ones(3)) ℱc ℝ^3 -> ℝ^3 @@ -46,14 +46,14 @@ end `IDCT(x::AbstractArray)` -Creates a `LinearOperator` which, when multiplied with an array `x::AbstractArray{N}`, returns the `N`-dimensional Discrete Cosine Transform of `x`. +Creates a `LinearOperator` which, when multiplied with an array `x::AbstractArray{N}`, returns the `N`-dimensional Discrete Cosine Transform of `x`. ```julia julia> IDCT(Complex{Float64},(10,10)) -ℱc⁻¹ ℂ^(10, 10) -> ℂ^(10, 10) +ℱc⁻¹ ℂ^(10, 10) -> ℂ^(10, 10) julia> IDCT(10,10) -ℱc⁻¹ ℝ^(10, 10) -> ℂ^(10, 10) +ℱc⁻¹ ℝ^(10, 10) -> ℂ^(10, 10) julia> A = IDCT(ones(3)) ℱc⁻¹ ℝ^3 -> ℝ^3 diff --git a/src/linearoperators/DiagOp.jl b/src/linearoperators/DiagOp.jl index 4b200c8..cb161d5 100644 --- a/src/linearoperators/DiagOp.jl +++ b/src/linearoperators/DiagOp.jl @@ -27,14 +27,14 @@ end # Constructors ###standard constructor Operator{N}(DomainType::Type, DomainDim::NTuple{N,Int}) -function DiagOp(DomainType::Type, DomainDim::NTuple{N,Int}, d::T) where {N, T <: AbstractArray} +function DiagOp(DomainType::Type, DomainDim::NTuple{N,Int}, d::T) where {N, T <: AbstractArray} size(d) != DomainDim && error("dimension of d must coincide with DomainDim") C = eltype(d) <: Complex ? complex(DomainType) : DomainType DiagOp{N, DomainType, C, T}(DomainDim, d) end ###standard constructor with Scalar -function DiagOp(DomainType::Type, DomainDim::NTuple{N,Int}, d::T) where {N, T <: Number} +function DiagOp(DomainType::Type, DomainDim::NTuple{N,Int}, d::T) where {N, T <: Number} C = eltype(d) <: Complex ? Complex{DomainType} : DomainType DiagOp{N, DomainType, C, T}(DomainDim, d) end diff --git a/src/linearoperators/Eye.jl b/src/linearoperators/Eye.jl index 65b1146..d238a59 100644 --- a/src/linearoperators/Eye.jl +++ b/src/linearoperators/Eye.jl @@ -26,7 +26,7 @@ end # Constructors ###standard constructor Operator{N}(DomainType::Type, DomainDim::NTuple{N,Int}) -Eye(DomainType::Type, DomainDim::NTuple{N,Int}) where {N} = Eye{DomainType,N}(DomainDim) +Eye(DomainType::Type, DomainDim::NTuple{N,Int}) where {N} = Eye{DomainType,N}(DomainDim) ### Eye(t::Type, dims::Vararg{Integer}) = Eye(t,dims) diff --git a/src/linearoperators/Filt.jl b/src/linearoperators/Filt.jl index 4455095..8515a86 100644 --- a/src/linearoperators/Filt.jl +++ b/src/linearoperators/Filt.jl @@ -5,7 +5,7 @@ export Filt `Filt(x::AbstractVector, b::AbstractVector, [a::AbstractVector,])` -Creates a `LinearOperator` which, when multiplied with an array `x::AbstractVector`, returns a vector `y` filtered by an IIR filter of coefficients `b` and `a`. If only `b` is provided a FIR is used to comute `y` instead. +Creates a `LinearOperator` which, when multiplied with an array `x::AbstractVector`, returns a vector `y` filtered by an IIR filter of coefficients `b` and `a`. If only `b` is provided a FIR is used to comute `y` instead. """ struct Filt{T,N} <: LinearOperator @@ -151,7 +151,7 @@ size(L::Filt) = L.dim_in, L.dim_in fun_name(L::Filt) = size(L.a,1) != 1 ? "IIR" : "FIR" -#TODO find out a way to verify this, +#TODO find out a way to verify this, # probably for IIR it means zeros inside unit circle is_full_row_rank(L::Filt) = true is_full_column_rank(L::Filt) = true diff --git a/src/linearoperators/FiniteDiff.jl b/src/linearoperators/FiniteDiff.jl index 313c6f8..d2f6cd5 100644 --- a/src/linearoperators/FiniteDiff.jl +++ b/src/linearoperators/FiniteDiff.jl @@ -7,7 +7,7 @@ export FiniteDiff `FiniteDiff(x::AbstractArray, direction = 1)` -Creates a `LinearOperator` which, when multiplied with an array `x::AbstractArray{N}`, returns the discretized gradient over the specified `direction` obtained using forward finite differences. +Creates a `LinearOperator` which, when multiplied with an array `x::AbstractArray{N}`, returns the discretized gradient over the specified `direction` obtained using forward finite differences. ```julia julia> FiniteDiff(Float64,(3,)) @@ -22,56 +22,48 @@ true ``` """ -struct FiniteDiff{T,N,D,C <: CartesianIndices{N}} <: LinearOperator +struct FiniteDiff{T,N,D} <: LinearOperator dim_in::NTuple{N,Int} - idx::C function FiniteDiff{T,N,D}(dim_in) where {T,N,D} D > N && error("direction is bigger the number of dimension $N") - idx = CartesianIndices(([i == D ? (2:d) : (1:d) for (i,d) in enumerate(dim_in)]...,)) - new{T,N,D,typeof(idx)}(dim_in,idx) + new{T,N,D}(dim_in) end end # Constructors #default constructor FiniteDiff(domainType::Type, dim_in::NTuple{N,Int}, dir::Int64 = 1) where {N} = -FiniteDiff{domainType,N,dir}(dim_in) + FiniteDiff{domainType,N,dir}(dim_in) FiniteDiff(dim_in::NTuple{N,Int}, dir::Int64 = 1) where {N} = -FiniteDiff(Float64, dim_in, dir) + FiniteDiff(Float64, dim_in, dir) FiniteDiff(x::AbstractArray{T,N}, dir::Int64 = 1) where {T,N} = FiniteDiff(eltype(x), size(x), dir) # Mappings -@generated function mul!(y::AbstractArray{T,N}, - L::FiniteDiff{T,N,D}, - b::AbstractArray{T,N}) where {T,N,D} - z = zeros(Int,N) - z[D] = 1 - idx = CartesianIndex(z...) - ex = quote - for I in L.idx - y[I-$idx] = b[I]-b[I-$idx] - end - return y - end +function mul!(y::AbstractArray{T,N}, + L::FiniteDiff{T,N,D}, + b::AbstractArray{T,N}) where {T,N,D} + idx_1 = CartesianIndices(([i == D ? (2:d) : (1:d) for (i,d) in enumerate(L.dim_in)]...,)) + idx_2 = CartesianIndices(([i == D ? (1:d-1) : (1:d) for (i,d) in enumerate(L.dim_in)]...,)) + y .= b[idx_1] .- b[idx_2] + return y end -@generated function mul!(y::AbstractArray{T,N}, - L::AdjointOperator{FiniteDiff{T,N,D,C}}, - b::AbstractArray{T,N}) where {T,N,D,C} - z = zeros(Int,N) - z[D] = 1 - idx = CartesianIndex(z...) - ex = quote - for I in CartesianIndices(size(y)) - y[I] = - I[$D] == 1 ? -b[I] : - I[$D] == size(y,$D) ? b[I-$idx] : -b[I]+b[I-$idx] - end - return y - end +function mul!(y::AbstractArray{T,N}, + L::AdjointOperator{FiniteDiff{T,N,D}}, + b::AbstractArray{T,N}) where {T,N,D} + dim_in = L.A.dim_in + idx_start = CartesianIndices(([i == D ? 1 : (1:d) for (i,d) in enumerate(dim_in)]...,)) + idx_between_1 = CartesianIndices(([i == D ? (1:d-2) : (1:d) for (i,d) in enumerate(dim_in)]...,)) + idx_between_2 = CartesianIndices(([i == D ? (2:d-1) : (1:d) for (i,d) in enumerate(dim_in)]...,)) + idx_end_1 = CartesianIndices(([i == D ? (d-1:d-1) : (1:d) for (i,d) in enumerate(dim_in)]...,)) + idx_end_2 = CartesianIndices(([i == D ? (d:d) : (1:d) for (i,d) in enumerate(dim_in)]...,)) + y[idx_start] .= -b[idx_start] + y[idx_between_2] .= b[idx_between_1] .- b[idx_between_2] + y[idx_end_2] .= b[idx_end_1] + return y end # Properties @@ -79,7 +71,7 @@ end domainType(L::FiniteDiff{T, N}) where {T, N} = T codomainType(L::FiniteDiff{T, N}) where {T, N} = T -function size(L::FiniteDiff{T,N,D}) where {T,N,D} +function size(L::FiniteDiff{T,N,D}) where {T,N,D} dim_out = [L.dim_in...] dim_out[D] = dim_out[D]-1 return ((dim_out...,), L.dim_in) @@ -92,5 +84,3 @@ fun_name(L::FiniteDiff{T,N,D}) where {T,N,D} = "δx$D" is_full_row_rank(L::FiniteDiff) = true - - diff --git a/src/linearoperators/GetIndex.jl b/src/linearoperators/GetIndex.jl index f111037..d3653db 100644 --- a/src/linearoperators/GetIndex.jl +++ b/src/linearoperators/GetIndex.jl @@ -11,7 +11,7 @@ Creates a `LinearOperator` which, when multiplied with `x`, returns `x[idx]`. julia> x = collect(linspace(1,10,10)); julia> G = GetIndex(Float64,(10,), 1:3) -↓ ℝ^10 -> ℝ^3 +↓ ℝ^10 -> ℝ^3 julia> G*x 3-element Array{Float64,1}: @@ -33,7 +33,7 @@ struct GetIndex{N,M,T<:Tuple} <: LinearOperator end # Constructors -# default +# default function GetIndex(domainType::Type,dim_in::NTuple{M,Int},idx::T) where {M,T<:Tuple} length(idx) > M && error("cannot slice object of dimension $dim_in with $idx") dim_out = get_dim_out(dim_in,idx...) @@ -51,13 +51,13 @@ GetIndex(x::AbstractArray, idx::Tuple) = GetIndex(eltype(x), size(x), idx) # Mappings -function mul!(y::Array{T1,N},L::GetIndex{N,M,T2},b::Array{T1,M}) where {T1,N,M,T2} +function mul!(y::AbstractArray{T1,N}, L::GetIndex{N,M,T2}, b::AbstractArray{T1,M}) where {T1,N,M,T2} y .= view(b,L.idx...) end -function mul!(y::Array{T1,M},L::AdjointOperator{GetIndex{N,M,T2}},b::AbstractArray{T1,N}) where {T1,N,M,T2} - fill!(y,0.) - setindex!(y,b,L.A.idx...) +function mul!(y::AbstractArray{T1,M}, L::AdjointOperator{GetIndex{N,M,T2}}, b::AbstractArray{T1,N}) where {T1,N,M,T2} + fill!(y, 0.) + setindex!(y, b, L.A.idx...) end # Properties @@ -82,7 +82,7 @@ function get_dim_out(dim, args...) if length(args) != 1 dim2 = () for i in eachindex(args) - if args[i] != Colon() + if args[i] != Colon() !(typeof(args[i]) <: Int) && ( dim2 = (dim2..., length(args[i]) ) ) else dim2 = (dim2..., dim[i]) diff --git a/src/linearoperators/IRDFT.jl b/src/linearoperators/IRDFT.jl index d5e0f98..8d1d007 100644 --- a/src/linearoperators/IRDFT.jl +++ b/src/linearoperators/IRDFT.jl @@ -9,7 +9,7 @@ Creates a `LinearOperator` which, when multiplied with a complex array `x`, retu ```julia julia> A = IRDFT(Complex{Float64},(10,),19) -ℱ⁻¹ ℂ^10 -> ℝ^19 +ℱ⁻¹ ℂ^10 -> ℝ^19 julia> A = IRDFT((5,10,8),19,2) ℱ⁻¹ ℂ^(5, 10, 8) -> ℝ^(5, 19, 8) @@ -23,7 +23,7 @@ struct IRDFT{T <:Number, T1<:AbstractFFTs.Plan, T2<:AbstractFFTs.Plan, T3<:NTuple{N,Any} - } <: LinearOperator + } <: LinearOperator dim_in::NTuple{N,Int} dim_out::NTuple{N,Int} A::T1 @@ -35,7 +35,7 @@ end #standard constructor function IRDFT(x::AbstractArray{Complex{T},N}, d::Int, dims::Int=1) where {T<:Number,N} - A = plan_irfft(x,d,dims) + A = plan_irfft(x,d,dims) dim_in = size(x) dim_out = () idx = () @@ -43,7 +43,7 @@ function IRDFT(x::AbstractArray{Complex{T},N}, d::Int, dims::Int=1) where {T<:Nu dim_out = i == dims ? (dim_out..., d) : (dim_out...,dim_in[i] ) idx = i == dims ? (idx... , 2:ceil(Int,d/2)) : (idx... ,Colon() ) end - At = plan_rfft(zeros(dim_out),dims) + At = plan_rfft(similar(x, T, dim_out),dims) IRDFT{T,N,dims,typeof(A),typeof(At),typeof(idx)}(dim_in,dim_out,A,At,idx) end @@ -63,7 +63,7 @@ function mul!(y::C2, L::AdjointOperator{IRDFT{T,N,D,T1,T2,T3}}, b::C1) where {N,T,D,T1,T2,T3,C1<:AbstractArray{T,N}, C2<:AbstractArray{Complex{T},N}} - + A = L.A mul!(y,A.At,b) y ./= size(b,D) diff --git a/src/linearoperators/LBFGS.jl b/src/linearoperators/LBFGS.jl index 4b6064a..2da4647 100644 --- a/src/linearoperators/LBFGS.jl +++ b/src/linearoperators/LBFGS.jl @@ -27,10 +27,10 @@ mutable struct LBFGS{R, T <: AbstractArray, M, I <: Integer} <: LinearOperator curridx::I s::T y::T - s_M::Array{T, 1} - y_M::Array{T, 1} - ys_M::Array{R, 1} - alphas::Array{R, 1} + s_M::AbstractArray{T, 1} + y_M::AbstractArray{T, 1} + ys_M::AbstractArray{R, 1} + alphas::AbstractArray{R, 1} H::R end @@ -42,8 +42,10 @@ function LBFGS(x::T, M::I) where {T <: AbstractArray, I <: Integer} s = zero(x) y = zero(x) R = real(eltype(x)) - ys_M = zeros(R, M) - alphas = zeros(R, M) + ys_M = similar(x, R, M) + fill!(ys_M, one(R)) + alphas = similar(x, R, M) + fill!(alphas, one(R)) LBFGS{R, T, M, I}(0, 0, s, y, s_M, y_M, ys_M, alphas, one(R)) end diff --git a/src/linearoperators/LMatrixOp.jl b/src/linearoperators/LMatrixOp.jl index 91d0f28..3ee58a6 100644 --- a/src/linearoperators/LMatrixOp.jl +++ b/src/linearoperators/LMatrixOp.jl @@ -10,7 +10,7 @@ Creates a `LinearOperator` which, when multiplied with a matrix `X::AbstractMatr ```julia julia> op = LMatrixOp(Float64,(3,4),ones(4)) -(⋅)b ℝ^(3, 4) -> ℝ^3 +(⋅)b ℝ^(3, 4) -> ℝ^3 julia> op = LMatrixOp(ones(4),3) (⋅)b ℝ^(3, 4) -> ℝ^3 @@ -24,7 +24,7 @@ julia> op*ones(3,4) ``` """ -struct LMatrixOp{T, A <: Union{AbstractVector,AbstractMatrix}, +struct LMatrixOp{T, A <: Union{AbstractVector,AbstractMatrix}, B <: AbstractMatrix} <: LinearOperator b::A bt::B @@ -34,23 +34,23 @@ end ##TODO decide what to do when domainType is given, with conversion one loses pointer to data... # Constructors function LMatrixOp(DomainType::Type, - DomainDim::Tuple{Int,Int}, b::A) where {A <: Union{AbstractVector,AbstractMatrix}} + DomainDim::Tuple{Int,Int}, b::A) where {A <: Union{AbstractVector,AbstractMatrix}} bt = b' LMatrixOp{DomainType, A, typeof(bt)}(b,bt,DomainDim[1]) end -LMatrixOp(b::A, n_row_in::Int) where {T,A<:Union{AbstractVector{T},AbstractMatrix{T}}} = -LMatrixOp(T,(n_row_in,size(b,1)),b) +LMatrixOp(b::A, n_row_in::Int) where {T,A<:Union{AbstractVector{T},AbstractMatrix{T}}} = +LMatrixOp(T,(n_row_in,size(b,1)),b) # Mappings -mul!(y::C, L::LMatrixOp{T,A,B}, X::AbstractMatrix{T} ) where {T,A,B,C<:Union{AbstractVector,AbstractMatrix}} = +mul!(y::C, L::LMatrixOp{T,A,B}, X::AbstractMatrix{T} ) where {T,A,B,C<:Union{AbstractVector,AbstractMatrix}} = mul!(y,X,L.b) -function mul!(y::AbstractMatrix{T}, L::AdjointOperator{LMatrixOp{T,A,B}}, Y::AbstractVector{T} ) where {T,A,B} +function mul!(y::AbstractMatrix{T}, L::AdjointOperator{LMatrixOp{T,A,B}}, Y::AbstractVector{T} ) where {T,A,B} y .= L.A.bt.*Y end -function mul!(y::AbstractMatrix{T}, L::AdjointOperator{LMatrixOp{T,A,B}}, Y::AbstractMatrix{T} ) where {T,A,B} +function mul!(y::AbstractMatrix{T}, L::AdjointOperator{LMatrixOp{T,A,B}}, Y::AbstractMatrix{T} ) where {T,A,B} mul!(y,Y,L.A.b') end @@ -65,5 +65,5 @@ size(L::LMatrixOp{T,A,B}) where {T,A <: AbstractMatrix,B <: AbstractMatrix} = (L #TODO -#is_full_row_rank(L::LMatrixOp) = +#is_full_row_rank(L::LMatrixOp) = #is_full_column_rank(L::MatrixOp) = diff --git a/src/linearoperators/MIMOFilt.jl b/src/linearoperators/MIMOFilt.jl index 5a00ba0..8f33462 100644 --- a/src/linearoperators/MIMOFilt.jl +++ b/src/linearoperators/MIMOFilt.jl @@ -5,14 +5,14 @@ export MIMOFilt `MIMOFilt(x::AbstractMatrix, b::Vector{AbstractVector}, [a::Vector{AbstractVector},])` -Creates a `LinearOperator` which, when multiplied with a matrix `X`, returns a matrix `Y`. Here a Multiple Input Multiple Output system is evaluated: the columns of `X` and `Y` represent the input signals and output signals respectively. +Creates a `LinearOperator` which, when multiplied with a matrix `X`, returns a matrix `Y`. Here a Multiple Input Multiple Output system is evaluated: the columns of `X` and `Y` represent the input signals and output signals respectively. ```math -\\mathbf{y}_i = \\sum_{j = 1}^{M} \\mathbf{h}_{i,j} * \\mathbf{x}_j +\\mathbf{y}_i = \\sum_{j = 1}^{M} \\mathbf{h}_{i,j} * \\mathbf{x}_j ``` where ``\\mathbf{y}_i`` and ``\\mathbf{x}_j`` are the ``i``-th and ``j``-th columns of the output `Y` and input `X` matrices respectively. -The filters ``\\mathbf{h}_{i,j}`` can be represented either by providing coefficients `B` and `A` (IIR) or `B` alone (FIR). These coefficients must be given in a `Vector` of `Vector`s. +The filters ``\\mathbf{h}_{i,j}`` can be represented either by providing coefficients `B` and `A` (IIR) or `B` alone (FIR). These coefficients must be given in a `Vector` of `Vector`s. For example for a `3` by `2` MIMO system (i.e. `size(X,2) == 3` inputs and `size(Y,2) == 2` outputs) `B` must be: @@ -30,7 +30,7 @@ julia> A = [[1.;1.;1.],[2.;2.;2.],[ 3.],[ 4.],[ 5.],[ 6.], #A = [ a11 , a12 , a13 , a21 , a22, a23 , ] julia> op = MIMOFilt(Float64, (m,n), B, A) -※ ℝ^(10, 3) -> ℝ^(10, 2) +※ ℝ^(10, 3) -> ℝ^(10, 2) julia> X = randn(m,n); #input signals @@ -160,7 +160,7 @@ codomainType(L::MIMOFilt{T, M}) where {T, M} = T size(L::MIMOFilt) = L.dim_out, L.dim_in -#TODO find out a way to verify this, +#TODO find out a way to verify this, # probably for IIR it means zeros inside unit circle is_full_row_rank(L::MIMOFilt) = true is_full_column_rank(L::MIMOFilt) = true diff --git a/src/linearoperators/MatrixOp.jl b/src/linearoperators/MatrixOp.jl index f8868ff..fba605b 100644 --- a/src/linearoperators/MatrixOp.jl +++ b/src/linearoperators/MatrixOp.jl @@ -13,7 +13,7 @@ The input `x` can be also a matrix: the number of columns must be given either i ```julia julia> MatrixOp(Float64,(10,),randn(20,10)) -▒ ℝ^10 -> ℝ^20 +▒ ℝ^10 -> ℝ^20 julia> MatrixOp(randn(20,10)) ▒ ℝ^10 -> ℝ^20 @@ -36,8 +36,8 @@ end ##TODO decide what to do when domainType is given, with conversion one loses pointer to data... ###standard constructor Operator{N}(DomainType::Type, DomainDim::NTuple{N,Int}) -function MatrixOp(DomainType::Type, DomainDim::NTuple{N,Int}, A::M) where {N, T, M <: AbstractMatrix{T}} - N > 2 && error("cannot multiply a Matrix by a n-dimensional Variable with n > 2") +function MatrixOp(DomainType::Type, DomainDim::NTuple{N,Int}, A::M) where {N, T, M <: AbstractMatrix{T}} + N > 2 && error("cannot multiply a Matrix by a n-dimensional Variable with n > 2") size(A,2) != DomainDim[1] && error("wrong input dimensions") if N == 1 MatrixOp{DomainType, T, M}(A, 1) @@ -63,8 +63,8 @@ mul!(y::AbstractArray, L::MatrixOp{D, T, M}, b::AbstractArray) where {D, T, M} = mul!(y::AbstractArray, L::AdjointOperator{MatrixOp{D, T, M}}, b::AbstractArray) where {D, T, M} = mul!(y, L.A.A', b) # Special Case, real b, complex matrix -function mul!(y::AbstractArray, L::AdjointOperator{MatrixOp{D, T, M}}, b::AbstractArray) where {D <: Real , T <: Complex, M} - yc = zeros(T,size(y)) +function mul!(y::AbstractArray, L::AdjointOperator{MatrixOp{D, T, M}}, b::AbstractArray) where {D <: Real , T <: Complex, M} + yc = similar(y,T,size(y)) mul!(yc, L.A.A', b) y .= real.(yc) end diff --git a/src/linearoperators/MyLinOp.jl b/src/linearoperators/MyLinOp.jl index 45a8851..7de0a94 100644 --- a/src/linearoperators/MyLinOp.jl +++ b/src/linearoperators/MyLinOp.jl @@ -38,8 +38,8 @@ MyLinOp{N,M, domainType, codomainType}(dim_out, dim_in, Fwd!, Adj! ) # Mappings -mul!(y::Array{C,N}, L::MyLinOp{N,M,C,D}, b::Array{D,M}) where {N,M,C,D} = L.Fwd!(y,b) -mul!(y::Array{C,N}, L::AdjointOperator{MyLinOp{N,M,C,D}}, b::Array{D,M}) where {N,M,C,D} = L.A.Adj!(y,b) +mul!(y::AbstractArray{C,N}, L::MyLinOp{N,M,C,D}, b::AbstractArray{D,M}) where {N,M,C,D} = L.Fwd!(y,b) +mul!(y::AbstractArray{C,N}, L::AdjointOperator{MyLinOp{N,M,C,D}}, b::AbstractArray{D,M}) where {N,M,C,D} = L.A.Adj!(y,b) # Properties diff --git a/src/linearoperators/RDFT.jl b/src/linearoperators/RDFT.jl index 97967e1..7d2daa8 100644 --- a/src/linearoperators/RDFT.jl +++ b/src/linearoperators/RDFT.jl @@ -7,7 +7,7 @@ export RDFT `RDFT(x::AbstractArray [,dims=1])` -Creates a `LinearOperator` which, when multiplied with a real array `x`, returns the DFT over the dimension `dims`, exploiting Hermitian symmetry. +Creates a `LinearOperator` which, when multiplied with a real array `x`, returns the DFT over the dimension `dims`, exploiting Hermitian symmetry. ```julia julia> RDFT(Float64,(10,10)) @@ -24,7 +24,7 @@ struct RDFT{T <:Number, T1<:AbstractFFTs.Plan, T2<:AbstractFFTs.Plan, T3<:AbstractArray{Complex{T},N} - } <: LinearOperator + } <: LinearOperator dim_in::NTuple{N,Int} dim_out::NTuple{N,Int} A::T1 @@ -39,9 +39,9 @@ end #standard constructor function RDFT(x::AbstractArray{T,N}, dims::Int=1) where {T<:Real,N} - A = plan_rfft(x,dims) - b2 = zeros(complex(T),size(x)) - y2 = zeros(complex(T),size(x)) + A = plan_rfft(x,dims) + b2 = similar(x,complex(T),size(x)) + y2 = similar(x,complex(T),size(x)) At = plan_bfft(y2,dims) dim_in = size(x) dim_out = () diff --git a/src/linearoperators/Variation.jl b/src/linearoperators/Variation.jl index 68de91d..fa3670c 100644 --- a/src/linearoperators/Variation.jl +++ b/src/linearoperators/Variation.jl @@ -7,7 +7,7 @@ export Variation `Variation(x::AbstractArray)` -Creates a `LinearOperator` which, when multiplied with an array `x::AbstractArray{N}`, returns a matrix with its `i`th column consisting of the vectorized discretized gradient over the `i`th `direction obtained using forward finite differences. +Creates a `LinearOperator` which, when multiplied with an array `x::AbstractArray{N}`, returns a matrix with its `i`th column consisting of the vectorized discretized gradient over the `i`th `direction obtained using forward finite differences. ```julia julia> Variation(Float64,(10,2)) @@ -32,7 +32,7 @@ end # Constructors #default constructor -function Variation(domainType::Type, dim_in::NTuple{N,Int}) where {N} +function Variation(domainType::Type, dim_in::NTuple{N,Int}) where {N} N == 1 && error("use FiniteDiff instead!") Variation{domainType,N}(dim_in) end @@ -44,20 +44,19 @@ Variation(x::AbstractArray) = Variation(eltype(x), size(x)) # Mappings -@generated function mul!(y::AbstractArray{T,2}, +@generated function mul!(y::AbstractArray{T,2}, A::Variation{T,N}, b::AbstractArray{T,N}) where {T,N} - ex = :() for i = 1:N z = zeros(Int,N) z[i] = 1 z = (z...,) - ex = :($ex; y[cnt,$i] = I[$i] == 1 ? b[I+CartesianIndex($z)]-b[I] : + ex = :($ex; y[cnt,$i] = I[$i] == 1 ? b[I+CartesianIndex($z)]-b[I] : b[I]-b[I-CartesianIndex($z)]) end - ex2 = quote + ex2 = quote cnt = 0 for I in CartesianIndices(size(b)) cnt += 1 @@ -67,7 +66,7 @@ Variation(x::AbstractArray) = Variation(eltype(x), size(x)) end end -@generated function mul!(y::AbstractArray{T,N}, +@generated function mul!(y::AbstractArray{T,N}, A::AdjointOperator{Variation{T,N}}, b::AbstractArray{T,2}) where {T,N} ex = :(y[I] = I[1] == 1 ? -(b[cnt,1] + b[cnt+1,1]) : @@ -77,8 +76,8 @@ end Nx = :(size(y,1)) for i = 2:N - ex = quote - $ex + ex = quote + $ex y[I] += I[$i] == 1 ? -(b[cnt,$i] + b[cnt+$Nx,$i]) : I[$i] == 2 ? b[cnt,$i] + b[cnt-$Nx,$i] - b[cnt+$Nx,$i] : I[$i] == size(y,$i) ? b[cnt,$i] : b[cnt, $i] - b[cnt+$Nx,$i] @@ -86,7 +85,7 @@ end Nx = :($Nx*size(y,$i)) end - ex2 = quote + ex2 = quote cnt = 0 for I in CartesianIndices(size(y)) cnt += 1 diff --git a/src/linearoperators/Xcorr.jl b/src/linearoperators/Xcorr.jl index 24691bd..fb54e2f 100644 --- a/src/linearoperators/Xcorr.jl +++ b/src/linearoperators/Xcorr.jl @@ -6,7 +6,7 @@ export Xcorr `Xcorr(x::AbstractVector, h::AbstractVector)` -Creates a `LinearOperator` which, when multiplied with an array `x::AbstractVector`, returns the cross correlation between `x` and `h`. Uses `xcross`. +Creates a `LinearOperator` which, when multiplied with an array `x::AbstractVector`, returns the cross correlation between `x` and `h`. Uses `xcross`. """ struct Xcorr{T,H <:AbstractVector{T}} <: LinearOperator @@ -15,7 +15,7 @@ struct Xcorr{T,H <:AbstractVector{T}} <: LinearOperator end # Constructors -function Xcorr(DomainType::Type, DomainDim::NTuple{N,Int}, h::H) where {H<:AbstractVector, N} +function Xcorr(DomainType::Type, DomainDim::NTuple{N,Int}, h::H) where {H<:AbstractVector, N} eltype(h) != DomainType && error("eltype(h) is $(eltype(h)), should be $(DomainType)") N != 1 && error("Xcorr treats only SISO, check Filt and MIMOFilt for MIMO") Xcorr{DomainType,H}(DomainDim,h) @@ -40,7 +40,7 @@ end domainType(L::Xcorr{T}) where {T} = T codomainType(L::Xcorr{T}) where {T} = T -#TODO find out a way to verify this, +#TODO find out a way to verify this, is_full_row_rank(L::Xcorr) = true is_full_column_rank(L::Xcorr) = true diff --git a/src/linearoperators/ZeroPad.jl b/src/linearoperators/ZeroPad.jl index 6bdd992..751a0ea 100644 --- a/src/linearoperators/ZeroPad.jl +++ b/src/linearoperators/ZeroPad.jl @@ -6,7 +6,7 @@ export ZeroPad `ZeroPad(x::AbstractArray, zp::Tuple)` -Create a `LinearOperator` which, when multiplied to an array `x` of size `dim_in`, returns an expanded array `y` of size `dim_in .+ zp` where `y[1:dim_in[1], 1:dim_in[2] ... ] = x` and zero elsewhere. +Create a `LinearOperator` which, when multiplied to an array `x` of size `dim_in`, returns an expanded array `y` of size `dim_in .+ zp` where `y[1:dim_in[1], 1:dim_in[2] ... ] = x` and zero elsewhere. ```julia julia> Z = ZeroPad((2,2),(0,2)) @@ -44,8 +44,8 @@ ZeroPad(x::AbstractArray, zp::Vararg{Int,N}) where {N} = ZeroPad( # builds - #for i1 =1:size(y,1), i2 =1:size(y,2) - # y[i1,i2] = i1 <= size(b,1) && i2 <= size(b,2) ? b[i1,i2] : 0. + #for i1 =1:size(y,1), i2 =1:size(y,2) + # y[i1,i2] = i1 <= size(b,1) && i2 <= size(b,2) ? b[i1,i2] : 0. #end ex = "for " @@ -53,21 +53,21 @@ ZeroPad(x::AbstractArray, zp::Vararg{Int,N}) where {N} = ZeroPad( ex = ex[1:end-1] #remove , - ex *= " y[" + ex *= " y[" for i = 1:N ex *= "i$i," end ex = ex[1:end-1] #remove , - ex *= "] = " + ex *= "] = " for i = 1:N ex *= " i$i <= size(b,$i) &&" end ex = ex[1:end-2] #remove && - ex *= " ? b[" + ex *= " ? b[" for i = 1:N ex *= "i$i," end ex = ex[1:end-1] #remove , - ex *= "] : 0. end" + ex *= "] : 0. end" ex = Meta.parse(ex) @@ -83,17 +83,17 @@ end ex = "for " for i = 1:N ex *= "i$i =1:size(y,$i)," end - ex *= " y[" + ex *= " y[" idx = "" for i = 1:N idx *= "i$i," end idx = idx[1:end-1] #remove , - ex *= idx + ex *= idx - ex *= "] = b[" - ex *= idx - ex *= "] end" + ex *= "] = b[" + ex *= idx + ex *= "] end" ex = Meta.parse(ex) diff --git a/src/linearoperators/Zeros.jl b/src/linearoperators/Zeros.jl index 382468a..4dba105 100644 --- a/src/linearoperators/Zeros.jl +++ b/src/linearoperators/Zeros.jl @@ -22,21 +22,21 @@ struct Zeros{C,N,D,M} <: LinearOperator end # Constructors -#default -Zeros(domainType::Type, dim_in::NTuple{M,Int}, +#default +Zeros(domainType::Type, dim_in::NTuple{M,Int}, codomainType::Type, dim_out::NTuple{N,Int}) where {N,M} = Zeros{codomainType,N,domainType,M}(dim_out,dim_in) -Zeros(domainType::Type, dim_in::NTuple{M,Int}, dim_out::NTuple{N,Int}) where {N,M} = +Zeros(domainType::Type, dim_in::NTuple{M,Int}, dim_out::NTuple{N,Int}) where {N,M} = Zeros{domainType,N,domainType,M}(dim_out,dim_in) -function Zeros(domainType::NTuple{NN,Type}, dim_in::NTuple{NN,Tuple}, +function Zeros(domainType::NTuple{NN,Type}, dim_in::NTuple{NN,Tuple}, codomainType::Type, dim_out::Tuple) where {NN} HCAT([Zeros(domainType[i], dim_in[i], codomainType, dim_out) for i =1:NN]...) end -function Zeros(domainType::Type, dim_in::Tuple, - codomainType::NTuple{NN,Type}, dim_out::NTuple{NN,Tuple}) where {NN} +function Zeros(domainType::Type, dim_in::Tuple, + codomainType::NTuple{NN,Type}, dim_out::NTuple{NN,Tuple}) where {NN} VCAT([Zeros(domainType, dim_in, codomainType[i], dim_out[i]) for i =1:NN]...) end diff --git a/src/nonlinearoperators/Atan.jl b/src/nonlinearoperators/Atan.jl index f1669db..783bca4 100644 --- a/src/nonlinearoperators/Atan.jl +++ b/src/nonlinearoperators/Atan.jl @@ -13,7 +13,7 @@ struct Atan{T,N} <: NonLinearOperator dim::NTuple{N,Int} end -function Atan(DomainType::Type, DomainDim::NTuple{N,Int}) where {N} +function Atan(DomainType::Type, DomainDim::NTuple{N,Int}) where {N} Atan{DomainType,N}(DomainDim) end @@ -24,8 +24,8 @@ function mul!(y::AbstractArray{T,N}, L::Atan{T,N}, x::AbstractArray{T,N}) where y .= atan.(x) end -function mul!(y::AbstractArray, - J::AdjointOperator{Jacobian{A,TT}}, +function mul!(y::AbstractArray, + J::AdjointOperator{Jacobian{A,TT}}, b::AbstractArray) where {T, N, A <: Atan{T,N}, TT <: AbstractArray{T,N}} L = J.A y .= conj.(1.0./(1.0.+ L.x.^2)).*b diff --git a/src/nonlinearoperators/Cos.jl b/src/nonlinearoperators/Cos.jl index f5fd9cf..ccf2140 100644 --- a/src/nonlinearoperators/Cos.jl +++ b/src/nonlinearoperators/Cos.jl @@ -13,7 +13,7 @@ struct Cos{T,N} <: NonLinearOperator dim::NTuple{N,Int} end -function Cos(DomainType::Type, DomainDim::NTuple{N,Int}) where {N} +function Cos(DomainType::Type, DomainDim::NTuple{N,Int}) where {N} Cos{DomainType,N}(DomainDim) end @@ -24,8 +24,8 @@ function mul!(y::AbstractArray{T,N}, L::Cos{T,N}, x::AbstractArray{T,N}) where { y .= cos.(x) end -function mul!(y::AbstractArray, - J::AdjointOperator{Jacobian{A,TT}}, +function mul!(y::AbstractArray, + J::AdjointOperator{Jacobian{A,TT}}, b::AbstractArray) where {T,N, A<: Cos{T,N}, TT <: AbstractArray{T,N}} L = J.A y .= -conj.(sin.(L.x)).*b diff --git a/src/nonlinearoperators/Exp.jl b/src/nonlinearoperators/Exp.jl index e85fa7b..3270d33 100644 --- a/src/nonlinearoperators/Exp.jl +++ b/src/nonlinearoperators/Exp.jl @@ -13,7 +13,7 @@ struct Exp{T,N} <: NonLinearOperator dim::NTuple{N,Int} end -function Exp(DomainType::Type, DomainDim::NTuple{N,Int}) where {N} +function Exp(DomainType::Type, DomainDim::NTuple{N,Int}) where {N} Exp{DomainType,N}(DomainDim) end @@ -24,8 +24,8 @@ function mul!(y::AbstractArray{T,N}, L::Exp{T,N}, x::AbstractArray{T,N}) where { y .= exp.(x) end -function mul!(y::AbstractArray, - J::AdjointOperator{Jacobian{A,TT}}, +function mul!(y::AbstractArray, + J::AdjointOperator{Jacobian{A,TT}}, b::AbstractArray) where {T,N, A<: Exp{T,N}, TT <: AbstractArray{T,N} } L = J.A y .= conj.(exp.(L.x)).*b diff --git a/src/nonlinearoperators/Pow.jl b/src/nonlinearoperators/Pow.jl index ea05751..6ad187d 100644 --- a/src/nonlinearoperators/Pow.jl +++ b/src/nonlinearoperators/Pow.jl @@ -11,7 +11,7 @@ struct Pow{T,N,I<:Real} <: NonLinearOperator p::I end -function Pow(DomainType::Type, DomainDim::NTuple{N,Int}, p::I) where {N, I <: Real} +function Pow(DomainType::Type, DomainDim::NTuple{N,Int}, p::I) where {N, I <: Real} Pow{DomainType, N, I}(DomainDim, p) end @@ -21,8 +21,8 @@ function mul!(y::AbstractArray{T,N}, L::Pow{T,N,I}, x::AbstractArray{T,N}) where y .= x.^L.p end -function mul!(y::AbstractArray, - J::AdjointOperator{Jacobian{Pow{T, N, I},TT}}, +function mul!(y::AbstractArray, + J::AdjointOperator{Jacobian{Pow{T, N, I},TT}}, b::AbstractArray) where {T, N, I, TT <: AbstractArray{T,N}} L = J.A y .= conj.(L.A.p.*(L.x).^(L.A.p-1)).*b diff --git a/src/nonlinearoperators/Sech.jl b/src/nonlinearoperators/Sech.jl index 2556ece..35929c6 100644 --- a/src/nonlinearoperators/Sech.jl +++ b/src/nonlinearoperators/Sech.jl @@ -13,7 +13,7 @@ struct Sech{T,N} <: NonLinearOperator dim::NTuple{N,Int} end -function Sech(DomainType::Type, DomainDim::NTuple{N,Int}) where {N} +function Sech(DomainType::Type, DomainDim::NTuple{N,Int}) where {N} Sech{DomainType,N}(DomainDim) end @@ -24,8 +24,8 @@ function mul!(y::AbstractArray{T,N}, L::Sech{T,N}, x::AbstractArray{T,N}) where y .= sech.(x) end -function mul!(y::AbstractArray, - J::AdjointOperator{Jacobian{A,TT}}, +function mul!(y::AbstractArray, + J::AdjointOperator{Jacobian{A,TT}}, b::AbstractArray) where {T,N, A<: Sech{T,N}, TT <: AbstractArray{T,N}} L = J.A y .= -conj.( tanh.(L.x) .* sech.(L.x) ).*b diff --git a/src/nonlinearoperators/Sigmoid.jl b/src/nonlinearoperators/Sigmoid.jl index 1d5beec..a1e823e 100644 --- a/src/nonlinearoperators/Sigmoid.jl +++ b/src/nonlinearoperators/Sigmoid.jl @@ -14,7 +14,7 @@ struct Sigmoid{T,N,G<:Real} <: NonLinearOperator gamma::G end -function Sigmoid(DomainType::Type, DomainDim::NTuple{N,Int}, gamma::G=1.) where {N, G <: Real} +function Sigmoid(DomainType::Type, DomainDim::NTuple{N,Int}, gamma::G=1.) where {N, G <: Real} Sigmoid{DomainType,N,G}(DomainDim,gamma) end @@ -25,12 +25,12 @@ function mul!(y::AbstractArray{T,N}, L::Sigmoid{T,N,G}, x::AbstractArray{T,N}) w end -function mul!(y::AbstractArray, - J::AdjointOperator{Jacobian{A,TT}}, +function mul!(y::AbstractArray, + J::AdjointOperator{Jacobian{A,TT}}, b::AbstractArray) where {T,N,G, A<: Sigmoid{T,N,G}, TT <: AbstractArray{T,N}} L = J.A y .= exp.(-L.A.gamma.*L.x) - y ./= (1 .+y).^2 + y ./= (1 .+y).^2 y .= conj.(L.A.gamma.*y) y .*= b end diff --git a/src/nonlinearoperators/Sin.jl b/src/nonlinearoperators/Sin.jl index 3c79628..34fd74c 100644 --- a/src/nonlinearoperators/Sin.jl +++ b/src/nonlinearoperators/Sin.jl @@ -13,7 +13,7 @@ struct Sin{T,N} <: NonLinearOperator dim::NTuple{N,Int} end -function Sin(DomainType::Type, DomainDim::NTuple{N,Int}) where {N} +function Sin(DomainType::Type, DomainDim::NTuple{N,Int}) where {N} Sin{DomainType,N}(DomainDim) end @@ -24,8 +24,8 @@ function mul!(y::AbstractArray{T,N}, L::Sin{T,N}, x::AbstractArray{T,N}) where { y .= sin.(x) end -function mul!(y::AbstractArray, - J::AdjointOperator{Jacobian{A,TT}}, +function mul!(y::AbstractArray, + J::AdjointOperator{Jacobian{A,TT}}, b::AbstractArray) where {T,N, A<: Sin{T,N}, TT <: AbstractArray{T,N}} L = J.A y .= conj.(cos.(L.x)).*b diff --git a/src/nonlinearoperators/SoftMax.jl b/src/nonlinearoperators/SoftMax.jl index f26408d..1057948 100644 --- a/src/nonlinearoperators/SoftMax.jl +++ b/src/nonlinearoperators/SoftMax.jl @@ -11,10 +11,14 @@ Creates the softmax non-linear operator with input dimensions `dim_in`. """ struct SoftMax{T,N} <: NonLinearOperator dim::NTuple{N,Int} - buf::Array{T,N} + buf::AbstractArray{T,N} end -function SoftMax(DomainType::Type, DomainDim::NTuple{N,Int}) where {N} +function SoftMax(x::AbstractArray{T,N}) where {T,N} + SoftMax{N,T}(size(x),similar(x)) +end + +function SoftMax(DomainType::Type, DomainDim::NTuple{N,Int}) where {N} SoftMax{DomainType,N}(DomainDim,zeros(DomainType,DomainDim)) end @@ -25,15 +29,15 @@ function mul!(y::AbstractArray{T,N}, L::SoftMax{T,N}, x::AbstractArray{T,N}) whe y ./= sum(y) end -function mul!(y::AbstractArray, - J::AdjointOperator{Jacobian{A,TT}}, +function mul!(y::AbstractArray, + J::AdjointOperator{Jacobian{A,TT}}, b::AbstractArray) where {T, N, A<: SoftMax{T,N}, TT <: AbstractArray{T,N} } L = J.A fill!(y,zero(T)) L.A.buf .= exp.(L.x.-maximum(L.x)) L.A.buf ./= sum(L.A.buf) for i in eachindex(y) - y[i] = -L.A.buf[i]*dot(L.A.buf,b) + y[i] = -L.A.buf[i]*dot(L.A.buf,b) y[i] += L.A.buf[i]*b[i] end return y diff --git a/src/nonlinearoperators/SoftPlus.jl b/src/nonlinearoperators/SoftPlus.jl index 0ee54b1..757996b 100644 --- a/src/nonlinearoperators/SoftPlus.jl +++ b/src/nonlinearoperators/SoftPlus.jl @@ -13,7 +13,7 @@ struct SoftPlus{T,N} <: NonLinearOperator dim::NTuple{N,Int} end -function SoftPlus(DomainType::Type, DomainDim::NTuple{N,Int}) where {N} +function SoftPlus(DomainType::Type, DomainDim::NTuple{N,Int}) where {N} SoftPlus{DomainType,N}(DomainDim) end @@ -23,8 +23,8 @@ function mul!(y::AbstractArray{T,N}, L::SoftPlus{T,N}, x::AbstractArray{T,N}) wh y .= log.(1 .+exp.(x)) end -function mul!(y::AbstractArray, - J::AdjointOperator{Jacobian{A,TT}}, +function mul!(y::AbstractArray, + J::AdjointOperator{Jacobian{A,TT}}, b::AbstractArray) where {T, N, A <: SoftPlus{T,N}, TT <: AbstractArray{T,N} } L = J.A y .= 1 ./(1 .+exp.(-L.x)).*b diff --git a/src/nonlinearoperators/Tanh.jl b/src/nonlinearoperators/Tanh.jl index ea5e99e..ba8e45d 100644 --- a/src/nonlinearoperators/Tanh.jl +++ b/src/nonlinearoperators/Tanh.jl @@ -13,7 +13,7 @@ struct Tanh{T,N} <: NonLinearOperator dim::NTuple{N,Int} end -function Tanh(DomainType::Type, DomainDim::NTuple{N,Int}) where {N} +function Tanh(DomainType::Type, DomainDim::NTuple{N,Int}) where {N} Tanh{DomainType,N}(DomainDim) end @@ -24,8 +24,8 @@ function mul!(y::AbstractArray{T,N}, L::Tanh{T,N}, x::AbstractArray{T,N}) where y .= tanh.(x) end -function mul!(y::AbstractArray, - J::AdjointOperator{Jacobian{A,TT}}, +function mul!(y::AbstractArray, + J::AdjointOperator{Jacobian{A,TT}}, b::AbstractArray) where {T,N, A<: Tanh{T,N}, TT <: AbstractArray{T,N}} L = J.A y .= conj.(sech.(L.x).^2).*b diff --git a/src/properties.jl b/src/properties.jl index 03e6740..f7e5586 100644 --- a/src/properties.jl +++ b/src/properties.jl @@ -1,10 +1,12 @@ -import Base: size, ndims -import LinearAlgebra: diag +import Base: size, ndims +import LinearAlgebra: diag -export ndoms, - domainType, +export ndoms, + domainType, codomainType, + domainStorageType, + codomainStorageType, is_linear, is_eye, is_null, @@ -16,7 +18,7 @@ export ndoms, is_full_row_rank, is_full_column_rank, is_sliced, - diag_AcA, + diag_AcA, diag_AAc, displacement, remove_displacement @@ -52,10 +54,63 @@ julia> codomainType(vcat(Eye(Complex{Float64},(10,)),DFT(Complex{Float64},10))) """ codomainType +""" +`inputStorageType(A::AbstractOperator)` + +Returns the type of the storage for the input of the operator. + +```julia +julia> inputStorageType(DFT(10)) +Array{Complex{Float64},1} + +julia> inputStorageType(hcat(Eye(Complex{Float64},(10,)),DFT(Complex{Float64},10))) +ArrayPartition{Complex{Float64},1,Array{Complex{Float64},1},Array{Complex{Float64},1}} +``` +""" +function domainStorageType(L::AbstractOperator) + dt = domainType(L) + return if dt isa Tuple + arrayTypes = Tuple{[Array{t, d} for (t, d) in zip(dt, length.(size(L,2)))]...} + ArrayPartition{promote_type(dt...), arrayTypes} + else + Array{dt, length(size(L,2))} + end +end + +""" +`outputStorageType(A::AbstractOperator)` + +Returns the type of the storage of for the output of the operator. + +```julia +julia> outputStorageType(DFT(10)) +Array{Complex{Float64},1} + +julia> outputStorageType(vcat(Eye(Complex{Float64},(10,)),DFT(Complex{Float64},10))) +ArrayPartition{Complex{Float64},1,Array{Complex{Float64},1},Array{Complex{Float64},1}} +``` +""" +function codomainStorageType(L::AbstractOperator) + dt = codomainType(L) + return if dt isa Tuple + arrayTypes = Tuple{[Array{t, d} for (t, d) in zip(dt, length.(size(L,1)))]...} + ArrayPartition{promote_type(dt...), arrayTypes} + else + Array{dt, length(size(L,1))} + end +end + +allocateInDomain(L::AbstractOperator, dims...=size(L,2)...) = allocate(domainStorageType(L), dims...) +allocateInCodomain(L::AbstractOperator, dims...=size(L,1)...) = allocate(codomainStorageType(L), dims...) + +allocate(::Type{T}, dims...) where {T <: AbstractArray} = T(undef, dims...) +allocate(::Type{ArrayPartition{T,S}}, dims...) where {T,S} = + ArrayPartition([allocate(s, d...) for (s,d) in zip(S.parameters, dims)]...) + """ `size(A::AbstractOperator, [dom,])` -Returns the size of an `AbstractOperator`. Type `size(A,1)` for the size of the codomain and `size(A,2)` for the size of the codomain. +Returns the size of an `AbstractOperator`. Type `size(A,1)` for the size of the codomain and `size(A,2)` for the size of the codomain. """ size(L::AbstractOperator, i::Int) = size(L)[i] @@ -286,19 +341,15 @@ julia> displacement(A) ``` """ -function displacement(S::AbstractOperator) - D = domainType(S) - if typeof(D) <: Tuple - x = ArrayPartition(zeros.(D, size(S, 2))...) - else - x = zeros(D, size(S, 2)) - end - d = S*x - if all(y -> y == d[1], d ) - return d[1] - else - return d - end +function displacement(S::AbstractOperator) + x = allocateInDomain(S) + fill!(x, 0) + d = S*x + if all(y -> y == d[1], d) + return d[1] + else + return d + end end """ @@ -319,13 +370,13 @@ end #printing function Base.show(io::IO, L::AbstractOperator) - print(io, fun_name(L)*" "*fun_space(L) ) + print(io, fun_name(L)*" "*fun_space(L)) end -function fun_space(L::AbstractOperator) +function fun_space(L::AbstractOperator) dom = fun_dom(L,2) codom = fun_dom(L,1) - return dom*"->"*codom + return dom*"->"*codom end function fun_dom(L::AbstractOperator,n::Int) diff --git a/src/syntax.jl b/src/syntax.jl index 228a78f..c649333 100644 --- a/src/syntax.jl +++ b/src/syntax.jl @@ -12,12 +12,7 @@ adjoint(L::T) where {T <: AbstractOperator} = AdjointOperator(L) ###### * ###### function (*)(L::AbstractOperator, b::AbstractArray) - C = codomainType(L) - if typeof(C) <: Tuple - y = ArrayPartition(zeros.(codomainType(L), size(L, 1))...) - else - y = zeros(codomainType(L), size(L, 1)) - end + y = allocateInCodomain(L) mul!(y, L, b) return y end @@ -86,8 +81,8 @@ function getindex(H::VCAT, idx::Union{AbstractArray,Int}) end end -getindex(H::A, idx::Union{AbstractArray,Int}) where {L <: HCAT, D, S, A<: AffineAdd{L,D,S}} = -AffineAdd(getindex(H.A, idx), H.d, S) +getindex(H::A, idx::Union{AbstractArray,Int}) where {L <: HCAT, D, S, A<: AffineAdd{L,D,S}} = +AffineAdd(getindex(H.A, idx), H.d, S) # get index of scale getindex(A::S,idx...) where {T, L, S <:Scale{T,L}} = Scale(A.coeff,A.coeff_conj,getindex(A.A,idx...))