fixed MatrixOp, Eye constructor from BlockVector

Author: nantonel

src/AbstractOperators.jl

@@ -42,7 +42,6 @@ include("linearoperators/Filt.jl")
 include("linearoperators/MIMOFilt.jl")
 include("linearoperators/Xcorr.jl")
 include("linearoperators/LBFGS.jl")
-# include("linearoperators/BlkDiagLBFGS.jl")
 
 # Calculus rules
 

src/calculus/DCAT.jl

@@ -181,10 +181,11 @@ function codomainType(H::DCAT)
 	return (codomain...)
 end
 
+is_eye(L::DCAT) = all(is_eye.(L.A))
 is_linear(L::DCAT) = all(is_linear.(L.A))
 is_diagonal(L::DCAT) = all(is_diagonal.(L.A))
 is_AcA_diagonal(L::DCAT) = all(is_AcA_diagonal.(L.A))
-is_Ac_diagonal(L::DCAT) = all(is_Ac_diagonal.(L.A))
+is_AAc_diagonal(L::DCAT) = all(is_AAc_diagonal.(L.A))
 is_orthogonal(L::DCAT) = all(is_orthogonal.(L.A))
 is_invertible(L::DCAT) = all(is_invertible.(L.A))
 is_full_row_rank(L::DCAT) = all(is_full_row_rank.(L.A))
@@ -208,3 +209,10 @@ function permute(H::DCAT{N,L,P1,P2}, p::AbstractVector{Int}) where {N,L,P1,P2}
 
 	DCAT(H.A,new_part,H.idxC)
 end
+
+# special cases
+# Eye constructor
+Eye(x::A) where {N, A <: NTuple{N,AbstractArray}} = DCAT(Eye.(x)...)
+diag(L::DCAT{N,NTuple{N,E}}) where {N, E <: Eye} = 1.
+diag_AAc(L::DCAT{N,NTuple{N,E}}) where {N, E <: Eye} = 1.
+diag_AcA(L::DCAT{N,NTuple{N,E}}) where {N, E <: Eye} = 1.

src/linearoperators/Eye.jl

@@ -33,6 +33,7 @@ Eye(DomainType::Type, DomainDim::NTuple{N,Int}) where {N} = Eye{DomainType,N}(Do
 Eye(t::Type, dims::Vararg{Integer}) = Eye(t,dims)
 Eye(dims::NTuple{N, Integer}) where {N} = Eye(Float64,dims)
 Eye(dims::Vararg{Integer}) = Eye(Float64,dims)
+Eye(x::A) where {A <: AbstractArray} = Eye(eltype(x), size(x))
 
 # Mappings
 

src/linearoperators/LBFGS.jl

@@ -1,4 +1,4 @@
-export LBFGS, update!
+export LBFGS, update!, reset!
 
 """
 `LBFGS(domainType::Type,dim_in::Tuple, M::Integer)`
@@ -18,6 +18,9 @@ julia> update!(L,x,x_prev,grad,grad_prev); # update memory
 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
@@ -81,6 +84,17 @@ function update!(L::LBFGS{R, T, M, I}, x::T, x_prev::T, gradx::T, gradx_prev::T)
 	return L
 end
 
+"""
+`reset!(L::LBFGS)`
+
+Cancels the memory of `L`.
+"""
+
+function reset!(L::LBFGS)
+	L.currmem, L.curridx = 0, 0
+	L.H = 1.0
+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)

src/linearoperators/MatrixOp.jl

@@ -28,7 +28,7 @@ julia> MatrixOp(randn(20,10),4)
 
 """
 
-struct MatrixOp{T, M <: AbstractMatrix{T}} <: LinearOperator
+struct MatrixOp{D, T, M <: AbstractMatrix{T}} <: LinearOperator
 	A::M
 	n_col_in::Integer
 end
@@ -37,35 +37,43 @@ end
 
 ##TODO decide what to do when domainType is given, with conversion one loses pointer to data...
 ###standard constructor Operator{N}(DomainType::Type, DomainDim::NTuple{N,Int})
-function MatrixOp(DomainType::Type, DomainDim::NTuple{N,Int}, A::M) where {N, M <: AbstractMatrix} 
+function MatrixOp(DomainType::Type, DomainDim::NTuple{N,Int}, A::M) where {N, T, M <: AbstractMatrix{T}} 
 	N > 2 && error("cannot multiply a Matrix by a n-dimensional Variable with n > 2") 
 	size(A,2) != DomainDim[1] && error("wrong input dimensions")
 	if N == 1
-		MatrixOp{DomainType, M}(A, 1)
+		MatrixOp{DomainType, T, M}(A, 1)
 	else
-		MatrixOp{DomainType, M}(A, DomainDim[2])
+		MatrixOp{DomainType, T, M}(A, DomainDim[2])
 	end
 end
 ###
 
-MatrixOp(A::M)                      where {M <: AbstractMatrix} = MatrixOp{eltype(A), M}(A, 1)
-MatrixOp(T::Type, A::M)             where {M <: AbstractMatrix} = MatrixOp{T, M}(A, 1)
-MatrixOp(A::M, n::Integer)          where {M <: AbstractMatrix} = MatrixOp{eltype(A), M}(A, n)
-MatrixOp(T::Type, A::M, n::Integer) where {M <: AbstractMatrix} = MatrixOp{T, M}(A, n)
+MatrixOp(A::M)                      where {M <: AbstractMatrix} = MatrixOp(eltype(A), (size(A,2),), A)
+MatrixOp(D::Type, A::M)             where {M <: AbstractMatrix} = MatrixOp(D, (size(A,2),), A)
+MatrixOp(A::M, n::Integer)          where {M <: AbstractMatrix} = MatrixOp(eltype(A), (size(A,2), n), A)
+MatrixOp(D::Type, A::M, n::Integer) where {M <: AbstractMatrix} = MatrixOp(D, (size(A,2), n), A)
 
 import Base: convert
-convert(::Type{LinearOperator}, L::M) where {T,M<:AbstractMatrix{T}} = MatrixOp{T,M}(L,1)
-convert(::Type{LinearOperator}, L::M, n::Integer) where {T,M<:AbstractMatrix{T}} = MatrixOp{T,M}(L, n)
+convert(::Type{LinearOperator}, L::M) where {T, M<:AbstractMatrix{T}} = MatrixOp{T, T, M}(L,1)
+convert(::Type{LinearOperator}, L::M, n::Integer) where {T, M<:AbstractMatrix{T}} = MatrixOp{T, T, M}(L, n)
+convert(::Type{LinearOperator}, dom::Type, dim_in::Tuple, L::AbstractMatrix) = MatrixOp(dom, dim_in, L)
 
 # Mappings
 
- A_mul_B!(y::AbstractArray, L::MatrixOp{M, T}, b::AbstractArray) where {M, T} = A_mul_B!(y, L.A, b)
-Ac_mul_B!(y::AbstractArray, L::MatrixOp{M, T}, b::AbstractArray) where {M, T} = Ac_mul_B!(y, L.A, b)
+ 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)
+
+# Special Case, real b, complex matrix
+function Ac_mul_B!(y::AbstractArray, L::MatrixOp{D, T}, b::AbstractArray) where {D <: Real , T <: Complex} 
+	yc = zeros(T,size(y))
+	Ac_mul_B!(yc, L.A, b)
+	y .= real.(yc)
+end
 
 # Properties
 
-  domainType(L::MatrixOp{T, M}) where {T, M} = T
-codomainType(L::MatrixOp{T, M}) where {T, M} = T
+  domainType(L::MatrixOp{D, T}) where {D, T} = D
+codomainType(L::MatrixOp{D, T}) where {D, T} = D <: Real && T <: Complex ? T : D
 
 function size(L::MatrixOp)
 	if L.n_col_in == 1

src/properties.jl

@@ -48,7 +48,6 @@ julia> codomainType(vcat(Eye(Complex{Float64},(10,)),DFT(Complex{Float64},10)))
 """
 codomainType
 
-
 """
 `size(A::AbstractOperator, [dom,])`
 
@@ -254,6 +253,15 @@ false
 """
 is_full_column_rank(L::AbstractOperator) = false
 
+
+import Base: convert
+function convert(::Type{T}, dom::Type, dim_in::Tuple, L::T) where {T <: AbstractOperator}
+	domainType(L) != dom && error("cannot convert operator with domain $(domainType(L)) to operator with domain $dom ")
+	size(L,1) != dim_in && error("cannot convert operator with size $(size(L,1)) to operator with domain $dim_in ")
+	return L
+end
+
+
 #printing
 function Base.show(io::IO, L::AbstractOperator)
 	print(io, fun_name(L)" "*fun_space(L) )

src/utilities/block.jl

@@ -14,6 +14,7 @@ export RealOrComplex,
        blockset!,
        blockvecdot,
        blockzeros,
+       blockones,
        blockaxpy!,
        blockiszero
 
@@ -53,8 +54,8 @@ blockcopy!(y::AbstractArray, x::AbstractArray) = copy!(y, x)
 blockset!(y::Tuple, x) = blockset!.(y, x)
 blockset!(y::AbstractArray, x) = (y .= x)
 
-blockvecdot(x::T, y::T) where {T <: Tuple} = sum(blockvecdot.(x,y))
-blockvecdot(x::AbstractArray{R}, y::AbstractArray{R}) where {R <: Number} = real(vecdot(x, y))
+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} = real(vecdot(x, y))
 # inner product must be always real see section 4.2 of TFOCS manual
 
 blockzeros(t::Tuple, s::Tuple) = blockzeros.(t, s)
@@ -64,6 +65,13 @@ blockzeros(n::NTuple{N, Integer} where {N}) = zeros(n)
 blockzeros(n::Integer) = zeros(n)
 blockzeros(a::AbstractArray) = zeros(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)
+
 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)
 

test/test_block.jl

@@ -66,6 +66,21 @@ yb = blockzeros(blockeltype(xb),blocksize(xb))
 @test y == zeros(2)
 @test yb == (zeros(2),zeros(3)+im*zeros(3),zeros(2,3))
 
+y = blockones(x)
+yb = blockones(xb)
+@test y == ones(2)
+@test yb == (ones(2),ones(3)+im*zeros(3),ones(2,3))
+
+y = blockones(blocksize(x))
+yb = blockones(blocksize(xb))
+@test y == ones(2)
+@test yb == (ones(2),ones(3)+im*zeros(3),ones(2,3))
+
+y = blockones(blockeltype(x),blocksize(x))
+yb = blockones(blockeltype(xb),blocksize(xb))
+@test y == ones(2)
+@test yb == (ones(2),ones(3)+im*zeros(3),ones(2,3))
+
 blockaxpy!(y,x,2,x2)
 blockaxpy!(yb,xb,2,x2b)
 

test/test_lbfgs.jl

@@ -76,6 +76,17 @@ for i = 1:5
 
 end
 
+@test blockones(size(H,1)) != H*blockones(size(H,1))
+@test blockones(size(HH,1)) != HH*blockones(size(HH,1))
+
+#testing reset
+
+reset!(H)
+reset!(HH)
+
+@test blockones(size(H,1)) == H*blockones(size(H,1))
+@test blockones(size(HH,1)) == HH*blockones(size(HH,1))
+
 end
 
 test_lbfgs()

test/test_linear_operators.jl

@@ -282,6 +282,7 @@ y1 = test_op(op, x1, randn(n), verb)
 op = Eye(Float64, (n,))
 op = Eye((n,))
 op = Eye(n)
+op = Eye(x1)
 
 #properties
 @test is_linear(op)           == true
@@ -472,21 +473,47 @@ op = GetIndex(Float64,(n,),(1:k,))
 
 ######## MatrixOp ############
 
+# real matrix, real input
 n,m = 5,4
 A = randn(n,m)
 op = MatrixOp(Float64,(m,),A)
 x1 = randn(m)
 y1 = test_op(op, x1, randn(n), verb)
 y2 = A*x1
 
+# real matrix, complex input
+n,m = 5,4
+A = randn(n,m)
+op = MatrixOp(Complex{Float64},(m,),A)
+x1 = randn(m)+im.*randn(m)
+y1 = test_op(op, x1, randn(n)+im*randn(n), verb)
+y2 = A*x1
+
+# complex matrix, complex input
+n,m = 5,4
+A = randn(n,m)+im*randn(n,m)
+op = MatrixOp(Complex{Float64},(m,),A)
+x1 = randn(m)+im.*randn(m)
+y1 = test_op(op, x1, randn(n)+im*randn(n), verb)
+y2 = A*x1
+
+# complex matrix, real input
+n,m = 5,4
+A = randn(n,m)+im*randn(n,m)
+op = MatrixOp(Float64,(m,),A)
+x1 = randn(m)
+y1 = test_op(op, x1, randn(n)+im*randn(n), verb)
+y2 = A*x1
+
 @test all(vecnorm.(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))
 x1 = randn(m,c)
-y1 = test_op(op, x1, randn(n,c), verb)
+y1 = test_op(op, x1, randn(n,c).+randn(n,c), verb)
 y2 = A*x1
 
 # other constructors
@@ -497,6 +524,7 @@ op = MatrixOp(Float64, A, c)
 
 op = convert(LinearOperator,A)
 op = convert(LinearOperator,A,c)
+op = convert(LinearOperator, Complex{Float64}, size(x1), A)
 
 ##properties
 @test is_linear(op)           == true