Merge pull request #252 from chriscoey/generic_reals

WIP: allow generic real number types (higher precision) where possible

Manifest.toml

@@ -29,9 +29,9 @@ version = "0.5.4"
 
 [[CSV]]
 deps = ["CategoricalArrays", "DataFrames", "Dates", "Mmap", "Parsers", "PooledArrays", "Profile", "Tables", "Unicode", "WeakRefStrings"]
-git-tree-sha1 = "f64241c9688ae3eb003bce26ffd9ed863cfb824c"
+git-tree-sha1 = "239240112fc5f18fb934942a5b6f207a8f9e45d2"
 uuid = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
-version = "0.5.2"
+version = "0.5.5"
 
 [[Calculus]]
 deps = ["Compat"]
@@ -77,9 +77,9 @@ version = "4.0.0"
 
 [[DataFrames]]
 deps = ["CategoricalArrays", "Compat", "IteratorInterfaceExtensions", "Missings", "PooledArrays", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "StatsBase", "TableTraits", "Tables", "Unicode"]
-git-tree-sha1 = "3891a62fd843662af9f78f25bdd415530b9b9c1e"
+git-tree-sha1 = "279baa6358fd5e944deccab88434f69c74cfc722"
 uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
-version = "0.18.2"
+version = "0.18.3"
 
 [[DataStructures]]
 deps = ["InteractiveUtils", "OrderedCollections", "Random", "Serialization", "Test"]
@@ -145,6 +145,12 @@ git-tree-sha1 = "429b70c7a08ec34a20e0864366c20dfb8966e5e5"
 uuid = "92c85e6c-cbff-5e0c-80f7-495c94daaecd"
 version = "0.5.1"
 
+[[GenericLinearAlgebra]]
+deps = ["LinearAlgebra", "Printf", "Random", "Test"]
+git-tree-sha1 = "ca235f9c4652b31525232a36d7832f5ee681d76a"
+uuid = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
+version = "0.1.0"
+
 [[InteractiveUtils]]
 deps = ["Markdown"]
 uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
@@ -230,9 +236,9 @@ version = "0.9.7"
 
 [[Parsers]]
 deps = ["Dates", "Test"]
-git-tree-sha1 = "162855122e7d2b7ffbcdd8d19d6b18472f2117bc"
+git-tree-sha1 = "eaed2db080700f1013f0fc05667ecb2a082265a1"
 uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
-version = "0.3.4"
+version = "0.3.5"
 
 [[Pkg]]
 deps = ["Dates", "LibGit2", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"]
@@ -364,10 +370,10 @@ deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
 uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
 
 [[SumOfSquares]]
-deps = ["JuMP", "LinearAlgebra", "MathOptInterface", "MultivariateMoments", "MultivariatePolynomials", "Pkg", "PolyJuMP", "Reexport", "SemialgebraicSets", "SparseArrays", "Test"]
-git-tree-sha1 = "2fde6e851a9b9b5d1f33b108d8fa71857346051b"
+deps = ["JuMP", "LinearAlgebra", "MathOptInterface", "MultivariateMoments", "MultivariatePolynomials", "Pkg", "PolyJuMP", "Reexport", "SemialgebraicSets", "SparseArrays"]
+git-tree-sha1 = "153b6aa62a5712644ae0786417fc249ea5794580"
 uuid = "4b9e565b-77fc-50a5-a571-1244f986bda1"
-version = "0.3.1"
+version = "0.3.2"
 
 [[TableTraits]]
 deps = ["IteratorInterfaceExtensions"]
@@ -377,9 +383,9 @@ version = "1.0.0"
 
 [[Tables]]
 deps = ["IteratorInterfaceExtensions", "LinearAlgebra", "Requires", "TableTraits", "Test"]
-git-tree-sha1 = "83b4a0261e5d01274f12b35d4c2212386fb15569"
+git-tree-sha1 = "c5d784c61e9d243a5a6a8458d19f535b70bdedeb"
 uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
-version = "0.2.3"
+version = "0.2.4"
 
 [[Test]]
 deps = ["Distributed", "InteractiveUtils", "Logging", "Random"]

Project.toml

@@ -13,6 +13,7 @@ DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 GSL = "92c85e6c-cbff-5e0c-80f7-495c94daaecd"
+GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"

examples/envelope/native.jl

@@ -15,9 +15,9 @@ using Test
 import Hypatia
 const HYP = Hypatia
 const CO = HYP.Cones
-const MU = HYP.ModelUtilities
 const MO = HYP.Models
 const SO = HYP.Solvers
+const MU = HYP.ModelUtilities
 
 function envelope(
     npoly::Int,
@@ -51,7 +51,7 @@ function envelope(
         h = zeros(npoly * U)
     end
 
-    cones = [CO.WSOSPolyInterp(U, [P0, PWts...], !primal_wsos) for k in 1:npoly]
+    cones = [CO.WSOSPolyInterp{Float64, Float64}(U, [P0, PWts...], !primal_wsos) for k in 1:npoly]
     cone_idxs = [(1 + (k - 1) * U):(k * U) for k in 1:npoly]
 
     return (c = c, A = A, b = b, G = G, h = h, cones = cones, cone_idxs = cone_idxs)
@@ -67,8 +67,8 @@ envelope6() = envelope(2, 30, 1, 30, primal_wsos = false)
 function test_envelope(instance::Function; options, rseed::Int = 1)
     Random.seed!(rseed)
     d = instance()
-    model = MO.PreprocessedLinearModel(d.c, d.A, d.b, d.G, d.h, d.cones, d.cone_idxs)
-    solver = SO.HSDSolver(model; options...)
+    model = MO.PreprocessedLinearModel{Float64}(d.c, d.A, d.b, d.G, d.h, d.cones, d.cone_idxs)
+    solver = SO.HSDSolver{Float64}(model; options...)
     SO.solve(solver)
     r = SO.get_certificates(solver, model, test = true, atol = 1e-4, rtol = 1e-4)
     @test r.status == :Optimal

examples/linearopt/native.jl

@@ -15,7 +15,6 @@ const HYP = Hypatia
 const CO = HYP.Cones
 const MO = HYP.Models
 const SO = HYP.Solvers
-const MU = HYP.ModelUtilities
 
 function linearopt(
     m::Int,
@@ -30,7 +29,7 @@ function linearopt(
     c = rand(0.0:9.0, n)
     G = Diagonal(-1.0I, n) # TODO uniformscaling
     h = zeros(n)
-    cones = [CO.Nonnegative(n)]
+    cones = [CO.Nonnegative{Float64}(n)]
     cone_idxs = [1:n]
     return (c = c, A = A, b = b, G = G, h = h, cones = cones, cone_idxs = cone_idxs)
 end
@@ -43,8 +42,8 @@ linearopt4() = linearopt(15, 20, nzfrac = 1 / 4)
 function test_linearopt(instance::Function; options, rseed::Int = 1)
     Random.seed!(rseed)
     d = instance()
-    model = MO.PreprocessedLinearModel(d.c, d.A, d.b, d.G, d.h, d.cones, d.cone_idxs)
-    solver = SO.HSDSolver(model; options...)
+    model = MO.PreprocessedLinearModel{Float64}(d.c, d.A, d.b, d.G, d.h, d.cones, d.cone_idxs)
+    solver = SO.HSDSolver{Float64}(model; options...)
     SO.solve(solver)
     r = SO.get_certificates(solver, model, test = true, atol = 1e-3, rtol = 1e-3)
     @test r.status == :Optimal

examples/polymin/JuMP.jl

@@ -7,8 +7,6 @@ see description in examples/polymin/native.jl
 import Random
 using LinearAlgebra
 using Test
-import MathOptInterface
-const MOI = MathOptInterface
 import JuMP
 import SumOfSquares
 import Hypatia

examples/polymin/native.jl

@@ -20,6 +20,7 @@ const HYP = Hypatia
 const CO = HYP.Cones
 const MO = HYP.Models
 const SO = HYP.Solvers
+const MU = HYP.ModelUtilities
 
 include(joinpath(@__DIR__, "data.jl"))
 
@@ -48,7 +49,7 @@ function polyminreal(
         G = Diagonal(-1.0I, U) # TODO use UniformScaling
         h = zeros(U)
     end
-    cones = [CO.WSOSPolyInterp(U, [P0, PWts...], !primal_wsos)]
+    cones = [CO.WSOSPolyInterp{Float64, Float64}(U, [P0, PWts...], !primal_wsos)]
     cone_idxs = [1:U]
 
     return (c = c, A = A, b = b, G = G, h = h, cones = cones, cone_idxs = cone_idxs, true_obj = true_obj)
@@ -75,7 +76,7 @@ polyminreal17() = polyminreal(:lotkavolterra, 3, primal_wsos = false)
 function polymincomplex(
     polyname::Symbol,
     halfdeg::Int;
-    primal_wsos = true,
+    primal_wsos::Bool = true,
     sample_factor::Int = 100,
     use_QR::Bool = false,
     )
@@ -142,7 +143,7 @@ function polymincomplex(
         G = Diagonal(-1.0I, U)
         h = zeros(U)
     end
-    cones = [CO.WSOSPolyInterp(U, P_data, !primal_wsos)]
+    cones = [CO.WSOSPolyInterp{Float64, ComplexF64}(U, P_data, !primal_wsos)]
     cone_idxs = [1:U]
 
     return (c = c, A = A, b = b, G = G, h = h, cones = cones, cone_idxs = cone_idxs, true_obj = true_obj)
@@ -166,8 +167,8 @@ polymincomplex14() = polymincomplex(:denseunit1d, 2, primal_wsos = false)
 function test_polymin(instance::Function; options, rseed::Int = 1)
     Random.seed!(rseed)
     d = instance()
-    model = MO.PreprocessedLinearModel(d.c, d.A, d.b, d.G, d.h, d.cones, d.cone_idxs)
-    solver = SO.HSDSolver(model; options...)
+    model = MO.PreprocessedLinearModel{Float64}(d.c, d.A, d.b, d.G, d.h, d.cones, d.cone_idxs)
+    solver = SO.HSDSolver{Float64}(model; options...)
     SO.solve(solver)
     r = SO.get_certificates(solver, model, test = true, atol = 1e-4, rtol = 1e-4)
     @test r.status == :Optimal

src/Cones/Cones.jl

@@ -7,58 +7,55 @@ functions and caches for cones
 module Cones
 
 using LinearAlgebra
+import LinearAlgebra.BlasFloat
 using ForwardDiff
 using DiffResults
-# using TimerOutputs
+import Hypatia.HypReal
+import Hypatia.HypRealOrComplex
+import Hypatia.hyp_AtA!
+import Hypatia.hyp_chol!
+import Hypatia.hyp_ldiv_chol_L!
+import Hypatia.rt2
+import Hypatia.rt2i
 
-abstract type Cone end
+abstract type Cone{T <: HypReal} end
 
 include("orthant.jl")
 include("epinorminf.jl")
 include("epinormeucl.jl")
 include("epipersquare.jl")
+include("semidefinite.jl")
 include("hypoperlog.jl")
 include("epiperpower.jl")
-include("epipersumexp.jl")
 include("hypogeomean.jl")
 include("epinormspectral.jl")
-include("semidefinite.jl")
 include("hypoperlogdet.jl")
+include("epipersumexp.jl")
 include("wsospolyinterp.jl")
 include("wsospolyinterpmat.jl")
 include("wsospolyinterpsoc.jl")
 
 use_dual(cone::Cone) = cone.use_dual
-load_point(cone::Cone, point::AbstractVector{Float64}) = (cone.point = point)
+load_point(cone::Cone{T}, point::AbstractVector{T}) where {T <: HypReal} = (cone.point = point)
 dimension(cone::Cone) = cone.dim
 
 function factorize_hess(cone::Cone)
-    @. cone.H2 = cone.H
-
-    # cone.F = bunchkaufman!(Symmetric(cone.H2, :U), true, check = false)
-    # return issuccess(cone.F)
-
-    cone.F = cholesky!(Symmetric(cone.H2, :U), Val(true), check = false)
+    copyto!(cone.H2, cone.H)
+    cone.F = hyp_chol!(Symmetric(cone.H2, :U))
     return isposdef(cone.F)
 end
 
 grad(cone::Cone) = cone.g
 hess(cone::Cone) = Symmetric(cone.H, :U)
-inv_hess(cone::Cone) = inv(cone.F)
+inv_hess(cone::Cone) = Symmetric(inv(cone.F), :U)
 hess_fact(cone::Cone) = cone.F
-# hessL(cone::Cone) = cone.F.L
-# inv_hessL(cone::Cone) = inv(cone.F.L)
-hess_prod!(prod::AbstractArray{Float64}, arr::AbstractArray{Float64}, cone::Cone) = mul!(prod, Symmetric(cone.H, :U), arr)
-inv_hess_prod!(prod::AbstractArray{Float64}, arr::AbstractArray{Float64}, cone::Cone) = ldiv!(prod, cone.F, arr)
-# hessL_prod!(prod::AbstractArray{Float64}, arr::AbstractArray{Float64}, cone::Cone) = mul!(prod, cone.F.L, arr)
-# inv_hessL_prod!(prod::AbstractArray{Float64}, arr::AbstractArray{Float64}, cone::Cone) = ldiv!(prod, cone.F.L, arr)
+hess_prod!(prod::AbstractVecOrMat{T}, arr::AbstractVecOrMat{T}, cone::Cone) where {T <: HypReal} = mul!(prod, Symmetric(cone.H, :U), arr)
+inv_hess_prod!(prod::AbstractVecOrMat{T}, arr::AbstractVecOrMat{T}, cone::Cone) where {T <: HypReal} = ldiv!(prod, cone.F, arr)
 
 # utilities for converting between smat and svec forms (lower triangle) for symmetric matrices
 # TODO only need to do lower triangle if use symmetric matrix types
-const rt2 = sqrt(2)
-const rt2i = inv(rt2)
 
-function smat_to_svec!(vec::AbstractVector{T}, mat::AbstractMatrix{T}) where {T <: Real}
+function smat_to_svec!(vec::AbstractVector{T}, mat::AbstractMatrix{T}) where {T <: HypReal}
     k = 1
     m = size(mat, 1)
     for i in 1:m, j in 1:i
@@ -72,7 +69,7 @@ function smat_to_svec!(vec::AbstractVector{T}, mat::AbstractMatrix{T}) where {T
     return vec
 end
 
-function svec_to_smat!(mat::AbstractMatrix{T}, vec::AbstractVector{T}) where {T <: Real}
+function svec_to_smat!(mat::AbstractMatrix{T}, vec::AbstractVector{T}) where {T <: HypReal}
     k = 1
     m = size(mat, 1)
     for i in 1:m, j in 1:i
@@ -86,7 +83,7 @@ function svec_to_smat!(mat::AbstractMatrix{T}, vec::AbstractVector{T}) where {T
     return mat
 end
 
-function smat_to_svec!(vec::AbstractVector{T}, mat::AbstractMatrix{Complex{T}}) where {T <: Real}
+function smat_to_svec!(vec::AbstractVector{T}, mat::AbstractMatrix{Complex{T}}) where {T <: HypReal}
     k = 1
     m = size(mat, 1)
     for i in 1:m, j in 1:i
@@ -104,7 +101,7 @@ function smat_to_svec!(vec::AbstractVector{T}, mat::AbstractMatrix{Complex{T}})
     return vec
 end
 
-function svec_to_smat!(mat::AbstractMatrix{Complex{T}}, vec::AbstractVector{T}) where {T <: Real}
+function svec_to_smat!(mat::AbstractMatrix{Complex{T}}, vec::AbstractVector{T}) where {T <: HypReal}
     k = 1
     m = size(mat, 1)
     for i in 1:m, j in 1:i

src/Cones/epinormeucl.jl

@@ -8,53 +8,53 @@ barrier from "Self-Scaled Barriers and Interior-Point Methods for Convex Program
 -log(u^2 - norm(w)^2)
 =#
 
-mutable struct EpiNormEucl <: Cone
+mutable struct EpiNormEucl{T <: HypReal} <: Cone{T}
     use_dual::Bool
     dim::Int
 
-    point::AbstractVector{Float64}
-    g::Vector{Float64}
-    H::Matrix{Float64}
-    Hi::Matrix{Float64}
+    point::AbstractVector{T}
+    g::Vector{T}
+    H::Matrix{T}
+    Hi::Matrix{T}
 
-    function EpiNormEucl(dim::Int, is_dual::Bool)
-        cone = new()
+    function EpiNormEucl{T}(dim::Int, is_dual::Bool) where {T <: HypReal}
+        cone = new{T}()
         cone.use_dual = is_dual
         cone.dim = dim
         return cone
     end
 end
 
-EpiNormEucl(dim::Int) = EpiNormEucl(dim, false)
+EpiNormEucl{T}(dim::Int) where {T <: HypReal} = EpiNormEucl{T}(dim, false)
 
-function setup_data(cone::EpiNormEucl)
+function setup_data(cone::EpiNormEucl{T}) where {T <: HypReal}
     dim = cone.dim
-    cone.g = Vector{Float64}(undef, dim)
-    cone.H = Matrix{Float64}(undef, dim, dim)
+    cone.g = Vector{T}(undef, dim)
+    cone.H = Matrix{T}(undef, dim, dim)
     cone.Hi = similar(cone.H)
     return
 end
 
 get_nu(cone::EpiNormEucl) = 1
 
-set_initial_point(arr::AbstractVector{Float64}, cone::EpiNormEucl) = (@. arr = 0.0; arr[1] = 1.0; arr)
+set_initial_point(arr::AbstractVector{T}, cone::EpiNormEucl{T}) where {T <: HypReal} = (@. arr = zero(T); arr[1] = one(T); arr)
 
-function check_in_cone(cone::EpiNormEucl)
+function check_in_cone(cone::EpiNormEucl{T}) where {T <: HypReal}
     u = cone.point[1]
     w = view(cone.point, 2:cone.dim)
-    if u <= 0.0
+    if u <= zero(T)
         return false
     end
     dist = abs2(u) - sum(abs2, w)
-    if dist <= 0.0
+    if dist <= zero(T)
         return false
     end
 
     @. cone.g = cone.point / dist
-    cone.g[1] *= -1.0
+    cone.g[1] = -cone.g[1]
 
     Hi = cone.Hi
-    mul!(Hi, cone.point, cone.point') # TODO syrk
+    mul!(Hi, cone.point, cone.point') # TODO use syrk
     @. Hi += Hi
     Hi[1, 1] -= dist
     for j in 2:cone.dim
@@ -64,7 +64,7 @@ function check_in_cone(cone::EpiNormEucl)
     H = cone.H
     @. H = Hi
     for j in 2:cone.dim
-        H[1, j] = H[j, 1] = -H[j, 1]
+        H[1, j] = H[j, 1] = -H[j, 1] # TODO only need upper tri
     end
     @. H *= abs2(inv(dist))
 
@@ -73,4 +73,4 @@ end
 
 inv_hess(cone::EpiNormEucl) = Symmetric(cone.Hi, :U)
 
-inv_hess_prod!(prod::AbstractArray{Float64}, arr::AbstractArray{Float64}, cone::EpiNormEucl) = mul!(prod, Symmetric(cone.Hi, :U), arr)
+inv_hess_prod!(prod::AbstractVecOrMat{T}, arr::AbstractVecOrMat{T}, cone::EpiNormEucl{T}) where {T <: HypReal} = mul!(prod, Symmetric(cone.Hi, :U), arr)