diff --git a/Project.toml b/Project.toml index 40e5d4f9f..c3193f319 100644 --- a/Project.toml +++ b/Project.toml @@ -34,6 +34,7 @@ DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" GSL = "92c85e6c-cbff-5e0c-80f7-495c94daaecd" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" PolyJuMP = "ddf597a6-d67e-5340-b84c-e37d84115374" @@ -42,4 +43,4 @@ SemialgebraicSets = "8e049039-38e8-557d-ae3a-bc521ccf6204" SumOfSquares = "4b9e565b-77fc-50a5-a571-1244f986bda1" [targets] -test = ["CSV", "DataFrames", "DataStructures", "DelimitedFiles", "Distributions", "DynamicPolynomials", "GSL", "JuMP", "PolyJuMP", "Random", "SemialgebraicSets", "SumOfSquares"] +test = ["CSV", "DataFrames", "DataStructures", "DelimitedFiles", "Distributions", "DynamicPolynomials", "ForwardDiff", "GSL", "JuMP", "PolyJuMP", "Random", "SemialgebraicSets", "SumOfSquares"] diff --git a/examples/common_JuMP.jl b/examples/common_JuMP.jl index cc78dc46b..a08e998bf 100644 --- a/examples/common_JuMP.jl +++ b/examples/common_JuMP.jl @@ -216,7 +216,7 @@ cone_from_hyp(cone::Cones.EpiPerSquare) = MOI.RotatedSecondOrderCone(Cones.dimen cone_from_hyp(cone::Cones.HypoPerLog) = (@assert Cones.dimension(cone) == 3; MOI.ExponentialCone()) cone_from_hyp(cone::Cones.EpiRelEntropy) = MOI.RelativeEntropyCone(Cones.dimension(cone)) cone_from_hyp(cone::Cones.HypoGeoMean) = MOI.GeometricMeanCone(Cones.dimension(cone)) -cone_from_hyp(cone::Cones.Power) = (@assert Cones.dimension(cone) == 3; MOI.PowerCone{Float64}(cone.alpha[1])) +cone_from_hyp(cone::Cones.GeneralizedPower) = (@assert Cones.dimension(cone) == 3; MOI.PowerCone{Float64}(cone.alpha[1])) cone_from_hyp(cone::Cones.EpiNormSpectral) = (Cones.use_dual_barrier(cone) ? MOI.NormNuclearCone : MOI.NormSpectralCone)(cone.n, cone.m) cone_from_hyp(cone::Cones.PosSemidefTri{T, R}) where {R <: Hypatia.RealOrComplex{T}} where {T <: Real} = MOI.PositiveSemidefiniteConeTriangle(cone.side) cone_from_hyp(cone::Cones.LinMatrixIneq{T}) where {T <: Real} = Hypatia.LinMatrixIneqCone{T}(cone.As) diff --git a/examples/matrixcompletion/native.jl b/examples/matrixcompletion/native.jl index 442bcf197..26c4d423e 100644 --- a/examples/matrixcompletion/native.jl +++ b/examples/matrixcompletion/native.jl @@ -195,21 +195,21 @@ function build(inst::MatrixCompletionNative{T}) where {T <: Real} G_geo_unknown[1, 1] = -1 G_geo_unknown[2, 2] = -1 G_geo_newvars[3, 1] = -1 - push!(cones, Cones.Power{T}(fill(inv(T(2)), 2), 1)) + push!(cones, Cones.GeneralizedPower{T}(fill(inv(T(2)), 2), 1)) offset = 4 # loop over new vars for i in 1:(num_unknown - 3) G_geo_newvars[offset + 2, i + 1] = -1 G_geo_newvars[offset + 1, i] = -1 G_geo_unknown[offset, i + 2] = -1 - push!(cones, Cones.Power{T}([inv(T(i + 2)), T(i + 1) / T(i + 2)], 1)) + push!(cones, Cones.GeneralizedPower{T}([inv(T(i + 2)), T(i + 1) / T(i + 2)], 1)) offset += 3 end # last row also special because hypograph variable is fixed G_geo_unknown[offset, num_unknown] = -1 G_geo_newvars[offset + 1, num_unknown - 2] = -1 - push!(cones, Cones.Power{T}([inv(T(num_unknown)), T(num_unknown - 1) / T(num_unknown)], 1)) + push!(cones, Cones.GeneralizedPower{T}([inv(T(num_unknown)), T(num_unknown - 1) / T(num_unknown)], 1)) h = vcat(h_norm, zeros(T, 3 * (num_unknown - 2)), T[0, 0, 1]) # if using extended with spectral objective G_geo needs to be prepadded with an epigraph variable diff --git a/examples/maxvolume/native.jl b/examples/maxvolume/native.jl index 903da1640..3db44c79f 100644 --- a/examples/maxvolume/native.jl +++ b/examples/maxvolume/native.jl @@ -49,14 +49,14 @@ function build(inst::MaxVolumeNative{T}) where {T <: Real} G_geo_orig[1, 1] = -1 G_geo_orig[2, 2] = -1 G_geo_newvars[3, 1] = -1 - push!(cones, Cones.Power{T}(fill(inv(T(2)), 2), 1)) + push!(cones, Cones.GeneralizedPower{T}(fill(inv(T(2)), 2), 1)) offset = 4 # loop over new vars for i in 1:(n - 3) G_geo_newvars[offset + 2, i + 1] = -1 G_geo_newvars[offset + 1, i] = -1 G_geo_orig[offset, i + 2] = -1 - push!(cones, Cones.Power{T}([inv(T(i + 2)), T(i + 1) / T(i + 2)], 1)) + push!(cones, Cones.GeneralizedPower{T}([inv(T(i + 2)), T(i + 1) / T(i + 2)], 1)) offset += 3 end # last row also special becuase hypograph variable is involved @@ -65,7 +65,7 @@ function build(inst::MaxVolumeNative{T}) where {T <: Real} G = [ vcat(zeros(T, len_power - 1), -one(T)) G_geo_orig G_geo_newvars ] - push!(cones, Cones.Power{T}([inv(T(n)), T(n - 1) / T(n)], 1)) + push!(cones, Cones.GeneralizedPower{T}([inv(T(n)), T(n - 1) / T(n)], 1)) h = zeros(T, 3 * (n - 1)) else @assert inst.use_epipersquare == true diff --git a/src/Cones/Cones.jl b/src/Cones/Cones.jl index 17e363f47..859031ca3 100644 --- a/src/Cones/Cones.jl +++ b/src/Cones/Cones.jl @@ -1,7 +1,5 @@ #= functions and caches for cones - -TODO maybe write a fallback dual feas check that checks if ray of dual point intersects dikin ellipsoid at primal point =# module Cones @@ -33,29 +31,29 @@ hessian_cache(T::Type{<:Real}) = DensePosDefCache{T}() abstract type Cone{T <: Real} end include("nonnegative.jl") +include("possemideftri.jl") +include("doublynonnegativetri.jl") +include("possemideftrisparse.jl") +include("linmatrixineq.jl") include("epinorminf.jl") include("epinormeucl.jl") -include("epiperentropy.jl") -include("epipertraceentropytri.jl") include("epipersquare.jl") -include("epirelentropy.jl") -include("hypoperlog.jl") -include("power.jl") -include("hypopowermean.jl") -include("hypogeomean.jl") include("epinormspectral.jl") -include("linmatrixineq.jl") -include("possemideftri.jl") -include("possemideftrisparse.jl") -include("doublynonnegativetri.jl") include("matrixepipersquare.jl") -include("hypoperlogdettri.jl") +include("generalizedpower.jl") +include("hypopowermean.jl") +include("hypogeomean.jl") include("hyporootdettri.jl") +include("hypoperlog.jl") +include("hypoperlogdettri.jl") +include("epiperentropy.jl") +include("epipertraceentropytri.jl") +include("epirelentropy.jl") include("epitracerelentropytri.jl") include("wsosinterpnonnegative.jl") +include("wsosinterppossemideftri.jl") include("wsosinterpepinormone.jl") include("wsosinterpepinormeucl.jl") -include("wsosinterppossemideftri.jl") use_dual_barrier(cone::Cone) = cone.use_dual_barrier dimension(cone::Cone) = cone.dim diff --git a/src/Cones/power.jl b/src/Cones/generalizedpower.jl similarity index 90% rename from src/Cones/power.jl rename to src/Cones/generalizedpower.jl index 38ade9715..1814d20b9 100644 --- a/src/Cones/power.jl +++ b/src/Cones/generalizedpower.jl @@ -7,7 +7,7 @@ barrier from "On self-concordant barriers for generalized power cones" by Roy & -log(prod_i((u_i)^(2 * alpha_i)) - norm_2(w)^2) - sum_i((1 - alpha_i)*log(u_i)) =# -mutable struct Power{T <: Real} <: Cone{T} +mutable struct GeneralizedPower{T <: Real} <: Cone{T} use_dual_barrier::Bool use_heuristic_neighborhood::Bool dim::Int @@ -37,7 +37,7 @@ mutable struct Power{T <: Real} <: Cone{T} auiproduuw::Vector{T} tempm::Vector{T} - function Power{T}( + function GeneralizedPower{T}( alpha::Vector{T}, n::Int; use_dual::Bool = false, @@ -60,10 +60,10 @@ mutable struct Power{T <: Real} <: Cone{T} end end -dimension(cone::Power) = length(cone.alpha) + cone.n +dimension(cone::GeneralizedPower) = length(cone.alpha) + cone.n # TODO only allocate the fields we use -function setup_extra_data(cone::Power{T}) where {T <: Real} +function setup_extra_data(cone::GeneralizedPower{T}) where {T <: Real} dim = cone.dim cone.hess = Symmetric(zeros(T, dim, dim), :U) cone.inv_hess = Symmetric(zeros(T, dim, dim), :U) @@ -75,16 +75,16 @@ function setup_extra_data(cone::Power{T}) where {T <: Real} return cone end -get_nu(cone::Power) = length(cone.alpha) + 1 +get_nu(cone::GeneralizedPower) = length(cone.alpha) + 1 -function set_initial_point(arr::AbstractVector, cone::Power) +function set_initial_point(arr::AbstractVector, cone::GeneralizedPower) m = length(cone.alpha) @. @views arr[1:m] = sqrt(1 + cone.alpha) @views arr[(m + 1):cone.dim] .= 0 return arr end -function update_feas(cone::Power{T}) where {T <: Real} +function update_feas(cone::GeneralizedPower{T}) where {T <: Real} @assert !cone.feas_updated m = length(cone.alpha) @views u = cone.point[1:m] @@ -101,7 +101,7 @@ function update_feas(cone::Power{T}) where {T <: Real} return cone.is_feas end -function is_dual_feas(cone::Power{T}) where {T <: Real} +function is_dual_feas(cone::GeneralizedPower{T}) where {T <: Real} alpha = cone.alpha m = length(cone.alpha) @views u = cone.dual_point[1:m] @@ -115,7 +115,7 @@ function is_dual_feas(cone::Power{T}) where {T <: Real} return false end -function update_grad(cone::Power) +function update_grad(cone::GeneralizedPower) @assert cone.is_feas m = length(cone.alpha) @views u = cone.point[1:m] @@ -132,7 +132,7 @@ function update_grad(cone::Power) return cone.grad end -function update_hess(cone::Power) +function update_hess(cone::GeneralizedPower) @assert cone.grad_updated m = length(cone.alpha) @views u = cone.point[1:m] @@ -167,7 +167,7 @@ function update_hess(cone::Power) return cone.hess end -function hess_prod!(prod::AbstractVecOrMat, arr::AbstractVecOrMat, cone::Power) +function hess_prod!(prod::AbstractVecOrMat, arr::AbstractVecOrMat, cone::GeneralizedPower) @assert cone.grad_updated m = length(cone.alpha) @views u = cone.point[1:m] @@ -194,7 +194,7 @@ function hess_prod!(prod::AbstractVecOrMat, arr::AbstractVecOrMat, cone::Power) return prod end -function correction(cone::Power, primal_dir::AbstractVector) +function correction(cone::GeneralizedPower, primal_dir::AbstractVector) m = length(cone.alpha) @views u = cone.point[1:m] @views w = cone.point[(m + 1):end] diff --git a/src/MathOptInterface/cones.jl b/src/MathOptInterface/cones.jl index 7e79dad04..0fe886407 100644 --- a/src/MathOptInterface/cones.jl +++ b/src/MathOptInterface/cones.jl @@ -7,21 +7,21 @@ cone_from_moi(::Type{<:Real}, cone::MOI.AbstractVectorSet) = error("MOI set $con # MOI predefined cones +cone_from_moi(::Type{T}, cone::MOI.PositiveSemidefiniteConeTriangle) where {T <: Real} = Cones.PosSemidefTri{T, T}(MOI.dimension(cone)) cone_from_moi(::Type{T}, cone::MOI.NormInfinityCone) where {T <: Real} = Cones.EpiNormInf{T, T}(MOI.dimension(cone)) cone_from_moi(::Type{T}, cone::MOI.NormOneCone) where {T <: Real} = Cones.EpiNormInf{T, T}(MOI.dimension(cone), use_dual = true) cone_from_moi(::Type{T}, cone::MOI.SecondOrderCone) where {T <: Real} = Cones.EpiNormEucl{T}(MOI.dimension(cone)) cone_from_moi(::Type{T}, cone::MOI.RotatedSecondOrderCone) where {T <: Real} = Cones.EpiPerSquare{T}(MOI.dimension(cone)) -cone_from_moi(::Type{T}, cone::MOI.ExponentialCone) where {T <: Real} = Cones.HypoPerLog{T}(3) -cone_from_moi(::Type{T}, cone::MOI.DualExponentialCone) where {T <: Real} = Cones.HypoPerLog{T}(3, use_dual = true) -cone_from_moi(::Type{T}, cone::MOI.PowerCone{T}) where {T <: Real} = Cones.Power{T}(T[cone.exponent, 1 - cone.exponent], 1) -cone_from_moi(::Type{T}, cone::MOI.DualPowerCone{T}) where {T <: Real} = Cones.Power{T}(T[cone.exponent, 1 - cone.exponent], 1, use_dual = true) -cone_from_moi(::Type{T}, cone::MOI.GeometricMeanCone) where {T <: Real} = (l = MOI.dimension(cone) - 1; Cones.HypoGeoMean{T}(1 + l)) -cone_from_moi(::Type{T}, cone::MOI.RelativeEntropyCone) where {T <: Real} = Cones.EpiRelEntropy{T}(MOI.dimension(cone)) cone_from_moi(::Type{T}, cone::MOI.NormSpectralCone) where {T <: Real} = Cones.EpiNormSpectral{T, T}(extrema((cone.row_dim, cone.column_dim))...) cone_from_moi(::Type{T}, cone::MOI.NormNuclearCone) where {T <: Real} = Cones.EpiNormSpectral{T, T}(extrema((cone.row_dim, cone.column_dim))..., use_dual = true) -cone_from_moi(::Type{T}, cone::MOI.PositiveSemidefiniteConeTriangle) where {T <: Real} = Cones.PosSemidefTri{T, T}(MOI.dimension(cone)) -cone_from_moi(::Type{T}, cone::MOI.LogDetConeTriangle) where {T <: Real} = Cones.HypoPerLogdetTri{T, T}(MOI.dimension(cone)) +cone_from_moi(::Type{T}, cone::MOI.PowerCone{T}) where {T <: Real} = Cones.GeneralizedPower{T}(T[cone.exponent, 1 - cone.exponent], 1) +cone_from_moi(::Type{T}, cone::MOI.DualPowerCone{T}) where {T <: Real} = Cones.GeneralizedPower{T}(T[cone.exponent, 1 - cone.exponent], 1, use_dual = true) +cone_from_moi(::Type{T}, cone::MOI.GeometricMeanCone) where {T <: Real} = (l = MOI.dimension(cone) - 1; Cones.HypoGeoMean{T}(1 + l)) cone_from_moi(::Type{T}, cone::MOI.RootDetConeTriangle) where {T <: Real} = Cones.HypoRootdetTri{T, T}(MOI.dimension(cone)) +cone_from_moi(::Type{T}, cone::MOI.ExponentialCone) where {T <: Real} = Cones.HypoPerLog{T}(3) +cone_from_moi(::Type{T}, cone::MOI.DualExponentialCone) where {T <: Real} = Cones.HypoPerLog{T}(3, use_dual = true) +cone_from_moi(::Type{T}, cone::MOI.LogDetConeTriangle) where {T <: Real} = Cones.HypoPerLogdetTri{T, T}(MOI.dimension(cone)) +cone_from_moi(::Type{T}, cone::MOI.RelativeEntropyCone) where {T <: Real} = Cones.EpiRelEntropy{T}(MOI.dimension(cone)) # transformations fallbacks needs_untransform(::MOI.AbstractVectorSet) = false @@ -125,6 +125,43 @@ end MOI.dimension(cone::NonnegativeCone) = cone.dim cone_from_moi(::Type{T}, cone::NonnegativeCone{T}) where {T <: Real} = Cones.Nonnegative{T}(cone.dim) +export PosSemidefTriCone +struct PosSemidefTriCone{T <: Real, R <: RealOrComplex{T}} <: MOI.AbstractVectorSet + dim::Int +end +MOI.dimension(cone::PosSemidefTriCone) = cone.dim +cone_from_moi(::Type{T}, cone::PosSemidefTriCone{T, R}) where {R <: RealOrComplex{T}} where {T <: Real} = Cones.PosSemidefTri{T, R}(cone.dim) + +export DoublyNonnegativeTriCone +struct DoublyNonnegativeTriCone{T <: Real} <: MOI.AbstractVectorSet + dim::Int + use_dual::Bool +end +DoublyNonnegativeTriCone{T}(dim::Int) where {T <: Real} = DoublyNonnegativeTriCone{T}(dim, false) +MOI.dimension(cone::DoublyNonnegativeTriCone where {T <: Real}) = cone.dim +cone_from_moi(::Type{T}, cone::DoublyNonnegativeTriCone{T}) where {T <: Real} = Cones.DoublyNonnegativeTri{T}(cone.dim, use_dual = cone.use_dual) + +export PosSemidefTriSparseCone +struct PosSemidefTriSparseCone{T <: Real, R <: RealOrComplex{T}} <: MOI.AbstractVectorSet + side::Int + row_idxs::Vector{Int} + col_idxs::Vector{Int} + use_dual::Bool +end +PosSemidefTriSparseCone{T, R}(side::Int, row_idxs::Vector{Int}, col_idxs::Vector{Int}) where {R <: RealOrComplex{T}} where {T <: Real} = PosSemidefTriSparseCone{T, R}(side, row_idxs, col_idxs, false) +MOI.dimension(cone::PosSemidefTriSparseCone{T, T} where {T <: Real}) = length(cone.row_idxs) +MOI.dimension(cone::PosSemidefTriSparseCone{T, Complex{T}} where {T <: Real}) = (2 * length(cone.row_idxs) - cone.side) +cone_from_moi(::Type{T}, cone::PosSemidefTriSparseCone{T, R}) where {R <: RealOrComplex{T}} where {T <: Real} = Cones.PosSemidefTriSparse{T, R}(cone.side, cone.row_idxs, cone.col_idxs, use_dual = cone.use_dual) + +export LinMatrixIneqCone +struct LinMatrixIneqCone{T <: Real} <: MOI.AbstractVectorSet + As::Vector + use_dual::Bool +end +LinMatrixIneqCone{T}(As::Vector) where {T <: Real} = LinMatrixIneqCone{T}(As, false) +MOI.dimension(cone::LinMatrixIneqCone) = length(cone.As) +cone_from_moi(::Type{T}, cone::LinMatrixIneqCone{T}) where {T <: Real} = Cones.LinMatrixIneq{T}(cone.As, use_dual = cone.use_dual) + export EpiNormInfinityCone struct EpiNormInfinityCone{T <: Real, R <: RealOrComplex{T}} <: MOI.AbstractVectorSet dim::Int @@ -148,61 +185,6 @@ end MOI.dimension(cone::EpiPerSquareCone) = cone.dim cone_from_moi(::Type{T}, cone::EpiPerSquareCone{T}) where {T <: Real} = Cones.EpiPerSquare{T}(cone.dim) -export PowerCone -struct PowerCone{T <: Real} <: MOI.AbstractVectorSet - alpha::Vector{T} - n::Int - use_dual::Bool -end -PowerCone{T}(alpha::Vector{T}, n::Int) where {T <: Real} = PowerCone{T}(alpha, n, false) -MOI.dimension(cone::PowerCone) = length(cone.alpha) + cone.n -cone_from_moi(::Type{T}, cone::PowerCone{T}) where {T <: Real} = Cones.Power{T}(cone.alpha, cone.n, use_dual = cone.use_dual) - -export HypoPerLogCone -struct HypoPerLogCone{T <: Real} <: MOI.AbstractVectorSet - dim::Int - use_dual::Bool -end -HypoPerLogCone{T}(dim::Int) where {T <: Real} = HypoPerLogCone{T}(dim, false) -MOI.dimension(cone::HypoPerLogCone) = cone.dim -cone_from_moi(::Type{T}, cone::HypoPerLogCone{T}) where {T <: Real} = Cones.HypoPerLog{T}(cone.dim, use_dual = cone.use_dual) - -export EpiRelEntropyCone -struct EpiRelEntropyCone{T <: Real} <: MOI.AbstractVectorSet - dim::Int - use_dual::Bool -end -EpiRelEntropyCone{T}(dim::Int) where {T <: Real} = EpiRelEntropyCone{T}(dim, false) -MOI.dimension(cone::EpiRelEntropyCone) = cone.dim -cone_from_moi(::Type{T}, cone::EpiRelEntropyCone{T}) where {T <: Real} = Cones.EpiRelEntropy{T}(cone.dim, use_dual = cone.use_dual) - -export EpiPerEntropyCone -struct EpiPerEntropyCone{T <: Real} <: MOI.AbstractVectorSet - dim::Int - use_dual::Bool -end -EpiPerEntropyCone{T}(dim::Int) where {T <: Real} = EpiPerEntropyCone{T}(dim, false) -MOI.dimension(cone::EpiPerEntropyCone) = cone.dim -cone_from_moi(::Type{T}, cone::EpiPerEntropyCone{T}) where {T <: Real} = Cones.EpiPerEntropy{T}(cone.dim, use_dual = cone.use_dual) - -export HypoGeoMeanCone -struct HypoGeoMeanCone{T <: Real} <: MOI.AbstractVectorSet - dim::Int - use_dual::Bool -end -HypoGeoMeanCone{T}(dim::Int) where {T <: Real} = HypoGeoMeanCone{T}(dim, false) -MOI.dimension(cone::HypoGeoMeanCone) = cone.dim -cone_from_moi(::Type{T}, cone::HypoGeoMeanCone{T}) where {T <: Real} = Cones.HypoGeoMean{T}(cone.dim, use_dual = cone.use_dual) - -export HypoPowerMeanCone -struct HypoPowerMeanCone{T <: Real} <: MOI.AbstractVectorSet - alpha::Vector{T} - use_dual::Bool -end -HypoPowerMeanCone{T}(alpha::Vector{T}) where {T <: Real} = HypoPowerMeanCone{T}(alpha, false) -MOI.dimension(cone::HypoPowerMeanCone) = 1 + length(cone.alpha) -cone_from_moi(::Type{T}, cone::HypoPowerMeanCone{T}) where {T <: Real} = Cones.HypoPowerMean{T}(cone.alpha, use_dual = cone.use_dual) - export EpiNormSpectralCone struct EpiNormSpectralCone{T <: Real, R <: RealOrComplex{T}} <: MOI.AbstractVectorSet n::Int @@ -225,42 +207,51 @@ MOI.dimension(cone::MatrixEpiPerSquareCone{T, T} where {T <: Real}) = Cones.svec MOI.dimension(cone::MatrixEpiPerSquareCone{T, Complex{T}} where {T <: Real}) = cone.n ^ 2 + 1 + 2 * cone.n * cone.m cone_from_moi(::Type{T}, cone::MatrixEpiPerSquareCone{T, R}) where {R <: RealOrComplex{T}} where {T <: Real} = Cones.MatrixEpiPerSquare{T, R}(cone.n, cone.m, use_dual = cone.use_dual) -export LinMatrixIneqCone -struct LinMatrixIneqCone{T <: Real} <: MOI.AbstractVectorSet - As::Vector +export GeneralizedPowerCone +struct GeneralizedPowerCone{T <: Real} <: MOI.AbstractVectorSet + alpha::Vector{T} + n::Int use_dual::Bool end -LinMatrixIneqCone{T}(As::Vector) where {T <: Real} = LinMatrixIneqCone{T}(As, false) -MOI.dimension(cone::LinMatrixIneqCone) = length(cone.As) -cone_from_moi(::Type{T}, cone::LinMatrixIneqCone{T}) where {T <: Real} = Cones.LinMatrixIneq{T}(cone.As, use_dual = cone.use_dual) +GeneralizedPowerCone{T}(alpha::Vector{T}, n::Int) where {T <: Real} = GeneralizedPowerCone{T}(alpha, n, false) +MOI.dimension(cone::GeneralizedPowerCone) = length(cone.alpha) + cone.n +cone_from_moi(::Type{T}, cone::GeneralizedPowerCone{T}) where {T <: Real} = Cones.GeneralizedPower{T}(cone.alpha, cone.n, use_dual = cone.use_dual) -export PosSemidefTriCone -struct PosSemidefTriCone{T <: Real, R <: RealOrComplex{T}} <: MOI.AbstractVectorSet +export HypoPowerMeanCone +struct HypoPowerMeanCone{T <: Real} <: MOI.AbstractVectorSet + alpha::Vector{T} + use_dual::Bool +end +HypoPowerMeanCone{T}(alpha::Vector{T}) where {T <: Real} = HypoPowerMeanCone{T}(alpha, false) +MOI.dimension(cone::HypoPowerMeanCone) = 1 + length(cone.alpha) +cone_from_moi(::Type{T}, cone::HypoPowerMeanCone{T}) where {T <: Real} = Cones.HypoPowerMean{T}(cone.alpha, use_dual = cone.use_dual) + +export HypoGeoMeanCone +struct HypoGeoMeanCone{T <: Real} <: MOI.AbstractVectorSet dim::Int + use_dual::Bool end -MOI.dimension(cone::PosSemidefTriCone) = cone.dim -cone_from_moi(::Type{T}, cone::PosSemidefTriCone{T, R}) where {R <: RealOrComplex{T}} where {T <: Real} = Cones.PosSemidefTri{T, R}(cone.dim) +HypoGeoMeanCone{T}(dim::Int) where {T <: Real} = HypoGeoMeanCone{T}(dim, false) +MOI.dimension(cone::HypoGeoMeanCone) = cone.dim +cone_from_moi(::Type{T}, cone::HypoGeoMeanCone{T}) where {T <: Real} = Cones.HypoGeoMean{T}(cone.dim, use_dual = cone.use_dual) -export PosSemidefTriSparseCone -struct PosSemidefTriSparseCone{T <: Real, R <: RealOrComplex{T}} <: MOI.AbstractVectorSet - side::Int - row_idxs::Vector{Int} - col_idxs::Vector{Int} +export HypoRootdetTriCone +struct HypoRootdetTriCone{T <: Real, R <: RealOrComplex{T}} <: MOI.AbstractVectorSet + dim::Int use_dual::Bool end -PosSemidefTriSparseCone{T, R}(side::Int, row_idxs::Vector{Int}, col_idxs::Vector{Int}) where {R <: RealOrComplex{T}} where {T <: Real} = PosSemidefTriSparseCone{T, R}(side, row_idxs, col_idxs, false) -MOI.dimension(cone::PosSemidefTriSparseCone{T, T} where {T <: Real}) = length(cone.row_idxs) -MOI.dimension(cone::PosSemidefTriSparseCone{T, Complex{T}} where {T <: Real}) = (2 * length(cone.row_idxs) - cone.side) -cone_from_moi(::Type{T}, cone::PosSemidefTriSparseCone{T, R}) where {R <: RealOrComplex{T}} where {T <: Real} = Cones.PosSemidefTriSparse{T, R}(cone.side, cone.row_idxs, cone.col_idxs, use_dual = cone.use_dual) +HypoRootdetTriCone{T, R}(dim::Int) where {R <: RealOrComplex{T}} where {T <: Real} = HypoRootdetTriCone{T, R}(dim, false) +MOI.dimension(cone::HypoRootdetTriCone where {T <: Real}) = cone.dim +cone_from_moi(::Type{T}, cone::HypoRootdetTriCone{T, R}) where {R <: RealOrComplex{T}} where {T <: Real} = Cones.HypoRootdetTri{T, R}(cone.dim, use_dual = cone.use_dual) -export DoublyNonnegativeTriCone -struct DoublyNonnegativeTriCone{T <: Real} <: MOI.AbstractVectorSet +export HypoPerLogCone +struct HypoPerLogCone{T <: Real} <: MOI.AbstractVectorSet dim::Int use_dual::Bool end -DoublyNonnegativeTriCone{T}(dim::Int) where {T <: Real} = DoublyNonnegativeTriCone{T}(dim, false) -MOI.dimension(cone::DoublyNonnegativeTriCone where {T <: Real}) = cone.dim -cone_from_moi(::Type{T}, cone::DoublyNonnegativeTriCone{T}) where {T <: Real} = Cones.DoublyNonnegativeTri{T}(cone.dim, use_dual = cone.use_dual) +HypoPerLogCone{T}(dim::Int) where {T <: Real} = HypoPerLogCone{T}(dim, false) +MOI.dimension(cone::HypoPerLogCone) = cone.dim +cone_from_moi(::Type{T}, cone::HypoPerLogCone{T}) where {T <: Real} = Cones.HypoPerLog{T}(cone.dim, use_dual = cone.use_dual) export HypoPerLogdetTriCone struct HypoPerLogdetTriCone{T <: Real, R <: RealOrComplex{T}} <: MOI.AbstractVectorSet @@ -271,23 +262,14 @@ HypoPerLogdetTriCone{T, R}(dim::Int) where {R <: RealOrComplex{T}} where {T <: R MOI.dimension(cone::HypoPerLogdetTriCone) = cone.dim cone_from_moi(::Type{T}, cone::HypoPerLogdetTriCone{T, R}) where {R <: RealOrComplex{T}} where {T <: Real} = Cones.HypoPerLogdetTri{T, R}(cone.dim, use_dual = cone.use_dual) -export HypoRootdetTriCone -struct HypoRootdetTriCone{T <: Real, R <: RealOrComplex{T}} <: MOI.AbstractVectorSet - dim::Int - use_dual::Bool -end -HypoRootdetTriCone{T, R}(dim::Int) where {R <: RealOrComplex{T}} where {T <: Real} = HypoRootdetTriCone{T, R}(dim, false) -MOI.dimension(cone::HypoRootdetTriCone where {T <: Real}) = cone.dim -cone_from_moi(::Type{T}, cone::HypoRootdetTriCone{T, R}) where {R <: RealOrComplex{T}} where {T <: Real} = Cones.HypoRootdetTri{T, R}(cone.dim, use_dual = cone.use_dual) - -export EpiTraceRelEntropyTriCone -struct EpiTraceRelEntropyTriCone{T <: Real} <: MOI.AbstractVectorSet +export EpiPerEntropyCone +struct EpiPerEntropyCone{T <: Real} <: MOI.AbstractVectorSet dim::Int use_dual::Bool end -EpiTraceRelEntropyTriCone{T}(dim::Int) where {T <: Real} = EpiTraceRelEntropyTriCone{T}(dim, false) -MOI.dimension(cone::EpiTraceRelEntropyTriCone where {T <: Real}) = cone.dim -cone_from_moi(::Type{T}, cone::EpiTraceRelEntropyTriCone{T}) where {T <: Real} = Cones.EpiTraceRelEntropyTri{T}(cone.dim, use_dual = cone.use_dual) +EpiPerEntropyCone{T}(dim::Int) where {T <: Real} = EpiPerEntropyCone{T}(dim, false) +MOI.dimension(cone::EpiPerEntropyCone) = cone.dim +cone_from_moi(::Type{T}, cone::EpiPerEntropyCone{T}) where {T <: Real} = Cones.EpiPerEntropy{T}(cone.dim, use_dual = cone.use_dual) export EpiPerTraceEntropyTriCone struct EpiPerTraceEntropyTriCone{T <: Real} <: MOI.AbstractVectorSet @@ -298,6 +280,24 @@ EpiPerTraceEntropyTriCone{T}(dim::Int) where {T <: Real} = EpiPerTraceEntropyTri MOI.dimension(cone::EpiPerTraceEntropyTriCone where {T <: Real}) = cone.dim cone_from_moi(::Type{T}, cone::EpiPerTraceEntropyTriCone{T}) where {T <: Real} = Cones.EpiPerTraceEntropyTri{T}(cone.dim, use_dual = cone.use_dual) +export EpiRelEntropyCone +struct EpiRelEntropyCone{T <: Real} <: MOI.AbstractVectorSet + dim::Int + use_dual::Bool +end +EpiRelEntropyCone{T}(dim::Int) where {T <: Real} = EpiRelEntropyCone{T}(dim, false) +MOI.dimension(cone::EpiRelEntropyCone) = cone.dim +cone_from_moi(::Type{T}, cone::EpiRelEntropyCone{T}) where {T <: Real} = Cones.EpiRelEntropy{T}(cone.dim, use_dual = cone.use_dual) + +export EpiTraceRelEntropyTriCone +struct EpiTraceRelEntropyTriCone{T <: Real} <: MOI.AbstractVectorSet + dim::Int + use_dual::Bool +end +EpiTraceRelEntropyTriCone{T}(dim::Int) where {T <: Real} = EpiTraceRelEntropyTriCone{T}(dim, false) +MOI.dimension(cone::EpiTraceRelEntropyTriCone where {T <: Real}) = cone.dim +cone_from_moi(::Type{T}, cone::EpiTraceRelEntropyTriCone{T}) where {T <: Real} = Cones.EpiTraceRelEntropyTri{T}(cone.dim, use_dual = cone.use_dual) + export WSOSInterpNonnegativeCone struct WSOSInterpNonnegativeCone{T <: Real, R <: RealOrComplex{T}} <: MOI.AbstractVectorSet U::Int @@ -319,63 +319,63 @@ WSOSInterpPosSemidefTriCone{T}(R::Int, U::Int, Ps::Vector{Matrix{T}}) where {T < MOI.dimension(cone::WSOSInterpPosSemidefTriCone) = cone.U * Cones.svec_length(cone.R) cone_from_moi(::Type{T}, cone::WSOSInterpPosSemidefTriCone{T}) where {T <: Real} = Cones.WSOSInterpPosSemidefTri{T}(cone.R, cone.U, cone.Ps, use_dual = cone.use_dual) -export WSOSInterpEpiNormEuclCone -struct WSOSInterpEpiNormEuclCone{T <: Real} <: MOI.AbstractVectorSet +export WSOSInterpEpiNormOneCone +struct WSOSInterpEpiNormOneCone{T <: Real} <: MOI.AbstractVectorSet R::Int U::Int Ps::Vector{Matrix{T}} use_dual::Bool end -WSOSInterpEpiNormEuclCone{T}(R::Int, U::Int, Ps::Vector{Matrix{T}}) where {T <: Real} = WSOSInterpEpiNormEuclCone{T}(R, U, Ps, false) -MOI.dimension(cone::WSOSInterpEpiNormEuclCone) = cone.U * cone.R -cone_from_moi(::Type{T}, cone::WSOSInterpEpiNormEuclCone{T}) where {T <: Real} = Cones.WSOSInterpEpiNormEucl{T}(cone.R, cone.U, cone.Ps, use_dual = cone.use_dual) +WSOSInterpEpiNormOneCone{T}(R::Int, U::Int, Ps::Vector{Matrix{T}}) where {T <: Real} = WSOSInterpEpiNormOneCone{T}(R, U, Ps, false) +MOI.dimension(cone::WSOSInterpEpiNormOneCone) = cone.U * cone.R +cone_from_moi(::Type{T}, cone::WSOSInterpEpiNormOneCone{T}) where {T <: Real} = Cones.WSOSInterpEpiNormOne{T}(cone.R, cone.U, cone.Ps, use_dual = cone.use_dual) -export WSOSInterpEpiNormOneCone -struct WSOSInterpEpiNormOneCone{T <: Real} <: MOI.AbstractVectorSet +export WSOSInterpEpiNormEuclCone +struct WSOSInterpEpiNormEuclCone{T <: Real} <: MOI.AbstractVectorSet R::Int U::Int Ps::Vector{Matrix{T}} use_dual::Bool end -WSOSInterpEpiNormOneCone{T}(R::Int, U::Int, Ps::Vector{Matrix{T}}) where {T <: Real} = WSOSInterpEpiNormOneCone{T}(R, U, Ps, false) -MOI.dimension(cone::WSOSInterpEpiNormOneCone) = cone.U * cone.R -cone_from_moi(::Type{T}, cone::WSOSInterpEpiNormOneCone{T}) where {T <: Real} = Cones.WSOSInterpEpiNormOne{T}(cone.R, cone.U, cone.Ps, use_dual = cone.use_dual) +WSOSInterpEpiNormEuclCone{T}(R::Int, U::Int, Ps::Vector{Matrix{T}}) where {T <: Real} = WSOSInterpEpiNormEuclCone{T}(R, U, Ps, false) +MOI.dimension(cone::WSOSInterpEpiNormEuclCone) = cone.U * cone.R +cone_from_moi(::Type{T}, cone::WSOSInterpEpiNormEuclCone{T}) where {T <: Real} = Cones.WSOSInterpEpiNormEucl{T}(cone.R, cone.U, cone.Ps, use_dual = cone.use_dual) # all cones const HypatiaCones{T <: Real} = Union{ NonnegativeCone{T}, + PosSemidefTriCone{T, T}, + PosSemidefTriCone{T, Complex{T}}, + DoublyNonnegativeTriCone{T}, + PosSemidefTriSparseCone{T, T}, + PosSemidefTriSparseCone{T, Complex{T}}, + LinMatrixIneqCone{T}, EpiNormInfinityCone{T, T}, EpiNormInfinityCone{T, Complex{T}}, EpiNormEuclCone{T}, EpiPerSquareCone{T}, - PowerCone{T}, - HypoPerLogCone{T}, - EpiPerEntropyCone{T}, - EpiRelEntropyCone{T}, - HypoGeoMeanCone{T}, - HypoPowerMeanCone{T}, EpiNormSpectralCone{T, T}, EpiNormSpectralCone{T, Complex{T}}, MatrixEpiPerSquareCone{T, T}, MatrixEpiPerSquareCone{T, Complex{T}}, - LinMatrixIneqCone{T}, - PosSemidefTriCone{T, T}, - PosSemidefTriCone{T, Complex{T}}, - PosSemidefTriSparseCone{T, T}, - PosSemidefTriSparseCone{T, Complex{T}}, - DoublyNonnegativeTriCone{T}, - HypoPerLogdetTriCone{T, T}, - HypoPerLogdetTriCone{T, Complex{T}}, + GeneralizedPowerCone{T}, + HypoPowerMeanCone{T}, + HypoGeoMeanCone{T}, HypoRootdetTriCone{T, T}, HypoRootdetTriCone{T, Complex{T}}, - EpiTraceRelEntropyTriCone{T}, + HypoPerLogCone{T}, + HypoPerLogdetTriCone{T, T}, + HypoPerLogdetTriCone{T, Complex{T}}, + EpiPerEntropyCone{T}, EpiPerTraceEntropyTriCone{T}, + EpiRelEntropyCone{T}, + EpiTraceRelEntropyTriCone{T}, WSOSInterpNonnegativeCone{T, T}, WSOSInterpNonnegativeCone{T, Complex{T}}, WSOSInterpPosSemidefTriCone{T}, - WSOSInterpEpiNormEuclCone{T}, WSOSInterpEpiNormOneCone{T}, + WSOSInterpEpiNormEuclCone{T}, } const SupportedCones{T <: Real} = Union{ @@ -383,21 +383,21 @@ const SupportedCones{T <: Real} = Union{ MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives, + MOI.PositiveSemidefiniteConeTriangle, MOI.NormInfinityCone, MOI.NormOneCone, MOI.SecondOrderCone, MOI.RotatedSecondOrderCone, - MOI.ExponentialCone, - MOI.DualExponentialCone, + MOI.NormSpectralCone, + MOI.NormNuclearCone, MOI.PowerCone{T}, MOI.DualPowerCone{T}, MOI.GeometricMeanCone, - MOI.RelativeEntropyCone, - MOI.NormSpectralCone, - MOI.NormNuclearCone, - MOI.PositiveSemidefiniteConeTriangle, - MOI.LogDetConeTriangle, MOI.RootDetConeTriangle, + MOI.ExponentialCone, + MOI.DualExponentialCone, + MOI.LogDetConeTriangle, + MOI.RelativeEntropyCone, } Base.copy(cone::HypatiaCones) = cone # NOTE maybe should deep copy the cone struct, but this is expensive diff --git a/test/barrier.jl b/test/barrier.jl deleted file mode 100644 index 3170df927..000000000 --- a/test/barrier.jl +++ /dev/null @@ -1,599 +0,0 @@ -#= -tests for primitive cone barrier oracles -=# - -using Test -import Random -import GenericLinearAlgebra.eigen -using LinearAlgebra -using SparseArrays -# import ForwardDiff # TODO not using barrier functions -import Hypatia -import Hypatia.ModelUtilities -import Hypatia.Cones - -function test_barrier_oracles( - cone::Cones.Cone{T}, - barrier::Function; - noise::T = T(0.1), - scale::T = T(1e-3), - tol::Real = 1000eps(T), - init_tol::Real = tol, - init_only::Bool = false, - ) where {T <: Real} - Random.seed!(1) - - Cones.setup_data(cone) - dim = Cones.dimension(cone) - point = zeros(T, dim) - dual_point = copy(point) - Cones.set_initial_point(point, cone) - Cones.set_initial_point(dual_point, cone) - load_reset_check(cone, point, dual_point) - @test cone.point == point - @test cone.dual_point == dual_point - - if isfinite(init_tol) - # tests for centrality of initial point - minus_grad = -Cones.grad(cone) - @test dot(point, minus_grad) ≈ norm(point) * norm(minus_grad) atol=init_tol rtol=init_tol - @test point ≈ minus_grad atol=init_tol rtol=init_tol - @test Cones.in_neighborhood(cone, one(T), one(T)) - end - init_only && return - - # perturb and scale the initial point and check feasible - perturb_scale(point, dual_point, noise, one(T)) - load_reset_check(cone, point, dual_point) - - # test gradient and Hessian oracles - test_grad_hess(cone, point, dual_point, tol = tol) - - # check gradient and Hessian agree with ForwardDiff - Cones.reset_data(cone) - @test Cones.is_feas(cone) - @test Cones.is_dual_feas(cone) - grad = Cones.grad(cone) - hess = Cones.hess(cone) - # if dim < 8 - # grad = Cones.grad(cone) - # fd_grad = ForwardDiff.gradient(barrier, point) - # @test grad ≈ fd_grad atol=tol rtol=tol - # hess = Cones.hess(cone) - # fd_hess = ForwardDiff.hessian(barrier, point) - # @test hess ≈ fd_hess atol=tol rtol=tol - # end - - if Cones.use_correction(cone) - # check correction satisfies log-homog property F'''(s)[s, s] = -2F''(s) * s = 2F'(s) - @test -Cones.correction(cone, point) ≈ grad atol=tol rtol=tol - # check correction term agrees with directional 3rd derivative - (primal_dir, dual_dir) = perturb_scale(zeros(T, dim), zeros(T, dim), noise, one(T)) - corr = Cones.correction(cone, primal_dir) - @test dot(corr, point) ≈ dot(primal_dir, hess * primal_dir) atol=tol rtol=tol - - # barrier_dir(point, t) = barrier(point + t * primal_dir) - # @test -2 * corr ≈ ForwardDiff.gradient(x -> ForwardDiff.derivative(s -> ForwardDiff.derivative(t -> barrier_dir(x, t), s), 0), point) atol=tol rtol=tol - end - - return -end - -function test_grad_hess(cone::Cones.Cone{T}, point::Vector{T}, dual_point::Vector{T}; tol::Real = 1000eps(T)) where {T <: Real} - # TODO not currently using dual_point - nu = Cones.get_nu(cone) - dim = length(point) - grad = Cones.grad(cone) - hess = Matrix(Cones.hess(cone)) - inv_hess = Matrix(Cones.inv_hess(cone)) - - @test dot(point, grad) ≈ -nu atol=tol rtol=tol - @test hess * inv_hess ≈ I atol=tol rtol=tol - - prod_mat = zeros(T, dim, dim) - @test Cones.hess_prod!(prod_mat, inv_hess, cone) ≈ I atol=tol rtol=tol - @test Cones.inv_hess_prod!(prod_mat, hess, cone) ≈ I atol=tol rtol=tol - - prod = zero(point) - @test hess * point ≈ -grad atol=tol rtol=tol - @test Cones.hess_prod!(prod, point, cone) ≈ -grad atol=tol rtol=tol - @test Cones.inv_hess_prod!(prod, grad, cone) ≈ -point atol=tol rtol=tol - - if Cones.use_sqrt_oracles(cone) - prod_mat2 = Matrix(Cones.hess_sqrt_prod!(prod_mat, inv_hess, cone)') - @test Cones.hess_sqrt_prod!(prod_mat, prod_mat2, cone) ≈ I atol=tol rtol=tol - Cones.inv_hess_sqrt_prod!(prod_mat2, Matrix(one(T) * I, dim, dim), cone) - @test prod_mat2' * prod_mat2 ≈ inv_hess atol=tol rtol=tol - end - - return -end - -function load_reset_check(cone::Cones.Cone{T}, point::Vector{T}, dual_point::Vector{T}) where {T <: Real} - Cones.load_point(cone, point) - Cones.load_dual_point(cone, dual_point) - Cones.reset_data(cone) - @test Cones.is_feas(cone) - @test Cones.is_dual_feas(cone) - return -end - -function perturb_scale(point::Vector{T}, dual_point::Vector{T}, noise::T, scale::T) where {T <: Real} - if !iszero(noise) - @. point += 2 * noise * rand(T) - noise - @. dual_point += 2 * noise * rand(T) - noise - end - if !isone(scale) - point .*= scale - end - return (point, dual_point) -end - -# primitive cone barrier tests - -function test_nonnegative_barrier(T::Type{<:Real}) - barrier = (s -> -sum(log, s)) - for dim in [1, 2, 3, 6] - test_barrier_oracles(Cones.Nonnegative{T}(dim), barrier) - end - return -end - -function test_epinorminf_barrier(T::Type{<:Real}) - for n in [1, 2, 3, 5] - # real epinorminf cone - function R_barrier(s) - (u, w) = (s[1], s[2:end]) - return -sum(log(abs2(u) - abs2(wj)) for wj in w) + (n - 1) * log(u) - end - test_barrier_oracles(Cones.EpiNormInf{T, T}(1 + n), R_barrier) - - # complex epinorminf cone - function C_barrier(s) - (u, wr) = (s[1], s[2:end]) - w = Cones.rvec_to_cvec!(zeros(Complex{eltype(s)}, n), wr) - return -sum(log(abs2(u) - abs2(wj)) for wj in w) + (n - 1) * log(u) - end - test_barrier_oracles(Cones.EpiNormInf{T, Complex{T}}(1 + 2n), C_barrier) - end - return -end - -function test_epinormeucl_barrier(T::Type{<:Real}) - function barrier(s) - (u, w) = (s[1], s[2:end]) - return -log(abs2(u) - sum(abs2, w)) - end - for dim in [2, 3, 4, 6] - test_barrier_oracles(Cones.EpiNormEucl{T}(dim), barrier) - end - return -end - -function test_epipersquare_barrier(T::Type{<:Real}) - function barrier(s) - (u, v, w) = (s[1], s[2], s[3:end]) - return -log(2 * u * v - sum(abs2, w)) - end - for dim in [3, 4, 5, 7] - test_barrier_oracles(Cones.EpiPerSquare{T}(dim), barrier) - end - return -end - -function test_hypoperlog_barrier(T::Type{<:Real}) - function barrier(s) - (u, v, w) = (s[1], s[2], s[3:end]) - return -log(v * sum(log(wj / v) for wj in w) - u) - sum(log, w) - log(v) - end - for dim in [3, 5, 10] - test_barrier_oracles(Cones.HypoPerLog{T}(dim), barrier, init_tol = 1e-5) - end - for dim in [15, 65, 75, 100] - test_barrier_oracles(Cones.HypoPerLog{T}(dim), barrier, init_tol = 1e-1, init_only = true) - end - return -end - -function test_epiperentropy_barrier(T::Type{<:Real}) - function barrier(s) - (u, v, w) = (s[1], s[2], s[3:end]) - return -log(u - sum(wi * log(wi / v) for wi in w)) - log(v) - sum(log(wi) for wi in w) - end - for w_dim in [1, 2, 3] - test_barrier_oracles(Cones.EpiPerEntropy{T}(2 + w_dim), barrier, init_tol = 1e-5) - end - for w_dim in [15, 65, 75, 100] - test_barrier_oracles(Cones.EpiPerEntropy{T}(2 + w_dim), barrier, init_tol = 1e-1, init_only = true) - end - return -end - -function test_epipertraceentropytri_barrier(T::Type{<:Real}) - for side in [1, 2, 3, 12, 20] - dim = 2 + Cones.svec_length(side) - function barrier(s) - (u, v, w) = (s[1], s[2], s[3:end]) - W = Hermitian(Cones.svec_to_smat!(zeros(eltype(s), side, side), w, sqrt(T(2))), :U) - return -log(u - dot(W, log(W / v))) - log(v) - logdet(W) - end - if side <= 3 - test_barrier_oracles(Cones.EpiPerTraceEntropyTri{T}(dim), barrier, init_tol = 1e-5) - else - test_barrier_oracles(Cones.EpiPerTraceEntropyTri{T}(dim), barrier, init_tol = 1e-1, init_only = true) - end - end - return -end - -function test_epirelentropy_barrier(T::Type{<:Real}) - function barrier(s) - (u, v, w) = (s[1], s[2:2:(end - 1)], s[3:2:end]) - return -log(u - sum(wi * log(wi / vi) for (vi, wi) in zip(v, w))) - sum(log(vi) + log(wi) for (vi, wi) in zip(v, w)) - end - for w_dim in [1, 2, 3] - test_barrier_oracles(Cones.EpiRelEntropy{T}(1 + 2 * w_dim), barrier, init_tol = 1e-5) - end - for w_dim in [15, 65, 75, 100] - test_barrier_oracles(Cones.EpiRelEntropy{T}(1 + 2 * w_dim), barrier, init_tol = 1e-1, init_only = true) - end - return -end - -function test_hypogeomean_barrier(T::Type{<:Real}) - for dim in [2, 3, 5, 8] - invn = inv(T(dim - 1)) - function barrier(s) - (u, w) = (s[1], s[2:end]) - return -log(exp(sum(invn * log(w[j]) for j in eachindex(w))) - u) - sum(log, w) - end - test_barrier_oracles(Cones.HypoGeoMean{T}(dim), barrier) - end - return -end - -function test_hypopowermean_barrier(T::Type{<:Real}) - Random.seed!(1) - for dim in [2, 3, 5, 15, 90, 120] - alpha = rand(T, dim - 1) .+ 1 - alpha ./= sum(alpha) - function barrier(s) - (u, w) = (s[1], s[2:end]) - return -log(exp(sum(alpha[j] * log(w[j]) for j in eachindex(w))) - u) - sum(log, w) - end - cone = Cones.HypoPowerMean{T}(alpha) - test_barrier_oracles(cone, barrier, init_tol = 1e-2, init_only = (dim > 5)) - # test initial point when all alphas are the same - cone = Cones.HypoPowerMean{T}(fill(inv(T(dim - 1)), dim - 1)) - test_barrier_oracles(cone, barrier, init_tol = sqrt(eps(T)), init_only = true) - end - return -end - -function test_power_barrier(T::Type{<:Real}) - Random.seed!(1) - for (m, n) in [(2, 1), (2, 2), (4, 1), (2, 4)] - alpha = rand(T, m) .+ 1 - alpha ./= sum(alpha) - function barrier(s) - (u, w) = (s[1:m], s[(m + 1):end]) - return -log(exp(2 * sum(alpha[j] * log(u[j]) for j in eachindex(u))) - sum(abs2, w)) - sum((1 - alpha[j]) * log(u[j]) for j in eachindex(u)) - end - test_barrier_oracles(Cones.Power{T}(alpha, n), barrier) - end - return -end - -function test_epinormspectral_barrier(T::Type{<:Real}) - for (n, m) in [(1, 1), (1, 2), (2, 2), (2, 4), (3, 4)] - # real epinormspectral barrier - function R_barrier(s) - (u, W) = (s[1], reshape(s[2:end], n, m)) - return -logdet(cholesky!(Hermitian(abs2(u) * I - W * W'))) + (n - 1) * log(u) - end - test_barrier_oracles(Cones.EpiNormSpectral{T, T}(n, m), R_barrier) - - # complex epinormspectral barrier - function C_barrier(s) - u = s[1] - W = Cones.rvec_to_cvec!(zeros(Complex{eltype(s)}, n, m), s[2:end]) - return -logdet(cholesky!(Hermitian(abs2(u) * I - W * W'))) + (n - 1) * log(u) - end - test_barrier_oracles(Cones.EpiNormSpectral{T, Complex{T}}(n, m), C_barrier) - end - return -end - -function test_epitracerelentropytri_barrier(T::Type{<:Real}) - rt2 = sqrt(T(2)) - for side in [1, 2, 3, 8, 12] - svec_dim = Cones.svec_length(side) - cone = Cones.EpiTraceRelEntropyTri{T}(2 * svec_dim + 1) - function barrier(s) - u = s[1] - u = s[1] - V = Hermitian(Cones.svec_to_smat!(similar(s, side, side), s[2:(svec_dim + 1)], rt2), :U) - W = Hermitian(Cones.svec_to_smat!(similar(s, side, side), s[(svec_dim + 2):end], rt2), :U) - return -log(u - tr(W * log(W) - W * log(V))) - logdet(V) - logdet(W) - end - if side <= 3 - test_barrier_oracles(cone, barrier, init_tol = 1e-5) - else - test_barrier_oracles(cone, barrier, init_tol = 1e-1, init_only = true) - end - end - return -end - -function test_linmatrixineq_barrier(T::Type{<:Real}) - if !(T <: LinearAlgebra.BlasReal) - return # TODO currently failing with BigFloat due to an apparent cholesky bug - end - Random.seed!(1) - Rs_list = [[T, T], [Complex{T}, Complex{T}], [T, Complex{T}, T], [Complex{T}, T, T]] - for side in [2, 3, 5], Rs in Rs_list - As = Vector{LinearAlgebra.HermOrSym{R, Matrix{R}} where {R <: Hypatia.RealOrComplex{T}}}(undef, length(Rs)) - A_1_half = rand(Rs[1], side, side) - As[1] = Hermitian(A_1_half * A_1_half' + I, :U) - for i in 2:length(Rs) - As[i] = Hermitian(rand(Rs[i], side, side), :U) - end - barrier(s) = -logdet(cholesky!(Hermitian(sum(s_i * A_i for (s_i, A_i) in zip(s, As)), :U))) - test_barrier_oracles(Cones.LinMatrixIneq{T}(As), barrier, noise = T(1e-2), init_tol = Inf) - end - return -end - -function test_possemideftri_barrier(T::Type{<:Real}) - for side in [1, 2, 3, 5] - # real PSD cone - function R_barrier(s) - S = zeros(eltype(s), side, side) - Cones.svec_to_smat!(S, s, sqrt(T(2))) - return -logdet(cholesky!(Symmetric(S, :U))) - end - dim = Cones.svec_length(side) - test_barrier_oracles(Cones.PosSemidefTri{T, T}(dim), R_barrier) - - # complex PSD cone - function C_barrier(s) - S = zeros(Complex{eltype(s)}, side, side) - Cones.svec_to_smat!(S, s, sqrt(T(2))) - return -logdet(cholesky!(Hermitian(S, :U))) - end - dim = side^2 - test_barrier_oracles(Cones.PosSemidefTri{T, Complex{T}}(dim), C_barrier) - end - return -end - -function test_possemideftrisparse_barrier(T::Type{<:Real}) - if !(T <: LinearAlgebra.BlasReal) - return # only works with BLAS real types - end - Random.seed!(1) - invrt2 = inv(sqrt(T(2))) - - for side in [1, 2, 5, 10, 25, 50] - # generate random sparsity pattern for lower triangle - sparsity = inv(sqrt(side)) - (row_idxs, col_idxs, _) = findnz(tril!(sprand(Bool, side, side, sparsity)) + I) - - # real sparse PSD cone - function R_barrier(s) - scal_s = copy(s) - for i in eachindex(s) - if row_idxs[i] != col_idxs[i] - scal_s[i] *= invrt2 - end - end - S = Matrix(sparse(row_idxs, col_idxs, scal_s, side, side)) - return -logdet(cholesky(Symmetric(S, :L))) - end - test_barrier_oracles(Cones.PosSemidefTriSparse{T, T}(side, row_idxs, col_idxs), R_barrier) - - # complex sparse PSD cone - function C_barrier(s) - scal_s = zeros(Complex{eltype(s)}, length(row_idxs)) - idx = 1 - for i in eachindex(scal_s) - if row_idxs[i] == col_idxs[i] - scal_s[i] = s[idx] - idx += 1 - else - scal_s[i] = invrt2 * Complex(s[idx], s[idx + 1]) - idx += 2 - end - end - S = Matrix(sparse(row_idxs, col_idxs, scal_s, side, side)) - return -logdet(cholesky!(Hermitian(S, :L))) - end - test_barrier_oracles(Cones.PosSemidefTriSparse{T, Complex{T}}(side, row_idxs, col_idxs), C_barrier) - end - return -end - -function test_doublynonnegativetri_barrier(T::Type{<:Real}) - for side in [1, 2, 3, 5] - function barrier(s) - S = zeros(eltype(s), side, side) - Cones.svec_to_smat!(S, s, sqrt(T(2))) - offdiags = vcat([div(i * (i - 1), 2) .+ (1:(i - 1)) for i in 2:side]...) - return -logdet(cholesky!(Hermitian(S, :U))) - mapreduce(log, +, s[offdiags]; init = zero(eltype(s))) - end - dim = Cones.svec_length(side) - test_barrier_oracles(Cones.DoublyNonnegativeTri{T}(dim), barrier, init_tol = sqrt(eps(T)), init_only = (side > 6)) - end - return -end - -function test_matrixepipersquare_barrier(T::Type{<:Real}) - for (n, m) in [(1, 1), (1, 2), (2, 2), (2, 4), (3, 4)] - # real matrixepipersquare barrier - per_idx = Cones.svec_length(n) + 1 - function R_barrier(s) - U = Cones.svec_to_smat!(zeros(eltype(s), n, n), s[1:(per_idx - 1)], sqrt(T(2))) - v = s[per_idx] - W = reshape(s[(per_idx + 1):end], n, m) - return -logdet(cholesky!(Symmetric(2 * v * U - W * W', :U))) + (n - 1) * log(v) - end - test_barrier_oracles(Cones.MatrixEpiPerSquare{T, T}(n, m), R_barrier) - - # complex matrixepipersquare barrier - per_idx = n ^ 2 + 1 - function C_barrier(s) - U = Cones.svec_to_smat!(zeros(Complex{eltype(s)}, n, n), s[1:(per_idx - 1)], sqrt(T(2))) - v = s[per_idx] - W = Cones.rvec_to_cvec!(zeros(Complex{eltype(s)}, n, m), s[(per_idx + 1):end]) - return -logdet(cholesky!(Hermitian(2 * v * U - W * W', :U))) + (n - 1) * log(v) - end - test_barrier_oracles(Cones.MatrixEpiPerSquare{T, Complex{T}}(n, m), C_barrier) - end - return -end - -function test_hypoperlogdettri_barrier(T::Type{<:Real}) - for side in [1, 2, 3, 5, 8] - # real logdet barrier - dim = 2 + Cones.svec_length(side) - cone = Cones.HypoPerLogdetTri{T, T}(dim) - function R_barrier(s) - (u, v, W) = (s[1], s[2], zeros(eltype(s), side, side)) - Cones.svec_to_smat!(W, s[3:end], sqrt(T(2))) - return -log(v * logdet(cholesky!(Symmetric(W / v, :U))) - u) - logdet(cholesky!(Symmetric(W, :U))) - log(v) - end - if side <= 5 - test_barrier_oracles(cone, R_barrier, init_tol = 1e-5) - else - test_barrier_oracles(cone, R_barrier, init_tol = 1e-1, init_only = true) - end - - # complex logdet barrier - dim = 2 + side^2 - cone = Cones.HypoPerLogdetTri{T, Complex{T}}(dim) - function C_barrier(s) - (u, v, W) = (s[1], s[2], zeros(Complex{eltype(s)}, side, side)) - Cones.svec_to_smat!(W, s[3:end], sqrt(T(2))) - return -log(v * logdet(cholesky!(Hermitian(W / v, :U))) - u) - logdet(cholesky!(Hermitian(W, :U))) - log(v) - end - if side <= 5 - test_barrier_oracles(cone, C_barrier, init_tol = 1e-5) - else - test_barrier_oracles(cone, C_barrier, init_tol = 1e-1, init_only = true) - end - end - return -end - -function test_hyporootdettri_barrier(T::Type{<:Real}) - for side in [1, 2, 3, 5, 8] - # real rootdet barrier - dim = 1 + Cones.svec_length(side) - cone = Cones.HypoRootdetTri{T, T}(dim) - function R_barrier(s) - (u, W) = (s[1], zeros(eltype(s), side, side)) - Cones.svec_to_smat!(W, s[2:end], sqrt(T(2))) - fact_W = cholesky!(Symmetric(W, :U)) - return -log(exp(logdet(fact_W) / side) - u) - logdet(fact_W) - end - test_barrier_oracles(cone, R_barrier) - - # complex rootdet barrier - dim = 1 + side^2 - cone = Cones.HypoRootdetTri{T, Complex{T}}(dim) - function C_barrier(s) - (u, W) = (s[1], zeros(Complex{eltype(s)}, side, side)) - Cones.svec_to_smat!(W, s[2:end], sqrt(T(2))) - fact_W = cholesky!(Hermitian(W, :U)) - return -log(exp(logdet(fact_W) / side) - u) - logdet(fact_W) - end - test_barrier_oracles(cone, C_barrier) - end - return -end - -function test_wsosinterpnonnegative_barrier(T::Type{<:Real}) - Random.seed!(1) - for (n, halfdeg) in [(1, 1), (1, 2), (1, 3), (2, 1), (2, 2), (3, 1)] - (U, _, Ps, _) = ModelUtilities.interpolate(ModelUtilities.Box{T}(-ones(T, n), ones(T, n)), halfdeg, sample = false) # use a unit box domain - barrier(s) = -sum(logdet(cholesky!(Symmetric(P' * Diagonal(s) * P))) for P in Ps) - cone = Cones.WSOSInterpNonnegative{T, T}(U, Ps) - test_barrier_oracles(cone, barrier, init_tol = Inf) # TODO center and test initial points - end - # TODO also test complex case Cones.WSOSInterpNonnegative{T, Complex{T}} - need complex MU interp functions first - return -end - -function test_wsosinterppossemideftri_barrier(T::Type{<:Real}) - Random.seed!(1) - rt2i = inv(sqrt(T(2))) - for (n, halfdeg, R) in [(1, 1, 1), (1, 1, 4), (2, 2, 1), (2, 1, 3), (3, 1, 2)] - (U, _, Ps, _) = ModelUtilities.interpolate(ModelUtilities.Box{T}(-ones(T, n), ones(T, n)), halfdeg, sample = false) - cone = Cones.WSOSInterpPosSemidefTri{T}(R, U, Ps) - function barrier(s) - bar = zero(eltype(s)) - for P in Ps - L = size(P, 2) - Lambda = zeros(eltype(s), R * L, R * L) - for i in 1:R, j in 1:i - Lambdaij = P' * Diagonal(s[Cones.block_idxs(U, Cones.svec_idx(i, j))]) * P - if i != j - Lambdaij .*= rt2i - end - Lambda[Cones.block_idxs(L, i), Cones.block_idxs(L, j)] = Lambdaij - end - bar -= logdet(cholesky!(Symmetric(Lambda, :L))) - end - return bar - end - test_barrier_oracles(cone, barrier, init_tol = Inf) - end - return -end - -function test_wsosinterpepinormeucl_barrier(T::Type{<:Real}) - Random.seed!(1) - for (n, halfdeg, R) in [(1, 1, 2), (1, 2, 4), (2, 2, 3), (3, 1, 2)] - (U, _, Ps, _) = ModelUtilities.interpolate(ModelUtilities.Box{T}(-ones(T, n), ones(T, n)), halfdeg, sample = false) - cone = Cones.WSOSInterpEpiNormEucl{T}(R, U, Ps) - function barrier(s) - bar = zero(eltype(s)) - for P in Ps - Lambda = P' * Diagonal(s[1:U]) * P - Lambda1fact = cholesky!(Symmetric(copy(Lambda), :L)) - for i in 2:R - Lambdai = P' * Diagonal(s[Cones.block_idxs(U, i)]) * P - ldiv!(Lambda1fact.L, Lambdai) - Lambda -= Lambdai' * Lambdai - end - bar -= logdet(cholesky!(Symmetric(Lambda))) + logdet(Lambda1fact) - end - return bar - end - test_barrier_oracles(cone, barrier, init_tol = Inf) - end - return -end - -function test_wsosinterpepinormone_barrier(T::Type{<:Real}) - Random.seed!(1) - for n in 1:3, halfdeg in 1:2, R in 2:3 - (U, _, Ps, _) = ModelUtilities.interpolate(ModelUtilities.Box{T}(-ones(T, n), ones(T, n)), halfdeg, sample = false) - cone = Cones.WSOSInterpEpiNormOne{T}(R, U, Ps) - function barrier(point) - bar = zero(eltype(point)) - for P in cone.Ps - lambda_1 = Symmetric(P' * Diagonal(point[1:U]) * P, :U) - fact_1 = cholesky(lambda_1) - for i in 2:R - lambda_i = Symmetric(P' * Diagonal(point[Cones.block_idxs(U, i)]) * P) - bar -= logdet(lambda_1 - lambda_i * (fact_1 \ lambda_i)) - end - bar -= logdet(fact_1) - end - return bar - end - test_barrier_oracles(cone, barrier, init_tol = Inf) - end - return -end diff --git a/test/cone.jl b/test/cone.jl new file mode 100644 index 000000000..58556552d --- /dev/null +++ b/test/cone.jl @@ -0,0 +1,806 @@ +#= +tests for primitive cone barrier oracles +=# + +using Test +import Random +import Random.randn +import DataStructures +import GenericLinearAlgebra.eigen +using LinearAlgebra +using SparseArrays +import ForwardDiff +import Hypatia +import Hypatia.ModelUtilities +import Hypatia.Cones +import Hypatia.RealOrComplex + +Random.randn(::Type{BigFloat}, dims::Integer...) = BigFloat.(randn(dims...)) +Random.randn(::Type{Complex{BigFloat}}, dims::Integer...) = Complex{BigFloat}.(randn(ComplexF64, dims...)) + +# sanity check oracles +function test_oracles( + cone::Cones.Cone{T}; + noise::T = T(1e-1), + scale::T = T(1e-2), + tol::Real = 1e3 * eps(T), + init_only::Bool = false, + init_tol::Real = tol, + ) where {T <: Real} + Random.seed!(1) + dim = Cones.dimension(cone) + Cones.setup_data(cone) + Cones.reset_data(cone) + + point = zeros(T, dim) + Cones.set_initial_point(point, cone) + Cones.load_point(cone, point) + @test Cones.is_feas(cone) + @test cone.point == point + + dual_point = -Cones.grad(cone) + Cones.load_dual_point(cone, dual_point) + @test Cones.is_dual_feas(cone) + @test cone.dual_point == dual_point + @test Cones.in_neighborhood(cone, one(T), one(T)) + + # test centrality of initial point + if isfinite(init_tol) + @test point ≈ dual_point atol=init_tol rtol=init_tol + end + init_only && return + + # perturb and scale the initial point + perturb_scale(point, noise, scale) + perturb_scale(dual_point, noise, inv(scale)) + + Cones.reset_data(cone) + Cones.load_point(cone, point) + @test Cones.is_feas(cone) + Cones.load_dual_point(cone, dual_point) + @test Cones.is_dual_feas(cone) + + # test gradient and Hessian oracles + nu = Cones.get_nu(cone) + grad = Cones.grad(cone) + @test dot(point, grad) ≈ -nu atol=tol rtol=tol + + hess = Matrix(Cones.hess(cone)) + inv_hess = Matrix(Cones.inv_hess(cone)) + @test hess * inv_hess ≈ I atol=tol rtol=tol + + prod_vec = zero(point) + @test hess * point ≈ -grad atol=tol rtol=tol + @test Cones.hess_prod!(prod_vec, point, cone) ≈ -grad atol=tol rtol=tol + @test Cones.inv_hess_prod!(prod_vec, grad, cone) ≈ -point atol=tol rtol=tol + + prod_mat = zeros(T, dim, dim) + @test Cones.hess_prod!(prod_mat, inv_hess, cone) ≈ I atol=tol rtol=tol + @test Cones.inv_hess_prod!(prod_mat, hess, cone) ≈ I atol=tol rtol=tol + + if Cones.use_sqrt_oracles(cone) + prod_mat2 = Matrix(Cones.hess_sqrt_prod!(prod_mat, inv_hess, cone)') + @test Cones.hess_sqrt_prod!(prod_mat, prod_mat2, cone) ≈ I atol=tol rtol=tol + Cones.inv_hess_sqrt_prod!(prod_mat2, Matrix(one(T) * I, dim, dim), cone) + @test prod_mat2' * prod_mat2 ≈ inv_hess atol=tol rtol=tol + end + + # test correction oracle + if Cones.use_correction(cone) + @test -Cones.correction(cone, point) ≈ grad atol=tol rtol=tol + + dir = perturb_scale(zeros(T, dim), noise, one(T)) + corr = Cones.correction(cone, dir) + @test dot(corr, point) ≈ dot(dir, hess * dir) atol=tol rtol=tol + end + + return +end + +# check some oracles agree with ForwardDiff +function test_barrier( + cone::Cones.Cone{T}, + barrier::Function; + noise::T = T(1e-1), + scale::T = T(1e-2), + tol::Real = 1e4 * eps(T), + ) where {T <: Real} + Random.seed!(1) + dim = Cones.dimension(cone) + Cones.setup_data(cone) + + point = zeros(T, dim) + Cones.set_initial_point(point, cone) + perturb_scale(point, noise, scale) + + Cones.reset_data(cone) + Cones.load_point(cone, point) + @test Cones.is_feas(cone) + + fd_grad = ForwardDiff.gradient(barrier, point) + @test Cones.grad(cone) ≈ fd_grad atol=tol rtol=tol + + dir = 10 * randn(T, dim) + barrier_dir(s, t) = barrier(s + t * dir) + + fd_hess_dir = ForwardDiff.gradient(s -> ForwardDiff.derivative(t -> barrier_dir(s, t), 0), point) + @test Cones.hess(cone) * dir ≈ fd_hess_dir atol=tol rtol=tol + + if Cones.use_correction(cone) + fd_third_dir = ForwardDiff.gradient(s2 -> ForwardDiff.derivative(s -> ForwardDiff.derivative(t -> barrier_dir(s2, t), s), 0), point) + @test -2 * Cones.correction(cone, dir) ≈ fd_third_dir atol=tol rtol=tol + end + + return +end + +# get memory allocations for oracles +function get_allocs( + cone::Cones.Cone{T}; + noise::T = T(1e-1), + scale::T = T(1e-2), + ) where {T <: Real} + Random.seed!(1) + allocs = DataStructures.OrderedDict{Symbol, Int}() + dim = Cones.dimension(cone) + allocs[:dimension] = dim + + allocs[:setup_data] = @allocated Cones.setup_data(cone) + Cones.reset_data(cone) + + point = zeros(T, dim) + Cones.set_initial_point(point, cone) + Cones.load_point(cone, point) + allocs[:is_feas] = @allocated Cones.is_feas(cone) + + dual_point = -Cones.grad(cone) + Cones.load_dual_point(cone, dual_point) + allocs[:is_dual_feas] = @allocated Cones.is_dual_feas(cone) + + allocs[:grad] = @allocated Cones.grad(cone) + allocs[:hess] = @allocated Cones.hess(cone) + allocs[:inv_hess] = @allocated Cones.inv_hess(cone) + + point1 = randn(T, dim) + point2 = zero(point1) + allocs[:hess_prod] = @allocated Cones.hess_prod!(point2, point1, cone) + allocs[:inv_hess_prod] = @allocated Cones.inv_hess_prod!(point2, point1, cone) + + if Cones.use_sqrt_oracles(cone) + allocs[:hess_sqrt_prod] = @allocated Cones.hess_sqrt_prod!(point2, point1, cone) + allocs[:inv_hess_sqrt_prod] = @allocated Cones.inv_hess_sqrt_prod!(point2, point1, cone) + end + + if Cones.use_correction(cone) + allocs[:correction] = @allocated Cones.correction(cone, point1) + end + + allocs[:in_neighborhood] = @allocated Cones.in_neighborhood(cone, one(T), one(T)) + + return allocs +end + +function perturb_scale( + point::Vector{T}, + noise::T, + scale::T, + ) where {T <: Real} + if !iszero(noise) + @. point += 2 * noise * rand(T) - noise + end + if !isone(scale) + point .*= scale + end + return point +end + +# cone utilities + +logdet_pd(W::Hermitian) = logdet(cholesky!(copy(W))) + +# TODO maybe move to ModelUtilities + +dim_vec(d::Int, ::Type{<:Real}) = d +dim_vec(d::Int, ::Type{<:Complex}) = 2 * d + +dim_herm(d::Int, ::Type{<:Real}) = Cones.svec_length(d) +dim_herm(d::Int, ::Type{<:Complex}) = d^2 + +function new_vec(w::Vector, dw::Int, T::Type{<:Real}) + @assert length(w) == dw + return w +end +function new_vec(w::Vector, dw::Int, R::Type{Complex{T}}) where {T <: Real} + @assert length(w) == 2 * dw + wR = zeros(Complex{eltype(w)}, dw) + Cones.rvec_to_cvec!(wR, w) + return wR +end + +function new_mat_herm(w::Vector, dW::Int, T::Type{<:Real}) + @assert length(w) == dim_herm(dW, T) + W = similar(w, dW, dW) + Cones.svec_to_smat!(W, w, sqrt(T(2))) + return Hermitian(W, :U) +end +function new_mat_herm(w::Vector, dW::Int, R::Type{Complex{T}}) where {T <: Real} + @assert length(w) == dim_herm(dW, R) + W = zeros(Complex{eltype(w)}, dW, dW) + Cones.svec_to_smat!(W, w, sqrt(T(2))) + return Hermitian(W, :U) +end + +function rand_sppsd_pattern(dW::Int) + sparsity = inv(sqrt(dW)) + (row_idxs, col_idxs, _) = findnz(tril!(sprand(Bool, dW, dW, sparsity)) + I) + return (row_idxs, col_idxs) +end + +function rand_herms(ds::Int, Rd::Vector, T::Type{<:Real}) + Ps = Vector{LinearAlgebra.HermOrSym{R, Matrix{R}} where {R <: RealOrComplex{T}}}(undef, length(Rd)) + A_1_half = randn(Rd[1], ds, ds) + Ps[1] = Hermitian(A_1_half * A_1_half' + I, :U) + for i in 2:length(Rd) + Ps[i] = Hermitian(randn(Rd[i], ds, ds), :U) + end + return Ps +end + +function rand_powers(T, d) + Random.seed!(1) + alpha = rand(T, d) .+ 1 + alpha ./= sum(alpha) + return alpha +end + +function rand_interp(num_vars::Int, halfdeg::Int, T::Type{<:Real}) + Random.seed!(1) + # use a unit box domain + domain = ModelUtilities.Box{T}(-ones(T, num_vars), ones(T, num_vars)) + (d, _, Ps, _) = ModelUtilities.interpolate(domain, halfdeg, sample = false) + return (d, Ps) +end + + +# cones + +# Nonnegative +function test_oracles(C::Type{<:Cones.Nonnegative}) + for d in [1, 2, 6] + test_oracles(C(d)) + end +end + +function test_barrier(C::Type{<:Cones.Nonnegative}) + barrier = (s -> -sum(log, s)) + test_barrier(C(3), barrier) +end + +get_allocs(C::Type{<:Cones.Nonnegative}) = get_allocs(C(9)) + + +# PosSemidefTri +function test_oracles(C::Type{Cones.PosSemidefTri{T, R}}) where {T, R} + for dW in [1, 2, 3, 5] + test_oracles(C(dim_herm(dW, R))) + end +end + +function test_barrier(C::Type{Cones.PosSemidefTri{T, R}}) where {T, R} + dW = 3 + barrier(s) = -logdet_pd(Hermitian(new_mat_herm(s, dW, R), :U)) + test_barrier(C(dim_herm(dW, R)), barrier) +end + +get_allocs(C::Type{Cones.PosSemidefTri{T, R}}) where {T, R} = get_allocs(C(dim_herm(4, R))) + + +# DoublyNonnegativeTri +function test_oracles(C::Type{Cones.DoublyNonnegativeTri{T}}) where T + for dW in [1, 2, 5] + test_oracles(C(Cones.svec_length(dW)), init_tol = sqrt(eps(T))) + end + for dW in [10, 20] + test_oracles(C(Cones.svec_length(dW)), init_tol = sqrt(eps(T)), init_only = true) + end +end + +function test_barrier(C::Type{Cones.DoublyNonnegativeTri{T}}) where T + dW = 3 + function barrier(s) + W = new_mat_herm(s, dW, T) + offdiags = vcat([div(i * (i - 1), 2) .+ (1:(i - 1)) for i in 2:dW]...) + return -logdet_pd(Hermitian(W, :U)) - sum(log, s[offdiags]) + end + test_barrier(C(Cones.svec_length(dW)), barrier) +end + +get_allocs(C::Type{<:Cones.DoublyNonnegativeTri}) = get_allocs(C(10)) + + +# PosSemidefTriSparse +function test_oracles(C::Type{<:Cones.PosSemidefTriSparse}) + for dW in [1, 2, 10, 25, 40] + test_oracles(C(dW, rand_sppsd_pattern(dW)...)) + end +end + +function test_barrier(C::Type{Cones.PosSemidefTriSparse{T, T}}) where T + dW = 25 + (row_idxs, col_idxs) = rand_sppsd_pattern(dW) + invrt2 = inv(sqrt(T(2))) + function barrier(s) + scal_s = copy(s) + scal_s[row_idxs .!= col_idxs] .*= invrt2 + W = Matrix(sparse(row_idxs, col_idxs, scal_s, dW, dW)) + return -logdet_pd(Hermitian(W, :L)) + end + test_barrier(C(dW, row_idxs, col_idxs), barrier) +end + +function test_barrier(C::Type{Cones.PosSemidefTriSparse{T, Complex{T}}}) where T + dW = 25 + (row_idxs, col_idxs) = rand_sppsd_pattern(dW) + invrt2 = inv(sqrt(T(2))) + function barrier(s) + scal_s = zeros(Complex{eltype(s)}, length(row_idxs)) + idx = 1 + for i in eachindex(scal_s) + if row_idxs[i] == col_idxs[i] + scal_s[i] = s[idx] + idx += 1 + else + scal_s[i] = invrt2 * Complex(s[idx], s[idx + 1]) + idx += 2 + end + end + W = Matrix(sparse(row_idxs, col_idxs, scal_s, dW, dW)) + return -logdet_pd(Hermitian(W, :L)) + end + test_barrier(C(dW, row_idxs, col_idxs), barrier) +end + +get_allocs(C::Type{<:Cones.PosSemidefTriSparse}) = get_allocs(C(15, rand_sppsd_pattern(15)...)) + + +# LinMatrixIneq +function test_oracles(C::Type{Cones.LinMatrixIneq{T}}) where T + Random.seed!(1) + Rd_list = [[T, T], [T, Complex{T}], [Complex{T}, T, T]] + for ds in [2, 3, 4], Rd in Rd_list + test_oracles(C(rand_herms(ds, Rd, T)), noise = T(1e-2), init_tol = Inf) + end +end + +function test_barrier(C::Type{Cones.LinMatrixIneq{T}}) where T + Random.seed!(1) + Ps = rand_herms(2, [T, Complex{T}], T) + barrier(s) = -logdet_pd(Hermitian(sum(s[i] * Ps[i] for i in eachindex(s)), :U)) + test_barrier(C(Ps), barrier) +end + +get_allocs(C::Type{Cones.LinMatrixIneq{T}}) where T = get_allocs(C(rand_herms(3, [T, Complex{T}, T, Complex{T}], T))) + + +# EpiNormInf +function test_oracles(C::Type{Cones.EpiNormInf{T, R}}) where {T, R} + for dw in [1, 2, 5] + test_oracles(C(1 + dim_vec(dw, R))) + end +end + +function test_barrier(C::Type{Cones.EpiNormInf{T, R}}) where {T, R} + dw = 2 + function barrier(s) + u = s[1] + w = new_vec(s[2:end], dw, R) + return -sum(log(abs2(u) - abs2(wi)) for wi in w) + (dw - 1) * log(u) + end + test_barrier(C(1 + dim_vec(dw, R)), barrier) +end + +get_allocs(C::Type{<:Cones.EpiNormInf}) = get_allocs(C(9)) + + +# EpiNormEucl +function test_oracles(C::Type{<:Cones.EpiNormEucl}) + for dw in [1, 2, 5] + test_oracles(C(1 + dw)) + end +end + +function test_barrier(C::Type{<:Cones.EpiNormEucl}) + function barrier(s) + (u, w) = (s[1], s[2:end]) + return -log(abs2(u) - sum(abs2, w)) + end + test_barrier(C(3), barrier) +end + +get_allocs(C::Type{<:Cones.EpiNormEucl}) = get_allocs(C(9)) + + +# EpiPerSquare +function test_oracles(C::Type{<:Cones.EpiPerSquare}) + for dw in [1, 2, 5] + test_oracles(C(2 + dw)) + end +end + +function test_barrier(C::Type{<:Cones.EpiPerSquare}) + function barrier(s) + (u, v, w) = (s[1], s[2], s[3:end]) + return -log(2 * u * v - sum(abs2, w)) + end + test_barrier(C(4), barrier) +end + +get_allocs(C::Type{<:Cones.EpiPerSquare}) = get_allocs(C(9)) + + +# EpiNormSpectral +function test_oracles(C::Type{<:Cones.EpiNormSpectral}) + for (dr, ds) in [(1, 1), (1, 2), (2, 2), (2, 4), (3, 4)] + test_oracles(C(dr, ds)) + end +end + +function test_barrier(C::Type{Cones.EpiNormSpectral{T, R}}) where {T, R} + (dr, ds) = (2, 3) + function barrier(s) + u = s[1] + W = reshape(new_vec(s[2:end], dr * ds, R), dr, ds) + return -logdet_pd(Hermitian(abs2(u) * I - W * W')) + (dr - 1) * log(u) + end + test_barrier(C(dr, ds), barrier) +end + +get_allocs(C::Type{<:Cones.EpiNormSpectral}) = get_allocs(C(2, 2)) + + +# MatrixEpiPerSquare +function test_oracles(C::Type{<:Cones.MatrixEpiPerSquare}) + for (dr, ds) in [(1, 1), (1, 2), (2, 2), (2, 4), (3, 4)] + test_oracles(C(dr, ds)) + end +end + +function test_barrier(C::Type{Cones.MatrixEpiPerSquare{T, R}}) where {T, R} + (dr, ds) = (2, 2) + du = dim_herm(dr, R) + function barrier(s) + U = new_mat_herm(s[1:du], dr, R) + v = s[du + 1] + W = reshape(new_vec(s[(du + 2):end], dr * ds, R), dr, ds) + return -logdet_pd(Hermitian(2 * v * U - W * W', :U)) + (dr - 1) * log(v) + end + test_barrier(C(dr, ds), barrier) +end + +get_allocs(C::Type{<:Cones.MatrixEpiPerSquare}) = get_allocs(C(2, 2)) + + +# GeneralizedPower +function test_oracles(C::Type{Cones.GeneralizedPower{T}}) where T + for (du, dw) in [(2, 1), (3, 2), (4, 1), (2, 4)] + test_oracles(C(rand_powers(T, du), dw)) + end +end + +function test_barrier(C::Type{Cones.GeneralizedPower{T}}) where T + (du, dw) = (2, 2) + alpha = rand_powers(T, du) + function barrier(s) + (u, w) = (s[1:du], s[(du + 1):end]) + return -log(exp(2 * sum(alpha[i] * log(u[i]) for i in eachindex(u))) - sum(abs2, w)) - sum((1 - alpha[i]) * log(u[i]) for i in eachindex(u)) + end + test_barrier(C(alpha, dw), barrier) +end + +get_allocs(C::Type{Cones.GeneralizedPower{T}}) where T = get_allocs(C(rand_powers(T, 4), 5)) + + +# HypoPowerMean +function test_oracles(C::Type{Cones.HypoPowerMean{T}}) where T + for dw in [1, 2, 5] + test_oracles(C(rand_powers(T, dw)), init_tol = 1e-2) + end + for dw in [15, 40, 100] + test_oracles(C(rand_powers(T, dw)), init_tol = 1e-1, init_only = true) + end +end + +function test_barrier(C::Type{Cones.HypoPowerMean{T}}) where T + alpha = rand_powers(T, 3) + function barrier(s) + (u, w) = (s[1], s[2:end]) + return -log(exp(sum(alpha[i] * log(w[i]) for i in eachindex(w))) - u) - sum(log, w) + end + test_barrier(C(alpha), barrier) +end + +get_allocs(C::Type{Cones.HypoPowerMean{T}}) where T = get_allocs(C(rand_powers(T, 8))) + + +# HypoGeoMean +function test_oracles(C::Type{<:Cones.HypoGeoMean}) + for dw in [1, 2, 5] + test_oracles(C(1 + dw)) + end +end + +function test_barrier(C::Type{<:Cones.HypoGeoMean}) + function barrier(s) + (u, w) = (s[1], s[2:end]) + sumlogw = sum(log, w) + return -log(exp(sumlogw / length(w)) - u) - sumlogw + end + test_barrier(C(3), barrier) +end + +get_allocs(C::Type{<:Cones.HypoGeoMean}) = get_allocs(C(9)) + + +# HypoRootdetTri +function test_oracles(C::Type{Cones.HypoRootdetTri{T, R}}) where {T, R} + for dW in [1, 2, 4] + test_oracles(C(1 + dim_herm(dW, R))) + end +end + +function test_barrier(C::Type{Cones.HypoRootdetTri{T, R}}) where {T, R} + dW = 3 + function barrier(s) + (u, w) = (s[1], s[2:end]) + logdet_W = logdet_pd(new_mat_herm(w, dW, R)) + return -log(exp(logdet_W / dW) - u) - logdet_W + end + test_barrier(C(1 + dim_herm(dW, R)), barrier) +end + +get_allocs(C::Type{Cones.HypoRootdetTri{T, R}}) where {T, R} = get_allocs(C(1 + dim_herm(3, R))) + + +# HypoPerLog +function test_oracles(C::Type{<:Cones.HypoPerLog}) + for dw in [1, 2, 5] + test_oracles(C(2 + dw), init_tol = 1e-5) + end + for dw in [15, 40, 100] + test_oracles(C(2 + dw), init_tol = 1e-1, init_only = true) + end +end + +function test_barrier(C::Type{<:Cones.HypoPerLog}) + function barrier(s) + (u, v, w) = (s[1], s[2], s[3:end]) + return -log(v * sum(log(wi / v) for wi in w) - u) - log(v) - sum(log, w) + end + test_barrier(C(4), barrier) +end + +get_allocs(C::Type{<:Cones.HypoPerLog}) = get_allocs(C(9)) + + +# HypoPerLogdetTri +function test_oracles(C::Type{Cones.HypoPerLogdetTri{T, R}}) where {T, R} + for dW in [1, 2, 4] + test_oracles(C(2 + dim_herm(dW, R)), init_tol = 1e-4) + end + for dW in [8, 12] + test_oracles(C(2 + dim_herm(dW, R)), init_tol = 1e-1, init_only = true) + end +end + +function test_barrier(C::Type{Cones.HypoPerLogdetTri{T, R}}) where {T, R} + dW = 3 + function barrier(s) + (u, v, w) = (s[1], s[2], s[3:end]) + W = new_mat_herm(w, dW, R) + return -log(v * logdet_pd(W / v) - u) - log(v) - logdet_pd(W) + end + test_barrier(C(2 + dim_herm(dW, R)), barrier) +end + +get_allocs(C::Type{Cones.HypoPerLogdetTri{T, R}}) where {T, R} = get_allocs(C(2 + dim_herm(3, R))) + + +# EpiPerEntropy +function test_oracles(C::Type{<:Cones.EpiPerEntropy}) + for dw in [1, 2, 5] + test_oracles(C(2 + dw), init_tol = 1e-5) + end + for dw in [15, 40, 100] + test_oracles(C(2 + dw), init_tol = 1e-1, init_only = true) + end +end + +function test_barrier(C::Type{<:Cones.EpiPerEntropy}) + function barrier(s) + (u, v, w) = (s[1], s[2], s[3:end]) + return -log(u - sum(wi * log(wi / v) for wi in w)) - log(v) - sum(log, w) + end + test_barrier(C(4), barrier) +end + +get_allocs(C::Type{<:Cones.EpiPerEntropy}) = get_allocs(C(9)) + + +# EpiPerTraceEntropyTri +function test_oracles(C::Type{<:Cones.EpiPerTraceEntropyTri}) + for dW in [1, 2, 4] + test_oracles(C(2 + Cones.svec_length(dW)), init_tol = 1e-4) + end + for dW in [8, 12] + test_oracles(C(2 + Cones.svec_length(dW)), init_tol = 1e-1, init_only = true) + end +end + +function test_barrier(C::Type{Cones.EpiPerTraceEntropyTri{T}}) where T + dW = 3 + function barrier(s) + (u, v, w) = (s[1], s[2], s[3:end]) + W = new_mat_herm(w, dW, T) + return -log(u - dot(W, log(W / v))) - log(v) - logdet_pd(W) + end + test_barrier(C(2 + Cones.svec_length(dW)), barrier) +end + +get_allocs(C::Type{<:Cones.EpiPerTraceEntropyTri}) = get_allocs(C(8)) + + +# EpiRelEntropy +function test_oracles(C::Type{<:Cones.EpiRelEntropy}) + for dw in [1, 2, 4] + test_oracles(C(1 + 2 * dw), init_tol = 1e-5) + end + for dw in [15, 40, 100] + test_oracles(C(1 + 2 * dw), init_tol = 1e-1, init_only = true) + end +end + +function test_barrier(C::Type{<:Cones.EpiRelEntropy}) + function barrier(s) + (u, v, w) = (s[1], s[2:2:(end - 1)], s[3:2:end]) + return -log(u - sum(wi * log(wi / vi) for (vi, wi) in zip(v, w))) - sum(log, v) - sum(log, w) + end + test_barrier(C(5), barrier) +end + +get_allocs(C::Type{<:Cones.EpiRelEntropy}) = get_allocs(C(9)) + + +# EpiTraceRelEntropyTri +function test_oracles(C::Type{<:Cones.EpiTraceRelEntropyTri}) + for dW in [1, 2, 4] + test_oracles(C(1 + 2 * Cones.svec_length(dW)), init_tol = 1e-4) + end + for dW in [6, 10] + test_oracles(C(1 + 2 * Cones.svec_length(dW)), init_tol = 1e-1, init_only = true) + end +end + +function test_barrier(C::Type{Cones.EpiTraceRelEntropyTri{T}}) where T + dW = 3 + dw = Cones.svec_length(dW) + function barrier(s) + (u, v, w) = (s[1], s[1 .+ (1:dw)], s[(2 + dw):end]) + V = new_mat_herm(v, dW, T) + W = new_mat_herm(w, dW, T) + return -log(u - dot(W, log(W) - log(V))) - logdet_pd(V) - logdet_pd(W) + end + test_barrier(C(1 + 2 * dw), barrier, tol = 1e8 * eps(T)) +end + +get_allocs(C::Type{<:Cones.EpiTraceRelEntropyTri}) = get_allocs(C(13)) + + +# WSOSInterpNonnegative +# TODO test complex case, but need complex interpolation +function test_oracles(C::Type{Cones.WSOSInterpNonnegative{T, R}}) where {T, R} + for (num_vars, halfdeg) in [(1, 1), (1, 3), (2, 1), (2, 2), (3, 1)] + (d, Ps) = rand_interp(num_vars, halfdeg, R) + test_oracles(C(d, Ps), init_tol = Inf) + end +end + +function test_barrier(C::Type{Cones.WSOSInterpNonnegative{T, R}}) where {T, R} + (d, Ps) = rand_interp(2, 1, R) + barrier(s) = -sum(logdet_pd(Hermitian(P' * Diagonal(s) * P)) for P in Ps) + test_barrier(C(d, Ps), barrier) +end + +get_allocs(C::Type{Cones.WSOSInterpNonnegative{T, R}}) where {T, R} = get_allocs(C(rand_interp(3, 1, R)...)) + + +# WSOSInterpPosSemidefTri +function test_oracles(C::Type{Cones.WSOSInterpPosSemidefTri{T}}) where T + for (num_vars, halfdeg, ds) in [(1, 1, 1), (1, 1, 4), (2, 2, 1), (2, 1, 3), (3, 1, 2)] + (d, Ps) = rand_interp(num_vars, halfdeg, T) + test_oracles(C(ds, d, Ps), init_tol = Inf) + end +end + +function test_barrier(C::Type{Cones.WSOSInterpPosSemidefTri{T}}) where T + (d, Ps) = rand_interp(1, 1, T) + ds = 3 + invrt2 = inv(sqrt(T(2))) + function barrier(s) + function ldlamP(P) + dt = size(P, 2) + lam = similar(s, ds * dt, ds * dt) + for i in 1:ds, j in 1:i + lamij = P' * Diagonal(s[Cones.block_idxs(d, Cones.svec_idx(i, j))]) * P + if i != j + lamij .*= invrt2 + end + lam[Cones.block_idxs(dt, i), Cones.block_idxs(dt, j)] = lamij + end + return -logdet_pd(Hermitian(lam, :L)) + end + return sum(ldlamP, Ps) + end + test_barrier(C(ds, d, Ps), barrier) +end + +get_allocs(C::Type{Cones.WSOSInterpPosSemidefTri{T}}) where T = get_allocs(C(2, rand_interp(2, 1, T)...)) + + +# WSOSInterpEpiNormEucl +function test_oracles(C::Type{Cones.WSOSInterpEpiNormEucl{T}}) where T + for (num_vars, halfdeg, ds) in [(1, 1, 1), (1, 2, 3), (2, 2, 2), (3, 1, 1)] + (d, Ps) = rand_interp(num_vars, halfdeg, T) + test_oracles(C(1 + ds, d, Ps), init_tol = Inf) + end +end + +function test_barrier(C::Type{Cones.WSOSInterpEpiNormEucl{T}}) where T + (d, Ps) = rand_interp(1, 1, T) + ds = 2 + invrt2 = inv(sqrt(T(2))) + function barrier(s) + function ldlamP(P) + lam = P' * Diagonal(s[1:d]) * P + lam1fact = cholesky(Hermitian(lam, :L)) + PL1 = lam1fact.L \ P' + for i in 1:ds + lamLi = PL1 * Diagonal(s[Cones.block_idxs(d, 1 + i)]) * P + lam -= lamLi' * lamLi + end + return -logdet(lam1fact) - logdet_pd(Hermitian(lam)) + end + return sum(ldlamP, Ps) + end + test_barrier(C(1 + ds, d, Ps), barrier) +end + +get_allocs(C::Type{Cones.WSOSInterpEpiNormEucl{T}}) where T = get_allocs(C(3, rand_interp(2, 1, T)...)) + + +# WSOSInterpEpiNormOne +function test_oracles(C::Type{Cones.WSOSInterpEpiNormOne{T}}) where T + for (num_vars, halfdeg, ds) in [(1, 1, 1), (1, 2, 3), (2, 2, 2), (3, 1, 1)] + (d, Ps) = rand_interp(num_vars, halfdeg, T) + test_oracles(C(1 + ds, d, Ps), init_tol = Inf) + end +end + +function test_barrier(C::Type{Cones.WSOSInterpEpiNormOne{T}}) where T + (d, Ps) = rand_interp(1, 1, T) + ds = 2 + invrt2 = inv(sqrt(T(2))) + function barrier(s) + function ldlamP(P) + lam = P' * Diagonal(s[1:d]) * P + lam1fact = cholesky(Hermitian(lam, :L)) + PL1 = lam1fact.L \ P' + lamLs = [PL1 * Diagonal(s[Cones.block_idxs(d, 1 + i)]) * P for i in 1:ds] + lams = [Hermitian(lam - lamLi' * lamLi) for lamLi in lamLs] + return -logdet(lam1fact) - sum(logdet_pd, lams) + end + return sum(ldlamP, Ps) + end + test_barrier(C(1 + ds, d, Ps), barrier) +end + +get_allocs(C::Type{Cones.WSOSInterpEpiNormOne{T}}) where T = get_allocs(C(3, rand_interp(2, 1, T)...)) diff --git a/test/moi.jl b/test/moi.jl index a4455493c..cdebb05f7 100644 --- a/test/moi.jl +++ b/test/moi.jl @@ -24,24 +24,23 @@ unit_exclude = [ conic_exclude = String[ # "lin", + # "sdp", # "norminf", # "normone", # "soc", # "rsoc", - # "exp", - # "dualexp", + # "normspec", + # "normnuc", # "pow", # "dualpow", # "geomean", - # "relentr", - # "normspec", - # "normnuc", - # "sdp", - # "logdet", # TODO currently failing due to barrier with high parameter # "rootdet", - # square logdet and rootdet cones not handled - "logdets", "rootdets", + # "exp", + # "dualexp", + # "logdet", + "logdets", + # "relentr", ] function test_moi(T::Type{<:Real}; solver_options...) diff --git a/test/moicones.jl b/test/moicones.jl index 2ced46e53..c814c3f79 100644 --- a/test/moicones.jl +++ b/test/moicones.jl @@ -14,6 +14,14 @@ import Hypatia.Cones function test_moi_cones(T::Type{<:Real}) # MOI predefined cones + @testset "PositiveSemidefiniteConeTriangle" begin + moi_cone = MOI.PositiveSemidefiniteConeTriangle(3) + hyp_cone = Hypatia.cone_from_moi(T, moi_cone) + @test hyp_cone isa Cones.PosSemidefTri{T, T} + @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 6 + @test !Cones.use_dual_barrier(hyp_cone) + end + @testset "NormInfinityCone" begin moi_cone = MOI.NormInfinityCone(3) hyp_cone = Hypatia.cone_from_moi(T, moi_cone) @@ -46,19 +54,19 @@ function test_moi_cones(T::Type{<:Real}) @test !Cones.use_dual_barrier(hyp_cone) end - @testset "ExponentialCone" begin - moi_cone = MOI.ExponentialCone() + @testset "NormSpectralCone" begin + moi_cone = MOI.NormSpectralCone(2, 3) hyp_cone = Hypatia.cone_from_moi(T, moi_cone) - @test hyp_cone isa Cones.HypoPerLog{T} - @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3 + @test hyp_cone isa Cones.EpiNormSpectral{T, T} + @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 7 @test !Cones.use_dual_barrier(hyp_cone) end - @testset "DualExponentialCone" begin - moi_cone = MOI.DualExponentialCone() + @testset "NormNuclearCone" begin + moi_cone = MOI.NormNuclearCone(2, 3) hyp_cone = Hypatia.cone_from_moi(T, moi_cone) - @test hyp_cone isa Cones.HypoPerLog{T} - @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3 + @test hyp_cone isa Cones.EpiNormSpectral{T, T} + @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 7 @test Cones.use_dual_barrier(hyp_cone) end @@ -66,7 +74,7 @@ function test_moi_cones(T::Type{<:Real}) iT5 = inv(T(5)) moi_cone = MOI.PowerCone{T}(iT5) hyp_cone = Hypatia.cone_from_moi(T, moi_cone) - @test hyp_cone isa Cones.Power{T} + @test hyp_cone isa Cones.GeneralizedPower{T} @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3 @test !Cones.use_dual_barrier(hyp_cone) @test hyp_cone.alpha == T[iT5, 1 - iT5] @@ -76,7 +84,7 @@ function test_moi_cones(T::Type{<:Real}) iT5 = inv(T(5)) moi_cone = MOI.DualPowerCone{T}(iT5) hyp_cone = Hypatia.cone_from_moi(T, moi_cone) - @test hyp_cone isa Cones.Power{T} + @test hyp_cone isa Cones.GeneralizedPower{T} @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3 @test Cones.use_dual_barrier(hyp_cone) @test hyp_cone.alpha == T[iT5, 1 - iT5] @@ -90,38 +98,30 @@ function test_moi_cones(T::Type{<:Real}) @test !Cones.use_dual_barrier(hyp_cone) end - @testset "RelativeEntropyCone" begin - moi_cone = MOI.RelativeEntropyCone(3) + @testset "RootDetConeTriangle" begin + moi_cone = MOI.RootDetConeTriangle(3) hyp_cone = Hypatia.cone_from_moi(T, moi_cone) - @test hyp_cone isa Cones.EpiRelEntropy{T} - @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3 + @test hyp_cone isa Cones.HypoRootdetTri{T, T} + @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 7 @test !Cones.use_dual_barrier(hyp_cone) end - @testset "NormSpectralCone" begin - moi_cone = MOI.NormSpectralCone(2, 3) + @testset "ExponentialCone" begin + moi_cone = MOI.ExponentialCone() hyp_cone = Hypatia.cone_from_moi(T, moi_cone) - @test hyp_cone isa Cones.EpiNormSpectral{T, T} - @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 7 + @test hyp_cone isa Cones.HypoPerLog{T} + @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3 @test !Cones.use_dual_barrier(hyp_cone) end - @testset "NormNuclearCone" begin - moi_cone = MOI.NormNuclearCone(2, 3) + @testset "DualExponentialCone" begin + moi_cone = MOI.DualExponentialCone() hyp_cone = Hypatia.cone_from_moi(T, moi_cone) - @test hyp_cone isa Cones.EpiNormSpectral{T, T} - @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 7 + @test hyp_cone isa Cones.HypoPerLog{T} + @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3 @test Cones.use_dual_barrier(hyp_cone) end - @testset "PositiveSemidefiniteConeTriangle" begin - moi_cone = MOI.PositiveSemidefiniteConeTriangle(3) - hyp_cone = Hypatia.cone_from_moi(T, moi_cone) - @test hyp_cone isa Cones.PosSemidefTri{T, T} - @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 6 - @test !Cones.use_dual_barrier(hyp_cone) - end - @testset "LogDetConeTriangle" begin moi_cone = MOI.LogDetConeTriangle(3) hyp_cone = Hypatia.cone_from_moi(T, moi_cone) @@ -130,11 +130,11 @@ function test_moi_cones(T::Type{<:Real}) @test !Cones.use_dual_barrier(hyp_cone) end - @testset "RootDetConeTriangle" begin - moi_cone = MOI.RootDetConeTriangle(3) + @testset "RelativeEntropyCone" begin + moi_cone = MOI.RelativeEntropyCone(3) hyp_cone = Hypatia.cone_from_moi(T, moi_cone) - @test hyp_cone isa Cones.HypoRootdetTri{T, T} - @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 7 + @test hyp_cone isa Cones.EpiRelEntropy{T} + @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3 @test !Cones.use_dual_barrier(hyp_cone) end @@ -147,78 +147,76 @@ function test_moi_cones(T::Type{<:Real}) @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3 end - @testset "EpiNormInfinity" begin - moi_cone = Hypatia.EpiNormInfinityCone{T, T}(3) + @testset "PosSemidefTri" begin + moi_cone = Hypatia.PosSemidefTriCone{T, T}(6) hyp_cone = Hypatia.cone_from_moi(T, moi_cone) - @test hyp_cone isa Cones.EpiNormInf{T, T} - @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3 + @test hyp_cone isa Cones.PosSemidefTri{T, T} + @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 6 - moi_cone = Hypatia.EpiNormInfinityCone{T, Complex{T}}(5) + moi_cone = Hypatia.PosSemidefTriCone{T, Complex{T}}(9) hyp_cone = Hypatia.cone_from_moi(T, moi_cone) - @test hyp_cone isa Cones.EpiNormInf{T, Complex{T}} - @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 5 + @test hyp_cone isa Cones.PosSemidefTri{T, Complex{T}} + @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 9 end - @testset "EpiNormEucl" begin - moi_cone = Hypatia.EpiNormEuclCone{T}(3) + @testset "DoublyNonnegativeTri" begin + moi_cone = Hypatia.DoublyNonnegativeTriCone{T}(3) hyp_cone = Hypatia.cone_from_moi(T, moi_cone) - @test hyp_cone isa Cones.EpiNormEucl{T} + @test hyp_cone isa Cones.DoublyNonnegativeTri{T} @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3 end - @testset "EpiPerSquare" begin - moi_cone = Hypatia.EpiPerSquareCone{T}(4) - hyp_cone = Hypatia.cone_from_moi(T, moi_cone) - @test hyp_cone isa Cones.EpiPerSquare{T} - @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 4 - end + @testset "PosSemidefTriSparse" begin + if T <: LinearAlgebra.BlasReal # only works with BLAS real types + Random.seed!(1) + side = 5 + (row_idxs, col_idxs, _) = SparseArrays.findnz(tril!(SparseArrays.sprand(Bool, side, side, 0.3)) + I) - @testset "HypoPerLog" begin - moi_cone = Hypatia.HypoPerLogCone{T}(4) - hyp_cone = Hypatia.cone_from_moi(T, moi_cone) - @test hyp_cone isa Cones.HypoPerLog{T} - @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 4 + moi_cone = Hypatia.PosSemidefTriSparseCone{T, T}(side, row_idxs, col_idxs) + hyp_cone = Hypatia.cone_from_moi(T, moi_cone) + @test hyp_cone isa Cones.PosSemidefTriSparse{T, T} + @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == length(row_idxs) + + moi_cone = Hypatia.PosSemidefTriSparseCone{T, Complex{T}}(side, row_idxs, col_idxs) + hyp_cone = Hypatia.cone_from_moi(T, moi_cone) + @test hyp_cone isa Cones.PosSemidefTriSparse{T, Complex{T}} + @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 2 * length(row_idxs) - side + end end - @testset "EpiRelEntropy" begin - moi_cone = Hypatia.EpiRelEntropyCone{T}(5) + @testset "LinMatrixIneq" begin + As = [Symmetric(Matrix(one(T) * I, 2, 2)), Hermitian(Complex{T}[1 0; 0 -1])] + moi_cone = Hypatia.LinMatrixIneqCone{T}(As) hyp_cone = Hypatia.cone_from_moi(T, moi_cone) - @test hyp_cone isa Cones.EpiRelEntropy{T} - @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 5 + @test hyp_cone isa Cones.LinMatrixIneq{T} + @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 2 + @test hyp_cone.As == As end - @testset "EpiPerEntropy" begin - moi_cone = Hypatia.EpiPerEntropyCone{T}(3) + @testset "EpiNormInfinity" begin + moi_cone = Hypatia.EpiNormInfinityCone{T, T}(3) hyp_cone = Hypatia.cone_from_moi(T, moi_cone) - @test hyp_cone isa Cones.EpiPerEntropy{T} + @test hyp_cone isa Cones.EpiNormInf{T, T} @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3 - end - @testset "HypoGeoMean" begin - moi_cone = Hypatia.HypoGeoMeanCone{T}(3) + moi_cone = Hypatia.EpiNormInfinityCone{T, Complex{T}}(5) hyp_cone = Hypatia.cone_from_moi(T, moi_cone) - @test hyp_cone isa Cones.HypoGeoMean{T} - @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3 + @test hyp_cone isa Cones.EpiNormInf{T, Complex{T}} + @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 5 end - @testset "HypoPowerMean" begin - alpha = rand(T, 2) - alpha ./= sum(alpha) - moi_cone = Hypatia.HypoPowerMeanCone{T}(alpha) + @testset "EpiNormEucl" begin + moi_cone = Hypatia.EpiNormEuclCone{T}(3) hyp_cone = Hypatia.cone_from_moi(T, moi_cone) - @test hyp_cone isa Cones.HypoPowerMean{T} + @test hyp_cone isa Cones.EpiNormEucl{T} @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3 - @test hyp_cone.alpha == alpha end - @testset "Power" begin - alpha = rand(T, 2) - alpha ./= sum(alpha) - moi_cone = Hypatia.PowerCone{T}(alpha, 3) + @testset "EpiPerSquare" begin + moi_cone = Hypatia.EpiPerSquareCone{T}(4) hyp_cone = Hypatia.cone_from_moi(T, moi_cone) - @test hyp_cone isa Cones.Power{T} - @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 5 - @test hyp_cone.alpha == alpha + @test hyp_cone isa Cones.EpiPerSquare{T} + @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 4 end @testset "EpiNormSpectral" begin @@ -245,43 +243,50 @@ function test_moi_cones(T::Type{<:Real}) @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 17 end - @testset "LinMatrixIneq" begin - As = [Symmetric(Matrix(one(T) * I, 2, 2)), Hermitian(Complex{T}[1 0; 0 -1])] - moi_cone = Hypatia.LinMatrixIneqCone{T}(As) + @testset "GeneralizedPower" begin + alpha = rand(T, 2) + alpha ./= sum(alpha) + moi_cone = Hypatia.GeneralizedPowerCone{T}(alpha, 3) hyp_cone = Hypatia.cone_from_moi(T, moi_cone) - @test hyp_cone isa Cones.LinMatrixIneq{T} - @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 2 - @test hyp_cone.As == As + @test hyp_cone isa Cones.GeneralizedPower{T} + @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 5 + @test hyp_cone.alpha == alpha end - @testset "PosSemidefTri" begin - moi_cone = Hypatia.PosSemidefTriCone{T, T}(6) + @testset "HypoPowerMean" begin + alpha = rand(T, 2) + alpha ./= sum(alpha) + moi_cone = Hypatia.HypoPowerMeanCone{T}(alpha) hyp_cone = Hypatia.cone_from_moi(T, moi_cone) - @test hyp_cone isa Cones.PosSemidefTri{T, T} - @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 6 + @test hyp_cone isa Cones.HypoPowerMean{T} + @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3 + @test hyp_cone.alpha == alpha + end - moi_cone = Hypatia.PosSemidefTriCone{T, Complex{T}}(9) + @testset "HypoGeoMean" begin + moi_cone = Hypatia.HypoGeoMeanCone{T}(3) hyp_cone = Hypatia.cone_from_moi(T, moi_cone) - @test hyp_cone isa Cones.PosSemidefTri{T, Complex{T}} - @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 9 + @test hyp_cone isa Cones.HypoGeoMean{T} + @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3 end - @testset "PosSemidefTriSparse" begin - if T <: LinearAlgebra.BlasReal # only works with BLAS real types - Random.seed!(1) - side = 5 - (row_idxs, col_idxs, _) = SparseArrays.findnz(tril!(SparseArrays.sprand(Bool, side, side, 0.3)) + I) + @testset "HypoRootdetTri" begin + moi_cone = Hypatia.HypoRootdetTriCone{T, T}(7) + hyp_cone = Hypatia.cone_from_moi(T, moi_cone) + @test hyp_cone isa Cones.HypoRootdetTri{T, T} + @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 7 - moi_cone = Hypatia.PosSemidefTriSparseCone{T, T}(side, row_idxs, col_idxs) - hyp_cone = Hypatia.cone_from_moi(T, moi_cone) - @test hyp_cone isa Cones.PosSemidefTriSparse{T, T} - @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == length(row_idxs) + moi_cone = Hypatia.HypoRootdetTriCone{T, Complex{T}}(10) + hyp_cone = Hypatia.cone_from_moi(T, moi_cone) + @test hyp_cone isa Cones.HypoRootdetTri{T, Complex{T}} + @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 10 + end - moi_cone = Hypatia.PosSemidefTriSparseCone{T, Complex{T}}(side, row_idxs, col_idxs) - hyp_cone = Hypatia.cone_from_moi(T, moi_cone) - @test hyp_cone isa Cones.PosSemidefTriSparse{T, Complex{T}} - @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 2 * length(row_idxs) - side - end + @testset "HypoPerLog" begin + moi_cone = Hypatia.HypoPerLogCone{T}(4) + hyp_cone = Hypatia.cone_from_moi(T, moi_cone) + @test hyp_cone isa Cones.HypoPerLog{T} + @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 4 end @testset "HypoPerLogdetTri" begin @@ -296,23 +301,25 @@ function test_moi_cones(T::Type{<:Real}) @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 11 end - @testset "HypoRootdetTri" begin - moi_cone = Hypatia.HypoRootdetTriCone{T, T}(7) + @testset "EpiPerEntropy" begin + moi_cone = Hypatia.EpiPerEntropyCone{T}(3) hyp_cone = Hypatia.cone_from_moi(T, moi_cone) - @test hyp_cone isa Cones.HypoRootdetTri{T, T} - @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 7 + @test hyp_cone isa Cones.EpiPerEntropy{T} + @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3 + end - moi_cone = Hypatia.HypoRootdetTriCone{T, Complex{T}}(10) + @testset "EpiPerTraceEntropyTriCone" begin + moi_cone = Hypatia.EpiPerTraceEntropyTriCone{T}(3) hyp_cone = Hypatia.cone_from_moi(T, moi_cone) - @test hyp_cone isa Cones.HypoRootdetTri{T, Complex{T}} - @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 10 + @test hyp_cone isa Cones.EpiPerTraceEntropyTri{T} + @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3 end - @testset "DoublyNonnegativeTri" begin - moi_cone = Hypatia.DoublyNonnegativeTriCone{T}(3) + @testset "EpiRelEntropy" begin + moi_cone = Hypatia.EpiRelEntropyCone{T}(5) hyp_cone = Hypatia.cone_from_moi(T, moi_cone) - @test hyp_cone isa Cones.DoublyNonnegativeTri{T} - @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3 + @test hyp_cone isa Cones.EpiRelEntropy{T} + @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 5 end @testset "EpiTraceRelEntropyTriCone" begin @@ -322,13 +329,6 @@ function test_moi_cones(T::Type{<:Real}) @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3 end - @testset "EpiPerTraceEntropyTriCone" begin - moi_cone = Hypatia.EpiPerTraceEntropyTriCone{T}(3) - hyp_cone = Hypatia.cone_from_moi(T, moi_cone) - @test hyp_cone isa Cones.EpiPerTraceEntropyTri{T} - @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3 - end - @testset "WSOSInterpNonnegative" begin Ps = [rand(T, 3, 2), rand(T, 3, 1)] moi_cone = Hypatia.WSOSInterpNonnegativeCone{T, T}(3, Ps) @@ -354,6 +354,15 @@ function test_moi_cones(T::Type{<:Real}) @test hyp_cone.Ps == Ps end + @testset "WSOSInterpEpiNormOne" begin + Ps = [rand(T, 3, 2), rand(T, 3, 1)] + moi_cone = Hypatia.WSOSInterpEpiNormOneCone{T}(2, 3, Ps) + hyp_cone = Hypatia.cone_from_moi(T, moi_cone) + @test hyp_cone isa Cones.WSOSInterpEpiNormOne{T} + @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 6 + @test hyp_cone.Ps == Ps + end + @testset "WSOSInterpEpiNormEucl" begin Ps = [rand(T, 3, 2), rand(T, 3, 1)] moi_cone = Hypatia.WSOSInterpEpiNormEuclCone{T}(2, 3, Ps) diff --git a/test/nativeinstances.jl b/test/nativeinstances.jl index d099680eb..acb7e4104 100644 --- a/test/nativeinstances.jl +++ b/test/nativeinstances.jl @@ -161,6 +161,84 @@ function inconsistent2(T; options...) @test r.status == Solvers.DualInconsistent end +function primalinfeas1(T; options...) + tol = test_tol(T) + c = T[1, 0] + A = T[1 1] + b = [-T(2)] + G = SparseMatrixCSC(-one(T) * I, 2, 2) + h = zeros(T, 2) + cones = Cone{T}[Cones.Nonnegative{T}(2)] + + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.PrimalInfeasible +end + +function primalinfeas2(T; options...) + tol = test_tol(T) + c = T[1, 1, 1] + A = zeros(T, 0, 3) + b = T[] + G = vcat(SparseMatrixCSC(-one(T) * I, 3, 3), Diagonal([one(T), one(T), -one(T)])) + h = vcat(zeros(T, 3), one(T), one(T), -T(2)) + cones = Cone{T}[Cones.EpiNormEucl{T}(3), Cones.Nonnegative{T}(3)] + + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.PrimalInfeasible +end + +function primalinfeas3(T; options...) + tol = test_tol(T) + c = zeros(T, 3) + A = SparseMatrixCSC(-one(T) * I, 3, 3) + b = [one(T), one(T), T(3)] + G = SparseMatrixCSC(-one(T) * I, 3, 3) + h = zeros(T, 3) + cones = Cone{T}[Cones.HypoPerLog{T}(3)] + + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.PrimalInfeasible +end + +function dualinfeas1(T; options...) + tol = test_tol(T) + c = T[-1, -1, 0] + A = zeros(T, 0, 3) + b = T[] + G = repeat(SparseMatrixCSC(-one(T) * I, 3, 3), 2, 1) + h = zeros(T, 6) + cones = Cone{T}[Cones.EpiNormInf{T, T}(3), Cones.EpiNormInf{T, T}(3, use_dual = true)] + + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.DualInfeasible +end + +function dualinfeas2(T; options...) + tol = test_tol(T) + c = T[-1, 0] + A = zeros(T, 0, 2) + b = T[] + G = T[-1 0; 0 0; 0 -1] + h = T[0, 1, 0] + cones = Cone{T}[Cones.EpiPerSquare{T}(3)] + + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.DualInfeasible +end + +function dualinfeas3(T; options...) + tol = test_tol(T) + c = T[0, 1, 1, 0] + A = zeros(T, 0, 4) + b = T[] + G = SparseMatrixCSC(-one(T) * I, 4, 4) + h = zeros(T, 4) + cones = Cone{T}[Cones.EpiPerSquare{T}(4)] + + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.DualInfeasible +end + function nonnegative1(T; options...) tol = test_tol(T) Random.seed!(1) @@ -223,229 +301,223 @@ function nonnegative4(T; options...) @test r.z ≈ [2, 2, 0] atol=tol rtol=tol end -function epinorminf1(T; options...) +function possemideftri1(T; options...) tol = test_tol(T) - Tirt2 = inv(sqrt(T(2))) - c = T[0, -1, -1] - A = T[1 0 0; 0 1 0] - b = [one(T), Tirt2] - G = SparseMatrixCSC(-one(T) * I, 3, 3) + c = T[0, -1, 0] + A = T[1 0 0; 0 0 1] + b = T[0.5, 1] + G = -one(T) * I h = zeros(T, 3) - cones = Cone{T}[Cones.EpiNormInf{T, T}(3)] + cones = Cone{T}[Cones.PosSemidefTri{T, T}(3)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - @test r.primal_obj ≈ -1 - Tirt2 atol=tol rtol=tol - @test r.x ≈ [1, Tirt2, 1] atol=tol rtol=tol - @test r.y ≈ [1, 1] atol=tol rtol=tol -end - -function epinorminf2(T; options...) - tol = 10 * test_tol(T) - l = 3 - L = 2l + 1 - c = collect(T, -l:l) - A = spzeros(T, 2, L) - A[1, 1] = A[1, L] = A[2, 1] = 1; A[2, L] = -1 - b = T[0, 0] - G = [spzeros(T, 1, L); sparse(one(T) * I, L, L); spzeros(T, 1, L); sparse(T(2) * I, L, L)] - h = zeros(T, 2L + 2); h[1] = 1; h[L + 2] = 1 - cones = Cone{T}[Cones.EpiNormInf{T, T}(L + 1, use_dual = true), Cones.EpiNormInf{T, T}(L + 1, use_dual = false)] - - r = build_solve_check(c, A, b, G, h, cones, tol; obj_offset = one(T), options...) - @test r.status == Solvers.Optimal - @test r.primal_obj ≈ -l + 2 atol=tol rtol=tol - @test r.x[2] ≈ 0.5 atol=tol rtol=tol - @test r.x[end - 1] ≈ -0.5 atol=tol rtol=tol - @test sum(abs, r.x) ≈ 1 atol=tol rtol=tol -end - -function epinorminf3(T; options...) - tol = test_tol(T) - c = T[1, 0, 0, 0, 0, 0] - A = zeros(T, 0, 6) - b = zeros(T, 0) - G = Diagonal(-one(T) * I, 6) - h = zeros(T, 6) - - for use_dual in (false, true) - cones = Cone{T}[Cones.EpiNormInf{T, T}(6, use_dual = use_dual)] - - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test r.primal_obj ≈ 0 atol=tol rtol=tol - @test norm(r.x) ≈ 0 atol=tol rtol=tol - end + @test r.primal_obj ≈ -one(T) atol=tol rtol=tol + @test r.x[2] ≈ one(T) atol=tol rtol=tol end -function epinorminf4(T; options...) +function possemideftri2(T; options...) tol = test_tol(T) - c = T[0, 1, -1] - A = T[1 0 0; 0 1 0] - b = T[1, -0.4] - G = -one(T) * I + c = T[0, -1, 0] + A = T[1 0 1] + b = T[0] + G = Diagonal(-one(T) * I, 3) h = zeros(T, 3) - cones = Cone{T}[Cones.EpiNormInf{T, T}(3, use_dual = true)] + cones = Cone{T}[Cones.PosSemidefTri{T, T}(3)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - @test r.primal_obj ≈ -1 atol=tol rtol=tol - @test r.x ≈ [1, -0.4, 0.6] atol=tol rtol=tol - @test r.y ≈ [1, 0] atol=tol rtol=tol + @test r.primal_obj ≈ 0 atol=tol rtol=tol + @test norm(r.x) ≈ 0 atol=tol rtol=tol end -function epinorminf5(T; options...) +function possemideftri3(T; options...) tol = test_tol(T) + rt2 = sqrt(T(2)) Random.seed!(1) - c = T[1, 0, 0, 0, 0, 0] - A = rand(T(-9):T(9), 6, 6) - b = vec(sum(A, dims = 2)) - G = rand(T, 6, 6) - h = vec(sum(G, dims = 2)) - cones = Cone{T}[Cones.EpiNormInf{T, T}(6, use_dual = true)] + c = T[1] + A = zeros(T, 0, 1) + b = T[] + rand_mat = Hermitian(rand(T, 2, 2), :U) + G = reshape(T[-1, 0, -1], 3, 1) + h = -Cones.smat_to_svec!(zeros(T, 3), rand_mat, rt2) + cones = Cone{T}[Cones.PosSemidefTri{T, T}(3)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - @test r.primal_obj ≈ 1 atol=tol rtol=tol + eig_max = maximum(eigvals(rand_mat)) + @test r.primal_obj ≈ eig_max atol=tol rtol=tol + @test r.x[1] ≈ eig_max atol=tol rtol=tol end -function epinorminf6(T; options...) +function possemideftri4(T; options...) tol = test_tol(T) - c = T[0, -1, -1, -1, -1] - A = T[1 0 0 0 0; 0 1 0 0 0; 0 0 0 1 0; 0 0 0 0 1] - b = T[2, 0, 1, 0] - G = SparseMatrixCSC(-one(T) * I, 5, 5) - h = zeros(T, 5) - cones = Cone{T}[Cones.EpiNormInf{T, Complex{T}}(5)] + rt2 = sqrt(T(2)) + Random.seed!(1) + s = 3 + rand_mat = Hermitian(rand(T, s, s), :U) + dim = sum(1:s) + c = -Cones.smat_to_svec!(zeros(T, dim), rand_mat, rt2) + A = reshape(Cones.smat_to_svec!(zeros(T, dim), Matrix{T}(I, s, s), rt2), 1, dim) + b = T[1] + G = Diagonal(-one(T) * I, dim) + h = zeros(T, dim) + cones = Cone{T}[Cones.PosSemidefTri{T, T}(dim)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - @test r.primal_obj ≈ -3 atol=tol rtol=tol - @test r.x ≈ [2, 0, 2, 1, 0] atol=tol rtol=tol + eig_max = maximum(eigvals(rand_mat)) + @test r.primal_obj ≈ -eig_max atol=tol rtol=tol end -function epinorminf7(T; options...) +function possemideftri5(T; options...) tol = test_tol(T) - c = T[1, 0, 0, 0, 0, 0, 0] - A = zeros(T, 0, 7) - b = zeros(T, 0) - G = -one(T) * I - h = zeros(T, 7) - - for use_dual in (false, true) - cones = Cone{T}[Cones.EpiNormInf{T, Complex{T}}(7, use_dual = use_dual)] - - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test r.primal_obj ≈ 0 atol=tol rtol=tol - @test norm(r.x) ≈ 0 atol=tol rtol=tol - end -end - -function epinorminf8(T; options...) - tol = eps(T) ^ 0.2 - c = T[1, -1, 1, 1] - A = T[1 0 0 0 ; 0 1 0 0; 0 0 1 0] - b = T[-0.4, 0.3, -0.3] - G = vcat(zeros(T, 1, 4), Diagonal(T[-1, -1, -1, -1])) - h = T[1, 0, 0, 0, 0] - cones = Cone{T}[Cones.EpiNormInf{T, Complex{T}}(5, use_dual = true)] + Trt2 = sqrt(T(2)) + Trt2i = inv(Trt2) + c = T[1, 0, 0, 1] + A = T[0 0 1 0] + b = T[1] + G = Diagonal(-one(T) * I, 4) + h = zeros(T, 4) + cones = Cone{T}[Cones.PosSemidefTri{T, Complex{T}}(4)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - @test r.primal_obj ≈ -1.4 atol=tol rtol=tol - @test r.x ≈ [-0.4, 0.3, -0.3, -0.4] atol=tol rtol=tol - @test r.y ≈ [0, 0.25, -0.25] atol=tol rtol=tol - @test r.z ≈ [1.25, 1, -0.75, 0.75, 1] atol=tol rtol=tol + @test r.primal_obj ≈ Trt2 atol=tol rtol=tol + @test r.x ≈ [Trt2i, 0, 1, Trt2i] atol=tol rtol=tol end -function epinormeucl1(T; options...) +function possemideftri6(T; options...) tol = test_tol(T) - Trt2 = sqrt(T(2)) - Tirt2 = inv(Trt2) - c = T[0, -1, -1] - A = T[10 0 0; 0 10 0] - b = T[10, 10Tirt2] - G = Matrix{T}(-I, 3, 3) - h = zeros(T, 3) - cones = Cone{T}[Cones.EpiNormEucl{T}(3)] + rt2 = sqrt(T(2)) + Random.seed!(1) + c = T[1] + A = zeros(T, 0, 1) + b = T[] + rand_mat = Hermitian(rand(Complex{T}, 2, 2), :U) + G = reshape(T[-1, 0, 0, -1], 4, 1) + h = -Cones.smat_to_svec!(zeros(T, 4), rand_mat, rt2) + cones = Cone{T}[Cones.PosSemidefTri{T, Complex{T}}(4)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - @test r.primal_obj ≈ -Trt2 atol=tol rtol=tol - @test r.x ≈ [1, Tirt2, Tirt2] atol=tol rtol=tol - @test r.y ≈ [Trt2 / 10, 0] atol=tol rtol=tol + eig_max = maximum(eigvals(rand_mat)) + @test r.primal_obj ≈ eig_max atol=tol rtol=tol + @test r.x[1] ≈ eig_max atol=tol rtol=tol end -function epinormeucl2(T; options...) +function possemideftri7(T; options...) tol = test_tol(T) - c = T[0, -1, -1] - A = T[1 0 0] - b = T[0] - G = -one(T) * I - h = zeros(T, 3) - cones = Cone{T}[Cones.EpiNormEucl{T}(3)] + rt2 = sqrt(T(2)) + Random.seed!(1) + side = 3 + rand_mat = Hermitian(rand(Complex{T}, side, side), :U) + dim = abs2(side) + c = -Cones.smat_to_svec!(zeros(T, dim), rand_mat, rt2) + A = reshape(Cones.smat_to_svec!(zeros(T, dim), Matrix{Complex{T}}(I, side, side), rt2), 1, dim) + b = T[1] + G = Diagonal(-one(T) * I, dim) + h = zeros(T, dim) + cones = Cone{T}[Cones.PosSemidefTri{T, Complex{T}}(dim)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - @test r.primal_obj ≈ 0 atol=tol rtol=tol - @test norm(r.x) ≈ 0 atol=tol rtol=tol + eig_max = maximum(eigvals(rand_mat)) + @test r.primal_obj ≈ -eig_max atol=tol rtol=tol end -function epinormeucl3(T; options...) +function possemideftri8(T; options...) tol = test_tol(T) - c = T[1, 0, 0] - A = T[0 1 0] - b = T[1] - G = Diagonal(-one(T) * I, 3) - h = zeros(T, 3) - cones = Cone{T}[Cones.EpiNormEucl{T}(3)] + rt2 = sqrt(T(2)) + c = T[1] + A = zeros(T, 0, 1) + b = T[] + G = zeros(T, 15, 1) + G[[1, 3, 6, 10, 15]] .= -1 + h = zeros(T, 15) + @. h[[7, 8, 9, 11, 12, 13]] = rt2 * [1, 1, 0, 1, -1, 1] + cones = Cone{T}[Cones.PosSemidefTri{T, T}(15)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - @test r.primal_obj ≈ 1 atol=tol rtol=tol - @test r.x ≈ [1, 1, 0] atol=tol rtol=tol + rt3 = sqrt(T(3)) + @test r.primal_obj ≈ rt3 atol=tol rtol=tol + @test r.s ≈ [rt3, 0, rt3, 0, 0, rt3, rt2, rt2, 0, rt3, rt2, -rt2, rt2, 0, rt3] atol=tol rtol=tol + inv6 = inv(T(6)) + rt2inv6 = rt2 / 6 + invrt6 = inv(rt2 * rt3) + @test r.z ≈ [inv6, -rt2inv6, inv6, rt2inv6, -rt2inv6, inv6, 0, 0, 0, 0, -invrt6, invrt6, -invrt6, 0, inv(T(2))] atol=tol rtol=tol end -function epipersquare1(T; options...) +function possemideftri9(T; options...) tol = test_tol(T) - c = T[0, 0, -1, -1] - A = T[1 0 0 0; 0 1 0 0] - b = T[0.5, 1] - G = Matrix{T}(-I, 4, 4) - h = zeros(T, 4) - cones = Cone{T}[Cones.EpiPerSquare{T}(4)] + rt2 = sqrt(T(2)) + inv2 = inv(T(2)) + c = vcat(one(T), zeros(T, 9)) + A = zeros(T, 0, 10) + b = T[] + G = zeros(T, 16, 10) + G[1, 2] = G[1, 4] = G[1, 7] = G[1, 8] = G[1, 10] = inv2 + G[1, 1] = G[2, 2] = G[4, 4] = G[7, 7] = G[11, 8] = G[16, 10] = -1 + G[3, 3] = G[5, 5] = G[6, 6] = G[15, 9] = -rt2 + h = zeros(T, 16) + @. h[[8, 9, 10, 12, 13, 14]] = rt2 * [1, 1, 0, 1, -1, 1] + cones = Cone{T}[Cones.Nonnegative{T}(1), Cones.PosSemidefTri{T, T}(15)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - @test r.primal_obj ≈ -sqrt(T(2)) atol=tol rtol=tol - @test r.x[3:4] ≈ [1, 1] / sqrt(T(2)) atol=tol rtol=tol + rt3 = sqrt(T(3)) + @test r.primal_obj ≈ rt2 + rt3 atol=tol rtol=tol + invrt2 = inv(rt2) + invrt3 = inv(rt3) + @test r.s ≈ [0, invrt2 + invrt3, 1 - rt2 / rt3, invrt2 + invrt3, rt2 * invrt3, -rt2 * invrt3, invrt3, rt2, rt2, 0, rt2, rt2, -rt2, rt2, 0, rt3] atol=tol rtol=tol + invrt6 = invrt2 * invrt3 + @test r.z ≈ [1, inv2, 0, inv2, 0, 0, inv2, -inv2, -inv2, 0, inv2, -invrt6, invrt6, -invrt6, 0, inv2] atol=tol rtol=tol end -function epipersquare2(T; options...) +function doublynonnegativetri1(T; options...) tol = test_tol(T) - Tirt2 = inv(sqrt(T(2))) - c = T[0, 0, -1] - A = T[1 0 0; 0 1 0] - b = T[Tirt2 / 2, Tirt2] + n = q = 3 + c = T[0, 1, 0] + A = T[1 0 0; 0 0 1] + b = ones(T, 2) + G = SparseMatrixCSC(-one(T) * I, q, n) + for use_dual in [false, true] + use_dual = false + h = (use_dual ? T[1, 0, 1] : zeros(T, q)) + cones = Cone{T}[Cones.DoublyNonnegativeTri{T}(q, use_dual = use_dual)] + + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.Optimal + @test r.primal_obj ≈ 0 atol=tol rtol=tol + @test r.x ≈ T[1, 0, 1] atol=tol rtol=tol + @test r.s ≈ T[1, 0, 1] atol=tol rtol=tol + end +end + +function doublynonnegativetri2(T; options...) + tol = test_tol(T) + c = T[0, -1, 0] + A = T[1 0 0; 0 0 1] + b = T[1.0, 1.5] G = -one(T) * I - h = zeros(T, 3) - cones = Cone{T}[Cones.EpiPerSquare{T}(3)] + h = [-inv(T(2)), zero(T), -inv(T(2))] + cones = Cone{T}[Cones.DoublyNonnegativeTri{T}(3)] - r = build_solve_check(c, A, b, G, h, cones, tol; obj_offset = -one(T), options...) + r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - @test r.primal_obj ≈ -Tirt2 - 1 atol=tol rtol=tol - @test r.x[2] ≈ Tirt2 atol=tol rtol=tol + @test r.primal_obj ≈ -one(T) atol=tol rtol=tol + @test r.x[2] ≈ one(T) atol=tol rtol=tol end -function epipersquare3(T; options...) +function doublynonnegativetri3(T; options...) tol = test_tol(T) - c = T[0, 1, -1, -1] - A = T[1 0 0 0] + c = ones(T, 3) + A = T[1 0 1] b = T[0] - G = SparseMatrixCSC(-one(T) * I, 4, 4) - h = zeros(T, 4) - cones = Cone{T}[Cones.EpiPerSquare{T}(4)] + G = Diagonal(-one(T) * I, 3) + h = zeros(T, 3) + cones = Cone{T}[Cones.PosSemidefTri{T, T}(3)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal @@ -453,808 +525,487 @@ function epipersquare3(T; options...) @test norm(r.x) ≈ 0 atol=tol rtol=tol end -function epipersquare4(T; options...) +function possemideftrisparse1(T; options...) + if !(T <: LinearAlgebra.BlasReal) + return # only works with BLAS real types + end tol = test_tol(T) - c = zeros(T, 7) - c[1] = -1 - A = zeros(T, 0, 7) - b = zeros(T, 0) - G = sparse( - [1, 1, 2, 2, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11], - [1, 7, 2, 3, 4, 2, 3, 5, 4, 7, 6, 5, 6, 7], - T[1, -0.5, 1, 1, 1, -1, -1, -1, -1, -0.5, -1, -1, -1, -1], - 11, 7) - h = zeros(T, 11) - h[2] = 3 - cones = Cone{T}[Cones.Nonnegative{T}(2), Cones.EpiPerSquare{T}(3), Cones.EpiPerSquare{T}(3), Cones.EpiPerSquare{T}(3)] + c = T[0, -1, 0] + A = T[1 0 0; 0 0 1] + b = T[0.5, 1] + G = -one(T) * I + h = zeros(T, 3) + row_idxs = [1, 2, 2] + col_idxs = [1, 1, 2] + cones = Cone{T}[Cones.PosSemidefTriSparse{T, T}(2, row_idxs, col_idxs)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - @test r.primal_obj ≈ -1 atol=tol rtol=tol - rt2 = sqrt(T(2)) - @test r.x ≈ [1, 1, 1, 1, rt2, rt2, 2] atol=tol rtol=tol - @test r.s ≈ [0, 0, 1, 1, rt2, 1, 1, rt2, rt2, rt2, 2] atol=tol rtol=tol - inv3 = inv(T(3)) - rt2inv3 = inv3 * rt2 - @test r.z ≈ [1, inv3, inv3, inv3, -rt2inv3, inv3, inv3, -rt2inv3, rt2inv3, rt2inv3, -2 * inv3] atol=tol rtol=tol + @test r.primal_obj ≈ -one(T) atol=tol rtol=tol + @test r.x[2] ≈ one(T) atol=tol rtol=tol end -# TODO add use_dual = true tests -function epirelentropy1(T; options...) +function possemideftrisparse2(T; options...) + if !(T <: LinearAlgebra.BlasReal) + return # only works with BLAS real types + end tol = test_tol(T) - Random.seed!(1) - for w_dim in [1, 2, 3] - dim = 1 + 2 * w_dim - c = T[1] - A = zeros(T, 0, 1) - b = zeros(T, 0) - G = zeros(T, dim, 1) - G[1, 1] = -1 - h = zeros(T, dim) - h[2:2:(end - 1)] .= 1 - w = rand(T, w_dim) .+ 1 - h[3:2:end] .= w - cones = Cone{T}[Cones.EpiRelEntropy{T}(dim)] + Trt2 = sqrt(T(2)) + Trt2i = inv(Trt2) + c = T[1, 0, 0, 1] + A = T[0 0 1 0] + b = T[1] + G = Diagonal(-one(T) * I, 4) + h = zeros(T, 4) + row_idxs = [1, 2, 2] + col_idxs = [1, 1, 2] + cones = Cone{T}[Cones.PosSemidefTriSparse{T, Complex{T}}(2, row_idxs, col_idxs)] - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test r.primal_obj ≈ sum(wi * log(wi) for wi in w) atol=tol rtol=tol - end + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.Optimal + @test r.primal_obj ≈ Trt2 atol=tol rtol=tol + @test r.x ≈ [Trt2i, 0, 1, Trt2i] atol=tol rtol=tol end -function epirelentropy2(T; options...) +function possemideftrisparse3(T; options...) + if !(T <: LinearAlgebra.BlasReal) + return # only works with BLAS real types + end tol = test_tol(T) - for w_dim in [1, 2, 4] - dim = 1 + 2 * w_dim - c = fill(-one(T), w_dim) - A = zeros(T, 0, w_dim) - b = zeros(T, 0) - G = zeros(T, dim, w_dim) - for i in 1:w_dim - G[2i + 1, i] = -1 + rt2 = sqrt(T(2)) + Random.seed!(1) + for is_complex in (false, true), side in [1, 2, 5, 20] + R = (is_complex ? Complex{T} : T) + rand_mat_L = tril!(sprand(R, side, side, inv(sqrt(side))) + I) + (row_idxs, col_idxs, vals) = findnz(rand_mat_L) + dim = (is_complex ? side + 2 * (length(row_idxs) - side) : length(row_idxs)) + c = zeros(T, dim) + A = zeros(T, 1, dim) + idx = 1 + for (i, v) in enumerate(vals) # scale + if row_idxs[i] == col_idxs[i] + c[idx] = -real(v) + A[idx] = 1 + else + c[idx] = -rt2 * real(v) + if is_complex + idx += 1 + c[idx] = -rt2 * imag(v) + end + end + idx += 1 end + b = T[1] + G = Diagonal(-one(T) * I, dim) h = zeros(T, dim) - h[2:2:(dim - 1)] .= 1 - cones = Cone{T}[Cones.EpiRelEntropy{T}(dim)] + cones = Cone{T}[Cones.PosSemidefTriSparse{T, R}(side, row_idxs, col_idxs, use_dual = true)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - @test r.primal_obj ≈ -w_dim atol=tol rtol=tol - @test r.x ≈ fill(1, w_dim) atol=tol rtol=tol + eig_max = maximum(eigvals(Hermitian(Matrix(rand_mat_L), :L))) + @test r.primal_obj ≈ -eig_max atol=tol rtol=tol end end -function epirelentropy3(T; options...) - tol = test_tol(T) - for w_dim in [2, 4] - dim = 1 + 2 * w_dim - c = fill(-one(T), w_dim) - A = ones(T, 1, w_dim) - b = T[dim] - G = zeros(T, dim, w_dim) - for i in 1:w_dim - G[2i, i] = -1 - end - h = zeros(T, dim) - h[3:2:end] .= 1 - cones = Cone{T}[Cones.EpiRelEntropy{T}(dim)] - - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test r.primal_obj ≈ -dim atol=tol rtol=tol - @test r.x ≈ fill(dim / w_dim, w_dim) atol=tol rtol=tol - end -end - -function epirelentropy4(T; options...) +function possemideftrisparse4(T; options...) + if !(T <: LinearAlgebra.BlasReal) + return # only works with BLAS real types + end tol = test_tol(T) + rt2 = sqrt(T(2)) c = T[1] A = zeros(T, 0, 1) - b = zeros(T, 0) - G = Matrix{T}(-I, 5, 1) - h = T[0, 1, 2, 5, 3] - cones = Cone{T}[Cones.EpiRelEntropy{T}(5)] + b = T[] + G = zeros(T, 10, 1) + G[[1, 2, 3, 6, 10]] .= -1 + h = zeros(T, 10) + @. h[[4, 5, 7, 8, 9]] = rt2 * [1, 1, 1, -1, 1] + row_idxs = [1, 2, 3, 4, 4, 4, 5, 5, 5, 5] + col_idxs = [1, 2, 3, 1, 2, 4, 1, 2, 3, 5] + cones = Cone{T}[Cones.PosSemidefTriSparse{T, T}(5, row_idxs, col_idxs)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - entr = 2 * log(T(2)) + 3 * log(3 / T(5)) - @test r.primal_obj ≈ entr atol=tol rtol=tol - @test r.s ≈ [entr, 1, 2, 5, 3] atol=tol rtol=tol - @test r.z ≈ [1, 2, log(inv(T(2))) - 1, 3 / T(5), log(5 / T(3)) - 1] atol=tol rtol=tol + rt3 = sqrt(T(3)) + @test r.primal_obj ≈ rt3 atol=tol rtol=tol + @test r.s ≈ [rt3, rt3, rt3, rt2, rt2, rt3, rt2, -rt2, rt2, rt3] atol=tol rtol=tol + inv6 = inv(T(6)) + rt2inv6 = rt2 / 6 + invrt6 = inv(rt2 * rt3) + @test r.z ≈ [inv6, inv6, inv6, 0, 0, 0, -invrt6, invrt6, -invrt6, inv(T(2))] atol=tol rtol=tol end -function epirelentropy5(T; options...) +function possemideftrisparse5(T; options...) + if !(T <: LinearAlgebra.BlasReal) + return # only works with BLAS real types + end tol = test_tol(T) - c = T[0, -1] - A = zeros(T, 0, 2) - b = zeros(T, 0) - G = vcat(zeros(T, 1, 2), repeat(T[0 0; -1 -1], 3), [-1, 0]') - h = T[0, 1, 0, 1, 0, 1, 0, 0] - cones = Cone{T}[Cones.EpiRelEntropy{T}(7), Cones.Nonnegative{T}(1)] + rt2 = sqrt(T(2)) + inv2 = inv(T(2)) + c = vcat(one(T), zeros(T, 5)) + A = zeros(T, 0, 6) + b = T[] + G = vcat(Matrix(-one(T) * I, 6, 6), zeros(T, 6, 6)) + G[1, 2:end] .= inv2 + h = vcat(zeros(T, 6), rt2 * [1, 1, 0, 1, -1, 1]) + row_idxs = [1, 2, 3, 4, 5, 4, 4, 4, 5, 5, 5] + col_idxs = [1, 2, 3, 4, 5, 1, 2, 3, 1, 2, 3] + cones = Cone{T}[Cones.Nonnegative{T}(1), Cones.PosSemidefTriSparse{T, T}(5, row_idxs, col_idxs, use_dual = true)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - @test r.primal_obj ≈ -1 atol=tol rtol=tol - @test r.s ≈ [0, 1, 1, 1, 1, 1, 1, 0] atol=tol rtol=tol - @test r.z ≈ inv(T(3)) * [1, 1, -1, 1, -1, 1, -1, 3] atol=tol rtol=tol + rt3 = sqrt(T(3)) + @test r.primal_obj ≈ rt2 + rt3 atol=tol rtol=tol + invrt2 = inv(rt2) + invrt3 = inv(rt3) + @test r.s ≈ [0, invrt2 + invrt3, invrt2 + invrt3, invrt3, rt2, rt3, rt2, rt2, 0, rt2, -rt2, rt2] atol=tol rtol=tol + invrt6 = invrt2 * invrt3 + @test r.z ≈ [1, inv2, inv2, inv2, inv2, inv2, -inv2, -inv2, 0, -invrt6, invrt6, -invrt6] atol=tol rtol=tol end -# TODO add use_dual = true tests -function epiperentropy1(T; options...) +function linmatrixineq1(T; options...) tol = test_tol(T) Random.seed!(1) - for w_dim in [1, 2, 3] - dim = 2 + w_dim + for side in [2, 4], R in [T, Complex{T}] c = T[1] A = zeros(T, 0, 1) - b = zeros(T, 0) - G = zeros(T, dim, 1) + b = T[] + G = zeros(T, 2, 1) G[1, 1] = -1 - h = zeros(T, dim) - h[2] = 1 - w = rand(T, w_dim) .+ 1 - h[3:end] .= w - cones = Cone{T}[Cones.EpiPerEntropy{T}(dim)] + h = T[0, 2] + A_1_half = rand(R, side, side) + A_1 = Hermitian(A_1_half * A_1_half' + 2I) + F = eigen(A_1) + val_1 = F.values[end] + vec_1 = F.vectors[:, end] + As = [A_1, Hermitian(-vec_1 * vec_1')] + cones = Cone{T}[Cones.LinMatrixIneq{T}(As)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - @test r.primal_obj ≈ sum(wi * log(wi) for wi in w) atol=tol rtol=tol + @test r.primal_obj ≈ 2 / val_1 atol=tol rtol=tol + @test r.s ≈ [2 / val_1, 2] atol=tol rtol=tol end end -function epiperentropy2(T; options...) +function linmatrixineq2(T; options...) tol = test_tol(T) - for w_dim in [1, 2, 4] - dim = 2 + w_dim - c = fill(-one(T), w_dim) - A = zeros(T, 0, w_dim) - b = zeros(T, 0) - G = vcat(zeros(T, 2, w_dim), Matrix{T}(-I, w_dim, w_dim)) + Random.seed!(1) + for Rs in [[T, T], [Complex{T}, Complex{T}], [T, Complex{T}, T], [Complex{T}, T, T]] + dim = length(Rs) + c = ones(T, dim - 1) + A = zeros(T, 0, dim - 1) + b = T[] + G = vcat(spzeros(T, 1, dim - 1), sparse(-one(T) * I, dim - 1, dim - 1)) h = zeros(T, dim) - h[2] = 1 - cones = Cone{T}[Cones.EpiPerEntropy{T}(dim)] + h[1] = 1 + As = Hermitian[] + for R in Rs + A_half = rand(R, 3, 3) + push!(As, Hermitian(A_half * A_half')) + end + As[1] += I + cones = Cone{T}[Cones.LinMatrixIneq{T}(As)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - @test r.primal_obj ≈ -w_dim atol=tol rtol=tol - @test r.x ≈ fill(1, w_dim) atol=tol rtol=tol + @test r.primal_obj < 0 end end -function epiperentropy3(T; options...) - tol = test_tol(T) - for w_dim in [2, 4] - dim = 2 + w_dim - c = T[-1] - A = ones(T, 1, 1) - b = T[dim] - G = zeros(T, dim, 1) - G[2, 1] = -1 - h = zeros(T, dim) - h[3:end] .= 1 - cones = Cone{T}[Cones.EpiPerEntropy{T}(dim)] +function linmatrixineq3(T; options...) + dense1 = [1 0; 0 1] + dense2 = [1 0; 0 -1] + sparse1 = sparse(dense1) + sparse2 = sparse(dense2) + diag1 = Diagonal([1, 1]) + diag2 = Diagonal([1, -1]) + # NOTE not all combinations work due to missing methods in LinearAlgebra + As_list = [ + [dense1, dense2], + # [dense1, sparse2], + [dense1, diag2], + # [sparse1, dense2], + [sparse1, sparse2], + # [sparse1, diag2], + [diag1, dense2], + # [diag1, sparse2], + [diag1, diag2], + [I, dense2], + # [I, sparse2], + [I, diag2], + ] + for As in As_list + tol = test_tol(T) + c = T[1] + A = zeros(T, 0, 1) + b = T[] + G = zeros(T, 2, 1) + G[1, 1] = -1 + h = T[0, -1] + cones = Cone{T}[Cones.LinMatrixIneq{T}(As)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - @test r.primal_obj ≈ -dim atol=tol rtol=tol - @test r.x ≈ [dim] atol=tol rtol=tol + @test r.primal_obj ≈ 1 atol=tol rtol=tol + @test r.s ≈ [1, -1] atol=tol rtol=tol end end -function epiperentropy4(T; options...) +function epinorminf1(T; options...) tol = test_tol(T) - c = T[1] - A = zeros(T, 0, 1) - b = zeros(T, 0) - G = Matrix{T}(-I, 4, 1) - h = T[0, 5, 2, 3] - cones = Cone{T}[Cones.EpiPerEntropy{T}(4)] + Tirt2 = inv(sqrt(T(2))) + c = T[0, -1, -1] + A = T[1 0 0; 0 1 0] + b = [one(T), Tirt2] + G = SparseMatrixCSC(-one(T) * I, 3, 3) + h = zeros(T, 3) + cones = Cone{T}[Cones.EpiNormInf{T, T}(3)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - entr = 2 * log(2 / T(5)) + 3 * log(3 / T(5)) - @test r.primal_obj ≈ entr atol=tol rtol=tol - @test r.s ≈ [entr, 5, 2, 3] atol=tol rtol=tol - @test r.z ≈ [1, 1, log(5 / T(2)) - 1, log(5 / T(3)) - 1] atol=tol rtol=tol + @test r.primal_obj ≈ -1 - Tirt2 atol=tol rtol=tol + @test r.x ≈ [1, Tirt2, 1] atol=tol rtol=tol + @test r.y ≈ [1, 1] atol=tol rtol=tol end -function epiperentropy5(T; options...) +function epinorminf2(T; options...) + tol = 10 * test_tol(T) + l = 3 + L = 2l + 1 + c = collect(T, -l:l) + A = spzeros(T, 2, L) + A[1, 1] = A[1, L] = A[2, 1] = 1; A[2, L] = -1 + b = T[0, 0] + G = [spzeros(T, 1, L); sparse(one(T) * I, L, L); spzeros(T, 1, L); sparse(T(2) * I, L, L)] + h = zeros(T, 2L + 2); h[1] = 1; h[L + 2] = 1 + cones = Cone{T}[Cones.EpiNormInf{T, T}(L + 1, use_dual = true), Cones.EpiNormInf{T, T}(L + 1, use_dual = false)] + + r = build_solve_check(c, A, b, G, h, cones, tol; obj_offset = one(T), options...) + @test r.status == Solvers.Optimal + @test r.primal_obj ≈ -l + 2 atol=tol rtol=tol + @test r.x[2] ≈ 0.5 atol=tol rtol=tol + @test r.x[end - 1] ≈ -0.5 atol=tol rtol=tol + @test sum(abs, r.x) ≈ 1 atol=tol rtol=tol +end + +function epinorminf3(T; options...) tol = test_tol(T) - c = T[0, -1] - A = zeros(T, 0, 2) + c = T[1, 0, 0, 0, 0, 0] + A = zeros(T, 0, 6) b = zeros(T, 0) - G = vcat(zeros(T, 2, 2), fill(-1, 3, 2), [-1, 0]') - h = T[0, 1, 0, 0, 0, 0] - cones = Cone{T}[Cones.EpiPerEntropy{T}(5), Cones.Nonnegative{T}(1)] + G = Diagonal(-one(T) * I, 6) + h = zeros(T, 6) + + for use_dual in (false, true) + cones = Cone{T}[Cones.EpiNormInf{T, T}(6, use_dual = use_dual)] + + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.Optimal + @test r.primal_obj ≈ 0 atol=tol rtol=tol + @test norm(r.x) ≈ 0 atol=tol rtol=tol + end +end + +function epinorminf4(T; options...) + tol = test_tol(T) + c = T[0, 1, -1] + A = T[1 0 0; 0 1 0] + b = T[1, -0.4] + G = -one(T) * I + h = zeros(T, 3) + cones = Cone{T}[Cones.EpiNormInf{T, T}(3, use_dual = true)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal @test r.primal_obj ≈ -1 atol=tol rtol=tol - @test r.s ≈ [0, 1, 1, 1, 1, 0] atol=tol rtol=tol - @test r.z ≈ inv(T(3)) * [1, 3, -1, -1, -1, 3] atol=tol rtol=tol + @test r.x ≈ [1, -0.4, 0.6] atol=tol rtol=tol + @test r.y ≈ [1, 0] atol=tol rtol=tol end -function epitracerelentropytri1(T; options...) +function epinorminf5(T; options...) tol = test_tol(T) Random.seed!(1) - rt2 = sqrt(T(2)) - side = 4 - svec_dim = Cones.svec_length(side) - cone_dim = 2 * svec_dim + 1 - c = T[1] - A = zeros(T, 0, 1) - b = T[] - W = rand(T, side, side) - W = W * W' - V = rand(T, side, side) - V = V * V' - h = vcat(zero(T), Cones.smat_to_svec!(zeros(T, svec_dim), V, rt2), Cones.smat_to_svec!(zeros(T, svec_dim), W, rt2)) - G = zeros(T, cone_dim, 1) - G[1, 1] = -1 - cones = Cone{T}[Cones.EpiTraceRelEntropyTri{T}(cone_dim)] + c = T[1, 0, 0, 0, 0, 0] + A = rand(T(-9):T(9), 6, 6) + b = vec(sum(A, dims = 2)) + G = rand(T, 6, 6) + h = vec(sum(G, dims = 2)) + cones = Cone{T}[Cones.EpiNormInf{T, T}(6, use_dual = true)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - @test r.primal_obj ≈ tr(W * log(W) - W * log(V)) atol=tol rtol=tol + @test r.primal_obj ≈ 1 atol=tol rtol=tol end -function epitracerelentropytri2(T; options...) +function epinorminf6(T; options...) tol = test_tol(T) - rt2 = sqrt(T(2)) - side = 3 - svec_dim = Cones.svec_length(side) - cone_dim = 2 * svec_dim + 1 - c = vcat(zeros(T, svec_dim), -ones(T, svec_dim)) - A = hcat(Matrix{T}(I, svec_dim, svec_dim), zeros(T, svec_dim, svec_dim)) - b = zeros(T, svec_dim) - b[[sum(1:i) for i in 1:side]] .= 1 - h = vcat(T(5), zeros(T, 2 * svec_dim)) - G = vcat(zeros(T, 1, 2 * svec_dim), ModelUtilities.vec_to_svec!(Diagonal(-one(T) * I, 2 * svec_dim))) - cones = Cone{T}[Cones.EpiTraceRelEntropyTri{T}(cone_dim)] + c = T[0, -1, -1, -1, -1] + A = T[1 0 0 0 0; 0 1 0 0 0; 0 0 0 1 0; 0 0 0 0 1] + b = T[2, 0, 1, 0] + G = SparseMatrixCSC(-one(T) * I, 5, 5) + h = zeros(T, 5) + cones = Cone{T}[Cones.EpiNormInf{T, Complex{T}}(5)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - W = Hermitian(Cones.svec_to_smat!(zeros(T, side, side), r.s[(svec_dim + 2):end], rt2), :U) - @test tr(W * log(W)) ≈ T(5) atol=tol rtol=tol + @test r.primal_obj ≈ -3 atol=tol rtol=tol + @test r.x ≈ [2, 0, 2, 1, 0] atol=tol rtol=tol end -function epitracerelentropytri3(T; options...) +function epinorminf7(T; options...) tol = test_tol(T) - rt2 = sqrt(T(2)) - side = 3 - svec_dim = Cones.svec_length(side) - cone_dim = 2 * svec_dim + 1 - c = vcat(zeros(T, svec_dim + 1), ones(T, svec_dim)) - A = hcat(one(T), zeros(T, 1, 2 * svec_dim)) - b = [zero(T)] - h = zeros(T, cone_dim) - G = Diagonal(-one(T) * I, cone_dim) - cones = Cone{T}[Cones.EpiTraceRelEntropyTri{T}(cone_dim)] + c = T[1, 0, 0, 0, 0, 0, 0] + A = zeros(T, 0, 7) + b = zeros(T, 0) + G = -one(T) * I + h = zeros(T, 7) - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test r.primal_obj ≈ zero(T) atol=tol rtol=tol - @test r.s[1] ≈ zero(T) atol=tol rtol=tol - @test r.s[(svec_dim + 2):end] ≈ zeros(T, svec_dim) atol=tol rtol=tol + for use_dual in (false, true) + cones = Cone{T}[Cones.EpiNormInf{T, Complex{T}}(7, use_dual = use_dual)] + + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.Optimal + @test r.primal_obj ≈ 0 atol=tol rtol=tol + @test norm(r.x) ≈ 0 atol=tol rtol=tol + end end -function epitracerelentropytri4(T; options...) - tol = test_tol(T) - Random.seed!(1) - rt2 = sqrt(T(2)) - side = 3 - svec_dim = Cones.svec_length(side) - cone_dim = 2 * svec_dim + 1 - c = vcat(zero(T), ones(T, svec_dim), zeros(T, svec_dim)) - A = hcat(one(T), zeros(T, 1, 2 * svec_dim)) - b = [zero(T)] - h = zeros(T, cone_dim) - G = Diagonal(-one(T) * I, cone_dim) - cones = Cone{T}[Cones.EpiTraceRelEntropyTri{T}(cone_dim)] +function epinorminf8(T; options...) + tol = eps(T) ^ 0.2 + c = T[1, -1, 1, 1] + A = T[1 0 0 0 ; 0 1 0 0; 0 0 1 0] + b = T[-0.4, 0.3, -0.3] + G = vcat(zeros(T, 1, 4), Diagonal(T[-1, -1, -1, -1])) + h = T[1, 0, 0, 0, 0] + cones = Cone{T}[Cones.EpiNormInf{T, Complex{T}}(5, use_dual = true)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - @test r.primal_obj ≈ zero(T) atol=tol rtol=tol - @test r.s ≈ zeros(T, cone_dim) atol=tol rtol=tol + @test r.primal_obj ≈ -1.4 atol=tol rtol=tol + @test r.x ≈ [-0.4, 0.3, -0.3, -0.4] atol=tol rtol=tol + @test r.y ≈ [0, 0.25, -0.25] atol=tol rtol=tol + @test r.z ≈ [1.25, 1, -0.75, 0.75, 1] atol=tol rtol=tol end -function epipertraceentropytri1(T; options...) +function epinormeucl1(T; options...) tol = test_tol(T) - Random.seed!(1) - rt2 = sqrt(T(2)) - side = 4 - svec_dim = Cones.svec_length(side) - cone_dim = svec_dim + 2 - c = T[1] - A = zeros(T, 0, 1) - b = T[] - W = rand(T, side, side) - W = W * W' - v = rand(T) - h = vcat(zero(T), v, Cones.smat_to_svec!(zeros(T, svec_dim), W, rt2)) - G = zeros(T, cone_dim, 1) - G[1, 1] = -1 - cones = Cone{T}[Cones.EpiPerTraceEntropyTri{T}(cone_dim)] + Trt2 = sqrt(T(2)) + Tirt2 = inv(Trt2) + c = T[0, -1, -1] + A = T[10 0 0; 0 10 0] + b = T[10, 10Tirt2] + G = Matrix{T}(-I, 3, 3) + h = zeros(T, 3) + cones = Cone{T}[Cones.EpiNormEucl{T}(3)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - # TODO need https://github.com/JuliaLinearAlgebra/GenericLinearAlgebra.jl/issues/51 to use log with BF - (vals_W, vecs_W) = eigen(Hermitian(W, :U)) - log_W = vecs_W * Diagonal(log.(vals_W)) * vecs_W' - @test r.primal_obj ≈ tr(W * log_W - W * log(v)) atol=tol rtol=tol + @test r.primal_obj ≈ -Trt2 atol=tol rtol=tol + @test r.x ≈ [1, Tirt2, Tirt2] atol=tol rtol=tol + @test r.y ≈ [Trt2 / 10, 0] atol=tol rtol=tol end -function epipertraceentropytri2(T; options...) +function epinormeucl2(T; options...) tol = test_tol(T) - rt2 = sqrt(T(2)) - side = 3 - svec_dim = Cones.svec_length(side) - cone_dim = svec_dim + 2 - c = vcat(zero(T), -ones(T, svec_dim)) - A = hcat(one(T), zeros(T, 1, svec_dim)) - b = T[1] - h = vcat(T(5), zeros(T, svec_dim + 1)) - G = vcat(zeros(T, 1, svec_dim + 1), ModelUtilities.vec_to_svec!(Diagonal(-one(T) * I, svec_dim + 1))) - cones = Cone{T}[Cones.EpiPerTraceEntropyTri{T}(cone_dim)] + c = T[0, -1, -1] + A = T[1 0 0] + b = T[0] + G = -one(T) * I + h = zeros(T, 3) + cones = Cone{T}[Cones.EpiNormEucl{T}(3)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - W = Hermitian(Cones.svec_to_smat!(zeros(T, side, side), r.s[3:end], rt2), :U) - @test r.s[2] ≈ 1 atol=tol rtol=tol - @test tr(W * log(W)) ≈ T(5) atol=tol rtol=tol + @test r.primal_obj ≈ 0 atol=tol rtol=tol + @test norm(r.x) ≈ 0 atol=tol rtol=tol end -function epipertraceentropytri3(T; options...) +function epinormeucl3(T; options...) tol = test_tol(T) - rt2 = sqrt(T(2)) - side = 3 - svec_dim = Cones.svec_length(side) - cone_dim = svec_dim + 2 - c = vcat(zeros(T, 2), ones(T, svec_dim)) - A = hcat(one(T), zeros(T, 1, svec_dim + 1)) - b = [zero(T)] - h = zeros(T, cone_dim) - G = Diagonal(-one(T) * I, cone_dim) - cones = Cone{T}[Cones.EpiPerTraceEntropyTri{T}(cone_dim)] + c = T[1, 0, 0] + A = T[0 1 0] + b = T[1] + G = Diagonal(-one(T) * I, 3) + h = zeros(T, 3) + cones = Cone{T}[Cones.EpiNormEucl{T}(3)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - @test r.primal_obj ≈ zero(T) atol=tol rtol=tol - @test r.s[1] ≈ zero(T) atol=tol rtol=tol - @test r.s[3:end] ≈ zeros(T, svec_dim) atol=tol rtol=tol + @test r.primal_obj ≈ 1 atol=tol rtol=tol + @test r.x ≈ [1, 1, 0] atol=tol rtol=tol end -function epipertraceentropytri4(T; options...) +function epipersquare1(T; options...) tol = test_tol(T) - rt2 = sqrt(T(2)) - side = 3 - svec_dim = Cones.svec_length(side) - cone_dim = svec_dim + 2 - c = vcat(zero(T), one(T), zeros(T, svec_dim)) - A = hcat(one(T), zeros(T, 1, svec_dim + 1)) - b = [zero(T)] - h = zeros(T, cone_dim) - G = Diagonal(-one(T) * I, cone_dim) - cones = Cone{T}[Cones.EpiPerTraceEntropyTri{T}(cone_dim)] + c = T[0, 0, -1, -1] + A = T[1 0 0 0; 0 1 0 0] + b = T[0.5, 1] + G = Matrix{T}(-I, 4, 4) + h = zeros(T, 4) + cones = Cone{T}[Cones.EpiPerSquare{T}(4)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - @test r.primal_obj ≈ zero(T) atol=tol rtol=tol - @test r.s ≈ zeros(T, cone_dim) atol=tol rtol=tol + @test r.primal_obj ≈ -sqrt(T(2)) atol=tol rtol=tol + @test r.x[3:4] ≈ [1, 1] / sqrt(T(2)) atol=tol rtol=tol end -function hypoperlog1(T; options...) +function epipersquare2(T; options...) tol = test_tol(T) - Texph = exp(T(0.5)) - c = T[1, 1, 1] - A = T[0 1 0; 1 0 0] - b = T[2, 1] + Tirt2 = inv(sqrt(T(2))) + c = T[0, 0, -1] + A = T[1 0 0; 0 1 0] + b = T[Tirt2 / 2, Tirt2] G = -one(T) * I h = zeros(T, 3) - cones = Cone{T}[Cones.HypoPerLog{T}(3)] + cones = Cone{T}[Cones.EpiPerSquare{T}(3)] - r = build_solve_check(c, A, b, G, h, cones, tol; options...) + r = build_solve_check(c, A, b, G, h, cones, tol; obj_offset = -one(T), options...) @test r.status == Solvers.Optimal - @test r.primal_obj ≈ 2 * Texph + 3 atol=tol rtol=tol - @test r.x ≈ [1, 2, 2 * Texph] atol=tol rtol=tol - @test r.y ≈ -[1 + Texph / 2, 1 + Texph] atol=tol rtol=tol - @test r.z ≈ c + A' * r.y atol=tol rtol=tol + @test r.primal_obj ≈ -Tirt2 - 1 atol=tol rtol=tol + @test r.x[2] ≈ Tirt2 atol=tol rtol=tol end -function hypoperlog2(T; options...) +function epipersquare3(T; options...) tol = test_tol(T) - c = T[-1, 0, 0] - A = T[0 1 0] + c = T[0, 1, -1, -1] + A = T[1 0 0 0] b = T[0] - G = Diagonal(-one(T) * I, 3) - h = zeros(T, 3) - cones = Cone{T}[Cones.HypoPerLog{T}(3)] + G = SparseMatrixCSC(-one(T) * I, 4, 4) + h = zeros(T, 4) + cones = Cone{T}[Cones.EpiPerSquare{T}(4)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal @test r.primal_obj ≈ 0 atol=tol rtol=tol + @test norm(r.x) ≈ 0 atol=tol rtol=tol end -function hypoperlog3(T; options...) +function epipersquare4(T; options...) tol = test_tol(T) - c = T[1, 1, 1] - A = zeros(T, 0, 3) + c = zeros(T, 7) + c[1] = -1 + A = zeros(T, 0, 7) b = zeros(T, 0) - G = sparse([1, 2, 3, 4], [1, 2, 3, 1], -ones(T, 4)) - h = zeros(T, 4) - cones = Cone{T}[Cones.HypoPerLog{T}(3), Cones.Nonnegative{T}(1)] - - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test r.primal_obj ≈ 0 atol=tol rtol=tol - @test norm(r.x) ≈ 0 atol=tol rtol=tol - @test isempty(r.y) -end - -function hypoperlog4(T; options...) - tol = test_tol(T) - Texp2 = exp(T(-2)) - c = T[0, 0, 1] - A = T[0 1 0; 1 0 0] - b = T[1, -1] - G = SparseMatrixCSC(-one(T) * I, 3, 3) - h = zeros(T, 3) - cones = Cone{T}[Cones.HypoPerLog{T}(3, use_dual = true)] - - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test r.primal_obj ≈ Texp2 atol=tol rtol=tol - @test r.x ≈ [-1, 1, Texp2] atol=tol rtol=tol -end - -function hypoperlog5(T; options...) - tol = test_tol(T) - Tlogq = log(T(0.25)) - c = T[-1, 0, 0] - A = T[0 1 1] - b = T[1] - G = sparse([1, 3, 4], [1, 2, 3], -ones(T, 3)) - h = T[0, 1, 0, 0] - cones = Cone{T}[Cones.HypoPerLog{T}(4)] - - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test r.primal_obj ≈ -Tlogq atol=tol rtol=tol - @test r.x ≈ [Tlogq, 0.5, 0.5] atol=tol rtol=tol - @test r.y ≈ [2] atol=tol rtol=tol -end - -function hypoperlog6(T; options...) - tol = test_tol(T) - c = T[-1, 0, 0] - A = zeros(T, 0, 3) - b = zeros(T, 0) - G = sparse([1, 3, 4], [1, 2, 3], -ones(T, 3)) - h = zeros(T, 4) - cones = Cone{T}[Cones.HypoPerLog{T}(4)] - - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test r.primal_obj ≈ 0 atol=tol rtol=tol - @test r.x[1] ≈ 0 atol=tol rtol=tol - @test isempty(r.y) -end - -function hypoperlog7(T; options...) - tol = test_tol(T) - c = zeros(T, 4) - c[1] = -2 - A = zeros(T, 0, 4) - b = zeros(T, 0) - G = sparse( - [1, 2, 2, 3, 4, 5, 6], - [2, 1, 3, 4, 4, 3, 2], - T[1, 1, -1, -1, -1, -1, -1], - 6, 4) - h = zeros(T, 6) - h[1] = 2 - cones = Cone{T}[Cones.Nonnegative{T}(3), Cones.HypoPerLog{T}(3)] - - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test r.primal_obj ≈ -4 atol=tol rtol=tol - @test r.x ≈ [2, 2, 2, 0] atol=tol rtol=tol - @test r.s ≈ [0, 0, 0, 0, 2, 2] atol=tol rtol=tol - @test r.z ≈ [2, 2, 2, -2, -2, 2] atol=tol rtol=tol -end - -function hypogeomean1(T; options...) - tol = test_tol(T) - c = T[-1, 0, 0] - A = T[0 0 1; 0 1 0] - b = T[0.5, 1] - G = -one(T) * I - h = zeros(T, 3) - - for use_dual in (false, true) - cones = Cone{T}[Cones.HypoGeoMean{T}(3, use_dual = use_dual)] - - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test r.primal_obj ≈ (use_dual ? 0 : -inv(sqrt(T(2)))) atol=tol rtol=tol - @test r.x[2:3] ≈ [1, 0.5] atol=tol rtol=tol - end -end - -function hypogeomean2(T; options...) - tol = test_tol(T) - l = 4 - c = vcat(zero(T), ones(T, l)) - A = T[one(T) zeros(T, 1, l)] - G = SparseMatrixCSC(-one(T) * I, l + 1, l + 1) - h = zeros(T, l + 1) - - for use_dual in (false, true) - b = use_dual ? [-one(T)] : [one(T)] - cones = Cone{T}[Cones.HypoGeoMean{T}(l + 1, use_dual = use_dual)] - - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test r.primal_obj ≈ (use_dual ? 1 : l) atol=tol rtol=tol - @test r.x[2:end] ≈ (use_dual ? fill(inv(T(l)), l) : ones(l)) atol=tol rtol=tol - end -end - -function hypogeomean3(T; options...) - tol = test_tol(T) - l = 4 - c = ones(T, l) - A = zeros(T, 0, l) - b = zeros(T, 0) - G = [zeros(T, 1, l); Matrix{T}(-I, l, l)] - h = zeros(T, l + 1) - - for use_dual in (false, true) - cones = Cone{T}[Cones.HypoGeoMean{T}(l + 1, use_dual = use_dual)] - - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test r.primal_obj ≈ 0 atol=tol rtol=tol - @test norm(r.x) ≈ 0 atol=tol rtol=tol - end -end - -function hypogeomean4(T; options...) - tol = test_tol(T) - c = T[-1, 0, 0, 0] - A = zeros(T, 0, 4) - b = zeros(T, 0) - G = vcat(Matrix{T}(-I, 4, 4), T[0, 1, 1, 1]') - h = T[0, 0, 0, 0, 3] - cones = Cone{T}[Cones.HypoGeoMean{T}(4), Cones.Nonnegative{T}(1)] - - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test r.primal_obj ≈ -1 atol=tol rtol=tol - @test r.x ≈ [1, 1, 1, 1] atol=tol rtol=tol - @test r.s ≈ [1, 1, 1, 1, 0] atol=tol rtol=tol - @test r.z ≈ vcat(-1, fill(inv(T(3)), 4)) atol=tol rtol=tol -end - -function hypogeomean5(T; options...) - tol = test_tol(T) - c = T[-2, 0] - A = zeros(T, 0, 2) - b = zeros(T, 0) - G = sparse([1, 2, 3], [1, 2, 2], T[-1, -1, 1], 3, 2) - h = T[0, 0, 2] - cones = Cone{T}[Cones.HypoGeoMean{T}(2), Cones.Nonnegative{T}(1)] - - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test r.primal_obj ≈ -4 atol=tol rtol=tol - @test r.x ≈ [2, 2] atol=tol rtol=tol - @test r.s ≈ [2, 2, 0] atol=tol rtol=tol - @test r.z ≈ [-2, 2, 2] atol=tol rtol=tol -end - -function hypogeomean6(T; options...) - tol = test_tol(T) - c = zeros(T, 10) - c[1] = -1 - A = hcat(zeros(T, 9), Matrix{T}(I, 9, 9)) - b = ones(T, 9) - G = -one(T) * I - h = zeros(T, 10) - cones = Cone{T}[Cones.HypoGeoMean{T}(10)] - - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test r.primal_obj ≈ -1 atol=tol rtol=tol - @test r.x ≈ ones(T, 10) atol=tol rtol=tol - @test r.z ≈ vcat(-one(T), fill(inv(T(9)), 9)) atol=tol rtol=tol - @test r.y ≈ fill(inv(T(9)), 9) atol=tol rtol=tol -end - -# TODO test with non-equal alphas -function hypopowermean1(T; options...) - tol = test_tol(T) - c = T[-1, 0, 0] - A = T[0 0 1; 0 1 0] - b = T[0.5, 1] - G = -one(T) * I - h = zeros(T, 3) - - for use_dual in (false, true) - cones = Cone{T}[Cones.HypoPowerMean{T}(ones(T, 2) / 2, use_dual = use_dual)] - - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test r.primal_obj ≈ (use_dual ? 0 : -inv(sqrt(T(2)))) atol=tol rtol=tol - @test r.x[2:3] ≈ [1, 0.5] atol=tol rtol=tol - end -end - -function hypopowermean2(T; options...) - tol = test_tol(T) - l = 4 - c = vcat(zero(T), ones(T, l)) - A = T[one(T) zeros(T, 1, l)] - G = SparseMatrixCSC(-one(T) * I, l + 1, l + 1) - h = zeros(T, l + 1) - - for use_dual in (false, true) - b = use_dual ? [-one(T)] : [one(T)] - cones = Cone{T}[Cones.HypoPowerMean{T}(fill(inv(T(l)), l), use_dual = use_dual)] - - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test r.primal_obj ≈ (use_dual ? 1 : l) atol=tol rtol=tol - @test r.x[2:end] ≈ (use_dual ? fill(inv(T(l)), l) : ones(l)) atol=tol rtol=tol - end -end - -function hypopowermean3(T; options...) - tol = test_tol(T) - l = 4 - c = ones(T, l) - A = zeros(T, 0, l) - b = zeros(T, 0) - G = [zeros(T, 1, l); Matrix{T}(-I, l, l)] - h = zeros(T, l + 1) - - for use_dual in (false, true) - cones = Cone{T}[Cones.HypoPowerMean{T}(fill(inv(T(l)), l), use_dual = use_dual)] - - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test r.primal_obj ≈ 0 atol=tol rtol=tol - @test norm(r.x) ≈ 0 atol=tol rtol=tol - end -end - -function hypopowermean4(T; options...) - tol = test_tol(T) - c = T[-1, 0, 0, 0] - A = zeros(T, 0, 4) - b = zeros(T, 0) - G = vcat(Matrix{T}(-I, 4, 4), T[0, 1, 1, 1]') - h = T[0, 0, 0, 0, 3] - cones = Cone{T}[Cones.HypoPowerMean{T}(fill(inv(T(3)), 3)), Cones.Nonnegative{T}(1)] - - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test r.primal_obj ≈ -1 atol=tol rtol=tol - @test r.x ≈ [1, 1, 1, 1] atol=tol rtol=tol - @test r.s ≈ [1, 1, 1, 1, 0] atol=tol rtol=tol - @test r.z ≈ vcat(-1, fill(inv(T(3)), 4)) atol=tol rtol=tol -end - -function hypopowermean5(T; options...) - tol = test_tol(T) - c = T[-2, 0] - A = zeros(T, 0, 2) - b = zeros(T, 0) - G = sparse([1, 2, 3], [1, 2, 2], T[-1, -1, 1], 3, 2) - h = T[0, 0, 2] - cones = Cone{T}[Cones.HypoPowerMean{T}([one(T)]), Cones.Nonnegative{T}(1)] - - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test r.primal_obj ≈ -4 atol=tol rtol=tol - @test r.x ≈ [2, 2] atol=tol rtol=tol - @test r.s ≈ [2, 2, 0] atol=tol rtol=tol - @test r.z ≈ [-2, 2, 2] atol=tol rtol=tol -end - -function hypopowermean6(T; options...) - tol = test_tol(T) - c = zeros(T, 10) - c[1] = -1 - A = hcat(zeros(T, 9), Matrix{T}(I, 9, 9)) - b = ones(T, 9) - G = -one(T) * I - h = zeros(T, 10) - cones = Cone{T}[Cones.HypoPowerMean{T}(fill(inv(T(9)), 9))] + G = sparse( + [1, 1, 2, 2, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11], + [1, 7, 2, 3, 4, 2, 3, 5, 4, 7, 6, 5, 6, 7], + T[1, -0.5, 1, 1, 1, -1, -1, -1, -1, -0.5, -1, -1, -1, -1], + 11, 7) + h = zeros(T, 11) + h[2] = 3 + cones = Cone{T}[Cones.Nonnegative{T}(2), Cones.EpiPerSquare{T}(3), Cones.EpiPerSquare{T}(3), Cones.EpiPerSquare{T}(3)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal @test r.primal_obj ≈ -1 atol=tol rtol=tol - @test r.x ≈ ones(T, 10) atol=tol rtol=tol - @test r.z ≈ vcat(-one(T), fill(inv(T(9)), 9)) atol=tol rtol=tol - @test r.y ≈ fill(inv(T(9)), 9) atol=tol rtol=tol -end - -function power1(T; options...) - tol = test_tol(T) - c = T[0, 0, 1] - A = T[1 0 0; 0 1 0] - b = T[0.5, 1] - G = Matrix{T}(-I, 3, 3) - h = zeros(T, 3) - - for use_dual in (false, true) - cones = Cone{T}[Cones.Power{T}(ones(T, 2) / 2, 1, use_dual = use_dual)] - - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test r.primal_obj ≈ (use_dual ? -sqrt(T(2)) : -inv(sqrt(T(2)))) atol=tol rtol=tol - @test r.x[3] ≈ (use_dual ? -sqrt(T(2)) : -inv(sqrt(T(2)))) atol=tol rtol=tol - @test r.x[1:2] ≈ [0.5, 1] atol=tol rtol=tol - end -end - -function power2(T; options...) - tol = test_tol(T) - c = T[0, 0, -1, -1] - A = T[0 1 0 0; 1 0 0 0] - b = T[0.5, 1] - G = -one(T) * I - h = zeros(T, 4) - - for use_dual in (false, true) - cones = Cone{T}[Cones.Power{T}(ones(T, 2) / 2, 2, use_dual = use_dual)] - - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test r.primal_obj ≈ (use_dual ? -T(2) : -1) atol=tol rtol=tol - @test norm(r.x[3:4]) ≈ (use_dual ? sqrt(T(2)) : inv(sqrt(T(2)))) atol=tol rtol=tol - @test r.x[3:4] ≈ (use_dual ? ones(T, 2) : fill(inv(T(2)), 2)) atol=tol rtol=tol - @test r.x[1:2] ≈ [1, 0.5] atol=tol rtol=tol - end -end - -function power3(T; options...) - tol = test_tol(T) - l = 4 - c = vcat(fill(T(10), l), zeros(T, 2)) - A = T[zeros(T, 1, l) one(T) zero(T); zeros(T, 1, l) zero(T) one(T)] - G = SparseMatrixCSC(-T(10) * I, l + 2, l + 2) - h = zeros(T, l + 2) - - for use_dual in (false, true) - b = [one(T), zero(T)] - cones = Cone{T}[Cones.Power{T}(fill(inv(T(l)), l), 2, use_dual = use_dual)] - - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test r.primal_obj ≈ (use_dual ? 10 : 10 * T(l)) atol=tol rtol=tol - @test r.x[1:l] ≈ (use_dual ? fill(inv(T(l)), l) : ones(l)) atol=tol rtol=tol - end -end - -function power4(T; options...) - tol = test_tol(T) - l = 4 - c = ones(T, l) - A = zeros(T, 0, l) - b = zeros(T, 0) - G = [zeros(T, 3, l); Matrix{T}(-I, l, l)] - h = zeros(T, l + 3) - - for use_dual in (false, true) - cones = Cone{T}[Cones.Power{T}(fill(inv(T(l)), l), 3, use_dual = use_dual)] - - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test r.primal_obj ≈ 0 atol=tol rtol=tol - @test norm(r.x) ≈ 0 atol=tol rtol=tol - end + rt2 = sqrt(T(2)) + @test r.x ≈ [1, 1, 1, 1, rt2, rt2, 2] atol=tol rtol=tol + @test r.s ≈ [0, 0, 1, 1, rt2, 1, 1, rt2, rt2, rt2, 2] atol=tol rtol=tol + inv3 = inv(T(3)) + rt2inv3 = inv3 * rt2 + @test r.z ≈ [1, inv3, inv3, inv3, -rt2inv3, inv3, inv3, -rt2inv3, rt2inv3, rt2inv3, -2 * inv3] atol=tol rtol=tol end function epinormspectral1(T; options...) @@ -1454,693 +1205,935 @@ function matrixepipersquare2(T; options...) end end -function matrixepipersquare3(T; options...) - tol = 5 * test_tol(T) - Random.seed!(1) - (Xn, Xm) = (3, 4) - for is_complex in (false, true) - R = (is_complex ? Complex{T} : T) - per_idx = (is_complex ? Xn ^ 2 + 1 : Cones.svec_length(Xn) + 1) - W_dim = (is_complex ? 2 : 1) * Xn * Xm - dim = per_idx + W_dim - c = ones(T, W_dim) - A = zeros(T, 0, W_dim) - b = T[] - G = vcat(zeros(T, per_idx, W_dim), Matrix{T}(-10I, W_dim, W_dim)) - h = zeros(T, dim) - U_half = rand(R, Xn, Xn) - U = Hermitian(U_half * U_half') - @views Cones.smat_to_svec!(h[1:(per_idx - 1)], U.data, sqrt(T(2))) +function matrixepipersquare3(T; options...) + tol = 5 * test_tol(T) + Random.seed!(1) + (Xn, Xm) = (3, 4) + for is_complex in (false, true) + R = (is_complex ? Complex{T} : T) + per_idx = (is_complex ? Xn ^ 2 + 1 : Cones.svec_length(Xn) + 1) + W_dim = (is_complex ? 2 : 1) * Xn * Xm + dim = per_idx + W_dim + c = ones(T, W_dim) + A = zeros(T, 0, W_dim) + b = T[] + G = vcat(zeros(T, per_idx, W_dim), Matrix{T}(-10I, W_dim, W_dim)) + h = zeros(T, dim) + U_half = rand(R, Xn, Xn) + U = Hermitian(U_half * U_half') + @views Cones.smat_to_svec!(h[1:(per_idx - 1)], U.data, sqrt(T(2))) + + for use_dual in (false, true) + cones = Cone{T}[Cones.MatrixEpiPerSquare{T, R}(Xn, Xm, use_dual = use_dual)] + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.Optimal + @test norm(r.x) ≈ 0 atol=2tol rtol=2tol + end + end +end + +function generalizedpower1(T; options...) + tol = test_tol(T) + c = T[0, 0, 1] + A = T[1 0 0; 0 1 0] + b = T[0.5, 1] + G = Matrix{T}(-I, 3, 3) + h = zeros(T, 3) + + for use_dual in (false, true) + cones = Cone{T}[Cones.GeneralizedPower{T}(ones(T, 2) / 2, 1, use_dual = use_dual)] + + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.Optimal + @test r.primal_obj ≈ (use_dual ? -sqrt(T(2)) : -inv(sqrt(T(2)))) atol=tol rtol=tol + @test r.x[3] ≈ (use_dual ? -sqrt(T(2)) : -inv(sqrt(T(2)))) atol=tol rtol=tol + @test r.x[1:2] ≈ [0.5, 1] atol=tol rtol=tol + end +end + +function generalizedpower2(T; options...) + tol = test_tol(T) + c = T[0, 0, -1, -1] + A = T[0 1 0 0; 1 0 0 0] + b = T[0.5, 1] + G = -one(T) * I + h = zeros(T, 4) + + for use_dual in (false, true) + cones = Cone{T}[Cones.GeneralizedPower{T}(ones(T, 2) / 2, 2, use_dual = use_dual)] + + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.Optimal + @test r.primal_obj ≈ (use_dual ? -T(2) : -1) atol=tol rtol=tol + @test norm(r.x[3:4]) ≈ (use_dual ? sqrt(T(2)) : inv(sqrt(T(2)))) atol=tol rtol=tol + @test r.x[3:4] ≈ (use_dual ? ones(T, 2) : fill(inv(T(2)), 2)) atol=tol rtol=tol + @test r.x[1:2] ≈ [1, 0.5] atol=tol rtol=tol + end +end + +function generalizedpower3(T; options...) + tol = test_tol(T) + l = 4 + c = vcat(fill(T(10), l), zeros(T, 2)) + A = T[zeros(T, 1, l) one(T) zero(T); zeros(T, 1, l) zero(T) one(T)] + G = SparseMatrixCSC(-T(10) * I, l + 2, l + 2) + h = zeros(T, l + 2) + + for use_dual in (false, true) + b = [one(T), zero(T)] + cones = Cone{T}[Cones.GeneralizedPower{T}(fill(inv(T(l)), l), 2, use_dual = use_dual)] + + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.Optimal + @test r.primal_obj ≈ (use_dual ? 10 : 10 * T(l)) atol=tol rtol=tol + @test r.x[1:l] ≈ (use_dual ? fill(inv(T(l)), l) : ones(l)) atol=tol rtol=tol + end +end + +function generalizedpower4(T; options...) + tol = test_tol(T) + l = 4 + c = ones(T, l) + A = zeros(T, 0, l) + b = zeros(T, 0) + G = [zeros(T, 3, l); Matrix{T}(-I, l, l)] + h = zeros(T, l + 3) + + for use_dual in (false, true) + cones = Cone{T}[Cones.GeneralizedPower{T}(fill(inv(T(l)), l), 3, use_dual = use_dual)] + + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.Optimal + @test r.primal_obj ≈ 0 atol=tol rtol=tol + @test norm(r.x) ≈ 0 atol=tol rtol=tol + end +end +# TODO test with non-equal alphas +function hypopowermean1(T; options...) + tol = test_tol(T) + c = T[-1, 0, 0] + A = T[0 0 1; 0 1 0] + b = T[0.5, 1] + G = -one(T) * I + h = zeros(T, 3) + + for use_dual in (false, true) + cones = Cone{T}[Cones.HypoPowerMean{T}(ones(T, 2) / 2, use_dual = use_dual)] + + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.Optimal + @test r.primal_obj ≈ (use_dual ? 0 : -inv(sqrt(T(2)))) atol=tol rtol=tol + @test r.x[2:3] ≈ [1, 0.5] atol=tol rtol=tol + end +end + +function hypopowermean2(T; options...) + tol = test_tol(T) + l = 4 + c = vcat(zero(T), ones(T, l)) + A = T[one(T) zeros(T, 1, l)] + G = SparseMatrixCSC(-one(T) * I, l + 1, l + 1) + h = zeros(T, l + 1) + + for use_dual in (false, true) + b = use_dual ? [-one(T)] : [one(T)] + cones = Cone{T}[Cones.HypoPowerMean{T}(fill(inv(T(l)), l), use_dual = use_dual)] + + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.Optimal + @test r.primal_obj ≈ (use_dual ? 1 : l) atol=tol rtol=tol + @test r.x[2:end] ≈ (use_dual ? fill(inv(T(l)), l) : ones(l)) atol=tol rtol=tol + end +end + +function hypopowermean3(T; options...) + tol = test_tol(T) + l = 4 + c = ones(T, l) + A = zeros(T, 0, l) + b = zeros(T, 0) + G = [zeros(T, 1, l); Matrix{T}(-I, l, l)] + h = zeros(T, l + 1) + + for use_dual in (false, true) + cones = Cone{T}[Cones.HypoPowerMean{T}(fill(inv(T(l)), l), use_dual = use_dual)] + + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.Optimal + @test r.primal_obj ≈ 0 atol=tol rtol=tol + @test norm(r.x) ≈ 0 atol=tol rtol=tol + end +end + +function hypopowermean4(T; options...) + tol = test_tol(T) + c = T[-1, 0, 0, 0] + A = zeros(T, 0, 4) + b = zeros(T, 0) + G = vcat(Matrix{T}(-I, 4, 4), T[0, 1, 1, 1]') + h = T[0, 0, 0, 0, 3] + cones = Cone{T}[Cones.HypoPowerMean{T}(fill(inv(T(3)), 3)), Cones.Nonnegative{T}(1)] + + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.Optimal + @test r.primal_obj ≈ -1 atol=tol rtol=tol + @test r.x ≈ [1, 1, 1, 1] atol=tol rtol=tol + @test r.s ≈ [1, 1, 1, 1, 0] atol=tol rtol=tol + @test r.z ≈ vcat(-1, fill(inv(T(3)), 4)) atol=tol rtol=tol +end + +function hypopowermean5(T; options...) + tol = test_tol(T) + c = T[-2, 0] + A = zeros(T, 0, 2) + b = zeros(T, 0) + G = sparse([1, 2, 3], [1, 2, 2], T[-1, -1, 1], 3, 2) + h = T[0, 0, 2] + cones = Cone{T}[Cones.HypoPowerMean{T}([one(T)]), Cones.Nonnegative{T}(1)] + + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.Optimal + @test r.primal_obj ≈ -4 atol=tol rtol=tol + @test r.x ≈ [2, 2] atol=tol rtol=tol + @test r.s ≈ [2, 2, 0] atol=tol rtol=tol + @test r.z ≈ [-2, 2, 2] atol=tol rtol=tol +end + +function hypopowermean6(T; options...) + tol = test_tol(T) + c = zeros(T, 10) + c[1] = -1 + A = hcat(zeros(T, 9), Matrix{T}(I, 9, 9)) + b = ones(T, 9) + G = -one(T) * I + h = zeros(T, 10) + cones = Cone{T}[Cones.HypoPowerMean{T}(fill(inv(T(9)), 9))] - for use_dual in (false, true) - cones = Cone{T}[Cones.MatrixEpiPerSquare{T, R}(Xn, Xm, use_dual = use_dual)] - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test norm(r.x) ≈ 0 atol=2tol rtol=2tol - end - end + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.Optimal + @test r.primal_obj ≈ -1 atol=tol rtol=tol + @test r.x ≈ ones(T, 10) atol=tol rtol=tol + @test r.z ≈ vcat(-one(T), fill(inv(T(9)), 9)) atol=tol rtol=tol + @test r.y ≈ fill(inv(T(9)), 9) atol=tol rtol=tol end -function linmatrixineq1(T; options...) +function hypogeomean1(T; options...) tol = test_tol(T) - Random.seed!(1) - for side in [2, 4], R in [T, Complex{T}] - c = T[1] - A = zeros(T, 0, 1) - b = T[] - G = zeros(T, 2, 1) - G[1, 1] = -1 - h = T[0, 2] - A_1_half = rand(R, side, side) - A_1 = Hermitian(A_1_half * A_1_half' + 2I) - F = eigen(A_1) - val_1 = F.values[end] - vec_1 = F.vectors[:, end] - As = [A_1, Hermitian(-vec_1 * vec_1')] - cones = Cone{T}[Cones.LinMatrixIneq{T}(As)] + c = T[-1, 0, 0] + A = T[0 0 1; 0 1 0] + b = T[0.5, 1] + G = -one(T) * I + h = zeros(T, 3) + + for use_dual in (false, true) + cones = Cone{T}[Cones.HypoGeoMean{T}(3, use_dual = use_dual)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - @test r.primal_obj ≈ 2 / val_1 atol=tol rtol=tol - @test r.s ≈ [2 / val_1, 2] atol=tol rtol=tol + @test r.primal_obj ≈ (use_dual ? 0 : -inv(sqrt(T(2)))) atol=tol rtol=tol + @test r.x[2:3] ≈ [1, 0.5] atol=tol rtol=tol end end -function linmatrixineq2(T; options...) +function hypogeomean2(T; options...) tol = test_tol(T) - Random.seed!(1) - for Rs in [[T, T], [Complex{T}, Complex{T}], [T, Complex{T}, T], [Complex{T}, T, T]] - dim = length(Rs) - c = ones(T, dim - 1) - A = zeros(T, 0, dim - 1) - b = T[] - G = vcat(spzeros(T, 1, dim - 1), sparse(-one(T) * I, dim - 1, dim - 1)) - h = zeros(T, dim) - h[1] = 1 - As = Hermitian[] - for R in Rs - A_half = rand(R, 3, 3) - push!(As, Hermitian(A_half * A_half')) - end - As[1] += I - cones = Cone{T}[Cones.LinMatrixIneq{T}(As)] + l = 4 + c = vcat(zero(T), ones(T, l)) + A = T[one(T) zeros(T, 1, l)] + G = SparseMatrixCSC(-one(T) * I, l + 1, l + 1) + h = zeros(T, l + 1) + + for use_dual in (false, true) + b = use_dual ? [-one(T)] : [one(T)] + cones = Cone{T}[Cones.HypoGeoMean{T}(l + 1, use_dual = use_dual)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - @test r.primal_obj < 0 + @test r.primal_obj ≈ (use_dual ? 1 : l) atol=tol rtol=tol + @test r.x[2:end] ≈ (use_dual ? fill(inv(T(l)), l) : ones(l)) atol=tol rtol=tol end end -function linmatrixineq3(T; options...) - dense1 = [1 0; 0 1] - dense2 = [1 0; 0 -1] - sparse1 = sparse(dense1) - sparse2 = sparse(dense2) - diag1 = Diagonal([1, 1]) - diag2 = Diagonal([1, -1]) - # NOTE not all combinations work due to missing methods in LinearAlgebra - As_list = [ - [dense1, dense2], - # [dense1, sparse2], - [dense1, diag2], - # [sparse1, dense2], - [sparse1, sparse2], - # [sparse1, diag2], - [diag1, dense2], - # [diag1, sparse2], - [diag1, diag2], - [I, dense2], - # [I, sparse2], - [I, diag2], - ] - for As in As_list - tol = test_tol(T) - c = T[1] - A = zeros(T, 0, 1) - b = T[] - G = zeros(T, 2, 1) - G[1, 1] = -1 - h = T[0, -1] - cones = Cone{T}[Cones.LinMatrixIneq{T}(As)] +function hypogeomean3(T; options...) + tol = test_tol(T) + l = 4 + c = ones(T, l) + A = zeros(T, 0, l) + b = zeros(T, 0) + G = [zeros(T, 1, l); Matrix{T}(-I, l, l)] + h = zeros(T, l + 1) + + for use_dual in (false, true) + cones = Cone{T}[Cones.HypoGeoMean{T}(l + 1, use_dual = use_dual)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - @test r.primal_obj ≈ 1 atol=tol rtol=tol - @test r.s ≈ [1, -1] atol=tol rtol=tol + @test r.primal_obj ≈ 0 atol=tol rtol=tol + @test norm(r.x) ≈ 0 atol=tol rtol=tol end end -function possemideftri1(T; options...) +function hypogeomean4(T; options...) tol = test_tol(T) - c = T[0, -1, 0] - A = T[1 0 0; 0 0 1] - b = T[0.5, 1] + c = T[-1, 0, 0, 0] + A = zeros(T, 0, 4) + b = zeros(T, 0) + G = vcat(Matrix{T}(-I, 4, 4), T[0, 1, 1, 1]') + h = T[0, 0, 0, 0, 3] + cones = Cone{T}[Cones.HypoGeoMean{T}(4), Cones.Nonnegative{T}(1)] + + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.Optimal + @test r.primal_obj ≈ -1 atol=tol rtol=tol + @test r.x ≈ [1, 1, 1, 1] atol=tol rtol=tol + @test r.s ≈ [1, 1, 1, 1, 0] atol=tol rtol=tol + @test r.z ≈ vcat(-1, fill(inv(T(3)), 4)) atol=tol rtol=tol +end + +function hypogeomean5(T; options...) + tol = test_tol(T) + c = T[-2, 0] + A = zeros(T, 0, 2) + b = zeros(T, 0) + G = sparse([1, 2, 3], [1, 2, 2], T[-1, -1, 1], 3, 2) + h = T[0, 0, 2] + cones = Cone{T}[Cones.HypoGeoMean{T}(2), Cones.Nonnegative{T}(1)] + + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.Optimal + @test r.primal_obj ≈ -4 atol=tol rtol=tol + @test r.x ≈ [2, 2] atol=tol rtol=tol + @test r.s ≈ [2, 2, 0] atol=tol rtol=tol + @test r.z ≈ [-2, 2, 2] atol=tol rtol=tol +end + +function hypogeomean6(T; options...) + tol = test_tol(T) + c = zeros(T, 10) + c[1] = -1 + A = hcat(zeros(T, 9), Matrix{T}(I, 9, 9)) + b = ones(T, 9) G = -one(T) * I - h = zeros(T, 3) - cones = Cone{T}[Cones.PosSemidefTri{T, T}(3)] + h = zeros(T, 10) + cones = Cone{T}[Cones.HypoGeoMean{T}(10)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - @test r.primal_obj ≈ -one(T) atol=tol rtol=tol - @test r.x[2] ≈ one(T) atol=tol rtol=tol + @test r.primal_obj ≈ -1 atol=tol rtol=tol + @test r.x ≈ ones(T, 10) atol=tol rtol=tol + @test r.z ≈ vcat(-one(T), fill(inv(T(9)), 9)) atol=tol rtol=tol + @test r.y ≈ fill(inv(T(9)), 9) atol=tol rtol=tol end -function possemideftri2(T; options...) +function hyporootdettri1(T; options...) tol = test_tol(T) - c = T[0, -1, 0] - A = T[1 0 1] - b = T[0] - G = Diagonal(-one(T) * I, 3) - h = zeros(T, 3) - cones = Cone{T}[Cones.PosSemidefTri{T, T}(3)] + rt2 = sqrt(T(2)) + Random.seed!(1) + side = 3 + for is_complex in (false, true) + dim = (is_complex ? 1 + side^2 : 1 + Cones.svec_length(side)) + R = (is_complex ? Complex{T} : T) + c = T[-1] + A = zeros(T, 0, 1) + b = T[] + G = Matrix{T}(-I, dim, 1) + mat_half = rand(R, side, side) + mat = mat_half * mat_half' + h = zeros(T, dim) + Cones.smat_to_svec!(view(h, 2:dim), mat, rt2) + cones = Cone{T}[Cones.HypoRootdetTri{T, R}(dim)] + + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.Optimal + @test r.x[1] ≈ -r.primal_obj atol=tol rtol=tol + sol_mat = zeros(R, side, side) + Cones.svec_to_smat!(sol_mat, r.s[2:end], rt2) + @test det(cholesky!(Hermitian(sol_mat, :U))) ^ inv(T(side)) ≈ r.s[1] atol=tol rtol=tol + Cones.svec_to_smat!(sol_mat, r.z[2:end] .* T(side), rt2) + @test det(cholesky!(Hermitian(sol_mat, :U))) ^ inv(T(side)) ≈ -r.z[1] atol=tol rtol=tol + end +end + +function hyporootdettri2(T; options...) + tol = test_tol(T) + rt2 = sqrt(T(2)) + Random.seed!(1) + side = 4 + for is_complex in (false, true) + dim = (is_complex ? 1 + side^2 : 1 + Cones.svec_length(side)) + R = (is_complex ? Complex{T} : T) + c = T[1] + A = zeros(T, 0, 1) + b = T[] + G = Matrix{T}(-I, dim, 1) + mat_half = rand(R, side, side) + mat = mat_half * mat_half' + h = zeros(T, dim) + Cones.smat_to_svec!(view(h, 2:dim), mat, rt2) + cones = Cone{T}[Cones.HypoRootdetTri{T, R}(dim, use_dual = true)] - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test r.primal_obj ≈ 0 atol=tol rtol=tol - @test norm(r.x) ≈ 0 atol=tol rtol=tol + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.Optimal + @test r.x[1] ≈ r.primal_obj atol=tol rtol=tol + sol_mat = zeros(R, side, side) + Cones.svec_to_smat!(sol_mat, r.s[2:end] .* T(side), rt2) + @test det(cholesky!(Hermitian(sol_mat, :U))) ^ inv(T(side)) ≈ -r.s[1] atol=tol rtol=tol + Cones.svec_to_smat!(sol_mat, r.z[2:end], rt2) + @test det(cholesky!(Hermitian(sol_mat, :U))) ^ inv(T(side)) ≈ r.z[1] atol=tol rtol=tol + end end -function possemideftri3(T; options...) - tol = test_tol(T) +function hyporootdettri3(T; options...) + # max u: u <= rootdet(W) where W is not full rank + tol = eps(T) ^ 0.15 rt2 = sqrt(T(2)) Random.seed!(1) - c = T[1] - A = zeros(T, 0, 1) - b = T[] - rand_mat = Hermitian(rand(T, 2, 2), :U) - G = reshape(T[-1, 0, -1], 3, 1) - h = -Cones.smat_to_svec!(zeros(T, 3), rand_mat, rt2) - cones = Cone{T}[Cones.PosSemidefTri{T, T}(3)] + side = 3 + for is_complex in (false, true) + dim = (is_complex ? 1 + side^2 : 1 + Cones.svec_length(side)) + R = (is_complex ? Complex{T} : T) + c = T[-1] + A = zeros(T, 0, 1) + b = T[] + G = SparseMatrixCSC(-one(T) * I, dim, 1) + mat_half = T(0.2) * rand(R, side, side - 1) + mat = mat_half * mat_half' + h = zeros(T, dim) + Cones.smat_to_svec!(view(h, 2:dim), mat, rt2) + cones = Cone{T}[Cones.HypoRootdetTri{T, R}(dim)] - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - eig_max = maximum(eigvals(rand_mat)) - @test r.primal_obj ≈ eig_max atol=tol rtol=tol - @test r.x[1] ≈ eig_max atol=tol rtol=tol + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.Optimal + @test r.primal_obj ≈ 0 atol=tol rtol=tol + @test r.x[1] ≈ zero(T) atol=tol rtol=tol + end end -function possemideftri4(T; options...) +function hyporootdettri4(T; options...) tol = test_tol(T) rt2 = sqrt(T(2)) - Random.seed!(1) - s = 3 - rand_mat = Hermitian(rand(T, s, s), :U) - dim = sum(1:s) - c = -Cones.smat_to_svec!(zeros(T, dim), rand_mat, rt2) - A = reshape(Cones.smat_to_svec!(zeros(T, dim), Matrix{T}(I, s, s), rt2), 1, dim) - b = T[1] - G = Diagonal(-one(T) * I, dim) - h = zeros(T, dim) - cones = Cone{T}[Cones.PosSemidefTri{T, T}(dim)] + c = T[-1, 0, 0, 0] + A = zeros(T, 0, 4) + b = T[] + G = zeros(T, 6, 4) + G[1, 1] = G[2, 2] = G[4, 4] = -1 + G[3, 3] = -rt2 + G[5, 2] = G[6, 4] = 1 + h = T[0, 0, 0, 0, 1, 1] + cones = Cone{T}[Cones.HypoRootdetTri{T, T}(4), Cones.Nonnegative{T}(2)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - eig_max = maximum(eigvals(rand_mat)) - @test r.primal_obj ≈ -eig_max atol=tol rtol=tol + @test r.primal_obj ≈ -1 atol=tol rtol=tol + @test r.x ≈ [1, 1, 0, 1] atol=tol rtol=tol + @test r.z ≈ [-1, 0.5, 0, 0.5, 0.5, 0.5] atol=tol rtol=tol end -function possemideftri5(T; options...) +function hypoperlog1(T; options...) tol = test_tol(T) - Trt2 = sqrt(T(2)) - Trt2i = inv(Trt2) - c = T[1, 0, 0, 1] - A = T[0 0 1 0] - b = T[1] - G = Diagonal(-one(T) * I, 4) - h = zeros(T, 4) - cones = Cone{T}[Cones.PosSemidefTri{T, Complex{T}}(4)] + Texph = exp(T(0.5)) + c = T[1, 1, 1] + A = T[0 1 0; 1 0 0] + b = T[2, 1] + G = -one(T) * I + h = zeros(T, 3) + cones = Cone{T}[Cones.HypoPerLog{T}(3)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - @test r.primal_obj ≈ Trt2 atol=tol rtol=tol - @test r.x ≈ [Trt2i, 0, 1, Trt2i] atol=tol rtol=tol + @test r.primal_obj ≈ 2 * Texph + 3 atol=tol rtol=tol + @test r.x ≈ [1, 2, 2 * Texph] atol=tol rtol=tol + @test r.y ≈ -[1 + Texph / 2, 1 + Texph] atol=tol rtol=tol + @test r.z ≈ c + A' * r.y atol=tol rtol=tol end -function possemideftri6(T; options...) +function hypoperlog2(T; options...) tol = test_tol(T) - rt2 = sqrt(T(2)) - Random.seed!(1) - c = T[1] - A = zeros(T, 0, 1) - b = T[] - rand_mat = Hermitian(rand(Complex{T}, 2, 2), :U) - G = reshape(T[-1, 0, 0, -1], 4, 1) - h = -Cones.smat_to_svec!(zeros(T, 4), rand_mat, rt2) - cones = Cone{T}[Cones.PosSemidefTri{T, Complex{T}}(4)] + c = T[-1, 0, 0] + A = T[0 1 0] + b = T[0] + G = Diagonal(-one(T) * I, 3) + h = zeros(T, 3) + cones = Cone{T}[Cones.HypoPerLog{T}(3)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - eig_max = maximum(eigvals(rand_mat)) - @test r.primal_obj ≈ eig_max atol=tol rtol=tol - @test r.x[1] ≈ eig_max atol=tol rtol=tol + @test r.primal_obj ≈ 0 atol=tol rtol=tol end -function possemideftri7(T; options...) +function hypoperlog3(T; options...) tol = test_tol(T) - rt2 = sqrt(T(2)) - Random.seed!(1) - side = 3 - rand_mat = Hermitian(rand(Complex{T}, side, side), :U) - dim = abs2(side) - c = -Cones.smat_to_svec!(zeros(T, dim), rand_mat, rt2) - A = reshape(Cones.smat_to_svec!(zeros(T, dim), Matrix{Complex{T}}(I, side, side), rt2), 1, dim) - b = T[1] - G = Diagonal(-one(T) * I, dim) - h = zeros(T, dim) - cones = Cone{T}[Cones.PosSemidefTri{T, Complex{T}}(dim)] + c = T[1, 1, 1] + A = zeros(T, 0, 3) + b = zeros(T, 0) + G = sparse([1, 2, 3, 4], [1, 2, 3, 1], -ones(T, 4)) + h = zeros(T, 4) + cones = Cone{T}[Cones.HypoPerLog{T}(3), Cones.Nonnegative{T}(1)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - eig_max = maximum(eigvals(rand_mat)) - @test r.primal_obj ≈ -eig_max atol=tol rtol=tol + @test r.primal_obj ≈ 0 atol=tol rtol=tol + @test norm(r.x) ≈ 0 atol=tol rtol=tol + @test isempty(r.y) end -function possemideftri8(T; options...) +function hypoperlog4(T; options...) tol = test_tol(T) - rt2 = sqrt(T(2)) - c = T[1] - A = zeros(T, 0, 1) - b = T[] - G = zeros(T, 15, 1) - G[[1, 3, 6, 10, 15]] .= -1 - h = zeros(T, 15) - @. h[[7, 8, 9, 11, 12, 13]] = rt2 * [1, 1, 0, 1, -1, 1] - cones = Cone{T}[Cones.PosSemidefTri{T, T}(15)] + Texp2 = exp(T(-2)) + c = T[0, 0, 1] + A = T[0 1 0; 1 0 0] + b = T[1, -1] + G = SparseMatrixCSC(-one(T) * I, 3, 3) + h = zeros(T, 3) + cones = Cone{T}[Cones.HypoPerLog{T}(3, use_dual = true)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - rt3 = sqrt(T(3)) - @test r.primal_obj ≈ rt3 atol=tol rtol=tol - @test r.s ≈ [rt3, 0, rt3, 0, 0, rt3, rt2, rt2, 0, rt3, rt2, -rt2, rt2, 0, rt3] atol=tol rtol=tol - inv6 = inv(T(6)) - rt2inv6 = rt2 / 6 - invrt6 = inv(rt2 * rt3) - @test r.z ≈ [inv6, -rt2inv6, inv6, rt2inv6, -rt2inv6, inv6, 0, 0, 0, 0, -invrt6, invrt6, -invrt6, 0, inv(T(2))] atol=tol rtol=tol + @test r.primal_obj ≈ Texp2 atol=tol rtol=tol + @test r.x ≈ [-1, 1, Texp2] atol=tol rtol=tol end -function possemideftri9(T; options...) +function hypoperlog5(T; options...) tol = test_tol(T) - rt2 = sqrt(T(2)) - inv2 = inv(T(2)) - c = vcat(one(T), zeros(T, 9)) - A = zeros(T, 0, 10) - b = T[] - G = zeros(T, 16, 10) - G[1, 2] = G[1, 4] = G[1, 7] = G[1, 8] = G[1, 10] = inv2 - G[1, 1] = G[2, 2] = G[4, 4] = G[7, 7] = G[11, 8] = G[16, 10] = -1 - G[3, 3] = G[5, 5] = G[6, 6] = G[15, 9] = -rt2 - h = zeros(T, 16) - @. h[[8, 9, 10, 12, 13, 14]] = rt2 * [1, 1, 0, 1, -1, 1] - cones = Cone{T}[Cones.Nonnegative{T}(1), Cones.PosSemidefTri{T, T}(15)] + Tlogq = log(T(0.25)) + c = T[-1, 0, 0] + A = T[0 1 1] + b = T[1] + G = sparse([1, 3, 4], [1, 2, 3], -ones(T, 3)) + h = T[0, 1, 0, 0] + cones = Cone{T}[Cones.HypoPerLog{T}(4)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - rt3 = sqrt(T(3)) - @test r.primal_obj ≈ rt2 + rt3 atol=tol rtol=tol - invrt2 = inv(rt2) - invrt3 = inv(rt3) - @test r.s ≈ [0, invrt2 + invrt3, 1 - rt2 / rt3, invrt2 + invrt3, rt2 * invrt3, -rt2 * invrt3, invrt3, rt2, rt2, 0, rt2, rt2, -rt2, rt2, 0, rt3] atol=tol rtol=tol - invrt6 = invrt2 * invrt3 - @test r.z ≈ [1, inv2, 0, inv2, 0, 0, inv2, -inv2, -inv2, 0, inv2, -invrt6, invrt6, -invrt6, 0, inv2] atol=tol rtol=tol + @test r.primal_obj ≈ -Tlogq atol=tol rtol=tol + @test r.x ≈ [Tlogq, 0.5, 0.5] atol=tol rtol=tol + @test r.y ≈ [2] atol=tol rtol=tol end -function possemideftrisparse1(T; options...) - if !(T <: LinearAlgebra.BlasReal) - return # only works with BLAS real types - end +function hypoperlog6(T; options...) tol = test_tol(T) - c = T[0, -1, 0] - A = T[1 0 0; 0 0 1] - b = T[0.5, 1] - G = -one(T) * I - h = zeros(T, 3) - row_idxs = [1, 2, 2] - col_idxs = [1, 1, 2] - cones = Cone{T}[Cones.PosSemidefTriSparse{T, T}(2, row_idxs, col_idxs)] + c = T[-1, 0, 0] + A = zeros(T, 0, 3) + b = zeros(T, 0) + G = sparse([1, 3, 4], [1, 2, 3], -ones(T, 3)) + h = zeros(T, 4) + cones = Cone{T}[Cones.HypoPerLog{T}(4)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - @test r.primal_obj ≈ -one(T) atol=tol rtol=tol - @test r.x[2] ≈ one(T) atol=tol rtol=tol + @test r.primal_obj ≈ 0 atol=tol rtol=tol + @test r.x[1] ≈ 0 atol=tol rtol=tol + @test isempty(r.y) end -function possemideftrisparse2(T; options...) - if !(T <: LinearAlgebra.BlasReal) - return # only works with BLAS real types - end +function hypoperlog7(T; options...) tol = test_tol(T) - Trt2 = sqrt(T(2)) - Trt2i = inv(Trt2) - c = T[1, 0, 0, 1] - A = T[0 0 1 0] - b = T[1] - G = Diagonal(-one(T) * I, 4) - h = zeros(T, 4) - row_idxs = [1, 2, 2] - col_idxs = [1, 1, 2] - cones = Cone{T}[Cones.PosSemidefTriSparse{T, Complex{T}}(2, row_idxs, col_idxs)] + c = zeros(T, 4) + c[1] = -2 + A = zeros(T, 0, 4) + b = zeros(T, 0) + G = sparse( + [1, 2, 2, 3, 4, 5, 6], + [2, 1, 3, 4, 4, 3, 2], + T[1, 1, -1, -1, -1, -1, -1], + 6, 4) + h = zeros(T, 6) + h[1] = 2 + cones = Cone{T}[Cones.Nonnegative{T}(3), Cones.HypoPerLog{T}(3)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - @test r.primal_obj ≈ Trt2 atol=tol rtol=tol - @test r.x ≈ [Trt2i, 0, 1, Trt2i] atol=tol rtol=tol + @test r.primal_obj ≈ -4 atol=tol rtol=tol + @test r.x ≈ [2, 2, 2, 0] atol=tol rtol=tol + @test r.s ≈ [0, 0, 0, 0, 2, 2] atol=tol rtol=tol + @test r.z ≈ [2, 2, 2, -2, -2, 2] atol=tol rtol=tol end -function possemideftrisparse3(T; options...) - if !(T <: LinearAlgebra.BlasReal) - return # only works with BLAS real types - end +function hypoperlogdettri1(T; options...) tol = test_tol(T) rt2 = sqrt(T(2)) Random.seed!(1) - for is_complex in (false, true), side in [1, 2, 5, 20] + side = 4 + for is_complex in (false, true) + dim = (is_complex ? 2 + side^2 : 2 + Cones.svec_length(side)) R = (is_complex ? Complex{T} : T) - rand_mat_L = tril!(sprand(R, side, side, inv(sqrt(side))) + I) - (row_idxs, col_idxs, vals) = findnz(rand_mat_L) - dim = (is_complex ? side + 2 * (length(row_idxs) - side) : length(row_idxs)) - c = zeros(T, dim) - A = zeros(T, 1, dim) - idx = 1 - for (i, v) in enumerate(vals) # scale - if row_idxs[i] == col_idxs[i] - c[idx] = -real(v) - A[idx] = 1 - else - c[idx] = -rt2 * real(v) - if is_complex - idx += 1 - c[idx] = -rt2 * imag(v) - end - end - idx += 1 - end + c = T[-1, 0] + A = T[0 1] b = T[1] - G = Diagonal(-one(T) * I, dim) + G = Matrix{T}(-I, dim, 2) + mat_half = rand(R, side, side) + mat = mat_half * mat_half' + I h = zeros(T, dim) - cones = Cone{T}[Cones.PosSemidefTriSparse{T, R}(side, row_idxs, col_idxs, use_dual = true)] + Cones.smat_to_svec!(view(h, 3:dim), mat, rt2) + cones = Cone{T}[Cones.HypoPerLogdetTri{T, R}(dim)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - eig_max = maximum(eigvals(Hermitian(Matrix(rand_mat_L), :L))) - @test r.primal_obj ≈ -eig_max atol=tol rtol=tol + @test r.x[1] ≈ -r.primal_obj atol=tol rtol=tol + @test r.x[2] ≈ 1 atol=tol rtol=tol + sol_mat = zeros(R, side, side) + Cones.svec_to_smat!(sol_mat, r.s[3:end] / r.s[2], rt2) + @test r.s[2] * logdet(cholesky!(Hermitian(sol_mat, :U))) ≈ r.s[1] atol=tol rtol=tol + Cones.svec_to_smat!(sol_mat, -r.z[3:end] / r.z[1], rt2) + @test r.z[1] * (logdet(cholesky!(Hermitian(sol_mat, :U))) + T(side)) ≈ r.z[2] atol=tol rtol=tol end end -function possemideftrisparse4(T; options...) - if !(T <: LinearAlgebra.BlasReal) - return # only works with BLAS real types - end +function hypoperlogdettri2(T; options...) tol = test_tol(T) rt2 = sqrt(T(2)) - c = T[1] - A = zeros(T, 0, 1) - b = T[] - G = zeros(T, 10, 1) - G[[1, 2, 3, 6, 10]] .= -1 - h = zeros(T, 10) - @. h[[4, 5, 7, 8, 9]] = rt2 * [1, 1, 1, -1, 1] - row_idxs = [1, 2, 3, 4, 4, 4, 5, 5, 5, 5] - col_idxs = [1, 2, 3, 1, 2, 4, 1, 2, 3, 5] - cones = Cone{T}[Cones.PosSemidefTriSparse{T, T}(5, row_idxs, col_idxs)] + Random.seed!(1) + side = 2 + for is_complex in (false, true) + dim = (is_complex ? 2 + side^2 : 2 + Cones.svec_length(side)) + R = (is_complex ? Complex{T} : T) + c = T[0, 1] + A = T[1 0] + b = T[-1] + G = Matrix{T}(-I, dim, 2) + mat_half = rand(R, side, side) + mat = mat_half * mat_half' + h = zeros(T, dim) + Cones.smat_to_svec!(view(h, 3:dim), mat, rt2) + cones = Cone{T}[Cones.HypoPerLogdetTri{T, R}(dim, use_dual = true)] - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - rt3 = sqrt(T(3)) - @test r.primal_obj ≈ rt3 atol=tol rtol=tol - @test r.s ≈ [rt3, rt3, rt3, rt2, rt2, rt3, rt2, -rt2, rt2, rt3] atol=tol rtol=tol - inv6 = inv(T(6)) - rt2inv6 = rt2 / 6 - invrt6 = inv(rt2 * rt3) - @test r.z ≈ [inv6, inv6, inv6, 0, 0, 0, -invrt6, invrt6, -invrt6, inv(T(2))] atol=tol rtol=tol + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.Optimal + @test r.x[2] ≈ r.primal_obj atol=tol rtol=tol + @test r.x[1] ≈ -1 atol=tol rtol=tol + sol_mat = zeros(R, side, side) + Cones.svec_to_smat!(sol_mat, -r.s[3:end] / r.s[1], rt2) + @test r.s[1] * (logdet(cholesky!(Hermitian(sol_mat, :U))) + T(side)) ≈ r.s[2] atol=tol rtol=tol + Cones.svec_to_smat!(sol_mat, r.z[3:end] / r.z[2], rt2) + @test r.z[2] * logdet(cholesky!(Hermitian(sol_mat, :U))) ≈ r.z[1] atol=tol rtol=tol + end end -function possemideftrisparse5(T; options...) - if !(T <: LinearAlgebra.BlasReal) - return # only works with BLAS real types +function hypoperlogdettri3(T; options...) + tol = test_tol(T) + rt2 = sqrt(T(2)) + Random.seed!(1) + side = 3 + for is_complex in (false, true) + dim = (is_complex ? 2 + side^2 : 2 + Cones.svec_length(side)) + R = (is_complex ? Complex{T} : T) + c = T[-1, 0] + A = T[0 1] + b = T[0] + G = SparseMatrixCSC(-one(T) * I, dim, 2) + mat_half = rand(R, side, side) + mat = mat_half * mat_half' + h = zeros(T, dim) + Cones.smat_to_svec!(view(h, 3:dim), mat, rt2) + cones = Cone{T}[Cones.HypoPerLogdetTri{T, R}(dim)] + + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.Optimal + @test r.x[1] ≈ -r.primal_obj atol=tol rtol=tol + @test norm(r.x) ≈ 0 atol=tol rtol=tol end +end + +function hypoperlogdettri4(T; options...) tol = test_tol(T) rt2 = sqrt(T(2)) - inv2 = inv(T(2)) - c = vcat(one(T), zeros(T, 5)) - A = zeros(T, 0, 6) - b = T[] - G = vcat(Matrix(-one(T) * I, 6, 6), zeros(T, 6, 6)) - G[1, 2:end] .= inv2 - h = vcat(zeros(T, 6), rt2 * [1, 1, 0, 1, -1, 1]) - row_idxs = [1, 2, 3, 4, 5, 4, 4, 4, 5, 5, 5] - col_idxs = [1, 2, 3, 4, 5, 1, 2, 3, 1, 2, 3] - cones = Cone{T}[Cones.Nonnegative{T}(1), Cones.PosSemidefTriSparse{T, T}(5, row_idxs, col_idxs, use_dual = true)] + c = T[-1, 0, 0, 0, 0] + A = zeros(T, 1, 5) + A[1, 2] = 1 + b = T[1] + G = zeros(T, 7, 5) + G[1, 1] = G[2, 2] = G[3, 3] = G[5, 5] = -1 + G[4, 4] = -rt2 + G[6, 3] = G[7, 5] = 1 + h = T[0, 0, 0, 0, 0, 1, 1] + cones = Cone{T}[Cones.HypoPerLogdetTri{T, T}(5), Cones.Nonnegative{T}(2)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - rt3 = sqrt(T(3)) - @test r.primal_obj ≈ rt2 + rt3 atol=tol rtol=tol - invrt2 = inv(rt2) - invrt3 = inv(rt3) - @test r.s ≈ [0, invrt2 + invrt3, invrt2 + invrt3, invrt3, rt2, rt3, rt2, rt2, 0, rt2, -rt2, rt2] atol=tol rtol=tol - invrt6 = invrt2 * invrt3 - @test r.z ≈ [1, inv2, inv2, inv2, inv2, inv2, -inv2, -inv2, 0, -invrt6, invrt6, -invrt6] atol=tol rtol=tol + @test r.primal_obj ≈ 0 atol=tol rtol=tol + @test r.x ≈ [0, 1, 1, 0, 1] atol=tol rtol=tol + @test r.y ≈ [-2] atol=tol rtol=tol + @test r.z ≈ [-1, -2, 1, 0, 1, 1, 1] atol=tol rtol=tol end -function doublynonnegativetri1(T; options...) +# TODO add use_dual = true tests +function epiperentropy1(T; options...) tol = test_tol(T) - n = q = 3 - c = T[0, 1, 0] - A = T[1 0 0; 0 0 1] - b = ones(T, 2) - G = SparseMatrixCSC(-one(T) * I, q, n) - for use_dual in [false, true] - use_dual = false - h = (use_dual ? T[1, 0, 1] : zeros(T, q)) - cones = Cone{T}[Cones.DoublyNonnegativeTri{T}(q, use_dual = use_dual)] + Random.seed!(1) + for w_dim in [1, 2, 3] + dim = 2 + w_dim + c = T[1] + A = zeros(T, 0, 1) + b = zeros(T, 0) + G = zeros(T, dim, 1) + G[1, 1] = -1 + h = zeros(T, dim) + h[2] = 1 + w = rand(T, w_dim) .+ 1 + h[3:end] .= w + cones = Cone{T}[Cones.EpiPerEntropy{T}(dim)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - @test r.primal_obj ≈ 0 atol=tol rtol=tol - @test r.x ≈ T[1, 0, 1] atol=tol rtol=tol - @test r.s ≈ T[1, 0, 1] atol=tol rtol=tol + @test r.primal_obj ≈ sum(wi * log(wi) for wi in w) atol=tol rtol=tol end end -function doublynonnegativetri2(T; options...) +function epiperentropy2(T; options...) tol = test_tol(T) - c = T[0, -1, 0] - A = T[1 0 0; 0 0 1] - b = T[1.0, 1.5] - G = -one(T) * I - h = [-inv(T(2)), zero(T), -inv(T(2))] - cones = Cone{T}[Cones.DoublyNonnegativeTri{T}(3)] + for w_dim in [1, 2, 4] + dim = 2 + w_dim + c = fill(-one(T), w_dim) + A = zeros(T, 0, w_dim) + b = zeros(T, 0) + G = vcat(zeros(T, 2, w_dim), Matrix{T}(-I, w_dim, w_dim)) + h = zeros(T, dim) + h[2] = 1 + cones = Cone{T}[Cones.EpiPerEntropy{T}(dim)] - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test r.primal_obj ≈ -one(T) atol=tol rtol=tol - @test r.x[2] ≈ one(T) atol=tol rtol=tol + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.Optimal + @test r.primal_obj ≈ -w_dim atol=tol rtol=tol + @test r.x ≈ fill(1, w_dim) atol=tol rtol=tol + end end -function doublynonnegativetri3(T; options...) +function epiperentropy3(T; options...) tol = test_tol(T) - c = ones(T, 3) - A = T[1 0 1] - b = T[0] - G = Diagonal(-one(T) * I, 3) - h = zeros(T, 3) - cones = Cone{T}[Cones.PosSemidefTri{T, T}(3)] + for w_dim in [2, 4] + dim = 2 + w_dim + c = T[-1] + A = ones(T, 1, 1) + b = T[dim] + G = zeros(T, dim, 1) + G[2, 1] = -1 + h = zeros(T, dim) + h[3:end] .= 1 + cones = Cone{T}[Cones.EpiPerEntropy{T}(dim)] + + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.Optimal + @test r.primal_obj ≈ -dim atol=tol rtol=tol + @test r.x ≈ [dim] atol=tol rtol=tol + end +end + +function epiperentropy4(T; options...) + tol = test_tol(T) + c = T[1] + A = zeros(T, 0, 1) + b = zeros(T, 0) + G = Matrix{T}(-I, 4, 1) + h = T[0, 5, 2, 3] + cones = Cone{T}[Cones.EpiPerEntropy{T}(4)] + + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.Optimal + entr = 2 * log(2 / T(5)) + 3 * log(3 / T(5)) + @test r.primal_obj ≈ entr atol=tol rtol=tol + @test r.s ≈ [entr, 5, 2, 3] atol=tol rtol=tol + @test r.z ≈ [1, 1, log(5 / T(2)) - 1, log(5 / T(3)) - 1] atol=tol rtol=tol +end + +function epiperentropy5(T; options...) + tol = test_tol(T) + c = T[0, -1] + A = zeros(T, 0, 2) + b = zeros(T, 0) + G = vcat(zeros(T, 2, 2), fill(-1, 3, 2), [-1, 0]') + h = T[0, 1, 0, 0, 0, 0] + cones = Cone{T}[Cones.EpiPerEntropy{T}(5), Cones.Nonnegative{T}(1)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - @test r.primal_obj ≈ 0 atol=tol rtol=tol - @test norm(r.x) ≈ 0 atol=tol rtol=tol + @test r.primal_obj ≈ -1 atol=tol rtol=tol + @test r.s ≈ [0, 1, 1, 1, 1, 0] atol=tol rtol=tol + @test r.z ≈ inv(T(3)) * [1, 3, -1, -1, -1, 3] atol=tol rtol=tol end -function hypoperlogdettri1(T; options...) +function epipertraceentropytri1(T; options...) tol = test_tol(T) - rt2 = sqrt(T(2)) Random.seed!(1) + rt2 = sqrt(T(2)) side = 4 - for is_complex in (false, true) - dim = (is_complex ? 2 + side^2 : 2 + Cones.svec_length(side)) - R = (is_complex ? Complex{T} : T) - c = T[-1, 0] - A = T[0 1] - b = T[1] - G = Matrix{T}(-I, dim, 2) - mat_half = rand(R, side, side) - mat = mat_half * mat_half' + I - h = zeros(T, dim) - Cones.smat_to_svec!(view(h, 3:dim), mat, rt2) - cones = Cone{T}[Cones.HypoPerLogdetTri{T, R}(dim)] + svec_dim = Cones.svec_length(side) + cone_dim = svec_dim + 2 + c = T[1] + A = zeros(T, 0, 1) + b = T[] + W = rand(T, side, side) + W = W * W' + v = rand(T) + h = vcat(zero(T), v, Cones.smat_to_svec!(zeros(T, svec_dim), W, rt2)) + G = zeros(T, cone_dim, 1) + G[1, 1] = -1 + cones = Cone{T}[Cones.EpiPerTraceEntropyTri{T}(cone_dim)] - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test r.x[1] ≈ -r.primal_obj atol=tol rtol=tol - @test r.x[2] ≈ 1 atol=tol rtol=tol - sol_mat = zeros(R, side, side) - Cones.svec_to_smat!(sol_mat, r.s[3:end] / r.s[2], rt2) - @test r.s[2] * logdet(cholesky!(Hermitian(sol_mat, :U))) ≈ r.s[1] atol=tol rtol=tol - Cones.svec_to_smat!(sol_mat, -r.z[3:end] / r.z[1], rt2) - @test r.z[1] * (logdet(cholesky!(Hermitian(sol_mat, :U))) + T(side)) ≈ r.z[2] atol=tol rtol=tol - end + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.Optimal + # TODO need https://github.com/JuliaLinearAlgebra/GenericLinearAlgebra.jl/issues/51 to use log with BF + (vals_W, vecs_W) = eigen(Hermitian(W, :U)) + log_W = vecs_W * Diagonal(log.(vals_W)) * vecs_W' + @test r.primal_obj ≈ tr(W * log_W - W * log(v)) atol=tol rtol=tol end -function hypoperlogdettri2(T; options...) +function epipertraceentropytri2(T; options...) tol = test_tol(T) rt2 = sqrt(T(2)) - Random.seed!(1) - side = 2 - for is_complex in (false, true) - dim = (is_complex ? 2 + side^2 : 2 + Cones.svec_length(side)) - R = (is_complex ? Complex{T} : T) - c = T[0, 1] - A = T[1 0] - b = T[-1] - G = Matrix{T}(-I, dim, 2) - mat_half = rand(R, side, side) - mat = mat_half * mat_half' - h = zeros(T, dim) - Cones.smat_to_svec!(view(h, 3:dim), mat, rt2) - cones = Cone{T}[Cones.HypoPerLogdetTri{T, R}(dim, use_dual = true)] + side = 3 + svec_dim = Cones.svec_length(side) + cone_dim = svec_dim + 2 + c = vcat(zero(T), -ones(T, svec_dim)) + A = hcat(one(T), zeros(T, 1, svec_dim)) + b = T[1] + h = vcat(T(5), zeros(T, svec_dim + 1)) + G = vcat(zeros(T, 1, svec_dim + 1), ModelUtilities.vec_to_svec!(Diagonal(-one(T) * I, svec_dim + 1))) + cones = Cone{T}[Cones.EpiPerTraceEntropyTri{T}(cone_dim)] - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test r.x[2] ≈ r.primal_obj atol=tol rtol=tol - @test r.x[1] ≈ -1 atol=tol rtol=tol - sol_mat = zeros(R, side, side) - Cones.svec_to_smat!(sol_mat, -r.s[3:end] / r.s[1], rt2) - @test r.s[1] * (logdet(cholesky!(Hermitian(sol_mat, :U))) + T(side)) ≈ r.s[2] atol=tol rtol=tol - Cones.svec_to_smat!(sol_mat, r.z[3:end] / r.z[2], rt2) - @test r.z[2] * logdet(cholesky!(Hermitian(sol_mat, :U))) ≈ r.z[1] atol=tol rtol=tol - end + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.Optimal + W = Hermitian(Cones.svec_to_smat!(zeros(T, side, side), r.s[3:end], rt2), :U) + @test r.s[2] ≈ 1 atol=tol rtol=tol + @test tr(W * log(W)) ≈ T(5) atol=tol rtol=tol end -function hypoperlogdettri3(T; options...) +function epipertraceentropytri3(T; options...) tol = test_tol(T) rt2 = sqrt(T(2)) - Random.seed!(1) side = 3 - for is_complex in (false, true) - dim = (is_complex ? 2 + side^2 : 2 + Cones.svec_length(side)) - R = (is_complex ? Complex{T} : T) - c = T[-1, 0] - A = T[0 1] - b = T[0] - G = SparseMatrixCSC(-one(T) * I, dim, 2) - mat_half = rand(R, side, side) - mat = mat_half * mat_half' - h = zeros(T, dim) - Cones.smat_to_svec!(view(h, 3:dim), mat, rt2) - cones = Cone{T}[Cones.HypoPerLogdetTri{T, R}(dim)] + svec_dim = Cones.svec_length(side) + cone_dim = svec_dim + 2 + c = vcat(zeros(T, 2), ones(T, svec_dim)) + A = hcat(one(T), zeros(T, 1, svec_dim + 1)) + b = [zero(T)] + h = zeros(T, cone_dim) + G = Diagonal(-one(T) * I, cone_dim) + cones = Cone{T}[Cones.EpiPerTraceEntropyTri{T}(cone_dim)] - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test r.x[1] ≈ -r.primal_obj atol=tol rtol=tol - @test norm(r.x) ≈ 0 atol=tol rtol=tol - end + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.Optimal + @test r.primal_obj ≈ zero(T) atol=tol rtol=tol + @test r.s[1] ≈ zero(T) atol=tol rtol=tol + @test r.s[3:end] ≈ zeros(T, svec_dim) atol=tol rtol=tol end -function hypoperlogdettri4(T; options...) +function epipertraceentropytri4(T; options...) tol = test_tol(T) rt2 = sqrt(T(2)) - c = T[-1, 0, 0, 0, 0] - A = zeros(T, 1, 5) - A[1, 2] = 1 - b = T[1] - G = zeros(T, 7, 5) - G[1, 1] = G[2, 2] = G[3, 3] = G[5, 5] = -1 - G[4, 4] = -rt2 - G[6, 3] = G[7, 5] = 1 - h = T[0, 0, 0, 0, 0, 1, 1] - cones = Cone{T}[Cones.HypoPerLogdetTri{T, T}(5), Cones.Nonnegative{T}(2)] + side = 3 + svec_dim = Cones.svec_length(side) + cone_dim = svec_dim + 2 + c = vcat(zero(T), one(T), zeros(T, svec_dim)) + A = hcat(one(T), zeros(T, 1, svec_dim + 1)) + b = [zero(T)] + h = zeros(T, cone_dim) + G = Diagonal(-one(T) * I, cone_dim) + cones = Cone{T}[Cones.EpiPerTraceEntropyTri{T}(cone_dim)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - @test r.primal_obj ≈ 0 atol=tol rtol=tol - @test r.x ≈ [0, 1, 1, 0, 1] atol=tol rtol=tol - @test r.y ≈ [-2] atol=tol rtol=tol - @test r.z ≈ [-1, -2, 1, 0, 1, 1, 1] atol=tol rtol=tol + @test r.primal_obj ≈ zero(T) atol=tol rtol=tol + @test r.s ≈ zeros(T, cone_dim) atol=tol rtol=tol end -function hyporootdettri1(T; options...) +# TODO add use_dual = true tests +function epirelentropy1(T; options...) tol = test_tol(T) - rt2 = sqrt(T(2)) Random.seed!(1) - side = 3 - for is_complex in (false, true) - dim = (is_complex ? 1 + side^2 : 1 + Cones.svec_length(side)) - R = (is_complex ? Complex{T} : T) - c = T[-1] + for w_dim in [1, 2, 3] + dim = 1 + 2 * w_dim + c = T[1] A = zeros(T, 0, 1) - b = T[] - G = Matrix{T}(-I, dim, 1) - mat_half = rand(R, side, side) - mat = mat_half * mat_half' + b = zeros(T, 0) + G = zeros(T, dim, 1) + G[1, 1] = -1 h = zeros(T, dim) - Cones.smat_to_svec!(view(h, 2:dim), mat, rt2) - cones = Cone{T}[Cones.HypoRootdetTri{T, R}(dim)] + h[2:2:(end - 1)] .= 1 + w = rand(T, w_dim) .+ 1 + h[3:2:end] .= w + cones = Cone{T}[Cones.EpiRelEntropy{T}(dim)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - @test r.x[1] ≈ -r.primal_obj atol=tol rtol=tol - sol_mat = zeros(R, side, side) - Cones.svec_to_smat!(sol_mat, r.s[2:end], rt2) - @test det(cholesky!(Hermitian(sol_mat, :U))) ^ inv(T(side)) ≈ r.s[1] atol=tol rtol=tol - Cones.svec_to_smat!(sol_mat, r.z[2:end] .* T(side), rt2) - @test det(cholesky!(Hermitian(sol_mat, :U))) ^ inv(T(side)) ≈ -r.z[1] atol=tol rtol=tol + @test r.primal_obj ≈ sum(wi * log(wi) for wi in w) atol=tol rtol=tol end end -function hyporootdettri2(T; options...) +function epirelentropy2(T; options...) tol = test_tol(T) - rt2 = sqrt(T(2)) - Random.seed!(1) - side = 4 - for is_complex in (false, true) - dim = (is_complex ? 1 + side^2 : 1 + Cones.svec_length(side)) - R = (is_complex ? Complex{T} : T) - c = T[1] - A = zeros(T, 0, 1) - b = T[] - G = Matrix{T}(-I, dim, 1) - mat_half = rand(R, side, side) - mat = mat_half * mat_half' + for w_dim in [1, 2, 4] + dim = 1 + 2 * w_dim + c = fill(-one(T), w_dim) + A = zeros(T, 0, w_dim) + b = zeros(T, 0) + G = zeros(T, dim, w_dim) + for i in 1:w_dim + G[2i + 1, i] = -1 + end h = zeros(T, dim) - Cones.smat_to_svec!(view(h, 2:dim), mat, rt2) - cones = Cone{T}[Cones.HypoRootdetTri{T, R}(dim, use_dual = true)] + h[2:2:(dim - 1)] .= 1 + cones = Cone{T}[Cones.EpiRelEntropy{T}(dim)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - @test r.x[1] ≈ r.primal_obj atol=tol rtol=tol - sol_mat = zeros(R, side, side) - Cones.svec_to_smat!(sol_mat, r.s[2:end] .* T(side), rt2) - @test det(cholesky!(Hermitian(sol_mat, :U))) ^ inv(T(side)) ≈ -r.s[1] atol=tol rtol=tol - Cones.svec_to_smat!(sol_mat, r.z[2:end], rt2) - @test det(cholesky!(Hermitian(sol_mat, :U))) ^ inv(T(side)) ≈ r.z[1] atol=tol rtol=tol + @test r.primal_obj ≈ -w_dim atol=tol rtol=tol + @test r.x ≈ fill(1, w_dim) atol=tol rtol=tol end end -function hyporootdettri3(T; options...) - # max u: u <= rootdet(W) where W is not full rank - tol = eps(T) ^ 0.15 - rt2 = sqrt(T(2)) - Random.seed!(1) - side = 3 - for is_complex in (false, true) - dim = (is_complex ? 1 + side^2 : 1 + Cones.svec_length(side)) - R = (is_complex ? Complex{T} : T) - c = T[-1] - A = zeros(T, 0, 1) - b = T[] - G = SparseMatrixCSC(-one(T) * I, dim, 1) - mat_half = T(0.2) * rand(R, side, side - 1) - mat = mat_half * mat_half' +function epirelentropy3(T; options...) + tol = test_tol(T) + for w_dim in [2, 4] + dim = 1 + 2 * w_dim + c = fill(-one(T), w_dim) + A = ones(T, 1, w_dim) + b = T[dim] + G = zeros(T, dim, w_dim) + for i in 1:w_dim + G[2i, i] = -1 + end h = zeros(T, dim) - Cones.smat_to_svec!(view(h, 2:dim), mat, rt2) - cones = Cone{T}[Cones.HypoRootdetTri{T, R}(dim)] + h[3:2:end] .= 1 + cones = Cone{T}[Cones.EpiRelEntropy{T}(dim)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal - @test r.primal_obj ≈ 0 atol=tol rtol=tol - @test r.x[1] ≈ zero(T) atol=tol rtol=tol + @test r.primal_obj ≈ -dim atol=tol rtol=tol + @test r.x ≈ fill(dim / w_dim, w_dim) atol=tol rtol=tol end end -function hyporootdettri4(T; options...) +function epirelentropy4(T; options...) tol = test_tol(T) - rt2 = sqrt(T(2)) - c = T[-1, 0, 0, 0] - A = zeros(T, 0, 4) - b = T[] - G = zeros(T, 6, 4) - G[1, 1] = G[2, 2] = G[4, 4] = -1 - G[3, 3] = -rt2 - G[5, 2] = G[6, 4] = 1 - h = T[0, 0, 0, 0, 1, 1] - cones = Cone{T}[Cones.HypoRootdetTri{T, T}(4), Cones.Nonnegative{T}(2)] + c = T[1] + A = zeros(T, 0, 1) + b = zeros(T, 0) + G = Matrix{T}(-I, 5, 1) + h = T[0, 1, 2, 5, 3] + cones = Cone{T}[Cones.EpiRelEntropy{T}(5)] + + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.Optimal + entr = 2 * log(T(2)) + 3 * log(3 / T(5)) + @test r.primal_obj ≈ entr atol=tol rtol=tol + @test r.s ≈ [entr, 1, 2, 5, 3] atol=tol rtol=tol + @test r.z ≈ [1, 2, log(inv(T(2))) - 1, 3 / T(5), log(5 / T(3)) - 1] atol=tol rtol=tol +end + +function epirelentropy5(T; options...) + tol = test_tol(T) + c = T[0, -1] + A = zeros(T, 0, 2) + b = zeros(T, 0) + G = vcat(zeros(T, 1, 2), repeat(T[0 0; -1 -1], 3), [-1, 0]') + h = T[0, 1, 0, 1, 0, 1, 0, 0] + cones = Cone{T}[Cones.EpiRelEntropy{T}(7), Cones.Nonnegative{T}(1)] r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.Optimal @test r.primal_obj ≈ -1 atol=tol rtol=tol - @test r.x ≈ [1, 1, 0, 1] atol=tol rtol=tol - @test r.z ≈ [-1, 0.5, 0, 0.5, 0.5, 0.5] atol=tol rtol=tol + @test r.s ≈ [0, 1, 1, 1, 1, 1, 1, 0] atol=tol rtol=tol + @test r.z ≈ inv(T(3)) * [1, 1, -1, 1, -1, 1, -1, 3] atol=tol rtol=tol end function epitracerelentropytri1(T; options...) @@ -2169,7 +2162,6 @@ end function epitracerelentropytri2(T; options...) tol = test_tol(T) - Random.seed!(1) rt2 = sqrt(T(2)) side = 3 svec_dim = Cones.svec_length(side) @@ -2190,7 +2182,6 @@ end function epitracerelentropytri3(T; options...) tol = test_tol(T) - Random.seed!(1) rt2 = sqrt(T(2)) side = 3 svec_dim = Cones.svec_length(side) @@ -2211,7 +2202,6 @@ end function epitracerelentropytri4(T; options...) tol = test_tol(T) - Random.seed!(1) rt2 = sqrt(T(2)) side = 3 svec_dim = Cones.svec_length(side) @@ -2285,6 +2275,72 @@ function wsosinterpnonnegative3(T; options...) @test r.primal_obj ≈ zero(T) atol=tol rtol=tol end +function wsosinterppossemideftri1(T; options...) + # convexity parameter for (x + 1) ^ 2 * (x - 1) ^ 2 + tol = test_tol(T) + (U, pts, Ps) = ModelUtilities.interpolate(ModelUtilities.Box{T}([-one(T)], [one(T)]), 1) + DynamicPolynomials.@polyvar x + fn = (x + 1) ^ 2 * (x - 1) ^ 2 + H = DynamicPolynomials.differentiate(fn, x, 2) + + c = T[-1] + A = zeros(T, 0, 1) + b = T[] + G = ones(T, U, 1) + h = T[H(pts[u, :]...) for u in 1:U] + cones = Cone{T}[Cones.WSOSInterpPosSemidefTri{T}(1, U, Ps)] + + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.Optimal + @test r.primal_obj ≈ T(4) atol=tol rtol=tol + @test r.x[1] ≈ -T(4) atol=tol rtol=tol +end + +function wsosinterppossemideftri2(T; options...) + # convexity parameter for x[1] ^ 4 - 3 * x[2] ^ 2 + tol = test_tol(T) + (U, pts, Ps) = ModelUtilities.interpolate(ModelUtilities.FreeDomain{T}(2), 1) + DynamicPolynomials.@polyvar x[1:2] + fn = x[1] ^ 4 - 3 * x[2] ^ 2 + H = DynamicPolynomials.differentiate(fn, x, 2) + + c = T[-1] + A = zeros(T, 0, 1) + b = T[] + G = vcat(ones(T, U, 1), zeros(T, U, 1), ones(T, U, 1)) + h = T[H[i, j](pts[u, :]...) for i in 1:2 for j in 1:i for u in 1:U] + ModelUtilities.vec_to_svec!(h, incr = U) + cones = Cone{T}[Cones.WSOSInterpPosSemidefTri{T}(2, U, Ps)] + + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.Optimal + @test r.primal_obj ≈ T(6) atol=tol rtol=tol + @test r.x[1] ≈ -T(6) atol=tol rtol=tol +end + +function wsosinterppossemideftri3(T; options...) + tol = test_tol(T) + (U, pts, Ps) = ModelUtilities.interpolate(ModelUtilities.FreeDomain{T}(1), 3) + DynamicPolynomials.@polyvar x + M = [ + (x + 2x^3) 1; + (-x^2 + 2) (3x^2 - x + 1); + ] + M = M' * M / 10 + + c = T[] + A = zeros(T, 0, 0) + b = T[] + h = T[M[i, j](pts[u, :]...) for i in 1:2 for j in 1:i for u in 1:U] + G = zeros(T, length(h), 0) + ModelUtilities.vec_to_svec!(h, incr = U) + cones = Cone{T}[Cones.WSOSInterpPosSemidefTri{T}(2, U, Ps)] + + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.Optimal + @test r.primal_obj ≈ 0 atol=tol rtol=tol +end + function wsosinterpepinormone1(T; options...) # min t(x) : t(x) >= abs(x ^ 2) on [-1, 1] where t(x) a constant tol = test_tol(T) @@ -2419,150 +2475,6 @@ function wsosinterpepinormeucl3(T; options...) @test abs2.(r.x) ≈ abs2.([fn_sol(pts[u, :]...) for u in 1:U]) atol=tol rtol=tol end -function wsosinterppossemideftri1(T; options...) - # convexity parameter for (x + 1) ^ 2 * (x - 1) ^ 2 - tol = test_tol(T) - (U, pts, Ps) = ModelUtilities.interpolate(ModelUtilities.Box{T}([-one(T)], [one(T)]), 1) - DynamicPolynomials.@polyvar x - fn = (x + 1) ^ 2 * (x - 1) ^ 2 - H = DynamicPolynomials.differentiate(fn, x, 2) - - c = T[-1] - A = zeros(T, 0, 1) - b = T[] - G = ones(T, U, 1) - h = T[H(pts[u, :]...) for u in 1:U] - cones = Cone{T}[Cones.WSOSInterpPosSemidefTri{T}(1, U, Ps)] - - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test r.primal_obj ≈ T(4) atol=tol rtol=tol - @test r.x[1] ≈ -T(4) atol=tol rtol=tol -end - -function wsosinterppossemideftri2(T; options...) - # convexity parameter for x[1] ^ 4 - 3 * x[2] ^ 2 - tol = test_tol(T) - (U, pts, Ps) = ModelUtilities.interpolate(ModelUtilities.FreeDomain{T}(2), 1) - DynamicPolynomials.@polyvar x[1:2] - fn = x[1] ^ 4 - 3 * x[2] ^ 2 - H = DynamicPolynomials.differentiate(fn, x, 2) - - c = T[-1] - A = zeros(T, 0, 1) - b = T[] - G = vcat(ones(T, U, 1), zeros(T, U, 1), ones(T, U, 1)) - h = T[H[i, j](pts[u, :]...) for i in 1:2 for j in 1:i for u in 1:U] - ModelUtilities.vec_to_svec!(h, incr = U) - cones = Cone{T}[Cones.WSOSInterpPosSemidefTri{T}(2, U, Ps)] - - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test r.primal_obj ≈ T(6) atol=tol rtol=tol - @test r.x[1] ≈ -T(6) atol=tol rtol=tol -end - -function wsosinterppossemideftri3(T; options...) - tol = test_tol(T) - (U, pts, Ps) = ModelUtilities.interpolate(ModelUtilities.FreeDomain{T}(1), 3) - DynamicPolynomials.@polyvar x - M = [ - (x + 2x^3) 1; - (-x^2 + 2) (3x^2 - x + 1); - ] - M = M' * M / 10 - - c = T[] - A = zeros(T, 0, 0) - b = T[] - h = T[M[i, j](pts[u, :]...) for i in 1:2 for j in 1:i for u in 1:U] - G = zeros(T, length(h), 0) - ModelUtilities.vec_to_svec!(h, incr = U) - cones = Cone{T}[Cones.WSOSInterpPosSemidefTri{T}(2, U, Ps)] - - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.Optimal - @test r.primal_obj ≈ 0 atol=tol rtol=tol -end - -function primalinfeas1(T; options...) - tol = test_tol(T) - c = T[1, 0] - A = T[1 1] - b = [-T(2)] - G = SparseMatrixCSC(-one(T) * I, 2, 2) - h = zeros(T, 2) - cones = Cone{T}[Cones.Nonnegative{T}(2)] - - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.PrimalInfeasible -end - -function primalinfeas2(T; options...) - tol = test_tol(T) - c = T[1, 1, 1] - A = zeros(T, 0, 3) - b = T[] - G = vcat(SparseMatrixCSC(-one(T) * I, 3, 3), Diagonal([one(T), one(T), -one(T)])) - h = vcat(zeros(T, 3), one(T), one(T), -T(2)) - cones = Cone{T}[Cones.EpiNormEucl{T}(3), Cones.Nonnegative{T}(3)] - - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.PrimalInfeasible -end - -function primalinfeas3(T; options...) - tol = test_tol(T) - c = zeros(T, 3) - A = SparseMatrixCSC(-one(T) * I, 3, 3) - b = [one(T), one(T), T(3)] - G = SparseMatrixCSC(-one(T) * I, 3, 3) - h = zeros(T, 3) - cones = Cone{T}[Cones.HypoPerLog{T}(3)] - - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.PrimalInfeasible -end - -function dualinfeas1(T; options...) - tol = test_tol(T) - c = T[-1, -1, 0] - A = zeros(T, 0, 3) - b = T[] - G = repeat(SparseMatrixCSC(-one(T) * I, 3, 3), 2, 1) - h = zeros(T, 6) - cones = Cone{T}[Cones.EpiNormInf{T, T}(3), Cones.EpiNormInf{T, T}(3, use_dual = true)] - - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.DualInfeasible -end - -function dualinfeas2(T; options...) - tol = test_tol(T) - c = T[-1, 0] - A = zeros(T, 0, 2) - b = T[] - G = T[-1 0; 0 0; 0 -1] - h = T[0, 1, 0] - cones = Cone{T}[Cones.EpiPerSquare{T}(3)] - - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.DualInfeasible -end - -function dualinfeas3(T; options...) - tol = test_tol(T) - c = T[0, 1, 1, 0] - A = zeros(T, 0, 4) - b = T[] - G = SparseMatrixCSC(-one(T) * I, 4, 4) - h = zeros(T, 4) - cones = Cone{T}[Cones.EpiPerSquare{T}(4)] - - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.DualInfeasible -end - function indirect1(T; options...) tol = 1e-3 Random.seed!(1) diff --git a/test/nativesets.jl b/test/nativesets.jl index 675ac8059..e586c5fe9 100644 --- a/test/nativesets.jl +++ b/test/nativesets.jl @@ -22,32 +22,32 @@ inst_infeas = [ inst_cones_few = [ "nonnegative1", + "possemideftri1", + "possemideftri5", + "doublynonnegativetri1", + "possemideftrisparse2", + "possemideftrisparse5", + "linmatrixineq1", "epinorminf4", "epinorminf6", "epinormeucl1", "epipersquare1", - "epiperentropy4", - "epirelentropy4", - "hypoperlog1", - "power1", - "hypogeomean1", - "hypopowermean1", "epinormspectral1", "matrixepipersquare2", - "linmatrixineq1", - "possemideftri1", - "possemideftri5", - "possemideftrisparse2", - "possemideftrisparse5", - "doublynonnegativetri1", - "hypoperlogdettri1", + "generalizedpower1", + "hypopowermean1", + "hypogeomean1", "hyporootdettri1", + "hypoperlog1", + "hypoperlogdettri1", + "epiperentropy4", "epipertraceentropytri1", + "epirelentropy4", "epitracerelentropytri1", "wsosinterpnonnegative1", "wsosinterppossemideftri1", - "wsosinterpepinormeucl2", "wsosinterpepinormone1", + "wsosinterpepinormeucl2", ] # superset of inst_cones_few @@ -56,6 +56,26 @@ inst_cones_many = [ "nonnegative2", "nonnegative3", "nonnegative4", + "possemideftri1", + "possemideftri2", + "possemideftri3", + "possemideftri4", + "possemideftri5", + "possemideftri6", + "possemideftri7", + "possemideftri8", + "possemideftri9", + "doublynonnegativetri1", + "doublynonnegativetri2", + "doublynonnegativetri3", + "possemideftrisparse1", + "possemideftrisparse2", + "possemideftrisparse3", + "possemideftrisparse4", + "possemideftrisparse5", + "linmatrixineq1", + "linmatrixineq2", + "linmatrixineq3", "epinorminf1", "epinorminf2", "epinorminf3", @@ -71,27 +91,17 @@ inst_cones_many = [ "epipersquare2", "epipersquare3", "epipersquare4", - "epiperentropy1", - "epiperentropy2", - "epiperentropy3", - "epiperentropy4", - "epiperentropy5", - "epirelentropy1", - "epirelentropy2", - "epirelentropy3", - "epirelentropy4", - "epirelentropy5", - "hypoperlog1", - "hypoperlog2", - "hypoperlog3", - "hypoperlog4", - "hypoperlog5", - "hypoperlog6", - "hypoperlog7", - "power1", - "power2", - "power3", - "power4", + "epinormspectral1", + "epinormspectral2", + "epinormspectral3", + "epinormspectral4", + "matrixepipersquare1", + "matrixepipersquare2", + "matrixepipersquare3", + "generalizedpower1", + "generalizedpower2", + "generalizedpower3", + "generalizedpower4", "hypopowermean1", "hypopowermean2", "hypopowermean3", @@ -104,61 +114,51 @@ inst_cones_many = [ "hypogeomean4", "hypogeomean5", "hypogeomean6", - "epinormspectral1", - "epinormspectral2", - "epinormspectral3", - "epinormspectral4", - "linmatrixineq1", - "linmatrixineq2", - "linmatrixineq3", - "possemideftri1", - "possemideftri2", - "possemideftri3", - "possemideftri4", - "possemideftri5", - "possemideftri6", - "possemideftri7", - "possemideftri8", - "possemideftri9", - "possemideftrisparse1", - "possemideftrisparse2", - "possemideftrisparse3", - "possemideftrisparse4", - "possemideftrisparse5", - "doublynonnegativetri1", - "doublynonnegativetri2", - "doublynonnegativetri3", - "matrixepipersquare1", - "matrixepipersquare2", - "matrixepipersquare3", - "epitracerelentropytri1", - "epitracerelentropytri2", - "epitracerelentropytri3", - "epitracerelentropytri4", - "epipertraceentropytri1", - "epipertraceentropytri2", - "epipertraceentropytri3", - "epipertraceentropytri4", - "hypoperlogdettri1", - "hypoperlogdettri2", - "hypoperlogdettri3", - "hypoperlogdettri4", "hyporootdettri1", "hyporootdettri2", "hyporootdettri3", "hyporootdettri4", + "hypoperlog1", + "hypoperlog2", + "hypoperlog3", + "hypoperlog4", + "hypoperlog5", + "hypoperlog6", + "hypoperlog7", + "hypoperlogdettri1", + "hypoperlogdettri2", + "hypoperlogdettri3", + "hypoperlogdettri4", + "epiperentropy1", + "epiperentropy2", + "epiperentropy3", + "epiperentropy4", + "epiperentropy5", + "epipertraceentropytri1", + "epipertraceentropytri2", + "epipertraceentropytri3", + "epipertraceentropytri4", + "epirelentropy1", + "epirelentropy2", + "epirelentropy3", + "epirelentropy4", + "epirelentropy5", + "epitracerelentropytri1", + "epitracerelentropytri2", + "epitracerelentropytri3", + "epitracerelentropytri4", "wsosinterpnonnegative1", "wsosinterpnonnegative2", "wsosinterpnonnegative3", + "wsosinterppossemideftri1", + "wsosinterppossemideftri2", + "wsosinterppossemideftri3", "wsosinterpepinormone1", "wsosinterpepinormone2", "wsosinterpepinormone3", "wsosinterpepinormeucl1", "wsosinterpepinormeucl2", "wsosinterpepinormeucl3", - "wsosinterppossemideftri1", - "wsosinterppossemideftri2", - "wsosinterppossemideftri3", ] inst_indirect = [ diff --git a/test/runbarriertests.jl b/test/runbarriertests.jl deleted file mode 100644 index 375d56c9d..000000000 --- a/test/runbarriertests.jl +++ /dev/null @@ -1,52 +0,0 @@ -#= -run barrier tests -=# - -using Test -using Printf -include(joinpath(@__DIR__, "barrier.jl")) - -barrier_test_names = [ - "nonnegative", - "epinorminf", - "epinormeucl", - "epiperentropy", - "epipersquare", - "epipertraceentropytri", - "epitracerelentropytri", - "epirelentropy", - "hypoperlog", - "power", - "hypopowermean", - "hypogeomean", - "epinormspectral", - "linmatrixineq", # NOTE failing on Julia v1.5.1 with ForwardDiff or BigFloat - "possemideftri", - "possemideftrisparse", - "doublynonnegativetri", - "matrixepipersquare", - "hypoperlogdettri", - "hyporootdettri", - "epitracerelentropytri", - "wsosinterpnonnegative", - "wsosinterpepinormone", - "wsosinterpepinormeucl", - "wsosinterppossemideftri", - ] - -real_types = [ - Float64, - Float32, - BigFloat, - ] - -@testset "barrier tests" begin -@testset "$name" for name in barrier_test_names -@testset "$T" for T in real_types - println("$name: $T ...") - test_time = @elapsed eval(Symbol("test_", name, "_barrier"))(T) - @printf("%8.2e seconds\n", test_time) -end -end -end -; diff --git a/test/runconetests.jl b/test/runconetests.jl new file mode 100644 index 000000000..0d30bef86 --- /dev/null +++ b/test/runconetests.jl @@ -0,0 +1,99 @@ +#= +run barrier tests +=# + +using Test +using Printf +import Hypatia.Cones +include(joinpath(@__DIR__, "cone.jl")) + +function cone_types(T::Type{<:Real}) + cones_T = [ + Cones.Nonnegative{T}, + Cones.PosSemidefTri{T, T}, + Cones.PosSemidefTri{T, Complex{T}}, + Cones.DoublyNonnegativeTri{T}, + Cones.LinMatrixIneq{T}, + Cones.EpiNormInf{T, T}, + Cones.EpiNormInf{T, Complex{T}}, + Cones.EpiNormEucl{T}, + Cones.EpiPerSquare{T}, + Cones.EpiNormSpectral{T, T}, + Cones.EpiNormSpectral{T, Complex{T}}, + Cones.MatrixEpiPerSquare{T, T}, + Cones.MatrixEpiPerSquare{T, Complex{T}}, + Cones.GeneralizedPower{T}, + Cones.HypoPowerMean{T}, + Cones.HypoGeoMean{T}, + Cones.HypoRootdetTri{T, T}, + Cones.HypoRootdetTri{T, Complex{T}}, + Cones.HypoPerLog{T}, + Cones.HypoPerLogdetTri{T, T}, + Cones.HypoPerLogdetTri{T, Complex{T}}, + Cones.EpiPerEntropy{T}, + Cones.EpiPerTraceEntropyTri{T}, + Cones.EpiRelEntropy{T}, + Cones.EpiTraceRelEntropyTri{T}, # TODO tighten tol on test_barrier + Cones.WSOSInterpNonnegative{T, T}, + # Cones.WSOSInterpNonnegative{T, Complex{T}}, # TODO + Cones.WSOSInterpPosSemidefTri{T}, + Cones.WSOSInterpEpiNormEucl{T}, + Cones.WSOSInterpEpiNormOne{T}, + ] + if T <: LinearAlgebra.BlasReal + append!(cones_T, [ + Cones.PosSemidefTriSparse{T, T}, + Cones.PosSemidefTriSparse{T, Complex{T}}, + ]) + end + return cones_T +end + +@testset "cone tests" begin + +println("starting oracle tests") +@testset "oracle tests" begin +real_types = [ + Float64, + Float32, + BigFloat, + ] +@testset "$cone" for T in real_types, cone in cone_types(T) + println("$cone") + test_time = @elapsed test_oracles(cone) + @printf("%8.2e seconds\n", test_time) +end +end + +println("\nstarting barrier tests") +@testset "barrier tests" begin +real_types = [ + Float64, + # Float32, + # BigFloat, + ] +@testset "$cone" for T in real_types, cone in cone_types(T) + println("$cone") + test_time = @elapsed test_barrier(cone) + @printf("%8.2e seconds\n", test_time) +end +end + +println("\nstarting allocation tests") +@testset "allocation tests" begin +real_types = [ + Float64, + # Float32, + # BigFloat, + ] +@testset "$cone" for T in real_types, cone in cone_types(T) + println("\n$cone") + test_time = @elapsed allocs = get_allocs(cone) + display(allocs) + @printf("%8.2e seconds\n", test_time) +end +println() +end + +end +; diff --git a/test/runtests.jl b/test/runtests.jl index 6d9367af6..436953d54 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,7 +8,7 @@ using Hypatia test_files = [ "modelutilities", - "barrier", + "cone", "native", "moi", "examples",