Basic finite difference functionality - performant derivatives of real-valued callables

Author: dextorious

src/DiffEqDiffTools.jl

@@ -2,6 +2,6 @@ __precompile__()
 
 module DiffEqDiffTools
 
-# package code goes here
+include("finitediff.jl")
 
 end # module

src/finitediff.jl

@@ -0,0 +1,128 @@
+#=
+Very heavily inspired by Calculus.jl, but with an emphasis on performance and DiffEq API convenience.
+=#
+
+#=
+Compute the finite difference interval epsilon.
+Reference: Numerical Recipes, chapter 5.7.
+=#
+@inline function compute_epsilon{T<:Real}(::Type{Val{:forward}}, x::T, eps_sqrt::T=sqrt(eps(T)))
+    eps_sqrt * max(one(T), abs(x))
+end
+
+@inline function compute_epsilon{T<:Real}(::Type{Val{:central}}, x::T, eps_cbrt::T=cbrt(eps(T)))
+    eps_cbrt * max(one(T), abs(x))
+end
+
+@inline function compute_epsilon{T<:Complex}(::Type{Val{:complex}}, x::T)
+    eps(real(x))
+end
+
+
+#=
+Compute the derivative df of a real-valued callable f on a collection of points x.
+Generic fallbacks for AbstractArrays that are not StridedArrays.
+=#
+function finite_difference{T<:Real}(f, x::AbstractArray{T}, ::Type{Val{:central}}, ::Union{Void,AbstractArray{T}}=nothing)
+    df = zeros(T, size(x))
+    finite_difference!(df, f, x, Val{:central})
+end
+
+function finite_difference!{T<:Real}(df::AbstractArray{T}, f, x::AbstractArray{T}, ::Type{Val{:central}}, ::Union{Void,AbstractArray{T}}=nothing)
+    eps_sqrt = sqrt(eps(T))
+    epsilon = compute_epsilon.(Val{:central}, x, eps_sqrt)
+    @. df = (f(x+epsilon) - f(x-epsilon)) / (2 * epsilon)
+end
+
+function finite_difference{T<:Real}(f, x::AbstractArray{T}, ::Type{Val{:forward}}, f_x::AbstractArray{T}=f.(x))
+    df = zeros(T, size(x))
+    finite_difference!(df, f, x, Val{:forward}, f_x)
+end
+
+function finite_difference!{T<:Real}(df::AbstractArray{T}, f, x::AbstractArray{T}, ::Type{Val{:forward}}, f_x::AbstractArray{T}=f.(x))
+    eps_cbrt = cbrt(eps(T))
+    epsilon = compute_epsilon.(Val{:forward}, x, eps_cbrt)
+    @. df = (f(x+epsilon) - f_x) / epsilon
+end
+
+
+#=
+Compute the derivative df of a real-valued callable f on a collection of points x.
+Optimized implementations for StridedArrays.
+=#
+function finite_difference{T<:Real}(f, x::StridedArray{T}, ::Type{Val{:central}}, ::Union{Void,StridedArray{T}}=nothing)
+    df = zeros(T, size(x))
+    finite_difference!(df, f, x, Val{:central})
+end
+
+function finite_difference!{T<:Real}(df::StridedArray{T}, f, x::StridedArray{T}, ::Type{Val{:central}}, ::Union{Void,StridedArray{T}}=nothing)
+    eps_sqrt = sqrt(eps(T))
+    @inbounds for i in 1 : length(x)
+        epsilon = compute_epsilon(Val{:central}, x[i], eps_sqrt)
+        epsilon_double_inv = one(T) / (2*epsilon)
+        x_plus, x_minus = x[i]+epsilon, x[i]-epsilon
+        df[i] = (f(x_plus) - f(x_minus)) * epsilon_double_inv
+    end
+    df
+end
+
+function finite_difference{T<:Real}(f, x::StridedArray{T}, ::Type{Val{:forward}}, fx::Union{Void,StridedArray{T}})
+    df = zeros(T, size(x))
+    if typeof(fx) == Void
+        finite_difference!(df, f, x, Val{:forward})
+    else
+        finite_difference!(df, f, x, Val{:forward}, fx)
+    end
+    df
+end
+
+function finite_difference!{T<:Real}(df::StridedArray{T}, f, x::StridedArray{T}, ::Type{Val{:forward}})
+    eps_cbrt = cbrt(eps(T))
+    @fastmath @inbounds for i in 1 : length(x)
+        epsilon = compute_epsilon(Val{:forward}, x[i], eps_cbrt)
+        epsilon_inv = one(T) / epsilon
+        x_plus = x[i] + epsilon
+        df[i] = (f(x_plus) - f(x[i])) * epsilon_inv
+    end
+    df
+end
+
+function finite_difference!{T<:Real}(df::StridedArray{T}, f, x::StridedArray{T}, ::Type{Val{:forward}}, fx::StridedArray{T})
+    eps_cbrt = cbrt(eps(T))
+    @fastmath @inbounds for i in 1 : length(x)
+        epsilon = compute_epsilon(Val{:forward}, x[i], eps_cbrt)
+        epsilon_inv = one(T) / epsilon
+        x_plus = x[i] + epsilon
+        df[i] = (f(x_plus) - fx[i]) * epsilon_inv
+    end
+    df
+end
+
+#=
+Compute the derivative df of a real-valued callable f on a collection of points x.
+Single point implementations.
+=#
+function finite_difference{T<:Real}(f, x::T, t::DataType, f_x::Union{Void,T}=nothing)
+    epsilon = compute_epsilon(t, x)
+    finite_difference_kernel(f, x, t, epsilon, f_x)
+end
+
+@inline function finite_difference_kernel{T<:Real}(f, x::T, ::Type{Val{:forward}}, epsilon::T, f_x::Union{Void,T})
+    if typeof(f_x) == Void
+        return (f(x+epsilon) - f(x)) / epsilon
+    else
+        return (f(x+epsilon) - f_x) / epsilon
+    end
+end
+
+@inline function finite_difference_kernel{T<:Real}(f, x::T, ::Type{Val{:central}}, epsilon::T, ::Union{Void,T}=nothing)
+    (f(x+epsilon) - f(x-epsilon)) / (2 * epsilon)
+end
+
+# TODO: derivatives for complex-valued callables
+
+
+#=
+Compute the Jacobian matrix of a real-valued callable f.
+=#
+# TODO

test/finitedifftests.jl

@@ -0,0 +1,8 @@
+x = collect(linspace(-2π, 2π, 100))
+y = sin.(x)
+df = zeros(100)
+df_ref = cos.(x)
+
+@test maximum(abs.(DiffEqDiffTools.finite_difference!(df, sin, x, Val{:central}) - df_ref)) < 1e-8
+@test maximum(abs.(DiffEqDiffTools.finite_difference!(df, sin, x, Val{:forward}) - df_ref)) < 1e-4
+@test maximum(abs.(DiffEqDiffTools.finite_difference!(df, sin, x, Val{:forward}, y) - df_ref)) < 1e-4

test/runtests.jl

@@ -1,5 +1,4 @@
 using DiffEqDiffTools
 using Base.Test
 
-# write your own tests here
-@test 1 == 2
+include("finitedifftests.jl")