diff --git a/.travis.yml b/.travis.yml index f4a0e2c..4c394b0 100644 --- a/.travis.yml +++ b/.travis.yml @@ -4,7 +4,8 @@ os: - linux - osx julia: - - 0.6 + - 0.7 + - 1.0 - nightly matrix: allow_failures: @@ -14,7 +15,7 @@ notifications: # script: # - julia -e 'Pkg.clone(pwd()); Pkg.build("AbstractOperators"); Pkg.test("AbstractOperators"; coverage=true)' after_success: - - julia -e 'cd(Pkg.dir("AbstractOperators")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(process_folder())' - - julia -e 'cd(Pkg.dir("AbstractOperators")); Pkg.add("Coverage"); using Coverage; Codecov.submit(Codecov.process_folder())' - - julia -e 'Pkg.add("Documenter")' - - julia -e 'cd(Pkg.dir("AbstractOperators")); include(joinpath("docs", "make.jl"))' + - julia -e 'using Pkg; cd(Pkg.dir("AbstractOperators")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(process_folder())' + - julia -e 'using Pkg; cd(Pkg.dir("AbstractOperators")); Pkg.add("Coverage"); using Coverage; Codecov.submit(Codecov.process_folder())' + - julia -e 'using Pkg; Pkg.add("Documenter")' + - julia -e 'using Pkg; cd(Pkg.dir("AbstractOperators")); include(joinpath("docs", "make.jl"))' diff --git a/README.md b/README.md index 838077f..dbd87a8 100644 --- a/README.md +++ b/README.md @@ -14,17 +14,15 @@ This is particularly useful in iterative algorithms and in first order large-sca ## Installation -To install the package, use the following in the Julia command line +To install the package, hit `]` from the Julia command line to enter the package manager, then ```julia -Pkg.add("AbstractOperators") +pkg> add AbstractOperators ``` -Remember to `Pkg.update()` to keep the package up to date. - ## Usage -With `using AbstractOperators` the package imports several methods like multiplication `*` and transposition `'` (and their in-place methods `A_mul_B!`, `Ac_mul_B!`). +With `using AbstractOperators` the package imports several methods like multiplication `*` and adjoint transposition `'` (and their in-place methods `mul!`). For example, one can create a 2-D Discrete Fourier Transform as follows: @@ -45,13 +43,13 @@ julia> y = A*x -0.905575+1.98446im 0.441199-0.913338im 0.315788+3.29666im 0.174273+0.318065im -0.905575-1.98446im 0.174273-0.318065im 0.315788-3.29666im 0.441199+0.913338im -julia> A_mul_B!(y,A,x) == A*x #in-place evaluation +julia> mul!(y, A, x) == A*x #in-place evaluation true julia> all(A'*y - *(size(x)...)*x .< 1e-12) true -julia> Ac_mul_B!(x,A,y) #in-place evaluation +julia> mul!(x, A',y) #in-place evaluation 3×4 Array{Float64,2}: -2.99091 9.45611 -19.799 1.6327 -11.1841 11.2365 -26.3614 11.7261 @@ -97,6 +95,11 @@ julia> V*ones(3,3) A list of the available `AbstractOperators` and calculus rules can be found in the [documentation](https://kul-forbes.github.io/AbstractOperators.jl/latest). +## Related packages + +* [ProximalOperators.jl](https://github.com/kul-forbes/ProximalOperators.jl) +* [ProximalAlgorithms.jl](https://github.com/kul-forbes/ProximalAlgorithms.jl) +* [StructuredOptimization.jl](https://github.com/kul-forbes/StructuredOptimization.jl) ## Credits diff --git a/REQUIRE b/REQUIRE index 137767a..ae922ae 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1 +1,4 @@ -julia 0.6 +julia 0.7 +AbstractFFTs v0.3.2 +FFTW 0.2.4 +DSP 0.5.1 diff --git a/appveyor.yml b/appveyor.yml index 6a3be44..77bae98 100644 --- a/appveyor.yml +++ b/appveyor.yml @@ -1,7 +1,17 @@ environment: matrix: - - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x64/0.6/julia-0.6-latest-win64.exe" -# - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x86/0.6/julia-0.6-latest-win32.exe" + - julia_version: 0.7 + - julia_version: 1 + - julia_version: nightly + +platform: +# - x86 # 32-bit + - x64 # 64-bit + +matrix: + allow_failures: + - julia_version: 1 + - julia_version: nightly branches: only: @@ -15,19 +25,18 @@ notifications: on_build_status_changed: false install: - - ps: "[System.Net.ServicePointManager]::SecurityProtocol = [System.Net.SecurityProtocolType]::Tls12" -# Download most recent Julia Windows binary - - ps: (new-object net.webclient).DownloadFile( - $env:JULIA_URL, - "C:\projects\julia-binary.exe") -# Run installer silently, output to C:\projects\julia - - C:\projects\julia-binary.exe /S /D=C:\projects\julia + - ps: iex ((new-object net.webclient).DownloadString("https://raw.githubusercontent.com/JuliaCI/Appveyor.jl/version-1/bin/install.ps1")) build_script: -# Need to convert from shallow to complete for Pkg.clone to work - - IF EXIST .git\shallow (git fetch --unshallow) - - C:\projects\julia\bin\julia -e "versioninfo(); - Pkg.clone(pwd(), \"AbstractOperators\"); Pkg.build(\"AbstractOperators\")" + - echo "%JL_BUILD_SCRIPT%" + - C:\julia\bin\julia -e "%JL_BUILD_SCRIPT%" test_script: - - C:\projects\julia\bin\julia --check-bounds=yes -e "Pkg.test(\"AbstractOperators\")" + - echo "%JL_TEST_SCRIPT%" + - C:\julia\bin\julia -e "%JL_TEST_SCRIPT%" + +# # Uncomment to support code coverage upload. Should only be enabled for packages +# # which would have coverage gaps without running on Windows +# on_success: +# - echo "%JL_CODECOV_SCRIPT%" +# - C:\julia\bin\julia -e "%JL_CODECOV_SCRIPT%" diff --git a/docs/src/index.md b/docs/src/index.md index ad28544..8c662de 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -7,17 +7,15 @@ This is particularly useful in iterative algorithms and in first order large-sca ## Installation -To install the package, use the following in the Julia command line +To install the package, hit `]` from the Julia command line to enter the package manager, then ```julia -Pkg.add("AbstractOperators") +pkg> add AbstractOperators ``` -Remember to `Pkg.update()` to keep the package up to date. - ## Usage -With `using AbstractOperators` the package imports several methods like multiplication `*` and transposition `'` (and their in-place methods `A_mul_B!`, `Ac_mul_B!`). +With `using AbstractOperators` the package imports several methods like multiplication `*` and adjoint transposition `'` (and their in-place methods `mul!`). For example, one can create a 2-D Discrete Fourier Transform as follows: @@ -38,13 +36,13 @@ julia> y = A*x -0.905575+1.98446im 0.441199-0.913338im 0.315788+3.29666im 0.174273+0.318065im -0.905575-1.98446im 0.174273-0.318065im 0.315788-3.29666im 0.441199+0.913338im -julia> A_mul_B!(y,A,x) == A*x #in-place evaluation +julia> mul!(y, A, x) == A*x #in-place evaluation true julia> all(A'*y - *(size(x)...)*x .< 1e-12) true -julia> Ac_mul_B!(x,A,y) #in-place evaluation +julia> mul!(x, A',y) #in-place evaluation 3×4 Array{Float64,2}: -2.99091 9.45611 -19.799 1.6327 -11.1841 11.2365 -26.3614 11.7261 @@ -90,6 +88,11 @@ julia> V*ones(3,3) A list of the available `AbstractOperators` and calculus rules can be found in the [documentation](https://kul-forbes.github.io/AbstractOperators.jl/latest). +## Related packages + +* [ProximalOperators.jl](https://github.com/kul-forbes/ProximalOperators.jl) +* [ProximalAlgorithms.jl](https://github.com/kul-forbes/ProximalAlgorithms.jl) +* [StructuredOptimization.jl](https://github.com/kul-forbes/StructuredOptimization.jl) ## Credits diff --git a/src/AbstractOperators.jl b/src/AbstractOperators.jl index b0517fb..4dd35c2 100644 --- a/src/AbstractOperators.jl +++ b/src/AbstractOperators.jl @@ -2,6 +2,8 @@ __precompile__() module AbstractOperators +using LinearAlgebra, AbstractFFTs, DSP, FFTW + # Block stuff include("utilities/block.jl") using AbstractOperators.BlockArrays @@ -11,7 +13,7 @@ abstract type AbstractOperator end abstract type LinearOperator <: AbstractOperator end abstract type NonLinearOperator <: AbstractOperator end -import Base: A_mul_B!, Ac_mul_B! +import LinearAlgebra: mul! export LinearOperator, NonLinearOperator, @@ -21,7 +23,9 @@ export LinearOperator, include("properties.jl") -# Linear operators +include("calculus/AdjointOperator.jl") + +## Linear operators include("linearoperators/MyLinOp.jl") include("linearoperators/Zeros.jl") @@ -30,7 +34,7 @@ include("linearoperators/Eye.jl") include("linearoperators/DiagOp.jl") include("linearoperators/GetIndex.jl") include("linearoperators/MatrixOp.jl") -include("linearoperators/MatrixMul.jl") +include("linearoperators/LMatrixOp.jl") include("linearoperators/DFT.jl") include("linearoperators/RDFT.jl") include("linearoperators/IRDFT.jl") @@ -53,7 +57,6 @@ include("calculus/Compose.jl") include("calculus/Reshape.jl") include("calculus/BroadCast.jl") include("calculus/Sum.jl") -include("calculus/Transpose.jl") include("calculus/AffineAdd.jl") include("calculus/Jacobian.jl") include("calculus/NonLinearCompose.jl") diff --git a/src/calculus/AdjointOperator.jl b/src/calculus/AdjointOperator.jl new file mode 100644 index 0000000..2c53237 --- /dev/null +++ b/src/calculus/AdjointOperator.jl @@ -0,0 +1,54 @@ +export AdjointOperator + +""" +`AdjointOperator(A::AbstractOperator)` + +Shorthand constructor: + +`'(A::AbstractOperator)` + +Returns the adjoint operator of `A`. + +```julia +julia> AdjointOperator(DFT(10)) +ℱᵃ ℂ^10 -> ℝ^10 + +julia> [DFT(10); DCT(10)]' +[ℱ;ℱc]ᵃ ℂ^10 ℝ^10 -> ℝ^10 +``` +""" +struct AdjointOperator{T <: AbstractOperator} <: AbstractOperator + A::T + 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 +end + +# Constructors + +AdjointOperator(L::AdjointOperator) = L.A + +# Properties + +size(L::AdjointOperator) = size(L.A,2), size(L.A,1) + +domainType(L::AdjointOperator) = codomainType(L.A) +codomainType(L::AdjointOperator) = domainType(L.A) + +fun_name(L::AdjointOperator) = fun_name(L.A)*"ᵃ" + +is_linear(L::AdjointOperator) = is_linear(L.A) +is_null(L::AdjointOperator) = is_null(L.A) +is_eye(L::AdjointOperator) = is_eye(L.A) +is_diagonal(L::AdjointOperator) = is_diagonal(L.A) +is_AcA_diagonal(L::AdjointOperator) = is_AAc_diagonal(L.A) +is_AAc_diagonal(L::AdjointOperator) = is_AcA_diagonal(L.A) +is_orthogonal(L::AdjointOperator) = is_orthogonal(L.A) +is_invertible(L::AdjointOperator) = is_invertible(L.A) +is_full_row_rank(L::AdjointOperator) = is_full_column_rank(L.A) +is_full_column_rank(L::AdjointOperator) = is_full_row_rank(L.A) + +diag(L::AdjointOperator) = diag(L.A) +diag_AcA(L::AdjointOperator) = diag_AAc(L.A) +diag_AAc(L::AdjointOperator) = diag_AcA(L.A) diff --git a/src/calculus/AffineAdd.jl b/src/calculus/AffineAdd.jl index 198a39b..8428559 100644 --- a/src/calculus/AffineAdd.jl +++ b/src/calculus/AffineAdd.jl @@ -46,17 +46,17 @@ end # Mappings # array -function A_mul_B!(y::DD, T::AffineAdd{L, D, true}, x) where {L <: AbstractOperator, DD, D} - A_mul_B!(y,T.A,x) +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 A_mul_B!(y::DD, T::AffineAdd{L, D, false}, x) where {L <: AbstractOperator, DD, D} - A_mul_B!(y,T.A,x) +function mul!(y::DD, T::AffineAdd{L, D, false}, x) where {L <: AbstractOperator, DD, D} + mul!(y,T.A,x) y .-= T.d end -Ac_mul_B!(y, T::AffineAdd{L, D}, x) where {L <: AbstractOperator, D} = Ac_mul_B!(y,T.A,x) +mul!(y, T::AdjointOperator{AffineAdd{L, D, S}}, x) where {L <: AbstractOperator, D, S} = mul!(y,T.A.A',x) # Properties @@ -93,7 +93,7 @@ function permute(T::AffineAdd{L,D,S}, p::AbstractVector{Int}) where {L,D,S} return AffineAdd(A,T.d,S) end -displacement(A::AffineAdd{L,D,true}) where {L,D} = A.d+displacement(A.A) -displacement(A::AffineAdd{L,D,false}) where {L,D} = -A.d+displacement(A.A) +displacement(A::AffineAdd{L,D,true}) where {L,D} = A.d .+ displacement(A.A) +displacement(A::AffineAdd{L,D,false}) where {L,D} = -A.d .+ displacement(A.A) remove_displacement(A::AffineAdd) = remove_displacement(A.A) diff --git a/src/calculus/BroadCast.jl b/src/calculus/BroadCast.jl index d4766f7..b1bac4f 100644 --- a/src/calculus/BroadCast.jl +++ b/src/calculus/BroadCast.jl @@ -27,7 +27,7 @@ struct BroadCast{N, D <: AbstractArray, M, C <: NTuple{M,Colon}, - I <: CartesianRange + I <: CartesianIndices } <: AbstractOperator A::L dim_out::NTuple{N,Int} @@ -44,12 +44,12 @@ struct BroadCast{N, Base.Broadcast.check_broadcast_shape(dim_out,size(A,1)) if size(A,1) != (1,) M = length(size(A,1)) - cols = ([Colon() for i = 1:M]...) - idxs = CartesianRange((dim_out[M+1:end]...)) + 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) else #singleton case M = 0 - idxs = CartesianRange((1,)) + idxs = CartesianIndices((1,)) new{N,L,T,D,M,NTuple{0,Colon},typeof(idxs)}(A,dim_out,bufC,bufD,(),idxs) end @@ -63,27 +63,29 @@ BroadCast(A, dim_out, zeros(codomainType(A),size(A,1)), zeros(domainType(A),size # Mappings -function A_mul_B!(y::CC, R::BroadCast{N,L,T,D,M}, b::DD) where {N,L,T,D,M,CC,DD} - A_mul_B!(R.bufC, R.A, b) +function mul!(y::CC, R::BroadCast{N,L,T,D,M}, b::DD) where {N,L,T,D,M,CC,DD} + mul!(R.bufC, R.A, b) y .= R.bufC end -function Ac_mul_B!(y::CC, R::BroadCast{N,L,T,D,M}, b::DD) where {N,L,T,D,M,CC,DD} +function mul!(y::CC, A::AdjointOperator{BroadCast{N,L,T,D,M,C,I}}, b::DD) where {N,L,T,D,M,C,I,CC,DD} + R = A.A fill!(y, 0.) for i in R.idxs - @views Ac_mul_B!(R.bufD, R.A, b[R.cols...,i.I...]) + @views mul!(R.bufD, R.A', b[R.cols...,i.I...]) y .+= R.bufD end return y end #singleton -function Ac_mul_B!(y::CC, R::BroadCast{N,L,T,D,0}, b::DD) where {N,L,T,D,CC,DD} +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) for bi in b bii[1] = bi - Ac_mul_B!(R.bufD, R.A, bii) + mul!(R.bufD, R.A', bii) y .+= R.bufD end return y @@ -91,17 +93,18 @@ 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 -function Ac_mul_B!(y::CC, R::BroadCast{2,L,T,D,2}, b::DD) where {L,T,D,CC,DD} +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.) for i = 1:size(b,1) - @views Ac_mul_B!(R.bufD, R.A, b[i,:]') + @views mul!(R.bufD, R.A, b[i,:]') y .+= R.bufD end return y end #singleton Eye -function Ac_mul_B!(y::CC, R::BroadCast{N,L,T,D,0}, b::DD) where {N,L<:Eye,T,D,CC,DD} +function mul!(y::CC, R::AdjointOperator{BroadCast{N,L,T,D,0,C,I}}, b::DD) where {N,L<:Eye,T,D,C,I,CC,DD} sum!(y,b) end @@ -115,5 +118,5 @@ codomainType( R::BroadCast) = codomainType(R.A) is_linear( R::BroadCast) = is_linear(R.A) is_null( R::BroadCast) = is_null(R.A) -fun_name(R::BroadCast) = "."fun_name(R.A) +fun_name(R::BroadCast) = "."*fun_name(R.A) remove_displacement(B::BroadCast) = BroadCast(remove_displacement(B.A), B.dim_out, B.bufC, B.bufD) diff --git a/src/calculus/Compose.jl b/src/calculus/Compose.jl index 701363b..f05b113 100644 --- a/src/calculus/Compose.jl +++ b/src/calculus/Compose.jl @@ -33,12 +33,9 @@ function Compose(L1::AbstractOperator, L2::AbstractOperator) if domainType(L1) != codomainType(L2) throw(DomainError()) end - Compose( L1, L2, Array{domainType(L1)}(size(L2,1)) ) + Compose( L1, L2, Array{domainType(L1)}(undef,size(L2,1)) ) end -Compose(A::NTuple{N,Any},buf::NTuple{M,Any}) where {N,M} = -Compose{N,M,typeof(A),typeof(buf)}(A,buf) - Compose(L1::AbstractOperator,L2::AbstractOperator,buf::AbstractArray) = Compose( (L2,L1), (buf,)) @@ -62,32 +59,32 @@ Compose(L1::Eye, L2::Eye) = L1 # Mappings -@generated function A_mul_B!(y::C, L::Compose{N,M,T1,T2},b::D) where {N,M,T1,T2,C,D} - ex = :(A_mul_B!(L.buf[1],L.A[1],b)) +@generated function mul!(y::C, L::Compose{N,M,T1,T2},b::D) where {N,M,T1,T2,C,D} + ex = :(mul!(L.buf[1],L.A[1],b)) for i = 2:M ex = quote $ex - A_mul_B!(L.buf[$i],L.A[$i], L.buf[$i-1]) + mul!(L.buf[$i],L.A[$i], L.buf[$i-1]) end end ex = quote $ex - A_mul_B!(y,L.A[N], L.buf[M]) + mul!(y,L.A[N], L.buf[M]) return y end end -@generated function Ac_mul_B!(y::D, L::Compose{N,M,T1,T2},b::C) where {N,M,T1,T2,C,D} - ex = :(Ac_mul_B!(L.buf[M],L.A[N],b)) +@generated function mul!(y::D, L::AdjointOperator{Compose{N,M,T1,T2}},b::C) where {N,M,T1,T2,C,D} + ex = :(mul!(L.A.buf[M],L.A.A[N]',b)) for i = M:-1:2 ex = quote $ex - Ac_mul_B!(L.buf[$i-1],L.A[$i], L.buf[$i]) + mul!(L.A.buf[$i-1],L.A.A[$i]', L.A.buf[$i]) end end ex = quote $ex - Ac_mul_B!(y,L.A[1], L.buf[1]) + mul!(y,L.A.A[1]', L.A.buf[1]) return y end end @@ -112,7 +109,6 @@ diag(L::Compose) = is_sliced(L) ? diag(L.A[2]) : prod(diag.(L.A)) diag_AAc(L::Compose) = is_AAc_diagonal(L) ? diag_AAc(L.A[2]) : error("is_AAc_diagonal( $(typeof(L) ) ) == false") # utils -import Base: permute function permute(C::Compose, p::AbstractVector{Int}) i = findfirst( x -> ndoms(x,2) > 1 , C.A) diff --git a/src/calculus/DCAT.jl b/src/calculus/DCAT.jl index dfaf947..6702f8c 100644 --- a/src/calculus/DCAT.jl +++ b/src/calculus/DCAT.jl @@ -3,17 +3,13 @@ export DCAT """ `DCAT(A::AbstractOperator...)` -Shorthand constructor: - -`blkdiag(α::Number,A::AbstractOperator)` - Block-diagonally concatenate `AbstractOperator`s. ```julia julia> D = DCAT(HCAT(Eye(2),Eye(2)),DFT(3)) [[I,I],0;0,ℱ] ℝ^2 ℝ^2 ℝ^4 -> ℝ^2 ℂ^3 -julia> blkdiag(Eye(10),Eye(10),FiniteDiff((4,4))) +julia> DCAT(Eye(10),Eye(10),FiniteDiff((4,4))) DCAT ℝ^10 ℝ^10 ℝ^(4, 4) -> ℝ^10 ℝ^10 ℝ^(3, 4) ``` @@ -52,13 +48,13 @@ function DCAT(A::Vararg{AbstractOperator}) K += 1 push!(idxs,K) else # stacked operator - idxs = push!(idxs,(collect(K+1:K+ndoms(A[i],d))...) ) + idxs = push!(idxs,(collect(K+1:K+ndoms(A[i],d))...,) ) for ii = 1:ndoms(A[i],d) K += 1 end end end - d == 1 ? (idxC = (idxs...)) : (idxD = (idxs...)) + d == 1 ? (idxC = (idxs...,)) : (idxD = (idxs...,)) end return DCAT(A,idxD,idxC) @@ -66,7 +62,7 @@ function DCAT(A::Vararg{AbstractOperator}) end # Mappings -@generated function A_mul_B!(y, H::DCAT{N,L,P1,P2}, b) where {N,L,P1,P2} +@generated function mul!(y, H::DCAT{N,L,P1,P2}, b) where {N,L,P1,P2} ex = :() @@ -74,33 +70,33 @@ end if fieldtype(P2,i) <: Int # flatten operator - # build A_mul_B!(y[H.idxC[i]], H.A[i], b) + # build mul!(y[H.idxC[i]], H.A[i], b) yy = :(y[H.idxC[$i]]) else - # staked operator - # build A_mul_B!(( y[H.idxC[i][1]], y[H.idxC[i][2]] ... ), H.A[i], b) + # stacked operator + # build mul!(( y[H.idxC[i][1]], y[H.idxC[i][2]] ... ), H.A[i], b) yy = "" for ii in eachindex(fieldnames(fieldtype(P2,i))) yy *= "y[H.idxC[$i][$ii]]," end - yy = parse(yy) + yy = Meta.parse(yy) end if fieldtype(P1,i) <: Int # flatten operator - # build Ac_mul_B!(H.buf, H.A[i], b[H.idxD[i]]) + # build mul!(H.buf, H.A[i], b[H.idxD[i]]) bb = :(b[H.idxD[$i]]) else - # staked operator - # build Ac_mul_B!(H.buf, H.A[i],( b[H.idxD[i][1]], b[H.idxD[i][2]] ... )) + # stacked operator + # build mul!(H.buf, H.A[i],( b[H.idxD[i][1]], b[H.idxD[i][2]] ... )) bb = "" for ii in eachindex(fieldnames(fieldtype(P1,i))) bb *= "b[H.idxD[$i][$ii]]," end - bb = parse(bb) + bb = Meta.parse(bb) end - ex = :($ex; A_mul_B!($yy,H.A[$i],$bb)) + ex = :($ex; mul!($yy,H.A[$i],$bb)) end ex = :($ex; return y) @@ -108,41 +104,41 @@ end end -@generated function Ac_mul_B!(y, H::DCAT{N,L,P1,P2}, b) where {N,L,P1,P2} +@generated function mul!(y, A::AdjointOperator{DCAT{N,L,P1,P2}}, b) where {N,L,P1,P2} - ex = :() + ex = :(H = A.A) for i = 1:N if fieldtype(P1,i) <: Int # flatten operator - # build Ac_mul_B!(y[H.idxD[i]], H.A[i], b) + # build mul!(y[H.idxD[i]], H.A[i]', b) yy = :(y[H.idxD[$i]]) else - # staked operator - # build Ac_mul_B!(( y[H.idxD[i][1]], y[H.idxD[i][2]] ... ), H.A[i], b) + # stacked operator + # build mul!(( y[H.idxD[i][1]], y[H.idxD[i][2]] ... ), H.A[i]', b) yy = "" for ii in eachindex(fieldnames(fieldtype(P1,i))) yy *= "y[H.idxD[$i][$ii]]," end - yy = parse(yy) + yy = Meta.parse(yy) end if fieldtype(P2,i) <: Int # flatten operator - # build Ac_mul_B!(H.buf, H.A[i], b[H.idxC[i]]) + # build mul!(H.buf, H.A[i]', b[H.idxC[i]]) bb = :(b[H.idxC[$i]]) else - # staked operator - # build Ac_mul_B!(H.buf, H.A[i],( b[H.idxC[i][1]], b[H.idxC[i][2]] ... )) + # stacked operator + # build mul!(H.buf, H.A[i]',( b[H.idxC[i][1]], b[H.idxC[i][2]] ... )) bb = "" for ii in eachindex(fieldnames(fieldtype(P2,i))) bb *= "b[H.idxC[$i][$ii]]," end - bb = parse(bb) + bb = Meta.parse(bb) end - ex = :($ex; Ac_mul_B!($yy,H.A[$i],$bb)) + ex = :($ex; mul!($yy,H.A[$i]',$bb)) end ex = :($ex; return y) @@ -159,26 +155,26 @@ function size(H::DCAT, i::Int) for s in size.(H.A,i) eltype(s) <: Int ? push!(sz,s) : push!(sz,s...) end - p = vcat([[idx... ] for idx in (i == 1? H.idxC : H.idxD) ]...) - ipermute!(sz,p) + p = vcat([[idx... ] for idx in (i == 1 ? H.idxC : H.idxD) ]...) + invpermute!(sz,p) - (sz...) + (sz...,) end -fun_name(L::DCAT) = length(L.A) == 2 ? "["fun_name(L.A[1])*",0;0,"*fun_name(L.A[2])*"]" : +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) domain = vcat([typeof(d)<:Tuple ? [d...] : d for d in domainType.(H.A)]...) p = vcat([[idx... ] for idx in H.idxD]...) - ipermute!(domain,p) - return (domain...) + invpermute!(domain,p) + return (domain...,) end 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]...) - ipermute!(codomain,p) - return (codomain...) + invpermute!(codomain,p) + return (codomain...,) end is_eye(L::DCAT) = all(is_eye.(L.A)) @@ -192,18 +188,16 @@ is_full_row_rank(L::DCAT) = all(is_full_row_rank.(L.A)) is_full_column_rank(L::DCAT) = all(is_full_column_rank.(L.A)) # utils -import Base: permute - function permute(H::DCAT{N,L,P1,P2}, p::AbstractVector{Int}) where {N,L,P1,P2} unfolded = vcat([[idx... ] for idx in H.idxD]...) - ipermute!(unfolded,p) + invpermute!(unfolded,p) new_part = () cnt = 0 for z in length.(H.idxD) - new_part = (new_part..., z == 1 ? unfolded[cnt+1] : (unfolded[cnt+1:z+cnt]...)) + new_part = (new_part..., z == 1 ? unfolded[cnt+1] : (unfolded[cnt+1:z+cnt]...,)) cnt += z end diff --git a/src/calculus/HCAT.jl b/src/calculus/HCAT.jl index 49848d3..80255cb 100644 --- a/src/calculus/HCAT.jl +++ b/src/calculus/HCAT.jl @@ -106,13 +106,13 @@ function HCAT(AA::NTuple{N,AbstractOperator}, buf::C) where {N,C} K += 1 push!(idxs,K) else # stacked operator - idxs = push!(idxs,(collect(K+1:K+ndoms(AA[i],2))...) ) + idxs = push!(idxs,(collect(K+1:K+ndoms(AA[i],2))...,) ) for ii = 1:ndoms(AA[i],2) K += 1 end end end - return HCAT(AA, (idxs...), buf, M) + return HCAT(AA, (idxs...,), buf, M) end end @@ -120,42 +120,42 @@ HCAT(A::AbstractOperator) = A # Mappings -@generated function A_mul_B!(y::C, H::HCAT{M,N,L,P,C}, b::DD) where {M,N,L,P,C,DD} +@generated function mul!(y::C, H::HCAT{M,N,L,P,C}, b::DD) where {M,N,L,P,C,DD} ex = :() if fieldtype(P,1) <: Int # flatten operator - # build A_mul_B!(y, H.A[1], b[H.idxs[1]]) + # build mul!(y, H.A[1], b[H.idxs[1]]) bb = :(b[H.idxs[1]]) else - # staked operator - # build A_mul_B!(y, H.A[1],( b[H.idxs[1][1]], b[H.idxs[1][2]] ... )) + # stacked operator + # build mul!(y, H.A[1],( b[H.idxs[1][1]], b[H.idxs[1][2]] ... )) bb = "" for ii in eachindex(fieldnames(fieldtype(P,1))) bb *= "b[H.idxs[1][$ii]]," end - bb = parse(bb) + bb = Meta.parse(bb) end - ex = :($ex; A_mul_B!(y,H.A[1],$bb)) # write on y + ex = :($ex; mul!(y,H.A[1],$bb)) # write on y for i = 2:N if fieldtype(P,i) <: Int # flatten operator - # build A_mul_B!(H.buf, H.A[i], b[H.idxs[i]]) + # build mul!(H.buf, H.A[i], b[H.idxs[i]]) bb = :(b[H.idxs[$i]]) else - # staked operator - # build A_mul_B!(H.buf, H.A[i],( b[H.idxs[i][1]], b[H.idxs[i][2]] ... )) + # stacked operator + # build mul!(H.buf, H.A[i],( b[H.idxs[i][1]], b[H.idxs[i][2]] ... )) bb = "" for ii in eachindex(fieldnames(fieldtype(P,i))) bb *= "b[H.idxs[$i][$ii]]," end - bb = parse(bb) + bb = Meta.parse(bb) end - ex = :($ex; A_mul_B!(H.buf,H.A[$i],$bb)) # write on H.buf + ex = :($ex; mul!(H.buf,H.A[$i],$bb)) # write on H.buf # sum H.buf with y if C <: AbstractArray @@ -172,27 +172,27 @@ HCAT(A::AbstractOperator) = A end -@generated function Ac_mul_B!(y::DD, H::HCAT{M,N,L,P,C}, b::C) where {M,N,L,P,C,DD} +@generated function mul!(y::DD, A::AdjointOperator{HCAT{M,N,L,P,C}}, b::C) where {M,N,L,P,C,DD} - ex = :() + ex = :(H = A.A) for i = 1:N if fieldtype(P,i) <: Int # flatten operator - # build Ac_mul_B!(y[H.idxs[i]], H.A[i], b) + # build mul!(y[H.idxs[i]], H.A[i]', b) yy = :(y[H.idxs[$i]]) else - # staked operator - # build Ac_mul_B!(( y[H.idxs[i][1]], y[H.idxs[i][2]] ... ), H.A[i], b) + # stacked operator + # build mul!(( y[H.idxs[i][1]], y[H.idxs[i][2]] ... ), H.A[i]', b) yy = "" for ii in eachindex(fieldnames(fieldtype(P,i))) yy *= "y[H.idxs[$i][$ii]]," end - yy = parse(yy) + yy = Meta.parse(yy) end - ex = :($ex; Ac_mul_B!($yy,H.A[$i],b)) + ex = :($ex; mul!($yy,H.A[$i]',b)) end ex = :($ex; return y) @@ -200,8 +200,8 @@ end end -# same as A_mul_B but skips `Zeros` -@generated function A_mul_B_skipZeros!(y::C, H::HCAT{M,N,L,P,C}, b::DD) where {M,N,L,P,C,DD} +# same as mul! but skips `Zeros` +@generated function mul_skipZeros!(y::C, H::HCAT{M,N,L,P,C}, b::DD) where {M,N,L,P,C,DD} ex = :() @@ -212,9 +212,9 @@ end for ii in eachindex(fieldnames(fieldtype(P,1))) bb *= "b[H.idxs[1][$ii]]," end - bb = parse(bb) + bb = Meta.parse(bb) end - ex = :($ex; A_mul_B!(y,H.A[1],$bb)) + ex = :($ex; mul!(y,H.A[1],$bb)) for i = 2:N if !(fieldtype(L,i) <: Zeros) @@ -226,10 +226,10 @@ end for ii in eachindex(fieldnames(fieldtype(P,i))) bb *= "b[H.idxs[$i][$ii]]," end - bb = parse(bb) + bb = Meta.parse(bb) end - ex = :($ex; A_mul_B!(H.buf,H.A[$i],$bb)) + ex = :($ex; mul!(H.buf,H.A[$i],$bb)) if C <: AbstractArray ex = :($ex; y .+= H.buf) @@ -246,10 +246,10 @@ end end -# same as Ac_mul_B but skips `Zeros` -@generated function Ac_mul_B_skipZeros!(y::DD, H::HCAT{M,N,L,P,C}, b::C) where {M,N,L,P,C,DD} +# same as mul! but skips `Zeros` +@generated function mul_skipZeros!(y::DD, A::AdjointOperator{HCAT{M,N,L,P,C}}, b::C) where {M,N,L,P,C,DD} - ex = :() + ex = :(H = A.A) for i = 1:N @@ -261,10 +261,10 @@ end for ii in eachindex(fieldnames(fieldtype(P,i))) yy *= "y[H.idxs[$i][$ii]]," end - yy = parse(yy) + yy = Meta.parse(yy) end - ex = :($ex; Ac_mul_B!($yy,H.A[$i],b)) + ex = :($ex; mul!($yy,H.A[$i]',b)) end end @@ -281,40 +281,38 @@ function size(H::HCAT) eltype(s) <: Int ? push!(size_in,s) : push!(size_in,s...) end p = vcat([[idx... ] for idx in H.idxs]...) - ipermute!(size_in,p) + invpermute!(size_in,p) - size(H.A[1],1), (size_in...) + size(H.A[1],1), (size_in...,) end -fun_name(L::HCAT) = length(L.A) == 2 ? "["fun_name(L.A[1])*","*fun_name(L.A[2])*"]" : "HCAT" +fun_name(L::HCAT) = length(L.A) == 2 ? "["*fun_name(L.A[1])*","*fun_name(L.A[2])*"]" : "HCAT" function domainType(H::HCAT) - domain = vcat([typeof(d)<:Tuple ? [d...] : d for d in domainType.(H.A)]...) + domain = vcat([typeof(d)<:Tuple ? [d...] : d for d in domainType.(H.A)]...) p = vcat([[idx... ] for idx in H.idxs]...) - ipermute!(domain,p) - return (domain...) + invpermute!(domain,p) + return (domain...,) end -codomainType(L::HCAT) = codomainType.(L.A[1]) +codomainType(L::HCAT) = codomainType.(Ref(L.A[1])) is_linear(L::HCAT) = all(is_linear.(L.A)) is_AAc_diagonal(L::HCAT) = all(is_AAc_diagonal.(L.A)) is_full_row_rank(L::HCAT) = any(is_full_row_rank.(L.A)) -diag_AAc(L::HCAT) = sum(diag_AAc.(L.A)) +diag_AAc(L::HCAT) = (+).(diag_AAc.(L.A)...) # utils -import Base: permute - function permute(H::HCAT{M,N,L,P,C}, p::AbstractVector{Int}) where {M,N,L,P,C} unfolded = vcat([[idx... ] for idx in H.idxs]...) - ipermute!(unfolded,p) + invpermute!(unfolded,p) new_part = () cnt = 0 for z in length.(H.idxs) - new_part = (new_part..., z == 1 ? unfolded[cnt+1] : (unfolded[cnt+1:z+cnt]...)) + new_part = (new_part..., z == 1 ? unfolded[cnt+1] : (unfolded[cnt+1:z+cnt]...,)) cnt += z end diff --git a/src/calculus/Hadamard.jl b/src/calculus/Hadamard.jl index 2031500..dc6889d 100644 --- a/src/calculus/Hadamard.jl +++ b/src/calculus/Hadamard.jl @@ -27,7 +27,6 @@ true ``` """ - struct Hadamard{M, C, V <: VCAT{M}} <: NonLinearOperator A::V buf::C @@ -64,8 +63,8 @@ function Hadamard(L1::AbstractOperator,L2::AbstractOperator) end # Mappings -function A_mul_B!(y, H::Hadamard{M,C,V}, b) where {M,C,V} - A_mul_B!(H.buf,H.A,b) +function mul!(y, H::Hadamard{M,C,V}, b) where {M,C,V} + mul!(H.buf,H.A,b) y .= H.buf[1] for i = 2:length(H.buf) @@ -77,12 +76,13 @@ end Jacobian(A::H, x::D) where {M, D<:Tuple, C, V, H <: Hadamard{M,C,V}} = HadamardJacobian(Jacobian(A.A,x),A.buf,A.buf2) -function Ac_mul_B!(y, J::HadamardJacobian{M,C,V}, b) where {M,C,V} +function mul!(y, A::AdjointOperator{HadamardJacobian{M,C,V}}, b) where {M,C,V} + J = A.A for i = 1:length(J.buf) c = (J.buf[1:i-1]...,J.buf[i+1:end]...,b) J.buf2[i] .= (.*)(c...) end - Ac_mul_B!(y, J.A, J.buf2) + mul!(y, J.A', J.buf2) end @@ -93,15 +93,13 @@ size(P::HadamardJacobian) = size(P.A[1],1), size(P.A[1],2) fun_name(L::Hadamard) = "⊙" fun_name(L::HadamardJacobian) = "J(⊙)" -domainType(L::Hadamard) = domainType.(L.A[1]) +domainType(L::Hadamard) = domainType.(Ref(L.A[1])) codomainType(L::Hadamard) = codomainType(L.A[1]) -domainType(L::HadamardJacobian) = domainType.(L.A[1]) +domainType(L::HadamardJacobian) = domainType.(Ref(L.A[1])) codomainType(L::HadamardJacobian) = codomainType(L.A[1]) # utils -import Base: permute - function permute(H::Hadamard, p::AbstractVector{Int}) A = VCAT([permute(a,p) for a in H.A.A]...) Hadamard(A,H.buf,H.buf2) diff --git a/src/calculus/Jacobian.jl b/src/calculus/Jacobian.jl index d4c86f6..9458b5a 100644 --- a/src/calculus/Jacobian.jl +++ b/src/calculus/Jacobian.jl @@ -36,7 +36,7 @@ function Jacobian(H::DCAT{N,L,P1,P2},x) where {N,L,P1,P2} if length(idx) == 1 A = (A...,jacobian(H.A[k],x[idx])) else - xx = ([x[i] for i in idx]...) + xx = ([x[i] for i in idx]...,) A = (A...,jacobian(H.A[k],xx)) end end @@ -49,15 +49,17 @@ function Jacobian(H::HCAT{M,N,L,P,C},x::D) where {M,N,L,P,C,D} if length(idx) == 1 A = (A...,jacobian(H.A[k],x[idx])) else - xx = ([x[i] for i in idx]...) + xx = ([x[i] for i in idx]...,) A = (A...,jacobian(H.A[k],xx)) end end HCAT(A,H.idxs,H.buf,M) end #Jacobian of VCAT -Jacobian(V::VCAT{M,N,L,P,C},x::D) where {M,N,L,P,C,D} = -VCAT(([Jacobian(a,x) for a in V.A]...), V.idxs, V.buf, M) +function Jacobian(V::VCAT{M,N,L,P,C},x::D) where {M,N,L,P,C,D} + JJ = ([Jacobian(a,x) for a in V.A]...,) + VCAT{M,N,typeof(JJ),P,C}(JJ, V.idxs, V.buf) +end #Jacobian of Compose function Jacobian(L::Compose, x::X) where {X<:AbstractArray} Compose(Jacobian.(L.A,(x,L.buf...)),L.buf) @@ -70,7 +72,7 @@ end Jacobian(R::Reshape{N,L},x::AbstractArray) where {N,L} = Reshape(Jacobian(R.A,x),R.dim_out) #Jacobian of Sum Jacobian(S::Sum{M,N,K,C,D},x::D) where {M,N,K,C,D} = -Sum(([Jacobian(a,x) for a in S.A]...),S.bufC,S.bufD,M,N) +Sum(([Jacobian(a,x) for a in S.A]...,),S.bufC,S.bufD,M,N) #Jacobian of Transpose Jacobian(T::Transpose{A}, x::AbstractArray) where {A <: AbstractOperator} = T #Jacobian of BroadCast diff --git a/src/calculus/NonLinearCompose.jl b/src/calculus/NonLinearCompose.jl index 69c68be..3e9a644 100644 --- a/src/calculus/NonLinearCompose.jl +++ b/src/calculus/NonLinearCompose.jl @@ -26,7 +26,6 @@ true ``` """ - struct NonLinearCompose{N, L1 <: HCAT{1}, L2 <: HCAT{1}, @@ -86,31 +85,33 @@ function Jacobian(P::NonLinearCompose{N,L,C,D},x::DD) where {M,N,L,C, end # Mappings -function A_mul_B!(y, P::NonLinearCompose{N,L,C,D}, b) where {N,L,C,D} - A_mul_B_skipZeros!(P.buf[1],P.A,b) - A_mul_B_skipZeros!(P.buf[2],P.B,b) - A_mul_B!(y,P.buf[1],P.buf[2]) +function mul!(y, P::NonLinearCompose{N,L1,L2,C,D}, b) where {N,L1,L2,C,D} + mul_skipZeros!(P.buf[1],P.A,b) + mul_skipZeros!(P.buf[2],P.B,b) + mul!(y,P.buf[1],P.buf[2]) end -function Ac_mul_B!(y, P::NonLinearComposeJac{N,L,C,D}, b) where {N,L,C,D} +function mul!(y, J::AdjointOperator{NonLinearComposeJac{N,L1,L2,C,D}}, b) where {N,L1,L2,C,D} - A_mul_Bc!(P.bufx[1],b,P.buf[2]) - Ac_mul_B_skipZeros!(y,P.A,P.bufx[1]) + P = J.A + mul!(P.bufx[1],b,P.buf[2]') + mul_skipZeros!(y,P.A',P.bufx[1]) - Ac_mul_B!(P.bufx[2],P.buf[1],b) - Ac_mul_B_skipZeros!(y,P.B,P.bufx[2]) + mul!(P.bufx[2],P.buf[1]',b) + mul_skipZeros!(y,P.B',P.bufx[2]) end # special case outer product -function Ac_mul_B!(y, P::NonLinearComposeJac{N,L,C,D}, b) where {N,L,C,D <: Tuple{AbstractVector,AbstractArray}} +function mul!(y, J::AdjointOperator{NonLinearComposeJac{N,L1,L2,C,D}}, b) where {N,L1,L2,C,D <: Tuple{AbstractVector,AbstractArray}} + P = J.A p = reshape(P.bufx[1], length(P.bufx[1]),1) - A_mul_Bc!(p,b,P.buf[2]) - Ac_mul_B_skipZeros!(y,P.A,P.bufx[1]) + mul!(p,b,P.buf[2]') + mul_skipZeros!(y,P.A',P.bufx[1]) - Ac_mul_B!(P.bufx[2],P.buf[1],b) - Ac_mul_B_skipZeros!(y,P.B,P.bufx[2]) + mul!(P.bufx[2],P.buf[1]',b) + mul_skipZeros!(y,P.B',P.bufx[2]) end @@ -127,18 +128,16 @@ function size(P::NonLinearComposeJac) size_out, size(P.A,2) end -fun_name(L::NonLinearCompose) = fun_name(L.A.A[1])"*"*fun_name(L.B.A[2]) -fun_name(L::NonLinearComposeJac) = fun_name(L.A.A[1])"*"*fun_name(L.B.A[2]) +fun_name(L::NonLinearCompose) = fun_name(L.A.A[1])*"*"*fun_name(L.B.A[2]) +fun_name(L::NonLinearComposeJac) = fun_name(L.A.A[1])*"*"*fun_name(L.B.A[2]) -domainType(L::NonLinearCompose) = domainType.(L.A) +domainType(L::NonLinearCompose) = domainType.(Ref(L.A)) codomainType(L::NonLinearCompose) = codomainType(L.A) -domainType(L::NonLinearComposeJac) = domainType.(L.A) +domainType(L::NonLinearComposeJac) = domainType.(Ref(L.A)) codomainType(L::NonLinearComposeJac) = codomainType(L.A) # utils -import Base: permute - function permute(P::NonLinearCompose{N,L,C,D}, p::AbstractVector{Int}) where {N,L,C,D} NonLinearCompose(permute(P.A,p),permute(P.B,p),P.buf,P.bufx) end diff --git a/src/calculus/Reshape.jl b/src/calculus/Reshape.jl index e73e6a7..b86913f 100644 --- a/src/calculus/Reshape.jl +++ b/src/calculus/Reshape.jl @@ -39,14 +39,15 @@ Reshape(A, dim_out) # Mappings -function A_mul_B!(y::C, R::Reshape{N,L}, b::D) where {N,L,C,D} +function mul!(y::C, R::Reshape{N,L}, b::D) where {N,L,C,D} y_res = reshape(y,size(R.A,1)) - A_mul_B!(y_res, R.A, b) + mul!(y_res, R.A, b) end -function Ac_mul_B!(y::D, R::Reshape{N,L}, b::C) where {N,L,C,D} +function mul!(y::D, A::AdjointOperator{Reshape{N,L}}, b::C) where {N,L,C,D} + R = A.A b_res = reshape(b,size(R.A,1)) - Ac_mul_B!(y, R.A, b_res) + mul!(y, R.A', b_res) end # Properties diff --git a/src/calculus/Scale.jl b/src/calculus/Scale.jl index 59c3df6..4e1876b 100644 --- a/src/calculus/Scale.jl +++ b/src/calculus/Scale.jl @@ -53,28 +53,30 @@ Scale(coeff::T,L::DiagOp) where {T<:RealOrComplex} = DiagOp(coeff*diag(L)) # Mappings -function A_mul_B!(y::C, L::Scale{T, A}, x::D) where {T, C <: AbstractArray, D, A <: AbstractOperator} - A_mul_B!(y, L.A, x) +function mul!(y::C, L::Scale{T, A}, x::D) where {T, C <: AbstractArray, D, A <: AbstractOperator} + mul!(y, L.A, x) y .*= L.coeff end -function A_mul_B!(y::C, L::Scale{T, A}, x::D) where {T, C <: Tuple, D, A <: AbstractOperator} - A_mul_B!(y, L.A, x) +function mul!(y::C, L::Scale{T, A}, x::D) where {T, C <: Tuple, D, A <: AbstractOperator} + mul!(y, L.A, x) for k in eachindex(y) y[k] .*= L.coeff end end -function Ac_mul_B!(y::D, L::Scale{T, A}, x::C) where {T, C, D <: AbstractArray, A <: AbstractOperator} - Ac_mul_B!(y, L.A, x) - y .*= L.coeff_conj +function mul!(y::D, S::AdjointOperator{Scale{T, A}}, x::C) where {T, C, D <: AbstractArray, A <: AbstractOperator} + L = S.A + mul!(y, L.A', x) + y .*= L.coeff_conj end -function Ac_mul_B!(y::D, L::Scale{T, A}, x::C) where {T, C, D <: Tuple, A <: AbstractOperator} - Ac_mul_B!(y, L.A, x) - for k in eachindex(y) - y[k] .*= L.coeff_conj - end +function mul!(y::D, S::AdjointOperator{Scale{T, A}}, x::C) where {T, C, D <: Tuple, A <: AbstractOperator} + L = S.A + mul!(y, L.A', x) + for k in eachindex(y) + y[k] .*= L.coeff_conj + end end # Properties diff --git a/src/calculus/Sum.jl b/src/calculus/Sum.jl index f0c3534..3f5c2ed 100644 --- a/src/calculus/Sum.jl +++ b/src/calculus/Sum.jl @@ -17,7 +17,7 @@ function Sum(A::L, bufC::C, bufD::D, M::Int, N::Int) where {C, D, K, L <: NTuple end if any([codomainType(A[1]) != codomainType(a) for a in A]) || any([codomainType(A[1]) != codomainType(a) for a in A]) - throw(DomainError()) + throw(DomainError(A,"cannot sum operator with different codomains")) end Sum{M, N, K, C, D, L}(A, bufC, bufD) end @@ -42,12 +42,12 @@ Sum((L1,L2.A...),L2.bufC,L2.bufD, M, N) # Mappings -@generated function A_mul_B!(y::C, S::Sum{M,N,K,C,D}, b::D) where {M,N,K,C,D} - ex = :(A_mul_B!(y,S.A[1],b)) +@generated function mul!(y::C, S::Sum{M,N,K,C,D}, b::D) where {M,N,K,C,D} + ex = :(mul!(y,S.A[1],b)) for i = 2:K ex = quote $ex - A_mul_B!(S.bufC,S.A[$i],b) + mul!(S.bufC,S.A[$i],b) end if C <: AbstractArray ex = :($ex; y .+= S.bufC) @@ -63,12 +63,12 @@ Sum((L1,L2.A...),L2.bufC,L2.bufD, M, N) end end -@generated function Ac_mul_B!(y::D, S::Sum{M,N,K,C,D}, b::C) where {M,N,K,C,D} - ex = :(Ac_mul_B!(y,S.A[1],b)) +@generated function mul!(y::D, A::AdjointOperator{Sum{M,N,K,C,D,L}}, b::C) where {M,N,K,C,D,L} + ex = :(S = A.A; mul!(y,S.A[1]',b)) for i = 2:K ex = quote $ex - Ac_mul_B!(S.bufD,S.A[$i],b) + mul!(S.bufD,S.A[$i]',b) end if D <: AbstractArray ex = :($ex; y .+= S.bufD) @@ -89,15 +89,15 @@ end size(L::Sum) = size(L.A[1]) domainType(S::Sum{M, N, K, C, D, L}) where {M,N,K,C,D<:AbstractArray,L} = domainType(S.A[1]) - domainType(S::Sum{M, N, K, C, D, L}) where {M,N,K,C,D<:Tuple ,L} = domainType.(S.A[1]) + domainType(S::Sum{M, N, K, C, D, L}) where {M,N,K,C,D<:Tuple ,L} = domainType.(Ref(S.A[1])) codomainType(S::Sum{M, N, K, C, D, L}) where {M,N,K,C<:AbstractArray,D,L} = codomainType(S.A[1]) -codomainType(S::Sum{M, N, K, C, D, L}) where {M,N,K,C<:Tuple ,D,L} = codomainType.(S.A[1]) +codomainType(S::Sum{M, N, K, C, D, L}) where {M,N,K,C<:Tuple ,D,L} = codomainType.(Ref(S.A[1])) fun_domain(S::Sum) = fun_domain(S.A[1]) fun_codomain(S::Sum) = fun_codomain(S.A[1]) fun_name(S::Sum) = -length(S.A) == 2 ? fun_name(S.A[1])"+"fun_name(S.A[2]) : "Σ" +length(S.A) == 2 ? fun_name(S.A[1])*"+"*fun_name(S.A[2]) : "Σ" is_linear(L::Sum) = all(is_linear.(L.A)) @@ -106,14 +106,11 @@ 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) = sum(diag.(L.A)) - +diag(L::Sum) = (+).(diag.(L.A)...,) # utils -import Base: permute - function permute(S::Sum{M,N}, p::AbstractVector{Int}) where {M,N} - AA = ([permute(A,p) for A in S.A]...) + AA = ([permute(A,p) for A in S.A]...,) return Sum(AA,S.bufC,S.bufD[p],M,N) end diff --git a/src/calculus/Transpose.jl b/src/calculus/Transpose.jl deleted file mode 100644 index 86444f4..0000000 --- a/src/calculus/Transpose.jl +++ /dev/null @@ -1,59 +0,0 @@ -export Transpose - -""" -`Transpose(A::AbstractOperator)` - -Shorthand constructor: - -`'(A::AbstractOperator)` - -Returns the adjoint operator of `A`. - -```julia -julia> Transpose(DFT(10)) -ℱᵃ ℂ^10 -> ℝ^10 - -julia> [DFT(10); DCT(10)]' -[ℱ;ℱc]ᵃ ℂ^10 ℝ^10 -> ℝ^10 -``` -""" -struct Transpose{T <: AbstractOperator} <: AbstractOperator - A::T - function Transpose(A::T) where {T<:AbstractOperator} - is_linear(A) == false && error("Cannot transpose a nonlinear operator. You might use `jacobian`") - new{T}(A) - end -end - -# Constructors - -Transpose(L::Transpose) = L.A - -# Mappings - - A_mul_B!(y, L::Transpose{T}, x) where {T<:AbstractOperator} = Ac_mul_B!(y, L.A, x) -Ac_mul_B!(y, L::Transpose{T}, x) where {T<:AbstractOperator} = A_mul_B!(y, L.A, x) - -# Properties - -size(L::Transpose) = size(L.A,2), size(L.A,1) - -domainType(L::Transpose) = codomainType(L.A) -codomainType(L::Transpose) = domainType(L.A) - -fun_name(L::Transpose) = fun_name(L.A)*"ᵃ" - -is_linear(L::Transpose) = is_linear(L.A) -is_null(L::Transpose) = is_null(L.A) -is_eye(L::Transpose) = is_eye(L.A) -is_diagonal(L::Transpose) = is_diagonal(L.A) -is_AcA_diagonal(L::Transpose) = is_AAc_diagonal(L.A) -is_AAc_diagonal(L::Transpose) = is_AcA_diagonal(L.A) -is_orthogonal(L::Transpose) = is_orthogonal(L.A) -is_invertible(L::Transpose) = is_invertible(L.A) -is_full_row_rank(L::Transpose) = is_full_column_rank(L.A) -is_full_column_rank(L::Transpose) = is_full_row_rank(L.A) - -diag(L::Transpose) = diag(L.A) -diag_AcA(L::Transpose) = diag_AAc(L.A) -diag_AAc(L::Transpose) = diag_AcA(L.A) diff --git a/src/calculus/VCAT.jl b/src/calculus/VCAT.jl index a06bbff..da48f3e 100644 --- a/src/calculus/VCAT.jl +++ b/src/calculus/VCAT.jl @@ -34,10 +34,10 @@ julia> V*ones(3) """ struct VCAT{M, # number of domains - N, # number of AbstractOperator - L <: NTuple{N,AbstractOperator}, - P <: NTuple{N,Union{Int,Tuple}}, - C <: Union{NTuple{M,AbstractArray}, AbstractArray}, + N, # number of AbstractOperator + L <: NTuple{N,AbstractOperator}, + P <: NTuple{N,Union{Int,Tuple}}, + C <: Union{NTuple{M,AbstractArray}, AbstractArray}, } <: AbstractOperator A::L # tuple of AbstractOperators idxs::P # indices @@ -62,6 +62,7 @@ function VCAT(A::L, idxs::P, buf::C, M::Int) where {N, if any([domainType(A[1]) != domainType(a) for a in A]) throw(error("operators must all share the same domainType!")) end + VCAT{M,N,L,P,C}(A, idxs, buf) end @@ -103,13 +104,13 @@ function VCAT(AA::NTuple{N,AbstractOperator}, buf::C) where {N,C} K += 1 push!(idxs,K) else # stacked operator - idxs = push!(idxs,(collect(K+1:K+ndoms(AA[i],1))...) ) + idxs = push!(idxs,(collect(K+1:K+ndoms(AA[i],1))...,) ) for ii = 1:ndoms(AA[i],1) K += 1 end end end - return VCAT(AA, (idxs...), buf, M) + return VCAT(AA, (idxs...,), buf, M) end end @@ -117,42 +118,42 @@ VCAT(A::AbstractOperator) = A # Mappings -@generated function Ac_mul_B!(y::C, H::VCAT{M,N,L,P,C}, b::DD) where {M,N,L,P,C,DD} +@generated function mul!(y::C, A::AdjointOperator{VCAT{M,N,L,P,C}}, b::DD) where {M,N,L,P,C,DD} - ex = :() + ex = :(H = A.A) if fieldtype(P,1) <: Int # flatten operator - # build Ac_mul_B!(y, H.A[1], b[H.idxs[1]]) + # build mul!(y, H.A[1]', b[H.idxs[1]]) bb = :(b[H.idxs[1]]) else - # staked operator - # build Ac_mul_B!(y, H.A[1],( b[H.idxs[1][1]], b[H.idxs[1][2]] ... )) + # stacked operator + # build mul!(y, H.A[1]',( b[H.idxs[1][1]], b[H.idxs[1][2]] ... )) bb = "" for ii in eachindex(fieldnames(fieldtype(P,1))) bb *= "b[H.idxs[1][$ii]]," end - bb = parse(bb) + bb = Meta.parse(bb) end - ex = :($ex; Ac_mul_B!(y,H.A[1],$bb)) # write on y + ex = :($ex; mul!(y,H.A[1]',$bb)) # write on y for i = 2:N if fieldtype(P,i) <: Int # flatten operator - # build Ac_mul_B!(H.buf, H.A[i], b[H.idxs[i]]) + # build mul!(H.buf, H.A[i], b[H.idxs[i]]) bb = :(b[H.idxs[$i]]) else - # staked operator - # build Ac_mul_B!(H.buf, H.A[i],( b[H.idxs[i][1]], b[H.idxs[i][2]] ... )) + # stacked operator + # build mul!(H.buf, H.A[i],( b[H.idxs[i][1]], b[H.idxs[i][2]] ... )) bb = "" for ii in eachindex(fieldnames(fieldtype(P,i))) bb *= "b[H.idxs[$i][$ii]]," end - bb = parse(bb) + bb = Meta.parse(bb) end - ex = :($ex; Ac_mul_B!(H.buf,H.A[$i],$bb)) # write on H.buf + ex = :($ex; mul!(H.buf,H.A[$i]',$bb)) # write on H.buf # sum H.buf with y if C <: AbstractArray @@ -169,7 +170,7 @@ VCAT(A::AbstractOperator) = A end -@generated function A_mul_B!(y::DD, H::VCAT{M,N,L,P,C}, b::C) where {M,N,L,P,C,DD} +@generated function mul!(y::DD, H::VCAT{M,N,L,P,C}, b::C) where {M,N,L,P,C,DD} ex = :() @@ -177,19 +178,19 @@ end if fieldtype(P,i) <: Int # flatten operator - # build A_mul_B!(y[H.idxs[i]], H.A[i], b) + # build mul!(y[H.idxs[i]], H.A[i], b) yy = :(y[H.idxs[$i]]) else - # staked operator - # build A_mul_B!(( y[H.idxs[i][1]], y[H.idxs[i][2]] ... ), H.A[i], b) + # stacked operator + # build mul!(( y[H.idxs[i][1]], y[H.idxs[i][2]] ... ), H.A[i], b) yy = "" for ii in eachindex(fieldnames(fieldtype(P,i))) yy *= "y[H.idxs[$i][$ii]]," end - yy = parse(yy) + yy = Meta.parse(yy) end - ex = :($ex; A_mul_B!($yy,H.A[$i],b)) + ex = :($ex; mul!($yy,H.A[$i],b)) end ex = :($ex; return y) @@ -197,8 +198,8 @@ end end -# same as Ac_mul_B but skips `Zeros` -@generated function Ac_mul_B_skipZeros!(y::C, H::VCAT{M,N,L,P,C}, b::DD) where {M,N,L,P,C,DD} +# same as mul but skips `Zeros` +@generated function mul_skipZeros!(y::C, H::VCAT{M,N,L,P,C}, b::DD) where {M,N,L,P,C,DD} ex = :() @@ -209,9 +210,9 @@ end for ii in eachindex(fieldnames(fieldtype(P,1))) bb *= "b[H.idxs[1][$ii]]," end - bb = parse(bb) + bb = Meta.parse(bb) end - ex = :($ex; Ac_mul_B!(y,H.A[1],$bb)) + ex = :($ex; mul!(y,H.A[1]',$bb)) for i = 2:N if !(fieldtype(L,i) <: Zeros) @@ -223,10 +224,10 @@ end for ii in eachindex(fieldnames(fieldtype(P,i))) bb *= "b[H.idxs[$i][$ii]]," end - bb = parse(bb) + bb = Meta.parse(bb) end - ex = :($ex; Ac_mul_B!(H.buf,H.A[$i],$bb)) + ex = :($ex; mul!(H.buf,H.A[$i]',$bb)) if C <: AbstractArray ex = :($ex; y .+= H.buf) @@ -243,8 +244,8 @@ end end -# same as A_mul_B but skips `Zeros` -@generated function A_mul_B_skipZeros!(y::DD, H::VCAT{M,N,L,P,C}, b::C) where {M,N,L,P,C,DD} +# same as mul but skips `Zeros` +@generated function mul_skipZeros!(y::DD, H::VCAT{M,N,L,P,C}, b::C) where {M,N,L,P,C,DD} ex = :() @@ -258,10 +259,10 @@ end for ii in eachindex(fieldnames(fieldtype(P,i))) yy *= "y[H.idxs[$i][$ii]]," end - yy = parse(yy) + yy = Meta.parse(yy) end - ex = :($ex; A_mul_B!($yy,H.A[$i],b)) + ex = :($ex; mul!($yy,H.A[$i],b)) end end @@ -278,40 +279,38 @@ function size(H::VCAT) eltype(s) <: Int ? push!(size_out,s) : push!(size_out,s...) end p = vcat([[idx... ] for idx in H.idxs]...) - ipermute!(size_out,p) + invpermute!(size_out,p) - (size_out...), size(H.A[1],2) + (size_out...,), size(H.A[1],2) end -fun_name(L::VCAT) = length(L.A) == 2 ? "["fun_name(L.A[1])*";"*fun_name(L.A[2])*"]" : "VCAT" +fun_name(L::VCAT) = length(L.A) == 2 ? "["*fun_name(L.A[1])*";"*fun_name(L.A[2])*"]" : "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]...) - ipermute!(codomain,p) - return (codomain...) + invpermute!(codomain,p) + return (codomain...,) end -domainType(L::VCAT) = domainType.(L.A[1]) +domainType(L::VCAT) = domainType.(Ref(L.A[1])) is_linear(L::VCAT) = all(is_linear.(L.A)) is_AcA_diagonal(L::VCAT) = all(is_AcA_diagonal.(L.A)) is_full_column_rank(L::VCAT) = any(is_full_column_rank.(L.A)) -diag_AcA(L::VCAT) = sum(diag_AcA.(L.A)) +diag_AcA(L::VCAT) = (+).(diag_AcA.(L.A)...,) # utils -import Base: permute - function permute(H::VCAT{M,N,L,P,C}, p::AbstractVector{Int}) where {M,N,L,P,C} unfolded = vcat([[idx... ] for idx in H.idxs]...) - ipermute!(unfolded,p) + invpermute!(unfolded,p) new_part = () cnt = 0 for z in length.(H.idxs) - new_part = (new_part..., z == 1 ? unfolded[cnt+1] : (unfolded[cnt+1:z+cnt]...)) + new_part = (new_part..., z == 1 ? unfolded[cnt+1] : (unfolded[cnt+1:z+cnt]...,)) cnt += z end diff --git a/src/linearoperators/Conv.jl b/src/linearoperators/Conv.jl index 0efaf37..e5dfd1f 100644 --- a/src/linearoperators/Conv.jl +++ b/src/linearoperators/Conv.jl @@ -8,7 +8,6 @@ export Conv 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}}, @@ -18,8 +17,8 @@ struct Conv{T, buf::H buf_c1::Hc buf_c2::Hc - R::Base.DFT.Plan - I::Base.DFT.Plan + R::AbstractFFTs.Plan + I::AbstractFFTs.Plan end # Constructors @@ -41,37 +40,37 @@ Conv(x::H, h::H) where {H} = Conv(eltype(x), size(x), h) # Mappings -function A_mul_B!(y::H, A::Conv{T,H}, b::H) where {T, 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) end - A_mul_B!(A.buf_c1, A.R, A.buf) + mul!(A.buf_c1, A.R, A.buf) for i in eachindex(A.buf) A.buf[i] = i <= length(b) ? b[i] : zero(T) end - A_mul_B!(A.buf_c2, A.R, A.buf) + mul!(A.buf_c2, A.R, A.buf) A.buf_c2 .*= A.buf_c1 - A_mul_B!(y,A.I,A.buf_c2) + mul!(y,A.I,A.buf_c2) end -function Ac_mul_B!(y::H, A::Conv{T,H}, b::H) where {T, H} - #y .= xcorr(b,A.h)[size(A,1)[1]:end-length(A.h)+1] #naive implementation - for i in eachindex(A.buf) - ii = length(A.buf)-i+1 - A.buf[ii] = i <= length(A.h) ? A.h[i] : zero(T) +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) end - A_mul_B!(A.buf_c1, A.R, A.buf) - for i in eachindex(A.buf) - A.buf[i] = b[i] + mul!(L.A.buf_c1, L.A.R, L.A.buf) + for i in eachindex(L.A.buf) + L.A.buf[i] = b[i] end - A_mul_B!(A.buf_c2, A.R, A.buf) - A.buf_c2 .*= A.buf_c1 - A_mul_B!(A.buf,A.I,A.buf_c2) - y[1] = A.buf[end] + mul!(L.A.buf_c2, L.A.R, L.A.buf) + L.A.buf_c2 .*= L.A.buf_c1 + mul!(L.A.buf,L.A.I,L.A.buf_c2) + y[1] = L.A.buf[end] for i = 2:length(y) - y[i] = A.buf[i-1] + y[i] = L.A.buf[i-1] end end diff --git a/src/linearoperators/DCT.jl b/src/linearoperators/DCT.jl index 3b7523e..5b0a455 100644 --- a/src/linearoperators/DCT.jl +++ b/src/linearoperators/DCT.jl @@ -30,11 +30,10 @@ julia> A*ones(3) ``` """ - struct DCT{N, C<:RealOrComplex, - T1<:Base.DFT.Plan, - T2<:Base.DFT.Plan} <: CosineTransform{N,C,T1,T2} + T1<:AbstractFFTs.Plan, + T2<:AbstractFFTs.Plan} <: CosineTransform{N,C,T1,T2} dim_in::NTuple{N,Int} A::T1 At::T2 @@ -68,11 +67,10 @@ julia> A*[1.;0.;0.] ``` """ - struct IDCT{N, C<:RealOrComplex, - T1<:Base.DFT.Plan, - T2<:Base.DFT.Plan} <: CosineTransform{N,C,T1,T2} + T1<:AbstractFFTs.Plan, + T2<:AbstractFFTs.Plan} <: CosineTransform{N,C,T1,T2} dim_in::NTuple{N,Int} A::T1 At::T2 @@ -106,20 +104,20 @@ end # Mappings -function A_mul_B!(y::AbstractArray{C,N},A::DCT{N,C,T1,T2},b::AbstractArray{C,N}) where {N,C,T1,T2} - A_mul_B!(y,A.A,b) +function mul!(y::AbstractArray{C,N},A::DCT{N,C,T1,T2},b::AbstractArray{C,N}) where {N,C,T1,T2} + mul!(y,A.A,b) end -function Ac_mul_B!(y::AbstractArray{C,N},A::DCT{N,C,T1,T2},b::AbstractArray{C,N}) where {N,C,T1,T2} - y .= A.At*b +function mul!(y::AbstractArray{C,N},A::AdjointOperator{DCT{N,C,T1,T2}},b::AbstractArray{C,N}) where {N,C,T1,T2} + y .= A.A.At*b end -function A_mul_B!(y::AbstractArray{C,N},A::IDCT{N,C,T1,T2},b::AbstractArray{C,N}) where {N,C,T1,T2} +function mul!(y::AbstractArray{C,N},A::IDCT{N,C,T1,T2},b::AbstractArray{C,N}) where {N,C,T1,T2} y .= A.A*b end -function Ac_mul_B!(y::AbstractArray{C,N},A::IDCT{N,C,T1,T2},b::AbstractArray{C,N}) where {N,C,T1,T2} - A_mul_B!(y,A.At,b) +function mul!(y::AbstractArray{C,N},A::AdjointOperator{IDCT{N,C,T1,T2}},b::AbstractArray{C,N}) where {N,C,T1,T2} + mul!(y,A.A.At,b) end # Properties diff --git a/src/linearoperators/DFT.jl b/src/linearoperators/DFT.jl index 4884703..6d1d96f 100644 --- a/src/linearoperators/DFT.jl +++ b/src/linearoperators/DFT.jl @@ -29,12 +29,11 @@ julia> A*ones(3) ``` """ - struct DFT{N, C<:RealOrComplex, D<:RealOrComplex, - T1<:Base.DFT.Plan, - T2<:Base.DFT.Plan} <: FourierTransform{N,C,D,T1,T2} + T1<:AbstractFFTs.Plan, + T2<:AbstractFFTs.Plan} <: FourierTransform{N,C,D,T1,T2} dim_in::NTuple{N,Int} A::T1 At::T2 @@ -68,12 +67,11 @@ julia> A*ones(3) ``` """ - struct IDFT{N, C<:RealOrComplex, D<:RealOrComplex, - T1<:Base.DFT.Plan, - T2<:Base.DFT.Plan} <: FourierTransform{N,C,D,T1,T2} + T1<:AbstractFFTs.Plan, + T2<:AbstractFFTs.Plan} <: FourierTransform{N,C,D,T1,T2} dim_in::NTuple{N,Int} A::T1 At::T2 @@ -116,57 +114,57 @@ IDFT(T::Type,dim_in::Vararg{Int}) = IDFT(T,dim_in) # Mappings -function A_mul_B!(y::AbstractArray{C,N}, - L::DFT{N,C,C,T1,T2}, - b::AbstractArray{C,N}) where {N,C<:Complex,T1,T2} - A_mul_B!(y,L.A,b) +function mul!(y::AbstractArray{C,N}, + L::DFT{N,C,C,T1,T2}, + b::AbstractArray{C,N}) where {N,C<:Complex,T1,T2} + mul!(y,L.A,b) end -function A_mul_B!(y::AbstractArray{C,N}, - L::DFT{N,C,D,T1,T2}, - b::AbstractArray{D,N}) where {N,C<:Complex,D<:Real,T1,T2} - A_mul_B!(y,L.A,complex(b)) +function mul!(y::AbstractArray{C,N}, + L::DFT{N,C,D,T1,T2}, + b::AbstractArray{D,N}) where {N,C<:Complex,D<:Real,T1,T2} + mul!(y,L.A,complex(b)) end -function Ac_mul_B!(y::AbstractArray{C,N}, - L::DFT{N,C,C,T1,T2}, - b::AbstractArray{C,N}) where {N,C<:Complex,T1,T2} - A_mul_B!(y,L.At,b) +function mul!(y::AbstractArray{C,N}, + L::AdjointOperator{DFT{N,C,C,T1,T2}}, + b::AbstractArray{C,N}) where {N,C<:Complex,T1,T2} + mul!(y,L.A.At,b) end -function Ac_mul_B!(y::AbstractArray{D,N}, - L::DFT{N,C,D,T1,T2}, - b::AbstractArray{C,N}) where {N,C<:Complex,D<:Real,T1,T2} +function mul!(y::AbstractArray{D,N}, + L::AdjointOperator{DFT{N,C,D,T1,T2}}, + b::AbstractArray{C,N}) where {N,C<:Complex,D<:Real,T1,T2} y2 = complex(y) - A_mul_B!(y2,L.At,b) + mul!(y2,L.A.At,b) y .= real.(y2) end -function A_mul_B!(y::AbstractArray{C,N}, - L::IDFT{N,C,C,T1,T2}, - b::AbstractArray{C,N}) where {N,C<:Complex,T1,T2} - A_mul_B!(y,L.A,b) +function mul!(y::AbstractArray{C,N}, + L::IDFT{N,C,C,T1,T2}, + b::AbstractArray{C,N}) where {N,C<:Complex,T1,T2} + mul!(y,L.A,b) end -function A_mul_B!(y::AbstractArray{C,N}, - L::IDFT{N,C,D,T1,T2}, - b::AbstractArray{D,N}) where {N,C<:Complex,D<:Real,T1,T2} - A_mul_B!(y,L.A,complex(b)) +function mul!(y::AbstractArray{C,N}, + L::IDFT{N,C,D,T1,T2}, + b::AbstractArray{D,N}) where {N,C<:Complex,D<:Real,T1,T2} + mul!(y,L.A,complex(b)) end -function Ac_mul_B!(y::AbstractArray{C,N}, - L::IDFT{N,C,C,T1,T2}, - b::AbstractArray{C,N}) where {N,C<:Complex,T1,T2} - A_mul_B!(y,L.At,b) +function mul!(y::AbstractArray{C,N}, + L::AdjointOperator{IDFT{N,C,C,T1,T2}}, + b::AbstractArray{C,N}) where {N,C<:Complex,T1,T2} + mul!(y,L.A.At,b) y ./= length(b) end -function Ac_mul_B!(y::AbstractArray{D,N}, - L::IDFT{N,C,D,T1,T2}, - b::AbstractArray{C,N}) where {N,C<:Complex,D<:Real,T1,T2} +function mul!(y::AbstractArray{D,N}, + L::AdjointOperator{IDFT{N,C,D,T1,T2}}, + b::AbstractArray{C,N}) where {N,C<:Complex,D<:Real,T1,T2} y2 = complex(y) - A_mul_B!(y2,L.At,b) + mul!(y2,L.A.At,b) y .= (/).(real.(y2), length(b)) end diff --git a/src/linearoperators/DiagOp.jl b/src/linearoperators/DiagOp.jl index 2a82ba1..4b200c8 100644 --- a/src/linearoperators/DiagOp.jl +++ b/src/linearoperators/DiagOp.jl @@ -19,7 +19,6 @@ julia> D*ones(2,2) ``` """ - struct DiagOp{N, D, C, T <: Union{AbstractArray{C,N}, Number}} <: LinearOperator dim_in::NTuple{N,Int} d::T @@ -46,16 +45,16 @@ DiagOp(DomainDim::NTuple{N,Int}, d::A) where {N, A <: Number} = DiagOp(Float64, # Mappings -function A_mul_B!(y::AbstractArray{C,N}, L::DiagOp{N, D, C, T}, b::AbstractArray{D,N}) where {N, D, C, T} +function mul!(y::AbstractArray{C,N}, L::DiagOp{N, D, C, T}, b::AbstractArray{D,N}) where {N, D, C, T} y .= L.d.*b end -function Ac_mul_B!(y::AbstractArray{D,N}, L::DiagOp{N, D, C, T}, b::AbstractArray{C,N}) where {N, D, C, T} - y .= conj.(L.d).*b +function mul!(y::AbstractArray{D,N}, L::AdjointOperator{DiagOp{N, D, C, T}}, b::AbstractArray{C,N}) where {N, D, C, T} + y .= conj.(L.A.d).*b end -function Ac_mul_B!(y::AbstractArray{D,N}, L::DiagOp{N, D, C, T}, b::AbstractArray{C,N}) where {N, D <: Real, C <: Complex, T} - y .= real.(conj.(L.d).*b) +function mul!(y::AbstractArray{D,N}, L::AdjointOperator{DiagOp{N, D, C, T}}, b::AbstractArray{C,N}) where {N, D <: Real, C <: Complex, T} + y .= real.(conj.(L.A.d).*b) end # Transformations (we'll see about this) diff --git a/src/linearoperators/Eye.jl b/src/linearoperators/Eye.jl index 2ccd216..65b1146 100644 --- a/src/linearoperators/Eye.jl +++ b/src/linearoperators/Eye.jl @@ -20,7 +20,6 @@ true ``` """ - struct Eye{T, N} <: LinearOperator dim::NTuple{N, Integer} end @@ -37,8 +36,8 @@ Eye(x::A) where {A <: AbstractArray} = Eye(eltype(x), size(x)) # Mappings -A_mul_B!(y::AbstractArray{T, N}, L::Eye{T, N}, b::AbstractArray{T, N}) where {T, N} = y .= b -Ac_mul_B!(y::AbstractArray{T, N}, L::Eye{T, N}, b::AbstractArray{T, N}) where {T, N} = A_mul_B!(y, L, b) +mul!(y::AbstractArray{T, N}, L::Eye{T, N}, b::AbstractArray{T, N}) where {T, N} = y .= b +mul!(y::AbstractArray{T, N}, L::AdjointOperator{Eye{T, N}}, b::AbstractArray{T, N}) where {T, N} = mul!(y, L.A, b) # Properties diag(L::Eye) = 1. diff --git a/src/linearoperators/Filt.jl b/src/linearoperators/Filt.jl index 9fc03fa..4455095 100644 --- a/src/linearoperators/Filt.jl +++ b/src/linearoperators/Filt.jl @@ -29,8 +29,8 @@ struct Filt{T,N} <: LinearOperator b ./= norml end # Pad the coefficients with zeros if needed - bs N && error("direction is bigger the number of dimension $N") - new{T,N,D}(dim_in) + idx = CartesianIndices(([i == D ? (2:d) : (1:d) for (i,d) in enumerate(dim_in)]...,)) + new{T,N,D,typeof(idx)}(dim_in,idx) end end @@ -44,32 +44,28 @@ FiniteDiff(x::AbstractArray{T,N}, dir::Int64 = 1) where {T,N} = FiniteDiff(elty # Mappings -@generated function A_mul_B!(y::AbstractArray{T,N}, - L::FiniteDiff{T,N,D}, - b::AbstractArray{T,N}) where {T,N,D} - o = ones(Int,N) - o[D] = 2 - I1 = CartesianIndex(o...) +@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 - I2 = CartesianIndex(size(b)) - for I in CartesianRange($I1,I2) + for I in L.idx y[I-$idx] = b[I]-b[I-$idx] end return y end end -@generated function Ac_mul_B!(y::AbstractArray{T,N}, - L::FiniteDiff{T,N,D}, - b::AbstractArray{T,N}) where {T,N,D} +@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 CartesianRange(size(y)) + 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] @@ -86,7 +82,7 @@ codomainType(L::FiniteDiff{T, N}) where {T, N} = T 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) + return ((dim_out...,), L.dim_in) end fun_name(L::FiniteDiff{T,N,1}) where {T,N} = "δx" diff --git a/src/linearoperators/GetIndex.jl b/src/linearoperators/GetIndex.jl index d3c9c5e..f111037 100644 --- a/src/linearoperators/GetIndex.jl +++ b/src/linearoperators/GetIndex.jl @@ -25,7 +25,6 @@ julia> GetIndex(randn(10,20,30),(1:2,1:4)) ``` """ - struct GetIndex{N,M,T<:Tuple} <: LinearOperator domainType::Type dim_out::NTuple{N,Int} @@ -52,13 +51,13 @@ GetIndex(x::AbstractArray, idx::Tuple) = GetIndex(eltype(x), size(x), idx) # Mappings -function A_mul_B!(y::Array{T1,N},L::GetIndex{N,M,T2},b::Array{T1,M}) where {T1,N,M,T2} +function mul!(y::Array{T1,N},L::GetIndex{N,M,T2},b::Array{T1,M}) where {T1,N,M,T2} y .= view(b,L.idx...) end -function Ac_mul_B!(y::Array{T1,M},L::GetIndex{N,M,T2},b::AbstractArray{T1,N}) where {T1,N,M,T2} +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.idx...) + setindex!(y,b,L.A.idx...) end # Properties diff --git a/src/linearoperators/IRDFT.jl b/src/linearoperators/IRDFT.jl index acdf1b6..d5e0f98 100644 --- a/src/linearoperators/IRDFT.jl +++ b/src/linearoperators/IRDFT.jl @@ -1,4 +1,3 @@ - export IRDFT """ @@ -18,12 +17,11 @@ julia> A = IRDFT((5,10,8),19,2) ``` """ - struct IRDFT{T <:Number, N, D, - T1<:Base.DFT.Plan, - T2<:Base.DFT.Plan, + T1<:AbstractFFTs.Plan, + T2<:AbstractFFTs.Plan, T3<:NTuple{N,Any} } <: LinearOperator dim_in::NTuple{N,Int} @@ -54,21 +52,22 @@ IRDFT(dim_in::NTuple{N,Int},d::Int,dims::Int=1) where {N} = IRDFT(zeros(Complex{ # Mappings -function A_mul_B!(y::C1, - L::IRDFT{T,N,D,T1,T2,T3}, - b::C2) where {N,T,D,T1,T2,T3,C1<:AbstractArray{T,N}, +function mul!(y::C1, + L::IRDFT{T,N,D,T1,T2,T3}, + b::C2) where {N,T,D,T1,T2,T3,C1<:AbstractArray{T,N}, C2<:AbstractArray{Complex{T},N}} - A_mul_B!(y,L.A,b) + mul!(y,L.A,b) end -function Ac_mul_B!(y::C2, - L::IRDFT{T,N,D,T1,T2,T3}, - b::C1) where {N,T,D,T1,T2,T3,C1<:AbstractArray{T,N}, +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_mul_B!(y,L.At,b) + A = L.A + mul!(y,A.At,b) y ./= size(b,D) - @views y[L.idx...] .*= 2 + @views y[A.idx...] .*= 2 return y end diff --git a/src/linearoperators/LBFGS.jl b/src/linearoperators/LBFGS.jl index 721bf6b..913be07 100644 --- a/src/linearoperators/LBFGS.jl +++ b/src/linearoperators/LBFGS.jl @@ -22,7 +22,6 @@ julia> d = L*grad; # compute new direction Use `reset!(L)` to cancel the memory of the operator. """ - mutable struct LBFGS{R, T <: BlockArray, M, I <: Integer} <: LinearOperator currmem::I curridx::I @@ -50,7 +49,7 @@ function LBFGS(domainType, dim_in, M::I) where {I <: Integer} end function LBFGS(dim_in, M::I) where {I <: Integer} - domainType = eltype(dim_in) <: Integer ? Float64 : ([Float64 for i in eachindex(dim_in)]...) + domainType = eltype(dim_in) <: Integer ? Float64 : ([Float64 for i in eachindex(dim_in)]...,) LBFGS(domainType, dim_in, M) end @@ -65,10 +64,11 @@ end See the documentation for `LBFGS`. """ - function update!(L::LBFGS{R, T, M, I}, x::T, x_prev::T, gradx::T, gradx_prev::T) where {R, T, M, I} - L.s .= x .- x_prev - L.y .= gradx .- gradx_prev + #L.s .= x .- x_prev + blockaxpy!(L.s, x, -1, x_prev) + #L.y .= gradx .- gradx_prev + blockaxpy!(L.y, gradx, -1, gradx_prev) ys = real(blockvecdot(L.s, L.y)) if ys > 0 L.curridx += 1 @@ -89,7 +89,6 @@ end Cancels the memory of `L`. """ - function reset!(L::LBFGS) L.currmem, L.curridx = 0, 0 L.H = 1.0 @@ -97,14 +96,16 @@ end # LBFGS operators are symmetric -Ac_mul_B!(x::T, L::LBFGS{R, T, M, I}, y::T) where {R, T, M, I} = A_mul_B!(x, L, y) +mul!(x::T, L::AdjointOperator{LBFGS{R, T, M, I}}, y::T) where {R, T, M, I} = mul!(x, L.A, y) # Two-loop recursion -function A_mul_B!(d::T, L::LBFGS{R, T, M, I}, gradx::T) where {R, T, M, I} - d .= gradx +function mul!(d::T, L::LBFGS{R, T, M, I}, gradx::T) where {R, T, M, I} + #d .= gradx + blockcopy!(d,gradx) idx = loop1!(d,L) - d .= (*).(L.H, d) + #d .= (*).(L.H, d) + blockscale!(d, L.H, d) d = loop2!(d,idx,L) end @@ -112,7 +113,8 @@ function loop1!(d::T, L::LBFGS{R, T, M, I}) where {R, T, M, I} idx = L.curridx for i = 1:L.currmem L.alphas[idx] = real(blockvecdot(L.s_M[idx], d))/L.ys_M[idx] - d .-= L.alphas[idx] .* L.y_M[idx] + #d .-= L.alphas[idx] .* L.y_M[idx] + blockcumscale!(d, -L.alphas[idx], L.y_M[idx]) idx -= 1 if idx == 0 idx = M end end @@ -124,7 +126,8 @@ function loop2!(d::T, idx::Int, L::LBFGS{R, T, M, I}) where {R, T, M, I} idx += 1 if idx > M idx = 1 end beta = real(blockvecdot(L.y_M[idx], d))/L.ys_M[idx] - d .+= (L.alphas[idx] - beta) .* L.s_M[idx] + #d .+= (L.alphas[idx] - beta) .* L.s_M[idx] + blockcumscale!(d, L.alphas[idx] - beta, L.s_M[idx]) end return d end diff --git a/src/linearoperators/MatrixMul.jl b/src/linearoperators/LMatrixOp.jl similarity index 72% rename from src/linearoperators/MatrixMul.jl rename to src/linearoperators/LMatrixOp.jl index 05452d0..91d0f28 100644 --- a/src/linearoperators/MatrixMul.jl +++ b/src/linearoperators/LMatrixOp.jl @@ -24,9 +24,8 @@ julia> op*ones(3,4) ``` """ - struct LMatrixOp{T, A <: Union{AbstractVector,AbstractMatrix}, - B <:Union{RowVector,AbstractMatrix}} <: LinearOperator + B <: AbstractMatrix} <: LinearOperator b::A bt::B n_row_in::Integer @@ -44,15 +43,15 @@ LMatrixOp(b::A, n_row_in::Int) where {T,A<:Union{AbstractVector{T},AbstractMatri LMatrixOp(T,(n_row_in,size(b,1)),b) # Mappings -A_mul_B!(y::C, L::LMatrixOp{T,A,B}, X::AbstractMatrix{T} ) where {T,A,B,C<:Union{AbstractVector,AbstractMatrix}} = -A_mul_B!(y,X,L.b) +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 Ac_mul_B!(y::AbstractMatrix{T}, L::LMatrixOp{T,A,B}, Y::AbstractVector{T} ) where {T,A,B} - y .= L.bt.*Y +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 Ac_mul_B!(y::AbstractMatrix{T}, L::LMatrixOp{T,A,B}, Y::AbstractMatrix{T} ) where {T,A,B} - A_mul_Bc!(y,Y,L.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 # Properties @@ -61,7 +60,7 @@ codomainType(L::LMatrixOp{T, A}) where {T, A} = T fun_name(L::LMatrixOp) = "(⋅)b" -size(L::LMatrixOp{T,A,B}) where {T,A <: AbstractVector,B <: RowVector} = (L.n_row_in,),(L.n_row_in, length(L.b)) +size(L::LMatrixOp{T,A,B}) where {T,A <: AbstractVector,B <: Adjoint} = (L.n_row_in,),(L.n_row_in, length(L.b)) size(L::LMatrixOp{T,A,B}) where {T,A <: AbstractMatrix,B <: AbstractMatrix} = (L.n_row_in,size(L.b,2)),(L.n_row_in, size(L.b,1)) #TODO diff --git a/src/linearoperators/MIMOFilt.jl b/src/linearoperators/MIMOFilt.jl index 139212d..5a00ba0 100644 --- a/src/linearoperators/MIMOFilt.jl +++ b/src/linearoperators/MIMOFilt.jl @@ -45,7 +45,6 @@ true ``` """ - struct MIMOFilt{T, A<:AbstractVector{T}} <: LinearOperator dim_out::Tuple{Int,Int} dim_in::Tuple{Int,Int} @@ -69,7 +68,7 @@ function MIMOFilt(domainType::Type, dim_in::NTuple{N,Int}, b::Vector, a::Vector) mod(length(b),dim_in[2]) !=0 && error("wrong number of filters") dim_out = (dim_in[1], div(length(b),dim_in[2]) ) - B,A,SI = Array{typeof(b[1]),1}(length(b)),Array{typeof(b[1]),1}(length(b)),Array{typeof(b[1]),1}(length(b)) + B,A,SI = similar(b),similar(b),similar(b) for i = 1:length(b) a[i][1] == 0 && error("filter vector a[$i][1] must be nonzero") @@ -90,8 +89,8 @@ function MIMOFilt(domainType::Type, dim_in::NTuple{N,Int}, b::Vector, a::Vector) end # Pad the coefficients with zeros if needed - bs MatrixOp(randn(20,10),4) ``` """ - struct MatrixOp{D, T, M <: AbstractMatrix{T}} <: LinearOperator A::M n_col_in::Integer @@ -60,13 +59,13 @@ convert(::Type{LinearOperator}, dom::Type, dim_in::Tuple, L::AbstractMatrix) = M # Mappings - A_mul_B!(y::AbstractArray, L::MatrixOp{D, T, M}, b::AbstractArray) where {D, T, M} = A_mul_B!(y, L.A, b) -Ac_mul_B!(y::AbstractArray, L::MatrixOp{D, T, M}, b::AbstractArray) where {D, T, M} = Ac_mul_B!(y, L.A, b) +mul!(y::AbstractArray, L::MatrixOp{D, T, M}, b::AbstractArray) where {D, T, M} = mul!(y, L.A, b) +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 Ac_mul_B!(y::AbstractArray, L::MatrixOp{D, T}, b::AbstractArray) where {D <: Real , T <: Complex} +function mul!(y::AbstractArray, L::AdjointOperator{MatrixOp{D, T, M}}, b::AbstractArray) where {D <: Real , T <: Complex, M} yc = zeros(T,size(y)) - Ac_mul_B!(yc, L.A, b) + mul!(yc, L.A.A', b) y .= real.(yc) end diff --git a/src/linearoperators/MyLinOp.jl b/src/linearoperators/MyLinOp.jl index da30127..45a8851 100644 --- a/src/linearoperators/MyLinOp.jl +++ b/src/linearoperators/MyLinOp.jl @@ -10,16 +10,15 @@ julia> n,m = 5,4; julia> A = randn(n,m); -julia> op = MyLinOp(Float64, (m,),(n,), (y,x) -> A_mul_B!(y,A,x), (y,x) -> Ac_mul_B!(y,A,x)) +julia> op = MyLinOp(Float64, (m,),(n,), (y,x) -> mul!(y,A,x), (y,x) -> mul!(y,A',x)) A ℝ^4 -> ℝ^5 -julia> op = MyLinOp(Float64, (m,), Float64, (n,), (y,x) -> A_mul_B!(y,A,x), (y,x) -> Ac_mul_B!(y,A,x)) +julia> op = MyLinOp(Float64, (m,), Float64, (n,), (y,x) -> mul!(y,A,x), (y,x) -> mul!(y,A',x)) A ℝ^4 -> ℝ^5 ``` """ - struct MyLinOp{N,M,C,D} <: LinearOperator dim_out::NTuple{N,Int} dim_in::NTuple{M,Int} @@ -39,14 +38,14 @@ MyLinOp{N,M, domainType, codomainType}(dim_out, dim_in, Fwd!, Adj! ) # Mappings - A_mul_B!(y::Array{C,N}, L::MyLinOp{N,M,C,D}, b::Array{D,M}) where {N,M,C,D} = L.Fwd!(y,b) -Ac_mul_B!(y::Array{C,N}, L::MyLinOp{N,M,C,D}, b::Array{D,M}) where {N,M,C,D} = L.Adj!(y,b) +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) # Properties size(L::MyLinOp) = (L.dim_out, L.dim_in) -codomainType{N,M,C,D}(L::MyLinOp{N,M,C,D}) = C - domainType{N,M,C,D}(L::MyLinOp{N,M,C,D}) = D +codomainType(L::MyLinOp{N,M,C,D}) where {N,M,C,D} = C + domainType(L::MyLinOp{N,M,C,D}) where {N,M,C,D} = D fun_name(L::MyLinOp) = "A" diff --git a/src/linearoperators/RDFT.jl b/src/linearoperators/RDFT.jl index 1eb14b4..97967e1 100644 --- a/src/linearoperators/RDFT.jl +++ b/src/linearoperators/RDFT.jl @@ -1,5 +1,4 @@ - -export RDFT, IRDFT +export RDFT """ `RDFT([domainType=Float64::Type,] dim_in::Tuple [,dims=1])` @@ -20,11 +19,10 @@ julia> RDFT((10,10,10),2) ``` """ - struct RDFT{T <:Number, N, - T1<:Base.DFT.Plan, - T2<:Base.DFT.Plan, + T1<:AbstractFFTs.Plan, + T2<:AbstractFFTs.Plan, T3<:AbstractArray{Complex{T},N} } <: LinearOperator dim_in::NTuple{N,Int} @@ -62,18 +60,19 @@ RDFT(T::Type,dim_in::Vararg{Int}) = RDFT(T,dim_in) # Mappings -function A_mul_B!(y::T3, - L::RDFT{T,N,T1,T2,T3}, - b::T4) where {N,T,T1,T2,T3,T4 <: AbstractArray{T,N}} - A_mul_B!(y,L.A,b) +function mul!(y::T3, + L::RDFT{T,N,T1,T2,T3}, + b::T4) where {N,T,T1,T2,T3,T4 <: AbstractArray{T,N}} + mul!(y,L.A,b) end -function Ac_mul_B!(y::T4, - L::RDFT{T,N,T1,T2,T3}, - b::T3) where {N,T,T1,T2,T3,T4 <: AbstractArray{T,N}} - A_mul_B!(L.b2,L.Zp,b) - A_mul_B!(L.y2,L.At,L.b2) - y .= real.(L.y2) +function mul!(y::T4, + L::AdjointOperator{RDFT{T,N,T1,T2,T3}}, + b::T3) where {N,T,T1,T2,T3,T4 <: AbstractArray{T,N}} + A = L.A + mul!(A.b2,A.Zp,b) + mul!(A.y2,A.At,A.b2) + y .= real.(A.y2) return y end diff --git a/src/linearoperators/Variation.jl b/src/linearoperators/Variation.jl index 7e5676f..68de91d 100644 --- a/src/linearoperators/Variation.jl +++ b/src/linearoperators/Variation.jl @@ -26,7 +26,6 @@ julia> Variation(ones(2,2))*[1. 2.; 1. 2.] ``` """ - struct Variation{T,N} <: LinearOperator dim_in::NTuple{N,Int} end @@ -45,22 +44,22 @@ Variation(x::AbstractArray) = Variation(eltype(x), size(x)) # Mappings -@generated function A_mul_B!(y::AbstractArray{T,2}, - A::Variation{T,N}, b::AbstractArray{T,N}) where {T,N} +@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...) + z = (z...,) ex = :($ex; y[cnt,$i] = I[$i] == 1 ? b[I+CartesianIndex($z)]-b[I] : b[I]-b[I-CartesianIndex($z)]) end ex2 = quote cnt = 0 - for I in CartesianRange(size(b)) + for I in CartesianIndices(size(b)) cnt += 1 $ex end @@ -68,8 +67,8 @@ Variation(x::AbstractArray) = Variation(eltype(x), size(x)) end end -@generated function Ac_mul_B!(y::AbstractArray{T,N}, - A::Variation{T,N}, b::AbstractArray{T,2}) where {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]) : I[1] == 2 ? b[cnt,1] + b[cnt-1,1] - b[cnt+1,1] : @@ -89,7 +88,7 @@ end ex2 = quote cnt = 0 - for I in CartesianRange(size(y)) + for I in CartesianIndices(size(y)) cnt += 1 $ex end diff --git a/src/linearoperators/Xcorr.jl b/src/linearoperators/Xcorr.jl index b276c74..042392c 100644 --- a/src/linearoperators/Xcorr.jl +++ b/src/linearoperators/Xcorr.jl @@ -1,4 +1,5 @@ export Xcorr +#TODO make more efficient """ `Xcorr([domainType=Float64::Type,] dim_in::Tuple, h::AbstractVector)` @@ -8,7 +9,6 @@ export Xcorr 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 dim_in::Tuple{Int} h::H @@ -24,11 +24,12 @@ Xcorr(x::H, h::H) where {H} = Xcorr(eltype(x), size(x), h) # Mappings -function A_mul_B!(y::H,A::Xcorr{T,H},b::H) where {T,H} +function mul!(y::H,A::Xcorr{T,H},b::H) where {T,H} y .= xcorr(b,A.h) end -function Ac_mul_B!(y::H,A::Xcorr{T,H},b::H) where {T,H} +function mul!(y::H,L::AdjointOperator{Xcorr{T,H}},b::H) where {T,H} + A = L.A l =floor(Int64,size(A,1)[1]/2) idx = l+1:l+length(y) y .= conv(b,A.h)[idx] diff --git a/src/linearoperators/ZeroPad.jl b/src/linearoperators/ZeroPad.jl index df79749..6bdd992 100644 --- a/src/linearoperators/ZeroPad.jl +++ b/src/linearoperators/ZeroPad.jl @@ -20,7 +20,6 @@ julia> Z*ones(2,2) ``` """ - struct ZeroPad{T,N} <: LinearOperator dim_in::NTuple{N,Int} zp::NTuple{N,Int} @@ -41,7 +40,7 @@ ZeroPad(x::AbstractArray, zp::NTuple{N,Int}) where {N} = ZeroPad( ZeroPad(x::AbstractArray, zp::Vararg{Int,N}) where {N} = ZeroPad(eltype(x), size(x), zp) # Mappings -@generated function A_mul_B!(y::AbstractArray{T,N},L::ZeroPad{T,N},b::AbstractArray{T,N}) where {T,N} +@generated function mul!(y::AbstractArray{T,N},L::ZeroPad{T,N},b::AbstractArray{T,N}) where {T,N} # builds @@ -70,11 +69,11 @@ ZeroPad(x::AbstractArray, zp::Vararg{Int,N}) where {N} = ZeroPad( ex = ex[1:end-1] #remove , ex *= "] : 0. end" - ex = parse(ex) + ex = Meta.parse(ex) end -@generated function Ac_mul_B!(y::AbstractArray{T,N},L::ZeroPad{T,N},b::AbstractArray{T,N}) where {T,N} +@generated function mul!(y::AbstractArray{T,N},L::AdjointOperator{ZeroPad{T,N}},b::AbstractArray{T,N}) where {T,N} #builds #for l = 1:size(y,1), m = 1:size(y,2) @@ -96,7 +95,7 @@ end ex *= idx ex *= "] end" - ex = parse(ex) + ex = Meta.parse(ex) end diff --git a/src/linearoperators/Zeros.jl b/src/linearoperators/Zeros.jl index 05df317..382468a 100644 --- a/src/linearoperators/Zeros.jl +++ b/src/linearoperators/Zeros.jl @@ -16,7 +16,6 @@ julia> Zeros([Eye(10,20) Eye(10,20)]) ``` """ - struct Zeros{C,N,D,M} <: LinearOperator dim_out::NTuple{N, Int} dim_in::NTuple{M, Int} @@ -45,8 +44,8 @@ Zeros(A::AbstractOperator) = Zeros(domainType(A),size(A,2),codomainType(A),size( # Mappings - A_mul_B!(y::AbstractArray{C,N}, A::Zeros{C,N,D,M}, b::AbstractArray{D,M}) where {C,N,D,M}= fill!(y,zero(C)) -Ac_mul_B!(y::AbstractArray{D,M}, A::Zeros{C,N,D,M}, b::AbstractArray{C,N}) where {C,N,D,M}= fill!(y,zero(D)) +mul!(y::AbstractArray{C,N}, A::Zeros{C,N,D,M}, b::AbstractArray{D,M}) where {C,N,D,M}= fill!(y,zero(C)) +mul!(y::AbstractArray{D,M}, A::AdjointOperator{Zeros{C,N,D,M}}, b::AbstractArray{C,N}) where {C,N,D,M}= fill!(y,zero(D)) # Properties diff --git a/src/nonlinearoperators/Atan.jl b/src/nonlinearoperators/Atan.jl index 360f2c2..f1669db 100644 --- a/src/nonlinearoperators/Atan.jl +++ b/src/nonlinearoperators/Atan.jl @@ -20,13 +20,14 @@ end Atan(DomainDim::NTuple{N,Int}) where {N} = Atan{Float64,N}(DomainDim) Atan(DomainDim::Vararg{Int}) = Atan{Float64,length(DomainDim)}(DomainDim) -function A_mul_B!(y::AbstractArray{T,N}, L::Atan{T,N}, x::AbstractArray{T,N}) where {T,N} +function mul!(y::AbstractArray{T,N}, L::Atan{T,N}, x::AbstractArray{T,N}) where {T,N} y .= atan.(x) end -function Ac_mul_B!(y::AbstractArray{T,N}, - L::Jacobian{A}, - b::AbstractArray{T,N}) where {T,N, A<: Atan{T,N}} +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 end diff --git a/src/nonlinearoperators/Cos.jl b/src/nonlinearoperators/Cos.jl index 5f96716..f5fd9cf 100644 --- a/src/nonlinearoperators/Cos.jl +++ b/src/nonlinearoperators/Cos.jl @@ -20,13 +20,14 @@ end Cos(DomainDim::NTuple{N,Int}) where {N} = Cos{Float64,N}(DomainDim) Cos(DomainDim::Vararg{Int}) = Cos{Float64,length(DomainDim)}(DomainDim) -function A_mul_B!(y::AbstractArray{T,N}, L::Cos{T,N}, x::AbstractArray{T,N}) where {T,N} +function mul!(y::AbstractArray{T,N}, L::Cos{T,N}, x::AbstractArray{T,N}) where {T,N} y .= cos.(x) end -function Ac_mul_B!(y::AbstractArray{T,N}, - L::Jacobian{A}, - b::AbstractArray{T,N}) where {T,N, A<: Cos{T,N}} +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 end diff --git a/src/nonlinearoperators/Exp.jl b/src/nonlinearoperators/Exp.jl index 166c3e9..e85fa7b 100644 --- a/src/nonlinearoperators/Exp.jl +++ b/src/nonlinearoperators/Exp.jl @@ -20,13 +20,14 @@ end Exp(DomainDim::NTuple{N,Int}) where {N} = Exp{Float64,N}(DomainDim) Exp(DomainDim::Vararg{Int}) = Exp{Float64,length(DomainDim)}(DomainDim) -function A_mul_B!(y::AbstractArray{T,N}, L::Exp{T,N}, x::AbstractArray{T,N}) where {T,N} +function mul!(y::AbstractArray{T,N}, L::Exp{T,N}, x::AbstractArray{T,N}) where {T,N} y .= exp.(x) end -function Ac_mul_B!(y::AbstractArray{T,N}, - L::Jacobian{A}, - b::AbstractArray{T,N}) where {T,N, A<: Exp{T,N}} +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 end diff --git a/src/nonlinearoperators/Pow.jl b/src/nonlinearoperators/Pow.jl index 6a26488..ea05751 100644 --- a/src/nonlinearoperators/Pow.jl +++ b/src/nonlinearoperators/Pow.jl @@ -17,13 +17,14 @@ end Pow(DomainDim::NTuple{N,Int}, p::I) where {N, I <: Real} = Pow{Float64,N,I}(DomainDim,p) -function A_mul_B!(y::AbstractArray{T,N}, L::Pow{T,N,I}, x::AbstractArray{T,N}) where {T,N,I} +function mul!(y::AbstractArray{T,N}, L::Pow{T,N,I}, x::AbstractArray{T,N}) where {T,N,I} y .= x.^L.p end -function Ac_mul_B!(y::AbstractArray{T,N}, - L::Jacobian{A}, - b::AbstractArray{T,N}) where {T, N, I, A<: Pow{T, N, I}} +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 end diff --git a/src/nonlinearoperators/Sigmoid.jl b/src/nonlinearoperators/Sigmoid.jl index c28d465..1d5beec 100644 --- a/src/nonlinearoperators/Sigmoid.jl +++ b/src/nonlinearoperators/Sigmoid.jl @@ -20,16 +20,17 @@ end Sigmoid(DomainDim::NTuple{N,Int}, gamma::G=1.) where {N,G} = Sigmoid{Float64,N,G}(DomainDim,gamma) -function A_mul_B!(y::AbstractArray{T,N}, L::Sigmoid{T,N,G}, x::AbstractArray{T,N}) where {T,N,G} - y .= (1.+exp.(-L.gamma.*x)).^(-1) +function mul!(y::AbstractArray{T,N}, L::Sigmoid{T,N,G}, x::AbstractArray{T,N}) where {T,N,G} + y .= (1 .+exp.(-L.gamma.*x)).^(-1) end -function Ac_mul_B!(y::AbstractArray{T,N}, - L::Jacobian{A}, - b::AbstractArray{T,N}) where {T,N,G, A<: Sigmoid{T,N,G}} +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 449eef7..3c79628 100644 --- a/src/nonlinearoperators/Sin.jl +++ b/src/nonlinearoperators/Sin.jl @@ -20,13 +20,14 @@ end Sin(DomainDim::NTuple{N,Int}) where {N} = Sin{Float64,N}(DomainDim) Sin(DomainDim::Vararg{Int}) = Sin{Float64,length(DomainDim)}(DomainDim) -function A_mul_B!(y::AbstractArray{T,N}, L::Sin{T,N}, x::AbstractArray{T,N}) where {T,N} +function mul!(y::AbstractArray{T,N}, L::Sin{T,N}, x::AbstractArray{T,N}) where {T,N} y .= sin.(x) end -function Ac_mul_B!(y::AbstractArray{T,N}, - L::Jacobian{A}, - b::AbstractArray{T,N}) where {T,N, A<: Sin{T,N}} +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 end diff --git a/src/nonlinearoperators/SoftMax.jl b/src/nonlinearoperators/SoftMax.jl index 22bcfd1..f26408d 100644 --- a/src/nonlinearoperators/SoftMax.jl +++ b/src/nonlinearoperators/SoftMax.jl @@ -20,19 +20,20 @@ end SoftMax(DomainDim::NTuple{N,Int}) where {N} = SoftMax(Float64,DomainDim) -function A_mul_B!(y::AbstractArray{T,N}, L::SoftMax{T,N}, x::AbstractArray{T,N}) where {T,N} +function mul!(y::AbstractArray{T,N}, L::SoftMax{T,N}, x::AbstractArray{T,N}) where {T,N} y .= exp.(x.-maximum(x)) y ./= sum(y) end -function Ac_mul_B!(y::AbstractArray{T,N}, - L::Jacobian{A}, - b::AbstractArray{T,N}) where {T,N, A<: SoftMax{T,N}} +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]*vecdot(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 409c65c..0ee54b1 100644 --- a/src/nonlinearoperators/SoftPlus.jl +++ b/src/nonlinearoperators/SoftPlus.jl @@ -19,14 +19,15 @@ end SoftPlus(DomainDim::NTuple{N,Int}) where {N} = SoftPlus{Float64,N}(DomainDim) -function A_mul_B!(y::AbstractArray{T,N}, L::SoftPlus{T,N}, x::AbstractArray{T,N}) where {T,N} - y .= log.(1.+exp.(x)) +function mul!(y::AbstractArray{T,N}, L::SoftPlus{T,N}, x::AbstractArray{T,N}) where {T,N} + y .= log.(1 .+exp.(x)) end -function Ac_mul_B!(y::AbstractArray{T,N}, - L::Jacobian{A}, - b::AbstractArray{T,N}) where {T,N, A <: SoftPlus{T,N}} - y .= 1./(1.+exp.(-L.x)).*b +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 end fun_name(L::SoftPlus) = "σ" diff --git a/src/nonlinearoperators/Tanh.jl b/src/nonlinearoperators/Tanh.jl index 18ffe88..ea5e99e 100644 --- a/src/nonlinearoperators/Tanh.jl +++ b/src/nonlinearoperators/Tanh.jl @@ -20,13 +20,14 @@ end Tanh(DomainDim::NTuple{N,Int}) where {N} = Tanh{Float64,N}(DomainDim) Tanh(DomainDim::Vararg{Int}) = Tanh{Float64,length(DomainDim)}(DomainDim) -function A_mul_B!(y::AbstractArray{T,N}, L::Tanh{T,N}, x::AbstractArray{T,N}) where {T,N} +function mul!(y::AbstractArray{T,N}, L::Tanh{T,N}, x::AbstractArray{T,N}) where {T,N} y .= tanh.(x) end -function Ac_mul_B!(y::AbstractArray{T,N}, - L::Jacobian{A}, - b::AbstractArray{T,N}) where {T,N, A<: Tanh{T,N}} +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 end diff --git a/src/properties.jl b/src/properties.jl index 9d0856f..bc87804 100644 --- a/src/properties.jl +++ b/src/properties.jl @@ -1,5 +1,6 @@ -import Base: size, ndims, diag +import Base: size, ndims +import LinearAlgebra: diag export ndoms, domainType, @@ -56,7 +57,6 @@ codomainType 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] """ @@ -65,12 +65,11 @@ size(L::AbstractOperator, i::Int) = size(L)[i] Returns a `Tuple` with the number of dimensions of the codomain and domain of an `AbstractOperator`. Type `ndims(A,1)` for the number of dimensions of the codomain and `ndims(A,2)` for the number of dimensions of the codomain. """ - ndims(L::AbstractOperator) = count_dims(size(L,1)), count_dims(size(L,2)) ndims(L::AbstractOperator, i::Int) = ndims(L)[i] -count_dims{N}(dims::NTuple{N,Int}) = N -count_dims{N}(dims::NTuple{N,Tuple}) = count_dims.(dims) +count_dims(dims::NTuple{N,Int}) where N = N +count_dims(dims::NTuple{N,Tuple}) where N = count_dims.(dims) """ `ndoms(L::AbstractOperator, [dom::Int]) -> (number of codomains, number of domains)` @@ -87,7 +86,7 @@ julia> ndoms(hcat(DFT(10,10),DFT(10,10))) julia> ndoms(hcat(DFT(10,10),DFT(10,10)),2) 2 -julia> ndoms(blkdiag(DFT(10,10),DFT(10,10)) +julia> ndoms(DCAT(DFT(10,10),DFT(10,10))) (2,2) ``` """ @@ -314,7 +313,7 @@ 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) @@ -324,7 +323,7 @@ function fun_space(L::AbstractOperator) end function fun_dom(L::AbstractOperator,n::Int) - dm = n == 2? domainType(L) : codomainType(L) + dm = n == 2 ? domainType(L) : codomainType(L) sz = size(L,n) return string_dom(dm,sz) end @@ -337,5 +336,5 @@ end function string_dom(dm::Tuple,sz::Tuple) s = string_dom.(dm,sz) - length(s) > 3 ? s[1]*"..."s[end] : *(s...) + length(s) > 3 ? s[1]*"..."*s[end] : *(s...) end diff --git a/src/syntax.jl b/src/syntax.jl index acc8584..7c367c0 100644 --- a/src/syntax.jl +++ b/src/syntax.jl @@ -1,32 +1,25 @@ -import Base: blkdiag, transpose, *, +, -, getindex, hcat, vcat, reshape +import Base: adjoint, *, +, -, getindex, hcat, vcat, reshape export jacobian -###### blkdiag ###### -blkdiag(L::Vararg{AbstractOperator}) = DCAT(L...) - ###### ' ###### -transpose{T <: AbstractOperator}(L::T) = Transpose(L) +adjoint(L::T) where {T <: AbstractOperator} = AdjointOperator(L) ######+,-###### -(+){T <: AbstractOperator}(L::T) = L -(-){T <: AbstractOperator}(L::T) = Scale(-1.0, L) +(+)(L::T) where {T <: AbstractOperator} = L +(-)(L::T) where {T <: AbstractOperator} = Scale(-1.0, L) (+)(L1::AbstractOperator, L2::AbstractOperator) = Sum(L1, L2 ) (-)(L1::AbstractOperator, L2::AbstractOperator) = Sum(L1, -L2 ) ###### * ###### -function (*){T <: BlockArray}(L::AbstractOperator, b::T) +function (*)(L::AbstractOperator, b::T) where {T <: BlockArray} y = blockzeros(codomainType(L), size(L, 1)) - A_mul_B!(y, L, b) + mul!(y, L, b) return y end -*{T<:Number}(coeff::T, L::AbstractOperator) = Scale(coeff,L) +*(coeff::T, L::AbstractOperator) where {T<:Number} = Scale(coeff,L) *(L1::AbstractOperator, L2::AbstractOperator) = Compose(L1,L2) -# redefine .* -Base.broadcast(::typeof(*), d::AbstractArray, L::AbstractOperator) = DiagOp(codomainType(L), size(d), d)*L -Base.broadcast(::typeof(*), d::AbstractArray, L::Scale) = DiagOp(L.coeff*d)*L.A - # getindex function getindex(A::AbstractOperator,idx...) if ndoms(A,2) == 1 @@ -40,7 +33,7 @@ function getindex(A::AbstractOperator,idx...) end #get index of HCAT returns HCAT (or Operator) -function getindex{M,N,L,P,C,A<:HCAT{M,N,L,P,C}}(H::A, idx::Union{AbstractArray,Int}) +function getindex(H::A, idx::Union{AbstractArray,Int}) where {M,N,L,P,C,A<:HCAT{M,N,L,P,C}} unfolded = vcat([[i... ] for i in H.idxs]...) if length(idx) == length(unfolded) @@ -64,7 +57,7 @@ end #get index of HCAT returns HCAT (or Operator) -function getindex{M,N,L,P,C,A<:VCAT{M,N,L,P,C}}(H::A, idx::Union{AbstractArray,Int}) +function getindex(H::A, idx::Union{AbstractArray,Int}) where {M,N,L,P,C,A<:VCAT{M,N,L,P,C}} unfolded = vcat([[i... ] for i in H.idxs]...) if length(idx) == length(unfolded) @@ -90,14 +83,14 @@ getindex(H::A, idx::Union{AbstractArray,Int}) where {L <: HCAT, D, S, A<: Affine AffineAdd(getindex(H.A, idx), H.d, S) # get index of scale -getindex{T, L, S <:Scale{T,L}}(A::S,idx...) = Scale(A.coeff,A.coeff_conj,getindex(A.A,idx...)) +getindex(A::S,idx...) where {T, L, S <:Scale{T,L}} = Scale(A.coeff,A.coeff_conj,getindex(A.A,idx...)) hcat(L::Vararg{AbstractOperator}) = HCAT(L...) vcat(L::Vararg{AbstractOperator}) = VCAT(L...) ###### reshape ###### -reshape{N,A<:AbstractOperator}(L::A, idx::NTuple{N,Int}) = Reshape(L,idx) -reshape{A<:AbstractOperator}(L::A, idx::Vararg{Int}) = Reshape(L,idx) +reshape(L::A, idx::NTuple{N,Int}) where {N,A<:AbstractOperator} = Reshape(L,idx) +reshape(L::A, idx::Vararg{Int}) where {A<:AbstractOperator} = Reshape(L,idx) ###### jacobian ###### jacobian = Jacobian diff --git a/src/utilities/block.jl b/src/utilities/block.jl index 41f88c9..083a156 100644 --- a/src/utilities/block.jl +++ b/src/utilities/block.jl @@ -1,6 +1,8 @@ # Define block-arrays module BlockArrays +using LinearAlgebra + export RealOrComplex, BlockArray, blocksize, @@ -16,6 +18,8 @@ export RealOrComplex, blockzeros, blockones, blockaxpy!, + blockscale!, + blockcumscale!, blockiszero @@ -37,7 +41,7 @@ blocklength(x::Tuple) = sum(blocklength.(x)) blocklength(x::AbstractArray) = length(x) blockvecnorm(x::Tuple) = sqrt(real(blockvecdot(x, x))) -blockvecnorm(x::AbstractArray{R}) where {R <: Number} = vecnorm(x) +blockvecnorm(x::AbstractArray{R}) where {R <: Number} = norm(x) blockmaxabs(x::Tuple) = maximum(blockmaxabs.(x)) blockmaxabs(x::AbstractArray{R}) where {R <: Number}= maximum(abs, x) @@ -49,30 +53,36 @@ blockcopy(x::Tuple) = blockcopy.(x) blockcopy(x::AbstractArray) = copy(x) blockcopy!(y::Tuple, x::Tuple) = blockcopy!.(y, x) -blockcopy!(y::AbstractArray, x::AbstractArray) = copy!(y, x) +blockcopy!(y::AbstractArray, x::AbstractArray) = copyto!(y, x) blockset!(y::Tuple, x) = blockset!.(y, x) blockset!(y::AbstractArray, x) = (y .= x) blockvecdot(x::T1, y::T2) where {T1 <: Tuple, T2 <: Tuple} = sum(blockvecdot.(x,y)) -blockvecdot(x::AbstractArray{R1}, y::AbstractArray{R2}) where {R1 <: Number, R2 <: Number} = vecdot(x, y) +blockvecdot(x::AbstractArray{R1}, y::AbstractArray{R2}) where {R1 <: Number, R2 <: Number} = dot(x, y) blockzeros(t::Tuple, s::Tuple) = blockzeros.(t, s) blockzeros(t::Type, n::NTuple{N, Integer} where {N}) = zeros(t, n) blockzeros(t::Tuple) = blockzeros.(t) blockzeros(n::NTuple{N, Integer} where {N}) = zeros(n) blockzeros(n::Integer) = zeros(n) -blockzeros(a::AbstractArray) = zeros(a) +blockzeros(a::AbstractArray) = fill!(similar(a),zero(eltype(a))) blockones(t::Tuple, s::Tuple) = blockones.(t, s) blockones(t::Type, n::NTuple{N, Integer} where {N}) = ones(t, n) blockones(t::Tuple) = blockones.(t) blockones(n::NTuple{N, Integer} where {N}) = ones(n) blockones(n::Integer) = ones(n) -blockones(a::AbstractArray) = ones(a) +blockones(a::AbstractArray) = fill!(a,one(eltype(a))) + +blockscale!(z::Tuple, alpha::Real, y::Tuple) = blockscale!.(z, alpha, y) +blockscale!(z::AbstractArray, alpha::Real, y::AbstractArray) = (z .= alpha.*y) + +blockcumscale!(z::Tuple, alpha::Real, y::Tuple) = blockcumscale!.(z, alpha, y) +blockcumscale!(z::AbstractArray, alpha::Real, y::AbstractArray) = (z .+= alpha.*y) -blockaxpy!(z::Tuple, x, alpha::Real, y::Tuple) = blockaxpy!.(z, x, alpha, y) -blockaxpy!(z::AbstractArray, x, alpha::Real, y::AbstractArray) = (z .= x .+ alpha.*y) +blockaxpy!(z::Tuple, x::Tuple, alpha::Real, y::Tuple) = blockaxpy!.(z, x, alpha, y) +blockaxpy!(z::AbstractArray, x::AbstractArray, alpha::Real, y::AbstractArray) = (z .= x .+ alpha.*y) blockiszero(x::AbstractArray) = iszero(x) blockiszero(x::Tuple) = all(iszero.(x)) diff --git a/test/runtests.jl b/test/runtests.jl index 55d9c01..ad484f8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,9 +1,12 @@ using AbstractOperators -using Base.Test -using Base.Profile +using LinearAlgebra, FFTW, DSP, SparseArrays +using Printf +using Random +using Test +using Profile include("utils.jl") -srand(0) +Random.seed!(0) verb = true diff --git a/test/test_block.jl b/test/test_block.jl index 656777f..9ed39c9 100644 --- a/test/test_block.jl +++ b/test/test_block.jl @@ -1,6 +1,8 @@ @printf("\nTesting BlockArrays\n") -using AbstractOperators:BlockArrays +using AbstractOperators.BlockArrays +using LinearAlgebra +using Printf T = Float64 x = [2.0; 3.0] @@ -18,8 +20,8 @@ x2b = (randn(2), randn(3)+im.*randn(3), randn(2,3)) @test blocklength(x) == 2 @test blocklength(xb) == 2+3+6 -@test blockvecnorm(x) == vecnorm(x) -@test blockvecnorm(xb) == sqrt(vecnorm(xb[1])^2+vecnorm(xb[2])^2+vecnorm(xb[3])^2) +@test blockvecnorm(x) == norm(x) +@test blockvecnorm(xb) == sqrt(norm(xb[1])^2+norm(xb[2])^2+norm(xb[3])^2) @test blockmaxabs(x) == 3.0 @test blockmaxabs(xb) == 6.0 @@ -48,8 +50,8 @@ blockset!(yb,xb) z = blockvecdot(x,x2) zb = blockvecdot(xb,x2b) -@test z == vecdot(x,x2) -@test zb == vecdot(xb[1],x2b[1])+vecdot(xb[2],x2b[2])+vecdot(xb[3],x2b[3]) +@test z == dot(x,x2) +@test zb == dot(xb[1],x2b[1])+dot(xb[2],x2b[2])+dot(xb[3],x2b[3]) y = blockzeros(x) yb = blockzeros(xb) @@ -81,11 +83,26 @@ yb = blockones(blockeltype(xb),blocksize(xb)) @test y == ones(2) @test yb == (ones(2),ones(3)+im*zeros(3),ones(2,3)) +blockscale!(y,2,x2) +blockscale!(yb,2,x2b) + +@test y == 2 .*x2 +@test yb == (2 .*x2b[1], 2 .*x2b[2], 2 .*x2b[3]) + +blockcopy!(y,x) +blockcopy!(yb,xb) + +blockcumscale!(y,2,x2) +blockcumscale!(yb,2,x2b) + +@test y == x .+ 2 .*x2 +@test yb == (xb[1] .+ 2 .*x2b[1], xb[2] .+ 2 .*x2b[2], xb[3] .+ 2 .*x2b[3]) + blockaxpy!(y,x,2,x2) blockaxpy!(yb,xb,2,x2b) -@test y == x .+ 2.*x2 -@test yb == (xb[1] .+ 2.*x2b[1], xb[2] .+ 2.*x2b[2], xb[3] .+ 2.*x2b[3]) +@test y == x .+ 2 .*x2 +@test yb == (xb[1] .+ 2 .*x2b[1], xb[2] .+ 2 .*x2b[2], xb[3] .+ 2 .*x2b[3]) x = (ones(Float64, 5), zeros(Float64, 2, 3)) @test blockiszero(x) == false diff --git a/test/test_lbfgs.jl b/test/test_lbfgs.jl index 25158fe..07d4c21 100644 --- a/test/test_lbfgs.jl +++ b/test/test_lbfgs.jl @@ -63,13 +63,13 @@ for i = 1:5 dir_ref = dirs_ref[:,i] gradm = -grad - @time A_mul_B!(dir, H, gradm) - @test vecnorm(dir-dir_ref, Inf)/(1+vecnorm(dir_ref, Inf)) <= 1e-15 + @time mul!(dir, H, gradm) + @test norm(dir-dir_ref, Inf)/(1+norm(dir_ref, Inf)) <= 1e-15 gradm2 = (-grad,-grad) - @time A_mul_B!(dirdir, HH, gradm2) - @test vecnorm(dirdir[1]-dir_ref, Inf)/(1+vecnorm(dir_ref, Inf)) <= 1e-15 - @test vecnorm(dirdir[2]-dir_ref, Inf)/(1+vecnorm(dir_ref, Inf)) <= 1e-15 + @time mul!(dirdir, HH, gradm2) + @test norm(dirdir[1]-dir_ref, Inf)/(1+norm(dir_ref, Inf)) <= 1e-15 + @test norm(dirdir[2]-dir_ref, Inf)/(1+norm(dir_ref, Inf)) <= 1e-15 x_old = x; grad_old = grad; @@ -81,8 +81,8 @@ end #testing reset -reset!(H) -reset!(HH) +AbstractOperators.reset!(H) +AbstractOperators.reset!(HH) @test blockones(size(H,1)) == H*blockones(size(H,1)) @test blockones(size(HH,1)) == HH*blockones(size(HH,1)) diff --git a/test/test_linear_operators.jl b/test/test_linear_operators.jl index 2ad418e..620edb1 100644 --- a/test/test_linear_operators.jl +++ b/test/test_linear_operators.jl @@ -1,6 +1,6 @@ @printf("\nTesting linear operators\n") -####### Conv ############ +######## Conv ############ n,m = 5, 6 h = randn(m) op = Conv(Float64,(n,),h) @@ -8,7 +8,7 @@ x1 = randn(n) y1 = test_op(op, x1, randn(n+m-1), verb) y2 = conv(x1,h) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) # other constructors op = Conv(x1,h) @@ -25,14 +25,14 @@ op = Conv(x1,h) @test is_full_row_rank(op) == true @test is_full_column_rank(op) == true -######### DCT ############ +########## DCT ############ n = 4 op = DCT(Float64,(n,)) x1 = randn(n) y1 = test_op(op, x1, randn(n), verb) y2 = dct(x1) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) # other constructors op = DCT((n,)) @@ -55,7 +55,7 @@ m = 10 op = DCT(n,m) x1 = randn(n,m) -@test vecnorm(op'*(op*x1) - x1) <= 1e-12 +@test norm(op'*(op*x1) - x1) <= 1e-12 @test diag_AAc(op) == 1. @test diag_AcA(op) == 1. @@ -66,7 +66,7 @@ x1 = randn(n) y1 = test_op(op, x1, randn(n), verb) y2 = idct(x1) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) # other constructors op = IDCT((n,)) @@ -89,30 +89,31 @@ m = 10 op = IDCT(n,m) x1 = randn(n,m) -@test vecnorm(op'*(op*x1) - x1) <= 1e-12 +@test norm(op'*(op*x1) - x1) <= 1e-12 @test diag_AAc(op) == 1. @test diag_AcA(op) == 1. ######## DFT ############ +# seems like there is an object called DFT in Base julia 0.7 (however in 1.0 was rm) n = 4 -op = DFT(Float64,(n,)) +op = AbstractOperators.DFT(Float64,(n,)) x1 = randn(n) y1 = test_op(op, x1, fft(randn(n)), verb) y2 = fft(x1) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) -op = DFT(Complex{Float64},(n,)) +op = AbstractOperators.DFT(Complex{Float64},(n,)) x1 = randn(n)+im*randn(n) y1 = test_op(op, x1, fft(randn(n)), verb) y2 = fft(x1) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) # other constructors -op = DFT((n,)) -op = DFT(n,n) -op = DFT(Complex{Float64}, n,n) +op = AbstractOperators.DFT((n,)) +op = AbstractOperators.DFT(n,n) +op = AbstractOperators.DFT(Complex{Float64}, n,n) #properties @test is_linear(op) == true @@ -127,11 +128,11 @@ op = DFT(Complex{Float64}, n,n) @test is_full_column_rank(op) == true m = 10 -op = DFT(n,m) +op = AbstractOperators.DFT(n,m) x1 = randn(n,m) y1 = op*x1 -@test vecnorm(op'*(op*x1) - diag_AcA(op)*x1) <= 1e-12 -@test vecnorm(op*(op'*y1) - diag_AAc(op)*y1) <= 1e-12 +@test norm(op'*(op*x1) - diag_AcA(op)*x1) <= 1e-12 +@test norm(op*(op'*y1) - diag_AAc(op)*y1) <= 1e-12 ######### IDFT ############ n = 4 @@ -140,14 +141,14 @@ x1 = randn(n) y1 = test_op(op, x1, fft(randn(n)), verb) y2 = ifft(x1) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) op = IDFT(Complex{Float64},(n,)) x1 = randn(n)+im*randn(n) y1 = test_op(op, x1, fft(randn(n)), verb) y2 = ifft(x1) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) # other constructors op = IDFT((n,)) @@ -170,8 +171,8 @@ m = 10 op = IDFT(n,m) x1 = randn(n,m) y1 = op*x1 -@test vecnorm(op'*(op*x1) - diag_AcA(op)*x1) <= 1e-12 -@test vecnorm(op*(op'*y1) - diag_AAc(op)*y1) <= 1e-12 +@test norm(op'*(op*x1) - diag_AcA(op)*x1) <= 1e-12 +@test norm(op*(op'*y1) - diag_AAc(op)*y1) <= 1e-12 ####### RDFT ############ n = 4 @@ -180,7 +181,7 @@ x1 = randn(n) y1 = test_op(op, x1, rfft(x1), verb) y2 = rfft(x1) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) n,m,l = 4,8,5 op = RDFT(Float64,(n,m,l),2) @@ -188,7 +189,7 @@ x1 = randn(n,m,l) y1 = test_op(op, x1, rfft(x1,2), verb) y2 = rfft(x1,2) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) # other constructors op = RDFT((n,)) @@ -213,7 +214,7 @@ x1 = rfft(randn(n)) y1 = test_op(op, x1,irfft(randn(div(n,2)+1),n), verb) y2 = irfft(x1,n) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) n = 11 op = IRDFT(Complex{Float64},(div(n,2)+1,),n) @@ -221,7 +222,7 @@ x1 = rfft(randn(n)) y1 = test_op(op, x1,irfft(randn(div(n,2)+1),n), verb) y2 = irfft(x1,n) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) n,m,l = 4,19,5 op = IRDFT(Complex{Float64},(n,div(m,2)+1,l),m,2) @@ -229,7 +230,7 @@ x1 = rfft(randn(n,m,l),2) y1 = test_op(op, x1, irfft(x1,m,2), verb) y2 = irfft(x1,m,2) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) n,m,l = 4,18,5 op = IRDFT(Complex{Float64},(n,div(m,2)+1,l),m,2) @@ -237,7 +238,7 @@ x1 = rfft(randn(n,m,l),2) y1 = test_op(op, x1, irfft(x1,m,2), verb) y2 = irfft(x1,m,2) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) n,m,l = 5,18,5 op = IRDFT(Complex{Float64},(div(n,2)+1,m,l),n,1) @@ -245,7 +246,7 @@ x1 = rfft(randn(n,m,l),1) y1 = test_op(op, x1, irfft(x1,n,1), verb) y2 = irfft(x1,n,1) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) ## other constructors op = IRDFT((10,),19) @@ -262,7 +263,7 @@ op = IRDFT((10,),19) @test is_full_row_rank(op) == true @test is_full_column_rank(op) == false -######### DiagOp ############ +########## DiagOp ############ n = 4 d = randn(n) op = DiagOp(Float64,(n,),d) @@ -270,16 +271,16 @@ x1 = randn(n) y1 = test_op(op, x1, randn(n), verb) y2 = d.*x1 -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) n = 4 d = randn(n)+im*randn(n) op = DiagOp(Float64,(n,),d) x1 = randn(n) -y1 = test_op(op, x1, randn(n)+im*randn(n), verb) +y1 = test_op(op, x1, randn(n).+im*randn(n), verb) y2 = d.*x1 -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) n = 4 d = pi @@ -288,7 +289,7 @@ x1 = randn(n) y1 = test_op(op, x1, randn(n), verb) y2 = d.*x1 -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) n = 4 d = im @@ -297,13 +298,13 @@ x1 = randn(n) y1 = test_op(op, x1, randn(n)+im*randn(n), verb) y2 = d.*x1 -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) # other constructors d = randn(4) op = DiagOp(d) -d = randn(4)+im +d = randn(4).+im op = DiagOp(d) n = 4 @@ -326,8 +327,8 @@ op = DiagOp((n,), d) @test is_full_column_rank(DiagOp([ones(5);0])) == false @test diag(op) == d -@test vecnorm(op'*(op*x1) - diag_AcA(op).*x1) <= 1e-12 -@test vecnorm(op*(op'*x1) - diag_AAc(op).*x1) <= 1e-12 +@test norm(op'*(op*x1) .- diag_AcA(op).*x1) <= 1e-12 +@test norm(op*(op'*x1) .- diag_AAc(op).*x1) <= 1e-12 n = 4 d = pi @@ -335,8 +336,8 @@ op = DiagOp((n,), d) x1 = randn(n) @test diag(op) == d -@test vecnorm(op'*(op*x1) - diag_AcA(op).*x1) <= 1e-12 -@test vecnorm(op*(op'*x1) - diag_AAc(op).*x1) <= 1e-12 +@test norm(op'*(op*x1) .- diag_AcA(op).*x1) <= 1e-12 +@test norm(op*(op'*x1) .- diag_AAc(op).*x1) <= 1e-12 ########## Eye ############ n = 4 @@ -344,7 +345,7 @@ op = Eye(Float64,(n,)) x1 = randn(n) y1 = test_op(op, x1, randn(n), verb) -@test all(vecnorm.(y1 .- x1) .<= 1e-12) +@test all(norm.(y1 .- x1) .<= 1e-12) # other constructors op = Eye(Float64, (n,)) @@ -376,7 +377,7 @@ x1 = randn(n) y1 = test_op(op, x1, randn(n), verb) y2 = filt(b, a, x1) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) h = randn(10) op = Filt(Float64,(n,m),h) @@ -384,7 +385,7 @@ x1 = randn(n,m) y1 = test_op(op, x1, randn(n,m), verb) y2 = filt(h, [1.], x1) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) # other constructors Filt(n, b, a) @@ -406,67 +407,70 @@ Filt(x1, b) @test is_full_row_rank(op) == true @test is_full_column_rank(op) == true -####### FiniteDiff ############ +###### FiniteDiff ############ n= 10 op = FiniteDiff(Float64,(n,)) x1 = randn(n) y1 = test_op(op, x1, randn(n-1), verb) -y1 = op*collect(linspace(0,1,n)) -@test all(vecnorm.(y1 .- 1/9) .<= 1e-12) +y1 = op*collect(range(0,stop=1,length=n)) +@test all(norm.(y1 .- 1/9) .<= 1e-12) + +I1, J1, V1 = SparseArrays.spdiagm_internal(0 => ones(n-1)) +I2, J2, V2 = SparseArrays.spdiagm_internal(1 => ones(n-1)) +B = -sparse(I1,J1,V1,n-1,n)+sparse(I2,J2,V2,n-1,n) -B = -spdiagm(ones(n-1),0,n-1,n)+spdiagm(ones(n-1),1,n-1,n) @test norm(B*x1-op*x1) <= 1e-8 n,m= 10,5 op = FiniteDiff(Float64,(n,m)) x1 = randn(n,m) y1 = test_op(op, x1, randn(n-1,m), verb) -y1 = op*repmat(collect(linspace(0,1,n)),1,m) -@test all(vecnorm.(y1 .- 1/9) .<= 1e-12) +y1 = op*repeat(collect(range(0,stop=1,length=n)),1,m) +@test all(norm.(y1 .- 1/9) .<= 1e-12) n,m= 10,5 op = FiniteDiff(Float64,(n,m),2) x1 = randn(n,m) y1 = test_op(op, x1, randn(n,m-1), verb) -y1 = op*repmat(collect(linspace(0,1,n)),1,m) -@test all(vecnorm.(y1) .<= 1e-12) +y1 = op*repeat(collect(range(0,stop=1,length=n)),1,m) +@test all(norm.(y1) .<= 1e-12) n,m,l= 10,5,7 op = FiniteDiff(Float64,(n,m,l)) x1 = randn(n,m,l) y1 = test_op(op, x1, randn(n-1,m,l), verb) -y1 = op*reshape(repmat(collect(linspace(0,1,n)),1,m*l),n,m,l) -@test all(vecnorm.(y1 .- 1/9) .<= 1e-12) +y1 = op*reshape(repeat(collect(range(0,stop=1,length=n)),1,m*l),n,m,l) +@test all(norm.(y1 .- 1/9) .<= 1e-12) n,m,l= 10,5,7 op = FiniteDiff(Float64,(n,m,l),2) x1 = randn(n,m,l) y1 = test_op(op, x1, randn(n,m-1,l), verb) -y1 = op*reshape(repmat(collect(linspace(0,1,n)),1,m*l),n,m,l) -@test all(vecnorm.(y1) .<= 1e-12) +y1 = op*reshape(repeat(collect(range(0,stop=1,length=n)),1,m*l),n,m,l) +@test all(norm.(y1) .<= 1e-12) n,m,l= 10,5,7 op = FiniteDiff(Float64,(n,m,l),3) x1 = randn(n,m,l) y1 = test_op(op, x1, randn(n,m,l-1), verb) -y1 = op*reshape(repmat(collect(linspace(0,1,n)),1,m*l),n,m,l) -@test all(vecnorm.(y1) .<= 1e-12) +y1 = op*reshape(repeat(collect(range(0,stop=1,length=n)),1,m*l),n,m,l) +@test all(norm.(y1) .<= 1e-12) n,m,l,i = 5,6,2,3 op = FiniteDiff(Float64,(n,m,l,i),1) x1 = randn(n,m,l,i) y1 = test_op(op, x1, randn(n-1,m,l,i), verb) -y1 = op*reshape(repmat(collect(linspace(0,1,n)),1,m*l*i),n,m,l,i) -@test all(vecnorm.(y1 .- 1/(n-1)) .<= 1e-12) +y1 = op*reshape(repeat(collect(range(0,stop=1,length=n)),1,m*l*i),n,m,l,i) +@test all(norm.(y1 .- 1/(n-1)) .<= 1e-12) n,m,l,i = 5,6,2,3 op = FiniteDiff(Float64,(n,m,l,i),4) x1 = randn(n,m,l,i) y1 = test_op(op, x1, randn(n,m,l,i-1), verb) -y1 = op*reshape(repmat(collect(linspace(0,1,n)),1,m*l*i),n,m,l,i) -@test vecnorm(y1) <= 1e-12 +y1 = op*reshape(repeat(collect(range(0,stop=1,length=n)),1,m*l*i),n,m,l,i) +@test norm(y1) <= 1e-12 -@test_throws ErrorException op = FiniteDiff(Float64,(n,m,l), 4) +@test_throws ErrorException FiniteDiff(Float64,(n,m,l), 4) ## other constructors FiniteDiff((n,m)) @@ -491,7 +495,7 @@ op = GetIndex(Float64,(n,),(1:k,)) x1 = randn(n) y1 = test_op(op, x1, randn(k), verb) -@test all(vecnorm.(y1 .- x1[1:k]) .<= 1e-12) +@test all(norm.(y1 .- x1[1:k]) .<= 1e-12) n,m = 5,4 k = 3 @@ -499,27 +503,27 @@ op = GetIndex(Float64,(n,m),(1:k,:)) x1 = randn(n,m) y1 = test_op(op, x1, randn(k,m), verb) -@test all(vecnorm.(y1 .- x1[1:k,:]) .<= 1e-12) +@test all(norm.(y1 .- x1[1:k,:]) .<= 1e-12) n,m = 5,4 op = GetIndex(Float64,(n,m),(:,2)) x1 = randn(n,m) y1 = test_op(op, x1, randn(n), verb) -@test all(vecnorm.(y1 .- x1[:,2]) .<= 1e-12) +@test all(norm.(y1 .- x1[:,2]) .<= 1e-12) n,m,l = 5,4,3 op = GetIndex(Float64,(n,m,l),(1:3,2,:)) x1 = randn(n,m,l) y1 = test_op(op, x1, randn(3,3), verb) -@test all(vecnorm.(y1 .- x1[1:3,2,:]) .<= 1e-12) +@test all(norm.(y1 .- x1[1:3,2,:]) .<= 1e-12) # other constructors GetIndex((n,m), (1:k,:)) GetIndex(x1, (1:k,:)) -@test_throws ErrorException op = GetIndex(Float64,(n,m),(1:k,:,:)) +@test_throws ErrorException GetIndex(Float64,(n,m),(1:k,:,:)) op = GetIndex(Float64,(n,m),(1:n,1:m)) @test typeof(op) <: Eye @@ -540,7 +544,7 @@ op = GetIndex(Float64,(n,),(1:k,)) @test diag_AAc(op) == 1 -######## MatrixOp ############ +####### MatrixOp ############ # real matrix, real input n,m = 5,4 @@ -574,13 +578,13 @@ x1 = randn(m) y1 = test_op(op, x1, randn(n)+im*randn(n), verb) y2 = A*x1 -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) # complex matrix, real matrix input c = 3 op = MatrixOp(Float64,(m,c),A) -@test_throws ErrorException op = MatrixOp(Float64,(m,c,3),A) -@test_throws MethodError op = MatrixOp(Float64,(m,c),randn(n,m,2)) +@test_throws ErrorException MatrixOp(Float64,(m,c,3),A) +@test_throws MethodError MatrixOp(Float64,(m,c),randn(n,m,2)) x1 = randn(m,c) y1 = test_op(op, x1, randn(n,c).+randn(n,c), verb) y2 = A*x1 @@ -604,11 +608,11 @@ op = convert(LinearOperator, Complex{Float64}, size(x1), A) @test is_AAc_diagonal(op) == false @test is_orthogonal(op) == false @test is_invertible(op) == false -@test is_full_row_rank(MatrixOp(randn(srand(0),3,4))) == true -@test is_full_column_rank(MatrixOp(randn(srand(0),3,4))) == false +@test is_full_row_rank(MatrixOp(randn(Random.seed!(0),3,4))) == true +@test is_full_column_rank(MatrixOp(randn(Random.seed!(0),3,4))) == false -###### LMatrixOp ############ +##### LMatrixOp ############ n,m = 5,6 b = randn(m) @@ -617,7 +621,7 @@ x1 = randn(n,m) y1 = test_op(op, x1, randn(n), verb) y2 = x1*b -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) n,m = 5,6 b = randn(m)+im*randn(m) @@ -626,7 +630,7 @@ x1 = randn(n,m)+im*randn(n,m) y1 = test_op(op, x1, randn(n)+im*randn(n), verb) y2 = x1*b -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) n,m,l = 5,6,7 b = randn(m,l) @@ -635,7 +639,7 @@ x1 = randn(n,m) y1 = test_op(op, x1, randn(n,l), verb) y2 = x1*b -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) n,m,l = 5,6,7 b = randn(m,l)+im*randn(m,l) @@ -644,7 +648,7 @@ x1 = randn(n,m)+im*randn(n,m) y1 = test_op(op, x1, randn(n,l)+im*randn(n,l), verb) y2 = x1*b -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) ## other constructors op = LMatrixOp(b,n) @@ -665,13 +669,13 @@ op = LMatrixOp(b,n) n,m = 5,4; A = randn(n,m); -op = MyLinOp(Float64, (m,),(n,), (y,x) -> A_mul_B!(y,A,x), (y,x) -> Ac_mul_B!(y,A,x)) +op = MyLinOp(Float64, (m,),(n,), (y,x) -> mul!(y,A,x), (y,x) -> mul!(y,A',x)) x1 = randn(m) y1 = test_op(op, x1, randn(n), verb) y2 = A*x1 # other constructors -op = MyLinOp(Float64, (m,), Float64, (n,), (y,x) -> A_mul_B!(y,A,x), (y,x) -> Ac_mul_B!(y,A,x)) +op = MyLinOp(Float64, (m,), Float64, (n,), (y,x) -> mul!(y,A,x), (y,x) -> mul!(y,A',x)) ######### MIMOFilt ############ @@ -684,7 +688,7 @@ x1 = randn(m,n) y1 = test_op(op, x1, randn(m,1), verb) y2 = filt(b[1],a[1],x1[:,1])+filt(b[2],a[2],x1[:,2]) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) m,n = 10,3; #time samples, number of inputs b = [[1.;0.;1.],[1.;0.;1.],[1.;0.;1.],[1.;0.;1.],[1.;0.;1.],[1.;0.;1.], ]; @@ -695,7 +699,7 @@ x1 = randn(m,n) y1 = test_op(op, x1, randn(m,2), verb) y2 = [filt(b[1],a[1],x1[:,1])+filt(b[2],a[2],x1[:,2])+filt(b[3],a[3],x1[:,3]) filt(b[4],a[4],x1[:,1])+filt(b[5],a[5],x1[:,2])+filt(b[6],a[6],x1[:,3])] -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) m,n = 10,3 b = [randn(10),randn(5),randn(10),randn(2),randn(10),randn(10)] @@ -706,7 +710,7 @@ x1 = randn(m,n) y1 = test_op(op, x1, randn(m,2), verb) y2 = [filt(b[1],a[1],x1[:,1])+filt(b[2],a[2],x1[:,2])+filt(b[3],a[3],x1[:,3]) filt(b[4],a[4],x1[:,1])+filt(b[5],a[5],x1[:,2])+filt(b[6],a[6],x1[:,3])] -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) ## other constructors MIMOFilt((10,3), b, a) @@ -715,16 +719,16 @@ MIMOFilt(x1, b, a) MIMOFilt(x1, b) #errors -@test_throws ErrorException op = MIMOFilt(Float64, (10,3,2) ,b,a) +@test_throws ErrorException MIMOFilt(Float64, (10,3,2) ,b,a) a2 = [[1.0f0],[1.0f0],[1.0f0],[1.0f0],[1.0f0],[1.0f0]] b2 = convert.(Array{Float32,1},b) -@test_throws ErrorException op = MIMOFilt(Float64, (m,n),b2,a2) -@test_throws ErrorException op = MIMOFilt(Float64, (m,n),b,a[1:end-1]) +@test_throws ErrorException MIMOFilt(Float64, (m,n),b2,a2) +@test_throws ErrorException MIMOFilt(Float64, (m,n),b,a[1:end-1]) push!(a2,[1.0f0]) push!(b2,randn(Float32,10)) -@test_throws ErrorException op = MIMOFilt(Float32, (m,n),b2,a2) +@test_throws ErrorException MIMOFilt(Float32, (m,n),b2,a2) a[1][1] = 0. -@test_throws ErrorException op = MIMOFilt(Float64, (m,n) ,b,a) +@test_throws ErrorException MIMOFilt(Float64, (m,n) ,b,a) b = [randn(10),randn(5),randn(10),randn(2),randn(10),randn(10)] a = [[1.],[1.],[1.],[1.],[1.],[1.]] @@ -749,31 +753,31 @@ op = Variation(Float64,(n,m)) x1 = randn(n,m) y1 = test_op(op, x1, randn(m*n,2), verb) -y1 = op*repmat(collect(linspace(0,1,n)),1,m) -@test all(vecnorm.(y1[:,1] .- 1/(n-1) ) .<= 1e-12) -@test all(vecnorm.(y1[:,2] ) .<= 1e-12) +y1 = op*repeat(collect(range(0,stop=1,length=n)),1,m) +@test all(norm.(y1[:,1] .- 1/(n-1) ) .<= 1e-12) +@test all(norm.(y1[:,2] ) .<= 1e-12) -Dx = spdiagm(ones(n),0,n,n)-spdiagm(ones(n-1),-1,n,n) +Dx = spdiagm(0 => ones(n), -1 => -ones(n-1)) Dx[1,1],Dx[1,2] = -1,1 -Dy = spdiagm(ones(m),0,m,m)-spdiagm(ones(m-1),-1,m,m) +Dy = spdiagm(0 => ones(m), -1 => -ones(m-1)) Dy[1,1],Dy[1,2] = -1,1 -Dxx = kron(eye(m),Dx) -Dyy = kron(Dy,eye(n)) +Dxx = kron(sparse(I,m,m),Dx) +Dyy = kron(Dy,sparse(I,n,n)) TV = [Dxx;Dyy] x1 = randn(n,m) -@test vecnorm(op*x1-reshape(TV*(x1[:]),n*m,2))<1e-12 +@test norm(op*x1-reshape(TV*(x1[:]),n*m,2))<1e-12 n,m,l = 10,5,3 op = Variation(Float64,(n,m,l)) x1 = randn(n,m,l) y1 = test_op(op, x1, randn(m*n*l,3), verb) -y1 = op*reshape(repmat(collect(linspace(0,1,n)),1,m*l),n,m,l) -@test all(vecnorm.(y1[:,1] .- 1/(n-1) ) .<= 1e-12) -@test all(vecnorm.(y1[:,2] ) .<= 1e-12) -@test all(vecnorm.(y1[:,3] ) .<= 1e-12) +y1 = op*reshape(repeat(collect(range(0,stop=1,length=n)),1,m*l),n,m,l) +@test all(norm.(y1[:,1] .- 1/(n-1) ) .<= 1e-12) +@test all(norm.(y1[:,2] ) .<= 1e-12) +@test all(norm.(y1[:,3] ) .<= 1e-12) ### other constructors Variation(Float64, n,m) @@ -782,7 +786,7 @@ Variation(n,m) Variation(x1) ##errors -@test_throws ErrorException op = Variation(Float64,(n,)) +@test_throws ErrorException Variation(Float64,(n,)) ###properties @test is_linear(op) == true @@ -804,7 +808,7 @@ x1 = randn(n) y1 = test_op(op, x1, randn(n+m), verb) y2 = xcorr(x1, h) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) # other constructors op = Xcorr(x1,h) @@ -826,7 +830,7 @@ z = (5,) op = ZeroPad(Float64,n,z) x1 = randn(n) y1 = test_op(op, x1, randn(n.+z), verb) -@test all(vecnorm.(y1 .- [x1;zeros(5)] ) .<= 1e-12) +@test all(norm.(y1 .- [x1;zeros(5)] ) .<= 1e-12) n = (3,2) z = (5,3) @@ -835,7 +839,7 @@ x1 = randn(n) y1 = test_op(op, x1, randn(n.+z), verb) y2 = zeros(n.+z) y2[1:n[1],1:n[2]] = x1 -@test all(vecnorm.(y1 .- y2 ) .<= 1e-12) +@test all(norm.(y1 .- y2 ) .<= 1e-12) n = (3,2,2) z = (5,3,1) @@ -844,7 +848,7 @@ x1 = randn(n) y1 = test_op(op, x1, randn(n.+z), verb) y2 = zeros(n.+z) y2[1:n[1],1:n[2],1:n[3]] = x1 -@test all(vecnorm.(y1 .- y2 ) .<= 1e-12) +@test all(norm.(y1 .- y2 ) .<= 1e-12) # other constructors ZeroPad(n, z...) @@ -854,8 +858,8 @@ ZeroPad(x1, z) ZeroPad(x1, z...) #errors -@test_throws ErrorException op = ZeroPad(Float64,n,(1,2)) -@test_throws ErrorException op = ZeroPad(Float64,n,(1,-2,3)) +@test_throws ErrorException ZeroPad(Float64,n,(1,2)) +@test_throws ErrorException ZeroPad(Float64,n,(1,-2,3)) #properties @test is_linear(op) == true diff --git a/test/test_linear_operators_calculus.jl b/test/test_linear_operators_calculus.jl index 80d2d1c..78cc116 100644 --- a/test/test_linear_operators_calculus.jl +++ b/test/test_linear_operators_calculus.jl @@ -13,7 +13,7 @@ opC = Compose(opA2,opA1) x = randn(m1) y1 = test_op(opC, x, randn(m3), verb) y2 = A2*A1*x -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) # test Compose longer m1, m2, m3, m4 = 4, 7, 3, 2 @@ -30,8 +30,8 @@ x = randn(m1) y1 = test_op(opC1, x, randn(m4), verb) y2 = test_op(opC2, x, randn(m4), verb) y3 = A3*A2*A1*x -@test all(vecnorm.(y1 .- y2) .<= 1e-12) -@test all(vecnorm.(y3 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y3 .- y2) .<= 1e-12) #test Compose special cases @test typeof(opA1*Eye(m1)) == typeof(opA1) @@ -86,13 +86,13 @@ opC = Compose(Compose(Compose(opA4,opA3),opA2),opA1) x = randn(m1) -@test norm( opC*x - (A4*(A3*( A2*( A1*x+d1 )+d2 )+d3)-d4) ) < 1e-9 -@test norm( displacement(opC) - ( A4*(A3*(A2*d1+d2)+d3)-d4) ) < 1e-9 +@test norm( opC*x - (A4*(A3*( A2*( A1*x+d1 ) .+ d2 ) .+ d3)-d4) ) < 1e-9 +@test norm( displacement(opC) - ( A4*(A3*(A2*d1 .+d2 ) .+ d3)-d4) ) < 1e-9 opA4 = MatrixOp(A4) opC = AffineAdd(Compose(Compose(Compose(opA4,opA3),opA2),opA1),d4) -@test norm( opC*x - (A4*(A3*( A2*( A1*x+d1 )+d2 )+d3)+d4) ) < 1e-9 -@test norm( displacement(opC) - ( A4*(A3*(A2*d1+d2)+d3)+d4) ) < 1e-9 +@test norm( opC*x - (A4*(A3*( A2*( A1*x+d1 ) .+ d2 ) .+ d3)+d4) ) < 1e-9 +@test norm( displacement(opC) - ( A4*(A3*(A2*d1 .+ d2) .+ d3)+d4) ) < 1e-9 @test norm( remove_displacement(opC)*x - (A4*(A3*( A2*( A1*x ) ) ) ) ) < 1e-9 @@ -110,7 +110,7 @@ x1 = randn(n1) x2 = randn(n2) y1 = test_op(opD, (x1, x2), (randn(m1),randn(m2)), verb) y2 = (A1*x1, A2*x2) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) # test DCAT longer @@ -127,7 +127,7 @@ x2 = randn(n2) x3 = randn(n3) y1 = test_op(opD, (x1, x2, x3), (randn(m1),randn(m2),randn(m3)), verb) y2 = (A1*x1, A2*x2, A3*x3) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) #properties @test is_linear(opD) == true @@ -184,15 +184,15 @@ x2 = randn(n2) x3 = randn(n3) y1 = opD*(x1,x2,x3) y2 = (A1*x1+d1, A2*x2+d2, A3*x3+d3) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) -@test all(vecnorm.(displacement(opD) .- (d1,d2,d3)) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) +@test all(norm.(displacement(opD) .- (d1,d2,d3)) .<= 1e-12) y1 = remove_displacement(opD)*(x1,x2,x3) y2 = (A1*x1, A2*x2, A3*x3) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) -######################## -### test HCAT ####### -######################## +####################### +## test HCAT ####### +####################### m, n1, n2 = 4, 7, 5 A1 = randn(m, n1) @@ -204,13 +204,13 @@ x1 = randn(n1) x2 = randn(n2) y1 = test_op(opH, (x1, x2), randn(m), verb) y2 = A1*x1 + A2*x2 -@test vecnorm(y1-y2) <= 1e-12 +@test norm(y1-y2) <= 1e-12 #permuatation p = [2;1] opHp = opH[p] y1 = test_op(opHp, (x2, x1), randn(m), verb) -@test vecnorm(y1-y2) <= 1e-12 +@test norm(y1-y2) <= 1e-12 # test HCAT longer @@ -227,37 +227,37 @@ x2 = randn(n2) x3 = randn(n3) y1 = test_op(opH, (x1, x2, x3), randn(m), verb) y2 = A1*x1 + A2*x2 + A3*x3 -@test vecnorm(y1-y2) <= 1e-12 +@test norm(y1-y2) <= 1e-12 # test HCAT of HCAT opHH = HCAT(opH, opA2, opA3) y1 = test_op(opHH, (x1, x2, x3, x2, x3), randn(m), verb) y2 = A1*x1 + A2*x2 + A3*x3 + A2*x2 + A3*x3 -@test vecnorm(y1-y2) <= 1e-12 +@test norm(y1-y2) <= 1e-12 opHH = HCAT(opH, opH, opA3) x = (x1, x2, x3, x1, x2, x3, x3) y1 = test_op(opHH, x, randn(m), verb) y2 = A1*x1 + A2*x2 + A3*x3 + A1*x1 + A2*x2 + A3*x3 + A3*x3 -@test vecnorm(y1-y2) <= 1e-12 +@test norm(y1-y2) <= 1e-12 opA3 = MatrixOp(randn(n1,n1)) @test_throws Exception HCAT(opA1,opA2,opA3) -opF = DFT(Complex{Float64},(m,)) +opF = AbstractOperators.DFT(Complex{Float64},(m,)) @test_throws Exception HCAT(opA1,opF,opA2) # test utilities # permutation p = randperm(ndoms(opHH,2)) -opHP = permute(opHH,p) +opHP = AbstractOperators.permute(opHH,p) xp = x[p] y1 = test_op(opHP, xp, randn(m), verb) pp = randperm(ndoms(opHH,2)) -opHPP = permute(opHH,pp) +opHPP = AbstractOperators.permute(opHH,pp) xpp = x[pp] y1 = test_op(opHPP, xpp, randn(m), verb) @@ -281,8 +281,8 @@ op = HCAT(opA1, opA2, opA3) @test is_full_row_rank(op) == true @test is_full_column_rank(op) == false -d = randn(n1)+im*randn(n1) -op = HCAT(DiagOp(d), DFT(Complex{Float64},n1)) +d = randn(n1).+im.*randn(n1) +op = HCAT(DiagOp(d), AbstractOperators.DFT(Complex{Float64},n1)) @test is_null(op) == false @test is_eye(op) == false @test is_diagonal(op) == false @@ -293,10 +293,10 @@ op = HCAT(DiagOp(d), DFT(Complex{Float64},n1)) @test is_full_row_rank(op) == true @test is_full_column_rank(op) == false -@test diag_AAc(op) == d.*conj(d)+n1 +@test diag_AAc(op) == d .* conj(d) .+ n1 -y1 = randn(n1)+im*randn(n1) -@test norm(op*(op'*y1)-diag_AAc(op).*y1) <1e-12 +y1 = randn(n1) .+ im .* randn(n1) +@test norm(op*(op'*y1).-diag_AAc(op).*y1) <1e-12 #test displacement @@ -312,14 +312,14 @@ x1 = randn(n1) x2 = randn(n2) y1 = opH*(x1,x2) y2 = A1*x1+d1 + A2*x2+d2 -@test vecnorm(y1-y2) <= 1e-12 +@test norm(y1-y2) <= 1e-12 y1 = remove_displacement(opH)*(x1,x2) y2 = A1*x1 + A2*x2 -@test vecnorm(y1-y2) <= 1e-12 +@test norm(y1-y2) <= 1e-12 -######################### -#### test Reshape ####### -######################### +######################## +### test Reshape ####### +######################## m, n = 8, 4 dim_out = (2, 2, 2) @@ -330,7 +330,7 @@ opR = Reshape(opA1, dim_out...) x1 = randn(n) y1 = test_op(opR, x1, randn(dim_out), verb) y2 = reshape(A1*x1, dim_out) -@test vecnorm(y1-y2) <= 1e-12 +@test norm(y1-y2) <= 1e-12 @test_throws Exception Reshape(opA1,(2,2,1)) @@ -354,10 +354,10 @@ opR = Reshape(opA1, dim_out) x1 = randn(n) y1 = opR*x1 y2 = reshape(A1*x1+d1, dim_out) -@test vecnorm(y1-y2) <= 1e-12 +@test norm(y1-y2) <= 1e-12 y1 = remove_displacement(opR)*x1 y2 = reshape(A1*x1, dim_out) -@test vecnorm(y1-y2) <= 1e-12 +@test norm(y1-y2) <= 1e-12 ####################### ## test Scale ####### @@ -371,20 +371,20 @@ opS = Scale(coeff, opA1) x1 = randn(n) y1 = test_op(opS, x1, randn(m), verb) y2 = coeff*A1*x1 -@test vecnorm(y1-y2) <= 1e-12 +@test norm(y1-y2) <= 1e-12 coeff2 = 3 opS2 = Scale(coeff2, opS) y1 = test_op(opS2, x1, randn(m), verb) y2 = coeff2*coeff*A1*x1 -@test vecnorm(y1-y2) <= 1e-12 +@test norm(y1-y2) <= 1e-12 -opF = DFT(m,n) +opF = AbstractOperators.DFT(m,n) opS = Scale(coeff, opF) x1 = randn(m,n) y1 = test_op(opS, x1, fft(randn(m,n)), verb) y2 = coeff*(fft(x1)) -@test vecnorm(y1-y2) <= 1e-12 +@test norm(y1-y2) <= 1e-12 opS = Scale(coeff, opA1) @test is_null(opS) == is_null(opA1) @@ -397,7 +397,7 @@ opS = Scale(coeff, opA1) @test is_full_row_rank(opS) == is_full_row_rank(opA1) @test is_full_column_rank(opS) == is_full_column_rank(opA1) -op = Scale(-4.0,DFT(10)) +op = Scale(-4.0,AbstractOperators.DFT(10)) @test is_AAc_diagonal(op) == true @test diag_AAc(op) == 16*10 @@ -409,14 +409,14 @@ op = Scale(-4.0,ZeroPad((10,), 20)) d = randn(10) op = Scale(3,DiagOp(d)) @test typeof(op) <: DiagOp -@test norm(diag(op) - 3.*d) < 1e-12 +@test norm(diag(op) - 3 .*d) < 1e-12 # Scale with imaginary coeff gives error m, n = 8, 4 coeff = im A1 = randn(m, n) opA1 = MatrixOp(A1) -@test_throws ErrorException opS = Scale(coeff, opA1) +@test_throws ErrorException Scale(coeff, opA1) ## testing displacement m, n = 8, 4 @@ -428,10 +428,10 @@ opS = Scale(coeff, opA1) x1 = randn(n) y1 = opS*x1 y2 = coeff*(A1*x1+d1) -@test vecnorm(y1-y2) <= 1e-12 +@test norm(y1-y2) <= 1e-12 y1 = remove_displacement(opS)*x1 y2 = coeff*(A1*x1) -@test vecnorm(y1-y2) <= 1e-12 +@test norm(y1-y2) <= 1e-12 ######################### #### test Sum ####### @@ -446,7 +446,7 @@ opS = Sum(opA1,opA2) x1 = randn(n) y1 = test_op(opS, x1, randn(m), verb) y2 = A1*x1+A2*x1 -@test vecnorm(y1-y2) <= 1e-12 +@test norm(y1-y2) <= 1e-12 #test Sum longer m,n = 5,7 @@ -460,11 +460,11 @@ opS = Sum(opA1,opA2,opA3) x1 = randn(n) y1 = test_op(opS, x1, randn(m), verb) y2 = A1*x1+A2*x1+A3*x1 -@test vecnorm(y1-y2) <= 1e-12 +@test norm(y1-y2) <= 1e-12 opA3 = MatrixOp(randn(m,m)) @test_throws Exception Sum(opA1,opA3) -opF = DFT(Float64,(m,)) +opF = AbstractOperators.DFT(Float64,(m,)) @test_throws Exception Sum(opF,opA3) @test is_null(opS) == false @@ -480,7 +480,7 @@ opF = DFT(Float64,(m,)) d = randn(10) op = Sum(Scale(-3.1,Eye(10)),DiagOp(d)) @test is_diagonal(op) == true -@test norm( diag(op) - (d-3.1) )<1e-12 +@test norm( diag(op) - (d .-3.1) )<1e-12 #test displacement of sum m,n = 5,7 @@ -495,25 +495,25 @@ opA2 = AffineAdd(MatrixOp(A2), d2) opA3 = AffineAdd(MatrixOp(A3), d3) opS = Sum(opA1,opA2,opA3) x1 = randn(n) -y2 = A1*x1+A2*x1+A3*x1+d1+d2+d3 -@test vecnorm(opS*x1-y2) <= 1e-12 -@test vecnorm(displacement(opS) - (d1+d2+d3)) <= 1e-12 +y2 = A1*x1+A2*x1+A3*x1+d1.+d2+d3 +@test norm(opS*x1-y2) <= 1e-12 +@test norm(displacement(opS) - (d1.+d2+d3)) <= 1e-12 y2 = A1*x1+A2*x1+A3*x1 -@test vecnorm(remove_displacement(opS)*x1-y2) <= 1e-12 +@test norm(remove_displacement(opS)*x1-y2) <= 1e-12 -########################### -###### test Transpose###### -########################### +################################### +###### test AdjointOperator ####### +################################### m,n = 5,7 A1 = randn(m,n) opA1 = MatrixOp(A1) opA1t = MatrixOp(A1') -opT = Transpose(opA1) +opT = AdjointOperator(opA1) x1 = randn(m) y1 = test_op(opT, x1, randn(n), verb) y2 = A1'*x1 -@test vecnorm(y1-y2) <= 1e-12 +@test norm(y1-y2) <= 1e-12 @test is_null(opT) == is_null(opA1t) @test is_eye(opT) == is_eye(opA1t) @@ -526,11 +526,11 @@ y2 = A1'*x1 @test is_full_column_rank(opT) == is_full_column_rank(opA1t) d = randn(3) -op = Transpose(DiagOp(d)) +op = AdjointOperator(DiagOp(d)) @test is_diagonal(op) == true @test diag(op) == d -op = Transpose(ZeroPad((10,),5)) +op = AdjointOperator(ZeroPad((10,),5)) @test is_AcA_diagonal(op) == false @test is_AAc_diagonal(op) == true @test diag_AAc(op) == 1 @@ -548,17 +548,17 @@ opV = VCAT(opA1, opA2) x1 = randn(n) y1 = test_op(opV, x1, (randn(m1), randn(m2)), verb) y2 = (A1*x1, A2*x1) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) m1, n = 4, 5 A1 = randn(m1, n)+im*randn(m1, n) opA1 = MatrixOp(A1) -opA2 = DFT(n)' +opA2 = AbstractOperators.DFT(n)' opV = VCAT(opA1, opA2) x1 = randn(n)+im*randn(n) y1 = test_op(opV, x1, (randn(m1)+im*randn(m1), randn(n)), verb) y2 = (A1*x1, opA2*x1) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) #test VCAT longer m1, m2, m3, n = 4, 7, 3, 5 @@ -572,22 +572,22 @@ opV = VCAT(opA1, opA2, opA3) x1 = randn(n) y1 = test_op(opV, x1, (randn(m1), randn(m2), randn(m3)), verb) y2 = (A1*x1, A2*x1, A3*x1) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) #test VCAT of VCAT opVV = VCAT(opV,opA3) y1 = test_op(opVV, x1, (randn(m1), randn(m2), randn(m3), randn(m3)), verb) y2 = (A1*x1, A2*x1, A3*x1, A3*x1) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) opVV = VCAT(opA1,opV,opA3) y1 = test_op(opVV, x1, (randn(m1), randn(m1), randn(m2), randn(m3), randn(m3)), verb) y2 = (A1*x1, A1*x1, A2*x1, A3*x1, A3*x1) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) opA3 = MatrixOp(randn(m1,m1)) @test_throws Exception VCAT(opA1,opA2,opA3) -opF = DFT(Complex{Float64},(n,)) +opF = AbstractOperators.DFT(Complex{Float64},(n,)) @test_throws Exception VCAT(opA1,opF,opA2) ###properties @@ -610,7 +610,7 @@ op = VCAT(opA1, opA2, opA3) @test is_full_row_rank(op) == false @test is_full_column_rank(op) == true -op = VCAT(DFT(Complex{Float64},10), Eye(Complex{Float64},10) ) +op = VCAT(AbstractOperators.DFT(Complex{Float64},10), Eye(Complex{Float64},10) ) @test is_AcA_diagonal(op) == true @test diag_AcA(op) == 11 @@ -626,10 +626,10 @@ opV = VCAT(AffineAdd(opA1,d1), AffineAdd(opA2,d2)) x1 = randn(n) y1 = opV*x1 y2 = (A1*x1+d1, A2*x1+d2) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) y1 = remove_displacement(opV)*x1 y2 = (A1*x1, A2*x1) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) ######################### #### test BroadCast ##### @@ -644,7 +644,7 @@ x1 = randn(n) y1 = test_op(opR, x1, randn(dim_out), verb) y2 = zeros(dim_out) y2 .= A1*x1 -@test vecnorm(y1-y2) <= 1e-12 +@test norm(y1-y2) <= 1e-12 m, n, l, k = 8, 4, 5, 7 dim_out = (m, n, l, k) @@ -654,7 +654,7 @@ x1 = randn(m,n) y1 = test_op(opR, x1, randn(dim_out), verb) y2 = zeros(dim_out) y2 .= x1 -@test vecnorm(y1-y2) <= 1e-12 +@test norm(y1-y2) <= 1e-12 @test_throws Exception BroadCast(opA1,(m,m)) @@ -666,7 +666,7 @@ x1 = randn(m) y1 = test_op(opR, x1, randn(dim_out), verb) y2 = zeros(dim_out) y2 .= x1 -@test vecnorm(y1-y2) <= 1e-12 +@test norm(y1-y2) <= 1e-12 #colum in - matrix out m, l = 4, 5 @@ -677,13 +677,13 @@ x1 = randn(1,l) y1 = test_op(opR, x1, randn(dim_out), verb) y2 = zeros(dim_out) y2 .= x1 -@test vecnorm(y1-y2) <= 1e-12 +@test norm(y1-y2) <= 1e-12 op = HCAT(Eye(m,l),opR) x1 = (randn(m,l),randn(1,l)) y1 = test_op(op, x1, randn(dim_out), verb) y2 = x1[1].+x1[2] -@test vecnorm(y1-y2) <= 1e-12 +@test norm(y1-y2) <= 1e-12 m, n, l = 2, 5, 8 dim_out = (m, n, l) @@ -693,7 +693,7 @@ x1 = randn(m) y1 = test_op(opR, x1, randn(dim_out), verb) y2 = zeros(dim_out) y2 .= x1 -@test vecnorm(y1-y2) <= 1e-12 +@test norm(y1-y2) <= 1e-12 m, n, l = 1, 5, 8 dim_out = (m, n, l) @@ -703,7 +703,7 @@ x1 = randn(m) y1 = test_op(opR, x1, randn(dim_out), verb) y2 = zeros(dim_out) y2 .= x1 -@test vecnorm(y1-y2) <= 1e-12 +@test norm(y1-y2) <= 1e-12 m, n, l = 1, 5, 8 dim_out = (m, n, l) @@ -713,7 +713,7 @@ x1 = randn(m) y1 = test_op(opR, x1, randn(dim_out), verb) y2 = zeros(dim_out) y2 .= 2.4*x1 -@test vecnorm(y1-y2) <= 1e-12 +@test norm(y1-y2) <= 1e-12 @test is_null(opR) == is_null(opA1) @test is_eye(opR) == false @@ -737,16 +737,16 @@ x1 = randn(n) y1 = opR*x1 y2 = zeros(dim_out) y2 .= A1*x1+d1 -@test vecnorm(y1-y2) <= 1e-12 +@test norm(y1-y2) <= 1e-12 x1 = randn(n) y1 = remove_displacement(opR)*x1 y2 = zeros(dim_out) y2 .= A1*x1 -@test vecnorm(y1-y2) <= 1e-12 +@test norm(y1-y2) <= 1e-12 -########################## -#### test AffineAdd ##### -########################## +########################### +##### test AffineAdd ##### +########################### n,m = 5,6 A = randn(n,m) @@ -757,9 +757,9 @@ T = AffineAdd(opA,d) println(T) x1 = randn(m) y1 = T*x1 -@test vecnorm(y1-(A*x1+d)) <1e-9 +@test norm(y1-(A*x1+d)) <1e-9 r = randn(n) -@test vecnorm(T'*r-(A'*r)) <1e-9 +@test norm(T'*r-(A'*r)) <1e-9 @test displacement(T) == d @test norm(remove_displacement(T)*x1-A*x1) <1e-9 @@ -770,9 +770,9 @@ T = AffineAdd(opA,d,false) println(T) x1 = randn(m) y1 = T*x1 -@test vecnorm(y1-(A*x1-d)) <1e-9 +@test norm(y1-(A*x1-d)) <1e-9 r = randn(n) -@test vecnorm(T'*r-(A'*r)) <1e-9 +@test norm(T'*r-(A'*r)) <1e-9 @test displacement(T) == -d @test norm(remove_displacement(T)*x1-A*x1) <1e-9 @@ -783,15 +783,15 @@ T = AffineAdd(opA,pi) println(T) x1 = randn(m) y1 = T*x1 -@test vecnorm(y1-(A*x1+pi)) <1e-9 +@test norm(y1-(A*x1.+pi)) <1e-9 r = randn(n) -@test vecnorm(T'*r-(A'*r)) < 1e-9 +@test norm(T'*r-(A'*r)) < 1e-9 @test displacement(T) .- pi < 1e-9 @test norm(remove_displacement(T)*x1-A*x1) <1e-9 @test_throws DimensionMismatch AffineAdd(MatrixOp(randn(2,5)),randn(5)) -@test_throws ErrorException AffineAdd(DFT(4),randn(4)) -AffineAdd(DFT(4),pi) +@test_throws ErrorException AffineAdd(AbstractOperators.DFT(4),randn(4)) +AffineAdd(AbstractOperators.DFT(4),pi) @test_throws ErrorException AffineAdd(Eye(4),im*pi) # with scalar and vector @@ -801,13 +801,13 @@ T = AffineAdd(AffineAdd(opA,pi),d,false) println(T) x1 = randn(m) y1 = T*x1 -@test vecnorm(y1-(A*x1+pi-d)) <1e-9 +@test norm(y1-(A*x1 .+ pi .- d)) <1e-9 r = randn(n) -@test vecnorm(T'*r-(A'*r)) < 1e-9 -@test vecnorm(displacement(T) .- (pi .-d )) < 1e-9 +@test norm(T'*r-(A'*r)) < 1e-9 +@test norm(displacement(T) .- (pi .-d )) < 1e-9 T2 = remove_displacement(T) -@test vecnorm(T2*x1-(A*x1)) <1e-9 +@test norm(T2*x1-(A*x1)) <1e-9 # permute AddAffine n,m = 5,6 @@ -817,13 +817,13 @@ opH = HCAT(Eye(n),MatrixOp(A)) x = (randn(n),randn(m)) opHT = AffineAdd(opH,d) -@test norm(opHT*x-(x[1]+A*x[2]+d)) < 1e-12 +@test norm(opHT*x-(x[1]+A*x[2].+d)) < 1e-12 p = [2;1] -@test norm(permute(opHT,p)*x[p]-(x[1]+A*x[2]+d)) < 1e-12 +@test norm(AbstractOperators.permute(opHT,p)*x[p]-(x[1]+A*x[2].+d)) < 1e-12 -############################## -######### test combin. ####### -############################## +############################# +######## test combin. ####### +############################# ## test Compose of HCAT m1, m2, m3, m4 = 4, 7, 3, 2 @@ -840,12 +840,12 @@ y1 = test_op(opC, (x1,x2), randn(m4), verb) y2 = A3*(A1*x1+A2*x2) -@test vecnorm(y1-y2) < 1e-9 +@test norm(y1-y2) < 1e-9 -opCp = permute(opC,[2,1]) +opCp = AbstractOperators.permute(opC,[2,1]) y1 = test_op(opCp, (x2,x1), randn(m4), verb) -@test vecnorm(y1-y2) < 1e-9 +@test norm(y1-y2) < 1e-9 ## test HCAT of Compose of HCAT m5 = 10 @@ -855,17 +855,17 @@ opHC = HCAT(opC,MatrixOp(A4)) x = (x1,x2,x3) y1 = test_op(opHC, x, randn(m4), verb) -@test vecnorm(y1-(y2+A4*x3)) < 1e-9 +@test norm(y1-(y2+A4*x3)) < 1e-9 p = randperm(ndoms(opHC,2)) -opHP = permute(opHC,p) +opHP = AbstractOperators.permute(opHC,p) xp = x[p] y1 = test_op(opHP, xp, randn(m4), verb) pp = randperm(ndoms(opHC,2)) -opHPP = permute(opHC,pp) +opHPP = AbstractOperators.permute(opHC,pp) xpp = x[pp] y1 = test_op(opHPP, xpp, randn(m4), verb) @@ -889,26 +889,26 @@ opV = VCAT(opH1,opH2) x1, x2 = randn(m1), randn(m2) y1 = test_op(opV, (x1,x2), (randn(n1),randn(n2)), verb) y2 = (A1*x1+A2*x2,A3*x1+A4*x2) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) # test VCAT of HCAT's with complex num m1, m2, n1 = 4, 7, 5 A1 = randn(n1, m1)+im*randn(n1, m1) opA1 = MatrixOp(A1) -opA2 = DFT(n1) +opA2 = AbstractOperators.DFT(n1) opH1 = HCAT(opA1,opA2) m1, m2, n2 = 4, 7, 5 A3 = randn(n2, m1)+im*randn(n2,m1) opA3 = MatrixOp(A3) -opA4 = DFT(n2) +opA4 = AbstractOperators.DFT(n2) opH2 = HCAT(opA3,opA4) opV = VCAT(opH1,opH2) x1, x2 = randn(m1)+im*randn(m1), randn(n2) y1 = test_op(opV, (x1,x2), (randn(n1)+im*randn(n1),randn(n2)+im*randn(n2)), verb) y2 = (A1*x1+fft(x2),A3*x1+fft(x2)) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) # test HCAT of VCAT's @@ -923,9 +923,9 @@ x2 = randn(n2) y1 = test_op(opV, (x1, x2), (randn(m1), randn(m2)), verb) y2 = (A*x1 + B*x2, C*x1 + D*x2) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) -## test Sum of HCAT's +# test Sum of HCAT's m, n1, n2, n3 = 4, 7, 5, 3 A1 = randn(m, n1) @@ -949,10 +949,10 @@ x3 = randn(n3) y1 = test_op(opS, (x1, x2, x3), randn(m), verb) y2 = A1*x1 + B1*x1 + A2*x2 + B2*x2 + A3*x3 + B3*x3 -@test vecnorm(y1-y2) <= 1e-12 +@test norm(y1-y2) <= 1e-12 p = [3;2;1] -opSp = permute(opS,p) +opSp = AbstractOperators.permute(opS,p) y1 = test_op(opSp, (x1, x2, x3)[p], randn(m), verb) # test Sum of VCAT's @@ -978,7 +978,7 @@ x = randn(n) y1 = test_op(opS, x, (randn(m1), randn(m2)), verb) y2 = (A1*x + B1*x +C1*x, A2*x + B2*x + C2*x) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) # test Scale of DCAT @@ -996,7 +996,7 @@ x2 = randn(n2) y = test_op(opS, (x1, x2), (randn(m1), randn(m2)), verb) z = (coeff*A1*x1, coeff*A2*x2) -@test all(vecnorm.(y .- z) .<= 1e-12) +@test all(norm.(y .- z) .<= 1e-12) # test Scale of VCAT @@ -1012,7 +1012,7 @@ x = randn(n) y = test_op(opS, x, (randn(m1), randn(m2)), verb) z = (coeff*A1*x, coeff*A2*x) -@test all(vecnorm.(y .- z) .<= 1e-12) +@test all(norm.(y .- z) .<= 1e-12) # test Scale of HCAT @@ -1029,7 +1029,7 @@ x2 = randn(n2) y = test_op(opS, (x1, x2), randn(m), verb) z = coeff*(A1*x1 + A2*x2) -@test all(vecnorm.(y .- z) .<= 1e-12) +@test all(norm.(y .- z) .<= 1e-12) # test DCAT of HCATs @@ -1075,7 +1075,7 @@ opSS = Scale(coeff,opS) x1 = randn(n) y1 = test_op(opSS, x1, randn(m), verb) y2 = coeff*(A1*x1+A2*x1) -@test vecnorm(y1-y2) <= 1e-12 +@test norm(y1-y2) <= 1e-12 # test Scale of Compose @@ -1091,5 +1091,5 @@ opS = Scale(coeff,opC) x = randn(m1) y1 = test_op(opS, x, randn(m3), verb) y2 = coeff*(A2*A1*x) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) diff --git a/test/test_nonlinear_operators.jl b/test/test_nonlinear_operators.jl index 8323bea..ea64163 100644 --- a/test/test_nonlinear_operators.jl +++ b/test/test_nonlinear_operators.jl @@ -30,7 +30,7 @@ op = SoftMax(Float64,(n,m,l)) y, grad = test_NLop(op,x,r,verb) -# SoftMax +# SoftPlus n = 10 x = randn(n) r = randn(n) @@ -52,7 +52,7 @@ op = Exp(Float64,(n,m,l)) y, grad = test_NLop(op,x,r,verb) -# Sin +## Sin n,m,l = 4,5,6 x = randn(n,m,l) r = randn(n,m,l) @@ -61,7 +61,7 @@ op = Sin(Float64,(n,m,l)) y, grad = test_NLop(op,x,r,verb) -# Sin + Cos n,m,l = 4,5,6 x = randn(n,m,l) r = randn(n,m,l) diff --git a/test/test_nonlinear_operators_calculus.jl b/test/test_nonlinear_operators_calculus.jl index 24b5520..a1164cf 100644 --- a/test/test_nonlinear_operators_calculus.jl +++ b/test/test_nonlinear_operators_calculus.jl @@ -9,7 +9,7 @@ op = 30*A y, grad = test_NLop(op,x,r,verb) Y = 30*(A*x) -@test vecnorm(Y-y) <1e-8 +@test norm(Y-y) <1e-8 m = 3 x = randn(m) @@ -20,7 +20,7 @@ op = -A y, grad = test_NLop(op,x,r,verb) Y = -A*x -@test vecnorm(Y-y) <1e-8 +@test norm(Y-y) <1e-8 #testing DCAT n,m = 4,3 @@ -33,8 +33,8 @@ op = DCAT(MatrixOp(A),B) y, grad = test_NLop(op,x,r,verb) Y = (A*x[1],B*x[2]) -@test vecnorm(Y[1]-y[1]) <1e-8 -@test vecnorm(Y[2]-y[2]) <1e-8 +@test norm(Y[1]-y[1]) <1e-8 +@test norm(Y[2]-y[2]) <1e-8 #testing HCAT n,m = 4,3 @@ -47,7 +47,7 @@ op = HCAT(MatrixOp(A),B) y, grad = test_NLop(op,x,r,verb) Y = A*x[1]+B*x[2] -@test vecnorm(Y-y) <1e-8 +@test norm(Y-y) <1e-8 m,n = 3,5 x = (randn(m),randn(n)) @@ -60,10 +60,10 @@ op = HCAT(A,B) y, grad = test_NLop(op,x,r,verb) Y = A*x[1]+M*x[2] -@test vecnorm(Y-y) <1e-8 +@test norm(Y-y) <1e-8 p = [2,1] -opP = permute(op,p) +opP = AbstractOperators.permute(op,p) J = Jacobian(opP,x[p])' println(size(J,1)) y, grad = test_NLop(opP,x[p],r,verb) @@ -79,8 +79,8 @@ op = VCAT(MatrixOp(A),B) y, grad = test_NLop(op,x,r,verb) Y = (A*x,B*x) -@test vecnorm(Y[1]-y[1]) <1e-8 -@test vecnorm(Y[2]-y[2]) <1e-8 +@test norm(Y[1]-y[1]) <1e-8 +@test norm(Y[2]-y[2]) <1e-8 #testing HCAT of VCAT n,m1,m2,m3 = 4,3,2,7 @@ -103,8 +103,8 @@ op = HCAT(op1,op2,op3) y, grad = test_NLop(op,x,r,verb) Y = (A1*x1+A2*x2+A3*x3,B1*x1+B2*x2+B3*x3) -@test vecnorm(Y[1]-y[1]) <1e-8 -@test vecnorm(Y[2]-y[2]) <1e-8 +@test norm(Y[1]-y[1]) <1e-8 +@test norm(Y[2]-y[2]) <1e-8 #testing VCAT of HCAT m1,m2,m3,n1,n2 = 3,4,5,6,7 @@ -127,8 +127,8 @@ op = VCAT(op1,op2) y, grad = test_NLop(op,x,r,verb) Y = (A1*x1+B1*x2+C1*x3,A2*x1+B2*x2+C2*x3) -@test vecnorm(Y[1]-y[1]) <1e-8 -@test vecnorm(Y[2]-y[2]) <1e-8 +@test norm(Y[1]-y[1]) <1e-8 +@test norm(Y[2]-y[2]) <1e-8 #testing Compose @@ -146,7 +146,7 @@ op = Compose(opA,Compose(opB,opC)) y, grad = test_NLop(op,x,r,verb) Y = A*(opB*(opC*x)) -@test vecnorm(Y-y) <1e-8 +@test norm(Y-y) <1e-8 ## NN m,n,l = 4,7,5 @@ -159,7 +159,6 @@ A1 = HCAT(LMatrixOp(b,n) ,Eye(n)) op = Compose(opS1,A1) y, grad = test_NLop(op,x,r,verb) - ###testing Reshape n = 4 x = randn(n) @@ -170,7 +169,7 @@ op = Reshape(opS,2,2) y, grad = test_NLop(op,x,r,verb) Y = reshape(opS*x,2,2) -@test vecnorm(Y-y) <1e-8 +@test norm(Y-y) <1e-8 ###testing BroadCast n,l = 4,7 @@ -182,7 +181,7 @@ op = BroadCast(opS,(n,l)) y, grad = test_NLop(op,x,r,verb) Y = (opS*x).*ones(n,l) -@test vecnorm(Y-y) <1e-8 +@test norm(Y-y) <1e-8 n,l = 1,7 x = randn(n) @@ -193,7 +192,7 @@ op = BroadCast(opS,(n,l)) y, grad = test_NLop(op,x,r,verb) Y = (opS*x).*ones(n,l) -@test vecnorm(Y-y) <1e-8 +@test norm(Y-y) <1e-8 ##testing Sum m = 5 @@ -207,9 +206,9 @@ op = Sum(opA,opB) y, grad = test_NLop(op,x,r,verb) Y = A*x+opB*x -@test vecnorm(Y-y) <1e-8 +@test norm(Y-y) <1e-8 -# testing NonLinearCompose +## testing NonLinearCompose ##with vectors inner product n,m = 3,4 @@ -223,7 +222,7 @@ y, grad = test_NLop(P,x,r,verb) Y = x[1]*(A*x[2]) @test norm(Y - y) <= 1e-12 -#with vectors outer product +##with vectors outer product n, m = 3, 5 x = (randn(n),randn(1,m)) r = randn(n,m) @@ -271,8 +270,8 @@ Y = A*x[1]*B*x[2] #further test on gradient with analytical formulas grad2 = (A'*r)*(B*x[2])', B'*(A*x[1])'*r -@test vecnorm(grad[1]-grad2[1]) <1e-7 -@test vecnorm(grad[2]-grad2[2]) <1e-7 +@test norm(grad[1]-grad2[1]) <1e-7 +@test norm(grad[2]-grad2[2]) <1e-7 #with complex matrices l,m1,m2,n1,n2 = 2,3,4,5,6 @@ -289,8 +288,8 @@ Y = A*x[1]*B*x[2] ##test on gradient with analytical formulas grad2 = (A'*r)*(B*x[2])', B'*(A*x[1])'*r -@test vecnorm(grad[1]-grad2[1]) <1e-7 -@test vecnorm(grad[2]-grad2[2]) <1e-7 +@test norm(grad[1]-grad2[1]) <1e-7 +@test norm(grad[2]-grad2[2]) <1e-7 #nested NonLinearOp l1,l2,m1,m2,n1,n2 = 2,3,4,5,6,7 @@ -309,12 +308,12 @@ Y = A*x[1]*B*x[2]*C*x[3] #further test on gradient with analytical formulas grad2 = A'*(r*(B*x[2]*C*x[3])'), B'*((r'*A*x[1])'*(C*x[3])'), C'*(B*x[2])'*(A*x[1])'*r -@test vecnorm(grad[1]-grad2[1]) <1e-7 -@test vecnorm(grad[2]-grad2[2]) <1e-7 -@test vecnorm(grad[3]-grad2[3]) <1e-7 +@test norm(grad[1]-grad2[1]) <1e-7 +@test norm(grad[2]-grad2[2]) <1e-7 +@test norm(grad[3]-grad2[3]) <1e-7 p = randperm(length(x)) -Pp = permute(P,p) +Pp = AbstractOperators.permute(P,p) y, grad = test_NLop(Pp,x[p],r,verb) ## DNN @@ -340,7 +339,7 @@ Y = opS3*(x[1]*(opS2*(x[2]*(opS1*(x[3]*b+x[4]))))) @test norm(Y - y) <= 1e-12 p = randperm(length(x)) -L3p = permute(L3,p) +L3p = AbstractOperators.permute(L3,p) y, grad = test_NLop(L3p,x[p],r,verb) # Hadamard @@ -356,7 +355,7 @@ y, grad = test_NLop(H,x,r,verb) @test norm(y-(op1.A*x[1]).*(op2.A*x[2])) < 1e-9 p = [2;1] -Hp = permute(H,p) +Hp = AbstractOperators.permute(H,p) y, grad = test_NLop(Hp,x[p],r,verb) @test norm(y-(op1.A*x[1]).*(op2.A*x[2])) < 1e-9 @@ -373,20 +372,20 @@ y, grad = test_NLop(H,x,r,verb) @test norm(y-(op1.A*x[1]).*(op2.A*x[2]).*(op3.A*x[3])) < 1e-9 p = [2;1;3] -Hp = permute(H,p) +Hp = AbstractOperators.permute(H,p) y, grad = test_NLop(Hp,x[p],r,verb) @test norm(y-(op1.A*x[1]).*(op2.A*x[2]).*(op3.A*x[3])) < 1e-9 l,m,n = 10,3,7 A = randn(n,m)+im*randn(n,m) op1 = MatrixOp(A) -op3 = DFT(n) +op3 = AbstractOperators.DFT(n) H = Hadamard(op1,op3) x = randn(m)+im*randn(m),randn(n) y = H*x -@test vecnorm(y - (A*x[1]).*(fft(x[2])) ) <1e-9 +@test norm(y - (A*x[1]).*(fft(x[2])) ) <1e-9 grad = Jacobian(H,x)'*y # TODO add test on gradient @@ -402,7 +401,7 @@ x = randn.(size(H,2)) y, grad = test_NLop(H,x,r,verb) p = [2;1;3] -Hp = permute(H,p) +Hp = AbstractOperators.permute(H,p) y, grad = test_NLop(Hp,x[p],r,verb) ## Hadamard of Hadamard with NonLinear operators @@ -425,7 +424,7 @@ T = AffineAdd(Exp(n),d,false) r = randn(n) x = randn(size(T,2)) y, grad = test_NLop(T,x,r,verb) -@test vecnorm( y - (exp.(x)-d) ) < 1e-8 +@test norm( y - (exp.(x)-d) ) < 1e-8 ## AffineAdd and Compose NonLinearOperator n = 10 @@ -436,7 +435,7 @@ T = Compose(AffineAdd(Sin(n),d2),AffineAdd(Eye(n),d1)) r = randn(n) x = randn(size(T,2)) y, grad = test_NLop(T,x,r,verb) -@test vecnorm( y - (sin.(x+d1)+d2) ) < 1e-8 +@test norm( y - (sin.(x+d1)+d2) ) < 1e-8 n = 10 d1 = randn(n) @@ -447,7 +446,7 @@ T = Compose( AffineAdd(Sin(n),d3), Compose( AffineAdd(Exp(n),d2,false),AffineAdd r = randn(n) x = randn(size(T,2)) y, grad = test_NLop(T,x,r,verb) -@test vecnorm( y -( sin.(exp.(x+d1)-d2)+d3 )) < 1e-8 +@test norm( y -( sin.(exp.(x+d1)-d2).+d3 )) < 1e-8 ## AffineAdd and NonLinearCompose l,m1,m2,n1,n2 = 2,3,4,5,6 diff --git a/test/test_syntax.jl b/test/test_syntax.jl index eecfd69..eae744b 100644 --- a/test/test_syntax.jl +++ b/test/test_syntax.jl @@ -1,24 +1,5 @@ @printf("\nTesting linear operators syntax\n") -##### blkdiag #### -n1,m1 = 2,3 -n2,m2 = 4,5 -n3,m3 = 6,7 -A = randn(n1,m1) -B = randn(n2,m2) -C = randn(n3,m3) -opA = MatrixOp(A) -opB = MatrixOp(B) -opC = MatrixOp(C) -x1 = randn(m1) -x2 = randn(m2) -x3 = randn(m3) - -opB = blkdiag(opA,opB,opC) -y1 = opB*(x1,x2,x3) -y2 = (A*x1,B*x2,C*x3) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) - ###### ' ###### n,m = 5,3 A = randn(n,m) @@ -124,20 +105,6 @@ y1 = opc*x1 y2 = x1 @test norm(y1-y2) < 1e-9 -###### .* ###### -d = randn(n) - -opd = d.*opA -y1 = opd*x1 -y2 = d.*(A*x1) -@test norm(y1-y2) < 1e-9 - -ops = alpha*opA -opd = d.*ops -y1 = opd*x1 -y2 = d.*(alpha*A*x1) -@test norm(y1-y2) < 1e-9 - ###### getindex ###### ops = opA[1:n-1] @@ -150,20 +117,20 @@ x3 = randn(n,m,l) ops = opF[1:n-1,:,2:l] y1 = ops*x3 y2 = dct(x3)[1:n-1,:,2:l] -@test vecnorm(y1-y2) < 1e-9 +@test norm(y1-y2) < 1e-9 opF = DCT(n,m,l) x3 = randn(n,m,l) -ops = opF[1:n-1,2:m] +ops = opF[1:n-1,2:m,1] y1 = ops*x3 -y2 = dct(x3)[1:n-1,2:m] -@test vecnorm(y1-y2) < 1e-9 +y2 = dct(x3)[1:n-1,2:m,1] +@test norm(y1-y2) < 1e-9 opV = Variation(n,m,l) ops = opV[1:4] y1 = ops*x3 y2 = (opV*x3)[1:4] -@test vecnorm(y1-y2) < 1e-9 +@test norm(y1-y2) < 1e-9 ops = (opB*opA)[1:l-1] y1 = ops*x1 @@ -191,11 +158,11 @@ opH = HCAT(opA,opB,opC) opH2 = opH[1:2] y1 = opH2*(x1,x2) y2 = A*x1+B*x2 -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) opH3 = opH[3] y1 = opH3*x3 y2 = C*x3 -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) opHperm = opH[[3,1,2]] @test norm(opH*(x1,x2,x3) - opHperm*(x3,x1,x2)) <1e-12 @@ -243,11 +210,11 @@ opV = VCAT(opA,opB,opC) opV2 = opV[1:2] y1 = opV2*x1 y2 = (A*x1,B*x1) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) opV3 = opV[3] y1 = opV3*x3 y2 = C*x3 -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) ###### hcat ###### @@ -261,12 +228,12 @@ x1 = randn(m1) x2 = randn(m2) y1 = opH*(x1,x2) y2 = [A B]*[x1;x2] -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) opHH = [opH opB] y1 = opHH*(x1, x2, x2) y2 = [A B B]*[x1;x2;x2] -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) ###### vcat ###### @@ -279,12 +246,12 @@ opH = [opA; opB] x1 = randn(n) y1 = opH*x1 y2 = (A*x1,B*x1) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) opVV = [opA; opH] y1 = opVV*x1 y2 = (A*x1, A*x1, B*x1) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) ###### reshape ###### n,m = 10,5 @@ -295,7 +262,7 @@ opR = reshape(opA,2,5) opR = reshape(opA,(2,5)) y1 = opR*x1 y2 = reshape(A*x1,2,5) -@test all(vecnorm.(y1 .- y2) .<= 1e-12) +@test all(norm.(y1 .- y2) .<= 1e-12) # testing ndims & ndoms L = Variation((3,4,5)) @@ -310,7 +277,7 @@ H = hcat(L,L) @test ndoms(H) == (1,2) @test ndoms(H,1) == 1 @test ndoms(H,2) == 2 -D = blkdiag(L,L) +D = DCAT(L,L) @test ndims(D) == ((2,2),(3,3)) @test ndoms(D) == (2,2) @@ -329,7 +296,7 @@ LL = convert(AbstractOperator, Float64, (10,), L) LL = convert(LinearOperator, Float64, (10,), L) @test LL == L -@test_throws MethodError LL = convert(NonLinearOperator, Float64, (10,), L) +@test_throws MethodError convert(NonLinearOperator, Float64, (10,), L) ### displacement ### L = Eye(10) diff --git a/test/utils.jl b/test/utils.jl index d1793f7..21042ae 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -6,16 +6,17 @@ function test_op(A::AbstractOperator, x, y, verb::Bool = false) Ax = A*x Ax2 = AbstractOperators.blocksimilar(Ax) verb && println("forward preallocated") - A_mul_B!(Ax2, A, x) #verify in-place linear operator works - verb && @time A_mul_B!(Ax2, A, x) + mul!(Ax2, A, x) #verify in-place linear operator works + verb && @time mul!(Ax2, A, x) @test AbstractOperators.blockvecnorm(Ax .- Ax2) <= 1e-8 Acy = A'*y Acy2 = AbstractOperators.blocksimilar(Acy) verb && println("adjoint preallocated") - Ac_mul_B!(Acy2, A, y) #verify in-place linear operator works - verb && @time Ac_mul_B!(Acy2, A, y) + At = AdjointOperator(A) + mul!(Acy2, At, y) #verify in-place linear operator works + verb && @time mul!(Acy2, At, y) @test AbstractOperators.blockvecnorm(Acy .- Acy2) <= 1e-8 @@ -34,8 +35,8 @@ function test_NLop(A::AbstractOperator, x, y, verb::Bool = false) Ax = A*x Ax2 = AbstractOperators.blocksimilar(Ax) verb && println("forward preallocated") - A_mul_B!(Ax2, A, x) #verify in-place linear operator works - verb && @time A_mul_B!(Ax2, A, x) + mul!(Ax2, A, x) #verify in-place linear operator works + verb && @time mul!(Ax2, A, x) @test_throws ErrorException A' @@ -45,12 +46,12 @@ function test_NLop(A::AbstractOperator, x, y, verb::Bool = false) verb && println(J) grad = J'*y - A_mul_B!(Ax2, A, x) #redo forward - verb && println("jacobian Ac_mul_B! preallocated") + mul!(Ax2, A, x) #redo forward + verb && println("adjoint jacobian mul! preallocated") grad2 = AbstractOperators.blocksimilar(grad) - Ac_mul_B!(grad2, J, y) #verify in-place linear operator works - verb && A_mul_B!(Ax2, A, x) #redo forward - verb && @time Ac_mul_B!(grad2, J, y) + mul!(grad2, J', y) #verify in-place linear operator works + verb && mul!(Ax2, A, x) #redo forward + verb && @time mul!(grad2, J', y) @test AbstractOperators.blockvecnorm(grad .- grad2) < 1e-8 @@ -66,10 +67,10 @@ end ############# Finite Diff for Jacobian tests -function gradient_fd{A<:AbstractOperator}(op::A, - y0::AbstractArray, - x0::AbstractArray, - r::AbstractArray) +function gradient_fd(op::A, + y0::AbstractArray, + x0::AbstractArray, + r::AbstractArray) where {A<:AbstractOperator} y = copy(y0) J = zeros(*(size(op,1)...),*(size(op,2)...)) @@ -77,20 +78,20 @@ function gradient_fd{A<:AbstractOperator}(op::A, for i = 1:size(J,2) x = copy(x0) x[i] = x[i]+h - A_mul_B!(y,op,x) + mul!(y,op,x) J[:,i] .= ((y.-y0)./h)[:] end return reshape(J'*r[:],size(op,2)) end -function gradient_fd{N, A<:AbstractOperator}(op::A, - y0::AbstractArray, - x0::NTuple{N,AbstractArray}, - r::AbstractArray) +function gradient_fd(op::A, + y0::AbstractArray, + x0::NTuple{N,AbstractArray}, + r::AbstractArray) where {N, A<:AbstractOperator} y = copy(y0) - grad = zeros.(x0) + grad = AbstractOperators.blockzeros(x0) J = [ zeros(*(size(op,1)...),*(sz2...)) for sz2 in size(op,2)] h = sqrt(eps()) @@ -98,7 +99,7 @@ function gradient_fd{N, A<:AbstractOperator}(op::A, for i = 1:size(J[ii],2) x = deepcopy(x0) x[ii][i] = x[ii][i]+h - A_mul_B!(y,op,x) + mul!(y,op,x) J[ii][:,i] .= ((y.-y0)./h)[:] end grad[ii] .= reshape(J[ii]'*r[:],size(op,2)[ii]) @@ -106,20 +107,20 @@ function gradient_fd{N, A<:AbstractOperator}(op::A, return grad end -function gradient_fd{N, A<:AbstractOperator}(op::A, - y0::NTuple{N,AbstractArray}, - x0::AbstractArray, - r::NTuple{N,AbstractArray}) +function gradient_fd(op::A, + y0::NTuple{N,AbstractArray}, + x0::AbstractArray, + r::NTuple{N,AbstractArray}) where {N, A<:AbstractOperator} - y = zeros.(y0) - grad = zeros(x0) + y = AbstractOperators.blockzeros(y0) + grad = AbstractOperators.blockzeros(x0) J = [ zeros(*(sz1...),*(size(op,2)...)) for sz1 in size(op,1)] h = sqrt(eps()) for i in eachindex(x0) x = deepcopy(x0) x[i] = x[i]+h - A_mul_B!(y,op,x) + mul!(y,op,x) for ii in eachindex(J) J[ii][:,i] .= ((y[ii].-y0[ii])./h)[:] end @@ -130,12 +131,12 @@ function gradient_fd{N, A<:AbstractOperator}(op::A, return grad end -function gradient_fd{N,M, A<:AbstractOperator}(op::A, - y0::NTuple{N,AbstractArray}, - x0::NTuple{M,AbstractArray}, - r::NTuple{N,AbstractArray}) - grad = zeros.(x0) - y = zeros.(y0) +function gradient_fd(op::A, + y0::NTuple{N,AbstractArray}, + x0::NTuple{M,AbstractArray}, + r::NTuple{N,AbstractArray}) where {N,M, A<:AbstractOperator} + grad = AbstractOperators.blockzeros(x0) + y = AbstractOperators.blockzeros(y0) J = [ zeros(*(size(op,1)[i]...),*(size(op,2)[ii]...)) for ii = 1:M, i = 1:N ] @@ -144,7 +145,7 @@ function gradient_fd{N,M, A<:AbstractOperator}(op::A, for iii in eachindex(x[i]) x = deepcopy(x0) x[i][iii] = x[i][iii]+h - A_mul_B!(y,op,x) + mul!(y,op,x) for ii = 1:N J[i,ii][:,iii] .= ((y[ii].-y0[ii])./h)[:]