Skip to content

Commit c010b53

Browse files
joaquimgjoaquimg
authored
Add parameter diff for quad terms (#342)
* Add parameter diff for quad terms * Update ParametricOptInterface to version 0.15.0 * cleanup deps * add tests * fix test
1 parent 1c98163 commit c010b53

File tree

3 files changed

+692
-1
lines changed

3 files changed

+692
-1
lines changed

‎Project.toml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -23,5 +23,5 @@ JuMP = "1"
2323
LazyArrays = "1, 2"
2424
MathOptInterface = "1.18"
2525
MathOptSetDistances = "0.2.9"
26-
ParametricOptInterface = "0.14.1, 0.15"
26+
ParametricOptInterface = "0.15.1"
2727
julia = "1.6"

‎src/parameters.jl‎

Lines changed: 211 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -266,6 +266,213 @@ function _quadratic_objective_set_forward!(model::POI.Optimizer{T}) where {T}
266266
return
267267
end
268268

269+
function _cubic_objective_set_forward!(model::POI.Optimizer{T}) where {T}
270+
pf = MOI.get(
271+
model,
272+
POI.ParametricObjectiveFunction{POI.ParametricCubicFunction{T}}(),
273+
)
274+
sensitivity_data = _get_sensitivity_data(model)
275+
276+
cte = zero(T)
277+
affine_terms = MOI.ScalarAffineTerm{T}[]
278+
quadratic_terms = MOI.ScalarQuadraticTerm{T}[]
279+
280+
# pvv terms: Δp * coeff → quadratic perturbation
281+
for term in POI.cubic_parameter_variable_variable_terms(pf)
282+
p = term.index_1
283+
v1 = term.index_2
284+
v2 = term.index_3
285+
Δp = get(sensitivity_data.parameter_input_forward, p, zero(T))
286+
if !iszero(Δp)
287+
qcoeff = Δp * term.coefficient
288+
if v1 == v2
289+
qcoeff *= 2 # MOI diagonal convention
290+
end
291+
push!(quadratic_terms, MOI.ScalarQuadraticTerm{T}(qcoeff, v1, v2))
292+
end
293+
end
294+
295+
# ppv terms: chain rule on two params → affine perturbation
296+
for term in POI.cubic_parameter_parameter_variable_terms(pf)
297+
p1 = term.index_1
298+
p2 = term.index_2
299+
v = term.index_3
300+
Δp1 = get(sensitivity_data.parameter_input_forward, p1, zero(T))
301+
Δp2 = get(sensitivity_data.parameter_input_forward, p2, zero(T))
302+
p1_val = MOI.get(model, MOI.VariablePrimal(), p1)
303+
p2_val = MOI.get(model, MOI.VariablePrimal(), p2)
304+
acoeff =
305+
Δp1 * term.coefficient * p2_val + Δp2 * term.coefficient * p1_val
306+
if !iszero(acoeff)
307+
push!(affine_terms, MOI.ScalarAffineTerm{T}(acoeff, v))
308+
end
309+
end
310+
311+
# ppp terms: chain rule on three params → constant perturbation
312+
for term in POI.cubic_parameter_parameter_parameter_terms(pf)
313+
p1 = term.index_1
314+
p2 = term.index_2
315+
p3 = term.index_3
316+
Δp1 = get(sensitivity_data.parameter_input_forward, p1, zero(T))
317+
Δp2 = get(sensitivity_data.parameter_input_forward, p2, zero(T))
318+
Δp3 = get(sensitivity_data.parameter_input_forward, p3, zero(T))
319+
p1_val = MOI.get(model, MOI.VariablePrimal(), p1)
320+
p2_val = MOI.get(model, MOI.VariablePrimal(), p2)
321+
p3_val = MOI.get(model, MOI.VariablePrimal(), p3)
322+
cte +=
323+
term.coefficient * (
324+
Δp1 * p2_val * p3_val +
325+
p1_val * Δp2 * p3_val +
326+
p1_val * p2_val * Δp3
327+
)
328+
end
329+
330+
# Degree-2: p terms (affine parameter → constant perturbation)
331+
for term in POI.cubic_affine_parameter_terms(pf)
332+
p = term.variable
333+
Δp = get(sensitivity_data.parameter_input_forward, p, zero(T))
334+
cte += Δp * term.coefficient
335+
end
336+
# Degree-2: pp terms (parameter-parameter → constant perturbation)
337+
for term in POI.cubic_parameter_parameter_terms(pf)
338+
p_1 = term.variable_1
339+
p_2 = term.variable_2
340+
Δp1 = get(sensitivity_data.parameter_input_forward, p_1, zero(T))
341+
Δp2 = get(sensitivity_data.parameter_input_forward, p_2, zero(T))
342+
p1_val = MOI.get(model, MOI.VariablePrimal(), p_1)
343+
p2_val = MOI.get(model, MOI.VariablePrimal(), p_2)
344+
cte += Δp1 * term.coefficient * p2_val / ifelse(p_1 === p_2, T(2), T(1))
345+
cte += Δp2 * term.coefficient * p1_val / ifelse(p_1 === p_2, T(2), T(1))
346+
end
347+
# Degree-2: pv terms (parameter-variable → affine perturbation)
348+
for term in POI.cubic_parameter_variable_terms(pf)
349+
p = term.variable_1
350+
Δp = get(sensitivity_data.parameter_input_forward, p, zero(T))
351+
if !iszero(Δp)
352+
push!(
353+
affine_terms,
354+
MOI.ScalarAffineTerm{T}(Δp * term.coefficient, term.variable_2),
355+
)
356+
end
357+
end
358+
359+
# Send perturbation to inner model
360+
if !isempty(quadratic_terms) || !isempty(affine_terms) || !iszero(cte)
361+
MOI.set(
362+
model.optimizer,
363+
ForwardObjectiveFunction(),
364+
MOI.ScalarQuadraticFunction{T}(quadratic_terms, affine_terms, cte),
365+
)
366+
end
367+
return
368+
end
369+
370+
function _cubic_objective_get_reverse!(model::POI.Optimizer{T}) where {T}
371+
pf = MOI.get(
372+
model,
373+
POI.ParametricObjectiveFunction{POI.ParametricCubicFunction{T}}(),
374+
)
375+
pvv_terms = POI.cubic_parameter_variable_variable_terms(pf)
376+
ppv_terms = POI.cubic_parameter_parameter_variable_terms(pf)
377+
ppp_terms = POI.cubic_parameter_parameter_parameter_terms(pf)
378+
p_terms = POI.cubic_affine_parameter_terms(pf)
379+
pp_terms = POI.cubic_parameter_parameter_terms(pf)
380+
pv_terms = POI.cubic_parameter_variable_terms(pf)
381+
if isempty(pvv_terms) &&
382+
isempty(ppv_terms) &&
383+
isempty(ppp_terms) &&
384+
isempty(p_terms) &&
385+
isempty(pp_terms) &&
386+
isempty(pv_terms)
387+
return
388+
end
389+
sensitivity_data = _get_sensitivity_data(model)
390+
grad_pf = MOI.get(model.optimizer, ReverseObjectiveFunction())
391+
d_cte = MOI.constant(grad_pf)
392+
393+
# pvv terms: ∇p += 2 * coeff * dQ[v1, v2]
394+
for term in pvv_terms
395+
p = term.index_1
396+
v1 = term.index_2
397+
v2 = term.index_3
398+
dQ_ij = quad_sym_half(grad_pf, v1, v2)
399+
value = get!(sensitivity_data.parameter_output_backward, p, zero(T))
400+
sensitivity_data.parameter_output_backward[p] =
401+
value + T(2) * term.coefficient * dQ_ij
402+
end
403+
404+
# ppv terms: chain rule on two params
405+
for term in ppv_terms
406+
p1 = term.index_1
407+
p2 = term.index_2
408+
v = term.index_3
409+
dq_v = JuMP.coefficient(grad_pf, v)
410+
p1_val = MOI.get(model, MOI.VariablePrimal(), p1)
411+
p2_val = MOI.get(model, MOI.VariablePrimal(), p2)
412+
val1 = get!(sensitivity_data.parameter_output_backward, p1, zero(T))
413+
val2 = get!(sensitivity_data.parameter_output_backward, p2, zero(T))
414+
sensitivity_data.parameter_output_backward[p1] =
415+
val1 + term.coefficient * p2_val * dq_v
416+
sensitivity_data.parameter_output_backward[p2] =
417+
val2 + term.coefficient * p1_val * dq_v
418+
end
419+
420+
# ppp terms: chain rule on three params
421+
for term in ppp_terms
422+
p1 = term.index_1
423+
p2 = term.index_2
424+
p3 = term.index_3
425+
p1_val = MOI.get(model, MOI.VariablePrimal(), p1)
426+
p2_val = MOI.get(model, MOI.VariablePrimal(), p2)
427+
p3_val = MOI.get(model, MOI.VariablePrimal(), p3)
428+
val1 = get!(sensitivity_data.parameter_output_backward, p1, zero(T))
429+
val2 = get!(sensitivity_data.parameter_output_backward, p2, zero(T))
430+
val3 = get!(sensitivity_data.parameter_output_backward, p3, zero(T))
431+
sensitivity_data.parameter_output_backward[p1] =
432+
val1 + term.coefficient * p2_val * p3_val * d_cte
433+
sensitivity_data.parameter_output_backward[p2] =
434+
val2 + term.coefficient * p1_val * p3_val * d_cte
435+
sensitivity_data.parameter_output_backward[p3] =
436+
val3 + term.coefficient * p1_val * p2_val * d_cte
437+
end
438+
439+
# Degree-2: p terms
440+
for term in p_terms
441+
p = term.variable
442+
value = get!(sensitivity_data.parameter_output_backward, p, zero(T))
443+
sensitivity_data.parameter_output_backward[p] =
444+
value + term.coefficient * d_cte
445+
end
446+
# Degree-2: pp terms
447+
for term in pp_terms
448+
p_1 = term.variable_1
449+
p_2 = term.variable_2
450+
value_1 = get!(sensitivity_data.parameter_output_backward, p_1, zero(T))
451+
value_2 = get!(sensitivity_data.parameter_output_backward, p_2, zero(T))
452+
sensitivity_data.parameter_output_backward[p_1] =
453+
value_1 +
454+
term.coefficient *
455+
d_cte *
456+
MOI.get(model, MOI.VariablePrimal(), p_2) /
457+
ifelse(p_1 === p_2, T(2), T(1))
458+
sensitivity_data.parameter_output_backward[p_2] =
459+
value_2 +
460+
term.coefficient *
461+
d_cte *
462+
MOI.get(model, MOI.VariablePrimal(), p_1) /
463+
ifelse(p_1 === p_2, T(2), T(1))
464+
end
465+
# Degree-2: pv terms
466+
for term in pv_terms
467+
p = term.variable_1
468+
v = term.variable_2
469+
value = get!(sensitivity_data.parameter_output_backward, p, zero(T))
470+
sensitivity_data.parameter_output_backward[p] =
471+
value + term.coefficient * JuMP.coefficient(grad_pf, v)
472+
end
473+
return
474+
end
475+
269476
function empty_input_sensitivities!(model::POI.Optimizer{T}) where {T}
270477
empty_input_sensitivities!(model.optimizer)
271478
model.ext[_SENSITIVITY_DATA] = SensitivityData{T}()
@@ -287,6 +494,8 @@ function forward_differentiate!(model::POI.Optimizer{T}) where {T}
287494
_affine_objective_set_forward!(model)
288495
elseif obj_type <: POI.ParametricQuadraticFunction
289496
_quadratic_objective_set_forward!(model)
497+
elseif obj_type <: POI.ParametricCubicFunction
498+
_cubic_objective_set_forward!(model)
290499
end
291500
forward_differentiate!(model.optimizer)
292501
return
@@ -511,6 +720,8 @@ function reverse_differentiate!(model::POI.Optimizer)
511720
_affine_objective_get_reverse!(model)
512721
elseif obj_type <: POI.ParametricQuadraticFunction
513722
_quadratic_objective_get_reverse!(model)
723+
elseif obj_type <: POI.ParametricCubicFunction
724+
_cubic_objective_get_reverse!(model)
514725
end
515726
return
516727
end

0 commit comments

Comments
 (0)