Remove optim in favor of independent golden section implementation (#84)

* remove Optim dependency

* test optimizer

* minor reorg

* Apply suggestions from code review

Co-authored-by: Tamas K. Papp <[email protected]>

* add docstring for optimize

Co-authored-by: Tamas K. Papp <[email protected]>

Author: piever

.gitignore

@@ -1 +1,2 @@
 examples/.ipynb_checkpoints/*
+Manifest.toml

Project.toml

@@ -8,12 +8,10 @@ Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
-Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 
 [compat]
 Interpolations = "≥ 0.9"
-Optim = "≥ 0.16"
 julia = "^1"
 
 [extras]

src/KernelDensity.jl

@@ -3,7 +3,6 @@ module KernelDensity
 using DocStringExtensions: TYPEDEF, FIELDS
 using StatsBase
 using Distributions
-using Optim
 using Interpolations
 
 import StatsBase: RealVector, RealMatrix

src/univariate.jl

@@ -173,6 +173,76 @@ function kde(data::RealVector; bandwidth=default_bandwidth(data), kernel=Normal,
     kde(data,dist;boundary=boundary,npoints=npoints,weights=weights)
 end
 
+"""
+    optimize(f, x_lower, x_upper; iterations=1000, rel_tol=nothing, abs_tol=nothing)
+
+Minimize the function `f` in the interval `x_lower..x_upper`, using the
+[golden-section search](https://en.wikipedia.org/wiki/Golden-section_search).
+Return an approximate minimum `x̃` or error if such approximate minimum cannot be found.
+
+This algorithm assumes that `-f` is unimodal on the interval `x_lower..x_upper`,
+that is to say, there exists a unique `x` in `x_lower..x_upper` such that `f` is
+decreasing on `x_lower..x` and increasing on `x..x_upper`.
+
+`rel_tol` and `abs_tol` determine the relative and absolute tolerance, that is
+to say, the returned value `x̃` should differ from the actual minimum `x` at most
+`abs_tol + rel_tol * abs(x̃)`.
+If not manually specified, `rel_tol` and `abs_tol` default to `sqrt(eps(T))` and
+`eps(T)` respectively, where `T` is the floating point type of `x_lower` and `x_upper`.
+
+`iterations` determines the maximum number of iterations allowed before convergence.
+
+This is a private, unexported function, used internally to select the optimal bandwidth
+automatically.
+"""
+function optimize(f, x_lower, x_upper; iterations=1000, rel_tol=nothing, abs_tol=nothing)
+
+    if x_lower > x_upper
+        error("x_lower must be less than x_upper")
+    end
+
+    T = promote_type(typeof(x_lower/1), typeof(x_upper/1))
+    rtol = something(rel_tol, sqrt(eps(T)))
+    atol = something(abs_tol, eps(T))
+    
+    function midpoint_and_convergence(lower, upper)
+        midpoint = (lower + upper) / 2
+        tol = atol + rtol * midpoint
+        midpoint, (upper - lower) <= 2tol
+    end
+    
+    invphi::T = 0.5 * (sqrt(5) - 1)
+    invphisq::T = 0.5 * (3 - sqrt(5))
+    
+    a::T, b::T = x_lower, x_upper
+    h = b - a
+    c = a + invphisq * h
+    d = a + invphi * h
+    
+    fc, fd = f(c), f(d)
+    
+    for _ in 1:1000
+        h *= invphi
+        if fc < fd
+            m, converged = midpoint_and_convergence(a, d)
+            converged && return m
+            b = d
+            d, fd = c, fc
+            c = a + invphisq * h
+            fc = f(c)
+        else
+            m, converged = midpoint_and_convergence(c, b)
+            converged && return m
+            a = c
+            c, fc = d, fd
+            d = a + invphi * h
+            fd = f(d)
+        end
+    end
+
+    error("Reached maximum number of iterations without convergence.")
+end
+
 # Select bandwidth using least-squares cross validation, from:
 #   Density Estimation for Statistics and Data Analysis
 #   B. W. Silverman (1986)
@@ -194,7 +264,7 @@ function kde_lscv(data::RealVector, midpoints::R;
     c = -twoπ/(step(k.x)*K)
     hlb, hub = bandwidth_range
 
-    opt = Optim.optimize(hlb, hub) do h
+    minimizer = optimize(hlb, hub) do h
         dist = kernel_dist(kernel, h)
         ψ = 0.0
         for j = 1:length(ft2)-1
@@ -204,7 +274,7 @@ function kde_lscv(data::RealVector, midpoints::R;
         ψ*step(k.x)/K + pdf(dist,0.0)/ndata
     end
 
-    dist = kernel_dist(kernel, Optim.minimizer(opt))
+    dist = kernel_dist(kernel, minimizer)
     for j = 0:length(ft)-1
         ft[j+1] *= cf(dist, j*c)
     end

test/univariate.jl

@@ -63,3 +63,9 @@ end
 k11 = kde([0.0, 1.], r, bandwidth=1, weights=[0,1])
 k12 = kde([1.], r, bandwidth=1)
 @test k11.density ≈ k12.density
+
+rel_tol = sqrt(eps(Float64))
+abs_tol = eps(Float64)
+
+minimizer = @inferred KernelDensity.optimize(x -> (x - 1)^2, 0, 2)
+@test abs(minimizer - 1) <= abs_tol + rel_tol * abs(minimizer)