Merge pull request #1 from joehuchette/port-from-tulip

Port over code from Tulip

Project.toml

@@ -3,5 +3,8 @@ uuid = "c46b0367-4903-4ff3-bf29-73e6bd648aaa"
 authors = ["Mathieu Tanneau and contributors"]
 version = "0.1.0"
 
+[deps]
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+
 [compat]
 julia = "1.5"

src/MathOptPresolve.jl

@@ -1,5 +1,10 @@
 module MathOptPresolve
 
-# Write your package code here.
+using SparseArrays
+
+include("status.jl")
+include("problem_data.jl")
+include("solution.jl")
+include("presolve.jl")
 
 end

src/lp/dominated_column.jl

@@ -0,0 +1,157 @@
+struct DominatedColumn{T} <: PresolveTransformation{T}
+    j::Int
+    x::T  # Primal value
+    cj::T  # Objective
+    col::Col{T}  # Column
+end
+
+function remove_dominated_column!(ps::PresolveData{T}, j::Int; tol::T=100*sqrt(eps(T))) where{T}
+    ps.colflag[j] || return nothing
+
+    # Compute implied bounds on reduced cost: `ls ≤ s ≤ us`
+    ls = us = zero(T)
+    col = ps.pb0.acols[j]
+    for (i, aij) in zip(col.nzind, col.nzval)
+        (ps.rowflag[i] && !iszero(aij)) || continue
+
+        ls += aij * ( (aij >= zero(T)) ? ps.ly[i] : ps.uy[i] )
+        us += aij * ( (aij >= zero(T)) ? ps.uy[i] : ps.ly[i] )
+    end
+
+    # Check if column is dominated
+    cj = ps.obj[j]
+    if cj - us > tol
+        # Reduced cost is always positive => fix to lower bound (or problem is unbounded)
+        lb = ps.lcol[j]
+        @debug "Fixing dominated column $j to its lower bound $lb"
+
+        if !isfinite(lb)
+            # Problem is dual infeasible
+            @debug "Column $j is (lower) unbounded"
+            ps.status = Trm_DualInfeasible
+            ps.updated = true
+
+            # Resize problem
+            compute_index_mapping!(ps)
+            resize!(ps.solution, ps.nrow, ps.ncol)
+            ps.solution.x .= zero(T)
+            ps.solution.y_lower .= zero(T)
+            ps.solution.y_upper .= zero(T)
+            ps.solution.s_lower .= zero(T)
+            ps.solution.s_upper .= zero(T)
+
+            # Unbounded ray: xj = -1
+            ps.solution.primal_status = Sln_InfeasibilityCertificate
+            ps.solution.dual_status = Sln_Unknown
+            ps.solution.is_primal_ray = true
+            ps.solution.is_dual_ray = false
+            ps.solution.z_primal = ps.solution.z_dual = -T(Inf)
+            j_ = ps.new_var_idx[j]
+            ps.solution.x[j_] = -one(T)
+
+            return nothing
+        end
+
+        # Update objective
+        ps.obj0 += cj * lb
+
+        # Extract column and update rows
+        col_ = Col{T}(Int[], T[])
+        for (i, aij) in zip(col.nzind, col.nzval)
+            ps.rowflag[i] || continue
+
+            push!(col_.nzind, i)
+            push!(col_.nzval, aij)
+
+            # Update bounds and non-zeros
+            ps.lrow[i] -= aij * lb
+            ps.urow[i] -= aij * lb
+            ps.nzrow[i] -= 1
+
+            ps.nzrow[i] == 1 && push!(ps.row_singletons, i)
+        end
+
+        # Remove variable from problem
+        push!(ps.ops, DominatedColumn(j, lb, cj, col_))
+        ps.colflag[j] = false
+        ps.ncol -= 1
+        ps.updated = true
+
+    elseif cj - ls < -tol
+        # Reduced cost is always negative => fix to upper bound (or problem is unbounded)
+        ub = ps.ucol[j]
+
+        if !isfinite(ub)
+            # Problem is unbounded
+            @debug "Column $j is (upper) unbounded"
+
+            ps.status = Trm_DualInfeasible
+            ps.updated = true
+
+            # Resize solution
+            compute_index_mapping!(ps)
+            resize!(ps.solution, ps.nrow, ps.ncol)
+            ps.solution.x .= zero(T)
+            ps.solution.y_lower .= zero(T)
+            ps.solution.y_upper .= zero(T)
+            ps.solution.s_lower .= zero(T)
+            ps.solution.s_upper .= zero(T)
+
+            # Unbounded ray: xj = -1
+            ps.solution.primal_status = Sln_InfeasibilityCertificate
+            ps.solution.dual_status = Sln_Unknown
+            ps.solution.is_primal_ray = true
+            ps.solution.is_dual_ray = false
+            ps.solution.z_primal = ps.solution.z_dual = -T(Inf)
+            j_ = ps.new_var_idx[j]
+            ps.solution.x[j_] = one(T)
+
+            return nothing
+        end
+
+        @debug "Fixing dominated column $j to its upper bound $ub"
+
+        # Update objective
+        ps.obj0 += cj * ub
+
+        # Extract column and update rows
+        col_ = Col{T}(Int[], T[])
+        for (i, aij) in zip(col.nzind, col.nzval)
+            ps.rowflag[i] || continue
+
+            push!(col_.nzind, i)
+            push!(col_.nzval, aij)
+
+            # Update bounds and non-zeros
+            ps.lrow[i] -= aij * ub
+            ps.urow[i] -= aij * ub
+            ps.nzrow[i] -= 1
+
+            ps.nzrow[i] == 1 && push!(ps.row_singletons, i)
+        end
+
+        # Remove variable from problem
+        push!(ps.ops, DominatedColumn(j, ub, cj, col_))
+        ps.colflag[j] = false
+        ps.ncol -= 1
+        ps.updated = true
+    end
+
+    return nothing
+end
+
+function postsolve!(sol::Solution{T}, op::DominatedColumn{T}) where{T}
+    # Primal value
+    sol.x[op.j] = op.x
+
+    # Reduced cost
+    s = sol.is_dual_ray ? zero(T) : op.cj
+    for (i, aij) in zip(op.col.nzind, op.col.nzval)
+        s -= aij * (sol.y_lower[i] - sol.y_upper[i])
+    end
+
+    sol.s_lower[op.j] = pos_part(s)
+    sol.s_upper[op.j] = neg_part(s)
+
+    return nothing
+end

src/lp/empty_column.jl

@@ -0,0 +1,105 @@
+struct EmptyColumn{T} <: PresolveTransformation{T}
+    j::Int  # variable index
+    x::T  # Variable value
+    s::T  # Reduced cost
+end
+
+function remove_empty_column!(ps::PresolveData{T}, j::Int) where{T}
+    # Sanity check
+    (ps.colflag[j] && (ps.nzcol[j] == 0)) || return nothing
+
+    # Remove column
+    lb, ub = ps.lcol[j], ps.ucol[j]
+    cj = ps.obj[j]
+    @debug "Removing empty column $j" cj lb ub
+
+    if cj > zero(T)
+        if isfinite(lb)
+            # Set variable to lower bound
+            # Update objective constant
+            ps.obj0 += lb * cj
+            push!(ps.ops, EmptyColumn(j, lb, cj))
+        else
+            # Problem is dual infeasible
+            @debug "Column $j is (lower) unbounded"
+            ps.status = Trm_DualInfeasible
+            ps.updated = true
+
+            # Resize problem
+            compute_index_mapping!(ps)
+            resize!(ps.solution, ps.nrow, ps.ncol)
+            ps.solution.x .= zero(T)
+            ps.solution.y_lower .= zero(T)
+            ps.solution.y_upper .= zero(T)
+            ps.solution.s_lower .= zero(T)
+            ps.solution.s_upper .= zero(T)
+
+            # Unbounded ray: xj = -1
+            ps.solution.primal_status = Sln_InfeasibilityCertificate
+            ps.solution.dual_status = Sln_Unknown
+            ps.solution.is_primal_ray = true
+            ps.solution.is_dual_ray = false
+            ps.solution.z_primal = ps.solution.z_dual = -T(Inf)
+            j_ = ps.new_var_idx[j]
+            ps.solution.x[j_] = -one(T)
+
+            return nothing
+        end
+    elseif cj < zero(T)
+        if isfinite(ub)
+            # Set variable to upper bound
+            # Update objective constant
+            ps.obj0 += ub * cj
+            push!(ps.ops, EmptyColumn(j, ub, cj))
+        else
+            # Problem is dual infeasible
+            @debug "Column $j is (upper) unbounded"
+            ps.status = Trm_DualInfeasible
+            ps.updated = true
+
+            # Resize problem
+            compute_index_mapping!(ps)
+            resize!(ps.solution, ps.nrow, ps.ncol)
+            ps.solution.x .= zero(T)
+            ps.solution.y_lower .= zero(T)
+            ps.solution.y_upper .= zero(T)
+            ps.solution.s_lower .= zero(T)
+            ps.solution.s_upper .= zero(T)
+
+            # Unbounded ray: xj = 1
+            ps.solution.primal_status = Sln_InfeasibilityCertificate
+            ps.solution.dual_status = Sln_Unknown
+            ps.solution.is_primal_ray = true
+            ps.solution.is_dual_ray = false
+            ps.solution.z_primal = ps.solution.z_dual = -T(Inf)
+            j_ = ps.new_var_idx[j]
+            ps.solution.x[j_] = one(T)
+
+            return
+        end
+    else
+        # Any feasible value works
+        if isfinite(lb)
+            push!(ps.ops, EmptyColumn(j, lb, zero(T)))
+        elseif isfinite(ub)
+            push!(ps.ops, EmptyColumn(j, ub, zero(T)))
+        else
+            # Free variable with zero coefficient and empty column
+            push!(ps.ops, EmptyColumn(j, zero(T), zero(T)))
+        end
+
+    end
+
+    # Book=keeping
+    ps.colflag[j] = false
+    ps.updated = true
+    ps.ncol -= 1
+    return nothing
+end
+
+function postsolve!(sol::Solution{T}, op::EmptyColumn{T}) where{T}
+    sol.x[op.j] = op.x
+    sol.s_lower[op.j] = pos_part(op.s)
+    sol.s_upper[op.j] = neg_part(op.s)
+    return nothing
+end

src/lp/empty_row.jl

@@ -0,0 +1,77 @@
+# TODO: this is redundant with forcing constraint checks
+#   => an empty row is automatically redundant or infeasible
+
+struct EmptyRow{T} <: PresolveTransformation{T}
+    i::Int  # row index
+    y::T  # dual multiplier
+end
+
+function remove_empty_row!(ps::PresolveData{T}, i::Int) where{T}
+    # Sanity checks
+    (ps.rowflag[i] && ps.nzrow[i] == 0) || return nothing
+
+    # Check bounds
+    lb, ub = ps.lrow[i], ps.urow[i]
+
+    if ub < zero(T)
+        # Infeasible
+        @debug "Row $i is primal infeasible"
+        ps.status = Trm_PrimalInfeasible
+        ps.updated = true
+
+        # Resize problem
+        compute_index_mapping!(ps)
+        resize!(ps.solution, ps.nrow, ps.ncol)
+        ps.solution.x .= zero(T)
+        ps.solution.y_lower .= zero(T)
+        ps.solution.y_upper .= zero(T)
+        ps.solution.s_lower .= zero(T)
+        ps.solution.s_upper .= zero(T)
+
+        # Farkas ray: y⁺_i = 1 (any > 0 value works)
+        ps.solution.primal_status = Sln_Unknown
+        ps.solution.dual_status = Sln_InfeasibilityCertificate
+        ps.solution.is_primal_ray = false
+        ps.solution.is_dual_ray = true
+        ps.solution.z_primal = ps.solution.z_dual = T(Inf)
+        i_ = ps.new_con_idx[i]
+        ps.solution.y_upper[i] = one(T)
+        return
+    elseif lb > zero(T)
+        @debug "Row $i is primal infeasible"
+        ps.status = Trm_PrimalInfeasible
+        ps.updated = true
+
+        # Resize problem
+        compute_index_mapping!(ps)
+        resize!(ps.solution, ps.nrow, ps.ncol)
+        ps.solution.x .= zero(T)
+        ps.solution.y_lower .= zero(T)
+        ps.solution.y_upper .= zero(T)
+        ps.solution.s_lower .= zero(T)
+        ps.solution.s_upper .= zero(T)
+
+        # Farkas ray: y⁺_i = 1 (any > 0 value works)
+        ps.solution.primal_status = Sln_Unknown
+        ps.solution.dual_status = Sln_InfeasibilityCertificate
+        ps.solution.is_primal_ray = false
+        ps.solution.is_dual_ray = true
+        ps.solution.z_primal = ps.solution.z_dual = T(Inf)
+        i_ = ps.new_con_idx[i]
+        ps.solution.y_lower[i_] = one(T)
+        return
+    else
+        push!(ps.ops, EmptyRow(i, zero(T)))
+    end
+
+    # Book-keeping
+    ps.updated = true
+    ps.rowflag[i] = false
+    ps.nrow -= 1
+end
+
+function postsolve!(sol::Solution{T}, op::EmptyRow{T}) where{T}
+    sol.y_lower[op.i] = pos_part(op.y)
+    sol.y_upper[op.i] = neg_part(op.y)
+    return nothing
+end