define QMSolution type (#92)

Author: geoffroyleconte

src/QuadraticModels.jl

@@ -30,7 +30,7 @@ import NLPModelsModifiers: SlackModel, slack_meta
 
 import Base.convert
 
-export AbstractQuadraticModel, QuadraticModel, presolve, postsolve, postsolve!
+export AbstractQuadraticModel, QuadraticModel, presolve, postsolve, postsolve!, QMSolution
 
 include("linalg_utils.jl")
 include("qpmodel.jl")

src/presolve/empty_rows.jl

@@ -22,4 +22,4 @@ function empty_rows!(
   return empty_row_pass
 end
 
-postsolve!(pt::OutputPoint, operation::EmptyRow) = nothing
+postsolve!(sol::QMSolution, operation::EmptyRow) = nothing

src/presolve/free_rows.jl

@@ -21,4 +21,4 @@ function free_rows!(
   return free_row_pass
 end
 
-postsolve!(pt::OutputPoint, operation::FreeRow) = nothing
+postsolve!(sol::QMSolution, operation::FreeRow) = nothing

src/presolve/linear_singleton_columns.jl

@@ -81,16 +81,16 @@ function free_linear_singleton_columns!(
   return free_lsc_pass, c0 + c0_offset
 end
 
-function postsolve!(pt::OutputPoint{T, S}, operation::FreeLinearSingletonColumn{T, S}) where {T, S}
-  x = pt.x
+function postsolve!(sol::QMSolution{T, S}, operation::FreeLinearSingletonColumn{T, S}) where {T, S}
+  x = sol.x
   j = operation.j
   # x[j] = (coival - Σₖ Aik x[k]) / Aij , where k ≂̸ j
   x[j] = operation.conival
   for (i, Aij) in zip(operation.arowi.nzind, operation.arowi.nzval)
     x[j] -= Aij * x[i]
   end
   x[j] /= operation.aij
-  pt.s_l[j] = zero(T)
-  pt.s_u[j] = zero(T)
-  pt.y[operation.i] = operation.yi
+  sol.s_l[j] = zero(T)
+  sol.s_u[j] = zero(T)
+  sol.y[operation.i] = operation.yi
 end

src/presolve/postsolve_utils.jl

@@ -1,22 +1,22 @@
-function restore_x!(kept_cols, xin::S, xout::S, nvar) where {S}
+function restore_x!(kept_cols, x_in::S, x::S, nvar) where {S}
   # put x and xps inside xout according to kept_cols
   cx = 0
   for i = 1:nvar
     if kept_cols[i]
       cx += 1
-      xout[i] = xin[cx]
+      x[i] = x_in[cx]
     end
   end
 end
 
-function restore_y!(kept_rows::Vector{Bool}, y::Vector{T}, yout::Vector{T}, ncon) where {T}
+function restore_y!(kept_rows::Vector{Bool}, y_in::Vector{T}, y::Vector{T}, ncon) where {T}
   c_y = 0
   for i = 1:ncon
     if !kept_rows[i]
-      yout[i] = zero(T)
+      y[i] = zero(T)
     else
       c_y += 1
-      yout[i] = y[c_y]
+      y[i] = y_in[c_y]
     end
   end
 end

src/presolve/presolve.jl

@@ -1,6 +1,13 @@
 abstract type PresolveOperation{T, S} end
 
-mutable struct OutputPoint{T, S}
+"""
+Type used to define a solution point when using [`postsolve`](@ref).
+
+    sol = QMSolution(x, y, s_l, s_u)
+
+where `s_l` and `s_u` should be of type `SparseVector`.
+"""
+mutable struct QMSolution{T, S}
   x::S
   y::S
   s_l::SparseVector{T, Int}
@@ -347,49 +354,46 @@ end
 function postsolve!(
   qm::QuadraticModel{T, S},
   psqm::PresolvedQuadraticModel{T, S},
-  pt_out::OutputPoint{T, S},
-  x_in::S,
-  y_in::S,
+  sol::QMSolution{T, S},
+  sol_in::QMSolution{T, S},
 ) where {T, S}
+  x_in, y_in = sol_in.x, sol_in.y
   n_operations = length(psqm.psd.operations)
-  nvar = length(pt_out.x)
-  restore_x!(psqm.psd.kept_cols, x_in, pt_out.x, nvar)
-  ncon = length(pt_out.y)
-  restore_y!(psqm.psd.kept_rows, y_in, pt_out.y, ncon)
+  nvar = length(sol.x)
+  restore_x!(psqm.psd.kept_cols, x_in, sol.x, nvar)
+  ncon = length(sol.y)
+  restore_y!(psqm.psd.kept_rows, y_in, sol.y, ncon)
   for i = n_operations:-1:1
     operation_i = psqm.psd.operations[i]
-    postsolve!(pt_out, operation_i)
+    postsolve!(sol, operation_i)
   end
 end
 
 """
-    pt = postsolve(qm::QuadraticModel{T, S}, psqm::PresolvedQuadraticModel{T, S}, 
-                   x_in::S, y_in::S,
-                   s_l_in::SparseVector{T, Int},
-                   s_u_in::SparseVector{T, Int}) where {T, S}
+    sol = postsolve(qm::QuadraticModel{T, S}, psqm::PresolvedQuadraticModel{T, S}, 
+                    sol_in::QMSolution{T, S}) where {T, S}
 
-Retrieve the solution `(x, y, s_l, s_u)` of the original QP `qm` given the solution of the presolved QP (`psqm`)
-`x_in, y_in, s_l_in, s_u_in`.
+Retrieve the solution `sol = (x, y, s_l, s_u)` of the original QP `qm` given the solution of the presolved QP (`psqm`)
+`sol_in` of type [`QMSolution`](@ref).
 """
 function postsolve(
   qm::QuadraticModel{T, S},
   psqm::PresolvedQuadraticModel{T, S},
-  x_in::S,
-  y_in::S,
-  s_l::SparseVector{T, Int},
-  s_u::SparseVector{T, Int},
+  sol_in::QMSolution{T, S},
 ) where {T, S}
-  x_out = similar(qm.meta.x0)
-  y_out = similar(qm.meta.y0)
+  x = similar(qm.meta.x0)
+  y = similar(qm.meta.y0)
+  s_l = sol_in.s_l
+  s_u = sol_in.s_u
 
   ilow, iupp = s_l.nzind, s_u.nzind
   restore_ilow_iupp!(ilow, iupp, psqm.psd.kept_cols)
-  pt_out = OutputPoint(
-    x_out,
-    y_out,
+  sol = QMSolution(
+    x,
+    y,
     SparseVector(qm.meta.nvar, ilow, s_l.nzval),
     SparseVector(qm.meta.nvar, iupp, s_u.nzval),
   )
-  postsolve!(qm, psqm, pt_out, x_in, y_in)
-  return pt_out
+  postsolve!(qm, psqm, sol, sol_in)
+  return sol
 end

src/presolve/remove_ifix.jl

@@ -87,15 +87,15 @@ function remove_ifix!(
   return c0ps, ifix_pass
 end
 
-function postsolve!(pt::OutputPoint{T, S}, operation::RemoveIfix{T, S}) where {T, S}
+function postsolve!(sol::QMSolution{T, S}, operation::RemoveIfix{T, S}) where {T, S}
   j = operation.j
-  pt.x[j] = operation.xj
-  ATyj = @views dot(operation.acolj.nzval, pt.y[operation.acolj.nzind])
-  Hxj = @views dot(operation.hcolj.nzval, pt.x[operation.hcolj.nzind])
+  sol.x[j] = operation.xj
+  ATyj = @views dot(operation.acolj.nzval, sol.y[operation.acolj.nzind])
+  Hxj = @views dot(operation.hcolj.nzval, sol.x[operation.hcolj.nzind])
   s = operation.cj + Hxj - ATyj
   if s > zero(T)
-    pt.s_l[j] = s
+    sol.s_l[j] = s
   else
-    pt.s_u[j] = -s
+    sol.s_u[j] = -s
   end
 end

src/presolve/singleton_rows.jl

@@ -63,19 +63,19 @@ function singleton_rows!(
   return singl_row_pass
 end
 
-function postsolve!(pt::OutputPoint{T, S}, operation::SingletonRow{T, S}) where {T, S}
+function postsolve!(sol::QMSolution{T, S}, operation::SingletonRow{T, S}) where {T, S}
   i, j = operation.i, operation.j
   Aij = operation.Aij
-  pt.y[i] = zero(T)
+  sol.y[i] = zero(T)
   if operation.tightened_lvar
-    pt.y[i] += pt.s_l[j] / Aij
-    pt.s_l[j] = zero(T)
+    sol.y[i] += sol.s_l[j] / Aij
+    sol.s_l[j] = zero(T)
   end
   if operation.tightened_uvar
-    pt.y[i] -= pt.s_u[j] / Aij
-    pt.s_u[j] = zero(T)
+    sol.y[i] -= sol.s_u[j] / Aij
+    sol.s_u[j] = zero(T)
   end
   if !operation.tightened_lvar && !operation.tightened_uvar
-    pt.y[i] = zero(T)
+    sol.y[i] = zero(T)
   end
 end

src/presolve/unconstrained_reductions.jl

@@ -41,4 +41,4 @@ function unconstrained_reductions!(
   return unbounded
 end
 
-postsolve!(pt::OutputPoint, operation::UnconstrainedReduction) = nothing
+postsolve!(sol::QMSolution, operation::UnconstrainedReduction) = nothing

test/test_presolve.jl

@@ -48,12 +48,9 @@
   @test psqp.meta.ifix == Int[]
   @test psqp.meta.nvar == 2
 
-  x_in = [4.0; 7.0]
-  y_in = [2.0; 2.0]
-  s_l = sparse([3.0; 2.0])
-  s_u = sparse([0.0; 0.0])
-  pt_out = postsolve(qp, psqp, x_in, y_in, s_l, s_u)
-  @test pt_out.x == [4.0; 2.0; 7.0]
+  sol_in = QMSolution([4.0; 7.0], [2.0; 2.0], sparse([3.0; 2.0]), sparse([0.0; 0.0]))
+  sol = postsolve(qp, psqp, sol_in)
+  @test sol.x == [4.0; 2.0; 7.0]
 
   # test that solves the problem
   qp2 = QuadraticModel(
@@ -113,12 +110,9 @@ end
   @test psqp.meta.lcon == psqp.meta.ucon == bps_true
   @test psqp.meta.ncon == 3
 
-  x_in = [4.0; 7.0; 4.0]
-  y_in = [2.0; 2.0; 4.0]
-  s_l = sparse([3.0; 4.0; 2.0])
-  s_u = sparse([0.0; 3.0; 0.0])
-  pt_out = postsolve(qp, psqp, x_in, y_in, s_l, s_u)
-  @test pt_out.y == [2.0; 0.0; 0.0; 2.0; 0.0; 4.0]
+  sol_in = QMSolution([4.0; 7.0; 4.0], [2.0; 2.0; 4.0], sparse([3.0; 4.0; 2.0]), sparse([0.0; 3.0; 0.0]))
+  sol = postsolve(qp, psqp, sol_in)
+  @test sol.y == [2.0; 0.0; 0.0; 2.0; 0.0; 4.0]
 end
 
 @testset "presolve singleton rows" begin
@@ -165,12 +159,9 @@ end
   @test Aps == sparse(Aps_true)
   @test psqp.meta.ncon == 2
 
-  x_in = [4.0; 7.0; 4.0]
-  y_in = [2.0; 4.0]
-  s_l = sparse([3.0; 4.0; 2.0])
-  s_u = sparse([0.0; 3.0; 0.0])
-  pt_out = postsolve(qp, psqp, x_in, y_in, s_l, s_u)
-  @test pt_out.y == [2.0; 0.0; 0.0; 1.0; 0.0; 4.0]
+  sol_in = QMSolution([4.0; 7.0; 4.0], [2.0; 4.0], sparse([3.0; 4.0; 2.0]), sparse([0.0; 3.0; 0.0]))
+  sol = postsolve(qp, psqp, sol_in)
+  @test sol.y == [2.0; 0.0; 0.0; 1.0; 0.0; 4.0]
 end
 
 @testset "presolve singleton rows and ifix" begin
@@ -239,10 +230,7 @@ end
   statsps = presolve(qp)
   psqp = statsps.solver_specific[:presolvedQM]
 
-  x_in = [4.0; 4.0]
-  y_in = [2.0; 4.0]
-  s_l = sparse([3.0; 2.0])
-  s_u = sparse([0.0; 0.0])
-  pt_out = postsolve(qp, psqp, x_in, y_in, s_l, s_u)
-  @test pt_out.x ≈ [4.0, -2.5, 4.0]
+  sol_in = QMSolution([4.0; 4.0], [2.0; 4.0], sparse([3.0; 2.0]), sparse([0.0; 0.0]))
+  sol = postsolve(qp, psqp, sol_in)
+  @test sol.x ≈ [4.0, -2.5, 4.0]
 end