fix gauss tests

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

test/gauss.jl

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

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")