Skip to content

Commit dd8156c

Browse files
eeshan9815eeshan9815
committed
Add slope function for multidim case with tests and benchmarks
1 parent 3cd778a commit dd8156c

File tree

4 files changed

+188
-35
lines changed

4 files changed

+188
-35
lines changed

‎perf/slopes.jl

Lines changed: 50 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -64,3 +64,53 @@ function benchmark_time()
6464
println(dfnew)
6565
dfnew
6666
end
67+
68+
struct SlopesMulti
69+
f::Function
70+
x::IntervalBox
71+
c::Vector
72+
sol::Matrix{Interval}
73+
end
74+
75+
function benchmark_multi()
76+
77+
rts = SlopesMulti[]
78+
f(x, y) = SVector(x^2 + y^2 - 1, y - 2x)
79+
f(X) = f(X...)
80+
X = (-6..6) × (-6..6)
81+
c = [0.0, 0.0]
82+
push!(rts, SlopesMulti(f, X, c, [-6..6 -6..6; -2.. -2 1..1]))
83+
84+
function g(x)
85+
(x1, x2, x3) = x
86+
SVector( x1^2 + x2^2 + x3^2 - 1,
87+
x1^2 + x3^2 - 0.25,
88+
x1^2 + x2^2 - 4x3
89+
)
90+
end
91+
92+
X = (-5..5)
93+
XX = IntervalBox(X, 3)
94+
cc = [0.0, 0.0, 0.0]
95+
push!(rts, SlopesMulti(g, XX, cc, [-5..5 -5..5 -5..5; -5..5 0..0 -5..5; -5..5 -5..5 -4.. -4]))
96+
97+
function h(x)
98+
(x1, x2, x3) = x
99+
SVector( x1 + x2^2 + x3^2 - 1,
100+
x1^2 + x3 - 0.25,
101+
x1^2 + x2 - 4x3
102+
)
103+
end
104+
105+
XXX = IntervalBox(1..5, 2..6, -3..7)
106+
ccc = [3.0, 4.0, 2.0]
107+
push!(rts, SlopesMulti(h, XXX, ccc, [1..1 6..10 -1..9; 4..8 0..0 1..1; 4..8 1..1 -4.. -4]))
108+
109+
for i in 1:length(rts)
110+
println("\nFunction $i")
111+
println("Slope Expansion: ")
112+
println(DataFrame(slope(rts[i].f, rts[i].x, rts[i].c)))
113+
println("\nJacobian: ")
114+
println(DataFrame(ForwardDiff.jacobian(rts[i].f, rts[i].x)))
115+
end
116+
end

‎perf/slopes_tabular.txt

Lines changed: 45 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -35,3 +35,48 @@ Time
3535
│ 7 │ f7 │ 3.0198e-5 │ 7.2014e-5 │
3636
│ 8 │ f8 │ 7.43125e-6 │ 8.046e-6 │
3737
│ 9 │ f9 │ 7.118e-6 │ 7.8705e-6 │
38+
39+
40+
Multi-Dimensional Case:
41+
Function 1
42+
Slope Expansion:
43+
│ Row │ x1 │ x2 │
44+
├─────┼──────────┼─────────┤
45+
│ 1 │ [-6, 6] │ [-6, 6] │
46+
│ 2 │ [-2, -2] │ [1, 1] │
47+
48+
Jacobian:
49+
│ Row │ x1 │ x2 │
50+
├─────┼───────────┼───────────┤
51+
│ 1 │ [-12, 12] │ [-12, 12] │
52+
│ 2 │ [-2, -2] │ [1, 1] │
53+
54+
Function 2
55+
Slope Expansion:
56+
│ Row │ x1 │ x2 │ x3 │
57+
├─────┼─────────┼─────────┼──────────┤
58+
│ 1 │ [-5, 5] │ [-5, 5] │ [-5, 5] │
59+
│ 2 │ [-5, 5] │ [0, 0] │ [-5, 5] │
60+
│ 3 │ [-5, 5] │ [-5, 5] │ [-4, -4] │
61+
62+
Jacobian:
63+
│ Row │ x1 │ x2 │ x3 │
64+
├─────┼───────────┼───────────┼───────────┤
65+
│ 1 │ [-10, 10] │ [-10, 10] │ [-10, 10] │
66+
│ 2 │ [-10, 10] │ [0, 0] │ [-10, 10] │
67+
│ 3 │ [-10, 10] │ [-10, 10] │ [-4, -4] │
68+
69+
Function 3
70+
Slope Expansion:
71+
│ Row │ x1 │ x2 │ x3 │
72+
├─────┼────────┼─────────┼──────────┤
73+
│ 1 │ [1, 1] │ [6, 10] │ [-1, 9] │
74+
│ 2 │ [4, 8] │ [0, 0] │ [1, 1] │
75+
│ 3 │ [4, 8] │ [1, 1] │ [-4, -4] │
76+
77+
Jacobian:
78+
│ Row │ x1 │ x2 │ x3 │
79+
├─────┼─────────┼─────────┼──────────┤
80+
│ 1 │ [1, 1] │ [4, 12] │ [-6, 14] │
81+
│ 2 │ [2, 10] │ [0, 0] │ [1, 1] │
82+
│ 3 │ [2, 10] │ [1, 1] │ [-4, -4] │

‎src/slopes.jl

Lines changed: 35 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,12 @@
11
# Reference : Dietmar Ratz - An Optimized Interval Slope Arithmetic and its Application
2-
using IntervalArithmetic, ForwardDiff
2+
using IntervalArithmetic, ForwardDiff, StaticArrays
33
import Base: +, -, *, /, ^, sqrt, exp, log, sin, cos, tan, asin, acos, atan
44
import IntervalArithmetic: mid, interval
55

6-
function slope(f::Function, x::Interval, c::Number)
6+
"""
7+
Expands the slope of `f` over `x` at the point `c` (default `c = mid(x)`)
8+
"""
9+
function slope(f::Function, x::Interval, c::Number = mid(x))
710
try
811
f(slope_var(x, c)).s
912
catch y
@@ -13,6 +16,25 @@ function slope(f::Function, x::Interval, c::Number)
1316
end
1417
end
1518

19+
function slope(f::Function, x::Union{IntervalBox, SVector}, c::AbstractVector = mid.(x)) # multidim
20+
try
21+
T = typeof(x[1].lo)
22+
n = length(x)
23+
s = Vector{Slope{T}}[]
24+
for i in 1:n
25+
arr = fill(Slope(zero(T)), n)
26+
arr[i] = slope_var(x[i], c[i])
27+
push!(s, arr)
28+
end
29+
return slope.(hcat(([(f(s[i])) for i=1:n])...))
30+
catch y
31+
if isa(y, MethodError)
32+
ForwardDiff.jacobian(f, x)
33+
end
34+
end
35+
36+
end
37+
1638
struct Slope{T}
1739
x::Interval{T}
1840
c::Interval{T}
@@ -76,12 +98,19 @@ end
7698

7799
+(v::Slope, u) = u + v
78100

79-
-(v::Slope, u) = u - v
80-
-(u::Slope) = u * -1.0
81-
82101
*(v::Slope, u) = u * v
83102

84-
/(v::Slope, u) = u / v
103+
function -(u::Slope, v)
104+
Slope(u.x - v, u.c - v, u.s)
105+
end
106+
107+
function -(u::Slope)
108+
Slope(-u.x, -u.c, -u.s)
109+
end
110+
111+
function /(u::Slope, v)
112+
Slope(u.x / v, u.c / v, u.s / v)
113+
end
85114

86115
function sqr(u::Slope)
87116
Slope(u.x ^ 2, u.c ^ 2, (u.x + u.c) * u.s)

‎test/slopes.jl

Lines changed: 58 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
using IntervalArithmetic, IntervalRootFinding
2-
using ForwardDiff
2+
using ForwardDiff, StaticArrays
33
using Base.Test
44

55
struct Slopes{T}
@@ -9,39 +9,68 @@ struct Slopes{T}
99
sol::Interval{T}
1010
end
1111

12-
@testset "Automatic slope expansion(Float64)" begin
12+
@testset "Automatic slope expansion" begin
13+
for T in [Float64, BigFloat]
14+
s = interval(T(0.75), T(1.75))
15+
rts = Slopes{T}[]
16+
push!(rts, Slopes(x->((x + sin(x)) * exp(-x^2)), s, mid(s), interval(T(-2.8), T(0.1))))
17+
push!(rts, Slopes(x->(x^4 - 10x^3 + 35x^2 - 50x + 24), s, mid(s), interval(T(-44), T(38.5))))
18+
push!(rts, Slopes(x->((log(x + 1.25) - 0.84x) ^ 2), s, mid(s), interval(T(-0.16), T(0.44))))
19+
push!(rts, Slopes(x->(0.02x^2 - 0.03exp(-(20(x - 0.875))^2)), s, mid(s), interval(T(0.03), T(0.33))))
20+
push!(rts, Slopes(x->(exp(x^2)), s, mid(s), interval(T(6.03), T(33.23))))
21+
push!(rts, Slopes(x->(x^4 - 12x^3 + 47x^2 - 60x - 20exp(-x)), s, mid(s), interval(T(-39), T(65.56))))
22+
push!(rts, Slopes(x->(x^6 - 15x^4 + 27x^2 + 250), s, mid(s), interval(T(-146.9), T(67.1))))
23+
push!(rts, Slopes(x->(atan(cos(tan(x)))), s, mid(s), interval(T(1), T(2))))
24+
push!(rts, Slopes(x->(asin(cos(acos(sin(x))))), s, mid(s), interval(T(1.36), T(∞))))
1325

14-
s = interval(0.75, 1.75)
15-
rts = Slopes{Float64}[]
16-
push!(rts, Slopes(x->((x + sin(x)) * exp(-x^2)), s, mid(s), interval(-2.8, 0.1)))
17-
push!(rts, Slopes(x->(x^4 - 10x^3 + 35x^2 - 50x + 24), s, mid(s), interval(-44, 38.5)))
18-
push!(rts, Slopes(x->((log(x + 1.25) - 0.84x) ^ 2), s, mid(s), interval(-0.16, 0.44)))
19-
push!(rts, Slopes(x->(0.02x^2 - 0.03exp(-(20(x - 0.875))^2)), s, mid(s), interval(0.03, 0.33)))
20-
push!(rts, Slopes(x->(exp(x^2)), s, mid(s), interval(6.03, 33.23)))
21-
push!(rts, Slopes(x->(x^4 - 12x^3 + 47x^2 - 60x - 20exp(-x)), s, mid(s), interval(-39, 65.56)))
22-
push!(rts, Slopes(x->(x^6 - 15x^4 + 27x^2 + 250), s, mid(s), interval(-146.9, 67.1)))
23-
push!(rts, Slopes(x->(atan(cos(tan(x)))), s, mid(s), interval(1, 2)))
24-
push!(rts, Slopes(x->(asin(cos(acos(sin(x))))), s, mid(s), interval(1.36, ∞)))
25-
26-
for i in 1:length(rts)
27-
@test slope(rts[i].f, rts[i].x, rts[i].c) ⊆ rts[i].sol
26+
for i in 1:length(rts)
27+
@test slope(rts[i].f, rts[i].x, rts[i].c) ⊆ rts[i].sol
28+
end
2829
end
2930
end
3031

31-
@testset "Automatic slope expansion(BigFloat)" begin
32-
s = interval(BigFloat(0.75), BigFloat(1.75))
33-
rts = Slopes{BigFloat}[]
34-
push!(rts, Slopes(x->((x + sin(x)) * exp(-x^2)), s, mid(s), interval(BigFloat(-2.8), BigFloat(0.1))))
35-
push!(rts, Slopes(x->(x^4 - 10x^3 + 35x^2 - 50x + 24), s, mid(s), interval(BigFloat(-44), BigFloat(38.5))))
36-
push!(rts, Slopes(x->((log(x + 1.25) - 0.84x) ^ 2), s, mid(s), interval(BigFloat(-0.16), BigFloat(0.44))))
37-
push!(rts, Slopes(x->(0.02x^2 - 0.03exp(-(20(x - 0.875))^2)), s, mid(s), interval(BigFloat(0.03), BigFloat(0.33))))
38-
push!(rts, Slopes(x->(exp(x^2)), s, mid(s), interval(BigFloat(6.03), BigFloat(33.23))))
39-
push!(rts, Slopes(x->(x^4 - 12x^3 + 47x^2 - 60x - 20exp(-x)), s, mid(s), interval(BigFloat(-39), BigFloat(65.56))))
40-
push!(rts, Slopes(x->(x^6 - 15x^4 + 27x^2 + 250), s, mid(s), interval(BigFloat(-146.9), BigFloat(67.1))))
41-
push!(rts, Slopes(x->(atan(cos(tan(x)))), s, mid(s), interval(BigFloat(1), BigFloat(2))))
42-
push!(rts, Slopes(x->(asin(cos(acos(sin(x))))), s, mid(s), interval(BigFloat(1.36), BigFloat(∞))))
32+
struct SlopesMulti
33+
f::Function
34+
x::IntervalBox
35+
c::Vector
36+
sol::Matrix{Interval}
37+
end
38+
39+
@testset "Multidim slope expansion" begin
40+
41+
rts = SlopesMulti[]
42+
f(x, y) = SVector(x^2 + y^2 - 1, y - 2x)
43+
f(X) = f(X...)
44+
X = (-6..6) × (-6..6)
45+
c = [0.0, 0.0]
46+
push!(rts, SlopesMulti(f, X, c, [-6..6 -6..6; -2.. -2 1..1]))
47+
48+
function g(x)
49+
(x1, x2, x3) = x
50+
SVector( x1^2 + x2^2 + x3^2 - 1,
51+
x1^2 + x3^2 - 0.25,
52+
x1^2 + x2^2 - 4x3
53+
)
54+
end
55+
56+
X = (-5..5)
57+
XX = IntervalBox(X, 3)
58+
cc = [0.0, 0.0, 0.0]
59+
push!(rts, SlopesMulti(g, XX, cc, [-5..5 -5..5 -5..5; -5..5 0..0 -5..5; -5..5 -5..5 -4.. -4]))
60+
function h(x)
61+
(x1, x2, x3) = x
62+
SVector( x1 + x2^2 + x3^2 - 1,
63+
x1^2 + x3 - 0.25,
64+
x1^2 + x2 - 4x3
65+
)
66+
end
67+
68+
XXX = IntervalBox(1..5, 2..6, -3..7)
69+
ccc = [3.0, 4.0, 2.0]
70+
push!(rts, SlopesMulti(h, XXX, ccc, [1..1 6..10 -1..9; 4..8 0..0 1..1; 4..8 1..1 -4.. -4]))
4371

4472
for i in 1:length(rts)
45-
@test slope(rts[i].f, rts[i].x, rts[i].c) ⊆ rts[i].sol
73+
@test slope(rts[i].f, rts[i].x, rts[i].c) == rts[i].sol
4674
end
75+
4776
end

0 commit comments

Comments
 (0)