Updated README. Added STEPBYSTEP.md and VERSIONS.md. Fixed ^() defini…

…tions.

README.md

@@ -16,73 +16,40 @@ 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 Julia version was derived/written by Rob J Goedman ([email protected]).
+The initial Julia-ish version was derived/written by Rob J Goedman ([email protected]).
+
+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.]()
+
+For a quick into, see [STEPBYSTEP.md](https://github.com/JuliaDiff/HyperDualNumbers.jl/blob/master/STEPBYSTEP.md)
+
 Latest tagged versions:
+
 *  v1.1.0 (Julia 0.5 & 0.6, Oct 2017)
 *  v2.0.0 (Julia v0.7-, Oct 2017)
 *  v3.0.1 (Julia v0.7 & Julia v1.0), Aug 2018, Pkg(3))
+*  v4.0.0 (Julia v1.0, Nov 2018)
+
+For details see [VERSION.md](https://github.com/JuliaDiff/HyperDualNumbers.jl/blob/master/VERSIONS.md)
 
-The Julia package is structured similar to the [JuliaDiff/DualNumbers](https://github.com/JuliaDiff/DualNumbers.jl) package, which aims for complete support for `Dual` types for numerical functions in Julia. Currently, basic mathematical operations and trigonometric functions are supported.
 
 The following functions are specific to hyperdual numbers:
+
 * `Hyper`,
 * `Hyper256`,
 * `Hyper128`,
+* `ishyper`,
+* `hyper_show`
+* `realpart`,
+* `ε₁part()`, replaces eps1,
+* `ε₂part()`, replaces eps2,
+* `ε₁ε₂part()`, replaces eps1eps2
+
+In the future it is my intention to deprecate:
+
 * `hyper`,
 * `hyper256`,
 * `hyper128`,
-* `realpart`,
 * `eps1`,
 * `eps2`,
-* `eps1eps2`,
-* `ishyper`,
-* `hyper_show`
-
-Several other functions have been extended to accept hyperdual numbers, e.g.:
-`+`, ..., `<`, ..., `abs`, `log`, `sin`, ..., `erf`, `sqrt`, etc., see the final part of hyperdual.jl.
-
-[JuliaDiff](http://www.juliadiff.org) is a great starting point to learn about Julia packages related to Automatic Differentiation.
-
-### A walk-through example
-
-The example below demonstrates basic usage of hyperdual numbers by employing them to 
-perform automatic differentiation. The code for this example can be found in 
-`test/runtests.jl`.
-
-First install the package by using the Julia package manager:
-
-    Pkg.add("HyperDualNumbers")
-
-Then make the package available via
-
-    using HyperDualNumbers
-
-Use the `hyper()` function to define a hyperdual number, e.g.:
-
-    hd0 = hyper()
-    hd1 = hyper(1.0)
-    hd2 = hyper(3.0, 1.0, 1.0, 0.0)
-    hd3 = hyper(3//2, 1//1, 1//1,0//1)
-
-Let's say we want to calculate the first and second derivative of
-
-    f(x) = e^x / (sqrt(sin(x)^3 + cos(x)^3))
-
-To calculate these derivatives at a location `x`, evaluate your function at `hyper(x, 1.0, 1.0, 0.0)`. For example:
-
-    t0 = hyper(1.5, 1.0, 1.0, 0.0)
-    y = f(t0)
-
-For this example, you'll get the result
-
-    4.497780053946162 + 4.053427893898621ϵ1 + 4.053427893898621ϵ2 + 9.463073681596601ϵ1ϵ2
-
-The first term is the function value, the coefficients of both `ϵ1` and `ϵ2` (which correspond to the second and third arguments of `hyper`) are equal to the first derivative, and the coefficient of `ϵ1ϵ2` is the second derivative.
-
-You can extract these coefficients from the hyperdual number using the functions `realpart()`, `eps1()` or `eps2()` and `eps1eps2()`:
+* `eps1eps2`
 
-    println("f(1.5) = ", f(1.5))
-    println("f(t0) = ", realpart(f(t0)))
-    println("f'(t0) = ", eps1(f(t0)))
-    println("f'(t0) = ", eps2(f(t0)))
-    println("f''(t0) = ", eps1eps2(f(t0)))

STEPBYSTEP.md

@@ -0,0 +1,42 @@
+### A walk-through example
+
+The example below demonstrates basic usage of hyperdual numbers by employing them to perform automatic differentiation. The code for this example can be found in 
+`test/runtests.jl`.
+
+First install the package by using the Julia package manager:
+
+    Pkg.add("HyperDualNumbers")
+
+Then make the package available via
+
+    using HyperDualNumbers
+
+Use the `Hyper()` function to define a hyperdual number, e.g.:
+
+    hd0 = Hyper()
+    hd1 = Hyper(1.0)
+    hd2 = Hyper(3.0, 1.0, 1.0, 0.0)
+    hd3 = Hyper(3//2, 1//1, 1//1,0//1)
+
+Let's say we want to calculate the first and second derivative of
+
+    f(x) = ℯ^x / (sqrt(sin(x)^3 + cos(x)^3))
+
+To calculate these derivatives at a location `x`, evaluate your function at `hyper(x, 1.0, 1.0, 0.0)`. For example:
+
+    t0 = Hyper(1.5, 1.0, 1.0, 0.0)
+    y = f(t0)
+
+For this example, you'll get the result
+
+    4.497780053946162 + 4.053427893898621ϵ1 + 4.053427893898621ϵ2 + 9.463073681596601ϵ1ϵ2
+
+The first term is the function value, the coefficients of both `ϵ1` and `ϵ2` (which correspond to the second and third arguments of `hyper`) are equal to the first derivative, and the coefficient of `ϵ1ϵ2` is the second derivative.
+
+You can extract these coefficients from the hyperdual number using the functions `realpart()`, `ε₁part()` or `ε₂part()` and `ε₁ε₂part()`:
+
+    println("f(1.5) = ", f(1.5))
+    println("f(t0) = ", realpart(f(t0)))
+    println("f'(t0) = ", ε₁part(f(t0)))
+    println("f'(t0) = ", ε₂part(f(t0)))
+    println("f''(t0) = ", ε₁ε₂part(f(t0)))

VERSIONS.md

@@ -0,0 +1,27 @@
+### Versions
+
+#### v4.0.0
+
+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.)
+
+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.)
+
+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?)
+
+4. I renamed a few things to match the DualNumbers.jl nomenclature (e.g., realpart) and tried to keep aliases to older definitions (e.g., eps1eps2 is replaced by a clearer ε₁ε₂part to match the realpart construct, but eps1eps2 still available for backward compatibility)
+
+5. I changed a few tests from a == b to isequal(a, b) on purpose because == only compared the real part (h/h == 1 would give true even if the non-real parts are wrong for example).
+
+##### Outstanding questions/remarks:
+
+1. (Unsure, but) no `abs`, `real`, `conj`, `angle` or `imag`for hyperdual complex because differentiable nowhere in ℂ (these functions don't satisfy the Cauchy–Riemann equations). I guess it might be usable for printing or other checking purposes, but then I guess the user should define his own functions, at its own risk, rather than having this package suggest a value that does not make mathematical sense, right? (I may be completely wrong about this)
+
+2. No flipsign so far but maybe I should make it like the `abs` above and then have `abs(x) = flipsign(x, x)`?
+
+3. No `conjhyper`, `abshyper`, or `abs2hyper` because I don't understand their purpose in DualNumbers
+
+4. TODO: sinpi and cospi (erased here) should be generated in Calculus

src/hyperdual.jl

@@ -27,6 +27,8 @@ const Hyper64  = Hyper{Float16}
 const HyperComplex512  = Hyper{ComplexF64}
 const HyperComplex256 = Hyper{ComplexF32}
 const HyperComplex128 = Hyper{ComplexF16}
+
+# In the future to be deprecated:
 const hyper256 = Hyper256
 const hyper128 = Hyper128
 const hyper64 = Hyper64
@@ -347,11 +349,21 @@ function Base.:^(h::Hyper, a::Number)
         a*x^(a - 1)*z,
         a^2*x^(a - 2)*y*z - a*x^(a - 2)*y*z + a*w*x^(a - 1))
 end
+
+# Below definition is necesssaty to resolve a conflict with the
+# definition in MathConstants.jl
+function Base.:^(x::Irrational{:ℯ}, h::Hyper)
+    a, b, c, d = value(h), ε₁part(h), ε₂part(h), ε₁ε₂part(h)
+    return Hyper(x^a,
+        b*x^a*log(x),
+        c*x^a*log(x),
+        b*c*x^a*log(x)^2 + d*x^a*log(x))
+end
 function Base.:^(x::Number, h::Hyper)
     a, b, c, d = value(h), ε₁part(h), ε₂part(h), ε₁ε₂part(h)
     return Hyper(x^a,
         b*x^a*log(x),
-        c*x^a*log(x)*z,
+        c*x^a*log(x),
         b*c*x^a*log(x)^2 + d*x^a*log(x))
 end