Skip to content

Commit

Permalink
Add primal and dual tolerances in presolve (#100)
Browse files Browse the repository at this point in the history 
* Add Presolve tolerances

* Bump v0.7.4 --> v0.7.5
  • Loading branch information
@mtanneau
mtanneau authored Jun 26, 2021
1 parent 27fee90 commit 070a9f4
Show file tree
Hide file tree
Showing 8 changed files with 100 additions and 39 deletions.
2 changes: 1 addition & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "Tulip"
uuid = "6dd1b50a-3aae-11e9-10b5-ef983d2400fa"
authors = ["Mathieu Tanneau <[email protected]>"]
version = "0.7.4"
version = "0.7.5"

[deps]
CodecBzip2 = "523fee87-0ab8-5b00-afb7-3ecf72e48cfd"
Expand Down
36 changes: 22 additions & 14 deletions src/Presolve/Presolve.jl
View file
Original file line number Diff line number Diff line change
@@ -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

"""
Expand Down Expand Up @@ -28,6 +31,7 @@ whose dual writes
mutable struct PresolveData{T}
updated::Bool
status::TerminationStatus
options::PresolveOptions{T}

# Original problem
pb0::ProblemData{T}
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand All @@ -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])
Expand All @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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)
Expand All @@ -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
Expand All @@ -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)
Expand All @@ -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
Expand Down Expand Up @@ -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)
Expand All @@ -699,4 +707,4 @@ function remove_dominated_columns!(ps::PresolveData{T}) where{T}
ps.status == Trm_Unknown || break
end
return nothing
end
end
11 changes: 6 additions & 5 deletions src/Presolve/empty_column.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
end
7 changes: 4 additions & 3 deletions src/Presolve/empty_row.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down
2 changes: 1 addition & 1 deletion src/Tulip.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@ using SparseArrays

using TimerOutputs

version() = v"0.7.4"
version() = v"0.7.5"

include("utils.jl")

Expand Down
14 changes: 7 additions & 7 deletions src/model.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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"
Expand All @@ -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

Expand All @@ -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)
Expand Down Expand Up @@ -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] * τ_
Expand Down Expand Up @@ -212,4 +212,4 @@ function _extract_solution!(sol::Solution{T},
end

return nothing
end
end
65 changes: 58 additions & 7 deletions test/Presolve/empty_row.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -134,16 +131,70 @@ 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)

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
end
2 changes: 1 addition & 1 deletion test/runtests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down

2 comments on commit 070a9f4

@mtanneau
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register()

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/39711

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.7.5 -m "<description of version>" 070a9f4dff3bad91bd62fcac11e0b9e84da94df2
git push origin v0.7.5

Please sign in to comment.