Commit c4696a8e authored by Alberto Ramos's avatar Alberto Ramos

Some optimization

parent ad26ac05
......@@ -196,8 +196,9 @@ function uwerror(a::uwreal, ws::wspace, wpm::Dict{Int64,Vector{Float64}})
nid = unique_ids!(a, ws)
if (length(a.cfd) == 0)
a.cfd = Vector{cfdata}(undef, nid)
for j in 1:nid
new = cfdata(0.0, 0.0, 0.0, 0, Vector{Float64}(), Vector{Float64}())
a.cfd[j] = cfdata()
idx = ws.map_ids[a.ids[j]]
nd = ws.fluc[idx].nd
if (nd != 1)
......@@ -215,7 +216,7 @@ function uwerror(a::uwreal, ws::wspace, wpm::Dict{Int64,Vector{Float64}})
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]
a.cfd[j].var = a.cfd[j].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]
......@@ -225,56 +226,60 @@ function uwerror(a::uwreal, ws::wspace, wpm::Dict{Int64,Vector{Float64}})
end
if (nd == 1)
new.var = new.var^2
a.cfd[j].var = a.cfd[j].var^2
else
new.gamm = zeros(nt)
a.cfd[j].gamm = zeros(nt)
for k in 1:nrep
ftemp[k] .= ftemp[k].*conj(ftemp[k])
FFTW.ifft!(ftemp[k])
for ig in 1:min(nt,length(ftemp[k]))
new.gamm[ig] = new.gamm[ig] + real(ftemp[k][ig])
a.cfd[j].gamm[ig] = a.cfd[j].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))
nrcnt = 0
for k in 1:nrep
if (div(ws.fluc[idx].ivrep[k], 2*ws.fluc[idx].ibn) > ig-1 )
nrcnt = nrcnt + 1
end
end
a.cfd[j].gamm[ig] = a.cfd[j].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
iw = wopt_ulli(nd_eff, DEFAULT_STAU, a.cfd[j].gamm)
a.cfd[j].iw = iw
gwin = view(a.cfd[j].gamm, 2:iw)
dbias = a.cfd[j].gamm[1] + 2.0*sum(gwin)
a.cfd[j].gamm .= a.cfd[j].gamm .+ dbias/nd_eff
new.drho = zeros(nt)
if (new.gamm[1] != 0.0)
a.cfd[j].drho = zeros(nt)
if (a.cfd[j].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
cont = -2.0*a.cfd[j].gamm[k]*a.cfd[j].gamm[i]/a.cfd[j].gamm[1]^2
else
cont = 0.0
end
if ((i+k-2) < nt)
cont = cont + new.gamm[i+k-1]/new.gamm[1]
cont = cont + a.cfd[j].gamm[i+k-1]/a.cfd[j].gamm[1]
end
if (abs(i-k) < nt)
cont = cont + new.gamm[abs(i-k)+1]/new.gamm[1]
cont = cont + a.cfd[j].gamm[abs(i-k)+1]/a.cfd[j].gamm[1]
end
new.drho[i] = new.drho[i] + cont^2
a.cfd[j].drho[i] = a.cfd[j].drho[i] + cont^2
end
new.drho[i] = sqrt(new.drho[i]/nd_eff)
a.cfd[j].drho[i] = sqrt(a.cfd[j].drho[i]/nd_eff)
end
else
new.var = new.var^2
a.cfd[j].var = a.cfd[j].var^2
end
end
push!(a.cfd, new)
end
end
......@@ -323,7 +328,8 @@ function uwerror(a::uwreal, ws::wspace, wpm::Dict{Int64,Vector{Float64}})
end
iw = a.cfd[j].iw
vti = 0.5 + sum(a.cfd[j].gamm[2:iw])/a.cfd[j].gamm[1]
gwin = view(a.cfd[j].gamm, 2:iw)
vti = 0.5 + sum(gwin)/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
......
......@@ -201,7 +201,8 @@ function details(a::uwreal, ws::wspace, io::IO=stdout, names::Dict{Int64, String
n = 0
for i in 1:length(a.cfd)
if (length(a.cfd[i].gamm) > 0)
idx = ws.map_ids[a.ids[i]]
if (ws.fluc[idx].nd > 1)
n = n + 1
end
end
......@@ -217,10 +218,10 @@ function details(a::uwreal, ws::wspace, io::IO=stdout, names::Dict{Int64, String
end
ip = sortperm(v, rev=true)
for i in 1:length(a.cfd)
idx = ws.map_ids[a.ids[ip[i]]]
nd = ws.fluc[idx].nd
sid = truncate_ascii(get(names, a.ids[ip[i]], string(a.ids[ip[i]])), ntrunc)
if (length(a.cfd[ip[i]].gamm) > 0)
idx = ws.map_ids[a.ids[i]]
nd = ws.fluc[idx].nd
if (nd > 1)
Printf.@printf(" # %45s %6.2f %10d\n",
sid, 100.0 .* a.cfd[ip[i]].var ./ a.err^2, nd)
else
......
......@@ -14,13 +14,12 @@ function cobs(avgs::Vector{Float64}, cov::Array{Float64, 2}, ids::Vector{Int64})
ch = LinearAlgebra.cholesky(cov)
n = length(avgs)
p = Vector{uwreal}()
p = Vector{uwreal}(undef, n)
for j in 1:n
new = uwreal([avgs[j], ch.L[j, 1]], ids[1])
p[j] = uwreal([avgs[j], ch.L[j, 1]], ids[1])
for i in 2:n
new = new + uwreal([0.0, ch.L[j,i]], ids[i])
p[j] = p[j] + uwreal([0.0, ch.L[j,i]], ids[i])
end
push!(p, new)
end
return p
end
......@@ -55,7 +54,7 @@ function addobs(a::Vector{uwreal}, der::Array{Float64, 2}, mvl::Vector{Float64})
n = max(n,length(a[i].der))
end
x = Vector{uwreal}()
x = Vector{uwreal}(undef, size(der, 1))
for i in 1:size(der, 1)
p = [false for i in 1:n]
d = zeros(Float64, n)
......@@ -69,7 +68,7 @@ function addobs(a::Vector{uwreal}, der::Array{Float64, 2}, mvl::Vector{Float64})
end
end
end
push!(x, uwreal(mvl[i], p, d))
x[i] = uwreal(mvl[i], p, d)
end
return x
......
......@@ -10,6 +10,13 @@
###
mutable struct hyperd
v::Float64
d1::Float64
d2::Float64
d3::Float64
end
mutable struct cfdata
var::Float64
taui::Float64
......@@ -19,6 +26,15 @@ mutable struct cfdata
gamm::Vector{Float64}
drho::Vector{Float64}
function cfdata()
x = new()
x.var = 0.0
x.taui = 0.0
x.dtaui = 0.0
x.iw = 0
return x
end
end
mutable struct uwreal
......
......@@ -108,17 +108,24 @@ function fit_error(chisq::Function,
n = length(xp) # Number of fit parameters
m = length(data) # Number of data
ccsq(x::Vector) = chisq(x[1:n], x[n+1:end])
xav = zeros(Float64, n+m)
xav = Vector{Float64}(undef, n+m)
for i in 1:n
xav[i] = xp[i]
end
for i in n+1:n+m
xav[i] = data[i-n].mean
end
hess = ForwardDiff.hessian(ccsq, xav)
function cls(x0)
x1 = view(x0, 1:n)
x2 = view(x0, n+1:n+m)
return chisq(x1, x2)
end
hess = Array{Float64}(undef, n+m, n+m)
ForwardDiff.hessian!(hess, cls, xav)
hinv = LinearAlgebra.pinv(hess[1:n,1:n])
grad = - hinv[1:n,1:n] * hess[1:n,n+1:n+m]
......@@ -144,7 +151,7 @@ function fit_error(chisq::Function,
else
Ww = W
end
cse = chiexp(hess, data, Ww)
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