@@ -42,6 +42,7 @@ mutable struct EqQPSolver <: MOI.AbstractOptimizer
4242 # Reverse input seeds
4343 rev_dx:: Dict{MOI.VariableIndex,Float64}
4444 rev_dy:: Dict{EQ_CI,Float64}
45+ rev_dobj:: Float64
4546 # Reverse results
4647 dc:: Vector{Float64}
4748 dQ:: Matrix{Float64}
@@ -53,6 +54,7 @@ mutable struct EqQPSolver <: MOI.AbstractOptimizer
5354 # Forward results
5455 dx_fwd:: Vector{Float64}
5556 dnu_fwd:: Vector{Float64}
57+ fwd_obj_sensitivity:: Float64
5658 diff_time:: Float64
5759
5860 function EqQPSolver ()
@@ -71,6 +73,7 @@ mutable struct EqQPSolver <: MOI.AbstractOptimizer
7173 MOI. OPTIMIZE_NOT_CALLED,
7274 Dict {MOI.VariableIndex,Float64} (),
7375 Dict {EQ_CI,Float64} (),
76+ 0.0 ,
7477 Float64[],
7578 zeros (0 , 0 ),
7679 zeros (0 , 0 ),
@@ -80,6 +83,7 @@ mutable struct EqQPSolver <: MOI.AbstractOptimizer
8083 Float64[],
8184 Float64[],
8285 0.0 ,
86+ 0.0 ,
8387 )
8488 end
8589end
@@ -103,6 +107,7 @@ function MOI.empty!(s::EqQPSolver)
103107 s. status = MOI. OPTIMIZE_NOT_CALLED
104108 empty! (s. rev_dx)
105109 empty! (s. rev_dy)
110+ s. rev_dobj = 0.0
106111 s. dc = Float64[]
107112 s. dQ = zeros (0 , 0 )
108113 s. dA = zeros (0 , 0 )
@@ -111,6 +116,7 @@ function MOI.empty!(s::EqQPSolver)
111116 empty! (s. fwd_constraints)
112117 s. dx_fwd = Float64[]
113118 s. dnu_fwd = Float64[]
119+ s. fwd_obj_sensitivity = 0.0
114120 s. diff_time = 0.0
115121 return
116122end
@@ -386,6 +392,11 @@ function MOI.set(
386392 return
387393end
388394
395+ function MOI. set (s:: EqQPSolver , :: DiffOpt.ReverseObjectiveSensitivity , value)
396+ s. rev_dobj = value
397+ return
398+ end
399+
389400# Trigger backward differentiation
390401function MOI. set (s:: EqQPSolver , :: DiffOpt.BackwardDifferentiate , :: Nothing )
391402 s. diff_time = @elapsed _backward_differentiate! (s)
@@ -401,6 +412,11 @@ function _backward_differentiate!(s::EqQPSolver)
401412 for (ci, val) in s. rev_dy
402413 rhs[n+ ci. value] = val
403414 end
415+ # If dobj seed is nonzero, add the gradient of the objective w.r.t. x
416+ # dobj * (Qx + c) contributes to the rhs for the primal variables
417+ if ! iszero (s. rev_dobj)
418+ rhs[1 : n] .+ = s. rev_dobj .* (s. Q * s. x .+ s. c)
419+ end
404420 adj = s. kkt_factor \ rhs
405421 dx_adj = adj[1 : n]
406422 dnu_adj = adj[(n+ 1 ): end ]
@@ -411,6 +427,7 @@ function _backward_differentiate!(s::EqQPSolver)
411427 # Clear seeds
412428 empty! (s. rev_dx)
413429 empty! (s. rev_dy)
430+ s. rev_dobj = 0.0
414431 return
415432end
416433
@@ -473,6 +490,9 @@ function _forward_differentiate!(s::EqQPSolver)
473490 sol = s. kkt_factor \ rhs
474491 s. dx_fwd = sol[1 : n]
475492 s. dnu_fwd = sol[(n+ 1 ): end ]
493+ # Compute forward objective sensitivity: d(obj)/dt = (Qx + c)'dx + x'Q*dx/2 + d_c'x
494+ # Actually: obj = 1/2 x'Qx + c'x, so dobj/dt = (Qx+c)'dx_fwd + d_c'x
495+ s. fwd_obj_sensitivity = dot (s. Q * s. x .+ s. c, s. dx_fwd) + dot (d_c, s. x)
476496 # Clear inputs
477497 s. fwd_objective = nothing
478498 empty! (s. fwd_constraints)
@@ -492,13 +512,17 @@ function MOI.get(s::EqQPSolver, ::DiffOpt.ForwardConstraintDual, ci::EQ_CI)
492512 return s. dnu_fwd[ci. value]
493513end
494514
515+ # ForwardObjectiveSensitivity
516+ MOI. get (s:: EqQPSolver , :: DiffOpt.ForwardObjectiveSensitivity ) = s. fwd_obj_sensitivity
517+
495518# DifferentiateTimeSec
496519MOI. get (s:: EqQPSolver , :: DiffOpt.DifferentiateTimeSec ) = s. diff_time
497520
498521# empty_input_sensitivities!
499522function DiffOpt. empty_input_sensitivities! (s:: EqQPSolver )
500523 empty! (s. rev_dx)
501524 empty! (s. rev_dy)
525+ s. rev_dobj = 0.0
502526 s. fwd_objective = nothing
503527 empty! (s. fwd_constraints)
504528 return
@@ -834,6 +858,90 @@ function test_three_variable_problem()
834858 end
835859end
836860
861+ # ── Test reverse with dobj seed (ReverseObjectiveSensitivity) ─────────────
862+
863+ function test_reverse_dobj_seed ()
864+ model, x1, x2, c1 = _setup_model ()
865+
866+ # Set dobj = 1.0: sensitivity of the objective value w.r.t. problem data
867+ MOI. set (model, DiffOpt. ReverseObjectiveSensitivity (), 1.0 )
868+ DiffOpt. reverse_differentiate! (model)
869+
870+ # x* = [0, 1], nu* = [-3], Q = I, c = [3, 2]
871+ # grad_obj w.r.t. x = Q*x + c = [3, 3]
872+ # KKT solve with rhs = [3, 3, 0] gives the adjoint
873+ K = [1.0 0.0 1.0 ; 0.0 1.0 1.0 ; 1.0 1.0 0.0 ]
874+ rhs = [3.0 , 3.0 , 0.0 ]
875+ adj = K \ rhs
876+ dc_expected = - adj[1 : 2 ]
877+
878+ dobj = MOI. get (model, DiffOpt. ReverseObjectiveFunction ())
879+ @test JuMP. coefficient (dobj, x1) ≈ dc_expected[1 ] atol = ATOL
880+ @test JuMP. coefficient (dobj, x2) ≈ dc_expected[2 ] atol = ATOL
881+ end
882+
883+ # ── Test forward objective sensitivity ────────────────────────────────────
884+
885+ function test_forward_objective_sensitivity ()
886+ model, x1, x2, c1 = _setup_model ()
887+
888+ # Perturb dc = [1, 0]
889+ MOI. set (model, DiffOpt. ForwardObjectiveFunction (), 1.0 * x1)
890+ DiffOpt. forward_differentiate! (model)
891+
892+ # x* = [0, 1], Q = I, c = [3, 2]
893+ # K * [dx; dnu] = [-dc; 0] = [-1; 0; 0]
894+ K = [1.0 0.0 1.0 ; 0.0 1.0 1.0 ; 1.0 1.0 0.0 ]
895+ fwd = K \ [- 1.0 , 0.0 , 0.0 ]
896+ # dobj/dt = (Qx + c)'dx + dc'x = [3,3]'*dx + [1,0]'*[0,1]
897+ expected = dot ([3.0 , 3.0 ], fwd[1 : 2 ]) + dot ([1.0 , 0.0 ], [0.0 , 1.0 ])
898+
899+ @test MOI. get (model, DiffOpt. ForwardObjectiveSensitivity ()) ≈ expected atol = ATOL
900+ end
901+
902+ # ── Test forward twice (re-differentiation in forward mode) ──────────────
903+
904+ function test_forward_twice ()
905+ model, x1, x2, c1 = _setup_model ()
906+
907+ K = [1.0 0.0 1.0 ; 0.0 1.0 1.0 ; 1.0 1.0 0.0 ]
908+
909+ # First: perturb dc = [1, 0]
910+ MOI. set (model, DiffOpt. ForwardObjectiveFunction (), 1.0 * x1)
911+ DiffOpt. forward_differentiate! (model)
912+
913+ fwd1 = K \ [- 1.0 , 0.0 , 0.0 ]
914+ @test MOI. get (model, DiffOpt. ForwardVariablePrimal (), x1) ≈ fwd1[1 ] atol = ATOL
915+
916+ # Second: perturb dc = [0, 1]
917+ DiffOpt. empty_input_sensitivities! (model)
918+ MOI. set (model, DiffOpt. ForwardObjectiveFunction (), 1.0 * x2)
919+ DiffOpt. forward_differentiate! (model)
920+
921+ fwd2 = K \ [0.0 , - 1.0 , 0.0 ]
922+ @test MOI. get (model, DiffOpt. ForwardVariablePrimal (), x1) ≈ fwd2[1 ] atol = ATOL
923+ @test MOI. get (model, DiffOpt. ForwardVariablePrimal (), x2) ≈ fwd2[2 ] atol = ATOL
924+ end
925+
926+ # ── Test ReverseConstraintSet (not applicable without parameters, but test the getter path) ──
927+
928+ function test_reverse_constraint_set ()
929+ model, x1, x2, c1 = _setup_model ()
930+
931+ MOI. set (model, DiffOpt. ReverseVariablePrimal (), x1, 1.0 )
932+ DiffOpt. reverse_differentiate! (model)
933+
934+ dc, db, dA = _expected_reverse ([1.0 , 0.0 ], [0.0 ])
935+
936+ # ReverseConstraintFunction returns dA and db
937+ dcon = MOI. get (model, DiffOpt. ReverseConstraintFunction (), c1)
938+ @test MOI. constant (dcon) ≈ - db[1 ] atol = ATOL
939+
940+ # standard_form should convert the lazy function
941+ sf = DiffOpt. standard_form (dcon)
942+ @test sf isa MOI. ScalarAffineFunction
943+ end
944+
837945# ── Test empty_input_sensitivities! ──────────────────────────────────────────
838946
839947function test_empty_input_sensitivities ()
0 commit comments