release v0.4.1 (#119)

* Constraint reduction for inequalities with infinite RHS (#115)

* internal changes for sdp efficiency and rust compatibility

* Constructor float type defaults (fixes #117)

CHANGELOG.md

@@ -6,6 +6,17 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
 
 Version numbering in this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).  We aim to keep the core solver functionality and minor releases in sync between the Rust/Python and Julia implementations.   Small fixes that affect one implementation only may result in the patch release versions differing.
 
+## [0.4.1] - 2023-08-03
+### Changed 
+
+Added optional feature to remove inequality constraints with very large upper bounds.   This feature is enabled by default but can be turned off using the `presolve_enable` setting.  
+
+Bug fix in equilibration for NN and zero cones.
+### Julia specific changes
+
+Internal implementation of composite cone logic updated to more closely match the rust version.
+
+Internal modifications to SDP cone implementation to reduce allocations.
 ## [0.4.0] - 2023-25-02
 
 ### Changed 
@@ -43,8 +54,8 @@ Version numbering in this project adheres to [Semantic Versioning](https://semve
 - Initial release
 
 
-
-[0.4.0]: https://github.com/pyo3/pyo3/compare/v0.4.0...v0.3.0
-[0.3.0]: https://github.com/pyo3/pyo3/compare/v0.3.0...v0.2.0
-[0.2.0]: https://github.com/pyo3/pyo3/compare/v0.2.0...v0.1.0
-[0.1.0]: https://github.com/PyO3/pyo3/tree/0.1.0
+[0.4.1]: https://github.com/oxfordcontrol/Clarabel.jl/compare/v0.4.0...v0.4.1
+[0.4.0]: https://github.com/oxfordcontrol/Clarabel.jl/compare/v0.4.0...v0.3.0
+[0.3.0]: https://github.com/oxfordcontrol/Clarabel.jl/compare/v0.3.0...v0.2.0
+[0.2.0]: https://github.com/oxfordcontrol/Clarabel.jl/compare/v0.2.0...v0.1.0
+[0.1.0]: https://github.com/oxfordcontrol/Clarabel.jl/tree/0.1.0

Project.toml

@@ -1,7 +1,7 @@
 name = "Clarabel"
 uuid = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e"
 authors = ["Paul Goulart <[email protected]>"]
-version = "0.4.0"
+version = "0.4.1"
 
 [deps]
 AMD = "14f7f29c-3bd6-536c-9a0b-7339e30b5a3e"

README.md

@@ -10,7 +10,7 @@ Interior Point Conic Optimization for Julia
   <a href="https://codecov.io/gh/oxfordcontrol/Clarabel.jl"><img src="https://codecov.io/gh/oxfordcontrol/Clarabel.jl/branch/main/graph/badge.svg"></a>
   <a href="https://oxfordcontrol.github.io/ClarabelDocs/stable"><img src="https://img.shields.io/badge/Documentation-stable-purple.svg"></a>
   <a href="https://opensource.org/licenses/Apache-2.0"><img src="https://img.shields.io/badge/License-Apache%202.0-blue.svg"></a>
-  <a href="https://github.com/oxfordcontrol/Clarabel.jl/releases"><img src="https://img.shields.io/badge/Release-v0.4.0-blue.svg"></a>
+  <a href="https://github.com/oxfordcontrol/Clarabel.jl/releases"><img src="https://img.shields.io/badge/Release-v0.4.1-blue.svg"></a>
 </p>
 
 <p align="center">

src/Clarabel.jl

@@ -6,6 +6,17 @@ module Clarabel
     const DefaultInt   = LinearAlgebra.BlasInt
     const IdentityMatrix = UniformScaling{Bool}
 
+    #internal constraint RHS limits.  This let block 
+    #hides the INFINITY field in the module and makes 
+    #it accessible only via the get/set provided
+    let 
+        _INFINITY_DEFAULT = 1e20
+        INFINITY = _INFINITY_DEFAULT
+        global default_infinity() = INFINITY = _INFINITY_DEFAULT;
+        global set_infinity(v::Float64) = INFINITY =  Float64(v)
+        global get_infinity() = INFINITY
+    end 
+
     #version / release info
     include("./version.jl")
 
@@ -20,6 +31,7 @@ module Clarabel
     include("./settings.jl")
     include("./conicvector.jl")
     include("./statuscodes.jl")
+    include("./presolver.jl")
     include("./types.jl")
     include("./variables.jl")
     include("./residuals.jl")

src/cones/compositecone_type.jl

@@ -8,10 +8,8 @@ struct CompositeCone{T} <: AbstractCone{T}
 
     cones::Vector{AbstractCone{T}}
 
-    #API type specs and count of each cone type
-    cone_specs::Vector{SupportedCone}
-    types::Vector{DataType}
-    type_counts::Dict{DataType,Int}
+    #count of each cone type
+    type_counts::Dict{Type,Int}
 
     #overall size of the composite cone
     numel::DefaultInt
@@ -24,38 +22,31 @@ struct CompositeCone{T} <: AbstractCone{T}
     # the flag for symmetric cone check
     _is_symmetric::Bool
 
-    function CompositeCone{T}(cone_specs::Vector{<:SupportedCone}) where {T}
-
-        #make copy to protect from user interference after setup,
-        #and explicitly cast as the abstract type.  This prevents
-        #errors arising from input vectors that are all the same cone,
-        #and therefore more concretely typed than we want
-        cone_specs = Vector{SupportedCone}(cone_specs)
+    function CompositeCone{T}(cone_specs::Vector{SupportedCone}) where {T}
 
         ncones = length(cone_specs)
         cones  = Vector{AbstractCone{T}}(undef,ncones)
-        types = Vector{DataType}(undef,ncones)
 
         #count the number of each cone type
-        type_counts = Dict{DataType,Int}()
-        for coneT in keys(ConeDict)
-            type_counts[coneT] = count(C->isa(C,coneT), cone_specs)
+        type_counts = Dict{Type,Int}()
+        for (key, val) in ConeDict
+            type_counts[val] = count(C->isa(C,key), cone_specs)
         end
 
         #assumed symmetric to start
         _is_symmetric = true
 
         #create cones with the given dims
-        for i in eachindex(cone_specs)
-            types[i] = typeof(cone_specs[i])
-            if types[i] == ExponentialConeT
-                cones[i] = ConeDict[typeof(cone_specs[i])]{T}()
+        for (i, coneT) in enumerate(cone_specs)
+            typeT = typeof(coneT)
+            if typeT == ExponentialConeT
+                cones[i] = ConeDict[typeT]{T}()
                 _is_symmetric = false
-            elseif types[i] == PowerConeT
-                cones[i] = ConeDict[typeof(cone_specs[i])]{T}(T(cone_specs[i].α))
+            elseif typeT == PowerConeT
+                cones[i] = ConeDict[typeT]{T}(T(cone_specs[i].α))
                 _is_symmetric = false
             else
-                cones[i] = ConeDict[typeof(cone_specs[i])]{T}(cone_specs[i].dim)
+                cones[i] = ConeDict[typeT]{T}(cone_specs[i].dim)
             end
         end
 
@@ -68,7 +59,7 @@ struct CompositeCone{T} <: AbstractCone{T}
         headidx = Vector{Int}(undef,length(cones))
         _make_headidx!(headidx,cones)
 
-        return new(cones,cone_specs,types,type_counts,numel,degree,headidx,_is_symmetric)
+        return new(cones,type_counts,numel,degree,headidx,_is_symmetric)
     end
 end
 

src/cones/cone_types.jl

@@ -47,7 +47,7 @@ struct NonnegativeCone{T} <: AbstractCone{T}
 
     function NonnegativeCone{T}(dim) where {T}
 
-        dim >= 1 || throw(DomainError(dim, "dimension must be positive"))
+        dim >= 0 || throw(DomainError(dim, "dimension must be nonnegative"))
         w = zeros(T,dim)
         λ = zeros(T,dim)
         return new(dim,w,λ)
@@ -112,9 +112,10 @@ mutable struct PSDConeWork{T}
     B::Matrix{T}
     Hs::Matrix{T}
 
-    #workspace for various internal use
+    #workspace for various internal uses
     workmat1::Matrix{T}
     workmat2::Matrix{T}
+    workmat3::Matrix{T}
     workvec::Vector{T}
 
     function PSDConeWork{T}(n::Int) where {T}
@@ -133,10 +134,11 @@ mutable struct PSDConeWork{T}
 
         workmat1 = zeros(T,n,n)
         workmat2 = zeros(T,n,n)
+        workmat3 = zeros(T,n,n)
         workvec  = zeros(T,(n*(n+1))>>1)
 
         return new(cholS,cholZ,SVD,λ,Λisqrt,R,Rinv,
-                   kronRR,B,Hs,workmat1,workmat2,workvec)
+                   kronRR,B,Hs,workmat1,workmat2,workmat3,workvec)
     end
 end
 
@@ -239,7 +241,7 @@ PowerCone(args...) = PowerCone{DefaultFloat}(args...)
 A Dict that maps the user-facing SupportedCone types to
 the types used internally in the solver.   See [SupportedCone](@ref)
 """
-const ConeDict = Dict(
+const ConeDict = Dict{DataType,Type}(
            ZeroConeT => ZeroCone,
     NonnegativeConeT => NonnegativeCone,
     SecondOrderConeT => SecondOrderCone,

src/cones/coneops_compositecone.jl

@@ -143,7 +143,7 @@ function update_scaling!(
 
     # update cone scalings by passing subview to each of
     # the appropriate cone types.
-    for (cone,type,si,zi) in zip(cones,cones.types,s.views,z.views)
+    for (cone,si,zi) in zip(cones,s.views,z.views)
         @conedispatch is_scaling_success = update_scaling!(cone,si,zi,μ,scaling_strategy)
         # YC: currently, only check whether SOC variables are in the interior;
         # we could extend the interior checkfor other cones
@@ -247,7 +247,7 @@ function step_length(
     s     = s.views
 
     # Force symmetric cones first.   
-    for (cone,type,dzi,dsi,zi,si) in zip(cones,cones.types,dz,ds,z,s)
+    for (cone,dzi,dsi,zi,si) in zip(cones,dz,ds,z,s)
 
         if !is_symmetric(cone) continue end
         @conedispatch (nextαz,nextαs) = step_length(cone,dzi,dsi,zi,si,settings,α)

src/cones/coneops_nncone.jl

@@ -9,7 +9,7 @@ function rectify_equilibration!(
 ) where{T}
 
     #allow elementwise equilibration scaling
-    δ .= e
+    δ .= one(T)
     return false
 end
 
@@ -19,7 +19,11 @@ function margins(
     pd::PrimalOrDualCone
 ) where{T}
 
-    α = minimum(z)
+    if(length(z) > 0)  
+        α = minimum(z)
+    else  #catch cone of dimension zero
+        α = floatmax(T)
+    end
 
     # total positive margin is the sum of the 
     # nonnegative elements in the vector, since 

src/cones/coneops_psdtrianglecone.jl

@@ -14,8 +14,9 @@ function margins(
 
     Z = K.work.workmat1
     _svec_to_mat!(Z,z,K)
-    α = eigvals(Symmetric(Z),1:1)[1]  #minimum eigenvalue
-    β = sum(eigvals(Symmetric(Z),0,floatmax(T)))  #sum of positive eigenvalues
+    e = eigvals!(Symmetric(Z))
+    α = minimum(e)  #minimum eigenvalue
+    β = reduce((x,y) -> y > 0 ? x + y : x, e, init = 0.) # = sum(e[e.>0]) (no alloc)
     (α,β)
 
 end
@@ -94,16 +95,21 @@ function update_scaling!(
 
     #R is the same size as L2'*L1,
     #so use it as temporary workspace
-    f.R  .= L2'*L1
+    mul!(f.R,L2',L1)   #R = L2'*L1
+
     f.SVD = svd(f.R)
 
     #assemble  λ (diagonal), R and Rinv.
     f.λ           .= f.SVD.S
     f.Λisqrt.diag .= inv.(sqrt.(f.λ))
 
-    # PJG : allocating 
-    f.R    .= L1*(f.SVD.V)*f.Λisqrt
-    f.Rinv .= f.Λisqrt*(f.SVD.U)'*L2'
+    #f.R = L1*(f.SVD.V)*f.Λisqrt
+    mul!(f.R,L1,f.SVD.V);
+    mul!(f.R,f.R,f.Λisqrt) #mul! can take Rinv twice because Λ is diagonal
+
+    #f.Rinv .= f.Λisqrt*(f.SVD.U)'*L2'
+    mul!(f.Rinv,f.SVD.U',L2')
+    mul!(f.Rinv,f.Λisqrt,f.Rinv) #mul! can take Rinv twice because Λ is diagonal
 
     return is_scaling_success = true
 end
@@ -258,14 +264,18 @@ function mul_W!(
 
   (X,Y) = (K.work.workmat1,K.work.workmat2)
   map((M,v)->_svec_to_mat!(M,v,K),(X,Y),(x,y))
+  tmp = K.work.workmat3
 
   R = K.work.R
 
-  # PJG : allocating 
   if is_transpose === :T
-      Y .+= α*(R*X*R')  #W^T*x
+      #Y .+= α*(R*X*R')  #W^T*x
+      mul!(tmp,X,R')
+      mul!(Y,R,tmp,α,one(T))
   else  # :N
-      Y .+= α*(R'*X*R)  #W*x
+      #Y .+= α*(R'*X*R)  #W*x
+      mul!(tmp,R',X)
+      mul!(Y,tmp,R,α,one(T))
   end
 
   _mat_to_svec!(y,Y,K)
@@ -287,14 +297,18 @@ function mul_Winv!(
 
     (X,Y) = (K.work.workmat1,K.work.workmat2)
     map((M,v)->_svec_to_mat!(M,v,K),(X,Y),(x,y))
+    tmp = K.work.workmat3
 
     Rinv = K.work.Rinv
 
-    # PJG : allocating 
     if is_transpose === :T
-        Y .+= α*(Rinv*X*Rinv')  #W^{-T}*x
+        #Y .+= α*(Rinv*X*Rinv')  #W^{-T}*x
+        mul!(tmp,X,Rinv')
+        mul!(Y,Rinv,tmp,α,one(T))
     else # :N
-        Y .+= α*(Rinv'*X*Rinv)  #W^{-1}*x
+        #Y .+= α*(Rinv'*X*Rinv)  #W^{-1}*x
+        mul!(tmp,Rinv',X)
+        mul!(Y,tmp,Rinv,α,one(T))
     end
 
     _mat_to_svec!(y,Y,K)
@@ -338,8 +352,13 @@ function circ_op!(
     (Y,Z) = (K.work.workmat1,K.work.workmat2)
     map((M,v)->_svec_to_mat!(M,v,K),(Y,Z),(y,z))
 
-    # PJG : allocating 
-    Y  .= (Y*Z + Z*Y)/2
+    tmp = K.work.workmat3;
+
+    #Y  .= (Y*Z + Z*Y)/2 
+    # NB: works b/c Y and Z are both symmetric
+    mul!(tmp,Y,Z)
+    Y .= tmp
+    symmetric_part!(Y)
     _mat_to_svec!(x,Y,K)
 
     return nothing
@@ -383,9 +402,7 @@ function _step_length_psd_component(
     # we only need to populate the upper 
     # triangle 
     lrscale!(Λisqrt.diag,Δ,Λisqrt.diag)
-    M = Symmetric(Δ)
-
-    γ = eigvals(M,1:1)[1] #minimum eigenvalue
+    γ = eigvals!(Symmetric(Δ),1:1)[1] #minimum eigenvalue
     if γ < 0
         return min(inv(-γ),αmax)
     else

src/cones/coneops_zerocone.jl

@@ -19,7 +19,7 @@ function rectify_equilibration!(
 ) where{T}
 
     #allow elementwise equilibration scaling
-    δ .= e
+    δ .= one(T)
     return false
 end
 

src/info_print.jl

@@ -17,18 +17,22 @@ function info_print_configuration(
 
     if(settings.verbose == false) return end
 
+    if(is_reduced(data.presolver))
+        @printf("\npresolve: removed %i constraints\n", count_reduced(data.presolver))
+    end 
+
     @printf("\nproblem:\n")
     @printf("  variables     = %i\n", data.n)
     @printf("  constraints   = %i\n", data.m)
     @printf("  nnz(P)        = %i\n", nnz(data.P))
     @printf("  nnz(A)        = %i\n", nnz(data.A))
     @printf("  cones (total) = %i\n", length(cones))
-    print_conedims_by_type(cones, ZeroConeT)
-    print_conedims_by_type(cones, NonnegativeConeT)
-    print_conedims_by_type(cones, SecondOrderConeT)
-    print_conedims_by_type(cones, PSDTriangleConeT)
-    print_conedims_by_type(cones, PowerConeT)
-    print_conedims_by_type(cones, ExponentialConeT)
+    print_conedims_by_type(cones, ZeroCone)
+    print_conedims_by_type(cones, NonnegativeCone)
+    print_conedims_by_type(cones, SecondOrderCone)
+    print_conedims_by_type(cones, PSDTriangleCone)
+    print_conedims_by_type(cones, PowerCone)
+    print_conedims_by_type(cones, ExponentialCone)
     print_settings(settings, T)
     @printf("\n")
 
@@ -176,8 +180,8 @@ function print_conedims_by_type(cones::CompositeCone{T}, type) where {T}
         return #don't report if none
     end
 
-    nvars = map(K->Clarabel.numel(K), cones[isa.(cones.cone_specs,type)])
-    name  = rpad(string(type)[1:end-5],11)  #drops "ConeT part"
+    nvars = map(K->Clarabel.numel(K), cones[isa.(cones,type)])
+    name  = rpad(string(nameof(type))[1:end-4],11)  #drops "Cone" part
     @printf("    : %s = %i, ", name, count)
 
     if count == 1