juliadiff/dualnumbers.jl

Julia package for representing dual numbers and for performing dual algebra

dualnumbers.jl's Issues

/

/

bitmix not defined

DualNumbers.jl/src/dual.jl

Line 178 in 4603cc1

Base.hash(z::Dual) = (x = hash(value(z)); epsilon(z)==0 ? x : bitmix(x,hash(epsilon(z))))

0.4: @compat value.(v) for v::Vector{Dual{Float64}} returns Vector{Any}

This code is returning the wrong type, and without the @compat and . a deprecated note (on master DualNumbers):

julia> using DualNumbers, Compat

julia> v=[dual(1.,2.)]
1-element Array{DualNumbers.Dual{Float64},1}:
 1.0+2.0ɛ

julia> @compat value.(v)
1-element Array{Any,1}:
 1.0

ambiguity warnings

Version 0.3.0-prerelease+1290 (2014-01-28 23:09 UTC)
 Commit 1234340* (1 days old master)
 x86_64-w64-mingw32

        julia> using DualNumbers
Warning: New definition
    *(Real,Dual{T<:Real}) at \Users\Kees\.julia\DualNumbers\src\dual.jl:128
is ambiguous with:
    *(Bool,T<:Number) at bool.jl:52.
To fix, define
    *(Bool,_<:Dual{T<:Real})
before the new definition.
Warning: New definition
    *(Dual{T<:Real},Real) at \Users\Kees\.julia\DualNumbers\src\dual.jl:129
is ambiguous with:
    *(Number,Bool) at bool.jl:55.
To fix, define
    *(Dual{T<:Real},Bool)
before the new definition.

Support for fft of DualComplex Numbers

using DualNumbers;
x=[DualComplex256(1,1) for i in 1:10]  
fft(x)

I get the following error:

ERROR: MethodError: no method matching plan_fft(::Array{DualNumbers.Dual{Complex
{Float64}},1}, ::UnitRange{Int64})
Closest candidates are:
  plan_fft(::Union{Base.ReshapedArray{T<:Union{Complex{Float32}, Complex{Float64
}},N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicati
veInverse{Int64},N} where N} where A<:Union{DenseArray, SubArray{T,N,P,I,true} w
here I<:Tuple{Union{Base.Slice, UnitRange},Vararg{Any,N} where N} where P where
N where T}, DenseArray{T<:Union{Complex{Float32}, Complex{Float64}},N}, SubArray
{T<:Union{Complex{Float32}, Complex{Float64}},N,A,I,L} where L} where I<:Tuple{V
ararg{Union{Base.AbstractCartesianIndex, Int64, Range{Int64}},N} where N} where
A<:Union{Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.Multiplicative
Inverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{DenseArra
y, SubArray{T,N,P,I,true} where I<:Tuple{Union{Base.Slice, UnitRange},Vararg{Any
,N} where N} where P where N where T} where N where T, DenseArray}, ::Any; flags
, timelimit) where {T<:Union{Complex{Float32}, Complex{Float64}}, N} at fft\FFTW
.jl:596
  plan_fft(::AbstractArray{#s262,N} where N where #s262<:Real, ::Any; kws...) at
 dft.jl:205
  plan_fft(::AbstractArray{#s250,N} where N where #s250<:(Complex{#s249} where #
s249<:Union{Integer, Rational}), ::Any; kws...) at dft.jl:207
  ...
Stacktrace:
 [1] #plan_fft#1(::Array{Any,1}, ::Function, ::Array{DualNumbers.Dual{Complex{Fl
oat64}},1}) at .\dft.jl:58
 [2] fft(::Array{DualNumbers.Dual{Complex{Float64}},1}) at .\dft.jl:56

Consider replacing SpecialFunctions.jl by pure Julia alternatives

I've noticed that DualNumbers.jl checks if the function is available in SpecialFunctions.jl to wrap the definition with the Dual type. More recently, initiatives such as Bessels.jl have gained traction. They are pure Julia implementations of special functions that are faster and differentiable. What are your thoughts on replacing SpecialFunctions.jl by Bessels.jl and similar efforts to clean up the dependency tree downstream?

DualNumbers.jl is an indirect dependency of Distributions.jl which is pretty much everywhere in the ecosystem. So installing Distributions.jl nowadays comes with the downside of installing external _jll libraries and a ton of dependencies that could be avoided with Bessels.jl and friends.

I've already migrated all my packages from SpecialFunctions.jl to Bessels.jl because they provide all the bessel functions and the gamma (non-exported) functions I need without a single dependency.

cc: @heltonmc @cgeoga @oscardssmith

Consider renaming realpart to value or primalpart or similar

This package works a treat for evaluating the derivative functions at complex locations:

julia> using DualNumbers

julia> f(x) = (1, -2 *x) .* exp(-x^2) # value and derivative
f (generic function with 1 method)

julia> z = 1.0 + im # some complex location to evaluate the derivative
1.0 + 1.0im

julia> f(z)
(-0.4161468365471424 - 0.9092974268256817im, -0.9863011805570786 + 2.650888526745648im)

julia> realpart(f(Dual(z, 1))[1])
-0.4161468365471424 - 0.9092974268256817im

julia> dualpart(f(Dual(z, 1))[1])
-0.9863011805570786 + 2.650888526745648im

In this instance, realpart returns a complex value, which is a little misleading. @goretkin suggested in Slack that primalpart could be a better adjective. Any thoughts?

nan(Dual{T}) has type T

 nan(Dual{Float32})
 > NaN32

But a Dual NaN does exist, kinda:

convert(Dual{Float32},NaN)
>NaN32 + 0.0f0du

The definition of nan(f)

Base.nan(f)

Returns NaN (not-a-number) of the floating point type "f" or of
the same floating point type as "f"

There are two options for what could be done.

Either nan{Dual{T}} should a NaN in the real part and zero in the du part -- it makes sense to match the behavior of convert, and for the nan operator to apply to the real part.

Or it should do as nan(Complex{Float32}) does and throw an error:

`nan` has no method matching nan(::Type{Complex{Float32}}) while loading In[1], in expression starting on line 1

My code assumes that nan{T<:Number} gives me a NaN of type T.
Or at very least makes it clear that it can't be done.

abs(Dual(-0.0,1.0))

Perhaps we should have
abs(Dual(-0.0,1.0)) == Dual(0.0,-1.0)?

julia> abs(Dual(0.0,1.0))
0.0 + 1.0du

julia> abs(Dual(-0.0,1.0))
-0.0 + 1.0du

julia> abs(Dual(-eps(),1.0))
2.220446049250313e-16 - 1.0du

div() and rem()

now div() and rem() are not supported:

julia> x = Dual(1.1, 2.2);

julia> y = Dual(3.3, 4.4);

julia> div(y, x)
ERROR: MethodError: no method matching div(::Dual128, ::Dual128, ::RoundingMode{:ToZero})

julia> rem(y, x)
ERROR: MethodError: no method matching rem(::Dual128, ::Dual128)

how could we handle these functions?
thanks

Sqrt at 0.0

julia> sqrt(Dual(0.0))
dual(0.0,NaN)

So this surprised me... whats the intution behind it?

enable travis

@scidom, I've pushed a basic travis configuration file. Could you enable travis? (You'll need to log into https://travis-ci.org/).

erf function broken in julia Nightly

Hi,
i am getting the following error when running some test on the latest julia nightly:
Error During Test
Test threw an exception of type MethodError
MethodError: no method matching erf(::DualNumbers.Dual{Float64})
You may have intended to import Base.erf
Closest candidates are:
erf(::Complex{Float32}) at /home/travis/.julia/v0.7/SpecialFunctions/src/erf.jl:26
erf(::Complex{Float64}) at /home/travis/.julia/v0.7/SpecialFunctions/src/erf.jl:25
erf(::BigFloat) at /home/travis/.julia/v0.7/SpecialFunctions/src/erf.jl:14

the julia version is the following:
Julia Version 0.7.0-DEV.1063
on the previous version of julia-nightly i had no problem

Printing method calls "other" component `du`, but its accessed through `epsilon`

So is the other component the dual component, or the epsilon component?

Should we export ϵ ≡ ε?

zeta function have no method with Dual type？

x = dual(2, 1)
y = zeta(x)

will have a error: 'zeta' has no method matching zeta(::Dual{Int64})

Calculus.symbolic_derivative_1arg_list masking some definitions

I might be confused, but

julia> @which abs2(Dual(1.0,1.0))
abs2(z::DualNumbers.Dual{T<:Real}) at /Users/goretkin/.julia/v0.4/DualNumbers/src/dual.jl:186

https://github.com/JuliaDiff/DualNumbers.jl/blob/v0.1.5/src/dual.jl#L186
which is masking a previous definition on
https://github.com/JuliaDiff/DualNumbers.jl/blob/v0.1.5/src/dual.jl#L122

Add norm for Vector{T<:Dual}

julia> norm([dual(1.,2.),dual(3.,2.)])
ERROR: InexactError()
 in float at float.jl:121
 in generic_vecnormInf at linalg/generic.jl:92
 in generic_vecnorm2 at linalg/generic.jl:114
 in vecnorm at linalg/generic.jl:175
 in norm at linalg/generic.jl:191

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

isless broken for mixed types

julia> dual(1) < dual(2.0)
ERROR: MethodError: no method matching isless(::Dual{Int64}, ::Dual{Float64})
Closest candidates are:
  isless(::Missing, ::Any) at missing.jl:66
  isless(::InfiniteArrays.OrientedInfinity{Bool}, ::Number) at /Users/sheehanolver/.julia/packages/InfiniteArrays/Z4yap/src/Infinity.jl:145
  isless(::Dual{T<:Real}, ::Dual{T<:Real}) where T<:Real at /Users/sheehanolver/.julia/packages/DualNumbers/QeItX/src/dual.jl:174
  ...
Stacktrace:
 [1] <(::Dual{Int64}, ::Dual{Float64}) at ./operators.jl:260
 [2] top-level scope at none:0

julia> dual(1.0) < dual(2.0)
true

isnan(Dual(1.0,NaN)) == false

julia> isnan(Dual(1.,NaN))
false

This didn't come up as an issue for me, but I'm curious if this is a deliberate choice.

compare

julia> isnan(1. + (NaN)im)
true

isapprox throws an error

I think real(::Type{Dual{T}}) needs to be overriden:

julia> dual(1.0,1.0) ≈ dual(1.0,1.0)
ERROR: MethodError: no method matching real(::Type{DualNumbers.Dual{Float64}})
Closest candidates are:
  real(::Complex{T<:Real}) at complex.jl:36
  real(::Real) at complex.jl:38
  real{T<:Real}(::Type{T<:Real}) at complex.jl:42
  ...
 in isapprox(::DualNumbers.Dual{Float64}, ::DualNumbers.Dual{Float64}) at ./floatfuncs.jl:251

`Dual <: Real` would be more practical

I'm aware that DualNumbers are not mathematically a kind of real number; they are more like first-degree approximations of hyperreal numbers (or surreal numbers). But in practice, they are useful in places where a Real number is expected. If JuliaStats/Distributions.jl#475 is going to become the new form of Distributions, parameterized by Real, then we definitely want to be able to find the log density gradients using dualNumbers, which means that Dual should be a subtype of Real, not of Number.

(This seems to be the general consensus at JuliaStats/Distributions.jl#430 , by the way; but I came to this issue separately, as somebody who hopes to a GSoC project on Mamba.)

Are Dual4 and GraDual trying to do the same thing?

@mlubin: I saw the work on Dual4 you started, and that looks VERY exciting to me. I assume you eventually want to use tuples (i.e. fixed size arrays) to store the epsilons, so that this would essentially be a DualN datatype, right?

I am just wondering how that relates to the GraDual type in ForwardDiff.jl? Is that essentially exactly the same thing? Except for the implementation details. Right now GraDual uses an array to store the epsilons, which probably makes it very slow. But once fixed size arrays are part of core julia, one would just have to change the implementation details and probably make this really fast?

I guess this issue is just to clarify whether these two efforts could be unified. In some way it seems to me that this kind of functionality actually fits better into the ForwardDiff.jl package, but not sure.

Also pinging @scidom.

@show imε output

I noticed, in Julia 1.0.0, the following:

julia> @show imε;
imε = imɛ

The output produced by the @show macro uses 2 different characters for the epsilon: first, U+03b5 and then the correct U+02b5 (in some fonts the difference between ε and ɛ is very noticeable):

julia> for i in "imε = imɛ"; display(i); end
'i': ASCII/Unicode U+0069 (category Ll: Letter, lowercase)
'm': ASCII/Unicode U+006d (category Ll: Letter, lowercase)
'ε': Unicode U+03b5 (category Ll: Letter, lowercase)
' ': ASCII/Unicode U+0020 (category Zs: Separator, space)
'=': ASCII/Unicode U+003d (category Sm: Symbol, math)
' ': ASCII/Unicode U+0020 (category Zs: Separator, space)
'i': ASCII/Unicode U+0069 (category Ll: Letter, lowercase)
'm': ASCII/Unicode U+006d (category Ll: Letter, lowercase)
'ɛ': Unicode U+025b (category Ll: Letter, lowercase)

Any ideas of why this is so and how to fix it?

max when values agree?

The partial derivative of max currently uses the Julia Base implementation:

max(x,y) = ifelse(y < x, x, y)

Unfortunately, this may cause incorrect results when the Dual values agree but the partials do not:

max(Dual(3,1), Dual(3,-1))

Here, the values are both 3 but the partials have opposite signs. I think the partial here may be undefined.

Maybe something like:

function Base.max(a::Dual, b::Dual)
    if value(a) > value(b)
        return a
    elseif value(a) < value(b)
        return b
    else # value(a) == value(b)
        return ifelse(epsilon(a) == epsilon(b), a, Dual(value(a), NaN))
    end
end

Tag a new version

The last tag was in Feb, and I need a new tag to tag SingularIntegralEquations.jl

Complex Differentiation

I'd really like to be able to differentiate complex functions with ForwardDiff. For almost all analytic functions, the derivative can be taken in the same way as in the real case, so Julia's dynamic type system should make complex diferentiation work out of the box. It seems to me that the only reason it doesn't is because many functions in ForwardDiff accept Real arguments, but if we simply changed these to accept Numbers then we should be most of the way there. Functions that are not analytic such as real, imag, abs, conj, abs2, reim should all return errors when their derivatives are taken at a complex input. For example, consider these two functions:

f(z) = z^2
g(z) = Complex(real(z)^2 - imag(z)^2, 2 * real(z) * imag(z))

These functions are equal and both analytic, but it is not obvious by looking at g that it is analytic. In my opinion, it would be acceptable if derivative(f, im) worked while derivative(g, im) threw an error.

I already started experimenting with this. I copied the definitions of several functions in ForwardDiff and just replaced Real with Complex, like so:

ForwardDiff.can_dual(::Type{Complex{T}}) where {T<:Real} = true

@inline function ForwardDiff.derivative(f::F, x::R) where {F,R<:Complex}
  T = typeof(Tag(f, R))
  return extract_derivative(T, f(Dual{T}(x, one(x))))
end
@inline function Base.:*(partials::Partials, x::Complex)
  return Partials(scale_tuple(partials.values, x))
end
@inline function Base.:/(partials::Partials, x::Complex)
  return Partials(div_tuple_by_scalar(partials.values, x))
end
@inline Base.:*(x::Complex, partials::Partials) = partials * x
@inline Base.:*(partials::Partials{0,V}, x::Complex) where {V} = Partials{0,promote_type(V, typeof(x))}(tuple())
@inline Base.:*(x::Complex, partials::Partials{0,V}) where {V} = Partials{0,promote_type(V, typeof(x))}(tuple())
@inline Base.:/(partials::Partials{0,V}, x::Complex) where {V} = Partials{0,promote_type(V, typeof(x))}(tuple())

This was enough for me to correctly differentiate sin. I expect that if we just do this throughout ForwardDiff (and replace the Real definitions with Number, instead of adding Complex definitions like I did here), then we'll be most of the way there. Here's my plan:

1. Find and replace Real with Number throughout the repository
2. Tweak things until all the test cases pass
3. Add test cases for complex functions

I'm tempted to actually do this and make a pull request as soon as I have some time to do so.

provide == as well as isequal for Julia 0.3

Currently, you are providing an isequal method, but not ==, for various types, e.g. you have isequal(z::Dual, x::Real) but not ==(z::Dual, x::Real). In Julia 0.3, you will also need to provide ==. See JuliaLang/julia#6833.

(In Julia 0.2, == called isequal by default, so providing isequal was enough. In Julia 0.3, however, these two functions are swapped: isequal calls == by default. So overriding isequal is not enough to make things like != or == work. In order to remain backward-compatible with Julia 0.2, you will need to provide both methods.)

split dual AD and mathematical types

There seems to be some conflict over whether we should follow the mathematical convention for dual numbers or use the definitions that allow dual numbers to be dropped in automatically to compute derivatives. I propose separating these into two separate types to settle the conflict.

As the primary use of DualNumbers, as far as I'm aware, is for computing derivatives, I'd argue that the AD type should keep the name Dual. This also avoids breakage downstream. Perhaps DualAlg could be used for the algebraic/mathematical type?

@scidom @andreasnoackjensen

Argument of a Dual Number

The following method always seems to return 0 + 0ɛ for positive duals.

angle{T<:Real}(z::DualNumbers.Dual{T})
angle{T<:Real}(z::Dual{T}) = z ≥ 0 ? zero(z) : one(z)*π

This is a lot different than the definition of argument that I'm familiar with. In the generalized complex numbers, parameterized by p and q:
z = x + iy (x, y ∈ R) where i^2 =iq + p (q, p ∈ R)

(so i, ɛ, λ, whatever we call the imaginary unit)

Ignoring q and setting it to zero, and setting p = 0 we get the Dual numbers, p = 1 gives Double numbers, and p = -1 giving the regular Complex numbers. There is a general definition for angle or argument, magnitude or norm, and a very weird concept of unit circle, that is only actually what you would call a circle for p = -1.

So the unit "circle" is everywhere ||z|| = 1, which for duals is just at 1 and -1.

So for duals the argument is y / x. I will link to a paper, and just paste an image of the relevant part here.

Here is a link to the paper this is from:
https://people.rit.edu/harkin/research/articles/generalized_complex_numbers.pdf

The paper also shows how to implement the sinp, cosp, and tanp trig functions for any value of p.

don't use imag() for second component?

@johnmyleswhite suggested and I agree that it's awkward to use imag to access the second component of the dual number. What about epsilon?

Package compatibility caps

Ref: https://discourse.julialang.org/t/package-compatibility-caps/15301

Dual numbers and Complex numbers

Many functions of a complex argument are not differentiable, and there's some math here that I don't quite understand, but if I define the DFT as

function dft{T<:Number}(x::Vector{T})
    N = length(x)

    #T should not have any imaginary part 

    C = zeros(T,N)
    S = zeros(T,N)

    for k=0:(N-1)
        for n=0:(N-1)
            cs = exp(-im*2pi* k*n/N)
            c,s = real(cs), imag(cs)
            @inbounds C[k+1] += c*x[n+1]
            @inbounds S[k+1] += s*x[n+1]
        end
    end
    return C,S 
end

which just returns the real and imaginary parts of the transform separately,

and also define some function that maps, for example, a 1D time-series x to a Real

function spec(x)
    C,S = dft(x)
    h = floor(length(x) / 2)
    w = cat(1, [1.0:h ], zeros(length(x)-int(h)) )
    sum( (C.^2 + S.^2) .* w ) #some weighting on different frequencies
end

then I can take the derivative of spec with respect to the input timeseries

function grad{T<:Real}(fun::Function,x::Vector{T})
    N = length(x)
    G = zeros(T,N)
    for i = 1:N
        xd = [Dual(x[j],1.0*(j==i)) for j=1:N]
        G[i] = epsilon(fun(xd))
    end
    return G
end

julia> grad(spec,rand(10,))
10-element Array{Float64,1}:
 12.1547 
 19.9224 
 -3.22121
 18.8699 
 21.0265 
 -9.73494
  1.77734
 18.524  
 24.9879 
  9.14384

hopefully that example made sense. If I wanted to use the complex version of the DFT and suitably modify spec ((C.^2 + S.^2) would be taking the squared-magnitude of those complex numbers), then as far as Real numbers go, spec would return the same output for the same inputs. But I don't see how I would put a Dual number through that function.

Use dualpart to determine value on branch cut?

I'm proposing that functions with branch cuts be overriden to use the dualpart to determine direction of approach, that is realpart(f(dual(a,b))) should always return lim_h->0 f(a+h*b).

For example,

    f(z)=sqrt(z)
    realpart(f(dual(-1.,im)))     # correctly returns 1.im
    realpart(f(dual(-1.,-im)))     # should return -1.im but returns 1.im

While one can accomplish this example using -1.+0.im and -1.-0.im,this breaks if the branch cut is rotated, where using dual numbers to determine the angle could still work:

    f(z)=sqrt((1.+1.im)*z)
    z=-1/(1.+1.im)
    realpart(f(dual(z,-1.im)))  # would return -1im which is lim_{h->0} f(z-h*im)

Dual division rules

Hi, nice to see a dual number library in Julia. I have some issues for your consideration.

Dual(0, 2) / Dual(0, 6)
NaN + NaN*ɛ

When both terms are purely dual, and the denominator is not 0 + 0ε causing a division by zero, I believe this is the correct way to perform division:

Solving for x, y in:

(0 + bε) / (0 + dε) = (x + yε)
(0 + bε) = (x + yε) * (0 + dε)
(0 + bε) = x * 0 + xdε + yε * 0 + ydε^2
(0 + bε) = xdε
x = bε/dε
x = b/d
y = anything, it is annihilated. 

Because of this, in my implementation of Dual numbers dividing by two purely dual numbers results in a real, and I set the dual part to 0ε, (but it could be anything).

By doing this I am able to perform a division, where the result multiplied by the denominator gives you back the numerator, and result multiplied by the numerator gives you back the denominator, keeping the inverse relationship between multiplication and division consistent.

(0 + 2ε) /  (0 + 6ε) = (1/3 + 0ε)
(1/3 + 0ε) * (0 + 6ε) = (0 + 2ε)
(1/3 + 0ε) * (0 + 2ε) = (0 + 6ε)

I realize it makes no sense for an operation on 2 infinitesimal numbers to result in a real, and the selection of 0 is arbitrary since it has an infinite amount of solutions, but when writing my implementation I searched around on math exchange etc, and found this method as accepted. You can also think of this as the ε cancelling out during the division, resulting in a real.

What do you think of this?

I have a few more suggestions, I will make separate issues for. Thanks :)

Here is my implementation, albeit in C++

  //// Division
  Dual<T> operator/=(const Dual<T> &other){
    if(_real == 0 && other._real == 0 && other._dual != 0){
      _real = _dual / other._dual;
      _dual = 0;
    }else{
      //  Otherwise we divide the same way as a regular complex number
      //  d * d.conjugate() is always Real(d)^2
      (*this) *= other.conjugate();
      auto reciprocal = 1 / (other._real * other._real);
      _real *= reciprocal;
      _dual *= reciprocal;
    }
    return *this;
  }

  Dual<T> operator/(const Dual<T> &other) const {
    auto copy = *this;
    return copy /= other;
  }

A bug of the ^ operator

Operation

Dual(1, 0) ^ Dual(-3, 0)

gives error:

DomainError:
Cannot raise an integer x to a negative power -n.
Make x a float by adding a zero decimal (e.g. 2.0^-n instead of 2^-n), or write 1/x^n, float(x)^-n, or (x//1)^-n.
in power_by_squaring at /Applications/Julia-0.4.5.app/Contents/Resources/julia/lib/julia/sys.dylib
in ^ at /Users/kai/.julia/v0.4/DualNumbers/src/dual.jl:235
in include_string at /Users/kai/.julia/v0.4/CodeTools/src/eval.jl:28
in include_string at /Users/kai/.julia/v0.4/CodeTools/src/eval.jl:32
[inlined code] from /Users/kai/.julia/v0.4/Atom/src/eval.jl:39
in anonymous at /Users/kai/.julia/v0.4/Atom/src/eval.jl:62
in withpath at /Users/kai/.julia/v0.4/Requires/src/require.jl:37
in withpath at /Users/kai/.julia/v0.4/Atom/src/eval.jl:53
[inlined code] from /Users/kai/.julia/v0.4/Atom/src/eval.jl:61
in anonymous at task.jl:58

Is this a bug?

Odd behavior: `sparse(M::Array{Dual{Float64}})` zeroes dual values

I have no idea why this is happening, and it took me a long time to figure out that this was my issue in my application... Anyway I thought I should report it: I think calling sparse on a full dual-valued matrix "kills" the dual parts of the elements with no real part. MWE:

julia> M = zeros(2,2) + ones(2,2) * ε
2×2 Array{Dual{Float64},2}:
 0.0+1.0ɛ  0.0+1.0ɛ
 0.0+1.0ɛ  0.0+1.0ɛ

julia> sparse(M)
2×2 SparseMatrixCSC{Dual{Float64},Int64} with 0 stored entries

What's going on here? Am I missing something?

Replace this package with ForwardDiff's `Dual`?

It's annoying to have to port functionality between these two packages. It would be better if there were only be one Dual implementation that development efforts can focus on, if we can get away with it. AFAIK, ForwardDiff's Dual number is strictly superior to the version implemented in this package - in addition to allowing nested dual (i.e. hyper-dual) numbers and N partials derivatives at once, it features the same performance as this package's Dual in the single-element case, and is just as easy to write new derivative definitions for.

I propose we pull out ForwardDiff's Dual type and place it here instead. Of course, we'd want to ensure that we include any derivative definitions here that are missing from ForwardDiff, and vice versa.

The only controversial part of this change is that ForwardDiff's Dual is a subtype of Real, whereas the Dual here is a subtype of Number. I would require that we keep the proposed Dual a subtype of Real, since the goal of a dual number implementation is to programmatically support the Real interface. I really don't think this will matter in any practical sense, but others may disagree.

implement sinpi/cospi for Dual

There are probably other special functions that need to be implemented.