From 446cc0fde14ba07054a5598c2259bb0cb8bceb74 Mon Sep 17 00:00:00 2001 From: mtanneau Date: Sat, 26 Jun 2021 16:56:56 -0400 Subject: [PATCH 1/2] Add Presolve tolerances --- src/Presolve/Presolve.jl | 36 ++++++++++++-------- src/Presolve/empty_column.jl | 11 +++--- src/Presolve/empty_row.jl | 7 ++-- src/model.jl | 14 ++++---- test/Presolve/empty_row.jl | 65 ++++++++++++++++++++++++++++++++---- test/runtests.jl | 2 +- 6 files changed, 98 insertions(+), 37 deletions(-) diff --git a/src/Presolve/Presolve.jl b/src/Presolve/Presolve.jl index 389d7077..3778a932 100644 --- a/src/Presolve/Presolve.jl +++ b/src/Presolve/Presolve.jl @@ -1,5 +1,8 @@ Base.@kwdef mutable struct PresolveOptions{T} Level::Int = 1 # Presolve level + + TolerancePFeas::T = sqrt(eps(T)) # Primal feasibility tolerance + ToleranceDFeas::T = sqrt(eps(T)) # Dual feasibility tolerance end """ @@ -28,6 +31,7 @@ whose dual writes mutable struct PresolveData{T} updated::Bool status::TerminationStatus + options::PresolveOptions{T} # Original problem pb0::ProblemData{T} @@ -86,11 +90,15 @@ mutable struct PresolveData{T} # TODO: set of transformations for pre-post crush ops::Vector{PresolveTransformation{T}} - function PresolveData(pb::ProblemData{T}) where{T} + function PresolveData( + pb::ProblemData{T}, + options::PresolveOptions{T}=PresolveOptions{T}() + ) where{T} ps = new{T}() ps.updated = false ps.status = Trm_Unknown + ps.options = options ps.pb0 = pb ps.pb_red = nothing @@ -167,7 +175,7 @@ end # Extract pre-solved problem data, to be passed to the IPM solver function extract_reduced_problem!(ps::PresolveData{T}) where{T} - + pb = ProblemData{T}() pb.ncon = sum(ps.rowflag) @@ -193,7 +201,7 @@ function extract_reduced_problem!(ps::PresolveData{T}) where{T} inew = 0 for (iold, row) in enumerate(ps.pb0.arows) ps.rowflag[iold] || continue - + inew += 1 # Compute new row rind = Vector{Int}(undef, ps.nzrow[iold]) @@ -219,11 +227,11 @@ function extract_reduced_problem!(ps::PresolveData{T}) where{T} jnew = 0 for (jold, col) in enumerate(ps.pb0.acols) ps.colflag[jold] || continue - + jnew += 1 # Compute new row cind = Vector{Int}(undef, ps.nzcol[jold]) - cval = Vector{T}(undef, ps.nzcol[jold]) + cval = Vector{T}(undef, ps.nzcol[jold]) k = 0 for (iold, aij) in zip(col.nzind, col.nzval) @@ -379,7 +387,7 @@ function presolve!(ps::PresolveData{T}) where{T} # Identify row singletons ps.row_singletons = [i for (i, nz) in enumerate(ps.nzrow) if ps.rowflag[i] && nz == 1] - + # II. Passes ps.updated = true npasses = 0 # TODO: maximum number of passes @@ -392,7 +400,7 @@ function presolve!(ps::PresolveData{T}) where{T} ps.status == Trm_Unknown || return ps.status remove_empty_columns!(ps) ps.status == Trm_Unknown || return ps.status - + # Remove all fixed variables # TODO: remove empty variables as well @@ -444,7 +452,7 @@ function presolve!(ps::PresolveData{T}) where{T} ps.solution.z_primal = ps.obj0 ps.solution.z_dual = ps.obj0 end - + # Old <-> new index mapping compute_index_mapping!(ps) @@ -500,7 +508,7 @@ function bounds_consistency_checks!(ps::PresolveData{T}) where{T} ps.status = Trm_PrimalInfeasible ps.updated = true - # Resize problem + # Resize problem compute_index_mapping!(ps) resize!(ps.solution, ps.nrow, ps.ncol) ps.solution.x .= zero(T) @@ -518,7 +526,7 @@ function bounds_consistency_checks!(ps::PresolveData{T}) where{T} i_ = ps.new_con_idx[i] ps.solution.y_lower[i_] = one(T) ps.solution.y_upper[i_] = one(T) - + return end end @@ -529,7 +537,7 @@ function bounds_consistency_checks!(ps::PresolveData{T}) where{T} ps.status = Trm_PrimalInfeasible ps.updated = true - # Resize problem + # Resize problem compute_index_mapping!(ps) resize!(ps.solution, ps.nrow, ps.ncol) ps.solution.x .= zero(T) @@ -547,7 +555,7 @@ function bounds_consistency_checks!(ps::PresolveData{T}) where{T} j_ = ps.new_var_idx[j] ps.solution.s_lower[j_] = one(T) ps.solution.s_upper[j_] = one(T) - + return end end @@ -677,7 +685,7 @@ function remove_dominated_columns!(ps::PresolveData{T}) where{T} @debug "Col $j forces y$i >= $y_" ps.ly[i] = max(ps.ly[i], y_) end - + elseif !isfinite(l) && isfinite(u) # Upper-bounded variable: `aij * yi ≥ cj` if aij > zero(T) @@ -699,4 +707,4 @@ function remove_dominated_columns!(ps::PresolveData{T}) where{T} ps.status == Trm_Unknown || break end return nothing -end \ No newline at end of file +end diff --git a/src/Presolve/empty_column.jl b/src/Presolve/empty_column.jl index ec89edf4..8f746d23 100644 --- a/src/Presolve/empty_column.jl +++ b/src/Presolve/empty_column.jl @@ -13,7 +13,9 @@ function remove_empty_column!(ps::PresolveData{T}, j::Int) where{T} cj = ps.obj[j] @debug "Removing empty column $j" cj lb ub - if cj > zero(T) + ϵ = ps.options.ToleranceDFeas + + if cj > ϵ if isfinite(lb) # Set variable to lower bound # Update objective constant @@ -45,7 +47,7 @@ function remove_empty_column!(ps::PresolveData{T}, j::Int) where{T} return nothing end - elseif cj < zero(T) + elseif cj < -ϵ if isfinite(ub) # Set variable to upper bound # Update objective constant @@ -74,7 +76,7 @@ function remove_empty_column!(ps::PresolveData{T}, j::Int) where{T} 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 @@ -87,7 +89,6 @@ function remove_empty_column!(ps::PresolveData{T}, j::Int) where{T} # Free variable with zero coefficient and empty column push!(ps.ops, EmptyColumn(j, zero(T), zero(T))) end - end # Book=keeping @@ -102,4 +103,4 @@ function postsolve!(sol::Solution{T}, op::EmptyColumn{T}) where{T} sol.s_lower[op.j] = pos_part(op.s) sol.s_upper[op.j] = neg_part(op.s) return nothing -end \ No newline at end of file +end diff --git a/src/Presolve/empty_row.jl b/src/Presolve/empty_row.jl index 9dd20602..94135aa9 100644 --- a/src/Presolve/empty_row.jl +++ b/src/Presolve/empty_row.jl @@ -12,8 +12,9 @@ function remove_empty_row!(ps::PresolveData{T}, i::Int) where{T} # Check bounds lb, ub = ps.lrow[i], ps.urow[i] + ϵ = ps.options.TolerancePFeas - if ub < zero(T) + if ub < -ϵ # Infeasible @debug "Row $i is primal infeasible" ps.status = Trm_PrimalInfeasible @@ -37,7 +38,7 @@ function remove_empty_row!(ps::PresolveData{T}, i::Int) where{T} i_ = ps.new_con_idx[i] ps.solution.y_upper[i_] = one(T) return - elseif lb > zero(T) + elseif lb > ϵ @debug "Row $i is primal infeasible" ps.status = Trm_PrimalInfeasible ps.updated = true @@ -50,7 +51,7 @@ function remove_empty_row!(ps::PresolveData{T}, i::Int) where{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 diff --git a/src/model.jl b/src/model.jl index 6ce3dc13..233aedc0 100644 --- a/src/model.jl +++ b/src/model.jl @@ -55,7 +55,7 @@ function Base.empty!(m::Model{T}) where{T} m.presolve_data = nothing m.solver = nothing m.solution = nothing - + return nothing end @@ -71,7 +71,7 @@ function optimize!(model::Model{T}) where{T} "Number of threads must be > 0 (is $(model.params.Threads))" ) BLAS.set_num_threads(model.params.Threads) - + # Print initial stats if model.params.OutputLevel > 0 @printf "\nProblem info\n" @@ -80,13 +80,13 @@ function optimize!(model::Model{T}) where{T} @printf " Variables : %d\n" model.pbdata.nvar @printf " Non-zeros : %d\n" sum(length.([col.nzind for col in model.pbdata.acols])) end - + pb_ = model.pbdata # Presolve # TODO: improve the if-else ps_options = model.params.Presolve if ps_options.Level > 0 - model.presolve_data = PresolveData(model.pbdata) + model.presolve_data = PresolveData(model.pbdata, ps_options) t_ = @elapsed st = presolve!(model.presolve_data) model.status = st @@ -104,7 +104,7 @@ function optimize!(model::Model{T}) where{T} # Check presolve status if st == Trm_Optimal || st == Trm_PrimalInfeasible || st == Trm_DualInfeasible || st == Trm_PrimalDualInfeasible model.params.OutputLevel > 0 && println("Presolve solved the problem.") - + # Perform post-solve sol0 = Solution{T}(model.pbdata.ncon, model.pbdata.nvar) postsolve!(sol0, model.presolve_data.solution, model.presolve_data) @@ -169,7 +169,7 @@ function _extract_solution!(sol::Solution{T}, sol.is_primal_ray = is_primal_ray sol.is_dual_ray = is_dual_ray τ_ = (is_primal_ray || is_dual_ray) ? one(T) : inv(ipm.pt.τ) - + @. sol.x = ipm.pt.x[1:n] * τ_ @. sol.s_lower = ipm.pt.zl[1:n] * τ_ @. sol.s_upper = ipm.pt.zu[1:n] * τ_ @@ -212,4 +212,4 @@ function _extract_solution!(sol::Solution{T}, end return nothing -end \ No newline at end of file +end diff --git a/test/Presolve/empty_row.jl b/test/Presolve/empty_row.jl index 8485341d..c5e2fb3a 100644 --- a/test/Presolve/empty_row.jl +++ b/test/Presolve/empty_row.jl @@ -57,9 +57,6 @@ function empty_row_tests(T::Type) # (current problem only has 1 row) @test sol.y_lower[1] > zero(T) - test_empty_row_1(T) - test_empty_row_2(T) - return end @@ -95,7 +92,7 @@ function test_empty_row_1(T::Type) sol = ps.solution @test sol.dual_status == Tulip.Sln_InfeasibilityCertificate @test sol.z_primal == sol.z_dual == T(Inf) - + # Check Farkas ray # (current problem only has 1 row) @test sol.y_lower[1] > zero(T) @@ -134,7 +131,7 @@ function test_empty_row_2(T::Type) sol = ps.solution @test sol.dual_status == Tulip.Sln_InfeasibilityCertificate @test sol.z_primal == sol.z_dual == T(Inf) - + # Check Farkas ray # (current problem only has 1 row) @test sol.y_upper[1] > zero(T) @@ -142,8 +139,62 @@ function test_empty_row_2(T::Type) return nothing end +function test_empty_row_tolerances(T::Type) + # Adapted from https://github.com/ds4dm/Tulip.jl/issues/98 + #= + min x + y + z + s.t. x + y + z == 1 + x == ¹/₃ + y == ¹/₃ + z == ¹/₃ + x, y, z, ≥ 0 + + In the absence of numerical tolerances, x, y, and z get eliminated, + but rouding errors cause the first constraint to be 0 == ϵ ≈ 1e-16, + thereby rendering the problem infeasible. + =# + pb = Tulip.ProblemData{T}() + + m, n = 4, 3 + A = sparse( + [1, 1, 1, 2, 3, 4], + [1, 2, 3, 1, 2, 3], + T[1, 1, 1, 1, 1, 1], + m, n + ) + c = ones(T, n) + + Tulip.load_problem!(pb, "test", + true, c, zero(T), + A, T[1, 1//3, 1//3, 1//3], T[1, 1//3, 1//3, 1//3], + zeros(T, n), fill(T(Inf), n), + ["row1", "row2", "row3", "row4"], ["x", "y", "z"] + ) + + ps = Tulip.PresolveData(pb) + Tulip.presolve!(ps) + + @test ps.status == Tulip.Trm_Optimal + @test ps.nrow == 0 + @test ps.ncol == 0 + + # Check solution status & objective value + sol = ps.solution + @test sol.primal_status == Tulip.Sln_Optimal + @test sol.dual_status == Tulip.Sln_Optimal + @test sol.z_primal ≈ 1 + @test sol.z_dual ≈ 1 + + return nothing +end + @testset "Empty row" begin for T in TvTYPES - @testset "$T" begin empty_row_tests(T) end + @testset "$T" begin + empty_row_tests(T) + test_empty_row_1(T) + test_empty_row_2(T) + test_empty_row_tolerances(T) + end end -end \ No newline at end of file +end diff --git a/test/runtests.jl b/test/runtests.jl index 1814f145..1521a7d6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,9 +10,9 @@ import Convex const TvTYPES = [Float32, Float64, BigFloat] -# Check That Tulip.version() matches what's in the Project.toml @testset "Tulip" begin +# Check That Tulip.version() matches what's in the Project.toml tlp_ver = Tulip.version() toml_ver = TOML.parsefile("../Project.toml")["version"] @test tlp_ver == VersionNumber(toml_ver) From 89bcc5cc6d3ed8f2a330b2e7350b0140cab67e00 Mon Sep 17 00:00:00 2001 From: mtanneau Date: Sat, 26 Jun 2021 17:06:30 -0400 Subject: [PATCH 2/2] Bump v0.7.4 --> v0.7.5 --- Project.toml | 2 +- src/Tulip.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 204c5fb1..57cdd75b 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Tulip" uuid = "6dd1b50a-3aae-11e9-10b5-ef983d2400fa" authors = ["Mathieu Tanneau "] -version = "0.7.4" +version = "0.7.5" [deps] CodecBzip2 = "523fee87-0ab8-5b00-afb7-3ecf72e48cfd" diff --git a/src/Tulip.jl b/src/Tulip.jl index ddef9ff3..c0a52ea8 100644 --- a/src/Tulip.jl +++ b/src/Tulip.jl @@ -7,7 +7,7 @@ using SparseArrays using TimerOutputs -version() = v"0.7.4" +version() = v"0.7.5" include("utils.jl")