switch back to SparseArrays to avoid GC corruption

src/Convex.jl

@@ -130,11 +130,26 @@ function vec_tril(M)
     return M[inds]
 end
 
-const SPARSE_VECTOR{T} = GBVector{T,T}
-const SPARSE_MATRIX{T} = GBMatrix{T,T}
-spzeros(T, d) = GBVector{T,T}(d)
-spzeros(T, n, m) = GBMatrix{T,T}(n, m)
-spidentity(T, d) = GBMatrix{T,T}(Diagonal(ones(T, d)))
+# const SPARSE_VECTOR{T} = GBVector{T,T}
+# const SPARSE_MATRIX{T} = GBMatrix{T,T}
+# spzeros(T, d) = GBVector{T,T}(d)
+# spzeros(T, n, m) = GBMatrix{T,T}(n, m)
+# spidentity(T, d) = GBMatrix{T,T}(Diagonal(ones(T, d)))
+# create_sparse(T, args...) = GBMatrix{T,T}(args...)
+
+const SPARSE_VECTOR{T} = Vector{T}
+const SPARSE_MATRIX{T} = SparseMatrixCSC{T,Int}
+spzeros(T, d) = zeros(T, d)
+spzeros(T, n, m) = SparseArrays.spzeros(T, n, m)
+spidentity(T, d) = sparse(one(T)*I, d, d)
+# function create_sparse(::Type{T}, I, J, V, args...) where T
+#     return SparseArrays.sparse(I, J, T.(V), args...)::SPARSE_MATRIX{T}
+# end
+function create_sparse(::Type{T}, args...) where {T}
+    local result::SPARSE_MATRIX{T}
+    result = SparseArrays.sparse(args...)
+    return result
+end
 
 include("Context.jl")
 ### modeling framework

src/SparseTape.jl

@@ -16,7 +16,7 @@ function SparseAffineOperation(
     A::AbstractMatrix{T},
     b::AbstractVector{T},
 ) where {T}
-    return SparseAffineOperation{T}(SPARSE_MATRIX{T}(A), GBVector{T, T}(b))
+    return SparseAffineOperation{T}(create_sparse(T,A), SPARSE_VECTOR{T}(b))
 end
 
 
@@ -30,11 +30,25 @@ end
 # function SparseAffineOperation(A::AbstractSparseMatrix, b)
 #     return SparseAffineOperation(SparseMatrixCSC(A), b)
 # end
-# SparseAffineOperation(A, b) = SparseAffineOperation(sparse(A), b)
+# SparseAffineOperation(A, b) = SparseAffineOperation(create_sparse(A), b)
 
 mutable struct SparseTape{T}
     operations::Vector{SparseAffineOperation{T}}
     variables::Vector{MOI.VariableIndex}
+    function SparseTape{T}(operations::Vector{SparseAffineOperation{T}}, variables::Vector{MOI.VariableIndex}) where {T}
+        # Is this necessary?
+        # if !issorted(variables; by = x->x.value)
+        #     p = sortperm(variables; by = x->x.value)
+        #     op = foldl(compose, operations)
+        #     matrix = op.matrix[:, p]
+        #     vector = op.vector
+        #     operations = [SparseAffineOperation(matrix, vector)]
+        #     variables = variables[p]
+        # end
+        new(operations, variables)
+    end
+
+    SparseTape(operations::Vector{SparseAffineOperation{T}}, variables::Vector{MOI.VariableIndex}) where {T} = SparseTape{T}(operations, variables)
 end
 
 MOI.output_dimension(v::SparseTape) = size(v.operations[1].matrix, 1)
@@ -66,8 +80,8 @@ Base.real(tape::SparseTape) = tape
 function Base.imag(c::SparseTape{T}) where {T}
     n = MOI.output_dimension(c)
     m = length(c.variables)
-    mat = GBMatrix{T, T}(n,m)
-    v = zeros(T, n)
+    mat = spzeros(T, n, m)
+    v = spzeros(T, n)
     op = SparseAffineOperation(mat, v)
 
     # Hack re-use variables from input

src/atoms/affine/diag.jl

@@ -78,7 +78,7 @@ function _conic_form!(context::Context{T}, x::DiagAtom) where {T}
         sz_diag = Base.min(num_rows + k, num_cols)
     end
 
-    select_diag = SPARSE_MATRIX{T}(sz_diag, length(x.children[1]))
+    select_diag = spzeros(T, sz_diag, length(x.children[1]))
     for i in 1:sz_diag
         select_diag[i, start_index] = 1
         start_index += num_rows + 1

src/atoms/affine/diagm.jl

@@ -66,8 +66,8 @@ function _conic_form!(context::Context{T}, x::DiagMatrixAtom) where {T}
     I = collect(1:sz+1:sz*sz)
     J = collect(1:sz)
     V = one(T)
-    coeff = GBMatrix{T, T}(I, J, V, sz * sz, sz)
-    # coeff = sparse(, 1:sz, one(T),
+    coeff = create_sparse(T, I, J, V, sz * sz, sz)
+    # coeff = create_sparse(, 1:sz, one(T),
 
     return operate(add_operation, T, sign(x), coeff, obj)
 end

src/atoms/affine/index.jl

@@ -60,9 +60,8 @@ function _conic_form!(context::Context{T}, x::IndexAtom) where {T}
 
     if x.inds === nothing
         sz = length(x.cols) * length(x.rows)
-        J = Array{Int}(undef, sz)
+        J = Vector{Int}(undef, sz)
         k = 1
-
         num_rows = x.children[1].size[1]
         for c in x.cols
             for r in x.rows
@@ -71,9 +70,9 @@ function _conic_form!(context::Context{T}, x::IndexAtom) where {T}
             end
         end
 
-        index_matrix = SPARSE_MATRIX{T}(collect(1:sz), J, one(T), m, n)
+        index_matrix = create_sparse(T, collect(1:sz), J, one(T), m, n)
     else
-        index_matrix = SPARSE_MATRIX{T}(collect(1:length(x.inds)), collect(x.inds), one(T), m, n)
+        index_matrix = create_sparse(T, collect(1:length(x.inds)), collect(x.inds), one(T), m, n)
     end
 
     return operate(add_operation, T, sign(x), index_matrix, obj)

src/atoms/affine/multiply_divide.jl

@@ -66,7 +66,7 @@ function real_convert(::Type{T}, x::Number) where {T}
     return T(x)
 end
 function real_convert(::Type{T}, x::AbstractMatrix) where {T}
-    return GBMatrix{T,T}(x)
+    return create_sparse(T, x)
 end
 
 function real_convert(::Type{T}, x::SPARSE_MATRIX{T}) where {T}
@@ -132,7 +132,7 @@ function _conic_form!(context::Context{T}, x::MultiplyAtom) where {T}
             add_operation,
             T,
             sign(x),
-            kron(_id(T, x.size[2]), const_multiplier),
+            kron(spidentity(T, x.size[2]), const_multiplier),
             objective,
         )
 
@@ -151,15 +151,12 @@ function _conic_form!(context::Context{T}, x::MultiplyAtom) where {T}
             add_operation,
             T,
             sign(x),
-            kron(transpose(const_multiplier), _id(T, x.size[1])),
+            kron(transpose(const_multiplier), spidentity(T, x.size[1])),
             objective,
         )
     end
 end
 
-# _id(T, n) = Diagonal(one(T)*I, n)
-_id(T, n) = spidentity(T, n)
-
 function *(x::AbstractExpr, y::AbstractExpr)
     if isequal(x, y) && x.size == (1, 1)
         return square(x)

src/complex_operate.jl

@@ -118,11 +118,13 @@ end
 ## Binary
 
 function complex_operate(::typeof(add_operation), ::Type{T}, A, c) where {T}
+    L = real_operate(add_operation, T, real(A), real(c))
+    R = real_operate(add_operation, T, imag(A), imag(c))
     re = real_operate(
         -,
         T,
-        real_operate(add_operation, T, real(A), real(c)),
-        real_operate(add_operation, T, imag(A), imag(c)),
+        L,
+        R,
     )
 
     im = real_operate(
@@ -131,7 +133,7 @@ function complex_operate(::typeof(add_operation), ::Type{T}, A, c) where {T}
         real_operate(add_operation, T, real(A), imag(c)),
         real_operate(add_operation, T, imag(A), real(c)),
     )
-    if re isa Vector && im isa Vector
+    if re isa AbstractVector && im isa AbstractVector
         return ComplexStructOfVec(re, im)
     else
         return ComplexTape(re, im)

src/constant.jl

@@ -90,6 +90,8 @@ function _conic_form!(context::Context, C::ComplexConstant)
     )
 end
 
+constant(x::Constant) = x
+constant(x::ComplexConstant) = x
 function constant(x)
     # Convert to matrix
     x = [x;;]

src/operate.jl

@@ -1,9 +1,11 @@
 complex_promote(tape::ComplexTape) = tape
 complex_promote(tape::SparseTape) = tape
 complex_promote(v::ComplexStructOfVec) = v
+complex_promote(x::Real) = complex(x)
+complex_promote(x::Complex) = x
 
 function complex_promote(v::AbstractVector{T}) where {T}
-    return ComplexStructOfVec(v, SPARSE_VECTOR{T}(length(v)))
+    return ComplexStructOfVec(v, spzeros(T, length(v)))
 end
 
 # Here we run `complex_promote` and dispatch to either `real_operate`
@@ -19,7 +21,7 @@ function operate(op::F, ::Type{T}, sign::Sign, args...) where {F,T}
             x, y = args
             if !(
                 typeof(x) in
-                (T, Complex{T}, SPARSE_MATRIX{T}, GBMatrix{Complex{T},Complex{T}})
+                (T, Complex{T}, SPARSE_MATRIX{T},  SPARSE_MATRIX{Complex{T}}, Matrix{T},Matrix{Complex{T}})
             )
                 error(
                     "Convex.jl internal error: unexpected type $(typeof(x)) for first argument of operation $op of with sign=$sign",
@@ -75,7 +77,7 @@ function operate(op::F, ::Type{T}, sign::Sign, args...) where {F,T}
                 @assert length(args) == 2
 
                 x, y = args
-                if !(typeof(x) in (T, SPARSE_MATRIX{T}))
+                if !(typeof(x) in (T, SPARSE_MATRIX{T}, Matrix{T}))
                     error(
                         "Convex.jl internal error: unexpected type $(typeof(x)) for first argument of operation $op of with sign=$sign",
                     )

src/problem_depot/problems/sdp.jl

@@ -582,7 +582,7 @@ end
 ) where {T,test}
     a = 2 + 4im
     x = ComplexVariable()
-    objective = norm2(a - x)
+    objective = norm2(a + (-x))
     c1 = real(x) >= 0
     p = minimize(objective, c1; numeric_type = T)
 

src/real_operate.jl

@@ -98,7 +98,6 @@ function real_operate(
     tape::SparseTape{T},
 ) where {T<:Real}
     d = MOI.output_dimension(tape)
-
     return add_operation(
         tape,
         SparseAffineOperation(-spidentity(T, d), spzeros(T, d)),
@@ -150,7 +149,7 @@ function real_operate(
     return SparseTape([SparseAffineOperation(A, b)], x)
 end
 
-blockdiag(x, y) = Base.blockdiag(x, y)
+blockdiag(xs...) = SparseArrays.blockdiag(xs...)::SPARSE_MATRIX
 
 function blockdiag(xs::GBMatrix{T,T}...) where {T}
     N = length(xs)
@@ -285,7 +284,7 @@ function real_operate(
     A::AbstractMatrix,
     tape::SparseTape{T},
 ) where {T<:Real}
-    return add_operation(tape, SparseAffineOperation(A, zeros(T, size(A, 1))))
+    return add_operation(tape, SparseAffineOperation(A, spzeros(T, size(A, 1))))
 end
 
 function real_operate(
@@ -306,7 +305,7 @@ function real_operate(
     d = MOI.output_dimension(tape)
     return add_operation(
         tape,
-        SparseAffineOperation(spidentity(T, d), spzeros(T, d)),
+        SparseAffineOperation(x*spidentity(T, d), spzeros(T, d)),
     )
 end
 

src/variable_template.jl

@@ -53,7 +53,7 @@ function to_tape(v::MOI.VectorOfVariables, ::Context{T}) where {T}
     var_inds = v.variables
     d = length(var_inds)
     return SparseTape(
-        [SparseAffineOperation(spidentity(T, d), SPARSE_VECTOR{T}(d))],
+        [SparseAffineOperation(spidentity(T, d), spzeros(T,d))],
         var_inds,
     )
 end