From 0b1a5ad3af8be2a4944f2ddd30e62bb2067be62e Mon Sep 17 00:00:00 2001 From: Gaspard Jankowiak Date: Tue, 11 Feb 2025 13:42:14 +0100 Subject: [PATCH] fix gauss tests --- test/erf.jl | 18 +++++++++++------- test/gauss.jl | 49 +++++++++++++++++++++++++++--------------------- test/runtests.jl | 6 +----- 3 files changed, 40 insertions(+), 33 deletions(-) diff --git a/test/erf.jl b/test/erf.jl index d12f98e..98e50f3 100644 --- a/test/erf.jl +++ b/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) @@ -130,6 +132,8 @@ 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 diff --git a/test/gauss.jl b/test/gauss.jl index b245694..65c9620 100644 --- a/test/gauss.jl +++ b/test/gauss.jl @@ -1,32 +1,39 @@ - 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 @@ -39,8 +46,8 @@ end 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) diff --git a/test/runtests.jl b/test/runtests.jl index 91cf2a7..09a9ec0 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,11 +3,7 @@ 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")