From 79a828fd86b2f252b1555e16f927498ae3fde1ad Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Tue, 22 May 2018 16:03:59 +0530 Subject: [PATCH 01/36] Add automatic slope evaluation --- src/slopes.jl | 150 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 150 insertions(+) create mode 100644 src/slopes.jl diff --git a/src/slopes.jl b/src/slopes.jl new file mode 100644 index 00000000..62798766 --- /dev/null +++ b/src/slopes.jl @@ -0,0 +1,150 @@ +# Reference : Dietmar Ratz - An Optimized Interval Slope Arithmetic and its Application +using IntervalArithmetic +import Base: +, -, *, /, ^, sqrt, exp, log + +function slope(f::Function, x::Interval, c::Real) + f(SlopeVar(x, c)).fs +end + +struct SlopeType + fx::Interval + fc::Interval + fs::Interval +end + +function SlopeConst(c::Union{Real, Interval}) + SlopeType(c, c, 0) +end + +function SlopeVar(v::Real) + SlopeType(v, v, 1) +end + +function SlopeVar(v::Interval, c::Real) + SlopeType(v, c, 1) +end + +function fxValue(u::SlopeType) + u.fx +end + +function fcValue(u::SlopeType) + u.fc +end + +function fsValue(u::SlopeType) + u.fs +end + +function +(u::SlopeType, v::SlopeType) + SlopeType(u.fx + v.fx, u.fc + v.fc, u.fs + v.fs) +end + +function -(u::SlopeType, v::SlopeType) + SlopeType(u.fx - v.fx, u.fc - v.fc, u.fs - v.fs) +end + +function *(u::SlopeType, v::SlopeType) + SlopeType(u.fx * v.fx, u.fc * v.fc, u.fs * v.fc + u.fx * v.fs) +end + +function /(u::SlopeType, v::SlopeType) + SlopeType(u.fx / v.fx, u.fc / v.fc, (u.fs - (u.fc / v.fc) * v.fs) / v.fx) +end + +function +(u::Union{Interval, Real}, v::SlopeType) + SlopeType(u + v.fx, u + v.fc, v.fs) +end + +function -(u::Union{Interval, Real}, v::SlopeType) + SlopeType(u - v.fx, u - v.fc, -v.fs) +end + +function *(u::Union{Interval, Real}, v::SlopeType) + SlopeType(u * v.fx, u * v.fc, u * v.fs) +end + +function /(u::Union{Interval, Real}, v::SlopeType) + SlopeType(u / v.fx, u / v.fc, -(u / v.fc) * (v.fs / v.fx)) +end + ++(v::SlopeType, u::Union{Interval, Real}) = u + v + +-(v::SlopeType, u::Union{Interval, Real}) = u - v + +*(v::SlopeType, u::Union{Interval, Real}) = u * v + +/(v::SlopeType, u::Union{Interval, Real}) = u / v + +function sqr(u::SlopeType) + SlopeType(u.fx ^ 2, u.fc ^ 2, (u.fx + u.fc) * u.fs) +end + +function ^(u::SlopeType, k::Integer) + if k == 0 + return SlopeConst(1) + elseif k == 1 + return u + elseif k == 2 + return sqr(u) + else + hxi = interval(u.fx.lo) ^ k + hxs = interval(u.fx.hi) ^ k + hx = hull(hxi, hxs) + + if (k % 2 == 0) && (0 ∈ u.fx) + hx = interval(0, hx.hi) + end + + hc = u.fc ^ k + + i = u.fx.lo - u.fc.lo + s = u.fx.hi - u.fc.hi + + if ((i == 0) || (s == 0) || (k % 2 == 1 && Interval(0) ⪽ u.fx)) + h1 = k * (u.fx ^ (k - 1)) + else + if k % 2 == 0 || u.fx.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 SlopeType(hx, hc, h1 * u.fs) + end +end + +function sqrt(u::SlopeType) + SlopeType(sqrt(u.fx), sqrt(u.fc), u.fs / (sqrt(u.fx) + sqrt(u.fc))) +end + +function exp(u::SlopeType) + hx = exp(u.fx) + hc = exp(u.fc) + + i = u.fx.lo - u.fc.lo + s = u.fx.hi - u.fc.hi + + if (i == 0 || s == 0) + h1 = hx + else + h1 = interval((hx.lo - hc.lo) / i, (hx.hi - hc.hi) / s) + end + + SlopeType(hx, hc, h1 * u.fs) +end + +function log(u::SlopeType) + hx = log(u.fx) + hc = log(u.fc) + + i = u.fx.lo - u.fc.lo + s = u.fx.hi - u.fc.hi + + if (i == 0 || s == 0) + h1 = 1 / u.fx + else + h1 = interval((hx.hi - hc.hi) / s, (hx.lo - hc.lo) / i) + end + SlopeType(hx, hc, h1 * u.fs) +end From c9f03059549084b3ca88bcb3c35aabe3a8f7a11f Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Tue, 22 May 2018 16:32:18 +0530 Subject: [PATCH 02/36] Add functions which return derivative for completeness (outer bound) --- src/IntervalRootFinding.jl | 4 +-- src/slopes.jl | 55 +++++++++++++++++++++++++++++++++++--- 2 files changed, 54 insertions(+), 5 deletions(-) diff --git a/src/IntervalRootFinding.jl b/src/IntervalRootFinding.jl index 0496a574..e9ffc462 100644 --- a/src/IntervalRootFinding.jl +++ b/src/IntervalRootFinding.jl @@ -18,7 +18,7 @@ export derivative, jacobian, # reexport derivative from ForwardDiff Root, is_unique, roots, find_roots, - bisect, newton1d + bisect, newton1d, slope export isunique, root_status @@ -58,7 +58,7 @@ include("complex.jl") include("contractors.jl") include("roots.jl") include("newton1d.jl") - +include("slopes.jl") gradient(f) = X -> ForwardDiff.gradient(f, SVector(X)) diff --git a/src/slopes.jl b/src/slopes.jl index 62798766..627011dc 100644 --- a/src/slopes.jl +++ b/src/slopes.jl @@ -1,9 +1,15 @@ # Reference : Dietmar Ratz - An Optimized Interval Slope Arithmetic and its Application -using IntervalArithmetic -import Base: +, -, *, /, ^, sqrt, exp, log +using IntervalArithmetic, ForwardDiff +import Base: +, -, *, /, ^, sqrt, exp, log, sin, cos, tan, asin, acos, atan function slope(f::Function, x::Interval, c::Real) - f(SlopeVar(x, c)).fs + try + f(SlopeVar(x, c)).fs + catch y + if isa(y, MethodError) + ForwardDiff.derivative(f, x) + end + end end struct SlopeType @@ -71,6 +77,7 @@ end +(v::SlopeType, u::Union{Interval, Real}) = u + v -(v::SlopeType, u::Union{Interval, Real}) = u - v +-(u::SlopeType) = u * -1 *(v::SlopeType, u::Union{Interval, Real}) = u * v @@ -148,3 +155,45 @@ function log(u::SlopeType) end SlopeType(hx, hc, h1 * u.fs) end + +function sin(u::SlopeType) # Using derivative to upper bound the slope expansion for now + hx = sin(u.fx) + hc = sin(u.fc) + hs = cos(u.fx) + SlopeType(hx, hc, hs) +end + +function cos(u::SlopeType) # Using derivative to upper bound the slope expansion for now + hx = cos(u.fx) + hc = cos(u.fc) + hs = -sin(u.fx) + SlopeType(hx, hc, hs) +end + +function tan(u::SlopeType) # Using derivative to upper bound the slope expansion for now + hx = tan(u.fx) + hc = tan(u.fc) + hs = (sec(u.fx)) ^ 2 + SlopeType(hx, hc, hs) +end + +function asin(u::SlopeType) + hx = asin(u.fx) + hc = asin(u.fc) + hs = 1 / sqrt(1 - (u.fx ^ 2)) + SlopeType(hx, hc, hs) +end + +function acos(u::SlopeType) + hx = acos(u.fx) + hc = acos(u.fc) + hs = -1 / sqrt(1 - (u.fx ^ 2)) + SlopeType(hx, hc, hs) +end + +function atan(u::SlopeType) + hx = atan(u.fx) + hc = atan(u.fc) + hs = 1 / 1 + (u.fx ^ 2) + SlopeType(hx, hc, hs) +end From 0c76318bc506e917167504aafdc36ef18e716eed Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Tue, 22 May 2018 21:08:03 +0530 Subject: [PATCH 03/36] Add benchmarking and data --- perf/slopes.jl | 57 +++++++++++++++++++++++++++++++++++++++++ perf/slopes_tabular.txt | 23 +++++++++++++++++ 2 files changed, 80 insertions(+) create mode 100644 perf/slopes.jl create mode 100644 perf/slopes_tabular.txt diff --git a/perf/slopes.jl b/perf/slopes.jl new file mode 100644 index 00000000..0db3e0aa --- /dev/null +++ b/perf/slopes.jl @@ -0,0 +1,57 @@ +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)) + df = DataFrame() + df[:Method] = ["Automatic Differentiation", "Slope Expansion"] + for n in 1:length(f) + + t1 = ForwardDiff.derivative(f[n], 0.75..1.75) + t2 = slope(f[n], 0.75..1.75, 1.25) + 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)) + df = DataFrame() + df[:Method] = ["Automatic Differentiation", "Slope Expansion"] + for n in 1:length(f) + + t1 = @belapsed ForwardDiff.derivative($f[$n], 0.75..1.75) + t2 = @belapsed slope($f[$n], 0.75..1.75, 1.25) + 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 new file mode 100644 index 00000000..defdf488 --- /dev/null +++ b/perf/slopes_tabular.txt @@ -0,0 +1,23 @@ +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] │ + +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 │ From f837cafaa498a4eefe9c5a382da239c869c9f5e2 Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Wed, 23 May 2018 19:28:47 +0530 Subject: [PATCH 04/36] Add tests --- test/runtests.jl | 1 + test/slopes.jl | 25 +++++++++++++++++++++++++ 2 files changed, 26 insertions(+) create mode 100644 test/slopes.jl diff --git a/test/runtests.jl b/test/runtests.jl index e905e897..62e2d775 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,3 +5,4 @@ using Base.Test #include("findroots.jl") include("roots.jl") include("newton1d.jl") +include("slopes.jl") diff --git a/test/slopes.jl b/test/slopes.jl new file mode 100644 index 00000000..08f571ea --- /dev/null +++ b/test/slopes.jl @@ -0,0 +1,25 @@ +using IntervalArithmetic, IntervalRootFinding +using ForwardDiff +using Base.Test + +struct Slopes{T} + f::Function + x::Interval{T} + c::T + sol::Interval{T} +end + +@testset "Automatic slope expansion" begin + rts = Slopes{Float64}[] + push!(rts, Slopes(x->((x + sin(x)) * exp(-x^2)), interval(0.75, 1.75), 1.25, interval(-2.8, 0.1))) + push!(rts, Slopes(x->(x^4 - 10x^3 + 35x^2 - 50x + 24), interval(0.75, 1.75), 1.25, interval(-44, 38.5))) + push!(rts, Slopes(x->((log(x + 1.25) - 0.84x) ^ 2), interval(0.75, 1.75), 1.25, interval(-0.16, 0.44))) + push!(rts, Slopes(x->(0.02x^2 - 0.03exp(-(20(x - 0.875))^2)), interval(0.75, 1.75), 1.25, interval(0.04, 0.33))) + push!(rts, Slopes(x->(exp(x^2)), interval(0.75, 1.75), 1.25, interval(6.03, 33.23))) + push!(rts, Slopes(x->(x^4 - 12x^3 + 47x^2 - 60x - 20exp(-x)), interval(0.75, 1.75), 1.25, interval(-39, 65.56))) + push!(rts, Slopes(x->(x^6 - 15x^4 + 27x^2 + 250), interval(0.75, 1.75), 1.25, interval(-146.9, 67.1))) + + for i in 1:length(rts) + @test slope(rts[i].f, rts[i].x, rts[i].c) ⊆ rts[i].sol + end +end From 7b67107c77009c1ee0b7c56d5ace96fe2876447d Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Thu, 24 May 2018 14:36:22 +0530 Subject: [PATCH 05/36] Incorporate review comments --- perf/slopes.jl | 7 +- src/slopes.jl | 197 ++++++++++++++++++++++++------------------------- test/slopes.jl | 16 ++-- 3 files changed, 112 insertions(+), 108 deletions(-) diff --git a/perf/slopes.jl b/perf/slopes.jl index 0db3e0aa..6340639f 100644 --- a/perf/slopes.jl +++ b/perf/slopes.jl @@ -9,12 +9,15 @@ function benchmark_results() 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)) + + 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], 0.75..1.75) - t2 = slope(f[n], 0.75..1.75, 1.25) + t1 = ForwardDiff.derivative(f[n], s) + t2 = slope(f[n], s, mid(s)) df[Symbol("f" * "$n")] = [t1, t2] end a = [] diff --git a/src/slopes.jl b/src/slopes.jl index 627011dc..cef00280 100644 --- a/src/slopes.jl +++ b/src/slopes.jl @@ -1,10 +1,11 @@ # Reference : Dietmar Ratz - An Optimized Interval Slope Arithmetic and its Application using IntervalArithmetic, ForwardDiff import Base: +, -, *, /, ^, sqrt, exp, log, sin, cos, tan, asin, acos, atan +import IntervalArithmetic: mid, interval function slope(f::Function, x::Interval, c::Real) try - f(SlopeVar(x, c)).fs + f(slope_var(x, c)).s catch y if isa(y, MethodError) ForwardDiff.derivative(f, x) @@ -12,125 +13,123 @@ function slope(f::Function, x::Interval, c::Real) end end -struct SlopeType - fx::Interval - fc::Interval - fs::Interval +struct Slope{T} + x::Interval{T} + c::Interval{T} + s::Interval{T} + Slope{T}(a, b, c) where T = new(a, b, c) + Slope{T}(c) where T = Slope{T}(c, c, 0) end -function SlopeConst(c::Union{Real, Interval}) - SlopeType(c, c, 0) +function slope_var(v::Real) + Slope{Float64}(v, v, 1) end -function SlopeVar(v::Real) - SlopeType(v, v, 1) +function slope_var(v::Interval, c::Real) + Slope{Float64}(v, c, 1) end -function SlopeVar(v::Interval, c::Real) - SlopeType(v, c, 1) +function interval(u::Slope) + u.x end -function fxValue(u::SlopeType) - u.fx +function mid(u::Slope) + u.c end -function fcValue(u::SlopeType) - u.fc +function slope(u::Slope) + u.s end -function fsValue(u::SlopeType) - u.fs +function +(u::Slope, v::Slope) + Slope{Float64}(u.x + v.x, u.c + v.c, u.s + v.s) end -function +(u::SlopeType, v::SlopeType) - SlopeType(u.fx + v.fx, u.fc + v.fc, u.fs + v.fs) +function -(u::Slope, v::Slope) + Slope{Float64}(u.x - v.x, u.c - v.c, u.s - v.s) end -function -(u::SlopeType, v::SlopeType) - SlopeType(u.fx - v.fx, u.fc - v.fc, u.fs - v.fs) +function *(u::Slope, v::Slope) + Slope{Float64}(u.x * v.x, u.c * v.c, u.s * v.c + u.x * v.s) end -function *(u::SlopeType, v::SlopeType) - SlopeType(u.fx * v.fx, u.fc * v.fc, u.fs * v.fc + u.fx * v.fs) +function /(u::Slope, v::Slope) + Slope{Float64}(u.x / v.x, u.c / v.c, (u.s - (u.c / v.c) * v.s) / v.x) end -function /(u::SlopeType, v::SlopeType) - SlopeType(u.fx / v.fx, u.fc / v.fc, (u.fs - (u.fc / v.fc) * v.fs) / v.fx) +function +(u::Union{Interval, Real}, v::Slope) + Slope{Float64}(u + v.x, u + v.c, v.s) end -function +(u::Union{Interval, Real}, v::SlopeType) - SlopeType(u + v.fx, u + v.fc, v.fs) +function -(u::Union{Interval, Real}, v::Slope) + Slope{Float64}(u - v.x, u - v.c, -v.s) end -function -(u::Union{Interval, Real}, v::SlopeType) - SlopeType(u - v.fx, u - v.fc, -v.fs) +function *(u::Union{Interval, Real}, v::Slope) + Slope{Float64}(u * v.x, u * v.c, u * v.s) end -function *(u::Union{Interval, Real}, v::SlopeType) - SlopeType(u * v.fx, u * v.fc, u * v.fs) +function /(u::Union{Interval, Real}, v::Slope) + Slope{Float64}(u / v.x, u / v.c, -(u / v.c) * (v.s / v.x)) end -function /(u::Union{Interval, Real}, v::SlopeType) - SlopeType(u / v.fx, u / v.fc, -(u / v.fc) * (v.fs / v.fx)) -end - -+(v::SlopeType, u::Union{Interval, Real}) = u + v ++(v::Slope, u::Union{Interval, Real}) = u + v --(v::SlopeType, u::Union{Interval, Real}) = u - v --(u::SlopeType) = u * -1 +-(v::Slope, u::Union{Interval, Real}) = u - v +-(u::Slope) = u * -1 -*(v::SlopeType, u::Union{Interval, Real}) = u * v +*(v::Slope, u::Union{Interval, Real}) = u * v -/(v::SlopeType, u::Union{Interval, Real}) = u / v +/(v::Slope, u::Union{Interval, Real}) = u / v -function sqr(u::SlopeType) - SlopeType(u.fx ^ 2, u.fc ^ 2, (u.fx + u.fc) * u.fs) +function sqr(u::Slope) + Slope{Float64}(u.x ^ 2, u.c ^ 2, (u.x + u.c) * u.s) end -function ^(u::SlopeType, k::Integer) +function ^(u::Slope, k::Integer) if k == 0 - return SlopeConst(1) + return Slope{Float64}(1) elseif k == 1 return u elseif k == 2 return sqr(u) else - hxi = interval(u.fx.lo) ^ k - hxs = interval(u.fx.hi) ^ k + hxi = interval(u.x.lo) ^ k + hxs = interval(u.x.hi) ^ k hx = hull(hxi, hxs) - if (k % 2 == 0) && (0 ∈ u.fx) + if (k % 2 == 0) && (0 ∈ u.x) hx = interval(0, hx.hi) end - hc = u.fc ^ k + hc = u.c ^ k - i = u.fx.lo - u.fc.lo - s = u.fx.hi - u.fc.hi + i = u.x.lo - u.c.lo + s = u.x.hi - u.c.hi - if ((i == 0) || (s == 0) || (k % 2 == 1 && Interval(0) ⪽ u.fx)) - h1 = k * (u.fx ^ (k - 1)) + if ((i == 0) || (s == 0) || (k % 2 == 1 && Interval(0) ⪽ u.x)) + h1 = k * (u.x ^ (k - 1)) else - if k % 2 == 0 || u.fx.lo >= 0 + 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 SlopeType(hx, hc, h1 * u.fs) + return Slope{Float64}(hx, hc, h1 * u.s) end end -function sqrt(u::SlopeType) - SlopeType(sqrt(u.fx), sqrt(u.fc), u.fs / (sqrt(u.fx) + sqrt(u.fc))) +function sqrt(u::Slope) + Slope{Float64}(sqrt(u.x), sqrt(u.c), u.s / (sqrt(u.x) + sqrt(u.c))) end -function exp(u::SlopeType) - hx = exp(u.fx) - hc = exp(u.fc) +function exp(u::Slope) + hx = exp(u.x) + hc = exp(u.c) - i = u.fx.lo - u.fc.lo - s = u.fx.hi - u.fc.hi + i = u.x.lo - u.c.lo + s = u.x.hi - u.c.hi if (i == 0 || s == 0) h1 = hx @@ -138,62 +137,62 @@ function exp(u::SlopeType) h1 = interval((hx.lo - hc.lo) / i, (hx.hi - hc.hi) / s) end - SlopeType(hx, hc, h1 * u.fs) + Slope{Float64}(hx, hc, h1 * u.s) end -function log(u::SlopeType) - hx = log(u.fx) - hc = log(u.fc) +function log(u::Slope) + hx = log(u.x) + hc = log(u.c) - i = u.fx.lo - u.fc.lo - s = u.fx.hi - u.fc.hi + i = u.x.lo - u.c.lo + s = u.x.hi - u.c.hi if (i == 0 || s == 0) - h1 = 1 / u.fx + h1 = 1 / u.x else h1 = interval((hx.hi - hc.hi) / s, (hx.lo - hc.lo) / i) end - SlopeType(hx, hc, h1 * u.fs) + Slope{Float64}(hx, hc, h1 * u.s) end -function sin(u::SlopeType) # Using derivative to upper bound the slope expansion for now - hx = sin(u.fx) - hc = sin(u.fc) - hs = cos(u.fx) - SlopeType(hx, hc, hs) +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{Float64}(hx, hc, hs) end -function cos(u::SlopeType) # Using derivative to upper bound the slope expansion for now - hx = cos(u.fx) - hc = cos(u.fc) - hs = -sin(u.fx) - SlopeType(hx, hc, hs) +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{Float64}(hx, hc, hs) end -function tan(u::SlopeType) # Using derivative to upper bound the slope expansion for now - hx = tan(u.fx) - hc = tan(u.fc) - hs = (sec(u.fx)) ^ 2 - SlopeType(hx, hc, hs) +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{Float64}(hx, hc, hs) end -function asin(u::SlopeType) - hx = asin(u.fx) - hc = asin(u.fc) - hs = 1 / sqrt(1 - (u.fx ^ 2)) - SlopeType(hx, hc, hs) +function asin(u::Slope) + hx = asin(u.x) + hc = asin(u.c) + hs = 1 / sqrt(1 - (u.x ^ 2)) + Slope{Float64}(hx, hc, hs) end -function acos(u::SlopeType) - hx = acos(u.fx) - hc = acos(u.fc) - hs = -1 / sqrt(1 - (u.fx ^ 2)) - SlopeType(hx, hc, hs) +function acos(u::Slope) + hx = acos(u.x) + hc = acos(u.c) + hs = -1 / sqrt(1 - (u.x ^ 2)) + Slope{Float64}(hx, hc, hs) end -function atan(u::SlopeType) - hx = atan(u.fx) - hc = atan(u.fc) - hs = 1 / 1 + (u.fx ^ 2) - SlopeType(hx, hc, hs) +function atan(u::Slope) + hx = atan(u.x) + hc = atan(u.c) + hs = 1 / 1 + (u.x ^ 2) + Slope{Float64}(hx, hc, hs) end diff --git a/test/slopes.jl b/test/slopes.jl index 08f571ea..7efd1018 100644 --- a/test/slopes.jl +++ b/test/slopes.jl @@ -10,14 +10,16 @@ struct Slopes{T} end @testset "Automatic slope expansion" begin + + s = interval(0.75, 1.75) rts = Slopes{Float64}[] - push!(rts, Slopes(x->((x + sin(x)) * exp(-x^2)), interval(0.75, 1.75), 1.25, interval(-2.8, 0.1))) - push!(rts, Slopes(x->(x^4 - 10x^3 + 35x^2 - 50x + 24), interval(0.75, 1.75), 1.25, interval(-44, 38.5))) - push!(rts, Slopes(x->((log(x + 1.25) - 0.84x) ^ 2), interval(0.75, 1.75), 1.25, interval(-0.16, 0.44))) - push!(rts, Slopes(x->(0.02x^2 - 0.03exp(-(20(x - 0.875))^2)), interval(0.75, 1.75), 1.25, interval(0.04, 0.33))) - push!(rts, Slopes(x->(exp(x^2)), interval(0.75, 1.75), 1.25, interval(6.03, 33.23))) - push!(rts, Slopes(x->(x^4 - 12x^3 + 47x^2 - 60x - 20exp(-x)), interval(0.75, 1.75), 1.25, interval(-39, 65.56))) - push!(rts, Slopes(x->(x^6 - 15x^4 + 27x^2 + 250), interval(0.75, 1.75), 1.25, interval(-146.9, 67.1))) + push!(rts, Slopes(x->((x + sin(x)) * exp(-x^2)), s, mid(s), interval(-2.8, 0.1))) + push!(rts, Slopes(x->(x^4 - 10x^3 + 35x^2 - 50x + 24), s, mid(s), interval(-44, 38.5))) + push!(rts, Slopes(x->((log(x + 1.25) - 0.84x) ^ 2), s, mid(s), interval(-0.16, 0.44))) + push!(rts, Slopes(x->(0.02x^2 - 0.03exp(-(20(x - 0.875))^2)), s, mid(s), interval(0.04, 0.33))) + push!(rts, Slopes(x->(exp(x^2)), s, mid(s), interval(6.03, 33.23))) + push!(rts, Slopes(x->(x^4 - 12x^3 + 47x^2 - 60x - 20exp(-x)), s, mid(s), interval(-39, 65.56))) + push!(rts, Slopes(x->(x^6 - 15x^4 + 27x^2 + 250), s, mid(s), interval(-146.9, 67.1))) for i in 1:length(rts) @test slope(rts[i].f, rts[i].x, rts[i].c) ⊆ rts[i].sol From 3cd778a69839d6c185f5b7b6d0340ece8cdfd421 Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Fri, 25 May 2018 15:53:37 +0530 Subject: [PATCH 06/36] Use concrete types and add BigFloat tests --- perf/slopes.jl | 10 ++++-- perf/slopes_tabular.txt | 16 ++++++++- src/slopes.jl | 75 +++++++++++++++++++++-------------------- test/slopes.jl | 24 +++++++++++-- 4 files changed, 83 insertions(+), 42 deletions(-) diff --git a/perf/slopes.jl b/perf/slopes.jl index 6340639f..ea3e949c 100644 --- a/perf/slopes.jl +++ b/perf/slopes.jl @@ -9,6 +9,8 @@ function benchmark_results() 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() @@ -40,12 +42,16 @@ function benchmark_time() 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], 0.75..1.75) - t2 = @belapsed slope($f[$n], 0.75..1.75, 1.25) + t1 = @belapsed ForwardDiff.derivative($f[$n], $s) + t2 = @belapsed slope($f[$n], $s, mid($s)) df[Symbol("f" * "$n")] = [t1, t2] end a = [] diff --git a/perf/slopes_tabular.txt b/perf/slopes_tabular.txt index defdf488..2d744244 100644 --- a/perf/slopes_tabular.txt +++ b/perf/slopes_tabular.txt @@ -9,9 +9,11 @@ Results │ 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 │ @@ -21,3 +23,15 @@ Time │ 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/slopes.jl b/src/slopes.jl index cef00280..96df9f4e 100644 --- a/src/slopes.jl +++ b/src/slopes.jl @@ -3,7 +3,7 @@ using IntervalArithmetic, ForwardDiff import Base: +, -, *, /, ^, sqrt, exp, log, sin, cos, tan, asin, acos, atan import IntervalArithmetic: mid, interval -function slope(f::Function, x::Interval, c::Real) +function slope(f::Function, x::Interval, c::Number) try f(slope_var(x, c)).s catch y @@ -17,16 +17,17 @@ struct Slope{T} x::Interval{T} c::Interval{T} s::Interval{T} - Slope{T}(a, b, c) where T = new(a, b, c) - Slope{T}(c) where T = Slope{T}(c, c, 0) end -function slope_var(v::Real) - Slope{Float64}(v, v, 1) +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::Real) - Slope{Float64}(v, c, 1) +function slope_var(v::Interval, c::Number) + Slope(v, c, 1) end function interval(u::Slope) @@ -42,53 +43,53 @@ function slope(u::Slope) end function +(u::Slope, v::Slope) - Slope{Float64}(u.x + v.x, u.c + v.c, u.s + v.s) + Slope(u.x + v.x, u.c + v.c, u.s + v.s) end function -(u::Slope, v::Slope) - Slope{Float64}(u.x - v.x, u.c - v.c, u.s - v.s) + Slope(u.x - v.x, u.c - v.c, u.s - v.s) end function *(u::Slope, v::Slope) - Slope{Float64}(u.x * v.x, u.c * v.c, u.s * v.c + u.x * v.s) + Slope(u.x * v.x, u.c * v.c, u.s * v.c + u.x * v.s) end function /(u::Slope, v::Slope) - Slope{Float64}(u.x / v.x, u.c / v.c, (u.s - (u.c / v.c) * v.s) / v.x) + Slope(u.x / v.x, u.c / v.c, (u.s - (u.c / v.c) * v.s) / v.x) end -function +(u::Union{Interval, Real}, v::Slope) - Slope{Float64}(u + v.x, u + v.c, v.s) +function +(u, v::Slope) + Slope(u + v.x, u + v.c, v.s) end -function -(u::Union{Interval, Real}, v::Slope) - Slope{Float64}(u - v.x, u - v.c, -v.s) +function -(u, v::Slope) + Slope(u - v.x, u - v.c, -v.s) end -function *(u::Union{Interval, Real}, v::Slope) - Slope{Float64}(u * v.x, u * v.c, u * v.s) +function *(u, v::Slope) + Slope(u * v.x, u * v.c, u * v.s) end -function /(u::Union{Interval, Real}, v::Slope) - Slope{Float64}(u / v.x, u / v.c, -(u / v.c) * (v.s / v.x)) +function /(u, v::Slope) + Slope(u / v.x, u / v.c, -(u / v.c) * (v.s / v.x)) end -+(v::Slope, u::Union{Interval, Real}) = u + v ++(v::Slope, u) = u + v --(v::Slope, u::Union{Interval, Real}) = u - v --(u::Slope) = u * -1 +-(v::Slope, u) = u - v +-(u::Slope) = u * -1.0 -*(v::Slope, u::Union{Interval, Real}) = u * v +*(v::Slope, u) = u * v -/(v::Slope, u::Union{Interval, Real}) = u / v +/(v::Slope, u) = u / v function sqr(u::Slope) - Slope{Float64}(u.x ^ 2, u.c ^ 2, (u.x + u.c) * u.s) + Slope(u.x ^ 2, u.c ^ 2, (u.x + u.c) * u.s) end function ^(u::Slope, k::Integer) if k == 0 - return Slope{Float64}(1) + return Slope(1) elseif k == 1 return u elseif k == 2 @@ -107,7 +108,7 @@ function ^(u::Slope, k::Integer) i = u.x.lo - u.c.lo s = u.x.hi - u.c.hi - if ((i == 0) || (s == 0) || (k % 2 == 1 && Interval(0) ⪽ u.x)) + 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 @@ -116,12 +117,12 @@ function ^(u::Slope, k::Integer) h1 = interval((hxs.lo - hc.hi) / s, (hxi.lo - hc.hi) / i) end end - return Slope{Float64}(hx, hc, h1 * u.s) + return Slope(hx, hc, h1 * u.s) end end function sqrt(u::Slope) - Slope{Float64}(sqrt(u.x), sqrt(u.c), u.s / (sqrt(u.x) + sqrt(u.c))) + Slope(sqrt(u.x), sqrt(u.c), u.s / (sqrt(u.x) + sqrt(u.c))) end function exp(u::Slope) @@ -137,7 +138,7 @@ function exp(u::Slope) h1 = interval((hx.lo - hc.lo) / i, (hx.hi - hc.hi) / s) end - Slope{Float64}(hx, hc, h1 * u.s) + Slope(hx, hc, h1 * u.s) end function log(u::Slope) @@ -152,47 +153,47 @@ function log(u::Slope) else h1 = interval((hx.hi - hc.hi) / s, (hx.lo - hc.lo) / i) end - Slope{Float64}(hx, hc, h1 * u.s) + 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{Float64}(hx, hc, hs) + 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{Float64}(hx, hc, hs) + 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{Float64}(hx, hc, hs) + 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{Float64}(hx, hc, hs) + 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{Float64}(hx, hc, hs) + Slope(hx, hc, hs) end function atan(u::Slope) hx = atan(u.x) hc = atan(u.c) hs = 1 / 1 + (u.x ^ 2) - Slope{Float64}(hx, hc, hs) + Slope(hx, hc, hs) end diff --git a/test/slopes.jl b/test/slopes.jl index 7efd1018..6730583f 100644 --- a/test/slopes.jl +++ b/test/slopes.jl @@ -9,17 +9,37 @@ struct Slopes{T} sol::Interval{T} end -@testset "Automatic slope expansion" begin +@testset "Automatic slope expansion(Float64)" begin s = interval(0.75, 1.75) rts = Slopes{Float64}[] push!(rts, Slopes(x->((x + sin(x)) * exp(-x^2)), s, mid(s), interval(-2.8, 0.1))) push!(rts, Slopes(x->(x^4 - 10x^3 + 35x^2 - 50x + 24), s, mid(s), interval(-44, 38.5))) push!(rts, Slopes(x->((log(x + 1.25) - 0.84x) ^ 2), s, mid(s), interval(-0.16, 0.44))) - push!(rts, Slopes(x->(0.02x^2 - 0.03exp(-(20(x - 0.875))^2)), s, mid(s), interval(0.04, 0.33))) + push!(rts, Slopes(x->(0.02x^2 - 0.03exp(-(20(x - 0.875))^2)), s, mid(s), interval(0.03, 0.33))) push!(rts, Slopes(x->(exp(x^2)), s, mid(s), interval(6.03, 33.23))) push!(rts, Slopes(x->(x^4 - 12x^3 + 47x^2 - 60x - 20exp(-x)), s, mid(s), interval(-39, 65.56))) push!(rts, Slopes(x->(x^6 - 15x^4 + 27x^2 + 250), s, mid(s), interval(-146.9, 67.1))) + push!(rts, Slopes(x->(atan(cos(tan(x)))), s, mid(s), interval(1, 2))) + push!(rts, Slopes(x->(asin(cos(acos(sin(x))))), s, mid(s), interval(1.36, ∞))) + + for i in 1:length(rts) + @test slope(rts[i].f, rts[i].x, rts[i].c) ⊆ rts[i].sol + end +end + +@testset "Automatic slope expansion(BigFloat)" begin + s = interval(BigFloat(0.75), BigFloat(1.75)) + rts = Slopes{BigFloat}[] + push!(rts, Slopes(x->((x + sin(x)) * exp(-x^2)), s, mid(s), interval(BigFloat(-2.8), BigFloat(0.1)))) + push!(rts, Slopes(x->(x^4 - 10x^3 + 35x^2 - 50x + 24), s, mid(s), interval(BigFloat(-44), BigFloat(38.5)))) + push!(rts, Slopes(x->((log(x + 1.25) - 0.84x) ^ 2), s, mid(s), interval(BigFloat(-0.16), BigFloat(0.44)))) + push!(rts, Slopes(x->(0.02x^2 - 0.03exp(-(20(x - 0.875))^2)), s, mid(s), interval(BigFloat(0.03), BigFloat(0.33)))) + push!(rts, Slopes(x->(exp(x^2)), s, mid(s), interval(BigFloat(6.03), BigFloat(33.23)))) + push!(rts, Slopes(x->(x^4 - 12x^3 + 47x^2 - 60x - 20exp(-x)), s, mid(s), interval(BigFloat(-39), BigFloat(65.56)))) + push!(rts, Slopes(x->(x^6 - 15x^4 + 27x^2 + 250), s, mid(s), interval(BigFloat(-146.9), BigFloat(67.1)))) + push!(rts, Slopes(x->(atan(cos(tan(x)))), s, mid(s), interval(BigFloat(1), BigFloat(2)))) + push!(rts, Slopes(x->(asin(cos(acos(sin(x))))), s, mid(s), interval(BigFloat(1.36), BigFloat(∞)))) for i in 1:length(rts) @test slope(rts[i].f, rts[i].x, rts[i].c) ⊆ rts[i].sol From 8500a9914e69b6d47287d762e98d895c9104f0ef Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Sun, 27 May 2018 16:14:24 +0530 Subject: [PATCH 07/36] Add slope function for multidim case with tests and benchmarks --- perf/slopes.jl | 50 +++++++++++++++++++++++++++++++++++++++++ perf/slopes_tabular.txt | 45 +++++++++++++++++++++++++++++++++++++ src/slopes.jl | 41 ++++++++++++++++++++++++++++----- test/slopes.jl | 48 ++++++++++++++++++++++++++++++++++++++- 4 files changed, 177 insertions(+), 7 deletions(-) diff --git a/perf/slopes.jl b/perf/slopes.jl index ea3e949c..642a8d84 100644 --- a/perf/slopes.jl +++ b/perf/slopes.jl @@ -64,3 +64,53 @@ function benchmark_time() 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 index 2d744244..c8ac3c28 100644 --- a/perf/slopes_tabular.txt +++ b/perf/slopes_tabular.txt @@ -35,3 +35,48 @@ Time │ 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/slopes.jl b/src/slopes.jl index 96df9f4e..28345edc 100644 --- a/src/slopes.jl +++ b/src/slopes.jl @@ -1,9 +1,12 @@ # Reference : Dietmar Ratz - An Optimized Interval Slope Arithmetic and its Application -using IntervalArithmetic, ForwardDiff +using IntervalArithmetic, ForwardDiff, StaticArrays import Base: +, -, *, /, ^, sqrt, exp, log, sin, cos, tan, asin, acos, atan import IntervalArithmetic: mid, interval -function slope(f::Function, x::Interval, c::Number) +""" +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 @@ -13,6 +16,25 @@ function slope(f::Function, x::Interval, c::Number) 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} @@ -76,12 +98,19 @@ end +(v::Slope, u) = u + v --(v::Slope, u) = u - v --(u::Slope) = u * -1.0 - *(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) diff --git a/test/slopes.jl b/test/slopes.jl index 6730583f..d4e54f02 100644 --- a/test/slopes.jl +++ b/test/slopes.jl @@ -1,5 +1,5 @@ using IntervalArithmetic, IntervalRootFinding -using ForwardDiff +using ForwardDiff, StaticArrays using Base.Test struct Slopes{T} @@ -45,3 +45,49 @@ end @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 From 6257f207126773d2e1db55a6b1ce7bbae6cfa871 Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Sun, 1 Apr 2018 14:14:06 +0530 Subject: [PATCH 08/36] Add debug mode and bug fixes --- src/newton1d.jl | 33 ++++++++++++++++++++++++++++++--- 1 file changed, 30 insertions(+), 3 deletions(-) diff --git a/src/newton1d.jl b/src/newton1d.jl index cac4be47..774cc264 100644 --- a/src/newton1d.jl +++ b/src/newton1d.jl @@ -5,7 +5,7 @@ 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=10eps(T), abstol=eps(T), debug=false) L = Interval{T}[] @@ -15,15 +15,24 @@ function newton1d{T}(f::Function, f′::Function, x::Interval{T}; while !isempty(L) X = pop!(L) + + debug && (print("Current interval popped: "); @show X) + m = mid(X) if (isempty(X)) continue end if 0 ∉ f′(X) + + debug && println("0 ∉ f′(X)") + while true m = mid(X) N = m - (f(Interval(m)) / f′(X)) + + debug && (print("Newton step: "); @show (X, X ∩ N)) + X = X ∩ N if isempty(X) @@ -31,18 +40,25 @@ function newton1d{T}(f::Function, f′::Function, x::Interval{T}; elseif 0 ∈ f(Interval(prevfloat(m), nextfloat(m))) push!(R, Root(X, :unique)) + + debug && @show "Root found", X + break end end else # 0 ∈ f'(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))) # 0 ∈ fⁱ(x) - # Step 7 + + debug && println("0 ∈ fⁱ(x)") if 0 ∉ f(Interval(X.lo)) expansion_pt = X.lo @@ -59,7 +75,10 @@ function newton1d{T}(f::Function, f′::Function, x::Interval{T}; continue else - push!(R, Root(X, :unique)) + push!(R, Root(X, :unknown)) + + debug && @show "Multiple root found", X + continue end end @@ -67,8 +86,13 @@ function newton1d{T}(f::Function, f′::Function, x::Interval{T}; else # 0 ∉ fⁱ(x) + debug && println("0 ∉ fⁱ(x)") + if (diam(X)/mag(X)) < reltol && diam(f(X)) < abstol push!(R, Root(X, :unknown)) + + debug && @show "Tolerance root found", X + continue end end @@ -98,6 +122,9 @@ function newton1d{T}(f::Function, f′::Function, x::Interval{T}; else N = expansion_pt - (f(Interval(expansion_pt))/f′(X)) + + debug && (print("Newton step: "); @show (X, X ∩ N)) + X = X ∩ N m = mid(X) From aa50d9183fed3cf016ecbbad18ee2d38b951fcd8 Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Sun, 1 Apr 2018 15:39:38 +0530 Subject: [PATCH 09/36] Add tests --- test/newton1d.jl | 19 +++++++++++++++++-- 1 file changed, 17 insertions(+), 2 deletions(-) diff --git a/test/newton1d.jl b/test/newton1d.jl index 74fadcf7..095f7805 100644 --- a/test/newton1d.jl +++ b/test/newton1d.jl @@ -22,17 +22,24 @@ three_halves_pi = 3*big_pi/2 f′(x) = 2*x*exp(x^2) + sin(x) f1(x) = x^4 - 10x^3 + 35x^2 - 50x + 24 f1′(x) = 4x^3 - 30x^2 + 70x - 50 - + f2(x) = 4567x^2 - 9134x + 4567 + f2′(x) = 9134x - 9134 + f3(x) = (x^2 - 2)^2 + f3′(x) = 4x * (x^2 - 2) 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) 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) end @test length(rts1) == 3 @@ -42,12 +49,20 @@ three_halves_pi = 3*big_pi/2 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), sqrt(2), sqrt(2)] + for i = 1:length(rts5) + @test L1[i] in rts5[i].interval && :unknown == rts5[i].status + end end end From 1412b1f040892a4c96046e15a48d65d9f4e30fc2 Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Tue, 29 May 2018 00:11:00 +0530 Subject: [PATCH 10/36] Major robustness improvements and hard test cases added --- src/newton1d.jl | 123 +++++++++++++++++++++++++++++++++-------------- test/newton1d.jl | 30 +++++++++--- 2 files changed, 111 insertions(+), 42 deletions(-) diff --git a/src/newton1d.jl b/src/newton1d.jl index 774cc264..bffb89b1 100644 --- a/src/newton1d.jl +++ b/src/newton1d.jl @@ -5,14 +5,16 @@ 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=10eps(T), abstol=eps(T), debug=false) + reltol=eps(T), abstol=eps(T), debug=false, debugroot=false) L = Interval{T}[] R = Root{Interval{T}}[] + reps = reps1 = 0 push!(L, x) - + initial_width = ∞ + X = emptyinterval(T) while !isempty(L) X = pop!(L) @@ -28,56 +30,111 @@ function newton1d{T}(f::Function, f′::Function, x::Interval{T}; debug && println("0 ∉ f′(X)") while true - m = mid(X) - N = m - (f(Interval(m)) / f′(X)) - - debug && (print("Newton step: "); @show (X, X ∩ N)) + m = mid(X) + N = m - (f(interval(m)) / f′(X)) + + debug && (print("Newton step1: "); @show (X, X ∩ N)) + if X == X ∩ N + reps1 += 1 + if reps1 > 20 + reps1 = 0 + break + end + end X = X ∩ N - if isempty(X) + if (isempty(X) || diam(X) == 0) break - elseif 0 ∈ f(Interval(prevfloat(m), nextfloat(m))) - push!(R, Root(X, :unique)) - - debug && @show "Root found", X + elseif 0 ∈ f(interval(prevfloat(mid(X)), nextfloat(mid(X)))) + n = fa = fb = 0 + root_exist = true + while (n < 4 && (fa == 0 || fb == 0)) + if fa == 0 + if 0 ∈ f(interval(prevfloat(X.lo), nextfloat(X.lo))) + fa = 1 + else + N = X.lo - (f(interval(X.lo)) / f′(X)) + X = X ∩ N + if (isempty(X) || diam(X) == 0) + root_exist = false + break + end + end + end + if fb == 0 + if 0 ∈ f(interval(prevfloat(X.hi), nextfloat(X.hi))) + fb = 1 + else + if 0 ∈ f(interval(prevfloat(mid(X)), nextfloat(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) || diam(X) == 0) + root_exist = false + break + end + end + end + N = mid(X) - (f(interval(mid(X))) / f′(X)) + X = X ∩ N + if (isempty(X) || diam(X) == 0) + root_exist = false + break + end + n += 1 + end + if root_exist + push!(R, Root(X, :unique)) + debugroot && @show "Root found", X + end break end end else - # 0 ∈ f'(X) - + if diam(X) == initial_width + reps += 1 + if reps > 10 + push!(R, Root(X, :unknown)) + debugroot && @show "Repititive 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(interval(prevfloat(mid(X)), nextfloat(mid(X)))) # 0 ∈ fⁱ(x) debug && println("0 ∈ fⁱ(x)") - if 0 ∉ f(Interval(X.lo)) + if 0 ∉ f(interval(prevfloat(X.lo), nextfloat(X.lo))) expansion_pt = X.lo - elseif 0 ∉ f(Interval(X.hi)) + elseif 0 ∉ f(interval(prevfloat(X.hi), nextfloat(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(interval(prevfloat(x1), nextfloat(x1))) || 0 ∉ f(interval(prevfloat(x2), nextfloat(x2))) + push!(L, interval(X.lo, m)) + push!(L, interval(m, X.hi)) continue else push!(R, Root(X, :unknown)) - debug && @show "Multiple root found", X + debugroot && @show "Multiple root found", X continue end @@ -91,7 +148,7 @@ function newton1d{T}(f::Function, f′::Function, x::Interval{T}; if (diam(X)/mag(X)) < reltol && diam(f(X)) < abstol push!(R, Root(X, :unknown)) - debug && @show "Tolerance root found", X + debugroot && @show "Tolerance root found", X continue end @@ -102,28 +159,27 @@ function newton1d{T}(f::Function, f′::Function, x::Interval{T}; expansion_pt = mid(X) end - initial_width = diam(X) - a = f(Interval(expansion_pt)) + a = f(interval(expansion_pt)) b = 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 step: "); @show (X, X ∩ N)) + debug && (print("Newton step2: "); @show (X, X ∩ N)) X = X ∩ N m = mid(X) @@ -133,11 +189,6 @@ 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) end end diff --git a/test/newton1d.jl b/test/newton1d.jl index 095f7805..099dce00 100644 --- a/test/newton1d.jl +++ b/test/newton1d.jl @@ -18,14 +18,18 @@ 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 + f2(x) = 4567x^2 - 9134x + 4567 #double root f2′(x) = 9134x - 9134 - f3(x) = (x^2 - 2)^2 + 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) @@ -33,6 +37,8 @@ three_halves_pi = 3*big_pi/2 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) @@ -40,12 +46,14 @@ three_halves_pi = 3*big_pi/2 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 @@ -60,9 +68,19 @@ three_halves_pi = 3*big_pi/2 @test length(rts4) == 1 @test 1 in rts4[1].interval && :unknown == rts4[1].status - L1 = [-sqrt(2), -sqrt(2), sqrt(2), sqrt(2)] + 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 From 9086afc8b6e975af7b22e42d288fc03806aaf328 Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Thu, 31 May 2018 18:05:38 +0530 Subject: [PATCH 11/36] Add slope interval newton method --- src/IntervalRootFinding.jl | 3 ++- src/newton1d.jl | 9 +++++-- test/newton1d.jl | 52 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 61 insertions(+), 3 deletions(-) diff --git a/src/IntervalRootFinding.jl b/src/IntervalRootFinding.jl index e9ffc462..afa9987e 100644 --- a/src/IntervalRootFinding.jl +++ b/src/IntervalRootFinding.jl @@ -18,7 +18,8 @@ export derivative, jacobian, # reexport derivative from ForwardDiff Root, is_unique, roots, find_roots, - bisect, newton1d, slope + bisect, newton1d, slope, + slope_newton1d export isunique, root_status diff --git a/src/newton1d.jl b/src/newton1d.jl index bffb89b1..85c643cb 100644 --- a/src/newton1d.jl +++ b/src/newton1d.jl @@ -25,7 +25,7 @@ function newton1d{T}(f::Function, f′::Function, x::Interval{T}; continue end - if 0 ∉ f′(X) + if 0 ∉ ForwardDiff.derivative(f, X) debug && println("0 ∉ f′(X)") @@ -161,7 +161,7 @@ function newton1d{T}(f::Function, f′::Function, x::Interval{T}; a = f(interval(expansion_pt)) - b = f′(X) + b = ForwardDiff.derivative(f, X) if 0 < b.hi && 0 > b.lo && 0 ∉ a if a.hi < 0 @@ -204,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/test/newton1d.jl b/test/newton1d.jl index 099dce00..2b286652 100644 --- a/test/newton1d.jl +++ b/test/newton1d.jl @@ -84,3 +84,55 @@ three_halves_pi = 3*big_pi/2 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 From 2bb8c8912eca974f37fba5545e874a9cfec0593d Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Fri, 1 Jun 2018 13:48:47 +0530 Subject: [PATCH 12/36] Add slope function for multidim case with tests and benchmarks --- test/slopes.jl | 85 +++++++++++++++++++++++++++++++++----------------- 1 file changed, 57 insertions(+), 28 deletions(-) diff --git a/test/slopes.jl b/test/slopes.jl index d4e54f02..2e26fc81 100644 --- a/test/slopes.jl +++ b/test/slopes.jl @@ -9,41 +9,70 @@ struct Slopes{T} sol::Interval{T} end -@testset "Automatic slope expansion(Float64)" begin - - s = interval(0.75, 1.75) - rts = Slopes{Float64}[] - push!(rts, Slopes(x->((x + sin(x)) * exp(-x^2)), s, mid(s), interval(-2.8, 0.1))) - push!(rts, Slopes(x->(x^4 - 10x^3 + 35x^2 - 50x + 24), s, mid(s), interval(-44, 38.5))) - push!(rts, Slopes(x->((log(x + 1.25) - 0.84x) ^ 2), s, mid(s), interval(-0.16, 0.44))) - push!(rts, Slopes(x->(0.02x^2 - 0.03exp(-(20(x - 0.875))^2)), s, mid(s), interval(0.03, 0.33))) - push!(rts, Slopes(x->(exp(x^2)), s, mid(s), interval(6.03, 33.23))) - push!(rts, Slopes(x->(x^4 - 12x^3 + 47x^2 - 60x - 20exp(-x)), s, mid(s), interval(-39, 65.56))) - push!(rts, Slopes(x->(x^6 - 15x^4 + 27x^2 + 250), s, mid(s), interval(-146.9, 67.1))) - push!(rts, Slopes(x->(atan(cos(tan(x)))), s, mid(s), interval(1, 2))) - push!(rts, Slopes(x->(asin(cos(acos(sin(x))))), s, mid(s), interval(1.36, ∞))) +@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 + for i in 1:length(rts) + @test slope(rts[i].f, rts[i].x, rts[i].c) ⊆ rts[i].sol + end end end -@testset "Automatic slope expansion(BigFloat)" begin - s = interval(BigFloat(0.75), BigFloat(1.75)) - rts = Slopes{BigFloat}[] - push!(rts, Slopes(x->((x + sin(x)) * exp(-x^2)), s, mid(s), interval(BigFloat(-2.8), BigFloat(0.1)))) - push!(rts, Slopes(x->(x^4 - 10x^3 + 35x^2 - 50x + 24), s, mid(s), interval(BigFloat(-44), BigFloat(38.5)))) - push!(rts, Slopes(x->((log(x + 1.25) - 0.84x) ^ 2), s, mid(s), interval(BigFloat(-0.16), BigFloat(0.44)))) - push!(rts, Slopes(x->(0.02x^2 - 0.03exp(-(20(x - 0.875))^2)), s, mid(s), interval(BigFloat(0.03), BigFloat(0.33)))) - push!(rts, Slopes(x->(exp(x^2)), s, mid(s), interval(BigFloat(6.03), BigFloat(33.23)))) - push!(rts, Slopes(x->(x^4 - 12x^3 + 47x^2 - 60x - 20exp(-x)), s, mid(s), interval(BigFloat(-39), BigFloat(65.56)))) - push!(rts, Slopes(x->(x^6 - 15x^4 + 27x^2 + 250), s, mid(s), interval(BigFloat(-146.9), BigFloat(67.1)))) - push!(rts, Slopes(x->(atan(cos(tan(x)))), s, mid(s), interval(BigFloat(1), BigFloat(2)))) - push!(rts, Slopes(x->(asin(cos(acos(sin(x))))), s, mid(s), interval(BigFloat(1.36), BigFloat(∞)))) +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 + @test slope(rts[i].f, rts[i].x, rts[i].c) == rts[i].sol end + end struct SlopesMulti From deb8878ab71ebe4b99355e8ddb19e1dd7741d1f9 Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Sun, 13 May 2018 10:28:51 +0530 Subject: [PATCH 13/36] Add Gauss-Seidel method --- src/linear_eq.jl | 37 +++++++++++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) create mode 100644 src/linear_eq.jl diff --git a/src/linear_eq.jl b/src/linear_eq.jl new file mode 100644 index 00000000..d0e2cbea --- /dev/null +++ b/src/linear_eq.jl @@ -0,0 +1,37 @@ +""" +Preconditions the matrix A and b with the inverse of mid(A) +""" +function preconditioner{T}(A::Matrix{Interval{T}}, b::Array{Interval{T}}) + + Aᶜ = mid.(A) + B = pinv(Aᶜ) + + return B*A, B*b + +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. +""" +function gauss_seidel_interval!{T}(x::Array{Interval{T}}, A::Matrix{Interval{T}}, b::Array{Interval{T}}; precondition=true, maxiter=100) + + precondition && ((M, r) = preconditioner(A, b)) + + n = size(A, 1) + + for iter in 1:maxiter + for i in 1:n + Y = r[i] + for j in 1:n + (i == j) || (Y -= M[i, j] * x[j]) + end + Y = extended_div(Y, M[i, i]) + x[i] = hull((x[i] ∩ Y[1]), x[i] ∩ Y[2]) + end + end + x +end From 41368ed3f69f6114131232f5061d093382e41660 Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Mon, 14 May 2018 19:32:29 +0530 Subject: [PATCH 14/36] Use StaticArrays and add benchmarking --- src/benchmark.jl | 17 ++++++++++++++ src/linear_eq.jl | 58 +++++++++++++++++++++++++++++++++++++++++++++--- 2 files changed, 72 insertions(+), 3 deletions(-) create mode 100644 src/benchmark.jl diff --git a/src/benchmark.jl b/src/benchmark.jl new file mode 100644 index 00000000..e78e9071 --- /dev/null +++ b/src/benchmark.jl @@ -0,0 +1,17 @@ +using IntervalArithmetic, StaticArrays, BenchmarkTools, Compat + +include("linear_eq.jl") + +function randVec(n::Int) + a = randn(n) + A = Interval.(a) + sA = MVector{n}(A) + return A, sA + end + +function randMat(n::Int) + a = randn(n, n) + A = Interval.(a) + sA = MMatrix{n, n}(A) + return A, sA + end diff --git a/src/linear_eq.jl b/src/linear_eq.jl index d0e2cbea..7963cce9 100644 --- a/src/linear_eq.jl +++ b/src/linear_eq.jl @@ -1,21 +1,29 @@ +using IntervalArithmetic """ Preconditions the matrix A and b with the inverse of mid(A) """ function preconditioner{T}(A::Matrix{Interval{T}}, b::Array{Interval{T}}) Aᶜ = mid.(A) - B = pinv(Aᶜ) + B = inv(Aᶜ) return B*A, B*b end +function gauss_seidel_interval{T}(A::Matrix{Interval{T}}, b::Array{Interval{T}}; precondition=true, maxiter=100) + n = size(A, 1) + x = fill(-1e16..1e16, n) + 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!{T}(x::Array{Interval{T}}, A::Matrix{Interval{T}}, b::Array{Interval{T}}; precondition=true, maxiter=100) @@ -29,8 +37,52 @@ function gauss_seidel_interval!{T}(x::Array{Interval{T}}, A::Matrix{Interval{T}} for j in 1:n (i == j) || (Y -= M[i, j] * x[j]) end - Y = extended_div(Y, M[i, i]) - x[i] = hull((x[i] ∩ Y[1]), x[i] ∩ Y[2]) + Z = extended_div(Y, M[i, i]) + x[i] = hull((x[i] ∩ Z[1]), x[i] ∩ Z[2]) + end + end + x +end + +function preconditioner_static{T, N}(A::MMatrix{N, N, Interval{T}}, b::MVector{N, Interval{T}}) + + Aᶜ = mid.(A) + B = inv(Aᶜ) + + return B*A, B*b + +end + + +function gauss_seidel_interval_static{T, N}(A::MMatrix{N, N, Interval{T}}, b::MVector{N, Interval{T}}; precondition=true, maxiter=100) + + n = size(A, 1) + x = @MVector fill(-1e16..1e16, n) + gauss_seidel_interval_static!(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_static!{T, N}(x::MVector{N, Interval{T}}, A::MMatrix{N, N, Interval{T}}, b::MVector{N, Interval{T}}; precondition=true, maxiter=100) + + precondition && ((M, r) = preconditioner_static(A, b)) + + n = size(A, 1) + + for iter in 1:maxiter + for i in 1:n + Y = r[i] + for j in 1:n + (i == j) || (Y -= M[i, j] * x[j]) + end + Z = extended_div(Y, M[i, i]) + x[i] = hull((x[i] ∩ Z[1]), x[i] ∩ Z[2]) end end x From d337490d06302fb4eb35acd038c5da85b00f581a Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Tue, 15 May 2018 15:05:16 +0530 Subject: [PATCH 15/36] Add benchmark script --- perf/linear_eq.jl | 30 ++++++++ perf/linear_eq_results.txt | 149 +++++++++++++++++++++++++++++++++++++ src/benchmark.jl | 17 ----- 3 files changed, 179 insertions(+), 17 deletions(-) create mode 100644 perf/linear_eq.jl create mode 100644 perf/linear_eq_results.txt delete mode 100644 src/benchmark.jl diff --git a/perf/linear_eq.jl b/perf/linear_eq.jl new file mode 100644 index 00000000..a5c32314 --- /dev/null +++ b/perf/linear_eq.jl @@ -0,0 +1,30 @@ +using IntervalArithmetic, StaticArrays, BenchmarkTools, Compat + +include("../src/linear_eq.jl") + +function randVec(n::Int) + a = randn(n) + A = Interval.(a) + sA = MVector{n}(A) + return A, sA +end + +function randMat(n::Int) + a = randn(n, n) + A = Interval.(a) + sA = MMatrix{n, n}(A) + return A, sA +end + +function benchmark(max=10) + for n in 1:max + A, sA = randMat(n) + b, sb = randVec(n) + println("For n = ", n) + t1 = @btime gauss_seidel_interval($A, $b) + println("Array: ", t1) + t2 = @btime gauss_seidel_interval_static($sA, $sb) + println("MArray: ", t2) + println() + end +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/src/benchmark.jl b/src/benchmark.jl deleted file mode 100644 index e78e9071..00000000 --- a/src/benchmark.jl +++ /dev/null @@ -1,17 +0,0 @@ -using IntervalArithmetic, StaticArrays, BenchmarkTools, Compat - -include("linear_eq.jl") - -function randVec(n::Int) - a = randn(n) - A = Interval.(a) - sA = MVector{n}(A) - return A, sA - end - -function randMat(n::Int) - a = randn(n, n) - A = Interval.(a) - sA = MMatrix{n, n}(A) - return A, sA - end From 2c08d6c7216f5c054211270a0e52dfcada5fb968 Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Tue, 15 May 2018 19:26:53 +0530 Subject: [PATCH 16/36] Add another version with StaticArrays --- perf/linear_eq.jl | 18 +++++++++++------- src/linear_eq.jl | 44 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 55 insertions(+), 7 deletions(-) diff --git a/perf/linear_eq.jl b/perf/linear_eq.jl index a5c32314..ffcb7e97 100644 --- a/perf/linear_eq.jl +++ b/perf/linear_eq.jl @@ -5,26 +5,30 @@ include("../src/linear_eq.jl") function randVec(n::Int) a = randn(n) A = Interval.(a) - sA = MVector{n}(A) - return A, sA + mA = MVector{n}(A) + sA = SVector{n}(A) + return A, mA, sA end function randMat(n::Int) a = randn(n, n) A = Interval.(a) - sA = MMatrix{n, n}(A) - return A, sA + mA = MMatrix{n, n}(A) + sA = SMatrix{n, n}(A) + return A, mA, sA end function benchmark(max=10) for n in 1:max - A, sA = randMat(n) - b, sb = randVec(n) + A, mA, sA = randMat(n) + b, mb, sb = randVec(n) println("For n = ", n) t1 = @btime gauss_seidel_interval($A, $b) println("Array: ", t1) - t2 = @btime gauss_seidel_interval_static($sA, $sb) + t2 = @btime gauss_seidel_interval_static($mA, $mb) println("MArray: ", t2) + t3 = @btime gauss_seidel_interval_static1($sA, $sb) + println("SArray: ", t3) println() end end diff --git a/src/linear_eq.jl b/src/linear_eq.jl index 7963cce9..d5b4d1d3 100644 --- a/src/linear_eq.jl +++ b/src/linear_eq.jl @@ -87,3 +87,47 @@ function gauss_seidel_interval_static!{T, N}(x::MVector{N, Interval{T}}, A::MMat end x end + +function preconditioner_static1{T, N}(A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}}) + + Aᶜ = mid.(A) + B = inv(Aᶜ) + + return B*A, B*b + +end + + +function gauss_seidel_interval_static1{T, N}(A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}}; precondition=true, maxiter=100) + + n = size(A, 1) + x = @MVector fill(-1e16..1e16, n) + gauss_seidel_interval_static1!(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_static1!{T, N}(x::MVector{N, Interval{T}}, A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}}; precondition=true, maxiter=100) + + precondition && ((M, r) = preconditioner_static1(A, b)) + + n = size(A, 1) + + for iter in 1:maxiter + for i in 1:n + Y = r[i] + for j in 1:n + (i == j) || (Y -= M[i, j] * x[j]) + end + Z = extended_div(Y, M[i, i]) + x[i] = hull((x[i] ∩ Z[1]), x[i] ∩ Z[2]) + end + end + x +end From a105c0ee8714eeb066030fb456ef5985a69222e9 Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Thu, 17 May 2018 12:18:19 +0530 Subject: [PATCH 17/36] contractor approach --- src/linear_eq.jl | 29 ++++++++++++++++++++++++++++- 1 file changed, 28 insertions(+), 1 deletion(-) diff --git a/src/linear_eq.jl b/src/linear_eq.jl index d5b4d1d3..ea74bf72 100644 --- a/src/linear_eq.jl +++ b/src/linear_eq.jl @@ -1,4 +1,4 @@ -using IntervalArithmetic +using IntervalArithmetic, StaticArrays """ Preconditions the matrix A and b with the inverse of mid(A) """ @@ -131,3 +131,30 @@ function gauss_seidel_interval_static1!{T, N}(x::MVector{N, Interval{T}}, A::SMa end x end + +function gauss_seidel_contractor{T}(A::Matrix{Interval{T}}, b::Array{Interval{T}}; precondition=true, maxiter=100) + + n = size(A, 1) + x = fill(-1e16..1e16, n) + x = gauss_seidel_contractor!(x, A, b, precondition=precondition, maxiter=maxiter) + return x +end + +function gauss_seidel_contractor!{T}(x::Array{Interval{T}}, A::Matrix{Interval{T}}, b::Array{Interval{T}}; precondition=true, maxiter=100) + + precondition && ((A, b) = preconditioner(A, b)) + + n = size(A, 1) + + diagA = Diagonal(A) + extdiagA = deepcopy(A) + for i in 1:n + extdiagA[i, i] = Interval(0) + end + inv_diagA = inv(diagA) + + for iter in 1:maxiter + x = x .∩ (inv_diagA * (b - extdiagA * x)) + end + x +end From 44f93abc32894671a0ac8ec2c070eafeca30e032 Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Thu, 17 May 2018 13:37:56 +0530 Subject: [PATCH 18/36] Add tests --- perf/linear_eq.jl | 2 -- src/IntervalRootFinding.jl | 6 +++++- src/linear_eq.jl | 1 - test/linear_eq.jl | 16 ++++++++++++++++ 4 files changed, 21 insertions(+), 4 deletions(-) create mode 100644 test/linear_eq.jl diff --git a/perf/linear_eq.jl b/perf/linear_eq.jl index ffcb7e97..5027346e 100644 --- a/perf/linear_eq.jl +++ b/perf/linear_eq.jl @@ -1,7 +1,5 @@ using IntervalArithmetic, StaticArrays, BenchmarkTools, Compat -include("../src/linear_eq.jl") - function randVec(n::Int) a = randn(n) A = Interval.(a) diff --git a/src/IntervalRootFinding.jl b/src/IntervalRootFinding.jl index afa9987e..eed4c58b 100644 --- a/src/IntervalRootFinding.jl +++ b/src/IntervalRootFinding.jl @@ -19,7 +19,10 @@ export Root, is_unique, roots, find_roots, bisect, newton1d, slope, - slope_newton1d + slope_newton1d, + gauss_seidel_interval, gauss_seidel_interval!, + gauss_seidel_contractor, gauss_seidel_contractor!, + gauss_seidel_interval_static1, gauss_seidel_interval_static1! export isunique, root_status @@ -59,6 +62,7 @@ include("complex.jl") include("contractors.jl") include("roots.jl") include("newton1d.jl") +include("linear_eq.jl") include("slopes.jl") diff --git a/src/linear_eq.jl b/src/linear_eq.jl index ea74bf72..a17b3e25 100644 --- a/src/linear_eq.jl +++ b/src/linear_eq.jl @@ -1,4 +1,3 @@ -using IntervalArithmetic, StaticArrays """ Preconditions the matrix A and b with the inverse of mid(A) """ diff --git a/test/linear_eq.jl b/test/linear_eq.jl new file mode 100644 index 00000000..5f3f5410 --- /dev/null +++ b/test/linear_eq.jl @@ -0,0 +1,16 @@ +using IntervalArithmetic, StaticArrays, IntervalRootFinding + +@testset "Testing Gauss Seidel Method" 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) + @test all(x[i] .⊆ f(A[i], b[i])) + end + end + end From 55b7999fa3ea3d72c6c64e223fb586faca1bd364 Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Thu, 17 May 2018 15:06:03 +0530 Subject: [PATCH 19/36] replace deepcopy with copy --- src/linear_eq.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/linear_eq.jl b/src/linear_eq.jl index a17b3e25..30581ef3 100644 --- a/src/linear_eq.jl +++ b/src/linear_eq.jl @@ -146,7 +146,7 @@ function gauss_seidel_contractor!{T}(x::Array{Interval{T}}, A::Matrix{Interval{T n = size(A, 1) diagA = Diagonal(A) - extdiagA = deepcopy(A) + extdiagA = copy(A) for i in 1:n extdiagA[i, i] = Interval(0) end From 2b95f04a634f04d48125284d4cca5ead1f05983c Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Thu, 17 May 2018 15:28:52 +0530 Subject: [PATCH 20/36] Add test dependency --- test/linear_eq.jl | 1 + test/runtests.jl | 1 + 2 files changed, 2 insertions(+) diff --git a/test/linear_eq.jl b/test/linear_eq.jl index 5f3f5410..c815209d 100644 --- a/test/linear_eq.jl +++ b/test/linear_eq.jl @@ -14,3 +14,4 @@ using IntervalArithmetic, StaticArrays, IntervalRootFinding end end end +end diff --git a/test/runtests.jl b/test/runtests.jl index 62e2d775..679cd1ad 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,4 +5,5 @@ using Base.Test #include("findroots.jl") include("roots.jl") include("newton1d.jl") +include("linear_eq.jl") include("slopes.jl") From db8710f8653a145131749cf63a8bddb7737d3a88 Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Fri, 18 May 2018 04:42:30 +0530 Subject: [PATCH 21/36] another method with a static vector x --- perf/linear_eq.jl | 4 +++- src/linear_eq.jl | 44 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 47 insertions(+), 1 deletion(-) diff --git a/perf/linear_eq.jl b/perf/linear_eq.jl index 5027346e..4518121e 100644 --- a/perf/linear_eq.jl +++ b/perf/linear_eq.jl @@ -26,7 +26,9 @@ function benchmark(max=10) t2 = @btime gauss_seidel_interval_static($mA, $mb) println("MArray: ", t2) t3 = @btime gauss_seidel_interval_static1($sA, $sb) - println("SArray: ", t3) + println("SArray1: ", t3) + t4 = @btime gauss_seidel_interval_static2($sA, $sb) + println("SArray2: ", t4) println() end end diff --git a/src/linear_eq.jl b/src/linear_eq.jl index 30581ef3..9bd2029e 100644 --- a/src/linear_eq.jl +++ b/src/linear_eq.jl @@ -131,6 +131,50 @@ function gauss_seidel_interval_static1!{T, N}(x::MVector{N, Interval{T}}, A::SMa x end +function preconditioner_static2{T, N}(A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}}) + + Aᶜ = mid.(A) + B = inv(Aᶜ) + + return B*A, B*b + +end + + +function gauss_seidel_interval_static2{T, N}(A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}}; precondition=true, maxiter=100) + + n = size(A, 1) + x = @SVector fill(-1e16..1e16, n) + x = gauss_seidel_interval_static2!(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_static2!{T, N}(x::SVector{N, Interval{T}}, A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}}; precondition=true, maxiter=100) + + precondition && ((M, r) = preconditioner_static2(A, b)) + + n = size(A, 1) + + for iter in 1:maxiter + for i in 1:n + Y = r[i] + for j in 1:n + (i == j) || (Y -= M[i, j] * x[j]) + end + Z = extended_div(Y, M[i, i]) + x = setindex(x, hull((x[i] ∩ Z[1]), x[i] ∩ Z[2]), i) + end + end + x +end + function gauss_seidel_contractor{T}(A::Matrix{Interval{T}}, b::Array{Interval{T}}; precondition=true, maxiter=100) n = size(A, 1) From 7ee7ca695ed121fc90a23e5fbf6148def184f0b4 Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Fri, 18 May 2018 09:15:20 +0530 Subject: [PATCH 22/36] bug fix --- src/linear_eq.jl | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/src/linear_eq.jl b/src/linear_eq.jl index 9bd2029e..b818a1f7 100644 --- a/src/linear_eq.jl +++ b/src/linear_eq.jl @@ -26,17 +26,17 @@ Eldon Hansen and G. William Walster : Global Optimization Using Interval Analysi """ function gauss_seidel_interval!{T}(x::Array{Interval{T}}, A::Matrix{Interval{T}}, b::Array{Interval{T}}; precondition=true, maxiter=100) - precondition && ((M, r) = preconditioner(A, b)) + precondition && ((A, b) = preconditioner(A, b)) n = size(A, 1) for iter in 1:maxiter for i in 1:n - Y = r[i] + Y = b[i] for j in 1:n - (i == j) || (Y -= M[i, j] * x[j]) + (i == j) || (Y -= A[i, j] * x[j]) end - Z = extended_div(Y, M[i, i]) + Z = extended_div(Y, A[i, i]) x[i] = hull((x[i] ∩ Z[1]), x[i] ∩ Z[2]) end end @@ -70,17 +70,17 @@ Eldon Hansen and G. William Walster : Global Optimization Using Interval Analysi """ function gauss_seidel_interval_static!{T, N}(x::MVector{N, Interval{T}}, A::MMatrix{N, N, Interval{T}}, b::MVector{N, Interval{T}}; precondition=true, maxiter=100) - precondition && ((M, r) = preconditioner_static(A, b)) + precondition && ((A, b) = preconditioner_static(A, b)) n = size(A, 1) for iter in 1:maxiter for i in 1:n - Y = r[i] + Y = b[i] for j in 1:n - (i == j) || (Y -= M[i, j] * x[j]) + (i == j) || (Y -= A[i, j] * x[j]) end - Z = extended_div(Y, M[i, i]) + Z = extended_div(Y, A[i, i]) x[i] = hull((x[i] ∩ Z[1]), x[i] ∩ Z[2]) end end @@ -114,17 +114,17 @@ Eldon Hansen and G. William Walster : Global Optimization Using Interval Analysi """ function gauss_seidel_interval_static1!{T, N}(x::MVector{N, Interval{T}}, A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}}; precondition=true, maxiter=100) - precondition && ((M, r) = preconditioner_static1(A, b)) + precondition && ((A, b) = preconditioner_static1(A, b)) n = size(A, 1) for iter in 1:maxiter for i in 1:n - Y = r[i] + Y = b[i] for j in 1:n - (i == j) || (Y -= M[i, j] * x[j]) + (i == j) || (Y -= A[i, j] * x[j]) end - Z = extended_div(Y, M[i, i]) + Z = extended_div(Y, A[i, i]) x[i] = hull((x[i] ∩ Z[1]), x[i] ∩ Z[2]) end end @@ -158,17 +158,17 @@ Eldon Hansen and G. William Walster : Global Optimization Using Interval Analysi """ function gauss_seidel_interval_static2!{T, N}(x::SVector{N, Interval{T}}, A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}}; precondition=true, maxiter=100) - precondition && ((M, r) = preconditioner_static2(A, b)) + precondition && ((A, b) = preconditioner_static2(A, b)) n = size(A, 1) for iter in 1:maxiter for i in 1:n - Y = r[i] + Y = b[i] for j in 1:n - (i == j) || (Y -= M[i, j] * x[j]) + (i == j) || (Y -= A[i, j] * x[j]) end - Z = extended_div(Y, M[i, i]) + Z = extended_div(Y, A[i, i]) x = setindex(x, hull((x[i] ∩ Z[1]), x[i] ∩ Z[2]), i) end end From b1cb7a5e924f00911e18dc8ed62199c0201c3688 Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Fri, 18 May 2018 10:20:52 +0530 Subject: [PATCH 23/36] Display benchmark data in a table --- perf/linear_eq.jl | 21 ++++++++++----------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/perf/linear_eq.jl b/perf/linear_eq.jl index 4518121e..a41cbd1d 100644 --- a/perf/linear_eq.jl +++ b/perf/linear_eq.jl @@ -1,4 +1,4 @@ -using IntervalArithmetic, StaticArrays, BenchmarkTools, Compat +using BenchmarkTools, Compat, DataFrames, IntervalRootFinding, IntervalArithmetic, StaticArrays function randVec(n::Int) a = randn(n) @@ -17,18 +17,17 @@ function randMat(n::Int) end function benchmark(max=10) + df = DataFrame() + df[:Method] = ["Array", "MArray", "SArray1", "SArray2", "Contractor"] for n in 1:max A, mA, sA = randMat(n) b, mb, sb = randVec(n) - println("For n = ", n) - t1 = @btime gauss_seidel_interval($A, $b) - println("Array: ", t1) - t2 = @btime gauss_seidel_interval_static($mA, $mb) - println("MArray: ", t2) - t3 = @btime gauss_seidel_interval_static1($sA, $sb) - println("SArray1: ", t3) - t4 = @btime gauss_seidel_interval_static2($sA, $sb) - println("SArray2: ", t4) - println() + t1 = @belapsed gauss_seidel_interval($A, $b) + t2 = @belapsed gauss_seidel_interval_static($mA, $mb) + t3 = @belapsed gauss_seidel_interval_static1($sA, $sb) + t4 = @belapsed gauss_seidel_interval_static2($sA, $sb) + t5 = @belapsed gauss_seidel_contractor($A, $b) + df[Symbol("n = $n")] = [t1, t2, t3, t4, t5] end + println(df) end From 168aff6a076dbdfc009c9da588be2eacd34038c1 Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Fri, 18 May 2018 10:52:20 +0530 Subject: [PATCH 24/36] add tabular results --- perf/linear_eq.jl | 1 + perf/linear_eq_tabular.txt | 8 ++++++++ 2 files changed, 9 insertions(+) create mode 100644 perf/linear_eq_tabular.txt diff --git a/perf/linear_eq.jl b/perf/linear_eq.jl index a41cbd1d..01773a75 100644 --- a/perf/linear_eq.jl +++ b/perf/linear_eq.jl @@ -30,4 +30,5 @@ function benchmark(max=10) df[Symbol("n = $n")] = [t1, t2, t3, t4, t5] end println(df) + df end diff --git a/perf/linear_eq_tabular.txt b/perf/linear_eq_tabular.txt new file mode 100644 index 00000000..f0d99d30 --- /dev/null +++ b/perf/linear_eq_tabular.txt @@ -0,0 +1,8 @@ + +│ Row │ Method │ n = 1 │ n = 2 │ n = 3 │ n = 4 │ n = 5 │ n = 6 │ n = 7 │ n = 8 │ n = 9 │ n = 10 │ n = 11 │ n = 12 │ n = 13 │ +├─────┼────────────┼───────────┼───────────┼───────────┼─────────────┼─────────────┼─────────────┼─────────────┼─────────────┼─────────────┼─────────────┼─────────────┼─────────────┼─────────────┤ +│ 1 │ Array │ 1.0224e-5 │ 2.4714e-5 │ 4.5112e-5 │ 7.3156e-5 │ 0.000107752 │ 0.000146873 │ 0.000193987 │ 0.000246073 │ 0.000315383 │ 0.000383108 │ 0.000466884 │ 0.0005546 │ 0.000656006 │ +│ 2 │ MArray │ 1.5073e-5 │ 3.9778e-5 │ 8.3171e-5 │ 0.000147665 │ 0.000238183 │ 0.000384167 │ 0.000641308 │ 0.00112101 │ 0.00178696 │ 0.00268557 │ 0.00495489 │ 0.0106303 │ 0.0149287 │ +│ 3 │ SArray1 │ 1.0526e-5 │ 2.4568e-5 │ 4.4919e-5 │ 6.9534e-5 │ 0.000109565 │ 0.000150989 │ 0.00019636 │ 0.00025249 │ 0.000315217 │ 0.000386076 │ 0.000474434 │ 0.000566991 │ 0.000673982 │ +│ 4 │ SArray2 │ 9.858e-6 │ 2.4191e-5 │ 4.4197e-5 │ 6.8149e-5 │ 0.000102997 │ 0.000143808 │ 0.000187444 │ 0.000239507 │ 0.000297437 │ 0.000368891 │ 0.000449486 │ 0.000537369 │ 0.000639082 │ +│ 5 │ Contractor │ 2.2636e-5 │ 4.0218e-5 │ 6.3581e-5 │ 9.4589e-5 │ 0.000127638 │ 0.000172273 │ 0.000222144 │ 0.000279826 │ 0.000347763 │ 0.000415261 │ 0.000500709 │ 0.000617629 │ 0.000719819 │ From 81731663ef1da3b9aaf4f76b7654210463879b0e Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Fri, 18 May 2018 13:06:38 +0530 Subject: [PATCH 25/36] make results more readable --- perf/linear_eq.jl | 10 ++++++++-- perf/linear_eq_tabular.txt | 30 ++++++++++++++++++++++-------- 2 files changed, 30 insertions(+), 10 deletions(-) diff --git a/perf/linear_eq.jl b/perf/linear_eq.jl index 01773a75..c44a493e 100644 --- a/perf/linear_eq.jl +++ b/perf/linear_eq.jl @@ -29,6 +29,12 @@ function benchmark(max=10) t5 = @belapsed gauss_seidel_contractor($A, $b) df[Symbol("n = $n")] = [t1, t2, t3, t4, t5] end - println(df) - df + a = [] + for i in 1:max + push!(a, Symbol("n = $i")) + end + df1 = stack(df, a) + dfnew = unstack(df1, :variable, :Method, :value) + println(dfnew) + dfnew end diff --git a/perf/linear_eq_tabular.txt b/perf/linear_eq_tabular.txt index f0d99d30..df8f118f 100644 --- a/perf/linear_eq_tabular.txt +++ b/perf/linear_eq_tabular.txt @@ -1,8 +1,22 @@ - -│ Row │ Method │ n = 1 │ n = 2 │ n = 3 │ n = 4 │ n = 5 │ n = 6 │ n = 7 │ n = 8 │ n = 9 │ n = 10 │ n = 11 │ n = 12 │ n = 13 │ -├─────┼────────────┼───────────┼───────────┼───────────┼─────────────┼─────────────┼─────────────┼─────────────┼─────────────┼─────────────┼─────────────┼─────────────┼─────────────┼─────────────┤ -│ 1 │ Array │ 1.0224e-5 │ 2.4714e-5 │ 4.5112e-5 │ 7.3156e-5 │ 0.000107752 │ 0.000146873 │ 0.000193987 │ 0.000246073 │ 0.000315383 │ 0.000383108 │ 0.000466884 │ 0.0005546 │ 0.000656006 │ -│ 2 │ MArray │ 1.5073e-5 │ 3.9778e-5 │ 8.3171e-5 │ 0.000147665 │ 0.000238183 │ 0.000384167 │ 0.000641308 │ 0.00112101 │ 0.00178696 │ 0.00268557 │ 0.00495489 │ 0.0106303 │ 0.0149287 │ -│ 3 │ SArray1 │ 1.0526e-5 │ 2.4568e-5 │ 4.4919e-5 │ 6.9534e-5 │ 0.000109565 │ 0.000150989 │ 0.00019636 │ 0.00025249 │ 0.000315217 │ 0.000386076 │ 0.000474434 │ 0.000566991 │ 0.000673982 │ -│ 4 │ SArray2 │ 9.858e-6 │ 2.4191e-5 │ 4.4197e-5 │ 6.8149e-5 │ 0.000102997 │ 0.000143808 │ 0.000187444 │ 0.000239507 │ 0.000297437 │ 0.000368891 │ 0.000449486 │ 0.000537369 │ 0.000639082 │ -│ 5 │ Contractor │ 2.2636e-5 │ 4.0218e-5 │ 6.3581e-5 │ 9.4589e-5 │ 0.000127638 │ 0.000172273 │ 0.000222144 │ 0.000279826 │ 0.000347763 │ 0.000415261 │ 0.000500709 │ 0.000617629 │ 0.000719819 │ +│ 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 │ From 138839784423bcdfdad0dc719eefbca2941ca620 Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Sat, 19 May 2018 19:03:00 +0530 Subject: [PATCH 26/36] generalized the various implementation --- perf/linear_eq.jl | 16 ++- perf/linear_eq_tabular.txt | 14 ++ src/IntervalRootFinding.jl | 6 +- src/linear_eq.jl | 288 +++++++++++++++++++------------------ 4 files changed, 178 insertions(+), 146 deletions(-) diff --git a/perf/linear_eq.jl b/perf/linear_eq.jl index c44a493e..e554f92b 100644 --- a/perf/linear_eq.jl +++ b/perf/linear_eq.jl @@ -18,23 +18,25 @@ end function benchmark(max=10) df = DataFrame() - df[:Method] = ["Array", "MArray", "SArray1", "SArray2", "Contractor"] + 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_static($mA, $mb) - t3 = @belapsed gauss_seidel_interval_static1($sA, $sb) - t4 = @belapsed gauss_seidel_interval_static2($sA, $sb) - t5 = @belapsed gauss_seidel_contractor($A, $b) - df[Symbol("n = $n")] = [t1, t2, t3, t4, t5] + 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("n = $i")) + 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_tabular.txt b/perf/linear_eq_tabular.txt index df8f118f..85a8d916 100644 --- a/perf/linear_eq_tabular.txt +++ b/perf/linear_eq_tabular.txt @@ -20,3 +20,17 @@ │ 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 │ diff --git a/src/IntervalRootFinding.jl b/src/IntervalRootFinding.jl index eed4c58b..35a4ba96 100644 --- a/src/IntervalRootFinding.jl +++ b/src/IntervalRootFinding.jl @@ -21,8 +21,10 @@ export bisect, newton1d, slope, slope_newton1d, gauss_seidel_interval, gauss_seidel_interval!, - gauss_seidel_contractor, gauss_seidel_contractor!, - gauss_seidel_interval_static1, gauss_seidel_interval_static1! + gauss_seidel_contractor, gauss_seidel_contractor! + # gauss_seidel_interval_static1, gauss_seidel_interval_static1! + # gauss_seidel_interval_static2, gauss_seidel_interval_static2! + # gauss_seidel_interval_static, gauss_seidel_interval_static! export isunique, root_status diff --git a/src/linear_eq.jl b/src/linear_eq.jl index b818a1f7..83e96ba7 100644 --- a/src/linear_eq.jl +++ b/src/linear_eq.jl @@ -1,7 +1,7 @@ """ Preconditions the matrix A and b with the inverse of mid(A) """ -function preconditioner{T}(A::Matrix{Interval{T}}, b::Array{Interval{T}}) +function preconditioner(A::AbstractMatrix, b::AbstractArray) Aᶜ = mid.(A) B = inv(Aᶜ) @@ -10,101 +10,18 @@ function preconditioner{T}(A::Matrix{Interval{T}}, b::Array{Interval{T}}) end -function gauss_seidel_interval{T}(A::Matrix{Interval{T}}, b::Array{Interval{T}}; precondition=true, maxiter=100) +function gauss_seidel_interval(A::AbstractMatrix, b::AbstractArray; precondition=true, maxiter=100) n = size(A, 1) x = fill(-1e16..1e16, n) - 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!{T}(x::Array{Interval{T}}, A::Matrix{Interval{T}}, b::Array{Interval{T}}; precondition=true, maxiter=100) - - precondition && ((A, b) = preconditioner(A, b)) - - n = size(A, 1) - - 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 preconditioner_static{T, N}(A::MMatrix{N, N, Interval{T}}, b::MVector{N, Interval{T}}) - - Aᶜ = mid.(A) - B = inv(Aᶜ) - - return B*A, B*b - -end - - -function gauss_seidel_interval_static{T, N}(A::MMatrix{N, N, Interval{T}}, b::MVector{N, Interval{T}}; precondition=true, maxiter=100) - n = size(A, 1) - x = @MVector fill(-1e16..1e16, n) - gauss_seidel_interval_static!(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_static!{T, N}(x::MVector{N, Interval{T}}, A::MMatrix{N, N, Interval{T}}, b::MVector{N, Interval{T}}; precondition=true, maxiter=100) - - precondition && ((A, b) = preconditioner_static(A, b)) - - n = size(A, 1) - - 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 + if (typeof(b) <: SArray || typeof(b) <: MArray) + x = MVector{n}(x) end - x -end - -function preconditioner_static1{T, N}(A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}}) - - Aᶜ = mid.(A) - B = inv(Aᶜ) - - return B*A, B*b - -end - - -function gauss_seidel_interval_static1{T, N}(A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}}; precondition=true, maxiter=100) - n = size(A, 1) - x = @MVector fill(-1e16..1e16, n) - gauss_seidel_interval_static1!(x, A, b, precondition=precondition, maxiter=maxiter) + 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 @@ -112,13 +29,13 @@ 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_static1!{T, N}(x::MVector{N, Interval{T}}, A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}}; precondition=true, maxiter=100) +function gauss_seidel_interval!(x::AbstractArray, A::AbstractMatrix, b::AbstractArray; precondition=true, maxiter=100) - precondition && ((A, b) = preconditioner_static1(A, b)) + precondition && ((A, b) = preconditioner(A, b)) n = size(A, 1) - for iter in 1:maxiter + @inbounds for iter in 1:maxiter for i in 1:n Y = b[i] for j in 1:n @@ -130,60 +47,153 @@ function gauss_seidel_interval_static1!{T, N}(x::MVector{N, Interval{T}}, A::SMa end x end - -function preconditioner_static2{T, N}(A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}}) - - Aᶜ = mid.(A) - B = inv(Aᶜ) - - return B*A, B*b - -end - - -function gauss_seidel_interval_static2{T, N}(A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}}; precondition=true, maxiter=100) - - n = size(A, 1) - x = @SVector fill(-1e16..1e16, n) - x = gauss_seidel_interval_static2!(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_static2!{T, N}(x::SVector{N, Interval{T}}, A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}}; precondition=true, maxiter=100) - - precondition && ((A, b) = preconditioner_static2(A, b)) +# +# function preconditioner_static{T, N}(A::MMatrix{N, N, Interval{T}}, b::MVector{N, Interval{T}}) +# +# Aᶜ = mid.(A) +# B = inv(Aᶜ) +# +# return B*A, B*b +# +# end +# +# +# function gauss_seidel_interval_static{T, N}(A::MMatrix{N, N, Interval{T}}, b::MVector{N, Interval{T}}; precondition=true, maxiter=100) +# +# n = size(A, 1) +# x = @MVector fill(-1e16..1e16, n) +# gauss_seidel_interval_static!(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_static!{T, N}(x::MVector{N, Interval{T}}, A::MMatrix{N, N, Interval{T}}, b::MVector{N, Interval{T}}; precondition=true, maxiter=100) +# +# precondition && ((A, b) = preconditioner_static(A, b)) +# +# n = size(A, 1) +# +# 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 preconditioner_static1{T, N}(A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}}) +# +# Aᶜ = mid.(A) +# B = inv(Aᶜ) +# +# return B*A, B*b +# +# end +# +# +# function gauss_seidel_interval_static1{T, N}(A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}}; precondition=true, maxiter=100) +# +# n = size(A, 1) +# x = @MVector fill(-1e16..1e16, n) +# gauss_seidel_interval_static1!(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_static1!{T, N}(x::MVector{N, Interval{T}}, A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}}; precondition=true, maxiter=100) +# +# precondition && ((A, b) = preconditioner_static1(A, b)) +# +# n = size(A, 1) +# +# 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 preconditioner_static2{T, N}(A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}}) +# +# Aᶜ = mid.(A) +# B = inv(Aᶜ) +# +# return B*A, B*b +# +# end +# +# +# function gauss_seidel_interval_static2{T, N}(A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}}; precondition=true, maxiter=100) +# +# n = size(A, 1) +# x = @SVector fill(-1e16..1e16, n) +# x = gauss_seidel_interval_static2!(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_static2!{T, N}(x::SVector{N, Interval{T}}, A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}}; precondition=true, maxiter=100) +# +# precondition && ((A, b) = preconditioner_static2(A, b)) +# +# n = size(A, 1) +# +# 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 = setindex(x, hull((x[i] ∩ Z[1]), x[i] ∩ Z[2]), i) +# end +# end +# x +# end + +function gauss_seidel_contractor(A::AbstractMatrix, b::AbstractArray; precondition=true, maxiter=100) n = size(A, 1) + x = fill(-1e16..1e16, n) - 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 = setindex(x, hull((x[i] ∩ Z[1]), x[i] ∩ Z[2]), i) - end + if (typeof(b) <: SArray || typeof(b) <: MArray) + x = MVector{n}(x) end - x -end - -function gauss_seidel_contractor{T}(A::Matrix{Interval{T}}, b::Array{Interval{T}}; precondition=true, maxiter=100) - n = size(A, 1) - x = fill(-1e16..1e16, n) x = gauss_seidel_contractor!(x, A, b, precondition=precondition, maxiter=maxiter) return x end -function gauss_seidel_contractor!{T}(x::Array{Interval{T}}, A::Matrix{Interval{T}}, b::Array{Interval{T}}; precondition=true, maxiter=100) +function gauss_seidel_contractor!(x::AbstractArray, A::AbstractMatrix, b::AbstractArray; precondition=true, maxiter=100) precondition && ((A, b) = preconditioner(A, b)) @@ -192,7 +202,11 @@ function gauss_seidel_contractor!{T}(x::Array{Interval{T}}, A::Matrix{Interval{T diagA = Diagonal(A) extdiagA = copy(A) for i in 1:n - extdiagA[i, i] = Interval(0) + if (typeof(b) <: SArray) + extdiagA = setindex(extdiagA, Interval(0), i, i) + else + extdiagA[i, i] = Interval(0) + end end inv_diagA = inv(diagA) From 7c9cbd91694790946e4782650d39041df0a13151 Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Sun, 20 May 2018 00:21:50 +0530 Subject: [PATCH 27/36] Add Gauss Elimination with preconditioning --- perf/linear_eq.jl | 21 +++++ src/linear_eq.jl | 216 ++++++++++++++++------------------------------ 2 files changed, 96 insertions(+), 141 deletions(-) diff --git a/perf/linear_eq.jl b/perf/linear_eq.jl index e554f92b..a418bfba 100644 --- a/perf/linear_eq.jl +++ b/perf/linear_eq.jl @@ -40,3 +40,24 @@ function benchmark(max=10) 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/src/linear_eq.jl b/src/linear_eq.jl index 83e96ba7..bc394492 100644 --- a/src/linear_eq.jl +++ b/src/linear_eq.jl @@ -1,3 +1,4 @@ +using IntervalArithmetic """ Preconditions the matrix A and b with the inverse of mid(A) """ @@ -15,10 +16,6 @@ function gauss_seidel_interval(A::AbstractMatrix, b::AbstractArray; precondition n = size(A, 1) x = fill(-1e16..1e16, n) - if (typeof(b) <: SArray || typeof(b) <: MArray) - x = MVector{n}(x) - end - gauss_seidel_interval!(x, A, b, precondition=precondition, maxiter=maxiter) return x end @@ -47,148 +44,11 @@ function gauss_seidel_interval!(x::AbstractArray, A::AbstractMatrix, b::Abstract end x end -# -# function preconditioner_static{T, N}(A::MMatrix{N, N, Interval{T}}, b::MVector{N, Interval{T}}) -# -# Aᶜ = mid.(A) -# B = inv(Aᶜ) -# -# return B*A, B*b -# -# end -# -# -# function gauss_seidel_interval_static{T, N}(A::MMatrix{N, N, Interval{T}}, b::MVector{N, Interval{T}}; precondition=true, maxiter=100) -# -# n = size(A, 1) -# x = @MVector fill(-1e16..1e16, n) -# gauss_seidel_interval_static!(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_static!{T, N}(x::MVector{N, Interval{T}}, A::MMatrix{N, N, Interval{T}}, b::MVector{N, Interval{T}}; precondition=true, maxiter=100) -# -# precondition && ((A, b) = preconditioner_static(A, b)) -# -# n = size(A, 1) -# -# 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 preconditioner_static1{T, N}(A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}}) -# -# Aᶜ = mid.(A) -# B = inv(Aᶜ) -# -# return B*A, B*b -# -# end -# -# -# function gauss_seidel_interval_static1{T, N}(A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}}; precondition=true, maxiter=100) -# -# n = size(A, 1) -# x = @MVector fill(-1e16..1e16, n) -# gauss_seidel_interval_static1!(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_static1!{T, N}(x::MVector{N, Interval{T}}, A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}}; precondition=true, maxiter=100) -# -# precondition && ((A, b) = preconditioner_static1(A, b)) -# -# n = size(A, 1) -# -# 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 preconditioner_static2{T, N}(A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}}) -# -# Aᶜ = mid.(A) -# B = inv(Aᶜ) -# -# return B*A, B*b -# -# end -# -# -# function gauss_seidel_interval_static2{T, N}(A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}}; precondition=true, maxiter=100) -# -# n = size(A, 1) -# x = @SVector fill(-1e16..1e16, n) -# x = gauss_seidel_interval_static2!(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_static2!{T, N}(x::SVector{N, Interval{T}}, A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}}; precondition=true, maxiter=100) -# -# precondition && ((A, b) = preconditioner_static2(A, b)) -# -# n = size(A, 1) -# -# 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 = setindex(x, hull((x[i] ∩ Z[1]), x[i] ∩ Z[2]), i) -# end -# end -# x -# end function gauss_seidel_contractor(A::AbstractMatrix, b::AbstractArray; precondition=true, maxiter=100) n = size(A, 1) x = fill(-1e16..1e16, n) - - if (typeof(b) <: SArray || typeof(b) <: MArray) - x = MVector{n}(x) - end - x = gauss_seidel_contractor!(x, A, b, precondition=precondition, maxiter=maxiter) return x end @@ -215,3 +75,77 @@ function gauss_seidel_contractor!(x::AbstractArray, A::AbstractMatrix, b::Abstra end x end + +function gauss_elimination_interval(A::AbstractMatrix, b::AbstractArray; precondition=true) + + n = size(A, 1) + x = fill(-1e16..1e16, n) + + x = gauss_seidel_interval!(x, A, b, precondition=precondition) + + 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_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_seidel_interval!(x, A, b, precondition=precondition) + + 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_elimination_interval1!(x::AbstractArray, a::AbstractMatrix, b::AbstractArray; precondition=true) + + precondition && ((a, b) = preconditioner(a, b)) + + a \ b +end From 1d8b6b19fe649848fd3c51ac6a2d263422e22f6a Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Sun, 20 May 2018 00:47:52 +0530 Subject: [PATCH 28/36] Add benchmarking data --- perf/linear_eq_tabular.txt | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/perf/linear_eq_tabular.txt b/perf/linear_eq_tabular.txt index 85a8d916..4922245f 100644 --- a/perf/linear_eq_tabular.txt +++ b/perf/linear_eq_tabular.txt @@ -1,3 +1,5 @@ +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 │ @@ -34,3 +36,29 @@ │ 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 + +20×3 DataFrames.DataFrame +│ Row │ n │ Base.\\ │ Gauss Elimination │ +├─────┼────┼─────────────┼───────────────────┤ +│ 1 │ 1 │ 1.0859e-5 │ 1.0402e-5 │ +│ 2 │ 10 │ 0.00037311 │ 0.000372278 │ +│ 3 │ 11 │ 0.000457845 │ 0.000457603 │ +│ 4 │ 12 │ 0.000543028 │ 0.00054428 │ +│ 5 │ 13 │ 0.000642694 │ 0.00064314 │ +│ 6 │ 14 │ 0.00076788 │ 0.000766503 │ +│ 7 │ 15 │ 0.000877223 │ 0.000874887 │ +│ 8 │ 16 │ 0.00102177 │ 0.00102244 │ +│ 9 │ 17 │ 0.00114972 │ 0.00115541 │ +│ 10 │ 18 │ 0.00150682 │ 0.00153407 │ +│ 11 │ 19 │ 0.00172544 │ 0.00170808 │ +│ 12 │ 2 │ 2.6519e-5 │ 2.5677e-5 │ +│ 13 │ 20 │ 0.00192433 │ 0.0019773 │ +│ 14 │ 3 │ 4.6299e-5 │ 4.6316e-5 │ +│ 15 │ 4 │ 7.2201e-5 │ 7.2106e-5 │ +│ 16 │ 5 │ 0.000101941 │ 0.000102027 │ +│ 17 │ 6 │ 0.000139823 │ 0.000139717 │ +│ 18 │ 7 │ 0.000189156 │ 0.000188991 │ +│ 19 │ 8 │ 0.000239304 │ 0.000238721 │ +│ 20 │ 9 │ 0.000305721 │ 0.000305388 │ From f4a6b53f9cd4ceb4cbc08156cddfaedeb76716b3 Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Sun, 20 May 2018 13:18:26 +0530 Subject: [PATCH 29/36] Add tests and fix typos --- perf/linear_eq_tabular.txt | 31 ++++++++----------------------- src/IntervalRootFinding.jl | 6 ++---- src/linear_eq.jl | 19 ++++++------------- test/linear_eq.jl | 4 ++-- 4 files changed, 18 insertions(+), 42 deletions(-) diff --git a/perf/linear_eq_tabular.txt b/perf/linear_eq_tabular.txt index 4922245f..69028385 100644 --- a/perf/linear_eq_tabular.txt +++ b/perf/linear_eq_tabular.txt @@ -39,26 +39,11 @@ Gauss Seidel Gauss Elimination -20×3 DataFrames.DataFrame -│ Row │ n │ Base.\\ │ Gauss Elimination │ -├─────┼────┼─────────────┼───────────────────┤ -│ 1 │ 1 │ 1.0859e-5 │ 1.0402e-5 │ -│ 2 │ 10 │ 0.00037311 │ 0.000372278 │ -│ 3 │ 11 │ 0.000457845 │ 0.000457603 │ -│ 4 │ 12 │ 0.000543028 │ 0.00054428 │ -│ 5 │ 13 │ 0.000642694 │ 0.00064314 │ -│ 6 │ 14 │ 0.00076788 │ 0.000766503 │ -│ 7 │ 15 │ 0.000877223 │ 0.000874887 │ -│ 8 │ 16 │ 0.00102177 │ 0.00102244 │ -│ 9 │ 17 │ 0.00114972 │ 0.00115541 │ -│ 10 │ 18 │ 0.00150682 │ 0.00153407 │ -│ 11 │ 19 │ 0.00172544 │ 0.00170808 │ -│ 12 │ 2 │ 2.6519e-5 │ 2.5677e-5 │ -│ 13 │ 20 │ 0.00192433 │ 0.0019773 │ -│ 14 │ 3 │ 4.6299e-5 │ 4.6316e-5 │ -│ 15 │ 4 │ 7.2201e-5 │ 7.2106e-5 │ -│ 16 │ 5 │ 0.000101941 │ 0.000102027 │ -│ 17 │ 6 │ 0.000139823 │ 0.000139717 │ -│ 18 │ 7 │ 0.000189156 │ 0.000188991 │ -│ 19 │ 8 │ 0.000239304 │ 0.000238721 │ -│ 20 │ 9 │ 0.000305721 │ 0.000305388 │ +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/src/IntervalRootFinding.jl b/src/IntervalRootFinding.jl index 35a4ba96..da64ac5b 100644 --- a/src/IntervalRootFinding.jl +++ b/src/IntervalRootFinding.jl @@ -21,10 +21,8 @@ export bisect, newton1d, slope, slope_newton1d, gauss_seidel_interval, gauss_seidel_interval!, - gauss_seidel_contractor, gauss_seidel_contractor! - # gauss_seidel_interval_static1, gauss_seidel_interval_static1! - # gauss_seidel_interval_static2, gauss_seidel_interval_static2! - # gauss_seidel_interval_static, gauss_seidel_interval_static! + gauss_seidel_contractor, gauss_seidel_contractor!, + gauss_elimination_interval, gauss_elimination_interval! export isunique, root_status diff --git a/src/linear_eq.jl b/src/linear_eq.jl index bc394492..f86e5bcc 100644 --- a/src/linear_eq.jl +++ b/src/linear_eq.jl @@ -1,4 +1,3 @@ -using IntervalArithmetic """ Preconditions the matrix A and b with the inverse of mid(A) """ @@ -81,16 +80,14 @@ function gauss_elimination_interval(A::AbstractMatrix, b::AbstractArray; precond n = size(A, 1) x = fill(-1e16..1e16, n) - x = gauss_seidel_interval!(x, A, b, precondition=precondition) + x = gauss_elimination_interval!(x, A, b, precondition=precondition) 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 +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) @@ -132,16 +129,12 @@ function gauss_elimination_interval1(A::AbstractMatrix, b::AbstractArray; precon n = size(A, 1) x = fill(-1e16..1e16, n) - x = gauss_seidel_interval!(x, A, b, precondition=precondition) + x = gauss_elimination_interval1!(x, A, b, precondition=precondition) 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 +Using `Base.\`` """ function gauss_elimination_interval1!(x::AbstractArray, a::AbstractMatrix, b::AbstractArray; precondition=true) diff --git a/test/linear_eq.jl b/test/linear_eq.jl index c815209d..f2972d4b 100644 --- a/test/linear_eq.jl +++ b/test/linear_eq.jl @@ -1,6 +1,6 @@ using IntervalArithmetic, StaticArrays, IntervalRootFinding -@testset "Testing Gauss Seidel Method" begin +@testset "Linear Equations" begin A = [[2..3 0..1;1..2 2..3], ] b = [[0..120, 60..240], ] @@ -9,7 +9,7 @@ using IntervalArithmetic, StaticArrays, IntervalRootFinding for i in 1:length(A) for precondition in (false, true) - for f in (gauss_seidel_interval, gauss_seidel_contractor) + for f in (gauss_seidel_interval, gauss_seidel_contractor, gauss_elimination_interval) @test all(x[i] .⊆ f(A[i], b[i])) end end From fa378a82494ecca56317fed2ddc3b75fb9f738ab Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Sun, 20 May 2018 13:51:51 +0530 Subject: [PATCH 30/36] Add examples --- examples/linear_eq.jl | 32 ++++++++++++++++++++++++++++++++ src/linear_eq.jl | 11 ++++++----- 2 files changed, 38 insertions(+), 5 deletions(-) create mode 100644 examples/linear_eq.jl 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/src/linear_eq.jl b/src/linear_eq.jl index f86e5bcc..294faf13 100644 --- a/src/linear_eq.jl +++ b/src/linear_eq.jl @@ -13,8 +13,8 @@ end function gauss_seidel_interval(A::AbstractMatrix, b::AbstractArray; precondition=true, maxiter=100) n = size(A, 1) - x = fill(-1e16..1e16, n) - + x = similar(b) + x .= -1e16..1e16 gauss_seidel_interval!(x, A, b, precondition=precondition, maxiter=maxiter) return x end @@ -47,7 +47,8 @@ end function gauss_seidel_contractor(A::AbstractMatrix, b::AbstractArray; precondition=true, maxiter=100) n = size(A, 1) - x = fill(-1e16..1e16, n) + x = similar(b) + x .= -1e16..1e16 x = gauss_seidel_contractor!(x, A, b, precondition=precondition, maxiter=maxiter) return x end @@ -78,8 +79,8 @@ end function gauss_elimination_interval(A::AbstractMatrix, b::AbstractArray; precondition=true) n = size(A, 1) - x = fill(-1e16..1e16, n) - + x = similar(b) + x .= -1e16..1e16 x = gauss_elimination_interval!(x, A, b, precondition=precondition) return x From 22470c5dd9aaada80c8a219a6cad0d9f7fef7187 Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Fri, 1 Jun 2018 14:55:06 +0530 Subject: [PATCH 31/36] Documentation and minor fixes --- src/IntervalRootFinding.jl | 2 +- src/newton1d.jl | 54 +++++++++++++++++++------------------- 2 files changed, 28 insertions(+), 28 deletions(-) diff --git a/src/IntervalRootFinding.jl b/src/IntervalRootFinding.jl index da64ac5b..5f329f81 100644 --- a/src/IntervalRootFinding.jl +++ b/src/IntervalRootFinding.jl @@ -27,7 +27,7 @@ export export isunique, root_status -import IntervalArithmetic.interval +import IntervalArithmetic: interval, wideinterval diff --git a/src/newton1d.jl b/src/newton1d.jl index 85c643cb..9b3dea4a 100644 --- a/src/newton1d.jl +++ b/src/newton1d.jl @@ -7,16 +7,16 @@ 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, 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) - while !isempty(L) - X = pop!(L) + 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) @@ -32,10 +32,10 @@ function newton1d{T}(f::Function, f′::Function, x::Interval{T}; 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 + if X == X ∩ N # Checking if Newton step was redundant reps1 += 1 if reps1 > 20 reps1 = 0 @@ -44,36 +44,36 @@ function newton1d{T}(f::Function, f′::Function, x::Interval{T}; end X = X ∩ N - if (isempty(X) || diam(X) == 0) + if (isempty(X)) # No root in X break - elseif 0 ∈ f(interval(prevfloat(mid(X)), nextfloat(mid(X)))) + 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)) + while (n < 4 && (fa == 0 || fb == 0)) # Narrowing the interval further if fa == 0 - if 0 ∈ f(interval(prevfloat(X.lo), nextfloat(X.lo))) + if 0 ∈ f(wideinterval(X.lo)) fa = 1 else N = X.lo - (f(interval(X.lo)) / f′(X)) X = X ∩ N - if (isempty(X) || diam(X) == 0) + if (isempty(X)) root_exist = false break end end end if fb == 0 - if 0 ∈ f(interval(prevfloat(X.hi), nextfloat(X.hi))) + if 0 ∈ f(wideinterval(X.hi)) fb = 1 else - if 0 ∈ f(interval(prevfloat(mid(X)), nextfloat(mid(X)))) + 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) || diam(X) == 0) + if (isempty(X)) root_exist = false break end @@ -81,7 +81,7 @@ function newton1d{T}(f::Function, f′::Function, x::Interval{T}; end N = mid(X) - (f(interval(mid(X))) / f′(X)) X = X ∩ N - if (isempty(X) || diam(X) == 0) + if (isempty(X)) root_exist = false break end @@ -89,19 +89,19 @@ function newton1d{T}(f::Function, f′::Function, x::Interval{T}; end if root_exist push!(R, Root(X, :unique)) - debugroot && @show "Root found", X + debugroot && @show "Root found", X # Storing determined unique root end break end end - else - if diam(X) == initial_width + 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 "Repititive root found", X + debugroot && @show "Repeated root found", X reps = 0 continue end @@ -112,21 +112,21 @@ function newton1d{T}(f::Function, f′::Function, x::Interval{T}; 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(prevfloat(mid(X)), nextfloat(mid(X)))) + if 0 ∈ f(wideinterval(mid(X))) # if root in X, narrow interval further # 0 ∈ fⁱ(x) debug && println("0 ∈ fⁱ(x)") - if 0 ∉ f(interval(prevfloat(X.lo), nextfloat(X.lo))) + if 0 ∉ f(wideinterval(X.lo)) expansion_pt = X.lo - elseif 0 ∉ f(interval(prevfloat(X.hi), nextfloat(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(prevfloat(x1), nextfloat(x1))) || 0 ∉ f(interval(prevfloat(x2), nextfloat(x2))) + if 0 ∉ f(wideinterval(x1)) || 0 ∉ f(wideinterval(x2)) push!(L, interval(X.lo, m)) push!(L, interval(m, X.hi)) continue @@ -145,7 +145,7 @@ function newton1d{T}(f::Function, f′::Function, x::Interval{T}; debug && println("0 ∉ fⁱ(x)") - if (diam(X)/mag(X)) < reltol && diam(f(X)) < abstol + 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 @@ -189,7 +189,7 @@ function newton1d{T}(f::Function, f′::Function, x::Interval{T}; end end - push!(L, X) + push!(L, X) # Pushing X into L to be processed again end end From b406d3b95e3157e74517db3eeb940021f4f32e9b Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Thu, 7 Jun 2018 16:49:41 +0530 Subject: [PATCH 32/36] Add linear hull for solving system of linear equations --- perf/linear_eq.jl | 5 +++-- perf/linear_eq_tabular.txt | 20 ++++++++++++-------- src/IntervalRootFinding.jl | 2 +- src/linear_eq.jl | 29 +++++++++++++++++++++++++++++ test/linear_eq.jl | 2 +- 5 files changed, 46 insertions(+), 12 deletions(-) diff --git a/perf/linear_eq.jl b/perf/linear_eq.jl index a418bfba..de582343 100644 --- a/perf/linear_eq.jl +++ b/perf/linear_eq.jl @@ -43,13 +43,14 @@ end function benchmark_elimination(max=10) df = DataFrame() - df[:Method] = ["Gauss Elimination", "Base.\\"] + 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) - df[Symbol("$n")] = [t1, t2] + t3 = @belapsed linear_hull($A, $b) + df[Symbol("$n")] = [t1, t2, t3] end a = [] for i in 1:max diff --git a/perf/linear_eq_tabular.txt b/perf/linear_eq_tabular.txt index 69028385..f06bfc27 100644 --- a/perf/linear_eq_tabular.txt +++ b/perf/linear_eq_tabular.txt @@ -39,11 +39,15 @@ Gauss Seidel 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 │ +│ 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/src/IntervalRootFinding.jl b/src/IntervalRootFinding.jl index 5f329f81..12fc586d 100644 --- a/src/IntervalRootFinding.jl +++ b/src/IntervalRootFinding.jl @@ -19,7 +19,7 @@ export Root, is_unique, roots, find_roots, bisect, newton1d, slope, - slope_newton1d, + slope_newton1d, linear_hull, gauss_seidel_interval, gauss_seidel_interval!, gauss_seidel_contractor, gauss_seidel_contractor!, gauss_elimination_interval, gauss_elimination_interval! diff --git a/src/linear_eq.jl b/src/linear_eq.jl index 294faf13..93d6c6a9 100644 --- a/src/linear_eq.jl +++ b/src/linear_eq.jl @@ -143,3 +143,32 @@ function gauss_elimination_interval1!(x::AbstractArray, a::AbstractMatrix, b::Ab a \ b end + +function linear_hull(M::AbstractMatrix, r::AbstractArray) + + n = size(M, 1) + + ((M, r) = preconditioner(M, r)) + M_lo = inf.(M) + M_hi = sup.(M) + if all(.≤(M_lo, zero(M_lo))) + return M \ r + end + P = inv(M_lo) + 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] + end + end + H = interval.(min.(H1, H2), max.(H1, H2)) + + return H + else + return M \ r + end +end diff --git a/test/linear_eq.jl b/test/linear_eq.jl index f2972d4b..a7ee7775 100644 --- a/test/linear_eq.jl +++ b/test/linear_eq.jl @@ -9,7 +9,7 @@ using IntervalArithmetic, StaticArrays, IntervalRootFinding for i in 1:length(A) for precondition in (false, true) - for f in (gauss_seidel_interval, gauss_seidel_contractor, gauss_elimination_interval) + 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 From 8697a5f1554f3ef2d9c91e9c38856fc4fc525ef7 Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Fri, 8 Jun 2018 13:26:55 +0530 Subject: [PATCH 33/36] Add doc string --- src/linear_eq.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/linear_eq.jl b/src/linear_eq.jl index 93d6c6a9..503431be 100644 --- a/src/linear_eq.jl +++ b/src/linear_eq.jl @@ -144,6 +144,10 @@ function gauss_elimination_interval1!(x::AbstractArray, a::AbstractMatrix, b::Ab 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) From 375deabbe94aef838f78b53904996fadbfb64660 Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Sat, 9 Jun 2018 15:35:14 +0530 Subject: [PATCH 34/36] Add support for StaticArrays --- src/linear_eq.jl | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/linear_eq.jl b/src/linear_eq.jl index 503431be..c5ec54ac 100644 --- a/src/linear_eq.jl +++ b/src/linear_eq.jl @@ -166,7 +166,12 @@ function linear_hull(M::AbstractMatrix, r::AbstractArray) H2 = C .* Z for i in 1:n if Z[i] < 0 - H2[i] = Z[i] + #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)) From b6de4a397338e4dcea0c8ffea285a87c8c86c1b6 Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Thu, 7 Jun 2018 13:55:19 +0530 Subject: [PATCH 35/36] Initial implementation of multidim interval Newton method --- src/newtonnd.jl | 112 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 112 insertions(+) create mode 100644 src/newtonnd.jl diff --git a/src/newtonnd.jl b/src/newtonnd.jl new file mode 100644 index 00000000..5a153080 --- /dev/null +++ b/src/newtonnd.jl @@ -0,0 +1,112 @@ +""" +Preconditions the matrix A and b with the inverse of mid(A) +""" +function preconditioner(A::AbstractMatrix, b::AbstractArray) + + Aᶜ = mid.(A) + B = pinv(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ᴵ) + 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 + + use_B = false + next_iter = false + + initial_width = diam(Xᴵ) + debug && (print("Current interval popped: "); println(Xᴵ)) + if use_B + # TODO Compute X using B in Step 19 + else + X = mid(Xᴵ) + end + + J = SMatrix{n}{n}(deriv(Xᴵ)) # either jacobian or calculated using slopes + + # Xᴵ = IntervalBox((X + linear_hull(J, -f(interval.(X)))) .∩ Xᴵ) + Xᴵ = IntervalBox((X + (J \ -f(interval.(X)))) .∩ Xᴵ) + + 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...) From 3a52b0d09f3190c5ca572e04aa8834b25f58b84a Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Sun, 10 Jun 2018 18:27:06 +0530 Subject: [PATCH 36/36] Inner iteration to improve choice of point of expansion --- src/linear_eq.jl | 5 +++++ src/newtonnd.jl | 30 ++++++++++++++++++++++-------- 2 files changed, 27 insertions(+), 8 deletions(-) diff --git a/src/linear_eq.jl b/src/linear_eq.jl index c5ec54ac..70bbd5e4 100644 --- a/src/linear_eq.jl +++ b/src/linear_eq.jl @@ -153,12 +153,17 @@ 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) diff --git a/src/newtonnd.jl b/src/newtonnd.jl index 5a153080..7d892f52 100644 --- a/src/newtonnd.jl +++ b/src/newtonnd.jl @@ -4,7 +4,7 @@ Preconditions the matrix A and b with the inverse of mid(A) function preconditioner(A::AbstractMatrix, b::AbstractArray) Aᶜ = mid.(A) - B = pinv(Aᶜ) # TODO If Aᶜ is singular + B = inv(Aᶜ) # TODO If Aᶜ is singular return B*A, B*b @@ -23,7 +23,10 @@ function newtonnd(f::Function, deriv::Function, x::IntervalBox{NUM, T}; reltol=e if !all(0 .∈ f(Xᴵ)) continue end + Xᴵ¹ = copy(Xᴵ) + use_B = false + debug && (print("Current interval popped: "); println(Xᴵ)) if (isempty(Xᴵ)) @@ -45,22 +48,33 @@ function newtonnd(f::Function, deriv::Function, x::IntervalBox{NUM, T}; reltol=e while true - use_B = false next_iter = false initial_width = diam(Xᴵ) debug && (print("Current interval popped: "); println(Xᴵ)) + X = mid(Xᴵ) if use_B - # TODO Compute X using B in Step 19 - else - X = mid(Xᴵ) + 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 - - # Xᴵ = IntervalBox((X + linear_hull(J, -f(interval.(X)))) .∩ Xᴵ) + 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