From b50939b96b9048e67096ac7aa25bb1828528df44 Mon Sep 17 00:00:00 2001
From: Chris Coey <coey.chris@gmail.com>
Date: Thu, 1 Apr 2021 00:37:08 +1000
Subject: [PATCH 1/6] modify matrixcompletion example

---
 examples/matrixcompletion/JuMP.jl   | 17 +++++++----------
 examples/matrixcompletion/native.jl |  8 +++-----
 test/runexamplestests.jl            |  2 +-
 3 files changed, 11 insertions(+), 16 deletions(-)

diff --git a/examples/matrixcompletion/JuMP.jl b/examples/matrixcompletion/JuMP.jl
index 028affd03..899e1f6bc 100644
--- a/examples/matrixcompletion/JuMP.jl
+++ b/examples/matrixcompletion/JuMP.jl
@@ -14,20 +14,17 @@ function build(inst::MatrixCompletionJuMP{T}) where {T <: Float64}
     @assert size_ratio >= 1
     num_cols = size_ratio * num_rows
 
-    (rows, cols, Avals) = findnz(sprand(num_rows, num_cols, 0.1))
-    mat_to_vec_idx(i::Int, j::Int) = (j - 1) * num_rows + i
-    is_unknown = fill(true, num_rows * num_cols)
-    for (i, j) in zip(rows, cols)
-        is_unknown[mat_to_vec_idx(i, j)] = false
-    end
+    (rows, cols, Avals) = findnz(sprandn(num_rows, num_cols, 0.8))
+    is_known = vec(Matrix(sparse(rows, cols, trues(length(Avals)), num_rows, num_cols)))
 
     model = JuMP.Model()
-    JuMP.@variable(model, X[1:num_rows, 1:num_cols])
+    JuMP.@variable(model, X[1:length(is_known)])
     JuMP.@variable(model, t)
     JuMP.@objective(model, Min, t)
-    JuMP.@constraint(model, vcat(t, vec(X)) in MOI.NormSpectralCone(num_rows, num_cols))
-    JuMP.@constraint(model, vcat(1, X[is_unknown]) in MOI.GeometricMeanCone(1 + sum(is_unknown)))
-    JuMP.@constraint(model, X[.!is_unknown] .== Avals)
+    JuMP.@constraint(model, vcat(t, X) in MOI.NormSpectralCone(num_rows, num_cols))
+    X_unknown = X[.!is_known]
+    JuMP.@constraint(model, vcat(1, X_unknown) in MOI.GeometricMeanCone(1 + length(X_unknown)))
+    JuMP.@constraint(model, X[is_known] .== Avals)
 
     return model
 end
diff --git a/examples/matrixcompletion/native.jl b/examples/matrixcompletion/native.jl
index 26c4d423e..6cca1f621 100644
--- a/examples/matrixcompletion/native.jl
+++ b/examples/matrixcompletion/native.jl
@@ -26,18 +26,16 @@ function build(inst::MatrixCompletionNative{T}) where {T <: Real}
     mn = m * n
     rt2 = sqrt(T(2))
 
-    num_known = round(Int, mn * 0.1)
+    num_known = round(Int, mn * 0.8)
     known_rows = rand(1:m, num_known)
     known_cols = rand(1:n, num_known)
-    known_vals = rand(T, num_known) .- T(0.5)
-
-    mat_to_vec_idx(i::Int, j::Int) = (j - 1) * m + i
+    known_vals = 2 * rand(T, num_known) .- 1
 
     is_known = fill(false, mn)
     # h for the rows that X (the matrix and not epigraph variable) participates in
     h_norm_x = zeros(T, m * n)
     for (k, (i, j)) in enumerate(zip(known_rows, known_cols))
-        known_idx = mat_to_vec_idx(i, j)
+        known_idx = (j - 1) * m + i
         # if not using the epinorminf cone, indices relate to X'
         h_norm_x[known_idx] = known_vals[k]
         is_known[known_idx] = true
diff --git a/test/runexamplestests.jl b/test/runexamplestests.jl
index 3022ea108..6652d6dca 100644
--- a/test/runexamplestests.jl
+++ b/test/runexamplestests.jl
@@ -17,7 +17,7 @@ default_options = (
     verbose = false,
     # verbose = true,
     default_tol_relax = 10,
-    stepper = Solvers.CombinedStepper{Float64}(),
+    # stepper = Solvers.CombinedStepper{Float64}(),
     # stepper = Solvers.PredOrCentStepper{Float64}(),
     iter_limit = 250,
     )

From 449f522cd2f8547f57b5684b42c7af6f1c5486c4 Mon Sep 17 00:00:00 2001
From: Chris Coey <coey.chris@gmail.com>
Date: Thu, 1 Apr 2021 15:01:17 +1000
Subject: [PATCH 2/6] update rel entr cone definition and oracles; improve cone
 oracles tests printing

---
 src/Cones/epiperentropy.jl    |   2 +-
 src/Cones/epirelentropy.jl    | 117 ++++++++++++++++------------------
 src/MathOptInterface/cones.jl |  52 +++++++--------
 test/cone.jl                  | 105 +++++++++++++++++-------------
 test/nativeinstances.jl       |  25 ++++----
 test/runconetests.jl          |   3 +-
 6 files changed, 158 insertions(+), 146 deletions(-)

diff --git a/src/Cones/epiperentropy.jl b/src/Cones/epiperentropy.jl
index 97fc24add..61880a9d9 100644
--- a/src/Cones/epiperentropy.jl
+++ b/src/Cones/epiperentropy.jl
@@ -78,7 +78,7 @@ function update_feas(cone::EpiPerEntropy{T}) where T
     v = cone.point[2]
     @views w = cone.point[3:cone.dim]
 
-    if (v > eps(T)) && all(wi -> wi > eps(T), w)
+    if (v > eps(T)) && all(>(eps(T)), w)
         @. cone.lwv = log(w / v)
         cone.z = u - dot(w, cone.lwv)
         cone.is_feas = (cone.z > eps(T))
diff --git a/src/Cones/epirelentropy.jl b/src/Cones/epirelentropy.jl
index fc393501a..2f8d9c987 100644
--- a/src/Cones/epirelentropy.jl
+++ b/src/Cones/epirelentropy.jl
@@ -1,6 +1,6 @@
 #=
 (closure of) epigraph of sum of perspectives of entropies (AKA vector relative entropy cone)
-(u in R, v in R_+^n, w in R_+^n) : u >= sum_i w_i*log(w_i/v_i) TODO update description here for non-contiguous v/w
+(u in R, v in R_+^n, w in R_+^n) : u >= sum_i w_i*log(w_i/v_i)
 
 barrier from "Primal-Dual Interior-Point Methods for Domain-Driven Formulations" by Karimi & Tuncel, 2019
 -log(u - sum_i w_i*log(w_i/v_i)) - sum_i (log(v_i) + log(w_i))
@@ -10,8 +10,8 @@ mutable struct EpiRelEntropy{T <: Real} <: Cone{T}
     use_dual_barrier::Bool
     dim::Int
     w_dim::Int
-    v_idxs
-    w_idxs
+    v_idxs::UnitRange{Int64}
+    w_idxs::UnitRange{Int64}
 
     point::Vector{T}
     dual_point::Vector{T}
@@ -24,8 +24,6 @@ mutable struct EpiRelEntropy{T <: Real} <: Cone{T}
     hess_updated::Bool
     inv_hess_updated::Bool
     inv_hess_aux_updated::Bool
-    inv_hess_sqrt_aux_updated::Bool
-    use_inv_hess_sqrt::Bool
     is_feas::Bool
     hess::Symmetric{T, Matrix{T}}
     inv_hess::Symmetric{T, SparseMatrixCSC{T, Int}}
@@ -40,12 +38,6 @@ mutable struct EpiRelEntropy{T <: Real} <: Cone{T}
     Hivv::Vector{T}
     Hivw::Vector{T}
     Hiww::Vector{T}
-    rtiuu::T
-    rtiuv::Vector{T}
-    rtiuw::Vector{T}
-    rtivv::Vector{T}
-    rtivw::Vector{T}
-    rtiww::Vector{T}
     temp1::Vector{T}
     temp2::Vector{T}
 
@@ -58,31 +50,18 @@ mutable struct EpiRelEntropy{T <: Real} <: Cone{T}
         cone = new{T}()
         cone.use_dual_barrier = use_dual
         cone.dim = dim
-        cone.w_dim = div(dim - 1, 2)
-        cone.v_idxs = 2:2:dim
-        cone.w_idxs = 3:2:dim
+        w_dim = cone.w_dim = div(dim - 1, 2)
+        cone.v_idxs = 1 .+ (1:w_dim)
+        cone.w_idxs = w_dim .+ cone.v_idxs
         return cone
     end
 end
 
 use_heuristic_neighborhood(cone::EpiRelEntropy) = false
 
-reset_data(cone::EpiRelEntropy) = (cone.feas_updated = cone.grad_updated = cone.hess_updated = cone.inv_hess_updated = cone.inv_hess_aux_updated = cone.inv_hess_sqrt_aux_updated = false)
-
-function use_sqrt_oracles(cone::EpiRelEntropy{T}) where T
-    cone.inv_hess_sqrt_aux_updated || update_inv_hess_sqrt_aux(cone)
-    cone.use_inv_hess_sqrt || return false
-    # check numerics
-    # TODO tune
-    tol = sqrt(sqrt(eps(T))) * dimension(cone)
-    nu = get_nu(cone)
-    vec1 = cone.vec1
-    hess_sqrt_prod!(vec1, cone.point, cone)
-    hess_viol = abs(1 - sum(abs2, vec1) / nu)
-    inv_hess_sqrt_prod!(vec1, grad(cone), cone)
-    inv_hess_viol = abs(1 - sum(abs2, vec1) / nu)
-    return (max(hess_viol, inv_hess_viol) < tol)
-end
+reset_data(cone::EpiRelEntropy) = (cone.feas_updated = cone.grad_updated = cone.hess_updated = cone.inv_hess_updated = cone.inv_hess_aux_updated = false)
+
+use_sqrt_oracles(cone::EpiRelEntropy) = false # TODO remove in favor of BHB oracles
 
 # TODO only allocate the fields we use
 function setup_extra_data(cone::EpiRelEntropy{T}) where {T <: Real}
@@ -95,14 +74,8 @@ function setup_extra_data(cone::EpiRelEntropy{T}) where {T <: Real}
     cone.Hivv = zeros(T, w_dim)
     cone.Hivw = zeros(T, w_dim)
     cone.Hiww = zeros(T, w_dim)
-    cone.rtiuv = zeros(T, w_dim)
-    cone.rtiuw = zeros(T, w_dim)
-    cone.rtivv = zeros(T, w_dim)
-    cone.rtivw = zeros(T, w_dim)
-    cone.rtiww = zeros(T, w_dim)
     cone.temp1 = zeros(T, w_dim)
     cone.temp2 = zeros(T, w_dim)
-    cone.use_inv_hess_sqrt = true
     return cone
 end
 
@@ -121,7 +94,7 @@ function update_feas(cone::EpiRelEntropy{T}) where T
     @views v = cone.point[cone.v_idxs]
     @views w = cone.point[cone.w_idxs]
 
-    if all(vi -> vi > eps(T), v) && all(wi -> wi > eps(T), w)
+    if all(>(eps(T)), v) && all(>(eps(T)), w)
         @. cone.lwv = log(w / v)
         cone.z = u - dot(w, cone.lwv)
         cone.is_feas = (cone.z > eps(T))
@@ -138,7 +111,7 @@ function is_dual_feas(cone::EpiRelEntropy{T}) where T
     @views v = cone.dual_point[cone.v_idxs]
     @views w = cone.dual_point[cone.w_idxs]
 
-    if all(vi -> vi > eps(T), v) && u > eps(T)
+    if all(>(eps(T)), v) && u > eps(T)
         return all(u * (1 + log(vi / u)) + wi > eps(T) for (vi, wi) in zip(v, w))
     end
 
@@ -249,19 +222,41 @@ end
 # updates for nonzero values in the inverse Hessian
 function update_inv_hess(cone::EpiRelEntropy{T}) where T
     cone.inv_hess_aux_updated || update_inv_hess_aux(cone)
+    dim = cone.dim
+    w_dim = cone.w_dim
 
     if !isdefined(cone, :inv_hess)
         # initialize sparse idxs for upper triangle of inverse Hessian
-        cone.inv_hess = Symmetric(sparse_arrow_block2(T, cone.w_dim), :U)
+        nnz_tri = 2 * dim - 1 + w_dim
+        I = Vector{Int}(undef, nnz_tri)
+        J = Vector{Int}(undef, nnz_tri)
+        idxs1 = 1:dim
+        @views I[idxs1] .= 1
+        @views J[idxs1] .= idxs1
+        idxs2 = (dim + 1):(2 * dim - 1)
+        @views I[idxs2] .= 2:dim
+        @views J[idxs2] .= 2:dim
+        idxs3 = (2 * dim):nnz_tri
+        @views I[idxs3] .= 1 .+ (1:w_dim)
+        @views J[idxs3] .= (1 + w_dim) .+ (1:w_dim)
+        V = ones(T, nnz_tri)
+        cone.inv_hess = Symmetric(sparse(I, J, V, dim, dim), :U)
     end
 
     # modify nonzeros of upper triangle of inverse Hessian
     nzval = cone.inv_hess.data.nzval
     nzval[1] = cone.Hiuu
-    nz_idx = 1
+    nz_idx = 2
+    @inbounds for i in 1:cone.w_dim
+        nzval[nz_idx] = cone.Hiuv[i]
+        nzval[nz_idx + 1] = cone.Hivv[i]
+        nz_idx += 2
+    end
     @inbounds for i in 1:cone.w_dim
-        @. nzval[nz_idx .+ (1:5)] = (cone.Hiuv[i], cone.Hivv[i], cone.Hiuw[i], cone.Hivw[i], cone.Hiww[i])
-        nz_idx += 5
+        nzval[nz_idx] = cone.Hiuw[i]
+        nzval[nz_idx + 1] = cone.Hivw[i]
+        nzval[nz_idx + 2] = cone.Hiww[i]
+        nz_idx += 3
     end
 
     cone.inv_hess_updated = true
@@ -298,26 +293,24 @@ end
 
 function inv_hess_prod!(prod::AbstractVecOrMat, arr::AbstractVecOrMat, cone::EpiRelEntropy)
     cone.inv_hess_aux_updated || update_inv_hess_aux(cone)
-    return arrow_prod(prod, arr, cone.Hiuu, cone.Hiuv, cone.Hiuw, cone.Hivv, cone.Hivw, cone.Hiww)
-end
-
-function update_inv_hess_sqrt_aux(cone::EpiRelEntropy)
-    cone.inv_hess_aux_updated || update_inv_hess_aux(cone)
-    @assert !cone.inv_hess_sqrt_aux_updated
-    cone.rtiuu = arrow_sqrt(cone.Hiuu, cone.Hiuv, cone.Hiuw, cone.Hivv, cone.Hivw, cone.Hiww, cone.rtiuv, cone.rtiuw, cone.rtivv, cone.rtivw, cone.rtiww)
-    cone.use_inv_hess_sqrt = !iszero(cone.rtiuu)
-    cone.inv_hess_sqrt_aux_updated = true
-    return
-end
 
-function hess_sqrt_prod!(prod::AbstractVecOrMat, arr::AbstractVecOrMat, cone::EpiRelEntropy)
-    @assert cone.inv_hess_sqrt_aux_updated && cone.use_inv_hess_sqrt
-    return inv_arrow_sqrt_prod(prod, arr, cone.rtiuu, cone.rtiuv, cone.rtiuw, cone.rtivv, cone.rtivw, cone.rtiww)
-end
+    @inbounds @views begin
+        u_arr = arr[1, :]
+        u_prod = prod[1, :]
+        v_arr = arr[cone.v_idxs, :]
+        v_prod = prod[cone.v_idxs, :]
+        w_arr = arr[cone.w_idxs, :]
+        w_prod = prod[cone.w_idxs, :]
+        @. u_prod = cone.Hiuu * u_arr
+        mul!(u_prod, v_arr', cone.Hiuv, true, true)
+        mul!(u_prod, w_arr', cone.Hiuw, true, true)
+        mul!(v_prod, cone.Hiuv, u_arr')
+        mul!(w_prod, cone.Hiuw, u_arr')
+        @. v_prod += cone.Hivv * v_arr + cone.Hivw * w_arr
+        @. w_prod += cone.Hivw * v_arr + cone.Hiww * w_arr
+    end
 
-function inv_hess_sqrt_prod!(prod::AbstractVecOrMat, arr::AbstractVecOrMat, cone::EpiRelEntropy)
-    @assert cone.inv_hess_sqrt_aux_updated && cone.use_inv_hess_sqrt
-    return arrow_sqrt_prod(prod, arr, cone.rtiuu, cone.rtiuv, cone.rtiuw, cone.rtivv, cone.rtivw, cone.rtiww)
+    return prod
 end
 
 function correction(cone::EpiRelEntropy, primal_dir::AbstractVector)
@@ -363,8 +356,8 @@ end
 # TODO remove this in favor of new hess_nz_count etc functions that directly use uu, uw, ww etc
 inv_hess_nz_count(cone::EpiRelEntropy) = 3 * cone.dim - 2 + 2 * cone.w_dim
 inv_hess_nz_count_tril(cone::EpiRelEntropy) = 2 * cone.dim - 1 + cone.w_dim
-inv_hess_nz_idxs_col(cone::EpiRelEntropy, j::Int) = (j == 1 ? (1:cone.dim) : (iseven(j) ? [1, j, j + 1] : [1, j - 1, j]))
-inv_hess_nz_idxs_col_tril(cone::EpiRelEntropy, j::Int) = (j == 1 ? (1:cone.dim) : (iseven(j) ? [j, j + 1] : [j]))
+inv_hess_nz_idxs_col(cone::EpiRelEntropy, j::Int) = (j == 1 ? (1:cone.dim) : (j <= (1 + cone.w_dim) ? [1, j, j + cone.w_dim] : [1, j - cone.w_dim, j]))
+inv_hess_nz_idxs_col_tril(cone::EpiRelEntropy, j::Int) = (j == 1 ? (1:cone.dim) : (j <= (1 + cone.w_dim) ? [j, j + cone.w_dim] : [j]))
 
 # see analysis in https://github.com/lkapelevich/HypatiaSupplements.jl/tree/master/centralpoints
 function get_central_ray_epirelentropy(w_dim::Int)
diff --git a/src/MathOptInterface/cones.jl b/src/MathOptInterface/cones.jl
index 0fe886407..56fb0cb9c 100644
--- a/src/MathOptInterface/cones.jl
+++ b/src/MathOptInterface/cones.jl
@@ -30,32 +30,32 @@ permute_affine(::MOI.AbstractVectorSet, idxs::AbstractVector) = idxs
 rescale_affine(::MOI.AbstractVectorSet, vals::AbstractVector) = vals
 rescale_affine(::MOI.AbstractVectorSet, vals::AbstractVector, ::AbstractVector) = vals
 
-# transformation for MOI relative entropy cone
-needs_untransform(::MOI.RelativeEntropyCone) = true
-
-function untransform_affine(::MOI.RelativeEntropyCone, vals::AbstractVector)
-    w_dim = div(length(vals) - 1, 2)
-    v_vals = vals[2:2:(end - 1)]
-    w_vals = vals[3:2:end]
-    vals[2:(w_dim + 1)] = v_vals
-    vals[(w_dim + 2):end] = w_vals
-    return vals
-end
-
-function permute_affine(cone::MOI.RelativeEntropyCone, idxs::AbstractVector)
-    dim = MOI.dimension(cone)
-    w_dim = div(dim - 1, 2)
-    new_idxs = collect(idxs)
-    for (i, idx) in enumerate(idxs)
-        idx <= 1 && continue
-        if idx <= 1 + w_dim
-            new_idxs[i] = 2 * (idx - 1)
-        else
-            new_idxs[i] = 2 * (idx - w_dim) - 1
-        end
-    end
-    return new_idxs
-end
+# # transformation for MOI relative entropy cone
+# needs_untransform(::MOI.RelativeEntropyCone) = true
+#
+# function untransform_affine(::MOI.RelativeEntropyCone, vals::AbstractVector)
+#     w_dim = div(length(vals) - 1, 2)
+#     v_vals = vals[2:2:(end - 1)]
+#     w_vals = vals[3:2:end]
+#     vals[2:(w_dim + 1)] = v_vals
+#     vals[(w_dim + 2):end] = w_vals
+#     return vals
+# end
+#
+# function permute_affine(cone::MOI.RelativeEntropyCone, idxs::AbstractVector)
+#     dim = MOI.dimension(cone)
+#     w_dim = div(dim - 1, 2)
+#     new_idxs = collect(idxs)
+#     for (i, idx) in enumerate(idxs)
+#         idx <= 1 && continue
+#         if idx <= 1 + w_dim
+#             new_idxs[i] = 2 * (idx - 1)
+#         else
+#             new_idxs[i] = 2 * (idx - w_dim) - 1
+#         end
+#     end
+#     return new_idxs
+# end
 
 # transformations (transposition of matrix) for MOI rectangular matrix cones with matrix of more rows than columns
 
diff --git a/test/cone.jl b/test/cone.jl
index 58556552d..14c9f809f 100644
--- a/test/cone.jl
+++ b/test/cone.jl
@@ -5,7 +5,6 @@ tests for primitive cone barrier oracles
 using Test
 import Random
 import Random.randn
-import DataStructures
 import GenericLinearAlgebra.eigen
 using LinearAlgebra
 using SparseArrays
@@ -134,50 +133,69 @@ function test_barrier(
     return
 end
 
-# get memory allocations for oracles
-function get_allocs(
+# show time and memory allocation for oracles
+function show_time_alloc(
     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
+    println("dimension: ", dim)
 
-    allocs[:setup_data] = @allocated Cones.setup_data(cone)
+    println("setup_data")
+    @time 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)
+    @assert 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)
+    @assert 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)
+    Cones.reset_data(cone)
+
+    Cones.load_point(cone, point)
+    println("is_feas")
+    @time Cones.is_feas(cone)
+
+    Cones.load_dual_point(cone, dual_point)
+    println("is_dual_feas")
+    @time Cones.is_dual_feas(cone)
+
+    println("grad")
+    @time Cones.grad(cone)
+    println("hess")
+    @time Cones.hess(cone)
+    println("inv_hess")
+    @time 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)
+    println("hess_prod")
+    @time Cones.hess_prod!(point2, point1, cone)
+    println("inv_hess_prod")
+    @time 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)
+        println("hess_sqrt_prod")
+        @time Cones.hess_sqrt_prod!(point2, point1, cone)
+        println("inv_hess_sqrt_prod")
+        @time Cones.inv_hess_sqrt_prod!(point2, point1, cone)
     end
 
     if Cones.use_correction(cone)
-        allocs[:correction] = @allocated Cones.correction(cone, point1)
+        println("correction")
+        @time Cones.correction(cone, point1)
     end
 
-    allocs[:in_neighborhood] = @allocated Cones.in_neighborhood(cone, one(T), one(T))
+    println("in_neighborhood")
+    @time Cones.in_neighborhood(cone, one(T), one(T))
 
-    return allocs
+    return
 end
 
 function perturb_scale(
@@ -276,7 +294,7 @@ function test_barrier(C::Type{<:Cones.Nonnegative})
     test_barrier(C(3), barrier)
 end
 
-get_allocs(C::Type{<:Cones.Nonnegative}) = get_allocs(C(9))
+show_time_alloc(C::Type{<:Cones.Nonnegative}) = show_time_alloc(C(9))
 
 
 # PosSemidefTri
@@ -292,7 +310,7 @@ function test_barrier(C::Type{Cones.PosSemidefTri{T, R}}) where {T, R}
     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)))
+show_time_alloc(C::Type{Cones.PosSemidefTri{T, R}}) where {T, R} = show_time_alloc(C(dim_herm(4, R)))
 
 
 # DoublyNonnegativeTri
@@ -315,7 +333,7 @@ function test_barrier(C::Type{Cones.DoublyNonnegativeTri{T}}) where T
     test_barrier(C(Cones.svec_length(dW)), barrier)
 end
 
-get_allocs(C::Type{<:Cones.DoublyNonnegativeTri}) = get_allocs(C(10))
+show_time_alloc(C::Type{<:Cones.DoublyNonnegativeTri}) = show_time_alloc(C(10))
 
 
 # PosSemidefTriSparse
@@ -360,7 +378,7 @@ function test_barrier(C::Type{Cones.PosSemidefTriSparse{T, Complex{T}}}) where T
     test_barrier(C(dW, row_idxs, col_idxs), barrier)
 end
 
-get_allocs(C::Type{<:Cones.PosSemidefTriSparse}) = get_allocs(C(15, rand_sppsd_pattern(15)...))
+show_time_alloc(C::Type{<:Cones.PosSemidefTriSparse}) = show_time_alloc(C(15, rand_sppsd_pattern(15)...))
 
 
 # LinMatrixIneq
@@ -379,7 +397,7 @@ function test_barrier(C::Type{Cones.LinMatrixIneq{T}}) where T
     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)))
+show_time_alloc(C::Type{Cones.LinMatrixIneq{T}}) where T = show_time_alloc(C(rand_herms(3, [T, Complex{T}, T, Complex{T}], T)))
 
 
 # EpiNormInf
@@ -399,7 +417,7 @@ function test_barrier(C::Type{Cones.EpiNormInf{T, R}}) where {T, R}
     test_barrier(C(1 + dim_vec(dw, R)), barrier)
 end
 
-get_allocs(C::Type{<:Cones.EpiNormInf}) = get_allocs(C(9))
+show_time_alloc(C::Type{<:Cones.EpiNormInf}) = show_time_alloc(C(9))
 
 
 # EpiNormEucl
@@ -417,7 +435,7 @@ function test_barrier(C::Type{<:Cones.EpiNormEucl})
     test_barrier(C(3), barrier)
 end
 
-get_allocs(C::Type{<:Cones.EpiNormEucl}) = get_allocs(C(9))
+show_time_alloc(C::Type{<:Cones.EpiNormEucl}) = show_time_alloc(C(9))
 
 
 # EpiPerSquare
@@ -435,7 +453,7 @@ function test_barrier(C::Type{<:Cones.EpiPerSquare})
     test_barrier(C(4), barrier)
 end
 
-get_allocs(C::Type{<:Cones.EpiPerSquare}) = get_allocs(C(9))
+show_time_alloc(C::Type{<:Cones.EpiPerSquare}) = show_time_alloc(C(9))
 
 
 # EpiNormSpectral
@@ -455,7 +473,7 @@ function test_barrier(C::Type{Cones.EpiNormSpectral{T, R}}) where {T, R}
     test_barrier(C(dr, ds), barrier)
 end
 
-get_allocs(C::Type{<:Cones.EpiNormSpectral}) = get_allocs(C(2, 2))
+show_time_alloc(C::Type{<:Cones.EpiNormSpectral}) = show_time_alloc(C(2, 2))
 
 
 # MatrixEpiPerSquare
@@ -477,7 +495,7 @@ function test_barrier(C::Type{Cones.MatrixEpiPerSquare{T, R}}) where {T, R}
     test_barrier(C(dr, ds), barrier)
 end
 
-get_allocs(C::Type{<:Cones.MatrixEpiPerSquare}) = get_allocs(C(2, 2))
+show_time_alloc(C::Type{<:Cones.MatrixEpiPerSquare}) = show_time_alloc(C(2, 2))
 
 
 # GeneralizedPower
@@ -497,7 +515,7 @@ function test_barrier(C::Type{Cones.GeneralizedPower{T}}) where T
     test_barrier(C(alpha, dw), barrier)
 end
 
-get_allocs(C::Type{Cones.GeneralizedPower{T}}) where T = get_allocs(C(rand_powers(T, 4), 5))
+show_time_alloc(C::Type{Cones.GeneralizedPower{T}}) where T = show_time_alloc(C(rand_powers(T, 4), 5))
 
 
 # HypoPowerMean
@@ -519,7 +537,7 @@ function test_barrier(C::Type{Cones.HypoPowerMean{T}}) where T
     test_barrier(C(alpha), barrier)
 end
 
-get_allocs(C::Type{Cones.HypoPowerMean{T}}) where T = get_allocs(C(rand_powers(T, 8)))
+show_time_alloc(C::Type{Cones.HypoPowerMean{T}}) where T = show_time_alloc(C(rand_powers(T, 8)))
 
 
 # HypoGeoMean
@@ -538,7 +556,7 @@ function test_barrier(C::Type{<:Cones.HypoGeoMean})
     test_barrier(C(3), barrier)
 end
 
-get_allocs(C::Type{<:Cones.HypoGeoMean}) = get_allocs(C(9))
+show_time_alloc(C::Type{<:Cones.HypoGeoMean}) = show_time_alloc(C(9))
 
 
 # HypoRootdetTri
@@ -558,7 +576,7 @@ function test_barrier(C::Type{Cones.HypoRootdetTri{T, R}}) where {T, R}
     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)))
+show_time_alloc(C::Type{Cones.HypoRootdetTri{T, R}}) where {T, R} = show_time_alloc(C(1 + dim_herm(3, R)))
 
 
 # HypoPerLog
@@ -579,7 +597,7 @@ function test_barrier(C::Type{<:Cones.HypoPerLog})
     test_barrier(C(4), barrier)
 end
 
-get_allocs(C::Type{<:Cones.HypoPerLog}) = get_allocs(C(9))
+show_time_alloc(C::Type{<:Cones.HypoPerLog}) = show_time_alloc(C(9))
 
 
 # HypoPerLogdetTri
@@ -602,7 +620,7 @@ function test_barrier(C::Type{Cones.HypoPerLogdetTri{T, R}}) where {T, R}
     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)))
+show_time_alloc(C::Type{Cones.HypoPerLogdetTri{T, R}}) where {T, R} = show_time_alloc(C(2 + dim_herm(3, R)))
 
 
 # EpiPerEntropy
@@ -623,7 +641,7 @@ function test_barrier(C::Type{<:Cones.EpiPerEntropy})
     test_barrier(C(4), barrier)
 end
 
-get_allocs(C::Type{<:Cones.EpiPerEntropy}) = get_allocs(C(9))
+show_time_alloc(C::Type{<:Cones.EpiPerEntropy}) = show_time_alloc(C(9))
 
 
 # EpiPerTraceEntropyTri
@@ -646,7 +664,7 @@ function test_barrier(C::Type{Cones.EpiPerTraceEntropyTri{T}}) where T
     test_barrier(C(2 + Cones.svec_length(dW)), barrier)
 end
 
-get_allocs(C::Type{<:Cones.EpiPerTraceEntropyTri}) = get_allocs(C(8))
+show_time_alloc(C::Type{<:Cones.EpiPerTraceEntropyTri}) = show_time_alloc(C(8))
 
 
 # EpiRelEntropy
@@ -660,14 +678,15 @@ function test_oracles(C::Type{<:Cones.EpiRelEntropy})
 end
 
 function test_barrier(C::Type{<:Cones.EpiRelEntropy})
+    dw = 2
     function barrier(s)
-        (u, v, w) = (s[1], s[2:2:(end - 1)], s[3:2:end])
+        (u, v, w) = (s[1], s[2:(1 + dw)], s[(2 + dw):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))
+show_time_alloc(C::Type{<:Cones.EpiRelEntropy}) = show_time_alloc(C(9))
 
 
 # EpiTraceRelEntropyTri
@@ -692,7 +711,7 @@ function test_barrier(C::Type{Cones.EpiTraceRelEntropyTri{T}}) where T
     test_barrier(C(1 + 2 * dw), barrier, tol = 1e8 * eps(T))
 end
 
-get_allocs(C::Type{<:Cones.EpiTraceRelEntropyTri}) = get_allocs(C(13))
+show_time_alloc(C::Type{<:Cones.EpiTraceRelEntropyTri}) = show_time_alloc(C(13))
 
 
 # WSOSInterpNonnegative
@@ -710,7 +729,7 @@ function test_barrier(C::Type{Cones.WSOSInterpNonnegative{T, R}}) where {T, R}
     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)...))
+show_time_alloc(C::Type{Cones.WSOSInterpNonnegative{T, R}}) where {T, R} = show_time_alloc(C(rand_interp(3, 1, R)...))
 
 
 # WSOSInterpPosSemidefTri
@@ -743,7 +762,7 @@ function test_barrier(C::Type{Cones.WSOSInterpPosSemidefTri{T}}) where T
     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)...))
+show_time_alloc(C::Type{Cones.WSOSInterpPosSemidefTri{T}}) where T = show_time_alloc(C(2, rand_interp(2, 1, T)...))
 
 
 # WSOSInterpEpiNormEucl
@@ -774,7 +793,7 @@ function test_barrier(C::Type{Cones.WSOSInterpEpiNormEucl{T}}) where T
     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)...))
+show_time_alloc(C::Type{Cones.WSOSInterpEpiNormEucl{T}}) where T = show_time_alloc(C(3, rand_interp(2, 1, T)...))
 
 
 # WSOSInterpEpiNormOne
@@ -803,4 +822,4 @@ function test_barrier(C::Type{Cones.WSOSInterpEpiNormOne{T}}) where T
     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)...))
+show_time_alloc(C::Type{Cones.WSOSInterpEpiNormOne{T}}) where T = show_time_alloc(C(3, rand_interp(2, 1, T)...))
diff --git a/test/nativeinstances.jl b/test/nativeinstances.jl
index acb7e4104..0bff7910b 100644
--- a/test/nativeinstances.jl
+++ b/test/nativeinstances.jl
@@ -2048,9 +2048,9 @@ function epirelentropy1(T; options...)
         G = zeros(T, dim, 1)
         G[1, 1] = -1
         h = zeros(T, dim)
-        h[2:2:(end - 1)] .= 1
+        h[2:(1 + w_dim)] .= 1
         w = rand(T, w_dim) .+ 1
-        h[3:2:end] .= w
+        h[(2 + w_dim):end] .= w
         cones = Cone{T}[Cones.EpiRelEntropy{T}(dim)]
 
         r = build_solve_check(c, A, b, G, h, cones, tol; options...)
@@ -2068,10 +2068,10 @@ function epirelentropy2(T; options...)
         b = zeros(T, 0)
         G = zeros(T, dim, w_dim)
         for i in 1:w_dim
-            G[2i + 1, i] = -1
+            G[1 + w_dim + i, i] = -1
         end
         h = zeros(T, dim)
-        h[2:2:(dim - 1)] .= 1
+        h[1 .+ (1:w_dim)] .= 1
         cones = Cone{T}[Cones.EpiRelEntropy{T}(dim)]
 
         r = build_solve_check(c, A, b, G, h, cones, tol; options...)
@@ -2090,10 +2090,10 @@ function epirelentropy3(T; options...)
         b = T[dim]
         G = zeros(T, dim, w_dim)
         for i in 1:w_dim
-            G[2i, i] = -1
+            G[1 + i, i] = -1
         end
         h = zeros(T, dim)
-        h[3:2:end] .= 1
+        h[(2 + w_dim):end] .= 1
         cones = Cone{T}[Cones.EpiRelEntropy{T}(dim)]
 
         r = build_solve_check(c, A, b, G, h, cones, tol; options...)
@@ -2109,15 +2109,15 @@ function epirelentropy4(T; options...)
     A = zeros(T, 0, 1)
     b = zeros(T, 0)
     G = Matrix{T}(-I, 5, 1)
-    h = T[0, 1, 2, 5, 3]
+    h = T[0, 1, 5, 2, 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
+    @test r.s ≈ [entr, 1, 5, 2, 3] atol=tol rtol=tol
+    @test r.z ≈ [1, 2, 3 / T(5), log(inv(T(2))) - 1, log(5 / T(3)) - 1] atol=tol rtol=tol
 end
 
 function epirelentropy5(T; options...)
@@ -2125,15 +2125,16 @@ function epirelentropy5(T; options...)
     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]
+    G = vcat(zeros(T, 4, 2), fill(-one(T), 3, 2), [-1, 0]')
+    h = zeros(T, 8)
+    h[2:4] .= 1
     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.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
+    @test r.z ≈ inv(T(3)) * [1, 1, 1, 1, -1, -1, -1, 3] atol=tol rtol=tol
 end
 
 function epitracerelentropytri1(T; options...)
diff --git a/test/runconetests.jl b/test/runconetests.jl
index 0d30bef86..11b155425 100644
--- a/test/runconetests.jl
+++ b/test/runconetests.jl
@@ -88,8 +88,7 @@ real_types = [
     ]
 @testset "$cone" for T in real_types, cone in cone_types(T)
     println("\n$cone")
-    test_time = @elapsed allocs = get_allocs(cone)
-    display(allocs)
+    test_time = @elapsed show_time_alloc(cone)
     @printf("%8.2e seconds\n", test_time)
 end
 println()

From 5d97fe774d4135658baed8bb106677219c30b8cc Mon Sep 17 00:00:00 2001
From: Chris Coey <coey.chris@gmail.com>
Date: Thu, 1 Apr 2021 15:39:53 +1000
Subject: [PATCH 3/6] update infinity norm cone oracles

---
 src/Cones/Cones.jl         | 210 ++++++++++--------------------------
 src/Cones/epinorminf.jl    | 215 ++++++++++++++++++++-----------------
 src/Cones/epirelentropy.jl |   1 +
 3 files changed, 170 insertions(+), 256 deletions(-)

diff --git a/src/Cones/Cones.jl b/src/Cones/Cones.jl
index 859031ca3..b866cf619 100644
--- a/src/Cones/Cones.jl
+++ b/src/Cones/Cones.jl
@@ -346,8 +346,14 @@ vec_copy_to!(v1::AbstractVecOrMat{Complex{T}}, v2::AbstractVecOrMat{T}) where {T
 # utilities for hessians for cones with PSD parts
 
 # TODO parallelize
-function symm_kron(H::AbstractMatrix{T}, mat::AbstractMatrix{T}, rt2::T; upper_only::Bool = true) where {T <: Real}
+function symm_kron(
+    H::AbstractMatrix{T},
+    mat::AbstractMatrix{T},
+    rt2::T;
+    upper_only::Bool = true,
+    ) where {T <: Real}
     side = size(mat, 1)
+
     col_idx = 1
     @inbounds for l in 1:side
         for k in 1:(l - 1)
@@ -364,6 +370,7 @@ function symm_kron(H::AbstractMatrix{T}, mat::AbstractMatrix{T}, rt2::T; upper_o
             end
             col_idx += 1
         end
+
         row_idx = 1
         for j in 1:side
             upper_only && row_idx > col_idx && continue
@@ -376,12 +383,18 @@ function symm_kron(H::AbstractMatrix{T}, mat::AbstractMatrix{T}, rt2::T; upper_o
         end
         col_idx += 1
     end
+
     return H
 end
 
 # TODO test output for non-Hermitian mat, the result may need transposing
-function symm_kron(H::AbstractMatrix{T}, mat::AbstractMatrix{Complex{T}}, rt2::T) where {T <: Real}
+function symm_kron(
+    H::AbstractMatrix{T},
+    mat::AbstractMatrix{Complex{T}},
+    rt2::T,
+    ) where {T <: Real}
     side = size(mat, 1)
+
     col_idx = 1
     for i in 1:side, j in 1:i
         row_idx = 1
@@ -428,15 +441,28 @@ function symm_kron(H::AbstractMatrix{T}, mat::AbstractMatrix{Complex{T}}, rt2::T
             col_idx += 2
         end
     end
+
     return H
 end
 
-function hess_element(H::Matrix{T}, r_idx::Int, c_idx::Int, term1::T, term2::T) where {T <: Real}
+function hess_element(
+    H::Matrix{T},
+    r_idx::Int,
+    c_idx::Int,
+    term1::T,
+    term2::T,
+    ) where {T <: Real}
     @inbounds H[r_idx, c_idx] = term1 + term2
     return
 end
 
-function hess_element(H::Matrix{T}, r_idx::Int, c_idx::Int, term1::Complex{T}, term2::Complex{T}) where {T <: Real}
+function hess_element(
+    H::Matrix{T},
+    r_idx::Int,
+    c_idx::Int,
+    term1::Complex{T},
+    term2::Complex{T},
+    ) where {T <: Real}
     @inbounds begin
         H[r_idx, c_idx] = real(term1) + real(term2)
         H[r_idx + 1, c_idx] = imag(term2) - imag(term1)
@@ -446,165 +472,28 @@ function hess_element(H::Matrix{T}, r_idx::Int, c_idx::Int, term1::Complex{T}, t
     return
 end
 
-# set up sparse arrow matrix data structure
-# TODO remove this in favor of new hess_nz_count etc functions that directly use uu, uw, ww etc
-function sparse_arrow(T::Type{<:Real}, w_dim::Int)
-    dim = w_dim + 1
-    nnz_tri = 2 * dim - 1
-    I = Vector{Int}(undef, nnz_tri)
-    J = Vector{Int}(undef, nnz_tri)
-    idxs1 = 1:dim
-    @views I[idxs1] .= 1
-    @views J[idxs1] .= idxs1
-    idxs2 = (dim + 1):(2 * dim - 1)
-    @views I[idxs2] .= 2:dim
-    @views J[idxs2] .= 2:dim
-    V = ones(T, nnz_tri)
-    return sparse(I, J, V, dim, dim)
-end
-
-# 2x2 block case
-function sparse_arrow_block2(T::Type{<:Real}, w_dim::Int)
-    dim = 2 * w_dim + 1
-    nnz_tri = 2 * dim - 1 + w_dim
-    I = Vector{Int}(undef, nnz_tri)
-    J = Vector{Int}(undef, nnz_tri)
-    idxs1 = 1:dim
-    @views I[idxs1] .= 1
-    @views J[idxs1] .= idxs1
-    idxs2 = (dim + 1):(2 * dim - 1)
-    @views I[idxs2] .= 2:dim
-    @views J[idxs2] .= 2:dim
-    idxs3 = (2 * dim):nnz_tri
-    @views I[idxs3] .= 2:2:dim
-    @views J[idxs3] .= 3:2:dim
-    V = ones(T, nnz_tri)
-    return sparse(I, J, V, dim, dim)
-end
-
-# lmul with arrow matrix
-function arrow_prod(prod::AbstractVecOrMat{T}, arr::AbstractVecOrMat{T}, uu::T, uw::Vector{T}, ww::Vector{T}) where {T <: Real}
-    @inbounds @views begin
-        arru = arr[1, :]
-        arrw = arr[2:end, :]
-        produ = prod[1, :]
-        prodw = prod[2:end, :]
-        copyto!(produ, arru)
-        mul!(produ, arrw', uw, true, uu)
-        mul!(prodw, uw, arru')
-        @. prodw += ww * arrw
-    end
-    return prod
-end
-
-# 2x2 block case
-function arrow_prod(prod::AbstractVecOrMat{T}, arr::AbstractVecOrMat{T}, uu::T, uv::Vector{T}, uw::Vector{T}, vv::Vector{T}, vw::Vector{T}, ww::Vector{T}) where {T <: Real}
-    @inbounds @views begin
-        arru = arr[1, :]
-        arrv = arr[2:2:end, :]
-        arrw = arr[3:2:end, :]
-        produ = prod[1, :]
-        prodv = prod[2:2:end, :]
-        prodw = prod[3:2:end, :]
-        @. produ = uu * arru
-        mul!(produ, arrv', uv, true, true)
-        mul!(produ, arrw', uw, true, true)
-        mul!(prodv, uv, arru')
-        mul!(prodw, uw, arru')
-        @. prodv += vv * arrv + vw * arrw
-        @. prodw += ww * arrw + vw * arrv
-    end
-    return prod
-end
-
-# factorize arrow matrix
-function arrow_sqrt(uu::T, uw::Vector{T}, ww::Vector{T}, rtuw::Vector{T}, rtww::Vector{T}) where {T <: Real}
-    tol = sqrt(eps(T))
-    any(<(tol), ww) && return zero(T)
-    @. rtww = sqrt(ww)
-    @. rtuw = uw / rtww
-    diff = uu - sum(abs2, rtuw)
-    (diff < tol) && return zero(T)
-    return sqrt(diff)
-end
-
-# 2x2 block case
-function arrow_sqrt(uu::T, uv::Vector{T}, uw::Vector{T}, vv::Vector{T}, vw::Vector{T}, ww::Vector{T}, rtuv::Vector{T}, rtuw::Vector{T}, rtvv::Vector{T}, rtvw::Vector{T}, rtww::Vector{T}) where {T <: Real}
-    tol = sqrt(eps(T))
-    any(<(tol), ww) && return zero(T)
-    @. rtww = sqrt(ww)
-    @. rtvw = vw / rtww
-    @. rtuw = uw / rtww
-    @. rtuv = vv - abs2(rtvw)
-    any(<(tol), rtuv) && return zero(T)
-    @. rtvv = sqrt(rtuv)
-    @. rtuv = (uv - rtvw * rtuw) / rtvv
-    diff = uu - sum(abs2, rtuv) - sum(abs2, rtuw)
-    (diff < tol) && return zero(T)
-    return sqrt(diff)
-end
-
-# lmul with lower Cholesky factor of arrow matrix
-function arrow_sqrt_prod(prod::AbstractVecOrMat{T}, arr::AbstractVecOrMat{T}, rtuu::T, rtuw::Vector{T}, rtww::Vector{T}) where {T <: Real}
-    @inbounds @views begin
-        arr1 = arr[1, :]
-        @. prod[1, :] = rtuu * arr1
-        @. prod[2:end, :] = rtuw * arr1' + rtww * arr[2:end, :]
-    end
-    return prod
-end
-
-# 2x2 block case
-function arrow_sqrt_prod(prod::AbstractVecOrMat{T}, arr::AbstractVecOrMat{T}, rtuu::T, rtuv::Vector{T}, rtuw::Vector{T}, rtvv::Vector{T}, rtvw::Vector{T}, rtww::Vector{T}) where {T <: Real}
-    @inbounds @views begin
-        arr1 = arr[1, :]
-        arrv = arr[2:2:end, :]
-        @. prod[1, :] = rtuu * arr1
-        @. prod[2:2:end, :] = rtuv * arr1' + rtvv * arrv
-        @. prod[3:2:end, :] = rtuw * arr1' + rtvw * arrv + rtww * arr[3:2:end, :]
-    end
-    return prod
-end
-
-# ldiv with upper Cholesky factor of arrow matrix
-function inv_arrow_sqrt_prod(prod::AbstractVecOrMat{T}, arr::AbstractVecOrMat{T}, rtuu::T, rtuw::Vector{T}, rtww::Vector{T}) where {T <: Real}
-    @inbounds @. @views prod[2:end, :] = arr[2:end, :] / rtww
-    @inbounds @views for j in 1:size(arr, 2)
-        prod[1, j] = (arr[1, j] - dot(prod[2:end, j], rtuw)) / rtuu
-    end
-    return prod
-end
-
-# 2x2 block case
-function inv_arrow_sqrt_prod(prod::AbstractVecOrMat{T}, arr::AbstractVecOrMat{T}, rtuu::T, rtuv::Vector{T}, rtuw::Vector{T}, rtvv::Vector{T}, rtvw::Vector{T}, rtww::Vector{T}) where {T <: Real}
-    @inbounds @views begin
-        prodw = prod[3:2:end, :]
-        @. prodw = arr[3:2:end, :] / rtww
-        @. prod[2:2:end, :] = (arr[2:2:end, :] - rtvw * prodw) / rtvv
-    end
-    @inbounds @views for j in 1:size(arr, 2)
-        prod[1, j] = (arr[1, j] - dot(prod[2:2:end, j], rtuv) - dot(prod[3:2:end, j], rtuw)) / rtuu
-    end
-    return prod
-end
-
 function grad_logm!(
-        mat::Matrix{T},
-        vecs::Matrix{T},
-        tempmat1::Matrix{T},
-        tempmat2::Matrix{T},
-        tempvec::Vector{T},
-        diff_mat::AbstractMatrix{T},
-        rt2::T,
-        ) where T
+    mat::Matrix{T},
+    vecs::Matrix{T},
+    tempmat1::Matrix{T},
+    tempmat2::Matrix{T},
+    tempvec::Vector{T},
+    diff_mat::AbstractMatrix{T},
+    rt2::T,
+    ) where {T <: Real}
     veckron = symm_kron(tempmat1, vecs, rt2, upper_only = false)
     smat_to_svec!(tempvec, diff_mat, one(T))
     mul!(tempmat2, veckron, Diagonal(tempvec))
     return mul!(mat, tempmat2, veckron')
 end
 
-function diff_mat!(mat::Matrix{T}, vals::Vector{T}, log_vals::Vector{T}) where T
+function diff_mat!(
+    mat::Matrix{T},
+    vals::Vector{T},
+    log_vals::Vector{T},
+    ) where {T <: Real}
     rteps = sqrt(eps(T))
+
     @inbounds for j in eachindex(vals)
         (vj, lvj) = (vals[j], log_vals[j])
         for i in 1:(j - 1)
@@ -613,14 +502,21 @@ function diff_mat!(mat::Matrix{T}, vals::Vector{T}, log_vals::Vector{T}) where T
         end
         mat[j, j] = inv(vj)
     end
+
     return mat
 end
 
-function diff_tensor!(diff_tensor::Array{T, 3}, diff_mat::AbstractMatrix{T}, vals::Vector{T}) where T
+function diff_tensor!(
+    diff_tensor::Array{T, 3},
+    diff_mat::AbstractMatrix{T},
+    vals::Vector{T},
+    ) where {T <: Real}
     rteps = sqrt(eps(T))
     d = size(diff_mat, 1)
+
     @inbounds for k in 1:d, j in 1:k, i in 1:j
         (vi, vj, vk) = (vals[i], vals[j], vals[k])
+
         if abs(vj - vk) < rteps
             if abs(vi - vj) < rteps
                 vijk = (vi + vj + vk) / 3
@@ -632,9 +528,11 @@ function diff_tensor!(diff_tensor::Array{T, 3}, diff_mat::AbstractMatrix{T}, val
         else
             t = (diff_mat[i, j] - diff_mat[i, k]) / (vj - vk)
         end
+
         diff_tensor[i, j, k] = diff_tensor[i, k, j] = diff_tensor[j, i, k] =
             diff_tensor[j, k, i] = diff_tensor[k, i, j] = diff_tensor[k, j, i] = t
     end
+
     return diff_tensor
 end
 
diff --git a/src/Cones/epinorminf.jl b/src/Cones/epinorminf.jl
index dba8f5efe..e004582a5 100644
--- a/src/Cones/epinorminf.jl
+++ b/src/Cones/epinorminf.jl
@@ -24,8 +24,6 @@ mutable struct EpiNormInf{T <: Real, R <: RealOrComplex{T}} <: Cone{T}
     hess_updated::Bool
     inv_hess_updated::Bool
     hess_aux_updated::Bool
-    hess_sqrt_aux_updated::Bool
-    use_hess_sqrt::Bool
     inv_hess_aux_updated::Bool
     is_feas::Bool
     hess::Symmetric{T, SparseMatrixCSC{T, Int}}
@@ -41,15 +39,9 @@ mutable struct EpiNormInf{T <: Real, R <: RealOrComplex{T}} <: Cone{T}
     Hrere::Vector{T}
     Hreim::Vector{T}
     Himim::Vector{T}
-    rtuim::Vector{T}
-    rtreim::Vector{T}
-    rtimim::Vector{T}
     Hiure::Vector{T}
     Hiuim::Vector{T}
     schur::T
-    rtuu::T
-    rture::Vector{T}
-    rtrere::Vector{T}
     idet::Vector{T}
 
     function EpiNormInf{T, R}(
@@ -68,22 +60,9 @@ end
 
 use_heuristic_neighborhood(cone::EpiNormInf) = false
 
-reset_data(cone::EpiNormInf) = (cone.feas_updated = cone.grad_updated = cone.hess_updated = cone.inv_hess_updated = cone.hess_aux_updated = cone.hess_sqrt_aux_updated = cone.inv_hess_aux_updated = false)
-
-function use_sqrt_oracles(cone::EpiNormInf{T}) where T
-    cone.hess_sqrt_aux_updated || update_hess_sqrt_aux(cone)
-    cone.use_hess_sqrt || return false
-    # check numerics
-    # TODO tune
-    tol = sqrt(sqrt(eps(T))) * dimension(cone)
-    nu = get_nu(cone)
-    vec1 = cone.vec1
-    hess_sqrt_prod!(vec1, cone.point, cone)
-    hess_viol = abs(1 - sum(abs2, vec1) / nu)
-    inv_hess_sqrt_prod!(vec1, grad(cone), cone)
-    inv_hess_viol = abs(1 - sum(abs2, vec1) / nu)
-    return (max(hess_viol, inv_hess_viol) < tol)
-end
+reset_data(cone::EpiNormInf) = (cone.feas_updated = cone.grad_updated = cone.hess_updated = cone.inv_hess_updated = cone.hess_aux_updated = cone.inv_hess_aux_updated = false)
+
+use_sqrt_oracles(cone::EpiNormInf) = false # TODO remove in favor of BHB oracles
 
 # TODO only allocate the fields we use
 function setup_extra_data(cone::EpiNormInf{T, R}) where {R <: RealOrComplex{T}} where {T <: Real}
@@ -95,19 +74,13 @@ function setup_extra_data(cone::EpiNormInf{T, R}) where {R <: RealOrComplex{T}}
     cone.Hure = zeros(T, n)
     cone.Hrere = zeros(T, n)
     cone.Hiure = zeros(T, n)
-    cone.rture = zeros(T, n)
-    cone.rtrere = zeros(T, n)
     if cone.is_complex
         cone.Huim = zeros(T, n)
         cone.Hreim = zeros(T, n)
         cone.Himim = zeros(T, n)
-        cone.rtuim = zeros(T, n)
-        cone.rtreim = zeros(T, n)
-        cone.rtimim = zeros(T, n)
         cone.Hiuim = zeros(T, n)
         cone.idet = zeros(T, n)
     end
-    cone.use_hess_sqrt = true
     return cone
 end
 
@@ -191,44 +164,113 @@ function update_hess_aux(cone::EpiNormInf{T}) where {T <: Real}
     return
 end
 
-function update_hess(cone::EpiNormInf{T}) where {T <: Real}
+function update_hess(cone::EpiNormInf{T, T}) where {T <: Real}
     cone.hess_aux_updated || update_hess_aux(cone)
+    dim = cone.dim
 
     if !isdefined(cone, :hess)
         # initialize sparse idxs for upper triangle of Hessian
-        spfun = (cone.is_complex ? sparse_arrow_block2 : sparse_arrow)
-        cone.hess = Symmetric(spfun(T, cone.n), :U)
+        nnz_tri = 2 * dim - 1
+        I = Vector{Int}(undef, nnz_tri)
+        J = Vector{Int}(undef, nnz_tri)
+        idxs1 = 1:dim
+        @views I[idxs1] .= 1
+        @views J[idxs1] .= idxs1
+        idxs2 = (dim + 1):(2 * dim - 1)
+        @views I[idxs2] .= 2:dim
+        @views J[idxs2] .= 2:dim
+        V = ones(T, nnz_tri)
+
+        cone.hess = Symmetric(sparse(I, J, V, dim, dim), :U)
     end
 
     # modify nonzeros of upper triangle of Hessian
     nzval = cone.hess.data.nzval
     nzval[1] = cone.Huu
-    if cone.is_complex
-        nz_idx = 1
-        @inbounds for i in 1:cone.n
-            @. nzval[nz_idx .+ (1:5)] = (cone.Hure[i], cone.Hrere[i], cone.Huim[i], cone.Hreim[i], cone.Himim[i])
-            nz_idx += 5
-        end
-    else
-        nz_idx = 2
-        @inbounds for i in 1:cone.n
-            nzval[nz_idx] = cone.Hure[i]
-            nzval[nz_idx + 1] = cone.Hrere[i]
-            nz_idx += 2
-        end
+    nz_idx = 2
+    @inbounds for i in 1:cone.n
+        nzval[nz_idx] = cone.Hure[i]
+        nzval[nz_idx + 1] = cone.Hrere[i]
+        nz_idx += 2
     end
 
     cone.hess_updated = true
     return cone.hess
 end
 
-function hess_prod!(prod::AbstractVecOrMat, arr::AbstractVecOrMat, cone::EpiNormInf)
+function update_hess(cone::EpiNormInf{T, Complex{T}}) where {T <: Real}
     cone.hess_aux_updated || update_hess_aux(cone)
-    if cone.is_complex
-        return arrow_prod(prod, arr, cone.Huu, cone.Hure, cone.Huim, cone.Hrere, cone.Hreim, cone.Himim)
-    else
-        return arrow_prod(prod, arr, cone.Huu, cone.Hure, cone.Hrere)
+    dim = cone.dim
+
+    if !isdefined(cone, :hess)
+        # initialize sparse idxs for upper triangle of Hessian
+        nnz_tri = 2 * dim - 1 + cone.n
+        I = Vector{Int}(undef, nnz_tri)
+        J = Vector{Int}(undef, nnz_tri)
+        idxs1 = 1:dim
+        @views I[idxs1] .= 1
+        @views J[idxs1] .= idxs1
+        idxs2 = (dim + 1):(2 * dim - 1)
+        @views I[idxs2] .= 2:dim
+        @views J[idxs2] .= 2:dim
+        idxs3 = (2 * dim):nnz_tri
+        @views I[idxs3] .= 2:2:dim
+        @views J[idxs3] .= 3:2:dim
+        V = ones(T, nnz_tri)
+
+        cone.hess = Symmetric(sparse(I, J, V, dim, dim), :U)
+    end
+
+    # modify nonzeros of upper triangle of Hessian
+    nzval = cone.hess.data.nzval
+    nzval[1] = cone.Huu
+    nz_idx = 1
+    @inbounds for i in 1:cone.n
+        @. nzval[nz_idx .+ (1:5)] = (cone.Hure[i], cone.Hrere[i], cone.Huim[i], cone.Hreim[i], cone.Himim[i])
+        nz_idx += 5
+    end
+
+    cone.hess_updated = true
+    return cone.hess
+end
+
+function hess_prod!(prod::AbstractVecOrMat{T}, arr::AbstractVecOrMat{T}, cone::EpiNormInf{T, T}) where {T <: Real}
+    cone.hess_aux_updated || update_hess_aux(cone)
+
+    @inbounds @views begin
+        u_arr = arr[1, :]
+        w_arr = arr[2:end, :]
+        u_prod = prod[1, :]
+        w_prod = prod[2:end, :]
+        copyto!(u_prod, u_arr)
+        mul!(u_prod, w_arr', cone.Hure, true, cone.Huu)
+        mul!(w_prod, cone.Hure, u_arr')
+        @. w_prod += cone.Hrere * w_arr
     end
+
+    return prod
+end
+
+function hess_prod!(prod::AbstractVecOrMat{T}, arr::AbstractVecOrMat{T}, cone::EpiNormInf{T, Complex{T}}) where {T <: Real}
+    cone.hess_aux_updated || update_hess_aux(cone)
+
+    @inbounds @views begin
+        u_arr = arr[1, :]
+        v_arr = arr[2:2:end, :]
+        w_arr = arr[3:2:end, :]
+        u_prod = prod[1, :]
+        v_prod = prod[2:2:end, :]
+        w_prod = prod[3:2:end, :]
+        @. u_prod = cone.Huu * u_arr
+        mul!(u_prod, v_arr', cone.Hure, true, true)
+        mul!(u_prod, w_arr', cone.Huim, true, true)
+        mul!(v_prod, cone.Hure, u_arr')
+        mul!(w_prod, cone.Huim, u_arr')
+        @. v_prod += cone.Hrere * v_arr + cone.Hreim * w_arr
+        @. w_prod += cone.Hreim * v_arr + cone.Himim * w_arr
+    end
+
+    return prod
 end
 
 function update_inv_hess_aux(cone::EpiNormInf{T}) where {T <: Real}
@@ -307,66 +349,39 @@ function update_inv_hess(cone::EpiNormInf{T}) where {T <: Real}
     return cone.inv_hess
 end
 
-function inv_hess_prod!(prod::AbstractVecOrMat, arr::AbstractVecOrMat, cone::EpiNormInf)
+function inv_hess_prod!(prod::AbstractVecOrMat{T}, arr::AbstractVecOrMat{T}, cone::EpiNormInf{T, T}) where {T <: Real}
     cone.inv_hess_aux_updated || update_inv_hess_aux(cone)
 
-    if cone.is_complex
-        @inbounds @views begin
-            copyto!(prod[1, :], arr[1, :])
-            mul!(prod[1, :], arr[2:2:end, :]', cone.Hiure, true, true)
-            mul!(prod[1, :], arr[3:2:end, :]', cone.Hiuim, true, true)
-            @. prod[2:2:end, :] = cone.Hiure * prod[1, :]'
-            @. prod[3:2:end, :] = cone.Hiuim * prod[1, :]'
-            prod ./= cone.schur
-        end
-
-        @inbounds @views for j in 1:cone.n
-            j2 = 2j
-            @. prod[j2, :] += (cone.Himim[j] * arr[j2, :] - cone.Hreim[j] * arr[j2 + 1, :]) / cone.idet[j]
-            @. prod[j2 + 1, :] += (cone.Hrere[j] * arr[j2 + 1, :] - cone.Hreim[j] * arr[j2, :]) / cone.idet[j]
-        end
-    else
-        @inbounds @views begin
-            copyto!(prod[1, :], arr[1, :])
-            mul!(prod[1, :], arr[2:end, :]', cone.Hiure, true, true)
-            @. prod[2:end, :] = cone.Hiure * prod[1, :]'
-            prod ./= cone.schur
-            @. prod[2:end, :] += arr[2:end, :] / cone.Hrere
-        end
+    @inbounds @views begin
+        copyto!(prod[1, :], arr[1, :])
+        mul!(prod[1, :], arr[2:end, :]', cone.Hiure, true, true)
+        @. prod[2:end, :] = cone.Hiure * prod[1, :]'
+        prod ./= cone.schur
+        @. prod[2:end, :] += arr[2:end, :] / cone.Hrere
     end
 
     return prod
 end
 
-function update_hess_sqrt_aux(cone::EpiNormInf)
-    cone.hess_aux_updated || update_hess_aux(cone)
-    @assert !cone.hess_sqrt_aux_updated
-    if cone.is_complex
-        cone.rtuu = arrow_sqrt(cone.Huu, cone.Hure, cone.Huim, cone.Hrere, cone.Hreim, cone.Himim, cone.rture, cone.rtuim, cone.rtrere, cone.rtreim, cone.rtimim)
-    else
-        cone.rtuu = arrow_sqrt(cone.Huu, cone.Hure, cone.Hrere, cone.rture, cone.rtrere)
-    end
-    cone.use_hess_sqrt = !iszero(cone.rtuu)
-    cone.hess_sqrt_aux_updated = true
-    return
-end
+function inv_hess_prod!(prod::AbstractVecOrMat{T}, arr::AbstractVecOrMat{T}, cone::EpiNormInf{T, Complex{T}}) where {T <: Real}
+    cone.inv_hess_aux_updated || update_inv_hess_aux(cone)
 
-function hess_sqrt_prod!(prod::AbstractVecOrMat, arr::AbstractVecOrMat, cone::EpiNormInf)
-    @assert cone.hess_sqrt_aux_updated && cone.use_hess_sqrt
-    if cone.is_complex
-        return arrow_sqrt_prod(prod, arr, cone.rtuu, cone.rture, cone.rtuim, cone.rtrere, cone.rtreim, cone.rtimim)
-    else
-        return arrow_sqrt_prod(prod, arr, cone.rtuu, cone.rture, cone.rtrere)
+    @inbounds @views begin
+        copyto!(prod[1, :], arr[1, :])
+        mul!(prod[1, :], arr[2:2:end, :]', cone.Hiure, true, true)
+        mul!(prod[1, :], arr[3:2:end, :]', cone.Hiuim, true, true)
+        @. prod[2:2:end, :] = cone.Hiure * prod[1, :]'
+        @. prod[3:2:end, :] = cone.Hiuim * prod[1, :]'
+        prod ./= cone.schur
     end
-end
 
-function inv_hess_sqrt_prod!(prod::AbstractVecOrMat, arr::AbstractVecOrMat, cone::EpiNormInf)
-    @assert cone.hess_sqrt_aux_updated && cone.use_hess_sqrt
-    if cone.is_complex
-        return inv_arrow_sqrt_prod(prod, arr, cone.rtuu, cone.rture, cone.rtuim, cone.rtrere, cone.rtreim, cone.rtimim)
-    else
-        return inv_arrow_sqrt_prod(prod, arr, cone.rtuu, cone.rture, cone.rtrere)
+    @inbounds @views for j in 1:cone.n
+        j2 = 2j
+        @. prod[j2, :] += (cone.Himim[j] * arr[j2, :] - cone.Hreim[j] * arr[j2 + 1, :]) / cone.idet[j]
+        @. prod[j2 + 1, :] += (cone.Hrere[j] * arr[j2 + 1, :] - cone.Hreim[j] * arr[j2, :]) / cone.idet[j]
     end
+
+    return prod
 end
 
 function correction(cone::EpiNormInf{T}, primal_dir::AbstractVector{T}) where {T <: Real}
diff --git a/src/Cones/epirelentropy.jl b/src/Cones/epirelentropy.jl
index 2f8d9c987..828ae8bd4 100644
--- a/src/Cones/epirelentropy.jl
+++ b/src/Cones/epirelentropy.jl
@@ -240,6 +240,7 @@ function update_inv_hess(cone::EpiRelEntropy{T}) where T
         @views I[idxs3] .= 1 .+ (1:w_dim)
         @views J[idxs3] .= (1 + w_dim) .+ (1:w_dim)
         V = ones(T, nnz_tri)
+
         cone.inv_hess = Symmetric(sparse(I, J, V, dim, dim), :U)
     end
 

From 6a5fcd5a72b92a2d40eb1f1e3849e282b6fdf259 Mon Sep 17 00:00:00 2001
From: Chris Coey <coey.chris@gmail.com>
Date: Thu, 1 Apr 2021 19:29:44 +1000
Subject: [PATCH 4/6] cleanup infinity norm oracles

---
 src/Cones/epinorminf.jl       | 43 ++++++++++++++++++-----------------
 src/MathOptInterface/cones.jl | 27 ----------------------
 test/runmoitests.jl           |  2 +-
 3 files changed, 23 insertions(+), 49 deletions(-)

diff --git a/src/Cones/epinorminf.jl b/src/Cones/epinorminf.jl
index e004582a5..ce5ca49a3 100644
--- a/src/Cones/epinorminf.jl
+++ b/src/Cones/epinorminf.jl
@@ -353,11 +353,12 @@ function inv_hess_prod!(prod::AbstractVecOrMat{T}, arr::AbstractVecOrMat{T}, con
     cone.inv_hess_aux_updated || update_inv_hess_aux(cone)
 
     @inbounds @views begin
-        copyto!(prod[1, :], arr[1, :])
-        mul!(prod[1, :], arr[2:end, :]', cone.Hiure, true, true)
-        @. prod[2:end, :] = cone.Hiure * prod[1, :]'
-        prod ./= cone.schur
-        @. prod[2:end, :] += arr[2:end, :] / cone.Hrere
+        w_arr = arr[2:end, :]
+        u_prod = prod[1, :]
+        copyto!(u_prod, arr[1, :])
+        mul!(u_prod, w_arr', cone.Hiure, true, true)
+        u_prod ./= cone.schur
+        @. prod[2:end, :] = cone.Hiure * u_prod' + w_arr / cone.Hrere
     end
 
     return prod
@@ -367,18 +368,15 @@ function inv_hess_prod!(prod::AbstractVecOrMat{T}, arr::AbstractVecOrMat{T}, con
     cone.inv_hess_aux_updated || update_inv_hess_aux(cone)
 
     @inbounds @views begin
-        copyto!(prod[1, :], arr[1, :])
-        mul!(prod[1, :], arr[2:2:end, :]', cone.Hiure, true, true)
-        mul!(prod[1, :], arr[3:2:end, :]', cone.Hiuim, true, true)
-        @. prod[2:2:end, :] = cone.Hiure * prod[1, :]'
-        @. prod[3:2:end, :] = cone.Hiuim * prod[1, :]'
-        prod ./= cone.schur
-    end
-
-    @inbounds @views for j in 1:cone.n
-        j2 = 2j
-        @. prod[j2, :] += (cone.Himim[j] * arr[j2, :] - cone.Hreim[j] * arr[j2 + 1, :]) / cone.idet[j]
-        @. prod[j2 + 1, :] += (cone.Hrere[j] * arr[j2 + 1, :] - cone.Hreim[j] * arr[j2, :]) / cone.idet[j]
+        u_prod = prod[1, :]
+        re_arr = arr[2:2:end, :]
+        im_arr = arr[3:2:end, :]
+        copyto!(u_prod, arr[1, :])
+        mul!(u_prod, re_arr', cone.Hiure, true, true)
+        mul!(u_prod, im_arr', cone.Hiuim, true, true)
+        u_prod ./= cone.schur
+        @. prod[2:2:end, :] = cone.Hiure * u_prod' + (cone.Himim * re_arr - cone.Hreim * im_arr) / cone.idet
+        @. prod[3:2:end, :] = cone.Hiuim * u_prod' + (cone.Hrere * im_arr - cone.Hreim * re_arr) / cone.idet
     end
 
     return prod
@@ -392,13 +390,15 @@ function correction(cone::EpiNormInf{T}, primal_dir::AbstractVector{T}) where {T
 
     u3 = T(1.5) / u
     udu = udir / u
-    corr1 = -udir * sum(z * (u3 - z) * z for z in cone.uden) * udir - udu * (cone.n - 1) / u * udu
+    corr[1] = -udir * sum(z * (u3 - z) * z for z in cone.uden) * udir - udu * (cone.n - 1) / u * udu
+
     @inbounds for i in 1:cone.n
         deni = -4 * cone.den[i]
         udeni = 2 * cone.uden[i]
         suuw = udir * (-1 + udeni * u)
         wi = cone.w[i]
         wdeni = 2 * cone.wden[i]
+
         if cone.is_complex
             (wdenire, wdeniim) = reim(wdeni)
             (wire, wiim) = reim(wi)
@@ -414,7 +414,8 @@ function correction(cone::EpiNormInf{T}, primal_dir::AbstractVector{T}) where {T
             imimwimimim = wdeniim * (2 + uimimim)
             imimwrereim = wdeniim * uimimre * dire
             imimwimimre = wdenire * uimimim * diim
-            corr1 += (2 * (uuwre * dire + uuwim * diim) + uimimrere * dire + uimimimim * diim + 2 * uimimimre * diim * dire) / deni
+
+            corr[1] += (2 * (uuwre * dire + uuwim * diim) + uimimrere * dire + uimimimim * diim + 2 * uimimimre * diim * dire) / deni
             corr[2i] = (udir * (2 * (uimimrere + uimimimre * diim) + uuwre) + (abs2(dire) * imimwrerere + diim * (2 * imimwrereim + imimwimimre))) / deni
             corr[2i + 1] = (udir * (2 * (uimimimim + uimimimre * dire) + uuwim) + (abs2(diim) * imimwimimim + dire * (2 * imimwimimre + imimwrereim))) / deni
         else
@@ -422,11 +423,11 @@ function correction(cone::EpiNormInf{T}, primal_dir::AbstractVector{T}) where {T
             uuw = suuw * wdeni
             uimim = 1 + wdeni * wi
             uimim2 = -udeni * uimim * di
-            corr1 += di * (2 * uuw + uimim2) / deni
+
+            corr[1] += di * (2 * uuw + uimim2) / deni
             corr[1 + i] = (udir * (uuw + 2 * uimim2) + di * wdeni * (2 + uimim) * di) / deni
         end
     end
-    corr[1] = corr1
 
     return corr
 end
diff --git a/src/MathOptInterface/cones.jl b/src/MathOptInterface/cones.jl
index 56fb0cb9c..07f2bf2f6 100644
--- a/src/MathOptInterface/cones.jl
+++ b/src/MathOptInterface/cones.jl
@@ -30,33 +30,6 @@ permute_affine(::MOI.AbstractVectorSet, idxs::AbstractVector) = idxs
 rescale_affine(::MOI.AbstractVectorSet, vals::AbstractVector) = vals
 rescale_affine(::MOI.AbstractVectorSet, vals::AbstractVector, ::AbstractVector) = vals
 
-# # transformation for MOI relative entropy cone
-# needs_untransform(::MOI.RelativeEntropyCone) = true
-#
-# function untransform_affine(::MOI.RelativeEntropyCone, vals::AbstractVector)
-#     w_dim = div(length(vals) - 1, 2)
-#     v_vals = vals[2:2:(end - 1)]
-#     w_vals = vals[3:2:end]
-#     vals[2:(w_dim + 1)] = v_vals
-#     vals[(w_dim + 2):end] = w_vals
-#     return vals
-# end
-#
-# function permute_affine(cone::MOI.RelativeEntropyCone, idxs::AbstractVector)
-#     dim = MOI.dimension(cone)
-#     w_dim = div(dim - 1, 2)
-#     new_idxs = collect(idxs)
-#     for (i, idx) in enumerate(idxs)
-#         idx <= 1 && continue
-#         if idx <= 1 + w_dim
-#             new_idxs[i] = 2 * (idx - 1)
-#         else
-#             new_idxs[i] = 2 * (idx - w_dim) - 1
-#         end
-#     end
-#     return new_idxs
-# end
-
 # transformations (transposition of matrix) for MOI rectangular matrix cones with matrix of more rows than columns
 
 NonSquareMatrixCone = Union{MOI.NormSpectralCone, MOI.NormNuclearCone}
diff --git a/test/runmoitests.jl b/test/runmoitests.jl
index b9cf02bcc..31dbf5a67 100644
--- a/test/runmoitests.jl
+++ b/test/runmoitests.jl
@@ -36,7 +36,7 @@ end
         default_options = (
             # verbose = true,
             verbose = false,
-            default_tol_relax = 3, # tune to make tests pass
+            default_tol_relax = 4,
             stepper = stepper{T}(),
             )
         test_info = "$T, $use_dense_model"

From 3f9523738235a76c674698be803749b627006116 Mon Sep 17 00:00:00 2001
From: Chris Coey <coey.chris@gmail.com>
Date: Thu, 1 Apr 2021 20:34:59 +1000
Subject: [PATCH 5/6] update alpha schedule

---
 src/Solvers/search.jl    | 2 +-
 test/runconetests.jl     | 2 +-
 test/runexamplestests.jl | 8 ++++----
 3 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/src/Solvers/search.jl b/src/Solvers/search.jl
index 8e75b296d..cd47153a8 100644
--- a/src/Solvers/search.jl
+++ b/src/Solvers/search.jl
@@ -21,7 +21,7 @@ mutable struct StepSearcher{T <: Real}
         step_searcher.cone_order = collect(1:length(cones))
         step_searcher.min_nbhd = T(0.01) # TODO tune
         step_searcher.max_nbhd = T(0.99) # TODO tune, maybe should be different for cones without third order correction
-        step_searcher.alpha_sched = T[0.9999, 0.999, 0.99, 0.97, 0.95, 0.9, 0.85, 0.8, 0.7, 0.6, 0.5, 0.3, 0.1] # TODO tune
+        step_searcher.alpha_sched = T[0.9999, 0.999, 0.99, 0.97, 0.95, 0.9, 0.85, 0.8, 0.7, 0.6, 0.5, 0.3, 0.1, 0.05, 0.01, 0.005, 0.001, 0.0005] # TODO tune
         step_searcher.prev_sched = 0
         return step_searcher
     end
diff --git a/test/runconetests.jl b/test/runconetests.jl
index 11b155425..6956922bf 100644
--- a/test/runconetests.jl
+++ b/test/runconetests.jl
@@ -79,7 +79,7 @@ real_types = [
 end
 end
 
-println("\nstarting allocation tests")
+println("\nstarting time/allocation measurements")
 @testset "allocation tests" begin
 real_types = [
     Float64,
diff --git a/test/runexamplestests.jl b/test/runexamplestests.jl
index 6652d6dca..c8cdca2db 100644
--- a/test/runexamplestests.jl
+++ b/test/runexamplestests.jl
@@ -14,8 +14,8 @@ script_verbose = false
 
 # default options to solvers
 default_options = (
-    verbose = false,
-    # verbose = true,
+    # verbose = false,
+    verbose = true,
     default_tol_relax = 10,
     # stepper = Solvers.CombinedStepper{Float64}(),
     # stepper = Solvers.PredOrCentStepper{Float64}(),
@@ -27,7 +27,7 @@ instance_sets = [
     ("minimal", Float64, 60),
     # ("minimal", Float32, 60),
     # ("minimal", BigFloat, 60),
-    # ("fast", Float64, 60),
+    ("fast", Float64, 60),
     # ("slow", Float64, 120),
     ]
 
@@ -100,7 +100,7 @@ for (inst_set, real_T, time_limit) in instance_sets
     new_default_options = (; default_options..., time_limit = time_limit)
     ex_type_T = ex_type{real_T}
 
-    println("\nstarting $ex_type_T $inst_set tests")
+    println("\nstarting $ex_type_T $inst_set tests\n")
     @testset "$ex_type_T $inst_set" begin
         run_instance_set(inst_subset, ex_type_T, info_perf, new_default_options, script_verbose, perf, results_path)
     end

From 355a749a092345631de5200802f01e2de215b8bf Mon Sep 17 00:00:00 2001
From: Chris Coey <coey.chris@gmail.com>
Date: Thu, 1 Apr 2021 21:12:13 +1000
Subject: [PATCH 6/6] fix output in examples tests

---
 test/runexamplestests.jl | 6 +++---
 test/runtests.jl         | 2 +-
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/test/runexamplestests.jl b/test/runexamplestests.jl
index c8cdca2db..19615b3a1 100644
--- a/test/runexamplestests.jl
+++ b/test/runexamplestests.jl
@@ -14,8 +14,8 @@ script_verbose = false
 
 # default options to solvers
 default_options = (
-    # verbose = false,
-    verbose = true,
+    verbose = false,
+    # verbose = true,
     default_tol_relax = 10,
     # stepper = Solvers.CombinedStepper{Float64}(),
     # stepper = Solvers.PredOrCentStepper{Float64}(),
@@ -27,7 +27,7 @@ instance_sets = [
     ("minimal", Float64, 60),
     # ("minimal", Float32, 60),
     # ("minimal", BigFloat, 60),
-    ("fast", Float64, 60),
+    # ("fast", Float64, 60),
     # ("slow", Float64, 120),
     ]
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 436953d54..a5c3f847c 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -24,9 +24,9 @@ println()
 all_test_time = @elapsed for t in test_files
     @info("starting $t tests")
     test_time = @elapsed include("run$(t)tests.jl")
+    flush(stdout); flush(stderr)
     @info("finished $t tests in $(@sprintf("%8.2e seconds", test_time))")
     println()
-    flush(stdout); flush(stderr)
 end
 @info("finished all tests in $(@sprintf("%8.2e seconds", all_test_time))")
 end