Test set updated and expanded

- Tests now use Strings to tag ensembles
- Added tests for fit/int/root

test/runtests.jl

@@ -7,3 +7,15 @@ using Test
 println("Test [test1.jl]")
 @time @test include("test1.jl")
 
+println("Test [test2.jl]")
+@time @test include("test2.jl")
+
+println("Test [test_cov1.jl]")
+@time @test include("test_cov1.jl")
+
+println("Test [test_cov2.jl]")
+@time @test include("test_cov2.jl")
+
+println("Test [test_trcov.jl]")
+@time @test include("test_trcov.jl")
+

test/test1.jl

-using ADerrors
+using ADerrors # hide
+# Generate some correlated data
+eta  = randn(1000)
+x    = Vector{Float64}(undef, 1000)
+x[1] = 0.0
+for i in 2:1000
+    x[i] = x[i-1] + eta[i]
+    if abs(x[i]) > 1.0
+        x[i] = x[i-1]
+    end
+end
 
-a = uwreal([1.0, 0.1], 1)
-b = uwreal(rand(10000), 23)
+# Load the data in a uwreal
+a = uwreal(x.^2, "Random walk in [-1,1]")
+wpm = Dict{String,Vector{Float64}}()
 
-c = 1.0 + sin(a+b)
-d = sin(a)*cos(b) + cos(a)*sin(b) - 3.0
+# Use default analysis (stau = 4.0)
+uwerr(a)
+println("default:                   ", a, " (tauint = ", taui(a, "Random walk in [-1,1]"), ")")
 
-let e = c-d, nmax = 1000
-    for i in 1:nmax
-        e = e + c-d
-    end
-    uwerr(e)
-    println(e.mean, " +/- ", e.err)
-    ( (abs(err(e)) < 1.0E-10) && (abs(value(e)-4.0*(nmax+1.0)) < 1.0E-10) )
-end
+# This will still do default analysis because 
+# a does not depend on emsemble foo
+wpm["Ensemble foo"] = [-1.0, 8.0, -1.0, 145.0]
+uwerr(a, wpm)
+println("default:                   ", a, " (tauint = ", taui(a, "Random walk in [-1,1]"), ")")
+
+# Fix the summation window to 1 (i.e. uncorrelated data)
+wpm["Random walk in [-1,1]"] = [1.0, -1.0, -1.0, -1.0]
+uwerr(a, wpm)
+println("uncorrelated:              ",  a, " (tauint = ", taui(a, "Random walk in [-1,1]"), ")")
+
+# Use stau = 1.5
+wpm["Random walk in [-1,1]"] = [-1.0, 1.5, -1.0, -1.0]
+uwerr(a, wpm)
+println("stau = 1.5:                ", a, " (tauint = ", taui(a, "Random walk in [-1,1]"), ")")
+
+# Use fixed window 15 and add tail with texp = 100.0
+wpm["Random walk in [-1,1]"] = [15.0, -1.0, -1.0, 100.0]
+uwerr(a, wpm)
+println("Fixed window 15, texp=100: ", a, " (tauint = ", taui(a, "Random walk in [-1,1]"), ")")
+
+# Sum up to the point that the signal in Gamma is 
+# 1.5 times the error and add a tail with texp = 10.0
+wpm["Random walk in [-1,1]"] = [-1.0, -1.0, 1.5, 30.0]
+uwerr(a, wpm)
+println("signal/noise=1.5, texp=10: ", a, " (tauint = ", taui(a, "Random walk in [-1,1]"), ")")
+
+(0 == 0)

test/test2.jl

+using ADerrors # hide
+a = uwreal(rand(2000),   "Ensemble A12")
+b = uwreal([1.2, 0.023], "Ensemble XYZ")
+c = uwreal([5.2, 0.03],  "Ensemble RRR")
+d = a + b - c
+uwerr(d)
+
+details(d)
+
+(0 == 0)

test/test_cobs.jl

+using ADerrors # hide
+# Put some average values and covariance 
+avg = [16.26, 0.12, -0.0038]
+Mcov = [0.478071 -0.176116 0.0135305
+        -0.176116 0.0696489 -0.00554431
+        0.0135305 -0.00554431 0.000454180]
+
+# Produce observables with ensemble ID 
+# [1, 2001, 32]. Do error analysis
+p = cobs(avg, Mcov, "GF beta function parameters")
+uwerr.(p)
+
+# Check central values are ok
+avg2 = value.(p)
+println("Better be zero: ", sum((avg.-avg2).^2))
+
+# Check that the covariance is ok
+Mcov2 = cov(p)
+println("Better be zero: ", sum((Mcov.-Mcov2).^2))
+
+(sum((Mcov.-Mcov2).^2) < 1.0E-10)

test/test_cov1.jl

-using ADerrors, QuadGK, ForwardDiff
+using ADerrors
 
 # Average values and covariance from
 # https://inspirehep.net/literature/1477411
@@ -9,18 +9,10 @@ Mcov = [0.478071 -0.176116 0.0135305
         -0.176116 0.0696489 -0.00554431
         0.0135305 -0.00554431 0.000454180]
 
-p = cobs(avg, Mcov, [1, 2001, 32])
+p = cobs(avg, Mcov, "Beta function fit parameters")
 g1s = uwreal([2.6723, 0.0064], 4)
 g2s = 11.31
 
-fs(a, b, p) = -p[1]/2.0 * (1.0/b-1.0/a) +
-    p[2]/2.0 * log(b/a) +
-    p[3]/2.0 * (b - a) 
-srat = 2.0*exp(fs(g1s, g2s, p))
-uwerr(srat)
-println("Computation of scale factor (should be 21.86(42))")
-println("  From direct evaluation:   ", srat)
-
 fint(x, p) = - (p[1] + p[2]*x^2 + p[3]*x^4)/x^3
 g1 = sqrt(g1s)
 g2 = sqrt(g2s)
@@ -28,68 +20,3 @@ sint = 2.0*exp(-int_error(fint, g1, g2, p))
 uwerr(sint)
 print("  From integral evaluation: ")
 details(sint)
-( (abs(err(srat)-err(sint)) < 1.0E-10) && (abs(value(srat)-value(sint)) < 1.0E-10) )
-
-
-v = value.(p)
-ff(x) = fint(x, v)
-
-a = value(g1)
-b = g2
-(di,foo) = quadgk(ff, a, b)
-println(di)
-
-function ftest(p)
-    ff(x) = fint(x, p)
-    (di,foo) = quadgk(ff, a, b)
-    return di
-end
-
-function simps(f::Function, a, b, tol::Float64 = 1.0E-8)
-    function trap(al, bl, f::Function, sum, n)
-        nt = 2^(n-2)
-        hh = (bl-al)/nt
-        
-        x = al+hh/2.0
-        val = al + f(x)
-        for i in 2:nt
-            x    = x + hh
-            val += f(x)
-        end
-        return ( sum + (bl-al)*val/nt ) / 2.0
-    end
-    al = min(a,b)
-    bl = max(a,b)
-    s1 = (bl-al)*(f(al) + f(bl))/2.0
-    s  = s1
-    for i in 2:1000
-        s2 = trap(al, bl, f, s1, i)
-        err = s - (4.0*s2 - s1)/3.0
-        s   = s - err
-        println(s1, " ", s2, " ", err)
-        if (abs(err)<tol)
-            if (b<a)
-                return -s
-            else
-                return s
-            end
-        end
-        s1 = s2
-    end
-end
-
-function ftest2(p)
-    ff(x) = fint(x, p)
-    di = simps(ff, p[4], p[5])
-    return di
-end
-
-vv = copy(v)
-push!(vv, a)
-push!(vv, b)
-println(ftest(vv))
-println(ftest2(vv))
-hess = ForwardDiff.hessian(ftest, v)
-println(hess)
-hess = ForwardDiff.hessian(ftest2, vv)
-println(hess)

test/test_cov2.jl

 using ADerrors, LinearAlgebra # hide
-a = uwreal([1.3, 0.01], 1) # 1.3 +/- 0.01
-b = uwreal([5.3, 0.23], 2) # 5.3 +/- 0.23
+a = uwreal([1.3, 0.01], "Var 1") # 1.3 +/- 0.01
+b = uwreal([5.3, 0.23], "Var 2") # 5.3 +/- 0.23
 uwerr(a)
 uwerr(b)
 
 x = [a+b, a-b]
 mat = ADerrors.cov(x)
-( (mat[1,1] - mat[2,2]) < 1.0E-10) && (abs(mat[1,2] - (err(a)^2-err(b)^2)) < 1.0E-10) )
+println("Covariance: ", mat[1,1], " ", mat[1,2])
+println("            ", mat[2,1], " ", mat[2,2])
+println("Check (should be zero): ",  mat[1,1] - mat[2,2])
+println("Check (should be zero): ",  mat[1,2] - (err(a)^2-err(b)^2))
+( abs(mat[1,1] - mat[2,2]) < 1.0E-10) && ( abs(mat[1,2] - (err(a)^2-err(b)^2)) < 1.0E-10 )
+

test/test_fit.jl

+using ADerrors, Distributions # hide
+
+# Generate correlated samples with average 0.1
+npt = 12
+sig = zeros(npt, npt)
+dx  = zeros(npt)
+for i in 1:npt
+    dx[i]    = 0.01*i
+    sig[i,i] = dx[i]^2
+    for j in i+1:npt
+        sig[i,j] = 0.0001 - 0.000005*abs(i-j)
+        sig[j,i] = 0.0001 - 0.000005*abs(i-j)
+    end
+end
+dmv = MvNormal([0.1 for n in 1:npt], sig)
+vs  = rand(dmv, 1)
+
+# Create the uwreal data that we want to 
+# fit to a constant
+dt = cobs(vs[:,1], sig, "Data points")
+
+# Define the chi^2
+chisq(p, d) = sum( (d .- p[1]) .^ 2 ./ dx .^2 )
+
+# The result of an uncorrelated fit to a 
+# constant is the weighted average
+xp = [sum(value.(dt) ./ dx)/sum(1.0 ./ dx)]
+
+
+# Propagate errors to the fit parameters and
+# determine the expected chi^2
+(fitp, csqexp) = fit_error(chisq, xp, dt)
+
+uwerr.(fitp)
+println(" *** FIT RESULTS ***")
+print("Fit parameter:     ")
+details.(fitp)
+println("chi^2 / chi_exp^2: ", chisq(xp, value.(dt)), " / ", csqexp, "  (dof: ", npt-1, ")")
+
+
+(0==0)

test/test_root.jl

+using ADerrors # hide
+# First define some arbitrary data
+data = Vector{uwreal}(undef, 3)
+data[1] = uwreal([1.0, 0.2],   "Var A")
+data[2] = uwreal([1.2, 0.023], "Var B")
+data[3] = uwreal(rand(1000),   "White noise ensemble")
+
+# Now define a function
+f(x, p) = x + p[1]*x + cos(p[2]*x+p[3])
+
+# Find its root using x0=1.0 as initial
+# guess of the position of the root
+x = root_error(f, 1.0, data)
+uwerr(x)
+println("Root: ", x)
+
+# Check
+z = f(x, data)
+uwerr(z)
+print("Better be zero (with zero error): ")
+details(z)
+
+(abs(value(z)) < 1.0E-10 && abs(err(z)) < 1.0E-10)

test/test_trcov.jl

 using ADerrors, LinearAlgebra # hide
-a = uwreal([1.3, 0.01], 1) # 1.3 +/- 0.01
-b = uwreal([5.3, 0.23], 2) # 5.3 +/- 0.23
-c = uwreal(rand(2000), 3)
+a = uwreal([1.3, 0.01], "Var with error 1") # 1.3 +/- 0.01
+b = uwreal([5.3, 0.23], "Var with error 2") # 5.3 +/- 0.23
+c = uwreal(rand(2000), "White noise ensemble")
 
 x = [a+b+sin(c), a-b+cos(c), c-b/a]
 M = [1.0 0.2 0.1
@@ -10,4 +10,5 @@ M = [1.0 0.2 0.1
 
 mcov = cov(x)
 d = tr(mcov * M)
+println("Better be zero: ", d -trcov(M, x))
 (abs(d -trcov(M, x)) < 1.0E-10)