diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index e321006..c6f57eb 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -24,6 +24,8 @@ jobs: exclude: - os: macOS-latest arch: x86 + - os: windows-latest + arch: x86 steps: - uses: actions/checkout@v6 - uses: julia-actions/setup-julia@v2 diff --git a/Project.toml b/Project.toml index 37ad2e5..8ab743c 100644 --- a/Project.toml +++ b/Project.toml @@ -5,6 +5,7 @@ version = "0.6.3" [deps] BranchAndPrune = "d3bc4f2e-91e6-11e9-365e-cd067da536ce" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" @@ -13,6 +14,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [compat] BranchAndPrune = "0.2.1" ForwardDiff = "0.10, 1" +InteractiveUtils = "1.10.0" IntervalArithmetic = "1.0.3" Reexport = "1" StaticArrays = "1" diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl index 02ca843..aafb7ef 100644 --- a/benchmark/benchmarks.jl +++ b/benchmark/benchmarks.jl @@ -20,7 +20,7 @@ S = SUITE["Smiley"] = BenchmarkGroup() for example in (SmileyExample22, SmileyExample52, SmileyExample54) #, SmileyExample55) s = S[example.title] = BenchmarkGroup() for contractor in (Newton, Krawczyk) - s[string(contractor)] = @benchmarkable roots($(example.f), $(example.region) ; contractor = $contractor, abstol = $tol) + s[string(contractor)] = @benchmarkable roots($(example.f), $(example.region) ; contractor = $contractor, abstol = $tol, infer_root_type = false) end end @@ -39,7 +39,7 @@ L = 5.0 X = SVector(interval(-L, (L+1)), interval(-L, (L+1))) for contractor in (Newton, Krawczyk) - S[string(contractor)] = @benchmarkable roots($(∇f), $X ; contractor = $contractor, abstol = 1e-5) + S[string(contractor)] = @benchmarkable roots($(∇f), $X ; contractor = $contractor, abstol = 1e-5, infer_root_type = false) end @@ -66,7 +66,7 @@ for (k, dr) in enumerate(dr_functions) if k != 8 # dr8 is excluded as it has too many roots for contractor in (Newton, Krawczyk) - s[string(contractor)] = @benchmarkable roots($dr, $X ; contractor = $contractor, abstol = $tol) + s[string(contractor)] = @benchmarkable roots($dr, $X ; contractor = $contractor, abstol = $tol, infer_root_type = false) end end end diff --git a/perf/linear_eq.jl b/perf/linear_eq.jl deleted file mode 100644 index ae25d51..0000000 --- a/perf/linear_eq.jl +++ /dev/null @@ -1,63 +0,0 @@ -using BenchmarkTools, Compat, DataFrames, IntervalRootFinding, IntervalArithmetic, StaticArrays - -function randVec(n::Int) - a = randn(n) - A = interval.(a) - mA = MVector{n}(A) - sA = SVector{n}(A) - return A, mA, sA -end - -function randMat(n::Int) - a = randn(n, n) - A = interval.(a) - mA = MMatrix{n, n}(A) - sA = SMatrix{n, n}(A) - return A, mA, sA -end - -function benchmark(max=10) - df = DataFrame() - df[:Method] = ["Array", "MArray", "SArray", "Contractor", "ContractorMArray", "ContractorSArray"] - for n in 1:max - A, mA, sA = randMat(n) - b, mb, sb = randVec(n) - t1 = @belapsed gauss_seidel_interval($A, $b) - t2 = @belapsed gauss_seidel_interval($mA, $mb) - t3 = @belapsed gauss_seidel_interval($sA, $sb) - t4 = @belapsed gauss_seidel_contractor($A, $b) - t5 = @belapsed gauss_seidel_contractor($mA, $mb) - t6 = @belapsed gauss_seidel_contractor($sA, $sb) - df[Symbol("$n")] = [t1, t2, t3, t4, t5, t6] - end - a = [] - for i in 1:max - push!(a, Symbol("$i")) - end - df1 = stack(df, a) - dfnew = unstack(df1, :variable, :Method, :value) - dfnew = rename(dfnew, :variable => :n) - println(dfnew) - dfnew -end - -function benchmark_elimination(max=10) - df = DataFrame() - df[:Method] = ["Gauss Elimination", "Base.\\"] - for n in 1:max - A, mA, sA = randMat(n) - b, mb, sb = randVec(n) - t1 = @belapsed gauss_elimination_interval($A, $b) - t2 = @belapsed gauss_elimination_interval1($A, $b) - df[Symbol("$n")] = [t1, t2] - end - a = [] - for i in 1:max - push!(a, Symbol("$i")) - end - df1 = stack(df, a) - dfnew = unstack(df1, :variable, :Method, :value) - dfnew = rename(dfnew, :variable => :n) - println(dfnew) - dfnew -end diff --git a/perf/linear_eq_results.txt b/perf/linear_eq_results.txt deleted file mode 100644 index ff08cef..0000000 --- a/perf/linear_eq_results.txt +++ /dev/null @@ -1,149 +0,0 @@ -For n = 1 - 11.382 μs (420 allocations: 16.83 KiB) -Array: IntervalArithmetic.Interval{Float64}[[-0.858728, -0.858727]] - 10.980 μs (421 allocations: 16.77 KiB) -MArray: IntervalArithmetic.Interval{Float64}[[-0.858728, -0.858727]] - -For n = 2 - 27.322 μs (822 allocations: 33.84 KiB) -Array: IntervalArithmetic.Interval{Float64}[[0.946446, 0.946447], [0.40483, 0.404831]] - 26.037 μs (821 allocations: 32.41 KiB) -MArray: IntervalArithmetic.Interval{Float64}[[0.946446, 0.946447], [0.40483, 0.404831]] - -For n = 3 - 50.094 μs (1222 allocations: 50.20 KiB) -Array: IntervalArithmetic.Interval{Float64}[[1.21978, 1.21979], [-2.29245, -2.29244], [-5.4913, -5.49129]] - 47.146 μs (1221 allocations: 48.06 KiB) -MArray: IntervalArithmetic.Interval{Float64}[[1.21978, 1.21979], [-2.29245, -2.29244], [-5.4913, -5.49129]] - -For n = 4 - 78.244 μs (1626 allocations: 66.92 KiB) -Array: IntervalArithmetic.Interval{Float64}[[-1.23305, -1.23304], [-0.267412, -0.267411], [-27.3277, -27.3276], [-8.3822, -8.38219]] - 74.046 μs (1621 allocations: 63.70 KiB) -MArray: IntervalArithmetic.Interval{Float64}[[-1.23305, -1.23304], [-0.267412, -0.267411], [-27.3277, -27.3276], [-8.3822, -8.38219]] - -For n = 5 - 109.184 μs (2026 allocations: 83.63 KiB) -Array: IntervalArithmetic.Interval{Float64}[[0.670434, 0.670435], [0.876143, 0.876144], [0.0960248, 0.0960249], [-0.109509, -0.109508], [0.375848, 0.375849]] - 112.329 μs (2032 allocations: 83.41 KiB) -MArray: IntervalArithmetic.Interval{Float64}[[0.670434, 0.670435], [0.876143, 0.876144], [0.0960248, 0.0960249], [-0.109509, -0.109508], [0.375848, 0.375849]] - -For n = 6 - 148.001 μs (2426 allocations: 100.09 KiB) -Array: IntervalArithmetic.Interval{Float64}[[-2.40426, -2.40425], [-2.04742, -2.04741], [-1.39128, -1.39127], [-1.28434, -1.28433], [3.76269, 3.7627], [-1.28934, -1.28933]] - 148.599 μs (2432 allocations: 99.73 KiB) -MArray: IntervalArithmetic.Interval{Float64}[[-2.40426, -2.40425], [-2.04742, -2.04741], [-1.39128, -1.39127], [-1.28434, -1.28433], [3.76269, 3.7627], [-1.28934, -1.28933]] - -For n = 7 - 191.907 μs (2826 allocations: 116.92 KiB) -Array: IntervalArithmetic.Interval{Float64}[[0.706613, 0.706614], [2.97926, 2.97927], [-2.55176, -2.55175], [-1.51898, -1.51897], [-0.987214, -0.987213], [-2.07295, -2.07294], [2.29913, 2.29914]] - 194.024 μs (2832 allocations: 116.34 KiB) -MArray: IntervalArithmetic.Interval{Float64}[[0.706613, 0.706614], [2.97926, 2.97927], [-2.55176, -2.55175], [-1.51898, -1.51897], [-0.987214, -0.987213], [-2.07295, -2.07294], [2.29913, 2.29914]] - -For n = 8 - 246.387 μs (3226 allocations: 133.58 KiB) -Array: IntervalArithmetic.Interval{Float64}[[-0.149055, -0.149054], [7.1198, 7.11981], [-8.30751, -8.3075], [-2.37305, -2.37304], [0.0150316, 0.0150317], [3.35756, 3.35757], [3.9637, 3.96371], [-4.7502, -4.75019]] - 248.388 μs (3232 allocations: 132.72 KiB) -MArray: IntervalArithmetic.Interval{Float64}[[-0.149055, -0.149054], [7.1198, 7.11981], [-8.30751, -8.3075], [-2.37305, -2.37304], [0.0150316, 0.0150317], [3.35756, 3.35757],[3.9637, 3.96371], [-4.7502, -4.75019]] - -For n = 9 - 313.552 μs (3626 allocations: 150.41 KiB) -Array: IntervalArithmetic.Interval{Float64}[[-0.319338, -0.319337], [1.05902, 1.05903], [-1.87228, -1.87227], [-0.878659, -0.878658], [1.2596, 1.25961], [-2.90817, -2.90816],[-2.77684, -2.77683], [-3.76999, -3.76998], [-0.355478, -0.355477]] - 309.842 μs (3632 allocations: 149.23 KiB) -MArray: IntervalArithmetic.Interval{Float64}[[-0.319338, -0.319337], [1.05902, 1.05903], [-1.87228, -1.87227], [-0.878659, -0.878658], [1.2596, 1.25961], [-2.90817, -2.90816], [-2.77684, -2.77683], [-3.76999, -3.76998], [-0.355478, -0.355477]] - -For n = 10 - 383.613 μs (4026 allocations: 167.34 KiB) -Array: IntervalArithmetic.Interval{Float64}[[2.98301, 2.98302], [0.995662, 0.995663], [0.150873, 0.150874], [1.03224, 1.03225], [-0.887891, -0.88789], [0.547691, 0.547692], [-0.530325, -0.530324], [5.39277, 5.39278], [0.230642, 0.230643], [-1.87602, -1.87601]] - 393.352 μs (4032 allocations: 165.84 KiB) -MArray: IntervalArithmetic.Interval{Float64}[[2.98301, 2.98302], [0.995662, 0.995663], [0.150873, 0.150874], [1.03224, 1.03225], [-0.887891, -0.88789], [0.547691, 0.547692], [-0.530325, -0.530324], [5.39277, 5.39278], [0.230642, 0.230643], [-1.87602, -1.87601]] - -For n = 11 - 460.411 μs (4426 allocations: 184.31 KiB) -Array: IntervalArithmetic.Interval{Float64}[[-1.49959, -1.49958], [0.280678, 0.280679], [-1.53757, -1.53756], [1.2517, 1.25171], [-0.605238, -0.605237], [0.825546, 0.825547],[0.755109, 0.75511], [0.18521, 0.185211], [-1.82804, -1.82803], [0.384323, 0.384324], [0.241048, 0.241049]] - 458.773 μs (4432 allocations: 182.59 KiB) -MArray: IntervalArithmetic.Interval{Float64}[[-1.49959, -1.49958], [0.280678, 0.280679], [-1.53757, -1.53756], [1.2517, 1.25171], [-0.605238, -0.605237], [0.825546, 0.825547], [0.755109, 0.75511], [0.18521, 0.185211], [-1.82804, -1.82803], [0.384323, 0.384324], [0.241048, 0.241049]] - -For n = 12 - 550.373 μs (4826 allocations: 201.45 KiB) -Array: IntervalArithmetic.Interval{Float64}[[-2.79797, -2.79796], [-1.73355, -1.73354], [-0.606508, -0.606507], [-2.44883, -2.44882], [1.71346, 1.71347], [-2.72265, -2.72264], [-1.55253, -1.55252], [0.289974, 0.289975], [0.252046, 0.252047], [0.762136, 0.762137], [1.5864, 1.58641], [-1.50617, -1.50616]] - 560.079 μs (4832 allocations: 199.20 KiB) -MArray: IntervalArithmetic.Interval{Float64}[[-2.79797, -2.79796], [-1.73355, -1.73354], [-0.606508, -0.606507], [-2.44883, -2.44882], [1.71346, 1.71347], [-2.72265, -2.72264], [-1.55253, -1.55252], [0.289974, 0.289975], [0.252046, 0.252047], [0.762136, 0.762137], [1.5864, 1.58641], [-1.50617, -1.50616]] - -For n = 13 - 656.259 μs (5226 allocations: 218.75 KiB) -Array: IntervalArithmetic.Interval{Float64}[[2.82409, 2.8241], [-1.48454, -1.48453], [-0.31499, -0.314989], [3.10565, 3.10566], [1.40989, 1.4099], [1.96286, 1.96287], [1.07219, 1.0722], [-4.11245, -4.11244], [-1.54954, -1.54953], [1.06494, 1.06495], [-3.6406, -3.64059], [1.50881, 1.50882], [1.19286, 1.19287]] - 676.472 μs (5232 allocations: 216.09 KiB) -MArray: IntervalArithmetic.Interval{Float64}[[2.82409, 2.8241], [-1.48454, -1.48453], [-0.31499, -0.314989], [3.10565, 3.10566], [1.40989, 1.4099], [1.96286, 1.96287], [1.07219, 1.0722], [-4.11245, -4.11244], [-1.54954, -1.54953], [1.06494, 1.06495], [-3.6406, -3.64059], [1.50881, 1.50882], [1.19286, 1.19287]] - -For n = 14 - 780.545 μs (5626 allocations: 236.19 KiB) -Array: IntervalArithmetic.Interval{Float64}[[0.78063, 0.780631], [0.0801395, 0.0801396], [-0.619735, -0.619734], [-0.408535, -0.408534], [-0.35918, -0.359179], [-0.279269, -0.279268], [-0.381064, -0.381063], [-0.583649, -0.583648], [0.15355, 0.153551], [0.303571, 0.303572], [-0.263495, -0.263494], [0.745922, 0.745923], [-1.25136, -1.25135], [-0.306904, -0.306903]] - 792.250 μs (5632 allocations: 233.17 KiB) -MArray: IntervalArithmetic.Interval{Float64}[[0.78063, 0.780631], [0.0801395, 0.0801396], [-0.619735, -0.619734], [-0.408535, -0.408534], [-0.35918, -0.359179], [-0.279269, -0.279268], [-0.381064, -0.381063], [-0.583649, -0.583648], [0.15355, 0.153551], [0.303571, 0.303572], [-0.263495, -0.263494], [0.745922, 0.745923], [-1.25136, -1.25135], [-0.306904, -0.306903]] - -For n = 15 - 896.447 μs (6026 allocations: 253.50 KiB) -Array: IntervalArithmetic.Interval{Float64}[[-1.31651, -1.3165], [1.50793, 1.50794], [-1.66271, -1.6627], [0.42862, 0.428621], [0.764462, 0.764463], [0.376219, 0.37622], [0.0408481, 0.0408482], [0.078441, 0.0784411], [-1.53466, -1.53465], [0.498318, 0.498319], [-0.257995, -0.257994], [-0.604852, -0.604851], [-0.386613, -0.386612], [-0.0438574, -0.0438573], [0.440631, 0.440632]] - 950.789 μs (6032 allocations: 250.02 KiB) -MArray: IntervalArithmetic.Interval{Float64}[[-1.31651, -1.3165], [1.50793, 1.50794], [-1.66271, -1.6627], [0.42862, 0.428621], [0.764462, 0.764463], [0.376219, 0.37622], [0.0408481, 0.0408482], [0.078441, 0.0784411], [-1.53466, -1.53465], [0.498318, 0.498319], [-0.257995, -0.257994], [-0.604852, -0.604851], [-0.386613, -0.386612], [-0.0438574, -0.0438573], [0.440631, 0.440632]] - -For n = 16 - 1.076 ms (6426 allocations: 270.61 KiB) -Array: IntervalArithmetic.Interval{Float64}[[-0.0621213, -0.0621212], [-0.436339, -0.436338], [0.680463, 0.680464], [0.55125, 0.551251], [0.979387, 0.979388], [-0.904311, -0.90431], [-1.34208, -1.34207], [2.02843, 2.02844], [1.99408, 1.99409], [-1.1726, -1.17259], [1.25116, 1.25117], [-0.313683, -0.313682], [1.09734, 1.09735], [0.232912, 0.232913], [-1.68377, -1.68376], [-0.425297, -0.425296]] - 1.116 ms (6432 allocations: 266.64 KiB) -MArray: IntervalArithmetic.Interval{Float64}[[-0.0621213, -0.0621212], [-0.436339, -0.436338], [0.680463, 0.680464], [0.55125, 0.551251], [0.979387, 0.979388], [-0.904311, -0.90431], [-1.34208, -1.34207], [2.02843, 2.02844], [1.99408, 1.99409], [-1.1726, -1.17259], [1.25116, 1.25117], [-0.313683, -0.313682], [1.09734, 1.09735], [0.232912, 0.232913], [-1.68377, -1.68376], [-0.425297, -0.425296]] - -For n = 17 - 1.205 ms (6826 allocations: 288.27 KiB) -Array: IntervalArithmetic.Interval{Float64}[[0.543059, 0.54306], [-1.29556, -1.29555], [-1.7007, -1.70069], [1.43802, 1.43803], [0.179184, 0.179185], [1.69012, 1.69013], [1.34535, 1.34536], [1.49156, 1.49157], [1.15319, 1.1532], [-0.629132, -0.629131], [-1.15364, -1.15363], [0.627152, 0.627153], [-1.82, -1.81999], [-2.69038, -2.69037], [2.13837, 2.13838], [-0.653879, -0.653878], [-0.585587, -0.585586]] - 1.284 ms (6832 allocations: 283.75 KiB) -MArray: IntervalArithmetic.Interval{Float64}[[0.543059, 0.54306], [-1.29556, -1.29555], [-1.7007, -1.70069], [1.43802, 1.43803], [0.179184, 0.179185], [1.69012, 1.69013], [1.34535, 1.34536], [1.49156, 1.49157], [1.15319, 1.1532], [-0.629132, -0.629131], [-1.15364, -1.15363], [0.627152, 0.627153], [-1.82, -1.81999], [-2.69038, -2.69037], [2.13837, 2.13838], [-0.653879, -0.653878], [-0.585587, -0.585586]] - -For n = 18 - 1.541 ms (7226 allocations: 305.64 KiB) -Array: IntervalArithmetic.Interval{Float64}[[-4.4038, -4.40379], [35.4039, 35.404], [23.2018, 23.2019], [29.5486, 29.5487], [-4.1716, -4.17159], [-10.229, -10.2289], [-4.36794, -4.36793], [-3.06734, -3.06733], [28.4815, 28.4816], [20.6252, 20.6253], [8.52634, 8.52635], [19.3466, 19.3467], [-17.4713, -17.4712], [-40.2325, -40.2324], [8.11921, 8.11922], [-23.8136, -23.8135], [17.8947, 17.8948], [55.5122, 55.5123]] - 1.666 ms (7232 allocations: 300.63 KiB) -MArray: IntervalArithmetic.Interval{Float64}[[-4.4038, -4.40379], [35.4039, 35.404], [23.2018, 23.2019], [29.5486, 29.5487], [-4.1716, -4.17159], [-10.229, -10.2289], [-4.36794, -4.36793], [-3.06734, -3.06733], [28.4815, 28.4816], [20.6252, 20.6253], [8.52634, 8.52635], [19.3466, 19.3467], [-17.4713, -17.4712], [-40.2325, -40.2324], [8.11921, 8.11922], [-23.8136, -23.8135], [17.8947, 17.8948], [55.5122, 55.5123]] - -For n = 19 - 1.769 ms (7626 allocations: 323.36 KiB) -Array: IntervalArithmetic.Interval{Float64}[[-0.16665, -0.166649], [-0.836062, -0.836061], [0.191566, 0.191567], [0.359444, 0.359445], [0.489007, 0.489008], [-0.125, -0.124999], [-0.397535, -0.397534], [-0.45875, -0.458749], [0.735181, 0.735182], [1.68832, 1.68833], [0.555231, 0.555232], [1.00734, 1.00735], [-0.433021, -0.43302], [0.742541, 0.742542], [-1.01507, -1.01506], [1.42938, 1.42939], [-0.306669, -0.306668], [-0.657369, -0.657368], [-0.139753, -0.139752]] - 1.930 ms (7632 allocations: 317.73 KiB) -MArray: IntervalArithmetic.Interval{Float64}[[-0.16665, -0.166649], [-0.836062, -0.836061], [0.191566, 0.191567], [0.359444, 0.359445], [0.489007, 0.489008], [-0.125, -0.124999], [-0.397535, -0.397534], [-0.45875, -0.458749], [0.735181, 0.735182], [1.68832, 1.68833], [0.555231, 0.555232], [1.00734, 1.00735], [-0.433021, -0.43302], [0.742541, 0.742542], [-1.01507, -1.01506], [1.42938, 1.42939], [-0.306669, -0.306668], [-0.657369, -0.657368], [-0.139753, -0.139752]] - -For n = 20 - 1.964 ms (8026 allocations: 340.89 KiB) -Array: IntervalArithmetic.Interval{Float64}[[1.35974, 1.35975], [0.0315064, 0.0315065], [1.62038, 1.62039], [1.49915, 1.49916], [0.846944, 0.846945], [1.45553, 1.45554], [4.53343, 4.53344], [1.90342, 1.90343], [-4.99214, -4.99213], [2.07785, 2.07786], [-0.444216, -0.444215], [0.736985, 0.736986], [1.07722, 1.07723], [2.75322, 2.75323], [1.59534, 1.59535], [2.78334, 2.78335], [-5.37279, -5.37278], [-0.843305, -0.843304], [2.58871, 2.58872], [0.805166, 0.805167]] - 2.225 ms (8032 allocations: 334.67 KiB) -MArray: IntervalArithmetic.Interval{Float64}[[1.35974, 1.35975], [0.0315064, 0.0315065], [1.62038, 1.62039], [1.49915, 1.49916], [0.846944, 0.846945], [1.45553, 1.45554], [4.53343, 4.53344], [1.90342, 1.90343], [-4.99214, -4.99213], [2.07785, 2.07786], [-0.444216, -0.444215], [0.736985, 0.736986], [1.07722, 1.07723], [2.75322, 2.75323], [1.59534, 1.59535], [2.78334, 2.78335], [-5.37279, -5.37278], [-0.843305, -0.843304], [2.58871, 2.58872], [0.805166, 0.805167]] - -For n = 21 - 2.279 ms (8426 allocations: 358.86 KiB) -Array: IntervalArithmetic.Interval{Float64}[[0.531469, 0.53147], [-1.55322, -1.55321], [0.610774, 0.610775], [-0.972668, -0.972667], [-0.995764, -0.995763], [1.53137, 1.53138], [0.226348, 0.226349], [-0.0360293, -0.0360292], [0.300785, 0.300786], [1.16974, 1.16975], [0.224403, 0.224404], [-1.59163, -1.59162], [-0.33666, -0.336659], [-0.365669, -0.365668], [0.671668, 0.671669], [0.166757, 0.166758], [-0.0187087, -0.0187086], [0.000503947, 0.000503948], [0.747613, 0.747614], [-0.436936, -0.436935], [-1.07276, -1.07275]] - 2.523 ms (8432 allocations: 351.97 KiB) -MArray: IntervalArithmetic.Interval{Float64}[[0.531469, 0.53147], [-1.55322, -1.55321], [0.610774, 0.610775], [-0.972668, -0.972667], [-0.995764, -0.995763], [1.53137, 1.53138], [0.226348, 0.226349], [-0.0360293, -0.0360292], [0.300785, 0.300786], [1.16974, 1.16975], [0.224403, 0.224404], [-1.59163, -1.59162], [-0.33666, -0.336659], [-0.365669, -0.365668], [0.671668, 0.671669], [0.166757, 0.166758], [-0.0187087, -0.0187086], [0.000503947, 0.000503948], [0.747613, 0.747614], [-0.436936, -0.436935], [-1.07276, -1.07275]] - -For n = 22 - 2.635 ms (8826 allocations: 376.55 KiB) -Array: IntervalArithmetic.Interval{Float64}[[-0.0576031, -0.057603], [-0.19045, -0.190449], [-0.145126, -0.145125], [0.0613842, 0.0613843], [-0.0716046, -0.0716045], [-0.020995, -0.0209949], [-0.399329, -0.399328], [-0.291664, -0.291663], [-0.131188, -0.131187], [0.287153, 0.287154], [-0.476477, -0.476476], [0.937626, 0.937627], [-0.38509, -0.385089], [0.0836971, 0.0836972], [0.176637, 0.176638], [0.887916, 0.887917], [0.0787311, 0.0787312], [-0.404081, -0.40408], [-0.319168, -0.319167], [-0.981298, -0.981297], [-0.0197936, -0.0197935], [0.589842, 0.589843]] - 2.997 ms (8832 allocations: 369.03 KiB) -MArray: IntervalArithmetic.Interval{Float64}[[-0.0576031, -0.057603], [-0.19045, -0.190449], [-0.145126, -0.145125], [0.0613842, 0.0613843], [-0.0716046, -0.0716045], [-0.020995, -0.0209949], [-0.399329, -0.399328], [-0.291664, -0.291663], [-0.131188, -0.131187], [0.287153, 0.287154], [-0.476477, -0.476476], [0.937626, 0.937627], [-0.38509, -0.385089], [0.0836971, 0.0836972], [0.176637, 0.176638], [0.887916, 0.887917], [0.0787311, 0.0787312], [-0.404081, -0.40408], [-0.319168, -0.319167], [-0.981298, -0.981297], [-0.0197936, -0.0197935], [0.589842, 0.589843]] - -For n = 23 - 2.810 ms (9226 allocations: 394.70 KiB) -Array: IntervalArithmetic.Interval{Float64}[[-1.89769, -1.89768], [-2.13486, -2.13485], [-0.0239536, -0.0239535], [2.26537, 2.26538], [0.142881, 0.142882], [-1.55981, -1.5598], [-1.5167, -1.51669], [-0.511803, -0.511802], [0.571971, 0.571972], [-0.00964798, -0.00964797], [-0.858116, -0.858115], [-0.239622, -0.239621], [-0.891741, -0.89174], [-2.30866, -2.30865], [-1.76962, -1.76961], [-2.61382, -2.61381], [-0.457754, -0.457753], [-1.03219, -1.03218], [-1.18658, -1.18657], [0.0264487, 0.0264488], [0.0913128, 0.0913129],[-0.0189498, -0.0189497], [-0.922417, -0.922416]] - 3.222 ms (9232 allocations: 386.45 KiB) -MArray: IntervalArithmetic.Interval{Float64}[[-1.89769, -1.89768], [-2.13486, -2.13485], [-0.0239536, -0.0239535], [2.26537, 2.26538], [0.142881, 0.142882], [-1.55981, -1.5598], [-1.5167, -1.51669], [-0.511803, -0.511802], [0.571971, 0.571972], [-0.00964798, -0.00964797], [-0.858116, -0.858115], [-0.239622, -0.239621], [-0.891741, -0.89174], [-2.30866, -2.30865], [-1.76962, -1.76961], [-2.61382, -2.61381], [-0.457754, -0.457753], [-1.03219, -1.03218], [-1.18658, -1.18657], [0.0264487, 0.0264488], [0.0913128, 0.0913129], [-0.0189498, -0.0189497], [-0.922417, -0.922416]] - -For n = 24 - 3.219 ms (9626 allocations: 412.64 KiB) -Array: IntervalArithmetic.Interval{Float64}[[-0.124356, -0.124355], [-0.00958074, -0.00958073], [0.539476, 0.539477], [-0.388642, -0.388641], [0.625504, 0.625505], [1.09303, 1.09304], [-1.11944, -1.11943], [0.638843, 0.638844], [-0.522925, -0.522924], [-1.37454, -1.37453], [0.642279, 0.64228], [0.764112, 0.764113], [-0.568954, -0.568953], [-0.558264, -0.558263], [0.58252, 0.582521], [0.735529, 0.73553], [-0.365248, -0.365247], [-1.32431, -1.3243], [-0.150569, -0.150568], [-0.399501, -0.3995], [1.37909, 1.3791], [0.00529648, 0.00529649], [0.21186, 0.211861], [0.6825, 0.682501]] - 3.707 ms (9632 allocations: 403.56 KiB) -MArray: IntervalArithmetic.Interval{Float64}[[-0.124356, -0.124355], [-0.00958074, -0.00958073], [0.539476, 0.539477], [-0.388642, -0.388641], [0.625504, 0.625505], [1.09303,1.09304], [-1.11944, -1.11943], [0.638843, 0.638844], [-0.522925, -0.522924], [-1.37454, -1.37453], [0.642279, 0.64228], [0.764112, 0.764113], [-0.568954, -0.568953], [-0.558264, -0.558263], [0.58252, 0.582521], [0.735529, 0.73553], [-0.365248, -0.365247], [-1.32431, -1.3243], [-0.150569, -0.150568], [-0.399501, -0.3995], [1.37909, 1.3791], [0.00529648, 0.00529649], [0.21186, 0.211861], [0.6825, 0.682501]] - -For n = 25 - 3.587 ms (10028 allocations: 431.00 KiB) -Array: IntervalArithmetic.Interval{Float64}[[-1.74061, -1.7406], [0.412884, 0.412885], [0.623495, 0.623496], [-1.42243, -1.42242], [-0.774224, -0.774223], [-0.0379006, -0.0379005], [0.735943, 0.735944], [0.0682942, 0.0682943], [0.992878, 0.992879], [-0.492704, -0.492703], [-1.38519, -1.38518], [0.387523, 0.387524], [-3.107, -3.10699], [0.0425826, 0.0425827], [-1.00909, -1.00908], [-0.164128, -0.164127], [0.0325966, 0.0325967], [-0.52808, -0.528079], [0.499164, 0.499165], [-0.00699219, -0.00699218], [0.034, 0.0340001], [0.878137, 0.878138], [0.396617, 0.396618], [-1.18615, -1.18614], [0.043228, 0.0432281]] - 4.199 ms (10032 allocations: 421.02 KiB) -MArray: IntervalArithmetic.Interval{Float64}[[-1.74061, -1.7406], [0.412884, 0.412885], [0.623495, 0.623496], [-1.42243, -1.42242], [-0.774224, -0.774223], [-0.0379006, -0.0379005], [0.735943, 0.735944], [0.0682942, 0.0682943], [0.992878, 0.992879], [-0.492704, -0.492703], [-1.38519, -1.38518], [0.387523, 0.387524], [-3.107, -3.10699], [0.0425826,0.0425827], [-1.00909, -1.00908], [-0.164128, -0.164127], [0.0325966, 0.0325967], [-0.52808, -0.528079], [0.499164, 0.499165], [-0.00699219, -0.00699218], [0.034, 0.0340001],[0.878137, 0.878138], [0.396617, 0.396618], [-1.18615, -1.18614], [0.043228, 0.0432281]] diff --git a/perf/linear_eq_tabular.txt b/perf/linear_eq_tabular.txt deleted file mode 100644 index 6902838..0000000 --- a/perf/linear_eq_tabular.txt +++ /dev/null @@ -1,49 +0,0 @@ -Gauss Seidel - -│ Row │ variable │ Array │ Contractor │ MArray │ SArray1 │ SArray2 │ -├─────┼──────────┼─────────────┼─────────────┼─────────────┼─────────────┼─────────────┤ -│ 1 │ n = 1 │ 1.1426e-5 │ 2.4173e-5 │ 1.5803e-5 │ 1.014e-5 │ 9.512e-6 │ -│ 2 │ n = 2 │ 2.4413e-5 │ 4.0504e-5 │ 3.9829e-5 │ 2.4083e-5 │ 2.4083e-5 │ -│ 3 │ n = 3 │ 4.5194e-5 │ 6.4343e-5 │ 9.1236e-5 │ 4.5533e-5 │ 5.1221e-5 │ -│ 4 │ n = 4 │ 7.2639e-5 │ 9.3671e-5 │ 0.00016058 │ 7.0972e-5 │ 6.9937e-5 │ -│ 5 │ n = 5 │ 0.000103281 │ 0.000127489 │ 0.000262814 │ 0.000104823 │ 0.000102729 │ -│ 6 │ n = 6 │ 0.000141315 │ 0.000169992 │ 0.000416266 │ 0.000144021 │ 0.000139086 │ -│ 7 │ n = 7 │ 0.000195001 │ 0.000217059 │ 0.000680848 │ 0.000198164 │ 0.000191145 │ -│ 8 │ n = 8 │ 0.000247193 │ 0.00027535 │ 0.00122559 │ 0.000251089 │ 0.000241972 │ -│ 9 │ n = 9 │ 0.000306262 │ 0.000338028 │ 0.0020223 │ 0.000310961 │ 0.000299153 │ -│ 10 │ n = 10 │ 0.00037781 │ 0.000414073 │ 0.00294134 │ 0.000386235 │ 0.000367335 │ -│ 11 │ n = 11 │ 0.000465017 │ 0.000489114 │ 0.00585246 │ 0.000473425 │ 0.000447176 │ -│ 12 │ n = 12 │ 0.00055403 │ 0.000587152 │ 0.0110894 │ 0.000573395 │ 0.00054274 │ -│ 13 │ n = 13 │ 0.000659994 │ 0.000692402 │ 0.0160662 │ 0.000686098 │ 0.000643391 │ -│ 14 │ n = 14 │ 0.000777718 │ 0.000793626 │ 0.0184595 │ 0.000800443 │ 0.000767576 │ -│ 15 │ n = 15 │ 0.000910174 │ 0.000952207 │ 0.0239422 │ 0.000972855 │ 0.000899402 │ -│ 16 │ n = 16 │ 0.00107785 │ 0.00111443 │ 0.0303295 │ 0.00113467 │ 0.0010867 │ -│ 17 │ n = 17 │ 0.00122326 │ 0.00125759 │ 0.0398127 │ 0.00134961 │ 0.00125675 │ -│ 18 │ n = 18 │ 0.00167176 │ 0.00159463 │ 0.043795 │ 0.00181409 │ 0.00162584 │ -│ 19 │ n = 19 │ 0.0018632 │ 0.00185617 │ 0.0549107 │ 0.00208529 │ 0.00196587 │ -│ 20 │ n = 20 │ 0.0020604 │ 0.00210595 │ 0.0635598 │ 0.00232704 │ 0.00212374 │ - -10×7 DataFrames.DataFrame -│ Row │ n │ Array │ Contractor │ ContractorMArray │ ContractorSArray │ MArray │ SArray │ -├─────┼────┼─────────────┼─────────────┼──────────────────┼──────────────────┼─────────────┼─────────────┤ -│ 1 │ 1 │ 1.0188e-5 │ 2.5276e-5 │ 0.000781355 │ 0.000764627 │ 1.6959e-5 │ 1.1216e-5 │ -│ 2 │ 10 │ 0.000377853 │ 0.000407817 │ 0.00122723 │ 0.00123259 │ 0.0029247 │ 0.000393896 │ -│ 3 │ 2 │ 2.609e-5 │ 4.3597e-5 │ 0.000802467 │ 0.000796962 │ 4.1039e-5 │ 2.6654e-5 │ -│ 4 │ 3 │ 4.7362e-5 │ 6.317e-5 │ 0.000832029 │ 0.000821289 │ 9.1695e-5 │ 4.7007e-5 │ -│ 5 │ 4 │ 7.0378e-5 │ 9.1502e-5 │ 0.00085429 │ 0.000844048 │ 0.000163604 │ 7.1166e-5 │ -│ 6 │ 5 │ 0.000104813 │ 0.000129112 │ 0.000898626 │ 0.000889969 │ 0.000268756 │ 0.000112476 │ -│ 7 │ 6 │ 0.000142316 │ 0.000168891 │ 0.000941422 │ 0.000922727 │ 0.000430442 │ 0.000146906 │ -│ 8 │ 7 │ 0.000192073 │ 0.000218861 │ 0.00108794 │ 0.000997821 │ 0.00073304 │ 0.000198962 │ -│ 9 │ 8 │ 0.000242832 │ 0.000273341 │ 0.00107159 │ 0.00106769 │ 0.00126407 │ 0.000250956 │ -│ 10 │ 9 │ 0.000308256 │ 0.000333722 │ 0.00116203 │ 0.00116428 │ 0.00196451 │ 0.000314499 │ - -Gauss Elimination - -5×3 DataFrames.DataFrame -│ Row │ n │ Base.\\ │ Gauss Elimination │ -├─────┼───┼────────────┼───────────────────┤ -│ 1 │ 1 │ 1.84e-6 │ 8.127e-6 │ -│ 2 │ 2 │ 3.1615e-6 │ 1.0029e-5 │ -│ 3 │ 3 │ 4.53986e-6 │ 1.1211e-5 │ -│ 4 │ 4 │ 6.8655e-6 │ 1.4342e-5 │ -│ 5 │ 5 │ 1.0773e-5 │ 1.7725e-5 │ diff --git a/perf/slopes.jl b/perf/slopes.jl deleted file mode 100644 index ea3e949..0000000 --- a/perf/slopes.jl +++ /dev/null @@ -1,66 +0,0 @@ -using BenchmarkTools, Compat, DataFrames, IntervalRootFinding, IntervalArithmetic, StaticArrays - -function benchmark_results() - f = [] # Example functions from Dietmar Ratz - An Optimized Interval Slope Arithmetic and its Application - push!(f, x->((x + sin(x)) * exp(-x^2))) - push!(f, x->(x^4 - 10x^3 + 35x^2 - 50x + 24)) - push!(f, x->((log(x + 1.25) - 0.84x) ^ 2)) - push!(f, x->(0.02x^2 - 0.03exp(-(20(x - 0.875))^2))) - push!(f, x->(exp(x^2))) - push!(f, x->(x^4 - 12x^3 + 47x^2 - 60x - 20exp(-x))) - push!(f, x->(x^6 - 15x^4 + 27x^2 + 250)) - push!(f, x->(atan(cos(tan(x))))) - push!(f, x->(asin(cos(acos(sin(x)))))) - - s = interval(0.75, 1.75) - df = DataFrame() - - df[:Method] = ["Automatic Differentiation", "Slope Expansion"] - for n in 1:length(f) - - t1 = ForwardDiff.derivative(f[n], s) - t2 = slope(f[n], s, mid(s)) - df[Symbol("f" * "$n")] = [t1, t2] - end - a = [] - for i in 1:length(f) - push!(a, Symbol("f" * "$i")) - end - df1 = stack(df, a) - dfnew = unstack(df1, :variable, :Method, :value) - dfnew = rename(dfnew, :variable => :Function) - println(dfnew) - dfnew -end - -function benchmark_time() - f = [] - push!(f, x->((x + sin(x)) * exp(-x^2))) - push!(f, x->(x^4 - 10x^3 + 35x^2 - 50x + 24)) - push!(f, x->((log(x + 1.25) - 0.84x) ^ 2)) - push!(f, x->(0.02x^2 - 0.03exp(-(20(x - 0.875))^2))) - push!(f, x->(exp(x^2))) - push!(f, x->(x^4 - 12x^3 + 47x^2 - 60x - 20exp(-x))) - push!(f, x->(x^6 - 15x^4 + 27x^2 + 250)) - push!(f, x->(atan(cos(tan(x))))) - push!(f, x->(asin(cos(acos(sin(x)))))) - - s = interval(0.75, 1.75) - df = DataFrame() - df[:Method] = ["Automatic Differentiation", "Slope Expansion"] - for n in 1:length(f) - - t1 = @belapsed ForwardDiff.derivative($f[$n], $s) - t2 = @belapsed slope($f[$n], $s, mid($s)) - df[Symbol("f" * "$n")] = [t1, t2] - end - a = [] - for i in 1:length(f) - push!(a, Symbol("f" * "$i")) - end - df1 = stack(df, a) - dfnew = unstack(df1, :variable, :Method, :value) - dfnew = rename(dfnew, :variable => :Function) - println(dfnew) - dfnew -end diff --git a/perf/slopes_tabular.txt b/perf/slopes_tabular.txt deleted file mode 100644 index 2d74424..0000000 --- a/perf/slopes_tabular.txt +++ /dev/null @@ -1,37 +0,0 @@ -Results - -│ Row │ Function │ Automatic Differentiation │ Slope Expansion │ -├─────┼──────────┼───────────────────────────┼───────────────────────┤ -│ 1 │ f1 │ [-5.44573, 0.886249] │ [-2.79917, 0.0521429] │ -│ 2 │ f2 │ [-87.6875, 77.0625] │ [-43.875, 38.25] │ -│ 3 │ f3 │ [-0.474861, 0.787211] │ [-0.159199, 0.432842] │ -│ 4 │ f4 │ [-2.97001, 21.0701] │ [0.04, 0.326667] │ -│ 5 │ f5 │ [2.63258, 74.8333] │ [6.03135, 33.2205] │ -│ 6 │ f6 │ [-94.5871, 115.135] │ [-38.9908, 65.5595] │ -│ 7 │ f7 │ [-279.639, 167.667] │ [-146.852, 67.0665] │ -│ 8 │ f8 │ [-∞, ∞] │ [1, 2] │ -│ 9 │ f9 │ [-∞, ∞] │ [1.3667, ∞] │ - -Time - -│ Row │ Function │ Automatic Differentiation │ Slope Expansion │ -├─────┼──────────┼───────────────────────────┼─────────────────┤ -│ 1 │ f1 │ 8.847e-6 │ 2.4192e-5 │ -│ 2 │ f2 │ 3.0109e-5 │ 7.6821e-5 │ -│ 3 │ f3 │ 6.1046e-6 │ 9.447e-6 │ -│ 4 │ f4 │ 1.6145e-5 │ 3.8827e-5 │ -│ 5 │ f5 │ 8.29233e-6 │ 2.2186e-5 │ -│ 6 │ f6 │ 3.0895e-5 │ 7.9524e-5 │ -│ 7 │ f7 │ 3.0847e-5 │ 7.6188e-5 │ - -│ Row │ Function │ Automatic Differentiation │ Slope Expansion │ -├─────┼──────────┼───────────────────────────┼─────────────────┤ -│ 1 │ f1 │ 8.34067e-6 │ 2.1679e-5 │ -│ 2 │ f2 │ 2.9766e-5 │ 7.194e-5 │ -│ 3 │ f3 │ 6.0278e-6 │ 7.827e-6 │ -│ 4 │ f4 │ 1.7114e-5 │ 3.5235e-5 │ -│ 5 │ f5 │ 8.36567e-6 │ 2.0747e-5 │ -│ 6 │ f6 │ 3.0061e-5 │ 7.273e-5 │ -│ 7 │ f7 │ 3.0198e-5 │ 7.2014e-5 │ -│ 8 │ f8 │ 7.43125e-6 │ 8.046e-6 │ -│ 9 │ f9 │ 7.118e-6 │ 7.8705e-6 │ diff --git a/src/IntervalRootFinding.jl b/src/IntervalRootFinding.jl index 3797cbf..c009c74 100644 --- a/src/IntervalRootFinding.jl +++ b/src/IntervalRootFinding.jl @@ -7,6 +7,7 @@ using IntervalArithmetic.Symbols using BranchAndPrune using ForwardDiff using ForwardDiff: derivative, jacobian +using InteractiveUtils using StaticArrays using LinearAlgebra diff --git a/src/roots.jl b/src/roots.jl index d05ade1..1762ff7 100644 --- a/src/roots.jl +++ b/src/roots.jl @@ -27,7 +27,7 @@ The returned `RootProblem` is an iterator that give access to the internal state of the search during the iteration, allowing to add callbacks and logging to the search. -Parameters +Keyword parameters ========== - `contractor`: Contractor used to determine the status of a region. Must be either `Newton`, `Krawczyk`, or `Bisection`. `Bisection` do not require @@ -46,6 +46,11 @@ Parameters Default: `0.0`. - `max_iteration`: The maximum number of iteration, which also corresponds to the maximum number of bisections allowed. Default: `100_000`. +- `infer_root_type`: When true, use the return type of the function as + type for the region in the returned roots, avoiding extra conversions + during the computation. + Otherwise, use the type of the provided region. + Default: `true`. Always `false` if the initial region is given as a `Root`. - `where_bisect`: Value used to bisect the region. It is used to avoid bisecting exactly on zero when starting with symmetrical regions, often leading to having a solution directly on the boundary of a region, @@ -55,7 +60,17 @@ Parameters further. Default: `[IntervalArithmetic.InconclusiveBooleanOperation]`. """ -RootProblem(f, region ; kwargs...) = RootProblem(f, Root(region, :unkown) ; kwargs...) +function RootProblem(f, region ; infer_root_type = true, kwargs...) + if infer_root_type + T = last(InteractiveUtils.@code_typed(f(region))) + if isconcretetype(T) + region = convert(T, region) + else + @warn "Could not infer the return type of the function (it may be type instable). Got $T" + end + end + RootProblem(f, Root(region, :unkown) ; kwargs...) +end function RootProblem( f, root::Root ; diff --git a/test/roots.jl b/test/roots.jl index c6878a1..6bbfa36 100644 --- a/test/roots.jl +++ b/test/roots.jl @@ -1,4 +1,3 @@ - using IntervalArithmetic, IntervalRootFinding, StaticArrays, ForwardDiff using Test @@ -245,46 +244,54 @@ end end end -@testset "Type stability" begin - f(x) = x^2 - 2 - df(x) = 2x +function f_vec(X) + x, y = X + return [x^2 - 2, y^2 - 2] +end - function A(X) - x, y = X - return [x^2 - 2, y^2 - 2] - end +function df_vec(X) + x, y = X + return [2x 0.0 ; 0.0 2y] +end - function dA(X) - x, y = X - return [2x 0.0 ; 0.0 2y] - end +function f_svec(X) + x, y = X + return SVector(x^2 - 2, y^2 - 2) +end - function S(X) - x, y = X - return SVector(x^2 - 2, y^2 - 2) - end +function df_svec(X) + x, y = X + return SMatrix{2, 2}(2x, 0.0, 0.0, 2y) +end - function dS(X) - x, y = X - return SMatrix{2, 2}(2x, 0.0, 0.0, 2y) - end - x = interval(0, 5) - a = [x, x] - s = SVector(x, x) +@testset "Type stability" begin + f(x) = x^2 - 2 + df(x) = 2x + + X = interval(0, 5) + x = Root(X, :unknown) + x_vec = Root([X, X], :unknown) + x_svec = Root(SVector(X, X), :unknown) for contractor in newtonlike_methods - @inferred roots(f, x ; contractor) - @inferred roots(f, x ; contractor, derivative = df) - - rA1 = @inferred roots(A, a ; contractor) - rA2 = @inferred roots(A, a ; contractor, derivative = dA) - @test eltype(rA1) <: Root{<:Vector} - @test eltype(rA2) <: Root{<:Vector} - - rS1 = @inferred roots(S, s ; contractor) - rS2 = @inferred roots(S, s ; contractor, derivative = dS) - @test eltype(rS1) <: Root{<:SVector} - @test eltype(rS2) <: Root{<:SVector} + @testset "$contractor" begin + @inferred roots(f, x ; contractor) + @inferred roots(f, x ; contractor, derivative = df) + + rts_vec1 = @inferred roots(f_vec, x_vec ; contractor) + rts_vec2 = @inferred roots(f_vec, x_vec ; contractor, derivative = df_vec) + @test eltype(rts_vec1) <: Root{<:Vector} + @test eltype(rts_vec2) <: Root{<:Vector} + + rts_svec1 = @inferred roots(f_svec, x_svec ; contractor) + rts_svec2 = @inferred roots(f_svec, x_svec ; contractor, derivative = df_svec) + @test eltype(rts_svec1) <: Root{<:SVector} + @test eltype(rts_svec2) <: Root{<:SVector} + + # Check that the root type is determined by the return type of the function + @test eltype(roots(f_vec, SVector(X, X) ; contractor)) <: Root{<:Vector} + @test eltype(roots(f_svec, [X, X] ; contractor)) <: Root{<:SVector} + end end end