JuliaMath / Calculus.jl

Calculus functions in Julia
Other
278 stars 76 forks source link

Calculus.gradient return an Array of Float64 instead of Array of Type T #141

Open StevenSiew opened 4 years ago

StevenSiew commented 4 years ago

The calculus module is great for performing finite differences for derivative, jacobian and hessian

But when I tried to use it with my own type of floating point numbers, it fails and I track down why it fails.

julia> using PDFPs

julia> Base.float(x::PDFP) = x

julia> using Calculus

julia> Calculus.derivative( x->x^2 , PDFP(7) )  # returns 14
PDFP(0, 1, [1, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])

julia> Calculus.gradient(x -> sin(x[1]) + cos(x[2]), [PDFP(0), PDFP(0)])
2-element Array{Float64,1}:
 0.99999999999
 0.0 

Now Calculus.gradient return an Array of Float64 instead of Array of PDFP

I have track down the problem to URL https : //github.com/JuliaMath/Calculus.jl/blob/master/src/finite_difference.jl

the problem is in the line "g = Vector{Float64}(undef, length(x))"

function finite_difference(f,
                           x::AbstractVector{T},
                           dtype::Symbol = :central) where T <: Number
    # Allocate memory for gradient
    g = Vector{Float64}(undef, length(x))

    # Mutate allocated gradient
    finite_difference!(f, float(x), g, dtype)

    # Return mutated gradient
    return g
end

I think it should have been

g = Vector{eltype(x)}(undef, length(x))