Fix LBFGS/BFGS callback receiving Dual numbers instead of scalar loss values

Fixes #1073

When using LBFGS or BFGS with bounds, Optim.jl wraps the optimizer in Fminbox,
which may use ForwardDiff internally for gradient computation. This resulted in
the callback receiving ForwardDiff.Dual numbers instead of scalar loss values,
causing incorrect (sometimes negative) values to be reported.

Changes:
- Use SciMLBase.value() to extract scalar values from potentially Dual-valued traces
- Updated all three _cb callback functions to use SciMLBase.value()
- Added comprehensive test case verifying callbacks receive correct scalar non-negative values
- No new dependencies required (SciMLBase already provides the functionality)

🤖 Generated with [Claude Code](https://claude.com/claude-code)

Co-Authored-By: Claude <[email protected]>

Author: ChrisRackauckas

lib/OptimizationOptimJL/Project.toml

@@ -3,7 +3,6 @@ uuid = "36348300-93cb-4f02-beb5-3c3902f8871e"
 authors = ["Vaibhav Dixit <[email protected]> and contributors"]
 version = "0.4.6"
 [deps]
-ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 OptimizationBase = "bca83a33-5cc9-4baa-983d-23429ab6bcbb"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -12,6 +11,7 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 
 [extras]
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
@@ -28,4 +28,4 @@ Reexport = "1.2"
 SciMLBase = "2.58"
 
 [targets]
-test = ["ModelingToolkit", "Random", "ReverseDiff", "Test", "Zygote"]
+test = ["ForwardDiff", "ModelingToolkit", "Random", "ReverseDiff", "Test", "Zygote"]

lib/OptimizationOptimJL/src/OptimizationOptimJL.jl

@@ -3,19 +3,9 @@ module OptimizationOptimJL
 using Reexport
 @reexport using Optim, OptimizationBase
 using SciMLBase, SparseArrays
-
-# Import ForwardDiff to handle Dual numbers in callbacks
-import ForwardDiff
-
 decompose_trace(trace::Optim.OptimizationTrace) = last(trace)
 decompose_trace(trace::Optim.OptimizationState) = trace
 
-# Extract scalar value from potentially Dual-valued trace values
-# This is needed because Optim.jl may use ForwardDiff internally for gradient computation,
-# resulting in Dual numbers in the trace, but callbacks should receive scalar values
-_scalar_value(x) = x  # Default case for regular numbers
-_scalar_value(x::ForwardDiff.Dual) = ForwardDiff.value(x)  # Extract value from Dual
-
 SciMLBase.allowsconstraints(::IPNewton) = true
 SciMLBase.allowsbounds(opt::Optim.AbstractOptimizer) = true
 SciMLBase.allowsbounds(opt::Optim.SimulatedAnnealing) = false
@@ -28,7 +18,8 @@ SciMLBase.requiresbounds(opt::Optim.SAMIN) = true
 end
 @static if isdefined(OptimizationBase, :supports_opt_cache_interface)
     OptimizationBase.supports_opt_cache_interface(opt::Optim.AbstractOptimizer) = true
-    OptimizationBase.supports_opt_cache_interface(opt::Union{Optim.Fminbox, Optim.SAMIN}) = true
+    OptimizationBase.supports_opt_cache_interface(opt::Union{
+        Optim.Fminbox, Optim.SAMIN}) = true
     OptimizationBase.supports_opt_cache_interface(opt::Optim.ConstrainedOptimizer) = true
 end
 function SciMLBase.requiresgradient(opt::Optim.AbstractOptimizer)
@@ -160,7 +151,8 @@ function SciMLBase.__solve(cache::OptimizationBase.OptimizationCache{
         metadata = trace_state.metadata
         θ = metadata[cache.opt isa Optim.NelderMead ? "centroid" : "x"]
         # Extract scalar value from potentially Dual-valued trace (issue #1073)
-        loss_val = _scalar_value(trace_state.value)
+        # Using SciMLBase.value to handle ForwardDiff.Dual numbers from Fminbox
+        loss_val = SciMLBase.value(trace_state.value)
         opt_state = OptimizationBase.OptimizationState(iter = trace_state.iteration,
             u = θ,
             p = cache.p,
@@ -283,7 +275,8 @@ function SciMLBase.__solve(cache::OptimizationBase.OptimizationCache{
             metadata["centroid"] :
             metadata["x"]
         # Extract scalar value from potentially Dual-valued trace (issue #1073)
-        loss_val = _scalar_value(trace_state.value)
+        # Using SciMLBase.value to handle ForwardDiff.Dual numbers from Fminbox
+        loss_val = SciMLBase.value(trace_state.value)
         opt_state = OptimizationBase.OptimizationState(iter = trace_state.iteration,
             u = θ,
             p = cache.p,
@@ -372,7 +365,8 @@ function SciMLBase.__solve(cache::OptimizationBase.OptimizationCache{
     function _cb(trace)
         metadata = decompose_trace(trace).metadata
         # Extract scalar value from potentially Dual-valued trace (issue #1073)
-        loss_val = _scalar_value(trace.value)
+        # Using SciMLBase.value to handle ForwardDiff.Dual numbers from Fminbox
+        loss_val = SciMLBase.value(trace.value)
         opt_state = OptimizationBase.OptimizationState(iter = trace.iteration,
             u = metadata["x"],
             p = cache.p,

lib/OptimizationOptimJL/test/runtests.jl

@@ -150,7 +150,8 @@ end
         G[1] = -2.0 * (1.0 - x[1]) - 400.0 * (x[2] - x[1]^2) * x[1]
         G[2] = 200.0 * (x[2] - x[1]^2)
     end
-    optprob = OptimizationFunction((x, p) -> -rosenbrock(x, p), OptimizationBase.AutoZygote(),
+    optprob = OptimizationFunction(
+        (x, p) -> -rosenbrock(x, p), OptimizationBase.AutoZygote(),
         grad = g!)
     prob = OptimizationProblem(optprob, x0, _p; sense = OptimizationBase.MaxSense)
     sol = solve(prob, BFGS())
@@ -171,7 +172,8 @@ end
     @test 10 * sol.objective < l1
 
     prob = OptimizationProblem(
-        optprob, x0, _p; sense = OptimizationBase.MaxSense, lb = [-1.0, -1.0], ub = [0.8, 0.8])
+        optprob, x0, _p; sense = OptimizationBase.MaxSense, lb = [-1.0, -1.0], ub = [
+            0.8, 0.8])
     sol = solve(prob, BFGS())
     @test 10 * sol.objective < l1
 
@@ -205,9 +207,8 @@ end
         # Create a non-negative loss function (sum of squares)
         loss_vals = Float64[]
         function test_callback(state, loss_val)
-            # Verify loss_val is a scalar, not a Dual number
-            @test loss_val isa Real
-            @test !(loss_val isa ForwardDiff.Dual)
+            # Verify loss_val is a scalar Float64, not a Dual number
+            @test loss_val isa Float64
             # For a sum-of-squares loss, values should be non-negative
             push!(loss_vals, loss_val)
             return false