From 4018990026f31acdfec74b2e22d1dca326c18288 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Wed, 30 Apr 2025 20:15:35 -0300 Subject: [PATCH 01/34] Move code to MOI and add JuMP layer --- src/ModelAnalyzer.jl | 1 + src/jump.jl | 29 +++ src/numerical.jl | 537 ++++++++++++++++++++++++++++++------------- test/numerical.jl | 7 +- test/runtests.jl | 3 +- 5 files changed, 419 insertions(+), 158 deletions(-) create mode 100644 src/jump.jl diff --git a/src/ModelAnalyzer.jl b/src/ModelAnalyzer.jl index 7f4040b..b7a8be4 100644 --- a/src/ModelAnalyzer.jl +++ b/src/ModelAnalyzer.jl @@ -109,4 +109,5 @@ include("numerical.jl") include("feasibility.jl") include("infeasibility.jl") +include("jump.jl") end diff --git a/src/jump.jl b/src/jump.jl new file mode 100644 index 0000000..6057cad --- /dev/null +++ b/src/jump.jl @@ -0,0 +1,29 @@ +using JuMP + +# struct JuMPData{T<:AbstractData} <: AbstractData +# data::T +# end + +struct JuMPIssue{T<:AbstractIssue} <: AbstractIssue + issue::T +end + +function analyze( + analyzer::T, + model::JuMP.Model; + kwargs..., +) where {T<:AbstractAnalyzer} + moi_model = JuMP.backend(model) + # Perform the analysis + result = analyze(analyzer, moi_model; kwargs...) + # return JuMPData(result) + return result +end + +# function _name(ref) +# name = JuMP.name(ref) +# if !isempty(name) +# return name +# end +# return "$(ref.index)" +# end \ No newline at end of file diff --git a/src/numerical.jl b/src/numerical.jl index 35d9677..414f742 100644 --- a/src/numerical.jl +++ b/src/numerical.jl @@ -6,15 +6,14 @@ module Numerical import ModelAnalyzer -import JuMP import LinearAlgebra -import JuMP.MOI as MOI import Printf +import MathOptInterface as MOI """ Analyzer() <: ModelAnalyzer.AbstractAnalyzer -The `Analyzer` type is used to analyze the coefficients of a JuMP model for +The `Analyzer` type is used to analyze the coefficients of a model for numerical issues. ## Example @@ -44,7 +43,7 @@ struct Analyzer <: ModelAnalyzer.AbstractAnalyzer end """ AbstractNumericalIssue <: AbstractNumericalIssue -Abstract type for numerical issues found during the analysis of a JuMP model. +Abstract type for numerical issues found during the analysis of a model. """ abstract type AbstractNumericalIssue <: ModelAnalyzer.AbstractIssue end @@ -60,7 +59,7 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.VariableNotInConstraints) ``` """ struct VariableNotInConstraints <: AbstractNumericalIssue - ref::JuMP.VariableRef + ref::MOI.VariableIndex end """ @@ -75,7 +74,7 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.EmptyConstraint) ``` """ struct EmptyConstraint <: AbstractNumericalIssue - ref::JuMP.ConstraintRef + ref::MOI.ConstraintIndex end """ @@ -91,7 +90,7 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.VariableBoundAsConstraint ``` """ struct VariableBoundAsConstraint <: AbstractNumericalIssue - ref::JuMP.ConstraintRef + ref::MOI.ConstraintIndex end """ @@ -107,7 +106,7 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.DenseConstraint) ``` """ struct DenseConstraint <: AbstractNumericalIssue - ref::JuMP.ConstraintRef + ref::MOI.ConstraintIndex nnz::Int end @@ -123,8 +122,8 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.SmallMatrixCoefficient) ``` """ struct SmallMatrixCoefficient <: AbstractNumericalIssue - ref::JuMP.ConstraintRef - variable::JuMP.VariableRef + ref::MOI.ConstraintIndex + variable::MOI.VariableIndex coefficient::Float64 end @@ -140,8 +139,8 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.LargeMatrixCoefficient) ``` """ struct LargeMatrixCoefficient <: AbstractNumericalIssue - ref::JuMP.ConstraintRef - variable::JuMP.VariableRef + ref::MOI.ConstraintIndex + variable::MOI.VariableIndex coefficient::Float64 end @@ -157,7 +156,7 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.SmallBoundCoefficient) ``` """ struct SmallBoundCoefficient <: AbstractNumericalIssue - variable::JuMP.VariableRef + variable::MOI.VariableIndex coefficient::Float64 end @@ -173,7 +172,7 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.LargeBoundCoefficient) ``` """ struct LargeBoundCoefficient <: AbstractNumericalIssue - variable::JuMP.VariableRef + variable::MOI.VariableIndex coefficient::Float64 end @@ -189,7 +188,7 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.SmallRHSCoefficient) ``` """ struct SmallRHSCoefficient <: AbstractNumericalIssue - ref::JuMP.ConstraintRef + ref::MOI.ConstraintIndex coefficient::Float64 end @@ -205,7 +204,7 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.LargeRHSCoefficient) ``` """ struct LargeRHSCoefficient <: AbstractNumericalIssue - ref::JuMP.ConstraintRef + ref::MOI.ConstraintIndex coefficient::Float64 end @@ -221,7 +220,7 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.SmallObjectiveCoefficient ``` """ struct SmallObjectiveCoefficient <: AbstractNumericalIssue - variable::JuMP.VariableRef + variable::MOI.VariableIndex coefficient::Float64 end @@ -237,7 +236,7 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.LargeObjectiveCoefficient ``` """ struct LargeObjectiveCoefficient <: AbstractNumericalIssue - variable::JuMP.VariableRef + variable::MOI.VariableIndex coefficient::Float64 end @@ -255,8 +254,8 @@ julia> ModelAnalyzer.summarize( ``` """ struct SmallObjectiveQuadraticCoefficient <: AbstractNumericalIssue - variable1::JuMP.VariableRef - variable2::JuMP.VariableRef + variable1::MOI.VariableIndex + variable2::MOI.VariableIndex coefficient::Float64 end @@ -274,8 +273,8 @@ julia> ModelAnalyzer.summarize( ``` """ struct LargeObjectiveQuadraticCoefficient <: AbstractNumericalIssue - variable1::JuMP.VariableRef - variable2::JuMP.VariableRef + variable1::MOI.VariableIndex + variable2::MOI.VariableIndex coefficient::Float64 end @@ -293,9 +292,9 @@ julia> ModelAnalyzer.summarize( ``` """ struct SmallMatrixQuadraticCoefficient <: AbstractNumericalIssue - ref::JuMP.ConstraintRef - variable1::JuMP.VariableRef - variable2::JuMP.VariableRef + ref::MOI.ConstraintIndex + variable1::MOI.VariableIndex + variable2::MOI.VariableIndex coefficient::Float64 end @@ -313,9 +312,9 @@ julia> ModelAnalyzer.summarize( ``` """ struct LargeMatrixQuadraticCoefficient <: AbstractNumericalIssue - ref::JuMP.ConstraintRef - variable1::JuMP.VariableRef - variable2::JuMP.VariableRef + ref::MOI.ConstraintIndex + variable1::MOI.VariableIndex + variable2::MOI.VariableIndex coefficient::Float64 end @@ -348,7 +347,7 @@ julia> ModelAnalyzer.summarize( ``` """ struct NonconvexQuadraticConstraint <: AbstractNumericalIssue - ref::JuMP.ConstraintRef + ref::MOI.ConstraintIndex end """ @@ -377,7 +376,7 @@ Base.@kwdef mutable struct Data <: ModelAnalyzer.AbstractData rhs_range::Vector{Float64} = sizehint!(Float64[1.0, 1.0], 2) objective_range::Vector{Float64} = sizehint!(Float64[1.0, 1.0], 2) # cache data - variables_in_constraints::Set{JuMP.VariableRef} = Set{JuMP.VariableRef}() + variables_in_constraints::Set{MOI.VariableIndex} = Set{MOI.VariableIndex}() # variables analysis variables_not_in_constraints::Vector{VariableNotInConstraints} = VariableNotInConstraints[] @@ -400,7 +399,7 @@ Base.@kwdef mutable struct Data <: ModelAnalyzer.AbstractData matrix_quadratic_large::Vector{LargeMatrixQuadraticCoefficient} = LargeMatrixQuadraticCoefficient[] # cache data - sense::JuMP.OptimizationSense = JuMP.FEASIBILITY_SENSE + sense::MOI.OptimizationSense = MOI.FEASIBILITY_SENSE # objective analysis objective_small::Vector{SmallObjectiveCoefficient} = SmallObjectiveCoefficient[] @@ -424,27 +423,15 @@ function _update_range(range::Vector{Float64}, value::Number) return 1 end -function _get_variable_data(data, variable, coefficient::Number) - if !(iszero(coefficient)) - _update_range(data.bounds_range, coefficient) - if abs(coefficient) < data.threshold_small - push!( - data.bounds_small, - SmallBoundCoefficient(variable, coefficient), - ) - elseif abs(coefficient) > data.threshold_large - push!( - data.bounds_large, - LargeBoundCoefficient(variable, coefficient), - ) - end - end +function _get_objective_data(data, func::MOI.VariableIndex) return end -function _get_objective_data(data, func::JuMP.GenericAffExpr) +function _get_objective_data(data, func::MOI.ScalarAffineFunction) nnz = 0 - for (variable, coefficient) in func.terms + for term in func.terms + variable = term.variable + coefficient = term.coefficient if iszero(coefficient) continue end @@ -464,10 +451,19 @@ function _get_objective_data(data, func::JuMP.GenericAffExpr) return end -function _get_objective_data(data, func::JuMP.GenericQuadExpr) - _get_objective_data(data, func.aff) +function _get_objective_data( + data, + func::MOI.ScalarQuadraticFunction{T}, +) where {T} + _get_objective_data( + data, + MOI.ScalarAffineFunction(func.affine_terms, func.constant), + ) nnz = 0 - for (v, coefficient) in func.terms + for term in func.quadratic_terms + coefficient = term.coefficient + v1 = term.variable_1 + v2 = term.variable_2 if iszero(coefficient) continue end @@ -475,21 +471,21 @@ function _get_objective_data(data, func::JuMP.GenericQuadExpr) if abs(coefficient) < data.threshold_small push!( data.objective_quadratic_small, - SmallObjectiveQuadraticCoefficient(v.a, v.b, coefficient), + SmallObjectiveQuadraticCoefficient(v1, v2, coefficient), ) elseif abs(coefficient) > data.threshold_large push!( data.objective_quadratic_large, - LargeObjectiveQuadraticCoefficient(v.a, v.b, coefficient), + LargeObjectiveQuadraticCoefficient(v1, v2, coefficient), ) end end data.has_quadratic_objective = true - if data.sense == JuMP.MAX_SENSE + if data.sense == MOI.MAX_SENSE if !_quadratic_vexity(func, -1) push!(data.nonconvex_objective, NonconvexQuadraticObjective()) end - elseif data.sense == JuMP.MIN_SENSE + elseif data.sense == MOI.MIN_SENSE if !_quadratic_vexity(func, 1) push!(data.nonconvex_objective, NonconvexQuadraticObjective()) end @@ -497,21 +493,26 @@ function _get_objective_data(data, func::JuMP.GenericQuadExpr) return end -function _quadratic_vexity(func::JuMP.GenericQuadExpr, sign::Int) - variables = JuMP.OrderedCollections.OrderedSet{JuMP.VariableRef}() - sizehint!(variables, 2 * length(func.terms)) - for v in keys(func.terms) - push!(variables, v.a) - push!(variables, v.b) +function _quadratic_vexity(func::MOI.ScalarQuadraticFunction, sign::Int) + variables = Set{MOI.VariableIndex}() + sizehint!(variables, 2 * length(func.quadratic_terms)) + for term in func.quadratic_terms + push!(variables, term.variable_1) + push!(variables, term.variable_2) end - var_map = Dict{JuMP.VariableRef,Int}() + var_map = Dict{MOI.VariableIndex,Int}() for (idx, var) in enumerate(variables) var_map[var] = idx end matrix = zeros(length(variables), length(variables)) - for (v, coefficient) in func.terms - matrix[var_map[v.a], var_map[v.b]] += sign * coefficient / 2 - matrix[var_map[v.b], var_map[v.a]] += sign * coefficient / 2 + for term in func.quadratic_terms + coefficient = term.coefficient + v1 = term.variable_1 + v2 = term.variable_2 + matrix[var_map[v1], var_map[v2]] += sign * coefficient / 2 + if v1 != v2 + matrix[var_map[v2], var_map[v1]] += sign * coefficient / 2 + end end ret = LinearAlgebra.cholesky!( LinearAlgebra.Symmetric(matrix), @@ -521,14 +522,41 @@ function _quadratic_vexity(func::JuMP.GenericQuadExpr, sign::Int) return LinearAlgebra.issuccess(ret) end +function _quadratic_vexity(func::MOI.VectorQuadraticFunction{T}, sign) where {T} + n = MOI.output_dimension(func) + quadratic_terms_vector = [MOI.ScalarQuadraticTerm{T}[] for i in 1:n] + for term in func.quadratic_terms + index = term.output_index + push!(quadratic_terms_vector[index], term.scalar_term) + end + for i in 1:n + if length(quadratic_terms_vector[i]) == 0 + continue + end + if !_quadratic_vexity( + MOI.ScalarQuadraticFunction{T}( + quadratic_terms_vector[i], + MOI.ScalarAffineTerm{T}[], + zero(T), + ), + sign, + ) + return false + end + end + return true +end + function _get_constraint_matrix_data( data, - ref::JuMP.ConstraintRef, - func::JuMP.GenericAffExpr; + ref::MOI.ConstraintIndex, + func::MOI.ScalarAffineFunction; ignore_extras = false, ) if length(func.terms) == 1 - if !ignore_extras && isapprox(first(values(func.terms)), 1.0) + coefficient = func.terms[1].coefficient + if !ignore_extras && isapprox(coefficient, 1.0) + # TODO: do this in the vector case push!(data.bound_rows, VariableBoundAsConstraint(ref)) data.matrix_nnz += 1 # in this case we do not count that the variable is in a constraint @@ -536,7 +564,9 @@ function _get_constraint_matrix_data( end end nnz = 0 - for (variable, coefficient) in func.terms + for term in func.terms + variable = term.variable + coefficient = term.coefficient if iszero(coefficient) continue end @@ -570,11 +600,14 @@ end function _get_constraint_matrix_data( data, - ref::JuMP.ConstraintRef, - func::JuMP.GenericQuadExpr, -) + ref::MOI.ConstraintIndex, + func::MOI.ScalarQuadraticFunction{T}, +) where {T} nnz = 0 - for (v, coefficient) in func.terms + for term in func.quadratic_terms + v1 = term.variable_1 + v2 = term.variable_2 + coefficient = term.coefficient if iszero(coefficient) continue end @@ -582,46 +615,113 @@ function _get_constraint_matrix_data( if abs(coefficient) < data.threshold_small push!( data.matrix_quadratic_small, - SmallMatrixQuadraticCoefficient(ref, v.a, v.b, coefficient), + SmallMatrixQuadraticCoefficient(ref, v1, v2, coefficient), ) elseif abs(coefficient) > data.threshold_large push!( data.matrix_quadratic_large, - LargeMatrixQuadraticCoefficient(ref, v.a, v.b, coefficient), + LargeMatrixQuadraticCoefficient(ref, v1, v2, coefficient), ) end - push!(data.variables_in_constraints, v.a) - push!(data.variables_in_constraints, v.b) + push!(data.variables_in_constraints, v1) + push!(data.variables_in_constraints, v2) end data.has_quadratic_constraints = true - _get_constraint_matrix_data(data, ref, func.aff, ignore_extras = nnz > 0) + _get_constraint_matrix_data( + data, + ref, + MOI.ScalarAffineFunction{T}(func.affine_terms, func.constant), + ignore_extras = nnz > 0, + ) return end function _get_constraint_matrix_data( data, - ref, - func::Vector{F}, -) where {F<:Union{JuMP.GenericAffExpr,JuMP.GenericQuadExpr}} - for f in func - _get_constraint_matrix_data(data, ref, f) + ref::MOI.ConstraintIndex, + func::MOI.VectorAffineFunction{T}, +) where {T} + for term in func.terms + variable = term.scalar_term.variable + coefficient = term.scalar_term.coefficient + # index = term.output_index + if iszero(coefficient) + continue + end + _update_range(data.matrix_range, coefficient) + if abs(coefficient) < data.threshold_small + push!( + data.matrix_small, + SmallMatrixCoefficient(ref, variable, coefficient), + ) + elseif abs(coefficient) > data.threshold_large + push!( + data.matrix_large, + LargeMatrixCoefficient(ref, variable, coefficient), + ) + end + push!(data.variables_in_constraints, variable) end - return true + return end -function _get_constraint_data( +function _get_constraint_matrix_data( data, - ref, - func::Vector{F}, - set, -) where {F<:Union{JuMP.GenericAffExpr,JuMP.GenericQuadExpr}} - for f in func - _get_constraint_data(data, ref, f, set) + ref::MOI.ConstraintIndex, + func::MOI.VectorQuadraticFunction{T}, +) where {T} + for term in func.quadratic_terms + v1 = term.scalar_term.variable_1 + v2 = term.scalar_term.variable_2 + coefficient = term.scalar_term.coefficient + if iszero(coefficient) + continue + end + _update_range(data.matrix_quadratic_range, coefficient) + if abs(coefficient) < data.threshold_small + push!( + data.matrix_quadratic_small, + SmallMatrixQuadraticCoefficient(ref, v1, v2, coefficient), + ) + elseif abs(coefficient) > data.threshold_large + push!( + data.matrix_quadratic_large, + LargeMatrixQuadraticCoefficient(ref, v1, v2, coefficient), + ) + end + push!(data.variables_in_constraints, v1) + push!(data.variables_in_constraints, v2) end - return true + _get_constraint_matrix_data( + data, + ref, + MOI.VectorAffineFunction{T}(func.affine_terms, func.constants), + # ignore_extras = nnz > 0, + ) + return +end + +function _get_constraint_matrix_data(data, ref::MOI.ConstraintIndex, func::MOI.VariableIndex) + # push!(data.variables_in_constraints, func) + return +end + +function _get_constraint_matrix_data(data, ref::MOI.ConstraintIndex, func::MOI.VectorOfVariables) + if length(func.variables) == 1 + return + end + for var in func.variables + push!(data.variables_in_constraints, var) + end + return end -function _get_constraint_data(data, ref, func::JuMP.GenericAffExpr, set) +function _get_constraint_data( + data, + ref, + func::Union{MOI.ScalarAffineFunction,MOI.ScalarQuadraticFunction}, + set, +) coefficient = func.constant if iszero(coefficient) return @@ -635,20 +735,58 @@ function _get_constraint_data(data, ref, func::JuMP.GenericAffExpr, set) return end -# this function is basically checking just the RHS just like the above -# skip additional checks for quadratics in non simples sets -function _get_constraint_data(data, ref, func::JuMP.GenericQuadExpr, set) - _get_constraint_data(data, ref, func.aff, set) +function _get_constraint_data( + data, + ref, + func::Union{MOI.VectorAffineFunction,MOI.VectorQuadraticFunction}, + set, +) + coefficients = func.constants + for i in eachindex(coefficients) + coefficient = coefficients[i] + if iszero(coefficient) + continue + end + _update_range(data.rhs_range, coefficient) + if abs(coefficient) < data.threshold_small + push!(data.rhs_small, SmallRHSCoefficient(ref, coefficient)) + elseif abs(coefficient) > data.threshold_large + push!(data.rhs_large, LargeRHSCoefficient(ref, coefficient)) + end + end + return +end + +function _get_constraint_data( + data, + ref, + func::MOI.ScalarQuadraticFunction{T}, + set::MOI.LessThan{T}, +) where {T} + _get_constraint_data( + data, + ref, + MOI.ScalarAffineFunction{T}(func.affine_terms, func.constant), + set, + ) + if !_quadratic_vexity(func, 1) + push!(data.nonconvex_rows, NonconvexQuadraticConstraint(ref)) + end return end function _get_constraint_data( data, ref, - func::JuMP.GenericQuadExpr, - set::Union{MOI.LessThan,MOI.Nonpositives}, -) - _get_constraint_data(data, ref, func.aff, set) + func::MOI.VectorQuadraticFunction{T}, + set::MOI.Nonpositives, +) where {T} + _get_constraint_data( + data, + ref, + MOI.VectorAffineFunction{T}(func.affine_terms, func.constants), + set, + ) if !_quadratic_vexity(func, 1) push!(data.nonconvex_rows, NonconvexQuadraticConstraint(ref)) end @@ -658,7 +796,7 @@ end function _get_constraint_data( data, ref, - func::JuMP.GenericAffExpr, + func::MOI.ScalarAffineFunction, set::MOI.LessThan, ) coefficient = set.upper - func.constant @@ -677,10 +815,33 @@ end function _get_constraint_data( data, ref, - func::JuMP.GenericQuadExpr, - set::Union{MOI.GreaterThan,MOI.Nonnegatives}, -) - _get_constraint_data(data, ref, func.aff, set) + func::MOI.ScalarQuadraticFunction{T}, + set::MOI.GreaterThan{T}, +) where {T} + _get_constraint_data( + data, + ref, + MOI.ScalarAffineFunction{T}(func.affine_terms, func.constant), + set, + ) + if !_quadratic_vexity(func, -1) + push!(data.nonconvex_rows, NonconvexQuadraticConstraint(ref)) + end + return +end + +function _get_constraint_data( + data, + ref, + func::MOI.VectorQuadraticFunction{T}, + set::MOI.Nonnegatives, +) where {T} + _get_constraint_data( + data, + ref, + MOI.VectorAffineFunction{T}(func.affine_terms, func.constants), + set, + ) if !_quadratic_vexity(func, -1) push!(data.nonconvex_rows, NonconvexQuadraticConstraint(ref)) end @@ -690,7 +851,7 @@ end function _get_constraint_data( data, ref, - func::JuMP.GenericAffExpr, + func::MOI.ScalarAffineFunction, set::MOI.GreaterThan, ) coefficient = set.lower - func.constant @@ -709,10 +870,31 @@ end function _get_constraint_data( data, ref, - func::JuMP.GenericQuadExpr, - set::Union{MOI.EqualTo,MOI.Interval,MOI.Zeros}, + func::MOI.ScalarQuadraticFunction, + set::Union{MOI.EqualTo,MOI.Interval}, ) - _get_constraint_data(data, ref, func.aff, set) + _get_constraint_data( + data, + ref, + MOI.ScalarAffineFunction(func.affine_terms, func.constant), + set, + ) + push!(data.nonconvex_rows, NonconvexQuadraticConstraint(ref)) + return +end + +function _get_constraint_data( + data, + ref, + func::MOI.VectorQuadraticFunction, + set::MOI.Zeros, +) + _get_constraint_data( + data, + ref, + MOI.VectorAffineFunction(func.affine_terms, func.constants), + set, + ) push!(data.nonconvex_rows, NonconvexQuadraticConstraint(ref)) return end @@ -720,7 +902,7 @@ end function _get_constraint_data( data, ref, - func::JuMP.GenericAffExpr, + func::MOI.ScalarAffineFunction, set::MOI.EqualTo, ) coefficient = set.value - func.constant @@ -739,7 +921,7 @@ end function _get_constraint_data( data, ref, - func::JuMP.GenericAffExpr, + func::MOI.ScalarAffineFunction, set::MOI.Interval, ) coefficient = set.upper - func.constant @@ -764,22 +946,83 @@ function _get_constraint_data( return end -function _get_constraint_matrix_data(data, func::Vector{JuMP.VariableRef}) - for var in func - push!(data.variables_in_constraints, var) +function _get_constraint_data( + data, + ref, + func::MOI.VariableIndex, + set::MOI.LessThan, +) + _get_variable_data(data, func, set.upper) + return +end + +function _get_constraint_data( + data, + ref, + func::MOI.VariableIndex, + set::MOI.GreaterThan, +) + _get_variable_data(data, func, set.lower) + return +end + +function _get_constraint_data( + data, + ref, + func::MOI.VariableIndex, + set::MOI.EqualTo, +) + _get_variable_data(data, func, set.value) + return +end + +function _get_constraint_data( + data, + ref, + func::MOI.VariableIndex, + set::MOI.Interval, +) + _get_variable_data(data, func, set.lower) + _get_variable_data(data, func, set.upper) + return +end + +function _get_constraint_data(data, ref, func::MOI.VariableIndex, set) + return +end + +function _get_variable_data(data, variable, coefficient::Number) + if !(iszero(coefficient)) + _update_range(data.bounds_range, coefficient) + if abs(coefficient) < data.threshold_small + push!( + data.bounds_small, + SmallBoundCoefficient(variable, coefficient), + ) + elseif abs(coefficient) > data.threshold_large + push!( + data.bounds_large, + LargeBoundCoefficient(variable, coefficient), + ) + end end return end +function _get_constraint_data(data, ref, func::MOI.VectorOfVariables, set) + return +end + """ analyze(model::Model; threshold_dense_fill_in = 0.10, threshold_dense_entries = 1000, threshold_small = 1e-5, threshold_large = 1e+5) -Analyze the coefficients of a JuMP model. +Analyze the coefficients of a model. """ function ModelAnalyzer.analyze( ::Analyzer, - model::JuMP.Model; + model::MOI.ModelLike, + ; threshold_dense_fill_in::Float64 = 0.10, threshold_dense_entries::Int = 1000, threshold_small::Float64 = 1e-5, @@ -792,47 +1035,37 @@ function ModelAnalyzer.analyze( data.threshold_large = threshold_large # initialize simples data - data.sense = JuMP.objective_sense(model) - data.number_of_variables = JuMP.num_variables(model) + data.sense = MOI.get(model, MOI.ObjectiveSense()) + data.number_of_variables = MOI.get(model, MOI.NumberOfVariables()) sizehint!(data.variables_in_constraints, data.number_of_variables) - data.number_of_constraints = - JuMP.num_constraints(model, count_variable_in_set_constraints = false) + # objective pass - _get_objective_data(data, JuMP.objective_function(model)) - # variables pass - for var in JuMP.all_variables(model) - if JuMP.has_lower_bound(var) - _get_variable_data(data, var, JuMP.lower_bound(var)) - end - if JuMP.has_upper_bound(var) - _get_variable_data(data, var, JuMP.upper_bound(var)) - end - if JuMP.is_fixed(var) - _get_variable_data(data, var, JuMP.fix_value(var)) - end - end + objective_type = MOI.get(model, MOI.ObjectiveFunctionType()) + obj_func = MOI.get(model, MOI.ObjectiveFunction{objective_type}()) + _get_objective_data(data, obj_func) + # constraints pass - for (F, S) in JuMP.list_of_constraint_types(model) - n = JuMP.num_constraints(model, F, S) + data.number_of_constraints = 0 + list_of_constraint_types = + MOI.get(model, MOI.ListOfConstraintTypesPresent()) + for (F, S) in list_of_constraint_types + list = MOI.get(model, MOI.ListOfConstraintIndices{F,S}()) + n = length(list) + data.number_of_constraints += n if n > 0 push!(data.constraint_info, (F, S, n)) end - F == JuMP.VariableRef && continue - if F == Vector{JuMP.VariableRef} - for con in JuMP.all_constraints(model, F, S) - con_obj = JuMP.constraint_object(con) - _get_constraint_matrix_data(data, con_obj.func) - end - continue - end - for con in JuMP.all_constraints(model, F, S) - con_obj = JuMP.constraint_object(con) - _get_constraint_matrix_data(data, con, con_obj.func) - _get_constraint_data(data, con, con_obj.func, con_obj.set) + for con in list + func = MOI.get(model, MOI.ConstraintFunction(), con) + set = MOI.get(model, MOI.ConstraintSet(), con) + _get_constraint_matrix_data(data, con, func) + _get_constraint_data(data, con, func, set) end end # second pass on variables after constraint pass - for var in JuMP.all_variables(model) + # variable index constraints are not counted in the constraints pass + list_of_variables = MOI.get(model, MOI.ListOfVariableIndices()) + for var in list_of_variables if !(var in data.variables_in_constraints) push!( data.variables_not_in_constraints, @@ -2085,11 +2318,7 @@ end # printing helpers function _name(ref) - name = JuMP.name(ref) - if !isempty(name) - return name - end - return "$(ref.index)" + return "$(ref)" end _print_value(x::Real) = Printf.@sprintf("%1.0e", x) diff --git a/test/numerical.jl b/test/numerical.jl index 1a6e98d..cdd2e3a 100644 --- a/test/numerical.jl +++ b/test/numerical.jl @@ -82,21 +82,22 @@ function test_constraint_bounds() @constraint(model, x == 4e+13) @constraint(model, x <= 1e-14) @constraint(model, x <= 1e+15) + @constraint(model, [x] in MOI.Nonnegatives(1)) @constraint(model, [x - 1e-16] in MOI.Nonnegatives(1)) @constraint(model, [x - 1e+17] in MOI.Nonnegatives(1)) data = ModelAnalyzer.analyze(ModelAnalyzer.Numerical.Analyzer(), model) list = ModelAnalyzer.list_of_issue_types(data) - @test length(list) == 4 + @test length(list) == 3 ret = ModelAnalyzer.list_of_issues( data, ModelAnalyzer.Numerical.VariableBoundAsConstraint, ) - @test length(ret) == 8 + @test length(ret) == 6 ret = ModelAnalyzer.list_of_issues( data, ModelAnalyzer.Numerical.VariableNotInConstraints, ) - @test length(ret) == 1 + @test length(ret) == 0 ret = ModelAnalyzer.list_of_issues( data, ModelAnalyzer.Numerical.SmallRHSCoefficient, diff --git a/test/runtests.jl b/test/runtests.jl index 12287cd..fb558e5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,7 +6,8 @@ using Test @testset "ModelAnalyzer" begin - for file in readdir(@__DIR__) + for file in ["numerical.jl"] + # for file in readdir(@__DIR__) if !endswith(file, ".jl") || file in ("runtests.jl",) continue end From 4f716b06a4183b4c876107d070f0c2d2418e86ee Mon Sep 17 00:00:00 2001 From: joaquimg Date: Thu, 1 May 2025 01:44:18 -0300 Subject: [PATCH 02/34] format --- src/jump.jl | 2 +- src/numerical.jl | 12 ++++++++++-- 2 files changed, 11 insertions(+), 3 deletions(-) diff --git a/src/jump.jl b/src/jump.jl index 6057cad..3720d3f 100644 --- a/src/jump.jl +++ b/src/jump.jl @@ -26,4 +26,4 @@ end # return name # end # return "$(ref.index)" -# end \ No newline at end of file +# end diff --git a/src/numerical.jl b/src/numerical.jl index 414f742..728ee83 100644 --- a/src/numerical.jl +++ b/src/numerical.jl @@ -701,12 +701,20 @@ function _get_constraint_matrix_data( return end -function _get_constraint_matrix_data(data, ref::MOI.ConstraintIndex, func::MOI.VariableIndex) +function _get_constraint_matrix_data( + data, + ref::MOI.ConstraintIndex, + func::MOI.VariableIndex, +) # push!(data.variables_in_constraints, func) return end -function _get_constraint_matrix_data(data, ref::MOI.ConstraintIndex, func::MOI.VectorOfVariables) +function _get_constraint_matrix_data( + data, + ref::MOI.ConstraintIndex, + func::MOI.VectorOfVariables, +) if length(func.variables) == 1 return end From 3185a6f38b384f92b19705aa8708dbd3e6f5f373 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Thu, 1 May 2025 02:35:17 -0300 Subject: [PATCH 03/34] Improve API --- src/ModelAnalyzer.jl | 90 ++++++++++- src/jump.jl | 82 ++++++++-- src/numerical.jl | 345 +++++++++++++++++++++++++++++++++---------- 3 files changed, 419 insertions(+), 98 deletions(-) diff --git a/src/ModelAnalyzer.jl b/src/ModelAnalyzer.jl index b7a8be4..7b5f496 100644 --- a/src/ModelAnalyzer.jl +++ b/src/ModelAnalyzer.jl @@ -25,7 +25,7 @@ See [`summarize`](@ref), [`list_of_issues`](@ref), and function analyze end """ - summarize([io::IO,] AbstractData; verbose = true, max_issues = typemax(Int), kwargs...) + summarize([io::IO,] AbstractData; name_source = nothing, verbose = true, max_issues = typemax(Int), kwargs...) Print a summary of the analysis results contained in `AbstractData` to the specified IO stream. If no IO stream is provided, it defaults to `stdout`. @@ -41,9 +41,11 @@ be a subtype of `AbstractIssue`). In the verbose case it will provide a text explaning the issue. In the non-verbose case it will provide just the issue name. - summarize([io::IO,] issue::AbstractIssue; verbose = true) + summarize([io::IO,] issue::AbstractIssue; name_source = nothing, verbose = true) This variant allows summarizing a single issue instance of type `AbstractIssue`. +The model tha led to the issue can be provided to `name_source`, it will be used +to generate the name of variables and constraints in the issue summary. """ function summarize end @@ -72,17 +74,23 @@ function summarize(io::IO, ::Type{T}; verbose = true) where {T<:AbstractIssue} end end -function summarize(io::IO, issue::AbstractIssue; verbose = true) +function summarize( + io::IO, + issue::AbstractIssue; + name_source = nothing, + verbose = true, +) if verbose - return _verbose_summarize(io, issue) + return _verbose_summarize(io, issue, name_source) else - return _summarize(io, issue) + return _summarize(io, issue, name_source) end end function summarize( io::IO, issues::Vector{T}; + name_source = nothing, verbose = true, max_issues = typemax(Int), ) where {T<:AbstractIssue} @@ -92,7 +100,7 @@ function summarize( print(io, "\n\n## List of issues\n\n") for issue in first(issues, max_issues) print(io, " * ") - summarize(io, issue, verbose = verbose) + summarize(io, issue, verbose = verbose, name_source = name_source) print(io, "\n") end return @@ -102,9 +110,79 @@ function summarize(data::AbstractData; kwargs...) return summarize(stdout, data; kwargs...) end +""" + value(issue::AbstractIssue) + +Return the value associated to a particular issue. The value is a number +with a different meaning depending on the type of issue. For example, for +some numerical issues, it can be the coefficient value. +""" +function value(issue::AbstractIssue, ::Nothing) + return value(issue) +end +function value(issue::AbstractIssue, ::MOI.ModelLike) + return value(issue) +end + +""" + variable(issue::AbstractIssue) + +Return the variable associated to a particular issue. +""" +function variable(issue::AbstractIssue, ::Nothing) + return variable(issue) +end +function variable(issue::AbstractIssue, ::MOI.ModelLike) + return variable(issue) +end + +""" + variables(issue::AbstractIssue) + +Return the variables associated to a particular issue. +""" +function variables(issue::AbstractIssue, ::Nothing) + return variables(issue) +end +function variables(issue::AbstractIssue, ::MOI.ModelLike) + return variables(issue) +end + +""" + constraint(issue::AbstractIssue) + +Return the constraint associated to a particular issue. +""" +function constraint(issue::AbstractIssue, ::Nothing) + return constraint(issue) +end +function constraint(issue::AbstractIssue, ::MOI.ModelLike) + return constraint(issue) +end + function _verbose_summarize end function _summarize end +function _name(ref::MOI.VariableIndex, model::MOI.ModelLike) + name = MOI.get(model, MOI.VariableName(), ref) + if !isempty(name) + return name + end + return "$ref" +end + +function _name(ref::MOI.ConstraintIndex, model::MOI.ModelLike) + name = MOI.get(model, MOI.ConstraintName(), ref) + if !isempty(name) + return name + end + return "$ref" +end + +function _name(ref, ::Nothing) + return "$ref" +end + include("numerical.jl") include("feasibility.jl") include("infeasibility.jl") diff --git a/src/jump.jl b/src/jump.jl index 3720d3f..437f6b5 100644 --- a/src/jump.jl +++ b/src/jump.jl @@ -1,29 +1,79 @@ using JuMP -# struct JuMPData{T<:AbstractData} <: AbstractData +# struct JuMPData{T<:AbstractData} <: ModelAnalyzer.AbstractData # data::T +# model::JuMP.Model # end -struct JuMPIssue{T<:AbstractIssue} <: AbstractIssue - issue::T -end +# struct JuMPIssue{T<:AbstractIssue} <: ModelAnalyzer.AbstractIssue +# issue::T +# model::JuMP.Model +# end -function analyze( +function ModelAnalyzer.analyze( analyzer::T, model::JuMP.Model; kwargs..., -) where {T<:AbstractAnalyzer} +) where {T<:ModelAnalyzer.AbstractAnalyzer} moi_model = JuMP.backend(model) - # Perform the analysis - result = analyze(analyzer, moi_model; kwargs...) - # return JuMPData(result) + result = ModelAnalyzer.analyze(analyzer, moi_model; kwargs...) + # return JuMPData(result, model) return result end -# function _name(ref) -# name = JuMP.name(ref) -# if !isempty(name) -# return name -# end -# return "$(ref.index)" -# end +function ModelAnalyzer._name(ref::MOI.VariableIndex, model::JuMP.Model) + jump_ref = JuMP.VariableRef(model, ref) + name = JuMP.name(jump_ref) + if !isempty(name) + return name + end + return "$jump_ref" +end + +function ModelAnalyzer._name(ref::MOI.ConstraintIndex, model::JuMP.Model) + jump_ref = JuMP.constraint_ref_with_index(model, ref) + name = JuMP.name(jump_ref) + if !isempty(name) + return name + end + return "$jump_ref" +end + +""" + variable(issue::ModelAnalyzer.AbstractIssue, model::JuMP.Model) + +Return the **JuMP** variable reference associated to a particular issue. +""" +function ModelAnalyzer.variable( + issue::ModelAnalyzer.AbstractIssue, + model::JuMP.Model, +) + ref = ModelAnalyzer.variable(issue) + return JuMP.VariableRef(model, ref) +end + +""" + variables(issue::ModelAnalyzer.AbstractIssue, model::JuMP.Model) + +Return the **JuMP** variable references associated to a particular issue. +""" +function ModelAnalyzer.variables( + issue::ModelAnalyzer.AbstractIssue, + model::JuMP.Model, +) + refs = ModelAnalyzer.variables(issue) + return JuMP.VariableRef.(model, refs) +end + +""" + constraint(issue::ModelAnalyzer.AbstractIssue, model::JuMP.Model) + +Return the **JuMP** constraint reference associated to a particular issue. +""" +function ModelAnalyzer.constraint( + issue::ModelAnalyzer.AbstractIssue, + model::JuMP.Model, +) + ref = ModelAnalyzer.constraint(issue) + return JuMP.constraint_ref_with_index(model, ref) +end diff --git a/src/numerical.jl b/src/numerical.jl index 728ee83..854a3ce 100644 --- a/src/numerical.jl +++ b/src/numerical.jl @@ -61,6 +61,7 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.VariableNotInConstraints) struct VariableNotInConstraints <: AbstractNumericalIssue ref::MOI.VariableIndex end +ModelAnalyzer.variable(issue::VariableNotInConstraints, model) = issue.ref """ EmptyConstraint <: AbstractNumericalIssue @@ -76,6 +77,7 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.EmptyConstraint) struct EmptyConstraint <: AbstractNumericalIssue ref::MOI.ConstraintIndex end +ModelAnalyzer.constraint(issue::EmptyConstraint, model) = issue.ref """ VariableBoundAsConstraint <: AbstractNumericalIssue @@ -92,6 +94,7 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.VariableBoundAsConstraint struct VariableBoundAsConstraint <: AbstractNumericalIssue ref::MOI.ConstraintIndex end +ModelAnalyzer.constraint(issue::EmptyConstraint, model) = issue.ref """ DenseConstraint <: AbstractNumericalIssue @@ -109,6 +112,8 @@ struct DenseConstraint <: AbstractNumericalIssue ref::MOI.ConstraintIndex nnz::Int end +ModelAnalyzer.constraint(issue::DenseConstraint, model) = issue.ref +ModelAnalyzer.value(issue::DenseConstraint) = issue.nnz """ SmallMatrixCoefficient <: AbstractNumericalIssue @@ -126,6 +131,9 @@ struct SmallMatrixCoefficient <: AbstractNumericalIssue variable::MOI.VariableIndex coefficient::Float64 end +ModelAnalyzer.variable(issue::SmallMatrixCoefficient, model) = issue.variable +ModelAnalyzer.constraint(issue::SmallMatrixCoefficient, model) = issue.ref +ModelAnalyzer.value(issue::SmallMatrixCoefficient) = issue.coefficient """ LargeMatrixCoefficient <: AbstractNumericalIssue @@ -143,6 +151,9 @@ struct LargeMatrixCoefficient <: AbstractNumericalIssue variable::MOI.VariableIndex coefficient::Float64 end +ModelAnalyzer.variable(issue::LargeMatrixCoefficient, model) = issue.variable +ModelAnalyzer.constraint(issue::LargeMatrixCoefficient, model) = issue.ref +ModelAnalyzer.value(issue::LargeMatrixCoefficient) = issue.coefficient """ SmallBoundCoefficient <: AbstractNumericalIssue @@ -159,6 +170,8 @@ struct SmallBoundCoefficient <: AbstractNumericalIssue variable::MOI.VariableIndex coefficient::Float64 end +ModelAnalyzer.variable(issue::SmallBoundCoefficient, model) = issue.variable +ModelAnalyzer.value(issue::SmallBoundCoefficient) = issue.coefficient """ LargeBoundCoefficient <: AbstractNumericalIssue @@ -175,6 +188,8 @@ struct LargeBoundCoefficient <: AbstractNumericalIssue variable::MOI.VariableIndex coefficient::Float64 end +ModelAnalyzer.variable(issue::LargeBoundCoefficient, model) = issue.variable +ModelAnalyzer.value(issue::LargeBoundCoefficient) = issue.coefficient """ SmallRHSCoefficient <: AbstractNumericalIssue @@ -191,6 +206,8 @@ struct SmallRHSCoefficient <: AbstractNumericalIssue ref::MOI.ConstraintIndex coefficient::Float64 end +ModelAnalyzer.constraint(issue::SmallRHSCoefficient, model) = issue.ref +ModelAnalyzer.value(issue::SmallRHSCoefficient) = issue.coefficient """ LargeRHSCoefficient <: AbstractNumericalIssue @@ -207,6 +224,8 @@ struct LargeRHSCoefficient <: AbstractNumericalIssue ref::MOI.ConstraintIndex coefficient::Float64 end +ModelAnalyzer.constraint(issue::LargeRHSCoefficient, model) = issue.ref +ModelAnalyzer.value(issue::LargeRHSCoefficient) = issue.coefficient """ SmallObjectiveCoefficient <: AbstractNumericalIssue @@ -223,6 +242,8 @@ struct SmallObjectiveCoefficient <: AbstractNumericalIssue variable::MOI.VariableIndex coefficient::Float64 end +ModelAnalyzer.variable(issue::SmallObjectiveCoefficient, model) = issue.variable +ModelAnalyzer.value(issue::SmallObjectiveCoefficient) = issue.coefficient """ LargeObjectiveCoefficient <: AbstractNumericalIssue @@ -239,6 +260,8 @@ struct LargeObjectiveCoefficient <: AbstractNumericalIssue variable::MOI.VariableIndex coefficient::Float64 end +ModelAnalyzer.variable(issue::LargeObjectiveCoefficient, model) = issue.variable +ModelAnalyzer.value(issue::LargeObjectiveCoefficient) = issue.coefficient """ SmallObjectiveQuadraticCoefficient <: AbstractNumericalIssue @@ -258,6 +281,15 @@ struct SmallObjectiveQuadraticCoefficient <: AbstractNumericalIssue variable2::MOI.VariableIndex coefficient::Float64 end +function ModelAnalyzer.variables( + issue::SmallObjectiveQuadraticCoefficient, + model, +) + return [issue.variable1, issue.variable2] +end +function ModelAnalyzer.value(issue::SmallObjectiveQuadraticCoefficient) + return issue.coefficient +end """ LargeObjectiveQuadraticCoefficient <: AbstractNumericalIssue @@ -277,6 +309,15 @@ struct LargeObjectiveQuadraticCoefficient <: AbstractNumericalIssue variable2::MOI.VariableIndex coefficient::Float64 end +function ModelAnalyzer.variables( + issue::LargeObjectiveQuadraticCoefficient, + model, +) + return [issue.variable1, issue.variable2] +end +function ModelAnalyzer.value(issue::LargeObjectiveQuadraticCoefficient) + return issue.coefficient +end """ SmallMatrixQuadraticCoefficient <: AbstractNumericalIssue @@ -297,6 +338,13 @@ struct SmallMatrixQuadraticCoefficient <: AbstractNumericalIssue variable2::MOI.VariableIndex coefficient::Float64 end +function ModelAnalyzer.variables(issue::SmallMatrixQuadraticCoefficient, model) + return [issue.variable1, issue.variable2] +end +function ModelAnalyzer.constraint(issue::SmallMatrixQuadraticCoefficient, model) + return issue.ref +end +ModelAnalyzer.value(issue::SmallMatrixQuadraticCoefficient) = issue.coefficient """ LargeMatrixQuadraticCoefficient <: AbstractNumericalIssue @@ -317,6 +365,13 @@ struct LargeMatrixQuadraticCoefficient <: AbstractNumericalIssue variable2::MOI.VariableIndex coefficient::Float64 end +function ModelAnalyzer.variables(issue::LargeMatrixQuadraticCoefficient, model) + return [issue.variable1, issue.variable2] +end +function ModelAnalyzer.constraint(issue::LargeMatrixQuadraticCoefficient, model) + return issue.ref +end +ModelAnalyzer.value(issue::LargeMatrixQuadraticCoefficient) = issue.coefficient """ NonconvexQuadraticObjective <: AbstractNumericalIssue @@ -349,6 +404,7 @@ julia> ModelAnalyzer.summarize( struct NonconvexQuadraticConstraint <: AbstractNumericalIssue ref::MOI.ConstraintIndex end +ModelAnalyzer.constraint(issue::NonconvexQuadraticConstraint, model) = issue.ref """ Data @@ -1768,77 +1824,153 @@ function ModelAnalyzer._verbose_summarize( ) end -function ModelAnalyzer._summarize(io::IO, issue::VariableNotInConstraints) - return print(io, _name(issue.ref)) +function ModelAnalyzer._summarize( + io::IO, + issue::VariableNotInConstraints, + name_source, +) + return print(io, ModelAnalyzer._name(issue.ref, name_source)) end -function ModelAnalyzer._summarize(io::IO, issue::EmptyConstraint) - return print(io, _name(issue.ref)) +function ModelAnalyzer._summarize(io::IO, issue::EmptyConstraint, name_source) + return print(io, ModelAnalyzer._name(issue.ref, name_source)) end -function ModelAnalyzer._summarize(io::IO, issue::VariableBoundAsConstraint) - return print(io, _name(issue.ref)) +function ModelAnalyzer._summarize( + io::IO, + issue::VariableBoundAsConstraint, + name_source, +) + return print(io, ModelAnalyzer._name(issue.ref, name_source)) end -function ModelAnalyzer._summarize(io::IO, issue::DenseConstraint) - return print(io, _name(issue.ref), " : ", issue.nnz) +function ModelAnalyzer._summarize(io::IO, issue::DenseConstraint, name_source) + return print( + io, + ModelAnalyzer._name(issue.ref, name_source), + " : ", + issue.nnz, + ) end -function ModelAnalyzer._summarize(io::IO, issue::SmallMatrixCoefficient) +function ModelAnalyzer._summarize( + io::IO, + issue::SmallMatrixCoefficient, + name_source, +) return print( io, - _name(issue.ref), + ModelAnalyzer._name(issue.ref, name_source), " -- ", - _name(issue.variable), + ModelAnalyzer._name(issue.variable, name_source), " : ", issue.coefficient, ) end -function ModelAnalyzer._summarize(io::IO, issue::LargeMatrixCoefficient) +function ModelAnalyzer._summarize( + io::IO, + issue::LargeMatrixCoefficient, + name_source, +) return print( io, - _name(issue.ref), + ModelAnalyzer._name(issue.ref, name_source), " -- ", - _name(issue.variable), + ModelAnalyzer._name(issue.variable, name_source), " : ", issue.coefficient, ) end -function ModelAnalyzer._summarize(io::IO, issue::SmallBoundCoefficient) - return print(io, _name(issue.variable), " : ", issue.coefficient) +function ModelAnalyzer._summarize( + io::IO, + issue::SmallBoundCoefficient, + name_source, +) + return print( + io, + ModelAnalyzer._name(issue.variable, name_source), + " : ", + issue.coefficient, + ) end -function ModelAnalyzer._summarize(io::IO, issue::LargeBoundCoefficient) - return print(io, _name(issue.variable), " : ", issue.coefficient) +function ModelAnalyzer._summarize( + io::IO, + issue::LargeBoundCoefficient, + name_source, +) + return print( + io, + ModelAnalyzer._name(issue.variable, name_source), + " : ", + issue.coefficient, + ) end -function ModelAnalyzer._summarize(io::IO, issue::SmallRHSCoefficient) - return print(io, _name(issue.ref), " : ", issue.coefficient) +function ModelAnalyzer._summarize( + io::IO, + issue::SmallRHSCoefficient, + name_source, +) + return print( + io, + ModelAnalyzer._name(issue.ref, name_source), + " : ", + issue.coefficient, + ) end -function ModelAnalyzer._summarize(io::IO, issue::LargeRHSCoefficient) - return print(io, _name(issue.ref), " : ", issue.coefficient) +function ModelAnalyzer._summarize( + io::IO, + issue::LargeRHSCoefficient, + name_source, +) + return print( + io, + ModelAnalyzer._name(issue.ref, name_source), + " : ", + issue.coefficient, + ) end -function ModelAnalyzer._summarize(io::IO, issue::SmallObjectiveCoefficient) - return print(io, _name(issue.variable), " : ", issue.coefficient) +function ModelAnalyzer._summarize( + io::IO, + issue::SmallObjectiveCoefficient, + name_source, +) + return print( + io, + ModelAnalyzer._name(issue.variable, name_source), + " : ", + issue.coefficient, + ) end -function ModelAnalyzer._summarize(io::IO, issue::LargeObjectiveCoefficient) - return print(io, _name(issue.variable), " : ", issue.coefficient) +function ModelAnalyzer._summarize( + io::IO, + issue::LargeObjectiveCoefficient, + name_source, +) + return print( + io, + ModelAnalyzer._name(issue.variable, name_source), + " : ", + issue.coefficient, + ) end function ModelAnalyzer._summarize( io::IO, issue::SmallObjectiveQuadraticCoefficient, + name_source, ) return print( io, - _name(issue.variable1), + ModelAnalyzer._name(issue.variable1, name_source), " -- ", - _name(issue.variable2), + ModelAnalyzer._name(issue.variable2, name_source), " : ", issue.coefficient, ) @@ -1847,12 +1979,13 @@ end function ModelAnalyzer._summarize( io::IO, issue::LargeObjectiveQuadraticCoefficient, + name_source, ) return print( io, - _name(issue.variable1), + ModelAnalyzer._name(issue.variable1, name_source), " -- ", - _name(issue.variable2), + ModelAnalyzer._name(issue.variable2, name_source), " : ", issue.coefficient, ) @@ -1861,14 +1994,15 @@ end function ModelAnalyzer._summarize( io::IO, issue::SmallMatrixQuadraticCoefficient, + name_source, ) return print( io, - _name(issue.ref), + ModelAnalyzer._name(issue.ref, name_source), " -- ", - _name(issue.variable1), + ModelAnalyzer._name(issue.variable1, name_source), " -- ", - _name(issue.variable2), + ModelAnalyzer._name(issue.variable2, name_source), " : ", issue.coefficient, ) @@ -1877,115 +2011,166 @@ end function ModelAnalyzer._summarize( io::IO, issue::LargeMatrixQuadraticCoefficient, + name_source, ) return print( io, - _name(issue.ref), + ModelAnalyzer._name(issue.ref, name_source), " -- ", - _name(issue.variable1), + ModelAnalyzer._name(issue.variable1, name_source), " -- ", - _name(issue.variable2), + ModelAnalyzer._name(issue.variable2, name_source), " : ", issue.coefficient, ) end -function ModelAnalyzer._summarize(io::IO, ::NonconvexQuadraticObjective) +function ModelAnalyzer._summarize( + io::IO, + ::NonconvexQuadraticObjective, + name_source, +) return print(io, "Objective is Nonconvex quadratic") end -function ModelAnalyzer._summarize(io::IO, issue::NonconvexQuadraticConstraint) - return print(io, _name(issue.ref)) +function ModelAnalyzer._summarize( + io::IO, + issue::NonconvexQuadraticConstraint, + name_source, +) + return print(io, ModelAnalyzer._name(issue.ref, name_source)) end function ModelAnalyzer._verbose_summarize( io::IO, issue::VariableNotInConstraints, + name_source, ) - return print(io, "Variable: ", _name(issue.ref)) + return print(io, "Variable: ", ModelAnalyzer._name(issue.ref, name_source)) end -function ModelAnalyzer._verbose_summarize(io::IO, issue::EmptyConstraint) - return print(io, "Constraint: ", _name(issue.ref)) +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::EmptyConstraint, + name_source, +) + return print( + io, + "Constraint: ", + ModelAnalyzer._name(issue.ref, name_source), + ) end function ModelAnalyzer._verbose_summarize( io::IO, issue::VariableBoundAsConstraint, + name_source, ) - return print(io, "Constraint: ", _name(issue.ref)) + return print( + io, + "Constraint: ", + ModelAnalyzer._name(issue.ref, name_source), + ) end -function ModelAnalyzer._verbose_summarize(io::IO, issue::DenseConstraint) +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::DenseConstraint, + name_source, +) return print( io, "Constraint: ", - _name(issue.ref), + ModelAnalyzer._name(issue.ref, name_source), " with ", issue.nnz, " non zero coefficients", ) end -function ModelAnalyzer._verbose_summarize(io::IO, issue::SmallMatrixCoefficient) +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::SmallMatrixCoefficient, + name_source, +) return print( io, "(Constraint -- Variable): (", - _name(issue.ref), + ModelAnalyzer._name(issue.ref, name_source), " -- ", - _name(issue.variable), + ModelAnalyzer._name(issue.variable, name_source), ") with coefficient ", issue.coefficient, ) end -function ModelAnalyzer._verbose_summarize(io::IO, issue::LargeMatrixCoefficient) +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::LargeMatrixCoefficient, + name_source, +) return print( io, "(Constraint -- Variable): (", - _name(issue.ref), + ModelAnalyzer._name(issue.ref, name_source), " -- ", - _name(issue.variable), + ModelAnalyzer._name(issue.variable, name_source), ") with coefficient ", issue.coefficient, ) end -function ModelAnalyzer._verbose_summarize(io::IO, issue::SmallBoundCoefficient) +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::SmallBoundCoefficient, + name_source, +) return print( io, "Variable: ", - _name(issue.variable), + ModelAnalyzer._name(issue.variable, name_source), " with bound ", issue.coefficient, ) end -function ModelAnalyzer._verbose_summarize(io::IO, issue::LargeBoundCoefficient) +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::LargeBoundCoefficient, + name_source, +) return print( io, "Variable: ", - _name(issue.variable), + ModelAnalyzer._name(issue.variable, name_source), " with bound ", issue.coefficient, ) end -function ModelAnalyzer._verbose_summarize(io::IO, issue::SmallRHSCoefficient) +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::SmallRHSCoefficient, + name_source, +) return print( io, "Constraint: ", - _name(issue.ref), + ModelAnalyzer._name(issue.ref, name_source), " with right-hand-side ", issue.coefficient, ) end -function ModelAnalyzer._verbose_summarize(io::IO, issue::LargeRHSCoefficient) +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::LargeRHSCoefficient, + name_source, +) return print( io, "Constraint: ", - _name(issue.ref), + ModelAnalyzer._name(issue.ref, name_source), " with right-hand-side ", issue.coefficient, ) @@ -1994,11 +2179,12 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::SmallObjectiveCoefficient, + name_source, ) return print( io, "Variable: ", - _name(issue.variable), + ModelAnalyzer._name(issue.variable, name_source), " with coefficient ", issue.coefficient, ) @@ -2007,11 +2193,12 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::LargeObjectiveCoefficient, + name_source, ) return print( io, "Variable: ", - _name(issue.variable), + ModelAnalyzer._name(issue.variable, name_source), " with coefficient ", issue.coefficient, ) @@ -2020,13 +2207,14 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::SmallObjectiveQuadraticCoefficient, + name_source, ) return print( io, "(Variable -- Variable): (", - _name(issue.variable1), + ModelAnalyzer._name(issue.variable1, name_source), " -- ", - _name(issue.variable2), + ModelAnalyzer._name(issue.variable2, name_source), ") with coefficient ", issue.coefficient, ) @@ -2035,13 +2223,14 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::LargeObjectiveQuadraticCoefficient, + name_source, ) return print( io, "(Variable -- Variable): (", - _name(issue.variable1), + ModelAnalyzer._name(issue.variable1, name_source), " -- ", - _name(issue.variable2), + ModelAnalyzer._name(issue.variable2, name_source), ") with coefficient ", issue.coefficient, ) @@ -2050,15 +2239,16 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::SmallMatrixQuadraticCoefficient, + name_source, ) return print( io, "(Constraint -- Variable -- Variable): (", - _name(issue.ref), + ModelAnalyzer._name(issue.ref, name_source), " -- ", - _name(issue.variable1), + ModelAnalyzer._name(issue.variable1, name_source), " -- ", - _name(issue.variable2), + ModelAnalyzer._name(issue.variable2, name_source), ") with coefficient ", issue.coefficient, ) @@ -2067,15 +2257,16 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::LargeMatrixQuadraticCoefficient, + name_source, ) return print( io, "(Constraint -- Variable -- Variable): (", - _name(issue.ref), + ModelAnalyzer._name(issue.ref, name_source), " -- ", - _name(issue.variable1), + ModelAnalyzer._name(issue.variable1, name_source), " -- ", - _name(issue.variable2), + ModelAnalyzer._name(issue.variable2, name_source), ") with coefficient ", issue.coefficient, ) @@ -2084,6 +2275,7 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::NonconvexQuadraticObjective, + name_source, ) return ModelAnalyzer._summarize(io, issue) end @@ -2091,8 +2283,13 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::NonconvexQuadraticConstraint, + name_source, ) - return print(io, "Constraint: ", _name(issue.ref)) + return print( + io, + "Constraint: ", + ModelAnalyzer._name(issue.ref, name_source), + ) end function ModelAnalyzer.list_of_issues( @@ -2325,10 +2522,6 @@ end # printing helpers -function _name(ref) - return "$(ref)" -end - _print_value(x::Real) = Printf.@sprintf("%1.0e", x) function _stringify_bounds(bounds::Vector{Float64}) From cea96e6ec1bb37393b2ecfd3d48d5ab2dfe8b231 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Thu, 1 May 2025 02:40:43 -0300 Subject: [PATCH 04/34] fixes --- src/ModelAnalyzer.jl | 2 ++ src/jump.jl | 2 +- src/numerical.jl | 4 ++-- 3 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/ModelAnalyzer.jl b/src/ModelAnalyzer.jl index 7b5f496..1c11246 100644 --- a/src/ModelAnalyzer.jl +++ b/src/ModelAnalyzer.jl @@ -5,6 +5,8 @@ module ModelAnalyzer +import MathOptInterface as MOI + abstract type AbstractIssue end abstract type AbstractData end diff --git a/src/jump.jl b/src/jump.jl index 437f6b5..88440d6 100644 --- a/src/jump.jl +++ b/src/jump.jl @@ -1,4 +1,4 @@ -using JuMP +import JuMP # struct JuMPData{T<:AbstractData} <: ModelAnalyzer.AbstractData # data::T diff --git a/src/numerical.jl b/src/numerical.jl index 854a3ce..d9f1fb4 100644 --- a/src/numerical.jl +++ b/src/numerical.jl @@ -94,7 +94,7 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.VariableBoundAsConstraint struct VariableBoundAsConstraint <: AbstractNumericalIssue ref::MOI.ConstraintIndex end -ModelAnalyzer.constraint(issue::EmptyConstraint, model) = issue.ref +ModelAnalyzer.constraint(issue::VariableBoundAsConstraint, model) = issue.ref """ DenseConstraint <: AbstractNumericalIssue @@ -2277,7 +2277,7 @@ function ModelAnalyzer._verbose_summarize( issue::NonconvexQuadraticObjective, name_source, ) - return ModelAnalyzer._summarize(io, issue) + return ModelAnalyzer._summarize(io, issue, name_source) end function ModelAnalyzer._verbose_summarize( From 7b0c0a2b4a67bc1e72404fbfb4baf9c923af147c Mon Sep 17 00:00:00 2001 From: joaquimg Date: Thu, 1 May 2025 12:18:12 -0300 Subject: [PATCH 05/34] name_source -> model --- src/ModelAnalyzer.jl | 18 ++-- src/numerical.jl | 227 +++++++++++++++++-------------------------- 2 files changed, 96 insertions(+), 149 deletions(-) diff --git a/src/ModelAnalyzer.jl b/src/ModelAnalyzer.jl index 1c11246..b7a6498 100644 --- a/src/ModelAnalyzer.jl +++ b/src/ModelAnalyzer.jl @@ -27,7 +27,7 @@ See [`summarize`](@ref), [`list_of_issues`](@ref), and function analyze end """ - summarize([io::IO,] AbstractData; name_source = nothing, verbose = true, max_issues = typemax(Int), kwargs...) + summarize([io::IO,] AbstractData; model = nothing, verbose = true, max_issues = typemax(Int), kwargs...) Print a summary of the analysis results contained in `AbstractData` to the specified IO stream. If no IO stream is provided, it defaults to `stdout`. @@ -43,11 +43,11 @@ be a subtype of `AbstractIssue`). In the verbose case it will provide a text explaning the issue. In the non-verbose case it will provide just the issue name. - summarize([io::IO,] issue::AbstractIssue; name_source = nothing, verbose = true) + summarize([io::IO,] issue::AbstractIssue; model = nothing, verbose = true) This variant allows summarizing a single issue instance of type `AbstractIssue`. -The model tha led to the issue can be provided to `name_source`, it will be used -to generate the name of variables and constraints in the issue summary. +The model that led to the issue can be provided to `model`, it will be +used to generate the name of variables and constraints in the issue summary. """ function summarize end @@ -79,20 +79,20 @@ end function summarize( io::IO, issue::AbstractIssue; - name_source = nothing, + model = nothing, verbose = true, ) if verbose - return _verbose_summarize(io, issue, name_source) + return _verbose_summarize(io, issue, model) else - return _summarize(io, issue, name_source) + return _summarize(io, issue, model) end end function summarize( io::IO, issues::Vector{T}; - name_source = nothing, + model = nothing, verbose = true, max_issues = typemax(Int), ) where {T<:AbstractIssue} @@ -102,7 +102,7 @@ function summarize( print(io, "\n\n## List of issues\n\n") for issue in first(issues, max_issues) print(io, " * ") - summarize(io, issue, verbose = verbose, name_source = name_source) + summarize(io, issue, verbose = verbose, model = model) print(io, "\n") end return diff --git a/src/numerical.jl b/src/numerical.jl index d9f1fb4..72457b9 100644 --- a/src/numerical.jl +++ b/src/numerical.jl @@ -1827,109 +1827,80 @@ end function ModelAnalyzer._summarize( io::IO, issue::VariableNotInConstraints, - name_source, + model, ) - return print(io, ModelAnalyzer._name(issue.ref, name_source)) + return print(io, ModelAnalyzer._name(issue.ref, model)) end -function ModelAnalyzer._summarize(io::IO, issue::EmptyConstraint, name_source) - return print(io, ModelAnalyzer._name(issue.ref, name_source)) +function ModelAnalyzer._summarize(io::IO, issue::EmptyConstraint, model) + return print(io, ModelAnalyzer._name(issue.ref, model)) end function ModelAnalyzer._summarize( io::IO, issue::VariableBoundAsConstraint, - name_source, + model, ) - return print(io, ModelAnalyzer._name(issue.ref, name_source)) + return print(io, ModelAnalyzer._name(issue.ref, model)) end -function ModelAnalyzer._summarize(io::IO, issue::DenseConstraint, name_source) - return print( - io, - ModelAnalyzer._name(issue.ref, name_source), - " : ", - issue.nnz, - ) +function ModelAnalyzer._summarize(io::IO, issue::DenseConstraint, model) + return print(io, ModelAnalyzer._name(issue.ref, model), " : ", issue.nnz) end -function ModelAnalyzer._summarize( - io::IO, - issue::SmallMatrixCoefficient, - name_source, -) +function ModelAnalyzer._summarize(io::IO, issue::SmallMatrixCoefficient, model) return print( io, - ModelAnalyzer._name(issue.ref, name_source), + ModelAnalyzer._name(issue.ref, model), " -- ", - ModelAnalyzer._name(issue.variable, name_source), + ModelAnalyzer._name(issue.variable, model), " : ", issue.coefficient, ) end -function ModelAnalyzer._summarize( - io::IO, - issue::LargeMatrixCoefficient, - name_source, -) +function ModelAnalyzer._summarize(io::IO, issue::LargeMatrixCoefficient, model) return print( io, - ModelAnalyzer._name(issue.ref, name_source), + ModelAnalyzer._name(issue.ref, model), " -- ", - ModelAnalyzer._name(issue.variable, name_source), + ModelAnalyzer._name(issue.variable, model), " : ", issue.coefficient, ) end -function ModelAnalyzer._summarize( - io::IO, - issue::SmallBoundCoefficient, - name_source, -) +function ModelAnalyzer._summarize(io::IO, issue::SmallBoundCoefficient, model) return print( io, - ModelAnalyzer._name(issue.variable, name_source), + ModelAnalyzer._name(issue.variable, model), " : ", issue.coefficient, ) end -function ModelAnalyzer._summarize( - io::IO, - issue::LargeBoundCoefficient, - name_source, -) +function ModelAnalyzer._summarize(io::IO, issue::LargeBoundCoefficient, model) return print( io, - ModelAnalyzer._name(issue.variable, name_source), + ModelAnalyzer._name(issue.variable, model), " : ", issue.coefficient, ) end -function ModelAnalyzer._summarize( - io::IO, - issue::SmallRHSCoefficient, - name_source, -) +function ModelAnalyzer._summarize(io::IO, issue::SmallRHSCoefficient, model) return print( io, - ModelAnalyzer._name(issue.ref, name_source), + ModelAnalyzer._name(issue.ref, model), " : ", issue.coefficient, ) end -function ModelAnalyzer._summarize( - io::IO, - issue::LargeRHSCoefficient, - name_source, -) +function ModelAnalyzer._summarize(io::IO, issue::LargeRHSCoefficient, model) return print( io, - ModelAnalyzer._name(issue.ref, name_source), + ModelAnalyzer._name(issue.ref, model), " : ", issue.coefficient, ) @@ -1938,11 +1909,11 @@ end function ModelAnalyzer._summarize( io::IO, issue::SmallObjectiveCoefficient, - name_source, + model, ) return print( io, - ModelAnalyzer._name(issue.variable, name_source), + ModelAnalyzer._name(issue.variable, model), " : ", issue.coefficient, ) @@ -1951,11 +1922,11 @@ end function ModelAnalyzer._summarize( io::IO, issue::LargeObjectiveCoefficient, - name_source, + model, ) return print( io, - ModelAnalyzer._name(issue.variable, name_source), + ModelAnalyzer._name(issue.variable, model), " : ", issue.coefficient, ) @@ -1964,13 +1935,13 @@ end function ModelAnalyzer._summarize( io::IO, issue::SmallObjectiveQuadraticCoefficient, - name_source, + model, ) return print( io, - ModelAnalyzer._name(issue.variable1, name_source), + ModelAnalyzer._name(issue.variable1, model), " -- ", - ModelAnalyzer._name(issue.variable2, name_source), + ModelAnalyzer._name(issue.variable2, model), " : ", issue.coefficient, ) @@ -1979,13 +1950,13 @@ end function ModelAnalyzer._summarize( io::IO, issue::LargeObjectiveQuadraticCoefficient, - name_source, + model, ) return print( io, - ModelAnalyzer._name(issue.variable1, name_source), + ModelAnalyzer._name(issue.variable1, model), " -- ", - ModelAnalyzer._name(issue.variable2, name_source), + ModelAnalyzer._name(issue.variable2, model), " : ", issue.coefficient, ) @@ -1994,15 +1965,15 @@ end function ModelAnalyzer._summarize( io::IO, issue::SmallMatrixQuadraticCoefficient, - name_source, + model, ) return print( io, - ModelAnalyzer._name(issue.ref, name_source), + ModelAnalyzer._name(issue.ref, model), " -- ", - ModelAnalyzer._name(issue.variable1, name_source), + ModelAnalyzer._name(issue.variable1, model), " -- ", - ModelAnalyzer._name(issue.variable2, name_source), + ModelAnalyzer._name(issue.variable2, model), " : ", issue.coefficient, ) @@ -2011,77 +1982,57 @@ end function ModelAnalyzer._summarize( io::IO, issue::LargeMatrixQuadraticCoefficient, - name_source, + model, ) return print( io, - ModelAnalyzer._name(issue.ref, name_source), + ModelAnalyzer._name(issue.ref, model), " -- ", - ModelAnalyzer._name(issue.variable1, name_source), + ModelAnalyzer._name(issue.variable1, model), " -- ", - ModelAnalyzer._name(issue.variable2, name_source), + ModelAnalyzer._name(issue.variable2, model), " : ", issue.coefficient, ) end -function ModelAnalyzer._summarize( - io::IO, - ::NonconvexQuadraticObjective, - name_source, -) +function ModelAnalyzer._summarize(io::IO, ::NonconvexQuadraticObjective, model) return print(io, "Objective is Nonconvex quadratic") end function ModelAnalyzer._summarize( io::IO, issue::NonconvexQuadraticConstraint, - name_source, + model, ) - return print(io, ModelAnalyzer._name(issue.ref, name_source)) + return print(io, ModelAnalyzer._name(issue.ref, model)) end function ModelAnalyzer._verbose_summarize( io::IO, issue::VariableNotInConstraints, - name_source, + model, ) - return print(io, "Variable: ", ModelAnalyzer._name(issue.ref, name_source)) + return print(io, "Variable: ", ModelAnalyzer._name(issue.ref, model)) end -function ModelAnalyzer._verbose_summarize( - io::IO, - issue::EmptyConstraint, - name_source, -) - return print( - io, - "Constraint: ", - ModelAnalyzer._name(issue.ref, name_source), - ) +function ModelAnalyzer._verbose_summarize(io::IO, issue::EmptyConstraint, model) + return print(io, "Constraint: ", ModelAnalyzer._name(issue.ref, model)) end function ModelAnalyzer._verbose_summarize( io::IO, issue::VariableBoundAsConstraint, - name_source, + model, ) - return print( - io, - "Constraint: ", - ModelAnalyzer._name(issue.ref, name_source), - ) + return print(io, "Constraint: ", ModelAnalyzer._name(issue.ref, model)) end -function ModelAnalyzer._verbose_summarize( - io::IO, - issue::DenseConstraint, - name_source, -) +function ModelAnalyzer._verbose_summarize(io::IO, issue::DenseConstraint, model) return print( io, "Constraint: ", - ModelAnalyzer._name(issue.ref, name_source), + ModelAnalyzer._name(issue.ref, model), " with ", issue.nnz, " non zero coefficients", @@ -2091,14 +2042,14 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::SmallMatrixCoefficient, - name_source, + model, ) return print( io, "(Constraint -- Variable): (", - ModelAnalyzer._name(issue.ref, name_source), + ModelAnalyzer._name(issue.ref, model), " -- ", - ModelAnalyzer._name(issue.variable, name_source), + ModelAnalyzer._name(issue.variable, model), ") with coefficient ", issue.coefficient, ) @@ -2107,14 +2058,14 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::LargeMatrixCoefficient, - name_source, + model, ) return print( io, "(Constraint -- Variable): (", - ModelAnalyzer._name(issue.ref, name_source), + ModelAnalyzer._name(issue.ref, model), " -- ", - ModelAnalyzer._name(issue.variable, name_source), + ModelAnalyzer._name(issue.variable, model), ") with coefficient ", issue.coefficient, ) @@ -2123,12 +2074,12 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::SmallBoundCoefficient, - name_source, + model, ) return print( io, "Variable: ", - ModelAnalyzer._name(issue.variable, name_source), + ModelAnalyzer._name(issue.variable, model), " with bound ", issue.coefficient, ) @@ -2137,12 +2088,12 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::LargeBoundCoefficient, - name_source, + model, ) return print( io, "Variable: ", - ModelAnalyzer._name(issue.variable, name_source), + ModelAnalyzer._name(issue.variable, model), " with bound ", issue.coefficient, ) @@ -2151,12 +2102,12 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::SmallRHSCoefficient, - name_source, + model, ) return print( io, "Constraint: ", - ModelAnalyzer._name(issue.ref, name_source), + ModelAnalyzer._name(issue.ref, model), " with right-hand-side ", issue.coefficient, ) @@ -2165,12 +2116,12 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::LargeRHSCoefficient, - name_source, + model, ) return print( io, "Constraint: ", - ModelAnalyzer._name(issue.ref, name_source), + ModelAnalyzer._name(issue.ref, model), " with right-hand-side ", issue.coefficient, ) @@ -2179,12 +2130,12 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::SmallObjectiveCoefficient, - name_source, + model, ) return print( io, "Variable: ", - ModelAnalyzer._name(issue.variable, name_source), + ModelAnalyzer._name(issue.variable, model), " with coefficient ", issue.coefficient, ) @@ -2193,12 +2144,12 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::LargeObjectiveCoefficient, - name_source, + model, ) return print( io, "Variable: ", - ModelAnalyzer._name(issue.variable, name_source), + ModelAnalyzer._name(issue.variable, model), " with coefficient ", issue.coefficient, ) @@ -2207,14 +2158,14 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::SmallObjectiveQuadraticCoefficient, - name_source, + model, ) return print( io, "(Variable -- Variable): (", - ModelAnalyzer._name(issue.variable1, name_source), + ModelAnalyzer._name(issue.variable1, model), " -- ", - ModelAnalyzer._name(issue.variable2, name_source), + ModelAnalyzer._name(issue.variable2, model), ") with coefficient ", issue.coefficient, ) @@ -2223,14 +2174,14 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::LargeObjectiveQuadraticCoefficient, - name_source, + model, ) return print( io, "(Variable -- Variable): (", - ModelAnalyzer._name(issue.variable1, name_source), + ModelAnalyzer._name(issue.variable1, model), " -- ", - ModelAnalyzer._name(issue.variable2, name_source), + ModelAnalyzer._name(issue.variable2, model), ") with coefficient ", issue.coefficient, ) @@ -2239,16 +2190,16 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::SmallMatrixQuadraticCoefficient, - name_source, + model, ) return print( io, "(Constraint -- Variable -- Variable): (", - ModelAnalyzer._name(issue.ref, name_source), + ModelAnalyzer._name(issue.ref, model), " -- ", - ModelAnalyzer._name(issue.variable1, name_source), + ModelAnalyzer._name(issue.variable1, model), " -- ", - ModelAnalyzer._name(issue.variable2, name_source), + ModelAnalyzer._name(issue.variable2, model), ") with coefficient ", issue.coefficient, ) @@ -2257,16 +2208,16 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::LargeMatrixQuadraticCoefficient, - name_source, + model, ) return print( io, "(Constraint -- Variable -- Variable): (", - ModelAnalyzer._name(issue.ref, name_source), + ModelAnalyzer._name(issue.ref, model), " -- ", - ModelAnalyzer._name(issue.variable1, name_source), + ModelAnalyzer._name(issue.variable1, model), " -- ", - ModelAnalyzer._name(issue.variable2, name_source), + ModelAnalyzer._name(issue.variable2, model), ") with coefficient ", issue.coefficient, ) @@ -2275,21 +2226,17 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::NonconvexQuadraticObjective, - name_source, + model, ) - return ModelAnalyzer._summarize(io, issue, name_source) + return ModelAnalyzer._summarize(io, issue, model) end function ModelAnalyzer._verbose_summarize( io::IO, issue::NonconvexQuadraticConstraint, - name_source, + model, ) - return print( - io, - "Constraint: ", - ModelAnalyzer._name(issue.ref, name_source), - ) + return print(io, "Constraint: ", ModelAnalyzer._name(issue.ref, model)) end function ModelAnalyzer.list_of_issues( From 2c79007973c713d24986eff563fa7359039b3fa9 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Thu, 1 May 2025 12:31:24 -0300 Subject: [PATCH 06/34] suggestions --- src/ModelAnalyzer.jl | 7 ++++++- src/feasibility.jl | 17 +++++++++------- src/infeasibility.jl | 8 ++++---- src/jump.jl | 46 +++++++++++++++++++++----------------------- 4 files changed, 42 insertions(+), 36 deletions(-) diff --git a/src/ModelAnalyzer.jl b/src/ModelAnalyzer.jl index b7a6498..3162277 100644 --- a/src/ModelAnalyzer.jl +++ b/src/ModelAnalyzer.jl @@ -14,7 +14,7 @@ abstract type AbstractData end abstract type AbstractAnalyzer end """ - analyze(analyzer::AbstractAnalyzer, model::JuMP.Model; kwargs...) + analyze(analyzer::AbstractAnalyzer, model::JuMP.GenericModel; kwargs...) Analyze a JuMP model using the specified analyzer. Depending on the analyzer, this keyword arguments might vary. @@ -122,6 +122,7 @@ some numerical issues, it can be the coefficient value. function value(issue::AbstractIssue, ::Nothing) return value(issue) end + function value(issue::AbstractIssue, ::MOI.ModelLike) return value(issue) end @@ -134,6 +135,7 @@ Return the variable associated to a particular issue. function variable(issue::AbstractIssue, ::Nothing) return variable(issue) end + function variable(issue::AbstractIssue, ::MOI.ModelLike) return variable(issue) end @@ -146,6 +148,7 @@ Return the variables associated to a particular issue. function variables(issue::AbstractIssue, ::Nothing) return variables(issue) end + function variables(issue::AbstractIssue, ::MOI.ModelLike) return variables(issue) end @@ -158,11 +161,13 @@ Return the constraint associated to a particular issue. function constraint(issue::AbstractIssue, ::Nothing) return constraint(issue) end + function constraint(issue::AbstractIssue, ::MOI.ModelLike) return constraint(issue) end function _verbose_summarize end + function _summarize end function _name(ref::MOI.VariableIndex, model::MOI.ModelLike) diff --git a/src/feasibility.jl b/src/feasibility.jl index cb38584..f163e49 100644 --- a/src/feasibility.jl +++ b/src/feasibility.jl @@ -85,7 +85,7 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Feasibility.DualViolation) ``` """ struct DualViolation <: AbstractFeasibilityIssue - ref::Union{JuMP.ConstraintRef,JuMP.VariableRef} + ref::Union{JuMP.ConstraintRef,JuMP.GenericVariableRef} violation::Float64 end @@ -651,7 +651,7 @@ function _name(ref::JuMP.ConstraintRef) return JuMP.name(ref) end -function _name(ref::JuMP.VariableRef) +function _name(ref::JuMP.GenericVariableRef) return JuMP.name(ref) end @@ -1062,12 +1062,12 @@ function _dual_point_to_dual_model_ref( if haskey(primal_con_dual_var, moi_con) vec_vars = primal_con_dual_var[moi_con] for (i, moi_var) in enumerate(vec_vars) - jump_var = JuMP.VariableRef(dual_model, moi_var) + jump_var = JuMP.GenericVariableRef{T}(dual_model, moi_var) dual_point[jump_var] = val[i] end elseif haskey(primal_con_dual_convar, moi_con) moi_convar = primal_con_dual_convar[moi_con] - jump_var = JuMP.VariableRef( + jump_var = JuMP.GenericVariableRef{T}( dual_model, MOI.VariableIndex(moi_convar.value), ) @@ -1105,7 +1105,10 @@ function _fix_ret( primal_model::JuMP.GenericModel{T}, dual_con_primal_all, ) where {T} - ret = Dict{Union{JuMP.ConstraintRef,JuMP.VariableRef},Union{T,Vector{T}}}() + ret = Dict{ + Union{JuMP.ConstraintRef,JuMP.GenericVariableRef{T}}, + Union{T,Vector{T}}, + }() for (jump_dual_con, val) in pre_ret # v is a variable in the dual jump model # we need the associated cosntraint in the primal jump model @@ -1115,7 +1118,7 @@ function _fix_ret( # variable in the dual model # constraint in the primal model jump_primal_var = - JuMP.VariableRef(primal_model, moi_primal_something) + JuMP.GenericVariableRef{T}(primal_model, moi_primal_something) # ret[jump_primal_var] = T[val] ret[jump_primal_var] = val else @@ -1144,7 +1147,7 @@ function _dualize2( if mode == JuMP.MANUAL error("Dualization does not support solvers in $(mode) mode") end - dual_model = JuMP.Model() + dual_model = JuMP.GenericModel() dual_problem = Dualization.DualProblem(JuMP.backend(dual_model)) Dualization.dualize(JuMP.backend(model), dual_problem; kwargs...) Dualization._fill_obj_dict_with_variables!(dual_model) diff --git a/src/infeasibility.jl b/src/infeasibility.jl index 079c8e1..4ffd6ee 100644 --- a/src/infeasibility.jl +++ b/src/infeasibility.jl @@ -50,7 +50,7 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Infeasibility.InfeasibleBounds) ```` """ struct InfeasibleBounds{T} <: AbstractInfeasibilitylIssue - variable::JuMP.VariableRef + variable::JuMP.GenericVariableRef{T} lb::T ub::T end @@ -71,7 +71,7 @@ julia> ModelAnalyzer.summarize( ``` """ struct InfeasibleIntegrality{T} <: AbstractInfeasibilitylIssue - variable::JuMP.VariableRef + variable::JuMP.GenericVariableRef{T} lb::T ub::T set::Union{MOI.Integer,MOI.ZeroOne}#, MOI.Semicontinuous{T}, MOI.Semiinteger{T}} @@ -147,7 +147,7 @@ function ModelAnalyzer.analyze( ) where {T} out = Data() - variables = Dict{JuMP.VariableRef,Interval{T}}() + variables = Dict{JuMP.GenericVariableRef{T},Interval{T}}() # first layer of infeasibility analysis is bounds consistency bounds_consistent = true @@ -196,7 +196,7 @@ function ModelAnalyzer.analyze( # second layer of infeasibility analysis is constraint range analysis range_consistent = true for (F, S) in JuMP.list_of_constraint_types(model) - F != JuMP.GenericAffExpr{T,JuMP.VariableRef} && continue + F != JuMP.GenericAffExpr{T,JuMP.GenericVariableRef{T}} && continue # TODO: handle quadratics !(S in (MOI.EqualTo{T}, MOI.LessThan{T}, MOI.GreaterThan{T})) && continue diff --git a/src/jump.jl b/src/jump.jl index 88440d6..9f60f9a 100644 --- a/src/jump.jl +++ b/src/jump.jl @@ -1,18 +1,13 @@ -import JuMP - -# struct JuMPData{T<:AbstractData} <: ModelAnalyzer.AbstractData -# data::T -# model::JuMP.Model -# end +# 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. -# struct JuMPIssue{T<:AbstractIssue} <: ModelAnalyzer.AbstractIssue -# issue::T -# model::JuMP.Model -# end +import JuMP function ModelAnalyzer.analyze( analyzer::T, - model::JuMP.Model; + model::JuMP.GenericModel; kwargs..., ) where {T<:ModelAnalyzer.AbstractAnalyzer} moi_model = JuMP.backend(model) @@ -21,8 +16,11 @@ function ModelAnalyzer.analyze( return result end -function ModelAnalyzer._name(ref::MOI.VariableIndex, model::JuMP.Model) - jump_ref = JuMP.VariableRef(model, ref) +function ModelAnalyzer._name( + ref::MOI.VariableIndex, + model::JuMP.GenericModel{T}, +) where {T} + jump_ref = JuMP.GenericVariableRef{T}(model, ref) name = JuMP.name(jump_ref) if !isempty(name) return name @@ -30,7 +28,7 @@ function ModelAnalyzer._name(ref::MOI.VariableIndex, model::JuMP.Model) return "$jump_ref" end -function ModelAnalyzer._name(ref::MOI.ConstraintIndex, model::JuMP.Model) +function ModelAnalyzer._name(ref::MOI.ConstraintIndex, model::JuMP.GenericModel) jump_ref = JuMP.constraint_ref_with_index(model, ref) name = JuMP.name(jump_ref) if !isempty(name) @@ -40,39 +38,39 @@ function ModelAnalyzer._name(ref::MOI.ConstraintIndex, model::JuMP.Model) end """ - variable(issue::ModelAnalyzer.AbstractIssue, model::JuMP.Model) + variable(issue::ModelAnalyzer.AbstractIssue, model::JuMP.GenericModel) Return the **JuMP** variable reference associated to a particular issue. """ function ModelAnalyzer.variable( issue::ModelAnalyzer.AbstractIssue, - model::JuMP.Model, -) + model::JuMP.GenericModel{T}, +) where {T} ref = ModelAnalyzer.variable(issue) - return JuMP.VariableRef(model, ref) + return JuMP.GenericVariableRef{T}(model, ref) end """ - variables(issue::ModelAnalyzer.AbstractIssue, model::JuMP.Model) + variables(issue::ModelAnalyzer.AbstractIssue, model::JuMP.GenericModel) Return the **JuMP** variable references associated to a particular issue. """ function ModelAnalyzer.variables( issue::ModelAnalyzer.AbstractIssue, - model::JuMP.Model, -) + model::JuMP.GenericModel{T}, +) where {T} refs = ModelAnalyzer.variables(issue) - return JuMP.VariableRef.(model, refs) + return JuMP.GenericVariableRef{T}.(model, refs) end """ - constraint(issue::ModelAnalyzer.AbstractIssue, model::JuMP.Model) + constraint(issue::ModelAnalyzer.AbstractIssue, model::JuMP.GenericModel) Return the **JuMP** constraint reference associated to a particular issue. """ function ModelAnalyzer.constraint( issue::ModelAnalyzer.AbstractIssue, - model::JuMP.Model, + model::JuMP.GenericModel, ) ref = ModelAnalyzer.constraint(issue) return JuMP.constraint_ref_with_index(model, ref) From 9241ee16b7d3aed016ba2e4d715da4c22270f770 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Thu, 1 May 2025 12:32:17 -0300 Subject: [PATCH 07/34] small fix --- src/jump.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/jump.jl b/src/jump.jl index 9f60f9a..f3a6e78 100644 --- a/src/jump.jl +++ b/src/jump.jl @@ -12,7 +12,6 @@ function ModelAnalyzer.analyze( ) where {T<:ModelAnalyzer.AbstractAnalyzer} moi_model = JuMP.backend(model) result = ModelAnalyzer.analyze(analyzer, moi_model; kwargs...) - # return JuMPData(result, model) return result end From ed420752dcf452c6d3ae0798ed74d3dae5421a73 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Wed, 7 May 2025 02:50:36 -0300 Subject: [PATCH 08/34] feasibility almost ready --- Project.toml | 9 +- docs/Project.toml | 6 +- docs/make.jl | 2 +- docs/src/feasibility.md | 3 +- .../ModelAnalyzerJuMPExt.jl | 6 + src/ModelAnalyzer.jl | 7 +- src/feasibility.jl | 880 ++++++++++-------- src/numerical.jl | 2 + test/feasibility.jl | 617 +++++++++--- test/runtests.jl | 5 +- 10 files changed, 975 insertions(+), 562 deletions(-) rename src/jump.jl => ext/ModelAnalyzerJuMPExt/ModelAnalyzerJuMPExt.jl (94%) diff --git a/Project.toml b/Project.toml index b858c6f..46d6f91 100644 --- a/Project.toml +++ b/Project.toml @@ -4,12 +4,17 @@ version = "0.1.0" [deps] Dualization = "191a621a-6537-11e9-281d-650236a99e60" -JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +[weakdeps] +JuMP = "4076af6c-e467-56ae-b986-b466b2749572" + +[extensions] +ModelAnalyzerJuMPExt = "JuMP" + [compat] -Dualization = "0.5.9" +Dualization = "0.6.0" JuMP = "1.24.0" MathOptInterface = "1.37.0" \ No newline at end of file diff --git a/docs/Project.toml b/docs/Project.toml index fd952a6..32b49dc 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,3 +1,7 @@ [deps] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -ModelAnalyzer = "d1179b25-476b-425c-b826-c7787f0fff83" \ No newline at end of file +ModelAnalyzer = "d1179b25-476b-425c-b826-c7787f0fff83" +JuMP = "4076af6c-e467-56ae-b986-b466b2749572" + +[compat] +JuMP = "1.24.0" \ No newline at end of file diff --git a/docs/make.jl b/docs/make.jl index b280ffc..bf2d820 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,4 +1,4 @@ -using Documenter, ModelAnalyzer +using Documenter, ModelAnalyzer, JuMP makedocs(; sitename = "ModelAnalyzer.jl documentation") diff --git a/docs/src/feasibility.md b/docs/src/feasibility.md index 8ce2b76..3fac140 100644 --- a/docs/src/feasibility.md +++ b/docs/src/feasibility.md @@ -17,7 +17,8 @@ Specifically, the possible issues are: ```@docs ModelAnalyzer.Feasibility.PrimalViolation -ModelAnalyzer.Feasibility.DualViolation +ModelAnalyzer.Feasibility.DualConstraintViolation +ModelAnalyzer.Feasibility.DualConstrainedVariableViolation ModelAnalyzer.Feasibility.ComplemetarityViolation ModelAnalyzer.Feasibility.DualObjectiveMismatch ModelAnalyzer.Feasibility.PrimalObjectiveMismatch diff --git a/src/jump.jl b/ext/ModelAnalyzerJuMPExt/ModelAnalyzerJuMPExt.jl similarity index 94% rename from src/jump.jl rename to ext/ModelAnalyzerJuMPExt/ModelAnalyzerJuMPExt.jl index f3a6e78..a5eb50a 100644 --- a/src/jump.jl +++ b/ext/ModelAnalyzerJuMPExt/ModelAnalyzerJuMPExt.jl @@ -3,7 +3,11 @@ # 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. +module ModelAnalyzerJuMPExt + +import ModelAnalyzer import JuMP +import MathOptInterface as MOI function ModelAnalyzer.analyze( analyzer::T, @@ -74,3 +78,5 @@ function ModelAnalyzer.constraint( ref = ModelAnalyzer.constraint(issue) return JuMP.constraint_ref_with_index(model, ref) end + +end # module ModelAnalyzerJuMPExt diff --git a/src/ModelAnalyzer.jl b/src/ModelAnalyzer.jl index c5fe576..a28aeb4 100644 --- a/src/ModelAnalyzer.jl +++ b/src/ModelAnalyzer.jl @@ -31,6 +31,8 @@ function analyze end Print a summary of the analysis results contained in `AbstractData` to the specified IO stream. If no IO stream is provided, it defaults to `stdout`. +The model that led to the issue can be provided to `model`, it will be +used to generate the name of variables and constraints in the issue summary. The `verbose` flag controls whether to print detailed information about each issue (if `true`) or a concise summary (if `false`). The `max_issues` argument controls the maximum number of issues to display in the summary. If there are @@ -202,7 +204,6 @@ end include("numerical.jl") include("feasibility.jl") -include("infeasibility.jl") +# include("infeasibility.jl") -include("jump.jl") -end +end # module ModelAnalyzer diff --git a/src/feasibility.jl b/src/feasibility.jl index 403b3be..280823a 100644 --- a/src/feasibility.jl +++ b/src/feasibility.jl @@ -3,24 +3,17 @@ # 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. -# TODO -# 1 - JuMP: primal_feasibility_report should have a typed error for not found stuff so we can capture -# 2 - Dualization: should consider and option to dont treat @variable(m, x >= 0) differently from @variable(m, x >= 1) -# 3 - Dualization: JuMP model dualization should hold the primal dual map, maybe a JuMP converted version -# 4 - Dualization: Primal dual map could work with a getindex for simpler usage - module Feasibility import ModelAnalyzer import Dualization -import JuMP -import JuMP.MOI as MOI +import MathOptInterface as MOI import Printf """ Analyzer() <: ModelAnalyzer.AbstractAnalyzer -The `Analyzer` type is used to perform feasibility analysis on a JuMP model. +The `Analyzer` type is used to perform feasibility analysis on a model. ## Example @@ -53,7 +46,7 @@ struct Analyzer <: ModelAnalyzer.AbstractAnalyzer end """ AbstractFeasibilityIssue <: AbstractNumericalIssue -Abstract type for feasibility issues found during the analysis of a JuMP model. +Abstract type for feasibility issues found during the analysis of a model. """ abstract type AbstractFeasibilityIssue <: ModelAnalyzer.AbstractIssue end @@ -69,23 +62,42 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Feasibility.PrimalViolation) ``` """ struct PrimalViolation <: AbstractFeasibilityIssue - ref::JuMP.ConstraintRef + ref::MOI.ConstraintIndex + violation::Float64 +end + +""" + DualConstraintViolation <: AbstractFeasibilityIssue + +The `DualConstraintViolation` issue is identified when a dual constraint has a +value that is not within the dual constraint's set. +This dual constraint corresponds to a primal variable. + +For more information, run: +```julia +julia> ModelAnalyzer.summarize(ModelAnalyzer.Feasibility.DualConstraintViolation) +``` +""" +struct DualConstraintViolation <: AbstractFeasibilityIssue + ref::MOI.VariableIndex violation::Float64 end """ - DualViolation <: AbstractFeasibilityIssue + DualConstrainedVariableViolation <: AbstractFeasibilityIssue -The `DualViolation` issue is identified when a constraint has a dual value +The `DualConstrainedVariableViolation` issue is identified when a dual +constraint, which is a constrained varaible constraint, has a value that is not within the dual constraint's set. +This dual constraint corresponds to a primal constraint. For more information, run: ```julia -julia> ModelAnalyzer.summarize(ModelAnalyzer.Feasibility.DualViolation) +julia> ModelAnalyzer.summarize(ModelAnalyzer.Feasibility.DualConstrainedVariableViolation) ``` """ -struct DualViolation <: AbstractFeasibilityIssue - ref::Union{JuMP.ConstraintRef,JuMP.GenericVariableRef} +struct DualConstrainedVariableViolation <: AbstractFeasibilityIssue + ref::MOI.ConstraintIndex violation::Float64 end @@ -103,7 +115,7 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Feasibility.ComplemetarityViolation ``` """ struct ComplemetarityViolation <: AbstractFeasibilityIssue - ref::JuMP.ConstraintRef + ref::MOI.ConstraintIndex violation::Float64 end @@ -179,7 +191,7 @@ end Data The `Data` structure holds the results of the feasibility analysis performed -by the `ModelAnalyzer.analyze` function for a JuMP model. It contains +by the `ModelAnalyzer.analyze` function for a model. It contains the configuration used for the analysis, the primal and dual points, and the lists of various feasibility issues found during the analysis. """ @@ -192,7 +204,9 @@ Base.@kwdef mutable struct Data <: ModelAnalyzer.AbstractData dual_check::Bool # analysis results primal::Vector{PrimalViolation} = PrimalViolation[] - dual::Vector{DualViolation} = DualViolation[] + dual::Vector{DualConstraintViolation} = DualConstraintViolation[] + dual_convar::Vector{DualConstrainedVariableViolation} = + DualConstrainedVariableViolation[] complementarity::Vector{ComplemetarityViolation} = ComplemetarityViolation[] # objective analysis dual_objective_mismatch::Vector{DualObjectiveMismatch} = @@ -208,8 +222,15 @@ function ModelAnalyzer._summarize(io::IO, ::Type{PrimalViolation}) return print(io, "# PrimalViolation") end -function ModelAnalyzer._summarize(io::IO, ::Type{DualViolation}) - return print(io, "# DualViolation") +function ModelAnalyzer._summarize(io::IO, ::Type{DualConstraintViolation}) + return print(io, "# DualConstraintViolation") +end + +function ModelAnalyzer._summarize( + io::IO, + ::Type{DualConstrainedVariableViolation}, +) + return print(io, "# DualConstrainedVariableViolation") end function ModelAnalyzer._summarize(io::IO, ::Type{ComplemetarityViolation}) @@ -271,15 +292,18 @@ function ModelAnalyzer._verbose_summarize(io::IO, ::Type{PrimalViolation}) ) end -function ModelAnalyzer._verbose_summarize(io::IO, ::Type{DualViolation}) +function ModelAnalyzer._verbose_summarize( + io::IO, + ::Type{DualConstraintViolation}, +) return print( io, """ - # DualViolation + # DualConstraintViolation ## What - A `DualViolation` issue is identified when a constraint has + A `DualConstraintViolation` issue is identified when a constraint has a dual value that is not within the dual constraint's set. ## Why @@ -309,6 +333,48 @@ function ModelAnalyzer._verbose_summarize(io::IO, ::Type{DualViolation}) ) end +function ModelAnalyzer._verbose_summarize( + io::IO, + ::Type{DualConstrainedVariableViolation}, +) + return print( + io, + """ + # DualConstrainedVariableViolation + + ## What + + A `DualConstrainedVariableViolation` issue is identified when a dual + constraint, which is a constrained varaible constraint, has a value + that is not within the dual constraint's set. + + ## Why + + This can happen due to a few reasons: + - The solver did not converge. + - The model is infeasible and the solver converged to an + infeasible point. + - The solver converged to a low accuracy solution, which might + happen due to transformations in the the model presolve or + due to numerical issues. + + ## How to fix + + Check the solver convergence log and the solver status. If the + solver did not converge, you might want to try alternative + solvers or adjust the solver options. If the solver converged + to an infeasible point, you might want to check the model + constraints and bounds. If the solver converged to a low + accuracy solution, you might want to adjust the solver options + or the model presolve. + + ## More information + + No extra information for this issue. + """, + ) +end + function ModelAnalyzer._verbose_summarize( io::IO, ::Type{ComplemetarityViolation}, @@ -489,78 +555,123 @@ function ModelAnalyzer._verbose_summarize( ) end -function ModelAnalyzer._summarize(io::IO, issue::PrimalViolation) - return print(io, _name(issue.ref), " : ", issue.violation) +function ModelAnalyzer._summarize(io::IO, issue::PrimalViolation, model) + return print( + io, + ModelAnalyzer._name(issue.ref, model), + " : ", + issue.violation, + ) +end + +function ModelAnalyzer._summarize(io::IO, issue::DualConstraintViolation, model) + return print( + io, + ModelAnalyzer._name(issue.ref, model), + " : ", + issue.violation, + ) end -function ModelAnalyzer._summarize(io::IO, issue::DualViolation) - return print(io, _name(issue.ref), " : ", issue.violation) +function ModelAnalyzer._summarize( + io::IO, + issue::DualConstrainedVariableViolation, + model, +) + return print( + io, + ModelAnalyzer._name(issue.ref, model), + " : ", + issue.violation, + ) end -function ModelAnalyzer._summarize(io::IO, issue::ComplemetarityViolation) - return print(io, _name(issue.ref), " : ", issue.violation) +function ModelAnalyzer._summarize(io::IO, issue::ComplemetarityViolation, model) + return print( + io, + ModelAnalyzer._name(issue.ref, model), + " : ", + issue.violation, + ) end -function ModelAnalyzer._summarize(io::IO, issue::DualObjectiveMismatch) - return ModelAnalyzer._verbose_summarize(io, issue) +function ModelAnalyzer._summarize(io::IO, issue::DualObjectiveMismatch, model) + return ModelAnalyzer._verbose_summarize(io, issue, model) end -function ModelAnalyzer._summarize(io::IO, issue::PrimalObjectiveMismatch) - return ModelAnalyzer._verbose_summarize(io, issue) +function ModelAnalyzer._summarize(io::IO, issue::PrimalObjectiveMismatch, model) + return ModelAnalyzer._verbose_summarize(io, issue, model) end -function ModelAnalyzer._summarize(io::IO, issue::PrimalDualMismatch) - return ModelAnalyzer._verbose_summarize(io, issue) +function ModelAnalyzer._summarize(io::IO, issue::PrimalDualMismatch, model) + return ModelAnalyzer._verbose_summarize(io, issue, model) end -function ModelAnalyzer._summarize(io::IO, issue::PrimalDualSolverMismatch) - return ModelAnalyzer._verbose_summarize(io, issue) +function ModelAnalyzer._summarize( + io::IO, + issue::PrimalDualSolverMismatch, + model, +) + return ModelAnalyzer._verbose_summarize(io, issue, model) end -function ModelAnalyzer._verbose_summarize(io::IO, issue::PrimalViolation) +function ModelAnalyzer._verbose_summarize(io::IO, issue::PrimalViolation, model) return print( io, "Constraint ", - _name(issue.ref), - " has violation ", + ModelAnalyzer._name(issue.ref, model), + " has primal violation ", issue.violation, ) end -function ModelAnalyzer._verbose_summarize(io::IO, issue::DualViolation) - if issue.ref isa JuMP.ConstraintRef - return print( - io, - "Constraint ", - _name(issue.ref), - " has violation ", - issue.violation, - ) - else - return print( - io, - "Variable ", - _name(issue.ref), - " has violation ", - issue.violation, - ) - end +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::DualConstraintViolation, + model, +) + return print( + io, + "Variables ", + ModelAnalyzer._name.(issue.ref, model), + " have dual violation ", + issue.violation, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::DualConstrainedVariableViolation, + model, +) + return print( + io, + "Constraint ", + ModelAnalyzer._name(issue.ref, model), + " has dual violation ", + issue.violation, + ) end function ModelAnalyzer._verbose_summarize( io::IO, issue::ComplemetarityViolation, + model, ) return print( io, "Constraint ", - _name(issue.ref), - " has violation ", + ModelAnalyzer._name(issue.ref, model), + " has complementarty violation ", issue.violation, ) end -function ModelAnalyzer._verbose_summarize(io::IO, issue::DualObjectiveMismatch) +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::DualObjectiveMismatch, + model, +) return print( io, "Dual objective mismatch: ", @@ -574,6 +685,7 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::PrimalObjectiveMismatch, + model, ) return print( io, @@ -585,7 +697,11 @@ function ModelAnalyzer._verbose_summarize( ) end -function ModelAnalyzer._verbose_summarize(io::IO, issue::PrimalDualMismatch) +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::PrimalDualMismatch, + model, +) return print( io, "Primal dual mismatch: ", @@ -599,6 +715,7 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::PrimalDualSolverMismatch, + model, ) return print( io, @@ -614,10 +731,20 @@ function ModelAnalyzer.list_of_issues(data::Data, ::Type{PrimalViolation}) return data.primal end -function ModelAnalyzer.list_of_issues(data::Data, ::Type{DualViolation}) +function ModelAnalyzer.list_of_issues( + data::Data, + ::Type{DualConstraintViolation}, +) return data.dual end +function ModelAnalyzer.list_of_issues( + data::Data, + ::Type{DualConstrainedVariableViolation}, +) + return data.dual_convar +end + function ModelAnalyzer.list_of_issues( data::Data, ::Type{ComplemetarityViolation}, @@ -647,19 +774,12 @@ function ModelAnalyzer.list_of_issues( return data.primal_dual_solver_mismatch end -function _name(ref::JuMP.ConstraintRef) - return JuMP.name(ref) -end - -function _name(ref::JuMP.GenericVariableRef) - return JuMP.name(ref) -end - function ModelAnalyzer.list_of_issue_types(data::Data) ret = Type[] for type in ( PrimalViolation, - DualViolation, + DualConstraintViolation, + DualConstrainedVariableViolation, ComplemetarityViolation, DualObjectiveMismatch, PrimalObjectiveMismatch, @@ -684,6 +804,7 @@ end function ModelAnalyzer.summarize( io::IO, data::Data; + model = nothing, verbose = true, max_issues = ModelAnalyzer.DEFAULT_MAX_ISSUES, configurations = true, @@ -701,6 +822,7 @@ function ModelAnalyzer.summarize( ModelAnalyzer.summarize( io, issues, + model = model, verbose = verbose, max_issues = max_issues, ) @@ -719,13 +841,24 @@ end function ModelAnalyzer.analyze( ::Analyzer, - model::JuMP.GenericModel; + model::MOI.ModelLike; primal_point = nothing, dual_point = nothing, atol::Float64 = 1e-6, skip_missing::Bool = false, dual_check = true, ) + can_dualize = false + if dual_check + can_dualize = _can_dualize(model) + if !can_dualize + println( + "The model cannot be dualized. Automatically setting `dual_check = false`.", + ) + dual_check = false + end + end + data = Data( primal_point = primal_point, dual_point = dual_point, @@ -735,10 +868,8 @@ function ModelAnalyzer.analyze( ) if data.primal_point === nothing - if !( - JuMP.primal_status(model) in - (MOI.FEASIBLE_POINT, MOI.NEARLY_FEASIBLE_POINT) - ) + primal_status = MOI.get(model, MOI.PrimalStatus()) + if !(primal_status in (MOI.FEASIBLE_POINT, MOI.NEARLY_FEASIBLE_POINT)) error( "No primal solution is available. You must provide a point at " * "which to check feasibility.", @@ -747,12 +878,9 @@ function ModelAnalyzer.analyze( data.primal_point = _last_primal_solution(model) end - can_dualize = _can_dualize(model) - if data.dual_point === nothing && can_dualize && dual_check - if !( - JuMP.dual_status(model) in - (MOI.FEASIBLE_POINT, MOI.NEARLY_FEASIBLE_POINT) - ) + if data.dual_point === nothing && dual_check + dual_status = MOI.get(model, MOI.DualStatus()) + if !(dual_status in (MOI.FEASIBLE_POINT, MOI.NEARLY_FEASIBLE_POINT)) error( "No dual solution is available. You must provide a point at " * "which to check feasibility. Or set dual_check = false.", @@ -761,20 +889,18 @@ function ModelAnalyzer.analyze( data.dual_point = _last_dual_solution(model) end - _dual_model = if can_dualize && dual_check - _dualize2(model) - else - nothing - end - _analyze_primal!(model, data) - if can_dualize && dual_check - _analyze_dual!(model, _dual_model, data) - end - if data.dual_point !== nothing + _dual_model = nothing + _map = nothing + if dual_check + dual_problem = + Dualization.dualize(model, consider_constrained_variables = false) + _dual_model = dual_problem.dual_model + _map = dual_problem.primal_dual_map + _analyze_dual!(model, _dual_model, _map, data) _analyze_complementarity!(model, data) end - _analyze_objectives!(model, _dual_model, data) + _analyze_objectives!(model, _dual_model, _map, data) sort!(data.primal, by = x -> abs(x.violation)) sort!(data.dual, by = x -> abs(x.violation)) sort!(data.complementarity, by = x -> abs(x.violation)) @@ -782,44 +908,203 @@ function ModelAnalyzer.analyze( end function _analyze_primal!(model, data) - dict = JuMP.primal_feasibility_report( - model, - data.primal_point; - atol = data.atol, - skip_missing = data.skip_missing, - ) - for (ref, violation) in dict - push!(data.primal, PrimalViolation(ref, violation)) + types = MOI.get(model, MOI.ListOfConstraintTypesPresent()) + for (F, S) in types + list = MOI.get(model, MOI.ListOfConstraintIndices{F,S}()) + for con in list + func = MOI.get(model, MOI.ConstraintFunction(), con) + failed = false + val = MOI.Utilities.eval_variables(model, func) do var_idx + if !haskey(data.primal_point, var_idx) + if data.skip_missing + failed = true + return NaN # nothing + else + error( + "Missing variable in primal point: $var_idx. " * + "Set skip_missing = true to ignore this error.", + ) + end + end + return data.primal_point[var_idx] + end + if failed + continue + end + set = MOI.get(model, MOI.ConstraintSet(), con) + dist = MOI.Utilities.distance_to_set(val, set) + if dist > data.atol + push!(data.primal, PrimalViolation(con, dist)) + end + end end return end -function _analyze_dual!(model, _dual_model, data) - dict = dual_feasibility_report( - model, - data.dual_point; - atol = data.atol, - skip_missing = data.skip_missing, - _dual_model = _dual_model, - ) - for (ref, violation) in dict - push!(data.dual, DualViolation(ref, violation)) +function _dual_point_to_dual_model_ref( + primal_model, + map::Dualization.PrimalDualMap, + dual_point, +) + new_dual_point = Dict{MOI.VariableIndex,Number}() + dual_var_to_primal_con = Dict{MOI.VariableIndex,MOI.ConstraintIndex}() + dual_con_to_primal_con = Dict{MOI.ConstraintIndex,MOI.ConstraintIndex}() + for (primal_con, val) in dual_point + dual_vars = Dualization._get_dual_variables(map, primal_con) + if length(dual_vars) != length(val) + error( + "The dual point entry for constraint $primal_con has " * + "length $(length(val)) but the dual variable " * + "length is $(length(dual_vars)).", + ) + end + for (idx, dual_var) in enumerate(dual_vars) + new_dual_point[dual_var] = val[idx] + dual_var_to_primal_con[dual_var] = primal_con + end + dual_con = Dualization._get_dual_constraint(map, primal_con) + if dual_con !== nothing + dual_con_to_primal_con[dual_con] = primal_con + # else + # if !(primal_con isa MOI.ConstraintIndex{MOI.VariableIndex,<:MOI.EqualTo} || + # primal_con isa MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.Zeros} + # SAF in EQ, etc... + #) + # error("Problem with dualization, see: $primal_con") + # end + end + end + primal_vars = MOI.get(primal_model, MOI.ListOfVariableIndices()) + dual_con_to_primal_vars = + Dict{MOI.ConstraintIndex,Vector{MOI.VariableIndex}}() + for primal_var in primal_vars + dual_con, idx = Dualization._get_dual_constraint(map, primal_var) + # TODO + idx = max(idx, 1) + if haskey(dual_con_to_primal_vars, dual_con) + vec = dual_con_to_primal_vars[dual_con] + if idx > length(vec) + resize!(vec, idx) + end + vec[idx] = primal_var + else + vec = Vector{MOI.VariableIndex}(undef, idx) + vec[idx] = primal_var + dual_con_to_primal_vars[dual_con] = vec + end + end + return new_dual_point, + dual_var_to_primal_con, + dual_con_to_primal_vars, + dual_con_to_primal_con +end + +function _analyze_dual!(model, dual_model, map, data) + dual_point, + dual_var_to_primal_con, + dual_con_to_primal_vars, + dual_con_to_primal_con = + _dual_point_to_dual_model_ref(model, map, data.dual_point) + types = MOI.get(dual_model, MOI.ListOfConstraintTypesPresent()) + for (F, S) in types + list = MOI.get(dual_model, MOI.ListOfConstraintIndices{F,S}()) + for con in list + func = MOI.get(dual_model, MOI.ConstraintFunction(), con) + failed = false + val = MOI.Utilities.eval_variables(dual_model, func) do var_idx + if !haskey(dual_point, var_idx) + if data.skip_missing + failed = true + return NaN # nothing + else + primal_con = dual_var_to_primal_con[var_idx] + error( + "Missing data for dual of constraint: $primal_con. " * + "Set skip_missing = true to ignore this error.", + ) + end + end + return dual_point[var_idx] + end + if failed + continue + end + set = MOI.get(dual_model, MOI.ConstraintSet(), con) + dist = MOI.Utilities.distance_to_set(val, set) + if dist > data.atol + if haskey(dual_con_to_primal_vars, con) + vars = dual_con_to_primal_vars[con] + if length(vars) != 1 + # TODO improve error + error( + "The dual constraint $con has " * + "length $(length(vars)) != 1", + ) + end + push!(data.dual, DualConstraintViolation(vars[], dist)) + else + con = dual_con_to_primal_con[con] + push!( + data.dual_convar, + DualConstrainedVariableViolation(con, dist), + ) + end + end + end end return end function _analyze_complementarity!(model, data) - constraint_list = - JuMP.all_constraints(model; include_variable_in_set_constraints = true) - for con in constraint_list - obj = JuMP.constraint_object(con) - func = obj.func - set = obj.set - func_val = - JuMP.value.(x -> data.primal_point[x], func) - _set_value(set) - comp_val = MOI.Utilities.set_dot(func_val, data.dual_point[con], set) - if abs(comp_val) > data.atol - push!(data.complementarity, ComplemetarityViolation(con, comp_val)) + types = MOI.get(model, MOI.ListOfConstraintTypesPresent()) + for (F, S) in types + list = MOI.get(model, MOI.ListOfConstraintIndices{F,S}()) + for con in list + func = MOI.get(model, MOI.ConstraintFunction(), con) + failed = false + val = MOI.Utilities.eval_variables(model, func) do var_idx + if !haskey(data.primal_point, var_idx) + if data.skip_missing + failed = true + return NaN # nothing + else + error( + "Missing variable in primal point: $var_idx. " * + "Set skip_missing = true to ignore this error.", + ) + end + end + return data.primal_point[var_idx] + end + set = MOI.get(model, MOI.ConstraintSet(), con) + val = val - _set_value(set) + if failed + continue + end + if !haskey(data.dual_point, con) + if data.skip_missing + continue + else + error( + "Missing dual value for constraint: $con. " * + "Set skip_missing = true to ignore this error.", + ) + end + end + if length(data.dual_point[con]) != length(val) + error( + "The dual point entry for constraint $con has " * + "length $(length(data.dual_point[con])) but the primal " * + "constraint length is $(length(val)) .", + ) + end + comp_val = MOI.Utilities.set_dot(val, data.dual_point[con], set) + if abs(comp_val) > data.atol + push!( + data.complementarity, + ComplemetarityViolation(con, comp_val), + ) + end end end return @@ -850,35 +1135,21 @@ function _set_value(set::MOI.EqualTo) return set.value end -function _analyze_objectives!( - model::JuMP.GenericModel{T}, - dual_model, - data, -) where {T} - if JuMP.primal_status(model) in - (MOI.FEASIBLE_POINT, MOI.NEARLY_FEASIBLE_POINT) - obj_val_solver = JuMP.objective_value(model) +function _analyze_objectives!(model::MOI.ModelLike, dual_model, map, data) + primal_status = MOI.get(model, MOI.PrimalStatus()) + dual_status = MOI.get(model, MOI.DualStatus()) + if primal_status in (MOI.FEASIBLE_POINT, MOI.NEARLY_FEASIBLE_POINT) + obj_val_solver = MOI.get(model, MOI.ObjectiveValue()) else obj_val_solver = nothing end - if JuMP.dual_status(model) in - (MOI.FEASIBLE_POINT, MOI.NEARLY_FEASIBLE_POINT) - dual_obj_val_solver = JuMP.dual_objective_value(model) + + if dual_status in (MOI.FEASIBLE_POINT, MOI.NEARLY_FEASIBLE_POINT) + dual_obj_val_solver = MOI.get(model, MOI.DualObjectiveValue()) else dual_obj_val_solver = nothing end - obj_func = JuMP.objective_function(model) - obj_val = JuMP.value(x -> data.primal_point[x], obj_func) - - if obj_val_solver !== nothing && - !isapprox(obj_val, obj_val_solver; atol = data.atol) - push!( - data.primal_objective_mismatch, - PrimalObjectiveMismatch(obj_val, obj_val_solver), - ) - end - if dual_obj_val_solver !== nothing && obj_val_solver !== nothing && !isapprox(obj_val_solver, dual_obj_val_solver; atol = data.atol) @@ -888,13 +1159,51 @@ function _analyze_objectives!( ) end - if dual_model !== nothing && data.dual_point !== nothing - dual_point_in_dual_model_ref = - _dual_point_to_dual_model_ref(dual_model, data.dual_point) - dual_obj_val = JuMP.value( - x -> dual_point_in_dual_model_ref[x], - JuMP.objective_function(dual_model), + obj_type = MOI.get(model, MOI.ObjectiveFunctionType()) + obj_func = MOI.get(model, MOI.ObjectiveFunction{obj_type}()) + obj_val = MOI.Utilities.eval_variables(model, obj_func) do var_idx + if !haskey(data.primal_point, var_idx) + if data.skip_missing + return NaN # nothing + else + error( + "Missing variable in primal point: $var_idx. " * + "Set skip_missing = true to ignore this error.", + ) + end + end + return data.primal_point[var_idx] + end + + if obj_val_solver !== nothing && + !isapprox(obj_val, obj_val_solver; atol = data.atol) + push!( + data.primal_objective_mismatch, + PrimalObjectiveMismatch(obj_val, obj_val_solver), ) + end + + if dual_model !== nothing && data.dual_point !== nothing + dual_point, dual_var_to_primal_con, _, _ = + _dual_point_to_dual_model_ref(model, map, data.dual_point) + + obj_type = MOI.get(dual_model, MOI.ObjectiveFunctionType()) + obj_func = MOI.get(dual_model, MOI.ObjectiveFunction{obj_type}()) + dual_obj_val = + MOI.Utilities.eval_variables(dual_model, obj_func) do var_idx + if !haskey(dual_point, var_idx) + if data.skip_missing + return NaN # nothing + else + primal_con = dual_var_to_primal_con[var_idx] + error( + "Missing data for dual of constraint: $primal_con. " * + "Set skip_missing = true to ignore this error.", + ) + end + end + return dual_point[var_idx] + end if dual_obj_val_solver !== nothing && !isapprox(dual_obj_val, dual_obj_val_solver; atol = data.atol) @@ -915,272 +1224,41 @@ function _analyze_objectives!( return end -### - -# unsafe as is its checked upstream -function _last_primal_solution(model::JuMP.GenericModel) - return Dict(v => JuMP.value(v) for v in JuMP.all_variables(model)) -end - -function _last_dual_solution(model::JuMP.GenericModel{T}) where {T} - if !JuMP.has_duals(model) - error( - "No dual solution is available. You must provide a point at " * - "which to check feasibility.", - ) - end - constraint_list = - JuMP.all_constraints(model; include_variable_in_set_constraints = true) - ret = Dict{JuMP.ConstraintRef,Vector{T}}() - for c in constraint_list - _dual = JuMP.dual(c) - if typeof(_dual) == Vector{T} - ret[c] = _dual - else - ret[c] = T[_dual] - end - end - return ret -end - -""" - dual_feasibility_report( - model::GenericModel{T}, - point::AbstractDict{GenericVariableRef{T},T} = _last_dual_solution(model), - atol::T = zero(T), - skip_missing::Bool = false, - )::Dict{Any,T} - -Given a dictionary `point`, which maps variables to dual values, return a -dictionary whose keys are the constraints with an infeasibility greater than the -supplied tolerance `atol`. The value corresponding to each key is the respective -infeasibility. Infeasibility is defined as the distance between the dual -value of the constraint (see `MOI.ConstraintDual`) and the nearest point by -Euclidean distance in the corresponding set. - -## Notes - - * If `skip_missing = true`, constraints containing variables that are not in - `point` will be ignored. - * If `skip_missing = false` and a partial dual solution is provided, an error - will be thrown. - * If no point is provided, the dual solution from the last time the model was - solved is used. - -## Example - -```jldoctest -julia> model = Model(); - -julia> @variable(model, 0.5 <= x <= 1); - -julia> dual_feasibility_report(model, Dict(x => 0.2)) -XXXX -``` -""" -function dual_feasibility_report( - model::JuMP.GenericModel{T}, - point::AbstractDict = _last_dual_solution(model); - atol::T = zero(T), - skip_missing::Bool = false, - _dual_model = nothing, # helps to avoid dualizing twice -) where {T} - if JuMP.num_nonlinear_constraints(model) > 0 - error( - "Nonlinear constraints are not supported. " * - "Use `dual_feasibility_report` instead.", - ) - end - if !skip_missing - constraint_list = JuMP.all_constraints( - model; - include_variable_in_set_constraints = true, - ) - for c in constraint_list - if !haskey(point, c) - error( - "point does not contain a dual for constraint $c. Provide " * - "a dual, or pass `skip_missing = true`.", - ) - end - end - end - dual_model = if _dual_model !== nothing - _dual_model - else - _dualize2(model) - end - dual_point = _dual_point_to_dual_model_ref(dual_model, point) - - dual_con_to_violation = JuMP.primal_feasibility_report( - dual_model, - dual_point; - atol = atol, - skip_missing = skip_missing, - ) - - # some dual model constraints are associated with primal model variables (primal_con_dual_var) - # if variable is free (almost a primal con = ConstraintIndex{MOI.VariableIndex, MOI.Reals}) - primal_var_dual_con = - dual_model.ext[:dualization_primal_dual_map].primal_var_dual_con - # if variable is bounded - primal_convar_dual_con = - dual_model.ext[:dualization_primal_dual_map].constrained_var_dual - # other dual model constraints (bounds) are associated with primal model constraints (non-bounds) - primal_con_dual_convar = - dual_model.ext[:dualization_primal_dual_map].primal_con_dual_con - - dual_con_primal_all = _build_dual_con_primal_all( - primal_var_dual_con, - primal_convar_dual_con, - primal_con_dual_convar, - ) - - ret = _fix_ret(dual_con_to_violation, model, dual_con_primal_all) - - return ret +function _last_primal_solution(model::MOI.ModelLike) + variables = MOI.get(model, MOI.ListOfVariableIndices()) + return Dict(v => MOI.get(model, MOI.VariablePrimal(), v) for v in variables) end -function _dual_point_to_dual_model_ref( - dual_model::JuMP.GenericModel{T}, - point, -) where {T} - - # point is a: - # dict mapping primal constraints to (dual) values - # we need to convert it to a: - # dict mapping the dual model variables to these (dual) values - - primal_con_dual_var = - dual_model.ext[:dualization_primal_dual_map].primal_con_dual_var - primal_con_dual_convar = - dual_model.ext[:dualization_primal_dual_map].primal_con_dual_con - - dual_point = Dict{JuMP.GenericVariableRef{T},T}() - for (jump_con, val) in point - moi_con = JuMP.index(jump_con) - if haskey(primal_con_dual_var, moi_con) - vec_vars = primal_con_dual_var[moi_con] - for (i, moi_var) in enumerate(vec_vars) - jump_var = JuMP.GenericVariableRef{T}(dual_model, moi_var) - dual_point[jump_var] = val[i] - end - elseif haskey(primal_con_dual_convar, moi_con) - moi_convar = primal_con_dual_convar[moi_con] - jump_var = JuMP.GenericVariableRef{T}( - dual_model, - MOI.VariableIndex(moi_convar.value), - ) - dual_point[jump_var] = val - else - # careful with the case where bounds do not become variables - # error("Constraint $jump_con is not associated with a variable in the dual model.") - end - end - return dual_point -end - -function _build_dual_con_primal_all( - primal_var_dual_con, - primal_convar_dual_con, - primal_con_dual_con, -) - # MOI.VariableIndex here represents MOI.ConstraintIndex{MOI.VariableIndex, MOI.Reals} - dual_con_primal_all = - Dict{MOI.ConstraintIndex,Union{MOI.ConstraintIndex,MOI.VariableIndex}}() - for (primal_var, dual_con) in primal_var_dual_con - dual_con_primal_all[dual_con] = primal_var - end - for (primal_con, dual_con) in primal_convar_dual_con - dual_con_primal_all[dual_con] = primal_con - end - for (primal_con, dual_con) in primal_con_dual_con - dual_con_primal_all[dual_con] = primal_con - end - return dual_con_primal_all -end - -function _fix_ret( - pre_ret, - primal_model::JuMP.GenericModel{T}, - dual_con_primal_all, -) where {T} - ret = Dict{ - Union{JuMP.ConstraintRef,JuMP.GenericVariableRef{T}}, - Union{T,Vector{T}}, - }() - for (jump_dual_con, val) in pre_ret - # v is a variable in the dual jump model - # we need the associated cosntraint in the primal jump model - moi_dual_con = JuMP.index(jump_dual_con) - moi_primal_something = dual_con_primal_all[moi_dual_con] - if moi_primal_something isa MOI.VariableIndex - # variable in the dual model - # constraint in the primal model - jump_primal_var = - JuMP.GenericVariableRef{T}(primal_model, moi_primal_something) - # ret[jump_primal_var] = T[val] - ret[jump_primal_var] = val - else - # constraint in the primal model - jump_primal_con = JuMP.constraint_ref_with_index( - primal_model, - moi_primal_something, - ) - # if val isa Vector - # ret[jump_primal_con] = val - # else - # ret[jump_primal_con] = T[val] - # end - ret[jump_primal_con] = val +function _last_dual_solution(model::MOI.ModelLike) + ret = Dict{MOI.ConstraintIndex,Union{Number,Vector{<:Number}}}() + types = MOI.get(model, MOI.ListOfConstraintTypesPresent()) + for (F, S) in types + list = MOI.get(model, MOI.ListOfConstraintIndices{F,S}()) + for con in list + val = MOI.get(model, MOI.ConstraintDual(), con) + ret[con] = val end end return ret end -function _dualize2( - model::JuMP.GenericModel, - optimizer_constructor = nothing; - kwargs..., -) - mode = JuMP.mode(model) - if mode == JuMP.MANUAL - error("Dualization does not support solvers in $(mode) mode") - end - dual_model = JuMP.GenericModel() - dual_problem = Dualization.DualProblem(JuMP.backend(dual_model)) - Dualization.dualize(JuMP.backend(model), dual_problem; kwargs...) - Dualization._fill_obj_dict_with_variables!(dual_model) - Dualization._fill_obj_dict_with_constraints!(dual_model) - if optimizer_constructor !== nothing - JuMP.set_optimizer(dual_model, optimizer_constructor) - end - dual_model.ext[:dualization_primal_dual_map] = dual_problem.primal_dual_map - return dual_model -end - -function _can_dualize(model::JuMP.GenericModel) - types = JuMP.list_of_constraint_types(model) +function _can_dualize(model::MOI.ModelLike) + types = MOI.get(model, MOI.ListOfConstraintTypesPresent()) - for (_F, S) in types - F = JuMP.moi_function_type(_F) + for (F, S) in types if !Dualization.supported_constraint(F, S) return false end end - _F = JuMP.objective_function_type(model) - F = JuMP.moi_function_type(_F) - - if !Dualization.supported_obj(F) - return false - end + F = MOI.get(model, MOI.ObjectiveFunctionType()) - if JuMP.num_nonlinear_constraints(model) > 0 + if !Dualization.supported_objective(F) return false end - if JuMP.objective_sense(model) == MOI.FEASIBILITY_SENSE + sense = MOI.get(model, MOI.ObjectiveSense()) + if sense == MOI.FEASIBILITY_SENSE return false end diff --git a/src/numerical.jl b/src/numerical.jl index a5379bc..66ff92b 100644 --- a/src/numerical.jl +++ b/src/numerical.jl @@ -2426,6 +2426,7 @@ end function ModelAnalyzer.summarize( io::IO, data::Data; + model = nothing, verbose = true, max_issues = ModelAnalyzer.DEFAULT_MAX_ISSUES, configurations = true, @@ -2451,6 +2452,7 @@ function ModelAnalyzer.summarize( ModelAnalyzer.summarize( io, issues, + model = model, verbose = verbose, max_issues = max_issues, ) diff --git a/test/feasibility.jl b/test/feasibility.jl index 93a46e4..a47ef41 100644 --- a/test/feasibility.jl +++ b/test/feasibility.jl @@ -5,7 +5,7 @@ module TestDualFeasibilityChecker -import ModelAnalyzer +using ModelAnalyzer using Test using JuMP import HiGHS @@ -25,34 +25,152 @@ function test_no_solution() model = Model() @variable(model, x, Bin) # do not support binary - @test_throws ErrorException ModelAnalyzer.Feasibility.dual_feasibility_report( + @test_throws "No primal" ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), model, - Dict(), ) # no dual solutions available - @test_throws ErrorException ModelAnalyzer.Feasibility.dual_feasibility_report( + # @test_throws ErrorException + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), model, + primal_point = Dict(JuMP.index(x) => 1.0), + dual_point = Dict(), ) - # non linear not accepted + list = ModelAnalyzer.list_of_issue_types(data) + @test isempty(list) + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(JuMP.index(x) => 1.0), + # dual_point = Dict(), + ) + list = ModelAnalyzer.list_of_issue_types(data) + @test isempty(list) + # non linear in primal accepted @constraint(model, c, x^4 >= 0) # this will make the model non-linear - @test_throws ErrorException ModelAnalyzer.Feasibility.dual_feasibility_report( + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), model, - Dict(), # to skip dual solutions error + primal_point = Dict(JuMP.index(x) => 1.0), ) + list = ModelAnalyzer.list_of_issue_types(data) + @test isempty(list) end function test_only_bounds() - # in this case the dual has no varaibles and has innocuous constraints - # this needs to be reviewed in Dualization.jl model = Model() @variable(model, x >= 0) @objective(model, Min, x) - report = ModelAnalyzer.Feasibility.dual_feasibility_report( + + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(JuMP.index(x) => 0.0), + dual_point = Dict(JuMP.index(LowerBoundRef(x)) => 1.0), + ) + list = ModelAnalyzer.list_of_issue_types(data) + @test isempty(list) + + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(JuMP.index(x) => -1.0), + dual_point = Dict(JuMP.index(LowerBoundRef(x)) => 1.0), + ) + list = ModelAnalyzer.list_of_issue_types(data) + @test list == Type[ + ModelAnalyzer.Feasibility.PrimalViolation, + ModelAnalyzer.Feasibility.ComplemetarityViolation, + ModelAnalyzer.Feasibility.PrimalDualMismatch, + ] + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.PrimalViolation, + ) + @test ret[] == ModelAnalyzer.Feasibility.PrimalViolation( + JuMP.index(LowerBoundRef(x)), + 1.0, + ) + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.ComplemetarityViolation, + ) + @test ret[] == ModelAnalyzer.Feasibility.ComplemetarityViolation( + JuMP.index(LowerBoundRef(x)), + -1.0, + ) + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.PrimalDualMismatch, + ) + @test ret[] == ModelAnalyzer.Feasibility.PrimalDualMismatch(-1.0, 0.0) + + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), model, - Dict(LowerBoundRef(x) => 1.0), + primal_point = Dict(JuMP.index(x) => 1.0), + dual_point = Dict(JuMP.index(LowerBoundRef(x)) => 1.0), + ) + list = ModelAnalyzer.list_of_issue_types(data) + @test list == [ + ModelAnalyzer.Feasibility.ComplemetarityViolation, + ModelAnalyzer.Feasibility.PrimalDualMismatch, + ] + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.ComplemetarityViolation, + ) + @test ret[] == ModelAnalyzer.Feasibility.ComplemetarityViolation( + JuMP.index(LowerBoundRef(x)), + 1.0, + ) + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.PrimalDualMismatch, ) - @test isempty(report) + @test ret[] == ModelAnalyzer.Feasibility.PrimalDualMismatch(1.0, 0.0) + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(JuMP.index(x) => 1.0), + dual_point = Dict(JuMP.index(LowerBoundRef(x)) => -1.0), + ) + list = ModelAnalyzer.list_of_issue_types(data) + @test list == [ + ModelAnalyzer.Feasibility.DualConstraintViolation, + ModelAnalyzer.Feasibility.DualConstrainedVariableViolation, + ModelAnalyzer.Feasibility.ComplemetarityViolation, + ModelAnalyzer.Feasibility.PrimalDualMismatch, + ] + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.DualConstraintViolation, + ) + @test ret[] == + ModelAnalyzer.Feasibility.DualConstraintViolation(JuMP.index(x), 2.0) + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.DualConstrainedVariableViolation, + ) + @test ret[] == ModelAnalyzer.Feasibility.DualConstrainedVariableViolation( + JuMP.index(LowerBoundRef(x)), + 1.0, + ) + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.ComplemetarityViolation, + ) + @test ret[] == ModelAnalyzer.Feasibility.ComplemetarityViolation( + JuMP.index(LowerBoundRef(x)), + -1.0, + ) + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.PrimalDualMismatch, + ) + @test ret[] == ModelAnalyzer.Feasibility.PrimalDualMismatch(1.0, 0.0) return end @@ -64,24 +182,110 @@ function test_no_lb() # the dual is: # Max 0 # Subject to - # y == 1 (as a constraint) # from x, a free "bounded" varaible - # y >= 0 (as a bound) # from c, a ">=" constraint + # y == 1 (as a constraint) + # y >= 0 (as a bound) # mayber force fail here # @test_throws ErrorException - report = - ModelAnalyzer.Feasibility.dual_feasibility_report(model, Dict(c => 1.0)) - @test isempty(report) - report = ModelAnalyzer.Feasibility.dual_feasibility_report( + + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(JuMP.index(x) => 0.0), + dual_point = Dict(JuMP.index(c) => 1.0), + ) + list = ModelAnalyzer.list_of_issue_types(data) + @test isempty(list) + + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), model, - Dict(c => [1.0]), + primal_point = Dict(JuMP.index(x) => -1.0), + dual_point = Dict(JuMP.index(c) => 1.0), + ) + list = ModelAnalyzer.list_of_issue_types(data) + @test list == Type[ + ModelAnalyzer.Feasibility.PrimalViolation, + ModelAnalyzer.Feasibility.ComplemetarityViolation, + ModelAnalyzer.Feasibility.PrimalDualMismatch, + ] + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.PrimalViolation, ) - @test isempty(report) - report = ModelAnalyzer.Feasibility.dual_feasibility_report( + @test ret[] == ModelAnalyzer.Feasibility.PrimalViolation(JuMP.index(c), 1.0) + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.ComplemetarityViolation, + ) + @test ret[] == + ModelAnalyzer.Feasibility.ComplemetarityViolation(JuMP.index(c), -1.0) + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.PrimalDualMismatch, + ) + @test ret[] == ModelAnalyzer.Feasibility.PrimalDualMismatch(-1.0, 0.0) + + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), model, - Dict(c => [3.3]), + primal_point = Dict(JuMP.index(x) => 1.0), + dual_point = Dict(JuMP.index(c) => 1.0), + ) + list = ModelAnalyzer.list_of_issue_types(data) + @test list == [ + ModelAnalyzer.Feasibility.ComplemetarityViolation, + ModelAnalyzer.Feasibility.PrimalDualMismatch, + ] + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.ComplemetarityViolation, + ) + @test ret[] == + ModelAnalyzer.Feasibility.ComplemetarityViolation(JuMP.index(c), 1.0) + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.PrimalDualMismatch, ) - @test report[x] == 2.3 - @test length(report) == 1 + @test ret[] == ModelAnalyzer.Feasibility.PrimalDualMismatch(1.0, 0.0) + + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(JuMP.index(x) => 1.0), + dual_point = Dict(JuMP.index(c) => -1.0), + ) + list = ModelAnalyzer.list_of_issue_types(data) + @test list == [ + ModelAnalyzer.Feasibility.DualConstraintViolation, + ModelAnalyzer.Feasibility.DualConstrainedVariableViolation, + ModelAnalyzer.Feasibility.ComplemetarityViolation, + ModelAnalyzer.Feasibility.PrimalDualMismatch, + ] + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.DualConstraintViolation, + ) + @test ret[] == + ModelAnalyzer.Feasibility.DualConstraintViolation(JuMP.index(x), 2.0) + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.DualConstrainedVariableViolation, + ) + @test ret[] == ModelAnalyzer.Feasibility.DualConstrainedVariableViolation( + JuMP.index(c), + 1.0, + ) + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.ComplemetarityViolation, + ) + @test ret[] == + ModelAnalyzer.Feasibility.ComplemetarityViolation(JuMP.index(c), -1.0) + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.PrimalDualMismatch, + ) + @test ret[] == ModelAnalyzer.Feasibility.PrimalDualMismatch(1.0, 0.0) end function test_lb0() @@ -92,25 +296,69 @@ function test_lb0() # the dual is: # Max 0.5 * y # Subject to - # - y >= -1 (as a constraint) # from x >= 0 (bound) - # y >= 0 (as a bound) # from c, a ">=" constraint - report = ModelAnalyzer.Feasibility.dual_feasibility_report( + # ... + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), model, - Dict(c => [1.0], LowerBoundRef(x) => [0.0]), + primal_point = Dict(JuMP.index(x) => 0.5), + dual_point = Dict( + JuMP.index(c) => 1.0, + JuMP.index(LowerBoundRef(x)) => 0.0, + ), ) - @test isempty(report) - report = ModelAnalyzer.Feasibility.dual_feasibility_report( + list = ModelAnalyzer.list_of_issue_types(data) + @test isempty(list) + + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), model, - Dict(c => [3.3], LowerBoundRef(x) => [0.0]), + primal_point = Dict(JuMP.index(x) => 0.5), + dual_point = Dict( + JuMP.index(c) => 3.3, + JuMP.index(LowerBoundRef(x)) => 0.0, + ), + ) + list = ModelAnalyzer.list_of_issue_types(data) + @test list == [ + ModelAnalyzer.Feasibility.DualConstraintViolation, + ModelAnalyzer.Feasibility.PrimalDualMismatch, + ] + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.DualConstraintViolation, ) - @test report[LowerBoundRef(x)] == 2.3 - @test length(report) == 1 - report = ModelAnalyzer.Feasibility.dual_feasibility_report( + @test ret[] == + ModelAnalyzer.Feasibility.DualConstraintViolation(JuMP.index(x), 2.3) + + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), model, - Dict(c => [-3.3], LowerBoundRef(x) => [0.0]), + primal_point = Dict(JuMP.index(x) => 0.5), + dual_point = Dict( + JuMP.index(c) => -3.3, + JuMP.index(LowerBoundRef(x)) => 0.0, + ), + ) + list = ModelAnalyzer.list_of_issue_types(data) + @test list == [ + ModelAnalyzer.Feasibility.DualConstraintViolation, + ModelAnalyzer.Feasibility.DualConstrainedVariableViolation, + ModelAnalyzer.Feasibility.PrimalDualMismatch, + ] + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.DualConstraintViolation, + ) + @test ret[] == + ModelAnalyzer.Feasibility.DualConstraintViolation(JuMP.index(x), 4.3) + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.DualConstrainedVariableViolation, + ) + @test ret[] == ModelAnalyzer.Feasibility.DualConstrainedVariableViolation( + JuMP.index(c), + 3.3, ) - @test report[c] == 3.3 - @test length(report) == 1 end function test_lb2() @@ -121,42 +369,108 @@ function test_lb2() # the dual is: # Max 0.5 * y + 2 * z # Subject to - # y + z == 1 (as a constraint) # from x, a free variable (bound is considered below) - # z >= 0 (as a bound) # from the "constraint" x >= 2 (bound in the above example) - # y >= 0 (as a bound) # from c, a ">=" constraint - report = ModelAnalyzer.Feasibility.dual_feasibility_report( + # ... + + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), model, - Dict(c => [1.0], LowerBoundRef(x) => [0.0]), + primal_point = Dict(JuMP.index(x) => 2.0), + dual_point = Dict( + JuMP.index(c) => 0.0, + JuMP.index(LowerBoundRef(x)) => 1.0, + ), ) - @test isempty(report) - report = ModelAnalyzer.Feasibility.dual_feasibility_report( + list = ModelAnalyzer.list_of_issue_types(data) + @test isempty(list) + + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), model, - Dict(c => [3.3], LowerBoundRef(x) => [0.0]), + primal_point = Dict(JuMP.index(x) => 2.0), + dual_point = Dict( + JuMP.index(c) => 3.3, + JuMP.index(LowerBoundRef(x)) => 0.0, + ), + ) + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.DualConstraintViolation, ) - @test report[x] == 2.3 - @test length(report) == 1 - report = ModelAnalyzer.Feasibility.dual_feasibility_report( + @test ret[] == + ModelAnalyzer.Feasibility.DualConstraintViolation(JuMP.index(x), 2.3) + + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), model, - Dict(c => [-3.3], LowerBoundRef(x) => [0.0]), + primal_point = Dict(JuMP.index(x) => 2.0), + dual_point = Dict( + JuMP.index(c) => -3.3, + JuMP.index(LowerBoundRef(x)) => 0.0, + ), + ) + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.DualConstraintViolation, ) - @test report[x] == 4.3 - @test report[c] == 3.3 - @test length(report) == 2 - report = ModelAnalyzer.Feasibility.dual_feasibility_report( + @test ret[] == + ModelAnalyzer.Feasibility.DualConstraintViolation(JuMP.index(x), 4.3) + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.DualConstrainedVariableViolation, + ) + @test ret[] == ModelAnalyzer.Feasibility.DualConstrainedVariableViolation( + JuMP.index(c), + 3.3, + ) + + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), model, - Dict(c => [-3.3], LowerBoundRef(x) => [-1.0]), + primal_point = Dict(JuMP.index(x) => 2.0), + dual_point = Dict( + JuMP.index(c) => -3.3, + JuMP.index(LowerBoundRef(x)) => -1.0, + ), ) - @test report[x] == 5.3 - @test report[c] == 3.3 - @test report[LowerBoundRef(x)] == 1.0 - @test length(report) == 3 - report = ModelAnalyzer.Feasibility.dual_feasibility_report( + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.DualConstraintViolation, + ) + @test ret[] == + ModelAnalyzer.Feasibility.DualConstraintViolation(JuMP.index(x), 5.3) + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.DualConstrainedVariableViolation, + ) + @test ret[1] == ModelAnalyzer.Feasibility.DualConstrainedVariableViolation( + JuMP.index(c), + 3.3, + ) + @test ret[2] == ModelAnalyzer.Feasibility.DualConstrainedVariableViolation( + JuMP.index(LowerBoundRef(x)), + 1.0, + ) + + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), model, - Dict(c => [-3.3]), + primal_point = Dict(JuMP.index(x) => 2.0), + dual_point = Dict( + JuMP.index(c) => -3.3, + # JuMP.index(LowerBoundRef(x)) => -1.0, + ), skip_missing = true, ) - @test report[c] == 3.3 - @test length(report) == 1 + + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.DualConstrainedVariableViolation, + ) + @test ret[] == ModelAnalyzer.Feasibility.DualConstrainedVariableViolation( + JuMP.index(c), + 3.3, + ) + return end function test_analyse_simple() @@ -202,9 +516,7 @@ function test_with_interval() @constraint(model, c, 2 * x in MOI.Interval(0.0, 3.0)) @objective(model, Min, x) optimize!(model) - @test !ModelAnalyzer.Feasibility._can_dualize(model) - # TODO this should eb a waning at least - data = ModelAnalyzer.analyze(ModelAnalyzer.Feasibility.Analyzer(), model) + @test !ModelAnalyzer.Feasibility._can_dualize(JuMP.backend(model)) return end @@ -261,94 +573,97 @@ function test_analyse_mip() return end -function test_analyse_no_opt() - model = Model(HiGHS.Optimizer) - set_silent(model) - @variable(model, x) - @constraint(model, c, x >= 0) - @objective(model, Min, x) - - # test no primal point - @test_throws ErrorException ModelAnalyzer.analyze( - ModelAnalyzer.Feasibility.Analyzer(), - model, - ) +#= - # test no dual point - @test_throws ErrorException ModelAnalyzer.analyze( - ModelAnalyzer.Feasibility.Analyzer(), - model, - primal_point = Dict(x => 1.0), - dual_check = true, - ) - - data = ModelAnalyzer.analyze( - ModelAnalyzer.Feasibility.Analyzer(), - model, - primal_point = Dict(x => 1.0), - dual_check = false, - ) - list = ModelAnalyzer.list_of_issue_types(data) - @test length(list) == 0 - - data = ModelAnalyzer.analyze( - ModelAnalyzer.Feasibility.Analyzer(), - model, - primal_point = Dict(x => -1.0), - dual_check = false, - ) - list = ModelAnalyzer.list_of_issue_types(data) - @test length(list) == 1 - ret = ModelAnalyzer.list_of_issues(data, list[1]) - @test ret[] == ModelAnalyzer.Feasibility.PrimalViolation(c, 1.0) - - data = ModelAnalyzer.analyze( - ModelAnalyzer.Feasibility.Analyzer(), - model, - primal_point = Dict(x => 1.0), - dual_point = Dict(c => 1.0), - ) - list = ModelAnalyzer.list_of_issue_types(data) - @test length(list) == 2 - ret = ModelAnalyzer.list_of_issues(data, list[1]) - @test ret[1] == ModelAnalyzer.Feasibility.ComplemetarityViolation(c, 1.0) - ret = ModelAnalyzer.list_of_issues(data, list[2]) - @test ret[1] == ModelAnalyzer.Feasibility.PrimalDualMismatch(1.0, 0.0) - - data = ModelAnalyzer.analyze( - ModelAnalyzer.Feasibility.Analyzer(), - model, - primal_point = Dict(x => 0.0), - dual_point = Dict(c => 1.0), - ) - list = ModelAnalyzer.list_of_issue_types(data) - @test length(list) == 0 - - data = ModelAnalyzer.analyze( - ModelAnalyzer.Feasibility.Analyzer(), - model, - primal_point = Dict(x => -1.0), - dual_point = Dict(c => 2.0), - ) - list = ModelAnalyzer.list_of_issue_types(data) - @test length(list) == 4 - ret = ModelAnalyzer.list_of_issues(data, list[1]) - @test ret[1] == ModelAnalyzer.Feasibility.PrimalViolation(c, 1.0) - ret = ModelAnalyzer.list_of_issues(data, list[2]) - @test ret[1] == ModelAnalyzer.Feasibility.DualViolation(x, 1.0) - ret = ModelAnalyzer.list_of_issues(data, list[3]) - @test ret[1] == ModelAnalyzer.Feasibility.ComplemetarityViolation(c, -2.0) - ret = ModelAnalyzer.list_of_issues(data, list[4]) - @test ret[1] == ModelAnalyzer.Feasibility.PrimalDualMismatch(-1.0, 0.0) - - buf = IOBuffer() - - ModelAnalyzer.summarize(buf, data) - - ModelAnalyzer.summarize(buf, data, verbose = false) - - return +function test_analyse_no_opt() +model = Model(HiGHS.Optimizer) +set_silent(model) +@variable(model, x) +@constraint(model, c, x >= 0) +@objective(model, Min, x) + +# test no primal point +@test_throws ErrorException ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, +) + +# test no dual point +@test_throws ErrorException ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(x => 1.0), + dual_check = true, +) + +data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(x => 1.0), + dual_check = false, +) +list = ModelAnalyzer.list_of_issue_types(data) +@test length(list) == 0 + +data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(x => -1.0), + dual_check = false, +) +list = ModelAnalyzer.list_of_issue_types(data) +@test length(list) == 1 +ret = ModelAnalyzer.list_of_issues(data, list[1]) +@test ret[] == ModelAnalyzer.Feasibility.PrimalViolation(c, 1.0) + +data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(x => 1.0), + dual_point = Dict(c => 1.0), +) +list = ModelAnalyzer.list_of_issue_types(data) +@test length(list) == 2 +ret = ModelAnalyzer.list_of_issues(data, list[1]) +@test ret[1] == ModelAnalyzer.Feasibility.ComplemetarityViolation(c, 1.0) +ret = ModelAnalyzer.list_of_issues(data, list[2]) +@test ret[1] == ModelAnalyzer.Feasibility.PrimalDualMismatch(1.0, 0.0) + +data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(x => 0.0), + dual_point = Dict(c => 1.0), +) +list = ModelAnalyzer.list_of_issue_types(data) +@test length(list) == 0 + +data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(x => -1.0), + dual_point = Dict(c => 2.0), +) +list = ModelAnalyzer.list_of_issue_types(data) +@test length(list) == 4 +ret = ModelAnalyzer.list_of_issues(data, list[1]) +@test ret[1] == ModelAnalyzer.Feasibility.PrimalViolation(c, 1.0) +ret = ModelAnalyzer.list_of_issues(data, list[2]) +@test ret[1] == ModelAnalyzer.Feasibility.DualConstraintViolation(x, 1.0) +ret = ModelAnalyzer.list_of_issues(data, list[3]) +@test ret[1] == ModelAnalyzer.Feasibility.ComplemetarityViolation(c, -2.0) +ret = ModelAnalyzer.list_of_issues(data, list[4]) +@test ret[1] == ModelAnalyzer.Feasibility.PrimalDualMismatch(-1.0, 0.0) + +buf = IOBuffer() + +ModelAnalyzer.summarize(buf, data) + +ModelAnalyzer.summarize(buf, data, verbose = false) + +return end +=# # these tests are harder to permorm with a real solver as they tipically # return coherent objectives diff --git a/test/runtests.jl b/test/runtests.jl index fb558e5..4331481 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,8 +6,9 @@ using Test @testset "ModelAnalyzer" begin - for file in ["numerical.jl"] - # for file in readdir(@__DIR__) + for file in ["feasibility.jl", "numerical.jl"] + # for file in ["feasibility.jl"] + # for file in readdir(@__DIR__) if !endswith(file, ".jl") || file in ("runtests.jl",) continue end From 87aacb238dcce495f1491d828f5d14ebc879b5f4 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Wed, 7 May 2025 10:39:18 -0300 Subject: [PATCH 09/34] Finish feasibility rewrite --- test/feasibility.jl | 179 ++++++++++++++++++++++---------------------- test/runtests.jl | 4 +- 2 files changed, 92 insertions(+), 91 deletions(-) diff --git a/test/feasibility.jl b/test/feasibility.jl index a47ef41..cbeda18 100644 --- a/test/feasibility.jl +++ b/test/feasibility.jl @@ -573,97 +573,98 @@ function test_analyse_mip() return end -#= - function test_analyse_no_opt() -model = Model(HiGHS.Optimizer) -set_silent(model) -@variable(model, x) -@constraint(model, c, x >= 0) -@objective(model, Min, x) - -# test no primal point -@test_throws ErrorException ModelAnalyzer.analyze( - ModelAnalyzer.Feasibility.Analyzer(), - model, -) - -# test no dual point -@test_throws ErrorException ModelAnalyzer.analyze( - ModelAnalyzer.Feasibility.Analyzer(), - model, - primal_point = Dict(x => 1.0), - dual_check = true, -) - -data = ModelAnalyzer.analyze( - ModelAnalyzer.Feasibility.Analyzer(), - model, - primal_point = Dict(x => 1.0), - dual_check = false, -) -list = ModelAnalyzer.list_of_issue_types(data) -@test length(list) == 0 - -data = ModelAnalyzer.analyze( - ModelAnalyzer.Feasibility.Analyzer(), - model, - primal_point = Dict(x => -1.0), - dual_check = false, -) -list = ModelAnalyzer.list_of_issue_types(data) -@test length(list) == 1 -ret = ModelAnalyzer.list_of_issues(data, list[1]) -@test ret[] == ModelAnalyzer.Feasibility.PrimalViolation(c, 1.0) - -data = ModelAnalyzer.analyze( - ModelAnalyzer.Feasibility.Analyzer(), - model, - primal_point = Dict(x => 1.0), - dual_point = Dict(c => 1.0), -) -list = ModelAnalyzer.list_of_issue_types(data) -@test length(list) == 2 -ret = ModelAnalyzer.list_of_issues(data, list[1]) -@test ret[1] == ModelAnalyzer.Feasibility.ComplemetarityViolation(c, 1.0) -ret = ModelAnalyzer.list_of_issues(data, list[2]) -@test ret[1] == ModelAnalyzer.Feasibility.PrimalDualMismatch(1.0, 0.0) - -data = ModelAnalyzer.analyze( - ModelAnalyzer.Feasibility.Analyzer(), - model, - primal_point = Dict(x => 0.0), - dual_point = Dict(c => 1.0), -) -list = ModelAnalyzer.list_of_issue_types(data) -@test length(list) == 0 - -data = ModelAnalyzer.analyze( - ModelAnalyzer.Feasibility.Analyzer(), - model, - primal_point = Dict(x => -1.0), - dual_point = Dict(c => 2.0), -) -list = ModelAnalyzer.list_of_issue_types(data) -@test length(list) == 4 -ret = ModelAnalyzer.list_of_issues(data, list[1]) -@test ret[1] == ModelAnalyzer.Feasibility.PrimalViolation(c, 1.0) -ret = ModelAnalyzer.list_of_issues(data, list[2]) -@test ret[1] == ModelAnalyzer.Feasibility.DualConstraintViolation(x, 1.0) -ret = ModelAnalyzer.list_of_issues(data, list[3]) -@test ret[1] == ModelAnalyzer.Feasibility.ComplemetarityViolation(c, -2.0) -ret = ModelAnalyzer.list_of_issues(data, list[4]) -@test ret[1] == ModelAnalyzer.Feasibility.PrimalDualMismatch(-1.0, 0.0) - -buf = IOBuffer() - -ModelAnalyzer.summarize(buf, data) - -ModelAnalyzer.summarize(buf, data, verbose = false) - -return + model = Model(HiGHS.Optimizer) + set_silent(model) + @variable(model, x) + @constraint(model, c, x >= 0) + @objective(model, Min, x) + + # test no primal point + @test_throws ErrorException ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + ) + + # test no dual point + @test_throws ErrorException ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(JuMP.index(x) => 1.0), + dual_check = true, + ) + + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(JuMP.index(x) => 1.0), + dual_check = false, + ) + list = ModelAnalyzer.list_of_issue_types(data) + @test length(list) == 0 + + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(JuMP.index(x) => -1.0), + dual_check = false, + ) + list = ModelAnalyzer.list_of_issue_types(data) + @test length(list) == 1 + ret = ModelAnalyzer.list_of_issues(data, list[1]) + @test ret[] == ModelAnalyzer.Feasibility.PrimalViolation(JuMP.index(c), 1.0) + + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(JuMP.index(x) => 1.0), + dual_point = Dict(JuMP.index(c) => 1.0), + ) + list = ModelAnalyzer.list_of_issue_types(data) + @test length(list) == 2 + ret = ModelAnalyzer.list_of_issues(data, list[1]) + @test ret[1] == + ModelAnalyzer.Feasibility.ComplemetarityViolation(JuMP.index(c), 1.0) + ret = ModelAnalyzer.list_of_issues(data, list[2]) + @test ret[1] == ModelAnalyzer.Feasibility.PrimalDualMismatch(1.0, 0.0) + + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(JuMP.index(x) => 0.0), + dual_point = Dict(JuMP.index(c) => 1.0), + ) + list = ModelAnalyzer.list_of_issue_types(data) + @test length(list) == 0 + + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(JuMP.index(x) => -1.0), + dual_point = Dict(JuMP.index(c) => 2.0), + ) + list = ModelAnalyzer.list_of_issue_types(data) + @test length(list) == 4 + ret = ModelAnalyzer.list_of_issues(data, list[1]) + @test ret[1] == + ModelAnalyzer.Feasibility.PrimalViolation(JuMP.index(c), 1.0) + ret = ModelAnalyzer.list_of_issues(data, list[2]) + @test ret[1] == + ModelAnalyzer.Feasibility.DualConstraintViolation(JuMP.index(x), 1.0) + ret = ModelAnalyzer.list_of_issues(data, list[3]) + @test ret[1] == + ModelAnalyzer.Feasibility.ComplemetarityViolation(JuMP.index(c), -2.0) + ret = ModelAnalyzer.list_of_issues(data, list[4]) + @test ret[1] == ModelAnalyzer.Feasibility.PrimalDualMismatch(-1.0, 0.0) + + buf = IOBuffer() + + ModelAnalyzer.summarize(buf, data) + + ModelAnalyzer.summarize(buf, data, verbose = false) + + return end -=# # these tests are harder to permorm with a real solver as they tipically # return coherent objectives diff --git a/test/runtests.jl b/test/runtests.jl index 4331481..dadf665 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,8 +7,8 @@ using Test @testset "ModelAnalyzer" begin for file in ["feasibility.jl", "numerical.jl"] - # for file in ["feasibility.jl"] - # for file in readdir(@__DIR__) + # for file in ["infeasibility.jl"] + # for file in readdir(@__DIR__) if !endswith(file, ".jl") || file in ("runtests.jl",) continue end From 85e00fb37520daafdbd1eaa3c410487100ff8c0b Mon Sep 17 00:00:00 2001 From: joaquimg Date: Thu, 8 May 2025 01:36:46 -0300 Subject: [PATCH 10/34] MOI infeasibility analysis --- src/ModelAnalyzer.jl | 2 +- src/_eval_variables.jl | 31 ++++ src/infeasibility.jl | 362 +++++++++++++++++++++++++++-------------- src/intervals.jl | 1 + test/infeasibility.jl | 37 +++-- test/runtests.jl | 4 +- 6 files changed, 298 insertions(+), 139 deletions(-) create mode 100644 src/_eval_variables.jl diff --git a/src/ModelAnalyzer.jl b/src/ModelAnalyzer.jl index f08659d..22b2637 100644 --- a/src/ModelAnalyzer.jl +++ b/src/ModelAnalyzer.jl @@ -216,6 +216,6 @@ end include("numerical.jl") include("feasibility.jl") -# include("infeasibility.jl") +include("infeasibility.jl") end # module ModelAnalyzer diff --git a/src/_eval_variables.jl b/src/_eval_variables.jl new file mode 100644 index 0000000..7c7d1d7 --- /dev/null +++ b/src/_eval_variables.jl @@ -0,0 +1,31 @@ +# 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, t::MOI.ScalarQuadraticTerm) + out = t.coefficient * value_fn(t.variable_1) * value_fn(t.variable_2) + return t.variable_1 == t.variable_2 ? out / 2 : out +end + +_eval_variables(value_fn::Function, f::MOI.VariableIndex) = value_fn(f) + +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 = Base.promote_op(*, 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.jl b/src/infeasibility.jl index aa3af17..9232cc6 100644 --- a/src/infeasibility.jl +++ b/src/infeasibility.jl @@ -6,15 +6,15 @@ module Infeasibility import ModelAnalyzer -import JuMP -import JuMP.MOI as MOI +import MathOptInterface as MOI include("intervals.jl") +include("_eval_variables.jl") """ Analyzer() <: ModelAnalyzer.AbstractAnalyzer -The `Analyzer` type is used to perform infeasibility analysis on a JuMP model. +The `Analyzer` type is used to perform infeasibility analysis on a model. ## Example ```julia @@ -33,7 +33,7 @@ struct Analyzer <: ModelAnalyzer.AbstractAnalyzer end """ AbstractInfeasibilitylIssue -Abstract type for infeasibility issues found during the analysis of a JuMP +Abstract type for infeasibility issues found during the analysis of a model. """ abstract type AbstractInfeasibilitylIssue <: ModelAnalyzer.AbstractIssue end @@ -50,7 +50,7 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Infeasibility.InfeasibleBounds) ```` """ struct InfeasibleBounds{T} <: AbstractInfeasibilitylIssue - variable::JuMP.GenericVariableRef{T} + variable::MOI.VariableIndex lb::T ub::T end @@ -71,7 +71,7 @@ julia> ModelAnalyzer.summarize( ``` """ struct InfeasibleIntegrality{T} <: AbstractInfeasibilitylIssue - variable::JuMP.GenericVariableRef{T} + variable::MOI.VariableIndex lb::T ub::T set::Union{MOI.Integer,MOI.ZeroOne}#, MOI.Semicontinuous{T}, MOI.Semiinteger{T}} @@ -95,7 +95,7 @@ julia> ModelAnalyzer.summarize( ``` """ struct InfeasibleConstraintRange{T} <: AbstractInfeasibilitylIssue - constraint::JuMP.ConstraintRef + constraint::MOI.ConstraintIndex lb::T ub::T set::Union{MOI.EqualTo{T},MOI.LessThan{T},MOI.GreaterThan{T}} @@ -118,7 +118,7 @@ julia> ModelAnalyzer.summarize( ``` """ struct IrreducibleInfeasibleSubset <: AbstractInfeasibilitylIssue - constraint::Vector{JuMP.ConstraintRef} + constraint::Vector{<:MOI.ConstraintIndex} end """ @@ -142,55 +142,119 @@ end function ModelAnalyzer.analyze( ::Analyzer, - model::JuMP.GenericModel{T}; + model::MOI.ModelLike; optimizer = nothing, -) where {T} +) out = Data() - variables = Dict{JuMP.GenericVariableRef{T},Interval{T}}() + 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 + + 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 + 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 - # first layer of infeasibility analysis is bounds consistency bounds_consistent = true - for var in JuMP.all_variables(model) - lb = if JuMP.has_lower_bound(var) - JuMP.lower_bound(var) - elseif JuMP.is_fixed(var) - JuMP.fix_value(var) - else - -Inf - end - ub = if JuMP.has_upper_bound(var) - JuMP.upper_bound(var) - elseif JuMP.is_fixed(var) - JuMP.fix_value(var) - else - Inf - end - if JuMP.is_integer(var) - if abs(ub - lb) < 1 && ceil(ub) == ceil(lb) - push!( - out.infeasible_integrality, - InfeasibleIntegrality(var, lb, ub, MOI.Integer()), - ) - bounds_consistent = false - end + + 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) + push!( + out.infeasible_integrality, + InfeasibleIntegrality(func, _lb, _ub, MOI.Integer()), + ) + bounds_consistent = false end - if JuMP.is_binary(var) - if lb > 0 && ub < 1 - push!( - out.infeasible_integrality, - InfeasibleIntegrality(var, lb, ub, MOI.ZeroOne()), - ) - 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 + push!( + out.infeasible_integrality, + InfeasibleIntegrality(func, _lb, _ub, MOI.ZeroOne()), + ) + bounds_consistent = false end - if lb > ub - push!(out.infeasible_bounds, InfeasibleBounds(var, lb, ub)) + 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) + variables[var] = Interval(_lb, _ub) end end + # check PSD diagonal >= 0 ? # other cones? if !bounds_consistent @@ -199,57 +263,94 @@ function ModelAnalyzer.analyze( # second layer of infeasibility analysis is constraint range analysis range_consistent = true - for (F, S) in JuMP.list_of_constraint_types(model) - F != JuMP.GenericAffExpr{T,JuMP.GenericVariableRef{T}} && continue - # TODO: handle quadratics - !(S in (MOI.EqualTo{T}, MOI.LessThan{T}, MOI.GreaterThan{T})) && + + 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 + if !haskey(variables, var_idx) + failed = true + return Interval(-Inf, Inf) + end + return variables[var_idx] + end + if failed continue - for con in JuMP.all_constraints(model, F, S) - con_obj = JuMP.constraint_object(con) - interval = JuMP.value(x -> variables[x], con_obj.func) - if con_obj.set isa MOI.EqualTo{T} - rhs = con_obj.set.value - if interval.lo > rhs || interval.hi < rhs - push!( - out.constraint_range, - InfeasibleConstraintRange( - con, - interval.lo, - interval.hi, - con_obj.set, - ), - ) - range_consistent = false - end - elseif con_obj.set isa MOI.LessThan{T} - rhs = con_obj.set.upper - if interval.lo > rhs - push!( - out.constraint_range, - InfeasibleConstraintRange( - con, - interval.lo, - interval.hi, - con_obj.set, - ), - ) - range_consistent = false - end - elseif con_obj.set isa MOI.GreaterThan{T} - rhs = con_obj.set.lower - if interval.hi < rhs - push!( - out.constraint_range, - InfeasibleConstraintRange( - con, - interval.lo, - interval.hi, - con_obj.set, - ), - ) - range_consistent = false - end + 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 + 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 + 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 @@ -263,15 +364,17 @@ function ModelAnalyzer.analyze( 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 + # 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 iis_elastic_filter(original_model::JuMP.GenericModel, optimizer) # if JuMP.termination_status(original_model) == MOI.OPTIMIZE_NOT_CALLED @@ -390,6 +493,8 @@ function iis_elastic_filter(original_model::JuMP.GenericModel, optimizer) return iis end +=# + # API function ModelAnalyzer._summarize(io::IO, ::Type{InfeasibleBounds{T}}) where {T} @@ -536,17 +641,29 @@ function ModelAnalyzer._verbose_summarize( ) end -function ModelAnalyzer._summarize(io::IO, issue::InfeasibleBounds{T}) where {T} - return print(io, _name(issue.variable), " : ", issue.lb, " !<= ", issue.ub) +function ModelAnalyzer._summarize( + io::IO, + issue::InfeasibleBounds{T}, + model, +) where {T} + return print( + io, + ModelAnalyzer._name(issue.variable, model), + " : ", + issue.lb, + " !<= ", + issue.ub, + ) end function ModelAnalyzer._summarize( io::IO, issue::InfeasibleIntegrality{T}, + model, ) where {T} return print( io, - _name(issue.variable), + ModelAnalyzer._name(issue.variable, model), " : [", issue.lb, "; ", @@ -559,10 +676,11 @@ end function ModelAnalyzer._summarize( io::IO, issue::InfeasibleConstraintRange{T}, + model, ) where {T} return print( io, - _name(issue.constraint), + ModelAnalyzer._name(issue.constraint, model), " : [", issue.lb, "; ", @@ -572,18 +690,27 @@ function ModelAnalyzer._summarize( ) end -function ModelAnalyzer._summarize(io::IO, issue::IrreducibleInfeasibleSubset) - return print(io, "IIS: ", join(map(_name, issue.constraint), ", ")) +function ModelAnalyzer._summarize( + io::IO, + issue::IrreducibleInfeasibleSubset, + model, +) + return print( + io, + "IIS: ", + join(map(x -> ModelAnalyzer._name(x, model), issue.constraint), ", "), + ) end function ModelAnalyzer._verbose_summarize( io::IO, issue::InfeasibleBounds{T}, + model, ) where {T} return print( io, "Variable: ", - _name(issue.variable), + ModelAnalyzer._name(issue.variable, model), " with lower bound ", issue.lb, " and upper bound ", @@ -594,11 +721,12 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::InfeasibleIntegrality{T}, + model, ) where {T} return print( io, "Variable: ", - _name(issue.variable), + ModelAnalyzer._name(issue.variable, model), " with lower bound ", issue.lb, " and upper bound ", @@ -611,11 +739,12 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::InfeasibleConstraintRange{T}, + model, ) where {T} return print( io, "Constraint: ", - _name(issue.constraint), + ModelAnalyzer._name(issue.constraint, model), " with computed lower bound ", issue.lb, " and computed upper bound ", @@ -628,11 +757,12 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::IrreducibleInfeasibleSubset, + model, ) return print( io, "Irreducible Infeasible Subset: ", - join(map(_name, issue.constraint), ", "), + join(map(x -> ModelAnalyzer._name(x, model), issue.constraint), ", "), ) end @@ -676,6 +806,7 @@ end function ModelAnalyzer.summarize( io::IO, data::Data; + model = nothing, verbose = true, max_issues = ModelAnalyzer.DEFAULT_MAX_ISSUES, ) @@ -687,6 +818,7 @@ function ModelAnalyzer.summarize( ModelAnalyzer.summarize( io, issues, + model = model, verbose = verbose, max_issues = max_issues, ) @@ -703,14 +835,4 @@ function Base.show(io::IO, data::Data) return print(io, "Infeasibility analysis found $n issues") end -# printing helpers - -function _name(ref) - name = JuMP.name(ref) - if !isempty(name) - return name - end - return "$(ref.index)" -end - end # module diff --git a/src/intervals.jl b/src/intervals.jl index 1449316..a4c3881 100644 --- a/src/intervals.jl +++ b/src/intervals.jl @@ -54,6 +54,7 @@ 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 diff --git a/test/infeasibility.jl b/test/infeasibility.jl index a50fee6..761e62f 100644 --- a/test/infeasibility.jl +++ b/test/infeasibility.jl @@ -32,8 +32,11 @@ function test_bounds() @test length(list) == 1 ret = ModelAnalyzer.list_of_issues(data, list[1]) @test length(ret) == 1 - @test ret[] == - ModelAnalyzer.Infeasibility.InfeasibleBounds{Float64}(y, 2.0, 1.0) + @test ret[] == ModelAnalyzer.Infeasibility.InfeasibleBounds{Float64}( + JuMP.index(y), + 2.0, + 1.0, + ) # buf = IOBuffer() ModelAnalyzer.summarize( @@ -74,7 +77,7 @@ function test_integrality() ret = ModelAnalyzer.list_of_issues(data, list[1]) @test length(ret) == 1 @test ret[] == ModelAnalyzer.Infeasibility.InfeasibleIntegrality{Float64}( - y, + JuMP.index(y), 2.2, 2.9, MOI.Integer(), @@ -121,7 +124,7 @@ function test_binary() ret = ModelAnalyzer.list_of_issues(data, list[1]) @test length(ret) == 1 @test ret[] == ModelAnalyzer.Infeasibility.InfeasibleIntegrality{Float64}( - x, + JuMP.index(x), 0.5, 0.8, MOI.ZeroOne(), @@ -142,7 +145,7 @@ function test_range() @test length(ret) == 1 @test ret[] == ModelAnalyzer.Infeasibility.InfeasibleConstraintRange{Float64}( - c, + JuMP.index(c), 11.0, 22.0, MOI.LessThan{Float64}(1.0), @@ -190,7 +193,7 @@ function test_range_neg() @test length(ret) == 1 @test ret[] == ModelAnalyzer.Infeasibility.InfeasibleConstraintRange{Float64}( - c, + JuMP.index(c), 11.0, 22.0, MOI.LessThan{Float64}(1.0), @@ -211,7 +214,7 @@ function test_range_equalto() @test length(ret) == 1 @test ret[] == ModelAnalyzer.Infeasibility.InfeasibleConstraintRange{Float64}( - c, + JuMP.index(c), 3.0, 3.0, MOI.EqualTo{Float64}(1.0), @@ -232,7 +235,7 @@ function test_range_equalto_2() @test length(ret) == 1 @test ret[] == ModelAnalyzer.Infeasibility.InfeasibleConstraintRange{Float64}( - c, + JuMP.index(c), 7.0, 7.0, MOI.EqualTo{Float64}(1.0), @@ -253,7 +256,7 @@ function test_range_greaterthan() @test length(ret) == 1 @test ret[] == ModelAnalyzer.Infeasibility.InfeasibleConstraintRange{Float64}( - c, + JuMP.index(c), 11.0, 22.0, MOI.GreaterThan{Float64}(100.0), @@ -274,7 +277,7 @@ function test_range_equalto_3() @test length(ret) == 1 @test ret[] == ModelAnalyzer.Infeasibility.InfeasibleConstraintRange{Float64}( - c, + JuMP.index(c), 11.0, 22.0, MOI.EqualTo{Float64}(100.0), @@ -282,6 +285,8 @@ function test_range_equalto_3() return end +#= + function test_iis_feasible() model = Model(HiGHS.Optimizer) set_silent(model) @@ -321,7 +326,7 @@ function test_iis() ret = ModelAnalyzer.list_of_issues(data, list[1]) @test length(ret) == 1 @test length(ret[].constraint) == 2 - @test Set([ret[].constraint[1], ret[].constraint[2]]) == Set([c2, c1]) + @test Set([ret[].constraint[1], ret[].constraint[2]]) == Set(JuMP.index.([c2, c1])) # buf = IOBuffer() ModelAnalyzer.summarize( @@ -382,7 +387,7 @@ function test_iis_multiple() @test length(ret) == 1 @test length(ret[].constraint) == 2 @test c2 in Set([ret[].constraint[1], ret[].constraint[2]]) - @test Set([ret[].constraint[1], ret[].constraint[2]]) ⊆ Set([c3, c2, c1]) + @test Set([ret[].constraint[1], ret[].constraint[2]]) ⊆ Set(JuMP.index.([c3, c2, c1])) return end @@ -408,7 +413,7 @@ function test_iis_interval_right() ret = ModelAnalyzer.list_of_issues(data, list[1]) @test length(ret) == 1 @test length(ret[].constraint) == 2 - @test Set([ret[].constraint[1], ret[].constraint[2]]) == Set([c2, c1]) + @test Set([ret[].constraint[1], ret[].constraint[2]]) == Set(JuMP.index.([c2, c1])) return end @@ -434,7 +439,7 @@ function test_iis_interval_left() ret = ModelAnalyzer.list_of_issues(data, list[1]) @test length(ret) == 1 @test length(ret[].constraint) == 2 - @test Set([ret[].constraint[1], ret[].constraint[2]]) == Set([c2, c1]) + @test Set([ret[].constraint[1], ret[].constraint[2]]) == Set(JuMP.index.([c2, c1])) return end @@ -463,10 +468,12 @@ function test_iis_spare() ret = ModelAnalyzer.list_of_issues(data, list[1]) @test length(ret) == 1 @test length(ret[].constraint) == 2 - @test Set([ret[].constraint[1], ret[].constraint[2]]) == Set([c2, c1]) + @test Set([ret[].constraint[1], ret[].constraint[2]]) == Set(JuMP.index.([c2, c1])) return end +=# + end # module TestInfeasibilityChecker.runtests() diff --git a/test/runtests.jl b/test/runtests.jl index dadf665..12287cd 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,9 +6,7 @@ using Test @testset "ModelAnalyzer" begin - for file in ["feasibility.jl", "numerical.jl"] - # for file in ["infeasibility.jl"] - # for file in readdir(@__DIR__) + for file in readdir(@__DIR__) if !endswith(file, ".jl") || file in ("runtests.jl",) continue end From 38ae419715347df8023bdef6ed0f47d9ac7e6d30 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Thu, 8 May 2025 23:07:08 -0300 Subject: [PATCH 11/34] use MA --- src/_eval_variables.jl | 7 +------ src/intervals.jl | 4 ++++ 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/src/_eval_variables.jl b/src/_eval_variables.jl index 7c7d1d7..d30cd88 100644 --- a/src/_eval_variables.jl +++ b/src/_eval_variables.jl @@ -9,11 +9,6 @@ function _eval_variables(value_fn::Function, t::MOI.ScalarAffineTerm) return t.coefficient * value_fn(t.variable) end -function _eval_variables(value_fn::Function, t::MOI.ScalarQuadraticTerm) - out = t.coefficient * value_fn(t.variable_1) * value_fn(t.variable_2) - return t.variable_1 == t.variable_2 ? out / 2 : out -end - _eval_variables(value_fn::Function, f::MOI.VariableIndex) = value_fn(f) function _eval_variables( @@ -22,7 +17,7 @@ function _eval_variables( ) where {T} # TODO: this conversion exists in JuMP, but not in MOI S = Base.promote_op(value_fn, MOI.VariableIndex) - U = Base.promote_op(*, T, S) + U = MOI.MA.promote_operation(*, T, S) out = convert(U, f.constant) for t in f.terms out += _eval_variables(value_fn, t) diff --git a/src/intervals.jl b/src/intervals.jl index a4c3881..f54eb84 100644 --- a/src/intervals.jl +++ b/src/intervals.jl @@ -63,6 +63,10 @@ function Interval(lo::T, hi::T) where {T<:Real} 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 From 2e39b0ccfcf5bf17f4f15655067f1988d4e44d6c Mon Sep 17 00:00:00 2001 From: joaquimg Date: Tue, 13 May 2025 02:21:25 -0300 Subject: [PATCH 12/34] MOI IIS --- src/infeasibility.jl | 158 ++++++++++++++++++++++++------------------ test/infeasibility.jl | 6 +- 2 files changed, 91 insertions(+), 73 deletions(-) diff --git a/src/infeasibility.jl b/src/infeasibility.jl index 9232cc6..ffe3678 100644 --- a/src/infeasibility.jl +++ b/src/infeasibility.jl @@ -364,25 +364,57 @@ function ModelAnalyzer.analyze( 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 + 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 + 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 iis_elastic_filter(original_model::JuMP.GenericModel, optimizer) +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))) + else + MOI.add_constraint(model, variable, MOI.LessThan{T}(zero(T))) + end + return +end - # if JuMP.termination_status(original_model) == MOI.OPTIMIZE_NOT_CALLED - # println("iis resolver cannot continue because model is not optimized") - # # JuMP.optimize!(original_model) - # end +function iis_elastic_filter(original_model::MOI.ModelLike, optimizer) + T = Float64 - status = JuMP.termination_status(original_model) + # handle optimize not called + status = MOI.get(original_model, MOI.TerminationStatus()) if status != MOI.INFEASIBLE println( "iis resolver cannot continue because model is found to be $(status) by the solver", @@ -390,16 +422,14 @@ function iis_elastic_filter(original_model::JuMP.GenericModel, optimizer) return nothing end - model, reference_map = JuMP.copy_model(original_model) - JuMP.set_optimizer(model, optimizer) - JuMP.set_silent(model) - # TODO handle ".ext" to avoid warning + model = MOI.instantiate(optimizer) + reference_map = MOI.copy_to(model, original_model) + MOI.set(model, MOI.Silent(), true) - constraint_to_affine = JuMP.relax_with_penalty!(model, default = 1.0) + constraint_to_affine = + MOI.modify(model, MOI.Utilities.PenaltyRelaxation(default = 1.0)) # might need to do something related to integers / binary - JuMP.optimize!(model) - max_iterations = length(constraint_to_affine) tolerance = 1e-5 @@ -407,47 +437,45 @@ function iis_elastic_filter(original_model::JuMP.GenericModel, optimizer) de_elastisized = [] for i in 1:max_iterations - if JuMP.termination_status(model) == MOI.INFEASIBLE + MOI.optimize!(model) + status = MOI.get(model, MOI.TerminationStatus()) + if status == MOI.INFEASIBLE break end for (con, func) in constraint_to_affine if length(func.terms) == 1 - var = collect(keys(func.terms))[1] - if JuMP.value(var) > tolerance - has_lower = JuMP.has_lower_bound(var) - JuMP.fix(var, 0.0; force = true) - # or delete(model, var) + 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)) end elseif length(func.terms) == 2 - var = collect(keys(func.terms)) - coef1 = func.terms[var[1]] - coef2 = func.terms[var[2]] - if JuMP.value(var[1]) > tolerance && - JuMP.value(var[2]) > tolerance + 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 JuMP.value(var[1]) > tolerance - has_lower = JuMP.has_lower_bound(var[1]) - JuMP.fix(var[1], 0.0; force = true) - # or delete(model, var[1]) + elseif value1 > tolerance + has_lower = _fix_to_zero(model, var1, T) delete!(constraint_to_affine, con) - constraint_to_affine[con] = coef2 * var[2] - push!(de_elastisized, (con, var[1], has_lower)) - elseif JuMP.value(var[2]) > tolerance - has_lower = JuMP.has_lower_bound(var[2]) - JuMP.fix(var[2], 0.0; force = true) - # or delete(model, var[2]) + 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 * var[1] - push!(de_elastisized, (con, var[2], has_lower)) + 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 - JuMP.optimize!(model) end end @@ -455,46 +483,38 @@ function iis_elastic_filter(original_model::JuMP.GenericModel, optimizer) # be careful with intervals # deletion filter - cadidates = JuMP.ConstraintRef[] + cadidates = MOI.ConstraintIndex[] for (con, var, has_lower) in de_elastisized - JuMP.unfix(var) - if has_lower - JuMP.set_lower_bound(var, 0.0) - else - JuMP.set_upper_bound(var, 0.0) - end - JuMP.optimize!(model) - if JuMP.termination_status(model) in - (MOI.INFEASIBLE, MOI.ALMOST_INFEASIBLE) + _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) # this constraint is not in IIS - elseif JuMP.termination_status(model) in - (MOI.OPTIMAL, MOI.ALMOST_OPTIMAL) + elseif status in (MOI.OPTIMAL, MOI.ALMOST_OPTIMAL) push!(cadidates, con) - JuMP.fix(var, 0.0, force = true) + _fix_to_zero(model, var, T) else - error( - "IIS failed due numerical instability, got status $(JuMP.termination_status(model))", - ) + error("IIS failed due numerical instability, got status $status") end end pre_iis = Set(cadidates) - iis = JuMP.ConstraintRef[] - for con in JuMP.all_constraints( - original_model, - include_variable_in_set_constraints = false, - ) - new_con = reference_map[con] - if new_con in pre_iis - push!(iis, con) + 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 end -=# - # API function ModelAnalyzer._summarize(io::IO, ::Type{InfeasibleBounds{T}}) where {T} diff --git a/test/infeasibility.jl b/test/infeasibility.jl index 761e62f..5e7e8b3 100644 --- a/test/infeasibility.jl +++ b/test/infeasibility.jl @@ -285,7 +285,6 @@ function test_range_equalto_3() return end -#= function test_iis_feasible() model = Model(HiGHS.Optimizer) @@ -363,6 +362,7 @@ function test_iis() return end + function test_iis_multiple() model = Model(HiGHS.Optimizer) set_silent(model) @@ -386,7 +386,7 @@ function test_iis_multiple() ret = ModelAnalyzer.list_of_issues(data, list[1]) @test length(ret) == 1 @test length(ret[].constraint) == 2 - @test c2 in Set([ret[].constraint[1], ret[].constraint[2]]) + @test JuMP.index(c2) in Set([ret[].constraint[1], ret[].constraint[2]]) @test Set([ret[].constraint[1], ret[].constraint[2]]) ⊆ Set(JuMP.index.([c3, c2, c1])) return end @@ -472,8 +472,6 @@ function test_iis_spare() return end -=# - end # module TestInfeasibilityChecker.runtests() From 997200a81933656d27b52bffcacca8ca6a5cd8a5 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Tue, 13 May 2025 02:21:52 -0300 Subject: [PATCH 13/34] format --- test/infeasibility.jl | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/test/infeasibility.jl b/test/infeasibility.jl index 5e7e8b3..176c795 100644 --- a/test/infeasibility.jl +++ b/test/infeasibility.jl @@ -285,7 +285,6 @@ function test_range_equalto_3() return end - function test_iis_feasible() model = Model(HiGHS.Optimizer) set_silent(model) @@ -325,7 +324,8 @@ function test_iis() ret = ModelAnalyzer.list_of_issues(data, list[1]) @test length(ret) == 1 @test length(ret[].constraint) == 2 - @test Set([ret[].constraint[1], ret[].constraint[2]]) == Set(JuMP.index.([c2, c1])) + @test Set([ret[].constraint[1], ret[].constraint[2]]) == + Set(JuMP.index.([c2, c1])) # buf = IOBuffer() ModelAnalyzer.summarize( @@ -362,7 +362,6 @@ function test_iis() return end - function test_iis_multiple() model = Model(HiGHS.Optimizer) set_silent(model) @@ -387,7 +386,8 @@ function test_iis_multiple() @test length(ret) == 1 @test length(ret[].constraint) == 2 @test JuMP.index(c2) in Set([ret[].constraint[1], ret[].constraint[2]]) - @test Set([ret[].constraint[1], ret[].constraint[2]]) ⊆ Set(JuMP.index.([c3, c2, c1])) + @test Set([ret[].constraint[1], ret[].constraint[2]]) ⊆ + Set(JuMP.index.([c3, c2, c1])) return end @@ -413,7 +413,8 @@ function test_iis_interval_right() ret = ModelAnalyzer.list_of_issues(data, list[1]) @test length(ret) == 1 @test length(ret[].constraint) == 2 - @test Set([ret[].constraint[1], ret[].constraint[2]]) == Set(JuMP.index.([c2, c1])) + @test Set([ret[].constraint[1], ret[].constraint[2]]) == + Set(JuMP.index.([c2, c1])) return end @@ -439,7 +440,8 @@ function test_iis_interval_left() ret = ModelAnalyzer.list_of_issues(data, list[1]) @test length(ret) == 1 @test length(ret[].constraint) == 2 - @test Set([ret[].constraint[1], ret[].constraint[2]]) == Set(JuMP.index.([c2, c1])) + @test Set([ret[].constraint[1], ret[].constraint[2]]) == + Set(JuMP.index.([c2, c1])) return end @@ -468,7 +470,8 @@ function test_iis_spare() ret = ModelAnalyzer.list_of_issues(data, list[1]) @test length(ret) == 1 @test length(ret[].constraint) == 2 - @test Set([ret[].constraint[1], ret[].constraint[2]]) == Set(JuMP.index.([c2, c1])) + @test Set([ret[].constraint[1], ret[].constraint[2]]) == + Set(JuMP.index.([c2, c1])) return end From e059c41699f629de645cdfb804fddda65db883cc Mon Sep 17 00:00:00 2001 From: joaquimg Date: Tue, 13 May 2025 14:06:02 -0300 Subject: [PATCH 14/34] add query functions in numerical --- src/numerical.jl | 77 ++++++++++++++++++++++++++++++++--------------- test/numerical.jl | 60 +++++++++++++++++++++++++++++------- 2 files changed, 102 insertions(+), 35 deletions(-) diff --git a/src/numerical.jl b/src/numerical.jl index 66ff92b..ee1008b 100644 --- a/src/numerical.jl +++ b/src/numerical.jl @@ -61,7 +61,8 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.VariableNotInConstraints) struct VariableNotInConstraints <: AbstractNumericalIssue ref::MOI.VariableIndex end -ModelAnalyzer.variable(issue::VariableNotInConstraints, model) = issue.ref + +ModelAnalyzer.variable(issue::VariableNotInConstraints) = issue.ref """ EmptyConstraint <: AbstractNumericalIssue @@ -77,7 +78,8 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.EmptyConstraint) struct EmptyConstraint <: AbstractNumericalIssue ref::MOI.ConstraintIndex end -ModelAnalyzer.constraint(issue::EmptyConstraint, model) = issue.ref + +ModelAnalyzer.constraint(issue::EmptyConstraint) = issue.ref """ VariableBoundAsConstraint <: AbstractNumericalIssue @@ -94,7 +96,8 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.VariableBoundAsConstraint struct VariableBoundAsConstraint <: AbstractNumericalIssue ref::MOI.ConstraintIndex end -ModelAnalyzer.constraint(issue::VariableBoundAsConstraint, model) = issue.ref + +ModelAnalyzer.constraint(issue::VariableBoundAsConstraint) = issue.ref """ DenseConstraint <: AbstractNumericalIssue @@ -112,7 +115,9 @@ struct DenseConstraint <: AbstractNumericalIssue ref::MOI.ConstraintIndex nnz::Int end -ModelAnalyzer.constraint(issue::DenseConstraint, model) = issue.ref + +ModelAnalyzer.constraint(issue::DenseConstraint) = issue.ref + ModelAnalyzer.value(issue::DenseConstraint) = issue.nnz """ @@ -131,8 +136,11 @@ struct SmallMatrixCoefficient <: AbstractNumericalIssue variable::MOI.VariableIndex coefficient::Float64 end -ModelAnalyzer.variable(issue::SmallMatrixCoefficient, model) = issue.variable -ModelAnalyzer.constraint(issue::SmallMatrixCoefficient, model) = issue.ref + +ModelAnalyzer.variable(issue::SmallMatrixCoefficient) = issue.variable + +ModelAnalyzer.constraint(issue::SmallMatrixCoefficient) = issue.ref + ModelAnalyzer.value(issue::SmallMatrixCoefficient) = issue.coefficient """ @@ -151,8 +159,11 @@ struct LargeMatrixCoefficient <: AbstractNumericalIssue variable::MOI.VariableIndex coefficient::Float64 end -ModelAnalyzer.variable(issue::LargeMatrixCoefficient, model) = issue.variable -ModelAnalyzer.constraint(issue::LargeMatrixCoefficient, model) = issue.ref + +ModelAnalyzer.variable(issue::LargeMatrixCoefficient) = issue.variable + +ModelAnalyzer.constraint(issue::LargeMatrixCoefficient) = issue.ref + ModelAnalyzer.value(issue::LargeMatrixCoefficient) = issue.coefficient """ @@ -170,7 +181,9 @@ struct SmallBoundCoefficient <: AbstractNumericalIssue variable::MOI.VariableIndex coefficient::Float64 end -ModelAnalyzer.variable(issue::SmallBoundCoefficient, model) = issue.variable + +ModelAnalyzer.variable(issue::SmallBoundCoefficient) = issue.variable + ModelAnalyzer.value(issue::SmallBoundCoefficient) = issue.coefficient """ @@ -188,7 +201,9 @@ struct LargeBoundCoefficient <: AbstractNumericalIssue variable::MOI.VariableIndex coefficient::Float64 end -ModelAnalyzer.variable(issue::LargeBoundCoefficient, model) = issue.variable + +ModelAnalyzer.variable(issue::LargeBoundCoefficient) = issue.variable + ModelAnalyzer.value(issue::LargeBoundCoefficient) = issue.coefficient """ @@ -206,7 +221,9 @@ struct SmallRHSCoefficient <: AbstractNumericalIssue ref::MOI.ConstraintIndex coefficient::Float64 end -ModelAnalyzer.constraint(issue::SmallRHSCoefficient, model) = issue.ref + +ModelAnalyzer.constraint(issue::SmallRHSCoefficient) = issue.ref + ModelAnalyzer.value(issue::SmallRHSCoefficient) = issue.coefficient """ @@ -224,7 +241,9 @@ struct LargeRHSCoefficient <: AbstractNumericalIssue ref::MOI.ConstraintIndex coefficient::Float64 end -ModelAnalyzer.constraint(issue::LargeRHSCoefficient, model) = issue.ref + +ModelAnalyzer.constraint(issue::LargeRHSCoefficient) = issue.ref + ModelAnalyzer.value(issue::LargeRHSCoefficient) = issue.coefficient """ @@ -242,7 +261,9 @@ struct SmallObjectiveCoefficient <: AbstractNumericalIssue variable::MOI.VariableIndex coefficient::Float64 end -ModelAnalyzer.variable(issue::SmallObjectiveCoefficient, model) = issue.variable + +ModelAnalyzer.variable(issue::SmallObjectiveCoefficient) = issue.variable + ModelAnalyzer.value(issue::SmallObjectiveCoefficient) = issue.coefficient """ @@ -260,7 +281,9 @@ struct LargeObjectiveCoefficient <: AbstractNumericalIssue variable::MOI.VariableIndex coefficient::Float64 end -ModelAnalyzer.variable(issue::LargeObjectiveCoefficient, model) = issue.variable + +ModelAnalyzer.variable(issue::LargeObjectiveCoefficient) = issue.variable + ModelAnalyzer.value(issue::LargeObjectiveCoefficient) = issue.coefficient """ @@ -281,12 +304,11 @@ struct SmallObjectiveQuadraticCoefficient <: AbstractNumericalIssue variable2::MOI.VariableIndex coefficient::Float64 end -function ModelAnalyzer.variables( - issue::SmallObjectiveQuadraticCoefficient, - model, -) + +function ModelAnalyzer.variables(issue::SmallObjectiveQuadraticCoefficient) return [issue.variable1, issue.variable2] end + function ModelAnalyzer.value(issue::SmallObjectiveQuadraticCoefficient) return issue.coefficient end @@ -309,12 +331,13 @@ struct LargeObjectiveQuadraticCoefficient <: AbstractNumericalIssue variable2::MOI.VariableIndex coefficient::Float64 end + function ModelAnalyzer.variables( issue::LargeObjectiveQuadraticCoefficient, - model, ) return [issue.variable1, issue.variable2] end + function ModelAnalyzer.value(issue::LargeObjectiveQuadraticCoefficient) return issue.coefficient end @@ -338,12 +361,15 @@ struct SmallMatrixQuadraticCoefficient <: AbstractNumericalIssue variable2::MOI.VariableIndex coefficient::Float64 end -function ModelAnalyzer.variables(issue::SmallMatrixQuadraticCoefficient, model) + +function ModelAnalyzer.variables(issue::SmallMatrixQuadraticCoefficient) return [issue.variable1, issue.variable2] end -function ModelAnalyzer.constraint(issue::SmallMatrixQuadraticCoefficient, model) + +function ModelAnalyzer.constraint(issue::SmallMatrixQuadraticCoefficient) return issue.ref end + ModelAnalyzer.value(issue::SmallMatrixQuadraticCoefficient) = issue.coefficient """ @@ -365,12 +391,15 @@ struct LargeMatrixQuadraticCoefficient <: AbstractNumericalIssue variable2::MOI.VariableIndex coefficient::Float64 end -function ModelAnalyzer.variables(issue::LargeMatrixQuadraticCoefficient, model) + +function ModelAnalyzer.variables(issue::LargeMatrixQuadraticCoefficient) return [issue.variable1, issue.variable2] end -function ModelAnalyzer.constraint(issue::LargeMatrixQuadraticCoefficient, model) + +function ModelAnalyzer.constraint(issue::LargeMatrixQuadraticCoefficient) return issue.ref end + ModelAnalyzer.value(issue::LargeMatrixQuadraticCoefficient) = issue.coefficient """ @@ -404,7 +433,7 @@ julia> ModelAnalyzer.summarize( struct NonconvexQuadraticConstraint <: AbstractNumericalIssue ref::MOI.ConstraintIndex end -ModelAnalyzer.constraint(issue::NonconvexQuadraticConstraint, model) = issue.ref +ModelAnalyzer.constraint(issue::NonconvexQuadraticConstraint) = issue.ref """ Data diff --git a/test/numerical.jl b/test/numerical.jl index ed8f1d9..62c1e85 100644 --- a/test/numerical.jl +++ b/test/numerical.jl @@ -205,6 +205,8 @@ function test_variable_not_in_constraints() ModelAnalyzer.Numerical.VariableNotInConstraints, ) @test length(ret) == 1 + @test ModelAnalyzer.variable(ret[]) == JuMP.index(x) + @test ModelAnalyzer.variable(ret[], model) == x # buf = IOBuffer() ModelAnalyzer.summarize( @@ -232,8 +234,8 @@ end function test_empty_constraint_model() model = Model() @variable(model, x) - @constraint(model, 2 * x == 5) - @constraint(model, 0.0 * x == 3) + @constraint(model, c1, 2 * x == 5) + @constraint(model, c2, 0.0 * x == 3) data = ModelAnalyzer.analyze(ModelAnalyzer.Numerical.Analyzer(), model) list = ModelAnalyzer.list_of_issue_types(data) @test length(list) == 1 @@ -242,6 +244,8 @@ function test_empty_constraint_model() ModelAnalyzer.Numerical.EmptyConstraint, ) @test length(ret) == 1 + @test ModelAnalyzer.constraint(ret[]) == JuMP.index(c2) + @test ModelAnalyzer.constraint(ret[], model) == c2 # buf = IOBuffer() ModelAnalyzer.summarize(buf, ModelAnalyzer.Numerical.EmptyConstraint) @@ -266,7 +270,7 @@ end function test_variable_bound_as_constraint() model = Model() @variable(model, x) - @constraint(model, x <= 2) + @constraint(model, c, x <= 2) @constraint(model, 3x <= 4) data = ModelAnalyzer.analyze(ModelAnalyzer.Numerical.Analyzer(), model) list = ModelAnalyzer.list_of_issue_types(data) @@ -276,6 +280,8 @@ function test_variable_bound_as_constraint() ModelAnalyzer.Numerical.VariableBoundAsConstraint, ) @test length(ret) == 1 + @test ModelAnalyzer.constraint(ret[]) == JuMP.index(c) + @test ModelAnalyzer.constraint(ret[], model) == c # buf = IOBuffer() ModelAnalyzer.summarize( @@ -303,7 +309,7 @@ end function test_dense_constraint() model = Model() @variable(model, x[1:10_000] <= 1) - @constraint(model, sum(x) <= 4) + @constraint(model, c, sum(x) <= 4) data = ModelAnalyzer.analyze(ModelAnalyzer.Numerical.Analyzer(), model) list = ModelAnalyzer.list_of_issue_types(data) @test length(list) == 1 @@ -312,6 +318,9 @@ function test_dense_constraint() ModelAnalyzer.Numerical.DenseConstraint, ) @test length(ret) == 1 + @test ModelAnalyzer.constraint(ret[]) == JuMP.index(c) + @test ModelAnalyzer.constraint(ret[], model) == c + @test ModelAnalyzer.value(ret[]) == 10_000 # buf = IOBuffer() ModelAnalyzer.summarize(buf, ModelAnalyzer.Numerical.DenseConstraint) @@ -337,7 +346,7 @@ end function test_small_matrix_coef() model = Model() @variable(model, x <= 1) - @constraint(model, 1e-9 * x <= 4) + @constraint(model, c, 1e-9 * x <= 4) data = ModelAnalyzer.analyze(ModelAnalyzer.Numerical.Analyzer(), model) list = ModelAnalyzer.list_of_issue_types(data) @test length(list) == 1 @@ -346,6 +355,9 @@ function test_small_matrix_coef() ModelAnalyzer.Numerical.SmallMatrixCoefficient, ) @test length(ret) == 1 + @test ModelAnalyzer.constraint(ret[], model) == c + @test ModelAnalyzer.variable(ret[], model) == x + @test ModelAnalyzer.value(ret[]) == 1e-9 # buf = IOBuffer() ModelAnalyzer.summarize(buf, ModelAnalyzer.Numerical.SmallMatrixCoefficient) @@ -373,7 +385,7 @@ end function test_large_matrix_coef() model = Model() @variable(model, x <= 1) - @constraint(model, 1e+9 * x <= 4) + @constraint(model, c, 1e+9 * x <= 4) data = ModelAnalyzer.analyze(ModelAnalyzer.Numerical.Analyzer(), model) list = ModelAnalyzer.list_of_issue_types(data) @test length(list) == 1 @@ -382,6 +394,9 @@ function test_large_matrix_coef() ModelAnalyzer.Numerical.LargeMatrixCoefficient, ) @test length(ret) == 1 + @test ModelAnalyzer.constraint(ret[], model) == c + @test ModelAnalyzer.variable(ret[], model) == x + @test ModelAnalyzer.value(ret[]) == 1e+9 # buf = IOBuffer() ModelAnalyzer.summarize(buf, ModelAnalyzer.Numerical.LargeMatrixCoefficient) @@ -418,6 +433,8 @@ function test_small_bound_coef() ModelAnalyzer.Numerical.SmallBoundCoefficient, ) @test length(ret) == 1 + @test ModelAnalyzer.variable(ret[], model) == x + @test ModelAnalyzer.value(ret[]) == 1e-9 # buf = IOBuffer() ModelAnalyzer.summarize(buf, ModelAnalyzer.Numerical.SmallBoundCoefficient) @@ -453,6 +470,8 @@ function test_large_bound_coef() ModelAnalyzer.Numerical.LargeBoundCoefficient, ) @test length(ret) == 1 + @test ModelAnalyzer.variable(ret[], model) == x + @test ModelAnalyzer.value(ret[]) == 1e+9 # buf = IOBuffer() ModelAnalyzer.summarize(buf, ModelAnalyzer.Numerical.LargeBoundCoefficient) @@ -479,7 +498,7 @@ end function test_small_rhs_coef() model = Model() @variable(model, x <= 1) - @constraint(model, 3 * x <= 1e-9) + @constraint(model, c, 3 * x <= 1e-9) data = ModelAnalyzer.analyze(ModelAnalyzer.Numerical.Analyzer(), model) list = ModelAnalyzer.list_of_issue_types(data) @test length(list) == 1 @@ -488,6 +507,8 @@ function test_small_rhs_coef() ModelAnalyzer.Numerical.SmallRHSCoefficient, ) @test length(ret) == 1 + @test ModelAnalyzer.constraint(ret[], model) == c + @test ModelAnalyzer.value(ret[]) == 1e-9 # buf = IOBuffer() ModelAnalyzer.summarize(buf, ModelAnalyzer.Numerical.SmallRHSCoefficient) @@ -514,7 +535,7 @@ end function test_large_rhs_coef() model = Model() @variable(model, x <= 1) - @constraint(model, 3 * x <= 1e+9) + @constraint(model, c, 3 * x <= 1e+9) data = ModelAnalyzer.analyze(ModelAnalyzer.Numerical.Analyzer(), model) list = ModelAnalyzer.list_of_issue_types(data) @test length(list) == 1 @@ -523,6 +544,8 @@ function test_large_rhs_coef() ModelAnalyzer.Numerical.LargeRHSCoefficient, ) @test length(ret) == 1 + @test ModelAnalyzer.constraint(ret[], model) == c + @test ModelAnalyzer.value(ret[]) == 1e+9 # buf = IOBuffer() ModelAnalyzer.summarize(buf, ModelAnalyzer.Numerical.LargeRHSCoefficient) @@ -559,6 +582,8 @@ function test_small_objective_coef() ModelAnalyzer.Numerical.SmallObjectiveCoefficient, ) @test length(ret) == 1 + @test ModelAnalyzer.variable(ret[], model) == x + @test ModelAnalyzer.value(ret[]) == 1e-9 # buf = IOBuffer() ModelAnalyzer.summarize( @@ -598,6 +623,8 @@ function test_large_objective_coef() ModelAnalyzer.Numerical.LargeObjectiveCoefficient, ) @test length(ret) == 1 + @test ModelAnalyzer.variable(ret[], model) == x + @test ModelAnalyzer.value(ret[]) == 1e+9 # buf = IOBuffer() ModelAnalyzer.summarize( @@ -637,6 +664,8 @@ function test_small_objective_coef_quad() ModelAnalyzer.Numerical.SmallObjectiveQuadraticCoefficient, ) @test length(ret) == 1 + @test ModelAnalyzer.variables(ret[], model) == [x, x] + @test_broken ModelAnalyzer.value(ret[]) == 1e-9 # 2e-9 TODO, what to return here # buf = IOBuffer() ModelAnalyzer.summarize( @@ -677,6 +706,8 @@ function test_large_objective_coef_quad() ModelAnalyzer.Numerical.LargeObjectiveQuadraticCoefficient, ) @test length(ret) == 1 + @test ModelAnalyzer.variables(ret[], model) == [x, x] + @test_broken ModelAnalyzer.value(ret[]) == 1e+9 # 2e+9 TODO, what to return here # buf = IOBuffer() ModelAnalyzer.summarize( @@ -707,7 +738,7 @@ end function test_small_matrix_coef_quad() model = Model() @variable(model, x <= 1) - @constraint(model, 1e-9 * x^2 + x <= 4) + @constraint(model, c, 1e-9 * x^2 + x <= 4) data = ModelAnalyzer.analyze(ModelAnalyzer.Numerical.Analyzer(), model) list = ModelAnalyzer.list_of_issue_types(data) @test length(list) == 1 @@ -716,6 +747,9 @@ function test_small_matrix_coef_quad() ModelAnalyzer.Numerical.SmallMatrixQuadraticCoefficient, ) @test length(ret) == 1 + @test ModelAnalyzer.variables(ret[], model) == [x, x] + @test ModelAnalyzer.constraint(ret[], model) == c + @test_broken ModelAnalyzer.value(ret[]) == 1e-9 # 2e-9 TODO, what to return here # buf = IOBuffer() ModelAnalyzer.summarize( @@ -746,7 +780,7 @@ end function test_large_matrix_coef_quad() model = Model() @variable(model, x <= 1) - @constraint(model, 1e+9 * x^2 <= 4) + @constraint(model, c, 1e+9 * x^2 <= 4) data = ModelAnalyzer.analyze(ModelAnalyzer.Numerical.Analyzer(), model) list = ModelAnalyzer.list_of_issue_types(data) @test length(list) == 1 @@ -755,6 +789,9 @@ function test_large_matrix_coef_quad() ModelAnalyzer.Numerical.LargeMatrixQuadraticCoefficient, ) @test length(ret) == 1 + @test ModelAnalyzer.variables(ret[], model) == [x, x] + @test ModelAnalyzer.constraint(ret[], model) == c + @test_broken ModelAnalyzer.value(ret[]) == 1e+9 # 2e+9 TODO, what to return here # buf = IOBuffer() ModelAnalyzer.summarize( @@ -861,7 +898,7 @@ end function test_constraint_nonconvex() model = Model() @variable(model, x <= 1) - @constraint(model, x^2 >= 4) + @constraint(model, c, x^2 >= 4) data = ModelAnalyzer.analyze(ModelAnalyzer.Numerical.Analyzer(), model) list = ModelAnalyzer.list_of_issue_types(data) @test length(list) == 1 @@ -870,6 +907,7 @@ function test_constraint_nonconvex() ModelAnalyzer.Numerical.NonconvexQuadraticConstraint, ) @test length(ret) == 1 + @test ModelAnalyzer.constraint(ret[], model) == c # buf = IOBuffer() ModelAnalyzer.summarize( From fc09db07c575274966adb35660d9419612c7c279 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Tue, 13 May 2025 14:06:27 -0300 Subject: [PATCH 15/34] format --- src/numerical.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/numerical.jl b/src/numerical.jl index ee1008b..8eb904a 100644 --- a/src/numerical.jl +++ b/src/numerical.jl @@ -332,9 +332,7 @@ struct LargeObjectiveQuadraticCoefficient <: AbstractNumericalIssue coefficient::Float64 end -function ModelAnalyzer.variables( - issue::LargeObjectiveQuadraticCoefficient, -) +function ModelAnalyzer.variables(issue::LargeObjectiveQuadraticCoefficient) return [issue.variable1, issue.variable2] end From f0d70fa57c448f104291f78f1b7ed807d6b05188 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Tue, 13 May 2025 15:23:23 -0300 Subject: [PATCH 16/34] Add query tools --- src/feasibility.jl | 24 +++++++++++++++++++++ test/feasibility.jl | 52 +++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 76 insertions(+) diff --git a/src/feasibility.jl b/src/feasibility.jl index 280823a..2a0a952 100644 --- a/src/feasibility.jl +++ b/src/feasibility.jl @@ -66,6 +66,10 @@ struct PrimalViolation <: AbstractFeasibilityIssue violation::Float64 end +ModelAnalyzer.constraint(issue::PrimalViolation) = issue.ref + +ModelAnalyzer.value(issue::PrimalViolation) = issue.violation + """ DualConstraintViolation <: AbstractFeasibilityIssue @@ -83,6 +87,10 @@ struct DualConstraintViolation <: AbstractFeasibilityIssue violation::Float64 end +ModelAnalyzer.variable(issue::DualConstraintViolation) = issue.ref + +ModelAnalyzer.value(issue::DualConstraintViolation) = issue.violation + """ DualConstrainedVariableViolation <: AbstractFeasibilityIssue @@ -101,6 +109,10 @@ struct DualConstrainedVariableViolation <: AbstractFeasibilityIssue violation::Float64 end +ModelAnalyzer.constraint(issue::DualConstrainedVariableViolation) = issue.ref + +ModelAnalyzer.value(issue::DualConstrainedVariableViolation) = issue.violation + """ ComplemetarityViolation <: AbstractFeasibilityIssue @@ -119,6 +131,10 @@ struct ComplemetarityViolation <: AbstractFeasibilityIssue violation::Float64 end +ModelAnalyzer.constraint(issue::ComplemetarityViolation) = issue.ref + +ModelAnalyzer.value(issue::ComplemetarityViolation) = issue.violation + """ DualObjectiveMismatch <: AbstractFeasibilityIssue @@ -136,6 +152,8 @@ struct DualObjectiveMismatch <: AbstractFeasibilityIssue obj_solver::Float64 end +# ModelAnalyzer.values(issue::DualObjectiveMismatch) = [issue.obj, issue.obj_solver] + """ PrimalObjectiveMismatch <: AbstractFeasibilityIssue @@ -153,6 +171,8 @@ struct PrimalObjectiveMismatch <: AbstractFeasibilityIssue obj_solver::Float64 end +# ModelAnalyzer.values(issue::PrimalObjectiveMismatch) = [issue.obj, issue.obj_solver] + """ PrimalDualMismatch <: AbstractFeasibilityIssue @@ -170,6 +190,8 @@ struct PrimalDualMismatch <: AbstractFeasibilityIssue dual::Float64 end +ModelAnalyzer.values(issue::PrimalDualMismatch) = [issue.primal, issue.dual] + """ PrimalDualSolverMismatch <: AbstractFeasibilityIssue @@ -187,6 +209,8 @@ struct PrimalDualSolverMismatch <: AbstractFeasibilityIssue dual::Float64 end +# ModelAnalyzer.values(issue::PrimalDualSolverMismatch) = [issue.primal, issue.dual] + """ Data diff --git a/test/feasibility.jl b/test/feasibility.jl index cbeda18..797527e 100644 --- a/test/feasibility.jl +++ b/test/feasibility.jl @@ -92,6 +92,8 @@ function test_only_bounds() JuMP.index(LowerBoundRef(x)), 1.0, ) + @test ModelAnalyzer.constraint(ret[], model) == LowerBoundRef(x) + @test ModelAnalyzer.value(ret[]) == 1.0 ret = ModelAnalyzer.list_of_issues( data, ModelAnalyzer.Feasibility.ComplemetarityViolation, @@ -100,11 +102,14 @@ function test_only_bounds() JuMP.index(LowerBoundRef(x)), -1.0, ) + @test ModelAnalyzer.constraint(ret[], model) == LowerBoundRef(x) + @test ModelAnalyzer.value(ret[]) == -1.0 ret = ModelAnalyzer.list_of_issues( data, ModelAnalyzer.Feasibility.PrimalDualMismatch, ) @test ret[] == ModelAnalyzer.Feasibility.PrimalDualMismatch(-1.0, 0.0) + # @test ModelAnalyzer.values(ret[]) == [-1.0, 0.0] data = ModelAnalyzer.analyze( ModelAnalyzer.Feasibility.Analyzer(), @@ -125,6 +130,8 @@ function test_only_bounds() JuMP.index(LowerBoundRef(x)), 1.0, ) + @test ModelAnalyzer.constraint(ret[], model) == LowerBoundRef(x) + @test ModelAnalyzer.value(ret[]) == 1.0 ret = ModelAnalyzer.list_of_issues( data, ModelAnalyzer.Feasibility.PrimalDualMismatch, @@ -150,6 +157,8 @@ function test_only_bounds() ) @test ret[] == ModelAnalyzer.Feasibility.DualConstraintViolation(JuMP.index(x), 2.0) + @test ModelAnalyzer.variable(ret[], model) == x + @test ModelAnalyzer.value(ret[]) == 2.0 ret = ModelAnalyzer.list_of_issues( data, ModelAnalyzer.Feasibility.DualConstrainedVariableViolation, @@ -158,6 +167,8 @@ function test_only_bounds() JuMP.index(LowerBoundRef(x)), 1.0, ) + @test ModelAnalyzer.constraint(ret[], model) == LowerBoundRef(x) + @test ModelAnalyzer.value(ret[]) == 1.0 ret = ModelAnalyzer.list_of_issues( data, ModelAnalyzer.Feasibility.ComplemetarityViolation, @@ -166,6 +177,8 @@ function test_only_bounds() JuMP.index(LowerBoundRef(x)), -1.0, ) + @test ModelAnalyzer.constraint(ret[], model) == LowerBoundRef(x) + @test ModelAnalyzer.value(ret[]) == -1.0 ret = ModelAnalyzer.list_of_issues( data, ModelAnalyzer.Feasibility.PrimalDualMismatch, @@ -213,12 +226,16 @@ function test_no_lb() ModelAnalyzer.Feasibility.PrimalViolation, ) @test ret[] == ModelAnalyzer.Feasibility.PrimalViolation(JuMP.index(c), 1.0) + @test ModelAnalyzer.constraint(ret[], model) == c + @test ModelAnalyzer.value(ret[]) == 1.0 ret = ModelAnalyzer.list_of_issues( data, ModelAnalyzer.Feasibility.ComplemetarityViolation, ) @test ret[] == ModelAnalyzer.Feasibility.ComplemetarityViolation(JuMP.index(c), -1.0) + @test ModelAnalyzer.constraint(ret[], model) == c + @test ModelAnalyzer.value(ret[]) == -1.0 ret = ModelAnalyzer.list_of_issues( data, ModelAnalyzer.Feasibility.PrimalDualMismatch, @@ -242,6 +259,8 @@ function test_no_lb() ) @test ret[] == ModelAnalyzer.Feasibility.ComplemetarityViolation(JuMP.index(c), 1.0) + @test ModelAnalyzer.constraint(ret[], model) == c + @test ModelAnalyzer.value(ret[]) == 1.0 ret = ModelAnalyzer.list_of_issues( data, ModelAnalyzer.Feasibility.PrimalDualMismatch, @@ -267,6 +286,8 @@ function test_no_lb() ) @test ret[] == ModelAnalyzer.Feasibility.DualConstraintViolation(JuMP.index(x), 2.0) + @test ModelAnalyzer.variable(ret[], model) == x + @test ModelAnalyzer.value(ret[]) == 2.0 ret = ModelAnalyzer.list_of_issues( data, ModelAnalyzer.Feasibility.DualConstrainedVariableViolation, @@ -275,12 +296,16 @@ function test_no_lb() JuMP.index(c), 1.0, ) + @test ModelAnalyzer.constraint(ret[], model) == c + @test ModelAnalyzer.value(ret[]) == 1.0 ret = ModelAnalyzer.list_of_issues( data, ModelAnalyzer.Feasibility.ComplemetarityViolation, ) @test ret[] == ModelAnalyzer.Feasibility.ComplemetarityViolation(JuMP.index(c), -1.0) + @test ModelAnalyzer.constraint(ret[], model) == c + @test ModelAnalyzer.value(ret[]) == -1.0 ret = ModelAnalyzer.list_of_issues( data, ModelAnalyzer.Feasibility.PrimalDualMismatch, @@ -329,6 +354,8 @@ function test_lb0() ) @test ret[] == ModelAnalyzer.Feasibility.DualConstraintViolation(JuMP.index(x), 2.3) + @test ModelAnalyzer.variable(ret[], model) == x + @test ModelAnalyzer.value(ret[]) == 2.3 data = ModelAnalyzer.analyze( ModelAnalyzer.Feasibility.Analyzer(), @@ -351,6 +378,8 @@ function test_lb0() ) @test ret[] == ModelAnalyzer.Feasibility.DualConstraintViolation(JuMP.index(x), 4.3) + @test ModelAnalyzer.variable(ret[], model) == x + @test ModelAnalyzer.value(ret[]) == 4.3 ret = ModelAnalyzer.list_of_issues( data, ModelAnalyzer.Feasibility.DualConstrainedVariableViolation, @@ -359,6 +388,9 @@ function test_lb0() JuMP.index(c), 3.3, ) + @test ModelAnalyzer.constraint(ret[], model) == c + @test ModelAnalyzer.value(ret[]) == 3.3 + return end function test_lb2() @@ -398,6 +430,8 @@ function test_lb2() ) @test ret[] == ModelAnalyzer.Feasibility.DualConstraintViolation(JuMP.index(x), 2.3) + @test ModelAnalyzer.variable(ret[], model) == x + @test ModelAnalyzer.value(ret[]) == 2.3 data = ModelAnalyzer.analyze( ModelAnalyzer.Feasibility.Analyzer(), @@ -422,6 +456,8 @@ function test_lb2() JuMP.index(c), 3.3, ) + @test ModelAnalyzer.constraint(ret[], model) == c + @test ModelAnalyzer.value(ret[]) == 3.3 data = ModelAnalyzer.analyze( ModelAnalyzer.Feasibility.Analyzer(), @@ -446,10 +482,14 @@ function test_lb2() JuMP.index(c), 3.3, ) + @test ModelAnalyzer.constraint(ret[1], model) == c + @test ModelAnalyzer.value(ret[1]) == 3.3 @test ret[2] == ModelAnalyzer.Feasibility.DualConstrainedVariableViolation( JuMP.index(LowerBoundRef(x)), 1.0, ) + @test ModelAnalyzer.constraint(ret[2], model) == LowerBoundRef(x) + @test ModelAnalyzer.value(ret[2]) == 1.0 data = ModelAnalyzer.analyze( ModelAnalyzer.Feasibility.Analyzer(), @@ -470,6 +510,8 @@ function test_lb2() JuMP.index(c), 3.3, ) + @test ModelAnalyzer.constraint(ret[], model) == c + @test ModelAnalyzer.value(ret[]) == 3.3 return end @@ -613,6 +655,8 @@ function test_analyse_no_opt() @test length(list) == 1 ret = ModelAnalyzer.list_of_issues(data, list[1]) @test ret[] == ModelAnalyzer.Feasibility.PrimalViolation(JuMP.index(c), 1.0) + @test ModelAnalyzer.constraint(ret[], model) == c + @test ModelAnalyzer.value(ret[]) == 1.0 data = ModelAnalyzer.analyze( ModelAnalyzer.Feasibility.Analyzer(), @@ -625,6 +669,8 @@ function test_analyse_no_opt() ret = ModelAnalyzer.list_of_issues(data, list[1]) @test ret[1] == ModelAnalyzer.Feasibility.ComplemetarityViolation(JuMP.index(c), 1.0) + @test ModelAnalyzer.constraint(ret[1], model) == c + @test ModelAnalyzer.value(ret[1]) == 1.0 ret = ModelAnalyzer.list_of_issues(data, list[2]) @test ret[1] == ModelAnalyzer.Feasibility.PrimalDualMismatch(1.0, 0.0) @@ -648,12 +694,18 @@ function test_analyse_no_opt() ret = ModelAnalyzer.list_of_issues(data, list[1]) @test ret[1] == ModelAnalyzer.Feasibility.PrimalViolation(JuMP.index(c), 1.0) + @test ModelAnalyzer.constraint(ret[1], model) == c + @test ModelAnalyzer.value(ret[1]) == 1.0 ret = ModelAnalyzer.list_of_issues(data, list[2]) @test ret[1] == ModelAnalyzer.Feasibility.DualConstraintViolation(JuMP.index(x), 1.0) + @test ModelAnalyzer.variable(ret[1], model) == x + @test ModelAnalyzer.value(ret[1]) == 1.0 ret = ModelAnalyzer.list_of_issues(data, list[3]) @test ret[1] == ModelAnalyzer.Feasibility.ComplemetarityViolation(JuMP.index(c), -2.0) + @test ModelAnalyzer.constraint(ret[1], model) == c + @test ModelAnalyzer.value(ret[1]) == -2.0 ret = ModelAnalyzer.list_of_issues(data, list[4]) @test ret[1] == ModelAnalyzer.Feasibility.PrimalDualMismatch(-1.0, 0.0) From 833e6c1a9a7a102c7644218fd618e4467c6f57c5 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Tue, 13 May 2025 17:50:11 -0300 Subject: [PATCH 17/34] Add query tools --- .../ModelAnalyzerJuMPExt.jl | 13 ++++++ src/ModelAnalyzer.jl | 39 +++++++++++++++++ src/infeasibility.jl | 18 ++++++++ test/infeasibility.jl | 42 +++++++++++++++++++ 4 files changed, 112 insertions(+) diff --git a/ext/ModelAnalyzerJuMPExt/ModelAnalyzerJuMPExt.jl b/ext/ModelAnalyzerJuMPExt/ModelAnalyzerJuMPExt.jl index a5eb50a..25eafae 100644 --- a/ext/ModelAnalyzerJuMPExt/ModelAnalyzerJuMPExt.jl +++ b/ext/ModelAnalyzerJuMPExt/ModelAnalyzerJuMPExt.jl @@ -79,4 +79,17 @@ function ModelAnalyzer.constraint( return JuMP.constraint_ref_with_index(model, ref) end +""" + constraintss(issue::ModelAnalyzer.AbstractIssue, model::JuMP.GenericModel) + +Return the **JuMP** constraints reference associated to a particular issue. +""" +function ModelAnalyzer.constraints( + issue::ModelAnalyzer.AbstractIssue, + model::JuMP.GenericModel, +) + ref = ModelAnalyzer.constraints(issue) + return JuMP.constraint_ref_with_index.(model, ref) +end + end # module ModelAnalyzerJuMPExt diff --git a/src/ModelAnalyzer.jl b/src/ModelAnalyzer.jl index 22b2637..a7302c1 100644 --- a/src/ModelAnalyzer.jl +++ b/src/ModelAnalyzer.jl @@ -151,6 +151,19 @@ function value(issue::AbstractIssue, ::MOI.ModelLike) return value(issue) end +""" + values(issue::AbstractIssue) + +Return the values associated to a particular issue. +""" +function values(issue::AbstractIssue, ::Nothing) + return values(issue) +end + +function values(issue::AbstractIssue, ::MOI.ModelLike) + return values(issue) +end + """ variable(issue::AbstractIssue) @@ -190,6 +203,32 @@ function constraint(issue::AbstractIssue, ::MOI.ModelLike) return constraint(issue) end +""" + constraints(issue::AbstractIssue) + +Return the constraints associated to a particular issue. +""" +function constraints(issue::AbstractIssue, ::Nothing) + return constraints(issue) +end + +function constraints(issue::AbstractIssue, ::MOI.ModelLike) + return constraints(issue) +end + +""" + set(issue::AbstractIssue) + +Return the set associated to a particular issue. +""" +function set(issue::AbstractIssue, ::Nothing) + return set(issue) +end + +function set(issue::AbstractIssue, ::MOI.ModelLike) + return set(issue) +end + function _verbose_summarize end function _summarize end diff --git a/src/infeasibility.jl b/src/infeasibility.jl index ffe3678..7bee7f1 100644 --- a/src/infeasibility.jl +++ b/src/infeasibility.jl @@ -55,6 +55,10 @@ struct InfeasibleBounds{T} <: AbstractInfeasibilitylIssue ub::T end +ModelAnalyzer.variable(issue::InfeasibleBounds) = issue.variable + +ModelAnalyzer.values(issue::InfeasibleBounds) = [issue.lb, issue.ub] + """ InfeasibleIntegrality{T} <: AbstractInfeasibilitylIssue @@ -77,6 +81,12 @@ struct InfeasibleIntegrality{T} <: AbstractInfeasibilitylIssue set::Union{MOI.Integer,MOI.ZeroOne}#, MOI.Semicontinuous{T}, MOI.Semiinteger{T}} end +ModelAnalyzer.variable(issue::InfeasibleIntegrality) = issue.variable + +ModelAnalyzer.values(issue::InfeasibleIntegrality) = [issue.lb, issue.ub] + +ModelAnalyzer.set(issue::InfeasibleIntegrality) = issue.set + """ InfeasibleConstraintRange{T} <: AbstractInfeasibilitylIssue @@ -101,6 +111,12 @@ struct InfeasibleConstraintRange{T} <: AbstractInfeasibilitylIssue set::Union{MOI.EqualTo{T},MOI.LessThan{T},MOI.GreaterThan{T}} end +ModelAnalyzer.constraint(issue::InfeasibleConstraintRange) = issue.constraint + +ModelAnalyzer.values(issue::InfeasibleConstraintRange) = [issue.lb, issue.ub] + +ModelAnalyzer.set(issue::InfeasibleConstraintRange) = issue.set + """ IrreducibleInfeasibleSubset <: AbstractInfeasibilitylIssue @@ -121,6 +137,8 @@ struct IrreducibleInfeasibleSubset <: AbstractInfeasibilitylIssue constraint::Vector{<:MOI.ConstraintIndex} end +ModelAnalyzer.constraints(issue::IrreducibleInfeasibleSubset) = issue.constraint + """ Data <: ModelAnalyzer.AbstractData diff --git a/test/infeasibility.jl b/test/infeasibility.jl index 176c795..ec96019 100644 --- a/test/infeasibility.jl +++ b/test/infeasibility.jl @@ -37,6 +37,8 @@ function test_bounds() 2.0, 1.0, ) + @test ModelAnalyzer.variable(ret[], model) == y + @test ModelAnalyzer.values(ret[]) == [2.0, 1.0] # buf = IOBuffer() ModelAnalyzer.summarize( @@ -82,6 +84,9 @@ function test_integrality() 2.9, MOI.Integer(), ) + @test ModelAnalyzer.variable(ret[], model) == y + @test ModelAnalyzer.values(ret[]) == [2.2, 2.9] + @test ModelAnalyzer.set(ret[]) == MOI.Integer() # buf = IOBuffer() ModelAnalyzer.summarize( @@ -129,6 +134,9 @@ function test_binary() 0.8, MOI.ZeroOne(), ) + @test ModelAnalyzer.variable(ret[], model) == x + @test ModelAnalyzer.values(ret[]) == [0.5, 0.8] + @test ModelAnalyzer.set(ret[]) == MOI.ZeroOne() return end @@ -150,6 +158,9 @@ function test_range() 22.0, MOI.LessThan{Float64}(1.0), ) + @test ModelAnalyzer.constraint(ret[], model) == c + @test ModelAnalyzer.values(ret[]) == [11.0, 22.0] + @test ModelAnalyzer.set(ret[]) == MOI.LessThan{Float64}(1.0) # buf = IOBuffer() ModelAnalyzer.summarize( @@ -198,6 +209,9 @@ function test_range_neg() 22.0, MOI.LessThan{Float64}(1.0), ) + @test ModelAnalyzer.constraint(ret[], model) == c + @test ModelAnalyzer.values(ret[]) == [11.0, 22.0] + @test ModelAnalyzer.set(ret[]) == MOI.LessThan{Float64}(1.0) return end @@ -219,6 +233,9 @@ function test_range_equalto() 3.0, MOI.EqualTo{Float64}(1.0), ) + @test ModelAnalyzer.constraint(ret[], model) == c + @test ModelAnalyzer.values(ret[]) == [3.0, 3.0] + @test ModelAnalyzer.set(ret[]) == MOI.EqualTo{Float64}(1.0) return end @@ -240,6 +257,9 @@ function test_range_equalto_2() 7.0, MOI.EqualTo{Float64}(1.0), ) + @test ModelAnalyzer.constraint(ret[], model) == c + @test ModelAnalyzer.values(ret[]) == [7.0, 7.0] + @test ModelAnalyzer.set(ret[]) == MOI.EqualTo{Float64}(1.0) return end @@ -261,6 +281,9 @@ function test_range_greaterthan() 22.0, MOI.GreaterThan{Float64}(100.0), ) + @test ModelAnalyzer.constraint(ret[], model) == c + @test ModelAnalyzer.values(ret[]) == [11.0, 22.0] + @test ModelAnalyzer.set(ret[]) == MOI.GreaterThan{Float64}(100.0) return end @@ -282,6 +305,9 @@ function test_range_equalto_3() 22.0, MOI.EqualTo{Float64}(100.0), ) + @test ModelAnalyzer.constraint(ret[], model) == c + @test ModelAnalyzer.values(ret[]) == [11.0, 22.0] + @test ModelAnalyzer.set(ret[]) == MOI.EqualTo{Float64}(100.0) return end @@ -326,6 +352,9 @@ function test_iis() @test length(ret[].constraint) == 2 @test Set([ret[].constraint[1], ret[].constraint[2]]) == Set(JuMP.index.([c2, c1])) + iis = ModelAnalyzer.constraints(ret[], model) + @test length(iis) == 2 + @test Set(iis) == Set([c2, c1]) # buf = IOBuffer() ModelAnalyzer.summarize( @@ -388,6 +417,10 @@ function test_iis_multiple() @test JuMP.index(c2) in Set([ret[].constraint[1], ret[].constraint[2]]) @test Set([ret[].constraint[1], ret[].constraint[2]]) ⊆ Set(JuMP.index.([c3, c2, c1])) + iis = ModelAnalyzer.constraints(ret[], model) + @test length(iis) == 2 + @test Set(iis) ⊆ Set([c3, c2, c1]) + @test c2 in iis return end @@ -415,6 +448,9 @@ function test_iis_interval_right() @test length(ret[].constraint) == 2 @test Set([ret[].constraint[1], ret[].constraint[2]]) == Set(JuMP.index.([c2, c1])) + iis = ModelAnalyzer.constraints(ret[], model) + @test length(iis) == 2 + @test Set(iis) == Set([c2, c1]) return end @@ -442,6 +478,9 @@ function test_iis_interval_left() @test length(ret[].constraint) == 2 @test Set([ret[].constraint[1], ret[].constraint[2]]) == Set(JuMP.index.([c2, c1])) + iis = ModelAnalyzer.constraints(ret[], model) + @test length(iis) == 2 + @test Set(iis) == Set([c2, c1]) return end @@ -472,6 +511,9 @@ function test_iis_spare() @test length(ret[].constraint) == 2 @test Set([ret[].constraint[1], ret[].constraint[2]]) == Set(JuMP.index.([c2, c1])) + iis = ModelAnalyzer.constraints(ret[], model) + @test length(iis) == 2 + @test Set(iis) == Set([c2, c1]) return end From a88f0e8a2bcc7e19892cd1cf64d6b516f695a226 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Wed, 14 May 2025 02:37:50 -0300 Subject: [PATCH 18/34] add tests --- src/_eval_variables.jl | 2 -- test/numerical.jl | 74 ++++++++++++++++++++++++++++++++++-------- 2 files changed, 61 insertions(+), 15 deletions(-) diff --git a/src/_eval_variables.jl b/src/_eval_variables.jl index d30cd88..de0bfab 100644 --- a/src/_eval_variables.jl +++ b/src/_eval_variables.jl @@ -9,8 +9,6 @@ function _eval_variables(value_fn::Function, t::MOI.ScalarAffineTerm) return t.coefficient * value_fn(t.variable) end -_eval_variables(value_fn::Function, f::MOI.VariableIndex) = value_fn(f) - function _eval_variables( value_fn::Function, f::MOI.ScalarAffineFunction{T}, diff --git a/test/numerical.jl b/test/numerical.jl index 62c1e85..26540cb 100644 --- a/test/numerical.jl +++ b/test/numerical.jl @@ -192,42 +192,84 @@ function test_constraint_bounds_quad() return end -function test_variable_not_in_constraints() +function test_no_names() model = Model() + set_string_names_on_creation(model, false) @variable(model, x) @variable(model, y) @constraint(model, 7y >= 3) + @constraint(model, z, 0.0 * y == 3) data = ModelAnalyzer.analyze(ModelAnalyzer.Numerical.Analyzer(), model) list = ModelAnalyzer.list_of_issue_types(data) - @test length(list) == 1 + @test length(list) == 2 ret = ModelAnalyzer.list_of_issues( data, ModelAnalyzer.Numerical.VariableNotInConstraints, ) @test length(ret) == 1 - @test ModelAnalyzer.variable(ret[]) == JuMP.index(x) - @test ModelAnalyzer.variable(ret[], model) == x # buf = IOBuffer() - ModelAnalyzer.summarize( - buf, - ModelAnalyzer.Numerical.VariableNotInConstraints, + ModelAnalyzer.summarize(buf, ret[1], verbose = true) + str = String(take!(buf)) + @test startswith(str, "Variable: ") + ModelAnalyzer.summarize(buf, ret[1], verbose = false) + str = String(take!(buf)) + # + ModelAnalyzer.summarize(buf, ret[1], verbose = true, model = model) + str = String(take!(buf)) + @test startswith(str, "Variable: ") + ModelAnalyzer.summarize(buf, ret[1], verbose = false, model = model) + str = String(take!(buf)) + # + # + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Numerical.EmptyConstraint, ) + @test length(ret) == 1 + # + ModelAnalyzer.summarize(buf, ret[1], verbose = true) str = String(take!(buf)) - @test startswith(str, "# `VariableNotInConstraints`") - ModelAnalyzer.summarize( - buf, + @test startswith(str, "Constraint: ") + ModelAnalyzer.summarize(buf, ret[1], verbose = false) + str = String(take!(buf)) + # + ModelAnalyzer.summarize(buf, ret[1], verbose = true, model = model) + str = String(take!(buf)) + @test startswith(str, "Constraint: ") + ModelAnalyzer.summarize(buf, ret[1], verbose = false, model = model) + str = String(take!(buf)) + return +end + +function test_variable_not_in_constraints() + model = Model() + @variable(model, x) + @variable(model, y) + @constraint(model, 7y >= 3) + data = ModelAnalyzer.analyze(ModelAnalyzer.Numerical.Analyzer(), model) + list = ModelAnalyzer.list_of_issue_types(data) + @test length(list) == 1 + ret = ModelAnalyzer.list_of_issues( + data, ModelAnalyzer.Numerical.VariableNotInConstraints, - verbose = false, ) - str = String(take!(buf)) - @test str == "# VariableNotInConstraints" + @test length(ret) == 1 + @test ModelAnalyzer.variable(ret[]) == JuMP.index(x) + @test ModelAnalyzer.variable(ret[], model) == x # + buf = IOBuffer() ModelAnalyzer.summarize(buf, ret[1], verbose = true) str = String(take!(buf)) @test startswith(str, "Variable: ") ModelAnalyzer.summarize(buf, ret[1], verbose = false) str = String(take!(buf)) + # + ModelAnalyzer.summarize(buf, ret[1], verbose = true, model = model) + str = String(take!(buf)) + @test startswith(str, "Variable: ") + ModelAnalyzer.summarize(buf, ret[1], verbose = false, model = model) + str = String(take!(buf)) return end @@ -264,6 +306,12 @@ function test_empty_constraint_model() @test startswith(str, "Constraint: ") ModelAnalyzer.summarize(buf, ret[1], verbose = false) str = String(take!(buf)) + # + ModelAnalyzer.summarize(buf, ret[1], verbose = true, model = model) + str = String(take!(buf)) + @test startswith(str, "Constraint: ") + ModelAnalyzer.summarize(buf, ret[1], verbose = false, model = model) + str = String(take!(buf)) return end From 1cd08a357315015f7d4d4cbe6a8f65c59ba7b179 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Thu, 22 May 2025 00:16:55 -0700 Subject: [PATCH 19/34] add tests --- test/numerical.jl | 89 ++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 88 insertions(+), 1 deletion(-) diff --git a/test/numerical.jl b/test/numerical.jl index 26540cb..78c070f 100644 --- a/test/numerical.jl +++ b/test/numerical.jl @@ -988,7 +988,94 @@ function test_empty_model() return end -# TODO, test SDP and empty model +function test_vector_functions() + model = Model() + @variable(model, x[1:3]) + @constraint(model, c1, [1e-9 * x[1]] in Nonnegatives()) + @constraint(model, c2, [1e+9 * x[1]] in Nonnegatives()) + @constraint(model, c3, [x[2], x[1]] in Nonnegatives()) + @constraint(model, c4, [-1e-9 * x[1] * x[1]] in Nonnegatives()) + @constraint(model, c5, [1e+9 * x[1] * x[1]] in Nonnegatives()) + data = ModelAnalyzer.analyze(ModelAnalyzer.Numerical.Analyzer(), model) + @show list = ModelAnalyzer.list_of_issue_types(data) + @test length(list) == 6 + # + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Numerical.SmallMatrixCoefficient, + ) + @test length(ret) == 1 + @test ModelAnalyzer.constraint(ret[], model) == c1 + @test ModelAnalyzer.variable(ret[], model) == x[1] + @test ModelAnalyzer.value(ret[]) == 1e-9 + # + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Numerical.LargeMatrixCoefficient, + ) + @test length(ret) == 1 + @test ModelAnalyzer.constraint(ret[], model) == c2 + @test ModelAnalyzer.variable(ret[], model) == x[1] + @test ModelAnalyzer.value(ret[]) == 1e+9 + # + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Numerical.SmallMatrixQuadraticCoefficient, + ) + @test length(ret) == 1 + @test ModelAnalyzer.constraint(ret[], model) == c4 + @test ModelAnalyzer.variables(ret[], model) == [x[1], x[1]] + @test_broken ModelAnalyzer.value(ret[]) == 1e-9 + # + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Numerical.LargeMatrixQuadraticCoefficient, + ) + @test length(ret) == 1 + @test ModelAnalyzer.constraint(ret[], model) == c5 + @test ModelAnalyzer.variables(ret[], model) == [x[1], x[1]] + @test_broken ModelAnalyzer.value(ret[]) == 1e+9 + # + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Numerical.VariableNotInConstraints, + ) + @test length(ret) == 1 + @test ModelAnalyzer.variable(ret[], model) == x[3] + return +end + +function test_variable_interval() + model = Model() + @variable(model, x in MOI.Interval(1e-9, 1e+9)) + @objective(model, Min, x) + data = ModelAnalyzer.analyze(ModelAnalyzer.Numerical.Analyzer(), model) + @show list = ModelAnalyzer.list_of_issue_types(data) + @test length(list) == 3 + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Numerical.SmallBoundCoefficient, + ) + @test length(ret) == 1 + @test ModelAnalyzer.variable(ret[], model) == x + @test ModelAnalyzer.value(ret[]) == 1e-9 + # + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Numerical.LargeBoundCoefficient, + ) + @test length(ret) == 1 + @test ModelAnalyzer.variable(ret[], model) == x + @test ModelAnalyzer.value(ret[]) == 1e+9 + # + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Numerical.VariableNotInConstraints, + ) + @test length(ret) == 1 + @test ModelAnalyzer.variable(ret[], model) == x + return +end function test_many() model = Model() From b5253c79c872625f546886b6d6300b82599578ac Mon Sep 17 00:00:00 2001 From: joaquimg Date: Thu, 22 May 2025 00:54:42 -0700 Subject: [PATCH 20/34] add tests --- test/numerical.jl | 52 +++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 50 insertions(+), 2 deletions(-) diff --git a/test/numerical.jl b/test/numerical.jl index 78c070f..fd49aca 100644 --- a/test/numerical.jl +++ b/test/numerical.jl @@ -988,6 +988,53 @@ function test_empty_model() return end +function test_nonconvex_zeros() + model = Model() + @variable(model, x[1:1]) + @constraint(model, c, [x[1] * x[1]] in Zeros()) + data = ModelAnalyzer.analyze(ModelAnalyzer.Numerical.Analyzer(), model) + list = ModelAnalyzer.list_of_issue_types(data) + @test length(list) == 1 + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Numerical.NonconvexQuadraticConstraint, + ) + @test length(ret) == 1 + @test ModelAnalyzer.constraint(ret[], model) == c + return +end + +function test_vi_in_nonstandard_set() + model = Model() + @variable(model, x[1:1]) + @constraint(model, c, x[1] in MOI.ZeroOne()) + @constraint(model, c2, 2x[1] == 0) + data = ModelAnalyzer.analyze(ModelAnalyzer.Numerical.Analyzer(), model) + list = ModelAnalyzer.list_of_issue_types(data) + @test length(list) == 0 + return +end + +function test_saf_in_nonstandard_set() + model = Model() + @variable(model, x[1:1]) + @constraint(model, c, 2x[1] in MOI.ZeroOne()) + data = ModelAnalyzer.analyze(ModelAnalyzer.Numerical.Analyzer(), model) + list = ModelAnalyzer.list_of_issue_types(data) + @test length(list) == 0 + return +end + +function test_vaf_in_nonstandard_set() + model = Model() + @variable(model, x[1:1]) + @constraint(model, c, [2x[1], x[1], x[1]] in SecondOrderCone()) + data = ModelAnalyzer.analyze(ModelAnalyzer.Numerical.Analyzer(), model) + list = ModelAnalyzer.list_of_issue_types(data) + @test length(list) == 0 + return +end + function test_vector_functions() model = Model() @variable(model, x[1:3]) @@ -996,8 +1043,9 @@ function test_vector_functions() @constraint(model, c3, [x[2], x[1]] in Nonnegatives()) @constraint(model, c4, [-1e-9 * x[1] * x[1]] in Nonnegatives()) @constraint(model, c5, [1e+9 * x[1] * x[1]] in Nonnegatives()) + @constraint(model, c6, [2 * x[1] * x[1]] in Zeros()) data = ModelAnalyzer.analyze(ModelAnalyzer.Numerical.Analyzer(), model) - @show list = ModelAnalyzer.list_of_issue_types(data) + list = ModelAnalyzer.list_of_issue_types(data) @test length(list) == 6 # ret = ModelAnalyzer.list_of_issues( @@ -1050,7 +1098,7 @@ function test_variable_interval() @variable(model, x in MOI.Interval(1e-9, 1e+9)) @objective(model, Min, x) data = ModelAnalyzer.analyze(ModelAnalyzer.Numerical.Analyzer(), model) - @show list = ModelAnalyzer.list_of_issue_types(data) + list = ModelAnalyzer.list_of_issue_types(data) @test length(list) == 3 ret = ModelAnalyzer.list_of_issues( data, From 4fc7565568f8ed9638794b1039c56805c1cb7ae2 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Thu, 22 May 2025 00:54:56 -0700 Subject: [PATCH 21/34] add comments --- src/feasibility.jl | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/src/feasibility.jl b/src/feasibility.jl index 2a0a952..b36a736 100644 --- a/src/feasibility.jl +++ b/src/feasibility.jl @@ -97,7 +97,12 @@ ModelAnalyzer.value(issue::DualConstraintViolation) = issue.violation The `DualConstrainedVariableViolation` issue is identified when a dual constraint, which is a constrained varaible constraint, has a value that is not within the dual constraint's set. -This dual constraint corresponds to a primal constraint. +During the dualization process, each primal constraint is mapped to a dual +variable, this dual variable is tipically a constrained variable with the +dual set of the primal constraint. If the primal constraint is a an equality +type constraint, the dual variable is a free variable, hence, not constrained +(dual) variable. +This dual constraint corresponds to a primal (non-equality) constraint. For more information, run: ```julia @@ -371,6 +376,11 @@ function ModelAnalyzer._verbose_summarize( A `DualConstrainedVariableViolation` issue is identified when a dual constraint, which is a constrained varaible constraint, has a value that is not within the dual constraint's set. + During the dualization process, each primal constraint is mapped to a dual + variable, this dual variable is tipically a constrained variable with the + dual set of the primal constraint. If the primal constraint is a an equality + type constraint, the dual variable is a free variable, hence, not constrained + (dual) variable. ## Why From 5ecd302c3bdedd27210abeec9934ffcfd6806887 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Thu, 22 May 2025 00:58:52 -0700 Subject: [PATCH 22/34] simplify api --- src/ModelAnalyzer.jl | 28 ---------------------------- 1 file changed, 28 deletions(-) diff --git a/src/ModelAnalyzer.jl b/src/ModelAnalyzer.jl index a7302c1..1197206 100644 --- a/src/ModelAnalyzer.jl +++ b/src/ModelAnalyzer.jl @@ -143,10 +143,6 @@ Return the value associated to a particular issue. The value is a number with a different meaning depending on the type of issue. For example, for some numerical issues, it can be the coefficient value. """ -function value(issue::AbstractIssue, ::Nothing) - return value(issue) -end - function value(issue::AbstractIssue, ::MOI.ModelLike) return value(issue) end @@ -156,10 +152,6 @@ end Return the values associated to a particular issue. """ -function values(issue::AbstractIssue, ::Nothing) - return values(issue) -end - function values(issue::AbstractIssue, ::MOI.ModelLike) return values(issue) end @@ -169,10 +161,6 @@ end Return the variable associated to a particular issue. """ -function variable(issue::AbstractIssue, ::Nothing) - return variable(issue) -end - function variable(issue::AbstractIssue, ::MOI.ModelLike) return variable(issue) end @@ -182,10 +170,6 @@ end Return the variables associated to a particular issue. """ -function variables(issue::AbstractIssue, ::Nothing) - return variables(issue) -end - function variables(issue::AbstractIssue, ::MOI.ModelLike) return variables(issue) end @@ -195,10 +179,6 @@ end Return the constraint associated to a particular issue. """ -function constraint(issue::AbstractIssue, ::Nothing) - return constraint(issue) -end - function constraint(issue::AbstractIssue, ::MOI.ModelLike) return constraint(issue) end @@ -208,10 +188,6 @@ end Return the constraints associated to a particular issue. """ -function constraints(issue::AbstractIssue, ::Nothing) - return constraints(issue) -end - function constraints(issue::AbstractIssue, ::MOI.ModelLike) return constraints(issue) end @@ -221,10 +197,6 @@ end Return the set associated to a particular issue. """ -function set(issue::AbstractIssue, ::Nothing) - return set(issue) -end - function set(issue::AbstractIssue, ::MOI.ModelLike) return set(issue) end From ceccb9ff9702256a6dcfd1ee858d5fb9eecc6c8c Mon Sep 17 00:00:00 2001 From: joaquimg Date: Thu, 22 May 2025 01:09:07 -0700 Subject: [PATCH 23/34] cleanup api --- src/ModelAnalyzer.jl | 12 +++--------- test/infeasibility.jl | 3 +++ test/numerical.jl | 4 ++++ 3 files changed, 10 insertions(+), 9 deletions(-) diff --git a/src/ModelAnalyzer.jl b/src/ModelAnalyzer.jl index 1197206..c913b18 100644 --- a/src/ModelAnalyzer.jl +++ b/src/ModelAnalyzer.jl @@ -143,18 +143,14 @@ Return the value associated to a particular issue. The value is a number with a different meaning depending on the type of issue. For example, for some numerical issues, it can be the coefficient value. """ -function value(issue::AbstractIssue, ::MOI.ModelLike) - return value(issue) -end +function value end """ values(issue::AbstractIssue) Return the values associated to a particular issue. """ -function values(issue::AbstractIssue, ::MOI.ModelLike) - return values(issue) -end +function values end """ variable(issue::AbstractIssue) @@ -197,9 +193,7 @@ end Return the set associated to a particular issue. """ -function set(issue::AbstractIssue, ::MOI.ModelLike) - return set(issue) -end +function set end function _verbose_summarize end diff --git a/test/infeasibility.jl b/test/infeasibility.jl index ec96019..7af0850 100644 --- a/test/infeasibility.jl +++ b/test/infeasibility.jl @@ -481,6 +481,9 @@ function test_iis_interval_left() iis = ModelAnalyzer.constraints(ret[], model) @test length(iis) == 2 @test Set(iis) == Set([c2, c1]) + iis = ModelAnalyzer.constraints(ret[], JuMP.backend(model)) + @test length(iis) == 2 + @test Set(iis) == Set(JuMP.index.([c2, c1])) return end diff --git a/test/numerical.jl b/test/numerical.jl index fd49aca..705b4f4 100644 --- a/test/numerical.jl +++ b/test/numerical.jl @@ -367,6 +367,7 @@ function test_dense_constraint() ) @test length(ret) == 1 @test ModelAnalyzer.constraint(ret[]) == JuMP.index(c) + @test ModelAnalyzer.constraint(ret[], JuMP.backend(model)) == JuMP.index(c) @test ModelAnalyzer.constraint(ret[], model) == c @test ModelAnalyzer.value(ret[]) == 10_000 # @@ -405,6 +406,7 @@ function test_small_matrix_coef() @test length(ret) == 1 @test ModelAnalyzer.constraint(ret[], model) == c @test ModelAnalyzer.variable(ret[], model) == x + @test ModelAnalyzer.variable(ret[], JuMP.backend(model)) == JuMP.index(x) @test ModelAnalyzer.value(ret[]) == 1e-9 # buf = IOBuffer() @@ -713,6 +715,8 @@ function test_small_objective_coef_quad() ) @test length(ret) == 1 @test ModelAnalyzer.variables(ret[], model) == [x, x] + @test ModelAnalyzer.variables(ret[], JuMP.backend(model)) == + JuMP.index.([x, x]) @test_broken ModelAnalyzer.value(ret[]) == 1e-9 # 2e-9 TODO, what to return here # buf = IOBuffer() From 97cde568ed7bece0cd659e08632fdc16185a8872 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Thu, 22 May 2025 01:33:13 -0700 Subject: [PATCH 24/34] add tests --- test/feasibility.jl | 47 ++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 46 insertions(+), 1 deletion(-) diff --git a/test/feasibility.jl b/test/feasibility.jl index 797527e..fcf1a10 100644 --- a/test/feasibility.jl +++ b/test/feasibility.jl @@ -109,7 +109,7 @@ function test_only_bounds() ModelAnalyzer.Feasibility.PrimalDualMismatch, ) @test ret[] == ModelAnalyzer.Feasibility.PrimalDualMismatch(-1.0, 0.0) - # @test ModelAnalyzer.values(ret[]) == [-1.0, 0.0] + @test ModelAnalyzer.values(ret[]) == [-1.0, 0.0] data = ModelAnalyzer.analyze( ModelAnalyzer.Feasibility.Analyzer(), @@ -718,6 +718,51 @@ function test_analyse_no_opt() return end +function test_dual_constrained_variable() + model = Model() + @variable(model, x >= 0) + @objective(model, Min, x) + + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(JuMP.index(x) => 1.0), + dual_point = Dict(JuMP.index(LowerBoundRef(x)) => -1.0), + ) + list = ModelAnalyzer.list_of_issue_types(data) + @test list == [ + ModelAnalyzer.Feasibility.DualConstraintViolation, + ModelAnalyzer.Feasibility.DualConstrainedVariableViolation, + ModelAnalyzer.Feasibility.ComplemetarityViolation, + ModelAnalyzer.Feasibility.PrimalDualMismatch, + ] + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.DualConstrainedVariableViolation, + ) + @test ret[] == ModelAnalyzer.Feasibility.DualConstrainedVariableViolation( + JuMP.index(LowerBoundRef(x)), + 1.0, + ) + @test ModelAnalyzer.constraint(ret[], model) == LowerBoundRef(x) + @test ModelAnalyzer.value(ret[]) == 1.0 + + buf = IOBuffer() + ModelAnalyzer.summarize( + buf, + ModelAnalyzer.Feasibility.DualConstrainedVariableViolation, + verbose = true, + ) + ModelAnalyzer.summarize( + buf, + ModelAnalyzer.Feasibility.DualConstrainedVariableViolation, + verbose = false, + ) + ModelAnalyzer.summarize(buf, ret[], verbose = true) + ModelAnalyzer.summarize(buf, ret[], verbose = false) + return +end + # these tests are harder to permorm with a real solver as they tipically # return coherent objectives function test_lowlevel_mismatch() From 0b3e9e482874d6c782c02e405febe34e92d5446c Mon Sep 17 00:00:00 2001 From: joaquimg Date: Thu, 22 May 2025 01:38:48 -0700 Subject: [PATCH 25/34] add tests --- test/numerical.jl | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/test/numerical.jl b/test/numerical.jl index 705b4f4..ef80e29 100644 --- a/test/numerical.jl +++ b/test/numerical.jl @@ -218,6 +218,9 @@ function test_no_names() ModelAnalyzer.summarize(buf, ret[1], verbose = true, model = model) str = String(take!(buf)) @test startswith(str, "Variable: ") + ModelAnalyzer.summarize(buf, ret[1], verbose = true, model = JuMP.backend(model)) + str = String(take!(buf)) + @test startswith(str, "Variable: ") ModelAnalyzer.summarize(buf, ret[1], verbose = false, model = model) str = String(take!(buf)) # @@ -237,6 +240,9 @@ function test_no_names() ModelAnalyzer.summarize(buf, ret[1], verbose = true, model = model) str = String(take!(buf)) @test startswith(str, "Constraint: ") + ModelAnalyzer.summarize(buf, ret[1], verbose = true, model = JuMP.backend(model)) + str = String(take!(buf)) + @test startswith(str, "Constraint: ") ModelAnalyzer.summarize(buf, ret[1], verbose = false, model = model) str = String(take!(buf)) return @@ -268,6 +274,9 @@ function test_variable_not_in_constraints() ModelAnalyzer.summarize(buf, ret[1], verbose = true, model = model) str = String(take!(buf)) @test startswith(str, "Variable: ") + ModelAnalyzer.summarize(buf, ret[1], verbose = true, model = JuMP.backend(model)) + str = String(take!(buf)) + @test startswith(str, "Variable: ") ModelAnalyzer.summarize(buf, ret[1], verbose = false, model = model) str = String(take!(buf)) return @@ -310,6 +319,9 @@ function test_empty_constraint_model() ModelAnalyzer.summarize(buf, ret[1], verbose = true, model = model) str = String(take!(buf)) @test startswith(str, "Constraint: ") + ModelAnalyzer.summarize(buf, ret[1], verbose = true, model = JuMP.backend(model)) + str = String(take!(buf)) + @test startswith(str, "Constraint: ") ModelAnalyzer.summarize(buf, ret[1], verbose = false, model = model) str = String(take!(buf)) return From e1288d6988d507b6e0308361cde1691113a63c77 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Thu, 22 May 2025 01:42:22 -0700 Subject: [PATCH 26/34] format --- test/numerical.jl | 28 ++++++++++++++++++++++++---- 1 file changed, 24 insertions(+), 4 deletions(-) diff --git a/test/numerical.jl b/test/numerical.jl index ef80e29..8a072c6 100644 --- a/test/numerical.jl +++ b/test/numerical.jl @@ -218,7 +218,12 @@ function test_no_names() ModelAnalyzer.summarize(buf, ret[1], verbose = true, model = model) str = String(take!(buf)) @test startswith(str, "Variable: ") - ModelAnalyzer.summarize(buf, ret[1], verbose = true, model = JuMP.backend(model)) + ModelAnalyzer.summarize( + buf, + ret[1], + verbose = true, + model = JuMP.backend(model), + ) str = String(take!(buf)) @test startswith(str, "Variable: ") ModelAnalyzer.summarize(buf, ret[1], verbose = false, model = model) @@ -240,7 +245,12 @@ function test_no_names() ModelAnalyzer.summarize(buf, ret[1], verbose = true, model = model) str = String(take!(buf)) @test startswith(str, "Constraint: ") - ModelAnalyzer.summarize(buf, ret[1], verbose = true, model = JuMP.backend(model)) + ModelAnalyzer.summarize( + buf, + ret[1], + verbose = true, + model = JuMP.backend(model), + ) str = String(take!(buf)) @test startswith(str, "Constraint: ") ModelAnalyzer.summarize(buf, ret[1], verbose = false, model = model) @@ -274,7 +284,12 @@ function test_variable_not_in_constraints() ModelAnalyzer.summarize(buf, ret[1], verbose = true, model = model) str = String(take!(buf)) @test startswith(str, "Variable: ") - ModelAnalyzer.summarize(buf, ret[1], verbose = true, model = JuMP.backend(model)) + ModelAnalyzer.summarize( + buf, + ret[1], + verbose = true, + model = JuMP.backend(model), + ) str = String(take!(buf)) @test startswith(str, "Variable: ") ModelAnalyzer.summarize(buf, ret[1], verbose = false, model = model) @@ -319,7 +334,12 @@ function test_empty_constraint_model() ModelAnalyzer.summarize(buf, ret[1], verbose = true, model = model) str = String(take!(buf)) @test startswith(str, "Constraint: ") - ModelAnalyzer.summarize(buf, ret[1], verbose = true, model = JuMP.backend(model)) + ModelAnalyzer.summarize( + buf, + ret[1], + verbose = true, + model = JuMP.backend(model), + ) str = String(take!(buf)) @test startswith(str, "Constraint: ") ModelAnalyzer.summarize(buf, ret[1], verbose = false, model = model) From e719b9643b9d06b1e1a85a26becb9878adfdff1b Mon Sep 17 00:00:00 2001 From: joaquimg Date: Fri, 23 May 2025 00:18:23 -0700 Subject: [PATCH 27/34] more tests --- src/feasibility.jl | 37 ++++++++++- test/feasibility.jl | 150 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 186 insertions(+), 1 deletion(-) diff --git a/src/feasibility.jl b/src/feasibility.jl index b36a736..97e409c 100644 --- a/src/feasibility.jl +++ b/src/feasibility.jl @@ -983,6 +983,39 @@ function _dual_point_to_dual_model_ref( new_dual_point = Dict{MOI.VariableIndex,Number}() dual_var_to_primal_con = Dict{MOI.VariableIndex,MOI.ConstraintIndex}() dual_con_to_primal_con = Dict{MOI.ConstraintIndex,MOI.ConstraintIndex}() + for (F, S) in MOI.get(primal_model, MOI.ListOfConstraintTypesPresent()) + list = MOI.get(primal_model, MOI.ListOfConstraintIndices{F,S}()) + for primal_con in list + dual_vars = Dualization._get_dual_variables(map, primal_con) + val = get(dual_point, primal_con, nothing) + if !isnothing(val) + if length(dual_vars) != length(val) + error( + "The dual point entry for constraint $primal_con has " * + "length $(length(val)) but the dual variable " * + "length is $(length(dual_vars)).", + ) + end + for (idx, dual_var) in enumerate(dual_vars) + new_dual_point[dual_var] = val[idx] + end + end + for dual_var in dual_vars + dual_var_to_primal_con[dual_var] = primal_con + end + dual_con = Dualization._get_dual_constraint(map, primal_con) + if dual_con !== nothing + dual_con_to_primal_con[dual_con] = primal_con + # else + # if !(primal_con isa MOI.ConstraintIndex{MOI.VariableIndex,<:MOI.EqualTo} || + # primal_con isa MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.Zeros} + # SAF in EQ, etc... + #) + # error("Problem with dualization, see: $primal_con") + # end + end + end + end for (primal_con, val) in dual_point dual_vars = Dualization._get_dual_variables(map, primal_con) if length(dual_vars) != length(val) @@ -1247,7 +1280,9 @@ function _analyze_objectives!(model::MOI.ModelLike, dual_model, map, data) ) end - if !isapprox(obj_val, dual_obj_val; atol = data.atol) + if !isapprox(obj_val, dual_obj_val; atol = data.atol) && + !isnan(dual_obj_val) && + !isnan(obj_val) push!( data.primal_dual_mismatch, PrimalDualMismatch(obj_val, dual_obj_val), diff --git a/test/feasibility.jl b/test/feasibility.jl index fcf1a10..4c98137 100644 --- a/test/feasibility.jl +++ b/test/feasibility.jl @@ -796,6 +796,156 @@ function test_lowlevel_mismatch() return end +function test_skip_missing_primal() + model = Model() + set_silent(model) + @variable(model, x >= 0) + @variable(model, y >= 0) + @constraint(model, c, x + y >= 0) + @objective(model, Min, x) + + @test_throws ErrorException ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(JuMP.index(x) => 1.0), + # dual_point = Dict(JuMP.index(c) => 1.0), + skip_missing = false, + dual_check = false, + ) + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(JuMP.index(x) => 1.0), + # dual_point = Dict(JuMP.index(c) => 1.0), + skip_missing = true, + dual_check = false, + ) + list = ModelAnalyzer.list_of_issue_types(data) + @test length(list) == 0 + return +end + +function test_skip_missing_primal_var_not_in_con() + model = Model() + set_silent(model) + @variable(model, x) + @variable(model, y) + @constraint(model, c, x >= 0) + @objective(model, Min, x + y) + + @test_throws ErrorException ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(JuMP.index(x) => 0.0), + # dual_point = Dict(JuMP.index(c) => 1.0), + skip_missing = false, + dual_check = false, + ) + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(JuMP.index(x) => 0.0), + # dual_point = Dict(JuMP.index(c) => 1.0), + skip_missing = true, + dual_check = false, + ) + list = ModelAnalyzer.list_of_issue_types(data) + @test length(list) == 0 + return +end + +function test_skip_missing_primal_empty_con() + model = Model() + set_silent(model) + @variable(model, x) + @constraint(model, c, 1 >= 0) + @constraint(model, c2, x >= 0) + @objective(model, Min, x) + + @test_throws ErrorException ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(JuMP.index(x) => 0.0), + dual_point = Dict(JuMP.index(c2) => 1.0), + skip_missing = false, + dual_check = true, + ) + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(JuMP.index(x) => 0.0), + dual_point = Dict(JuMP.index(c2) => 1.0), + skip_missing = true, + dual_check = true, + ) + @show list = ModelAnalyzer.list_of_issue_types(data) + @test length(list) == 0 + return +end + +function test_skip_missing_dual() + model = Model() + set_silent(model) + @variable(model, x) + @constraint(model, c1, x >= 0) + @constraint(model, c2, x >= 2) + @objective(model, Min, x) + + @test_throws ErrorException ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(JuMP.index(x) => 2.0), + dual_point = Dict(JuMP.index(c1) => 1.0), + skip_missing = false, + dual_check = true, + ) + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(JuMP.index(x) => 2.0), + dual_point = Dict(JuMP.index(c1) => 0.0), + skip_missing = true, + dual_check = true, + ) + @show list = ModelAnalyzer.list_of_issue_types(data) + @test length(list) == 0 + return +end + +function test_dual_bad_size() + model = Model() + set_silent(model) + @variable(model, x) + @constraint(model, c1, x >= 0) + @objective(model, Min, x) + + @test_throws ErrorException ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(JuMP.index(x) => 2.0), + dual_point = Dict(JuMP.index(c1) => [1.0, 2.0]), + ) + return +end + +function test_dual_vector() + model = Model() + set_silent(model) + @variable(model, x) + @constraint(model, c1, [x, 2x - 1] in Nonnegatives()) + @objective(model, Min, x) + + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(JuMP.index(x) => 0.5), + dual_point = Dict(JuMP.index(c1) => [0.0, 0.5]), + ) + @show list = ModelAnalyzer.list_of_issue_types(data) + @test length(list) == 0 + return +end + end # module TestDualFeasibilityChecker.runtests() From 7e36757f038257275664a9113a11522349dbff8d Mon Sep 17 00:00:00 2001 From: joaquimg Date: Fri, 23 May 2025 00:24:06 -0700 Subject: [PATCH 28/34] more tests --- test/feasibility.jl | 36 ++++++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/test/feasibility.jl b/test/feasibility.jl index 4c98137..bdd6774 100644 --- a/test/feasibility.jl +++ b/test/feasibility.jl @@ -946,6 +946,42 @@ function test_dual_vector() return end +function test_nl_con() + model = Model() + set_silent(model) + @variable(model, x) + @constraint(model, c1, x^3 == 0) + @objective(model, Min, x) + + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(JuMP.index(x) => 0.0), + dual_point = Dict(JuMP.index(c1) => 0.0), + ) + @show list = ModelAnalyzer.list_of_issue_types(data) + @test length(list) == 0 + return +end + +function test_nl_obj() + model = Model() + set_silent(model) + @variable(model, x) + @constraint(model, c1, x == 0) + @objective(model, Min, x^3) + + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(JuMP.index(x) => 0.0), + dual_point = Dict(JuMP.index(c1) => 0.0), + ) + @show list = ModelAnalyzer.list_of_issue_types(data) + @test length(list) == 0 + return +end + end # module TestDualFeasibilityChecker.runtests() From 90903ee59e3769414e60759e0d78ed1e2adb7375 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Fri, 23 May 2025 00:39:59 -0700 Subject: [PATCH 29/34] fix and test --- src/feasibility.jl | 25 ------------------------- test/feasibility.jl | 10 ++++++++++ 2 files changed, 10 insertions(+), 25 deletions(-) diff --git a/src/feasibility.jl b/src/feasibility.jl index 97e409c..01c7222 100644 --- a/src/feasibility.jl +++ b/src/feasibility.jl @@ -1016,31 +1016,6 @@ function _dual_point_to_dual_model_ref( end end end - for (primal_con, val) in dual_point - dual_vars = Dualization._get_dual_variables(map, primal_con) - if length(dual_vars) != length(val) - error( - "The dual point entry for constraint $primal_con has " * - "length $(length(val)) but the dual variable " * - "length is $(length(dual_vars)).", - ) - end - for (idx, dual_var) in enumerate(dual_vars) - new_dual_point[dual_var] = val[idx] - dual_var_to_primal_con[dual_var] = primal_con - end - dual_con = Dualization._get_dual_constraint(map, primal_con) - if dual_con !== nothing - dual_con_to_primal_con[dual_con] = primal_con - # else - # if !(primal_con isa MOI.ConstraintIndex{MOI.VariableIndex,<:MOI.EqualTo} || - # primal_con isa MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.Zeros} - # SAF in EQ, etc... - #) - # error("Problem with dualization, see: $primal_con") - # end - end - end primal_vars = MOI.get(primal_model, MOI.ListOfVariableIndices()) dual_con_to_primal_vars = Dict{MOI.ConstraintIndex,Vector{MOI.VariableIndex}}() diff --git a/test/feasibility.jl b/test/feasibility.jl index bdd6774..1cfb60f 100644 --- a/test/feasibility.jl +++ b/test/feasibility.jl @@ -822,6 +822,16 @@ function test_skip_missing_primal() ) list = ModelAnalyzer.list_of_issue_types(data) @test length(list) == 0 + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(JuMP.index(y) => 1.0), + dual_point = Dict(JuMP.index(c) => 1.0), + skip_missing = true, + dual_check = true, + ) + @show list = ModelAnalyzer.list_of_issue_types(data) + @test length(list) == 0 return end From 1426502b374781a824a10cd036d85820de4a1d39 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Fri, 23 May 2025 00:54:26 -0700 Subject: [PATCH 30/34] add tests and cleanup --- src/infeasibility.jl | 56 ++++++++++++++++++++++++------------------- test/feasibility.jl | 12 +++++----- test/infeasibility.jl | 17 +++++++++++++ 3 files changed, 54 insertions(+), 31 deletions(-) diff --git a/src/infeasibility.jl b/src/infeasibility.jl index 7bee7f1..21a2852 100644 --- a/src/infeasibility.jl +++ b/src/infeasibility.jl @@ -293,15 +293,16 @@ function ModelAnalyzer.analyze( func = MOI.get(model, MOI.ConstraintFunction(), con) failed = false interval = _eval_variables(func) do var_idx - if !haskey(variables, var_idx) - failed = true - return Interval(-Inf, Inf) - end + # 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 + # if failed + # continue + # end rhs = set.value if interval.lo > rhs || interval.hi < rhs push!( @@ -323,15 +324,16 @@ function ModelAnalyzer.analyze( func = MOI.get(model, MOI.ConstraintFunction(), con) failed = false interval = _eval_variables(func) do var_idx - if !haskey(variables, var_idx) - failed = true - return Interval(-Inf, Inf) - end + # 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 + # if failed + # continue + # end rhs = set.upper if interval.lo > rhs push!( @@ -353,15 +355,16 @@ function ModelAnalyzer.analyze( func = MOI.get(model, MOI.ConstraintFunction(), con) failed = false interval = _eval_variables(func) do var_idx - if !haskey(variables, var_idx) - failed = true - return Interval(-Inf, Inf) - end + # 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 + # if failed + # continue + # end rhs = set.lower if interval.hi < rhs push!( @@ -401,8 +404,9 @@ function _fix_to_zero(model, variable::MOI.VariableIndex, ::Type{T}) where {T} if MOI.is_valid(model, lb_idx) MOI.delete(model, lb_idx) has_lower = true - elseif MOI.is_valid(model, ub_idx) - MOI.delete(model, ub_idx) + # 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 @@ -422,8 +426,9 @@ function _set_bound_zero( MOI.delete(model, eq_idx) if has_lower MOI.add_constraint(model, variable, MOI.GreaterThan{T}(zero(T))) - else - MOI.add_constraint(model, variable, MOI.LessThan{T}(zero(T))) + # MOI.PenaltyRelaxation only creates variables with LB + # else + # MOI.add_constraint(model, variable, MOI.LessThan{T}(zero(T))) end return end @@ -479,6 +484,7 @@ function iis_elastic_filter(original_model::MOI.ModelLike, optimizer) 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 diff --git a/test/feasibility.jl b/test/feasibility.jl index 1cfb60f..0d65111 100644 --- a/test/feasibility.jl +++ b/test/feasibility.jl @@ -830,7 +830,7 @@ function test_skip_missing_primal() skip_missing = true, dual_check = true, ) - @show list = ModelAnalyzer.list_of_issue_types(data) + list = ModelAnalyzer.list_of_issue_types(data) @test length(list) == 0 return end @@ -888,7 +888,7 @@ function test_skip_missing_primal_empty_con() skip_missing = true, dual_check = true, ) - @show list = ModelAnalyzer.list_of_issue_types(data) + list = ModelAnalyzer.list_of_issue_types(data) @test length(list) == 0 return end @@ -917,7 +917,7 @@ function test_skip_missing_dual() skip_missing = true, dual_check = true, ) - @show list = ModelAnalyzer.list_of_issue_types(data) + list = ModelAnalyzer.list_of_issue_types(data) @test length(list) == 0 return end @@ -951,7 +951,7 @@ function test_dual_vector() primal_point = Dict(JuMP.index(x) => 0.5), dual_point = Dict(JuMP.index(c1) => [0.0, 0.5]), ) - @show list = ModelAnalyzer.list_of_issue_types(data) + list = ModelAnalyzer.list_of_issue_types(data) @test length(list) == 0 return end @@ -969,7 +969,7 @@ function test_nl_con() primal_point = Dict(JuMP.index(x) => 0.0), dual_point = Dict(JuMP.index(c1) => 0.0), ) - @show list = ModelAnalyzer.list_of_issue_types(data) + list = ModelAnalyzer.list_of_issue_types(data) @test length(list) == 0 return end @@ -987,7 +987,7 @@ function test_nl_obj() primal_point = Dict(JuMP.index(x) => 0.0), dual_point = Dict(JuMP.index(c1) => 0.0), ) - @show list = ModelAnalyzer.list_of_issue_types(data) + list = ModelAnalyzer.list_of_issue_types(data) @test length(list) == 0 return end diff --git a/test/infeasibility.jl b/test/infeasibility.jl index 7af0850..cde2582 100644 --- a/test/infeasibility.jl +++ b/test/infeasibility.jl @@ -311,6 +311,23 @@ function test_range_equalto_3() return end +function test_interval() + model = Model(HiGHS.Optimizer) + set_silent(model) + @variable(model, x in MOI.Interval(0, 10)) + @variable(model, 0 <= y <= 20) + @constraint(model, c1, x + y <= 1) + @objective(model, Max, x + y) + optimize!(model) + data = ModelAnalyzer.analyze( + ModelAnalyzer.Infeasibility.Analyzer(), + model, + optimizer = HiGHS.Optimizer, + ) + list = ModelAnalyzer.list_of_issue_types(data) + @test length(list) == 0 +end + function test_iis_feasible() model = Model(HiGHS.Optimizer) set_silent(model) From 6f641ad45af5b9564654c53429606c2a189dcb29 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Fri, 23 May 2025 01:27:05 -0700 Subject: [PATCH 31/34] format --- src/infeasibility.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/infeasibility.jl b/src/infeasibility.jl index 21a2852..2dab25d 100644 --- a/src/infeasibility.jl +++ b/src/infeasibility.jl @@ -404,9 +404,9 @@ function _fix_to_zero(model, variable::MOI.VariableIndex, ::Type{T}) where {T} 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) + # 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 @@ -426,9 +426,9 @@ function _set_bound_zero( 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))) + # MOI.PenaltyRelaxation only creates variables with LB + # else + # MOI.add_constraint(model, variable, MOI.LessThan{T}(zero(T))) end return end From e9035d73bffc0b45f8320d13f3c98e449759478d Mon Sep 17 00:00:00 2001 From: joaquimg Date: Sat, 24 May 2025 01:03:24 -0700 Subject: [PATCH 32/34] add features and tests --- src/feasibility.jl | 49 ++++++++++++++++++++---------- test/feasibility.jl | 74 +++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 107 insertions(+), 16 deletions(-) diff --git a/src/feasibility.jl b/src/feasibility.jl index 01c7222..654f22f 100644 --- a/src/feasibility.jl +++ b/src/feasibility.jl @@ -23,6 +23,8 @@ julia> data = ModelAnalyzer.analyze( model; primal_point::Union{Nothing, Dict} = nothing, dual_point::Union{Nothing, Dict} = nothing, + primal_objective::Union{Nothing, Float64} = nothing, + dual_objective::Union{Nothing, Float64} = nothing, atol::Float64 = 1e-6, skip_missing::Bool = false, dual_check = true, @@ -35,11 +37,18 @@ The additional parameters: - `dual_point`: The dual solution point to use for feasibility checking. If `nothing` and the model can be dualized, it will use the current dual solution from the model. +- `primal_objective`: The primal objective value considered as the solver + objective. If `nothing`, it will use the current primal objective value from + the model (solver). +- `dual_objective`: The dual objective value considered as the solver + objective. If `nothing`, it will use the current dual objective value from + the model (solver). - `atol`: The absolute tolerance for feasibility checking. - `skip_missing`: If `true`, constraints with missing variables in the provided point will be ignored. - `dual_check`: If `true`, it will perform dual feasibility checking. Disabling - the dual check will also disable complementarity checking. + the dual check will also disable complementarity checking and dual objective + checks. """ struct Analyzer <: ModelAnalyzer.AbstractAnalyzer end @@ -228,6 +237,8 @@ Base.@kwdef mutable struct Data <: ModelAnalyzer.AbstractData # analysis configuration primal_point::Union{Nothing,AbstractDict} dual_point::Union{Nothing,AbstractDict} + primal_objective::Union{Nothing,Float64} + dual_objective::Union{Nothing,Float64} atol::Float64 skip_missing::Bool dual_check::Bool @@ -878,6 +889,8 @@ function ModelAnalyzer.analyze( model::MOI.ModelLike; primal_point = nothing, dual_point = nothing, + primal_objective::Union{Nothing,Float64} = nothing, + dual_objective::Union{Nothing,Float64} = nothing, atol::Float64 = 1e-6, skip_missing::Bool = false, dual_check = true, @@ -896,6 +909,8 @@ function ModelAnalyzer.analyze( data = Data( primal_point = primal_point, dual_point = dual_point, + primal_objective = primal_objective, + dual_objective = dual_objective, atol = atol, skip_missing = skip_missing, dual_check = dual_check, @@ -1021,14 +1036,16 @@ function _dual_point_to_dual_model_ref( Dict{MOI.ConstraintIndex,Vector{MOI.VariableIndex}}() for primal_var in primal_vars dual_con, idx = Dualization._get_dual_constraint(map, primal_var) - # TODO + @assert idx == -1 idx = max(idx, 1) if haskey(dual_con_to_primal_vars, dual_con) - vec = dual_con_to_primal_vars[dual_con] - if idx > length(vec) - resize!(vec, idx) - end - vec[idx] = primal_var + # TODO: this should never be hit because there will be no primal + # constrained variables. + # vec = dual_con_to_primal_vars[dual_con] + # if idx > length(vec) + # resize!(vec, idx) + # end + # vec[idx] = primal_var else vec = Vector{MOI.VariableIndex}(undef, idx) vec[idx] = primal_var @@ -1076,13 +1093,9 @@ function _analyze_dual!(model, dual_model, map, data) if dist > data.atol if haskey(dual_con_to_primal_vars, con) vars = dual_con_to_primal_vars[con] - if length(vars) != 1 - # TODO improve error - error( - "The dual constraint $con has " * - "length $(length(vars)) != 1", - ) - end + # TODO: this must be true because we dont consider + # constrained variables in the dualization. + @assert length(vars) == 1 push!(data.dual, DualConstraintViolation(vars[], dist)) else con = dual_con_to_primal_con[con] @@ -1180,13 +1193,17 @@ end function _analyze_objectives!(model::MOI.ModelLike, dual_model, map, data) primal_status = MOI.get(model, MOI.PrimalStatus()) dual_status = MOI.get(model, MOI.DualStatus()) - if primal_status in (MOI.FEASIBLE_POINT, MOI.NEARLY_FEASIBLE_POINT) + if data.primal_objective !== nothing + obj_val_solver = data.primal_objective + elseif primal_status in (MOI.FEASIBLE_POINT, MOI.NEARLY_FEASIBLE_POINT) obj_val_solver = MOI.get(model, MOI.ObjectiveValue()) else obj_val_solver = nothing end - if dual_status in (MOI.FEASIBLE_POINT, MOI.NEARLY_FEASIBLE_POINT) + if data.dual_objective !== nothing + dual_obj_val_solver = data.dual_objective + elseif dual_status in (MOI.FEASIBLE_POINT, MOI.NEARLY_FEASIBLE_POINT) dual_obj_val_solver = MOI.get(model, MOI.DualObjectiveValue()) else dual_obj_val_solver = nothing diff --git a/test/feasibility.jl b/test/feasibility.jl index 0d65111..793e56d 100644 --- a/test/feasibility.jl +++ b/test/feasibility.jl @@ -956,6 +956,80 @@ function test_dual_vector() return end +function test_feasibility_sense() + model = Model() + set_silent(model) + @variable(model, x) + @constraint(model, c1, x >= 0) + + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(JuMP.index(x) => 0.0), + # dual_point = Dict(JuMP.index(c1) => [0.0, 0.5]), + ) + list = ModelAnalyzer.list_of_issue_types(data) + @test length(list) == 0 + return +end + +function test_objective() + model = Model() + set_silent(model) + @variable(model, x) + @constraint(model, c, x >= 0) + @objective(model, Min, x) + + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(JuMP.index(x) => 0.0), + dual_point = Dict(JuMP.index(c) => 1.0), + primal_objective = 0.0, + dual_objective = 0.0, + ) + list = ModelAnalyzer.list_of_issue_types(data) + @test length(list) == 0 + + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(JuMP.index(x) => 0.0), + dual_point = Dict(JuMP.index(c) => 1.0), + primal_objective = 1.0, + dual_objective = 0.0, + ) + list = ModelAnalyzer.list_of_issue_types(data) + @test length(list) == 2 + @test ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.PrimalObjectiveMismatch, + )[] == ModelAnalyzer.Feasibility.PrimalObjectiveMismatch(0.0, 1.0) + @test ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.PrimalDualSolverMismatch, + )[] == ModelAnalyzer.Feasibility.PrimalDualSolverMismatch(1.0, 0.0) + data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + model, + primal_point = Dict(JuMP.index(x) => 0.0), + dual_point = Dict(JuMP.index(c) => 1.0), + primal_objective = 0.0, + dual_objective = 1.0, + ) + list = ModelAnalyzer.list_of_issue_types(data) + @test length(list) == 2 + @test ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.DualObjectiveMismatch, + )[] == ModelAnalyzer.Feasibility.DualObjectiveMismatch(0.0, 1.0) + @test ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.PrimalDualSolverMismatch, + )[] == ModelAnalyzer.Feasibility.PrimalDualSolverMismatch(0.0, 1.0) + return +end + function test_nl_con() model = Model() set_silent(model) From ad37c029b28b3bdfa8b3e798c7c21b38225db9e1 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Sat, 24 May 2025 01:17:46 -0700 Subject: [PATCH 33/34] add tests --- test/numerical.jl | 37 +++++++++++++++++++++++++++++++++++-- 1 file changed, 35 insertions(+), 2 deletions(-) diff --git a/test/numerical.jl b/test/numerical.jl index 8a072c6..3d5d6df 100644 --- a/test/numerical.jl +++ b/test/numerical.jl @@ -192,6 +192,22 @@ function test_constraint_bounds_quad() return end +function test_constraint_bounds_quad() + model = Model() + @variable(model, x) + @constraint(model, c, [-x^2 - 3] in MOI.Nonpositives(1)) + data = ModelAnalyzer.analyze(ModelAnalyzer.Numerical.Analyzer(), model) + list = ModelAnalyzer.list_of_issue_types(data) + @test length(list) == 1 + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Numerical.NonconvexQuadraticConstraint, + ) + @test length(ret) == 1 + @test ModelAnalyzer.constraint(ret[]) == JuMP.index(c) + return +end + function test_no_names() model = Model() set_string_names_on_creation(model, false) @@ -1044,10 +1060,27 @@ function test_vi_in_nonstandard_set() model = Model() @variable(model, x[1:1]) @constraint(model, c, x[1] in MOI.ZeroOne()) - @constraint(model, c2, 2x[1] == 0) + @constraint(model, c1, 3x[1] + 1e-9 in MOI.ZeroOne()) + @constraint(model, c2, 4x[1] - 1e+9 in MOI.ZeroOne()) + @constraint(model, c3, 2x[1] == 0) data = ModelAnalyzer.analyze(ModelAnalyzer.Numerical.Analyzer(), model) list = ModelAnalyzer.list_of_issue_types(data) - @test length(list) == 0 + @test length(list) == 2 + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Numerical.SmallRHSCoefficient, + ) + @test length(ret) == 1 + @test ModelAnalyzer.constraint(ret[], model) == c1 + @test ModelAnalyzer.value(ret[]) == 1e-9 + # + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Numerical.LargeRHSCoefficient, + ) + @test length(ret) == 1 + @test ModelAnalyzer.constraint(ret[], model) == c2 + @test ModelAnalyzer.value(ret[]) == -1e+9 return end From 2fc4f3fc9f40ae564b9aad0db8f2d5af232339f8 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Sat, 24 May 2025 01:32:51 -0700 Subject: [PATCH 34/34] add docs --- docs/src/analyzer.md | 12 ++++++++++++ docs/src/index.md | 14 ++++++++++++++ 2 files changed, 26 insertions(+) diff --git a/docs/src/analyzer.md b/docs/src/analyzer.md index ab8c463..e2854f1 100644 --- a/docs/src/analyzer.md +++ b/docs/src/analyzer.md @@ -22,3 +22,15 @@ and summarize them individually. The following functions are useful for this: ModelAnalyzer.list_of_issue_types ModelAnalyzer.list_of_issues ``` + +It is possible to extract data from the issues with the methods: + +```@docs +ModelAnalyzer.variables +ModelAnalyzer.variable +ModelAnalyzer.constraints +ModelAnalyzer.constraint +ModelAnalyzer.set +ModelAnalyzer.values +ModelAnalyzer.value +``` diff --git a/docs/src/index.md b/docs/src/index.md index b18086b..001b9c1 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -129,3 +129,17 @@ ModelAnalyzer.summarize(issues) # individual issues can also be summarized ModelAnalyzer.summarize(issues[1]) ``` + +### Non JuMP (or MOI) models + +If you dont have a JuMP (or MOI) model, you can still use this package reading from a file. + +```julia +model = Model(); +@variable(model, x >= 0); +@objective(model, Min, 2 * x + 1); +filename = joinpath(mktempdir(), "model.mps"); +write_to_file(model, filename; generic_names = true) +new_model = read_from_file(filename; use_nlp_block = false) +print(new_model) +```