Commit 7c0e8824 by Javier

### Merge branch 'javier'

parents b62250a8 55a5dc8b
This diff is collapsed.
 ... ... @@ -7,7 +7,6 @@ version = "0.1.0" ADerrors = "5e92007d-7bf1-471c-8ceb-4591b8b567a9" BDIO = "375f315e-f2c4-11e9-2ef9-134f02f79e27" LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" LeastSquaresOptim = "0fc2ff8b-aaa3-5acd-a817-1944a5e08891" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LsqFit = "2fda8390-95c7-5789-9bda-21331edee243" Optim = "429524aa-4258-5aef-a3af-852621145aeb" ... ...
 Home · juobs Documentation
Home · juobs Documentation
 ... ... @@ -78,4 +78,4 @@ evals = getall_eigvals(matrices, 5) #where t_0=5 Julia>source
juobs.getall_eigvecsFunction
getall_eigvecs(a::Vector{Matrix}, delta_t; iter=30 )

This function solves a GEVP problem, returning the eigenvectors, for a list of matrices.

$C(t_i)v_i = λ_i C(t_i-\delta_t) v_i$, with i=1:lenght(a)

Here delta_t is the time shift within the two matrices of the problem, and is kept fixed. It takes as input:

• a::Vector{Matrix} : a vector of matrices

• delta_t::Int64 : the fixed time shift t-t_0

• iter=30 : the number of iterations of the qr algorithm used to extract the eigenvalues

It returns:

• res = Vector{Matrix{uwreal}}

where each res[i] is a matrix with the eigenvectors as columns Examples:

mat_array = get_matrix(diag, upper_diag)
evecs = getall_eigvecs(mat_array, 5)
source
evecs = getall_eigvecs(mat_array, 5)source
This diff is collapsed.
 ... ... @@ -24,4 +24,4 @@ truncate_data!(Y, nc) dat = read_mesons([path1, path2], "G5", "G5") Y = read_ms.([path1_ms, path2_ms]) truncate_data!(dat, [nc1, nc2]) truncate_data!(Y, [nc1, nc2])source truncate_data!(Y, [nc1, nc2])source
 Search · juobs Documentation

Search · juobs Documentation

 ... ... @@ -494,9 +494,9 @@ function get_model(x, p, n) return s end @doc raw""" comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Union{Matrix{Float64}, Nothing}=nothing, npol::Int64=2, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Union{Matrix{Float64}, Nothing}=nothing, npol::Int64=2, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, info::Bool=false) comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Union{Vector{Matrix{Float64}}, Nothing}=nothing, npol::Int64=2, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Union{Vector{Matrix{Float64}}, Nothing}=nothing, npol::Int64=2, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, info::Bool=false) Computes t0 using the energy density of the action Ysl(Yang-Mills action). t0 is computed in the plateau plat. ... ... @@ -504,13 +504,16 @@ A polynomial interpolation in t is performed to find t0, where npol is the The flag pl allows to show the plot. The flag info provides extra output that contains information about the primary observables. The function returns the primary observables  and  (it returns the observable if rw=nothing) @example #Single replica Y = read_ms(path) rw = read_ms(path_rw) t0 = comp_t0(Y, [38, 58], L=32) t0_r = comp_t0(Y, [38, 58], L=32, rw=rw) t0, Yobs = comp_t0(Y, [38, 58], L=32, info=true) t0_r, WYobs, Wobs = comp_t0(Y, [38, 58], L=32, rw=rw, info=true) #Two replicas Y1 = read_ms(path1) ... ... @@ -525,7 +528,7 @@ t0_r = comp_t0(Y, [38, 58], L=32, rw=[rw1, rw2], pl=true) """ function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Union{Matrix{Float64}, Nothing}=nothing, npol::Int64=2, ws::ADerrors.wspace=ADerrors.wsg, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, info::Bool=false) Ysl = Y.obs t = Y.t ... ... @@ -551,16 +554,27 @@ function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, end end Ysl = isnothing(rw) ? Ysl : apply_rw(Ysl, rw) xmax = size(Ysl, 2) nt0 = t0_guess(t, Ysl, plat, L) dt0 = iseven(npol) ? Int64(npol / 2) : Int64((npol+1)/ 2) Y_aux = Matrix{uwreal}(undef, xmax, 2*dt0+1) if !isnothing(rw) Ysl_r, W = apply_rw(Ysl, rw) W_obs = uwreal(W, id) WY_aux = Matrix{uwreal}(undef, xmax, 2*dt0+1) end for i = 1:xmax k = 1 for j = nt0-dt0:nt0+dt0 if isnothing(rw) Y_aux[i, k] = uwreal(Ysl[:, i, j], id) else WY_aux[i, k] = uwreal(Ysl_r[:, i, j], id) Y_aux[i, k] = WY_aux[i, k] / W_obs end k = k + 1 end end ... ... @@ -601,12 +615,18 @@ function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, title(string(L"$t/a^2 =$", t[nt0])) display(gcf()) end if info && !isnothing(rw) return (t0, WY_aux, W_obs) elseif info && isnothing(rw) return (t0, Y_aux) else return t0 end end function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Union{Vector{Matrix{Float64}}, Nothing}=nothing, npol::Int64=2, ws::ADerrors.wspace=ADerrors.wsg, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, info::Bool=false) nr = length(Y) Ysl = getfield.(Y, :obs) ... ... @@ -641,18 +661,32 @@ function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false end end Ysl = isnothing(rw) ? Ysl : apply_rw(Ysl, rw) tmp = Ysl[1] [tmp = cat(tmp, Ysl[k], dims=1) for k = 2:nr] [tmp = cat(tmp, Ysl[k], dims=1) for k=2:nr] nt0 = t0_guess(t, tmp, plat, L) xmax = size(tmp, 2) dt0 = iseven(npol) ? Int64(npol / 2) : Int64((npol+1) / 2) Y_aux = Matrix{uwreal}(undef, xmax, 2*dt0+1) if !isnothing(rw) Ysl_r, W = apply_rw(Ysl, rw) tmp_r = Ysl_r[1] tmp_W = W[1] [tmp_r = cat(tmp_r, Ysl_r[k], dims=1) for k = 2:nr] [tmp_W = cat(tmp_W, W[k], dims=1) for k = 2:nr] W_obs = uwreal(tmp_W, id, replica) WY_aux = Matrix{uwreal}(undef, xmax, 2*dt0+1) end for i = 1:xmax k = 1 for j = nt0-dt0:nt0+dt0 if isnothing(rw) Y_aux[i, k] = uwreal(tmp[:, i, j], id, replica) else WY_aux[i, k] = uwreal(tmp_r[:, i, j], id, replica) Y_aux[i, k] = WY_aux[i, k] / W_obs end k = k + 1 end end ... ... @@ -693,7 +727,13 @@ function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false title(string(L"$t/a^2 =$", t[nt0])) display(gcf()) end if info && !isnothing(rw) return (t0, WY_aux, W_obs) elseif info && isnothing(rw) return (t0, Y_aux) else return t0 end end function t0_guess(t::Vector{Float64}, Ysl::Array{Float64, 3}, plat::Vector{Int64}, L::Int64) ... ...
 ... ... @@ -101,7 +101,8 @@ read_mesons(path, "G5", "G5", id="H100") read_mesons(path, "G5_d2", "G5_d2", legacy=true)  """ function read_mesons(path::String, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{String, Nothing}=nothing, legacy::Bool=false) function read_mesons(path::String, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{String, Nothing}=nothing, legacy::Bool=false, nnoise_trunc::Union{Int64, Nothing}=nothing) t1 = isnothing(g1) ? nothing : findfirst(x-> x==g1, gamma_name) - 1 t2 = isnothing(g2) ? nothing : findfirst(x-> x==g2, gamma_name) - 1 if isnothing(id) ... ... @@ -119,6 +120,8 @@ function read_mesons(path::String, g1::Union{String, Nothing}=nothing, g2::Union tvals = g_header.tvals nnoise = g_header.nnoise nnoise_trunc = isnothing(nnoise_trunc) ? nnoise : min(nnoise, nnoise_trunc) fsize = filesize(path) datsize = 4 + sum(getfield.(c_header, :dsize)) * tvals * nnoise #data_size / ncnfg ... ... @@ -144,13 +147,13 @@ function read_mesons(path::String, g1::Union{String, Nothing}=nothing, g2::Union tmp = Array{Float64}(undef, tvals*nnoise) read!(data, tmp) tmp2 = reshape(tmp, (nnoise, tvals)) tmp2 = mean(tmp2, dims=1) tmp2 = mean(tmp2[1:nnoise_trunc, :], dims=1) data_re[c, icfg, :] = tmp2[1, :] elseif c_header[k].is_real == 0 tmp = Array{Float64}(undef, 2*tvals*nnoise) read!(data, tmp) tmp2 = reshape(tmp, (2, nnoise, tvals)) tmp2 = mean(tmp2, dims=2) tmp2 = mean(tmp2[:, 1:nnoise_trunc, :], dims=2) data_re[c, icfg, :] = tmp2[1, 1, :] data_im[c, icfg, :] = tmp2[2, 1, :] ... ... @@ -172,8 +175,102 @@ function read_mesons(path::String, g1::Union{String, Nothing}=nothing, g2::Union return res end function read_mesons(path::Vector{String}, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{String, Nothing}=nothing, legacy::Bool=false) res = read_mesons.(path, g1, g2, id=id, legacy=legacy) function read_mesons(path::Vector{String}, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{String, Nothing}=nothing, legacy::Bool=false, nnoise_trunc::Union{Int64, Nothing}=nothing) res = read_mesons.(path, g1, g2, id=id, legacy=legacy, nnoise_trunc=nnoise_trunc) nrep = length(res) ncorr = length(res[1]) cdata = Vector{Vector{CData}}(undef, ncorr) for icorr = 1:ncorr cdata[icorr] = Vector{CData}(undef, nrep) for r = 1:nrep cdata[icorr][r] = res[r][icorr] end end return cdata end function read_mesons_correction(path::String, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{String, Nothing}=nothing, legacy::Bool=false, nnoise_trunc::Union{Int64, Nothing}=nothing) t1 = isnothing(g1) ? nothing : findfirst(x-> x==g1, gamma_name) - 1 t2 = isnothing(g2) ? nothing : findfirst(x-> x==g2, gamma_name) - 1 if isnothing(id) bname = basename(path) m = findfirst(r"[A-Z][0-9]{3}r[0-9]{3}", bname) id = bname[m[1:4]] #id = parse(Int64, bname[m[2:4]]) end data = open(path, "r") g_header = read_GHeader(path) c_header = read_CHeader(path, legacy=legacy) ncorr = g_header.ncorr tvals = g_header.tvals nnoise = g_header.nnoise nnoise_trunc = isnothing(nnoise_trunc) ? nnoise : min(nnoise, nnoise_trunc) fsize = filesize(path) datsize = 4 + sum(getfield.(c_header, :dsize)) * tvals * nnoise #data_size / ncnfg ncfg = div(fsize - g_header.hsize - sum(getfield.(c_header, :hsize)), datsize) #(total size - header_size) / data_size corr_match = findall(x-> (x.type1==t1 || isnothing(t1)) && (x.type2==t2 || isnothing(t2)), c_header) seek(data, g_header.hsize + sum(c.hsize for c in c_header)) res = Array{CData}(undef, div(length(corr_match), 2)) # Modification: total length is divided by 2 data_re = zeros(div(length(corr_match), 2), ncfg, tvals) # Modification: total length is divided by 2 data_im = zeros(div(length(corr_match), 2), ncfg, tvals) # Modification: total length is divided by 2 vcfg = Array{Int32}(undef, ncfg) for icfg = 1:ncfg vcfg[icfg] = read(data, Int32) c = 1 sgn = +1 # sign it changes at ncorr / 2. O_exact - O_sloppy for k = 1:ncorr if k in corr_match if c_header[k].is_real == 1 tmp = Array{Float64}(undef, tvals*nnoise) read!(data, tmp) tmp2 = reshape(tmp, (nnoise, tvals)) tmp2 = mean(tmp2[1:nnoise_trunc, :], dims=1) data_re[c, icfg, :] = data_re[c, icfg, :] + sgn * tmp2[1, :] elseif c_header[k].is_real == 0 tmp = Array{Float64}(undef, 2*tvals*nnoise) read!(data, tmp) tmp2 = reshape(tmp, (2, nnoise, tvals)) tmp2 = mean(tmp2[:, 1:nnoise_trunc, :], dims=2) data_re[c, icfg, :] = data_re[c, icfg, :] + sgn * tmp2[1, 1, :] data_im[c, icfg, :] = data_im[c, icfg, :] + sgn * tmp2[2, 1, :] end c += 1 else seek(data, position(data) + c_header[k].dsize*tvals*nnoise) end if k == div(ncorr, 2) c = 1 sgn = -1 end end end for c = 1:div(length(corr_match), 2) res[c] = CData(c_header[corr_match[c]], vcfg, data_re[c, :, :], data_im[c, :, :], id) end close(data) return res end function read_mesons_correction(path::Vector{String}, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{String, Nothing}=nothing, legacy::Bool=false, nnoise_trunc::Union{Int64, Nothing}=nothing) res = read_mesons_correction.(path, g1, g2, id=id, legacy=legacy, nnoise_trunc=nnoise_trunc) nrep = length(res) ncorr = length(res[1]) ... ...