From 7d0c09332eea8bff8f1e86a22369f77d40d0cf9b Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Fri, 1 Jun 2018 13:48:47 +0530 Subject: [PATCH] Add slope function for multidim case with tests and benchmarks --- perf/slopes.jl | 50 +++++++++++++++++++++++++++++++++++++++++ perf/slopes_tabular.txt | 45 +++++++++++++++++++++++++++++++++++++ src/slopes.jl | 14 +++++++++++- test/slopes.jl | 48 ++++++++++++++++++++++++++++++++++++++- 4 files changed, 155 insertions(+), 2 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 82119666..630d67f2 100644 --- a/src/slopes.jl +++ b/src/slopes.jl @@ -1,5 +1,5 @@ # 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 @@ -7,6 +7,18 @@ function slope(f::Function, x::Interval, c::Number) f(slope_var(x, c)).s end +function slope(f::Function, x::Union{IntervalBox, SVector}, c::AbstractVector = mid.(x)) # multidim + 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])...)) +end + struct Slope{T} x::Interval{T} # Interval on which slope is evaluated c::Interval{T} # Point about which slope is evaluated (Interval to get bounded rounding errors) diff --git a/test/slopes.jl b/test/slopes.jl index 8513899e..c296bbac 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} @@ -28,3 +28,49 @@ end end end end + +struct SlopesMulti + f::Function + x::IntervalBox + c::Vector + sol::Matrix{Interval} +end + +@testset "Multidim slope expansion" begin + + rts = SlopesMulti[] + f(x, y) = SVector(x^2 + y^2 - 1, y - 2x) + f(X) = f(X...) + X = (-6..6) × (-6..6) + c = [0.0, 0.0] + push!(rts, SlopesMulti(f, X, c, [-6..6 -6..6; -2.. -2 1..1])) + + function g(x) + (x1, x2, x3) = x + SVector( x1^2 + x2^2 + x3^2 - 1, + x1^2 + x3^2 - 0.25, + x1^2 + x2^2 - 4x3 + ) + end + + X = (-5..5) + XX = IntervalBox(X, 3) + cc = [0.0, 0.0, 0.0] + push!(rts, SlopesMulti(g, XX, cc, [-5..5 -5..5 -5..5; -5..5 0..0 -5..5; -5..5 -5..5 -4.. -4])) + function h(x) + (x1, x2, x3) = x + SVector( x1 + x2^2 + x3^2 - 1, + x1^2 + x3 - 0.25, + x1^2 + x2 - 4x3 + ) + end + + XXX = IntervalBox(1..5, 2..6, -3..7) + ccc = [3.0, 4.0, 2.0] + push!(rts, SlopesMulti(h, XXX, ccc, [1..1 6..10 -1..9; 4..8 0..0 1..1; 4..8 1..1 -4.. -4])) + + for i in 1:length(rts) + @test slope(rts[i].f, rts[i].x, rts[i].c) == rts[i].sol + end + +end