diff --git a/ext/ModelAnalyzerJuMPExt/ModelAnalyzerJuMPExt.jl b/ext/ModelAnalyzerJuMPExt/ModelAnalyzerJuMPExt.jl index 18d84c9..ce830a6 100644 --- a/ext/ModelAnalyzerJuMPExt/ModelAnalyzerJuMPExt.jl +++ b/ext/ModelAnalyzerJuMPExt/ModelAnalyzerJuMPExt.jl @@ -92,4 +92,4 @@ function ModelAnalyzer.constraints( return JuMP.constraint_ref_with_index.(model, ref) end -end # module ModelAnalyzerJuMPExt +end # module ModelAnalyzerJuMPExt diff --git a/src/Feasibility/Feasibility.jl b/src/Feasibility/Feasibility.jl new file mode 100644 index 0000000..d21c0bc --- /dev/null +++ b/src/Feasibility/Feasibility.jl @@ -0,0 +1,17 @@ +# 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 Feasibility + +import Dualization +import MathOptInterface as MOI +import ModelAnalyzer +import Printf + +include("structs.jl") +include("analyze.jl") +include("summarize.jl") + +end # module Feasibility diff --git a/src/Feasibility/analyze.jl b/src/Feasibility/analyze.jl new file mode 100644 index 0000000..cebe6c2 --- /dev/null +++ b/src/Feasibility/analyze.jl @@ -0,0 +1,446 @@ +# 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 ModelAnalyzer.analyze( + ::Analyzer, + 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 + 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.", + ) + end + data.primal_point = _last_primal_solution(model) + end + + 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.", + ) + end + data.dual_point = _last_dual_solution(model) + end + + _analyze_primal!(model, data) + _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, _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)) + return data +end + +function _analyze_primal!(model, data) + 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 _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) + 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 +end + +# not needed because it would have stoped in dualization before +# function _set_value(set::MOI.AbstractScalarSet) +# return 0.0 +# end +# function _set_value(set::MOI.Interval) +# error("Interval sets are not supported.") +# return (set.lower, set.upper) +# end + +function _set_value(set::MOI.AbstractVectorSet) + return zeros(MOI.dimension(set)) +end + +function _set_value(set::MOI.LessThan) + return set.upper +end + +function _set_value(set::MOI.GreaterThan) + return set.lower +end + +function _set_value(set::MOI.EqualTo) + return set.value +end + +function _analyze_objectives!(model::MOI.ModelLike, dual_model, map, data) + primal_status = MOI.get(model, MOI.PrimalStatus()) + dual_status = MOI.get(model, MOI.DualStatus()) + if 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 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 + + if dual_obj_val_solver !== nothing && + obj_val_solver !== nothing && + !isapprox(obj_val_solver, dual_obj_val_solver; atol = data.atol) + push!( + data.primal_dual_solver_mismatch, + PrimalDualSolverMismatch(obj_val_solver, dual_obj_val_solver), + ) + end + + 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) + push!( + data.dual_objective_mismatch, + DualObjectiveMismatch(dual_obj_val, dual_obj_val_solver), + ) + end + + 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), + ) + end + end + + return +end + +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 _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 _can_dualize(model::MOI.ModelLike) + types = MOI.get(model, MOI.ListOfConstraintTypesPresent()) + + for (F, S) in types + if !Dualization.supported_constraint(F, S) + return false + end + end + + F = MOI.get(model, MOI.ObjectiveFunctionType()) + + if !Dualization.supported_objective(F) + return false + end + + sense = MOI.get(model, MOI.ObjectiveSense()) + if sense == MOI.FEASIBILITY_SENSE + return false + end + + return true +end diff --git a/src/Feasibility/structs.jl b/src/Feasibility/structs.jl new file mode 100644 index 0000000..e8225ea --- /dev/null +++ b/src/Feasibility/structs.jl @@ -0,0 +1,252 @@ +# 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. + +""" + Analyzer() <: ModelAnalyzer.AbstractAnalyzer + +The `Analyzer` type is used to perform feasibility analysis on a model. + +## Example + +```julia +julia> data = ModelAnalyzer.analyze( + ModelAnalyzer.Feasibility.Analyzer(), + 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, +); +``` + +The additional parameters: +- `primal_point`: The primal solution point to use for feasibility checking. + If `nothing`, it will use the current primal solution from optimized model. +- `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 and dual objective + checks. +""" +struct Analyzer <: ModelAnalyzer.AbstractAnalyzer end + +""" + AbstractFeasibilityIssue <: AbstractNumericalIssue + +Abstract type for feasibility issues found during the analysis of a model. +""" +abstract type AbstractFeasibilityIssue <: ModelAnalyzer.AbstractIssue end + +""" + PrimalViolation <: AbstractFeasibilityIssue + +The `PrimalViolation` issue is identified when a primal constraint has a +left-hand-side value that is not within the constraint's set. + +For more information, run: +```julia +julia> ModelAnalyzer.summarize(ModelAnalyzer.Feasibility.PrimalViolation) +``` +""" +struct PrimalViolation <: AbstractFeasibilityIssue + ref::MOI.ConstraintIndex + violation::Float64 +end + +ModelAnalyzer.constraint(issue::PrimalViolation) = issue.ref + +ModelAnalyzer.value(issue::PrimalViolation) = issue.violation + +""" + DualConstraintViolation <: AbstractFeasibilityIssue + +The `DualConstraintViolation` issue is identified when a dual constraint has a +value that is not within the dual constraint's set. +This dual constraint corresponds to a primal variable. + +For more information, run: +```julia +julia> ModelAnalyzer.summarize(ModelAnalyzer.Feasibility.DualConstraintViolation) +``` +""" +struct DualConstraintViolation <: AbstractFeasibilityIssue + ref::MOI.VariableIndex + violation::Float64 +end + +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.DualConstrainedVariableViolation) +``` +""" +struct DualConstrainedVariableViolation <: AbstractFeasibilityIssue + ref::MOI.ConstraintIndex + violation::Float64 +end + +ModelAnalyzer.constraint(issue::DualConstrainedVariableViolation) = issue.ref + +ModelAnalyzer.value(issue::DualConstrainedVariableViolation) = issue.violation + +""" + ComplemetarityViolation <: AbstractFeasibilityIssue + +The `ComplemetarityViolation` issue is identified when a pair of primal +constraint and dual variable has a nonzero complementarity value, i.e., the +inner product of the primal constraint's slack and the dual variable's +violation is not zero. + +For more information, run: +```julia +julia> ModelAnalyzer.summarize(ModelAnalyzer.Feasibility.ComplemetarityViolation) +``` +""" +struct ComplemetarityViolation <: AbstractFeasibilityIssue + ref::MOI.ConstraintIndex + violation::Float64 +end + +ModelAnalyzer.constraint(issue::ComplemetarityViolation) = issue.ref + +ModelAnalyzer.value(issue::ComplemetarityViolation) = issue.violation + +""" + DualObjectiveMismatch <: AbstractFeasibilityIssue + +The `DualObjectiveMismatch` issue is identified when the dual objective value +computed from problem data and the dual solution does not match the solver's +dual objective value. + +For more information, run: +```julia +julia> ModelAnalyzer.summarize(ModelAnalyzer.Feasibility.DualObjectiveMismatch) +``` +""" +struct DualObjectiveMismatch <: AbstractFeasibilityIssue + obj::Float64 + obj_solver::Float64 +end + +# ModelAnalyzer.values(issue::DualObjectiveMismatch) = [issue.obj, issue.obj_solver] + +""" + PrimalObjectiveMismatch <: AbstractFeasibilityIssue + +The `PrimalObjectiveMismatch` issue is identified when the primal objective +value computed from problem data and the primal solution does not match +the solver's primal objective value. + +For more information, run: +```julia +julia> ModelAnalyzer.summarize(ModelAnalyzer.Feasibility.PrimalObjectiveMismatch) +``` +""" +struct PrimalObjectiveMismatch <: AbstractFeasibilityIssue + obj::Float64 + obj_solver::Float64 +end + +# ModelAnalyzer.values(issue::PrimalObjectiveMismatch) = [issue.obj, issue.obj_solver] + +""" + PrimalDualMismatch <: AbstractFeasibilityIssue + +The `PrimalDualMismatch` issue is identified when the primal objective value +computed from problem data and the primal solution does not match the dual +objective value computed from problem data and the dual solution. + +For more information, run: +```julia +julia> ModelAnalyzer.summarize(ModelAnalyzer.Feasibility.PrimalDualMismatch) +``` +""" +struct PrimalDualMismatch <: AbstractFeasibilityIssue + primal::Float64 + dual::Float64 +end + +ModelAnalyzer.values(issue::PrimalDualMismatch) = [issue.primal, issue.dual] + +""" + PrimalDualSolverMismatch <: AbstractFeasibilityIssue + +The `PrimalDualSolverMismatch` issue is identified when the primal objective +value reported by the solver does not match the dual objective value reported +by the solver. + +For more information, run: +```julia +julia> ModelAnalyzer.summarize(ModelAnalyzer.Feasibility.PrimalDualSolverMismatch) +``` +""" +struct PrimalDualSolverMismatch <: AbstractFeasibilityIssue + primal::Float64 + 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 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. +""" +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{DualConstraintViolation} = DualConstraintViolation[] + dual_convar::Vector{DualConstrainedVariableViolation} = + DualConstrainedVariableViolation[] + complementarity::Vector{ComplemetarityViolation} = ComplemetarityViolation[] + # objective analysis + dual_objective_mismatch::Vector{DualObjectiveMismatch} = + DualObjectiveMismatch[] + primal_objective_mismatch::Vector{PrimalObjectiveMismatch} = + PrimalObjectiveMismatch[] + primal_dual_mismatch::Vector{PrimalDualMismatch} = PrimalDualMismatch[] + primal_dual_solver_mismatch::Vector{PrimalDualSolverMismatch} = + PrimalDualSolverMismatch[] +end diff --git a/src/Feasibility/summarize.jl b/src/Feasibility/summarize.jl new file mode 100644 index 0000000..8ff6917 --- /dev/null +++ b/src/Feasibility/summarize.jl @@ -0,0 +1,630 @@ +# 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 ModelAnalyzer._summarize(io::IO, ::Type{PrimalViolation}) + return print(io, "# PrimalViolation") +end + +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}) + return print(io, "# ComplemetarityViolation") +end + +function ModelAnalyzer._summarize(io::IO, ::Type{DualObjectiveMismatch}) + return print(io, "# DualObjectiveMismatch") +end + +function ModelAnalyzer._summarize(io::IO, ::Type{PrimalObjectiveMismatch}) + return print(io, "# PrimalObjectiveMismatch") +end + +function ModelAnalyzer._summarize(io::IO, ::Type{PrimalDualMismatch}) + return print(io, "# PrimalDualMismatch") +end + +function ModelAnalyzer._summarize(io::IO, ::Type{PrimalDualSolverMismatch}) + return print(io, "# PrimalDualSolverMismatch") +end + +function ModelAnalyzer._verbose_summarize(io::IO, ::Type{PrimalViolation}) + return print( + io, + """ + # PrimalViolation + + ## What + + A `PrimalViolation` issue is identified when a constraint has + function , i.e., a left-hand-side value, that is not within + the constraint's set. + + ## Why + + This can happen due to a few reasons: + - The solver did not converge. + - The model is infeasible and the solver converged to an + infeasible point. + - The solver converged to a low accuracy solution, which might + happen due to transformations in the the model presolve or + due to numerical issues. + + ## How to fix + + Check the solver convergence log and the solver status. If the + solver did not converge, you might want to try alternative + solvers or adjust the solver options. If the solver converged + to an infeasible point, you might want to check the model + constraints and bounds. If the solver converged to a low + accuracy solution, you might want to adjust the solver options + or the model presolve. + + ## More information + + No extra information for this issue. + """, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + ::Type{DualConstraintViolation}, +) + return print( + io, + """ + # DualConstraintViolation + + ## What + + A `DualConstraintViolation` issue is identified when a constraint has + a dual value that is not within the dual constraint's set. + + ## Why + + This can happen due to a few reasons: + - The solver did not converge. + - The model is infeasible and the solver converged to an + infeasible point. + - The solver converged to a low accuracy solution, which might + happen due to transformations in the the model presolve or + due to numerical issues. + + ## How to fix + + Check the solver convergence log and the solver status. If the + solver did not converge, you might want to try alternative + solvers or adjust the solver options. If the solver converged + to an infeasible point, you might want to check the model + constraints and bounds. If the solver converged to a low + accuracy solution, you might want to adjust the solver options + or the model presolve. + + ## More information + + No extra information for this issue. + """, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + ::Type{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}, +) + return print( + io, + """ + # ComplemetarityViolation + + ## What + + A `ComplemetarityViolation` issue is identified when a pair of + primal constraint and dual varaible has a nonzero + complementarity value, i.e., the inner product of the primal + constraint's slack and the dual variable's violation is + not zero. + + ## 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{DualObjectiveMismatch}) + return print( + io, + """ + # DualObjectiveMismatch + + ## What + + A `DualObjectiveMismatch` issue is identified when the dual + objective value computed from problema data and the dual + solution does not match the solver's dual objective + value. + + ## Why + + This can happen due to: + - The solver performed presolve transformations and the + reported dual objective is reported from the transformed + problem. + - Bad problem numerical conditioning, very large and very + small coefficients might be present in the model. + + ## How to fix + + Check the solver convergence log and the solver status. + Consider reviewing the coefficients of the objective function. + Consider reviewing the options set in the solver. + + ## More information + + No extra information for this issue. + """, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + ::Type{PrimalObjectiveMismatch}, +) + return print( + io, + """ + # PrimalObjectiveMismatch + + ## What + + A `PrimalObjectiveMismatch` issue is identified when the primal + objective value computed from problema data and the primal + solution does not match the solver's primal objective + value. + + ## Why + + This can happen due to: + - The solver performed presolve transformations and the + reported primal objective is reported from the transformed + problem. + - Bad problem numerical conditioning, very large and very + small coefficients might be present in the model. + + ## How to fix + + Check the solver convergence log and the solver status. + Consider reviewing the coefficients of the objective function. + Consider reviewing the options set in the solver. + + ## More information + + No extra information for this issue. + """, + ) +end + +function ModelAnalyzer._verbose_summarize(io::IO, ::Type{PrimalDualMismatch}) + return print( + io, + """ + # PrimalDualMismatch + + ## What + + A `PrimalDualMismatch` issue is identified when the primal + objective value computed from problema data and the primal + solution does not match the dual objective value computed + from problem data and the dual solution. + + ## Why + + This can happen due to: + - The solver did not converge. + - Bad problem numerical conditioning, very large and very + small coefficients might be present in the model. + + ## How to fix + + Check the solver convergence log and the solver status. + Consider reviewing the coefficients of the model. + + ## More information + + No extra information for this issue. + """, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + ::Type{PrimalDualSolverMismatch}, +) + return print( + io, + """ + # PrimalDualSolverMismatch + + ## What + + A `PrimalDualSolverMismatch` issue is identified when the primal + objective value reported by the solver does not match the dual + objective value reported by the solver. + + ## Why + + This can happen due to: + - The solver did not converge. + + ## How to fix + + Check the solver convergence log and the solver status. + + ## More information + + No extra information for this issue. + """, + ) +end + +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::DualConstrainedVariableViolation, + model, +) + return print( + io, + ModelAnalyzer._name(issue.ref, model), + " : ", + issue.violation, + ) +end + +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, model) + return ModelAnalyzer._verbose_summarize(io, issue, model) +end + +function ModelAnalyzer._summarize(io::IO, issue::PrimalObjectiveMismatch, model) + return ModelAnalyzer._verbose_summarize(io, issue, model) +end + +function ModelAnalyzer._summarize(io::IO, issue::PrimalDualMismatch, model) + return ModelAnalyzer._verbose_summarize(io, issue, model) +end + +function ModelAnalyzer._summarize( + io::IO, + issue::PrimalDualSolverMismatch, + model, +) + return ModelAnalyzer._verbose_summarize(io, issue, model) +end + +function ModelAnalyzer._verbose_summarize(io::IO, issue::PrimalViolation, model) + return print( + io, + "Constraint ", + ModelAnalyzer._name(issue.ref, model), + " has primal 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 ", + ModelAnalyzer._name(issue.ref, model), + " has complementarty violation ", + issue.violation, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::DualObjectiveMismatch, + model, +) + return print( + io, + "Dual objective mismatch: ", + issue.obj, + " (computed) vs ", + issue.obj_solver, + " (reported by solver)\n", + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::PrimalObjectiveMismatch, + model, +) + return print( + io, + "Primal objective mismatch: ", + issue.obj, + " (computed) vs ", + issue.obj_solver, + " (reported by solver)\n", + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::PrimalDualMismatch, + model, +) + return print( + io, + "Primal dual mismatch: ", + issue.primal, + " (computed primal) vs ", + issue.dual, + " (computed dual)\n", + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::PrimalDualSolverMismatch, + model, +) + return print( + io, + "Solver reported objective mismatch: ", + issue.primal, + " (reported primal) vs ", + issue.dual, + " (reported dual)\n", + ) +end + +function ModelAnalyzer.list_of_issues(data::Data, ::Type{PrimalViolation}) + return data.primal +end + +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}, +) + return data.complementarity +end + +function ModelAnalyzer.list_of_issues(data::Data, ::Type{DualObjectiveMismatch}) + return data.dual_objective_mismatch +end + +function ModelAnalyzer.list_of_issues( + data::Data, + ::Type{PrimalObjectiveMismatch}, +) + return data.primal_objective_mismatch +end + +function ModelAnalyzer.list_of_issues(data::Data, ::Type{PrimalDualMismatch}) + return data.primal_dual_mismatch +end + +function ModelAnalyzer.list_of_issues( + data::Data, + ::Type{PrimalDualSolverMismatch}, +) + return data.primal_dual_solver_mismatch +end + +function ModelAnalyzer.list_of_issue_types(data::Data) + ret = Type[] + for type in ( + PrimalViolation, + DualConstraintViolation, + DualConstrainedVariableViolation, + ComplemetarityViolation, + DualObjectiveMismatch, + PrimalObjectiveMismatch, + PrimalDualMismatch, + PrimalDualSolverMismatch, + ) + if !isempty(ModelAnalyzer.list_of_issues(data, type)) + push!(ret, type) + end + end + return ret +end + +function summarize_configurations(io::IO, data::Data) + print(io, "## Configuration\n\n") + # print(io, " - point: ", data.point, "\n") + print(io, " atol: ", data.atol, "\n") + print(io, " skip_missing: ", data.skip_missing, "\n") + return +end + +function ModelAnalyzer.summarize( + io::IO, + data::Data; + model = nothing, + verbose = true, + max_issues = ModelAnalyzer.DEFAULT_MAX_ISSUES, + configurations = true, +) + print(io, "## Feasibility Analysis\n\n") + if configurations + summarize_configurations(io, data) + print(io, "\n") + end + # add maximum primal, dual and compl + # add sum of primal, dual and compl + for issue_type in ModelAnalyzer.list_of_issue_types(data) + issues = ModelAnalyzer.list_of_issues(data, issue_type) + print(io, "\n\n") + ModelAnalyzer.summarize( + io, + issues, + model = model, + verbose = verbose, + max_issues = max_issues, + ) + end + return +end + +function Base.show(io::IO, data::Data) + n = sum( + length(ModelAnalyzer.list_of_issues(data, T)) for + T in ModelAnalyzer.list_of_issue_types(data); + init = 0, + ) + return print(io, "Feasibility analysis found $n issues") +end diff --git a/src/Infeasibility/Infeasibility.jl b/src/Infeasibility/Infeasibility.jl new file mode 100644 index 0000000..d3e2f7d --- /dev/null +++ b/src/Infeasibility/Infeasibility.jl @@ -0,0 +1,18 @@ +# 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 Infeasibility + +import MathOptInterface as MOI +import ModelAnalyzer + +include("intervals.jl") +include("_eval_variables.jl") + +include("structs.jl") +include("analyze.jl") +include("summarize.jl") + +end # module Infeasibility diff --git a/src/_eval_variables.jl b/src/Infeasibility/_eval_variables.jl similarity index 100% rename from src/_eval_variables.jl rename to src/Infeasibility/_eval_variables.jl diff --git a/src/infeasibility.jl b/src/Infeasibility/analyze.jl similarity index 55% rename from src/infeasibility.jl rename to src/Infeasibility/analyze.jl index a93edf6..7c94ddb 100644 --- a/src/infeasibility.jl +++ b/src/Infeasibility/analyze.jl @@ -3,161 +3,6 @@ # 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 Infeasibility - -import ModelAnalyzer -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 model. - -## Example -```julia -julia> data = ModelAnalyzer.analyze( - Analyzer(), - model, - optimizer = nothing,, -) -``` - -The additional keyword argument `optimizer` is used to specify the optimizer to -use for the IIS resolver. -""" -struct Analyzer <: ModelAnalyzer.AbstractAnalyzer end - -""" - AbstractInfeasibilitylIssue - -Abstract type for infeasibility issues found during the analysis of a -model. -""" -abstract type AbstractInfeasibilitylIssue <: ModelAnalyzer.AbstractIssue end - -""" - InfeasibleBounds{T} <: AbstractInfeasibilitylIssue - -The `InfeasibleBounds` issue is identified when a variable has a lower bound -that is greater than its upper bound. - -For more information, run: -```julia -julia> ModelAnalyzer.summarize(ModelAnalyzer.Infeasibility.InfeasibleBounds) -```` -""" -struct InfeasibleBounds{T} <: AbstractInfeasibilitylIssue - 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 - -The `InfeasibleIntegrality` issue is identified when a variable has an -integrality constraint (like `MOI.Integer` or `MOI.ZeroOne`) that is not -consistent with its bounds. That is, the bounds do not allow for any -integer value to be feasible. - -For more information, run: -```julia -julia> ModelAnalyzer.summarize( - ModelAnalyzer.Infeasibility.InfeasibleIntegrality -) -``` -""" -struct InfeasibleIntegrality{T} <: AbstractInfeasibilitylIssue - 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 - -The `InfeasibleConstraintRange` issue is identified when a constraint cannot -be satisfied given the variable bounds. This analysis only considers one -constraint at a time and all variable bounds of variables involved in the -constraint. -This issue can only be found is all variable bounds are consistent, that is, -no issues of type `InfeasibleBounds` were found in the first layer of analysis. - -For more information, run: -```julia -julia> ModelAnalyzer.summarize( - ModelAnalyzer.Infeasibility.InfeasibleConstraintRange -) -``` -""" -struct InfeasibleConstraintRange{T} <: AbstractInfeasibilitylIssue - 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 - -The `IrreducibleInfeasibleSubset` issue is identified when a subset of -constraints cannot be satisfied simultaneously. This is typically found -by the IIS resolver after the first two layers of infeasibility analysis -have been completed with no issues, that is, no issues of any other type -were found. - -For more information, run: -```julia -julia> ModelAnalyzer.summarize( - ModelAnalyzer.Infeasibility.IrreducibleInfeasibleSubset -) -``` -""" -struct IrreducibleInfeasibleSubset <: AbstractInfeasibilitylIssue - constraint::Vector{<:MOI.ConstraintIndex} -end - -ModelAnalyzer.constraints(issue::IrreducibleInfeasibleSubset) = issue.constraint - -""" - Data <: ModelAnalyzer.AbstractData - -The `Data` type is used to store the results of the infeasibility analysis. -This type contains vectors of the various infeasibility issues found during -the analysis, including `InfeasibleBounds`, `InfeasibleIntegrality`, -`InfeasibleConstraintRange`, and `IrreducibleInfeasibleSubset`. -""" -Base.@kwdef mutable struct Data <: ModelAnalyzer.AbstractData - infeasible_bounds::Vector{InfeasibleBounds} = InfeasibleBounds[] - infeasible_integrality::Vector{InfeasibleIntegrality} = - InfeasibleIntegrality[] - - constraint_range::Vector{InfeasibleConstraintRange} = - InfeasibleConstraintRange[] - - iis::Vector{IrreducibleInfeasibleSubset} = IrreducibleInfeasibleSubset[] -end - function ModelAnalyzer.analyze( ::Analyzer, model::MOI.ModelLike; @@ -604,328 +449,3 @@ function iis_elastic_filter(original_model::MOI.ModelLike, optimizer) return iis end - -# API - -function ModelAnalyzer._summarize(io::IO, ::Type{<:InfeasibleBounds}) - return print(io, "# InfeasibleBounds") -end - -function ModelAnalyzer._summarize(io::IO, ::Type{<:InfeasibleIntegrality}) - return print(io, "# InfeasibleIntegrality") -end - -function ModelAnalyzer._summarize(io::IO, ::Type{<:InfeasibleConstraintRange}) - return print(io, "# InfeasibleConstraintRange") -end - -function ModelAnalyzer._summarize(io::IO, ::Type{<:IrreducibleInfeasibleSubset}) - return print(io, "# IrreducibleInfeasibleSubset") -end - -function ModelAnalyzer._verbose_summarize(io::IO, ::Type{<:InfeasibleBounds}) - return print( - io, - """ - # `InfeasibleBounds` - - ## What - - A `InfeasibleBounds` issue is identified when a variable has an - lower bound that is greater than its upper bound. - - ## Why - - This can be a sign of a mistake in the model formulation. This error - will lead to infeasibility in the optimization problem. - - ## How to fix - - Fix one of both of the bounds. - - ## More information - - No extra information for this issue. - """, - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - ::Type{<:InfeasibleIntegrality}, -) - return print( - io, - """ - # `InfeasibleIntegrality` - - ## What - - A `InfeasibleIntegrality` issue is identified when a variable has an - and integrality constraint and the bounds do not allow for any integer - value to be feasible. - - ## Why - - This can be a sign of a mistake in the model formulation. This error - will lead to infeasibility in the optimization problem. - - ## How to fix - - Fix one of both of the bounds or remove the integrality constraint. - - ## More information - - No extra information for this issue. - """, - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - ::Type{<:InfeasibleConstraintRange}, -) - return print( - io, - """ - # `InfeasibleConstraintRange` - - ## What - - A `InfeasibleConstraintRange` issue is identified when given the variable bounds - a constraint cannot be satisfied. This analysis only considers one contraint at - a time and all variable bounds of variables involved in the constraint. - - ## Why - - This can be a sign of a mistake in the model formulation. This error - will lead to infeasibility in the optimization problem. - - ## How to fix - - Fix the bounds of variables or the constraint. - - ## More information - - No extra information for this issue. - """, - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - ::Type{<:IrreducibleInfeasibleSubset}, -) - return print( - io, - """ - # `IrreducibleInfeasibleSubset` - - ## What - - An `IrreducibleInfeasibleSubset` issue is identified when a subset of constraints - cannot be satisfied simultaneously. - - ## Why - - This can be a sign of a mistake in the model formulation. This error - will lead to infeasibility in the optimization problem. - - ## How to fix - - Fix the constraints in question. - - ## More information - - No extra information for this issue. - """, - ) -end - -function ModelAnalyzer._summarize(io::IO, issue::InfeasibleBounds, model) - return print( - io, - ModelAnalyzer._name(issue.variable, model), - " : ", - issue.lb, - " !<= ", - issue.ub, - ) -end - -function ModelAnalyzer._summarize(io::IO, issue::InfeasibleIntegrality, model) - return print( - io, - ModelAnalyzer._name(issue.variable, model), - " : [", - issue.lb, - "; ", - issue.ub, - "], ", - issue.set, - ) -end - -function ModelAnalyzer._summarize( - io::IO, - issue::InfeasibleConstraintRange, - model, -) - return print( - io, - ModelAnalyzer._name(issue.constraint, model), - " : [", - issue.lb, - "; ", - issue.ub, - "], !in ", - issue.set, - ) -end - -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, - model, -) - return print( - io, - "Variable: ", - ModelAnalyzer._name(issue.variable, model), - " with lower bound ", - issue.lb, - " and upper bound ", - issue.ub, - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - issue::InfeasibleIntegrality, - model, -) - return print( - io, - "Variable: ", - ModelAnalyzer._name(issue.variable, model), - " with lower bound ", - issue.lb, - " and upper bound ", - issue.ub, - " and integrality constraint: ", - issue.set, - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - issue::InfeasibleConstraintRange, - model, -) - return print( - io, - "Constraint: ", - ModelAnalyzer._name(issue.constraint, model), - " with computed lower bound ", - issue.lb, - " and computed upper bound ", - issue.ub, - " and set: ", - issue.set, - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - issue::IrreducibleInfeasibleSubset, - model, -) - return print( - io, - "Irreducible Infeasible Subset: ", - join(map(x -> ModelAnalyzer._name(x, model), issue.constraint), ", "), - ) -end - -function ModelAnalyzer.list_of_issues(data::Data, ::Type{InfeasibleBounds}) - return data.infeasible_bounds -end - -function ModelAnalyzer.list_of_issues(data::Data, ::Type{InfeasibleIntegrality}) - return data.infeasible_integrality -end - -function ModelAnalyzer.list_of_issues( - data::Data, - ::Type{InfeasibleConstraintRange}, -) - return data.constraint_range -end - -function ModelAnalyzer.list_of_issues( - data::Data, - ::Type{IrreducibleInfeasibleSubset}, -) - return data.iis -end - -function ModelAnalyzer.list_of_issue_types(data::Data) - ret = Type[] - for type in ( - InfeasibleBounds, - InfeasibleIntegrality, - InfeasibleConstraintRange, - IrreducibleInfeasibleSubset, - ) - if !isempty(ModelAnalyzer.list_of_issues(data, type)) - push!(ret, type) - end - end - return ret -end - -function ModelAnalyzer.summarize( - io::IO, - data::Data; - model = nothing, - verbose = true, - max_issues = ModelAnalyzer.DEFAULT_MAX_ISSUES, -) - print(io, "## Infeasibility Analysis\n\n") - - for issue_type in ModelAnalyzer.list_of_issue_types(data) - issues = ModelAnalyzer.list_of_issues(data, issue_type) - print(io, "\n\n") - ModelAnalyzer.summarize( - io, - issues, - model = model, - verbose = verbose, - max_issues = max_issues, - ) - end - return -end - -function Base.show(io::IO, data::Data) - n = sum( - length(ModelAnalyzer.list_of_issues(data, T)) for - T in ModelAnalyzer.list_of_issue_types(data); - init = 0, - ) - return print(io, "Infeasibility analysis found $n issues") -end - -end # module diff --git a/src/intervals.jl b/src/Infeasibility/intervals.jl similarity index 100% rename from src/intervals.jl rename to src/Infeasibility/intervals.jl diff --git a/src/Infeasibility/structs.jl b/src/Infeasibility/structs.jl new file mode 100644 index 0000000..f24a2c3 --- /dev/null +++ b/src/Infeasibility/structs.jl @@ -0,0 +1,151 @@ +# 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. + +""" + Analyzer() <: ModelAnalyzer.AbstractAnalyzer + +The `Analyzer` type is used to perform infeasibility analysis on a model. + +## Example +```julia +julia> data = ModelAnalyzer.analyze( + Analyzer(), + model, + optimizer = nothing,, +) +``` + +The additional keyword argument `optimizer` is used to specify the optimizer to +use for the IIS resolver. +""" +struct Analyzer <: ModelAnalyzer.AbstractAnalyzer end + +""" + AbstractInfeasibilitylIssue + +Abstract type for infeasibility issues found during the analysis of a +model. +""" +abstract type AbstractInfeasibilitylIssue <: ModelAnalyzer.AbstractIssue end + +""" + InfeasibleBounds{T} <: AbstractInfeasibilitylIssue + +The `InfeasibleBounds` issue is identified when a variable has a lower bound +that is greater than its upper bound. + +For more information, run: +```julia +julia> ModelAnalyzer.summarize(ModelAnalyzer.Infeasibility.InfeasibleBounds) +```` +""" +struct InfeasibleBounds{T} <: AbstractInfeasibilitylIssue + 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 + +The `InfeasibleIntegrality` issue is identified when a variable has an +integrality constraint (like `MOI.Integer` or `MOI.ZeroOne`) that is not +consistent with its bounds. That is, the bounds do not allow for any +integer value to be feasible. + +For more information, run: +```julia +julia> ModelAnalyzer.summarize( + ModelAnalyzer.Infeasibility.InfeasibleIntegrality +) +``` +""" +struct InfeasibleIntegrality{T} <: AbstractInfeasibilitylIssue + 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 + +The `InfeasibleConstraintRange` issue is identified when a constraint cannot +be satisfied given the variable bounds. This analysis only considers one +constraint at a time and all variable bounds of variables involved in the +constraint. +This issue can only be found is all variable bounds are consistent, that is, +no issues of type `InfeasibleBounds` were found in the first layer of analysis. + +For more information, run: +```julia +julia> ModelAnalyzer.summarize( + ModelAnalyzer.Infeasibility.InfeasibleConstraintRange +) +``` +""" +struct InfeasibleConstraintRange{T} <: AbstractInfeasibilitylIssue + 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 + +The `IrreducibleInfeasibleSubset` issue is identified when a subset of +constraints cannot be satisfied simultaneously. This is typically found +by the IIS resolver after the first two layers of infeasibility analysis +have been completed with no issues, that is, no issues of any other type +were found. + +For more information, run: +```julia +julia> ModelAnalyzer.summarize( + ModelAnalyzer.Infeasibility.IrreducibleInfeasibleSubset +) +``` +""" +struct IrreducibleInfeasibleSubset <: AbstractInfeasibilitylIssue + constraint::Vector{<:MOI.ConstraintIndex} +end + +ModelAnalyzer.constraints(issue::IrreducibleInfeasibleSubset) = issue.constraint + +""" + Data <: ModelAnalyzer.AbstractData + +The `Data` type is used to store the results of the infeasibility analysis. +This type contains vectors of the various infeasibility issues found during +the analysis, including `InfeasibleBounds`, `InfeasibleIntegrality`, +`InfeasibleConstraintRange`, and `IrreducibleInfeasibleSubset`. +""" +Base.@kwdef mutable struct Data <: ModelAnalyzer.AbstractData + infeasible_bounds::Vector{InfeasibleBounds} = InfeasibleBounds[] + infeasible_integrality::Vector{InfeasibleIntegrality} = + InfeasibleIntegrality[] + + constraint_range::Vector{InfeasibleConstraintRange} = + InfeasibleConstraintRange[] + + iis::Vector{IrreducibleInfeasibleSubset} = IrreducibleInfeasibleSubset[] +end diff --git a/src/Infeasibility/summarize.jl b/src/Infeasibility/summarize.jl new file mode 100644 index 0000000..54cd369 --- /dev/null +++ b/src/Infeasibility/summarize.jl @@ -0,0 +1,325 @@ +# 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 ModelAnalyzer._summarize(io::IO, ::Type{<:InfeasibleBounds}) + return print(io, "# InfeasibleBounds") +end + +function ModelAnalyzer._summarize(io::IO, ::Type{<:InfeasibleIntegrality}) + return print(io, "# InfeasibleIntegrality") +end + +function ModelAnalyzer._summarize(io::IO, ::Type{<:InfeasibleConstraintRange}) + return print(io, "# InfeasibleConstraintRange") +end + +function ModelAnalyzer._summarize(io::IO, ::Type{<:IrreducibleInfeasibleSubset}) + return print(io, "# IrreducibleInfeasibleSubset") +end + +function ModelAnalyzer._verbose_summarize(io::IO, ::Type{<:InfeasibleBounds}) + return print( + io, + """ + # `InfeasibleBounds` + + ## What + + A `InfeasibleBounds` issue is identified when a variable has an + lower bound that is greater than its upper bound. + + ## Why + + This can be a sign of a mistake in the model formulation. This error + will lead to infeasibility in the optimization problem. + + ## How to fix + + Fix one of both of the bounds. + + ## More information + + No extra information for this issue. + """, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + ::Type{<:InfeasibleIntegrality}, +) + return print( + io, + """ + # `InfeasibleIntegrality` + + ## What + + A `InfeasibleIntegrality` issue is identified when a variable has an + and integrality constraint and the bounds do not allow for any integer + value to be feasible. + + ## Why + + This can be a sign of a mistake in the model formulation. This error + will lead to infeasibility in the optimization problem. + + ## How to fix + + Fix one of both of the bounds or remove the integrality constraint. + + ## More information + + No extra information for this issue. + """, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + ::Type{<:InfeasibleConstraintRange}, +) + return print( + io, + """ + # `InfeasibleConstraintRange` + + ## What + + A `InfeasibleConstraintRange` issue is identified when given the variable bounds + a constraint cannot be satisfied. This analysis only considers one contraint at + a time and all variable bounds of variables involved in the constraint. + + ## Why + + This can be a sign of a mistake in the model formulation. This error + will lead to infeasibility in the optimization problem. + + ## How to fix + + Fix the bounds of variables or the constraint. + + ## More information + + No extra information for this issue. + """, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + ::Type{<:IrreducibleInfeasibleSubset}, +) + return print( + io, + """ + # `IrreducibleInfeasibleSubset` + + ## What + + An `IrreducibleInfeasibleSubset` issue is identified when a subset of constraints + cannot be satisfied simultaneously. + + ## Why + + This can be a sign of a mistake in the model formulation. This error + will lead to infeasibility in the optimization problem. + + ## How to fix + + Fix the constraints in question. + + ## More information + + No extra information for this issue. + """, + ) +end + +function ModelAnalyzer._summarize(io::IO, issue::InfeasibleBounds, model) + return print( + io, + ModelAnalyzer._name(issue.variable, model), + " : ", + issue.lb, + " !<= ", + issue.ub, + ) +end + +function ModelAnalyzer._summarize(io::IO, issue::InfeasibleIntegrality, model) + return print( + io, + ModelAnalyzer._name(issue.variable, model), + " : [", + issue.lb, + "; ", + issue.ub, + "], ", + issue.set, + ) +end + +function ModelAnalyzer._summarize( + io::IO, + issue::InfeasibleConstraintRange, + model, +) + return print( + io, + ModelAnalyzer._name(issue.constraint, model), + " : [", + issue.lb, + "; ", + issue.ub, + "], !in ", + issue.set, + ) +end + +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, + model, +) + return print( + io, + "Variable: ", + ModelAnalyzer._name(issue.variable, model), + " with lower bound ", + issue.lb, + " and upper bound ", + issue.ub, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::InfeasibleIntegrality, + model, +) + return print( + io, + "Variable: ", + ModelAnalyzer._name(issue.variable, model), + " with lower bound ", + issue.lb, + " and upper bound ", + issue.ub, + " and integrality constraint: ", + issue.set, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::InfeasibleConstraintRange, + model, +) + return print( + io, + "Constraint: ", + ModelAnalyzer._name(issue.constraint, model), + " with computed lower bound ", + issue.lb, + " and computed upper bound ", + issue.ub, + " and set: ", + issue.set, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::IrreducibleInfeasibleSubset, + model, +) + return print( + io, + "Irreducible Infeasible Subset: ", + join(map(x -> ModelAnalyzer._name(x, model), issue.constraint), ", "), + ) +end + +function ModelAnalyzer.list_of_issues(data::Data, ::Type{InfeasibleBounds}) + return data.infeasible_bounds +end + +function ModelAnalyzer.list_of_issues(data::Data, ::Type{InfeasibleIntegrality}) + return data.infeasible_integrality +end + +function ModelAnalyzer.list_of_issues( + data::Data, + ::Type{InfeasibleConstraintRange}, +) + return data.constraint_range +end + +function ModelAnalyzer.list_of_issues( + data::Data, + ::Type{IrreducibleInfeasibleSubset}, +) + return data.iis +end + +function ModelAnalyzer.list_of_issue_types(data::Data) + ret = Type[] + for type in ( + InfeasibleBounds, + InfeasibleIntegrality, + InfeasibleConstraintRange, + IrreducibleInfeasibleSubset, + ) + if !isempty(ModelAnalyzer.list_of_issues(data, type)) + push!(ret, type) + end + end + return ret +end + +function ModelAnalyzer.summarize( + io::IO, + data::Data; + model = nothing, + verbose = true, + max_issues = ModelAnalyzer.DEFAULT_MAX_ISSUES, +) + print(io, "## Infeasibility Analysis\n\n") + + for issue_type in ModelAnalyzer.list_of_issue_types(data) + issues = ModelAnalyzer.list_of_issues(data, issue_type) + print(io, "\n\n") + ModelAnalyzer.summarize( + io, + issues, + model = model, + verbose = verbose, + max_issues = max_issues, + ) + end + return +end + +function Base.show(io::IO, data::Data) + n = sum( + length(ModelAnalyzer.list_of_issues(data, T)) for + T in ModelAnalyzer.list_of_issue_types(data); + init = 0, + ) + return print(io, "Infeasibility analysis found $n issues") +end diff --git a/src/ModelAnalyzer.jl b/src/ModelAnalyzer.jl index c913b18..55eb526 100644 --- a/src/ModelAnalyzer.jl +++ b/src/ModelAnalyzer.jl @@ -219,8 +219,8 @@ function _name(ref, ::Nothing) return "$ref" end -include("numerical.jl") -include("feasibility.jl") -include("infeasibility.jl") +include("Numerical/Numerical.jl") +include("Feasibility/Feasibility.jl") +include("Infeasibility/Infeasibility.jl") -end # module ModelAnalyzer +end # module ModelAnalyzer diff --git a/src/Numerical/Numerical.jl b/src/Numerical/Numerical.jl new file mode 100644 index 0000000..6ac45cc --- /dev/null +++ b/src/Numerical/Numerical.jl @@ -0,0 +1,17 @@ +# 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 Numerical + +import LinearAlgebra +import MathOptInterface as MOI +import ModelAnalyzer +import Printf + +include("structs.jl") +include("analyze.jl") +include("summarize.jl") + +end # module Numerical diff --git a/src/Numerical/analyze.jl b/src/Numerical/analyze.jl new file mode 100644 index 0000000..ae2c539 --- /dev/null +++ b/src/Numerical/analyze.jl @@ -0,0 +1,689 @@ +# 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. + +""" + 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 model. + +""" +function ModelAnalyzer.analyze( + ::Analyzer, + model::MOI.ModelLike, + ; + threshold_dense_fill_in::Float64 = 0.10, + threshold_dense_entries::Int = 1000, + threshold_small::Float64 = 1e-5, + threshold_large::Float64 = 1e+5, +) + data = Data() + data.threshold_dense_fill_in = threshold_dense_fill_in + data.threshold_dense_entries = threshold_dense_entries + data.threshold_small = threshold_small + data.threshold_large = threshold_large + + # initialize simples data + 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) + + # objective pass + objective_type = MOI.get(model, MOI.ObjectiveFunctionType()) + obj_func = MOI.get(model, MOI.ObjectiveFunction{objective_type}()) + _get_objective_data(data, obj_func) + + # constraints pass + 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 + 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 + # 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, + VariableNotInConstraints(var), + ) + end + end + sort!(data.dense_rows, by = x -> x.nnz, rev = true) + sort!(data.matrix_small, by = x -> abs(x.coefficient)) + sort!(data.matrix_large, by = x -> abs(x.coefficient), rev = true) + sort!(data.bounds_small, by = x -> abs(x.coefficient)) + sort!(data.bounds_large, by = x -> abs(x.coefficient), rev = true) + sort!(data.rhs_small, by = x -> abs(x.coefficient)) + sort!(data.rhs_large, by = x -> abs(x.coefficient), rev = true) + sort!(data.matrix_quadratic_small, by = x -> abs(x.coefficient)) + sort!(data.matrix_quadratic_large, by = x -> abs(x.coefficient), rev = true) + # objective + sort!(data.objective_small, by = x -> abs(x.coefficient)) + sort!(data.objective_large, by = x -> abs(x.coefficient), rev = true) + sort!(data.objective_quadratic_small, by = x -> abs(x.coefficient)) + sort!( + data.objective_quadratic_large, + by = x -> abs(x.coefficient), + rev = true, + ) + return data +end + +function _update_range(range::Vector{Float64}, value::Number) + range[1] = min(range[1], abs(value)) + range[2] = max(range[2], abs(value)) + return 1 +end + +function _get_objective_data(data, func::MOI.VariableIndex) + return +end + +function _get_objective_data(data, func::MOI.ScalarAffineFunction) + nnz = 0 + for term in func.terms + variable = term.variable + coefficient = term.coefficient + if iszero(coefficient) + continue + end + nnz += _update_range(data.objective_range, coefficient) + if abs(coefficient) < data.threshold_small + push!( + data.objective_small, + SmallObjectiveCoefficient(variable, coefficient), + ) + elseif abs(coefficient) > data.threshold_large + push!( + data.objective_large, + LargeObjectiveCoefficient(variable, coefficient), + ) + end + end + return +end + +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 term in func.quadratic_terms + coefficient = term.coefficient + v1 = term.variable_1 + v2 = term.variable_2 + if iszero(coefficient) + continue + end + nnz += _update_range(data.objective_quadratic_range, coefficient) + if abs(coefficient) < data.threshold_small + push!( + data.objective_quadratic_small, + SmallObjectiveQuadraticCoefficient(v1, v2, coefficient), + ) + elseif abs(coefficient) > data.threshold_large + push!( + data.objective_quadratic_large, + LargeObjectiveQuadraticCoefficient(v1, v2, coefficient), + ) + end + end + data.has_quadratic_objective = true + if data.sense == MOI.MAX_SENSE + if !_quadratic_vexity(func, -1) + push!(data.nonconvex_objective, NonconvexQuadraticObjective()) + end + elseif data.sense == MOI.MIN_SENSE + if !_quadratic_vexity(func, 1) + push!(data.nonconvex_objective, NonconvexQuadraticObjective()) + end + end + return +end + +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{MOI.VariableIndex,Int}() + for (idx, var) in enumerate(variables) + var_map[var] = idx + end + matrix = zeros(length(variables), length(variables)) + 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), + LinearAlgebra.RowMaximum(), + check = false, + ) + 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::MOI.ConstraintIndex, + func::MOI.ScalarAffineFunction; + ignore_extras = false, +) + if length(func.terms) == 1 + 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 + return + end + end + nnz = 0 + for term in func.terms + variable = term.variable + coefficient = term.coefficient + if iszero(coefficient) + continue + end + nnz += _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 + if nnz == 0 + if !ignore_extras + push!(data.empty_rows, EmptyConstraint(ref)) + end + return + end + if nnz / data.number_of_variables > data.threshold_dense_fill_in && + nnz > data.threshold_dense_entries + push!(data.dense_rows, DenseConstraint(ref, nnz)) + end + data.matrix_nnz += nnz + return +end + +function _get_constraint_matrix_data( + data, + ref::MOI.ConstraintIndex, + func::MOI.ScalarQuadraticFunction{T}, +) where {T} + nnz = 0 + for term in func.quadratic_terms + v1 = term.variable_1 + v2 = term.variable_2 + coefficient = term.coefficient + if iszero(coefficient) + continue + end + nnz += _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 + data.has_quadratic_constraints = true + _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::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 +end + +function _get_constraint_matrix_data( + data, + 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 + _get_constraint_matrix_data( + data, + ref, + MOI.VectorAffineFunction{T}(func.affine_terms, func.constants), + # ignore_extras = nnz > 0, + ) + return +end + +function _get_constraint_matrix_data( + data, + ref::MOI.ConstraintIndex, + func::MOI.VariableIndex, +) + # push!(data.variables_in_constraints, func) + return +end + +function _get_constraint_matrix_data( + data, + ref::MOI.ConstraintIndex, + func::MOI.VectorOfVariables, +) + if length(func.variables) == 1 + return + end + for var in func.variables + push!(data.variables_in_constraints, var) + end + return +end + +function _get_constraint_data( + data, + ref, + func::Union{MOI.ScalarAffineFunction,MOI.ScalarQuadraticFunction}, + set, +) + coefficient = func.constant + if iszero(coefficient) + return + 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 + return +end + +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::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 + return +end + +function _get_constraint_data( + data, + ref, + func::MOI.ScalarAffineFunction, + set::MOI.LessThan, +) + coefficient = set.upper - func.constant + if iszero(coefficient) + return + 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 + return +end + +function _get_constraint_data( + data, + ref, + 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 + return +end + +function _get_constraint_data( + data, + ref, + func::MOI.ScalarAffineFunction, + set::MOI.GreaterThan, +) + coefficient = set.lower - func.constant + if iszero(coefficient) + return + 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 + return +end + +function _get_constraint_data( + data, + ref, + func::MOI.ScalarQuadraticFunction, + set::Union{MOI.EqualTo,MOI.Interval}, +) + _get_constraint_data( + data, + ref, + MOI.ScalarAffineFunction(func.affine_terms, func.constant), + set, + ) + push!(data.nonconvex_rows, NonconvexQuadraticConstraint(ref)) + return +end + +function _get_constraint_data( + data, + ref, + func::MOI.VectorQuadraticFunction, + set::MOI.Zeros, +) + _get_constraint_data( + data, + ref, + MOI.VectorAffineFunction(func.affine_terms, func.constants), + set, + ) + push!(data.nonconvex_rows, NonconvexQuadraticConstraint(ref)) + return +end + +function _get_constraint_data( + data, + ref, + func::MOI.ScalarAffineFunction, + set::MOI.EqualTo, +) + coefficient = set.value - func.constant + if iszero(coefficient) + return + 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 + return +end + +function _get_constraint_data( + data, + ref, + func::MOI.ScalarAffineFunction, + set::MOI.Interval, +) + coefficient = set.upper - func.constant + if !(iszero(coefficient)) + _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 + coefficient = set.lower - func.constant + if iszero(coefficient) + return + 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 + return +end + +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 diff --git a/src/Numerical/structs.jl b/src/Numerical/structs.jl new file mode 100644 index 0000000..a6e3882 --- /dev/null +++ b/src/Numerical/structs.jl @@ -0,0 +1,494 @@ +# 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. + +""" + Analyzer() <: ModelAnalyzer.AbstractAnalyzer + +The `Analyzer` type is used to analyze the coefficients of a model for +numerical issues. + +## Example + +```julia +julia> data = ModelAnalyzer.analyze( + ModelAnalyzer.Numerical.Analyzer(), + model; + threshold_dense_fill_in = 0.10, + threshold_dense_entries = 1000, + threshold_small = 1e-5, + threshold_large = 1e+5, +) +``` + +The additional parameters: +- `threshold_dense_fill_in`: The threshold for the fraction of non-zero entries + in a constraint to be considered dense. +- `threshold_dense_entries`: The minimum number of non-zero entries for a + constraint to be considered dense. +- `threshold_small`: The threshold for small coefficients in the model. +- `threshold_large`: The threshold for large coefficients in the model. + +""" +struct Analyzer <: ModelAnalyzer.AbstractAnalyzer end + +""" + AbstractNumericalIssue <: AbstractNumericalIssue + +Abstract type for numerical issues found during the analysis of a model. +""" +abstract type AbstractNumericalIssue <: ModelAnalyzer.AbstractIssue end + +""" + VariableNotInConstraints <: AbstractNumericalIssue + +The `VariableNotInConstraints` issue is identified when a variable appears in no +constraints. + +For more information, run: +```julia +julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.VariableNotInConstraints) +``` +""" +struct VariableNotInConstraints <: AbstractNumericalIssue + ref::MOI.VariableIndex +end + +ModelAnalyzer.variable(issue::VariableNotInConstraints) = issue.ref + +""" + EmptyConstraint <: AbstractNumericalIssue + +The `EmptyConstraint` issue is identified when a constraint has no coefficients +different from zero. + +For more information, run: +```julia +julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.EmptyConstraint) +``` +""" +struct EmptyConstraint <: AbstractNumericalIssue + ref::MOI.ConstraintIndex +end + +ModelAnalyzer.constraint(issue::EmptyConstraint) = issue.ref + +""" + VariableBoundAsConstraint <: AbstractNumericalIssue + +The `VariableBoundAsConstraint` issue is identified when a constraint is +equivalent to a variable bound, that is, the constraint has only one non-zero +coefficient, and this coefficient is equal to one. + +For more information, run: +```julia +julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.VariableBoundAsConstraint) +``` +""" +struct VariableBoundAsConstraint <: AbstractNumericalIssue + ref::MOI.ConstraintIndex +end + +ModelAnalyzer.constraint(issue::VariableBoundAsConstraint) = issue.ref + +""" + DenseConstraint <: AbstractNumericalIssue + +The `DenseConstraint` issue is identified when a constraint has a fraction of +non-zero entries greater than `threshold_dense_fill_in` and the number of +non-zero entries is greater than `threshold_dense_entries`. + +For more information, run: +```julia +julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.DenseConstraint) +``` +""" +struct DenseConstraint <: AbstractNumericalIssue + ref::MOI.ConstraintIndex + nnz::Int +end + +ModelAnalyzer.constraint(issue::DenseConstraint) = issue.ref + +ModelAnalyzer.value(issue::DenseConstraint) = issue.nnz + +""" + SmallMatrixCoefficient <: AbstractNumericalIssue + +The `SmallMatrixCoefficient` issue is identified when a matrix coefficient in a +constraint is smaller than `threshold_small`. + +For more information, run: +```julia +julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.SmallMatrixCoefficient) +``` +""" +struct SmallMatrixCoefficient <: AbstractNumericalIssue + 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 + +The `LargeMatrixCoefficient` issue is identified when a matrix coefficient in a +constraint is larger than `threshold_large`. + +For more information, run: +```julia +julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.LargeMatrixCoefficient) +``` +""" +struct LargeMatrixCoefficient <: AbstractNumericalIssue + 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 + +The `SmallBoundCoefficient` issue is identified when a variable's bound +(coefficient) is smaller than `threshold_small`. + +For more information, run: +```julia +julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.SmallBoundCoefficient) +``` +""" +struct SmallBoundCoefficient <: AbstractNumericalIssue + variable::MOI.VariableIndex + coefficient::Float64 +end + +ModelAnalyzer.variable(issue::SmallBoundCoefficient) = issue.variable + +ModelAnalyzer.value(issue::SmallBoundCoefficient) = issue.coefficient + +""" + LargeBoundCoefficient <: AbstractNumericalIssue + +The `LargeBoundCoefficient` issue is identified when a variable's bound +(coefficient) is larger than `threshold_large`. + +For more information, run: +```julia +julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.LargeBoundCoefficient) +``` +""" +struct LargeBoundCoefficient <: AbstractNumericalIssue + variable::MOI.VariableIndex + coefficient::Float64 +end + +ModelAnalyzer.variable(issue::LargeBoundCoefficient) = issue.variable + +ModelAnalyzer.value(issue::LargeBoundCoefficient) = issue.coefficient + +""" + SmallRHSCoefficient <: AbstractNumericalIssue + +The `SmallRHSCoefficient` issue is identified when the right-hand-side (RHS) +coefficient of a constraint is smaller than `threshold_small`. + +For more information, run: +```julia +julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.SmallRHSCoefficient) +``` +""" +struct SmallRHSCoefficient <: AbstractNumericalIssue + ref::MOI.ConstraintIndex + coefficient::Float64 +end + +ModelAnalyzer.constraint(issue::SmallRHSCoefficient) = issue.ref + +ModelAnalyzer.value(issue::SmallRHSCoefficient) = issue.coefficient + +""" + LargeRHSCoefficient <: AbstractNumericalIssue + +The `LargeRHSCoefficient` issue is identified when the right-hand-side (RHS) +coefficient of a constraint is larger than `threshold_large`. + +For more information, run: +```julia +julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.LargeRHSCoefficient) +``` +""" +struct LargeRHSCoefficient <: AbstractNumericalIssue + ref::MOI.ConstraintIndex + coefficient::Float64 +end + +ModelAnalyzer.constraint(issue::LargeRHSCoefficient) = issue.ref + +ModelAnalyzer.value(issue::LargeRHSCoefficient) = issue.coefficient + +""" + SmallObjectiveCoefficient <: AbstractNumericalIssue + +The `SmallObjectiveCoefficient` issue is identified when a coefficient in the +objective function is smaller than `threshold_small`. + +For more information, run: +```julia +julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.SmallObjectiveCoefficient) +``` +""" +struct SmallObjectiveCoefficient <: AbstractNumericalIssue + variable::MOI.VariableIndex + coefficient::Float64 +end + +ModelAnalyzer.variable(issue::SmallObjectiveCoefficient) = issue.variable + +ModelAnalyzer.value(issue::SmallObjectiveCoefficient) = issue.coefficient + +""" + LargeObjectiveCoefficient <: AbstractNumericalIssue + +The `LargeObjectiveCoefficient` issue is identified when a coefficient in the +objective function is larger than `threshold_large`. + +For more information, run: +```julia +julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.LargeObjectiveCoefficient) +``` +""" +struct LargeObjectiveCoefficient <: AbstractNumericalIssue + variable::MOI.VariableIndex + coefficient::Float64 +end + +ModelAnalyzer.variable(issue::LargeObjectiveCoefficient) = issue.variable + +ModelAnalyzer.value(issue::LargeObjectiveCoefficient) = issue.coefficient + +""" + SmallObjectiveQuadraticCoefficient <: AbstractNumericalIssue + +The `SmallObjectiveQuadraticCoefficient` issue is identified when a quadratic +coefficient in the objective function is smaller than `threshold_small`. + +For more information, run: +```julia +julia> ModelAnalyzer.summarize( + ModelAnalyzer.Numerical.SmallObjectiveQuadraticCoefficient +) +``` +""" +struct SmallObjectiveQuadraticCoefficient <: AbstractNumericalIssue + 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 + +The `LargeObjectiveQuadraticCoefficient` issue is identified when a quadratic +coefficient in the objective function is larger than `threshold_large`. + +For more information, run: +```julia +julia> ModelAnalyzer.summarize( + ModelAnalyzer.Numerical.LargeObjectiveQuadraticCoefficient +) +``` +""" +struct LargeObjectiveQuadraticCoefficient <: AbstractNumericalIssue + 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 + +The `SmallMatrixQuadraticCoefficient` issue is identified when a quadratic +coefficient in a constraint is smaller than `threshold_small`. + +For more information, run: +```julia +julia> ModelAnalyzer.summarize( + ModelAnalyzer.Numerical.SmallMatrixQuadraticCoefficient +) +``` +""" +struct SmallMatrixQuadraticCoefficient <: AbstractNumericalIssue + 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 + +The `LargeMatrixQuadraticCoefficient` issue is identified when a quadratic +coefficient in a constraint is larger than `threshold_large`. + +For more information, run: +```julia +julia> ModelAnalyzer.summarize( + ModelAnalyzer.Numerical.LargeMatrixQuadraticCoefficient +) +``` +""" +struct LargeMatrixQuadraticCoefficient <: AbstractNumericalIssue + 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 + +The `NonconvexQuadraticObjective` issue is identified when a quadratic +objective function is non-convex. + +For more information, run: +```julia +julia> ModelAnalyzer.summarize( + ModelAnalyzer.Numerical.NonconvexQuadraticObjective +) +``` +""" +struct NonconvexQuadraticObjective <: AbstractNumericalIssue end + +""" + NonconvexQuadraticConstraint + +The `NonconvexQuadraticConstraint` issue is identified when a quadratic +constraint is non-convex. + +For more information, run: +```julia +julia> ModelAnalyzer.summarize( + ModelAnalyzer.Numerical.NonconvexQuadraticConstraint +) +``` +""" +struct NonconvexQuadraticConstraint <: AbstractNumericalIssue + ref::MOI.ConstraintIndex +end +ModelAnalyzer.constraint(issue::NonconvexQuadraticConstraint) = issue.ref + +""" + Data + +The `Data` structure holds the results of the analysis performed by the +`ModelAnalyzer.Numerical.Analyzer`. It contains various thresholds and the +information about the model's variables, constraints, and objective function. +""" +Base.@kwdef mutable struct Data <: ModelAnalyzer.AbstractData + # analysis configuration + threshold_dense_fill_in::Float64 = 0.10 + threshold_dense_entries::Int = 1000 + threshold_small::Float64 = 1e-5 + threshold_large::Float64 = 1e+5 + # main numbers + number_of_variables::Int = 0 + number_of_constraints::Int = 0 + constraint_info::Vector{Tuple{DataType,DataType,Int}} = + Tuple{DataType,DataType,Int}[] + # objective_info::Any + matrix_nnz::Int = 0 + # ranges + matrix_range::Vector{Float64} = sizehint!(Float64[1.0, 1.0], 2) + bounds_range::Vector{Float64} = sizehint!(Float64[1.0, 1.0], 2) + 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{MOI.VariableIndex} = Set{MOI.VariableIndex}() + # variables analysis + variables_not_in_constraints::Vector{VariableNotInConstraints} = + VariableNotInConstraints[] + bounds_small::Vector{SmallBoundCoefficient} = SmallBoundCoefficient[] + bounds_large::Vector{LargeBoundCoefficient} = LargeBoundCoefficient[] + # constraints analysis + empty_rows::Vector{EmptyConstraint} = EmptyConstraint[] + bound_rows::Vector{VariableBoundAsConstraint} = VariableBoundAsConstraint[] + dense_rows::Vector{DenseConstraint} = DenseConstraint[] + matrix_small::Vector{SmallMatrixCoefficient} = SmallMatrixCoefficient[] + matrix_large::Vector{LargeMatrixCoefficient} = LargeMatrixCoefficient[] + rhs_small::Vector{SmallRHSCoefficient} = SmallRHSCoefficient[] + rhs_large::Vector{LargeRHSCoefficient} = LargeRHSCoefficient[] + # quadratic constraints analysis + has_quadratic_constraints::Bool = false + nonconvex_rows::Vector{NonconvexQuadraticConstraint} = + NonconvexQuadraticConstraint[] + matrix_quadratic_small::Vector{SmallMatrixQuadraticCoefficient} = + SmallMatrixQuadraticCoefficient[] + matrix_quadratic_large::Vector{LargeMatrixQuadraticCoefficient} = + LargeMatrixQuadraticCoefficient[] + # cache data + sense::MOI.OptimizationSense = MOI.FEASIBILITY_SENSE + # objective analysis + objective_small::Vector{SmallObjectiveCoefficient} = + SmallObjectiveCoefficient[] + objective_large::Vector{LargeObjectiveCoefficient} = + LargeObjectiveCoefficient[] + # quadratic objective analysis + has_quadratic_objective::Bool = false + objective_quadratic_range::Vector{Float64} = sizehint!(Float64[1.0, 1.0], 2) + matrix_quadratic_range::Vector{Float64} = sizehint!(Float64[1.0, 1.0], 2) + nonconvex_objective::Vector{NonconvexQuadraticObjective} = + NonconvexQuadraticObjective[] + objective_quadratic_small::Vector{SmallObjectiveQuadraticCoefficient} = + SmallObjectiveQuadraticCoefficient[] + objective_quadratic_large::Vector{LargeObjectiveQuadraticCoefficient} = + LargeObjectiveQuadraticCoefficient[] +end diff --git a/src/Numerical/summarize.jl b/src/Numerical/summarize.jl new file mode 100644 index 0000000..02ff118 --- /dev/null +++ b/src/Numerical/summarize.jl @@ -0,0 +1,1323 @@ +# 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 ModelAnalyzer._summarize(io::IO, ::Type{VariableNotInConstraints}) + return print(io, "# VariableNotInConstraints") +end + +function ModelAnalyzer._summarize(io::IO, ::Type{EmptyConstraint}) + return print(io, "# EmptyConstraint") +end + +function ModelAnalyzer._summarize(io::IO, ::Type{VariableBoundAsConstraint}) + return print(io, "# VariableBoundAsConstraint") +end + +function ModelAnalyzer._summarize(io::IO, ::Type{DenseConstraint}) + return print(io, "# DenseConstraint") +end + +function ModelAnalyzer._summarize(io::IO, ::Type{SmallMatrixCoefficient}) + return print(io, "# SmallMatrixCoefficient") +end + +function ModelAnalyzer._summarize(io::IO, ::Type{LargeMatrixCoefficient}) + return print(io, "# LargeMatrixCoefficient") +end + +function ModelAnalyzer._summarize(io::IO, ::Type{SmallBoundCoefficient}) + return print(io, "# SmallBoundCoefficient") +end + +function ModelAnalyzer._summarize(io::IO, ::Type{LargeBoundCoefficient}) + return print(io, "# LargeBoundCoefficient") +end + +function ModelAnalyzer._summarize(io::IO, ::Type{SmallRHSCoefficient}) + return print(io, "# SmallRHSCoefficient") +end + +function ModelAnalyzer._summarize(io::IO, ::Type{LargeRHSCoefficient}) + return print(io, "# LargeRHSCoefficient") +end + +function ModelAnalyzer._summarize(io::IO, ::Type{SmallObjectiveCoefficient}) + return print(io, "# SmallObjectiveCoefficient") +end + +function ModelAnalyzer._summarize(io::IO, ::Type{LargeObjectiveCoefficient}) + return print(io, "# LargeObjectiveCoefficient") +end + +function ModelAnalyzer._summarize( + io::IO, + ::Type{SmallObjectiveQuadraticCoefficient}, +) + return print(io, "# SmallObjectiveQuadraticCoefficient") +end + +function ModelAnalyzer._summarize( + io::IO, + ::Type{LargeObjectiveQuadraticCoefficient}, +) + return print(io, "# LargeObjectiveQuadraticCoefficient") +end + +function ModelAnalyzer._summarize( + io::IO, + ::Type{SmallMatrixQuadraticCoefficient}, +) + return print(io, "# SmallMatrixQuadraticCoefficient") +end + +function ModelAnalyzer._summarize( + io::IO, + ::Type{LargeMatrixQuadraticCoefficient}, +) + return print(io, "# LargeMatrixQuadraticCoefficient") +end + +function ModelAnalyzer._summarize(io::IO, ::Type{NonconvexQuadraticObjective}) + return print(io, "# NonconvexQuadraticObjective") +end + +function ModelAnalyzer._summarize(io::IO, ::Type{NonconvexQuadraticConstraint}) + return print(io, "# NonconvexQuadraticConstraint") +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + ::Type{VariableNotInConstraints}, +) + return print( + io, + """ + # `VariableNotInConstraints` + + ## What + + A `VariableNotInConstraints` issue is identified when a variable appears + in no constraints. If a variable only appears alone in a constraint and + it has a coefficient of 1 it is considered a + `VariableNotInConstraints`, because this emulates a bound. + + ## Why + + This can be a sign of a mistake in the model formulation. If a variable + is not used in any constraints, it is not affecting the solution of the + problem. Moreover, it might be leading to an unbounded problem. + + ## How to fix + + If the variable is not needed, remove it from the model. If the variable + is needed, check that it is correctly used in the constraints. + + ## More information + + No extra information for this issue. + """, + ) +end + +function ModelAnalyzer._verbose_summarize(io::IO, ::Type{EmptyConstraint}) + return print( + io, + """ + # `EmptyConstraint` + + ## What + + An `EmptyConstraint` issue is identified when a constraint has no + coefficients different from zero. + + ## Why + + This can be a sign of a mistake in the model formulation. An empty + constraint is not affecting the solution of the problem. Moreover, it + might be leading to an infeasible problem since the \"left-hand-side\" + of the constraint is always zero. + + ## How to fix + + Remove the empty constraint from the model. + + ## More information + + No extra information for this issue. + """, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + ::Type{VariableBoundAsConstraint}, +) + return print( + io, + """ + # `VariableBoundAsConstraint` + + ## What + + A `VariableBoundAsConstraint` issue is identified when a constraint is + equivalent to a variable bound, that is, the constraint has only one + non-zero coefficient, and this coefficient is equal to one. + + ## Why + + This can be a sign of a mistake in the model formulation. Variable + bounds are frequently used by solver in special ways that can lead to + better performance. + + ## How to fix + + Remove the constraint and use the variable bound directly. + + ## More information + + - https://support.gurobi.com/hc/en-us/community/posts/24066470832529/comments/24183896218385 + """, + ) +end + +function ModelAnalyzer._verbose_summarize(io::IO, ::Type{DenseConstraint}) + return print( + io, + """ + # `DenseConstraint` + + ## What + + A `DenseConstraint` issue is identified when a constraint has a high + number of non-zero coefficients. + + ## Why + + Dense constraints can lead to performance issues in the solution + process. Very few dense constraints might not be a problem. + + ## How to fix + + Check if the constraint can be simplified. A common + case that can be avoided is when there is an expression + `e = c1 * x1 + c2 * x2 + ... + cn * xn` where `c1, c2, ..., cn` are + constants and `x1, x2, ..., xn` are variables, and this expression is + used in many constraints. In this case, it is better to create a new + variable `y = e` and use `y` in the constraints. + + ## More information + + No extra information for this issue. + """, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + ::Type{SmallMatrixCoefficient}, +) + return print( + io, + """ + # `SmallMatrixCoefficient` + + ## What + + A `SmallMatrixCoefficient` issue is identified when a constraint has a + coefficient with a small absolute value. + + ## Why + + Small coefficients can lead to numerical instability in the solution + process. + + ## How to fix + + Check if the coefficient is correct. Check if the units of variables and + coefficients are correct. Check if the number makes is + reasonable given that solver have tolerances. Sometimes these + coefficients can be replaced by zeros. + + ## More information + + - https://jump.dev/JuMP.jl/stable/tutorials/getting_started/tolerances/ + """, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + ::Type{LargeMatrixCoefficient}, +) + return print( + io, + """ + # `LargeMatrixCoefficient` + + ## What + + A `LargeMatrixCoefficient` issue is identified when a constraint has a + coefficient with a large absolute value. + + ## Why + + Large coefficients can lead to numerical instability in the solution + process. + + ## How to fix + + Check if the coefficient is correct. Check if the units of variables and + coefficients are correct. Check if the number makes is + reasonable given that solver have tolerances. Sometimes these + coefficients can be replaced by zeros. + + ## More information + + - https://jump.dev/JuMP.jl/stable/tutorials/getting_started/tolerances/ + """, + ) +end + +function ModelAnalyzer._verbose_summarize(io::IO, ::Type{SmallBoundCoefficient}) + return print( + io, + """ + # `SmallBoundCoefficient` + + ## What + + A `SmallBoundCoefficient` issue is identified when a variable has a + bound with a small absolute value. + + ## Why + + Small bounds can lead to numerical instability in the solution process. + + ## How to fix + + Check if the bound is correct. Check if the units of variables and + coefficients are correct. Check if the number makes is + reasonable given that solver have tolerances. Sometimes these + bounds can be replaced by zeros. + + ## More information + + - https://jump.dev/JuMP.jl/stable/tutorials/getting_started/tolerances/ + """, + ) +end + +function ModelAnalyzer._verbose_summarize(io::IO, ::Type{LargeBoundCoefficient}) + return print( + io, + """ + # `LargeBoundCoefficient` + + ## What + + A `LargeBoundCoefficient` issue is identified when a variable has a + bound with a large absolute value. + + ## Why + + Large bounds can lead to numerical instability in the solution process. + + ## How to fix + + Check if the bound is correct. Check if the units of variables and + coefficients are correct. Check if the number makes is + reasonable given that solver have tolerances. Sometimes these + bounds can be replaced by zeros. + + ## More information + + - https://jump.dev/JuMP.jl/stable/tutorials/getting_started/tolerances/ + """, + ) +end + +function ModelAnalyzer._verbose_summarize(io::IO, ::Type{SmallRHSCoefficient}) + return print( + io, + """ + # `SmallRHSCoefficient` + + ## What + + A `SmallRHSCoefficient` issue is identified when a constraint has a + right-hand-side with a small absolute value. + + ## Why + + Small right-hand-sides can lead to numerical instability in the solution + process. + + ## How to fix + + Check if the right-hand-side is correct. Check if the units of variables + and coefficients are correct. Check if the number makes is + reasonable given that solver have tolerances. Sometimes these + right-hand-sides can be replaced by zeros. + + ## More information + + - https://jump.dev/JuMP.jl/stable/tutorials/getting_started/tolerances/ + """, + ) +end + +function ModelAnalyzer._verbose_summarize(io::IO, ::Type{LargeRHSCoefficient}) + return print( + io, + """ + # `LargeRHSCoefficient` + + ## What + + A `LargeRHSCoefficient` issue is identified when a constraint has a + right-hand-side with a large absolute value. + + ## Why + + Large right-hand-sides can lead to numerical instability in the solution + process. + + ## How to fix + + Check if the right-hand-side is correct. Check if the units of variables + and coefficients are correct. Check if the number makes is + reasonable given that solver have tolerances. Sometimes these + right-hand-sides can be replaced by zeros. + + ## More information + + - https://jump.dev/JuMP.jl/stable/tutorials/getting_started/tolerances/ + """, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + ::Type{SmallObjectiveCoefficient}, +) + return print( + io, + """ + # `SmallObjectiveCoefficient` + + ## What + + A `SmallObjectiveCoefficient` issue is identified when the objective + function has a coefficient with a small absolute value. + + ## Why + + Small coefficients can lead to numerical instability in the solution + process. + + ## How to fix + + Check if the coefficient is correct. Check if the units of variables and + coefficients are correct. Check if the number makes is + reasonable given that solver have tolerances. Sometimes these + coefficients can be replaced by zeros. + + ## More information + + - https://jump.dev/JuMP.jl/stable/tutorials/getting_started/tolerances/ + """, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + ::Type{LargeObjectiveCoefficient}, +) + return print( + io, + """ + # `LargeObjectiveCoefficient` + + ## What + + A `LargeObjectiveCoefficient` issue is identified when the objective + function has a coefficient with a large absolute value. + + ## Why + + Large coefficients can lead to numerical instability in the solution + process. + + ## How to fix + + Check if the coefficient is correct. Check if the units of variables and + coefficients are correct. Check if the number makes is + reasonable given that solver have tolerances. Sometimes these + coefficients can be replaced by zeros. + + ## More information + + - https://jump.dev/JuMP.jl/stable/tutorials/getting_started/tolerances/ + """, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + ::Type{SmallObjectiveQuadraticCoefficient}, +) + return print( + io, + """ + # `SmallObjectiveQuadraticCoefficient` + + ## What + + A `SmallObjectiveQuadraticCoefficient` issue is identified when the + objective function has a quadratic coefficient with a small absolute value. + + ## Why + + Small coefficients can lead to numerical instability in the solution + process. + + ## How to fix + + Check if the coefficient is correct. Check if the units of variables and + coefficients are correct. Check if the number makes is + reasonable given that solver have tolerances. Sometimes these + coefficients can be replaced by zeros. + + ## More information + + - https://jump.dev/JuMP.jl/stable/tutorials/getting_started/tolerances/ + """, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + ::Type{LargeObjectiveQuadraticCoefficient}, +) + return print( + io, + """ + # `LargeObjectiveQuadraticCoefficient` + + ## What + + A `LargeObjectiveQuadraticCoefficient` issue is identified when the + objective function has a quadratic coefficient with a large absolute value. + + ## Why + + Large coefficients can lead to numerical instability in the solution + process. + + ## How to fix + + Check if the coefficient is correct. Check if the units of variables and + coefficients are correct. Check if the number makes is + reasonable given that solver have tolerances. Sometimes these + coefficients can be replaced by zeros. + + ## More information + + - https://jump.dev/JuMP.jl/stable/tutorials/getting_started/tolerances/ + """, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + ::Type{SmallMatrixQuadraticCoefficient}, +) + return print( + io, + """ + # `SmallMatrixQuadraticCoefficient` + + ## What + + A `SmallMatrixQuadraticCoefficient` issue is identified when a quadratic + constraint has a coefficient with a small absolute value. + + ## Why + + Small coefficients can lead to numerical instability in the solution + process. + + ## How to fix + + Check if the coefficient is correct. Check if the units of variables and + coefficients are correct. Check if the number makes is + reasonable given that solver have tolerances. Sometimes these + coefficients can be replaced by zeros. + + ## More information + + - https://jump.dev/JuMP.jl/stable/tutorials/getting_started/tolerances/ + """, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + ::Type{NonconvexQuadraticObjective}, +) + return print( + io, + """ + # `NonconvexQuadraticObjective` + + ## What + + A `NonconvexQuadraticObjective` issue is identified when a quadratic + objective is nonconvex, that is, the quadratic matrix is not positive + semidefinite for minimization or the quadratic matrix is not negative + semidefinite for maximization. + + ## Why + + Nonconvex objectives are not expected by many solver and can lead to + wrong solutions or even convergence issues. + + ## How to fix + + Check if the objective is correct. Coefficient signs might have been + inverted. This also occurs if user fix a variable to emulate a + parameter, in this case some solvers will not be able to solve the + model properly, other tools such as ParametricOptInteface.jl might be + more suitable than fixing variables. + + ## More information + + No extra information for this issue. + """, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + ::Type{NonconvexQuadraticConstraint}, +) + return print( + io, + """ + # `NonconvexQuadraticConstraint` + + ## What + + A `NonconvexQuadraticConstraint` issue is identified when a quadratic + constraint is nonconvex, that is, the quadratic matrix is not positive + semidefinite. + + ## Why + + Nonconvex constraints are not expected by many solver and can lead to + wrong solutions or even convergence issues. + + ## How to fix + + Check if the constraint is correct. Coefficient signs might have been + inverted. This also occurs if user fix a variable to emulate a + parameter, in this case some solvers will not be able to solve the + model properly, other tools such as ParametricOptInteface.jl might be + more suitable than fixing variables. + + ## More information + + No extra information for this issue. + """, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + ::Type{LargeMatrixQuadraticCoefficient}, +) + return print( + io, + """ + # `LargeMatrixQuadraticCoefficient` + + ## What + + A `LargeMatrixQuadraticCoefficient` issue is identified when a quadratic + constraint has a coefficient with a large absolute value. + + ## Why + + Large coefficients can lead to numerical instability in the solution + process. + + ## How to fix + + Check if the coefficient is correct. Check if the units of variables and + coefficients are correct. Check if the number makes is + reasonable given that solver have tolerances. Sometimes these + coefficients can be replaced by zeros. + + ## More information + + - https://jump.dev/JuMP.jl/stable/tutorials/getting_started/tolerances/ + """, + ) +end + +function ModelAnalyzer._summarize( + io::IO, + issue::VariableNotInConstraints, + model, +) + return print(io, ModelAnalyzer._name(issue.ref, model)) +end + +function ModelAnalyzer._summarize(io::IO, issue::EmptyConstraint, model) + return print(io, ModelAnalyzer._name(issue.ref, model)) +end + +function ModelAnalyzer._summarize( + io::IO, + issue::VariableBoundAsConstraint, + model, +) + return print(io, ModelAnalyzer._name(issue.ref, model)) +end + +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, model) + return print( + io, + ModelAnalyzer._name(issue.ref, model), + " -- ", + ModelAnalyzer._name(issue.variable, model), + " : ", + issue.coefficient, + ) +end + +function ModelAnalyzer._summarize(io::IO, issue::LargeMatrixCoefficient, model) + return print( + io, + ModelAnalyzer._name(issue.ref, model), + " -- ", + ModelAnalyzer._name(issue.variable, model), + " : ", + issue.coefficient, + ) +end + +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, model) + return print( + io, + ModelAnalyzer._name(issue.variable, model), + " : ", + issue.coefficient, + ) +end + +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, model) + return print( + io, + ModelAnalyzer._name(issue.ref, model), + " : ", + issue.coefficient, + ) +end + +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, + model, +) + return print( + io, + ModelAnalyzer._name(issue.variable, model), + " : ", + issue.coefficient, + ) +end + +function ModelAnalyzer._summarize( + io::IO, + issue::SmallObjectiveQuadraticCoefficient, + model, +) + return print( + io, + ModelAnalyzer._name(issue.variable1, model), + " -- ", + ModelAnalyzer._name(issue.variable2, model), + " : ", + issue.coefficient, + ) +end + +function ModelAnalyzer._summarize( + io::IO, + issue::LargeObjectiveQuadraticCoefficient, + model, +) + return print( + io, + ModelAnalyzer._name(issue.variable1, model), + " -- ", + ModelAnalyzer._name(issue.variable2, model), + " : ", + issue.coefficient, + ) +end + +function ModelAnalyzer._summarize( + io::IO, + issue::SmallMatrixQuadraticCoefficient, + model, +) + return print( + io, + ModelAnalyzer._name(issue.ref, model), + " -- ", + ModelAnalyzer._name(issue.variable1, model), + " -- ", + ModelAnalyzer._name(issue.variable2, model), + " : ", + issue.coefficient, + ) +end + +function ModelAnalyzer._summarize( + io::IO, + issue::LargeMatrixQuadraticCoefficient, + model, +) + return print( + io, + ModelAnalyzer._name(issue.ref, model), + " -- ", + ModelAnalyzer._name(issue.variable1, model), + " -- ", + ModelAnalyzer._name(issue.variable2, model), + " : ", + issue.coefficient, + ) +end + +function ModelAnalyzer._summarize(io::IO, ::NonconvexQuadraticObjective, model) + return print(io, "Objective is Nonconvex quadratic") +end + +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: ", ModelAnalyzer._name(issue.ref, model)) +end + +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: ", ModelAnalyzer._name(issue.ref, model)) +end + +function ModelAnalyzer._verbose_summarize(io::IO, issue::DenseConstraint, model) + return print( + io, + "Constraint: ", + ModelAnalyzer._name(issue.ref, model), + " with ", + issue.nnz, + " non zero coefficients", + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::SmallMatrixCoefficient, + model, +) + return print( + io, + "(Constraint -- Variable): (", + ModelAnalyzer._name(issue.ref, model), + " -- ", + ModelAnalyzer._name(issue.variable, model), + ") with coefficient ", + issue.coefficient, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::LargeMatrixCoefficient, + model, +) + return print( + io, + "(Constraint -- Variable): (", + ModelAnalyzer._name(issue.ref, model), + " -- ", + ModelAnalyzer._name(issue.variable, model), + ") with coefficient ", + issue.coefficient, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::SmallBoundCoefficient, + model, +) + return print( + io, + "Variable: ", + ModelAnalyzer._name(issue.variable, model), + " with bound ", + issue.coefficient, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::LargeBoundCoefficient, + model, +) + return print( + io, + "Variable: ", + ModelAnalyzer._name(issue.variable, model), + " with bound ", + issue.coefficient, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::SmallRHSCoefficient, + model, +) + return print( + io, + "Constraint: ", + ModelAnalyzer._name(issue.ref, model), + " with right-hand-side ", + issue.coefficient, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::LargeRHSCoefficient, + model, +) + return print( + io, + "Constraint: ", + ModelAnalyzer._name(issue.ref, model), + " with right-hand-side ", + issue.coefficient, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::SmallObjectiveCoefficient, + model, +) + return print( + io, + "Variable: ", + ModelAnalyzer._name(issue.variable, model), + " with coefficient ", + issue.coefficient, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::LargeObjectiveCoefficient, + model, +) + return print( + io, + "Variable: ", + ModelAnalyzer._name(issue.variable, model), + " with coefficient ", + issue.coefficient, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::SmallObjectiveQuadraticCoefficient, + model, +) + return print( + io, + "(Variable -- Variable): (", + ModelAnalyzer._name(issue.variable1, model), + " -- ", + ModelAnalyzer._name(issue.variable2, model), + ") with coefficient ", + issue.coefficient, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::LargeObjectiveQuadraticCoefficient, + model, +) + return print( + io, + "(Variable -- Variable): (", + ModelAnalyzer._name(issue.variable1, model), + " -- ", + ModelAnalyzer._name(issue.variable2, model), + ") with coefficient ", + issue.coefficient, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::SmallMatrixQuadraticCoefficient, + model, +) + return print( + io, + "(Constraint -- Variable -- Variable): (", + ModelAnalyzer._name(issue.ref, model), + " -- ", + ModelAnalyzer._name(issue.variable1, model), + " -- ", + ModelAnalyzer._name(issue.variable2, model), + ") with coefficient ", + issue.coefficient, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::LargeMatrixQuadraticCoefficient, + model, +) + return print( + io, + "(Constraint -- Variable -- Variable): (", + ModelAnalyzer._name(issue.ref, model), + " -- ", + ModelAnalyzer._name(issue.variable1, model), + " -- ", + ModelAnalyzer._name(issue.variable2, model), + ") with coefficient ", + issue.coefficient, + ) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::NonconvexQuadraticObjective, + model, +) + return ModelAnalyzer._summarize(io, issue, model) +end + +function ModelAnalyzer._verbose_summarize( + io::IO, + issue::NonconvexQuadraticConstraint, + model, +) + return print(io, "Constraint: ", ModelAnalyzer._name(issue.ref, model)) +end + +function ModelAnalyzer.list_of_issues( + data::Data, + ::Type{VariableNotInConstraints}, +) + return data.variables_not_in_constraints +end + +function ModelAnalyzer.list_of_issues(data::Data, ::Type{EmptyConstraint}) + return data.empty_rows +end + +function ModelAnalyzer.list_of_issues( + data::Data, + ::Type{VariableBoundAsConstraint}, +) + return data.bound_rows +end + +function ModelAnalyzer.list_of_issues(data::Data, ::Type{DenseConstraint}) + return data.dense_rows +end + +function ModelAnalyzer.list_of_issues( + data::Data, + ::Type{SmallMatrixCoefficient}, +) + return data.matrix_small +end + +function ModelAnalyzer.list_of_issues( + data::Data, + ::Type{LargeMatrixCoefficient}, +) + return data.matrix_large +end + +function ModelAnalyzer.list_of_issues(data::Data, ::Type{SmallBoundCoefficient}) + return data.bounds_small +end + +function ModelAnalyzer.list_of_issues(data::Data, ::Type{LargeBoundCoefficient}) + return data.bounds_large +end + +function ModelAnalyzer.list_of_issues(data::Data, ::Type{SmallRHSCoefficient}) + return data.rhs_small +end + +function ModelAnalyzer.list_of_issues(data::Data, ::Type{LargeRHSCoefficient}) + return data.rhs_large +end + +function ModelAnalyzer.list_of_issues( + data::Data, + ::Type{SmallObjectiveCoefficient}, +) + return data.objective_small +end + +function ModelAnalyzer.list_of_issues( + data::Data, + ::Type{LargeObjectiveCoefficient}, +) + return data.objective_large +end + +function ModelAnalyzer.list_of_issues( + data::Data, + ::Type{SmallObjectiveQuadraticCoefficient}, +) + return data.objective_quadratic_small +end + +function ModelAnalyzer.list_of_issues( + data::Data, + ::Type{LargeObjectiveQuadraticCoefficient}, +) + return data.objective_quadratic_large +end + +function ModelAnalyzer.list_of_issues( + data::Data, + ::Type{SmallMatrixQuadraticCoefficient}, +) + return data.matrix_quadratic_small +end + +function ModelAnalyzer.list_of_issues( + data::Data, + ::Type{LargeMatrixQuadraticCoefficient}, +) + return data.matrix_quadratic_large +end + +function ModelAnalyzer.list_of_issues( + data::Data, + ::Type{NonconvexQuadraticObjective}, +) + return data.nonconvex_objective +end + +function ModelAnalyzer.list_of_issues( + data::Data, + ::Type{NonconvexQuadraticConstraint}, +) + return data.nonconvex_rows +end + +function ModelAnalyzer.list_of_issue_types(data::Data) + ret = Type[] + for type in ( + VariableNotInConstraints, + EmptyConstraint, + VariableBoundAsConstraint, + DenseConstraint, + SmallMatrixCoefficient, + LargeMatrixCoefficient, + SmallBoundCoefficient, + LargeBoundCoefficient, + SmallRHSCoefficient, + LargeRHSCoefficient, + SmallObjectiveCoefficient, + LargeObjectiveCoefficient, + SmallObjectiveQuadraticCoefficient, + LargeObjectiveQuadraticCoefficient, + SmallMatrixQuadraticCoefficient, + LargeMatrixQuadraticCoefficient, + NonconvexQuadraticConstraint, + NonconvexQuadraticObjective, + ) + if !isempty(ModelAnalyzer.list_of_issues(data, type)) + push!(ret, type) + end + end + return ret +end + +function summarize_configurations(io::IO, data::Data) + print(io, "## Configuration\n\n") + print(io, " Dense fill-in threshold: ", data.threshold_dense_fill_in, "\n") + print(io, " Dense entries threshold: ", data.threshold_dense_entries, "\n") + print(io, " Small coefficient threshold: ", data.threshold_small, "\n") + print(io, " Large coefficient threshold: ", data.threshold_large, "\n") + return +end + +function summarize_dimensions(io::IO, data::Data) + print(io, "## Dimensions\n\n") + print(io, " Number of variables: ", data.number_of_variables, "\n") + print(io, " Number of constraints: ", data.number_of_constraints, "\n") + print(io, " Number of nonzeros in matrix: ", data.matrix_nnz, "\n") + # types + println(io, " Constraint types:") + for (F, S, n) in data.constraint_info + println(io, " * ", F, "-", S, ": ", n) + end + return +end + +function summarize_ranges(io::IO, data::Data) + print(io, "## Coefficient ranges\n\n") + print(io, " Matrix: ", _stringify_bounds(data.matrix_range), "\n") + print(io, " Objective: ", _stringify_bounds(data.objective_range), "\n") + print(io, " Bounds: ", _stringify_bounds(data.bounds_range), "\n") + print(io, " RHS: ", _stringify_bounds(data.rhs_range), "\n") + if data.has_quadratic_objective + print( + io, + " Objective quadratic: ", + _stringify_bounds(data.objective_quadratic_range), + "\n", + ) + end + if data.has_quadratic_constraints + print( + io, + " Matrix quadratic: ", + _stringify_bounds(data.matrix_quadratic_range), + "\n", + ) + end + return +end + +function ModelAnalyzer.summarize( + io::IO, + data::Data; + model = nothing, + verbose = true, + max_issues = ModelAnalyzer.DEFAULT_MAX_ISSUES, + configurations = true, + dimensions = true, + ranges = true, +) + print(io, "## Numerical Analysis\n\n") + if configurations + summarize_configurations(io, data) + print(io, "\n") + end + if dimensions + summarize_dimensions(io, data) + print(io, "\n") + end + if ranges + summarize_ranges(io, data) + print(io, "\n") + end + for issue_type in ModelAnalyzer.list_of_issue_types(data) + issues = ModelAnalyzer.list_of_issues(data, issue_type) + print(io, "\n\n") + ModelAnalyzer.summarize( + io, + issues, + model = model, + verbose = verbose, + max_issues = max_issues, + ) + end + return +end + +function Base.show(io::IO, data::Data) + n = sum( + length(ModelAnalyzer.list_of_issues(data, T)) for + T in ModelAnalyzer.list_of_issue_types(data); + init = 0, + ) + return print(io, "Numerical analysis found $n issues") +end + +# printing helpers + +_print_value(x::Real) = Printf.@sprintf("%1.0e", x) + +function _stringify_bounds(bounds::Vector{Float64}) + lower = bounds[1] < Inf ? _print_value(bounds[1]) : "0e+00" + upper = bounds[2] > -Inf ? _print_value(bounds[2]) : "0e+00" + return string("[", lower, ", ", upper, "]") +end diff --git a/src/feasibility.jl b/src/feasibility.jl deleted file mode 100644 index 654f22f..0000000 --- a/src/feasibility.jl +++ /dev/null @@ -1,1329 +0,0 @@ -# Copyright (c) 2025: Joaquim Garcia, Oscar Dowson and contributors -# -# Use of this source code is governed by an MIT-style license that can be found -# in the LICENSE.md file or at https://opensource.org/licenses/MIT. - -module Feasibility - -import ModelAnalyzer -import Dualization -import MathOptInterface as MOI -import Printf - -""" - Analyzer() <: ModelAnalyzer.AbstractAnalyzer - -The `Analyzer` type is used to perform feasibility analysis on a model. - -## Example - -```julia -julia> data = ModelAnalyzer.analyze( - ModelAnalyzer.Feasibility.Analyzer(), - 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, -); -``` - -The additional parameters: -- `primal_point`: The primal solution point to use for feasibility checking. - If `nothing`, it will use the current primal solution from optimized model. -- `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 and dual objective - checks. -""" -struct Analyzer <: ModelAnalyzer.AbstractAnalyzer end - -""" - AbstractFeasibilityIssue <: AbstractNumericalIssue - -Abstract type for feasibility issues found during the analysis of a model. -""" -abstract type AbstractFeasibilityIssue <: ModelAnalyzer.AbstractIssue end - -""" - PrimalViolation <: AbstractFeasibilityIssue - -The `PrimalViolation` issue is identified when a primal constraint has a -left-hand-side value that is not within the constraint's set. - -For more information, run: -```julia -julia> ModelAnalyzer.summarize(ModelAnalyzer.Feasibility.PrimalViolation) -``` -""" -struct PrimalViolation <: AbstractFeasibilityIssue - ref::MOI.ConstraintIndex - violation::Float64 -end - -ModelAnalyzer.constraint(issue::PrimalViolation) = issue.ref - -ModelAnalyzer.value(issue::PrimalViolation) = issue.violation - -""" - DualConstraintViolation <: AbstractFeasibilityIssue - -The `DualConstraintViolation` issue is identified when a dual constraint has a -value that is not within the dual constraint's set. -This dual constraint corresponds to a primal variable. - -For more information, run: -```julia -julia> ModelAnalyzer.summarize(ModelAnalyzer.Feasibility.DualConstraintViolation) -``` -""" -struct DualConstraintViolation <: AbstractFeasibilityIssue - ref::MOI.VariableIndex - violation::Float64 -end - -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.DualConstrainedVariableViolation) -``` -""" -struct DualConstrainedVariableViolation <: AbstractFeasibilityIssue - ref::MOI.ConstraintIndex - violation::Float64 -end - -ModelAnalyzer.constraint(issue::DualConstrainedVariableViolation) = issue.ref - -ModelAnalyzer.value(issue::DualConstrainedVariableViolation) = issue.violation - -""" - ComplemetarityViolation <: AbstractFeasibilityIssue - -The `ComplemetarityViolation` issue is identified when a pair of primal -constraint and dual variable has a nonzero complementarity value, i.e., the -inner product of the primal constraint's slack and the dual variable's -violation is not zero. - -For more information, run: -```julia -julia> ModelAnalyzer.summarize(ModelAnalyzer.Feasibility.ComplemetarityViolation) -``` -""" -struct ComplemetarityViolation <: AbstractFeasibilityIssue - ref::MOI.ConstraintIndex - violation::Float64 -end - -ModelAnalyzer.constraint(issue::ComplemetarityViolation) = issue.ref - -ModelAnalyzer.value(issue::ComplemetarityViolation) = issue.violation - -""" - DualObjectiveMismatch <: AbstractFeasibilityIssue - -The `DualObjectiveMismatch` issue is identified when the dual objective value -computed from problem data and the dual solution does not match the solver's -dual objective value. - -For more information, run: -```julia -julia> ModelAnalyzer.summarize(ModelAnalyzer.Feasibility.DualObjectiveMismatch) -``` -""" -struct DualObjectiveMismatch <: AbstractFeasibilityIssue - obj::Float64 - obj_solver::Float64 -end - -# ModelAnalyzer.values(issue::DualObjectiveMismatch) = [issue.obj, issue.obj_solver] - -""" - PrimalObjectiveMismatch <: AbstractFeasibilityIssue - -The `PrimalObjectiveMismatch` issue is identified when the primal objective -value computed from problem data and the primal solution does not match -the solver's primal objective value. - -For more information, run: -```julia -julia> ModelAnalyzer.summarize(ModelAnalyzer.Feasibility.PrimalObjectiveMismatch) -``` -""" -struct PrimalObjectiveMismatch <: AbstractFeasibilityIssue - obj::Float64 - obj_solver::Float64 -end - -# ModelAnalyzer.values(issue::PrimalObjectiveMismatch) = [issue.obj, issue.obj_solver] - -""" - PrimalDualMismatch <: AbstractFeasibilityIssue - -The `PrimalDualMismatch` issue is identified when the primal objective value -computed from problem data and the primal solution does not match the dual -objective value computed from problem data and the dual solution. - -For more information, run: -```julia -julia> ModelAnalyzer.summarize(ModelAnalyzer.Feasibility.PrimalDualMismatch) -``` -""" -struct PrimalDualMismatch <: AbstractFeasibilityIssue - primal::Float64 - dual::Float64 -end - -ModelAnalyzer.values(issue::PrimalDualMismatch) = [issue.primal, issue.dual] - -""" - PrimalDualSolverMismatch <: AbstractFeasibilityIssue - -The `PrimalDualSolverMismatch` issue is identified when the primal objective -value reported by the solver does not match the dual objective value reported -by the solver. - -For more information, run: -```julia -julia> ModelAnalyzer.summarize(ModelAnalyzer.Feasibility.PrimalDualSolverMismatch) -``` -""" -struct PrimalDualSolverMismatch <: AbstractFeasibilityIssue - primal::Float64 - 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 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. -""" -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{DualConstraintViolation} = DualConstraintViolation[] - dual_convar::Vector{DualConstrainedVariableViolation} = - DualConstrainedVariableViolation[] - complementarity::Vector{ComplemetarityViolation} = ComplemetarityViolation[] - # objective analysis - dual_objective_mismatch::Vector{DualObjectiveMismatch} = - DualObjectiveMismatch[] - primal_objective_mismatch::Vector{PrimalObjectiveMismatch} = - PrimalObjectiveMismatch[] - primal_dual_mismatch::Vector{PrimalDualMismatch} = PrimalDualMismatch[] - primal_dual_solver_mismatch::Vector{PrimalDualSolverMismatch} = - PrimalDualSolverMismatch[] -end - -function ModelAnalyzer._summarize(io::IO, ::Type{PrimalViolation}) - return print(io, "# PrimalViolation") -end - -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}) - return print(io, "# ComplemetarityViolation") -end - -function ModelAnalyzer._summarize(io::IO, ::Type{DualObjectiveMismatch}) - return print(io, "# DualObjectiveMismatch") -end - -function ModelAnalyzer._summarize(io::IO, ::Type{PrimalObjectiveMismatch}) - return print(io, "# PrimalObjectiveMismatch") -end - -function ModelAnalyzer._summarize(io::IO, ::Type{PrimalDualMismatch}) - return print(io, "# PrimalDualMismatch") -end - -function ModelAnalyzer._summarize(io::IO, ::Type{PrimalDualSolverMismatch}) - return print(io, "# PrimalDualSolverMismatch") -end - -function ModelAnalyzer._verbose_summarize(io::IO, ::Type{PrimalViolation}) - return print( - io, - """ - # PrimalViolation - - ## What - - A `PrimalViolation` issue is identified when a constraint has - function , i.e., a left-hand-side value, that is not within - the constraint's set. - - ## Why - - This can happen due to a few reasons: - - The solver did not converge. - - The model is infeasible and the solver converged to an - infeasible point. - - The solver converged to a low accuracy solution, which might - happen due to transformations in the the model presolve or - due to numerical issues. - - ## How to fix - - Check the solver convergence log and the solver status. If the - solver did not converge, you might want to try alternative - solvers or adjust the solver options. If the solver converged - to an infeasible point, you might want to check the model - constraints and bounds. If the solver converged to a low - accuracy solution, you might want to adjust the solver options - or the model presolve. - - ## More information - - No extra information for this issue. - """, - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - ::Type{DualConstraintViolation}, -) - return print( - io, - """ - # DualConstraintViolation - - ## What - - A `DualConstraintViolation` issue is identified when a constraint has - a dual value that is not within the dual constraint's set. - - ## Why - - This can happen due to a few reasons: - - The solver did not converge. - - The model is infeasible and the solver converged to an - infeasible point. - - The solver converged to a low accuracy solution, which might - happen due to transformations in the the model presolve or - due to numerical issues. - - ## How to fix - - Check the solver convergence log and the solver status. If the - solver did not converge, you might want to try alternative - solvers or adjust the solver options. If the solver converged - to an infeasible point, you might want to check the model - constraints and bounds. If the solver converged to a low - accuracy solution, you might want to adjust the solver options - or the model presolve. - - ## More information - - No extra information for this issue. - """, - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - ::Type{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}, -) - return print( - io, - """ - # ComplemetarityViolation - - ## What - - A `ComplemetarityViolation` issue is identified when a pair of - primal constraint and dual varaible has a nonzero - complementarity value, i.e., the inner product of the primal - constraint's slack and the dual variable's violation is - not zero. - - ## 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{DualObjectiveMismatch}) - return print( - io, - """ - # DualObjectiveMismatch - - ## What - - A `DualObjectiveMismatch` issue is identified when the dual - objective value computed from problema data and the dual - solution does not match the solver's dual objective - value. - - ## Why - - This can happen due to: - - The solver performed presolve transformations and the - reported dual objective is reported from the transformed - problem. - - Bad problem numerical conditioning, very large and very - small coefficients might be present in the model. - - ## How to fix - - Check the solver convergence log and the solver status. - Consider reviewing the coefficients of the objective function. - Consider reviewing the options set in the solver. - - ## More information - - No extra information for this issue. - """, - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - ::Type{PrimalObjectiveMismatch}, -) - return print( - io, - """ - # PrimalObjectiveMismatch - - ## What - - A `PrimalObjectiveMismatch` issue is identified when the primal - objective value computed from problema data and the primal - solution does not match the solver's primal objective - value. - - ## Why - - This can happen due to: - - The solver performed presolve transformations and the - reported primal objective is reported from the transformed - problem. - - Bad problem numerical conditioning, very large and very - small coefficients might be present in the model. - - ## How to fix - - Check the solver convergence log and the solver status. - Consider reviewing the coefficients of the objective function. - Consider reviewing the options set in the solver. - - ## More information - - No extra information for this issue. - """, - ) -end - -function ModelAnalyzer._verbose_summarize(io::IO, ::Type{PrimalDualMismatch}) - return print( - io, - """ - # PrimalDualMismatch - - ## What - - A `PrimalDualMismatch` issue is identified when the primal - objective value computed from problema data and the primal - solution does not match the dual objective value computed - from problem data and the dual solution. - - ## Why - - This can happen due to: - - The solver did not converge. - - Bad problem numerical conditioning, very large and very - small coefficients might be present in the model. - - ## How to fix - - Check the solver convergence log and the solver status. - Consider reviewing the coefficients of the model. - - ## More information - - No extra information for this issue. - """, - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - ::Type{PrimalDualSolverMismatch}, -) - return print( - io, - """ - # PrimalDualSolverMismatch - - ## What - - A `PrimalDualSolverMismatch` issue is identified when the primal - objective value reported by the solver does not match the dual - objective value reported by the solver. - - ## Why - - This can happen due to: - - The solver did not converge. - - ## How to fix - - Check the solver convergence log and the solver status. - - ## More information - - No extra information for this issue. - """, - ) -end - -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::DualConstrainedVariableViolation, - model, -) - return print( - io, - ModelAnalyzer._name(issue.ref, model), - " : ", - issue.violation, - ) -end - -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, model) - return ModelAnalyzer._verbose_summarize(io, issue, model) -end - -function ModelAnalyzer._summarize(io::IO, issue::PrimalObjectiveMismatch, model) - return ModelAnalyzer._verbose_summarize(io, issue, model) -end - -function ModelAnalyzer._summarize(io::IO, issue::PrimalDualMismatch, model) - return ModelAnalyzer._verbose_summarize(io, issue, model) -end - -function ModelAnalyzer._summarize( - io::IO, - issue::PrimalDualSolverMismatch, - model, -) - return ModelAnalyzer._verbose_summarize(io, issue, model) -end - -function ModelAnalyzer._verbose_summarize(io::IO, issue::PrimalViolation, model) - return print( - io, - "Constraint ", - ModelAnalyzer._name(issue.ref, model), - " has primal 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 ", - ModelAnalyzer._name(issue.ref, model), - " has complementarty violation ", - issue.violation, - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - issue::DualObjectiveMismatch, - model, -) - return print( - io, - "Dual objective mismatch: ", - issue.obj, - " (computed) vs ", - issue.obj_solver, - " (reported by solver)\n", - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - issue::PrimalObjectiveMismatch, - model, -) - return print( - io, - "Primal objective mismatch: ", - issue.obj, - " (computed) vs ", - issue.obj_solver, - " (reported by solver)\n", - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - issue::PrimalDualMismatch, - model, -) - return print( - io, - "Primal dual mismatch: ", - issue.primal, - " (computed primal) vs ", - issue.dual, - " (computed dual)\n", - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - issue::PrimalDualSolverMismatch, - model, -) - return print( - io, - "Solver reported objective mismatch: ", - issue.primal, - " (reported primal) vs ", - issue.dual, - " (reported dual)\n", - ) -end - -function ModelAnalyzer.list_of_issues(data::Data, ::Type{PrimalViolation}) - return data.primal -end - -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}, -) - return data.complementarity -end - -function ModelAnalyzer.list_of_issues(data::Data, ::Type{DualObjectiveMismatch}) - return data.dual_objective_mismatch -end - -function ModelAnalyzer.list_of_issues( - data::Data, - ::Type{PrimalObjectiveMismatch}, -) - return data.primal_objective_mismatch -end - -function ModelAnalyzer.list_of_issues(data::Data, ::Type{PrimalDualMismatch}) - return data.primal_dual_mismatch -end - -function ModelAnalyzer.list_of_issues( - data::Data, - ::Type{PrimalDualSolverMismatch}, -) - return data.primal_dual_solver_mismatch -end - -function ModelAnalyzer.list_of_issue_types(data::Data) - ret = Type[] - for type in ( - PrimalViolation, - DualConstraintViolation, - DualConstrainedVariableViolation, - ComplemetarityViolation, - DualObjectiveMismatch, - PrimalObjectiveMismatch, - PrimalDualMismatch, - PrimalDualSolverMismatch, - ) - if !isempty(ModelAnalyzer.list_of_issues(data, type)) - push!(ret, type) - end - end - return ret -end - -function summarize_configurations(io::IO, data::Data) - print(io, "## Configuration\n\n") - # print(io, " - point: ", data.point, "\n") - print(io, " atol: ", data.atol, "\n") - print(io, " skip_missing: ", data.skip_missing, "\n") - return -end - -function ModelAnalyzer.summarize( - io::IO, - data::Data; - model = nothing, - verbose = true, - max_issues = ModelAnalyzer.DEFAULT_MAX_ISSUES, - configurations = true, -) - print(io, "## Feasibility Analysis\n\n") - if configurations - summarize_configurations(io, data) - print(io, "\n") - end - # add maximum primal, dual and compl - # add sum of primal, dual and compl - for issue_type in ModelAnalyzer.list_of_issue_types(data) - issues = ModelAnalyzer.list_of_issues(data, issue_type) - print(io, "\n\n") - ModelAnalyzer.summarize( - io, - issues, - model = model, - verbose = verbose, - max_issues = max_issues, - ) - end - return -end - -function Base.show(io::IO, data::Data) - n = sum( - length(ModelAnalyzer.list_of_issues(data, T)) for - T in ModelAnalyzer.list_of_issue_types(data); - init = 0, - ) - return print(io, "Feasibility analysis found $n issues") -end - -function ModelAnalyzer.analyze( - ::Analyzer, - 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 - 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.", - ) - end - data.primal_point = _last_primal_solution(model) - end - - 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.", - ) - end - data.dual_point = _last_dual_solution(model) - end - - _analyze_primal!(model, data) - _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, _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)) - return data -end - -function _analyze_primal!(model, data) - 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 _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) - 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 -end - -# not needed because it would have stoped in dualization before -# function _set_value(set::MOI.AbstractScalarSet) -# return 0.0 -# end -# function _set_value(set::MOI.Interval) -# error("Interval sets are not supported.") -# return (set.lower, set.upper) -# end - -function _set_value(set::MOI.AbstractVectorSet) - return zeros(MOI.dimension(set)) -end - -function _set_value(set::MOI.LessThan) - return set.upper -end - -function _set_value(set::MOI.GreaterThan) - return set.lower -end - -function _set_value(set::MOI.EqualTo) - return set.value -end - -function _analyze_objectives!(model::MOI.ModelLike, dual_model, map, data) - primal_status = MOI.get(model, MOI.PrimalStatus()) - dual_status = MOI.get(model, MOI.DualStatus()) - if 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 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 - - if dual_obj_val_solver !== nothing && - obj_val_solver !== nothing && - !isapprox(obj_val_solver, dual_obj_val_solver; atol = data.atol) - push!( - data.primal_dual_solver_mismatch, - PrimalDualSolverMismatch(obj_val_solver, dual_obj_val_solver), - ) - end - - 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) - push!( - data.dual_objective_mismatch, - DualObjectiveMismatch(dual_obj_val, dual_obj_val_solver), - ) - end - - 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), - ) - end - end - - return -end - -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 _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 _can_dualize(model::MOI.ModelLike) - types = MOI.get(model, MOI.ListOfConstraintTypesPresent()) - - for (F, S) in types - if !Dualization.supported_constraint(F, S) - return false - end - end - - F = MOI.get(model, MOI.ObjectiveFunctionType()) - - if !Dualization.supported_objective(F) - return false - end - - sense = MOI.get(model, MOI.ObjectiveSense()) - if sense == MOI.FEASIBILITY_SENSE - return false - end - - return true -end - -end # module diff --git a/src/numerical.jl b/src/numerical.jl deleted file mode 100644 index 8eb904a..0000000 --- a/src/numerical.jl +++ /dev/null @@ -1,2509 +0,0 @@ -# Copyright (c) 2025: Joaquim Garcia, Oscar Dowson and contributors -# -# Use of this source code is governed by an MIT-style license that can be found -# in the LICENSE.md file or at https://opensource.org/licenses/MIT. - -module Numerical - -import ModelAnalyzer -import LinearAlgebra -import Printf -import MathOptInterface as MOI - -""" - Analyzer() <: ModelAnalyzer.AbstractAnalyzer - -The `Analyzer` type is used to analyze the coefficients of a model for -numerical issues. - -## Example - -```julia -julia> data = ModelAnalyzer.analyze( - ModelAnalyzer.Numerical.Analyzer(), - model; - threshold_dense_fill_in = 0.10, - threshold_dense_entries = 1000, - threshold_small = 1e-5, - threshold_large = 1e+5, -) -``` - -The additional parameters: -- `threshold_dense_fill_in`: The threshold for the fraction of non-zero entries - in a constraint to be considered dense. -- `threshold_dense_entries`: The minimum number of non-zero entries for a - constraint to be considered dense. -- `threshold_small`: The threshold for small coefficients in the model. -- `threshold_large`: The threshold for large coefficients in the model. - -""" -struct Analyzer <: ModelAnalyzer.AbstractAnalyzer end - -""" - AbstractNumericalIssue <: AbstractNumericalIssue - -Abstract type for numerical issues found during the analysis of a model. -""" -abstract type AbstractNumericalIssue <: ModelAnalyzer.AbstractIssue end - -""" - VariableNotInConstraints <: AbstractNumericalIssue - -The `VariableNotInConstraints` issue is identified when a variable appears in no -constraints. - -For more information, run: -```julia -julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.VariableNotInConstraints) -``` -""" -struct VariableNotInConstraints <: AbstractNumericalIssue - ref::MOI.VariableIndex -end - -ModelAnalyzer.variable(issue::VariableNotInConstraints) = issue.ref - -""" - EmptyConstraint <: AbstractNumericalIssue - -The `EmptyConstraint` issue is identified when a constraint has no coefficients -different from zero. - -For more information, run: -```julia -julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.EmptyConstraint) -``` -""" -struct EmptyConstraint <: AbstractNumericalIssue - ref::MOI.ConstraintIndex -end - -ModelAnalyzer.constraint(issue::EmptyConstraint) = issue.ref - -""" - VariableBoundAsConstraint <: AbstractNumericalIssue - -The `VariableBoundAsConstraint` issue is identified when a constraint is -equivalent to a variable bound, that is, the constraint has only one non-zero -coefficient, and this coefficient is equal to one. - -For more information, run: -```julia -julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.VariableBoundAsConstraint) -``` -""" -struct VariableBoundAsConstraint <: AbstractNumericalIssue - ref::MOI.ConstraintIndex -end - -ModelAnalyzer.constraint(issue::VariableBoundAsConstraint) = issue.ref - -""" - DenseConstraint <: AbstractNumericalIssue - -The `DenseConstraint` issue is identified when a constraint has a fraction of -non-zero entries greater than `threshold_dense_fill_in` and the number of -non-zero entries is greater than `threshold_dense_entries`. - -For more information, run: -```julia -julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.DenseConstraint) -``` -""" -struct DenseConstraint <: AbstractNumericalIssue - ref::MOI.ConstraintIndex - nnz::Int -end - -ModelAnalyzer.constraint(issue::DenseConstraint) = issue.ref - -ModelAnalyzer.value(issue::DenseConstraint) = issue.nnz - -""" - SmallMatrixCoefficient <: AbstractNumericalIssue - -The `SmallMatrixCoefficient` issue is identified when a matrix coefficient in a -constraint is smaller than `threshold_small`. - -For more information, run: -```julia -julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.SmallMatrixCoefficient) -``` -""" -struct SmallMatrixCoefficient <: AbstractNumericalIssue - 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 - -The `LargeMatrixCoefficient` issue is identified when a matrix coefficient in a -constraint is larger than `threshold_large`. - -For more information, run: -```julia -julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.LargeMatrixCoefficient) -``` -""" -struct LargeMatrixCoefficient <: AbstractNumericalIssue - 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 - -The `SmallBoundCoefficient` issue is identified when a variable's bound -(coefficient) is smaller than `threshold_small`. - -For more information, run: -```julia -julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.SmallBoundCoefficient) -``` -""" -struct SmallBoundCoefficient <: AbstractNumericalIssue - variable::MOI.VariableIndex - coefficient::Float64 -end - -ModelAnalyzer.variable(issue::SmallBoundCoefficient) = issue.variable - -ModelAnalyzer.value(issue::SmallBoundCoefficient) = issue.coefficient - -""" - LargeBoundCoefficient <: AbstractNumericalIssue - -The `LargeBoundCoefficient` issue is identified when a variable's bound -(coefficient) is larger than `threshold_large`. - -For more information, run: -```julia -julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.LargeBoundCoefficient) -``` -""" -struct LargeBoundCoefficient <: AbstractNumericalIssue - variable::MOI.VariableIndex - coefficient::Float64 -end - -ModelAnalyzer.variable(issue::LargeBoundCoefficient) = issue.variable - -ModelAnalyzer.value(issue::LargeBoundCoefficient) = issue.coefficient - -""" - SmallRHSCoefficient <: AbstractNumericalIssue - -The `SmallRHSCoefficient` issue is identified when the right-hand-side (RHS) -coefficient of a constraint is smaller than `threshold_small`. - -For more information, run: -```julia -julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.SmallRHSCoefficient) -``` -""" -struct SmallRHSCoefficient <: AbstractNumericalIssue - ref::MOI.ConstraintIndex - coefficient::Float64 -end - -ModelAnalyzer.constraint(issue::SmallRHSCoefficient) = issue.ref - -ModelAnalyzer.value(issue::SmallRHSCoefficient) = issue.coefficient - -""" - LargeRHSCoefficient <: AbstractNumericalIssue - -The `LargeRHSCoefficient` issue is identified when the right-hand-side (RHS) -coefficient of a constraint is larger than `threshold_large`. - -For more information, run: -```julia -julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.LargeRHSCoefficient) -``` -""" -struct LargeRHSCoefficient <: AbstractNumericalIssue - ref::MOI.ConstraintIndex - coefficient::Float64 -end - -ModelAnalyzer.constraint(issue::LargeRHSCoefficient) = issue.ref - -ModelAnalyzer.value(issue::LargeRHSCoefficient) = issue.coefficient - -""" - SmallObjectiveCoefficient <: AbstractNumericalIssue - -The `SmallObjectiveCoefficient` issue is identified when a coefficient in the -objective function is smaller than `threshold_small`. - -For more information, run: -```julia -julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.SmallObjectiveCoefficient) -``` -""" -struct SmallObjectiveCoefficient <: AbstractNumericalIssue - variable::MOI.VariableIndex - coefficient::Float64 -end - -ModelAnalyzer.variable(issue::SmallObjectiveCoefficient) = issue.variable - -ModelAnalyzer.value(issue::SmallObjectiveCoefficient) = issue.coefficient - -""" - LargeObjectiveCoefficient <: AbstractNumericalIssue - -The `LargeObjectiveCoefficient` issue is identified when a coefficient in the -objective function is larger than `threshold_large`. - -For more information, run: -```julia -julia> ModelAnalyzer.summarize(ModelAnalyzer.Numerical.LargeObjectiveCoefficient) -``` -""" -struct LargeObjectiveCoefficient <: AbstractNumericalIssue - variable::MOI.VariableIndex - coefficient::Float64 -end - -ModelAnalyzer.variable(issue::LargeObjectiveCoefficient) = issue.variable - -ModelAnalyzer.value(issue::LargeObjectiveCoefficient) = issue.coefficient - -""" - SmallObjectiveQuadraticCoefficient <: AbstractNumericalIssue - -The `SmallObjectiveQuadraticCoefficient` issue is identified when a quadratic -coefficient in the objective function is smaller than `threshold_small`. - -For more information, run: -```julia -julia> ModelAnalyzer.summarize( - ModelAnalyzer.Numerical.SmallObjectiveQuadraticCoefficient -) -``` -""" -struct SmallObjectiveQuadraticCoefficient <: AbstractNumericalIssue - 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 - -The `LargeObjectiveQuadraticCoefficient` issue is identified when a quadratic -coefficient in the objective function is larger than `threshold_large`. - -For more information, run: -```julia -julia> ModelAnalyzer.summarize( - ModelAnalyzer.Numerical.LargeObjectiveQuadraticCoefficient -) -``` -""" -struct LargeObjectiveQuadraticCoefficient <: AbstractNumericalIssue - 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 - -The `SmallMatrixQuadraticCoefficient` issue is identified when a quadratic -coefficient in a constraint is smaller than `threshold_small`. - -For more information, run: -```julia -julia> ModelAnalyzer.summarize( - ModelAnalyzer.Numerical.SmallMatrixQuadraticCoefficient -) -``` -""" -struct SmallMatrixQuadraticCoefficient <: AbstractNumericalIssue - 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 - -The `LargeMatrixQuadraticCoefficient` issue is identified when a quadratic -coefficient in a constraint is larger than `threshold_large`. - -For more information, run: -```julia -julia> ModelAnalyzer.summarize( - ModelAnalyzer.Numerical.LargeMatrixQuadraticCoefficient -) -``` -""" -struct LargeMatrixQuadraticCoefficient <: AbstractNumericalIssue - 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 - -The `NonconvexQuadraticObjective` issue is identified when a quadratic -objective function is non-convex. - -For more information, run: -```julia -julia> ModelAnalyzer.summarize( - ModelAnalyzer.Numerical.NonconvexQuadraticObjective -) -``` -""" -struct NonconvexQuadraticObjective <: AbstractNumericalIssue end - -""" - NonconvexQuadraticConstraint - -The `NonconvexQuadraticConstraint` issue is identified when a quadratic -constraint is non-convex. - -For more information, run: -```julia -julia> ModelAnalyzer.summarize( - ModelAnalyzer.Numerical.NonconvexQuadraticConstraint -) -``` -""" -struct NonconvexQuadraticConstraint <: AbstractNumericalIssue - ref::MOI.ConstraintIndex -end -ModelAnalyzer.constraint(issue::NonconvexQuadraticConstraint) = issue.ref - -""" - Data - -The `Data` structure holds the results of the analysis performed by the -`ModelAnalyzer.Numerical.Analyzer`. It contains various thresholds and the -information about the model's variables, constraints, and objective function. -""" -Base.@kwdef mutable struct Data <: ModelAnalyzer.AbstractData - # analysis configuration - threshold_dense_fill_in::Float64 = 0.10 - threshold_dense_entries::Int = 1000 - threshold_small::Float64 = 1e-5 - threshold_large::Float64 = 1e+5 - # main numbers - number_of_variables::Int = 0 - number_of_constraints::Int = 0 - constraint_info::Vector{Tuple{DataType,DataType,Int}} = - Tuple{DataType,DataType,Int}[] - # objective_info::Any - matrix_nnz::Int = 0 - # ranges - matrix_range::Vector{Float64} = sizehint!(Float64[1.0, 1.0], 2) - bounds_range::Vector{Float64} = sizehint!(Float64[1.0, 1.0], 2) - 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{MOI.VariableIndex} = Set{MOI.VariableIndex}() - # variables analysis - variables_not_in_constraints::Vector{VariableNotInConstraints} = - VariableNotInConstraints[] - bounds_small::Vector{SmallBoundCoefficient} = SmallBoundCoefficient[] - bounds_large::Vector{LargeBoundCoefficient} = LargeBoundCoefficient[] - # constraints analysis - empty_rows::Vector{EmptyConstraint} = EmptyConstraint[] - bound_rows::Vector{VariableBoundAsConstraint} = VariableBoundAsConstraint[] - dense_rows::Vector{DenseConstraint} = DenseConstraint[] - matrix_small::Vector{SmallMatrixCoefficient} = SmallMatrixCoefficient[] - matrix_large::Vector{LargeMatrixCoefficient} = LargeMatrixCoefficient[] - rhs_small::Vector{SmallRHSCoefficient} = SmallRHSCoefficient[] - rhs_large::Vector{LargeRHSCoefficient} = LargeRHSCoefficient[] - # quadratic constraints analysis - has_quadratic_constraints::Bool = false - nonconvex_rows::Vector{NonconvexQuadraticConstraint} = - NonconvexQuadraticConstraint[] - matrix_quadratic_small::Vector{SmallMatrixQuadraticCoefficient} = - SmallMatrixQuadraticCoefficient[] - matrix_quadratic_large::Vector{LargeMatrixQuadraticCoefficient} = - LargeMatrixQuadraticCoefficient[] - # cache data - sense::MOI.OptimizationSense = MOI.FEASIBILITY_SENSE - # objective analysis - objective_small::Vector{SmallObjectiveCoefficient} = - SmallObjectiveCoefficient[] - objective_large::Vector{LargeObjectiveCoefficient} = - LargeObjectiveCoefficient[] - # quadratic objective analysis - has_quadratic_objective::Bool = false - objective_quadratic_range::Vector{Float64} = sizehint!(Float64[1.0, 1.0], 2) - matrix_quadratic_range::Vector{Float64} = sizehint!(Float64[1.0, 1.0], 2) - nonconvex_objective::Vector{NonconvexQuadraticObjective} = - NonconvexQuadraticObjective[] - objective_quadratic_small::Vector{SmallObjectiveQuadraticCoefficient} = - SmallObjectiveQuadraticCoefficient[] - objective_quadratic_large::Vector{LargeObjectiveQuadraticCoefficient} = - LargeObjectiveQuadraticCoefficient[] -end - -function _update_range(range::Vector{Float64}, value::Number) - range[1] = min(range[1], abs(value)) - range[2] = max(range[2], abs(value)) - return 1 -end - -function _get_objective_data(data, func::MOI.VariableIndex) - return -end - -function _get_objective_data(data, func::MOI.ScalarAffineFunction) - nnz = 0 - for term in func.terms - variable = term.variable - coefficient = term.coefficient - if iszero(coefficient) - continue - end - nnz += _update_range(data.objective_range, coefficient) - if abs(coefficient) < data.threshold_small - push!( - data.objective_small, - SmallObjectiveCoefficient(variable, coefficient), - ) - elseif abs(coefficient) > data.threshold_large - push!( - data.objective_large, - LargeObjectiveCoefficient(variable, coefficient), - ) - end - end - return -end - -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 term in func.quadratic_terms - coefficient = term.coefficient - v1 = term.variable_1 - v2 = term.variable_2 - if iszero(coefficient) - continue - end - nnz += _update_range(data.objective_quadratic_range, coefficient) - if abs(coefficient) < data.threshold_small - push!( - data.objective_quadratic_small, - SmallObjectiveQuadraticCoefficient(v1, v2, coefficient), - ) - elseif abs(coefficient) > data.threshold_large - push!( - data.objective_quadratic_large, - LargeObjectiveQuadraticCoefficient(v1, v2, coefficient), - ) - end - end - data.has_quadratic_objective = true - if data.sense == MOI.MAX_SENSE - if !_quadratic_vexity(func, -1) - push!(data.nonconvex_objective, NonconvexQuadraticObjective()) - end - elseif data.sense == MOI.MIN_SENSE - if !_quadratic_vexity(func, 1) - push!(data.nonconvex_objective, NonconvexQuadraticObjective()) - end - end - return -end - -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{MOI.VariableIndex,Int}() - for (idx, var) in enumerate(variables) - var_map[var] = idx - end - matrix = zeros(length(variables), length(variables)) - 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), - LinearAlgebra.RowMaximum(), - check = false, - ) - 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::MOI.ConstraintIndex, - func::MOI.ScalarAffineFunction; - ignore_extras = false, -) - if length(func.terms) == 1 - 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 - return - end - end - nnz = 0 - for term in func.terms - variable = term.variable - coefficient = term.coefficient - if iszero(coefficient) - continue - end - nnz += _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 - if nnz == 0 - if !ignore_extras - push!(data.empty_rows, EmptyConstraint(ref)) - end - return - end - if nnz / data.number_of_variables > data.threshold_dense_fill_in && - nnz > data.threshold_dense_entries - push!(data.dense_rows, DenseConstraint(ref, nnz)) - end - data.matrix_nnz += nnz - return -end - -function _get_constraint_matrix_data( - data, - ref::MOI.ConstraintIndex, - func::MOI.ScalarQuadraticFunction{T}, -) where {T} - nnz = 0 - for term in func.quadratic_terms - v1 = term.variable_1 - v2 = term.variable_2 - coefficient = term.coefficient - if iszero(coefficient) - continue - end - nnz += _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 - data.has_quadratic_constraints = true - _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::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 -end - -function _get_constraint_matrix_data( - data, - 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 - _get_constraint_matrix_data( - data, - ref, - MOI.VectorAffineFunction{T}(func.affine_terms, func.constants), - # ignore_extras = nnz > 0, - ) - return -end - -function _get_constraint_matrix_data( - data, - ref::MOI.ConstraintIndex, - func::MOI.VariableIndex, -) - # push!(data.variables_in_constraints, func) - return -end - -function _get_constraint_matrix_data( - data, - ref::MOI.ConstraintIndex, - func::MOI.VectorOfVariables, -) - if length(func.variables) == 1 - return - end - for var in func.variables - push!(data.variables_in_constraints, var) - end - return -end - -function _get_constraint_data( - data, - ref, - func::Union{MOI.ScalarAffineFunction,MOI.ScalarQuadraticFunction}, - set, -) - coefficient = func.constant - if iszero(coefficient) - return - 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 - return -end - -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::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 - return -end - -function _get_constraint_data( - data, - ref, - func::MOI.ScalarAffineFunction, - set::MOI.LessThan, -) - coefficient = set.upper - func.constant - if iszero(coefficient) - return - 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 - return -end - -function _get_constraint_data( - data, - ref, - 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 - return -end - -function _get_constraint_data( - data, - ref, - func::MOI.ScalarAffineFunction, - set::MOI.GreaterThan, -) - coefficient = set.lower - func.constant - if iszero(coefficient) - return - 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 - return -end - -function _get_constraint_data( - data, - ref, - func::MOI.ScalarQuadraticFunction, - set::Union{MOI.EqualTo,MOI.Interval}, -) - _get_constraint_data( - data, - ref, - MOI.ScalarAffineFunction(func.affine_terms, func.constant), - set, - ) - push!(data.nonconvex_rows, NonconvexQuadraticConstraint(ref)) - return -end - -function _get_constraint_data( - data, - ref, - func::MOI.VectorQuadraticFunction, - set::MOI.Zeros, -) - _get_constraint_data( - data, - ref, - MOI.VectorAffineFunction(func.affine_terms, func.constants), - set, - ) - push!(data.nonconvex_rows, NonconvexQuadraticConstraint(ref)) - return -end - -function _get_constraint_data( - data, - ref, - func::MOI.ScalarAffineFunction, - set::MOI.EqualTo, -) - coefficient = set.value - func.constant - if iszero(coefficient) - return - 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 - return -end - -function _get_constraint_data( - data, - ref, - func::MOI.ScalarAffineFunction, - set::MOI.Interval, -) - coefficient = set.upper - func.constant - if !(iszero(coefficient)) - _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 - coefficient = set.lower - func.constant - if iszero(coefficient) - return - 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 - return -end - -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 model. - -""" -function ModelAnalyzer.analyze( - ::Analyzer, - model::MOI.ModelLike, - ; - threshold_dense_fill_in::Float64 = 0.10, - threshold_dense_entries::Int = 1000, - threshold_small::Float64 = 1e-5, - threshold_large::Float64 = 1e+5, -) - data = Data() - data.threshold_dense_fill_in = threshold_dense_fill_in - data.threshold_dense_entries = threshold_dense_entries - data.threshold_small = threshold_small - data.threshold_large = threshold_large - - # initialize simples data - 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) - - # objective pass - objective_type = MOI.get(model, MOI.ObjectiveFunctionType()) - obj_func = MOI.get(model, MOI.ObjectiveFunction{objective_type}()) - _get_objective_data(data, obj_func) - - # constraints pass - 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 - 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 - # 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, - VariableNotInConstraints(var), - ) - end - end - sort!(data.dense_rows, by = x -> x.nnz, rev = true) - sort!(data.matrix_small, by = x -> abs(x.coefficient)) - sort!(data.matrix_large, by = x -> abs(x.coefficient), rev = true) - sort!(data.bounds_small, by = x -> abs(x.coefficient)) - sort!(data.bounds_large, by = x -> abs(x.coefficient), rev = true) - sort!(data.rhs_small, by = x -> abs(x.coefficient)) - sort!(data.rhs_large, by = x -> abs(x.coefficient), rev = true) - sort!(data.matrix_quadratic_small, by = x -> abs(x.coefficient)) - sort!(data.matrix_quadratic_large, by = x -> abs(x.coefficient), rev = true) - # objective - sort!(data.objective_small, by = x -> abs(x.coefficient)) - sort!(data.objective_large, by = x -> abs(x.coefficient), rev = true) - sort!(data.objective_quadratic_small, by = x -> abs(x.coefficient)) - sort!( - data.objective_quadratic_large, - by = x -> abs(x.coefficient), - rev = true, - ) - return data -end - -# API - -function ModelAnalyzer._summarize(io::IO, ::Type{VariableNotInConstraints}) - return print(io, "# VariableNotInConstraints") -end - -function ModelAnalyzer._summarize(io::IO, ::Type{EmptyConstraint}) - return print(io, "# EmptyConstraint") -end - -function ModelAnalyzer._summarize(io::IO, ::Type{VariableBoundAsConstraint}) - return print(io, "# VariableBoundAsConstraint") -end - -function ModelAnalyzer._summarize(io::IO, ::Type{DenseConstraint}) - return print(io, "# DenseConstraint") -end - -function ModelAnalyzer._summarize(io::IO, ::Type{SmallMatrixCoefficient}) - return print(io, "# SmallMatrixCoefficient") -end - -function ModelAnalyzer._summarize(io::IO, ::Type{LargeMatrixCoefficient}) - return print(io, "# LargeMatrixCoefficient") -end - -function ModelAnalyzer._summarize(io::IO, ::Type{SmallBoundCoefficient}) - return print(io, "# SmallBoundCoefficient") -end - -function ModelAnalyzer._summarize(io::IO, ::Type{LargeBoundCoefficient}) - return print(io, "# LargeBoundCoefficient") -end - -function ModelAnalyzer._summarize(io::IO, ::Type{SmallRHSCoefficient}) - return print(io, "# SmallRHSCoefficient") -end - -function ModelAnalyzer._summarize(io::IO, ::Type{LargeRHSCoefficient}) - return print(io, "# LargeRHSCoefficient") -end - -function ModelAnalyzer._summarize(io::IO, ::Type{SmallObjectiveCoefficient}) - return print(io, "# SmallObjectiveCoefficient") -end - -function ModelAnalyzer._summarize(io::IO, ::Type{LargeObjectiveCoefficient}) - return print(io, "# LargeObjectiveCoefficient") -end - -function ModelAnalyzer._summarize( - io::IO, - ::Type{SmallObjectiveQuadraticCoefficient}, -) - return print(io, "# SmallObjectiveQuadraticCoefficient") -end - -function ModelAnalyzer._summarize( - io::IO, - ::Type{LargeObjectiveQuadraticCoefficient}, -) - return print(io, "# LargeObjectiveQuadraticCoefficient") -end - -function ModelAnalyzer._summarize( - io::IO, - ::Type{SmallMatrixQuadraticCoefficient}, -) - return print(io, "# SmallMatrixQuadraticCoefficient") -end - -function ModelAnalyzer._summarize( - io::IO, - ::Type{LargeMatrixQuadraticCoefficient}, -) - return print(io, "# LargeMatrixQuadraticCoefficient") -end - -function ModelAnalyzer._summarize(io::IO, ::Type{NonconvexQuadraticObjective}) - return print(io, "# NonconvexQuadraticObjective") -end - -function ModelAnalyzer._summarize(io::IO, ::Type{NonconvexQuadraticConstraint}) - return print(io, "# NonconvexQuadraticConstraint") -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - ::Type{VariableNotInConstraints}, -) - return print( - io, - """ - # `VariableNotInConstraints` - - ## What - - A `VariableNotInConstraints` issue is identified when a variable appears - in no constraints. If a variable only appears alone in a constraint and - it has a coefficient of 1 it is considered a - `VariableNotInConstraints`, because this emulates a bound. - - ## Why - - This can be a sign of a mistake in the model formulation. If a variable - is not used in any constraints, it is not affecting the solution of the - problem. Moreover, it might be leading to an unbounded problem. - - ## How to fix - - If the variable is not needed, remove it from the model. If the variable - is needed, check that it is correctly used in the constraints. - - ## More information - - No extra information for this issue. - """, - ) -end - -function ModelAnalyzer._verbose_summarize(io::IO, ::Type{EmptyConstraint}) - return print( - io, - """ - # `EmptyConstraint` - - ## What - - An `EmptyConstraint` issue is identified when a constraint has no - coefficients different from zero. - - ## Why - - This can be a sign of a mistake in the model formulation. An empty - constraint is not affecting the solution of the problem. Moreover, it - might be leading to an infeasible problem since the \"left-hand-side\" - of the constraint is always zero. - - ## How to fix - - Remove the empty constraint from the model. - - ## More information - - No extra information for this issue. - """, - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - ::Type{VariableBoundAsConstraint}, -) - return print( - io, - """ - # `VariableBoundAsConstraint` - - ## What - - A `VariableBoundAsConstraint` issue is identified when a constraint is - equivalent to a variable bound, that is, the constraint has only one - non-zero coefficient, and this coefficient is equal to one. - - ## Why - - This can be a sign of a mistake in the model formulation. Variable - bounds are frequently used by solver in special ways that can lead to - better performance. - - ## How to fix - - Remove the constraint and use the variable bound directly. - - ## More information - - - https://support.gurobi.com/hc/en-us/community/posts/24066470832529/comments/24183896218385 - """, - ) -end - -function ModelAnalyzer._verbose_summarize(io::IO, ::Type{DenseConstraint}) - return print( - io, - """ - # `DenseConstraint` - - ## What - - A `DenseConstraint` issue is identified when a constraint has a high - number of non-zero coefficients. - - ## Why - - Dense constraints can lead to performance issues in the solution - process. Very few dense constraints might not be a problem. - - ## How to fix - - Check if the constraint can be simplified. A common - case that can be avoided is when there is an expression - `e = c1 * x1 + c2 * x2 + ... + cn * xn` where `c1, c2, ..., cn` are - constants and `x1, x2, ..., xn` are variables, and this expression is - used in many constraints. In this case, it is better to create a new - variable `y = e` and use `y` in the constraints. - - ## More information - - No extra information for this issue. - """, - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - ::Type{SmallMatrixCoefficient}, -) - return print( - io, - """ - # `SmallMatrixCoefficient` - - ## What - - A `SmallMatrixCoefficient` issue is identified when a constraint has a - coefficient with a small absolute value. - - ## Why - - Small coefficients can lead to numerical instability in the solution - process. - - ## How to fix - - Check if the coefficient is correct. Check if the units of variables and - coefficients are correct. Check if the number makes is - reasonable given that solver have tolerances. Sometimes these - coefficients can be replaced by zeros. - - ## More information - - - https://jump.dev/JuMP.jl/stable/tutorials/getting_started/tolerances/ - """, - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - ::Type{LargeMatrixCoefficient}, -) - return print( - io, - """ - # `LargeMatrixCoefficient` - - ## What - - A `LargeMatrixCoefficient` issue is identified when a constraint has a - coefficient with a large absolute value. - - ## Why - - Large coefficients can lead to numerical instability in the solution - process. - - ## How to fix - - Check if the coefficient is correct. Check if the units of variables and - coefficients are correct. Check if the number makes is - reasonable given that solver have tolerances. Sometimes these - coefficients can be replaced by zeros. - - ## More information - - - https://jump.dev/JuMP.jl/stable/tutorials/getting_started/tolerances/ - """, - ) -end - -function ModelAnalyzer._verbose_summarize(io::IO, ::Type{SmallBoundCoefficient}) - return print( - io, - """ - # `SmallBoundCoefficient` - - ## What - - A `SmallBoundCoefficient` issue is identified when a variable has a - bound with a small absolute value. - - ## Why - - Small bounds can lead to numerical instability in the solution process. - - ## How to fix - - Check if the bound is correct. Check if the units of variables and - coefficients are correct. Check if the number makes is - reasonable given that solver have tolerances. Sometimes these - bounds can be replaced by zeros. - - ## More information - - - https://jump.dev/JuMP.jl/stable/tutorials/getting_started/tolerances/ - """, - ) -end - -function ModelAnalyzer._verbose_summarize(io::IO, ::Type{LargeBoundCoefficient}) - return print( - io, - """ - # `LargeBoundCoefficient` - - ## What - - A `LargeBoundCoefficient` issue is identified when a variable has a - bound with a large absolute value. - - ## Why - - Large bounds can lead to numerical instability in the solution process. - - ## How to fix - - Check if the bound is correct. Check if the units of variables and - coefficients are correct. Check if the number makes is - reasonable given that solver have tolerances. Sometimes these - bounds can be replaced by zeros. - - ## More information - - - https://jump.dev/JuMP.jl/stable/tutorials/getting_started/tolerances/ - """, - ) -end - -function ModelAnalyzer._verbose_summarize(io::IO, ::Type{SmallRHSCoefficient}) - return print( - io, - """ - # `SmallRHSCoefficient` - - ## What - - A `SmallRHSCoefficient` issue is identified when a constraint has a - right-hand-side with a small absolute value. - - ## Why - - Small right-hand-sides can lead to numerical instability in the solution - process. - - ## How to fix - - Check if the right-hand-side is correct. Check if the units of variables - and coefficients are correct. Check if the number makes is - reasonable given that solver have tolerances. Sometimes these - right-hand-sides can be replaced by zeros. - - ## More information - - - https://jump.dev/JuMP.jl/stable/tutorials/getting_started/tolerances/ - """, - ) -end - -function ModelAnalyzer._verbose_summarize(io::IO, ::Type{LargeRHSCoefficient}) - return print( - io, - """ - # `LargeRHSCoefficient` - - ## What - - A `LargeRHSCoefficient` issue is identified when a constraint has a - right-hand-side with a large absolute value. - - ## Why - - Large right-hand-sides can lead to numerical instability in the solution - process. - - ## How to fix - - Check if the right-hand-side is correct. Check if the units of variables - and coefficients are correct. Check if the number makes is - reasonable given that solver have tolerances. Sometimes these - right-hand-sides can be replaced by zeros. - - ## More information - - - https://jump.dev/JuMP.jl/stable/tutorials/getting_started/tolerances/ - """, - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - ::Type{SmallObjectiveCoefficient}, -) - return print( - io, - """ - # `SmallObjectiveCoefficient` - - ## What - - A `SmallObjectiveCoefficient` issue is identified when the objective - function has a coefficient with a small absolute value. - - ## Why - - Small coefficients can lead to numerical instability in the solution - process. - - ## How to fix - - Check if the coefficient is correct. Check if the units of variables and - coefficients are correct. Check if the number makes is - reasonable given that solver have tolerances. Sometimes these - coefficients can be replaced by zeros. - - ## More information - - - https://jump.dev/JuMP.jl/stable/tutorials/getting_started/tolerances/ - """, - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - ::Type{LargeObjectiveCoefficient}, -) - return print( - io, - """ - # `LargeObjectiveCoefficient` - - ## What - - A `LargeObjectiveCoefficient` issue is identified when the objective - function has a coefficient with a large absolute value. - - ## Why - - Large coefficients can lead to numerical instability in the solution - process. - - ## How to fix - - Check if the coefficient is correct. Check if the units of variables and - coefficients are correct. Check if the number makes is - reasonable given that solver have tolerances. Sometimes these - coefficients can be replaced by zeros. - - ## More information - - - https://jump.dev/JuMP.jl/stable/tutorials/getting_started/tolerances/ - """, - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - ::Type{SmallObjectiveQuadraticCoefficient}, -) - return print( - io, - """ - # `SmallObjectiveQuadraticCoefficient` - - ## What - - A `SmallObjectiveQuadraticCoefficient` issue is identified when the - objective function has a quadratic coefficient with a small absolute value. - - ## Why - - Small coefficients can lead to numerical instability in the solution - process. - - ## How to fix - - Check if the coefficient is correct. Check if the units of variables and - coefficients are correct. Check if the number makes is - reasonable given that solver have tolerances. Sometimes these - coefficients can be replaced by zeros. - - ## More information - - - https://jump.dev/JuMP.jl/stable/tutorials/getting_started/tolerances/ - """, - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - ::Type{LargeObjectiveQuadraticCoefficient}, -) - return print( - io, - """ - # `LargeObjectiveQuadraticCoefficient` - - ## What - - A `LargeObjectiveQuadraticCoefficient` issue is identified when the - objective function has a quadratic coefficient with a large absolute value. - - ## Why - - Large coefficients can lead to numerical instability in the solution - process. - - ## How to fix - - Check if the coefficient is correct. Check if the units of variables and - coefficients are correct. Check if the number makes is - reasonable given that solver have tolerances. Sometimes these - coefficients can be replaced by zeros. - - ## More information - - - https://jump.dev/JuMP.jl/stable/tutorials/getting_started/tolerances/ - """, - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - ::Type{SmallMatrixQuadraticCoefficient}, -) - return print( - io, - """ - # `SmallMatrixQuadraticCoefficient` - - ## What - - A `SmallMatrixQuadraticCoefficient` issue is identified when a quadratic - constraint has a coefficient with a small absolute value. - - ## Why - - Small coefficients can lead to numerical instability in the solution - process. - - ## How to fix - - Check if the coefficient is correct. Check if the units of variables and - coefficients are correct. Check if the number makes is - reasonable given that solver have tolerances. Sometimes these - coefficients can be replaced by zeros. - - ## More information - - - https://jump.dev/JuMP.jl/stable/tutorials/getting_started/tolerances/ - """, - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - ::Type{NonconvexQuadraticObjective}, -) - return print( - io, - """ - # `NonconvexQuadraticObjective` - - ## What - - A `NonconvexQuadraticObjective` issue is identified when a quadratic - objective is nonconvex, that is, the quadratic matrix is not positive - semidefinite for minimization or the quadratic matrix is not negative - semidefinite for maximization. - - ## Why - - Nonconvex objectives are not expected by many solver and can lead to - wrong solutions or even convergence issues. - - ## How to fix - - Check if the objective is correct. Coefficient signs might have been - inverted. This also occurs if user fix a variable to emulate a - parameter, in this case some solvers will not be able to solve the - model properly, other tools such as ParametricOptInteface.jl might be - more suitable than fixing variables. - - ## More information - - No extra information for this issue. - """, - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - ::Type{NonconvexQuadraticConstraint}, -) - return print( - io, - """ - # `NonconvexQuadraticConstraint` - - ## What - - A `NonconvexQuadraticConstraint` issue is identified when a quadratic - constraint is nonconvex, that is, the quadratic matrix is not positive - semidefinite. - - ## Why - - Nonconvex constraints are not expected by many solver and can lead to - wrong solutions or even convergence issues. - - ## How to fix - - Check if the constraint is correct. Coefficient signs might have been - inverted. This also occurs if user fix a variable to emulate a - parameter, in this case some solvers will not be able to solve the - model properly, other tools such as ParametricOptInteface.jl might be - more suitable than fixing variables. - - ## More information - - No extra information for this issue. - """, - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - ::Type{LargeMatrixQuadraticCoefficient}, -) - return print( - io, - """ - # `LargeMatrixQuadraticCoefficient` - - ## What - - A `LargeMatrixQuadraticCoefficient` issue is identified when a quadratic - constraint has a coefficient with a large absolute value. - - ## Why - - Large coefficients can lead to numerical instability in the solution - process. - - ## How to fix - - Check if the coefficient is correct. Check if the units of variables and - coefficients are correct. Check if the number makes is - reasonable given that solver have tolerances. Sometimes these - coefficients can be replaced by zeros. - - ## More information - - - https://jump.dev/JuMP.jl/stable/tutorials/getting_started/tolerances/ - """, - ) -end - -function ModelAnalyzer._summarize( - io::IO, - issue::VariableNotInConstraints, - model, -) - return print(io, ModelAnalyzer._name(issue.ref, model)) -end - -function ModelAnalyzer._summarize(io::IO, issue::EmptyConstraint, model) - return print(io, ModelAnalyzer._name(issue.ref, model)) -end - -function ModelAnalyzer._summarize( - io::IO, - issue::VariableBoundAsConstraint, - model, -) - return print(io, ModelAnalyzer._name(issue.ref, model)) -end - -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, model) - return print( - io, - ModelAnalyzer._name(issue.ref, model), - " -- ", - ModelAnalyzer._name(issue.variable, model), - " : ", - issue.coefficient, - ) -end - -function ModelAnalyzer._summarize(io::IO, issue::LargeMatrixCoefficient, model) - return print( - io, - ModelAnalyzer._name(issue.ref, model), - " -- ", - ModelAnalyzer._name(issue.variable, model), - " : ", - issue.coefficient, - ) -end - -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, model) - return print( - io, - ModelAnalyzer._name(issue.variable, model), - " : ", - issue.coefficient, - ) -end - -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, model) - return print( - io, - ModelAnalyzer._name(issue.ref, model), - " : ", - issue.coefficient, - ) -end - -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, - model, -) - return print( - io, - ModelAnalyzer._name(issue.variable, model), - " : ", - issue.coefficient, - ) -end - -function ModelAnalyzer._summarize( - io::IO, - issue::SmallObjectiveQuadraticCoefficient, - model, -) - return print( - io, - ModelAnalyzer._name(issue.variable1, model), - " -- ", - ModelAnalyzer._name(issue.variable2, model), - " : ", - issue.coefficient, - ) -end - -function ModelAnalyzer._summarize( - io::IO, - issue::LargeObjectiveQuadraticCoefficient, - model, -) - return print( - io, - ModelAnalyzer._name(issue.variable1, model), - " -- ", - ModelAnalyzer._name(issue.variable2, model), - " : ", - issue.coefficient, - ) -end - -function ModelAnalyzer._summarize( - io::IO, - issue::SmallMatrixQuadraticCoefficient, - model, -) - return print( - io, - ModelAnalyzer._name(issue.ref, model), - " -- ", - ModelAnalyzer._name(issue.variable1, model), - " -- ", - ModelAnalyzer._name(issue.variable2, model), - " : ", - issue.coefficient, - ) -end - -function ModelAnalyzer._summarize( - io::IO, - issue::LargeMatrixQuadraticCoefficient, - model, -) - return print( - io, - ModelAnalyzer._name(issue.ref, model), - " -- ", - ModelAnalyzer._name(issue.variable1, model), - " -- ", - ModelAnalyzer._name(issue.variable2, model), - " : ", - issue.coefficient, - ) -end - -function ModelAnalyzer._summarize(io::IO, ::NonconvexQuadraticObjective, model) - return print(io, "Objective is Nonconvex quadratic") -end - -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: ", ModelAnalyzer._name(issue.ref, model)) -end - -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: ", ModelAnalyzer._name(issue.ref, model)) -end - -function ModelAnalyzer._verbose_summarize(io::IO, issue::DenseConstraint, model) - return print( - io, - "Constraint: ", - ModelAnalyzer._name(issue.ref, model), - " with ", - issue.nnz, - " non zero coefficients", - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - issue::SmallMatrixCoefficient, - model, -) - return print( - io, - "(Constraint -- Variable): (", - ModelAnalyzer._name(issue.ref, model), - " -- ", - ModelAnalyzer._name(issue.variable, model), - ") with coefficient ", - issue.coefficient, - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - issue::LargeMatrixCoefficient, - model, -) - return print( - io, - "(Constraint -- Variable): (", - ModelAnalyzer._name(issue.ref, model), - " -- ", - ModelAnalyzer._name(issue.variable, model), - ") with coefficient ", - issue.coefficient, - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - issue::SmallBoundCoefficient, - model, -) - return print( - io, - "Variable: ", - ModelAnalyzer._name(issue.variable, model), - " with bound ", - issue.coefficient, - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - issue::LargeBoundCoefficient, - model, -) - return print( - io, - "Variable: ", - ModelAnalyzer._name(issue.variable, model), - " with bound ", - issue.coefficient, - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - issue::SmallRHSCoefficient, - model, -) - return print( - io, - "Constraint: ", - ModelAnalyzer._name(issue.ref, model), - " with right-hand-side ", - issue.coefficient, - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - issue::LargeRHSCoefficient, - model, -) - return print( - io, - "Constraint: ", - ModelAnalyzer._name(issue.ref, model), - " with right-hand-side ", - issue.coefficient, - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - issue::SmallObjectiveCoefficient, - model, -) - return print( - io, - "Variable: ", - ModelAnalyzer._name(issue.variable, model), - " with coefficient ", - issue.coefficient, - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - issue::LargeObjectiveCoefficient, - model, -) - return print( - io, - "Variable: ", - ModelAnalyzer._name(issue.variable, model), - " with coefficient ", - issue.coefficient, - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - issue::SmallObjectiveQuadraticCoefficient, - model, -) - return print( - io, - "(Variable -- Variable): (", - ModelAnalyzer._name(issue.variable1, model), - " -- ", - ModelAnalyzer._name(issue.variable2, model), - ") with coefficient ", - issue.coefficient, - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - issue::LargeObjectiveQuadraticCoefficient, - model, -) - return print( - io, - "(Variable -- Variable): (", - ModelAnalyzer._name(issue.variable1, model), - " -- ", - ModelAnalyzer._name(issue.variable2, model), - ") with coefficient ", - issue.coefficient, - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - issue::SmallMatrixQuadraticCoefficient, - model, -) - return print( - io, - "(Constraint -- Variable -- Variable): (", - ModelAnalyzer._name(issue.ref, model), - " -- ", - ModelAnalyzer._name(issue.variable1, model), - " -- ", - ModelAnalyzer._name(issue.variable2, model), - ") with coefficient ", - issue.coefficient, - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - issue::LargeMatrixQuadraticCoefficient, - model, -) - return print( - io, - "(Constraint -- Variable -- Variable): (", - ModelAnalyzer._name(issue.ref, model), - " -- ", - ModelAnalyzer._name(issue.variable1, model), - " -- ", - ModelAnalyzer._name(issue.variable2, model), - ") with coefficient ", - issue.coefficient, - ) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - issue::NonconvexQuadraticObjective, - model, -) - return ModelAnalyzer._summarize(io, issue, model) -end - -function ModelAnalyzer._verbose_summarize( - io::IO, - issue::NonconvexQuadraticConstraint, - model, -) - return print(io, "Constraint: ", ModelAnalyzer._name(issue.ref, model)) -end - -function ModelAnalyzer.list_of_issues( - data::Data, - ::Type{VariableNotInConstraints}, -) - return data.variables_not_in_constraints -end - -function ModelAnalyzer.list_of_issues(data::Data, ::Type{EmptyConstraint}) - return data.empty_rows -end - -function ModelAnalyzer.list_of_issues( - data::Data, - ::Type{VariableBoundAsConstraint}, -) - return data.bound_rows -end - -function ModelAnalyzer.list_of_issues(data::Data, ::Type{DenseConstraint}) - return data.dense_rows -end - -function ModelAnalyzer.list_of_issues( - data::Data, - ::Type{SmallMatrixCoefficient}, -) - return data.matrix_small -end - -function ModelAnalyzer.list_of_issues( - data::Data, - ::Type{LargeMatrixCoefficient}, -) - return data.matrix_large -end - -function ModelAnalyzer.list_of_issues(data::Data, ::Type{SmallBoundCoefficient}) - return data.bounds_small -end - -function ModelAnalyzer.list_of_issues(data::Data, ::Type{LargeBoundCoefficient}) - return data.bounds_large -end - -function ModelAnalyzer.list_of_issues(data::Data, ::Type{SmallRHSCoefficient}) - return data.rhs_small -end - -function ModelAnalyzer.list_of_issues(data::Data, ::Type{LargeRHSCoefficient}) - return data.rhs_large -end - -function ModelAnalyzer.list_of_issues( - data::Data, - ::Type{SmallObjectiveCoefficient}, -) - return data.objective_small -end - -function ModelAnalyzer.list_of_issues( - data::Data, - ::Type{LargeObjectiveCoefficient}, -) - return data.objective_large -end - -function ModelAnalyzer.list_of_issues( - data::Data, - ::Type{SmallObjectiveQuadraticCoefficient}, -) - return data.objective_quadratic_small -end - -function ModelAnalyzer.list_of_issues( - data::Data, - ::Type{LargeObjectiveQuadraticCoefficient}, -) - return data.objective_quadratic_large -end - -function ModelAnalyzer.list_of_issues( - data::Data, - ::Type{SmallMatrixQuadraticCoefficient}, -) - return data.matrix_quadratic_small -end - -function ModelAnalyzer.list_of_issues( - data::Data, - ::Type{LargeMatrixQuadraticCoefficient}, -) - return data.matrix_quadratic_large -end - -function ModelAnalyzer.list_of_issues( - data::Data, - ::Type{NonconvexQuadraticObjective}, -) - return data.nonconvex_objective -end - -function ModelAnalyzer.list_of_issues( - data::Data, - ::Type{NonconvexQuadraticConstraint}, -) - return data.nonconvex_rows -end - -function ModelAnalyzer.list_of_issue_types(data::Data) - ret = Type[] - for type in ( - VariableNotInConstraints, - EmptyConstraint, - VariableBoundAsConstraint, - DenseConstraint, - SmallMatrixCoefficient, - LargeMatrixCoefficient, - SmallBoundCoefficient, - LargeBoundCoefficient, - SmallRHSCoefficient, - LargeRHSCoefficient, - SmallObjectiveCoefficient, - LargeObjectiveCoefficient, - SmallObjectiveQuadraticCoefficient, - LargeObjectiveQuadraticCoefficient, - SmallMatrixQuadraticCoefficient, - LargeMatrixQuadraticCoefficient, - NonconvexQuadraticConstraint, - NonconvexQuadraticObjective, - ) - if !isempty(ModelAnalyzer.list_of_issues(data, type)) - push!(ret, type) - end - end - return ret -end - -function summarize_configurations(io::IO, data::Data) - print(io, "## Configuration\n\n") - print(io, " Dense fill-in threshold: ", data.threshold_dense_fill_in, "\n") - print(io, " Dense entries threshold: ", data.threshold_dense_entries, "\n") - print(io, " Small coefficient threshold: ", data.threshold_small, "\n") - print(io, " Large coefficient threshold: ", data.threshold_large, "\n") - return -end - -function summarize_dimensions(io::IO, data::Data) - print(io, "## Dimensions\n\n") - print(io, " Number of variables: ", data.number_of_variables, "\n") - print(io, " Number of constraints: ", data.number_of_constraints, "\n") - print(io, " Number of nonzeros in matrix: ", data.matrix_nnz, "\n") - # types - println(io, " Constraint types:") - for (F, S, n) in data.constraint_info - println(io, " * ", F, "-", S, ": ", n) - end - return -end - -function summarize_ranges(io::IO, data::Data) - print(io, "## Coefficient ranges\n\n") - print(io, " Matrix: ", _stringify_bounds(data.matrix_range), "\n") - print(io, " Objective: ", _stringify_bounds(data.objective_range), "\n") - print(io, " Bounds: ", _stringify_bounds(data.bounds_range), "\n") - print(io, " RHS: ", _stringify_bounds(data.rhs_range), "\n") - if data.has_quadratic_objective - print( - io, - " Objective quadratic: ", - _stringify_bounds(data.objective_quadratic_range), - "\n", - ) - end - if data.has_quadratic_constraints - print( - io, - " Matrix quadratic: ", - _stringify_bounds(data.matrix_quadratic_range), - "\n", - ) - end - return -end - -function ModelAnalyzer.summarize( - io::IO, - data::Data; - model = nothing, - verbose = true, - max_issues = ModelAnalyzer.DEFAULT_MAX_ISSUES, - configurations = true, - dimensions = true, - ranges = true, -) - print(io, "## Numerical Analysis\n\n") - if configurations - summarize_configurations(io, data) - print(io, "\n") - end - if dimensions - summarize_dimensions(io, data) - print(io, "\n") - end - if ranges - summarize_ranges(io, data) - print(io, "\n") - end - for issue_type in ModelAnalyzer.list_of_issue_types(data) - issues = ModelAnalyzer.list_of_issues(data, issue_type) - print(io, "\n\n") - ModelAnalyzer.summarize( - io, - issues, - model = model, - verbose = verbose, - max_issues = max_issues, - ) - end - return -end - -function Base.show(io::IO, data::Data) - n = sum( - length(ModelAnalyzer.list_of_issues(data, T)) for - T in ModelAnalyzer.list_of_issue_types(data); - init = 0, - ) - return print(io, "Numerical analysis found $n issues") -end - -# printing helpers - -_print_value(x::Real) = Printf.@sprintf("%1.0e", x) - -function _stringify_bounds(bounds::Vector{Float64}) - lower = bounds[1] < Inf ? _print_value(bounds[1]) : "0e+00" - upper = bounds[2] > -Inf ? _print_value(bounds[2]) : "0e+00" - return string("[", lower, ", ", upper, "]") -end - -end diff --git a/test/runtests.jl b/test/runtests.jl index 12287cd..3c8b08f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,11 +7,10 @@ using Test @testset "ModelAnalyzer" begin for file in readdir(@__DIR__) - if !endswith(file, ".jl") || file in ("runtests.jl",) - continue - end - @testset "$file" begin - include(file) + if startswith(file, "test_") && endswith(file, ".jl") + @testset "$file" begin + include(file) + end end end end diff --git a/test/feasibility.jl b/test/test_Feasibility.jl similarity index 99% rename from test/feasibility.jl rename to test/test_Feasibility.jl index 793e56d..3948fab 100644 --- a/test/feasibility.jl +++ b/test/test_Feasibility.jl @@ -3,7 +3,7 @@ # 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 TestDualFeasibilityChecker +module TestFeasibility using ModelAnalyzer using Test @@ -30,7 +30,7 @@ function test_no_solution() model, ) # no dual solutions available - # @test_throws ErrorException + # @test_throws ErrorException data = ModelAnalyzer.analyze( ModelAnalyzer.Feasibility.Analyzer(), model, @@ -195,7 +195,7 @@ function test_no_lb() # the dual is: # Max 0 # Subject to - # y == 1 (as a constraint) + # y == 1 (as a constraint) # y >= 0 (as a bound) # mayber force fail here # @test_throws ErrorException @@ -1066,6 +1066,6 @@ function test_nl_obj() return end -end # module +end # module TestFeasibility -TestDualFeasibilityChecker.runtests() +TestFeasibility.runtests() diff --git a/test/infeasibility.jl b/test/test_Infeasibility.jl similarity index 99% rename from test/infeasibility.jl rename to test/test_Infeasibility.jl index d9d7dee..423ed7f 100644 --- a/test/infeasibility.jl +++ b/test/test_Infeasibility.jl @@ -3,12 +3,12 @@ # 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 TestInfeasibilityChecker +module TestInfeasibility -import ModelAnalyzer -using Test using JuMP +using Test import HiGHS +import ModelAnalyzer function runtests() for name in names(@__MODULE__; all = true) @@ -575,6 +575,6 @@ function test_iis_spare() return end -end # module +end # module TestInfeasibility -TestInfeasibilityChecker.runtests() +TestInfeasibility.runtests() diff --git a/test/numerical.jl b/test/test_Numerical.jl similarity index 99% rename from test/numerical.jl rename to test/test_Numerical.jl index 3d5d6df..5ad2e2f 100644 --- a/test/numerical.jl +++ b/test/test_Numerical.jl @@ -1320,6 +1320,6 @@ function test_more_than_max_issues() return end -end # module +end # module TestNumerical TestNumerical.runtests()