Pg/one big module

[goulart-paul]

Implements the following changes:

• Changes name from QOCS->COSMO everywhere
• Flattens everything into a single module COSMO
• All convex set types are now lightweight structures with no function pointers.
• Implements a SplitVector type that contains an array of 'views' into an array of AbstractConvexSets.
• Implements a container-type convex set of type "CompositeConvexSet" that acts as a wrapper around an array of AbstractConvexSet constraints, and dispatches calls to project! etc in the subviews of a SplitVector.
• Changed many of the base types and function signatures to enforce a common (but parametric) floating point type across solver components. We should be able to a make a 64 bit solver, or 32, but not mix types.
• Various code efficiency changes in different places.
• Sporadic enforcement of style guidelines, particularly or lower case function names.

[nrontsis]

Thanks a lot Paul for all the work.

The issues with Diagonal*CSC matrix multiplication have been resolved as of Julia 1.0.1. Thus, I believe that the specialised function rmul! and lmul! in algebra.jl are unnecessary. Running

using LinearAlgebra, SparseArrays
using BenchmarkTools
include("src/algebra.jl")

n = Int(1e4)
X = sprandn(n, n, 0.001)
d = randn(n)
D = Diagonal(d)


# Our versions
Z = copy(X) # Make a copy just for sanity purposes
@btime lmul!(D, Z)
@btime rmul!(Z, D)

# Julia's 1.0.1 version
@btime Y = D*X
@btime Y = X*D

gives

  1.198 ms (0 allocations: 0 bytes)
  733.798 μs (0 allocations: 0 bytes)
  357.161 μs (7 allocations: 1.60 MiB)
  337.446 μs (7 allocations: 1.60 MiB)

i.e. Julia's built-in functions are faster, even though they allocate memory at every run.

As a matter of fact I think that we are overly concerned about memory allocation. I don't think that we should be, unless there is a timing impact (say >10%) that supports that. For example, why should we store the diagonal scaling matrices when the performance of e.g.:

@btime Y = D*X
@btime Y = Diagonal(d)*X

is identical?

Julia's garbage collector is very efficient in most cases, so I don't see the benefit of aiming at minimising memory allocation at the cost of code readability.

[nrontsis]

I think that there still remain places where Float64 is hardcoded. @goulart-paul would you be happy if I were to try changing these by pushing into your branch?

Also, I would suggest replacing SplitView withAbstractArray/AbstractVector in the various function definitions, which is more general and standard. Plus, it would allow easier calls of any function in independent bechmarking/testing scripts, e.g. a call like the following

using LinearAlgebra, SparseArrays
include("src/COSMO.jl")
include("src/projections.jl")

n=1000
X = randn(n,n)
cone = PsdCone{Float64}(n*n)
project!(X[:], cone) 

which currently throws the error:

ERROR: LoadError: MethodError: no method matching project!(::Array{Float64,1}, ::PsdCone{Float64})
Closest candidates are:
  project!(!Matched::SubArray{T,1,Array{T,1},Tuple{UnitRange{Int64}},true}, ::PsdCone{T}) where T at /Users/nrontsis/Downloads/QOCS.jl/src/convexset.jl:133

[goulart-paul]

I agree that we should consider the 1.0.1 operations, which only appeared recently. I rewrote the algebra.jl code slightly to try to improve things via @inbounds but it didn't really help.

I don't understand why your example is so fast in the second case. If I replace this

@btime Y = D*X
@btime Y = X*D

with this

@btime $X .= $D*$X
@btime $X .= $X*$D

then the times are about the same between our version and the current 1.0.1.

Regarding memory allocation : I think my strategy here is to be careful about it where it doesn't seem to matter or if it makes an improvement, since I have one eye always on an embedded version. If your example is correct though I'm happy to not worry about it someplace like this case.

[goulart-paul]

I think that there still remain places where Float64 is hardcoded. @goulart-paul would you be happy if I were to try changing these by pushing into your branch?
Yes, we should certainly do this. I was just trying not to break too many things in one commit.

The only point where I would hesitate to do this would be for things like reported timings and solve accuracy to the user. Maybe then it makes sense to have a second parametric type defined, or just to use Float64 always.

Also, I would suggest replacing SplitView withAbstractArray/AbstractVector in the various function definitions, which is more general and standard. Plus, it would allow easier calls of any function in independent bechmarking/testing scripts

Yes, I agree with this. It was that way originally but I was getting weird results from --track-allocation=user when doing composite projection, and I thought maybe it was a dispatch problem so tried to make everything very strongly typed to debug it.

I think making all of the individual projection operations (except for CompositeConvexSet) work with an AbstractVector input is fine as long as it doesn't kill the performance, which I don't think it will.

Note however that the SplitVector{T} type should still contain an array of SplitView{T}, not AbstractArray{T}, since I think it is faster to have a single clear type within a container.

[nrontsis]

I guess my benchmark was flawed as running:

@btime lmul!(D, Z)

ends up replacing Z with D^k*Z (where k equal to the number of calls by btime), so Z ends up having huge number & Infs, which makes things slower.

I believe that what you say is correct, i.e. that your algebra.jl functions and Julia's 1.0.1 are similar in performance. However, I would still suggest to rely on Julia's for reasons of simplicity and maintenance.

[nrontsis]

The only point where I would hesitate to do this would be for things like reported timings and solve accuracy to the user. Maybe then it makes sense to have a second parametric type defined, or just to use Float64 always.

Maybe for the final reporting of the accuracy we can convert the matrices to Float64 when we do the final residuals calculation. For the timings, we can always save them in Float64, if necessary.

Maybe then it makes sense to have a second parametric type defined, or just to use Float64 always.

Do you mean not have the option of Float32 in the solver at all? I think using Float32 is very relevant e.g. when we use indirect solvers (CG and Block-Lanczos) internally, as, hopefully, the accuracy of the ADMM operations would be dominated by the limited accuracy of the indirect solvers and not by machine precision.

[goulart-paul]

I believe that what you say is correct, i.e. that your algebra.jl functions and Julia's 1.0.1 are similar in performance. However, I would still suggest to rely on Julia's for reasons of simplicity and maintenance.

Yes, I'm fine with this also. If the performance is about the same and we don't allocate though, then what we should really do is submit a fix to the main julia repo.

[goulart-paul]

Maybe then it makes sense to have a second parametric type defined, or just to use Float64 always.

Do you mean not have the option of Float32 in the solver at all? I think using Float32 is very relevant e.g. when we use indirect solvers (CG and Block-Lanczos) internally, as, hopefully, the accuracy of the ADMM operations would be dominated by the limited accuracy of the indirect solvers and not by machine precision.

No, I meant something different. I have tried to write the code such that it an use any floating point type, but that consistency of type is enforced across the iterates, problem data, residuals, working data etc. Everything is the same type T. This means that in the end you would get a completely FloatXX based solver depending on the inputs.

However, if T=Float16 we may want some components to be some other T2 = Float64 or whatever, specifically for things like timings or other reporting. We could then enforce all the timings to be the same type as each other, but not necessarily the same as the data.

Note that it's not clear to me whether what I have done is the right approach. You could make a case that some of the problem data should be allowed to be some small size Int, or even Bool, for things like sparsity constraints, so long as the KKT factors we form from it ends up as some variety of float.

[nrontsis]

I believe that what you say is correct, i.e. that your algebra.jl functions and Julia's 1.0.1 are similar in performance. However, I would still suggest to rely on Julia's for reasons of simplicity and maintenance.

Yes, I'm fine with this also. If the performance is about the same and we don't allocate though, then what we should really do is submit a fix to the main julia repo.

You can avoid the memory allocation by using mul! instead of *.

[goulart-paul]

I believe that what you say is correct, i.e. that your algebra.jl functions and Julia's 1.0.1 are similar in performance. However, I would still suggest to rely on Julia's for reasons of simplicity and maintenance.

Yes, I'm fine with this also. If the performance is about the same and we don't allocate though, then what we should really do is submit a fix to the main julia repo.

You can avoid the memory allocation by using mul! instead of *.

In that case there is no reason not to use the Julia base version, unless it is really bad with identity matrices or something. I'm not committed to our internal version and am happy to dump it. I can also imagine that having it there might cause namespace conflicts as well.

[migarstka]

Changes approved. I will merge this into devel and create a separate issue about the open lmul! rmul! question.