diff --git a/examples/linear_eq.jl b/examples/linear_eq.jl new file mode 100644 index 00000000..918bcb69 --- /dev/null +++ b/examples/linear_eq.jl @@ -0,0 +1,32 @@ +# Examples from Luc Jaulin, Michel Kieffer, Olivier Didrit and Eric Walter - Applied Interval Analysis + +using IntervalArithmetic, IntervalRootFinding, StaticArrays + +A = [4..5 -1..1 1.5..2.5; -0.5..0.5 -7.. -5 1..2; -1.5.. -0.5 -0.7.. -0.5 2..3] +sA = SMatrix{3}{3}(A) +mA = MMatrix{3}{3}(A) + +b = [3..4, 0..2, 3..4] +sb = SVector{3}(b) +mb = MVector{3}(b) + +p = fill(-1e16..1e16, 3) + +rts = gauss_seidel_interval!(p, A, b, precondition=true) # Gauss-Seidel Method; precondition=true by default +rts = gauss_seidel_interval!(p, sA, sb, precondition=true) # Gauss-Seidel Method; precondition=true by default +rts = gauss_seidel_interval!(p, mA, mb, precondition=true) # Gauss-Seidel Method; precondition=true by default + +rts = gauss_seidel_interval(A, b, precondition=true) # Gauss-Seidel Method; precondition=true by default +rts = gauss_seidel_interval(sA, sb, precondition=true) # Gauss-Seidel Method; precondition=true by default +rts = gauss_seidel_interval(mA, mb, precondition=true) # Gauss-Seidel Method; precondition=true by default + +rts = gauss_seidel_contractor!(p, A, b, precondition=true) # Gauss-Seidel Method (Vectorized); precondition=true by default +rts = gauss_seidel_contractor!(p, sA, sb, precondition=true) # Gauss-Seidel Method (Vectorized); precondition=true by default +rts = gauss_seidel_contractor!(p, mA, mb, precondition=true) # Gauss-Seidel Method (Vectorized); precondition=true by default + +rts = gauss_seidel_contractor(A, b, precondition=true) # Gauss-Seidel Method (Vectorized); precondition=true by default +rts = gauss_seidel_contractor(sA, sb, precondition=true) # Gauss-Seidel Method (Vectorized); precondition=true by default +rts = gauss_seidel_contractor(mA, mb, precondition=true) # Gauss-Seidel Method (Vectorized); precondition=true by default + +rts = gauss_elimination_interval!(p, A, b, precondition=true) # Gaussian Elimination; precondition=true by default +rts = gauss_elimination_interval(A, b, precondition=true) # Gaussian Elimination; precondition=true by default diff --git a/perf/linear_eq.jl b/perf/linear_eq.jl new file mode 100644 index 00000000..de582343 --- /dev/null +++ b/perf/linear_eq.jl @@ -0,0 +1,64 @@ +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.\\", "Linear Hull"] + 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) + t3 = @belapsed linear_hull($A, $b) + df[Symbol("$n")] = [t1, t2, t3] + 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 new file mode 100644 index 00000000..ff08cefe --- /dev/null +++ b/perf/linear_eq_results.txt @@ -0,0 +1,149 @@ +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 new file mode 100644 index 00000000..f06bfc27 --- /dev/null +++ b/perf/linear_eq_tabular.txt @@ -0,0 +1,53 @@ +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 + +│ Row │ n │ Base.\\ │ Gauss Elimination │ Linear Hull │ +├─────┼────┼────────────┼───────────────────┼─────────────┤ +│ 1 │ 1 │ 1.9415e-6 │ 8.40867e-6 │ 3.94457e-6 │ +│ 2 │ 10 │ 7.1184e-5 │ 7.8812e-5 │ 4.6999e-5 │ +│ 3 │ 2 │ 3.34387e-6 │ 1.012e-5 │ 5.30567e-6 │ +│ 4 │ 3 │ 4.548e-6 │ 1.0955e-5 │ 6.5832e-6 │ +│ 5 │ 4 │ 7.256e-6 │ 1.3845e-5 │ 1.1115e-5 │ +│ 6 │ 5 │ 1.0878e-5 │ 1.7847e-5 │ 1.4564e-5 │ +│ 7 │ 6 │ 1.5747e-5 │ 2.1526e-5 │ 1.9468e-5 │ +│ 8 │ 7 │ 2.296e-5 │ 3.109e-5 │ 2.7597e-5 │ +│ 9 │ 8 │ 3.4512e-5 │ 4.0778e-5 │ 3.7466e-5 │ +│ 10 │ 9 │ 4.8627e-5 │ 5.7028e-5 │ 5.3766e-5 │ diff --git a/perf/slopes.jl b/perf/slopes.jl new file mode 100644 index 00000000..642a8d84 --- /dev/null +++ b/perf/slopes.jl @@ -0,0 +1,116 @@ +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 + +struct SlopesMulti + f::Function + x::IntervalBox + c::Vector + sol::Matrix{Interval} +end + +function benchmark_multi() + + rts = SlopesMulti[] + f(x, y) = SVector(x^2 + y^2 - 1, y - 2x) + f(X) = f(X...) + X = (-6..6) × (-6..6) + c = [0.0, 0.0] + push!(rts, SlopesMulti(f, X, c, [-6..6 -6..6; -2.. -2 1..1])) + + function g(x) + (x1, x2, x3) = x + SVector( x1^2 + x2^2 + x3^2 - 1, + x1^2 + x3^2 - 0.25, + x1^2 + x2^2 - 4x3 + ) + end + + X = (-5..5) + XX = IntervalBox(X, 3) + cc = [0.0, 0.0, 0.0] + push!(rts, SlopesMulti(g, XX, cc, [-5..5 -5..5 -5..5; -5..5 0..0 -5..5; -5..5 -5..5 -4.. -4])) + + function h(x) + (x1, x2, x3) = x + SVector( x1 + x2^2 + x3^2 - 1, + x1^2 + x3 - 0.25, + x1^2 + x2 - 4x3 + ) + end + + XXX = IntervalBox(1..5, 2..6, -3..7) + ccc = [3.0, 4.0, 2.0] + push!(rts, SlopesMulti(h, XXX, ccc, [1..1 6..10 -1..9; 4..8 0..0 1..1; 4..8 1..1 -4.. -4])) + + for i in 1:length(rts) + println("\nFunction $i") + println("Slope Expansion: ") + println(DataFrame(slope(rts[i].f, rts[i].x, rts[i].c))) + println("\nJacobian: ") + println(DataFrame(ForwardDiff.jacobian(rts[i].f, rts[i].x))) + end +end diff --git a/perf/slopes_tabular.txt b/perf/slopes_tabular.txt new file mode 100644 index 00000000..c8ac3c28 --- /dev/null +++ b/perf/slopes_tabular.txt @@ -0,0 +1,82 @@ +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 │ + + +Multi-Dimensional Case: +Function 1 +Slope Expansion: +│ Row │ x1 │ x2 │ +├─────┼──────────┼─────────┤ +│ 1 │ [-6, 6] │ [-6, 6] │ +│ 2 │ [-2, -2] │ [1, 1] │ + +Jacobian: +│ Row │ x1 │ x2 │ +├─────┼───────────┼───────────┤ +│ 1 │ [-12, 12] │ [-12, 12] │ +│ 2 │ [-2, -2] │ [1, 1] │ + +Function 2 +Slope Expansion: +│ Row │ x1 │ x2 │ x3 │ +├─────┼─────────┼─────────┼──────────┤ +│ 1 │ [-5, 5] │ [-5, 5] │ [-5, 5] │ +│ 2 │ [-5, 5] │ [0, 0] │ [-5, 5] │ +│ 3 │ [-5, 5] │ [-5, 5] │ [-4, -4] │ + +Jacobian: +│ Row │ x1 │ x2 │ x3 │ +├─────┼───────────┼───────────┼───────────┤ +│ 1 │ [-10, 10] │ [-10, 10] │ [-10, 10] │ +│ 2 │ [-10, 10] │ [0, 0] │ [-10, 10] │ +│ 3 │ [-10, 10] │ [-10, 10] │ [-4, -4] │ + +Function 3 +Slope Expansion: +│ Row │ x1 │ x2 │ x3 │ +├─────┼────────┼─────────┼──────────┤ +│ 1 │ [1, 1] │ [6, 10] │ [-1, 9] │ +│ 2 │ [4, 8] │ [0, 0] │ [1, 1] │ +│ 3 │ [4, 8] │ [1, 1] │ [-4, -4] │ + +Jacobian: +│ Row │ x1 │ x2 │ x3 │ +├─────┼─────────┼─────────┼──────────┤ +│ 1 │ [1, 1] │ [4, 12] │ [-6, 14] │ +│ 2 │ [2, 10] │ [0, 0] │ [1, 1] │ +│ 3 │ [2, 10] │ [1, 1] │ [-4, -4] │ diff --git a/src/IntervalRootFinding.jl b/src/IntervalRootFinding.jl index 0496a574..12fc586d 100644 --- a/src/IntervalRootFinding.jl +++ b/src/IntervalRootFinding.jl @@ -18,12 +18,16 @@ export derivative, jacobian, # reexport derivative from ForwardDiff Root, is_unique, roots, find_roots, - bisect, newton1d + bisect, newton1d, slope, + slope_newton1d, linear_hull, + gauss_seidel_interval, gauss_seidel_interval!, + gauss_seidel_contractor, gauss_seidel_contractor!, + gauss_elimination_interval, gauss_elimination_interval! export isunique, root_status -import IntervalArithmetic.interval +import IntervalArithmetic: interval, wideinterval @@ -58,7 +62,8 @@ include("complex.jl") include("contractors.jl") include("roots.jl") include("newton1d.jl") - +include("linear_eq.jl") +include("slopes.jl") gradient(f) = X -> ForwardDiff.gradient(f, SVector(X)) diff --git a/src/linear_eq.jl b/src/linear_eq.jl new file mode 100644 index 00000000..70bbd5e4 --- /dev/null +++ b/src/linear_eq.jl @@ -0,0 +1,188 @@ +""" +Preconditions the matrix A and b with the inverse of mid(A) +""" +function preconditioner(A::AbstractMatrix, b::AbstractArray) + + Aᶜ = mid.(A) + B = inv(Aᶜ) + + return B*A, B*b + +end + +function gauss_seidel_interval(A::AbstractMatrix, b::AbstractArray; precondition=true, maxiter=100) + + n = size(A, 1) + x = similar(b) + x .= -1e16..1e16 + gauss_seidel_interval!(x, A, b, precondition=precondition, maxiter=maxiter) + return x +end +""" +Iteratively solves the system of interval linear +equations and returns the solution set. Uses the +Gauss-Seidel method (Hansen-Sengupta version) to solve the system. +Keyword `precondition` to turn preconditioning off. +Eldon Hansen and G. William Walster : Global Optimization Using Interval Analysis - Chapter 5 - Page 115 +""" +function gauss_seidel_interval!(x::AbstractArray, A::AbstractMatrix, b::AbstractArray; precondition=true, maxiter=100) + + precondition && ((A, b) = preconditioner(A, b)) + + n = size(A, 1) + + @inbounds for iter in 1:maxiter + for i in 1:n + Y = b[i] + for j in 1:n + (i == j) || (Y -= A[i, j] * x[j]) + end + Z = extended_div(Y, A[i, i]) + x[i] = hull((x[i] ∩ Z[1]), x[i] ∩ Z[2]) + end + end + x +end + +function gauss_seidel_contractor(A::AbstractMatrix, b::AbstractArray; precondition=true, maxiter=100) + + n = size(A, 1) + x = similar(b) + x .= -1e16..1e16 + x = gauss_seidel_contractor!(x, A, b, precondition=precondition, maxiter=maxiter) + return x +end + +function gauss_seidel_contractor!(x::AbstractArray, A::AbstractMatrix, b::AbstractArray; precondition=true, maxiter=100) + + precondition && ((A, b) = preconditioner(A, b)) + + n = size(A, 1) + + diagA = Diagonal(A) + extdiagA = copy(A) + for i in 1:n + if (typeof(b) <: SArray) + extdiagA = setindex(extdiagA, Interval(0), i, i) + else + extdiagA[i, i] = Interval(0) + end + end + inv_diagA = inv(diagA) + + for iter in 1:maxiter + x = x .∩ (inv_diagA * (b - extdiagA * x)) + end + x +end + +function gauss_elimination_interval(A::AbstractMatrix, b::AbstractArray; precondition=true) + + n = size(A, 1) + x = similar(b) + x .= -1e16..1e16 + x = gauss_elimination_interval!(x, A, b, precondition=precondition) + + return x +end +""" +Solves the system of linear equations using Gaussian Elimination, +with (or without) preconditioning. (kwarg - `precondition`) +Luc Jaulin, Michel Kieffer, Olivier Didrit and Eric Walter - Applied Interval Analysis - Page 72 +""" +function gauss_elimination_interval!(x::AbstractArray, a::AbstractMatrix, b::AbstractArray; precondition=true) + + precondition && ((a, b) = preconditioner(a, b)) + + n = size(a, 1) + + p = zeros(x) + + for i in 1:(n-1) + if 0 ∈ a[i, i] # diagonal matrix is not invertible + p .= entireinterval(b[1]) + return p .∩ x + end + + for j in (i+1):n + α = a[j, i] / a[i, i] + b[j] -= α * b[i] + + for k in (i+1):n + a[j, k] -= α * a[i, k] + end + end + end + + for i in n:-1:1 + sum = 0 + for j in (i+1):n + sum += a[i, j] * p[j] + end + p[i] = (b[i] - sum) / a[i, i] + end + + p .∩ x +end + +function gauss_elimination_interval1(A::AbstractMatrix, b::AbstractArray; precondition=true) + + n = size(A, 1) + x = fill(-1e16..1e16, n) + + x = gauss_elimination_interval1!(x, A, b, precondition=precondition) + + return x +end +""" +Using `Base.\`` +""" +function gauss_elimination_interval1!(x::AbstractArray, a::AbstractMatrix, b::AbstractArray; precondition=true) + + precondition && ((a, b) = preconditioner(a, b)) + + a \ b +end + +""" +Method to compute a tight hull of the solution set of a regular system of interval linear equations. +Reference - Eldon Hansen and G. William Walster : Global Optimization Using Interval Analysis - Chapter 5 - Section 5.8 +""" +function linear_hull(M::AbstractMatrix, r::AbstractArray) + + n = size(M, 1) + + ((M, r) = preconditioner(M, r)) + M = interval.((2*eye(n) - sup.(M)), sup.(M)) + M_lo = inf.(M) + M_hi = sup.(M) + P = eye(n) + if all(.≤(M_lo, zero(M_lo))) + return M \ r + end + P = inv(M_lo) + if any(isinf.(P)) || any(isnan.(P)) + return gauss_seidel_interval(M, r, precondition=false, maxiter=2) + end + if all(.≤(eye(n), (P))) + H1 = P * sup.(r) + C = 1 ./ (2diag(P) - 1) + Z = ((2mid.(r)) .* diag(P)) - H1 + H2 = C .* Z + for i in 1:n + if Z[i] < 0 + #H2[i] = Z[i] + if (typeof(M) <: SArray) + H2 = setindex(H2, Z[i], i) + else + H2[i] = Z[i] + end + end + end + H = interval.(min.(H1, H2), max.(H1, H2)) + + return H + else + return M \ r + end +end diff --git a/src/newton1d.jl b/src/newton1d.jl index cac4be47..9b3dea4a 100644 --- a/src/newton1d.jl +++ b/src/newton1d.jl @@ -5,61 +5,137 @@ Optional keyword arguments give the tolerances `reltol` and `abstol`. and a `debug` boolean argument that prints out diagnostic information.""" function newton1d{T}(f::Function, f′::Function, x::Interval{T}; - reltol=eps(T), abstol=eps(T), debug=false) + reltol=eps(T), abstol=eps(T), debug=false, debugroot=false) - L = Interval{T}[] + L = Interval{T}[] # Array to hold the intervals still to be processed - R = Root{Interval{T}}[] + R = Root{Interval{T}}[] # Array to hold the `root` objects obtained + reps = reps1 = 0 - push!(L, x) + push!(L, x) # Initialize + initial_width = ∞ + X = emptyinterval(T) # Initialize + while !isempty(L) # Until all intervals have been processed + X = pop!(L) # Process next interval + + debug && (print("Current interval popped: "); @show X) - while !isempty(L) - X = pop!(L) m = mid(X) if (isempty(X)) continue end - if 0 ∉ f′(X) + if 0 ∉ ForwardDiff.derivative(f, X) + + debug && println("0 ∉ f′(X)") + while true + m = mid(X) - N = m - (f(Interval(m)) / f′(X)) + N = m - (f(interval(m)) / f′(X)) # Newton step + + debug && (print("Newton step1: "); @show (X, X ∩ N)) + if X == X ∩ N # Checking if Newton step was redundant + reps1 += 1 + if reps1 > 20 + reps1 = 0 + break + end + end X = X ∩ N - if isempty(X) + if (isempty(X)) # No root in X break - elseif 0 ∈ f(Interval(prevfloat(m), nextfloat(m))) - push!(R, Root(X, :unique)) + elseif 0 ∈ f(wideinterval(mid(X))) # Root guaranteed to be in X + n = fa = fb = 0 + root_exist = true + while (n < 4 && (fa == 0 || fb == 0)) # Narrowing the interval further + if fa == 0 + if 0 ∈ f(wideinterval(X.lo)) + fa = 1 + else + N = X.lo - (f(interval(X.lo)) / f′(X)) + X = X ∩ N + if (isempty(X)) + root_exist = false + break + end + end + end + if fb == 0 + if 0 ∈ f(wideinterval(X.hi)) + fb = 1 + else + if 0 ∈ f(wideinterval(mid(X))) + N = X.hi - (f(interval(X.hi)) / f′(X)) + else + N = mid(X) - (f(interval(mid(X))) / f′(X)) + end + X = X ∩ N + if (isempty(X)) + root_exist = false + break + end + end + end + N = mid(X) - (f(interval(mid(X))) / f′(X)) + X = X ∩ N + if (isempty(X)) + root_exist = false + break + end + n += 1 + end + if root_exist + push!(R, Root(X, :unique)) + debugroot && @show "Root found", X # Storing determined unique root + end + break end end - else - # 0 ∈ f'(X) + else # 0 ∈ f′(X) + if diam(X) == initial_width # if no improvement occuring for a number of iterations + reps += 1 + if reps > 10 + push!(R, Root(X, :unknown)) + debugroot && @show "Repeated root found", X + reps = 0 + continue + end + end + initial_width = diam(X) + debug && println("0 ∈ f'(X)") + expansion_pt = Inf # expansion point for the newton step might be m, X.lo or X.hi according to some conditions - if 0 ∈ f(Interval(mid(X))) + if 0 ∈ f(wideinterval(mid(X))) # if root in X, narrow interval further # 0 ∈ fⁱ(x) - # Step 7 - if 0 ∉ f(Interval(X.lo)) + debug && println("0 ∈ fⁱ(x)") + + if 0 ∉ f(wideinterval(X.lo)) expansion_pt = X.lo - elseif 0 ∉ f(Interval(X.hi)) + elseif 0 ∉ f(wideinterval(X.hi)) expansion_pt = X.hi else - x1 = mid(Interval(X.lo, mid(X))) - x2 = mid(Interval(mid(X), X.hi)) - if 0 ∉ f(Interval(x1)) || 0 ∉ f(Interval(x2)) - push!(L, Interval(X.lo, m)) - push!(L, Interval(m, X.hi)) + x1 = mid(interval(X.lo, mid(X))) + x2 = mid(interval(mid(X), X.hi)) + if 0 ∉ f(wideinterval(x1)) || 0 ∉ f(wideinterval(x2)) + push!(L, interval(X.lo, m)) + push!(L, interval(m, X.hi)) continue else - push!(R, Root(X, :unique)) + push!(R, Root(X, :unknown)) + + debugroot && @show "Multiple root found", X + continue end end @@ -67,8 +143,13 @@ function newton1d{T}(f::Function, f′::Function, x::Interval{T}; else # 0 ∉ fⁱ(x) - if (diam(X)/mag(X)) < reltol && diam(f(X)) < abstol + debug && println("0 ∉ fⁱ(x)") + + if (diam(X)/mag(X)) < reltol && diam(f(X)) < abstol # checking if X is still within tolerances push!(R, Root(X, :unknown)) + + debugroot && @show "Tolerance root found", X + continue end end @@ -78,26 +159,28 @@ function newton1d{T}(f::Function, f′::Function, x::Interval{T}; expansion_pt = mid(X) end - initial_width = diam(X) - a = f(Interval(expansion_pt)) - b = f′(X) + a = f(interval(expansion_pt)) + b = ForwardDiff.derivative(f, X) if 0 < b.hi && 0 > b.lo && 0 ∉ a if a.hi < 0 - push!(L, X ∩ (expansion_pt - Interval(-Inf, a.hi / b.hi))) - push!(L, X ∩ (expansion_pt - Interval(a.hi / b.lo, Inf))) + push!(L, X ∩ (expansion_pt - interval(-Inf, a.hi / b.hi))) + push!(L, X ∩ (expansion_pt - interval(a.hi / b.lo, Inf))) elseif a.lo > 0 - push!(L, X ∩ (expansion_pt - Interval(-Inf, a.lo / b.lo))) - push!(L, X ∩ (expansion_pt - Interval(a.lo / b.hi, Inf))) + push!(L, X ∩ (expansion_pt - interval(-Inf, a.lo / b.lo))) + push!(L, X ∩ (expansion_pt - interval(a.lo / b.hi, Inf))) end continue else - N = expansion_pt - (f(Interval(expansion_pt))/f′(X)) + N = expansion_pt - (f(interval(expansion_pt))/f′(X)) + + debug && (print("Newton step2: "); @show (X, X ∩ N)) + X = X ∩ N m = mid(X) @@ -106,12 +189,7 @@ function newton1d{T}(f::Function, f′::Function, x::Interval{T}; end end - if diam(X) > initial_width/2 - push!(L, Interval(m, X.hi)) - X = Interval(X.lo, m) - end - - push!(L, X) + push!(L, X) # Pushing X into L to be processed again end end @@ -126,3 +204,8 @@ and a `debug` boolean argument that prints out diagnostic information.""" newton1d{T}(f::Function, x::Interval{T}; args...) = newton1d(f, x->D(f,x), x; args...) + +function slope_newton1d{T}(f::Function, x::Interval{T}; + reltol=eps(T), abstol=eps(T), debug=false, debugroot=false) + newton1d(f, x->slope(f, x), x, reltol=reltol, abstol=abstol, debug=debug, debugroot=debugroot) +end diff --git a/src/newtonnd.jl b/src/newtonnd.jl new file mode 100644 index 00000000..7d892f52 --- /dev/null +++ b/src/newtonnd.jl @@ -0,0 +1,126 @@ +""" +Preconditions the matrix A and b with the inverse of mid(A) +""" +function preconditioner(A::AbstractMatrix, b::AbstractArray) + + Aᶜ = mid.(A) + B = inv(Aᶜ) # TODO If Aᶜ is singular + + return B*A, B*b + +end + +function newtonnd(f::Function, deriv::Function, x::IntervalBox{NUM, T}; reltol=eps(T), abstol=eps(T), debug=false, debugroot=false) where {T<:AbstractFloat} where {NUM} # TODO Incorporate Hull and Box consistencies + + L = IntervalBox{NUM, T}[] # Array to hold the interval boxes still to be processed + + R = Root{IntervalBox{NUM, T}}[] # Array to hold the `root` objects obtained + + push!(L, x) # Initialize + n = size(X, 1) + while !isempty(L) # Until all interval boxes have been processed + Xᴵ = pop!(L) # Process next interval box + if !all(0 .∈ f(Xᴵ)) + continue + end + + Xᴵ¹ = copy(Xᴵ) + use_B = false + + debug && (print("Current interval popped: "); println(Xᴵ)) + + if (isempty(Xᴵ)) + continue + end + if diam(Xᴵ) < reltol + if max((abs.(IntervalBox(f(Xᴵ))))...) < abstol + + debugroot && @show "Tolerance root found", Xᴵ + + push!(R, Root(Xᴵ, :unknown)) + continue + else + continue + end + end + + next_iter = false + + while true + + next_iter = false + + initial_width = diam(Xᴵ) + debug && (print("Current interval popped: "); println(Xᴵ)) + X = mid(Xᴵ) + if use_B + for i in 1:3 + z = X .- B * f(X) + if all(z .∈ Xᴵ) + if max(abs.(f(z))...) ≥ max(abs.(f(X))...) + break + end + X = z + else + break + end + end + if any(X .∉ Xᴵ) + X = mid.(Xᴵ) + end + end + + J = SMatrix{n}{n}(deriv(Xᴵ)) # either jacobian or calculated using slopes + B, r = preconditioner(J, -f(interval.(X))) + # Xᴵ = IntervalBox((X + linear_hull(B, r, precondition=false)) .∩ Xᴵ) + Xᴵ = IntervalBox((X + (J \ -f(interval.(X)))) .∩ Xᴵ) + use_B = true + if (isempty(Xᴵ)) + next_iter = true + break + end + + if diam(Xᴵ) < reltol + if max((abs.(IntervalBox(f(Xᴵ))))...) < abstol + + debugroot && @show "Tolerance root found", Xᴵ + next_iter = true + push!(R, Root(Xᴵ, :unknown)) + break + else + next_iter = true + break + end + end + + if all(0 .∈ f(interval.(X))) && initial_width > 0.9diam(Xᴵ) + next_iter = true + push!(R, Root(Xᴵ, :unique)) + break + end + + if diam(Xᴵ) > initial_width / 8 + break + end + + #Criterion C + + + end + + if next_iter + continue + end + + if 0.25 * diam(Xᴵ¹) ≤ max((diam.(Xᴵ¹) .- diam.(Xᴵ))...) + push!(L, Xᴵ) + else + push!(L, bisect(Xᴵ)...) + end + end + + R + +end + +newtonnd(f::Function, x::IntervalBox{NUM, T}; args...) where {T<:AbstractFloat} where {NUM} = newtonnd(f, x->ForwardDiff.jacobian(f,x), x; args...) diff --git a/src/slopes.jl b/src/slopes.jl new file mode 100644 index 00000000..28345edc --- /dev/null +++ b/src/slopes.jl @@ -0,0 +1,228 @@ +# Reference : Dietmar Ratz - An Optimized Interval Slope Arithmetic and its Application +using IntervalArithmetic, ForwardDiff, StaticArrays +import Base: +, -, *, /, ^, sqrt, exp, log, sin, cos, tan, asin, acos, atan +import IntervalArithmetic: mid, interval + +""" +Expands the slope of `f` over `x` at the point `c` (default `c = mid(x)`) +""" +function slope(f::Function, x::Interval, c::Number = mid(x)) + try + f(slope_var(x, c)).s + catch y + if isa(y, MethodError) + ForwardDiff.derivative(f, x) + end + end +end + +function slope(f::Function, x::Union{IntervalBox, SVector}, c::AbstractVector = mid.(x)) # multidim + try + T = typeof(x[1].lo) + n = length(x) + s = Vector{Slope{T}}[] + for i in 1:n + arr = fill(Slope(zero(T)), n) + arr[i] = slope_var(x[i], c[i]) + push!(s, arr) + end + return slope.(hcat(([(f(s[i])) for i=1:n])...)) + catch y + if isa(y, MethodError) + ForwardDiff.jacobian(f, x) + end + end + +end + +struct Slope{T} + x::Interval{T} + c::Interval{T} + s::Interval{T} +end + +Slope(c) = Slope(c, c, 0) +Slope(a, b, c) = Slope(promote(convert(Interval, a), b, c)...) + +function slope_var(v::Number) + Slope(v, v, 1) +end + +function slope_var(v::Interval, c::Number) + Slope(v, c, 1) +end + +function interval(u::Slope) + u.x +end + +function mid(u::Slope) + u.c +end + +function slope(u::Slope) + u.s +end + +function +(u::Slope, v::Slope) + Slope(u.x + v.x, u.c + v.c, u.s + v.s) +end + +function -(u::Slope, v::Slope) + Slope(u.x - v.x, u.c - v.c, u.s - v.s) +end + +function *(u::Slope, v::Slope) + Slope(u.x * v.x, u.c * v.c, u.s * v.c + u.x * v.s) +end + +function /(u::Slope, v::Slope) + Slope(u.x / v.x, u.c / v.c, (u.s - (u.c / v.c) * v.s) / v.x) +end + +function +(u, v::Slope) + Slope(u + v.x, u + v.c, v.s) +end + +function -(u, v::Slope) + Slope(u - v.x, u - v.c, -v.s) +end + +function *(u, v::Slope) + Slope(u * v.x, u * v.c, u * v.s) +end + +function /(u, v::Slope) + Slope(u / v.x, u / v.c, -(u / v.c) * (v.s / v.x)) +end + ++(v::Slope, u) = u + v + +*(v::Slope, u) = u * v + +function -(u::Slope, v) + Slope(u.x - v, u.c - v, u.s) +end + +function -(u::Slope) + Slope(-u.x, -u.c, -u.s) +end + +function /(u::Slope, v) + Slope(u.x / v, u.c / v, u.s / v) +end + +function sqr(u::Slope) + Slope(u.x ^ 2, u.c ^ 2, (u.x + u.c) * u.s) +end + +function ^(u::Slope, k::Integer) + if k == 0 + return Slope(1) + elseif k == 1 + return u + elseif k == 2 + return sqr(u) + else + hxi = interval(u.x.lo) ^ k + hxs = interval(u.x.hi) ^ k + hx = hull(hxi, hxs) + + if (k % 2 == 0) && (0 ∈ u.x) + hx = interval(0, hx.hi) + end + + hc = u.c ^ k + + i = u.x.lo - u.c.lo + s = u.x.hi - u.c.hi + + if ((i == 0) || (s == 0) || (k % 2 == 1 && zero(u.x) ⪽ u.x)) + h1 = k * (u.x ^ (k - 1)) + else + if k % 2 == 0 || u.x.lo >= 0 + h1 = interval((hxi.hi - hc.lo) / i, (hxs.hi - hc.lo) / s) + else + h1 = interval((hxs.lo - hc.hi) / s, (hxi.lo - hc.hi) / i) + end + end + return Slope(hx, hc, h1 * u.s) + end +end + +function sqrt(u::Slope) + Slope(sqrt(u.x), sqrt(u.c), u.s / (sqrt(u.x) + sqrt(u.c))) +end + +function exp(u::Slope) + hx = exp(u.x) + hc = exp(u.c) + + i = u.x.lo - u.c.lo + s = u.x.hi - u.c.hi + + if (i == 0 || s == 0) + h1 = hx + else + h1 = interval((hx.lo - hc.lo) / i, (hx.hi - hc.hi) / s) + end + + Slope(hx, hc, h1 * u.s) +end + +function log(u::Slope) + hx = log(u.x) + hc = log(u.c) + + i = u.x.lo - u.c.lo + s = u.x.hi - u.c.hi + + if (i == 0 || s == 0) + h1 = 1 / u.x + else + h1 = interval((hx.hi - hc.hi) / s, (hx.lo - hc.lo) / i) + end + Slope(hx, hc, h1 * u.s) +end + +function sin(u::Slope) # Using derivative to upper bound the slope expansion for now + hx = sin(u.x) + hc = sin(u.c) + hs = cos(u.x) + Slope(hx, hc, hs) +end + +function cos(u::Slope) # Using derivative to upper bound the slope expansion for now + hx = cos(u.x) + hc = cos(u.c) + hs = -sin(u.x) + Slope(hx, hc, hs) +end + +function tan(u::Slope) # Using derivative to upper bound the slope expansion for now + hx = tan(u.x) + hc = tan(u.c) + hs = (sec(u.x)) ^ 2 + Slope(hx, hc, hs) +end + +function asin(u::Slope) + hx = asin(u.x) + hc = asin(u.c) + hs = 1 / sqrt(1 - (u.x ^ 2)) + Slope(hx, hc, hs) +end + +function acos(u::Slope) + hx = acos(u.x) + hc = acos(u.c) + hs = -1 / sqrt(1 - (u.x ^ 2)) + Slope(hx, hc, hs) +end + +function atan(u::Slope) + hx = atan(u.x) + hc = atan(u.c) + hs = 1 / 1 + (u.x ^ 2) + Slope(hx, hc, hs) +end diff --git a/test/linear_eq.jl b/test/linear_eq.jl new file mode 100644 index 00000000..a7ee7775 --- /dev/null +++ b/test/linear_eq.jl @@ -0,0 +1,17 @@ +using IntervalArithmetic, StaticArrays, IntervalRootFinding + +@testset "Linear Equations" begin + + A = [[2..3 0..1;1..2 2..3], ] + b = [[0..120, 60..240], ] + x = [[-120..90, -60..240], ] + + + for i in 1:length(A) + for precondition in (false, true) + for f in (gauss_seidel_interval, gauss_seidel_contractor, gauss_elimination_interval, linear_hull) + @test all(x[i] .⊆ f(A[i], b[i])) + end + end + end +end diff --git a/test/newton1d.jl b/test/newton1d.jl index 74fadcf7..2b286652 100644 --- a/test/newton1d.jl +++ b/test/newton1d.jl @@ -18,36 +18,121 @@ three_halves_pi = 3*big_pi/2 @testset "Testing newton1d" begin - f(x) = exp(x^2) - cos(x) + f(x) = exp(x^2) - cos(x) #double root f′(x) = 2*x*exp(x^2) + sin(x) - f1(x) = x^4 - 10x^3 + 35x^2 - 50x + 24 + f1(x) = x^4 - 10x^3 + 35x^2 - 50x + 24 #four unique roots f1′(x) = 4x^3 - 30x^2 + 70x - 50 - + f2(x) = 4567x^2 - 9134x + 4567 #double root + f2′(x) = 9134x - 9134 + f3(x) = (x^2 - 2)^2 #two double roots + f3′(x) = 4x * (x^2 - 2) + f4(x) = sin(x) - x #triple root at 0 + f4′(x) = cos(x) - 1 + f5(x) = (x^2 - 1)^4 * (x^2 - 2)^4 #two quadruple roots + f5′(x) = 8x * (-3 + 2x^2) * (2 - 3x^2 + x^4)^3 for autodiff in (false, true) if autodiff rts1 = newton1d(sin, -5..5) rts2 = newton1d(f, -∞..∞) rts3 = newton1d(f1, -10..10) + rts4 = newton1d(f2, -10..11) + rts5 = newton1d(f3, -10..10) + rts6 = newton1d(f4, -10..10) + rts7 = newton1d(f5, -10..10) else rts1 = newton1d(sin, cos, -5..5) rts2 = newton1d(f, f′, -∞..∞) rts3 = newton1d(f1, f1′, -10..10) + rts4 = newton1d(f2, f2′, -10..11) + rts5 = newton1d(f3, f3′, -10..10) + rts6 = newton1d(f4, f4′, -10..10) + rts7 = newton1d(f5, f5′, -10..10, reltol=0) end @test length(rts1) == 3 L = [ -pi_interval(Float64), 0..0, pi_interval(Float64)] for i = 1:length(rts1) - @test L[i] == rts1[i].interval && :unique == rts1[i].status + @test L[i] in rts1[i].interval && :unique == rts1[i].status end @test length(rts2) == 1 - @test (0..0) == rts2[1].interval && :unique == rts2[1].status + @test (0..0) == rts2[1].interval && :unknown == rts2[1].status @test length(rts3) == 4 L = [1, 2, 3, 4] for i = 1:length(rts3) @test L[i] in rts3[i].interval && :unique == rts3[i].status end + + @test length(rts4) == 1 + @test 1 in rts4[1].interval && :unknown == rts4[1].status + + L1 = [-sqrt(2), sqrt(2)] + for i = 1:length(rts5) + @test L1[i] in rts5[i].interval && :unknown == rts5[i].status + end + + @test length(rts6) == 1 + @test 0 in rts6[1].interval && :unknown == rts6[1].status + + @test length(rts7) == 4 + L = [-sqrt(2), -1, 1, sqrt(2)] + for i = 1:length(rts7) + @test L[i] in rts7[i].interval && :unknown == rts7[i].status + end + end end + +@testset "Testing slope newton1d" begin + + f(x) = exp(x^2) - cos(x) #double root + f1(x) = x^4 - 10x^3 + 35x^2 - 50x + 24 #four unique roots + f2(x) = 4567x^2 - 9134x + 4567 #double root + f3(x) = (x^2 - 2)^2 #two double roots + f4(x) = sin(x) - x #triple root at 0 + f5(x) = (x^2 - 1)^4 * (x^2 - 2)^4 #two quadruple roots + + rts1 = slope_newton1d(sin, -5..5) + rts2 = slope_newton1d(f, -∞..∞) + rts3 = slope_newton1d(f1, -10..10) + rts4 = slope_newton1d(f2, -10..11) + rts5 = slope_newton1d(f3, -10..10) + rts6 = slope_newton1d(f4, -10..10) + rts7 = slope_newton1d(f5, -10..10) + + @test length(rts1) == 3 + L = [ -pi_interval(Float64), 0..0, pi_interval(Float64)] + for i = 1:length(rts1) + @test L[i] in rts1[i].interval && :unique == rts1[i].status + end + + @test length(rts2) == 1 + @test (0..0) == rts2[1].interval && :unknown == rts2[1].status + + @test length(rts3) == 4 + L = [1, 2, 3, 4] + for i = 1:length(rts3) + @test L[i] in rts3[i].interval && :unique == rts3[i].status + end + + @test length(rts4) == 1 + @test 1 in rts4[1].interval && :unknown == rts4[1].status + + L1 = [-sqrt(2), sqrt(2)] + for i = 1:length(rts5) + @test L1[i] in rts5[i].interval && :unknown == rts5[i].status + end + + + @test length(rts6) == 1 + @test 0 in rts6[1].interval && :unknown == rts6[1].status + + @test length(rts7) == 4 + L = [-sqrt(2), -1, 1, sqrt(2)] + for i = 1:length(rts7) + @test L[i] in rts7[i].interval && :unknown == rts7[i].status + end + +end diff --git a/test/runtests.jl b/test/runtests.jl index e905e897..679cd1ad 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,3 +5,5 @@ using Base.Test #include("findroots.jl") include("roots.jl") include("newton1d.jl") +include("linear_eq.jl") +include("slopes.jl") diff --git a/test/slopes.jl b/test/slopes.jl new file mode 100644 index 00000000..2e26fc81 --- /dev/null +++ b/test/slopes.jl @@ -0,0 +1,122 @@ +using IntervalArithmetic, IntervalRootFinding +using ForwardDiff, StaticArrays +using Base.Test + +struct Slopes{T} + f::Function + x::Interval{T} + c::T + sol::Interval{T} +end + +@testset "Automatic slope expansion" begin + for T in [Float64, BigFloat] + s = interval(T(0.75), T(1.75)) + rts = Slopes{T}[] + push!(rts, Slopes(x->((x + sin(x)) * exp(-x^2)), s, mid(s), interval(T(-2.8), T(0.1)))) + push!(rts, Slopes(x->(x^4 - 10x^3 + 35x^2 - 50x + 24), s, mid(s), interval(T(-44), T(38.5)))) + push!(rts, Slopes(x->((log(x + 1.25) - 0.84x) ^ 2), s, mid(s), interval(T(-0.16), T(0.44)))) + push!(rts, Slopes(x->(0.02x^2 - 0.03exp(-(20(x - 0.875))^2)), s, mid(s), interval(T(0.03), T(0.33)))) + push!(rts, Slopes(x->(exp(x^2)), s, mid(s), interval(T(6.03), T(33.23)))) + push!(rts, Slopes(x->(x^4 - 12x^3 + 47x^2 - 60x - 20exp(-x)), s, mid(s), interval(T(-39), T(65.56)))) + push!(rts, Slopes(x->(x^6 - 15x^4 + 27x^2 + 250), s, mid(s), interval(T(-146.9), T(67.1)))) + push!(rts, Slopes(x->(atan(cos(tan(x)))), s, mid(s), interval(T(1), T(2)))) + push!(rts, Slopes(x->(asin(cos(acos(sin(x))))), s, mid(s), interval(T(1.36), T(∞)))) + + for i in 1:length(rts) + @test slope(rts[i].f, rts[i].x, rts[i].c) ⊆ rts[i].sol + end + end +end + +struct SlopesMulti + f::Function + x::IntervalBox + c::Vector + sol::Matrix{Interval} +end + +@testset "Multidim slope expansion" begin + + rts = SlopesMulti[] + f(x, y) = SVector(x^2 + y^2 - 1, y - 2x) + f(X) = f(X...) + X = (-6..6) × (-6..6) + c = [0.0, 0.0] + push!(rts, SlopesMulti(f, X, c, [-6..6 -6..6; -2.. -2 1..1])) + + function g(x) + (x1, x2, x3) = x + SVector( x1^2 + x2^2 + x3^2 - 1, + x1^2 + x3^2 - 0.25, + x1^2 + x2^2 - 4x3 + ) + end + + X = (-5..5) + XX = IntervalBox(X, 3) + cc = [0.0, 0.0, 0.0] + push!(rts, SlopesMulti(g, XX, cc, [-5..5 -5..5 -5..5; -5..5 0..0 -5..5; -5..5 -5..5 -4.. -4])) + function h(x) + (x1, x2, x3) = x + SVector( x1 + x2^2 + x3^2 - 1, + x1^2 + x3 - 0.25, + x1^2 + x2 - 4x3 + ) + end + + XXX = IntervalBox(1..5, 2..6, -3..7) + ccc = [3.0, 4.0, 2.0] + push!(rts, SlopesMulti(h, XXX, ccc, [1..1 6..10 -1..9; 4..8 0..0 1..1; 4..8 1..1 -4.. -4])) + + for i in 1:length(rts) + @test slope(rts[i].f, rts[i].x, rts[i].c) == rts[i].sol + end + +end + +struct SlopesMulti + f::Function + x::IntervalBox + c::Vector + sol::Matrix{Interval} +end + +@testset "Multidim slope expansion" begin + + rts = SlopesMulti[] + f(x, y) = SVector(x^2 + y^2 - 1, y - 2x) + f(X) = f(X...) + X = (-6..6) × (-6..6) + c = [0.0, 0.0] + push!(rts, SlopesMulti(f, X, c, [-6..6 -6..6; -2.. -2 1..1])) + + function g(x) + (x1, x2, x3) = x + SVector( x1^2 + x2^2 + x3^2 - 1, + x1^2 + x3^2 - 0.25, + x1^2 + x2^2 - 4x3 + ) + end + + X = (-5..5) + XX = IntervalBox(X, 3) + cc = [0.0, 0.0, 0.0] + push!(rts, SlopesMulti(g, XX, cc, [-5..5 -5..5 -5..5; -5..5 0..0 -5..5; -5..5 -5..5 -4.. -4])) + function h(x) + (x1, x2, x3) = x + SVector( x1 + x2^2 + x3^2 - 1, + x1^2 + x3 - 0.25, + x1^2 + x2 - 4x3 + ) + end + + XXX = IntervalBox(1..5, 2..6, -3..7) + ccc = [3.0, 4.0, 2.0] + push!(rts, SlopesMulti(h, XXX, ccc, [1..1 6..10 -1..9; 4..8 0..0 1..1; 4..8 1..1 -4.. -4])) + + for i in 1:length(rts) + @test slope(rts[i].f, rts[i].x, rts[i].c) == rts[i].sol + end + +end