bump version

Project.toml

@@ -1,21 +1,24 @@
 name = "TaylorTest"
 uuid = "967b12e5-70be-421a-a124-e33167727e0a"
 authors = ["Gaspard Jankowiak <gaspard.jankowiak@uni-graz.at>"]
-version = "0.1.0"
+version = "0.2.1"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
 
 [compat]
 LinearAlgebra = "1.11.0"
 SpecialFunctions = "2.5.0"
 TensorOperations = "5.1.3"
 Test = "1.11.0"
+UnicodePlots = "3.7.2"
 
 [extras]
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
 
 [targets]
-test = ["SpecialFunctions"]
+test = ["SpecialFunctions", "TensorOperations"]

README.md

@@ -1,3 +1,35 @@
 # TaylorTest
 
 Simple package to check derivatives
+
+# Usage
+
+```
+  check(f, Jf, x[, constant_components]; taylortestplot=false, taylortestdirection=nothing, f_kwargs...)
+```
+
+Returns true if `Jf` approximates the derivative/gradient/Jacobian of `f` at point `x` (along a random direction unless specified using `taylortestdirection`).
+`f_kwargs` are keywords arguments to be passed to `f` and `Jf`.
+`constant_components` is an optional `Vector{Int}` corresponding to components of the direction which should be set to zero,
+effectively ignoring the dependency of `f` on these components.
+If `taylortestplot` is `true`, a log-log plot of the error against the perturbation size will be shown.
+
+
+```
+  check!(f!, Jf!, x, size_f_x, size_Jf_x, [, constant_components]; taylortestplot=false, taylortestdirection=nothing, f_kwargs...)
+```
+
+Like `check` but handling non-allocating functions. Output size for both `f!` and the Jacobian `Jf!` must be provided (as `Tuple`s) via `size_f_x` and `size_Jf_x`.
+
+## Examples (see `test` directory for more)
+```julia
+import TaylorTest
+
+f = x -> cos(x)
+Jf = x -> -sin(x)
+
+TaylorTest.check(f, Jf, rand())
+
+# [ Info: Approximation order ~ 1.0
+# true
+```

src/TaylorTest.jl

@@ -1,16 +1,42 @@
 module TaylorTest
 
+export check, check!
+
 import LinearAlgebra: norm, dot
 import TensorOperations: @tensor
+import UnicodePlots
 
-function check(f, Jf, x, constant_components::Vector{Int}=Int[]; f_kwargs...)
-  direction = 2 * (rand(size(x)...) .- 0.5)
-  for cst in constant_components
-    direction[cst] = 0
-  end
+"""
+```
+  check(f, Jf, x[, constant_components]; taylortestplot=false, taylortestdirection=nothing, f_kwargs...)
+```
 
-  if direction isa Array
-    direction ./= norm(direction)
+Returns true if `Jf` approximates the derivative/gradient/Jacobian of `f` at point `x` (along a random direction unless specified using `taylortestdirection`).
+`f_kwargs` are keywords arguments to be passed to `f` and `Jf`.
+`constant_components` is an optional `Vector{Int}` corresponding to components of the direction which should be set to zero,
+effectively ignoring the dependency of `f` on these components.
+If `taylortestplot` is `true`, a log-log plot of the error against the perturbation size will be shown.
+
+See also: `check!`
+
+# Examples
+```julia-repl
+julia> f = x -> cos(x); Jf = x -> -sin(x); check(f, Jf, rand())
+true
+```
+"""
+function check(f, Jf, x, constant_components::Vector{Int}=Int[]; taylortestplot::Bool=false, taylortestdirection=nothing, f_kwargs...)
+  if isnothing(taylortestdirection)
+    direction = 2 * (rand(size(x)...) .- 0.5)
+    for cst in constant_components
+      direction[cst] = 0
+    end
+
+    if direction isa Array
+      direction ./= norm(direction)
+    end
+  else
+    direction = taylortestdirection
   end
 
   ε_array = 10.0 .^ (-5:1:-1)
@@ -28,15 +54,19 @@ function check(f, Jf, x, constant_components::Vector{Int}=Int[]; f_kwargs...)
     error_array = [norm(vec(f(x + ε * direction; f_kwargs...) - f_x) - ε * Jf_x * direction) for ε in ε_array]
   else
     @tensor begin
-      Jf_x_direction[i,j] := Jf_x[i,j,k] * direction[k]
+      Jf_x_direction[i, j] := Jf_x[i, j, k] * direction[k]
     end
 
     error_array = [norm(f(x + ε * direction; f_kwargs...) - f_x - ε * Jf_x_direction) for ε in ε_array]
   end
 
+  if taylortestplot
+    println(UnicodePlots.lineplot(ε_array, error_array; xscale=:log10, yscale=:log10, xlabel="ε", title="|| f(x + ε·d) - f(x) - ε Jf(x)[d] ||"))
+  end
+
   m = maximum(error_array)
   if m < 1e-8
-    @warn "f looks linear!"
+    @warn "f looks constant or linear!"
     return true
   end
 
@@ -45,4 +75,32 @@ function check(f, Jf, x, constant_components::Vector{Int}=Int[]; f_kwargs...)
   return isapprox(order, 1; atol=0.5)
 end
 
+"""
+```
+  check!(f!, Jf!, x, size_f_x, size_Jf_x, [, constant_components]; taylortestplot=false, taylortestdirection=nothing, f_kwargs...)
+```
+
+Like `check` but handling non-allocating functions. Output size for both `f!` and the Jacobian `Jf!` must be provided (as `Tuple`s) via `size_f_x` and `size_Jf_x`.
+
+See also: `check`
+"""
+function check!(f!, Jf!, x, size_f_x::Tuple, size_Jf_x::Tuple, constant_components::Vector{Int}=Int[];
+  taylortestdirection=nothing, taylortestplot::Bool=false, f_kwargs...)
+
+  f = x -> begin
+    # must be allocated everytime to avoid aliasing
+    f_x = zeros(size_f_x...)
+    f!(f_x, x; f_kwargs...)
+    f_x
+  end
+  Jf = x -> begin
+    # must be allocated everytime to avoid aliasing
+    Jf_x = zeros(size_Jf_x...)
+    Jf!(Jf_x, x; f_kwargs...)
+    Jf_x
+  end
+
+  return check(f, Jf, x, constant_components, taylortestdirection=taylortestdirection, taylortestplot=taylortestplot)
+end
+
 end # module TaylorTest

test/erf.jl

@@ -62,12 +62,12 @@ function compute_H_f2(p::Vector{Float64}; x::Float64=0.0)
   return H_f2
 end
 
-@testset "Erf aux functions" begin
+@testset "Erf auxiliary functions" begin
   L = 10 * (rand() - 0.5)
   μ = 10 * (rand() - 0.5)
   σ = 10rand()
 
-  x = σ * (2 * rand() - 0.5) + μ
+  x = σ * (2 * rand() - 1) + μ
   p = [L, μ, σ]
 
   @test TaylorTest.check(compute_f1, compute_J_f1, p; x=x)
@@ -98,7 +98,8 @@ end
 
     J_f1 = compute_J_f1(p; x=x)
     J_f2 = compute_J_f2(p; x=x)
-    return J_f1 .* gauss([sqrt(2), 0, 1]; x=f1) - J_f2 .* gauss([sqrt(2), 0, 1]; x=f2)
+
+    return J_f1 .* gauss([1, 0, 1 / sqrt(2)]; x=f1) - J_f2 .* gauss([1, 0, 1 / sqrt(2)]; x=f2)
   end
 
   function H_diff_of_erf(p::Vector{Float64}; x::Float64=0.0)
@@ -106,7 +107,7 @@ end
     L, μ, σ = p
     sq2σ = sqrt(2) * σ
 
-    # d²/dξ² ( erf( f1(x) ) - erf( f2(x) ) )/2 =   d²/dξ² f1(x) gauss([sqrt(2), 0, 1]; x=f1(x)) - d²/dξ² f2(x) gauss([sqrt(2), 0, 1]; x=f2(x))
+    # d²/dξ² ( erf( f1(x) ) - erf( f2(x) ) )/2 =   d²/dξ² f1(x) gauss([1, 0, 1/sqrt(2)]; x=f1(x)) - d²/dξ² f2(x) gauss([1, 0, 1/sqrt(2)]; x=f2(x))
     # + (d/dξ f1(x))^2 gauss([sqrt(2), 0, 1]; x=f1(x)) - d²/dξ² f2(x) gauss([sqrt(2), 0, 1]; x=f2(x))
 
     f1 = compute_f1(p; x=x)
@@ -118,9 +119,10 @@ end
     H_f1 = compute_H_f1(p; x=x)
     H_f2 = compute_H_f2(p; x=x)
 
-    return (H_f1 .* gauss([sqrt(2), 0, 1]; x=f1) - H_f2 .* gauss([sqrt(2), 0, 1]; x=f2)
-            -
-            J_f1' * J_f1 .* f1 .* gauss([2 * sqrt(2), 0, 1]; x=f1) + J_f2' * J_f2 .* f2 .* gauss([2 * sqrt(2), 0, 1]; x=f2))
+    g_f1 = gauss([1, 0, 1 / sqrt(2)]; x=f1)
+    g_f2 = gauss([1, 0, 1 / sqrt(2)]; x=f2)
+
+    return H_f1 .* g_f1 - H_f2 .* g_f2 - 2J_f1' * J_f1 .* f1 .* g_f1 + 2J_f2' * J_f2 .* f2 .* g_f2
   end
 
   L = 10 * (rand() - 0.5)

test/gauss.jl

@@ -1,50 +1,74 @@
-  function J_gauss(p::Vector{Float64}; x::Float64=0.0)
-    λ, μ, σ = p
-    c = [1 / λ (2 * (x - μ) / σ^2) (-1 + 2 * ((x - μ) / σ)^2) / σ]
-    return c .* gauss(p; x=x)
-  end
+function gauss(p::Vector{Float64}; x::Float64=0.0)
+  λ, μ, σ = p
+  return λ / sqrt(2π * σ^2) * exp(-(x - μ)^2 / 2σ^2)
+end
 
-  function H_gauss(p::Vector{Float64}; x::Float64=0.0)
-    λ, μ, σ = p
-    c = zeros(3, 3)
+function J_gauss(p::Vector{Float64}; x::Float64=0.0)
+  λ, μ, σ = p
+  c = [(1 / λ) ((x - μ) / σ^2) ((x - μ)^2 - σ^2) / σ^3]
+  return c .* gauss(p; x=x)
+end
 
-    # c[1,1] = 0
-    c[1, 2] = c[2, 1] = (2x) / (λ * σ^2) - (2 * μ) / (λ * σ^2)
-    c[1, 3] = c[3, 1] = (2 * (x - μ)^2 - σ^2) / (λ * σ^3)
+function H_gauss(p::Vector{Float64}; x::Float64=0.0)
+  λ, μ, σ = p
+  c = zeros(3, 3)
 
-    c[2, 2] = (4 * (x - μ)^2 - 2 * σ^2) / σ^4
-    c[2, 3] = c[3, 2] = (2λ * (x - μ) * (2μ^2 - 3σ^2 + 2x^2 - 4μ * x)) / (λ * σ^5)
+  # c[1,1] = 0
+  c[1, 2] = c[2, 1] = (x - μ) / (λ * σ^2)
+  c[1, 3] = c[3, 1] = -((μ + σ - x) * (-μ + σ + x)) / (λ * σ^3)
 
-    c[3, 3] = (2 * (σ^4 - 5σ^2 * (x - μ)^2 + 2 * (x - μ)^4)) / σ^6
-    return c .* gauss(p; x=x)
-  end
+  c[2, 2] = ((x - μ)^2 - σ^2) / σ^4
+  c[2, 3] = c[3, 2] = -((μ - x) * ((μ - x)^2 - 3σ^2)) / σ^5
+
+  c[3, 3] = (-5 * σ^2 * (μ - x)^2 + (μ - x)^4 + 2 * σ^4) / σ^6
+  return c .* gauss(p; x=x)
+end
 
 @testset "Gauss function" begin
-  λ = 10 * (rand() - 0.5)
+  λ = 10 * rand()
   μ = 10 * (rand() - 0.5)
   σ = 10rand()
-
-  x = σ * (2 * rand() - 0.5) + μ
   p = [λ, μ, σ]
 
+  x = σ * (2 * rand() - 1) + μ
+
+  @assert abs(gauss(p; x=x)) > 1e-7 "gauss(p; x=x) is too small: $(gauss(p; x=x))"
+
   @test TaylorTest.check(gauss, J_gauss, p; x=x)
   @test TaylorTest.check(J_gauss, H_gauss, p; x=x)
 end
 
 @testset "Composed Gauss function" begin
-  p = 1 .+ 5 * rand(6)
+  p = 1 .+ 5rand(6)
 
-  x = p[3] + p[4] + 2p[5] + sqrt(p[6]) * (2rand() - 1)
-
-  p_aux = _p -> [_p[1]/_p[2], _p[3] + _p[4] + 2_p[5], _p[6]]
+  p_aux = _p -> [_p[1] / _p[2], _p[3] + _p[4] + 2_p[5], _p[6]]
 
   J_p_aux = _p -> [1/_p[2] -_p[1]/_p[2]^2 0 0 0 0;
-                           0            0 1 1 2 0;
-                           0            0 0 0 0 1]
+    0 0 1 1 2 0;
+    0 0 0 0 0 1]
 
-  g = (_p; x=x) -> gauss(p_aux(_p); x=x)
-  Jg = (_p; x=x) -> J_gauss(p_aux(_p); x=x) * J_p_aux(p)
-  Hg = (_p; x=x) -> J_gauss(p_aux(_p); x=x) * J_p_aux(p)
+  H_p_aux = _p -> begin
+    h = zeros(3, 6, 6)
+    h[1, 1, 2] = h[1, 2, 1] = -1 / _p[2]^2
+    h[1, 2, 2] = 2_p[1] / _p[2]^3
+    h
+  end
+
+  g = (_p; x = 0.0) -> gauss(p_aux(_p); x=x)
+  Jg = (_p; x = 0.0) -> J_gauss(p_aux(_p); x=x) * J_p_aux(_p)
+  Hg = (_p; x = 0.0) -> begin
+    p = p_aux(_p)
+    _J_gauss = vec(J_gauss(p; x=x))
+    _H_p = H_p_aux(_p)
+    _J_p = J_p_aux(_p)
+    r = _J_p' * H_gauss(p; x=x) * _J_p
+    r .+= TO.tensorcontract([2, 3], _J_gauss, [1], _H_p, [1, 2, 3])
+    r
+  end
+
+  _p = p_aux(p)
+  x = _p[2] + _p[2] * 2 * (rand() - 0.5)
 
   @test TaylorTest.check(g, Jg, p; x=x)
+  @test TaylorTest.check(Jg, Hg, p; x=x)
 end

test/nonallocating.jl

@@ -0,0 +1,23 @@
+@testset "Non-allocating functions" begin
+  f! = (f_x, x) -> (f_x[1] = x[1]^2 - 2x[2]^2; f_x)
+  Jf! = (Jf_x, x) -> (Jf_x[1] = 2x[1]; Jf_x[2] = -4x[2]; Jf_x)
+  Hf! = (Hf_x, x) -> begin
+    Hf_x[1, 1] = 2.0
+    Hf_x[1, 2] = Hf_x[2, 1] = 0.0
+    Hf_x[2, 2] = -4.0
+    Hf_x
+  end
+
+  size_f_x = (1,)
+  size_Jf_x = (1, 2)
+  size_Hf_x = (2, 2)
+
+  f_x = zeros(size_f_x...)
+  Jf_x = zeros(size_Jf_x...)
+  Hf_x = zeros(size_Hf_x...)
+
+  x = rand(2)
+
+  @test TaylorTest.check!(f!, Jf!, x, size_f_x, size_Jf_x)
+  @test TaylorTest.check!(Jf!, Hf!, x, size_Jf_x, size_Hf_x)
+end

test/runtests.jl

@@ -3,13 +3,10 @@ import Test: @test, @testset
 import TaylorTest
 
 import SpecialFunctions: erf
-
-function gauss(p::Vector{Float64}; x::Float64=0.0)
-  λ, μ, σ = p
-  return λ / sqrt(2π * σ^2) * exp(-((x - μ) / σ)^2)
-end
+import TensorOperations as TO
 
 include("trig_functions.jl")
 include("gauss.jl")
 include("erf.jl")
 include("tensors.jl")
+include("nonallocating.jl")

test/tensors.jl

@@ -26,3 +26,61 @@ end
   @test TaylorTest.check(compute_Jf, compute_Hf, x)
 
 end
+
+TP = TO.tensorproduct
+
+function prod(A, B)
+  return TP([1, 2, 3], A, [1, 2], B, [3])
+end
+
+function revprod(A, B)
+  return TP([1, 3, 2], B, [3], A, [1, 2])
+end
+
+function symprod(A, B)
+  return prod(A, B) + revprod(A, B)
+end
+
+@testset "Product of scalar and vector functions" begin
+  v = x -> [exp(2 * x[1]), cos(x[1] * x[2]), 1.0]
+  Jv = x -> [
+    2*exp(2 * x[1]) 0;
+    -x[2]*sin(x[1] * x[2]) -x[1]*sin(x[1] * x[2]);
+    0 0
+  ]
+  Hv = x -> begin
+    h = zeros(3, 2, 2)
+    h[1, :, :] = [
+      4*exp(2 * x[1]) 0;
+      0 0
+    ]
+
+    h[2, :, :] = [
+      -x[2]^2*cos(x[1] * x[2]) -sin(x[1] * x[2])-x[1]*x[2]*cos(x[1] * x[2]);
+      -sin(x[1] * x[2])-x[1]*x[2]*cos(x[1] * x[2]) -x[1]^2*cos(x[1] * x[2])
+    ]
+    h
+  end
+
+  φ = x -> x[1] * x[2]^2
+  ∇φ = x -> [x[2]^2; 2x[1] * x[2]]
+  Hφ = x -> [
+    0 2x[2];
+    2x[2] 2*x[1]
+  ]
+
+  f = x -> v(x) * φ(x)
+  Jf = x -> φ(x) * Jv(x) + v(x) * ∇φ(x)'
+  Hf = x -> φ(x) * Hv(x) + symprod(Jv(x), ∇φ(x)) + TP(v(x), [1], Hφ(x), [2, 3])
+
+  x = 2rand(2) .- 0.5
+
+  @test TaylorTest.check(v, Jv, x)
+  @test TaylorTest.check(Jv, Hv, x)
+
+  @test TaylorTest.check(φ, ∇φ, x)
+  @test TaylorTest.check(∇φ, Hφ, x)
+
+  @test TaylorTest.check(f, Jf, x)
+  @test TaylorTest.check(Jf, Hf, x)
+end