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/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/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/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) +``` diff --git a/ext/ModelAnalyzerJuMPExt/ModelAnalyzerJuMPExt.jl b/ext/ModelAnalyzerJuMPExt/ModelAnalyzerJuMPExt.jl new file mode 100644 index 0000000..25eafae --- /dev/null +++ b/ext/ModelAnalyzerJuMPExt/ModelAnalyzerJuMPExt.jl @@ -0,0 +1,95 @@ +# 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. + +module ModelAnalyzerJuMPExt + +import ModelAnalyzer +import JuMP +import MathOptInterface as MOI + +function ModelAnalyzer.analyze( + analyzer::T, + model::JuMP.GenericModel; + kwargs..., +) where {T<:ModelAnalyzer.AbstractAnalyzer} + moi_model = JuMP.backend(model) + result = ModelAnalyzer.analyze(analyzer, moi_model; kwargs...) + return result +end + +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 + end + return "$jump_ref" +end + +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) + return name + end + return "$jump_ref" +end + +""" + 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.GenericModel{T}, +) where {T} + ref = ModelAnalyzer.variable(issue) + return JuMP.GenericVariableRef{T}(model, ref) +end + +""" + 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.GenericModel{T}, +) where {T} + refs = ModelAnalyzer.variables(issue) + return JuMP.GenericVariableRef{T}.(model, refs) +end + +""" + 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.GenericModel, +) + ref = ModelAnalyzer.constraint(issue) + 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 3c2697f..c913b18 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 @@ -12,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. @@ -25,10 +27,12 @@ See [`summarize`](@ref), [`list_of_issues`](@ref), and function analyze end """ - summarize([io::IO,] AbstractData; verbose = true, max_issues = 10, kwargs...) + summarize([io::IO,] AbstractData; model = nothing, verbose = true, max_issues = 10, 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`. +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 @@ -41,9 +45,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; model = nothing, verbose = true) This variant allows summarizing a single issue instance of type `AbstractIssue`. +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 @@ -76,11 +82,16 @@ function summarize(::Type{T}; kwargs...) where {T<:AbstractIssue} return summarize(stdout, T; kwargs...) end -function summarize(io::IO, issue::AbstractIssue; verbose = true) +function summarize( + io::IO, + issue::AbstractIssue; + model = nothing, + verbose = true, +) if verbose - return _verbose_summarize(io, issue) + return _verbose_summarize(io, issue, model) else - return _summarize(io, issue) + return _summarize(io, issue, model) end end @@ -93,6 +104,7 @@ const DEFAULT_MAX_ISSUES = 10 function summarize( io::IO, issues::Vector{T}; + model = nothing, verbose = true, max_issues = DEFAULT_MAX_ISSUES, ) where {T<:AbstractIssue} @@ -110,7 +122,7 @@ function summarize( end for issue in first(issues, max_issues) print(io, " * ") - summarize(io, issue, verbose = verbose) + summarize(io, issue, verbose = verbose, model = model) print(io, "\n") end return @@ -124,11 +136,91 @@ 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 end + +""" + values(issue::AbstractIssue) + +Return the values associated to a particular issue. +""" +function values end + +""" + variable(issue::AbstractIssue) + +Return the variable associated to a particular issue. +""" +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, ::MOI.ModelLike) + return variables(issue) +end + +""" + constraint(issue::AbstractIssue) + +Return the constraint associated to a particular issue. +""" +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, ::MOI.ModelLike) + return constraints(issue) +end + +""" + set(issue::AbstractIssue) + +Return the set associated to a particular issue. +""" +function set 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") -end +end # module ModelAnalyzer diff --git a/src/_eval_variables.jl b/src/_eval_variables.jl new file mode 100644 index 0000000..de0bfab --- /dev/null +++ b/src/_eval_variables.jl @@ -0,0 +1,24 @@ +# Copyright (c) 2025: Joaquim Garcia, Oscar Dowson and contributors +# +# Use of this source code is governed by an MIT-style license that can be found +# in the LICENSE.md file or at https://opensource.org/licenses/MIT. + +function _eval_variables end + +function _eval_variables(value_fn::Function, t::MOI.ScalarAffineTerm) + return t.coefficient * value_fn(t.variable) +end + +function _eval_variables( + value_fn::Function, + f::MOI.ScalarAffineFunction{T}, +) where {T} + # TODO: this conversion exists in JuMP, but not in MOI + S = Base.promote_op(value_fn, MOI.VariableIndex) + U = MOI.MA.promote_operation(*, T, S) + out = convert(U, f.constant) + for t in f.terms + out += _eval_variables(value_fn, t) + end + return out +end diff --git a/src/feasibility.jl b/src/feasibility.jl index 089eaf2..654f22f 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 @@ -30,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, @@ -42,18 +37,25 @@ 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 """ 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,26 +71,62 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Feasibility.PrimalViolation) ``` """ struct PrimalViolation <: AbstractFeasibilityIssue - ref::JuMP.ConstraintRef + ref::MOI.ConstraintIndex violation::Float64 end +ModelAnalyzer.constraint(issue::PrimalViolation) = issue.ref + +ModelAnalyzer.value(issue::PrimalViolation) = issue.violation + """ - DualViolation <: AbstractFeasibilityIssue + DualConstraintViolation <: AbstractFeasibilityIssue -The `DualViolation` issue is identified when a constraint has a dual value +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 + +ModelAnalyzer.variable(issue::DualConstraintViolation) = issue.ref + +ModelAnalyzer.value(issue::DualConstraintViolation) = issue.violation + +""" + DualConstrainedVariableViolation <: AbstractFeasibilityIssue + +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. +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 -julia> ModelAnalyzer.summarize(ModelAnalyzer.Feasibility.DualViolation) +julia> ModelAnalyzer.summarize(ModelAnalyzer.Feasibility.DualConstrainedVariableViolation) ``` """ -struct DualViolation <: AbstractFeasibilityIssue - ref::Union{JuMP.ConstraintRef,JuMP.VariableRef} +struct DualConstrainedVariableViolation <: AbstractFeasibilityIssue + ref::MOI.ConstraintIndex violation::Float64 end +ModelAnalyzer.constraint(issue::DualConstrainedVariableViolation) = issue.ref + +ModelAnalyzer.value(issue::DualConstrainedVariableViolation) = issue.violation + """ ComplemetarityViolation <: AbstractFeasibilityIssue @@ -103,10 +141,14 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Feasibility.ComplemetarityViolation ``` """ struct ComplemetarityViolation <: AbstractFeasibilityIssue - ref::JuMP.ConstraintRef + ref::MOI.ConstraintIndex violation::Float64 end +ModelAnalyzer.constraint(issue::ComplemetarityViolation) = issue.ref + +ModelAnalyzer.value(issue::ComplemetarityViolation) = issue.violation + """ DualObjectiveMismatch <: AbstractFeasibilityIssue @@ -124,6 +166,8 @@ struct DualObjectiveMismatch <: AbstractFeasibilityIssue obj_solver::Float64 end +# ModelAnalyzer.values(issue::DualObjectiveMismatch) = [issue.obj, issue.obj_solver] + """ PrimalObjectiveMismatch <: AbstractFeasibilityIssue @@ -141,6 +185,8 @@ struct PrimalObjectiveMismatch <: AbstractFeasibilityIssue obj_solver::Float64 end +# ModelAnalyzer.values(issue::PrimalObjectiveMismatch) = [issue.obj, issue.obj_solver] + """ PrimalDualMismatch <: AbstractFeasibilityIssue @@ -158,6 +204,8 @@ struct PrimalDualMismatch <: AbstractFeasibilityIssue dual::Float64 end +ModelAnalyzer.values(issue::PrimalDualMismatch) = [issue.primal, issue.dual] + """ PrimalDualSolverMismatch <: AbstractFeasibilityIssue @@ -175,11 +223,13 @@ struct PrimalDualSolverMismatch <: AbstractFeasibilityIssue dual::Float64 end +# ModelAnalyzer.values(issue::PrimalDualSolverMismatch) = [issue.primal, issue.dual] + """ 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. """ @@ -187,12 +237,16 @@ 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 # 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 +262,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 +332,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 +373,53 @@ 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. + 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 + + 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 +600,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 +730,7 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::PrimalObjectiveMismatch, + model, ) return print( io, @@ -585,7 +742,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 +760,7 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::PrimalDualSolverMismatch, + model, ) return print( io, @@ -614,10 +776,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 +819,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.VariableRef) - 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 +849,7 @@ end function ModelAnalyzer.summarize( io::IO, data::Data; + model = nothing, verbose = true, max_issues = ModelAnalyzer.DEFAULT_MAX_ISSUES, configurations = true, @@ -701,6 +867,7 @@ function ModelAnalyzer.summarize( ModelAnalyzer.summarize( io, issues, + model = model, verbose = verbose, max_issues = max_issues, ) @@ -719,26 +886,39 @@ end function ModelAnalyzer.analyze( ::Analyzer, - model::JuMP.GenericModel; + 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, ) + 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, + primal_objective = primal_objective, + dual_objective = dual_objective, atol = atol, skip_missing = skip_missing, dual_check = dual_check, ) 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 +927,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 +938,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 +957,209 @@ 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 (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 + 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) + @assert idx == -1 + idx = max(idx, 1) + if haskey(dual_con_to_primal_vars, dual_con) + # 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 + 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] + # 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] + 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 +1190,25 @@ 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 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 JuMP.dual_status(model) in - (MOI.FEASIBLE_POINT, MOI.NEARLY_FEASIBLE_POINT) - dual_obj_val_solver = JuMP.dual_objective_value(model) + + 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 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 +1218,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) @@ -904,7 +1272,9 @@ function _analyze_objectives!( ) 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), @@ -915,269 +1285,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 -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.VariableRef(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( - 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 +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 _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.VariableRef},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.VariableRef(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.Model() - 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::MOI.ModelLike) + types = MOI.get(model, MOI.ListOfConstraintTypesPresent()) -function _can_dualize(model::JuMP.GenericModel) - types = JuMP.list_of_constraint_types(model) - - 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/infeasibility.jl b/src/infeasibility.jl index 2680f35..2dab25d 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,11 +50,15 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Infeasibility.InfeasibleBounds) ```` """ struct InfeasibleBounds{T} <: AbstractInfeasibilitylIssue - variable::JuMP.VariableRef + variable::MOI.VariableIndex lb::T ub::T end +ModelAnalyzer.variable(issue::InfeasibleBounds) = issue.variable + +ModelAnalyzer.values(issue::InfeasibleBounds) = [issue.lb, issue.ub] + """ InfeasibleIntegrality{T} <: AbstractInfeasibilitylIssue @@ -71,12 +75,18 @@ julia> ModelAnalyzer.summarize( ``` """ struct InfeasibleIntegrality{T} <: AbstractInfeasibilitylIssue - variable::JuMP.VariableRef + variable::MOI.VariableIndex lb::T ub::T 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 @@ -95,12 +105,18 @@ 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}} end +ModelAnalyzer.constraint(issue::InfeasibleConstraintRange) = issue.constraint + +ModelAnalyzer.values(issue::InfeasibleConstraintRange) = [issue.lb, issue.ub] + +ModelAnalyzer.set(issue::InfeasibleConstraintRange) = issue.set + """ IrreducibleInfeasibleSubset <: AbstractInfeasibilitylIssue @@ -118,9 +134,11 @@ julia> ModelAnalyzer.summarize( ``` """ struct IrreducibleInfeasibleSubset <: AbstractInfeasibilitylIssue - constraint::Vector{JuMP.ConstraintRef} + constraint::Vector{<:MOI.ConstraintIndex} end +ModelAnalyzer.constraints(issue::IrreducibleInfeasibleSubset) = issue.constraint + """ Data <: ModelAnalyzer.AbstractData @@ -142,55 +160,119 @@ end function ModelAnalyzer.analyze( ::Analyzer, - model::JuMP.GenericModel{T}; + model::MOI.ModelLike; optimizer = nothing, -) where {T} +) out = Data() - variables = Dict{JuMP.VariableRef,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 +281,97 @@ 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 - # TODO: handle quadratics - !(S in (MOI.EqualTo{T}, MOI.LessThan{T}, MOI.GreaterThan{T})) && - 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 + + for con in MOI.get( + model, + MOI.ListOfConstraintIndices{ + MOI.ScalarAffineFunction{T}, + MOI.EqualTo{T}, + }(), + ) + set = MOI.get(model, MOI.ConstraintSet(), con) + func = MOI.get(model, MOI.ConstraintFunction(), con) + failed = false + interval = _eval_variables(func) do var_idx + # this only fails if we allow continuing after bounds issues + # if !haskey(variables, var_idx) + # failed = true + # return Interval(-Inf, Inf) + # end + return variables[var_idx] + end + # if failed + # continue + # end + rhs = set.value + if interval.lo > rhs || interval.hi < rhs + push!( + out.constraint_range, + InfeasibleConstraintRange(con, interval.lo, interval.hi, set), + ) + range_consistent = false + end + end + + for con in MOI.get( + model, + MOI.ListOfConstraintIndices{ + MOI.ScalarAffineFunction{T}, + MOI.LessThan{T}, + }(), + ) + set = MOI.get(model, MOI.ConstraintSet(), con) + func = MOI.get(model, MOI.ConstraintFunction(), con) + failed = false + interval = _eval_variables(func) do var_idx + # this only fails if we allow continuing after bounds issues + # if !haskey(variables, var_idx) + # failed = true + # return Interval(-Inf, Inf) + # end + return variables[var_idx] + end + # if failed + # continue + # end + rhs = set.upper + if interval.lo > rhs + push!( + out.constraint_range, + InfeasibleConstraintRange(con, interval.lo, interval.hi, set), + ) + range_consistent = false + end + end + + for con in MOI.get( + model, + MOI.ListOfConstraintIndices{ + MOI.ScalarAffineFunction{T}, + MOI.GreaterThan{T}, + }(), + ) + set = MOI.get(model, MOI.ConstraintSet(), con) + func = MOI.get(model, MOI.ConstraintFunction(), con) + failed = false + interval = _eval_variables(func) do var_idx + # this only fails if we allow continuing after bounds issues + # if !haskey(variables, var_idx) + # failed = true + # return Interval(-Inf, Inf) + # end + return variables[var_idx] + end + # if failed + # continue + # end + rhs = set.lower + if interval.hi < rhs + push!( + out.constraint_range, + InfeasibleConstraintRange(con, interval.lo, interval.hi, set), + ) + range_consistent = false end end @@ -272,14 +394,50 @@ function ModelAnalyzer.analyze( return out end -function iis_elastic_filter(original_model::JuMP.GenericModel, optimizer) +function _fix_to_zero(model, variable::MOI.VariableIndex, ::Type{T}) where {T} + ub_idx = + MOI.ConstraintIndex{MOI.VariableIndex,MOI.LessThan{T}}(variable.value) + lb_idx = MOI.ConstraintIndex{MOI.VariableIndex,MOI.GreaterThan{T}}( + variable.value, + ) + has_lower = false + if MOI.is_valid(model, lb_idx) + MOI.delete(model, lb_idx) + has_lower = true + # MOI.PenaltyRelaxation only creates variables with LB + # elseif MOI.is_valid(model, ub_idx) + # MOI.delete(model, ub_idx) + else + error("Variable is not bounded") + end + MOI.add_constraint(model, variable, MOI.EqualTo{T}(zero(T))) + return has_lower +end - # 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 _set_bound_zero( + model, + variable::MOI.VariableIndex, + has_lower::Bool, + ::Type{T}, +) where {T} + eq_idx = + MOI.ConstraintIndex{MOI.VariableIndex,MOI.EqualTo{T}}(variable.value) + @assert MOI.is_valid(model, eq_idx) + MOI.delete(model, eq_idx) + if has_lower + MOI.add_constraint(model, variable, MOI.GreaterThan{T}(zero(T))) + # MOI.PenaltyRelaxation only creates variables with LB + # else + # MOI.add_constraint(model, variable, MOI.LessThan{T}(zero(T))) + end + return +end + +function iis_elastic_filter(original_model::MOI.ModelLike, optimizer) + T = Float64 - 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", @@ -287,16 +445,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 @@ -304,47 +460,46 @@ 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 + # TODO: coef is alwayas 1.0 + 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 @@ -352,38 +507,32 @@ 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 @@ -536,17 +685,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 +720,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 +734,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 +765,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 +783,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 +801,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 +850,7 @@ end function ModelAnalyzer.summarize( io::IO, data::Data; + model = nothing, verbose = true, max_issues = ModelAnalyzer.DEFAULT_MAX_ISSUES, ) @@ -687,6 +862,7 @@ function ModelAnalyzer.summarize( ModelAnalyzer.summarize( io, issues, + model = model, verbose = verbose, max_issues = max_issues, ) @@ -703,14 +879,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..f54eb84 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 @@ -62,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 diff --git a/src/numerical.jl b/src/numerical.jl index 89cf660..8eb904a 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,9 +59,11 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.VariableNotInConstraints) ``` """ struct VariableNotInConstraints <: AbstractNumericalIssue - ref::JuMP.VariableRef + ref::MOI.VariableIndex end +ModelAnalyzer.variable(issue::VariableNotInConstraints) = issue.ref + """ EmptyConstraint <: AbstractNumericalIssue @@ -75,9 +76,11 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.EmptyConstraint) ``` """ struct EmptyConstraint <: AbstractNumericalIssue - ref::JuMP.ConstraintRef + ref::MOI.ConstraintIndex end +ModelAnalyzer.constraint(issue::EmptyConstraint) = issue.ref + """ VariableBoundAsConstraint <: AbstractNumericalIssue @@ -91,9 +94,11 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.VariableBoundAsConstraint ``` """ struct VariableBoundAsConstraint <: AbstractNumericalIssue - ref::JuMP.ConstraintRef + ref::MOI.ConstraintIndex end +ModelAnalyzer.constraint(issue::VariableBoundAsConstraint) = issue.ref + """ DenseConstraint <: AbstractNumericalIssue @@ -107,10 +112,14 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.DenseConstraint) ``` """ struct DenseConstraint <: AbstractNumericalIssue - ref::JuMP.ConstraintRef + ref::MOI.ConstraintIndex nnz::Int end +ModelAnalyzer.constraint(issue::DenseConstraint) = issue.ref + +ModelAnalyzer.value(issue::DenseConstraint) = issue.nnz + """ SmallMatrixCoefficient <: AbstractNumericalIssue @@ -123,11 +132,17 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.SmallMatrixCoefficient) ``` """ struct SmallMatrixCoefficient <: AbstractNumericalIssue - ref::JuMP.ConstraintRef - variable::JuMP.VariableRef + ref::MOI.ConstraintIndex + variable::MOI.VariableIndex coefficient::Float64 end +ModelAnalyzer.variable(issue::SmallMatrixCoefficient) = issue.variable + +ModelAnalyzer.constraint(issue::SmallMatrixCoefficient) = issue.ref + +ModelAnalyzer.value(issue::SmallMatrixCoefficient) = issue.coefficient + """ LargeMatrixCoefficient <: AbstractNumericalIssue @@ -140,11 +155,17 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.LargeMatrixCoefficient) ``` """ struct LargeMatrixCoefficient <: AbstractNumericalIssue - ref::JuMP.ConstraintRef - variable::JuMP.VariableRef + ref::MOI.ConstraintIndex + variable::MOI.VariableIndex coefficient::Float64 end +ModelAnalyzer.variable(issue::LargeMatrixCoefficient) = issue.variable + +ModelAnalyzer.constraint(issue::LargeMatrixCoefficient) = issue.ref + +ModelAnalyzer.value(issue::LargeMatrixCoefficient) = issue.coefficient + """ SmallBoundCoefficient <: AbstractNumericalIssue @@ -157,10 +178,14 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.SmallBoundCoefficient) ``` """ struct SmallBoundCoefficient <: AbstractNumericalIssue - variable::JuMP.VariableRef + variable::MOI.VariableIndex coefficient::Float64 end +ModelAnalyzer.variable(issue::SmallBoundCoefficient) = issue.variable + +ModelAnalyzer.value(issue::SmallBoundCoefficient) = issue.coefficient + """ LargeBoundCoefficient <: AbstractNumericalIssue @@ -173,10 +198,14 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.LargeBoundCoefficient) ``` """ struct LargeBoundCoefficient <: AbstractNumericalIssue - variable::JuMP.VariableRef + variable::MOI.VariableIndex coefficient::Float64 end +ModelAnalyzer.variable(issue::LargeBoundCoefficient) = issue.variable + +ModelAnalyzer.value(issue::LargeBoundCoefficient) = issue.coefficient + """ SmallRHSCoefficient <: AbstractNumericalIssue @@ -189,10 +218,14 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.SmallRHSCoefficient) ``` """ struct SmallRHSCoefficient <: AbstractNumericalIssue - ref::JuMP.ConstraintRef + ref::MOI.ConstraintIndex coefficient::Float64 end +ModelAnalyzer.constraint(issue::SmallRHSCoefficient) = issue.ref + +ModelAnalyzer.value(issue::SmallRHSCoefficient) = issue.coefficient + """ LargeRHSCoefficient <: AbstractNumericalIssue @@ -205,10 +238,14 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.LargeRHSCoefficient) ``` """ struct LargeRHSCoefficient <: AbstractNumericalIssue - ref::JuMP.ConstraintRef + ref::MOI.ConstraintIndex coefficient::Float64 end +ModelAnalyzer.constraint(issue::LargeRHSCoefficient) = issue.ref + +ModelAnalyzer.value(issue::LargeRHSCoefficient) = issue.coefficient + """ SmallObjectiveCoefficient <: AbstractNumericalIssue @@ -221,10 +258,14 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.SmallObjectiveCoefficient ``` """ struct SmallObjectiveCoefficient <: AbstractNumericalIssue - variable::JuMP.VariableRef + variable::MOI.VariableIndex coefficient::Float64 end +ModelAnalyzer.variable(issue::SmallObjectiveCoefficient) = issue.variable + +ModelAnalyzer.value(issue::SmallObjectiveCoefficient) = issue.coefficient + """ LargeObjectiveCoefficient <: AbstractNumericalIssue @@ -237,10 +278,14 @@ julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.LargeObjectiveCoefficient ``` """ struct LargeObjectiveCoefficient <: AbstractNumericalIssue - variable::JuMP.VariableRef + variable::MOI.VariableIndex coefficient::Float64 end +ModelAnalyzer.variable(issue::LargeObjectiveCoefficient) = issue.variable + +ModelAnalyzer.value(issue::LargeObjectiveCoefficient) = issue.coefficient + """ SmallObjectiveQuadraticCoefficient <: AbstractNumericalIssue @@ -255,11 +300,19 @@ julia> ModelAnalyzer.summarize( ``` """ struct SmallObjectiveQuadraticCoefficient <: AbstractNumericalIssue - variable1::JuMP.VariableRef - variable2::JuMP.VariableRef + variable1::MOI.VariableIndex + variable2::MOI.VariableIndex coefficient::Float64 end +function ModelAnalyzer.variables(issue::SmallObjectiveQuadraticCoefficient) + return [issue.variable1, issue.variable2] +end + +function ModelAnalyzer.value(issue::SmallObjectiveQuadraticCoefficient) + return issue.coefficient +end + """ LargeObjectiveQuadraticCoefficient <: AbstractNumericalIssue @@ -274,11 +327,19 @@ julia> ModelAnalyzer.summarize( ``` """ struct LargeObjectiveQuadraticCoefficient <: AbstractNumericalIssue - variable1::JuMP.VariableRef - variable2::JuMP.VariableRef + variable1::MOI.VariableIndex + variable2::MOI.VariableIndex coefficient::Float64 end +function ModelAnalyzer.variables(issue::LargeObjectiveQuadraticCoefficient) + return [issue.variable1, issue.variable2] +end + +function ModelAnalyzer.value(issue::LargeObjectiveQuadraticCoefficient) + return issue.coefficient +end + """ SmallMatrixQuadraticCoefficient <: AbstractNumericalIssue @@ -293,12 +354,22 @@ 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 +function ModelAnalyzer.variables(issue::SmallMatrixQuadraticCoefficient) + return [issue.variable1, issue.variable2] +end + +function ModelAnalyzer.constraint(issue::SmallMatrixQuadraticCoefficient) + return issue.ref +end + +ModelAnalyzer.value(issue::SmallMatrixQuadraticCoefficient) = issue.coefficient + """ LargeMatrixQuadraticCoefficient <: AbstractNumericalIssue @@ -313,12 +384,22 @@ 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 +function ModelAnalyzer.variables(issue::LargeMatrixQuadraticCoefficient) + return [issue.variable1, issue.variable2] +end + +function ModelAnalyzer.constraint(issue::LargeMatrixQuadraticCoefficient) + return issue.ref +end + +ModelAnalyzer.value(issue::LargeMatrixQuadraticCoefficient) = issue.coefficient + """ NonconvexQuadraticObjective <: AbstractNumericalIssue @@ -348,8 +429,9 @@ julia> ModelAnalyzer.summarize( ``` """ struct NonconvexQuadraticConstraint <: AbstractNumericalIssue - ref::JuMP.ConstraintRef + ref::MOI.ConstraintIndex end +ModelAnalyzer.constraint(issue::NonconvexQuadraticConstraint) = issue.ref """ Data @@ -377,7 +459,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 +482,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 +506,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 +534,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 +554,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 +576,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 +605,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 +647,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 +683,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 +698,121 @@ 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_data(data, ref, func::JuMP.GenericAffExpr, set) +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::Union{MOI.ScalarAffineFunction,MOI.ScalarQuadraticFunction}, + set, +) coefficient = func.constant if iszero(coefficient) return @@ -635,20 +826,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 +887,7 @@ end function _get_constraint_data( data, ref, - func::JuMP.GenericAffExpr, + func::MOI.ScalarAffineFunction, set::MOI.LessThan, ) coefficient = set.upper - func.constant @@ -677,10 +906,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 +942,7 @@ end function _get_constraint_data( data, ref, - func::JuMP.GenericAffExpr, + func::MOI.ScalarAffineFunction, set::MOI.GreaterThan, ) coefficient = set.lower - func.constant @@ -709,10 +961,15 @@ 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 @@ -720,7 +977,23 @@ end function _get_constraint_data( data, ref, - func::JuMP.GenericAffExpr, + 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 + +function _get_constraint_data( + data, + ref, + func::MOI.ScalarAffineFunction, set::MOI.EqualTo, ) coefficient = set.value - func.constant @@ -739,7 +1012,7 @@ end function _get_constraint_data( data, ref, - func::JuMP.GenericAffExpr, + func::MOI.ScalarAffineFunction, set::MOI.Interval, ) coefficient = set.upper - func.constant @@ -764,22 +1037,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 +1126,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, @@ -1527,77 +1851,124 @@ 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, + model, +) + return print(io, ModelAnalyzer._name(issue.ref, model)) end -function ModelAnalyzer._summarize(io::IO, issue::EmptyConstraint) - return print(io, _name(issue.ref)) +function ModelAnalyzer._summarize(io::IO, issue::EmptyConstraint, model) + return print(io, ModelAnalyzer._name(issue.ref, model)) end -function ModelAnalyzer._summarize(io::IO, issue::VariableBoundAsConstraint) - return print(io, _name(issue.ref)) +function ModelAnalyzer._summarize( + io::IO, + issue::VariableBoundAsConstraint, + model, +) + return print(io, ModelAnalyzer._name(issue.ref, model)) end -function ModelAnalyzer._summarize(io::IO, issue::DenseConstraint) - return print(io, _name(issue.ref), " : ", 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) +function ModelAnalyzer._summarize(io::IO, issue::SmallMatrixCoefficient, model) return print( io, - _name(issue.ref), + ModelAnalyzer._name(issue.ref, model), " -- ", - _name(issue.variable), + ModelAnalyzer._name(issue.variable, model), " : ", issue.coefficient, ) end -function ModelAnalyzer._summarize(io::IO, issue::LargeMatrixCoefficient) +function ModelAnalyzer._summarize(io::IO, issue::LargeMatrixCoefficient, model) return print( io, - _name(issue.ref), + ModelAnalyzer._name(issue.ref, model), " -- ", - _name(issue.variable), + ModelAnalyzer._name(issue.variable, model), " : ", 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, model) + return print( + io, + ModelAnalyzer._name(issue.variable, model), + " : ", + 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, model) + return print( + io, + ModelAnalyzer._name(issue.variable, model), + " : ", + 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, model) + return print( + io, + ModelAnalyzer._name(issue.ref, model), + " : ", + 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, model) + return print( + io, + ModelAnalyzer._name(issue.ref, model), + " : ", + 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, + model, +) + return print( + io, + ModelAnalyzer._name(issue.variable, model), + " : ", + 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, + model, +) + return print( + io, + ModelAnalyzer._name(issue.variable, model), + " : ", + issue.coefficient, + ) end function ModelAnalyzer._summarize( io::IO, issue::SmallObjectiveQuadraticCoefficient, + model, ) return print( io, - _name(issue.variable1), + ModelAnalyzer._name(issue.variable1, model), " -- ", - _name(issue.variable2), + ModelAnalyzer._name(issue.variable2, model), " : ", issue.coefficient, ) @@ -1606,12 +1977,13 @@ end function ModelAnalyzer._summarize( io::IO, issue::LargeObjectiveQuadraticCoefficient, + model, ) return print( io, - _name(issue.variable1), + ModelAnalyzer._name(issue.variable1, model), " -- ", - _name(issue.variable2), + ModelAnalyzer._name(issue.variable2, model), " : ", issue.coefficient, ) @@ -1620,14 +1992,15 @@ end function ModelAnalyzer._summarize( io::IO, issue::SmallMatrixQuadraticCoefficient, + model, ) return print( io, - _name(issue.ref), + ModelAnalyzer._name(issue.ref, model), " -- ", - _name(issue.variable1), + ModelAnalyzer._name(issue.variable1, model), " -- ", - _name(issue.variable2), + ModelAnalyzer._name(issue.variable2, model), " : ", issue.coefficient, ) @@ -1636,115 +2009,146 @@ end function ModelAnalyzer._summarize( io::IO, issue::LargeMatrixQuadraticCoefficient, + model, ) return print( io, - _name(issue.ref), + ModelAnalyzer._name(issue.ref, model), " -- ", - _name(issue.variable1), + ModelAnalyzer._name(issue.variable1, model), " -- ", - _name(issue.variable2), + ModelAnalyzer._name(issue.variable2, model), " : ", issue.coefficient, ) end -function ModelAnalyzer._summarize(io::IO, ::NonconvexQuadraticObjective) +function ModelAnalyzer._summarize(io::IO, ::NonconvexQuadraticObjective, model) 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, + model, +) + return print(io, ModelAnalyzer._name(issue.ref, model)) end function ModelAnalyzer._verbose_summarize( io::IO, issue::VariableNotInConstraints, + model, ) - return print(io, "Variable: ", _name(issue.ref)) + return print(io, "Variable: ", ModelAnalyzer._name(issue.ref, model)) end -function ModelAnalyzer._verbose_summarize(io::IO, issue::EmptyConstraint) - return print(io, "Constraint: ", _name(issue.ref)) +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, + model, ) - return print(io, "Constraint: ", _name(issue.ref)) + return print(io, "Constraint: ", ModelAnalyzer._name(issue.ref, model)) end -function ModelAnalyzer._verbose_summarize(io::IO, issue::DenseConstraint) +function ModelAnalyzer._verbose_summarize(io::IO, issue::DenseConstraint, model) return print( io, "Constraint: ", - _name(issue.ref), + ModelAnalyzer._name(issue.ref, model), " with ", issue.nnz, " non zero coefficients", ) end -function ModelAnalyzer._verbose_summarize(io::IO, issue::SmallMatrixCoefficient) +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::SmallMatrixCoefficient, + model, +) return print( io, "(Constraint -- Variable): (", - _name(issue.ref), + ModelAnalyzer._name(issue.ref, model), " -- ", - _name(issue.variable), + ModelAnalyzer._name(issue.variable, model), ") with coefficient ", issue.coefficient, ) end -function ModelAnalyzer._verbose_summarize(io::IO, issue::LargeMatrixCoefficient) +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::LargeMatrixCoefficient, + model, +) return print( io, "(Constraint -- Variable): (", - _name(issue.ref), + ModelAnalyzer._name(issue.ref, model), " -- ", - _name(issue.variable), + ModelAnalyzer._name(issue.variable, model), ") with coefficient ", issue.coefficient, ) end -function ModelAnalyzer._verbose_summarize(io::IO, issue::SmallBoundCoefficient) +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::SmallBoundCoefficient, + model, +) return print( io, "Variable: ", - _name(issue.variable), + ModelAnalyzer._name(issue.variable, model), " with bound ", issue.coefficient, ) end -function ModelAnalyzer._verbose_summarize(io::IO, issue::LargeBoundCoefficient) +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::LargeBoundCoefficient, + model, +) return print( io, "Variable: ", - _name(issue.variable), + ModelAnalyzer._name(issue.variable, model), " with bound ", issue.coefficient, ) end -function ModelAnalyzer._verbose_summarize(io::IO, issue::SmallRHSCoefficient) +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::SmallRHSCoefficient, + model, +) return print( io, "Constraint: ", - _name(issue.ref), + ModelAnalyzer._name(issue.ref, model), " with right-hand-side ", issue.coefficient, ) end -function ModelAnalyzer._verbose_summarize(io::IO, issue::LargeRHSCoefficient) +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::LargeRHSCoefficient, + model, +) return print( io, "Constraint: ", - _name(issue.ref), + ModelAnalyzer._name(issue.ref, model), " with right-hand-side ", issue.coefficient, ) @@ -1753,11 +2157,12 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::SmallObjectiveCoefficient, + model, ) return print( io, "Variable: ", - _name(issue.variable), + ModelAnalyzer._name(issue.variable, model), " with coefficient ", issue.coefficient, ) @@ -1766,11 +2171,12 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::LargeObjectiveCoefficient, + model, ) return print( io, "Variable: ", - _name(issue.variable), + ModelAnalyzer._name(issue.variable, model), " with coefficient ", issue.coefficient, ) @@ -1779,13 +2185,14 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::SmallObjectiveQuadraticCoefficient, + model, ) return print( io, "(Variable -- Variable): (", - _name(issue.variable1), + ModelAnalyzer._name(issue.variable1, model), " -- ", - _name(issue.variable2), + ModelAnalyzer._name(issue.variable2, model), ") with coefficient ", issue.coefficient, ) @@ -1794,13 +2201,14 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::LargeObjectiveQuadraticCoefficient, + model, ) return print( io, "(Variable -- Variable): (", - _name(issue.variable1), + ModelAnalyzer._name(issue.variable1, model), " -- ", - _name(issue.variable2), + ModelAnalyzer._name(issue.variable2, model), ") with coefficient ", issue.coefficient, ) @@ -1809,15 +2217,16 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::SmallMatrixQuadraticCoefficient, + model, ) return print( io, "(Constraint -- Variable -- Variable): (", - _name(issue.ref), + ModelAnalyzer._name(issue.ref, model), " -- ", - _name(issue.variable1), + ModelAnalyzer._name(issue.variable1, model), " -- ", - _name(issue.variable2), + ModelAnalyzer._name(issue.variable2, model), ") with coefficient ", issue.coefficient, ) @@ -1826,15 +2235,16 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::LargeMatrixQuadraticCoefficient, + model, ) return print( io, "(Constraint -- Variable -- Variable): (", - _name(issue.ref), + ModelAnalyzer._name(issue.ref, model), " -- ", - _name(issue.variable1), + ModelAnalyzer._name(issue.variable1, model), " -- ", - _name(issue.variable2), + ModelAnalyzer._name(issue.variable2, model), ") with coefficient ", issue.coefficient, ) @@ -1843,15 +2253,17 @@ end function ModelAnalyzer._verbose_summarize( io::IO, issue::NonconvexQuadraticObjective, + model, ) - return ModelAnalyzer._summarize(io, issue) + return ModelAnalyzer._summarize(io, issue, model) end function ModelAnalyzer._verbose_summarize( io::IO, issue::NonconvexQuadraticConstraint, + model, ) - return print(io, "Constraint: ", _name(issue.ref)) + return print(io, "Constraint: ", ModelAnalyzer._name(issue.ref, model)) end function ModelAnalyzer.list_of_issues( @@ -2041,6 +2453,7 @@ end function ModelAnalyzer.summarize( io::IO, data::Data; + model = nothing, verbose = true, max_issues = ModelAnalyzer.DEFAULT_MAX_ISSUES, configurations = true, @@ -2066,6 +2479,7 @@ function ModelAnalyzer.summarize( ModelAnalyzer.summarize( io, issues, + model = model, verbose = verbose, max_issues = max_issues, ) @@ -2084,14 +2498,6 @@ end # printing helpers -function _name(ref) - name = JuMP.name(ref) - if !isempty(name) - return name - end - return "$(ref.index)" -end - _print_value(x::Real) = Printf.@sprintf("%1.0e", x) function _stringify_bounds(bounds::Vector{Float64}) diff --git a/test/feasibility.jl b/test/feasibility.jl index 93a46e4..793e56d 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,165 @@ 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, - Dict(LowerBoundRef(x) => 1.0), + primal_point = Dict(JuMP.index(x) => -1.0), + dual_point = Dict(JuMP.index(LowerBoundRef(x)) => 1.0), ) - @test isempty(report) + 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, + ) + @test ModelAnalyzer.constraint(ret[], model) == LowerBoundRef(x) + @test ModelAnalyzer.value(ret[]) == 1.0 + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.ComplemetarityViolation, + ) + @test ret[] == ModelAnalyzer.Feasibility.ComplemetarityViolation( + 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(), + 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.ComplemetarityViolation, + ModelAnalyzer.Feasibility.PrimalDualMismatch, + ] + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.ComplemetarityViolation, + ) + @test ret[] == ModelAnalyzer.Feasibility.ComplemetarityViolation( + 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) + + 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) + @test ModelAnalyzer.variable(ret[], model) == x + @test ModelAnalyzer.value(ret[]) == 2.0 + 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 + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.ComplemetarityViolation, + ) + @test ret[] == ModelAnalyzer.Feasibility.ComplemetarityViolation( + 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) return end @@ -64,24 +195,122 @@ 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, - Dict(c => [1.0]), + primal_point = Dict(JuMP.index(x) => 0.0), + dual_point = Dict(JuMP.index(c) => 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, + 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 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, + ) + @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), ) - @test report[x] == 2.3 - @test length(report) == 1 + 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) + @test ModelAnalyzer.constraint(ret[], model) == c + @test ModelAnalyzer.value(ret[]) == 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, + 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) + @test ModelAnalyzer.variable(ret[], model) == x + @test ModelAnalyzer.value(ret[]) == 2.0 + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.DualConstrainedVariableViolation, + ) + @test ret[] == ModelAnalyzer.Feasibility.DualConstrainedVariableViolation( + 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, + ) + @test ret[] == ModelAnalyzer.Feasibility.PrimalDualMismatch(1.0, 0.0) end function test_lb0() @@ -92,25 +321,76 @@ 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, + ), ) - @test report[LowerBoundRef(x)] == 2.3 - @test length(report) == 1 - report = ModelAnalyzer.Feasibility.dual_feasibility_report( + 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 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(), 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, + ), ) - @test report[c] == 3.3 - @test length(report) == 1 + 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) + @test ModelAnalyzer.variable(ret[], model) == x + @test ModelAnalyzer.value(ret[]) == 4.3 + ret = ModelAnalyzer.list_of_issues( + data, + ModelAnalyzer.Feasibility.DualConstrainedVariableViolation, + ) + @test ret[] == ModelAnalyzer.Feasibility.DualConstrainedVariableViolation( + JuMP.index(c), + 3.3, + ) + @test ModelAnalyzer.constraint(ret[], model) == c + @test ModelAnalyzer.value(ret[]) == 3.3 + return end function test_lb2() @@ -121,42 +401,118 @@ 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) + @test ModelAnalyzer.variable(ret[], model) == x + @test ModelAnalyzer.value(ret[]) == 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, + ), ) - @test report[x] == 4.3 - @test report[c] == 3.3 - @test length(report) == 2 - report = ModelAnalyzer.Feasibility.dual_feasibility_report( + 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 ModelAnalyzer.constraint(ret[], model) == c + @test ModelAnalyzer.value(ret[]) == 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, + ), + ) + 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 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 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( + @test ModelAnalyzer.constraint(ret[2], model) == LowerBoundRef(x) + @test ModelAnalyzer.value(ret[2]) == 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, + ) + @test ModelAnalyzer.constraint(ret[], model) == c + @test ModelAnalyzer.value(ret[]) == 3.3 + return end function test_analyse_simple() @@ -202,9 +558,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 @@ -278,14 +632,14 @@ function test_analyse_no_opt() @test_throws ErrorException ModelAnalyzer.analyze( ModelAnalyzer.Feasibility.Analyzer(), model, - primal_point = Dict(x => 1.0), + primal_point = Dict(JuMP.index(x) => 1.0), dual_check = true, ) data = ModelAnalyzer.analyze( ModelAnalyzer.Feasibility.Analyzer(), model, - primal_point = Dict(x => 1.0), + primal_point = Dict(JuMP.index(x) => 1.0), dual_check = false, ) list = ModelAnalyzer.list_of_issue_types(data) @@ -294,32 +648,37 @@ function test_analyse_no_opt() data = ModelAnalyzer.analyze( ModelAnalyzer.Feasibility.Analyzer(), model, - primal_point = Dict(x => -1.0), + 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(c, 1.0) + @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(), model, - primal_point = Dict(x => 1.0), - dual_point = 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 length(list) == 2 ret = ModelAnalyzer.list_of_issues(data, list[1]) - @test ret[1] == ModelAnalyzer.Feasibility.ComplemetarityViolation(c, 1.0) + @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) data = ModelAnalyzer.analyze( ModelAnalyzer.Feasibility.Analyzer(), model, - primal_point = Dict(x => 0.0), - dual_point = Dict(c => 1.0), + 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 @@ -327,17 +686,26 @@ function test_analyse_no_opt() data = ModelAnalyzer.analyze( ModelAnalyzer.Feasibility.Analyzer(), model, - primal_point = Dict(x => -1.0), - dual_point = Dict(c => 2.0), + 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(c, 1.0) + @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.DualViolation(x, 1.0) + @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(c, -2.0) + @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) @@ -350,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() @@ -383,6 +796,276 @@ 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 + 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, + ) + 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, + ) + 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, + ) + 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]), + ) + list = ModelAnalyzer.list_of_issue_types(data) + @test length(list) == 0 + 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) + @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), + ) + 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), + ) + list = ModelAnalyzer.list_of_issue_types(data) + @test length(list) == 0 + return +end + end # module TestDualFeasibilityChecker.runtests() diff --git a/test/infeasibility.jl b/test/infeasibility.jl index a50fee6..cde2582 100644 --- a/test/infeasibility.jl +++ b/test/infeasibility.jl @@ -32,8 +32,13 @@ 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, + ) + @test ModelAnalyzer.variable(ret[], model) == y + @test ModelAnalyzer.values(ret[]) == [2.0, 1.0] # buf = IOBuffer() ModelAnalyzer.summarize( @@ -74,11 +79,14 @@ 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(), ) + @test ModelAnalyzer.variable(ret[], model) == y + @test ModelAnalyzer.values(ret[]) == [2.2, 2.9] + @test ModelAnalyzer.set(ret[]) == MOI.Integer() # buf = IOBuffer() ModelAnalyzer.summarize( @@ -121,11 +129,14 @@ 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(), ) + @test ModelAnalyzer.variable(ret[], model) == x + @test ModelAnalyzer.values(ret[]) == [0.5, 0.8] + @test ModelAnalyzer.set(ret[]) == MOI.ZeroOne() return end @@ -142,11 +153,14 @@ 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), ) + @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( @@ -190,11 +204,14 @@ 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), ) + @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 @@ -211,11 +228,14 @@ 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), ) + @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 @@ -232,11 +252,14 @@ 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), ) + @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 @@ -253,11 +276,14 @@ 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), ) + @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 @@ -274,14 +300,34 @@ 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), ) + @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 +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) @@ -321,7 +367,11 @@ 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])) + iis = ModelAnalyzer.constraints(ret[], model) + @test length(iis) == 2 + @test Set(iis) == Set([c2, c1]) # buf = IOBuffer() ModelAnalyzer.summarize( @@ -381,8 +431,13 @@ 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 Set([ret[].constraint[1], ret[].constraint[2]]) ⊆ Set([c3, c2, c1]) + @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 @@ -408,7 +463,11 @@ 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])) + iis = ModelAnalyzer.constraints(ret[], model) + @test length(iis) == 2 + @test Set(iis) == Set([c2, c1]) return end @@ -434,7 +493,14 @@ 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])) + 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 @@ -463,7 +529,11 @@ 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])) + iis = ModelAnalyzer.constraints(ret[], model) + @test length(iis) == 2 + @test Set(iis) == Set([c2, c1]) return end diff --git a/test/numerical.jl b/test/numerical.jl index a3c7462..3d5d6df 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, @@ -191,14 +192,32 @@ function test_constraint_bounds_quad() return end -function test_variable_not_in_constraints() +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) @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, @@ -206,33 +225,99 @@ function test_variable_not_in_constraints() @test length(ret) == 1 # 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, - ModelAnalyzer.Numerical.VariableNotInConstraints, + 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)) + # + # + 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`") + @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, - ModelAnalyzer.Numerical.VariableNotInConstraints, - verbose = false, + ret[1], + verbose = true, + model = JuMP.backend(model), ) str = String(take!(buf)) - @test str == "# VariableNotInConstraints" + @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, + ) + @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 = 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 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 @@ -241,6 +326,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) @@ -259,13 +346,27 @@ 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 = 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 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) @@ -275,6 +376,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( @@ -302,7 +405,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 @@ -311,6 +414,10 @@ function test_dense_constraint() ModelAnalyzer.Numerical.DenseConstraint, ) @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 # buf = IOBuffer() ModelAnalyzer.summarize(buf, ModelAnalyzer.Numerical.DenseConstraint) @@ -336,7 +443,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 @@ -345,6 +452,10 @@ 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.variable(ret[], JuMP.backend(model)) == JuMP.index(x) + @test ModelAnalyzer.value(ret[]) == 1e-9 # buf = IOBuffer() ModelAnalyzer.summarize(buf, ModelAnalyzer.Numerical.SmallMatrixCoefficient) @@ -372,7 +483,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 @@ -381,6 +492,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) @@ -417,6 +531,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) @@ -452,6 +568,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) @@ -478,7 +596,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 @@ -487,6 +605,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) @@ -513,7 +633,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 @@ -522,6 +642,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) @@ -558,6 +680,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( @@ -597,6 +721,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( @@ -636,6 +762,10 @@ function test_small_objective_coef_quad() ModelAnalyzer.Numerical.SmallObjectiveQuadraticCoefficient, ) @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() ModelAnalyzer.summarize( @@ -676,6 +806,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( @@ -706,7 +838,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 @@ -715,6 +847,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( @@ -745,7 +880,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 @@ -754,6 +889,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( @@ -860,7 +998,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 @@ -869,6 +1007,7 @@ function test_constraint_nonconvex() ModelAnalyzer.Numerical.NonconvexQuadraticConstraint, ) @test length(ret) == 1 + @test ModelAnalyzer.constraint(ret[], model) == c # buf = IOBuffer() ModelAnalyzer.summarize( @@ -901,7 +1040,159 @@ function test_empty_model() return end -# TODO, test SDP and empty model +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, 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) == 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 + +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]) + @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()) + @constraint(model, c6, [2 * x[1] * x[1]] in Zeros()) + data = ModelAnalyzer.analyze(ModelAnalyzer.Numerical.Analyzer(), model) + 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) + 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()