Skip to content

Commit d018188

Browse files
joaquimgjoaquimg
committed
add quadratic
1 parent 50112cc commit d018188

File tree

2 files changed

+228
-6
lines changed

2 files changed

+228
-6
lines changed

‎src/ModelAnalyzer.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,7 @@ using JuMP
99
import LinearAlgebra
1010
import Printf
1111
import MathOptInterface as MOI
12+
import OrderedCollections
1213

1314
include("coefficients.jl")
1415

‎src/coefficients.jl‎

Lines changed: 227 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -29,8 +29,6 @@ Base.@kwdef mutable struct CoefficientsData
2929
bound_rows::Vector{ConstraintRef} = ConstraintRef[]
3030
dense_rows::Vector{Tuple{ConstraintRef,Int}} = Tuple{ConstraintRef,Int}[]
3131

32-
nonconvex_rows::Vector{ConstraintRef} = ConstraintRef[]
33-
3432
matrix_small::Vector{Tuple{ConstraintRef,VariableRef,Float64}} =
3533
Tuple{ConstraintRef,VariableRef,Float64}[]
3634
matrix_large::Vector{Tuple{ConstraintRef,VariableRef,Float64}} =
@@ -53,10 +51,14 @@ Base.@kwdef mutable struct CoefficientsData
5351

5452
has_quadratic_objective::Bool = false
5553
has_quadratic_constraints::Bool = false
54+
sense::JuMP.OptimizationSense = JuMP.FEASIBILITY_SENSE
5655

5756
objective_quadratic_range = sizehint!(Float64[1.0, 1.0], 2)
5857
matrix_quadratic_range = sizehint!(Float64[1.0, 1.0], 2)
5958

59+
nonconvex_objective::Bool = false
60+
nonconvex_rows::Vector{ConstraintRef} = ConstraintRef[]
61+
6062
matrix_quadratic_small::Vector{
6163
Tuple{ConstraintRef,VariableRef,VariableRef,Float64},
6264
} = Tuple{ConstraintRef,VariableRef,VariableRef,Float64}[]
@@ -67,6 +69,7 @@ Base.@kwdef mutable struct CoefficientsData
6769
Tuple{VariableRef,VariableRef,Float64}[]
6870
objective_quadratic_large::Vector{Tuple{VariableRef,VariableRef,Float64}} =
6971
Tuple{VariableRef,VariableRef,Float64}[]
72+
7073
end
7174

7275
function _update_range(range::Vector{Float64}, value::Number)
@@ -112,10 +115,29 @@ function _get_constraint_data(
112115
return
113116
end
114117

115-
# function _update_range(range::Vector, func::JuMP.GenericAffExpr)
116-
# _update_range(range, func.constant)
117-
# return true
118-
# end
118+
function _get_constraint_data(
119+
data,
120+
ref::ConstraintRef,
121+
func::JuMP.GenericQuadExpr,
122+
)
123+
_get_constraint_data(data, ref, func.aff)
124+
nnz = 0
125+
for ((v1, v2), coefficient) in func.terms
126+
if coefficient ≈ 0.0
127+
continue
128+
end
129+
nnz += _update_range(data.matrix_quadratic_range, coefficient)
130+
if abs(coefficient) < data.threshold_small
131+
push!(data.matrix_quadratic_small, (ref, v1, v2, coefficient))
132+
elseif abs(coefficient) > data.threshold_large
133+
push!(data.matrix_quadratic_large, (ref, v1, v2, coefficient))
134+
end
135+
push!(data.variables_in_constraints, v1)
136+
push!(data.variables_in_constraints, v2)
137+
end
138+
data.has_quadratic_constraints = true
139+
return
140+
end
119141

120142
function _get_variable_data(data, variable, coefficient::Number)
121143
if !(coefficient ≈ 0.0)
@@ -145,13 +167,67 @@ function _get_objective_data(data, func::JuMP.GenericAffExpr)
145167
return
146168
end
147169

170+
function _get_objective_data(data, func::JuMP.GenericQuadExpr)
171+
_get_objective_data(data, func.aff)
172+
nnz = 0
173+
for ((v1, v2), coefficient) in func.terms
174+
if coefficient ≈ 0.0
175+
continue
176+
end
177+
nnz += _update_range(data.objective_quadratic_range, coefficient)
178+
if abs(coefficient) < data.threshold_small
179+
push!(data.objective_quadratic_small, (v1, v2, coefficient))
180+
elseif abs(coefficient) > data.threshold_large
181+
push!(data.objective_quadratic_large, (v1, v2, coefficient))
182+
end
183+
end
184+
data.has_quadratic_objective = true
185+
if data.sense == JuMP.MAX_SENSE
186+
data.nonconvex_objective = !_quadratic_vexity(func, -1)
187+
elseif data.sense == JuMP.MIN_SENSE
188+
data.nonconvex_objective = !_quadratic_vexity(func, 1)
189+
end
190+
return
191+
end
192+
193+
function _quadratic_vexity(func::JuMP.GenericQuadExpr, sign::Int)
194+
variables = OrderedCollections.OrderedSet{VariableRef}()
195+
sizehint!(variables, 2 * length(func.terms))
196+
for (v1, v2) in keys(func.terms)
197+
push!(variables, v1)
198+
push!(variables, v2)
199+
end
200+
var_map = Dict{VariableRef,Int}()
201+
for (idx, var) in enumerate(variables)
202+
var_map[var] = idx
203+
end
204+
matrix = zeros(length(variables), length(variables))
205+
for ((v1, v2), coefficient) in func.terms
206+
matrix[var_map[v1], var_map[v2]] += sign * coefficient / 2
207+
matrix[var_map[v2], var_map[v1]] += sign * coefficient / 2
208+
end
209+
ret = LinearAlgebra.cholesky!(
210+
LinearAlgebra.Symmetric(matrix),
211+
LinearAlgebra.RowMaximum(),
212+
check = false,
213+
)
214+
return LinearAlgebra.issuccess(ret)
215+
end
216+
148217
function _get_constraint_data(data, func::Vector{JuMP.GenericAffExpr}, set)
149218
for f in func
150219
_get_constraint_data(data, ref, f, set)
151220
end
152221
return true
153222
end
154223

224+
function _get_constraint_data(data, func::Vector{JuMP.GenericQuadExpr}, set)
225+
for f in func
226+
_get_constraint_data(data, ref, f, set)
227+
end
228+
return true
229+
end
230+
155231
function _get_constraint_data(data, ref, func::JuMP.GenericAffExpr, set)
156232
coefficient = func.constant
157233
if coefficient ≈ 0.0
@@ -166,6 +242,25 @@ function _get_constraint_data(data, ref, func::JuMP.GenericAffExpr, set)
166242
return
167243
end
168244

245+
function _get_constraint_data(data, ref, func::JuMP.GenericQuadExpr, set)
246+
_get_constraint_data(data, ref, func.aff, set)
247+
# skip additional checks for quadratics in non-scalar simples sets
248+
return
249+
end
250+
251+
function _get_constraint_data(
252+
data,
253+
ref,
254+
func::JuMP.GenericQuadExpr,
255+
set::MOI.LessThan,
256+
)
257+
_get_constraint_data(data, ref, func.aff, set)
258+
if !_quadratic_vexity(func, 1)
259+
push!(data.nonconvex_rows, ref)
260+
end
261+
return
262+
end
263+
169264
function _get_constraint_data(
170265
data,
171266
ref,
@@ -185,6 +280,19 @@ function _get_constraint_data(
185280
return
186281
end
187282

283+
function _get_constraint_data(
284+
data,
285+
ref,
286+
func::JuMP.GenericQuadExpr,
287+
set::MOI.GreaterThan,
288+
)
289+
_get_constraint_data(data, ref, func.aff, set)
290+
if !_quadratic_vexity(func, -1)
291+
push!(data.nonconvex_rows, ref)
292+
end
293+
return
294+
end
295+
188296
function _get_constraint_data(
189297
data,
190298
ref,
@@ -204,6 +312,17 @@ function _get_constraint_data(
204312
return
205313
end
206314

315+
function _get_constraint_data(
316+
data,
317+
ref,
318+
func::JuMP.GenericQuadExpr,
319+
set::Union{MOI.EqualTo,MOI.Interval},
320+
)
321+
_get_constraint_data(data, ref, func.aff, set)
322+
push!(data.nonconvex_rows, ref)
323+
return
324+
end
325+
207326
function _get_constraint_data(
208327
data,
209328
ref,
@@ -263,6 +382,7 @@ _update_range(data, func, set) = false
263382

264383
function coefficient_analysis(model::JuMP.Model)
265384
data = CoefficientsData()
385+
data.sense = JuMP.objective_sense(model)
266386
data.number_of_variables = JuMP.num_variables(model)
267387
sizehint!(data.variables_in_constraints, data.number_of_variables)
268388
data.number_of_constraints =
@@ -369,12 +489,40 @@ function _print_numerical_stability_report(
369489
println(io, " Small coefficients: ", data.threshold_small)
370490
println(io, " Large coefficients: ", data.threshold_large)
371491

492+
if data.has_quadratic_objective
493+
println(io, "\n Objective is quadratic:")
494+
if data.nonconvex_objective
495+
println(io, " Objective is nonconvex (numerically)")
496+
else
497+
println(io, " Objective is convex (numerically)")
498+
end
499+
println(io, "")
500+
end
501+
372502
println(io, " Coefficient ranges:")
373503
warnings = Tuple{String,String}[]
374504
_print_coefficients(io, "matrix", data, data.matrix_range, warnings)
375505
_print_coefficients(io, "objective", data, data.objective_range, warnings)
376506
_print_coefficients(io, "bounds", data, data.bounds_range, warnings)
377507
_print_coefficients(io, "rhs", data, data.rhs_range, warnings)
508+
if data.has_quadratic_objective
509+
_print_coefficients(
510+
io,
511+
"objective q",
512+
data,
513+
data.objective_quadratic_range,
514+
warnings,
515+
)
516+
end
517+
if data.has_quadratic_constraints
518+
_print_coefficients(
519+
io,
520+
"matrix q",
521+
data,
522+
data.matrix_quadratic_range,
523+
warnings,
524+
)
525+
end
378526

379527
# rows that should be bounds
380528
println(
@@ -385,6 +533,9 @@ function _print_numerical_stability_report(
385533
println(io, " Bound rows: ", length(data.bound_rows))
386534
println(io, " Dense constraints: ", length(data.dense_rows))
387535
println(io, " Empty constraints: ", length(data.empty_rows))
536+
if data.has_quadratic_objective
537+
println(io, " Nonconvex constraints: ", length(data.nonconvex_rows))
538+
end
388539
println(io, " Coefficients:")
389540
println(io, " matrix small: ", length(data.matrix_small))
390541
println(io, " matrix large: ", length(data.matrix_large))
@@ -448,6 +599,76 @@ function _print_numerical_stability_report(
448599
for (var, coeff) in first(data.objective_large, max_list)
449600
println(io, " * ", var, _name_string(var, names), ": ", coeff)
450601
end
602+
if data.has_quadratic_objective
603+
println(io, "\n Small objective quadratic coefficients:")
604+
for (v1, v2, coeff) in
605+
first(data.objective_quadratic_small, max_list)
606+
println(
607+
io,
608+
" * ",
609+
v1,
610+
_name_string(v1, names),
611+
"-",
612+
v2,
613+
_name_string(v2, names),
614+
": ",
615+
coeff,
616+
)
617+
end
618+
println(io, "\n Large objective quadratic coefficients:")
619+
for (v1, v2, coeff) in
620+
first(data.objective_quadratic_large, max_list)
621+
println(
622+
io,
623+
" * ",
624+
v1,
625+
_name_string(v1, names),
626+
"-",
627+
v2,
628+
_name_string(v2, names),
629+
": ",
630+
coeff,
631+
)
632+
end
633+
end
634+
if data.has_quadratic_constraints
635+
println(io, "\n Small matrix quadratic coefficients:")
636+
for (ref, v1, v2, coeff) in
637+
first(data.matrix_quadratic_small, max_list)
638+
println(
639+
io,
640+
" * ",
641+
ref,
642+
_name_string(ref, names),
643+
"-",
644+
v1,
645+
_name_string(v1, names),
646+
"-",
647+
v2,
648+
_name_string(v2, names),
649+
": ",
650+
coeff,
651+
)
652+
end
653+
println(io, "\n Large matrix quadratic coefficients:")
654+
for (ref, v1, v2, coeff) in
655+
first(data.matrix_quadratic_large, max_list)
656+
println(
657+
io,
658+
" * ",
659+
ref,
660+
_name_string(ref, names),
661+
"-",
662+
v1,
663+
_name_string(v1, names),
664+
"-",
665+
v2,
666+
_name_string(v2, names),
667+
": ",
668+
coeff,
669+
)
670+
end
671+
end
451672
end
452673

453674
if warn && !isempty(warnings)

0 commit comments

Comments
 (0)