Automatic evaluation of multidimensional slope expansions

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

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] │

src/slopes.jl

@@ -1,12 +1,24 @@
 # 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)
     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)

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