Skip to content

Commit

Permalink
Updated eps2, etc. Updated VERSIONS and README. Removed execution of …
Browse files Browse the repository at this point in the history
…test_power.jl
  • Loading branch information
goedman committed Nov 21, 2018
1 parent 55bc5e3 commit 72f8311
Show file tree
Hide file tree
Showing 7 changed files with 17 additions and 15 deletions.
2 changes: 1 addition & 1 deletion README.md
Expand Up @@ -16,7 +16,7 @@ Hyper-dual numbers can be used to compute first and second derivatives numerical

[The Development of Hyper-Dual Numbers for Exact Second Derivative Calculations](https://adl.stanford.edu/hyperdual/Fike_AIAA-2011-886.pdf)

The initial Julia-ish version was derived/written by Rob J Goedman (goedman@icloud.com).
The initial Julia-ish versions (up to v3.0.1) were derived/written by Rob J Goedman (goedman@icloud.com).

HyperDualNumbers.jl v4.0.0 has been completely redone by Benoit Pasquier to make it `Julia` and much better follows the structure of the [JuliaDiff/DualNumbers](https://github.com/JuliaDiff/DualNumbers.jl) package.]()

Expand Down
4 changes: 2 additions & 2 deletions VERSIONS.md
Expand Up @@ -6,9 +6,9 @@ A complete rewrite of the implementation by Benoit Pasquier.

##### Notes:

1. Redefined abs in such a way that the 2nd derivative of abs(-x^2) is correct. (Additionally, the previous abs(x) = sqrt(x^2) was quite slower — compare e.g., with Base.abs.)
1. v4.0.0 redefined abs in such a way that the 2nd derivative of abs(-x^2) is correct. (Additionally, the previous abs(x) = sqrt(x^2) was quite slower — compare e.g., with Base.abs.)

2. Did not include some functions (sometimes on purpose, sometimes because I did not know what to do with them). But I thought I should raise the issue that some functions should not be available for differentiation in ℂ because they are differentiable nowhere. Namely functions that inject into ℝ, like abs or angle, but also some quite common functions like conj. (The functions that are ℂ-differentiable should be in the code though.)
2. v4.0.0 did not include some functions (sometimes on purpose, sometimes because I did not know what to do with them). But I thought I should raise the issue that some functions should not be available for differentiation in ℂ because they are differentiable nowhere. Namely functions that inject into ℝ, like abs or angle, but also some quite common functions like conj. (The functions that are ℂ-differentiable should be in the code though.)

3. I chose to write my own list of first and second derivatives (using sagemath) to use for overloading most of the functions (like cos, sin, etc.). This was to match DualNumbers.jl's approach (which uses a list in Calculus.jl. (I saw that ForwardDiff.jl uses DiffRules.jl instead — so maybe this is were this list should go to be used by all other packages?)

Expand Down
2 changes: 1 addition & 1 deletion test/runtests.jl
Expand Up @@ -5,7 +5,7 @@ my_tests = [
"test_TimHoly_example.jl",
"test_Paper_example.jl",
"test_Erf_example.jl",
"test_power.jl",
#"test_power.jl",
"tests_for_coverage.jl",
"test_function_list.jl"
]
Expand Down
8 changes: 4 additions & 4 deletions test/test_Erf_example.jl
Expand Up @@ -35,9 +35,9 @@ S0 = Hyper(100.0, 1.0, 1.0, 0)
K=80.0;T=2.0;r=0.01;sigma=0.2;
Price=blsprice(S0,K,r,T,sigma);
@test (abs(realpart(Price)-blsprice(realpart(S0),K,r,T,sigma))<1e-15)
@test (abs(eps1(Price)-blsdelta(realpart(S0),K,r,T,sigma))<1e-15)
@test (abs(eps1eps2(Price)-blsgamma(realpart(S0),K,r,T,sigma))<1e-15)
@test (abs(ε₁part(Price)-blsdelta(realpart(S0),K,r,T,sigma))<1e-15)
@test (abs(ε₁ε₂part(Price)-blsgamma(realpart(S0),K,r,T,sigma))<1e-15)
println("$(realpart(Price)) = $( blsprice(realpart(S0),K,r,T,sigma))")
println("$(eps1(Price)) = $( blsdelta(realpart(S0),K,r,T,sigma))")
println("$(eps1eps2(Price)) = $( blsgamma(realpart(S0),K,r,T,sigma))")
println("$(ε₁part(Price)) = $( blsdelta(realpart(S0),K,r,T,sigma))")
println("$(ε₁ε₂part(Price)) = $( blsgamma(realpart(S0),K,r,T,sigma))")
println("Test Erf Passed")
4 changes: 2 additions & 2 deletions test/test_Paper_example.jl
Expand Up @@ -29,6 +29,6 @@ println()
println("f(t0) = ", f(t0))

@test abs(f(1.5) - realpart(f(t0))) < 5eps(1.0)
@test abs(eps1(f(t0))-4.053427893898621) < 5eps(1.0) && eps1(f(t0)) == eps2(f(t0))
@test abs(eps1eps2(f(t0))-9.463073681596601) < 5eps(1.0)
@test abs(ε₁part(f(t0))-4.053427893898621) < 5eps(1.0) && ε₁part(f(t0)) == ε₂part(f(t0))
@test abs(ε₁ε₂part(f(t0))-9.463073681596601) < 5eps(1.0)

8 changes: 4 additions & 4 deletions test/test_basics.jl
Expand Up @@ -16,7 +16,7 @@ hd9 = hyper(-1.0, -2.0, -3.0, -4.0)
hdNaN = hyper(0/0)

# Addition and subtraction
@test eps1(hd1) == 0.0
@test ε₁part(hd1) == 0.0
@test isequal(hd1 + hd2, hyper(2.0, 2.0, 3.0, 4.0))
@test isequal(-hd2, hyper(-1.0, -2.0, -3.0, -4.0))
@test isequal(+hd3, hd2 + hyper(0.0, 1.0, 0.0, 0.0))
Expand All @@ -29,15 +29,15 @@ hdNaN = hyper(0/0)
# NaN tests
@test isnan(hdNaN) == isnan(hyper(NaN, 0.0, 0.0, 0.0))
hdNaN128 = hyper128(NaN, 4.0, 3.0, 4.0)
@test isnan(hdNaN + hd4) == isnan(hdNaN128) && eps1eps2(hdNaN128) == 4.0
@test isnan(hdNaN + hd4) == isnan(hdNaN128) && ε₁ε₂part(hdNaN128) == 4.0

println("\nExamples of show() for hyperdual numbers with NaN:\n")
println("hdNaN = $(hdNaN)")
println("hdNaN + hd4 = $(hdNaN + hd4)")

@test isnan(hdNaN) == true
@test isnan(hdNaN+hd4) == true
@test eps1eps2(hdNaN+hd4) == 4.0
@test ε₁ε₂part(hdNaN+hd4) == 4.0

# Using and mixing 64 & 32 bits
@test isequal(hd6+hd7, hyper(2.0, 11.0, 6.0, 8.0))
Expand All @@ -61,7 +61,7 @@ println("\nTesting includes Tim Holy's division performance improvement.")
println("Testing includes ngedwin98's fixes for asin, acos and atan.")
@test isequal(asin(sin(hd8)), hd8)
@test isequal(atan(tan(hd8)), hd8)
@test eps1(acos(hd8)) == -eps1(asin(hd8)) && eps1eps2(acos(hd8)) == -eps1eps2(asin(hd8))
@test ε₁part(acos(hd8)) == -ε₁part(asin(hd8)) && ε₁ε₂part(acos(hd8)) == -ε₁ε₂part(asin(hd8))

# Mixing types
@test isequal(hd5*hd7, hd7^2)
Expand Down
4 changes: 3 additions & 1 deletion test/test_power.jl
Expand Up @@ -12,6 +12,7 @@ r1 = f3(t0)
@test round(ε₂part(r1), digits=4) == 790.7755
@test round(ε₁ε₂part(r1), digits=10) == 6883.7763738258
r1 |> display
println()

# test 2
f4(x) = -4x^2 + 18
Expand All @@ -22,4 +23,5 @@ r2 = f6(t0)
@test round(ε₁part(r2), digits=4) == 0.6355
@test round(ε₂part(r2), digits=4) == 0.6355
@test round(ε₁ε₂part(r2), digits=10) == -503.1625169931
r2 |> display
r2 |> display
println()

0 comments on commit 72f8311

Please sign in to comment.