add docs, fix issues ifix

Author: geoffroyleconte

src/QuadraticModels.jl

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

src/presolve/empty_rows.jl

@@ -5,10 +5,19 @@ function row_cnt!(Arows, row_cnt::Vector{Int})
   end
 end
 
-removed_empty_rows(row_cnt::Vector{Int}) = findall(isequal(0), row_cnt)
+find_empty_rows(row_cnt::Vector{Int}) = findall(isequal(0), row_cnt)
 
+"""
+    new_ncon = empty_rows!(Arows, lcon, ucon, ncon, row_cnt, empty_rows, Arows_s)
+
+Removes the empty rows of A, and the corresponding elements in lcon and ucon that are in `empty_rows`.
+`Arows_s` is a view of `Arows` sorted in ascending order.
+`row_cnt` is a vector of the number of elements per row.
+
+Returns the new number of constraints `new_ncon` and updates in-place `Arows`, `lcon`, `ucon`.
+"""
 function empty_rows!(Arows, lcon::Vector{T}, ucon::Vector{T}, ncon, row_cnt::Vector{Int}, 
-                     rows_rm::Vector{Int}, Arows_s) where {T}
+                     empty_rows::Vector{Int}, Arows_s) where {T}
   new_ncon = 0
   for i=1:ncon
     if row_cnt[i] == 0
@@ -23,24 +32,12 @@ function empty_rows!(Arows, lcon::Vector{T}, ucon::Vector{T}, ncon, row_cnt::Vec
   resize!(ucon, new_ncon)
 
   c_rm = 1
-  nrm = length(rows_rm)
+  nrm = length(empty_rows)
   for k=1:length(Arows)
-    while c_rm ≤ nrm && Arows_s[k] ≥ rows_rm[c_rm]
+    while c_rm ≤ nrm && Arows_s[k] ≥ empty_rows[c_rm]
       c_rm += 1 
     end
     Arows_s[k] -= c_rm - 1
   end
   return new_ncon
 end
-
-function restore_y!(y::Vector{T}, yout::Vector{T}, row_cnt, ncon) where {T}
-  c_y = 0
-  for i = 1:ncon
-    if row_cnt[i] == 0 || row_cnt[i] == 1
-      yout[i] = zero(T)
-    else
-      c_y += 1
-      yout[i] = y[c_y]
-    end
-  end
-end

src/presolve/postsolve_utils.jl

@@ -1,3 +1,30 @@
+function restore_ifix!(ifix, xrm, x, xout)
+  # put x and xrm inside xout
+  cfix, cx = 1, 1
+  nfix = length(ifix)
+  for i = 1:length(xout)
+    if cfix <= nfix && i == ifix[cfix]
+      xout[i] = xrm[cfix]
+      cfix += 1
+    else
+      xout[i] = x[cx]
+      cx += 1
+    end
+  end
+end
+
+function restore_y!(y::Vector{T}, yout::Vector{T}, kept_rows::Vector{Bool}, ncon) where {T}
+  c_y = 0
+  for i = 1:ncon
+    if !kept_rows[i]
+      yout[i] = zero(T)
+    else
+      c_y += 1
+      yout[i] = y[c_y]
+    end
+  end
+end
+
 function restore_ilow_iupp!(ilow, iupp, ifix)
   c_fix = 1
   nfix = length(ifix)

src/presolve/presolve.jl

@@ -6,7 +6,7 @@ include("postsolve_utils.jl")
 mutable struct PresolvedData{T, S}
   ifix::Vector{Int}
   xrm::S
-  row_cnt::Vector{Int}
+  kept_rows::Vector{Bool}
   nconps::Int
 end
 
@@ -17,6 +17,23 @@ mutable struct PresolvedQuadraticModel{T, S, M1, M2} <: AbstractQuadraticModel{T
   psd::PresolvedData{T, S}
 end
 
+function update_kept_rows!(kept_rows::Vector{Bool}, vec_to_rm::Vector{Int})
+  ncon = length(kept_rows)
+  n_rm = length(vec_to_rm)
+  offset = 0
+  c_v = 1
+  for i=1:ncon
+    if !kept_rows[i]
+      offset += 1
+    else
+      if c_v ≤ n_rm && vec_to_rm[c_v] + offset == i
+        kept_rows[i] = false
+        c_v += 1
+      end
+    end
+  end
+end
+
 """
     stats_ps = presolve(qm::QuadraticModel{T, S}; kwargs...)
 
@@ -25,13 +42,15 @@ Apply a presolve routine to `qm` and returns a
 from the package [`SolverCore.jl`](https://github.com/JuliaSmoothOptimizers/SolverCore.jl).
 The presolve operations currently implemented are:
 
+- [`empty_rows!`](@ref) : remove empty rows
+- [`singleton_rows!`](@ref) : remove singleton rows
 - [`remove_ifix!`](@ref) : remove fixed variables
 
 The `PresolvedQuadraticModel{T, S} <: AbstractQuadraticModel{T, S}` is located in the `solver_specific` field:
 
     psqm = stats_ps.solver_specific[:presolvedQM]
 
-and should be used to call [`postsolve!`](@ref).
+and should be used to call [`postsolve`](@ref).
 
 If the presolved problem has 0 variables, `stats_ps.solution` contains a solution of the primal problem,
 `stats_ps.multipliers` is a zero `SparseVector`, and, if we define
@@ -50,31 +69,39 @@ function presolve(
   lvar, uvar = psqm.meta.lvar, psqm.meta.uvar
   lcon, ucon = psqm.meta.lcon, psqm.meta.ucon
   nvar, ncon = psqm.meta.nvar, psqm.meta.ncon
+  # copy if same vector
+  lcon === ucon && (lcon = copy(lcon))
+  lvar === uvar && (lvar = copy(lvar))
 
   # empty rows
   row_cnt = zeros(Int, ncon)
+  kept_rows = fill(true, ncon)
   row_cnt!(psdata.A.rows, row_cnt) # number of coefficients per row
-  rows_rm = removed_empty_rows(row_cnt) # indices of the empty rows
-  if length(rows_rm) > 0
+  empty_rows = find_empty_rows(row_cnt) # indices of the empty rows
+  if length(empty_rows) > 0
+    empty_row_pass = true
+    update_kept_rows!(kept_rows, empty_rows)
     Arows_sortperm = sortperm(psdata.A.rows) # permute rows 
     Arows_s = @views psdata.A.rows[Arows_sortperm]
-    nconps = empty_rows!(psdata.A.rows, lcon, ucon, ncon, row_cnt, rows_rm, Arows_s)
+    nconps = empty_rows!(psdata.A.rows, lcon, ucon, ncon, row_cnt, empty_rows, Arows_s)
   else
+    empty_row_pass = false
     nconps = ncon
   end
 
   # remove singleton rows
-  if nconps != ncon
-    row_cnt2 = Vector{Int}(undef, nconps)
-  else
-    row_cnt2 = row_cnt
+  if empty_row_pass
+    resize!(row_cnt, nconps)
+    row_cnt .= 0
+    row_cnt!(psdata.A.rows, row_cnt) # number of coefficients per rows
   end
-  row_cnt2 .= 0
-  row_cnt!(psdata.A.rows, row_cnt2) # number of coefficients per rows
-  singl_rows = removed_singleton_rows(row_cnt2) # indices of the empty rows
+  singl_rows = find_singleton_rows(row_cnt) # indices of the singleton rows
   if length(singl_rows) > 0
-    nconps = singleton_rows!(psdata.A.rows, psdata.A.cols, psdata.A.vals, lcon, ucon, lvar, uvar, nvar, nconps, row_cnt2, singl_rows)
+    singl_row_pass = true
+    update_kept_rows!(kept_rows, singl_rows)
+    nconps = singleton_rows!(psdata.A.rows, psdata.A.cols, psdata.A.vals, lcon, ucon, lvar, uvar, nvar, nconps, row_cnt, singl_rows)
   else
+    singl_row_pass = false
     nconps = nconps
   end
 
@@ -146,7 +173,7 @@ function presolve(
       minimize = qm.meta.minimize,
       kwargs...,
     )
-    psd = PresolvedData{T, S}(ifix, xrm, row_cnt, nconps)
+    psd = PresolvedData{T, S}(ifix, xrm, kept_rows, nconps)
     ps = PresolvedQuadraticModel(psmeta, Counters(), psdata, psd)
     return GenericExecutionStats(
       :unknown,
@@ -158,13 +185,6 @@ function presolve(
   end
 end
 
-"""
-    postsolve!(qm::QuadraticModel{T, S}, psqm::PresolvedQuadraticModel{T, S}, 
-               x_in::S, x_out::S) where {T, S}
-
-Retrieve the solution `x_out` of the original QP `qm` given the solution of the presolved QP (`psqm`)
-`x_in`.
-"""
 function postsolve!(
   qm::QuadraticModel{T, S},
   psqm::PresolvedQuadraticModel{T, S},
@@ -182,11 +202,34 @@ function postsolve!(
     x_out .= @views x_in[1:(qm.meta.nvar)]
   end
   ncon = length(y_out)
-  restore_y!(y_in, y_out, psqm.psd.row_cnt, ncon)
+  restore_y!(y_in, y_out, psqm.psd.kept_rows, ncon)
 
   ilow, iupp = s_l.nzind, s_u.nzind
   restore_ilow_iupp!(ilow, iupp, ifix)
   s_l_out = SparseVector(qm.meta.nvar, ilow, s_l.nzval)
   s_u_out = SparseVector(qm.meta.nvar, iupp, s_u.nzval)
   return s_l_out, s_u_out
 end
+
+"""
+    x, y, s_l, s_u = 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}
+
+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`.
+"""
+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},
+) where {T, S}
+  x_out = similar(qm.meta.x0)
+  y_out = similar(qm.meta.y0)
+  s_l_out, s_u_out = postsolve!(qm, psqm, x_in, x_out, y_in, y_out, s_l, s_u)
+  return x_out, y_out, s_l_out, s_u_out
+end

src/presolve/remove_ifix.jl

@@ -95,8 +95,9 @@ function remove_ifix!(
       Ai, Aj, Ax = Arows[k], Acols[k], Avals[k]
       if Aj == newcurrentifix
         Arm_tmp += 1
-        lcon[Ai] -= Ax * xifix
-        ucon[Ai] -= Ax * xifix
+        con_offset = Ax * xifix
+        lcon[Ai] -= con_offset
+        ucon[Ai] -= con_offset
       else
         Arows[Awritepos] = Ai
         Acols[Awritepos] = (Aj < newcurrentifix) ? Aj : Aj - 1
@@ -138,18 +139,3 @@ function remove_ifix!(
 
   return xrm, c0ps, nvarrm
 end
-
-function restore_ifix!(ifix, xrm, x, xout)
-  # put x and xrm inside xout
-  cfix, cx = 1, 1
-  nfix = length(ifix)
-  for i = 1:length(xout)
-    if cfix <= nfix && i == ifix[cfix]
-      xout[i] = xrm[cfix]
-      cfix += 1
-    else
-      xout[i] = x[cx]
-      cx += 1
-    end
-  end
-end

src/presolve/singleton_rows.jl

@@ -1,5 +1,14 @@
-removed_singleton_rows(row_cnt::Vector{Int}) = findall(isequal(1), row_cnt)
+find_singleton_rows(row_cnt::Vector{Int}) = findall(isequal(1), row_cnt)
 
+"""
+    new_ncon = singleton_rows!(Arows, Acols, Avals, lcon, ucon,
+                               lvar, uvar, nvar, ncon, row_cnt, singl_rows)
+
+Removes the singleton rows of A, and the corresponding elements in lcon and ucon that are in `singl_rows`.
+`row_cnt` is a vector of the number of elements per row.
+
+Returns the new number of constraints `new_ncon` and updates in-place `Arows`, `Acols`, `Avals`, `lcon`, `ucon`, `lvar`, `uvar`.
+"""
 function singleton_rows!(Arows, Acols, Avals, lcon::Vector{T}, ucon::Vector{T}, lvar, uvar,
                          nvar, ncon, row_cnt::Vector{Int}, 
                          singl_rows::Vector{Int}) where {T}
@@ -21,7 +30,6 @@ function singleton_rows!(Arows, Acols, Avals, lcon::Vector{T}, ucon::Vector{T},
     while k <= Annz - idxsingl + 1
       Ai, Aj, Ax = Arows[k], Acols[k], Avals[k]
       if Ai == newcurrentisingl
-        # oldAi = Ai + idxsingl - 1
         if Ax > zero(T)
           lvar[Aj] = max(lvar[Aj], lcon[currentisingl] / Ax)
           uvar[Aj] = min(uvar[Aj], ucon[currentisingl] / Ax)
@@ -48,15 +56,7 @@ function singleton_rows!(Arows, Acols, Avals, lcon::Vector{T}, ucon::Vector{T},
     resize!(Avals, Annz)
   end
 
-  new_ncon = 0
-  for i=1:ncon
-    if row_cnt[i] != 1
-      new_ncon += 1
-      lcon[new_ncon] = lcon[i]
-      ucon[new_ncon] = ucon[i]
-    end
-  end
-  resize!(lcon, new_ncon)
-  resize!(ucon, new_ncon)
-  return new_ncon
+  deleteat!(lcon, singl_rows)
+  deleteat!(ucon, singl_rows)
+  return ncon - nsingl
 end

test/test_presolve.jl

@@ -53,10 +53,10 @@
   @test psqp.meta.nvar == 2
 
   x_in = [4.0; 7.0]
-  x_out = zeros(3)
   y_in = [2.0; 2.0]
-  y_out = zeros(2)
-  postsolve!(qp, psqp, x_in, x_out, y_in, y_out)
+  s_l = sparse([3.0; 2.0])
+  s_u = sparse([0.0; 0.0])
+  x_out, y_out, s_l, s_u = postsolve(qp, psqp, x_in, y_in, s_l, s_u)
   @test x_out == [4.0; 2.0; 7.0]
 
   # test that solves the problem
@@ -118,10 +118,10 @@ end
   @test psqp.meta.ncon == 3
 
   x_in = [4.0; 7.0; 4.0]
-  x_out = zeros(3)
   y_in = [2.0; 2.0; 4.0]
-  y_out = zeros(6)
-  postsolve!(qp, psqp, x_in, x_out, y_in, y_out)
+  s_l = sparse([3.0; 4.0; 2.0])
+  s_u = sparse([0.0; 3.0; 0.0])
+  x_out, y_out, s_l, s_u = postsolve(qp, psqp, x_in, y_in, s_l, s_u)
   @test y_out == [2.0; 0.0; 0.0; 2.0; 0.0; 4.0]
 end
 
@@ -170,10 +170,10 @@ end
   @test psqp.meta.ncon == 2
 
   x_in = [4.0; 7.0; 4.0]
-  x_out = zeros(3)
   y_in = [2.0; 4.0]
-  y_out = zeros(6)
-  postsolve!(qp, psqp, x_in, x_out, y_in, y_out)
+  s_l = sparse([3.0; 4.0; 2.0])
+  s_u = sparse([0.0; 3.0; 0.0])
+  x_out, y_out, s_l, s_u = postsolve(qp, psqp, x_in, y_in, s_l, s_u)
   @test y_out == [2.0; 0.0; 0.0; 0.0; 0.0; 4.0]
 end
 
@@ -221,10 +221,10 @@ end
   @test psqp.meta.ncon == 2
 
   x_in = [7.0; 4.0]
-  x_out = zeros(3)
   y_in = [2.0; 4.0]
-  y_out = zeros(6)
-  postsolve!(qp, psqp, x_in, x_out, y_in, y_out)
+  s_l = sparse([4.0; 2.0])
+  s_u = sparse([3.0; 0.0])
+  x_out, y_out, s_l, s_u = postsolve(qp, psqp, x_in, y_in, s_l, s_u)
   @test y_out == [2.0; 0.0; 0.0; 0.0; 0.0; 4.0]
   @test x_out == [4.0/3.0; 7.0; 4.0]
 end