Implement Polyhedra.isredundant

Author: blegat

.github/workflows/ci.yml

@@ -12,6 +12,9 @@ jobs:
       fail-fast: false
       matrix:
         include:
+          - version: '1.10'
+            os: ubuntu-latest
+            arch: x64
           - version: '1'
             os: ubuntu-latest
             arch: x64

Project.toml

@@ -1,7 +1,7 @@
 name = "LRSLib"
 uuid = "262c1cb6-76e2-5873-868b-19ece3183cc5"
 repo = "https://github.com/JuliaPolyhedra/LRSLib.jl.git"
-version = "0.8.2"
+version = "0.9.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -11,7 +11,7 @@ lrslib_jll = "3873f7d0-7b7c-52c3-bdf4-8ab39b8f337a"
 
 [compat]
 Polyhedra = "0.7, 0.8"
-julia = "1.6"
+julia = "1.10"
 lrslib_jll = "= 0.3.3"
 
 [extras]

src/lp.jl

@@ -172,8 +172,7 @@ function print_lpoutput(m::HMatrix, print_solution::Bool=false)
     for i in 0:P.d-1
         C_i = unsafe_load(P.C, i+1)
         Col_i = unsafe_load(P.Col, i+1)
-        idx = _unsafe_load_inequality(m, C_i)
-        print(" y_$(idx)=")
+        print(" y_$(_unsafe_load_inequality(m, C_i).value)=")
         temp1 = extractbigint(unsafe_load(Q.Lcm, Col_i+1)) * (-1)
         temp1 *= extractbigint(unsafe_load(unsafe_load(P.A, 1), Col_i+1))
         temp2 = extractbigint(unsafe_load(Q.Gcd, Col_i+1))

src/matrix.jl

@@ -177,7 +177,20 @@ Base.get(rep::RepMatrix, idx::Polyhedra.Index{Rational{BigInt}}) = Polyhedra.val
 
 function _unsafe_load_inequality(m::RepMatrix, idx)
     Q = unsafe_load(m.Q)
-    return unsafe_load(Q.inequality, idx - Q.lastdv + 1)
+    _to_index(m, unsafe_load(Q.inequality, idx - Q.lastdv + 1))
+end
+
+# For H-rep, LRS internal index `0` corresponds to objective and `1` correspond to the first halfspace/hyperplane
+# so it's the same 1-based indexing than `Polyhedra.HIndex`
+_to_index(::HMatrix{T}, i) = Polyhedra.HIndex{T,HalfSpace{T,Vector{T}}}(i)
+
+# For V-rep, LRS internal index `0` first halfspace/hyperplane
+# so it's the 0-based indexing instead of `Polyhedra.HIndex` so we need to do `+1`
+_to_index(::VMatrix{T}, i) = Polyhedra.VIndex{T,Vector{T}}(i + 1)
+
+function _unsafe_load_inequality(m::VMatrix, Polyhedra.VIndex)
+    Q = unsafe_load(m.Q)
+    return unsafe_load(Q.inequality, idx - Q.lastdv)
 end
 
 # H-representation

src/nash.jl

@@ -125,8 +125,8 @@ function nash2_main(hr1::HMatrix, hr2::HMatrix, linindex::Vector{Clong})
                 unsafe_load(unsafe_load(P1).B, i+1),
             )
             if ((unsafe_load(Q1).nlinearity == 0) ||
-                (j < unsafe_load(unsafe_load(Q1).linearity, 1)))
-                unsafe_store!(linearity, j, nlinearity+1)
+                (j.value < unsafe_load(unsafe_load(Q1).linearity, 1)))
+                unsafe_store!(linearity, j.value, nlinearity+1)
                 nlinearity += 1
             end
         end

src/redund.jl

@@ -1,3 +1,15 @@
+function _inequality_indices(m::HMatrix)
+    # FIXME Not sure what `lastdv` is doing, isn't it just always zero ?
+    lastdv = unsafe_load(m.Q).lastdv
+    return (lastdv + 1):(m_A + d)
+end
+
+function _inequality_indices(m::VMatrix)
+    # FIXME Not sure what `lastdv` is doing, isn't it just always zero ?
+    lastdv = unsafe_load(m.Q).lastdv
+    return (lastdv + 2):(m_A + d + 1)
+end
+
 function redund(m::RepMatrix)
     # if non-negative flag is set, non-negative constraints are not input
     # explicitly, and are not checked for redundancy
@@ -25,19 +37,19 @@ function redund(m::RepMatrix)
     # rows 0..lastdv are cost, decision variables, or linearities
     # other rows need to be tested
 
-    redset = BitSet()
+    redset = [] # FIXME concrete type
     for index in (lastdv + 1):(m_A + d)
-        ineq = _unsafe_load_inequality(m, index) # the input inequality number corr. to this index
 
         status = checkindex(m, index)
         if :redundant == checkindex(m, index)
-            push!(redset, ineq)
+            # the input inequality number corr. to this index
+            push!(redset, _unsafe_load_inequality(m, index))
         end
     end
     redset
 end
 
-function redundi(m::RepMatrix, ineq::Int)
+function Polyhedra.isredundant(m::RepMatrix, ineq::Polyhedra.Index)
     if m.status == :AtNoBasis
         getfirstbasis(m)
     end