bump version

Project.toml

@@ -1,18 +1,20 @@
 name = "TaylorTest"
 uuid = "967b12e5-70be-421a-a124-e33167727e0a"
 authors = ["Gaspard Jankowiak <gaspard.jankowiak@uni-graz.at>"]
-version = "0.1.2"
+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"

README.md

@@ -5,15 +5,23 @@ Simple package to check derivatives
 # Usage
 
 ```
-  check(f, Jf, x[, constant_components]; f_kwargs...)
+  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).
+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.
 
-## Examples
+
+```
+  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
 

src/TaylorTest.jl

@@ -1,17 +1,23 @@
 module TaylorTest
 
-export check
+export check, check!
 
 import LinearAlgebra: norm, dot
 import TensorOperations: @tensor
+import UnicodePlots
 
 """
-  `check(f, Jf, x[, constant_components]; f_kwargs...)`
+```
+  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).
+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
@@ -19,14 +25,18 @@ julia> f = x -> cos(x); Jf = x -> -sin(x); check(f, Jf, rand())
 true
 ```
 """
-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
+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)
+    if direction isa Array
+      direction ./= norm(direction)
+    end
+  else
+    direction = taylortestdirection
   end
 
   ε_array = 10.0 .^ (-5:1:-1)
@@ -50,6 +60,10 @@ function check(f, Jf, x, constant_components::Vector{Int}=Int[]; f_kwargs...)
     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 constant or linear!"
@@ -61,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

@@ -132,8 +132,6 @@ end
   x = σ * (2 * rand() - 0.5) + μ
   p = [L, μ, σ]
 
-  @show diff_of_erf(p; x=x)
-
   @test TaylorTest.check(diff_of_erf, J_diff_of_erf, p; x=x)
   @test TaylorTest.check(J_diff_of_erf, H_diff_of_erf, p, [2, 3]; x=x)
 end

test/gauss.jl

@@ -1,11 +1,11 @@
 function gauss(p::Vector{Float64}; x::Float64=0.0)
   λ, μ, σ = p
-  return λ / sqrt(2π * σ^2) * exp(-(x - μ)^2/2σ^2)
+  return λ / sqrt(2π * σ^2) * exp(-(x - μ)^2 / 2σ^2)
 end
 
 function J_gauss(p::Vector{Float64}; x::Float64=0.0)
   λ, μ, σ = p
-  c = [(1/λ) ((x - μ) / σ^2) ((x - μ)^2 - σ^2)/σ^3]
+  c = [(1 / λ) ((x - μ) / σ^2) ((x - μ)^2 - σ^2) / σ^3]
   return c .* gauss(p; x=x)
 end
 
@@ -14,13 +14,13 @@ function H_gauss(p::Vector{Float64}; x::Float64=0.0)
   c = zeros(3, 3)
 
   # c[1,1] = 0
-  c[1, 2] = c[2, 1] = (x-μ) / (λ * σ^2)
-  c[1, 3] = c[3, 1] = -((μ + σ - x)*(-μ + σ + x))/(λ*σ^3)
+  c[1, 2] = c[2, 1] = (x - μ) / (λ * σ^2)
+  c[1, 3] = c[3, 1] = -((μ + σ - x) * (-μ + σ + x)) / (λ * σ^3)
 
   c[2, 2] = ((x - μ)^2 - σ^2) / σ^4
-  c[2, 3] = c[3, 2] = -((μ - x)*((μ - x)^2 - 3σ^2))/σ^5
+  c[2, 3] = c[3, 2] = -((μ - x) * ((μ - x)^2 - 3σ^2)) / σ^5
 
-  c[3, 3] = (-5*σ^2*(μ - x)^2 + (μ - x)^4 + 2*σ^4)/σ^6
+  c[3, 3] = (-5 * σ^2 * (μ - x)^2 + (μ - x)^4 + 2 * σ^4) / σ^6
   return c .* gauss(p; x=x)
 end
 
@@ -39,19 +39,36 @@ end
 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

@@ -9,3 +9,4 @@ include("trig_functions.jl")
 include("gauss.jl")
 include("erf.jl")
 include("tensors.jl")
+include("nonallocating.jl")

test/tensors.jl

@@ -29,19 +29,6 @@ end
 
 TP = TO.tensorproduct
 
-function check_symmetry(a)
-  R = true
-  for i in axes(a, 1)
-    r = a[i,:,:]' == a[i,:,:]
-    R &= r
-    if !r
-      n = maximum(abs, a[i,:,:]' - a[i,:,:])
-      @warn "Tensor not symmetric at i=$i: $n"
-    end
-  end
-  return R
-end
-
 function prod(A, B)
   return TP([1, 2, 3], A, [1, 2], B, [3])
 end
@@ -76,14 +63,14 @@ end
   end
 
   φ = x -> x[1] * x[2]^2
-  ∇φ = x -> [x[2]^2; 2x[1]*x[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)'
+  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