slimgroup · PetersBas · Jan 23, 2022 · Nov 5, 2020 · Nov 6, 2020 · Nov 10, 2020
diff --git 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,
@@ -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/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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -4,12 +4,11 @@ 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)))
+    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
@@ -14,24 +14,14 @@ 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
     Iz = SparseMatrixCSC{TF}(LinearAlgebra.I,n2,n2) #z
     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
@@ -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,.5f0*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))

diff --git 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
-Original file line number
+Diff line change
@@ Expand Up @@
         end
       end
-    return R,offset
+      return R, offset
     end