diff --git a/src/MathOptIIS.jl b/src/MathOptIIS.jl index 0faa6fb..bbaffda 100644 --- a/src/MathOptIIS.jl +++ b/src/MathOptIIS.jl @@ -7,12 +7,748 @@ module MathOptIIS import MathOptInterface as MOI -include("iis.jl") -include("bound.jl") -include("range.jl") -include("solver.jl") +include("interval.jl") -@deprecate ConflictCount MOI.ConflictCount -@deprecate ConstraintConflictStatus MOI.ConstraintConflictStatus +# ============================================================================== +# User-facing API +# ============================================================================== + +abstract type AbstractAdditionalData end + +struct InfeasibilityData + # IIS constraints set + constraints::Vector{MOI.ConstraintIndex} + # variable-set constraints only for NoData IIS (from the iis solver) + # this will be an empty vector for most types of IIS + maybe_constraints::Vector{MOI.ConstraintIndex} + # indicates if the IIS is irreducible + irreducible::Bool + # additional data + metadata::AbstractAdditionalData + + function InfeasibilityData( + constraints::Vector{<:MOI.ConstraintIndex}, + irreducible::Bool, + metadata::AbstractAdditionalData; + maybe_constraints::Vector{MOI.ConstraintIndex} = MOI.ConstraintIndex[], + ) + return new(constraints, maybe_constraints, irreducible, metadata) + end +end + +struct NoData <: AbstractAdditionalData end + +""" + Optimizer() + +Create a new optimizer object. +""" +mutable struct Optimizer <: MOI.AbstractOptimizer + infeasible_model::Union{MOI.ModelLike,Nothing} + inner_optimizer::Any + # iis attributes + time_limit::Float64 + verbose::Bool + # result data + start_time::Float64 + status::MOI.ConflictStatusCode + results::Vector{InfeasibilityData} + + function Optimizer() + return new( + nothing, + nothing, + Inf, + false, + NaN, + MOI.COMPUTE_CONFLICT_NOT_CALLED, + InfeasibilityData[], + ) + end +end + +# MathOptIIS.InfeasibleModel + +""" + InfeasibleModel() <: MOI.AbstractModelAttribute + +A model attribute for passing the infeasible model to the IIS solver. + +## Example + +```julia +import MathOptIIS +solver = MathOptIIS.Optimizer() +MOI.set(solver, MathOptIIS.InfeasibleModel(), model) +MOI.set(solver, MathOptIIS.InnerOptimizer(), Optimizer) +MOI.compute_conflict!(solver) +``` +""" +struct InfeasibleModel <: MOI.AbstractModelAttribute end + +function MOI.set(optimizer::Optimizer, ::InfeasibleModel, model::MOI.ModelLike) + optimizer.infeasible_model = model + empty!(optimizer.results) + optimizer.status = MOI.COMPUTE_CONFLICT_NOT_CALLED + return +end + +# MathOptIIS.InnerOptimizer + +""" + InnerOptimizer() <: MOI.AbstractOptimizerAttribute + +A optimizer attribute for passing the inner optimizer to the IIS solver. + +## Example + +```julia +import MathOptIIS +solver = MathOptIIS.Optimizer() +MOI.set(solver, MathOptIIS.InfeasibleModel(), model) +MOI.set(solver, MathOptIIS.InnerOptimizer(), Optimizer) +MOI.compute_conflict!(solver) +``` +""" +struct InnerOptimizer <: MOI.AbstractOptimizerAttribute end + +function MOI.set(optimizer::Optimizer, ::InnerOptimizer, inner_optimizer) + optimizer.inner_optimizer = inner_optimizer + return +end + +# MOI.TimeLimitSec + +function MOI.set(optimizer::Optimizer, ::MOI.TimeLimitSec, value::Float64) + optimizer.time_limit = value + return +end + +MOI.get(optimizer::Optimizer, ::MOI.TimeLimitSec) = optimizer.time_limit + +# MOI.Silent + +function MOI.set(optimizer::Optimizer, ::MOI.Silent, value::Bool) + optimizer.verbose = !value + return +end + +MOI.get(optimizer::Optimizer, ::MOI.Silent) = !optimizer.verbose + +# MOI.ConflictStatus + +MOI.get(optimizer::Optimizer, ::MOI.ConflictStatus) = optimizer.status + +# MOI.ConflictCount + +MOI.get(optimizer::Optimizer, ::MOI.ConflictCount) = length(optimizer.results) + +# MOI.ConstraintConflictStatus + +function MOI.get( + optimizer::Optimizer, + attr::MOI.ConstraintConflictStatus, + con::MOI.ConstraintIndex, +) + if !(1 <= attr.conflict_index <= length(optimizer.results)) + return MOI.NOT_IN_CONFLICT + elseif con in optimizer.results[attr.conflict_index].constraints + return MOI.IN_CONFLICT + elseif con in optimizer.results[attr.conflict_index].maybe_constraints + return MOI.MAYBE_IN_CONFLICT + end + return MOI.NOT_IN_CONFLICT +end + +# MathOptIIS.ListOfConstraintIndicesInConflict + +struct ListOfConstraintIndicesInConflict <: MOI.AbstractModelAttribute + conflict_index::Int + + ListOfConstraintIndicesInConflict(conflict_index = 1) = new(conflict_index) +end + +function MOI.get(optimizer::Optimizer, attr::ListOfConstraintIndicesInConflict) + if !(1 <= attr.conflict_index <= length(optimizer.results)) + return MOI.ConstraintIndex[] + end + return optimizer.results[attr.conflict_index].constraints +end + +# MOI.compute_conflict! + +function _in_time(optimizer::Optimizer) + @assert optimizer.start_time != NaN + return time() - optimizer.start_time < optimizer.time_limit +end + +function MOI.compute_conflict!(optimizer::Optimizer) + optimizer.status = MOI.NO_CONFLICT_FOUND + empty!(optimizer.results) + optimizer.start_time = time() + if optimizer.verbose + println("[MathOptIIS] starting compute_conflict!") + end + if _feasibility_check(optimizer, optimizer.infeasible_model) + optimizer.status = MOI.NO_CONFLICT_EXISTS + return optimizer.results + end + # Step 1: check for inconsistent variable bounds + variable_info = _bound_infeasibility!(optimizer) + if !isempty(optimizer.results) + optimizer.status = MOI.CONFLICT_FOUND + return + elseif !_in_time(optimizer) + return + end + # Step 2: check for inconsistent constraints based on variable bounds + _range_infeasibility!(optimizer, variable_info) + if !isempty(optimizer.results) + optimizer.status = MOI.CONFLICT_FOUND + return + elseif !_in_time(optimizer) + return + end + # Step 3: elastic filter + _elastic_filter(optimizer, variable_info) + if !isempty(optimizer.results) + optimizer.status = MOI.CONFLICT_FOUND + end + if optimizer.verbose + println( + "[MathOptIIS] elastic filter found $(length(optimizer.results)) infeasible subsets", + ) + end + return +end + +# ============================================================================== +# Step 0: feasibility check +# ============================================================================== + +function _feasibility_check( + optimizer::Optimizer, + infeasible_model::MOI.ModelLike, +) + termination_status = MOI.get(infeasible_model, MOI.TerminationStatus()) + if optimizer.verbose + println("[MathOptIIS] model termination status: $(termination_status)") + end + if termination_status in + (MOI.OTHER_ERROR, MOI.INVALID_MODEL, MOI.OPTIMIZE_NOT_CALLED) + return false # because we can assert it is feasible + end + primal_status = MOI.get(infeasible_model, MOI.PrimalStatus()) + if optimizer.verbose + println("[MathOptIIS] model primal status: $(primal_status)") + end + if primal_status in (MOI.FEASIBLE_POINT, MOI.NEARLY_FEASIBLE_POINT) && !( + termination_status in + (MOI.INFEASIBLE, MOI.ALMOST_INFEASIBLE, MOI.LOCALLY_INFEASIBLE) + ) + return true + end + return false +end + +# ============================================================================== +# Step 1: check variable bounds and integrality restrictions for violations +# ============================================================================== + +struct BoundsData{T} <: AbstractAdditionalData + lower_bound::T + upper_bound::T +end + +struct IntegralityData{T} <: AbstractAdditionalData + lower_bound::T + upper_bound::T + set::Union{MOI.Integer,MOI.ZeroOne} +end + +function _bound_infeasibility!(optimizer::Optimizer) + return _bound_infeasibility!(optimizer, optimizer.infeasible_model, Float64) +end + +mutable struct _VariableInfo{T} + lower::Union{Nothing,MOI.AbstractScalarSet} + upper::Union{Nothing,MOI.AbstractScalarSet} + integer::Bool + zero_one::Bool + + _VariableInfo{T}() where {T} = new{T}(nothing, nothing, false, false) +end + +_update_info!(info::_VariableInfo, s::MOI.LessThan) = (info.upper = s) + +_update_info!(info::_VariableInfo, s::MOI.GreaterThan) = (info.lower = s) + +function _update_info!(info::_VariableInfo, s::Union{MOI.EqualTo,MOI.Interval}) + info.lower = info.upper = s + return +end + +_update_info!(info::_VariableInfo, ::MOI.Integer) = (info.integer = true) + +_update_info!(info::_VariableInfo, ::MOI.ZeroOne) = (info.zero_one = true) + +function _update_info!(info, model, ::Type{S}) where {S} + for ci in MOI.get(model, MOI.ListOfConstraintIndices{MOI.VariableIndex,S}()) + f = MOI.get(model, MOI.ConstraintFunction(), ci) + s = MOI.get(model, MOI.ConstraintSet(), ci) + _update_info!(info[f], s) + end + return +end + +function _ci(x::MOI.VariableIndex, ::S) where {S} + return MOI.ConstraintIndex{MOI.VariableIndex,S}(x.value) +end + +_lower(::Type{T}, ::Nothing) where {T} = typemin(T) +_lower(::Type{T}, s::Union{MOI.GreaterThan,MOI.Interval}) where {T} = s.lower +_lower(::Type{T}, s::MOI.EqualTo) where {T} = s.value + +_upper(::Type{T}, ::Nothing) where {T} = typemax(T) +_upper(::Type{T}, s::Union{MOI.LessThan,MOI.Interval}) where {T} = s.upper +_upper(::Type{T}, s::MOI.EqualTo) where {T} = s.value + +function _check_conflict(x::MOI.VariableIndex, info::_VariableInfo{T}) where {T} + lb, ub = _lower(T, info.lower), _upper(T, info.upper) + if ub < lb + return InfeasibilityData( + MOI.ConstraintIndex[_ci(x, info.lower), _ci(x, info.upper)], + true, + BoundsData{T}(lb, ub), + ) + elseif info.integer + # I don't think this `if` is correct? + if abs(ub - lb) < 1 && ceil(Int, ub) == ceil(Int, lb) + metadata = IntegralityData{T}(lb, ub, MOI.Integer()) + con = MOI.ConstraintIndex{MOI.VariableIndex,MOI.Integer}(x.value) + c = MOI.ConstraintIndex[con, _ci(x, info.lower), _ci(x, info.upper)] + return InfeasibilityData(c, true, metadata) + end + elseif info.zero_one + con = MOI.ConstraintIndex{MOI.VariableIndex,MOI.ZeroOne}(x.value) + if 0 < lb && ub < 1 + c = MOI.ConstraintIndex[con, _ci(x, info.lower), _ci(x, info.upper)] + metadata = IntegralityData{T}(lb, ub, MOI.ZeroOne()) + return InfeasibilityData(c, true, metadata) + elseif 1 < lb + c = MOI.ConstraintIndex[con, _ci(x, info.lower)] + metadata = IntegralityData{T}(lb, typemax(T), MOI.ZeroOne()) + return InfeasibilityData(c, true, metadata) + elseif ub < 0 + c = MOI.ConstraintIndex[con, _ci(x, info.upper)] + metadata = IntegralityData{T}(typemin(T), ub, MOI.ZeroOne()) + return InfeasibilityData(c, true, metadata) + end + end + return +end + +function _bound_infeasibility!( + optimizer::Optimizer, + infeasible_model::MOI.ModelLike, + ::Type{T}, +) where {T} + if optimizer.verbose + println("[MathOptIIS] starting bound analysis") + end + variable_info = Dict( + x => _VariableInfo{T}() for + x in MOI.get(infeasible_model, MOI.ListOfVariableIndices()) + ) + _update_info!(variable_info, infeasible_model, MOI.LessThan{T}) + _update_info!(variable_info, infeasible_model, MOI.GreaterThan{T}) + _update_info!(variable_info, infeasible_model, MOI.EqualTo{T}) + _update_info!(variable_info, infeasible_model, MOI.Interval{T}) + _update_info!(variable_info, infeasible_model, MOI.Integer) + _update_info!(variable_info, infeasible_model, MOI.ZeroOne) + results = InfeasibilityData[] + for (x, info) in variable_info + if (conflict = _check_conflict(x, info)) !== nothing + push!(results, conflict) + end + end + append!(optimizer.results, results) + if optimizer.verbose + println( + "[MathOptIIS] bound analysis found $(length(results)) infeasible subsets", + ) + end + return variable_info +end + +# ============================================================================== +# Step 2: propagate variable bounds through functions to detect obvious errors +# ============================================================================== + +struct RangeData{T} <: AbstractAdditionalData + lower_bound::T + upper_bound::T + set::MOI.AbstractSet +end + +_supports_interval(::Type{T}) where {T} = false + +function _range_infeasibility!( + optimizer::Optimizer, + variable_info::Dict{MOI.VariableIndex,_VariableInfo{T}}, +) where {T} + if optimizer.verbose + println("[MathOptIIS] starting range analysis") + end + variables = Dict{MOI.VariableIndex,Interval{T}}( + x => Interval(_lower(T, info.lower), _upper(T, info.upper)) for + (x, info) in variable_info + ) + results = InfeasibilityData[] + for (F, S) in + MOI.get(optimizer.infeasible_model, MOI.ListOfConstraintTypesPresent()) + if _supports_interval(F) && _supports_interval(S) + _range_infeasibility!( + optimizer, + optimizer.infeasible_model, + variable_info, + variables, + F, + S, + results, + ) + end + end + append!(optimizer.results, results) + if optimizer.verbose + println( + "[MathOptIIS] range analysis found $(length(results)) infeasible subsets", + ) + end + return +end + +function _range_infeasibility!( + optimizer::Optimizer, + infeasible_model::MOI.ModelLike, + variable_info::Dict{MOI.VariableIndex,_VariableInfo{T}}, + variables::Dict{MOI.VariableIndex,Interval{T}}, + ::Type{F}, + ::Type{S}, + results::Vector{InfeasibilityData}, +) where {T,F,S} + for con in MOI.get(infeasible_model, MOI.ListOfConstraintIndices{F,S}()) + if !_in_time(optimizer) + return + end + func = MOI.get(infeasible_model, MOI.ConstraintFunction(), con) + cons = Set{MOI.ConstraintIndex}() + interval = _compute_interval(variables, func, variable_info, cons) + set = MOI.get(infeasible_model, MOI.ConstraintSet(), con)::S + if !_valid_range(set, interval) + push!(cons, con) + metadata = RangeData{T}(interval.lo, interval.hi, set) + push!(results, InfeasibilityData(collect(cons), true, metadata)) + end + end + return +end + +_supports_interval(::Type{MOI.ScalarAffineFunction{T}}) where {T} = true + +function _compute_interval( + variables::Dict{MOI.VariableIndex,Interval{T}}, + f::MOI.ScalarAffineFunction, + variable_info::Dict{MOI.VariableIndex,_VariableInfo{T}}, + cons::Set{MOI.ConstraintIndex}, +) where {T} + out = convert(Interval{T}, f.constant) + for t in f.terms + out += t.coefficient * variables[t.variable] + if (s = variable_info[t.variable].lower) !== nothing + push!(cons, _ci(t.variable, s)) + end + if (s = variable_info[t.variable].upper) !== nothing + push!(cons, _ci(t.variable, s)) + end + end + return out +end + +_supports_interval(::Type{MOI.EqualTo{T}}) where {T} = true + +function _valid_range(set::MOI.EqualTo, interval) + return interval.lo <= set.value <= interval.hi +end + +_supports_interval(::Type{MOI.LessThan{T}}) where {T} = true + +_valid_range(set::MOI.LessThan, interval) = set.upper >= interval.lo + +_supports_interval(::Type{MOI.GreaterThan{T}}) where {T} = true + +_valid_range(set::MOI.GreaterThan, interval) = interval.hi >= set.lower + +_supports_interval(::Type{MOI.Interval{T}}) where {T} = true + +function _valid_range(set::MOI.Interval, interval) + return set.upper >= interval.lo && set.lower <= interval.hi +end + +# ============================================================================== +# Step 3: compute a proper IIS using an elastic filter +# ============================================================================== + +function _fix_slack(model, func::MOI.ScalarAffineFunction{T}) where {T} + for term in func.terms + MOI.add_constraint(model, term.variable, MOI.LessThan(zero(T))) + end + return +end + +function _unfix_slack(model, func::MOI.ScalarAffineFunction{T}) where {T} + for term in func.terms + x = term.variable + ci = MOI.ConstraintIndex{MOI.VariableIndex,MOI.LessThan{T}}(x.value) + @assert MOI.is_valid(model, ci) + MOI.delete(model, ci) + end + return +end + +function _elastic_filter( + optimizer::Optimizer, + variable_info::Dict{MOI.VariableIndex,_VariableInfo{T}}, +) where {T} + if optimizer.verbose + println("[MathOptIIS] starting elastic filter") + end + if optimizer.inner_optimizer === nothing + println( + "[MathOptIIS] elastic filter cannot continue because no optimizer was provided", + ) + return + end + model = MOI.instantiate(optimizer.inner_optimizer) + MOI.set(model, MOI.Silent(), true) + index_map = MOI.copy_to(model, optimizer.infeasible_model) + new_to_old_index_map = MOI.IndexMap() + for (k, v) in index_map + new_to_old_index_map[v] = k + end + # ========================================================================== + # Step 1: relax integrality and solve. + # ========================================================================== + relax_info = _relax_integrality(model, variable_info) + MOI.optimize!(model) + if !_is_feasible(model) + # The relaxed problem is infeasible. This is great news, because we can + # continue to find an IIS from the continuous relaxation. + if MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE + if optimizer.verbose + println( + "[MathOptIIS] using INFEASIBILITY_CERTIFICATE to construct infeasible subset", + ) + end + # If there is an infeasibility certificate, the non-zero rows an IIS + _iis_from_certificate( + optimizer, + model, + new_to_old_index_map, + relax_info, + ) + return + end + else + # The relaxed problem is feasible. This is significantly more difficult + # to deal with because we need to add back in the integrality + # restrictions. + for x in relax_info.integer + MOI.add_constraint(model, x, MOI.Integer()) + end + for (ci, x) in relax_info.binary + MOI.delete(model, ci) + MOI.add_constraint(model, x, MOI.ZeroOne()) + end + end + # ========================================================================== + # Step 2: relax constraints and solve elastic filter + # + # TODO(odow): we can improve this step by adding support for new constraints + # in MOI.Utilities.PenaltyRelaxation + # ========================================================================== + constraint_to_affine = MOI.modify( + model, + MOI.Utilities.PenaltyRelaxation(; default = one(T), warn = false), + ) + slack_obj = zero(MOI.ScalarAffineFunction{T}) + for v in values(constraint_to_affine) + MOI.Utilities.operate!(+, T, slack_obj, v) + end + MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE) + MOI.set(model, MOI.ObjectiveFunction{typeof(slack_obj)}(), slack_obj) + iis_candidate = Set{MOI.ConstraintIndex}() + constraints_to_check = Set(keys(constraint_to_affine)) + while !isempty(constraints_to_check) + if !_in_time(optimizer) + # We've run out of time. Throw all remaining constraints in as well + append!(iis_candidate, constraints_to_check) + break + end + MOI.optimize!(model) + if !_is_feasible(model) + break + end + for con in constraints_to_check + func = constraint_to_affine[con] + for t in func.terms + if MOI.get(model, MOI.VariablePrimal(), t.variable) > 0 + _fix_slack(model, func) + push!(iis_candidate, con) + delete!(constraints_to_check, con) + break + end + end + end + end + # ========================================================================== + # Step 3: did the elastic filter work? + # ========================================================================== + if isempty(iis_candidate) + # Uh'oh. The penalty relaxed model is infeasible, even with our slacks. + # This means that there must be other constraints that are problematic. + return nothing + end + # ========================================================================== + # Step 4: a deletion filter + # + # TODO(odow): we can improve this step by (un)fixing in batches. Instead of + # looping through the candidates one-by-one, we could use a bisection search + # to reduce how many searches we need to loop through. + # ========================================================================== + for con in copy(iis_candidate) + if !_in_time(optimizer) + break + end + _unfix_slack(model, constraint_to_affine[con]) + MOI.optimize!(model) + if _is_feasible(model) + _fix_slack(model, constraint_to_affine[con]) + else + delete!(iis_candidate, con) + end + end + iis = MOI.ConstraintIndex[new_to_old_index_map[c] for c in iis_candidate] + maybe_constraints = + _get_variables_in_constraints(optimizer.infeasible_model, iis) + push!( + optimizer.results, + InfeasibilityData(iis, true, NoData(); maybe_constraints), + ) + return +end + +function _is_feasible(model::MOI.ModelLike) + status = MOI.get(model, MOI.TerminationStatus()) + if status in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED) + return true + elseif status in (MOI.INFEASIBLE,) + return false + else + error("Unsupported termination status in IIS algorithm: $status") + end +end + +function _iis_from_certificate( + optimizer, + model::MOI.ModelLike, + index_map, + relax_info, +) + iis = MOI.ConstraintIndex[] + for (F, S) in MOI.get(model, MOI.ListOfConstraintTypesPresent()) + for ci in MOI.get(model, MOI.ListOfConstraintIndices{F,S}()) + if iszero(MOI.get(model, MOI.ConstraintDual(), ci)) + continue + elseif haskey(relax_info.binary, ci) + x = relax_info.binary[ci] + # 1.0 * x in [0, 1] appears, but the constraint is really ZeroOne + push!(iis, _ci(index_map[x], MOI.ZeroOne())) + else + push!(iis, index_map[ci]) + end + end + end + push!(optimizer.results, InfeasibilityData(iis, true, NoData())) + return +end + +function _relax_integrality( + model::MOI.ModelLike, + variable_info::Dict{MOI.VariableIndex,_VariableInfo{T}}, +) where {T} + F = MOI.VariableIndex + ret = (; + integer = MOI.VariableIndex[], + binary = Dict{ + MOI.ConstraintIndex{MOI.ScalarAffineFunction{T},MOI.Interval{T}}, + MOI.VariableIndex, + }(), + ) + for (x, info) in variable_info + if info.integer + MOI.delete(model, MOI.ConstraintIndex{F,MOI.Integer}(x.value)) + push!(ret.integer, x) + elseif info.zero_one + MOI.delete(model, MOI.ConstraintIndex{F,MOI.ZeroOne}(x.value)) + zero_one_set = MOI.Interval(zero(T), one(T)) + ci = MOI.add_constraint(model, one(T) * x, zero_one_set) + ret.binary[ci] = x + end + end + return ret +end + +function _get_variables_in_constraints( + model::MOI.ModelLike, + con::Vector{MOI.ConstraintIndex}, +) + variables = Set{MOI.VariableIndex}() + for c in con + _get_variables_in_constraints!(model, c, variables) + end + variable_constraints = MOI.ConstraintIndex[] + for (F, S) in MOI.get(model, MOI.ListOfConstraintTypesPresent()) + if !(F <: MOI.VariableIndex) + continue + end + for con in MOI.get(model, MOI.ListOfConstraintIndices{F,S}()) + if MOI.get(model, MOI.ConstraintFunction(), con) in variables + push!(variable_constraints, con) + end + end + end + return variable_constraints +end + +function _get_variables_in_constraints!( + model::MOI.ModelLike, + con::MOI.ConstraintIndex{<:MOI.ScalarAffineFunction}, + variables::Set{MOI.VariableIndex}, +) + f = MOI.get(model, MOI.ConstraintFunction(), con) + for term in f.terms + push!(variables, term.variable) + end + return +end + +function _get_variables_in_constraints!( + ::MOI.ModelLike, + ::MOI.ConstraintIndex, + ::Set{MOI.VariableIndex}, +) + return # skip +end end # module MathOptIIS diff --git a/src/bound.jl b/src/bound.jl deleted file mode 100644 index 8497bbf..0000000 --- a/src/bound.jl +++ /dev/null @@ -1,149 +0,0 @@ -# Copyright (c) 2025: Joaquim Dias Garcia, Oscar Dowson and contributors -# -# Use of this source code is governed by an MIT-style license that can be found -# in the LICENSE.md file or at https://opensource.org/licenses/MIT. - -function _bound_infeasibility!(optimizer::Optimizer, ::Type{T}) where {T} - return _bound_infeasibility!(optimizer.original_model, optimizer.results, T) -end - -function _bound_infeasibility!( - original_model::MOI.ModelLike, - results::Vector{InfeasibilityData}, - ::Type{T}, -) where {T} - variables = Dict{MOI.VariableIndex,Interval{T}}() - variable_indices = MOI.get(original_model, MOI.ListOfVariableIndices()) - lb = Dict{MOI.VariableIndex,T}() - lb_con = Dict{MOI.VariableIndex,MOI.ConstraintIndex}() - ub = Dict{MOI.VariableIndex,T}() - ub_con = Dict{MOI.VariableIndex,MOI.ConstraintIndex}() - for con in MOI.get( - original_model, - MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.EqualTo{T}}(), - ) - set = MOI.get(original_model, MOI.ConstraintSet(), con) - func = MOI.get(original_model, MOI.ConstraintFunction(), con) - lb[func] = set.value - lb_con[func] = con - ub[func] = set.value - ub_con[func] = con - end - for con in MOI.get( - original_model, - MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.LessThan{T}}(), - ) - set = MOI.get(original_model, MOI.ConstraintSet(), con) - func = MOI.get(original_model, MOI.ConstraintFunction(), con) - # lb[func] = typemin(T) - ub[func] = set.upper - ub_con[func] = con - end - for con in MOI.get( - original_model, - MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.GreaterThan{T}}(), - ) - set = MOI.get(original_model, MOI.ConstraintSet(), con) - func = MOI.get(original_model, MOI.ConstraintFunction(), con) - lb[func] = set.lower - lb_con[func] = con - # ub[func] = typemax(T) - end - for con in MOI.get( - original_model, - MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.Interval{T}}(), - ) - set = MOI.get(original_model, MOI.ConstraintSet(), con) - func = MOI.get(original_model, MOI.ConstraintFunction(), con) - lb[func] = set.lower - lb_con[func] = con - ub[func] = set.upper - ub_con[func] = con - end - # for con in MOI.get(original_model, MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.SemiContinuous{T}}()) - # set = MOI.get(original_model, MOI.ConstraintSet(), con) - # func = MOI.get(original_model, MOI.ConstraintFunction(), con) - # lb[func] = 0 # set.lower - # ub[func] = set.upper - # end - # for con in MOI.get(original_model, MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.SemiInteger{T}}()) - # set = MOI.get(original_model, MOI.ConstraintSet(), con) - # func = MOI.get(original_model, MOI.ConstraintFunction(), con) - # lb[func] = 0 #set.lower - # ub[func] = set.upper - # end - bounds_consistent = true - for con in MOI.get( - original_model, - MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.Integer}(), - ) - func = MOI.get(original_model, MOI.ConstraintFunction(), con) - _lb, _ub = get(lb, func, typemin(T)), get(ub, func, typemax(T)) - if abs(_ub - _lb) < 1 && ceil(_ub) == ceil(_lb) - push!( - results, - InfeasibilityData( - MOI.ConstraintIndex[con, lb_con[func], ub_con[func]], - true, - IntegralityData(_lb, _ub, MOI.Integer()), - ), - ) - bounds_consistent = false - end - end - for con in MOI.get( - original_model, - MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.ZeroOne}(), - ) - func = MOI.get(original_model, MOI.ConstraintFunction(), con) - _lb, _ub = get(lb, func, typemin(T)), get(ub, func, typemax(T)) - if 0 < _lb && _ub < 1 - push!( - results, - InfeasibilityData( - MOI.ConstraintIndex[con, lb_con[func], ub_con[func]], - true, - IntegralityData(_lb, _ub, MOI.ZeroOne()), - ), - ) - bounds_consistent = false - elseif _lb > 1 - push!( - results, - InfeasibilityData( - MOI.ConstraintIndex[con, lb_con[func]], - true, - IntegralityData(_lb, typemax(T), MOI.ZeroOne()), - ), - ) - bounds_consistent = false - elseif _ub < 0 - push!( - results, - InfeasibilityData( - MOI.ConstraintIndex[con, ub_con[func]], - true, - IntegralityData(typemin(T), _ub, MOI.ZeroOne()), - ), - ) - bounds_consistent = false - end - end - for var in variable_indices - _lb, _ub = get(lb, var, typemin(T)), get(ub, var, typemax(T)) - if _lb > _ub - push!( - results, - InfeasibilityData( - MOI.ConstraintIndex[lb_con[var], ub_con[var]], - true, - BoundsData(_lb, _ub), - ), - ) - bounds_consistent = false - else - variables[var] = Interval(_lb, _ub) - end - end - return bounds_consistent, variables, lb_con, ub_con -end diff --git a/src/iis.jl b/src/iis.jl deleted file mode 100644 index 75254ee..0000000 --- a/src/iis.jl +++ /dev/null @@ -1,380 +0,0 @@ -# Copyright (c) 2025: Joaquim Dias Garcia, Oscar Dowson and contributors -# -# Use of this source code is governed by an MIT-style license that can be found -# in the LICENSE.md file or at https://opensource.org/licenses/MIT. - -abstract type AbstractAdditionalData end - -struct InfeasibilityData - # IIS constraints set - constraints::Vector{MOI.ConstraintIndex} - # variable-set constraints only for NoData IIS (from the iis solver) - # this will be an empty vector for most types of IIS - maybe_constraints::Vector{MOI.ConstraintIndex} - # indicates if the IIS is irreducible - irreducible::Bool - # additional data - metadata::AbstractAdditionalData - function InfeasibilityData( - constraints::Vector{<:MOI.ConstraintIndex}, - irreducible::Bool, - metadata::AbstractAdditionalData; - maybe_constraints::Vector{MOI.ConstraintIndex} = MOI.ConstraintIndex[], - ) - return new(constraints, maybe_constraints, irreducible, metadata) - end -end - -struct BoundsData <: AbstractAdditionalData - lower_bound::Float64 - upper_bound::Float64 -end - -struct IntegralityData <: AbstractAdditionalData - lower_bound::Float64 - upper_bound::Float64 - set::Union{MOI.Integer,MOI.ZeroOne}#, MOI.Semicontinuous{T}, MOI.Semiinteger{T}} -end - -struct RangeData <: AbstractAdditionalData - lower_bound::Float64 - upper_bound::Float64 - set::Union{<:MOI.EqualTo,<:MOI.LessThan,<:MOI.GreaterThan} -end - -struct NoData <: AbstractAdditionalData end - -Base.@kwdef mutable struct Optimizer - original_model::Union{MOI.ModelLike,Nothing} = nothing - # - # iterative solver data - optimizer::Any = nothing # MOI.ModelLike - optimizer_attributes::Dict{Any,Any} = Dict{Any,Any}() - # - # iis attributes - time_limit::Float64 = Inf - verbose::Bool = false - skip_feasibility_check::Bool = false - stop_if_infeasible_bounds::Bool = true - stop_if_infeasible_ranges::Bool = true - deletion_filter::Bool = true - elastic_filter_tolerance::Float64 = 1e-5 - ignore_integrality::Bool = false - # - # result data - start_time::Float64 = NaN - status::MOI.ConflictStatusCode = MOI.COMPUTE_CONFLICT_NOT_CALLED - results::Vector{InfeasibilityData} = InfeasibilityData[] -end - -struct InfeasibleModel end - -function MOI.set(optimizer::Optimizer, ::InfeasibleModel, model::MOI.ModelLike) - optimizer.original_model = model - empty!(optimizer.results) - optimizer.status = MOI.COMPUTE_CONFLICT_NOT_CALLED - return -end - -# function MOI.get(optimizer::Optimizer, ::InfeasibleModel) -# return optimizer.original_model -# end - -struct InnerOptimizer end - -function MOI.set(optimizer::Optimizer, ::InnerOptimizer, solver) - optimizer.optimizer = solver - return -end - -# function MOI.get(optimizer::Optimizer, ::InnerOptimizer) -# return optimizer.optimizer -# end - -struct InnerOptimizerAttribute - attr::MOI.AbstractOptimizerAttribute -end - -function MOI.set(optimizer::Optimizer, attr::InnerOptimizerAttribute, value) - optimizer.optimizer_attributes[attr.attr] = value - return -end - -# function MOI.get( -# optimizer::Optimizer, -# attr::InnerOptimizerAttribute, -# ) -# return optimizer.optimizer_attributes[attr.attr] -# end - -function MOI.set(optimizer::Optimizer, ::MOI.TimeLimitSec, value::Float64) - optimizer.time_limit = value - return -end - -function MOI.get(optimizer::Optimizer, ::MOI.TimeLimitSec) - return optimizer.time_limit -end - -function MOI.set(optimizer::Optimizer, ::MOI.Silent, value::Bool) - optimizer.verbose = !value - return -end - -function MOI.get(optimizer::Optimizer, ::MOI.Silent) - return !optimizer.verbose -end - -struct SkipFeasibilityCheck <: MOI.AbstractOptimizerAttribute end - -function MOI.set(optimizer::Optimizer, ::SkipFeasibilityCheck, value::Bool) - optimizer.skip_feasibility_check = value - return -end - -function MOI.get(optimizer::Optimizer, ::SkipFeasibilityCheck) - return optimizer.skip_feasibility_check -end - -struct StopIfInfeasibleBounds end - -function MOI.set(optimizer::Optimizer, ::StopIfInfeasibleBounds, value::Bool) - optimizer.stop_if_infeasible_bounds = value - return -end - -function MOI.get(optimizer::Optimizer, ::StopIfInfeasibleBounds) - return optimizer.stop_if_infeasible_bounds -end - -struct StopIfInfeasibleRanges end - -function MOI.set(optimizer::Optimizer, ::StopIfInfeasibleRanges, value::Bool) - optimizer.stop_if_infeasible_ranges = value - return -end - -function MOI.get(optimizer::Optimizer, ::StopIfInfeasibleRanges) - return optimizer.stop_if_infeasible_ranges -end - -struct DeletionFilter end - -function MOI.set(optimizer::Optimizer, ::DeletionFilter, value::Bool) - optimizer.deletion_filter = value - return -end - -function MOI.get(optimizer::Optimizer, ::DeletionFilter) - return optimizer.deletion_filter -end - -struct ElasticFilterTolerance end - -function MOI.set(optimizer::Optimizer, ::ElasticFilterTolerance, value::Float64) - optimizer.elastic_filter_tolerance = value - return -end - -function MOI.get(optimizer::Optimizer, ::ElasticFilterTolerance) - return optimizer.elastic_filter_tolerance -end - -struct ElasticFilterIgnoreIntegrality end - -function MOI.set( - optimizer::Optimizer, - ::ElasticFilterIgnoreIntegrality, - value::Bool, -) - optimizer.ignore_integrality = value - return -end - -function MOI.get(optimizer::Optimizer, ::ElasticFilterIgnoreIntegrality) - return optimizer.ignore_integrality -end - -MOI.get(optimizer::Optimizer, ::MOI.ConflictStatus) = optimizer.status - -MOI.get(optimizer::Optimizer, ::MOI.ConflictCount) = length(optimizer.results) - -function MOI.get( - optimizer::Optimizer, - attr::MOI.ConstraintConflictStatus, - con::MOI.ConstraintIndex, -) - if attr.conflict_index > length(optimizer.results) - return MOI.NOT_IN_CONFLICT # or error - end - if con in optimizer.results[attr.conflict_index].constraints - return MOI.IN_CONFLICT - elseif con in optimizer.results[attr.conflict_index].maybe_constraints - return MOI.MAYBE_IN_CONFLICT - end - return MOI.NOT_IN_CONFLICT -end - -struct ListOfConstraintIndicesInConflict <: MOI.AbstractModelAttribute - conflict_index::Int - ListOfConstraintIndicesInConflict(conflict_index = 1) = new(conflict_index) -end - -function MOI.get(optimizer::Optimizer, attr::ListOfConstraintIndicesInConflict) - if attr.conflict_index > length(optimizer.results) - return MOI.ConstraintIndex[] - end - return optimizer.results[attr.conflict_index].constraints -end - -function _in_time(optimizer::Optimizer) - @assert optimizer.start_time != NaN - return time() - optimizer.start_time < optimizer.time_limit -end - -function MOI.compute_conflict!(optimizer::Optimizer) - optimizer.status = MOI.NO_CONFLICT_FOUND - empty!(optimizer.results) - optimizer.start_time = time() - if optimizer.verbose - println("Starting MathOptIIS IIS search.") - end - T = Float64 - is_feasible = _feasibility_check(optimizer, optimizer.original_model) - if is_feasible && !optimizer.skip_feasibility_check - optimizer.status = MOI.NO_CONFLICT_EXISTS - return optimizer.results - end - if optimizer.verbose - println("Starting bound analysis.") - end - bounds_consistent, variables, lb_con, ub_con = - _bound_infeasibility!(optimizer, T) - bound_infeasibilities = length(optimizer.results) - if optimizer.verbose - println( - "Complete bound analysis found $bound_infeasibilities infeasibilities.", - ) - end - if !isempty(optimizer.results) - optimizer.status = MOI.CONFLICT_FOUND - end - # check PSD diagonal >= 0 ? - # other cones? - if (!bounds_consistent && optimizer.stop_if_infeasible_bounds) || - !_in_time(optimizer) - return - end - # second layer of infeasibility analysis is constraint range analysis - if optimizer.verbose - println("Starting range analysis.") - end - range_consistent = - _range_infeasibility!(optimizer, T, variables, lb_con, ub_con) - range_infeasibilities = length(optimizer.results) - bound_infeasibilities - if optimizer.verbose - println( - "Complete range analysis found $range_infeasibilities infeasibilities.", - ) - end - if !isempty(optimizer.results) - optimizer.status = MOI.CONFLICT_FOUND - end - if (!range_consistent && optimizer.stop_if_infeasible_ranges) || - !_in_time(optimizer) - return - end - # third layer is an IIS resolver - if optimizer.optimizer === nothing - println( - "IIS resolver cannot continue because no optimizer was provided", - ) - return - end - if optimizer.verbose - println("Starting elastic filter solver.") - end - iis = _elastic_filter(optimizer) - # for now, only one iis is computed - if iis !== nothing - maybe_constraints = - _get_variables_in_constraints(optimizer.original_model, iis) - push!( - optimizer.results, - InfeasibilityData(iis, true, NoData(); maybe_constraints), - ) - optimizer.status = MOI.CONFLICT_FOUND - end - if optimizer.verbose - iis_infeasibilities = iis === nothing ? 0 : 1 - println( - "Complete elastic filter solver found $iis_infeasibilities infeasibilities.", - ) - end - - return -end - -function _get_variables_in_constraints( - model::MOI.ModelLike, - con::Vector{MOI.ConstraintIndex}, -) - variables = Set{MOI.VariableIndex}() - for c in con - _get_variables_in_constraints!(model, c, variables) - end - variable_constraints = MOI.ConstraintIndex[] - for (F, S) in MOI.get(model, MOI.ListOfConstraintTypesPresent()) - if !(F <: MOI.VariableIndex) - continue - end - for con in MOI.get(model, MOI.ListOfConstraintIndices{F,S}()) - if MOI.get(model, MOI.ConstraintFunction(), con) in variables - push!(variable_constraints, con) - end - end - end - return variable_constraints -end - -function _get_variables_in_constraints!( - model::MOI.ModelLike, - con::MOI.ConstraintIndex{<:MOI.ScalarAffineFunction}, - variables::Set{MOI.VariableIndex}, -) - f = MOI.get(model, MOI.ConstraintFunction(), con) - for term in f.terms - push!(variables, term.variable) - end - return -end - -function _get_variables_in_constraints!( - ::MOI.ModelLike, - ::MOI.ConstraintIndex, - ::Set{MOI.VariableIndex}, -) - return # skip -end - -function _feasibility_check(optimizer::Optimizer, original_model::MOI.ModelLike) - termination_status = MOI.get(original_model, MOI.TerminationStatus()) - if optimizer.verbose - println("Original model termination status: $(termination_status)") - end - if termination_status in - (MOI.OTHER_ERROR, MOI.INVALID_MODEL, MOI.OPTIMIZE_NOT_CALLED) - return false # because we can assert it is feasible - end - primal_status = MOI.get(original_model, MOI.PrimalStatus()) - if optimizer.verbose - println("Original model primal status: $(primal_status)") - end - if primal_status in (MOI.FEASIBLE_POINT, MOI.NEARLY_FEASIBLE_POINT) && !( - termination_status in - (MOI.INFEASIBLE, MOI.ALMOST_INFEASIBLE, MOI.LOCALLY_INFEASIBLE) - ) - return true - end - return false -end diff --git a/src/interval.jl b/src/interval.jl new file mode 100644 index 0000000..33f4210 --- /dev/null +++ b/src/interval.jl @@ -0,0 +1,30 @@ +# Copyright (c) 2025: Joaquim Dias Garcia, Oscar Dowson and contributors +# Copyright (c) 2014-2021: David P. Sanders & Luis Benet +# +# Use of this source code is governed by an MIT-style license that can be found +# in the LICENSE.md file or at https://opensource.org/licenses/MIT. + +# This type and the associated functions were inspired by IntervalArithmetic.jl + +struct Interval{T<:Real} + lo::T + hi::T + + function Interval(lo::T, hi::T) where {T<:Real} + @assert lo <= hi + return new{T}(lo, hi) + end +end + +Base.convert(::Type{Interval{T}}, x::T) where {T<:Real} = Interval(x, x) + +function Base.:+(a::Interval{T}, b::Interval{T}) where {T<:Real} + return Interval(a.lo + b.lo, a.hi + b.hi) +end + +function Base.:*(x::T, a::Interval{T}) where {T<:Real} + if x >= zero(T) + return Interval(a.lo * x, a.hi * x) + end + return Interval(a.hi * x, a.lo * x) +end diff --git a/src/range.jl b/src/range.jl deleted file mode 100644 index f6a309b..0000000 --- a/src/range.jl +++ /dev/null @@ -1,141 +0,0 @@ -# Copyright (c) 2025: Joaquim Dias Garcia, Oscar Dowson and contributors -# -# Use of this source code is governed by an MIT-style license that can be found -# in the LICENSE.md file or at https://opensource.org/licenses/MIT. - -# This type and the associated function were inspired by IntervalArithmetic.jl -# Copyright (c) 2014-2021: David P. Sanders & Luis Benet - -struct Interval{T<:Real} - lo::T - hi::T - - function Interval(lo::T, hi::T) where {T<:Real} - @assert lo <= hi - return new{T}(lo, hi) - end -end - -Base.convert(::Type{Interval{T}}, x::T) where {T<:Real} = Interval(x, x) - -Base.zero(::Type{Interval{T}}) where {T<:Real} = Interval(zero(T), zero(T)) - -Base.iszero(a::Interval) = iszero(a.hi) && iszero(a.lo) - -function Base.:+(a::Interval{T}, b::Interval{T}) where {T<:Real} - return Interval(a.lo + b.lo, a.hi + b.hi) -end - -function Base.:*(x::T, a::Interval{T}) where {T<:Real} - if iszero(a) || iszero(x) - return Interval(zero(T), zero(T)) - elseif x >= zero(T) - return Interval(a.lo * x, a.hi * x) - end - return Interval(a.hi * x, a.lo * x) -end - -# Back to functions written for MathOptIIS.jl - -_supports_range(::Type{MOI.ScalarAffineFunction{T}}) where {T} = true -_supports_range(::Type{MOI.EqualTo{T}}) where {T} = true -_supports_range(::Type{MOI.GreaterThan{T}}) where {T} = true -_supports_range(::Type{MOI.LessThan{T}}) where {T} = true -_supports_range(::Type{T}) where {T} = false - -function _range_infeasibility!( - optimizer::Optimizer, - ::Type{T}, - variables::Dict{MOI.VariableIndex,Interval{T}}, - lb_con::Dict{MOI.VariableIndex,MOI.ConstraintIndex}, - ub_con::Dict{MOI.VariableIndex,MOI.ConstraintIndex}, -) where {T} - range_consistent = true - for (F, S) in - MOI.get(optimizer.original_model, MOI.ListOfConstraintTypesPresent()) - if !_supports_range(F) || !_supports_range(S) - continue - end - range_consistent &= _range_infeasibility!( - optimizer, - optimizer.original_model, - T, - variables, - lb_con, - ub_con, - F, - S, - ) - end - return range_consistent -end - -function _range_infeasibility!( - optimizer::Optimizer, - original_model::MOI.ModelLike, - ::Type{T}, - variables::Dict{MOI.VariableIndex,Interval{T}}, - lb_con::Dict{MOI.VariableIndex,MOI.ConstraintIndex}, - ub_con::Dict{MOI.VariableIndex,MOI.ConstraintIndex}, - ::Type{F}, - ::Type{S}, -) where {T,F,S} - range_consistent = true - for con in MOI.get(original_model, MOI.ListOfConstraintIndices{F,S}()) - if !_in_time(optimizer) - return range_consistent - end - func = MOI.get(original_model, MOI.ConstraintFunction(), con) - interval = _eval_variables(variables, func) - if interval === nothing - continue - end - set = MOI.get(original_model, MOI.ConstraintSet(), con)::S - if !_invalid_range(set, interval) - continue - end - cons = Set{MOI.ConstraintIndex}() - push!(cons, con) - for t in func.terms - if (c = get(lb_con, t.variable, nothing)) !== nothing - push!(cons, c) - end - if (c = get(ub_con, t.variable, nothing)) !== nothing - push!(cons, c) - end - end - push!( - optimizer.results, - InfeasibilityData( - collect(cons), - true, # strictly speaking, we might need the proper "sides" - RangeData(interval.lo, interval.hi, set), - ), - ) - range_consistent = false - end - return range_consistent -end - -function _eval_variables( - map::AbstractDict{MOI.VariableIndex,U}, - f::MOI.ScalarAffineFunction, -) where {U} - out = convert(U, f.constant) - for t in f.terms - v = get(map, t.variable, nothing) - if v === nothing - return - end - out += t.coefficient * v - end - return out -end - -function _invalid_range(set::MOI.EqualTo, interval) - return !(interval.lo <= set.value <= interval.hi) -end - -_invalid_range(set::MOI.LessThan, interval) = set.upper < interval.lo - -_invalid_range(set::MOI.GreaterThan, interval) = interval.hi < set.lower diff --git a/src/solver.jl b/src/solver.jl deleted file mode 100644 index 952a11e..0000000 --- a/src/solver.jl +++ /dev/null @@ -1,229 +0,0 @@ -# Copyright (c) 2025: Joaquim Dias Garcia, Oscar Dowson and contributors -# -# Use of this source code is governed by an MIT-style license that can be found -# in the LICENSE.md file or at https://opensource.org/licenses/MIT. - -function _fix_slack(model, variable::MOI.VariableIndex, ::Type{T}) where {T} - lb_idx = MOI.ConstraintIndex{MOI.VariableIndex,MOI.GreaterThan{T}}( - variable.value, - ) - @assert MOI.is_valid(model, lb_idx) - MOI.delete(model, lb_idx) - MOI.add_constraint(model, variable, MOI.EqualTo{T}(zero(T))) - return -end - -function _unfix_slack(model, variable::MOI.VariableIndex, ::Type{T}) where {T} - eq_idx = - MOI.ConstraintIndex{MOI.VariableIndex,MOI.EqualTo{T}}(variable.value) - @assert MOI.is_valid(model, eq_idx) - MOI.delete(model, eq_idx) - MOI.add_constraint(model, variable, MOI.GreaterThan{T}(zero(T))) - return -end - -function _elastic_filter(optimizer::Optimizer) - T = Float64 - - model = MOI.instantiate(optimizer.optimizer) - MOI.set(model, MOI.Silent(), true) - for (k, v) in optimizer.optimizer_attributes - MOI.set(model, k, v) - end - reference_map = MOI.copy_to(model, optimizer.original_model) - - if optimizer.ignore_integrality - _cp_constraint_types = - MOI.get(model, MOI.ListOfConstraintTypesPresent()) - # _removed_constraints = Tuple[] - for (F, S) in _cp_constraint_types - con_list = [] - if S in ( - MOI.ZeroOne, - MOI.Integer, - MOI.Semicontinuous{T}, - MOI.Semiinteger{T}, - MOI.SOS1{T}, - MOI.SOS2{T}, - ) || S <: MOI.Indicator - con_list = MOI.get(model, MOI.ListOfConstraintIndices{F,S}()) - end - for con in con_list - MOI.delete(model, con) - # # save removed constraints - # func = MOI.get(model, MOI.ConstraintFunction(), con) - # set = MOI.get(model, MOI.ConstraintSet(), con) - # push!(_removed_constraints, (func, set)) - end - end - end - - obj_sense = MOI.get(model, MOI.ObjectiveSense()) - base_obj_type = MOI.get(model, MOI.ObjectiveFunctionType()) - base_obj_func = MOI.get(model, MOI.ObjectiveFunction{base_obj_type}()) - - constraint_to_affine = MOI.modify( - model, - MOI.Utilities.PenaltyRelaxation(; default = one(T), warn = false), - ) - # all slack variables added are of type ">= 0" - # might need to do something related to integers / binary - relaxed_obj_type = MOI.get(model, MOI.ObjectiveFunctionType()) - relaxed_obj_func = MOI.get(model, MOI.ObjectiveFunction{relaxed_obj_type}()) - pure_relaxed_obj_func = relaxed_obj_func - base_obj_func - - MOI.set( - model, - MOI.ObjectiveFunction{relaxed_obj_type}(), - pure_relaxed_obj_func, - ) - - max_iterations = length(constraint_to_affine) - - tolerance = optimizer.elastic_filter_tolerance - - de_elastisized = [] - - changed_obj = false - - # all (affine, non-bound) constraints are relaxed at this point - # we will try to set positive slacks to zero until the model infeasible - # the constraints of the fixed slacks are a IIS candidate - - for i in 1:max_iterations - if !_in_time(optimizer) - return nothing - end - MOI.optimize!(model) - status = MOI.get(model, MOI.TerminationStatus()) - if status in ( # possibily primal unbounded statuses - MOI.INFEASIBLE_OR_UNBOUNDED, - MOI.DUAL_INFEASIBLE, - MOI.ALMOST_DUAL_INFEASIBLE, - ) - # - end - if status in - (MOI.INFEASIBLE, MOI.ALMOST_INFEASIBLE, MOI.LOCALLY_INFEASIBLE) - break - end - for (con, func) in constraint_to_affine - @assert length(func.terms) <= 2 - if length(func.terms) == 1 - var = func.terms[1].variable - value = MOI.get(model, MOI.VariablePrimal(), var) - if value > tolerance - _fix_slack(model, var, T) - delete!(constraint_to_affine, con) - push!(de_elastisized, (con, var)) - end - elseif length(func.terms) == 2 - var1 = func.terms[1].variable - var2 = func.terms[2].variable - value1 = MOI.get(model, MOI.VariablePrimal(), var1) - value2 = MOI.get(model, MOI.VariablePrimal(), var2) - if value1 > tolerance && value2 > tolerance - error("IIS failed due numerical instability") - elseif value1 > tolerance - # TODO: coef is always 1.0 - _fix_slack(model, var1, T) - delete!(constraint_to_affine, con) - constraint_to_affine[con] = one(T) * var2 - push!(de_elastisized, (con, var1)) - elseif value2 > tolerance - _fix_slack(model, var2, T) - delete!(constraint_to_affine, con) - constraint_to_affine[con] = one(T) * var1 - push!(de_elastisized, (con, var2)) - end - end - end - end - - # consider deleting all no iis constraints - # be careful with intervals - - obj_type = MOI.get(model, MOI.ObjectiveFunctionType()) - obj_func = MOI.get(model, MOI.ObjectiveFunction{obj_type}()) - obj_sense = MOI.get(model, MOI.ObjectiveSense()) - - candidates = MOI.ConstraintIndex[] - if !optimizer.deletion_filter - for (con, var) in de_elastisized - push!(candidates, con) - end - empty!(de_elastisized) - end - - # deletion filter - for (con, var) in de_elastisized - if !_in_time(optimizer) - return nothing - end - _unfix_slack(model, var, T) - MOI.optimize!(model) - status = MOI.get(model, MOI.TerminationStatus()) - if status in - (MOI.INFEASIBLE, MOI.ALMOST_INFEASIBLE, MOI.LOCALLY_INFEASIBLE) - # this constraint is not in IIS - # hence it remains unfixed - elseif status in ( - MOI.OPTIMAL, - MOI.ALMOST_OPTIMAL, - MOI.LOCALLY_SOLVED, - MOI.ALMOST_LOCALLY_SOLVED, - ) - push!(candidates, con) - _fix_slack(model, var, T) - elseif status in ( # possibily primal unbounded statuses - MOI.INFEASIBLE_OR_UNBOUNDED, - MOI.DUAL_INFEASIBLE, - MOI.ALMOST_DUAL_INFEASIBLE, - ) - MOI.set(model, MOI.ObjectiveSense(), MOI.FEASIBILITY_SENSE) - MOI.optimize!(model) - primal_status = MOI.get(model, MOI.PrimalStatus()) - # the unbounded case: - if primal_status in (MOI.FEASIBLE_POINT, MOI.NEARLY_FEASIBLE_POINT) - # this constraint is not in IIS - push!(candidates, con) - _fix_slack(model, var, T) - MOI.set(model, MOI.ObjectiveSense(), obj_sense) - MOI.set( - model, - MOI.ObjectiveFunction{relaxed_obj_type}(), - pure_relaxed_obj_func, - ) - # the both primal and dual infeasible case: - # else - # nothing - end - else - error("IIS failed due numerical instability, got status $status") - end - end - - if isempty(candidates) - return nothing - end - - pre_iis = Set(candidates) - iis = MOI.ConstraintIndex[] - for (F, S) in - MOI.get(optimizer.original_model, MOI.ListOfConstraintTypesPresent()) - if F == MOI.VariableIndex - continue - end - for con in MOI.get( - optimizer.original_model, - MOI.ListOfConstraintIndices{F,S}(), - ) - new_con = reference_map[con] - if new_con in pre_iis - push!(iis, con) - end - end - end - - return iis -end diff --git a/test/runtests.jl b/test/runtests.jl index df156e1..f88e1c0 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -206,36 +206,6 @@ function test_range_and_bound() data[].constraints, [index(LowerBoundRef(z)), index(UpperBoundRef(z))], ) - @test MOI.get(solver, MOIIS.StopIfInfeasibleBounds()) == true - MOI.set(solver, MOIIS.StopIfInfeasibleBounds(), false) - @test MOI.get(solver, MOIIS.StopIfInfeasibleBounds()) == false - MOI.compute_conflict!(solver) - data = solver.results - @test length(data) == 2 - @test _isequal_unordered( - data[1].constraints, - [index(LowerBoundRef(z)), index(UpperBoundRef(z))], - ) - @test _isequal_unordered( - data[2].constraints, - [ - index(c), - index(UpperBoundRef(x)), - index(LowerBoundRef(x)), - index(UpperBoundRef(y)), - index(LowerBoundRef(y)), - ], - ) - @test data[2].irreducible - @test data[2].metadata == - MOIIS.RangeData(11.0, 22.0, MOI.LessThan{Float64}(1.0)) - @test MOI.get(solver, MOIIS.StopIfInfeasibleRanges()) == true - MOI.set(solver, MOIIS.StopIfInfeasibleRanges(), false) - @test MOI.get(solver, MOIIS.StopIfInfeasibleRanges()) == false - MOI.set(solver, MOIIS.InnerOptimizer(), HiGHS.Optimizer) - MOI.compute_conflict!(solver) - data = solver.results - @test length(data) == 2 return end @@ -255,17 +225,6 @@ function test_range_and_bound_2() data[].constraints, [index(LowerBoundRef(z)), index(UpperBoundRef(z))], ) - @test MOI.get(solver, MOIIS.StopIfInfeasibleBounds()) == true - MOI.set(solver, MOIIS.StopIfInfeasibleBounds(), false) - @test MOI.get(solver, MOIIS.StopIfInfeasibleBounds()) == false - MOI.compute_conflict!(solver) - data = solver.results - # the result is only one conflic again because the range fail cant be computed - @test length(data) == 1 - @test _isequal_unordered( - data[1].constraints, - [index(LowerBoundRef(z)), index(UpperBoundRef(z))], - ) return end @@ -407,13 +366,6 @@ function test_interval() MOI.compute_conflict!(solver) data = solver.results @test length(data) == 0 - @test !MOI.get(solver, MOIIS.SkipFeasibilityCheck()) - MOI.set(solver, MOIIS.SkipFeasibilityCheck(), true) - @test MOI.get(solver, MOIIS.SkipFeasibilityCheck()) - MOI.compute_conflict!(solver) - data = solver.results - @test length(data) == 0 - # TODO check status return end @@ -433,15 +385,6 @@ function test_pass_attribute() @test MOI.get(solver, MOI.Silent()) == true MOI.set(solver, MOI.Silent(), false) @test MOI.get(solver, MOI.Silent()) == false - @test MOI.get(solver, MOIIS.ElasticFilterTolerance()) == 1e-5 - MOI.set(solver, MOIIS.ElasticFilterTolerance(), 1e-3) - @test MOI.get(solver, MOIIS.ElasticFilterTolerance()) == 1e-3 - MOI.compute_conflict!(solver) - data = solver.results - @test length(data) == 0 - MOI.set(solver, MOIIS.SkipFeasibilityCheck(), true) - @test MOI.get(solver, MOIIS.SkipFeasibilityCheck()) - MOI.compute_conflict!(solver) data = solver.results @test length(data) == 0 return @@ -492,64 +435,6 @@ function test_iis() return end -function test_iis_no_deletion_filter() - model = Model(HiGHS.Optimizer) - set_silent(model) - @variable(model, 0 <= x <= 10) - @variable(model, 0 <= y <= 20) - @constraint(model, c1, x + y <= 1) - @constraint(model, c2, x + y >= 2) - @objective(model, Max, x + y) - optimize!(model) - solver = MOIIS.Optimizer() - MOI.set(solver, MOIIS.InfeasibleModel(), backend(model)) - MOI.compute_conflict!(solver) - data = solver.results - @test length(data) == 0 - solver = MOIIS.Optimizer() - MOI.set(solver, MOIIS.InfeasibleModel(), backend(model)) - MOI.set(solver, MOIIS.InnerOptimizer(), HiGHS.Optimizer) - @test MOI.get(solver, MOIIS.DeletionFilter()) == true - MOI.set(solver, MOIIS.DeletionFilter(), false) - @test MOI.get(solver, MOIIS.DeletionFilter()) == false - MOI.compute_conflict!(solver) - data = solver.results - @test length(data) == 1 - @test data[].irreducible - @test data[].metadata == MOIIS.NoData() - @test _isequal_unordered(data[].constraints, [index(c2), index(c1)]) - return -end - -function test_iis_ignore_integrality() - model = Model(HiGHS.Optimizer) - set_silent(model) - @variable(model, 0 <= x <= 10) - @variable(model, 0 <= y <= 20, Bin) - @constraint(model, c1, x + y <= 1) - @constraint(model, c2, x + y >= 2) - @objective(model, Max, x + y) - optimize!(model) - solver = MOIIS.Optimizer() - MOI.set(solver, MOIIS.InfeasibleModel(), backend(model)) - MOI.compute_conflict!(solver) - data = solver.results - @test length(data) == 0 - solver = MOIIS.Optimizer() - MOI.set(solver, MOIIS.InfeasibleModel(), backend(model)) - MOI.set(solver, MOIIS.InnerOptimizer(), HiGHS.Optimizer) - @test MOI.get(solver, MOIIS.ElasticFilterIgnoreIntegrality()) == false - MOI.set(solver, MOIIS.ElasticFilterIgnoreIntegrality(), true) - @test MOI.get(solver, MOIIS.ElasticFilterIgnoreIntegrality()) == true - MOI.compute_conflict!(solver) - data = solver.results - @test length(data) == 1 - @test data[].irreducible - @test data[].metadata == MOIIS.NoData() - @test _isequal_unordered(data[].constraints, [index(c2), index(c1)]) - return -end - function test_pass_attribute_inner() model = Model(HiGHS.Optimizer) set_silent(model) @@ -561,8 +446,14 @@ function test_pass_attribute_inner() optimize!(model) solver = MOIIS.Optimizer() MOI.set(solver, MOIIS.InfeasibleModel(), backend(model)) - MOI.set(solver, MOIIS.InnerOptimizer(), HiGHS.Optimizer) - MOI.set(solver, MOIIS.InnerOptimizerAttribute(MOI.TimeLimitSec()), 10.0) + MOI.set( + solver, + MOIIS.InnerOptimizer(), + MOI.OptimizerWithAttributes( + HiGHS.Optimizer, + MOI.TimeLimitSec() => 10.0, + ), + ) MOI.compute_conflict!(solver) data = solver.results @test length(data) == 1 @@ -693,22 +584,22 @@ function test_iis_spare() solver, MOI.ConstraintConflictStatus(), index(LowerBoundRef(x)), - ) == MOI.MAYBE_IN_CONFLICT + ) == MOI.NOT_IN_CONFLICT @test MOI.get( solver, MOI.ConstraintConflictStatus(), index(LowerBoundRef(y)), - ) == MOI.MAYBE_IN_CONFLICT + ) == MOI.NOT_IN_CONFLICT @test MOI.get( solver, MOI.ConstraintConflictStatus(), index(UpperBoundRef(x)), - ) == MOI.MAYBE_IN_CONFLICT + ) == MOI.NOT_IN_CONFLICT @test MOI.get( solver, MOI.ConstraintConflictStatus(), index(UpperBoundRef(y)), - ) == MOI.MAYBE_IN_CONFLICT + ) == MOI.NOT_IN_CONFLICT @test MOI.get( solver, MOI.ConstraintConflictStatus(), @@ -728,8 +619,6 @@ function test_iis_binary() @variable(model, x, Bin) @constraint(model, c1, x == 1 / 2) optimize!(model) - @show termination_status(model) - @show primal_status(model) solver = MOIIS.Optimizer() MOI.set(solver, MOIIS.InfeasibleModel(), backend(model)) MOI.set(solver, MOIIS.InnerOptimizer(), HiGHS.Optimizer) @@ -749,10 +638,245 @@ function test_iis_binary() return end -function test_deprecated() - @test (@test_deprecated MOIIS.ConflictCount()) == MOI.ConflictCount() - @test (@test_deprecated MOIIS.ConstraintConflictStatus(2)) == - MOI.ConstraintConflictStatus(2) +function test_verbose() + model = Model(HiGHS.Optimizer) + @variable(model, 0 <= x <= 10) + @variable(model, 0 <= y) + @constraint(model, 2 * y <= 40) + @constraint(model, c, x + y >= 35) + @objective(model, Max, x + y) + optimize!(model) + solver = MOIIS.Optimizer() + MOI.set(solver, MOI.Silent(), false) + MOI.set(solver, MOIIS.InfeasibleModel(), backend(model)) + MOI.set(solver, MOIIS.InnerOptimizer(), HiGHS.Optimizer) + dir = mktempdir() + open(joinpath(dir, "log.stdout"), "w") do io + return redirect_stdout(io) do + return MOI.compute_conflict!(solver) + end + end + contents = """ + [MathOptIIS] starting compute_conflict! + [MathOptIIS] model termination status: INFEASIBLE + [MathOptIIS] model primal status: NO_SOLUTION + [MathOptIIS] starting bound analysis + [MathOptIIS] bound analysis found 0 infeasible subsets + [MathOptIIS] starting range analysis + [MathOptIIS] range analysis found 0 infeasible subsets + [MathOptIIS] starting elastic filter + [MathOptIIS] using INFEASIBILITY_CERTIFICATE to construct infeasible subset + [MathOptIIS] elastic filter found 1 infeasible subsets + """ + @test read(joinpath(dir, "log.stdout"), String) == contents + return +end + +function _copy_conflict(model::MOI.ModelLike, solver::MOIIS.Optimizer) + filter_fn(::Any) = true + function filter_fn(cref::MOI.ConstraintIndex) + status = MOI.get(solver, MOI.ConstraintConflictStatus(), cref) + return status != MOI.NOT_IN_CONFLICT + end + new_model = MOI.Utilities.Model{Float64}() + filtered_src = MOI.Utilities.ModelFilter(filter_fn, model) + MOI.copy_to(new_model, filtered_src) + MOI.set(new_model, MOI.ObjectiveSense(), MOI.FEASIBILITY_SENSE) + return new_model +end + +function _test_compute_conflict(input, output) + model = HiGHS.Optimizer() + MOI.set(model, MOI.Silent(), true) + MOI.Utilities.loadfromstring!(model, input) + MOI.optimize!(model) + solver = MOIIS.Optimizer() + MOI.set(solver, MOIIS.InfeasibleModel(), model) + MOI.set(solver, MOIIS.InnerOptimizer(), HiGHS.Optimizer) + MOI.compute_conflict!(solver) + iis = _copy_conflict(model, solver) + target = MOI.Utilities.Model{Float64}() + MOI.Utilities.loadfromstring!(target, output) + A, B = sprint(print, iis), sprint(print, target) + if A != B + @info "IIS" + print(A) + @info "Target" + println(B) + end + @test A == B + return +end + +function test_relax_integrality_integer() + _test_compute_conflict( + """ + variables: x, y + maxobjective: 1.0 * x + y + 2.0 * y <= 40.0 + 1.0 * x + y >= 35.0 + x >= 0.0 + x <= 10.0 + x in Integer() + y >= 0.0 + """, + """ + variables: x, y + 2.0 * y <= 40.0 + 1.0 * x + y >= 35.0 + x <= 10.0 + """, + ) + return +end + +function test_relax_integrality_integer_continuous_infeasible() + _test_compute_conflict( + """ + variables: x + maxobjective: 1.0 * x + x in Integer() + 1.1 * x == 1.0 + """, + """ + variables: x + x in Integer() + 1.1 * x == 1.0 + """, + ) + return +end + +function test_relax_integrality_zero_one() + _test_compute_conflict( + """ + variables: x, y + x in ZeroOne() + x == 1.0 + 1.0 * y <= 0.5 + 1.0 * x + -1.0 * y == 0.0 + """, + """ + variables: x, y + x == 1.0 + 1.0 * y <= 0.5 + 1.0 * x + -1.0 * y == 0.0 + """, + ) + _test_compute_conflict( + """ + variables: x, y + x in ZeroOne() + x in Interval(0.3, 1.2) + 1.0 * y <= 0.5 + 1.0 * x + -1.0 * y == 0.0 + """, + """ + variables: x, y + x in ZeroOne() + x in Interval(0.3, 1.2) + 1.0 * y <= 0.5 + 1.0 * x + -1.0 * y == 0.0 + """, + ) + _test_compute_conflict( + """ + variables: x, y + x in ZeroOne() + x in Interval(0.8, 1.2) + 1.0 * y <= 0.5 + 1.0 * x + -1.0 * y == 0.0 + """, + """ + variables: x, y + x in Interval(0.8, 1.2) + 1.0 * y <= 0.5 + 1.0 * x + -1.0 * y == 0.0 + """, + ) + _test_compute_conflict( + """ + variables: x, y + x in ZeroOne() + x >= 0.8 + 1.0 * y <= 0.5 + 1.0 * x + -1.0 * y == 0.0 + """, + """ + variables: x, y + x >= 0.8 + 1.0 * y <= 0.5 + 1.0 * x + -1.0 * y == 0.0 + """, + ) + _test_compute_conflict( + """ + variables: x, y + x in ZeroOne() + x <= 0.4 + 1.0 * y >= 0.5 + 1.0 * x + -1.0 * y == 0.0 + """, + """ + variables: x, y + x <= 0.4 + 1.0 * y >= 0.5 + 1.0 * x + -1.0 * y == 0.0 + """, + ) + _test_compute_conflict( + """ + variables: x, y + x in ZeroOne() + 1.0 * y >= 1.5 + 1.0 * x + -1.0 * y == 0.0 + """, + """ + variables: x, y + x in ZeroOne() + 1.0 * y >= 1.5 + 1.0 * x + -1.0 * y == 0.0 + """, + ) + _test_compute_conflict( + """ + variables: x, y + x in ZeroOne() + 1.0 * y <= -0.1 + 1.0 * x + -1.0 * y == 0.0 + """, + """ + variables: x, y + x in ZeroOne() + 1.0 * y <= -0.1 + 1.0 * x + -1.0 * y == 0.0 + """, + ) + return +end + +function test_unsupported_termination_status() + model = Model(HiGHS.Optimizer) + @variable(model, 0 <= x <= 10) + @variable(model, 0 <= y) + @constraint(model, 2 * y <= 40) + @constraint(model, c, x + y >= 35) + @objective(model, Max, x + y) + optimize!(model) + solver = MOIIS.Optimizer() + MOI.set(solver, MOIIS.InfeasibleModel(), backend(model)) + opt = MOI.OptimizerWithAttributes( + HiGHS.Optimizer, + "simplex_iteration_limit" => 0, + ) + MOI.set(solver, MOIIS.InnerOptimizer(), opt) + status = MOI.ITERATION_LIMIT + @test_throws( + ErrorException( + "Unsupported termination status in IIS algorithm: $status", + ), + MOI.compute_conflict!(solver), + ) return end