exanauts · kibaekkim · Oct 7, 2020 · Oct 8, 2020 · Oct 8, 2020 · Oct 8, 2020
diff --git a/netlib/getall.sh b/netlib/getall.sh
@@ -0,0 +1,112 @@
+#!/bin/bash
+
+instances=(
+    "25fv47"
+    "80bau3b"
+    "adlittle"
+    "afiro"
+    "agg"
+    "agg2"
+    "agg3"
+    "ascii"
+    "bandm"
+    "beaconfd"
+    "blend"
+    "bnl1"
+    "bnl2"
+    "boeing1"
+    "boeing2"
+    "bore3d"
+    "brandy"
+    "capri"
+    "changes"
+    "cycle"
+    "czprob"
+    "d2q06c"
+    "d6cube"
+    "degen2"
+    "degen3"
+    "dfl001"
+    "e226"
+    "etamacro"
+    "fffff800"
+    "finnis"
+    "fit1d"
+    "fit1p"
+    "fit2d"
+    "fit2p"
+    "forplan"
+    "ganges"
+    "gfrd-pnc"
+    "greenbea"
+    "greenbeb"
+    "grow15"
+    "grow22"
+    "grow7"
+    "israel"
+    "kb2"
+    "kennington"
+    "lotfi"
+    "maros"
+    "maros-r7"
+    "minos"
+    "modszk1"
+    "nesm"
+    "perold"
+    "pilot"
+    "pilot.ja"
+    "pilot.we"
+    "pilot4"
+    "pilot87"
+    "pilotnov"
+    "recipe"
+    "sc105"
+    "sc205"
+    "sc50a"
+    "sc50b"
+    "scagr25"
+    "scagr7"
+    "scfxm1"
+    "scfxm2"
+    "scfxm3"
+    "scorpion"
+    "scrs8"
+    "scsd1"
+    "scsd6"
+    "scsd8"
+    "sctap1"
+    "sctap2"
+    "sctap3"
+    "seba"
+    "share1b"
+    "share2b"
+    "shell"
+    "ship04l"
+    "ship04s"
+    "ship08l"
+    "ship08s"
+    "ship12l"
+    "ship12s"
+    "sierra"
+    "stair"
+    "standata"
+    "standgub"
+    "standmps"
+    "stocfor1"
+    "stocfor2"
+    "stocfor3"
+    "stocfor3.old"
+    "truss"
+    "tuff"
+    "vtp.base"
+    "wood1p"
+    "woodw"
+)
+
+for i in "${instances[@]}"
+do
+    wget https://www.netlib.org/lp/data/$i
+    ./emps $i > $i.mps
+    rm $i
+done
+
diff --git a/src/Constant.jl b/src/Constant.jl
diff --git a/src/LP.jl b/src/LP.jl
@@ -43,7 +43,7 @@ function canonical_form(standard::MatOI.LPForm{T, AT, VT}) where {T, AT, VT}
     @assert length(xl) == lp_ncols
     @assert length(xu) == lp_ncols
     for i = 1:nineq
-        if standard.c_lb[ineq[i]] > -INF
+        if isFinite(standard.c_lb[ineq[i]])
             xu[ncols(standard) + i] = standard.c_ub[ineq[i]] - standard.c_lb[ineq[i]]
         end
     end
@@ -54,10 +54,10 @@ function canonical_form(standard::MatOI.LPForm{T, AT, VT}) where {T, AT, VT}
     # create the submatrix for slack
     I = Int64[]; J = Int64[]; V = Float64[]; j = 1;
     for i in ineq
-        if standard.c_lb[i] > -INF
+        if isFinite(standard.c_lb[i])
             push!(I, i); push!(J, j); push!(V, -1.0);
             j += 1
-        elseif standard.c_ub[i] < INF
+        elseif isFinite(standard.c_ub[i])
             b[i] = standard.c_ub[i]
             push!(I, i); push!(J, j); push!(V, 1.0);
             j += 1

diff --git a/src/Parameters.jl b/src/Parameters.jl
@@ -0,0 +1,10 @@
+"""
+    Parameters
+"""
+Base.@kwdef mutable struct Parameters
+    use_invB::Bool = true
+    max_iter::Int = 1000000000
+    max_history::Int = 100
+    print_iter_freq::Int = 10
+    pivot_rule::Type = DantzigRule
+end
diff --git a/src/PhaseOne/Artificial.jl b/src/PhaseOne/Artificial.jl
@@ -143,7 +143,7 @@ function run(prob::MatOI.LPSolverForm, Tv::Type; kwargs...)::Simplex.SpxData
 
     # load the problem
     spx = Simplex.SpxData(p1lp, Tv)
-    spx.pivot_rule = deepcopy(pivot_rule)
+    spx.params.pivot_rule = pivot_rule
 
     # set basis
     basic = collect((ncols(p1lp)-nrows(p1lp)+1):ncols(p1lp))

diff --git a/src/PhaseOne/Cplex.jl b/src/PhaseOne/Cplex.jl
@@ -75,7 +75,7 @@ function run(lp::MatOI.LPSolverForm, Tv::Type; kwargs...)::Simplex.SpxData
 
     # load the problem
     spx = Simplex.SpxData(cpxlp, Tv)
-    spx.pivot_rule = deepcopy(pivot_rule)
+    spx.params.pivot_rule = deepcopy(pivot_rule)
 
     # set basis
     Simplex.set_basis(spx, B)

diff --git a/src/PivotRules.jl b/src/PivotRules.jl
@@ -1,16 +1,17 @@
-# Bland's rule
-function pivot_Bland(spx::SpxData)
-    spx.enter = -1
+"""
+    BlandRule
 
-    # compute reduced cost
-    compute_reduced_cost(spx)
-    # @show spx.r
-    # @show spx.nonbasic
+Empty struct for Bland's pivot rule that can avoid cycling in pivoting
+"""
+struct BlandRule{T} <: AbstractPivotRule{T} end
 
+"Find an entering variable based on Bland's rule"
+function find_entering_variable!(::Type{BlandRule}, spx::SpxData)
+    spx.enter = -1
     min_enter = ncols(spx.lpdata) + 1
     for j = 1:length(spx.r)
         if spx.nonbasic[j] < min_enter
-            if spx.basis_status[spx.nonbasic[j]] == BASIS_AT_UPPER
+            if spx.basis_status[spx.nonbasic[j]] == Basis_At_Upper
                 if spx.r[j] > 1e-10
                     spx.enter = spx.nonbasic[j]
                     spx.enter_pos = j
@@ -28,77 +29,91 @@ function pivot_Bland(spx::SpxData)
     # @show spx.enter, spx.r[spx.enter_pos]
 end
 
-# Steepest-edge rule
 """
-Steepest-edge rule:
+    SteepestEdge
+
+Mutable struct for steepest-edge pivoting rule:
 
 This is based on the following references.
 1. Bieling et al. An efficient GPU implementation of the revised simplex method. 2010
 2. Forest and Goldfarb. Steepest-edge simplex algorithms for linear programming. 1992
 3. Goldfarb and Reid. A practicable steepest-edge simplex algorithm. 1977
 """
-function pivot_steepest_edge(spx::SpxData)
-
-    compute_reduced_cost(spx)
-    # @show spx.r
-    # @show spx.nonbasic
 
-    compute_steepest_edge_weights(spx)
-
-    # compute the local slope
-    # @show length(spx.r), length(spx.gamma)
-    spx.s .= (spx.r).^2 ./ spx.gamma
-
-    spx.enter = -1
-    max_s = 0.0
-    for j = 1:length(spx.r)
-        if spx.s[j] > max_s
-            if spx.basis_status[spx.nonbasic[j]] == BASIS_AT_UPPER && sign(spx.r[j]) > 0
-                spx.enter = spx.nonbasic[j]
-                spx.enter_pos = j
-                max_s = spx.s[j]
-            elseif spx.basis_status[spx.nonbasic[j]] == BASIS_AT_LOWER && sign(spx.r[j]) < 0
-                spx.enter = spx.nonbasic[j]
-                spx.enter_pos = j
-                max_s = spx.s[j]
-            end
-        end
+mutable struct SteepestEdgeRule{T} <: AbstractPivotRule{T}
+    gamma::T # steepest-edge weights
+    v::T     # B^{-T} d
+    s::T     # = r.^2 ./ gamma
+    g::T     # invB * A_j
+    function SteepestEdgeRule(T::Type, m::Int, n::Int)
+        rule = new{T}()
+        rule.gamma = T{Float64}(undef, n - m)
+        rule.v = T{Float64}(undef, m)
+        rule.s = T{Float64}(undef, n - m)
+        rule.g = T{Float64}(undef, m)
+        return rule
     end
-    # @show spx.enter, s[spx.enter_pos]
 end
 
-function compute_steepest_edge_weights(spx::SpxData)
+"Find an entering variable based on steepest edge pivot rule"
+function find_entering_variable!(::Type{SteepestEdgeRule}, spx::SpxData)
+
+    pd = spx.pivot_data
     @timeit TO "compute steepest edge weights" begin
         if spx.iter == 1
             for j = 1:length(spx.nonbasic)
-                spx.g .= spx.invB * spx.lpdata.A[:,spx.nonbasic[j]]
-                spx.gamma[j] = spx.g' * spx.g
+                pd.g .= spx.invB * spx.lpdata.A[:,spx.nonbasic[j]]
+                pd.gamma[j] = pd.g' * pd.g
             end
         else
             # Here, spx.enter and spx.leave were obtained from the previous iteration.
             @timeit TO "gamma_j" begin
-                spx.gamma .+= 2 .* (spx.lpdata.A[:,spx.nonbasic]' * spx.invB[:,spx.leave]) .* (spx.lpdata.A[:,spx.nonbasic]' * spx.v) .+ (spx.lpdata.A[:,spx.nonbasic]' * spx.invB[:,spx.leave]).^2 .* (1 + spx.d' * spx.d)
+                pd.gamma .+= 2 .* (spx.lpdata.A[:,spx.nonbasic]' * spx.invB[:,spx.leave]) .* (spx.lpdata.A[:,spx.nonbasic]' * pd.v) .+ (spx.lpdata.A[:,spx.nonbasic]' * spx.invB[:,spx.leave]).^2 .* (1 + spx.d' * spx.d)
             end
 
             @timeit TO "gamm_e" begin
-                spx.g .= spx.invB * spx.lpdata.A[:,spx.nonbasic[spx.enter_pos]]
-                spx.gamma[spx.enter_pos] = spx.g' * spx.g
+                pd.g .= spx.invB * spx.lpdata.A[:,spx.nonbasic[spx.enter_pos]]
+                pd.gamma[spx.enter_pos] = pd.g' * pd.g
             end
         end
     end
-end
 
-# Dantzig's rule
-function pivot_Dantzig(spx::SpxData)
+    # compute the local slope
+    # @show length(spx.r), length(pd.gamma)
+    pd.s .= (spx.r).^2 ./ pd.gamma
+
     spx.enter = -1
+    max_s = 0.0
+    for j = 1:length(spx.r)
+        if pd.s[j] > max_s
+            if spx.basis_status[spx.nonbasic[j]] == Basis_At_Upper && sign(spx.r[j]) > 0
+                spx.enter = spx.nonbasic[j]
+                spx.enter_pos = j
+                max_s = pd.s[j]
+            elseif spx.basis_status[spx.nonbasic[j]] == Basis_At_Lower && sign(spx.r[j]) < 0
+                spx.enter = spx.nonbasic[j]
+                spx.enter_pos = j
+                max_s = pd.s[j]
+            end
+        end
+    end
+    # @show spx.enter, s[spx.enter_pos]
+end
+
+"""
+    DantzigRule
 
-    # compute reduced cost
-    compute_reduced_cost(spx)
+Struct for Dantzig's pivot rule
+"""
+struct DantzigRule{T} <: AbstractPivotRule{T} end
 
+"Find an entering variable based on Dantzig's pivot rule"
+function find_entering_variable!(::Type{DantzigRule}, spx::SpxData)
+    spx.enter = -1
     max_r = 0.0
     for j = 1:length(spx.r)
         if abs(spx.r[j]) > max_r
-            if spx.basis_status[spx.nonbasic[j]] == BASIS_AT_UPPER
+            if spx.basis_status[spx.nonbasic[j]] == Basis_At_Upper
                 if spx.r[j] > 1e-6
                     spx.enter = spx.nonbasic[j]
                     spx.enter_pos = j
@@ -115,30 +130,3 @@ function pivot_Dantzig(spx::SpxData)
     end
     # @show spx.enter, spx.r[spx.enter_pos]
 end
-
-# Pivot to remove artificial variables from basis (based on Bland's rule)
-function pivot_to_remove_artificials(spx::SpxData)
-    spx.enter = -1
-
-    # compute reduced cost
-    compute_reduced_cost(spx)
-
-    min_enter = spx.start_artvars
-    for j = 1:length(spx.r)
-        if spx.nonbasic[j] < min_enter
-            if spx.basis_status[spx.nonbasic[j]] == BASIS_AT_UPPER
-                if spx.r[j] > -1e-10
-                    spx.enter = spx.nonbasic[j]
-                    spx.enter_pos = j
-                    min_enter = spx.enter
-                end
-            else
-                if spx.r[j] < +1e-10
-                    spx.enter = spx.nonbasic[j]
-                    spx.enter_pos = j
-                    min_enter = spx.enter
-                end
-            end
-        end
-    end
-end