Commit 7dc28722 authored by Alberto Ramos's avatar Alberto Ramos

First version of ADerrorsv2 in Julia

Working and crosschecked:

  - Computation of Gamma(t), drho(t) without replica
  - Indexing and error propagation
parents
# This file is machine-generated - editing it directly is not advised
[[AbstractFFTs]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "051c95d6836228d120f5f4b984dd5aba1624f716"
uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
version = "0.5.0"
[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
[[CommonSubexpressions]]
deps = ["Test"]
git-tree-sha1 = "efdaf19ab11c7889334ca247ff4c9f7c322817b0"
uuid = "bbf7d656-a473-5ed7-a52c-81e309532950"
version = "0.2.0"
[[CompilerSupportLibraries_jll]]
deps = ["Libdl", "Pkg"]
git-tree-sha1 = "7c4f882c41faa72118841185afc58a2eb00ef612"
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
version = "0.3.3+0"
[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
[[DiffResults]]
deps = ["StaticArrays"]
git-tree-sha1 = "da24935df8e0c6cf28de340b958f6aac88eaa0cc"
uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
version = "1.0.2"
[[DiffRules]]
deps = ["NaNMath", "Random", "SpecialFunctions"]
git-tree-sha1 = "eb0c34204c8410888844ada5359ac8b96292cfd1"
uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
version = "1.0.1"
[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
[[FFTW]]
deps = ["AbstractFFTs", "FFTW_jll", "IntelOpenMP_jll", "Libdl", "LinearAlgebra", "MKL_jll", "Reexport"]
git-tree-sha1 = "14536c95939aadcee44014728a459d2fe3ca9acf"
uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
version = "1.2.2"
[[FFTW_jll]]
deps = ["Libdl", "Pkg"]
git-tree-sha1 = "6c975cd606128d45d1df432fb812d6eb10fee00b"
uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a"
version = "3.3.9+5"
[[ForwardDiff]]
deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "NaNMath", "Random", "SpecialFunctions", "StaticArrays"]
git-tree-sha1 = "869540e4367122fbffaace383a5bdc34d6e5e5ac"
uuid = "f6369f11-7733-5829-9624-2563aa707210"
version = "0.10.10"
[[IntelOpenMP_jll]]
deps = ["Libdl", "Pkg"]
git-tree-sha1 = "fb8e1c7a5594ba56f9011310790e03b5384998d6"
uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0"
version = "2018.0.3+0"
[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
[[LibGit2]]
deps = ["Printf"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
[[MKL_jll]]
deps = ["IntelOpenMP_jll", "Libdl", "Pkg"]
git-tree-sha1 = "0ce9a7fa68c70cf83c49d05d2c04d91b47404b08"
uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7"
version = "2020.1.216+0"
[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
[[NaNMath]]
git-tree-sha1 = "928b8ca9b2791081dc71a51c55347c27c618760f"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "0.3.3"
[[OpenSpecFun_jll]]
deps = ["CompilerSupportLibraries_jll", "Libdl", "Pkg"]
git-tree-sha1 = "d51c416559217d974a1113522d5919235ae67a87"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.3+3"
[[Pkg]]
deps = ["Dates", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
[[Reexport]]
deps = ["Pkg"]
git-tree-sha1 = "7b1d07f411bc8ddb7977ec7f377b97b158514fe0"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "0.2.0"
[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
[[SpecialFunctions]]
deps = ["OpenSpecFun_jll"]
git-tree-sha1 = "d8d8b8a9f4119829410ecd706da4cc8594a1e020"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "0.10.3"
[[StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "5c06c0aeb81bef54aed4b3f446847905eb6cbda0"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "0.12.3"
[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
[[Test]]
deps = ["Distributed", "InteractiveUtils", "Logging", "Random"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
name = "ADerrors"
uuid = "5e92007d-7bf1-471c-8ceb-4591b8b567a9"
authors = ["Alberto Ramos <alberto.ramos@ific.uv.es>"]
version = "0.1.0"
[deps]
FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
###
### "THE BEER-WARE LICENSE":
### Alberto Ramos wrote this file. As long as you retain this
### notice you can do whatever you want with this stuff. If we meet some
### day, and you think this stuff is worth it, you can buy me a beer in
### return. <alberto.ramos@cern.ch>
###
### file: ADerrors.jl
### created: Wed Jun 17 13:00:26 2020
###
module ADerrors
import ForwardDiff, Statistics, FFTW#, BDIO
# Include data types
include("ADerrorsTypes.jl")
# Include computation of autoCF
include("ADerrorsCF.jl")
# Math operations
include("ADerrorsMath.jl")
# I/O
include("ADerrorsIO.jl")
export err, value, derror, taui, dtaui, window
export uwreal, uwerr
end # module
###
### "THE BEER-WARE LICENSE":
### Alberto Ramos wrote this file. As long as you retain this
### notice you can do whatever you want with this stuff. If we meet some
### day, and you think this stuff is worth it, you can buy me a beer in
### return. <alberto.ramos@cern.ch>
###
### file: ADerrorsCF.jl
### created: Wed Jun 17 13:19:13 2020
###
const MIN_LENGTH = 500
const MIN_REP_LENGTH = 4
const nprm=50
const iprm = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41,
43, 47, 53, 59, 61, 67, 71, 73,79,83,89,97,101,
103,107,109,113,127,131,137,139,149,151,157,163,
167,173,179,181,191,193,197,199,211,223,227,229]
const DEFAULT_STAU = 4.0
const DO_BIN = false
function get_nbin_vec(nd::Array{Int64, 1})::Int64
nbin::Int64 = 1
if (!DO_BIN)
return nbin
end
if (any(nd .== 1))
return nbin
end
nl = deepcopy(nd)
for i = 1:nprm
while (all(nl .% iprm[i] .== 0))
nl = map(x->div(x,iprm[i]), nl)
if (all(nl .< MIN_LENGTH))
return nbin
end
nbin = nbin * iprm[i]
end
end
return nbin
end
function bin_data(vin::Array{Float64, 1}, nbin::Int64)
vout = similar(Array{Float64, 1}, div(length(vin), nbin))
is::Int64 = 1
for i = 1:length(vout)
ie = is + nbin - 1
vout[i] = Statistics.mean(vin[is:ie])
is = ie+1
end
return vout
end
function uwcls(data::Vector{Float64}, id::Int64, ws::wspace)
if (length(data) == 2)
ws.nob += 1
new = fbd(1, 1, [data[2]], [1], [1], [1], Dict{Int64,Vector{Complex{Float64}}}())
push!(ws.fluc, new)
push!(ws.map_nob, id)
if (!haskey(ws.map_ids, id))
ws.map_ids[id] = ws.nob
end
p = [false for n in 1:ws.nob]
p[end] = true
d = [0.0 for n in 1:ws.nob]
d[end] = 1.0
return uwreal(data[1], 0.0, 0.0,
p, d, Vector{Int64}(), Vector{cfdata}())
else
ws.nob += 1
avg = Statistics.mean(data)
push!(ws.map_nob, id)
if (!haskey(ws.map_ids, id))
ws.map_ids[id] = ws.nob
end
fseries = Dict{Int64,Vector{Complex{Float64}}}()
nbin = get_nbin_vec([length(data)])
nbdt = div(length(data), nbin)
datapad = [bin_data(data .- avg, nbin);
zeros(Float64, nextpow(2,2*nbdt+1)-nbdt) ]
fseries[1] = FFTW.fft(datapad)
new = fbd(length(data), nbin,
data .- avg,
[length(data)],
[1], [length(data)],
fseries)
push!(ws.fluc, new)
p = [false for n in 1:ws.nob]
p[end] = true
d = [0.0 for n in 1:ws.nob]
d[end] = 1.0
return uwreal(avg, 0.0, 0.0,
p, d, Vector{Int64}(), Vector{cfdata}())
end
end
function unique_ids!(a::uwreal, ws::wspace)
if (length(a.ids) == 0)
for i in 1:length(a.prop)
if (a.prop[i])
if (any(a.ids[1:end] == ws.map_nob[i]))
continue
end
push!(a.ids, ws.map_nob[i])
end
end
end
return length(a.ids)
end
function wopt_ulli(nd::Int64, stau::Float64, gmm::Vector{Float64})
tiw = 0.5
if (gmm[1] == 0.0)
return 1
else
for i in 2:length(gmm)
tiw = tiw + gmm[i]/gmm[1]
if (tiw <= 0.5)
return i
else
tau = stau/log((2.0*tiw+1.0)/(2.0*tiw-1.0))
gw = exp(-(i-1.0)/tau) - tau/sqrt((i-1.0)*nd)
if (gw < 0.0)
return i
end
end
end
end
end
function uwerror(a::uwreal, ws::wspace, wpm::Dict{Int64,Vector{Float64}})
nid = unique_ids!(a, ws)
if (length(a.cfd) == 0)
for j in 1:nid
new = cfdata(0.0, 0.0, 0.0, 0, Vector{Float64}(), Vector{Float64}())
idx = ws.map_ids[a.ids[j]]
nd = ws.fluc[idx].nd
if (nd != 1)
nd_eff = div(nd, ws.fluc[idx].ibn)
(nt,ip) = findmax(ws.fluc[idx].ivrep)
nt = div(nt, 2*ws.fluc[idx].ibn)
nrep = length(ws.fluc[idx].ivrep)
ftemp = Dict{Int64,Vector{Complex{Float64}}}()
for i in 1:nrep
ns = length(ws.fluc[idx].fourier[i])
ftemp[i] = zeros(Complex{Float64}, ns)
end
end
for i in 1:length(a.prop)
if (a.prop[i] && (ws.map_nob[i] == a.ids[j]))
if (nd == 1)
new.var = new.var + a.der[i]*ws.fluc[i].delta[1]
else
for k in 1:nrep
ftemp[k] = ftemp[k] + a.der[i]*ws.fluc[i].fourier[k]
end
end
end
end
if (nd == 1)
new.var = new.var^2
else
new.gamm = zeros(nt)
for k in 1:nrep
ftemp[k] .= ftemp[k].*conj(ftemp[k])
FFTW.ifft!(ftemp[k])
for ig in 1:nt
new.gamm[ig] = real(ftemp[k][ig])
end
end
for ig in 1:nt
nrcnt = count(map(x -> div(x, 2*ws.fluc[idx].ibn), ws.fluc[idx].ivrep) .> ig-1)
new.gamm[ig] = new.gamm[ig] / (nd_eff - nrcnt*(ig-1))
end
iw = wopt_ulli(nd_eff, DEFAULT_STAU, new.gamm)
new.iw = iw
dbias = new.gamm[1] + 2.0*sum(new.gamm[2:iw])
new.gamm .= new.gamm .+ dbias/nd_eff
new.drho = zeros(nt)
if (new.gamm[1] != 0.0)
for i in 1:nt
is = max(1, i-iw-2) + 1
ie = i + iw-1
for k in is:ie
if (k < nt+1)
cont = -2.0*new.gamm[k]*new.gamm[i]/new.gamm[1]^2
else
cont = 0.0
end
if ((i+k-2) < nt)
cont = cont + new.gamm[i+k-1]/new.gamm[1]
end
if (abs(i-k) < nt)
cont = cont + new.gamm[abs(i-k)+1]/new.gamm[1]
end
new.drho[i] = new.drho[i] + cont^2
end
new.drho[i] = sqrt(new.drho[i]/nd_eff)
end
else
new.var = new.var^2
end
end
push!(a.cfd, new)
end
end
a.err = 0.0
a.derr = 0.0
for j in 1:nid
idx = ws.map_ids[a.ids[j]]
if (ws.fluc[idx].nd == 1)
a.cfd[j].taui = 0.5
a.cfd[k].dtuai = 0.0
vti = 0.0
else
ibn = ws.fluc[idx].ibn
nd_eff = div(ws.fluc[idx].nd, ibn)
(nt,ip) = findmax(ws.fluc[idx].ivrep)
nt = div(nt, 2*ibn)
if (a.cfd[j].gamm[1] == 0.0)
a.cfd[j].taui = 0.0
a.cfd[j].dtaui = 0.0
vti = 0.0
else
wp = zeros(4)
wp = get(wpm, a.ids[j], [-1.0,-1.0,-1.0,-1.0])
if (wp[1] > 0.0)
a.cfd[j].iw = round(wp[1])
elseif (wp[2] > 0.0)
a.cfd[j].iw = wopt_ulli(nd_eff, wp[2], a.cfd[j].gamm)
end
if (wp[3] > 0.0)
for k in 2:nt
if (a.cfd[j].drho[k]*wp[3] > a.cfd[j].gamm[k]/a.cfd[j].gamm[1])
continue
end
end
a.cfd[j].iw = k-1
end
if (wp[4] > 0.0)
texp = wp[4]/ibn
else
texp = 0.0
end
iw = a.cfd[j].iw
vti = 0.5 + sum(a.cfd[j].gamm[2:iw])/a.cfd[j].gamm[1]
a.cfd[j].dtaui = sqrt(vti^2 * (4.0*iw-2.0*vti+2.0)/nd_eff)
a.cfd[j].taui = vti + texp*a.cfd[j].gamm[iw+1]/a.cfd[j].gamm[1]
end
a.cfd[j].var = a.cfd[j].gamm[1] * 2.0*a.cfd[j].taui/nd_eff
end
a.err = a.err + a.cfd[j].var
if (iw > 1)
a.derr = a.derr + vti*(a.cfd[j].iw-0.5)/nd_eff
end
end
a.err = sqrt(a.err)
a.derr = sqrt(a.derr)
end
wsg = ADerrors.wspace(similar(Vector{ADerrors.fbd}, 0),
0,
similar(Vector{Int64}, 0),
Dict{Int64, Int64}())
empt = Dict{Int64,Vector{Float64}}()
uwreal(data::Vector{Float64}, id::Int64) = ADerrors.uwcls(data::Vector{Float64}, id::Int64, wsg)
uwerr(a::uwreal) = ADerrors.uwerror(a::uwreal, wsg, empt)
uwerr(a::uwreal, wpm::Dict{Int64,Vector{Float64}}) = ADerrors.uwerror(a::uwreal, wsg, wpm::Dict{Int64,Vector{Float64}})
neid(a::uwreal) = ADerrors.unique_ids!(a::uwreal, wsg)
err(a::uwreal) = a.err
value(a::uwreal) = a.mean
derror(a::uwreal) = a.derr
taui(a::uwreal, i::Int64) = a.cfd[i].taui
dtaui(a::uwreal, i::Int64) = a.cfd[i].taui
window(a::uwreal, i::Int64) = a.cfd[i].iw
###
### "THE BEER-WARE LICENSE":
### Alberto Ramos wrote this file. As long as you retain this
### notice you can do whatever you want with this stuff. If we meet some
### day, and you think this stuff is worth it, you can buy me a beer in
### return. <alberto.ramos@cern.ch>
###
### file: ADerrorsMath.jl
### created: Thu Jun 18 08:13:30 2020
###
for op in (:sin, :cos, :tan, :log, :exp, :sqrt, :sind, :cosd, :tand, :sinpi, :cospi, :sinh, :cosh, :tanh, :asin, :acos, :atan, :asind, :acosd, :atand, :sec, :csc, :cot, :secd, :cscd, :cotd, :asec, :acsc, :acot, :asecd, :acscd, :acotd, :sech, :csch, :coth, :asinh, :acosh, :atanh, :asech, :acsch, :acoth, :sinc, :cosc, :deg2rad, :rad2deg, :log2, :log10, :log1p, :exp2, :exp10, :expm1, :-)
@eval function Base.$op(a::uwreal)
return uwreal(Base.$op(a.mean),
a.prop, ForwardDiff.derivative($op, a.mean)*a.der)
end
end
for op in (:+, :-, :*, :/, :^, :atan, :hypot)
@eval function Base.$op(a::uwreal, b::uwreal)
function fvec(x::Vector)
return Base.$op(x[1], x[2])
end
grad = ForwardDiff.gradient(fvec, [a.mean, b.mean])
if (length(a.der) > length(b.der))
p = similar(a.prop)
d = similar(a.der)
for i in 1:length(b.der)
d[i] = grad[1]*a.der[i] + grad[2]*b.der[i]
p[i] = a.prop[i] || b.prop[i]
end
for i in length(b.der)+1:length(a.der)
d[i] = grad[1]*a.der[i]
p[i] = a.prop[i]
end
else
p = similar(b.prop)
d = similar(b.der)
for i in 1:length(a.der)
d[i] = grad[1]*a.der[i] + grad[2]*b.der[i]
p[i] = a.prop[i] || b.prop[i]
end
for i in length(a.der)+1:length(b.der)
d[i] = grad[2]*b.der[i]
p[i] = b.prop[i]
end
end
return uwreal(Base.$op(a.mean,b.mean), p, d)
end
@eval function Base.$op(a::uwreal, b::Number)
function fvec(x::Vector)
return Base.$op(x[1], b)
end
grad = ForwardDiff.gradient(fvec,[a.mean])
return uwreal(Base.$op(a.mean,b), a.prop, grad[1]*a.der)
end
@eval function Base.$op(b::Number, a::uwreal)
function fvec(x::Vector)
return Base.$op(b, x[1])
end
grad = ForwardDiff.gradient(fvec,[a.mean])
return uwreal(Base.$op(b,a.mean), a.prop, grad[1]*a.der)
end
end
for op in (:<, :>, :≤, :≥, :≠)
@eval Base.$op(a::uwreal, b::uwreal) = Base.$op(a.mean, b.mean)
@eval Base.$op(a::uwreal, b::Number) = Base.$op(a.mean, b)
@eval Base.$op(a::Number, b::uwreal) = Base.$op(a, b.mean)
end
###
### "THE BEER-WARE LICENSE":
### Alberto Ramos wrote this file. As long as you retain this
### notice you can do whatever you want with this stuff. If we meet some
### day, and you think this stuff is worth it, you can buy me a beer in
### return. <alberto.ramos@cern.ch>
###
### file: ADerrorsTypes.jl
### created: Wed Jun 17 13:00:32 2020
###
mutable struct cfdata
var::Float64
taui::Float64
dtaui::Float64
iw::Int64
gamm::Vector{Float64}
drho::Vector{Float64}
end
mutable struct uwreal
mean::Float64
err::Float64
derr::Float64
prop::Array{Bool, 1}
der::Array{Float64, 1}
ids::Array{Int64, 1}
cfd::Vector{cfdata}
end
mutable struct fbd
nd::Int64
ibn::Int64
delta::Array{Float64, 1}
ivrep::Array{Int64, 1}
is::Array{Int64, 1}
ie::Array{Int64, 1}
fourier::Dict{Int64,Vector{Complex{Float64}}}
end
mutable struct wspace
fluc::Array{fbd, 1}
nob::Int64
map_nob::Array{Int64, 1} # The id that corresponds to each ob
map_ids::Dict{Int64, Int64} # For each id, a fluc index
end
uwreal(v::Float64, n::Int) = uwreal(v, 0.0, 0.0,
[false for i in 1:n], [0.0 for i in 1:n],
Vector{Int64}(), Vector{cfdata}())
uwreal(v::Float64, prop::Vector{Bool}, der::Vector{Float64}) = uwreal(v, 0.0, 0.0,
prop, der,
Vector{Int64}(), Vector{cfdata}())
#/usr/bin/env julia14
# Start test script
using ADerrors
using Test
println("Test [test1.jl]")
@time @test include("test1.jl")
a = uwreal([1.0, 0.1], 1)
b = uwreal([0.5, 0.023], 23)
c = 1.0 + sin(a+b)
d = sin(a)*cos(b) + cos(a)*sin(b) - 3.0
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
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment