Merge pull request #114 from oxfordcontrol/releases/v0.4.0-dev

v0.4.0 release

CHANGELOG.md

@@ -0,0 +1,50 @@
+# Changelog
+
+Changes for the Julia version of Clarabel are documented in this file.   For the Rust version, see [here](https://github.com/oxfordcontrol/clarabel.rs/CHANGELOG.md).
+
+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.0] - 2023-25-02
+
+### Changed 
+
+- Internal fixes relating to initialization of iterates in symmetric cone problems.
+
+- Numerical stability improvements for second order cone constraints. 
+
+### Julia-specific changes
+
+- Modification of the internal calls to QDLDL.jl to allow for direct assignment of parameters in AMD ordering.   This release requires QDLDL.jl v0.4.0.
+
+- Makes Pardiso an optional dependency via Requires.jl.  To use Pardiso/MKL it is not necessary to import the Pardiso package directly before calling any part of the solver.  Fixes [#108](https://github.com/oxfordcontrol/Clarabel.jl/issues/108)
+
+
+## [0.3.0] - 2022-09-13
+
+### Changed 
+
+- Implements support for exponential and power cones
+
+- Numerical stability improvements
+
+- Various bug fixes
+
+## [0.2.0] - 2022-07-31
+
+- Companion rust/python implementation released starting from this version.
+
+- Ported all documentation to the common site [here](https://github.com/oxfordcontrol/ClarabelDocs)
+
+
+## [0.1.0] - 2022-07-04
+
+- 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

Project.toml

@@ -1,18 +1,18 @@
 name = "Clarabel"
 uuid = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e"
 authors = ["Paul Goulart <[email protected]>"]
-version = "0.3.0"
+version = "0.4.0"
 
 [deps]
 AMD = "14f7f29c-3bd6-536c-9a0b-7339e30b5a3e"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
-Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 QDLDL = "bfc457fd-c171-5ab7-bd9e-d5dbfc242d63"
+Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
@@ -23,9 +23,15 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 AMD = "0.4, 0.5"
 DataFrames = "1"
 MathOptInterface = "1.2"
-Pardiso = "0.5"
-PrettyTables = "0.12, 1"
-QDLDL = "0.3"
+PrettyTables = "1, 2"
+QDLDL = "0.4"
+Requires = "1"
 StaticArrays = "1"
 TimerOutputs = "0.5"
 julia = "1.2"
+
+[extras]
+Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2"
+
+[targets]
+test = ["Pardiso"]

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.3.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>
 </p>
 
 <p align="center">

src/Clarabel.jl

@@ -1,6 +1,7 @@
+__precompile__()
 module Clarabel
 
-    using SparseArrays, LinearAlgebra, Printf
+    using SparseArrays, LinearAlgebra, Printf, Requires
     const DefaultFloat = Float64
     const DefaultInt   = LinearAlgebra.BlasInt
     const IdentityMatrix = UniformScaling{Bool}
@@ -55,10 +56,22 @@ module Clarabel
     include("./utils/mathutils.jl")
     include("./utils/csc_assembly.jl")
 
+    #optional dependencies.  
+    #NB: This __init__ function and its @require statements 
+    #should be removed upon update of this package for use 
+    #with Julia v1.10+, after which weakdeps / external 
+    #dependencies will be natively supported 
+    function __init__()
+        @require Pardiso="46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2" begin
+            include("./kktsolvers/direct-ldl/directldl_mklpardiso.jl")  
+        end 
+    end
+
     #MathOptInterface for JuMP/Convex.jl
     module MOImodule
          include("./MOI_wrapper/MOI_wrapper.jl")
     end
     const Optimizer{T} = Clarabel.MOImodule.Optimizer{T}
 
 end
+

src/MOI_wrapper/MOI_wrapper.jl

@@ -90,6 +90,7 @@ mutable struct Optimizer{T} <: MOI.AbstractOptimizer
     solver_info::Union{Nothing,Clarabel.DefaultInfo{T}}
     solver_solution::Union{Nothing,Clarabel.DefaultSolution{T}}
     solver_nvars::Union{Nothing,Int}
+    use_quad_obj::Bool
     sense::MOI.OptimizationSense
     objconstant::T
     rowranges::Dict{Int, UnitRange{Int}}
@@ -101,10 +102,11 @@ mutable struct Optimizer{T} <: MOI.AbstractOptimizer
         solver_info     = nothing
         solver_solution = nothing
         solver_nvars    = nothing
+        use_quad_obj    = true
         sense = MOI.MIN_SENSE
         objconstant = zero(T)
         rowranges = Dict{Int, UnitRange{Int}}()
-        optimizer = new(solver_module,solver,solver_settings,solver_info,solver_solution,solver_nvars,sense,objconstant,rowranges)
+        optimizer = new(solver_module,solver,solver_settings,solver_info,solver_solution,solver_nvars,use_quad_obj,sense,objconstant,rowranges)
         for (key, value) in user_settings
             MOI.set(optimizer, MOI.RawOptimizerAttribute(string(key)), value)
         end
@@ -223,10 +225,27 @@ MOI.set(opt::Optimizer, ::MOI.Silent, v::Bool) = (opt.solver_settings.verbose =
 
 
 MOI.supports(::Optimizer, ::MOI.RawOptimizerAttribute) = true
-MOI.get(opt::Optimizer, param::MOI.RawOptimizerAttribute) =
-    getproperty(opt.solver_settings, Symbol(param.name))
-MOI.set(opt::Optimizer, param::MOI.RawOptimizerAttribute, value) =
-    setproperty!(opt.solver_settings, Symbol(param.name), value)
+
+function MOI.get(opt::Optimizer, param::MOI.RawOptimizerAttribute)
+
+    #catch wrapper level attributes.  #otherwise pass to solver
+    if(param.name ==  "use_quad_obj")
+        return opt.use_quad_obj
+    else
+        return getproperty(opt.solver_settings, Symbol(param.name))
+    end
+end
+
+
+function MOI.set(opt::Optimizer, param::MOI.RawOptimizerAttribute, value) 
+
+    #catch wrapper level attributes.  #otherwise pass to solver
+    if(param.name ==  "use_quad_obj")
+        opt.use_quad_obj = value
+    else
+        setproperty!(opt.solver_settings, Symbol(param.name), value)
+    end
+end
 
 MOI.supports(::Optimizer, ::MOI.VariablePrimal) = true
 function MOI.get(opt::Optimizer, a::MOI.VariablePrimal, vi::MOI.VariableIndex)
@@ -314,13 +333,23 @@ end
 # supported objective functions
 #------------------------------
 
-MOI.supports(
-    ::Optimizer{T},
+function MOI.supports(
+    opt::Optimizer{T},
     ::MOI.ObjectiveFunction{<:Union{
        MOI.ScalarAffineFunction{T},
+    }}
+) where {T} 
+    true
+end
+
+function MOI.supports(
+    opt::Optimizer{T},
+    ::MOI.ObjectiveFunction{<:Union{
        MOI.ScalarQuadraticFunction{T},
     }}
-) where {T} = true
+) where {T}
+    opt.use_quad_obj
+end
 
 
 #------------------------------

src/cones/cone_types.jl

@@ -4,6 +4,17 @@ using StaticArrays
 # -------------------------------------
 abstract type AbstractCone{T <: AbstractFloat} end
 
+#NB: this enum can't use Primal and Dual as its markers, 
+#since Dual is already used in the solve strategies.  Julia
+#won't allow a repeated use of the same marker in different 
+#enums, which is dumb.
+
+# marker for primal / dual distinctions 
+@enum PrimalOrDualCone begin 
+    PrimalCone
+    DualCone 
+end 
+
 # -------------------------------------
 # Zero Cone
 # -------------------------------------

src/cones/coneops_compositecone.jl

@@ -79,15 +79,33 @@ function rectify_equilibration!(
     return any_changed
 end
 
-# place a vector to some nearby point in the cone
-function shift_to_cone!(
+function margins(
     cones::CompositeCone{T},
-    z::ConicVector{T}
+    z::ConicVector{T},
+    pd::PrimalOrDualCone,
+) where {T}
+    α = typemax(T)
+    β = zero(T)
+    for (cone,zi) in zip(cones,z.views)
+        @conedispatch (αi,βi) = margins(cone,zi,pd)
+        α = min(α,αi)
+        β += βi
+    end
+
+    return (α,β)
+end
+
+function scaled_unit_shift!(
+    cones::CompositeCone{T},
+    z::ConicVector{T},
+    α::T,
+    pd::PrimalOrDualCone
 ) where {T}
 
     for (cone,zi) in zip(cones,z.views)
-        @conedispatch shift_to_cone!(cone,zi)
+        @conedispatch scaled_unit_shift!(cone,zi,α,pd)
     end
+
     return nothing
 end
 
@@ -125,11 +143,15 @@ function update_scaling!(
 
     # update cone scalings by passing subview to each of
     # the appropriate cone types.
-    for (cone,si,zi) in zip(cones,s.views,z.views)
-        @conedispatch update_scaling!(cone,si,zi,μ,scaling_strategy)
+    for (cone,type,si,zi) in zip(cones,cones.types,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
+        if !is_scaling_success
+            return is_scaling_success = false
+        end
     end
-
-    return nothing
+    return is_scaling_success = true
 end
 
 # The Hs block for each cone.
@@ -196,11 +218,12 @@ function Δs_from_Δz_offset!(
     cones::CompositeCone{T},
     out::ConicVector{T},
     ds::ConicVector{T},
-    work::ConicVector{T}
+    work::ConicVector{T},
+    z::ConicVector{T}
 ) where {T}
 
-    for (cone,outi,dsi,worki) in zip(cones,out.views,ds.views,work.views)
-        @conedispatch Δs_from_Δz_offset!(cone,outi,dsi,worki) 
+    for (cone,outi,dsi,worki,zi) in zip(cones,out.views,ds.views,work.views,z.views)
+        @conedispatch Δs_from_Δz_offset!(cone,outi,dsi,worki,zi) 
     end
 
     return nothing
@@ -233,13 +256,12 @@ function step_length(
 
     #if we have any nonsymmetric cones, then back off from full steps slightly
     #so that centrality checks and logarithms don't fail right at the boundaries
-    #PJG: is this still necessary?
     if(!is_symmetric(cones))
-        α = min(α,0.99)
+        α = min(α,settings.max_step_fraction)
     end
 
     # Force asymmetric cones last.  
-    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_defaults.jl

@@ -3,7 +3,7 @@ import Statistics: mean
 ## -------------------------------------------
 # Default implementations for cone operations
 # --------------------------------------------
-
+ 
 # degree of the cone is the same as dimension
 # and numel by default.   Degree is different
 # for the zero cone and SOC (0 and 1, respectively)
@@ -41,11 +41,29 @@ end
 # defined.   To define a new cone, you must
 # define implementations for each function below.
 
-# functions relating to unit vectors and cone initialization
+# returns (α,β) such that:
+# z - α⋅e is just on the cone boundary, with value
+# α >=0 indicates z \in cone, i.e. negative margin ===
+# outside of the cone.
+#
+# β is the sum of the margins that are positive.   For most 
+# cones this will just be β = max(0.,α), but for cones that 
+# are composites (e.g. the R_n^+), it is the sum of all of 
+# the positive margin terms.
+function margins(
+    K::AbstractCone{T},
+    z::AbstractVector{T},
+    pd::PrimalOrDualCone
+) where{T}
 
-function shift_to_cone!(
+    error("Incomplete cone operation specification: ",typeof(K))
+
+end
+
+function scaled_unit_shift!(
     K::AbstractCone{T},
-    z::AbstractVector{T}
+    z::AbstractVector{T},
+    α::T
 ) where{T}
 
     error("Incomplete cone operation specification: ",typeof(K))
@@ -80,6 +98,8 @@ function update_scaling!(
     scaling_strategy::ScalingStrategy
 ) where {T}
 
+    #NB: should return bool:  `true` on success.
+
     error("Incomplete cone operation specification: ",typeof(K))
 
 end
@@ -183,6 +203,7 @@ function Δs_from_Δz_offset!(
     out::AbstractVector{T},
     ds::AbstractVector{T},
     work::AbstractVector{T},
+    z::AbstractVector{T}
 ) where {T}
 
     error("Incomplete cone operation specification: ",typeof(K))
@@ -222,16 +243,6 @@ end
 # operations supported by symmetric cones only 
 # ---------------------------------------------
 
-# Add the scaled identity element e
-function add_scaled_e!(
-    K::AbstractCone{T},
-    x::AbstractVector{T},α::T
-) where {T}
-
-    error("Incomplete cone operation specification: ",typeof(K))
-
-end
-
 # implements y = αWx + βy
 function mul_W!(
     K::AbstractCone{T},