diff --git a/Project.toml b/Project.toml index b1d640a..29b0e2c 100644 --- a/Project.toml +++ b/Project.toml @@ -5,6 +5,7 @@ version = "0.1.0" [deps] Dualization = "191a621a-6537-11e9-281d-650236a99e60" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +MathOptIIS = "8c4f8055-bd93-4160-a86b-a0c04941dbff" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" @@ -17,5 +18,6 @@ MathOptAnalyzerJuMPExt = "JuMP" [compat] Dualization = "0.6.0" JuMP = "1.24.0" +MathOptIIS = "0.1.0" MathOptInterface = "1.37.0" julia = "1.10" diff --git a/src/Infeasibility/Infeasibility.jl b/src/Infeasibility/Infeasibility.jl index 9034d40..a1870df 100644 --- a/src/Infeasibility/Infeasibility.jl +++ b/src/Infeasibility/Infeasibility.jl @@ -8,9 +8,6 @@ module Infeasibility import MathOptInterface as MOI import MathOptAnalyzer -include("intervals.jl") -include("_eval_variables.jl") - include("structs.jl") include("analyze.jl") include("summarize.jl") diff --git a/src/Infeasibility/_eval_variables.jl b/src/Infeasibility/_eval_variables.jl deleted file mode 100644 index de0bfab..0000000 --- a/src/Infeasibility/_eval_variables.jl +++ /dev/null @@ -1,24 +0,0 @@ -# Copyright (c) 2025: Joaquim 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 _eval_variables end - -function _eval_variables(value_fn::Function, t::MOI.ScalarAffineTerm) - return t.coefficient * value_fn(t.variable) -end - -function _eval_variables( - value_fn::Function, - f::MOI.ScalarAffineFunction{T}, -) where {T} - # TODO: this conversion exists in JuMP, but not in MOI - S = Base.promote_op(value_fn, MOI.VariableIndex) - U = MOI.MA.promote_operation(*, T, S) - out = convert(U, f.constant) - for t in f.terms - out += _eval_variables(value_fn, t) - end - return out -end diff --git a/src/Infeasibility/analyze.jl b/src/Infeasibility/analyze.jl index 4262cb2..ce00656 100644 --- a/src/Infeasibility/analyze.jl +++ b/src/Infeasibility/analyze.jl @@ -3,449 +3,71 @@ # 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. +import MathOptIIS as MOIIS + function MathOptAnalyzer.analyze( ::Analyzer, model::MOI.ModelLike; optimizer = nothing, ) out = Data() - T = Float64 - variables = Dict{MOI.VariableIndex,Interval{T}}() - - variable_indices = MOI.get(model, MOI.ListOfVariableIndices()) - - lb = Dict{MOI.VariableIndex,T}() - ub = Dict{MOI.VariableIndex,T}() - - for con in MOI.get( - model, - MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.EqualTo{T}}(), - ) - set = MOI.get(model, MOI.ConstraintSet(), con) - func = MOI.get(model, MOI.ConstraintFunction(), con) - lb[func] = set.value - ub[func] = set.value - end + solver = MOIIS.Optimizer() + MOI.set(solver, MOIIS.InfeasibleModel(), model) - for con in MOI.get( - model, - MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.LessThan{T}}(), - ) - set = MOI.get(model, MOI.ConstraintSet(), con) - func = MOI.get(model, MOI.ConstraintFunction(), con) - # lb[func] = -Inf - ub[func] = set.upper + if optimizer !== nothing + MOI.set(solver, MOIIS.InnerOptimizer(), optimizer) end - for con in MOI.get( - model, - MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.GreaterThan{T}}(), - ) - set = MOI.get(model, MOI.ConstraintSet(), con) - func = MOI.get(model, MOI.ConstraintFunction(), con) - lb[func] = set.lower - # ub[func] = Inf - end - - for con in MOI.get( - model, - MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.Interval{T}}(), - ) - set = MOI.get(model, MOI.ConstraintSet(), con) - func = MOI.get(model, MOI.ConstraintFunction(), con) - lb[func] = set.lower - ub[func] = set.upper - end - - # for con in MOI.get(model, MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.SemiContinuous{T}}()) - # set = MOI.get(model, MOI.ConstraintSet(), con) - # func = MOI.get(model, MOI.ConstraintFunction(), con) - # lb[func] = 0 # set.lower - # ub[func] = set.upper - # end - - # for con in MOI.get(model, MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.SemiInteger{T}}()) - # set = MOI.get(model, MOI.ConstraintSet(), con) - # func = MOI.get(model, MOI.ConstraintFunction(), con) - # lb[func] = 0 #set.lower - # ub[func] = set.upper - # end + MOI.compute_conflict!(solver) - bounds_consistent = true + data = solver.results - for con in MOI.get( - model, - MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.Integer}(), - ) - func = MOI.get(model, MOI.ConstraintFunction(), con) - _lb = get(lb, func, -Inf) - _ub = get(ub, func, Inf) - if abs(_ub - _lb) < 1 && ceil(_ub) == ceil(_lb) + for iis in data + meta = iis.metadata + if typeof(meta) <: MOIIS.BoundsData + constraints = iis.constraints + @assert length(constraints) == 2 + func = MOI.get(model, MOI.ConstraintFunction(), constraints[1]) push!( - out.infeasible_integrality, - InfeasibleIntegrality(func, _lb, _ub, MOI.Integer()), + out.infeasible_bounds, + InfeasibleBounds{T}(func, meta.lower_bound, meta.upper_bound), ) - bounds_consistent = false - end - end - - for con in MOI.get( - model, - MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.ZeroOne}(), - ) - func = MOI.get(model, MOI.ConstraintFunction(), con) - _lb = get(lb, func, -Inf) - _ub = get(ub, func, Inf) - if _lb > 0 && _ub < 1 + elseif typeof(meta) <: MOIIS.IntegralityData + constraints = iis.constraints + @assert length(constraints) >= 2 + func = MOI.get(model, MOI.ConstraintFunction(), constraints[1]) push!( out.infeasible_integrality, - InfeasibleIntegrality(func, _lb, _ub, MOI.ZeroOne()), - ) - bounds_consistent = false - end - end - - for var in variable_indices - _lb = get(lb, var, -Inf) - _ub = get(ub, var, Inf) - if _lb > _ub - push!(out.infeasible_bounds, InfeasibleBounds(var, _lb, _ub)) - bounds_consistent = false - else - variables[var] = Interval(_lb, _ub) - end - end - - # check PSD diagonal >= 0 ? - # other cones? - if !bounds_consistent - return out - end - - # second layer of infeasibility analysis is constraint range analysis - range_consistent = true - - for con in MOI.get( - model, - MOI.ListOfConstraintIndices{ - MOI.ScalarAffineFunction{T}, - MOI.EqualTo{T}, - }(), - ) - set = MOI.get(model, MOI.ConstraintSet(), con) - func = MOI.get(model, MOI.ConstraintFunction(), con) - failed = false - interval = _eval_variables(func) do var_idx - # this only fails if we allow continuing after bounds issues - # if !haskey(variables, var_idx) - # failed = true - # return Interval(-Inf, Inf) - # end - return variables[var_idx] - end - # if failed - # continue - # end - rhs = set.value - if interval.lo > rhs || interval.hi < rhs - push!( - out.constraint_range, - InfeasibleConstraintRange(con, interval.lo, interval.hi, set), - ) - range_consistent = false - end - end - - for con in MOI.get( - model, - MOI.ListOfConstraintIndices{ - MOI.ScalarAffineFunction{T}, - MOI.LessThan{T}, - }(), - ) - set = MOI.get(model, MOI.ConstraintSet(), con) - func = MOI.get(model, MOI.ConstraintFunction(), con) - failed = false - interval = _eval_variables(func) do var_idx - # this only fails if we allow continuing after bounds issues - # if !haskey(variables, var_idx) - # failed = true - # return Interval(-Inf, Inf) - # end - return variables[var_idx] - end - # if failed - # continue - # end - rhs = set.upper - if interval.lo > rhs - push!( - out.constraint_range, - InfeasibleConstraintRange(con, interval.lo, interval.hi, set), - ) - range_consistent = false - end - end - - for con in MOI.get( - model, - MOI.ListOfConstraintIndices{ - MOI.ScalarAffineFunction{T}, - MOI.GreaterThan{T}, - }(), - ) - set = MOI.get(model, MOI.ConstraintSet(), con) - func = MOI.get(model, MOI.ConstraintFunction(), con) - failed = false - interval = _eval_variables(func) do var_idx - # this only fails if we allow continuing after bounds issues - # if !haskey(variables, var_idx) - # failed = true - # return Interval(-Inf, Inf) - # end - return variables[var_idx] - end - # if failed - # continue - # end - rhs = set.lower - if interval.hi < rhs - push!( - out.constraint_range, - InfeasibleConstraintRange(con, interval.lo, interval.hi, set), - ) - range_consistent = false - end - end - - if !range_consistent - return out - end - - # check if there is a optimizer - # third layer is an IIS resolver - if optimizer === nothing - println("iis resolver cannot continue because no optimizer is provided") - return out - end - iis = iis_elastic_filter(model, optimizer) - # for now, only one iis is computed - if iis !== nothing - push!(out.iis, IrreducibleInfeasibleSubset(iis)) - end - - return out -end - -function _fix_to_zero(model, variable::MOI.VariableIndex, ::Type{T}) where {T} - ub_idx = - MOI.ConstraintIndex{MOI.VariableIndex,MOI.LessThan{T}}(variable.value) - lb_idx = MOI.ConstraintIndex{MOI.VariableIndex,MOI.GreaterThan{T}}( - variable.value, - ) - has_lower = false - if MOI.is_valid(model, lb_idx) - MOI.delete(model, lb_idx) - has_lower = true - # MOI.PenaltyRelaxation only creates variables with LB - # elseif MOI.is_valid(model, ub_idx) - # MOI.delete(model, ub_idx) - else - error("Variable is not bounded") - end - MOI.add_constraint(model, variable, MOI.EqualTo{T}(zero(T))) - return has_lower -end - -function _set_bound_zero( - model, - variable::MOI.VariableIndex, - has_lower::Bool, - ::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) - if has_lower - MOI.add_constraint(model, variable, MOI.GreaterThan{T}(zero(T))) - # MOI.PenaltyRelaxation only creates variables with LB - # else - # MOI.add_constraint(model, variable, MOI.LessThan{T}(zero(T))) - end - return -end - -function iis_elastic_filter(original_model::MOI.ModelLike, optimizer) - T = Float64 - - # handle optimize not called - status = MOI.get(original_model, MOI.TerminationStatus()) - if !( - status in - (MOI.INFEASIBLE, MOI.ALMOST_INFEASIBLE, MOI.ALMOST_INFEASIBLE) - ) - println( - "iis resolver cannot continue because model is found to be $(status) by the solver", - ) - return nothing - end - - model = MOI.instantiate(optimizer) - reference_map = MOI.copy_to(model, original_model) - MOI.set(model, MOI.Silent(), true) - - 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 = 1.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 - - max_iterations = length(constraint_to_affine) - - tolerance = 1e-5 - - de_elastisized = [] - - changed_obj = false - - for i in 1:max_iterations - MOI.optimize!(model) - status = MOI.get(model, MOI.TerminationStatus()) - if status in ( # possibily primal unbounded - MOI.INFEASIBLE_OR_UNBOUNDED, - MOI.DUAL_INFEASIBLE, - MOI.ALMOST_DUAL_INFEASIBLE, - ) - #try with a pure relaxation objective - MOI.set( - model, - MOI.ObjectiveFunction{relaxed_obj_type}(), - pure_relaxed_obj_func, + InfeasibleIntegrality{T}( + func, + meta.lower_bound, + meta.upper_bound, + meta.set, + ), ) - changed_obj = true - MOI.optimize!(model) - end - if status in - (MOI.INFEASIBLE, MOI.ALMOST_INFEASIBLE, MOI.ALMOST_INFEASIBLE) - break - end - for (con, func) in constraint_to_affine - if length(func.terms) == 1 - var = func.terms[1].variable - value = MOI.get(model, MOI.VariablePrimal(), var) - if value > tolerance - has_lower = _fix_to_zero(model, var, T) - delete!(constraint_to_affine, con) - push!(de_elastisized, (con, var, has_lower)) + elseif typeof(meta) <: MOIIS.RangeData + constraints = iis.constraints + @assert length(constraints) >= 1 + # main_con = nothing + for con in constraints + if !(typeof(con) <: MOI.ConstraintIndex{MOI.VariableIndex}) + push!( + out.constraint_range, + InfeasibleConstraintRange{T}( + con, + meta.lower_bound, + meta.upper_bound, + meta.set, + ), + ) + break end - elseif length(func.terms) == 2 - var1 = func.terms[1].variable - coef1 = func.terms[1].coefficient - var2 = func.terms[2].variable - coef2 = func.terms[2].coefficient - 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 alwayas 1.0 - has_lower = _fix_to_zero(model, var1, T) - delete!(constraint_to_affine, con) - constraint_to_affine[con] = coef2 * var2 - push!(de_elastisized, (con, var1, has_lower)) - elseif value2 > tolerance - has_lower = _fix_to_zero(model, var2, T) - delete!(constraint_to_affine, con) - constraint_to_affine[con] = coef1 * var1 - push!(de_elastisized, (con, var2, has_lower)) - end - else - println( - "$con and relaxing function with more than two terms: $func", - ) - end - end - end - - if changed_obj - MOI.set( - model, - MOI.ObjectiveFunction{relaxed_obj_type}(), - relaxed_obj_func, - ) - 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()) - - # deletion filter - cadidates = MOI.ConstraintIndex[] - for (con, var, has_lower) in de_elastisized - _set_bound_zero(model, var, has_lower, T) - MOI.optimize!(model) - status = MOI.get(model, MOI.TerminationStatus()) - if status in - (MOI.INFEASIBLE, MOI.ALMOST_INFEASIBLE, MOI.ALMOST_INFEASIBLE) - # this constraint is not in IIS - elseif status in ( - MOI.OPTIMAL, - MOI.ALMOST_OPTIMAL, - MOI.LOCALLY_SOLVED, - MOI.ALMOST_LOCALLY_SOLVED, - ) - push!(cadidates, con) - _fix_to_zero(model, var, T) - elseif status in ( - MOI.INFEASIBLE_OR_UNBOUNDED, - MOI.DUAL_INFEASIBLE, - MOI.ALMOST_DUAL_INFEASIBLE, # possibily primal unbounded - ) - MOI.set(model, MOI.ObjectiveSense(), MOI.FEASIBILITY_SENSE) - MOI.optimize!(model) - primal_status = MOI.get(model, MOI.PrimalStatus()) - if primal_status in (MOI.FEASIBLE_POINT, MOI.NEARLY_FEASIBLE_POINT) - # this constraint is not in IIS - push!(cadidates, con) - _fix_to_zero(model, var, T) - MOI.set(model, MOI.ObjectiveSense(), obj_sense) - MOI.set(model, MOI.ObjectiveFunction{obj_type}(), obj_func) - else - error( - "IIS failed due numerical instability, got status $status,", - "then, for MOI.FEASIBILITY_SENSE objective, got primal status $primal_status", - ) end else - error("IIS failed due numerical instability, got status $status") + push!(out.iis, IrreducibleInfeasibleSubset(iis.constraints)) end end - - pre_iis = Set(cadidates) - iis = MOI.ConstraintIndex[] - for (F, S) in MOI.get(original_model, MOI.ListOfConstraintTypesPresent()) - if F == MOI.VariableIndex - continue - end - for con in MOI.get(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 + return out end diff --git a/src/Infeasibility/intervals.jl b/src/Infeasibility/intervals.jl deleted file mode 100644 index f54eb84..0000000 --- a/src/Infeasibility/intervals.jl +++ /dev/null @@ -1,112 +0,0 @@ -# Copyright (c) 2025: Joaquim 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 file was originally written by Joaquim Garcia and contributors in -# BilevelJuMP.jl, which is licensed under the MIT "Expat" License: - -# The BilevelJuMP.jl package is licensed under the MIT "Expat" License: - -# Copyright (c) 2019 Joaquim Dias Garcia, and contributors - -# Permission is hereby granted, free of charge, to any person obtaining a copy -# of this software and associated documentation files (the "Software"), to deal -# in the Software without restriction, including without limitation the rights -# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -# copies of the Software, and to permit persons to whom the Software is -# furnished to do so, subject to the following conditions: - -# The above copyright notice and this permission notice shall be included in all -# copies or substantial portions of the Software. - -# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -# SOFTWARE. - -# Function in this file are heavily inspired in IntervalArithmetic.jl, -# which is licensed under the MIT "Expat" License: -# -# Copyright (c) 2014-2021: David P. Sanders & Luis Benet -# Permission is hereby granted, free of charge, to any person obtaining a copy -# of this software and associated documentation files (the "Software"), to deal -# in the Software without restriction, including without limitation the rights -# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -# copies of the Software, and to permit persons to whom the Software is -# furnished to do so, subject to the following conditions: -# -# The above copyright notice and this permission notice shall be included in all -# copies or substantial portions of the Software. -# -# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -# SOFTWARE. - -struct Interval{T} - lo::T - hi::T -end - -function Interval(lo::T, hi::T) where {T<:Real} - # if hi < lo <= hi + eps(T) - # lo = hi - # end - @assert lo <= hi - return Interval{T}(lo, hi) -end - -function Base.zero(::Type{Interval{T}}) where {T<:Real} - return Interval(zero(T), zero(T)) -end - -function Base.iszero(a::Interval) - return iszero(a.hi) && iszero(a.lo) -end - -# this code is only used for interval += scalar_coef * interval -# so only bivariate + and * are needes - -# Base.:+(a::Interval) = a -# Base.:-(a::Interval) = Interval(-a.hi, -a.lo) - -# function Base.:+(a::Interval{T}, b::T) where {T<:Real} -# return Interval(a.lo + b, a.hi + b) -# end -# Base.:+(b::T, a::Interval{T}) where {T<:Real} = a + b - -# function Base.:-(a::Interval{T}, b::T) where {T<:Real} -# return Interval(a.lo - b, a.hi - b) -# end -# function Base.:-(b::T, a::Interval{T}) where {T<:Real} -# return Interval(b - a.hi, b - a.lo) -# end - -function Base.:+(a::Interval{T}, b::Interval{T}) where {T<:Real} - return Interval(a.lo + b.lo, a.hi + b.hi) -end - -# function Base.:-(a::Interval{T}, b::Interval{T}) where {T<:Real} -# return Interval(a.lo - b.hi, a.hi - b.lo) -# end - -## Multiplication -function Base.:*(x::T, a::Interval{T}) where {T<:Real} - (iszero(a) || iszero(x)) && return Interval(zero(T), zero(T)) - if x ≥ zero(T) - return Interval(a.lo * x, a.hi * x) - else - return Interval(a.hi * x, a.lo * x) - end -end - -# Base.:*(a::Interval{T}, x::T) where {T<:Real} = x * a - -Base.convert(::Type{Interval{T}}, x::T) where {T<:Real} = Interval(x, x)