diff --git a/docs/make.jl b/docs/make.jl index 233cc25..62c5b76 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -16,6 +16,6 @@ makedocs( ) deploydocs( - repo = "github.com/kul-forbes/AbstractOperators.jl.git", + repo = "github.com/kul-forbes/StructuredOptimization.jl.git", target = "build", ) diff --git a/src/solvers/build_solve.jl b/src/solvers/build_solve.jl index f9718b5..ac7a0ae 100644 --- a/src/solvers/build_solve.jl +++ b/src/solvers/build_solve.jl @@ -23,36 +23,36 @@ julia> build(p, PG()); """ function build(terms::Tuple, solver::ForwardBackwardSolver) - x = extract_variables(terms) - # Separate smooth and nonsmooth - smooth, nonsmooth = split_smooth(terms) - # Separate quadratic and nonquadratic - quadratic, smooth = split_quadratic(smooth) - kwargs = Array{Any, 1}() - if is_proximable(nonsmooth) - g = extract_proximable(x, nonsmooth) - append!(kwargs, [(:g, g)]) - if !isempty(quadratic) - fq = extract_functions(quadratic) - Aq = extract_operators(x, quadratic) - append!(kwargs, [(:fq, fq)]) - append!(kwargs, [(:Aq, Aq)]) - end - if !isempty(smooth) - if is_linear(smooth) - fs = extract_functions(smooth) - As = extract_operators(x, smooth) - append!(kwargs, [(:As, As)]) - else - fs = extract_functions_nodisp(smooth) - As = extract_affines(x, smooth) - fs = PrecomposeNonlinear(fs, As) - end - append!(kwargs, [(:fs, fs)]) - end - return build_iterator(x, solver; kwargs...) - end - error("Sorry, I cannot solve this problem") + x = extract_variables(terms) + # Separate smooth and nonsmooth + smooth, nonsmooth = split_smooth(terms) + # Separate quadratic and nonquadratic + quadratic, smooth = split_quadratic(smooth) + kwargs = Array{Any, 1}() + if is_proximable(nonsmooth) + g = extract_proximable(x, nonsmooth) + append!(kwargs, [(:g, g)]) + if !isempty(quadratic) + fq = extract_functions(quadratic) + Aq = extract_operators(x, quadratic) + append!(kwargs, [(:fq, fq)]) + append!(kwargs, [(:Aq, Aq)]) + end + if !isempty(smooth) + if is_linear(smooth) + fs = extract_functions(smooth) + As = extract_operators(x, smooth) + append!(kwargs, [(:As, As)]) + else + fs = extract_functions_nodisp(smooth) + As = extract_affines(x, smooth) + fs = PrecomposeNonlinear(fs, As) + end + append!(kwargs, [(:fs, fs)]) + end + return build_iterator(x, solver; kwargs...) + end + error("Sorry, I cannot solve this problem") end ################################################################################ @@ -83,10 +83,10 @@ julia> solve!(x_solver); """ function solve!(x_and_iter::Tuple{Tuple{Vararg{Variable}}, ProximalAlgorithms.ProximalAlgorithm}) - x, iterator = x_and_iter - it, x_star = ProximalAlgorithms.run!(iterator) - ~x .= x_star - return it, iterator + x, iterator = x_and_iter + it, x_star = ProximalAlgorithms.run!(iterator) + ~x .= x_star + return it, iterator end @@ -109,15 +109,15 @@ Variable(Float64, (4,)) julia> A, b = randn(10,4), randn(10); julia> solve(p,PG()); - it | gamma | fpr | - ------|------------|------------| - 1 | 7.6375e-02 | 1.8690e+00 | - 12 | 7.6375e-02 | 9.7599e-05 | +it | gamma | fpr | +------|------------|------------| +1 | 7.6375e-02 | 1.8690e+00 | +12 | 7.6375e-02 | 9.7599e-05 | ``` """ function solve(terms::Tuple, solver::ForwardBackwardSolver) - built_slv = build(terms, solver) - return solve!(built_slv) + built_slv = build(terms, solver) + return solve!(built_slv) end diff --git a/src/solvers/terms_extract.jl b/src/solvers/terms_extract.jl index b482557..389dea6 100644 --- a/src/solvers/terms_extract.jl +++ b/src/solvers/terms_extract.jl @@ -1,33 +1,33 @@ # returns all variables of a cost function, in terms of appearance -extract_variables(t::Term) = variables(t) - -function extract_variables(t::NTuple{N,Term}) where {N} - x = variables.(t) - xAll = x[1] - for i = 2:length(x) - for xi in x[i] - if (xi in xAll) == false - xAll = (xAll...,xi) - end - end - end - return xAll +extract_variables(t::TermOrExpr) = variables(t) + +function extract_variables(t::NTuple{N,TermOrExpr}) where {N} + x = variables.(t) + xAll = x[1] + for i = 2:length(x) + for xi in x[i] + if (xi in xAll) == false + xAll = (xAll...,xi) + end + end + end + return xAll end # extract functions from terms function extract_functions(t::Term) - f = displacement(t) == 0 ? t.f : PrecomposeDiagonal(t.f, 1.0, displacement(t)) #for now I keep this - f = t.lambda == 1. ? f : Postcompose(f, t.lambda) #for now I keep this - #TODO change this - return f + f = displacement(t) == 0 ? t.f : PrecomposeDiagonal(t.f, 1.0, displacement(t)) #for now I keep this + f = t.lambda == 1. ? f : Postcompose(f, t.lambda) #for now I keep this + #TODO change this + return f end extract_functions(t::NTuple{N,Term}) where {N} = SeparableSum(extract_functions.(t)) extract_functions(t::Tuple{Term}) = extract_functions(t[1]) # extract functions from terms without displacement function extract_functions_nodisp(t::Term) - f = t.lambda == 1. ? t.f : Postcompose(t.f, t.lambda) - return f + f = t.lambda == 1. ? t.f : Postcompose(t.f, t.lambda) + return f end extract_functions_nodisp(t::NTuple{N,Term}) where {N} = SeparableSum(extract_functions_nodisp.(t)) extract_functions_nodisp(t::Tuple{Term}) = extract_functions_nodisp(t[1]) @@ -37,29 +37,28 @@ extract_functions_nodisp(t::Tuple{Term}) = extract_functions_nodisp(t[1]) # returns all operators with an order dictated by xAll #single term, single variable -extract_operators(xAll::Tuple{Variable}, t::Term) = operator(t) - -extract_operators(xAll::NTuple{N,Variable}, t::Term) where {N} = extract_operators(xAll, (t,)) +extract_operators(xAll::Tuple{Variable}, t::TermOrExpr) = operator(t) +extract_operators(xAll::NTuple{N,Variable}, t::TermOrExpr) where {N} = extract_operators(xAll, (t,)) #multiple terms, multiple variables -function extract_operators(xAll::NTuple{N,Variable}, t::NTuple{M,Term}) where {N,M} - ops = () - for ti in t - tex = expand(xAll,ti) - ops = (ops...,sort_and_extract_operators(xAll,tex)) - end - return vcat(ops...) +function extract_operators(xAll::NTuple{N,Variable}, t::NTuple{M,TermOrExpr}) where {N,M} + ops = () + for ti in t + tex = expand(xAll,ti) + ops = (ops...,sort_and_extract_operators(xAll,tex)) + end + return vcat(ops...) end -sort_and_extract_operators(xAll::Tuple{Variable}, t::Term) = operator(t) +sort_and_extract_operators(xAll::Tuple{Variable}, t::TermOrExpr) = operator(t) -function sort_and_extract_operators(xAll::NTuple{N,Variable}, t::Term) where {N} - p = zeros(Int,N) - xL = variables(t) - for i in eachindex(xAll) - p[i] = findfirst( xi -> xi == xAll[i], xL) - end - return operator(t)[p] +function sort_and_extract_operators(xAll::NTuple{N,Variable}, t::TermOrExpr) where {N} + p = zeros(Int,N) + xL = variables(t) + for i in eachindex(xAll) + p[i] = findfirst( xi -> xi == xAll[i], xL) + end + return operator(t)[p] end # extract affines from terms @@ -67,100 +66,114 @@ end # returns all affines with an order dictated by xAll #single term, single variable -extract_affines(xAll::Tuple{Variable}, t::Term) = affine(t) +extract_affines(xAll::Tuple{Variable}, t::TermOrExpr) = affine(t) -extract_affines(xAll::NTuple{N,Variable}, t::Term) where {N} = extract_affines(xAll, (t,)) +extract_affines(xAll::NTuple{N,Variable}, t::TermOrExpr) where {N} = extract_affines(xAll, (t,)) #multiple terms, multiple variables -function extract_affines(xAll::NTuple{N,Variable}, t::NTuple{M,Term}) where {N,M} - ops = () - for ti in t - tex = expand(xAll,ti) - ops = (ops...,sort_and_extract_affines(xAll,tex)) - end - return vcat(ops...) +function extract_affines(xAll::NTuple{N,Variable}, t::NTuple{M,TermOrExpr}) where {N,M} + ops = () + for ti in t + tex = expand(xAll,ti) + ops = (ops...,sort_and_extract_affines(xAll,tex)) + end + return vcat(ops...) end -sort_and_extract_affines(xAll::Tuple{Variable}, t::Term) = affine(t) +sort_and_extract_affines(xAll::Tuple{Variable}, t::TermOrExpr) = affine(t) -function sort_and_extract_affines(xAll::NTuple{N,Variable}, t::Term) where {N} - p = zeros(Int,N) - xL = variables(t) - for i in eachindex(xAll) - p[i] = findfirst( xi -> xi == xAll[i], xL) - end - return affine(t)[p] +function sort_and_extract_affines(xAll::NTuple{N,Variable}, t::TermOrExpr) where {N} + p = zeros(Int,N) + xL = variables(t) + for i in eachindex(xAll) + p[i] = findfirst( xi -> xi == xAll[i], xL) + end + return affine(t)[p] end # expand term domain dimensions -function expand(xAll::NTuple{N,Variable}, t::Term{T1,T2,T3}) where {N,T1,T2,T3} - xt = variables(t) - C = codomainType(operator(t)) - size_out = size(operator(t),1) - ex = t.A - - for x in xAll - if !( x in variables(t) ) - ex += Zeros(eltype(~x),size(x),C,size_out)*x - end - end - return Term(t.lambda, t.f, ex) +function expand(xAll::NTuple{N,Variable}, t::Term) where {N} + xt = variables(t) + C = codomainType(operator(t)) + size_out = size(operator(t),1) + ex = t.A + + for x in xAll + if !( x in variables(t) ) + ex += Zeros(eltype(~x),size(x),C,size_out)*x + end + end + return Term(t.lambda, t.f, ex) +end + +function expand(xAll::NTuple{N,Variable}, ex::AbstractExpression) where {N} + ex = convert(Expression,ex) + xt = variables(ex) + C = codomainType(operator(ex)) + size_out = size(operator(ex),1) + + for x in xAll + if !( x in variables(ex) ) + ex += Zeros(eltype(~x),size(x),C,size_out)*x + end + end + return ex end # extract function and merge operator function extract_merge_functions(t::Term) - if is_sliced(t) - if typeof(operator(t)) <: Compose - op = operator(t).A[2] - else - op = Eye(size(operator(t),1)...) - end + if is_sliced(t) + if typeof(operator(t)) <: Compose + op = operator(t).A[2] else - op = operator(t) - end - if is_eye(op) - f = displacement(t) == 0 ? t.f : PrecomposeDiagonal(t.f, 1.0, displacement(t)) - elseif is_diagonal(op) - f = PrecomposeDiagonal(t.f, diag(op), displacement(t)) - elseif is_AAc_diagonal(op) - f = Precompose(t.f, op, diag_AAc(op), displacement(t)) + op = Eye(size(operator(t),1)...) end - f = t.lambda == 1. ? f : Postcompose(f, t.lambda) #for now I keep this - #TODO change this - return f + else + op = operator(t) + end + if is_eye(op) + f = displacement(t) == 0 ? t.f : PrecomposeDiagonal(t.f, 1.0, displacement(t)) + elseif is_diagonal(op) + f = PrecomposeDiagonal(t.f, diag(op), displacement(t)) + elseif is_AAc_diagonal(op) + f = Precompose(t.f, op, diag_AAc(op), displacement(t)) + end + f = t.lambda == 1. ? f : Postcompose(f, t.lambda) #for now I keep this + #TODO change this + return f end function extract_proximable(xAll::NTuple{N,Variable}, t::NTuple{M,Term}) where {N,M} - fs = () - for x in xAll - tx = () #terms containing x - for ti in t - if x in variables(ti) - tx = (tx...,ti) #collect terms containing x - end - end - if isempty(tx) - fx = IndFree() - elseif length(tx) == 1 #only one term per variable - fx = extract_proximable(x,tx[1]) - else - #multiple terms per variable - #currently this happens only with GetIndex - fxi,idxs = (),() - for ti in tx - fxi = (fxi..., extract_merge_functions(ti)) - idx = typeof(operator(ti)) <: Compose ? operator(ti).A[1].idx : operator(ti).idx - idxs = (idxs..., idx ) - end - fx = SlicedSeparableSum(fxi,idxs) - end - fs = (fs...,fx) - end - if length(fs) > 1 - return SeparableSum(fs) ##probably change constructor in Prox? - else - return fs[1] - end + fs = () + for x in xAll + tx = () #terms containing x + for ti in t + if x in variables(ti) + tx = (tx...,ti) #collect terms containing x + end + end + if isempty(tx) + fx = IndFree() + elseif length(tx) == 1 #only one term per variable + fx = extract_proximable(x,tx[1]) + else + #multiple terms per variable + #currently this happens only with GetIndex + fxi,idxs = (),() + for ti in tx + fxi = (fxi..., extract_merge_functions(ti)) + idx = typeof(operator(ti)) <: Compose ? operator(ti).A[1].idx : operator(ti).idx + idxs = (idxs..., idx ) + end + fx = SlicedSeparableSum(fxi,idxs) + end + fs = (fs...,fx) + end + if length(fs) > 1 + return SeparableSum(fs) ##probably change constructor in Prox? + else + return fs[1] + end end extract_proximable(xAll::Variable, t::Term) = extract_merge_functions(t) diff --git a/src/syntax/expressions/abstractOperator_bind.jl b/src/syntax/expressions/abstractOperator_bind.jl index 2ad7861..ed8977d 100644 --- a/src/syntax/expressions/abstractOperator_bind.jl +++ b/src/syntax/expressions/abstractOperator_bind.jl @@ -18,9 +18,9 @@ julia> reshape(A*x-b,2,5) """ function reshape(a::AbstractExpression, dims...) - A = convert(Expression,a) - op = Reshape(A.L, dims...) - return Expression{length(A.x)}(A.x,op) + A = convert(Expression,a) + op = Reshape(A.L, dims...) + return Expression{length(A.x)}(A.x,op) end #Reshape @@ -34,19 +34,19 @@ imported = [ ] importedFFTW = [ - :fft :(AbstractOperators.DFT); - :rfft :RDFT; - :irfft :IRDFT; - :ifft :IDFT; - :dct :DCT; - :idct :IDCT; - ] + :fft :(AbstractOperators.DFT); + :rfft :RDFT; + :irfft :IRDFT; + :ifft :IDFT; + :dct :DCT; + :idct :IDCT; + ] importedDSP = [ - :conv :Conv; - :xcorr :Xcorr; - :filt :Filt; - ] + :conv :Conv; + :xcorr :Xcorr; + :filt :Filt; + ] exported = [ :finitediff :FiniteDiff; @@ -60,41 +60,41 @@ exported = [ #importing functions from Base for f in imported[:,1] - @eval begin - import Base: $f - end + @eval begin + import Base: $f + end end #importing functions from FFTW for f in importedFFTW[:,1] - @eval begin - import FFTW: $f - export $f - end + @eval begin + import FFTW: $f + export $f + end end #importing functions from DSP for f in importedDSP[:,1] - @eval begin - import DSP: $f - export $f - end + @eval begin + import DSP: $f + export $f + end end #exporting functions for f in exported[:,1] - @eval begin - export $f - end + @eval begin + export $f + end end fun = [imported; importedFFTW; importedDSP; exported] for i = 1:size(fun,1) - f,fAbsOp = fun[i,1],fun[i,2] - @eval begin - function $f(a::AbstractExpression, args...) - A = convert(Expression,a) - op = $fAbsOp(codomainType(operator(A)),size(operator(A),1), args...) - return op*A - end - end + f,fAbsOp = fun[i,1],fun[i,2] + @eval begin + function $f(a::AbstractExpression, args...) + A = convert(Expression,a) + op = $fAbsOp(codomainType(operator(A)),size(operator(A),1), args...) + return op*A + end + end end ## docs diff --git a/src/syntax/expressions/addition.jl b/src/syntax/expressions/addition.jl index f786fbc..4756265 100644 --- a/src/syntax/expressions/addition.jl +++ b/src/syntax/expressions/addition.jl @@ -45,110 +45,110 @@ julia> ex3.+z """ function (+)(a::AbstractExpression, b::AbstractExpression) - A = convert(Expression,a) - B = convert(Expression,b) - if variables(A) == variables(B) + A = convert(Expression,a) + B = convert(Expression,b) + if variables(A) == variables(B) return Expression{length(A.x)}(A.x,affine(A)+affine(B)) - else + else opA = affine(A) xA = variables(A) opB = affine(B) xB = variables(B) xNew, opNew = Usum_op(xA,xB,opA,opB,true) return Expression{length(xNew)}(xNew,opNew) - end + end end # sum expressions function (-)(a::AbstractExpression, b::AbstractExpression) - A = convert(Expression,a) - B = convert(Expression,b) - if variables(A) == variables(B) + A = convert(Expression,a) + B = convert(Expression,b) + if variables(A) == variables(B) return Expression{length(A.x)}(A.x,affine(A)-affine(B)) - else + else opA = affine(A) xA = variables(A) opB = affine(B) xB = variables(B) xNew, opNew = Usum_op(xA,xB,opA,opB,false) return Expression{length(xNew)}(xNew,opNew) - end + end end #unsigned sum affines with single variables function Usum_op(xA::Tuple{Variable}, - xB::Tuple{Variable}, - A::AbstractOperator, - B::AbstractOperator,sign::Bool) - xNew = (xA...,xB...) - opNew = sign ? hcat(A,B) : hcat(A,-B) - return xNew, opNew + xB::Tuple{Variable}, + A::AbstractOperator, + B::AbstractOperator,sign::Bool) + xNew = (xA...,xB...) + opNew = sign ? hcat(A,B) : hcat(A,-B) + return xNew, opNew end #unsigned sum: HCAT + AbstractOperator function Usum_op(xA::NTuple{N,Variable}, - xB::Tuple{Variable}, - A::L1, - B::AbstractOperator,sign::Bool) where {N, M, L1<:HCAT{N}} - if xB[1] in xA + xB::Tuple{Variable}, + A::L1, + B::AbstractOperator,sign::Bool) where {N, M, L1<:HCAT{N}} + if xB[1] in xA idx = findfirst(xA.==Ref(xB[1])) S = sign ? A[idx]+B : A[idx]-B xNew = xA opNew = hcat(A[1:idx-1],S,A[idx+1:N] ) - else + else xNew = (xA...,xB...) opNew = sign ? hcat(A,B) : hcat(A,-B) - end - return xNew, opNew + end + return xNew, opNew end #unsigned sum: AbstractOperator+HCAT function Usum_op(xA::Tuple{Variable}, - xB::NTuple{N,Variable}, - A::AbstractOperator, - B::L2,sign::Bool) where {N, M, L2<:HCAT{N}} + xB::NTuple{N,Variable}, + A::AbstractOperator, + B::L2,sign::Bool) where {N, M, L2<:HCAT{N}} if xA[1] in xB idx = findfirst(xA.==Ref(xB[1])) S = sign ? A+B[idx] : B[idx]-A xNew = xB opNew = sign ? hcat(B[1:idx-1],S,B[idx+1:N] ) : -hcat(B[1:idx-1],S,B[idx+1:N] ) - else + else xNew = (xA...,xB...) opNew = sign ? hcat(A,B) : hcat(A,-B) - end + end - return xNew, opNew + return xNew, opNew end #unsigned sum: HCAT+HCAT function Usum_op(xA::NTuple{NA,Variable}, - xB::NTuple{NB,Variable}, - A::L1, - B::L2,sign::Bool) where {NA,NB,M, - L1<:HCAT{NB}, - L2<:HCAT{NB} } - xNew = xA - opNew = A - for i in eachindex(xB) + xB::NTuple{NB,Variable}, + A::L1, + B::L2,sign::Bool) where {NA,NB,M, + L1<:HCAT{NB}, + L2<:HCAT{NB} } + xNew = xA + opNew = A + for i in eachindex(xB) xNew, opNew = Usum_op(xNew, (xB[i],), opNew, B[i], sign) - end -return xNew,opNew + end + return xNew,opNew end #unsigned sum: multivar AbstractOperator + AbstractOperator function Usum_op(xA::NTuple{N,Variable}, - xB::Tuple{Variable}, - A::AbstractOperator, - B::AbstractOperator,sign::Bool) where {N} - if xB[1] in xA + xB::Tuple{Variable}, + A::AbstractOperator, + B::AbstractOperator,sign::Bool) where {N} + if xB[1] in xA Z = Zeros(A) #this will be an HCAT xNew, opNew = Usum_op(xA,xB,Z,B,sign) opNew += A - else + else xNew = (xA...,xB...) opNew = sign ? hcat(A,B) : hcat(A,-B) - end - return xNew, opNew + end + return xNew, opNew end """ @@ -184,19 +184,19 @@ julia> ex + b """ function (+)(a::AbstractExpression, b::Union{AbstractArray,Number}) - A = convert(Expression,a) + A = convert(Expression,a) return Expression{length(A.x)}(A.x,AffineAdd(affine(A),b)) end (+)(a::Union{AbstractArray,Number}, b::AbstractExpression) = b+a function (-)(a::AbstractExpression, b::Union{AbstractArray,Number}) - A = convert(Expression,a) + A = convert(Expression,a) return Expression{length(A.x)}(A.x,AffineAdd(affine(A),b,false)) end function (-)(a::Union{AbstractArray,Number}, b::AbstractExpression) - B = convert(Expression,b) + B = convert(Expression,b) return Expression{length(B.x)}(B.x,-AffineAdd(affine(B),a)) end # sum with array/scalar @@ -204,9 +204,9 @@ end #broadcasted + - function Broadcast.broadcasted(::typeof(+),a::AbstractExpression, b::AbstractExpression) - A = convert(Expression,a) - B = convert(Expression,b) - if size(affine(A),1) != size(affine(B),1) + A = convert(Expression,a) + B = convert(Expression,b) + if size(affine(A),1) != size(affine(B),1) if prod(size(affine(A),1)) > prod(size(affine(B),1)) B = Expression{length(B.x)}(variables(B), BroadCast(affine(B),size(affine(A),1))) @@ -215,14 +215,14 @@ function Broadcast.broadcasted(::typeof(+),a::AbstractExpression, b::AbstractExp BroadCast(affine(A),size(affine(B),1))) end return A+B - end - return A+B + end + return A+B end function Broadcast.broadcasted(::typeof(-),a::AbstractExpression, b::AbstractExpression) - A = convert(Expression,a) - B = convert(Expression,b) - if size(affine(A),1) != size(affine(B),1) + A = convert(Expression,a) + B = convert(Expression,b) + if size(affine(A),1) != size(affine(B),1) if prod(size(affine(A),1)) > prod(size(affine(B),1)) B = Expression{length(B.x)}(variables(B), BroadCast(affine(B),size(affine(A),1))) @@ -231,6 +231,6 @@ function Broadcast.broadcasted(::typeof(-),a::AbstractExpression, b::AbstractExp BroadCast(affine(A),size(affine(B),1))) end return A-B - end - return A-B + end + return A-B end diff --git a/src/syntax/expressions/expression.jl b/src/syntax/expressions/expression.jl index 14ad47c..08d1f53 100644 --- a/src/syntax/expressions/expression.jl +++ b/src/syntax/expressions/expression.jl @@ -1,28 +1,36 @@ struct Expression{N,A<:AbstractOperator} <: AbstractExpression - x::NTuple{N,Variable} - L::A - function Expression{N}(x::NTuple{N,Variable}, L::A) where {N,A<:AbstractOperator} - - # checks on L - ndoms(L,1) > 1 && throw(ArgumentError( - "cannot create expression with LinearOperator with `ndoms(L,1) > 1`")) + x::NTuple{N,Variable} + L::A + function Expression{N}(x::NTuple{N,Variable}, L::A) where {N,A<:AbstractOperator} + # checks on L + ndoms(L,1) > 1 && throw(ArgumentError( + "Cannot create expression with LinearOperator with `ndoms(L,1) > 1`" + )) + #checks on x + szL = size(L,2) + szx = size.(x) + check_sz = length(szx) == 1 ? szx[1] != szL : szx != szL + check_sz && throw(ArgumentError( + "Size of the operator domain $(size(L, 2)) must match size of the variable $(size.(x))" + )) + dmL = domainType(L) + dmx = eltype.(x) + check_dm = length(dmx) == 1 ? dmx[1] != dmL : dmx != dmL + check_dm && throw(ArgumentError( + "Type of the operator domain $(domainType(L)) must match type of the variable $(eltype.(x))" + )) + new{N,A}(x,L) + end +end - #checks on x - szL = size(L,2) - szx = size.(x) - check_sz = length(szx) == 1 ? szx[1] != szL : szx != szL - check_sz && throw(ArgumentError( - "Size of the operator domain $(size(L, 2)) must match size of the variable $(size.(x))")) +struct AdjointExpression{E <: AbstractExpression} <: AbstractExpression + ex::E +end - dmL = domainType(L) - dmx = eltype.(x) - check_dm = length(dmx) == 1 ? dmx[1] != dmL : dmx != dmL - check_dm && throw(ArgumentError( - "Type of the operator domain $(domainType(L)) must match type of the variable $(eltype.(x))")) +import Base: adjoint - new{N,A}(x,L) - end -end +adjoint(ex::AbstractExpression) = AdjointExpression(convert(Expression,ex)) +adjoint(ex::AdjointExpression) = ex.ex include("utils.jl") include("multiplication.jl") diff --git a/src/syntax/expressions/multiplication.jl b/src/syntax/expressions/multiplication.jl index db850a8..bbab925 100644 --- a/src/syntax/expressions/multiplication.jl +++ b/src/syntax/expressions/multiplication.jl @@ -26,8 +26,8 @@ julia> affine(ex2) """ function (*)(L::AbstractOperator, a::AbstractExpression) - A = convert(Expression,a) - Expression{length(A.x)}(A.x,L*affine(A)) + A = convert(Expression,a) + Expression{length(A.x)}(A.x,L*affine(A)) end """ @@ -48,53 +48,53 @@ julia> A*x Other types of multiplications are also possible: * left array multiplication - ```julia - julia> X = Variable(10,5) - Variable(Float64, (10, 5)) - - julia> X*randn(5,6) - - ``` +```julia +julia> X = Variable(10,5) +Variable(Float64, (10, 5)) + +julia> X*randn(5,6) + +``` * scalar multiplication: - ```julia - julia> π*X - - ``` +```julia +julia> π*X + +``` * elementwise multiplication: - ```julia - julia> randn(10,5).*X - - ``` +```julia +julia> randn(10,5).*X + +``` """ function (*)(m::T, a::Union{AbstractVector,AbstractMatrix}) where {T<:AbstractExpression} - M = convert(Expression,m) - op = LMatrixOp(codomainType(affine(M)),size(affine(M),1),a) - return op*M + M = convert(Expression,m) + op = LMatrixOp(codomainType(affine(M)),size(affine(M),1),a) + return op*M end #LMatrixOp function (*)(M::AbstractMatrix, a::T) where {T<:AbstractExpression} - A = convert(Expression,a) - op = MatrixOp(codomainType(affine(A)),size(affine(A),1),M) - return op*A + A = convert(Expression,a) + op = MatrixOp(codomainType(affine(A)),size(affine(A),1),M) + return op*A end #MatrixOp function Broadcast.broadcasted(::typeof(*), d::D, a::T) where {D <: Union{Number,AbstractArray}, T<:AbstractExpression} - A = convert(Expression,a) - op = DiagOp(codomainType(affine(A)),size(affine(A),1),d) - return op*A + A = convert(Expression,a) + op = DiagOp(codomainType(affine(A)),size(affine(A),1),d) + return op*A end Broadcast.broadcasted(::typeof(*), a::T, d::D) where {D <: Union{Number,AbstractArray}, T<:AbstractExpression} = d.*a #DiagOp function (*)(coeff::T1, a::T) where {T1<:Number, T<:AbstractExpression} - A = convert(Expression,a) - return Expression{length(A.x)}(A.x,coeff*affine(A)) + A = convert(Expression,a) + return Expression{length(A.x)}(A.x,coeff*affine(A)) end (*)(a::T, coeff::T1) where {T1<:Number, T<:AbstractExpression} = coeff*a ##Scale @@ -126,27 +126,49 @@ Elementwise multiplication between `AbstractExpression` (i.e. Hadamard product). """ function (*)(ex1::AbstractExpression, ex2::AbstractExpression) - ex1 = convert(Expression,ex1) - ex2 = convert(Expression,ex2) - if any([x in variables(ex2) for x in variables(ex1)]) - error("cannot muliply expressions containing the same variable") - end - op = NonLinearCompose(affine(ex1),affine(ex2)) - x = (variables(ex1)...,variables(ex2)...) - exp3 = Expression{length(x)}(x,op) - return exp3 + ex1 = convert(Expression,ex1) + ex2 = convert(Expression,ex2) + x = extract_variables((ex1,ex2)) + A = extract_affines(x, ex1) + B = extract_affines(x, ex2) + op = Ax_mul_Bx(A,B) + exp3 = Expression{length(x)}(x,op) + return exp3 +end +# Ax_mul_Bx + +function (*)(ex1::AdjointExpression, ex2::AbstractExpression) + ex1 = ex1.ex + ex2 = convert(Expression,ex2) + x = extract_variables((ex1,ex2)) + A = extract_affines(x, ex1) + B = extract_affines(x, ex2) + op = Axt_mul_Bx(A,B) + exp3 = Expression{length(x)}(x,op) + return exp3 +end +# Axt_mul_Bx + +function (*)(ex1::AbstractExpression, ex2::AdjointExpression) + ex1 = convert(Expression,ex1) + ex2 = ex2.ex + x = extract_variables((ex1,ex2)) + A = extract_affines(x, ex1) + B = extract_affines(x, ex2) + op = Ax_mul_Bxt(A,B) + exp3 = Expression{length(x)}(x,op) + return exp3 end -# NonLinearCompose +# Ax_mul_Bxt function Broadcast.broadcasted(::typeof(*), ex1::AbstractExpression, ex2::AbstractExpression) - ex1 = convert(Expression,ex1) - ex2 = convert(Expression,ex2) - if any([x in variables(ex2) for x in variables(ex1)]) - error("cannot muliply expressions containing the same variable") - end - op = Hadamard(affine(ex1),affine(ex2)) - x = (variables(ex1)...,variables(ex2)...) - exp3 = Expression{length(x)}(x,op) - return exp3 + ex1 = convert(Expression,ex1) + ex2 = convert(Expression,ex2) + x = extract_variables((ex1,ex2)) + A = extract_affines(x, ex1) + B = extract_affines(x, ex2) + op = HadamardProd(A,B) + exp3 = Expression{length(x)}(x,op) + return exp3 end # Hadamard diff --git a/src/syntax/expressions/utils.jl b/src/syntax/expressions/utils.jl index da6957f..9e82cb1 100644 --- a/src/syntax/expressions/utils.jl +++ b/src/syntax/expressions/utils.jl @@ -72,9 +72,9 @@ julia> ex = fft(x)+[1.+im*2.;0.;3.+im*4]; julia> displacement(ex) 3-element Array{Complex{Float64},1}: - 1.0+2.0im - 0.0+0.0im - 3.0+4.0im +1.0+2.0im +0.0+0.0im +3.0+4.0im ``` """ diff --git a/src/syntax/syntax.jl b/src/syntax/syntax.jl index 51fa9e7..514514b 100644 --- a/src/syntax/syntax.jl +++ b/src/syntax/syntax.jl @@ -4,3 +4,5 @@ include("variable.jl") include("expressions/expression.jl") include("terms/term.jl") include("problem.jl") + +const TermOrExpr = Union{Term,AbstractExpression} diff --git a/test/test_expressions.jl b/test/test_expressions.jl index f082c24..890786d 100644 --- a/test/test_expressions.jl +++ b/test/test_expressions.jl @@ -1,4 +1,9 @@ println("\nTesting linear expressions\n") +### AdjointExpression + +x1 = Variable(randn(2)) +@test typeof(x1') <: StructuredOptimization.AdjointExpression +@test typeof((x1')') <: StructuredOptimization.Expression #### * #### n, m1, m2, k = 3, 4, 5, 6 @@ -43,7 +48,7 @@ b2 = randn(n,n) opA1 = MatrixOp(A1,n) opA2 = MatrixOp(A2,n) x1, x2 = Variable(randn(m1,n)), Variable(randn(m2,n)) -# multiply Expressions (NonLinearCompose) +# multiply Expressions (Ax_mul_Bx) ex = (opA1*x1)*(opA2*x2) @test variables(ex) == (x1,x2) @test norm(affine(ex)*(~variables(ex)) - (A1*(~x1))*(A2*(~x2))) < 1e-12 @@ -56,7 +61,23 @@ ex = (opA1*x1)*(opA2*x2+b2) ex = (opA1*x1+b1)*(opA2*x2+b2) @test variables(ex) == (x1,x2) @test norm(affine(ex)*(~variables(ex)) - (A1*(~x1)+b1)*(A2*(~x2)+b2)) < 1e-12 -@test_throws ErrorException (opA1*x1)*(opA1*x1) +ex = (opA1*x1-b1)*(opA1*x1+b1) +@test variables(ex) == (x1,) +@test norm(affine(ex)*(~variables(ex)) - (A1*(~x1)-b1)*(A1*(~x1)+b1)) < 1e-12 +# multiply Expressions (Axt_mul_Bx) +ex = (opA1*x1)'*(opA2*x2) +@test variables(ex) == (x1,x2) +@test norm(affine(ex)*(~variables(ex)) - (A1*(~x1))'*(A2*(~x2))) < 1e-12 +ex = (opA1*x1-b1)'*(opA1*x1+b1) +@test variables(ex) == (x1,) +@test norm(affine(ex)*(~variables(ex)) - (A1*(~x1)-b1)'*(A1*(~x1)+b1)) < 1e-12 +# multiply Expressions (Ax_mul_Bxt) +ex = (opA1*x1)*(opA2*x2)' +@test variables(ex) == (x1,x2) +@test norm(affine(ex)*(~variables(ex)) - (A1*(~x1))*(A2*(~x2))') < 1e-12 +ex = (opA1*x1-b1)*(opA1*x1+b1)' +@test variables(ex) == (x1,) +@test norm(affine(ex)*(~variables(ex)) - (A1*(~x1)-b1)*(A1*(~x1)+b1)') < 1e-12 n, m1, m2, k = 3, 4, 5, 6 A1 = randn(n, m1) @@ -79,7 +100,9 @@ ex = (opA1*x1).*(opA2*x2+b2) ex = (opA1*x1+b1).*(opA2*x2+b2) @test variables(ex) == (x1,x2) @test norm(affine(ex)*(~variables(ex)) - (A1*(~x1)+b1).*(A2*(~x2)+b2)) < 1e-12 -@test_throws ErrorException (opA1*x1)*(opA1*x1) +ex = (opA1*x1-b1).*(opA1*x1+b1) +@test variables(ex) == (x1,) +@test norm(affine(ex)*(~variables(ex)) - (A1*(~x1)-b1).*(A1*(~x1)+b1)) < 1e-12 ##### reshape #### m,n = 8,10 @@ -292,3 +315,4 @@ ex3 = ex1-ex2 @test_throws DimensionMismatch MatrixOp(randn(10,20))*Variable(20)+randn(11) @test_throws ErrorException MatrixOp(randn(10,20))*Variable(20)+(3+im) + diff --git a/test/test_proxstuff.jl b/test/test_proxstuff.jl index 2e8e5e9..c4ce361 100644 --- a/test/test_proxstuff.jl +++ b/test/test_proxstuff.jl @@ -20,12 +20,16 @@ grad_f_x_ref = 3.0 * ( expx ./ (1 .+ expx).^2 ) .* (1.0 ./ (1.0 .+ expmx) - b) #with vectors l,m1,m2,n1,n2 = 2,3,4,5,6 x = ArrayPartition(randn(m1,m2),randn(n1,n2)) -A = randn(l,m1) -B = randn(m2,n1) +A = MatrixOp(randn(l,m1),m2) +B = MatrixOp(randn(m2,n1),n2) r = randn(l,n2) b = randn(l,n2) -G = AffineAdd(NonLinearCompose( MatrixOp(A,m2), MatrixOp(B,n2) ),b,false) +G = AffineAdd(Ax_mul_Bx( + HCAT(A,Zeros(codomainType(B), size(B,2), size(A,1) )), + HCAT(Zeros(codomainType(A), size(A,2), size(B,1) ),B) + ), + b,false) g = SqrNormL2(3.0) f = StructuredOptimization.PrecomposeNonlinear(g, G)