Hopefully more robust convergence estimation

Author: jagot

Project.toml

@@ -37,8 +37,9 @@ ArnoldiMethod = "ec485272-7323-5ecc-a04f-4719b315124d"
 PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+RollingFunctions = "b0e4dd01-7b14-53d8-9b45-175a3e362653"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
 
 [targets]
-test = ["Test", "ArnoldiMethod", "PrettyTables", "Random", "ProgressMeter", "UnicodePlots"]
+test = ["Test", "ArnoldiMethod", "PrettyTables", "Random", "RollingFunctions", "ProgressMeter", "UnicodePlots"]

test/derivative_accuracy_utils.jl

@@ -5,6 +5,8 @@ using PrettyTables
 using ProgressMeter
 using UnicodePlots
 
+using RollingFunctions
+
 function test_derivatives(R::AbstractQuasiMatrix,
                           f::Function, g::Function, h::Function)
     r = axes(R,1)
@@ -38,18 +40,31 @@ function test_derivatives(R::AbstractQuasiMatrix,
     R,fv,gv,hv,hv′,δgv,δhv,δhv′
 end
 
+function estimate_convergence_rate(ρ, ϵ)
+    i = findall(!iszero, ϵ)
+    x = log10.(ρ[i])
+    y = log10.(abs.(ϵ[i]))
+
+    f = (x,y) -> ([x ones(length(x))] \ y)[1]
+
+    # Find all slopes for consecutive groups of three errors
+    slopes = rolling(f, x, y, 3)
+    e = findfirst(s -> s < 0, slopes)
+    if e !== nothing && e > 5
+        slopes = slopes[1:max(1,e-1)]
+    end
+
+    maximum(slopes)
+end
+
 function error_slopes(hs, ϵ, names, verbosity=0)
     loghs = log10.(hs)
     logϵ = log10.(abs.(ϵ))
     n = size(logϵ,2)
 
     slopes = zeros(n)
     for j = 1:n
-        # To avoid the effect of round-off errors on the order
-        # estimation.
-        i = argmin(abs.(ϵ[:,j]))
-
-        slopes[j] = ([loghs[1:i] ones(i)] \ logϵ[1:i,j])[1]
+        slopes[j] = estimate_convergence_rate(hs, ϵ[:,j])
     end
 
     if verbosity > 0

test/fedvr/derivatives.jl

@@ -163,12 +163,12 @@ end
     # The convergence rate should be order - 1 (polynomial order =
     # polynomial degree + 1), but since the error slope fitting is a
     # bit error prone, we require that it is greater than order - 2.5.
-    @test_broken all(slopes[:,1] .> oo .- 2.5)
+    @test all(slopes[:,1] .> oo .- 2.5)
     @test all(slopes[:,1] .> oo .- 3.5)
     # Since the approximation to h is calculated by computing ∇∇f, we
     # lose one order extra, compared to ∇²f.
-    @test_broken all(slopes[:,2] .> oo .- 3.5)
+    @test all(slopes[:,2] .> oo .- 3.5)
     @test all(slopes[:,2] .> oo .- 4.5)
-    @test_broken all(slopes[:,3] .> oo .- 2.5)
+    @test all(slopes[:,3] .> oo .- 2.5)
     @test all(slopes[:,3] .> oo .- 4.0)
 end