From 11966a5fdb1792a08a5c7f859140b94a8d7fdf9d Mon Sep 17 00:00:00 2001 From: Mathias Louboutin Date: Thu, 5 Nov 2020 09:12:13 -0500 Subject: [PATCH 1/3] Fix JOLI support in Minkowski sets --- src/PARSDMM.jl | 4 +- src/PARSDMM_initialize.jl | 4 +- src/PARSDMM_precompute_distribute.jl | 2 +- ...PARSDMM_precompute_distribute_Minkowski.jl | 42 ++++-- src/Q_update!.jl | 2 +- src/SetIntersectionProjection.jl | 3 +- src/compute_relative_feasibility.jl | 8 +- src/mat2CDS.jl | 2 +- src/projectors/project_l1_Duchi!.jl | 137 +++++------------- src/rhs_compose.jl | 2 +- src/setup_constraints.jl | 2 +- src/setup_multi_level_PARSDMM.jl | 4 +- src/update_y_l.jl | 2 +- src/update_y_l_parallel.jl | 2 +- 14 files changed, 81 insertions(+), 135 deletions(-) diff --git a/src/PARSDMM.jl b/src/PARSDMM.jl index 458d73b..2ed9ff3 100755 --- a/src/PARSDMM.jl +++ b/src/PARSDMM.jl @@ -23,8 +23,8 @@ output: - log_PARSDMM - constains log information about various quantities per iteration """ function PARSDMM(m ::Vector{TF}, - AtA ::Union{Vector{SparseMatrixCSC{TF,TI}},Vector{Array{TF,2}}}, - TD_OP ::Union{Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joLinearFunction{TF,TF}}},DistributedArrays.DArray{Union{JOLI.joLinearFunction{TF,TF}, SparseMatrixCSC{TF,TI}},1,Array{Union{JOLI.joLinearFunction{TF,TF}, SparseMatrixCSC{TF,TI}},1}} }, + AtA ::Union{Vector{Array{TF, 2}}, Vector{SparseMatrixCSC{TF, TI}}, Vector{joAbstractLinearOperator{TF, TF}}, Vector{Union{Array{TF, 2}, SparseMatrixCSC{TF, TI}, joAbstractLinearOperator{TF, TF}}}}, + TD_OP ::Union{Vector{Union{SparseMatrixCSC{TF,TI},joAbstractLinearOperator{TF,TF}}},DistributedArrays.DArray{Union{JOLI.joAbstractLinearOperator{TF,TF}, SparseMatrixCSC{TF,TI}},1,Array{Union{joAbstractLinearOperator{TF,TF}, SparseMatrixCSC{TF,TI}},1}} }, set_Prop, P_sub ::Union{Vector{Any},DistributedArrays.DArray{Any,1,Array{Any,1}}}, comp_grid, diff --git a/src/PARSDMM_initialize.jl b/src/PARSDMM_initialize.jl index bf2971a..889f450 100755 --- a/src/PARSDMM_initialize.jl +++ b/src/PARSDMM_initialize.jl @@ -7,8 +7,8 @@ function PARSDMM_initialize( x ::Vector{TF}, l ::Union{Vector{Vector{TF}},DistributedArrays.DArray{Array{TF,1},1,Array{Array{TF,1},1}},Array{Any,1}}, y ::Union{Vector{Vector{TF}},DistributedArrays.DArray{Array{TF,1},1,Array{Array{TF,1},1}},Array{Any,1}}, - AtA ::Union{Vector{SparseMatrixCSC{TF,TI}},Vector{Array{TF,2}}}, - TD_OP ::Union{Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joLinearFunction{TF,TF}}},DistributedArrays.DArray{Union{JOLI.joLinearFunction{TF,TF}, SparseMatrixCSC{TF,TI}},1,Array{Union{JOLI.joLinearFunction{TF,TF}, SparseMatrixCSC{TF,TI}},1}} }, + AtA ::Union{Vector{Array{TF, 2}}, Vector{SparseMatrixCSC{TF, TI}}, Vector{joAbstractLinearOperator{TF, TF}}, Vector{Union{Array{TF, 2}, SparseMatrixCSC{TF, TI}, joAbstractLinearOperator{TF, TF}}}}, + TD_OP ::Union{Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joAbstractLinearOperator{TF,TF}}},DistributedArrays.DArray{Union{JOLI.joAbstractLinearOperator{TF,TF}, SparseMatrixCSC{TF,TI}},1,Array{Union{JOLI.joAbstractLinearOperator{TF,TF}, SparseMatrixCSC{TF,TI}},1}} }, set_Prop, P_sub ::Union{Vector{Any},DistributedArrays.DArray{Any,1,Array{Any,1}}}, comp_grid, diff --git a/src/PARSDMM_precompute_distribute.jl b/src/PARSDMM_precompute_distribute.jl index d77370f..872f3c7 100644 --- a/src/PARSDMM_precompute_distribute.jl +++ b/src/PARSDMM_precompute_distribute.jl @@ -4,7 +4,7 @@ export PARSDMM_precompute_distribute Precomputes and distributes some quantities that serve as input for PARSDMM.jl """ function PARSDMM_precompute_distribute( - TD_OP ::Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joLinearFunction{TF,TF}}}, + TD_OP ::Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joAbstractLinearOperator{TF,TF}}}, set_Prop, comp_grid, options diff --git a/src/PARSDMM_precompute_distribute_Minkowski.jl b/src/PARSDMM_precompute_distribute_Minkowski.jl index 854bb33..e105dfc 100644 --- a/src/PARSDMM_precompute_distribute_Minkowski.jl +++ b/src/PARSDMM_precompute_distribute_Minkowski.jl @@ -1,9 +1,9 @@ export PARSDMM_precompute_distribute_Minkowski function PARSDMM_precompute_distribute_Minkowski( - TD_OP_c1 ::Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joLinearFunction{TF,TF}}}, - TD_OP_c2 ::Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joLinearFunction{TF,TF}}}, - TD_OP_sum ::Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joLinearFunction{TF,TF}}}, + TD_OP_c1::Vector{Union{SparseMatrixCSC{TF,TI}, joAbstractLinearOperator{TF,TF}}}, + TD_OP_c2::Vector{Union{SparseMatrixCSC{TF,TI}, joAbstractLinearOperator{TF,TF}}}, + TD_OP_sum::Vector{Union{SparseMatrixCSC{TF,TI}, joAbstractLinearOperator{TF,TF}}}, set_Prop_c1, set_Prop_c2, set_Prop_sum, @@ -22,9 +22,12 @@ else s=p+q+r+1 end -AtA=Vector{SparseMatrixCSC{TF,TI}}(undef,s) -ZERO_MAT = spzeros(TF,N,N) -Id_MAT = SparseMatrixCSC{TF}(LinearAlgebra.I,N,N) +joli_op = false + +AtA = Vector{Union{Array{TF, 2}, SparseMatrixCSC{TF, TI}, joAbstractLinearOperator{TF, TF}}}(undef, s) +ZERO_MAT_J = joZeros(N, N; DDT=TF, RDT=TF) +ZERO_MAT = spzeros(TF, N, N) +Id_MAT = SparseMatrixCSC{TF}(LinearAlgebra.I, N, N) for i=1:p #first component if set_Prop_c1.dense[i]==true @@ -35,7 +38,9 @@ for i=1:p #first component error("provided a dense non orthogoal transform-domain operator") end else - AtA[i] = [ TD_OP_c1[i]'*TD_OP_c1[i] ZERO_MAT ; ZERO_MAT ZERO_MAT ] + typeof(TD_OP_c1[i]) <: joAbstractLinearOperator ? ZM = ZERO_MAT_J : ZM = ZERO_MAT + joli_op = joli_op || (typeof(TD_OP_c1[i]) <: joAbstractLinearOperator) + AtA[i] = [ TD_OP_c1[i]'*TD_OP_c1[i] ZM ; ZM ZM ] #set_Prop_c1.AtA_offsets[i] #do not need to add additional offset info end end @@ -48,7 +53,9 @@ for i=1:q #second component error("provided a dense non orthogoal transform-domain operator") end else - AtA[p+i] = [ ZERO_MAT ZERO_MAT ; ZERO_MAT TD_OP_c2[i]'*TD_OP_c2[i] ] + typeof(TD_OP_c2[i]) <: joAbstractLinearOperator ? ZM = ZERO_MAT_J : ZM = ZERO_MAT + joli_op = joli_op || (typeof(TD_OP_c1[i]) <: joAbstractLinearOperator) + AtA[p+i] = [ ZM ZM ; ZM TD_OP_c2[i]'*TD_OP_c2[i] ] #do not need to add additional offset info end end @@ -69,10 +76,12 @@ end #modify the transform-domain operators to include two blocks per operator, i.e., # TD_OP = [A] -> TD_OP = [A 0] or [0 A] or [A A]. This is one row of \tilde{A} for i=1:length(TD_OP_c1) - TD_OP_c1[i] = [ TD_OP_c1[i] spzeros(TF,size(TD_OP_c1[i],1),N) ]; + typeof(TD_OP_c1[i])<:joAbstractLinearOperator ? ZM = joZeros(size(TD_OP_c1[i],1), N; DDT=TF, RDT=TF) : ZM = spzeros(TF, size(TD_OP_c1[i],1), N) + TD_OP_c1[i] = [ TD_OP_c1[i] ZM ]; end for i=1:length(TD_OP_c2) - TD_OP_c2[i] = [ spzeros(TF,size(TD_OP_c2[i],1),N) TD_OP_c2[i] ]; + typeof(TD_OP_c2[i])<:joAbstractLinearOperator ? ZM = joZeros(size(TD_OP_c2[i],1), N; DDT=TF, RDT=TF) : ZM = spzeros(TF, size(TD_OP_c2[i],1), N) + TD_OP_c2[i] = [ ZM TD_OP_c2[i] ]; end for i=1:length(TD_OP_sum) TD_OP_sum[i] = [ TD_OP_sum[i] TD_OP_sum[i] ]; @@ -111,21 +120,22 @@ append!(set_Prop.ncvx,set_Prop_sum.ncvx) append!(set_Prop.tag,set_Prop_sum.tag) #if all AtA are banded -> convert to compressed diagonal storage (CDS/DIA) format -if sum(set_Prop.banded[1:s].=true)==s +if sum(set_Prop.banded[1:s].=true)==s && ~joli_op for i=1:s + (AtA[i],set_Prop.AtA_offsets[i]) = mat2CDS(AtA[i]) set_Prop.AtA_offsets[i]=convert(Vector{TI},set_Prop.AtA_offsets[i]) end - AtA=convert(Vector{Array{TF,2}},AtA); - set_Prop.AtA_offsets=set_Prop.AtA_offsets[1:s] + AtA=convert(Vector{Array{TF,2}}, AtA); + set_Prop.AtA_offsets=set_Prop.AtA_offsets[1:s] end #allocate arrays of vectors -y = Vector{Vector{TF}}(undef,s); -l = Vector{Vector{TF}}(undef,s); +y = Vector{Vector{TF}}(undef,s); +l = Vector{Vector{TF}}(undef,s); #create one TD_OP that contains all transform-domain operators -TD_OP = Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joLinearFunction{TF,TF}}}(undef,s) +TD_OP = Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joAbstractLinearOperator{TF,TF}}}(undef,s) TD_OP[1:p] = TD_OP_c1 TD_OP[1+p:p+q] = TD_OP_c2 TD_OP[1+p+q:s] = TD_OP_sum diff --git a/src/Q_update!.jl b/src/Q_update!.jl index 1775fc4..45a2bcc 100755 --- a/src/Q_update!.jl +++ b/src/Q_update!.jl @@ -6,7 +6,7 @@ if any of the rho[i] change """ function Q_update!( Q ::Union{Array{TF,2},SparseMatrixCSC{TF,TI}}, - AtA ::Union{Vector{Array{TF,2}},Vector{SparseMatrixCSC{TF,TI}}}, + AtA ::Union{Vector{Array{TF, 2}}, Vector{SparseMatrixCSC{TF, TI}}, Vector{joAbstractLinearOperator{TF, TF}}, Vector{Union{Array{TF, 2}, SparseMatrixCSC{TF, TI}, joAbstractLinearOperator{TF, TF}}}}, set_Prop, rho ::Vector{TF}, ind_updated ::Vector{TI}, diff --git a/src/SetIntersectionProjection.jl b/src/SetIntersectionProjection.jl index 678b7a3..848bcc6 100644 --- a/src/SetIntersectionProjection.jl +++ b/src/SetIntersectionProjection.jl @@ -74,6 +74,7 @@ include("projectors/project_histogram_relaxed.jl"); #other proximal maps include("prox_l2s!.jl") +include("prox_l1!.jl") #scripts that are required to run examples include("constraint_learning_by_observation.jl") @@ -141,7 +142,7 @@ mutable struct set_definitions min ::Union{Vector,Real} max ::Union{Vector,Real} app_mode ::Tuple{String,String} #("tensor/"matrix"/"slice"/"fiber" , "x","y","z") , x,y,z apply only to slice and fibers of a tensor or matrix - custom_TD_OP ::Tuple{Union{Array{Any},Array{Float64,2},Array{Float32,2}},Bool} #custom matrix for subspace, and boolean to indicate orthogonality + custom_TD_OP ::Tuple{Union{Array{Any},Array{Float64,2},Array{Float32,2}, joAbstractLinearOperator},Bool} #custom matrix for subspace, and boolean to indicate orthogonality end end # module diff --git a/src/compute_relative_feasibility.jl b/src/compute_relative_feasibility.jl index 60d3ce9..1d62de2 100644 --- a/src/compute_relative_feasibility.jl +++ b/src/compute_relative_feasibility.jl @@ -4,11 +4,9 @@ export compute_relative_feasibility Compute transform-domain relative feasibility w.r.t. a set r_feas = ||P(A*x)-A*x||/||A*x|| """ -function compute_relative_feasibility(m::Vector{TF}, - feasibility_initial::Vector{TF}, - TD_OP::Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joLinearFunction{TF,TF}}}, - P_sub - ) where {TF<:Real,TI<:Integer} +function compute_relative_feasibility(m::Vector{TF}, feasibility_initial::Vector{TF}, + TD_OP::Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joAbstractLinearOperator{TF,TF}}}, + P_sub) where {TF<:Real,TI<:Integer} feasibility_initial[1]=norm(P_sub[1](TD_OP[1]*m) .- TD_OP[1]*m) ./ (norm(TD_OP[1]*m)+(100*eps(TF))) #gc() check if we need this in julia >0.7 diff --git a/src/mat2CDS.jl b/src/mat2CDS.jl index 7cbe492..beba1e6 100755 --- a/src/mat2CDS.jl +++ b/src/mat2CDS.jl @@ -28,5 +28,5 @@ function mat2CDS(A::SparseMatrixCSC{TF,TI}) where {TF<:Real,TI<:Integer} end end -return R,offset + return R, offset end diff --git a/src/projectors/project_l1_Duchi!.jl b/src/projectors/project_l1_Duchi!.jl index 9ba677c..787cd5d 100644 --- a/src/projectors/project_l1_Duchi!.jl +++ b/src/projectors/project_l1_Duchi!.jl @@ -1,109 +1,46 @@ export project_l1_Duchi! -function project_l1_Duchi!(v::Union{Vector{TF},Vector{Complex{TF}}}, b::TF) where {TF<:Real} -# % PROJECTONTOL1BALL Projects point onto L1 ball of specified radius. -# % -# % w = ProjectOntoL1Ball(v, b) returns the vector w which is the solution -# % to the following constrained minimization problem: -# % -# % min ||w - v||_2 -# % s.t. ||w||_1 <= b. -# % -# % That is, performs Euclidean projection of v to the 1-norm ball of radius -# % b. -# % -# % Author: John Duchi (jduchi@cs.berkeley.edu) -# Translated (with some modification) to Julia 1.1 by Bas Peters -if (b <= TF(0)) - error("Radius of L1 ball is negative"); -end -if norm(v, 1) <= b - return v -end -lv=length(v) -u = Vector{TF}(undef,lv) -sv= Vector{TF}(undef,lv) -#u = sort(abs(v),'descend'); -#u = sort(abs.(v), rev=true) #faster than the line below -#u = sort(v, by=abs , rev=true) +""" + project_l1_Duchi!(v::Union{Vector{TF},Vector{Complex{TF}}}, b::TF) + +Projects ector onto L1 ball of specified radius. -#use RadixSort for Float32 (short keywords) -u=copy(v) -if TF==Float32 - u = sort!(abs.(u), rev=true,alg=RadixSort) -else - u = sort!(abs.(u), rev=true,alg=QuickSort) -end +w = ProjectOntoL1Ball(v, b) returns the vector w which is the solution + to the following constrained minimization problem: -#sv = cumsum(u); -sv = cumsum(u) + min ||w - v||_2 + s.t. ||w||_1 <= b. -#rho = find(u > (sv - b) ./ (1:length(u))', 1, 'last'); -#tmp = (u .> ((sv.-b)./ LinSpace(1,lv,lv)))#::BitVector + That is, performs Euclidean projection of v to the 1-norm ball of radius + b. + +Author: John Duchi (jduchi@cs.berkeley.edu) +Translated (with some modification) to Julia 1.1 by Bas Peters +""" +function project_l1_Duchi!(v::Union{Vector{TF},Vector{Complex{TF}}}, b::TF) where {TF<:Real} + b <= TF(0) && error("Radius of L1 ball is negative") + norm(v, 1) <= b && return v -rho = max(1,min(lv,findlast(u .> ((sv.-b)./ (1.0:1.0:lv))))) -#convert(TF,rho) why was this here... -#rho = findlast(tmp)::Int64 + lv = length(v) + u = Vector{TF}(undef,lv) + sv = Vector{TF}(undef,lv) -#theta = max(0, (sv(rho) - b) / rho); -theta = max.(TF(0) , (sv[rho] .- b) ./ rho)::TF + #use RadixSort for Float32 (short keywords) + copyto!(u, v) + if TF==Float32 + u = sort!(abs.(u), rev=true,alg=RadixSort) + else + u = sort!(abs.(u), rev=true,alg=QuickSort) + end -#w = sign(v) .* max(abs(v) - theta, 0); -v .= sign.(v) .* max.(abs.(v) .- theta, TF(0)) + cumsum!(sv, u) -return v -end #end project_l1_Duchi -# -# function project_l1_Duchi!{TF<:Real}(v::Union{Array{TF,3},Array{Complex{TF},3}}, b::TF,mode::String) -# -# if (b < TF(0)) -# error("Radius of L1 ball is negative"); -# end -# (n1,n2,n3)=size(v) -# -# ##Slice based projection -# if mode == "x_slice" -# u = Vector{TF}(n2*n3) -# sv= Vector{TF}(n2*n3) -# v = reshape(v,n1,n2*n3) -# Threads.@threads for i=1:n1 -# u = sort(abs.(v[i,:]), rev=true) -# sv = cumsum(u) -# rho = findlast(u .> ((sv.-b)./ (1.0:1.0:(n2*n3)))) -# theta = max.(TF(0) , (sv[rho] .- b) ./ rho)::TF -# @inbounds v[i,:] .= sign.(v[i,:]) .* max.(abs.(v[:,i]).-theta, TF(0)) -# end -# v = reshape(v,n1,n2,n3) -# elseif mode == "z_slice" -# u = Vector{TF}(n1*n2) -# sv= Vector{TF}(n1*n2) -# v = reshape(v,n1*n2,n3) -# Threads.@threads for i=1:n3 -# u = sort(abs.(v[:,i]), rev=true) -# sv = cumsum(u) -# rho = findlast(u .> ((sv.-b)./ (1.0:1.0:(n1*n2)))) -# theta = max.(TF(0) , (sv[rho] .- b) ./ rho)::TF -# @inbounds v[:,i] .= sign.(v[:,i]) .* max.(abs.(v[:,i]).-theta, TF(0)) -# end -# v = reshape(v,n1,n2,n3) -# elseif mode == "y_slice" -# permutedims(v,[2;1;3]); -# v = reshape(v,n2,n1*n3) -# u = Vector{TF}(n1*n3) -# sv= Vector{TF}(n1*n3) -# Threads.@threads for i=1:size(v,1) -# u = sort(abs.(v[i,:]), rev=true) -# sv = cumsum(u) -# rho = findlast(u .> ((sv.-b)./ (1.0:1.0:(n1*n3)))) -# theta = max.(TF(0) , (sv[rho] .- b) ./ rho)::TF -# @inbounds v[i,:] .= sign.(v[i,:]) .* max.(abs.(v[:,i]).-theta, TF(0)) -# end -# v = reshape(v,n2,n1,n3); -# permutedims(v,[2;1;3]); -# else -# error("provided incorred mode for l1_projection") -# end #end if -# v=vec(v) -# -# return v -# end #end project_l1_Duchi + # Thresholding level + rho = max(1, min(lv, findlast(u .> ((sv.-b)./ (1.0:1.0:lv))))) + theta = max.(TF(0) , (sv[rho] .- b) ./ rho)::TF + + # Projection as soft thresholding + v .= sign.(v) .* max.(abs.(v) .- theta, TF(0)) + + return v +end \ No newline at end of file diff --git a/src/rhs_compose.jl b/src/rhs_compose.jl index cce7283..3f0677c 100644 --- a/src/rhs_compose.jl +++ b/src/rhs_compose.jl @@ -8,7 +8,7 @@ function rhs_compose( l ::Union{ Vector{Vector{TF}},DistributedArrays.DArray{Array{TF,1},1,Array{Array{TF,1},1}} }, y ::Union{ Vector{Vector{TF}},DistributedArrays.DArray{Array{TF,1},1,Array{Array{TF,1},1}} }, rho ::Union{ Vector{TF}, DistributedArrays.DArray{TF,1,Array{TF,1}} }, - TD_OP ::Union{Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joLinearFunction{TF,TF}}},DistributedArrays.DArray{Union{JOLI.joLinearFunction{TF,TF}, SparseMatrixCSC{TF,TI}},1,Array{Union{JOLI.joLinearFunction{TF,TF}, SparseMatrixCSC{TF,TI}},1}} }, + TD_OP ::Union{Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joAbstractLinearOperator{TF,TF}}},DistributedArrays.DArray{Union{JOLI.joAbstractLinearOperator{TF,TF}, SparseMatrixCSC{TF,TI}},1,Array{Union{JOLI.joAbstractLinearOperator{TF,TF}, SparseMatrixCSC{TF,TI}},1}} }, p ::Integer, Blas_active ::Bool, parallel ::Bool diff --git a/src/setup_constraints.jl b/src/setup_constraints.jl index 1f034fe..091b51a 100755 --- a/src/setup_constraints.jl +++ b/src/setup_constraints.jl @@ -44,7 +44,7 @@ end #allocate P_sub = Vector{Any}(undef,nr_constraints) -TD_OP = Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joLinearFunction{TF,TF}}}(undef,nr_constraints) +TD_OP = Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joAbstractLinearOperator{TF,TF}}}(undef,nr_constraints) AtA = Vector{SparseMatrixCSC{TF,TI}}(undef,nr_constraints) #initialize storage for set properties diff --git a/src/setup_multi_level_PARSDMM.jl b/src/setup_multi_level_PARSDMM.jl index be87221..a86614d 100644 --- a/src/setup_multi_level_PARSDMM.jl +++ b/src/setup_multi_level_PARSDMM.jl @@ -34,10 +34,10 @@ end #quantities for PARSDMM for each level m_levels = Vector{Vector{TF}}(undef,n_levels) if options.parallel==false - TD_OP_levels = Vector{Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joLinearFunction{TF,TF}}}}(undef,n_levels) + TD_OP_levels = Vector{Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joAbstractLinearOperator{TF,TF}}}}(undef,n_levels) else #TD_OP_levels = Vector{DistributedArrays.DArray{SparseMatrixCSC{TF,TI},1,Array{SparseMatrixCSC{TF,TI},1}}}(n_levels) - TD_OP_levels = Vector{DistributedArrays.DArray{Union{JOLI.joLinearFunction{TF,TF}, SparseMatrixCSC{TF,TI}},1,Array{Union{JOLI.joLinearFunction{TF,TF}, SparseMatrixCSC{TF,TI}},1}}}(undef,n_levels) + TD_OP_levels = Vector{DistributedArrays.DArray{Union{JOLI.joAbstractLinearOperator{TF,TF}, SparseMatrixCSC{TF,TI}},1,Array{Union{JOLI.joAbstractLinearOperator{TF,TF}, SparseMatrixCSC{TF,TI}},1}}}(undef,n_levels) end AtA_levels = Vector{Vector{SparseMatrixCSC{TF,TI}}}(undef,n_levels) P_sub_levels = Vector{Vector{Any}}(undef,n_levels) diff --git a/src/update_y_l.jl b/src/update_y_l.jl index c9e93ed..6cffe29 100644 --- a/src/update_y_l.jl +++ b/src/update_y_l.jl @@ -15,7 +15,7 @@ function update_y_l( rho ::Vector{TF}, gamma ::Vector{TF}, prox ::Vector{Any}, - TD_OP ::Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joLinearFunction{TF,TF}}}, + TD_OP ::Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joAbstractLinearOperator{TF,TF}}}, log_PARSDMM, P_sub, counter ::Integer, diff --git a/src/update_y_l_parallel.jl b/src/update_y_l_parallel.jl index fb6ac32..5002545 100644 --- a/src/update_y_l_parallel.jl +++ b/src/update_y_l_parallel.jl @@ -14,7 +14,7 @@ function update_y_l_parallel( rho ::Vector{TF}, gamma ::Vector{TF}, prox ::Vector{Any}, - TD_OP ::Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joLinearFunction{TF,TF}}}, + TD_OP ::Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joAbstractLinearOperator{TF,TF}}}, P_sub, x_hat ::Vector{Vector{TF}}, r_pri ::Vector{Vector{TF}}, From 55b12bec979266c41236ceb5bbaae82b542bda89 Mon Sep 17 00:00:00 2001 From: mloubout Date: Fri, 6 Nov 2020 15:10:22 -0500 Subject: [PATCH 2/3] prox_l1 --- src/get_projector.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/get_projector.jl b/src/get_projector.jl index f1fad4f..d01e3a4 100644 --- a/src/get_projector.jl +++ b/src/get_projector.jl @@ -18,6 +18,14 @@ function get_projector(constraint,comp_grid,special_operator_list::Array{String, end end + if constraint.set_type == "prox_l1" + if constraint.TD_OP in special_operator_list + P = x -> copyto!(x,A'*prox_l1!(A*x,constraint.max)) + else + P = x -> prox_l1!(x,constraint.max) + end + end + if constraint.set_type == "l1" if constraint.TD_OP in special_operator_list P = x -> copyto!(x,A'*project_l1_Duchi!(A*x,constraint.max)) From d3a5fa67b60ed31970f860af9d3d295fd8736729 Mon Sep 17 00:00:00 2001 From: Mathias Louboutin Date: Mon, 9 Nov 2020 20:58:24 -0500 Subject: [PATCH 3/3] fix custom TD_OP --- src/PARSDMM.jl | 7 +++--- src/compute_relative_feasibility.jl | 3 ++- src/get_discrete_Grad.jl | 34 +++++------------------------ src/get_projector.jl | 2 +- src/projectors/project_l1_Duchi!.jl | 10 ++++----- src/prox_l2s!.jl | 7 +----- src/setup_constraints.jl | 8 +++---- src/update_y_l.jl | 2 +- src/update_y_l_parallel.jl | 2 +- 9 files changed, 24 insertions(+), 51 deletions(-) diff --git a/src/PARSDMM.jl b/src/PARSDMM.jl index 2ed9ff3..a7f7a9a 100755 --- a/src/PARSDMM.jl +++ b/src/PARSDMM.jl @@ -48,8 +48,8 @@ if isreal(m)==false || isreal(x)==false || isreal(l)==false || isreal(y)==false end # initialize -pp=length(TD_OP); -if feasibility_only==false; pp=pp-1; end; +pp = length(TD_OP); +feasibility_only == false && (pp=pp-1) (ind_ref,N,TD_OP,AtA,p,rho_update_frequency,adjust_gamma,adjust_rho,adjust_feasibility_rho,gamma_ini,rho,gamma,y,y_0,y_old,l,l_0,l_old, l_hat_0,x_0,x_old,r_dual,rhs,s,s_0,Q,prox,log_PARSDMM,l_hat,x_hat,r_pri,d_l_hat,d_H_hat,d_l, @@ -88,10 +88,9 @@ counter=2 x_solve_tol_ref=TF(1.0) #scalar required to determine tolerance for x-minimization, initial value doesn't matter -log_PARSDMM.T_ini=0.0#toq();; +log_PARSDMM.T_ini=0.0 for i=1:maxit #main loop - #form right hand side for x-minimization #tic(); rhs = rhs_compose(rhs,l,y,rho,TD_OP,p,Blas_active,parallel) diff --git a/src/compute_relative_feasibility.jl b/src/compute_relative_feasibility.jl index 1d62de2..1fc1db7 100644 --- a/src/compute_relative_feasibility.jl +++ b/src/compute_relative_feasibility.jl @@ -8,6 +8,7 @@ function compute_relative_feasibility(m::Vector{TF}, feasibility_initial::Vector TD_OP::Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joAbstractLinearOperator{TF,TF}}}, P_sub) where {TF<:Real,TI<:Integer} -feasibility_initial[1]=norm(P_sub[1](TD_OP[1]*m) .- TD_OP[1]*m) ./ (norm(TD_OP[1]*m)+(100*eps(TF))) + pm = TD_OP[1]*m + feasibility_initial[1]=norm(P_sub[1](pm) .- pm) ./ (norm(pm)+(100*eps(TF))) #gc() check if we need this in julia >0.7 end diff --git a/src/get_discrete_Grad.jl b/src/get_discrete_Grad.jl index 70e7a16..cdfec76 100644 --- a/src/get_discrete_Grad.jl +++ b/src/get_discrete_Grad.jl @@ -14,13 +14,7 @@ export get_discrete_Grad - TD_OP : transform domain operator as a sparse matrix """ function get_discrete_Grad(n1,n2,h1::TF,h2::TF,TD_type::String) where {TF<:Real} - -# if TF==Float64 -# TI=Int64 -# else -# TI=Int32 -# end -TI=Int64 + TI=Int64 #define difference matrix D acting on vectorized model using Kronecker products Ix = SparseMatrixCSC{TF}(LinearAlgebra.I,n1,n1) #x @@ -28,10 +22,6 @@ TI=Int64 Dx = spdiagm(0 => ones(TF,n1-1)*-1, 1 => ones(TF,n1-1)*1); Dx = Dx[1:end-1,:] ./ h1 Dz = spdiagm(0 => ones(TF,n2-1)*-1, 1 => ones(TF,n2-1)*1); Dz = Dz[1:end-1,:] ./ h2 - # Ix = convert(SparseMatrixCSC{TF,TI},Ix) - # Iz = convert(SparseMatrixCSC{TF,TI},Iz) - # Dx = convert(SparseMatrixCSC{TF,TI},Dx) - # Dz = convert(SparseMatrixCSC{TF,TI},Dz) if TD_type=="D_z" D_OP = kron(Dz,Ix) #D2z @@ -42,10 +32,11 @@ TI=Int64 D2x = kron(Iz,Dx) D_OP = vcat(D2z,D2x) #D2D end - -return D_OP + + return D_OP end + """ 3D version input: n1 : number of grid points in the first dimension @@ -58,13 +49,7 @@ end output : TD_OP : transform domain operator as a sparse matrix """ function get_discrete_Grad(n1,n2,n3,h1::TF,h2::TF,h3::TF,TD_type::String) where {TF<:Real} - -# if TF==Float64 -# TI=Int64 -# else -# TI=Int32 -# end -TI = Int64 + TI = Int64 #define difference matrix D acting on vectorized model using Kronecker products Ix = SparseMatrixCSC{TF}(LinearAlgebra.I,n1,n1) #x @@ -74,13 +59,6 @@ TI = Int64 Dy = spdiagm(0 => ones(TF,n2-1)*-1, 1 => ones(TF,n2-1)*1); Dy = Dy[1:end-1,:] ./ h2 Dz = spdiagm(0 => ones(TF,n3-1)*-1, 1 => ones(TF,n3-1)*1); Dz = Dz[1:end-1,:] ./ h3 - # Ix = convert(SparseMatrixCSC{TF,TI},Ix) - # Iy = convert(SparseMatrixCSC{TF,TI},Iy) - # Iz = convert(SparseMatrixCSC{TF,TI},Iz) - # Dx = convert(SparseMatrixCSC{TF,TI},Dx) - # Dy = convert(SparseMatrixCSC{TF,TI},Dy) - # Dz = convert(SparseMatrixCSC{TF,TI},Dz) - if TD_type=="D_z" D_OP = kron(Dz,Iy,Ix) elseif TD_type=="D_y" @@ -94,5 +72,5 @@ TI = Int64 D_OP = vcat(D3z,D3y,D3x) end -return D_OP + return D_OP end diff --git a/src/get_projector.jl b/src/get_projector.jl index d01e3a4..0bade18 100644 --- a/src/get_projector.jl +++ b/src/get_projector.jl @@ -20,7 +20,7 @@ function get_projector(constraint,comp_grid,special_operator_list::Array{String, if constraint.set_type == "prox_l1" if constraint.TD_OP in special_operator_list - P = x -> copyto!(x,A'*prox_l1!(A*x,constraint.max)) + P = x -> copyto!(x,.5f0*A'*prox_l1!(A*x,constraint.max)) else P = x -> prox_l1!(x,constraint.max) end diff --git a/src/projectors/project_l1_Duchi!.jl b/src/projectors/project_l1_Duchi!.jl index 787cd5d..690de78 100644 --- a/src/projectors/project_l1_Duchi!.jl +++ b/src/projectors/project_l1_Duchi!.jl @@ -22,15 +22,15 @@ function project_l1_Duchi!(v::Union{Vector{TF},Vector{Complex{TF}}}, b::TF) wher norm(v, 1) <= b && return v lv = length(v) - u = Vector{TF}(undef,lv) - sv = Vector{TF}(undef,lv) + u = similar(v) + sv = Vector{TF}(undef, lv) #use RadixSort for Float32 (short keywords) copyto!(u, v) if TF==Float32 - u = sort!(abs.(u), rev=true,alg=RadixSort) + u = sort!(abs.(u), rev=true, alg=RadixSort) else - u = sort!(abs.(u), rev=true,alg=QuickSort) + u = sort!(abs.(u), rev=true, alg=QuickSort) end cumsum!(sv, u) @@ -43,4 +43,4 @@ function project_l1_Duchi!(v::Union{Vector{TF},Vector{Complex{TF}}}, b::TF) wher v .= sign.(v) .* max.(abs.(v) .- theta, TF(0)) return v -end \ No newline at end of file +end diff --git a/src/prox_l2s!.jl b/src/prox_l2s!.jl index 8f5675d..77d0f31 100644 --- a/src/prox_l2s!.jl +++ b/src/prox_l2s!.jl @@ -1,11 +1,6 @@ export prox_l2s! -function prox_l2s!( - x ::Vector{TF}, - rho ::TF, - m ::Vector{TF} - ) where {TF<:Real} - +function prox_l2s!(x::Vector{TF}, rho::TF, m::Vector{TF}) where {TF<:Real} x .= (x .* rho .+ m) ./ (rho .+ 1.0); return x end diff --git a/src/setup_constraints.jl b/src/setup_constraints.jl index 091b51a..32cabff 100755 --- a/src/setup_constraints.jl +++ b/src/setup_constraints.jl @@ -44,14 +44,14 @@ end #allocate P_sub = Vector{Any}(undef,nr_constraints) -TD_OP = Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joAbstractLinearOperator{TF,TF}}}(undef,nr_constraints) +TD_OP = Vector{Union{Array{TF, 2}, SparseMatrixCSC{TF,TI},JOLI.joAbstractLinearOperator{TF,TF}}}(undef,nr_constraints) AtA = Vector{SparseMatrixCSC{TF,TI}}(undef,nr_constraints) #initialize storage for set properties set_Prop=set_properties(zeros(nr_constraints),zeros(nr_constraints),zeros(nr_constraints),Vector{Tuple{TI,TI}}(undef,nr_constraints),Vector{Tuple{String,String,String,String}}(undef,nr_constraints),zeros(nr_constraints),Vector{Vector{TI}}(undef,nr_constraints)) #other initialization -special_operator_list = ["DFT", "DCT", "wavelet"] #complex valued operators that are orthogonal +special_operator_list = ["DFT", "DCT", "wavelet", "curvelet"] #complex valued operators that are orthogonal N = prod(comp_grid.n) for i = 1:nr_constraints @@ -66,13 +66,13 @@ for i = 1:nr_constraints error("l1 and l2 constraints only available for matrix or tensor mode, currently") end - (A,AtA_diag,dense,TD_n,banded) = get_TD_operator(comp_grid,constraint[i].TD_OP,TF) + ~(constraint[i].custom_TD_OP[1] == []) && (A = constraint[i].custom_TD_OP[1]) P_sub[i] = get_projector(constraint[i],comp_grid,special_operator_list,A,TD_n,TF) if constraint[i].TD_OP in special_operator_list - TD_OP[i] = SparseMatrixCSC{TF}(LinearAlgebra.I,N,N) + TD_OP[i] = SparseMatrixCSC{TF}(LinearAlgebra.I,N,N) else TD_OP[i] = A end diff --git a/src/update_y_l.jl b/src/update_y_l.jl index 6cffe29..86636c4 100644 --- a/src/update_y_l.jl +++ b/src/update_y_l.jl @@ -15,7 +15,7 @@ function update_y_l( rho ::Vector{TF}, gamma ::Vector{TF}, prox ::Vector{Any}, - TD_OP ::Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joAbstractLinearOperator{TF,TF}}}, + TD_OP ::Vector{<:Union{SparseMatrixCSC{TF,TI}, joAbstractLinearOperator{TF,TF}}}, log_PARSDMM, P_sub, counter ::Integer, diff --git a/src/update_y_l_parallel.jl b/src/update_y_l_parallel.jl index 5002545..6ead3a9 100644 --- a/src/update_y_l_parallel.jl +++ b/src/update_y_l_parallel.jl @@ -14,7 +14,7 @@ function update_y_l_parallel( rho ::Vector{TF}, gamma ::Vector{TF}, prox ::Vector{Any}, - TD_OP ::Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joAbstractLinearOperator{TF,TF}}}, + TD_OP ::Vector{<:Union{SparseMatrixCSC{TF,TI},joAbstractLinearOperator{TF,TF}}}, P_sub, x_hat ::Vector{Vector{TF}}, r_pri ::Vector{Vector{TF}},