Merge pull request #95 from rmsrosa/convergence_recipe

Δt-convergence plot recipe for WorkPrecisionSet

Author: ChrisRackauckas

.gitignore

@@ -2,3 +2,7 @@
 *.jl.*.cov
 *.jl.mem
 Manifest.toml
+
+# vscode
+.vscode
+.vscode/*

Project.toml

@@ -35,7 +35,8 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SDEProblemLibrary = "c72e72a9-a271-4b2b-8966-303ed956772e"
 StochasticDelayDiffEq = "29a0d76e-afc8-11e9-03a4-eda52ae4b960"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["DelayDiffEq", "DDEProblemLibrary", "ODEProblemLibrary", "SDEProblemLibrary", "OrdinaryDiffEq", "ParameterizedFunctions", "Random", "StochasticDiffEq", "StochasticDelayDiffEq", "Test"]
+test = ["DelayDiffEq", "DDEProblemLibrary", "ODEProblemLibrary", "SDEProblemLibrary", "OrdinaryDiffEq", "ParameterizedFunctions", "Random", "StochasticDiffEq", "StochasticDelayDiffEq", "Plots", "Test"]

src/plotrecipes.jl

@@ -37,22 +37,69 @@ end
     wp.errors, wp.times
 end
 
-@recipe function f(wp_set::WorkPrecisionSet)
-    seriestype --> :path
-    linewidth --> 3
-    yguide --> "Time (s)"
-    xguide --> "Error"
-    xscale --> :log10
-    yscale --> :log10
-    marker --> :auto
-    errors = Vector{Any}(undef, 0)
-    times = Vector{Any}(undef, 0)
-    for i in 1:length(wp_set)
-        push!(errors, wp_set[i].errors)
-        push!(times, wp_set[i].times)
+@recipe function f(wp_set::WorkPrecisionSet; view = :benchmark, color = nothing)
+    if view == :benchmark
+        seriestype --> :path
+        linewidth --> 3
+        yguide --> "Time (s)"
+        xguide --> "Error"
+        xscale --> :log10
+        yscale --> :log10
+        marker --> :auto
+        errors = Vector{Any}(undef, 0)
+        times = Vector{Any}(undef, 0)
+        for i in 1:length(wp_set)
+            push!(errors, wp_set[i].errors)
+            push!(times, wp_set[i].times)
+        end
+        label --> reshape(wp_set.names, 1, length(wp_set))
+        return errors, times
+    elseif view == :dt_convergence
+        idts = filter(i -> haskey(wp_set.setups[i], :dts), 1:length(wp_set))
+        length(idts) > 0 ||
+            throw(ArgumentError("Convergence with respect to Δt requires runs with fixed time steps"))
+        dts = Vector{Any}(undef, 0)
+        errors = Vector{Any}(undef, 0)
+        ps = Vector{Any}(undef, 0)
+        convs = Vector{Any}(undef, 0)
+        for i in idts
+            push!(dts, wp_set.setups[i][:dts])
+            push!(errors, wp_set[i].errors)
+            lc, p = [one.(dts[end]) log.(dts[end])] \ log.(errors[end])
+            push!(ps, p)
+            push!(convs, exp(lc) * dts[end] .^ p)
+        end
+        names = wp_set.names[idts] .*
+                map(p -> " (Δtᵖ order p=$(round(p, sigdigits=2)))", ps)
+        if color === nothing
+            color = reshape(1:length(idts), 1, :)
+        end
+        @series begin
+            seriestype --> :path
+            linestyle --> :dash
+            linewidth --> 1
+            color --> color
+            xscale --> :log10
+            yscale --> :log10
+            markersize --> 0
+            label --> nothing
+            return dts, convs
+        end
+        @series begin
+            seriestype --> :path
+            linewidth --> 3
+            color --> color
+            yguide --> "Error"
+            xguide --> "Δt"
+            xscale --> :log10
+            yscale --> :log10
+            marker --> :auto
+            label --> reshape(names, 1, length(idts))
+            return dts, errors
+        end
+    else
+        throw(ArgumentError("view argument `$view` not implemented"))
     end
-    label --> reshape(wp_set.names, 1, length(wp_set))
-    errors, times
 end
 
 @recipe function f(tab::ODERKTableau; dx = 1 / 100, dy = 1 / 100, order_star = false,

test/plotrecipes_tests.jl

@@ -0,0 +1,63 @@
+using Test
+using OrdinaryDiffEq, DiffEqDevTools, Plots
+
+using Random
+Random.seed!(123)
+
+gr()
+
+# Linear ODE system
+f_rode_lin = function (du, u, p, t)
+    @inbounds for i in eachindex(u)
+        du[i] = cos(t) * u[i]
+    end
+end
+
+f_rode_lin_analytic = function (u₀, p, t)
+    u₀ * exp(sin(t))
+end
+
+tspan = (0.0, 10.0)
+prob = ODEProblem(ODEFunction(f_rode_lin, analytic = f_rode_lin_analytic), rand(10, 10),
+                  tspan)
+
+abstols = 1.0 ./ 10.0 .^ (3:8)
+reltols = 1.0 ./ 10.0 .^ (0:5)
+
+setups = [Dict(:alg => DP5())
+          Dict(:alg => Tsit5())]
+wp_names = ["DP5", "Tsit5"]
+wp = WorkPrecisionSet(prob, abstols, reltols, setups, names = wp_names,
+                      save_everystep = false,
+                      numruns = 100)
+
+plt = @test_nowarn plot(wp)
+@test plt[1][1][:x] ≈ wp[1].errors
+@test plt[1][2][:x] ≈ wp[2].errors
+@test plt[1][1][:label] == wp_names[1]
+@test plt[1][2][:label] == wp_names[2]
+@test_nowarn plot(wp, color = [1 2])
+@test_nowarn plot(wp, color = :blue)
+@test_throws ArgumentError plot(wp, view = :dt_convergence)
+
+dts = 1.0 ./ 2.0 .^ ((1:length(reltols)) .+ 1)
+setups = [Dict(:alg => Euler(), :dts => dts)
+          Dict(:alg => Heun(), :dts => dts)
+          Dict(:alg => Tsit5(), :dts => dts, :adaptive => false)
+          Dict(:alg => Tsit5())]
+wp_names = ["Euler", "Heun", "Tsit5 fixed step", "Tsit5 adaptive"]
+wp = WorkPrecisionSet(prob, abstols, reltols, setups, names = wp_names,
+                      save_everystep = false,
+                      numruns = 100)
+
+plt = @test_nowarn plot(wp)
+@test all(plt[1][i][:x] ≈ wp[i].errors for i in 1:4)
+@test all(plt[1][i][:label] == wp_names[i] for i in 1:4)
+
+plt = @test_nowarn plot(wp, view = :dt_convergence, legend = :bottomright)
+@test all(plt[1][i][:x] == plt[1][i + 3][:x] == dts == wp.setups[i][:dts] for i in 1:3)
+@test all(plt[1][i + 3][:y] ≈ wp[i].errors for i in 1:3)
+@test all(startswith(plt[1][i + 3][:label], wp_names[i]) for i in 1:3)
+@test_throws BoundsError plt[1][7]
+@test_nowarn plot(wp, view = :dt_convergence, color = [:red :orange :green])
+@test_nowarn plot(wp, view = :dt_convergence, color = :lightblue, title = "Δt Convergence")

test/runtests.jl

@@ -8,3 +8,4 @@ using Test
 @time @testset "ODE Tableau Convergence Tests" begin include("ode_tableau_convergence_tests.jl") end ## Windows 32-bit fails on Butcher62 convergence test
 @time @testset "Analyticless Stochastic WP" begin include("analyticless_stochastic_wp.jl") end
 @time @testset "Stability Region Tests" begin include("stability_region_test.jl") end
+@time @testset "Plot Recipes" begin include("plotrecipes_tests.jl") end