@PetersBas curvelets and any generic JOLI operator should work with that, I also added custom_operator to work with any constraints but not sure if works directly for everything.

Any news on this? Let me know if anything I can help with to make it happen @mloubout

Well it's kinda working but never really needed any of this so left it as is. Don't think Bad is maintaining it anymore neither

I plan on continuing to maintain this, will update to newer julia version sometime soon. What specific functionality do you need @ziyiyin97 ? Then I'll will test & incorporate the changes in this pull request.

I am looking into projection to norm(A*x,1)<sigma where A is a JOLI operator. Thanks @PetersBas and let me know if anything I can help with

Glad to hear, sorry I spoke for yourself.

@ziyiyin97 which JOLI operatorr(s) do you need? It's best if I include them myself, so I can make them a pre-set choice that you can easily select. That's much better than using the (limited) support for your own custom operators. The software does not treat all operators equally: for computational and memory efficiency, the algorithm needs to know which operators are sparse, dense, map to complex numbers, which operators take a long time to compute if formed on the fly instead of explicitly precomputing the matrix.

I can have a look at the JOLI operators that you need, so after including them with the right properties attached, everyone can use those efficiently.

I currently have this thing in mind: I have an array of images x1, x2, ..., xn and I want to do TV on its time-derivative (if you treat n as number of time steps), or do TV on x2-x1, x3-x1, ..., xn-x1. So basically the operator would be joKron(TV, joEye(n-1)) * D where TV is the total variance on an image, D is either the time-derivative, or to compute the differences from each time step of x to the first time step. It should also be extended to array of 3D cubes as well.

Although I suspect you could work out a very efficient algorithm with this specific structure of operator, I will be very grateful if we include a general code that works for any JOLI operator (i.e. we assume to only have forward and adjoint of it). One PR per JOLI operator seems a lot of work. @PetersBas

These constraints can be implemented with the current code, no additional functionality needs to be implemented. I'll add an example.

OK thanks @PetersBas for the quick reply! I am sorry if I miss some details in your software. Me carelessly claiming it's not working while it actually works doesn't sound good : )

There might be some other complicated constraints, including:

1. subsampling + curvelet + l1, i.e. in the previous array of images example, impose l1 on each curvelet transformed images, while they contain different sigma
2. there is mirror-extended curvelet which is now under gatech github (hasn't be fully released yet), equipped with forward and adjoint. That's still one of the fundamental reasons that I look for an implementation where the projection can generally work with an "assumed unknown" JOLI operator @mloubout correct me if I am wrong

Concerning the second point above I think a generic JOLI input would be quite good. In theory, JOLI allows you to build any linear operator for any need that one may wanna use for constraints. So while it may incur some performance hit since not natively implemented in the constraints. being able to play with custom linear constraints easily would be a great addition IMO. And as @ziyiyin97 said, in some case the JOLI operator may be under a license that doesn't permit it to be added here.

Yes, I agree the additions in this pr would be helpful. Will merge.

Thanks for the update but is this PR working now? I got

#bounds:
m_min     = 1500.0
m_max     = 4500.0
set_type  = "l1"
TD_OP     = "identity"
app_mode  = ("matrix","")
custom_TD_OP = (joMatrix(randn(options.FL, comp_grid.n)),false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

options.parallel       = false
(P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL)
(TD_OP,AtA,l,y)        = PARSDMM_precompute_distribute(TD_OP,set_Prop,comp_grid,options)

and

julia> (TD_OP,AtA,l,y)        = PARSDMM_precompute_distribute(TD_OP,set_Prop,comp_grid,options)
ERROR: MethodError: no method matching PARSDMM_precompute_distribute(::Vector{Union{joAbstractLinearOperator{Float32, Float32}, SparseArrays.SparseMatrixCSC{Float32, Int64}, Matrix{Float32}}}, ::set_properties, ::compgrid, ::PARSDMM_options)
Closest candidates are:
  PARSDMM_precompute_distribute(::Array{Union{joAbstractLinearOperator{TF, TF}, SparseArrays.SparseMatrixCSC{TF, TI}}, 1}, ::Any, ::Any, ::Any) where {TF<:Real, TI<:Integer} at /Users/francisyin/.julia/dev/SetIntersectionProjection/src/PARSDMM_precompute_distribute.jl:6
Stacktrace:
 [1] top-level scope
   @ REPL[39]:1
 [2] top-level scope
   @ ~/.julia/packages/CUDA/0IDh2/src/initialization.jl:52

Maybe a new type should be defined, like

const OpType{TF, TI}=Union{Matrix{TF, TI}, SparseArrays.SparseMatrixCSC{TF, TI}, joAbstractOperator{TF, TI}}

and some work is needed for some operations in PARSDMM, e.g. skipping mat2CDS. Am I right?

made many more changes, JOLI operators can now be used, see SetIntersectionProjection.jl/examples/ConstraintSetupExamples.jl for examples.