New attempt with Optim.jl v2.2

Project.toml

@@ -40,7 +40,7 @@ Makie = "0.24"
 Meshes = "0.57"
 NearestNeighbors = "0.4"
 OhMyThreads = "0.5 - 0.8"
-Optim = "2.0"
+Optim = "2.2"
 PrecompileTools = "1.2"
 Random = "1.10"
 Setfield = "1.0"

src/fitting.jl

@@ -132,8 +132,7 @@ _ustrip(u, x::Quantity) = ustrip(u, x)
 _weights(f, x, n) = isnothing(f) ? n / sum(n) : map(xᵢ -> ustrip(f(xᵢ)), x)
 
 function _optimize(J, L, λ, l, u, θₒ)
-  # https://github.com/JuliaNLSolvers/Optim.jl/issues/1228
-  s = Optim.optimize(θ -> J(θ) + λ * L(θ), l, u, θₒ, Fminbox(), Optim.Options(show_warnings=false))
+  s = Optim.optimize(θ -> J(θ) + λ * L(θ), l, u, θₒ, LBFGSB())
   ϵ = Optim.minimum(s)
   θ = Optim.minimizer(s)
   θ, ϵ

src/fitting/transiograms.jl

@@ -44,7 +44,7 @@ function _fit(
   # objective function
   function J(θ)
     b = ball(θ[1])
-    p = _proportions(θ[2:end])
+    p = _normalize(θ[2:end])
     τ = T(b, proportions=p)
     mat(i) = getindex.(Y, i)
     err(i) = sum(abs2, τ(x′[i]) - mat(i))
@@ -59,28 +59,27 @@ function _fit(
 
   # maximum range and proportions
   xmax = maximum(x′)
-  ymax = _ones
   rmax = isnothing(maxrange′) ? xmax : maxrange′
-  pmax = isnothing(maxproportions) ? ymax : maxproportions
-
-  # initial guess
-  rₒ = isnothing(range′) ? rmax / 3 : range′
-  pₒ = isnothing(proportions) ? 0.95 .* pmax : proportions
-  θₒ = [rₒ, pₒ...]
+  pmax = isnothing(maxproportions) ? _ones : maxproportions
 
   # box constraints
   δ = 1e-8
   rₗ, rᵤ = isnothing(range′) ? (zero(rmax), rmax) : (range′ - δ, range′ + δ)
-  pₗ, pᵤ = isnothing(proportions) ? (_zeros, pmax) : (proportions .- δ, proportions .+ δ)
+  pₗ, pᵤ = isnothing(proportions) ? (_zeros .+ δ, pmax) : (proportions .- δ, proportions .+ δ)
   l = [rₗ, pₗ...]
   u = [rᵤ, pᵤ...]
 
+  # initial guess
+  rₒ = isnothing(range′) ? rmax / 3 : range′
+  pₒ = isnothing(proportions) ? (pₗ .+ pᵤ) ./ 2 : proportions
+  θₒ = [rₒ, pₒ...]
+
   # solve optimization problem
   θ, ϵ = _optimize(J, L, λ, l, u, θₒ)
 
   # optimal transiogram (with units)
   b = ball(θ[1] * ux)
-  p = _proportions(θ[2:end])
+  p = _normalize(θ[2:end])
   τ = T(b, proportions=p)
 
   τ, ϵ
@@ -132,8 +131,8 @@ function _fit(
   # objective function
   function J(θ)
     b = ball(θ[1])
-    l = _lengths(θ[2:(k + 1)])
-    p = _proportions(θ[(k + 2):end])
+    l = _tuple(θ[2:(k + 1)])
+    p = _normalize(θ[(k + 2):end])
     τ = T(b, lengths=l, proportions=p)
     mat(i) = getindex.(Y, i)
     err(i) = sum(abs2, τ(x′[i]) - mat(i))
@@ -148,32 +147,31 @@ function _fit(
 
   # maximum range, proportions and lengths
   xmax = maximum(x′)
-  ymax = _ones
   rmax = isnothing(maxrange′) ? xmax : maxrange′
   lmax = isnothing(maxlengths′) ? ntuple(i -> xmax, k) : maxlengths′
-  pmax = isnothing(maxproportions) ? ymax : maxproportions
-
-  # initial guess
-  rₒ = isnothing(range′) ? rmax / 3 : range′
-  lₒ = isnothing(lengths′) ? lmax ./ 3 : lengths′
-  pₒ = isnothing(proportions) ? 0.95 .* pmax : proportions
-  θₒ = [rₒ, lₒ..., pₒ...]
+  pmax = isnothing(maxproportions) ? _ones : maxproportions
 
   # box constraints
   δ = 1e-8
   rₗ, rᵤ = isnothing(range′) ? (zero(rmax), rmax) : (range′ - δ, range′ + δ)
-  lₗ, lᵤ = isnothing(lengths′) ? (_zeros, lmax) : (lengths′ .- δ, lengths′ .+ δ)
-  pₗ, pᵤ = isnothing(proportions) ? (_zeros, pmax) : (proportions .- δ, proportions .+ δ)
+  lₗ, lᵤ = isnothing(lengths′) ? (_zeros .+ δ, lmax) : (lengths′ .- δ, lengths′ .+ δ)
+  pₗ, pᵤ = isnothing(proportions) ? (_zeros .+ δ, pmax) : (proportions .- δ, proportions .+ δ)
   l = [rₗ, lₗ..., pₗ...]
   u = [rᵤ, lᵤ..., pᵤ...]
 
+  # initial guess
+  rₒ = isnothing(range′) ? rmax / 3 : range′
+  lₒ = isnothing(lengths′) ? lmax ./ 3 : lengths′
+  pₒ = isnothing(proportions) ? (pₗ .+ pᵤ) ./ 2 : proportions
+  θₒ = [rₒ, lₒ..., pₒ...]
+
   # solve optimization problem
   θ, ϵ = _optimize(J, L, λ, l, u, θₒ)
 
   # optimal transiogram (with units)
   b = ball(θ[1] * ux)
-  l = _lengths(θ[2:(k + 1)] * ux)
-  p = _proportions(θ[(k + 2):end])
+  l = _tuple(θ[2:(k + 1)] * ux)
+  p = _normalize(θ[(k + 2):end])
   τ = T(b, lengths=l, proportions=p)
 
   τ, ϵ
@@ -205,10 +203,10 @@ end
 # HELPER FUNCTIONS
 # -----------------
 
-function _proportions(θ)
-  p = clamp.(θ, 0, 1)
-  s = sum(abs, p)
-  ntuple(i -> p[i] / s, length(p))
+function _normalize(p)
+  c = clamp.(p, 0, 1)
+  s = sum(abs, c)
+  ntuple(i -> c[i] / s, length(c))
 end
 
-_lengths(θ) = ntuple(i -> θ[i], length(θ))
+_tuple(θ) = ntuple(i -> θ[i], length(θ))