Skip to content

Commit 616e5c3

Browse files
joaquimgjoaquimg
committed
Numerical Analyzer
1 parent 2b868d6 commit 616e5c3

File tree

3 files changed

+292
-0
lines changed

3 files changed

+292
-0
lines changed

‎src/ModelAnalyzer.jl‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4,4 +4,12 @@
44
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.
55

66
module ModelAnalyzer
7+
8+
using JuMP
9+
import LinearAlgebra
10+
import Printf
11+
import MathOptInterface as MOI
12+
13+
include("coefficients.jl")
14+
715
end

‎src/coefficients.jl‎

Lines changed: 271 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,271 @@
1+
# Copyright (c) 2025: Joaquim Garcia, Oscar Dowson and contributors
2+
#
3+
# Use of this source code is governed by an MIT-style license that can be found
4+
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.
5+
6+
Base.@kwdef mutable struct CoefficientsData
7+
8+
number_of_variables::Int = 0
9+
number_of_constraints::Int = 0
10+
11+
constraint_info::Vector{Tuple{DataType, DataType, Int}} = Tuple{DataType, DataType, Int}[]
12+
13+
matrix_nnz::Int = 0
14+
15+
matrix_range::Vector{Float64} = sizehint!(Float64[1.0, 1.0], 2)
16+
bounds_range::Vector{Float64} = sizehint!(Float64[1.0, 1.0], 2)
17+
rhs_range::Vector{Float64} = sizehint!(Float64[1.0, 1.0], 2)
18+
objective_range::Vector{Float64} = sizehint!(Float64[1.0, 1.0], 2)
19+
20+
threshold_dense_row::Float64 = 0.10
21+
threshold_small_coefficient::Float64 = 1e-5
22+
threshold_large_coefficient::Float64 = 1e+5
23+
24+
bound_rows::Vector{ConstraintRef} = ConstraintRef[]
25+
dense_rows::Vector{Tuple{ConstraintRef, Int}} = Tuple{ConstraintRef, Int}[]
26+
small_coefficients::Vector{Tuple{ConstraintRef, VariableRef, Float64}} = Tuple{ConstraintRef, VariableRef, Float64}[]
27+
large_coefficients::Vector{Tuple{ConstraintRef, VariableRef, Float64}} = Tuple{ConstraintRef, VariableRef, Float64}[]
28+
29+
end
30+
31+
function _update_range(range::Vector{Float64}, value::Number)
32+
if !(value ≈ 0.0)
33+
range[1] = min(range[1], abs(value))
34+
range[2] = max(range[2], abs(value))
35+
return true
36+
end
37+
return false
38+
end
39+
40+
function _get_data(data, ref::ConstraintRef, func::JuMP.GenericAffExpr)
41+
if length(func.terms) == 1
42+
if first(values(func.terms)) ≈ 1.0
43+
push!(data.bound_rows, ref)
44+
data.matrix_nnz += 1
45+
return
46+
end
47+
end
48+
nnz = 0
49+
for (variable, coefficient) in func.terms
50+
nnz += _update_range(data.matrix_range, coefficient)
51+
if abs(coefficient) < data.threshold_small_coefficient
52+
push!(data.small_coefficients, (ref, variable, coefficient))
53+
elseif abs(coefficient) > data.threshold_large_coefficient
54+
push!(data.large_coefficients, (ref, variable, coefficient))
55+
end
56+
end
57+
if nnz / data.number_of_variables > data.threshold_dense_row && nnz > 100
58+
push!(data.dense_rows, (ref, nnz))
59+
end
60+
data.matrix_nnz += nnz
61+
return
62+
end
63+
64+
function _update_range(range::Vector, func::JuMP.GenericAffExpr)
65+
_update_range(range, func.constant)
66+
return true
67+
end
68+
69+
function _get_data(data, func::Vector{JuMP.GenericAffExpr}, set)
70+
for f in func
71+
_update_range(data, f, set)
72+
end
73+
return true
74+
end
75+
76+
function _get_data(data, func::JuMP.GenericAffExpr, set)
77+
_update_range(data.rhs_range, func.constant)
78+
return true
79+
end
80+
81+
function _get_data(data, func::JuMP.GenericAffExpr, set::MOI.LessThan)
82+
_update_range(data.rhs_range, set.upper - func.constant)
83+
return true
84+
end
85+
86+
function _get_data(data, func::JuMP.GenericAffExpr, set::MOI.GreaterThan)
87+
_update_range(data.rhs_range, set.lower - func.constant)
88+
return true
89+
end
90+
91+
function _get_data(data, func::JuMP.GenericAffExpr, set::MOI.EqualTo)
92+
_update_range(data.rhs_range, set.value - func.constant)
93+
return true
94+
end
95+
96+
function _get_data(data, func::JuMP.GenericAffExpr, set::MOI.Interval)
97+
_update_range(data.rhs_range, set.upper - func.constant)
98+
_update_range(data.rhs_range, set.lower - func.constant)
99+
return true
100+
end
101+
102+
# Default fallback for unsupported constraints.
103+
_update_range(data, func, set) = false
104+
105+
function coefficient_analysis(model::JuMP.Model)
106+
data = CoefficientsData()
107+
data.number_of_variables = JuMP.num_variables(model)
108+
data.number_of_constraints = JuMP.num_constraints(model, count_variable_in_set_constraints = false)
109+
_update_range(data.objective_range, JuMP.objective_function(model))
110+
for var in JuMP.all_variables(model)
111+
if JuMP.has_lower_bound(var)
112+
_update_range(data.bounds_range, JuMP.lower_bound(var))
113+
end
114+
if JuMP.has_upper_bound(var)
115+
_update_range(data.bounds_range, JuMP.upper_bound(var))
116+
end
117+
end
118+
for (F, S) in JuMP.list_of_constraint_types(model)
119+
n = JuMP.num_constraints(model, F, S)
120+
if n > 0
121+
push!(data.constraint_info, (F, S, n))
122+
end
123+
F == JuMP.VariableRef && continue
124+
F == Vector{JuMP.VariableRef} && continue
125+
for con in JuMP.all_constraints(model, F, S)
126+
con_obj = JuMP.constraint_object(con)
127+
_get_data(data, con, con_obj.func)
128+
_get_data(data, con_obj.func, con_obj.set)
129+
end
130+
end
131+
sort!(data.dense_rows, by = x -> x[2], rev = true)
132+
sort!(data.small_coefficients, by = x -> abs(x[3]))
133+
sort!(data.large_coefficients, by = x -> abs(x[3]), rev = true)
134+
return data
135+
end
136+
137+
# printing
138+
139+
_print_value(x::Real) = Printf.@sprintf("%1.0e", x)
140+
141+
function _stringify_bounds(bounds::Vector{Float64})
142+
lower = bounds[1] < Inf ? _print_value(bounds[1]) : "0e+00"
143+
upper = bounds[2] > -Inf ? _print_value(bounds[2]) : "0e+00"
144+
return string("[", lower, ", ", upper, "]")
145+
end
146+
147+
function _print_coefficients(
148+
io::IO,
149+
name::String,
150+
data,
151+
range,
152+
warnings::Vector{Tuple{String,String}},
153+
)
154+
println(
155+
io,
156+
" ",
157+
rpad(string(name, " range"), 17),
158+
_stringify_bounds(range),
159+
)
160+
if range[1] < data.threshold_small_coefficient
161+
push!(warnings, (name, "small"))
162+
end
163+
if range[2] > data.threshold_large_coefficient
164+
push!(warnings, (name, "large"))
165+
end
166+
return
167+
end
168+
169+
function _print_numerical_stability_report(
170+
io::IO,
171+
data::CoefficientsData;
172+
warn::Bool = true,
173+
verbose::Bool = true,
174+
max_list::Int = 10,
175+
)
176+
println(io, "Numerical stability report:")
177+
println(io, " Number of variables: ", data.number_of_variables)
178+
println(io, " Number of constraints: ", data.number_of_constraints)
179+
println(io, " Number of nonzeros in matrix: ", data.matrix_nnz)
180+
println(io, " Threshold for dense rows: ", data.threshold_dense_row)
181+
println(io, " Threshold for small coefficients: ", data.threshold_small_coefficient)
182+
println(io, " Threshold for large coefficients: ", data.threshold_large_coefficient)
183+
184+
println(io, " Coefficient ranges:")
185+
warnings = Tuple{String, String}[]
186+
_print_coefficients(io, "matrix", data, data.matrix_range, warnings)
187+
_print_coefficients(io, "objective", data, data.objective_range, warnings)
188+
_print_coefficients(io, "bounds", data, data.bounds_range, warnings)
189+
_print_coefficients(io, "rhs", data, data.rhs_range, warnings)
190+
191+
# types
192+
println(io, " Constraint types:")
193+
for (F, S, n) in data.constraint_info
194+
println(io, " * ", F, "-", S, ": ", n)
195+
end
196+
197+
# rows that should be bounds
198+
println(io, " Bound rows: ", length(data.bound_rows))
199+
if verbose
200+
c = 0
201+
for ref in data.bound_rows
202+
println(io, " * ", ref)
203+
c += 1
204+
if c >= max_list
205+
break
206+
end
207+
end
208+
end
209+
210+
println(io, " Dense constraints: ", length(data.dense_rows))
211+
println(io, " Small coefficients: ", length(data.small_coefficients))
212+
println(io, " Large coefficients: ", length(data.large_coefficients))
213+
214+
if verbose
215+
println(io, "")
216+
println(io, " Dense constraints:")
217+
c = 0
218+
for (ref, nnz) in data.dense_rows
219+
println(io, " * ", ref, ": ", nnz)
220+
c += 1
221+
if c >= max_list
222+
break
223+
end
224+
end
225+
println(io, "")
226+
println(io, " Small coefficients:")
227+
c = 0
228+
for (ref, var, coeff) in data.small_coefficients
229+
println(io, " * ", ref, ": ", var, " -> ", coeff)
230+
c += 1
231+
if c >= max_list
232+
break
233+
end
234+
end
235+
println(io, "")
236+
println(io, " Large coefficients:")
237+
c = 0
238+
for (ref, var, coeff) in data.large_coefficients
239+
println(io, " * ", ref, ": ", var, " -> ", coeff)
240+
c += 1
241+
if c >= max_list
242+
break
243+
end
244+
end
245+
end
246+
247+
if warn && !isempty(warnings)
248+
println(io, "\nWARNING: numerical stability issues detected")
249+
for (name, sense) in warnings
250+
println(io, " - $(name) range contains $(sense) coefficients")
251+
end
252+
println(
253+
io,
254+
"Very large or small absolute values of coefficients\n",
255+
"can cause numerical stability issues. Consider\n",
256+
"reformulating the model.",
257+
)
258+
end
259+
return
260+
end
261+
262+
function Base.show(io::IO, data::CoefficientsData; verbose::Bool = false)
263+
_print_numerical_stability_report(io, data, warn = true, verbose = verbose)
264+
return
265+
end
266+
267+
# TODO add names in the output
268+
# TODO add lists for rhs, bounds and obj coefs
269+
# TODO analyse quadratics
270+
# check variable that are not in constraints
271+
# check start poitn in bounds

‎test/runtests.jl‎

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,4 +7,17 @@ import ModelAnalyzer
77
import MathOptInterface as MOI
88
using Test
99
using JuMP
10+
import HiGHS
1011

12+
model = Model()
13+
@variable(model, x <= 1e9)
14+
@variable(model, y >= 1e-9)
15+
@constraint(model, x + y <= 1e8)
16+
@constraint(model, x + y + 1e7 <= 2)
17+
@constraint(model, 1e6*x + 1e-5*y >= 2)
18+
@constraint(model, x <= 100)
19+
@objective(model, Max, x + y)
20+
21+
data = ModelAnalyzer.coefficient_analysis(model)
22+
23+
show(data)

0 commit comments

Comments
 (0)