Commit c8b1159d by Alejandro Saez

### Merge branch 'javier' into alejandro

parents 7f6a1be6 5eb61b12
 Home · juobs Documentation
Home · juobs Documentation
 ... @@ -78,4 +78,4 @@ evals = getall_eigvals(matrices, 5) #where t_0=5 ... @@ -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)                          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
 ... @@ -18,12 +18,19 @@ ca = -0.006033 * g2 *( 1 + exp(p0 + p1/g2)) ... @@ -18,12 +18,19 @@ ca = -0.006033 * g2 *( 1 + exp(p0 + p1/g2)) m12 = mpcac(corr_a0p, corr_pp, [50, 60], pl=false, ca=ca)source
juobs.dec_constFunction
dec_const(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, y0::Int64; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         m12 = mpcac(corr_a0p, corr_pp, [50, 60], pl=false, ca=ca)
source
juobs.dec_constFunction
dec_const(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, y0::Int64; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)

dec_const(a0p::Corr, pp::Corr, plat::Vector{Int64}, m::uwreal; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)

Computes the bare decay constant using $A_0P$ and $PP$ correlators . The decay constant is computed in the plateau plat. Correlator can be passed as an Corr struct or Vector{uwreal}. If it is passed as a uwreal vector, effective mass m and source position y0 must be specified.

The flags pl and data allow to show the plots and return data as an extra result. The ca variable allows to compute dec_const using the improved axial current.

The method assumes that the source is close to the boundary. It takes the following ratio to cancel boundary effects. $R = \frac{f_A(x_0, y_0)}{\sqrt{f_P(T-y_0, y_0)}} * e^{m (x_0 - T/2)}$

data_pp = read_mesons(path, "G5", "G5", legacy=true)                                                                                                                                                                                                                       dec_const(a0p::Corr, pp::Corr, plat::Vector{Int64}, m::uwreal; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)

dec_const(a0pL::Vector{uwreal}, a0pR::Vector{uwreal}, ppL::Vector{uwreal}, ppR::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, y0::Int64; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)

dec_const(a0pL::Corr, a0pR::Corr, ppL::Corr, ppR::Corr, plat::Vector{Int64}, m::uwreal; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)

Computes the bare decay constant using $A_0P$ and $PP$ correlators . The decay constant is computed in the plateau plat. Correlator can be passed as an Corr struct or Vector{uwreal}. If it is passed as a uwreal vector, effective mass m and source position y0 must be specified.

The flags pl and data allow to show the plots and return data as an extra result. The ca variable allows to compute dec_const using the improved axial current.

The method assumes that the source is close to the boundary. It takes the following ratio to cancel boundary effects. $R = \frac{f_A(x_0, y_0)}{\sqrt{f_P(T-y_0, y_0)}} * e^{m (x_0 - T/2)}$

If left and right correlators are included in the input. The result is computed with the following ratio $R = \sqrt{f_A(x_0, y_0) * f_A(x_0, T - 1 - y_0) / f_P(T - 1 - y_0, y_0)}$

data_pp = read_mesons(path, "G5", "G5", legacy=true)
data_a0p = read_mesons(path, "G5", "G0G5", legacy=true)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       data_a0p = read_mesons(path, "G5", "G0G5", legacy=true)

corr_pp = corr_obs.(data_pp, L=32)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                corr_pp = corr_obs.(data_pp, L=32)
corr_a0p = corr_obs.(data_a0p, L=32)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                              corr_a0p = corr_obs.(data_a0p, L=32)

corr_a0pL, corr_a0pR = [corr_a0p[1], corr_a0p[2]]
corr_ppL, corr_ppR = [corr_pp[1], corr_pp[2]]

m = meff(corr_pp[1], [50, 60], pl=false)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          m = meff(corr_pp[1], [50, 60], pl=false)

beta = 3.46                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       beta = 3.46
...  @@ -32,12 +39,21 @@ p1 = -13.9847                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            ...  @@ -32,12 +39,21 @@ p1 = -13.9847
g2 = 6 / beta                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     g2 = 6 / beta
ca = -0.006033 * g2 *( 1 + exp(p0 + p1/g2))                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       ca = -0.006033 * g2 *( 1 + exp(p0 + p1/g2))

f = dec_const(corr_a0p[1], corr_pp[1], [50, 60], m, pl=true, ca=ca)
source
juobs.dec_const_pcvcFunction
dec_const_pcvc(corr::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu::Vector{Float64}, y0::Int64 ; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         f = dec_const(corr_a0p[1], corr_pp[1], [50, 60], m, pl=true, ca=ca)

f_ratio = dec_const(corr_a0pL, corr_a0pR, corr_ppL, corr_ppR, [50, 60], m, pl=true, ca=ca)
source
juobs.dec_const_pcvcFunction
dec_const_pcvc(corr::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu::Vector{Float64}, y0::Int64 ; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)

dec_const_pcvc(corr::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)

Computes decay constant using the PCVC relation for twisted mass fermions. The decay constant is computed in the plateau plat. Correlator can be passed as an Corr struct or Vector{uwreal}. If it is passed as a uwreal vector, vector of twisted masses mu and source position y0 must be specified.

The flags pl and data allow to show the plots and return data as an extra result.

The method assumes that the source is in the bulk.

data = read_mesons(path, "G5", "G5")                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    dec_const_pcvc(corr::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)

dec_const_pcvc(ppL::Vector{uwreal}, ppR::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu::Vector{Float64}, y0::Int64 ; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)

dec_const_pcvc(corrL::Corr, corrR::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)

Computes decay constant using the PCVC relation for twisted mass fermions. The decay constant is computed in the plateau plat. Correlator can be passed as an Corr struct or Vector{uwreal}. If it is passed as a uwreal vector, vector of twisted masses mu and source position y0 must be specified.

The flags pl and data allow to show the plots and return data as an extra result.

The method extract the matrix element assuming that the source is in the bulk. ** **If left and right correlators are included in the input. The result is computed with a ratio that cancels boundary effects: $R = \sqrt{f_P(x_0, y_0) * f_P(x_0, T - 1 - y_0) / f_P(T - 1 - y_0, y_0)}$

data = read_mesons(path, "G5", "G5")
corr_pp = corr_obs.(data, L=32)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                   corr_pp = corr_obs.(data, L=32)
m = meff(corr_pp[1], [50, 60], pl=false)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          m = meff(corr_pp[1], [50, 60], pl=false)
f = dec_const_pcvc(corr_pp[1], [50, 60], m, pl=false)
source
juobs.comp_t0Function
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)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                      f = dec_const_pcvc(corr_pp[1], [50, 60], m, pl=false)

#left and right correlators
f_ratio = dec_const_pcvc(ppL, ppR, [50, 60], m)
source
juobs.comp_t0Function
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::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)

Computes t0 using the energy density of the action Ysl(Yang-Mills action). t0 is computed in the plateau plat. A polynomial interpolation in t is performed to find t0, where npol is the degree of the polynomial (linear fit by default)

The flag pl allows to show the plot.

#Single replica                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  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)

Computes t0 using the energy density of the action Ysl(Yang-Mills action). t0 is computed in the plateau plat. A polynomial interpolation in t is performed to find t0, where npol is the degree of the polynomial (linear fit by default)

The flag pl allows to show the plot.

#Single replica
...  @@ -54,4 +70,4 @@ rw2 = read_ms(path_rw2)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    ...  @@ -54,4 +70,4 @@ rw2 = read_ms(path_rw2)

t0 = comp_t0([Y1, Y2], [38, 58], L=32, pl=true)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                   t0 = comp_t0([Y1, Y2], [38, 58], L=32, pl=true)
t0_r = comp_t0(Y, [38, 58], L=32, rw=[rw1, rw2], pl=true)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         t0_r = comp_t0(Y, [38, 58], L=32, rw=[rw1, rw2], pl=true)

source
source
 ... @@ -24,4 +24,4 @@ truncate_data!(Y, nc) ... @@ -24,4 +24,4 @@ truncate_data!(Y, nc) dat = read_mesons([path1, path2], "G5", "G5") dat = read_mesons([path1, path2], "G5", "G5") Y = read_ms.([path1_ms, path2_ms]) Y = read_ms.([path1_ms, path2_ms]) truncate_data!(dat, [nc1, nc2]) truncate_data!(dat, [nc1, nc2]) truncate_data!(Y, [nc1, nc2])source truncate_data!(Y, [nc1, nc2])source
 Search · juobs Documentation

Search · juobs Documentation

 ... @@ -50,4 +50,4 @@ m2_phys = fitp[1] + fitp[2] * phi2_phys

Given a model function with a number param of parameters and an array of uwreal, this function fit ydata with the given model and print fit information The method return an array upar with the best fit parameters with their errors. The flag wpm is an optional array of Float64 of lenght 4. The first three paramenters specify the criteria to determine the summation windows:

• vp[1]: The autocorrelation function is summed up to $t = round(vp[1])$.

• vp[2]: The sumation window is determined using U. Wolff poposal with $S_\tau = wpm[2]$

• vp[3]: The autocorrelation function $\Gamma(t)$ is summed up a point where its error $\delta\Gamma(t)$ is a factor vp[3] times larger than the signal.

An additional fourth parameter vp[4], tells ADerrors to add a tail to the error with $\tau_{exp} = wpm[4]$. Negative values of wpm[1:4] are ignored and only one component of wpm[1:3] needs to be positive. If the flag covaris set to true, fit_routine takes into account covariances between x and y for each data point.

@. model(x,p) = p[1] + p[2] * exp(-(p[3]-p[1])*x)       fit_routine(model::Function, xdata::Array{uwreal}, ydata::Array{uwreal}, param::Int64=3; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, covar::Bool=false)

Given a model function with a number param of parameters and an array of uwreal, this function fit ydata with the given model and print fit information The method return an array upar with the best fit parameters with their errors. The flag wpm is an optional array of Float64 of lenght 4. The first three paramenters specify the criteria to determine the summation windows:

• vp[1]: The autocorrelation function is summed up to $t = round(vp[1])$.

• vp[2]: The sumation window is determined using U. Wolff poposal with $S_\tau = wpm[2]$

• vp[3]: The autocorrelation function $\Gamma(t)$ is summed up a point where its error $\delta\Gamma(t)$ is a factor vp[3] times larger than the signal.

An additional fourth parameter vp[4], tells ADerrors to add a tail to the error with $\tau_{exp} = wpm[4]$. Negative values of wpm[1:4] are ignored and only one component of wpm[1:3] needs to be positive. If the flag covaris set to true, fit_routine takes into account covariances between x and y for each data point.

@. model(x,p) = p[1] + p[2] * exp(-(p[3]-p[1])*x)
@. model2(x,p) = p[1] + p[2] * x[:, 1] + (p[3] + p[4] * x[:, 1]) * x[:, 2]                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                 @. model2(x,p) = p[1] + p[2] * x[:, 1] + (p[3] + p[4] * x[:, 1]) * x[:, 2]
fit_routine(model, xdata, ydata, param=3)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  fit_routine(model, xdata, ydata, param=3)
fit_routine(model, xdata, ydata, param=3, covar=true)
source fit_routine(model, xdata, ydata, param=3, covar=true)source
 ... @@ -126,7 +126,7 @@ function energies(evals::Union{Vector{Vector{uwreal}},Array{Array{uwreal}} }; wp ... @@ -126,7 +126,7 @@ function energies(evals::Union{Vector{Vector{uwreal}},Array{Array{uwreal}} }; wp ratio = evals[t][i] / evals[t+1][i] ratio = evals[t][i] / evals[t+1][i] aux_en[t] = 0.5*log(ratio * ratio) aux_en[t] = 0.5*log(ratio * ratio) end end isnothing(wpm) ? uwerr.(aux_en) : uwerr.(aux_en, wpm) isnothing(wpm) ? uwerr.(aux_en) : [uwerr(a, wpm) for a in aux_en] eff_en[i] = copy(aux_en) eff_en[i] = copy(aux_en) end end return eff_en return eff_en ... ...
 ... @@ -313,18 +313,27 @@ end ... @@ -313,18 +313,27 @@ end dec_const_pcvc(corr::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) dec_const_pcvc(corr::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) dec_const_pcvc(ppL::Vector{uwreal}, ppR::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu::Vector{Float64}, y0::Int64 ; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) dec_const_pcvc(corrL::Corr, corrR::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) Computes decay constant using the PCVC relation for twisted mass fermions. The decay constant is computed in the plateau plat. Computes decay constant using the PCVC relation for twisted mass fermions. The decay constant is computed in the plateau plat. Correlator can be passed as an Corr struct or Vector{uwreal}. If it is passed as a uwreal vector, vector of twisted masses mu and source position y0 Correlator can be passed as an Corr struct or Vector{uwreal}. If it is passed as a uwreal vector, vector of twisted masses mu and source position y0 must be specified. must be specified. The flags pl and data allow to show the plots and return data as an extra result. The flags pl and data allow to show the plots and return data as an extra result. **The method assumes that the source is in the bulk.** **The method extract the matrix element assuming that the source is in the bulk. ** **If left and right correlators are included in the input. The result is computed with a ratio that cancels boundary effects:**  R = \sqrt{f_P(x_0, y_0) * f_P(x_0, T - 1 - y_0) / f_P(T - 1 - y_0, y_0)} @example @example data = read_mesons(path, "G5", "G5") data = read_mesons(path, "G5", "G5") corr_pp = corr_obs.(data, L=32) corr_pp = corr_obs.(data, L=32) m = meff(corr_pp[1], [50, 60], pl=false) m = meff(corr_pp[1], [50, 60], pl=false) f = dec_const_pcvc(corr_pp[1], [50, 60], m, pl=false) f = dec_const_pcvc(corr_pp[1], [50, 60], m, pl=false) #left and right correlators f_ratio = dec_const_pcvc(ppL, ppR, [50, 60], m)   """ """ function dec_const_pcvc(corr::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu::Vector{Float64}, y0::Int64 ; pl::Bool=true, data::Bool=false, function dec_const_pcvc(corr::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu::Vector{Float64}, y0::Int64 ; pl::Bool=true, data::Bool=false, ... @@ -340,23 +349,18 @@ function dec_const_pcvc(corr::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu ... @@ -340,23 +349,18 @@ function dec_const_pcvc(corr::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu R_av = plat_av(R, plat, wpm) R_av = plat_av(R, plat, wpm) f = sqrt(2) * (mu[1] + mu[2]) *R_av / m^1.5 f = sqrt(2) * (mu[1] + mu[2]) *R_av / m^1.5 if pl if pl if isnothing(wpm) isnothing(wpm) ? uwerr(f) : uwerr(f, wpm) uwerr(f) isnothing(wpm) ? uwerr(R_av) : uwerr(R_av, wpm) uwerr(R_av) uwerr.(R) else uwerr(f, wpm) uwerr(R_av, wpm) uwerr.(R, wpm) end v = value(R_av) v = value(R_av) e = err(R_av) e = err(R_av) figure() figure() fill_between(plat[1]:plat[2], v-e, v+e, color="green", alpha=0.75) lbl = string(L"$af =$", sprint(show, f)) errorbar(1:length(R), value.(R), err.(R), fmt="x", color="black") fill_between(plat[1]:plat[2], v-e, v+e, color="green", alpha=0.75, label=L"$R$") errorbar(1:length(R), value.(R), err.(R), fmt="x", color="black", label=lbl) ylabel(L"$R$") ylabel(L"$R$") xlabel(L"$x_0$") xlabel(L"$x_0$") legend() title(string(L"$\mu_1 =$", mu[1], L" $\mu_2 =$", mu[2])) title(string(L"$\mu_1 =$", mu[1], L" $\mu_2 =$", mu[2])) display(gcf()) display(gcf()) end end ... @@ -373,6 +377,50 @@ function dec_const_pcvc(corr::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=tru ... @@ -373,6 +377,50 @@ function dec_const_pcvc(corr::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=tru return dec_const_pcvc(corr.obs, plat, m, corr.mu, corr.y0, pl=pl, data=data, wpm=wpm) return dec_const_pcvc(corr.obs, plat, m, corr.mu, corr.y0, pl=pl, data=data, wpm=wpm) end end function dec_const_pcvc(ppL::Vector{uwreal}, ppR::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu::Vector{Float64}, y0::Int64 ; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) T = length(ppL) f1 = [ppL[T - y0] for k = 1:T] R = ((ppL .* ppR ./ f1).^2).^(1/4) R = R[2:end-1] R_av = plat_av(R, plat, wpm) f = sqrt(2) * (mu[1] + mu[2]) * R_av / m^1.5 if pl isnothing(wpm) ? uwerr(f) : uwerr(f, wpm) isnothing(wpm) ? uwerr(R_av) : uwerr(R_av, wpm) v = value(R_av) e = err(R_av) figure() lbl = string(L"$af =$", sprint(show, f)) fill_between(plat[1]:plat[2], v-e, v+e, color="green", alpha=0.75, label=L"$R$") errorbar(1:length(R), value.(R), err.(R), fmt="x", color="black", label=lbl) ylabel(L"$R$") xlabel(L"$x_0$") legend() title(string(L"$\mu_1 =$", mu[1], L" $\mu_2 =$", mu[2])) display(gcf()) end if !data return f else return (f, uwdot(sqrt(2) * (mu[1] + mu[2]) / m^1.5, R)) end end function dec_const_pcvc(corrL::Corr, corrR::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) T = length(corrL.obs) if (corrL.kappa == corrR.kappa) && (corrL.mu == corrR.mu) && (corrL.y0 == T - 1 - corrR.y0) return dec_const_pcvc(corrL.obs, corrR.obs, plat, m, corrL.mu, corrL.y0, pl=pl, data=data, wpm=wpm) else error("y0, kappa or mu values does not match") end end #t0 #t0 function get_model(x, p, n) function get_model(x, p, n) s = 0.0 s = 0.0 ... @@ -413,7 +461,7 @@ t0_r = comp_t0(Y, [38, 58], L=32, rw=[rw1, rw2], pl=true) ... @@ -413,7 +461,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, 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, 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, rw_info::Bool=false) Ysl = Y.obs Ysl = Y.obs t = Y.t t = Y.t ... @@ -421,32 +469,45 @@ function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, ... @@ -421,32 +469,45 @@ function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, replica = size.([Ysl], 1) replica = size.([Ysl], 1) #Truncation #Truncation n_ws = findfirst(x-> x == ws.str2id[id], ws.map_nob) if id in keys(ADerrors.wsg.str2id) if !isnothing(n_ws) n_ws = findfirst(x-> x == ws.str2id[id], ws.map_nob) ivrep_ws = ws.fluc[n_ws].ivrep if !isnothing(n_ws) ivrep_ws = ws.fluc[n_ws].ivrep if length(ivrep_ws) != 1 if length(ivrep_ws) != 1 error("Different number of replicas") error("Different number of replicas") end end if replica[1] > ivrep_ws[1] if replica[1] > ivrep_ws[1] println("Automatic truncation in Ysl ", ivrep_ws[1], " / ", replica[1], ". R = 1") println("Automatic truncation in Ysl ", ivrep_ws[1], " / ", replica[1], ". R = 1") Ysl = Ysl[1:ivrep_ws[1], :, :] Ysl = Ysl[1:ivrep_ws[1], :, :] elseif replica[1] < ivrep_ws[1] elseif replica[1] < ivrep_ws[1] error("Automatic truncation failed. R = 1\nTry using truncate_data!") error("Automatic truncation failed. R = 1\nTry using truncate_data!") end end end end end Ysl = isnothing(rw) ? Ysl : apply_rw(Ysl, rw) xmax = size(Ysl, 2) xmax = size(Ysl, 2) nt0 = t0_guess(t, Ysl, plat, L) nt0 = t0_guess(t, Ysl, plat, L) dt0 = iseven(npol) ? Int64(npol / 2) : Int64((npol+1)/ 2) dt0 = iseven(npol) ? Int64(npol / 2) : Int64((npol+1)/ 2) Y_aux = Matrix{uwreal}(undef, xmax, 2*dt0+1) 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 for i = 1:xmax k = 1 k = 1 for j = nt0-dt0:nt0+dt0 for j = nt0-dt0:nt0+dt0 Y_aux[i, k] = uwreal(Ysl[:, i, j], id) 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 k = k + 1 end end end end ... @@ -464,7 +525,7 @@ function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, ... @@ -464,7 +525,7 @@ function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, uwerr.(t2E) uwerr.(t2E) else else uwerr(t0, wpm) uwerr(t0, wpm) uwerr.(t2E, wpm) [uwerr(t2E_aux, wpm) for t2E_aux in t2E] end end v = value.(t2E) v = value.(t2E) ... @@ -487,7 +548,11 @@ function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, ... @@ -487,7 +548,11 @@ function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, title(string(L"$t/a^2 =$", t[nt0])) title(string(L"$t/a^2 =$", t[nt0])) display(gcf()) display(gcf()) end end return t0 if rw_info && !isnothing(rw) (t0, WY_aux, W_obs) else return t0 end end end function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false, function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false, ... @@ -504,38 +569,55 @@ function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false ... @@ -504,38 +569,55 @@ function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false if !all(id .== id[1]) if !all(id .== id[1]) error("IDs are not equal") error("IDs are not equal") end end id = id[1] #Truncation #Truncation n_ws = findfirst(x-> x == ws.str2id[id[1]], ws.map_nob) if id in keys(ADerrors.wsg.str2id) if !isnothing(n_ws) n_ws = findfirst(x-> x == ws.str2id[id], ws.map_nob) ivrep_ws = ws.fluc[n_ws].ivrep if !isnothing(n_ws) ivrep_ws = ws.fluc[n_ws].ivrep if length(replica) != length(ivrep_ws) if length(replica) != length(ivrep_ws) error("Different number of replicas") error("Different number of replicas") end end for k = 1:length(replica) for k = 1:length(replica) if replica[k] > ivrep_ws[k] if replica[k] > ivrep_ws[k] println("Automatic truncation in Ysl ", ivrep_ws[k], " / ", replica[k], ". R = ", k) println("Automatic truncation in Ysl ", ivrep_ws[k], " / ", replica[k], ". R = ", k) Ysl[k] = Ysl[k][1:ivrep_ws[k], :, :] Ysl[k] = Ysl[k][1:ivrep_ws[k], :, :] elseif replica[k] < ivrep_ws[k] elseif replica[k] < ivrep_ws[k] error("Automatic truncation failed. R = ", replica[k], "\nTry using truncate_data!") error("Automatic truncation failed. R = ", replica[k], "\nTry using truncate_data!") end end end replica = size.(Ysl, 1) end end replica = size.(Ysl, 1) end end Ysl = isnothing(rw) ? Ysl : apply_rw(Ysl, rw) tmp = Ysl[1] 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) nt0 = t0_guess(t, tmp, plat, L) xmax = size(tmp, 2) xmax = size(tmp, 2) dt0 = iseven(npol) ? Int64(npol / 2) : Int64((npol+1) / 2) dt0 = iseven(npol) ? Int64(npol / 2) : Int64((npol+1) / 2) Y_aux = Matrix{uwreal}(undef, xmax, 2*dt0+1) 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 for i = 1:xmax k = 1 k = 1 for j = nt0-dt0:nt0+dt0 for j = nt0-dt0:nt0+dt0 Y_aux[i, k] = uwreal(tmp[:, i, j], id[1], replica) 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 k = k + 1 end end end end ... @@ -553,7 +635,7 @@ function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false ... @@ -553,7 +635,7 @@ function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false uwerr.(t2E) uwerr.(t2E) else else uwerr(t0, wpm) uwerr(t0, wpm) uwerr.(t2E, wpm) [uwerr(t2E_aux, wpm) for t2E_aux in t2E] end end v = value.(t2E) v = value.(t2E) ... @@ -576,11 +658,19 @@ function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false ... @@ -576,11 +658,19 @@ function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false title(string(L"$t/a^2 =$", t[nt0])) title(string(L"$t/a^2 =$", t[nt0])) display(gcf()) display(gcf()) end end return t0 if rw_info && !isnothing(rw) (t0, WY_aux, W_obs) else return t0 end end end function t0_guess(t::Vector{Float64}, Ysl::Array{Float64, 3}, plat::Vector{Int64}, L::Int64) function t0_guess(t::Vector{Float64}, Ysl::Array{Float64, 3}, plat::Vector{Int64}, L::Int64) t2E_ax = t.^2 .* mean(mean(Ysl[:, plat[1]:plat[2], :], dims=2), dims=1)[1, 1, :] / L^3 t2E_ax = t.^2 .* mean(mean(Ysl[:, plat[1]:plat[2], :], dims=2), dims=1)[1, 1, :] / L^3 #look for values t2E: t2E > 0.3 and 0.0 < t2E_ax < 0.3 if !(any(t2E_ax .> 0.3) && any((t2E_ax .< 0.3) .& (t2E_ax .> 0.0))) error("Error finding solutions. Check volume normalization!") end t0_aux = minimum(abs.(t2E_ax .- 0.3)) t0_aux = minimum(abs.(t2E_ax .- 0.3)) nt0 = findfirst(x-> abs(x - 0.3) == t0_aux, t2E_ax) #index closest value to t0 nt0 = findfirst(x-> abs(x - 0.3) == t0_aux, t2E_ax) #index closest value to t0 return nt0 return nt0 ... ...
<
 ... @@ -4,8 +4,8 @@ function apply_rw(data::Array{Float64}, W::Matrix{Float64}) ... @@ -4,8 +4,8 @@ function apply_rw(data::Array{Float64}, W::Matrix{Float64}) W1 = W[1, 1:nc] W1 = W[1, 1:nc] W2 = W[2, 1:nc] W2 = W[2, 1:nc] data_r = data .* W1 .* W2 / mean(W1 .* W2) data_r = data .* W1 .* W2 return data_r return (data_r, W1 .* W2) end end