From 76136dbc0599a024afe6c425016c587cb6f544b1 Mon Sep 17 00:00:00 2001 From: Gaspard Jankowiak Date: Mon, 10 Feb 2025 12:08:26 +0100 Subject: [PATCH] import --- Project.toml | 21 +++++++ src/TaylorTest.jl | 48 +++++++++++++++ test/erf.jl | 137 +++++++++++++++++++++++++++++++++++++++++ test/gauss.jl | 50 +++++++++++++++ test/runtests.jl | 15 +++++ test/tensors.jl | 28 +++++++++ test/trig_functions.jl | 4 ++ 7 files changed, 303 insertions(+) create mode 100644 Project.toml create mode 100644 src/TaylorTest.jl create mode 100644 test/erf.jl create mode 100644 test/gauss.jl create mode 100644 test/runtests.jl create mode 100644 test/tensors.jl create mode 100644 test/trig_functions.jl diff --git a/Project.toml b/Project.toml new file mode 100644 index 0000000..5b354b6 --- /dev/null +++ b/Project.toml @@ -0,0 +1,21 @@ +name = "TaylorTest" +uuid = "967b12e5-70be-421a-a124-e33167727e0a" +authors = ["Gaspard Jankowiak "] +version = "0.1.0" + +[deps] +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[compat] +LinearAlgebra = "1.11.0" +SpecialFunctions = "2.5.0" +TensorOperations = "5.1.3" +Test = "1.11.0" + +[extras] +SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" + +[targets] +test = ["SpecialFunctions"] diff --git a/src/TaylorTest.jl b/src/TaylorTest.jl new file mode 100644 index 0000000..807b9b9 --- /dev/null +++ b/src/TaylorTest.jl @@ -0,0 +1,48 @@ +module TaylorTest + +import LinearAlgebra: norm, dot +import TensorOperations: @tensor + +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 + + if direction isa Array + direction ./= norm(direction) + end + + ε_array = 10.0 .^ (-5:1:-1) + + n = size(ε_array) + f_x = f(x; f_kwargs...) + Jf_x = Jf(x; f_kwargs...) + + if ndims(f_x) == 0 + # if f is a scalar function, avoid potential ∇f vs Jac(f) by using dot + error_array = [norm(f(x + ε * direction; f_kwargs...) - (f_x + ε * dot(Jf_x, direction))) for ε in ε_array] + elseif ndims(f_x) == 1 || prod(size(f_x)) == maximum(size(f_x)) + # f is essentially a vector, potentially horizontal + # f simply use matrix multiplication + 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] + end + + error_array = [norm(f(x + ε * direction; f_kwargs...) - f_x - ε * Jf_x_direction) for ε in ε_array] + end + + m = maximum(error_array) + if m < 1e-8 + @warn "f looks linear!" + return true + end + + order = trunc(([ones(n) log.(ε_array)]\log.(error_array))[2] - 1; digits=2) + @info "Approximation order ~ $order" + return isapprox(order, 1; atol=0.5) +end + +end # module TaylorTest diff --git a/test/erf.jl b/test/erf.jl new file mode 100644 index 0000000..d12f98e --- /dev/null +++ b/test/erf.jl @@ -0,0 +1,137 @@ +function compute_f1(p::Vector{Float64}; x::Float64=0.0) + L, μ, σ = p + + sq2σ = sqrt(2) * σ + f1 = (x - (μ - 0.5L)) / sq2σ + + return f1 +end + +function compute_J_f1(p::Vector{Float64}; x::Float64=0.0) + L, μ, σ = p + + sq2σ = sqrt(2) * σ + J_f1 = [0.5 -1 -(x - (μ - 0.5L)) / σ] / sq2σ + return J_f1 +end + +function compute_H_f1(p::Vector{Float64}; x::Float64=0.0) + L, μ, σ = p + + sq2σ = sqrt(2) * σ + + J_f1 = [0.5 -1 -(x - (μ - 0.5L)) / σ] / sq2σ + + H_f1 = zeros(3, 3) + H_f1[1, 3] = H_f1[3, 1] = -0.5 / (σ * sq2σ) + H_f1[2, 3] = H_f1[3, 2] = 1 / (σ * sq2σ) + H_f1[3, 3] = -2 * J_f1[3] / σ + + return H_f1 +end + +function compute_f2(p::Vector{Float64}; x::Float64=0.0) + L, μ, σ = p + + sq2σ = sqrt(2) * σ + f2 = (x - (μ + 0.5L)) / sq2σ + + return f2 +end + +function compute_J_f2(p::Vector{Float64}; x::Float64=0.0) + L, μ, σ = p + + sq2σ = sqrt(2) * σ + J_f2 = [-0.5 -1 -(x - (μ + 0.5L)) / σ] / sq2σ + return J_f2 +end + +function compute_H_f2(p::Vector{Float64}; x::Float64=0.0) + L, μ, σ = p + + sq2σ = sqrt(2) * σ + + J_f2 = [-0.5 -1 -(x - (μ + 0.5L)) / σ] / sq2σ + + H_f2 = zeros(3, 3) + H_f2[1, 3] = H_f2[3, 1] = 0.5 / (σ * sq2σ) + H_f2[2, 3] = H_f2[3, 2] = 1 / (σ * sq2σ) + H_f2[3, 3] = -2 * J_f2[3] / σ + + return H_f2 +end + +@testset "Erf aux functions" begin + L = 10 * (rand() - 0.5) + μ = 10 * (rand() - 0.5) + σ = 10rand() + + x = σ * (2 * rand() - 0.5) + μ + p = [L, μ, σ] + + @test TaylorTest.check(compute_f1, compute_J_f1, p; x=x) + @test TaylorTest.check(compute_J_f1, compute_H_f1, p; x=x) + @test TaylorTest.check(compute_f2, compute_J_f2, p; x=x) + @test TaylorTest.check(compute_J_f2, compute_H_f2, p; x=x) +end + +@testset "Difference of error functions" begin + function diff_of_erf(p::Vector{Float64}; x::Float64=0.0) + L, μ, σ = p + + f1 = compute_f1(p; x=x) + f2 = compute_f2(p; x=x) + + return 0.5 * erf(f2, f1) + end + + function J_diff_of_erf(p::Vector{Float64}; x::Float64=0.0) + # d/dx erf(x)/2 = gauss([sqrt(2), 0, 1]; x=x) + 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)) + + f1 = compute_f1(p; x=x) + f2 = compute_f2(p; x=x) + + 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) + end + + function H_diff_of_erf(p::Vector{Float64}; x::Float64=0.0) + # d/dx erf(x)/2 = gauss([sqrt(2), 0, 1]; x=x) + 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ξ 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) + f2 = compute_f2(p; x=x) + + J_f1 = compute_J_f1(p; x=x) + J_f2 = compute_J_f2(p; x=x) + + 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)) + end + + L = 10 * (rand() - 0.5) + μ = 10 * (rand() - 0.5) + σ = 10rand() + + x = σ * (2 * rand() - 0.5) + μ + p = [L, μ, σ] + + @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 + + diff --git a/test/gauss.jl b/test/gauss.jl new file mode 100644 index 0000000..b245694 --- /dev/null +++ b/test/gauss.jl @@ -0,0 +1,50 @@ + 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 H_gauss(p::Vector{Float64}; x::Float64=0.0) + λ, μ, σ = p + c = zeros(3, 3) + + # 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) + + 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[3, 3] = (2 * (σ^4 - 5σ^2 * (x - μ)^2 + 2 * (x - μ)^4)) / σ^6 + return c .* gauss(p; x=x) + end + +@testset "Gauss function" begin + λ = 10 * (rand() - 0.5) + μ = 10 * (rand() - 0.5) + σ = 10rand() + + x = σ * (2 * rand() - 0.5) + μ + p = [λ, μ, σ] + + @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) + + 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]] + + 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] + + 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) + + @test TaylorTest.check(g, Jg, p; x=x) +end diff --git a/test/runtests.jl b/test/runtests.jl new file mode 100644 index 0000000..91cf2a7 --- /dev/null +++ b/test/runtests.jl @@ -0,0 +1,15 @@ +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 + +include("trig_functions.jl") +include("gauss.jl") +include("erf.jl") +include("tensors.jl") diff --git a/test/tensors.jl b/test/tensors.jl new file mode 100644 index 0000000..ae4f8b7 --- /dev/null +++ b/test/tensors.jl @@ -0,0 +1,28 @@ +compute_f = x -> [x[1]^2; x[2]^2 + x[3]; -x[3]^2 + x[4] + x[5]^3] + +compute_Jf = x -> [ + 2x[1] 0 0 0 0; + 0 2x[2] 1 0 0; + 0 0 -2x[3] 1 3*x[5]^2 +] + +function compute_Hf(x) + Hf = zeros(3, 5, 5) + + Hf[1, 1, 1] = 2 + Hf[2, 2, 2] = 2 + Hf[3, 3, 3] = -2 + Hf[3, 5, 5] = 6 * x[5] + + return Hf +end + +@testset "Tensors" begin + m, n = 3, 5 + + x = 2rand(n) .- 0.5 + + @test TaylorTest.check(compute_f, compute_Jf, x) + @test TaylorTest.check(compute_Jf, compute_Hf, x) + +end diff --git a/test/trig_functions.jl b/test/trig_functions.jl new file mode 100644 index 0000000..d754e15 --- /dev/null +++ b/test/trig_functions.jl @@ -0,0 +1,4 @@ +@testset "trigonometric functions" begin + @test TaylorTest.check(x -> cos(x), x -> -sin(x), 2π * rand()) + @test TaylorTest.check(x -> sin(x), x -> cos(x), 2π * rand()) +end