Fix bounds errors in default x_abschange and x_relchange (#1222)

* Fix bounds errors in default `x_abschange` and `x_relchange`

* Add code comment

Author: devmotion

src/Optim.jl

@@ -146,7 +146,6 @@ include("multivariate/precon.jl") # preconditioning functionality
 
 # utilities
 include("utilities/generic.jl") # generic utilities
-include("utilities/maxdiff.jl") # find largest difference
 include("utilities/update.jl")  # trace code
 
 # Unconstrained optimization

src/utilities/assess_convergence.jl

@@ -4,9 +4,26 @@ f_relchange(state::AbstractOptimizerState) = f_relchange(state.f_x, state.f_x_pr
 f_relchange(f_x, f_x_previous) = abs(f_x - f_x_previous) / abs(f_x)
 
 x_abschange(state::AbstractOptimizerState) = x_abschange(state.x, state.x_previous)
-x_abschange(x, x_previous) = maxdiff(x, x_previous)
+x_abschange(x::AbstractArray{<:Number}, x_previous::AbstractArray{<:Number}) = Linfdist(x, x_previous)
 x_relchange(state::AbstractOptimizerState) = x_relchange(state.x, state.x_previous)
-x_relchange(x, x_previous) = maxdiff(x, x_previous) / Base.maximum(abs, x) # Base.maximum !== maximum
+x_relchange(x::AbstractArray{<:Number}, x_previous::AbstractArray{<:Number}) = Linfdist(x, x_previous) / Base.maximum(abs, x) # Base.maximum !== maximum
+
+# Copied and adapted from https://github.com/JuliaStats/StatsBase.jl/blob/d70c4a203177c79d851c1cdca450bc6dbd2a4683/src/deviation.jl#L100-L110
+#    Linfdist(a, b)
+#
+# Compute the L∞ distance, also called the Chebyshev distance, between
+# two arrays: ``\\max_{1≤i≤n} |a_i - b_i|``.
+# Efficient equivalent of `maximum(abs, a - b)`.
+function Linfdist(x::AbstractArray{<:Number}, y::AbstractArray{<:Number})
+    n = length(x)
+    length(y) == n || throw(DimensionMismatch("Inconsistent array lengths."))
+    if iszero(n)
+        return zero(abs(zero(eltype(x)) - zero(eltype(y))))
+    else
+        broadcasted = Broadcast.broadcasted((xi, yi) -> abs(xi - yi), vec(x), vec(y))
+        return Base.maximum(Broadcast.instantiate(broadcasted)) # Base.maximum !== maximum
+    end
+end
 
 g_residual(state::NelderMeadState) = state.nm_x
 g_residual(state::ZerothOrderState) = oftype(state.f_x, NaN)

src/utilities/maxdiff.jl

@@ -1,17 +1,16 @@
-mutable struct DummyState <: Optim.AbstractOptimizerState
-    x::Any
-    x_previous::Any
-    f_x::Any
-    f_x_previous::Any
-    g_x::Any
+mutable struct DummyState{TX,TF,TG} <: Optim.AbstractOptimizerState
+    x::TX
+    x_previous::TX
+    f_x::TF
+    f_x_previous::TF
+    g_x::TG
 end
 
-mutable struct DummyStateZeroth <: Optim.ZerothOrderState
-    x::Any
-    x_previous::Any
-    f_x::Any
-    f_x_previous::Any
-    g_x::Any
+mutable struct DummyStateZeroth{TX,TF} <: Optim.ZerothOrderState
+    x::TX
+    x_previous::TX
+    f_x::TF
+    f_x_previous::TF
 end
 
 mutable struct DummyOptions
@@ -81,7 +80,7 @@ mutable struct DummyMethodZeroth <: Optim.ZerothOrderOptimizer end
     @test Optim.initial_convergence(ds, opt) == (true, false)
 
     # Zeroth order methods have no gradient -> returns false by default
-    ds = DummyStateZeroth(x1, x0, f1, f0, g)
+    ds = DummyStateZeroth(x1, x0, f1, f0)
     dm = DummyMethodZeroth()
 
     x = ones(2)
@@ -92,4 +91,18 @@ mutable struct DummyMethodZeroth <: Optim.ZerothOrderOptimizer end
 
     # should check all other methods as well
 
+    for T in (Float32, Float64)
+        ds = DummyState(T[-1.3, 2.5, -4.1], T[-1.1, 2.8, -4.0], f_x, f0, zeros(3))
+        @test @inferred(Optim.x_abschange(ds))::T ≈ 0.3
+        @test iszero((s -> @allocated(Optim.x_abschange(s)))(ds))
+        @test @inferred(Optim.x_relchange(ds))::T ≈ 0.3 / 4.1
+        @test iszero((s -> @allocated(Optim.x_relchange(s)))(ds))
+
+        # Special case: Empty state
+        ds = DummyState(T[], T[], f_x, f0, empty(g_x))
+        @test iszero(@inferred(Optim.x_abschange(ds))::T)
+        @test iszero((s -> @allocated(Optim.x_abschange(s)))(ds))
+        @test isnan(@inferred(Optim.x_relchange(ds))::T)
+        @test iszero((s -> @allocated(Optim.x_relchange(s)))(ds))
+    end
 end