minor modification on local branch (macbook)

src/gevp.jl

@@ -12,7 +12,7 @@ rcParams["axes.labelsize"] =16
 plt.rc("text", usetex=true)
 ##
 #============= SET UP VARIABLES ===========#
-const path_data = "/Volumes/Macintosh HD/Lattice/data"
+const path_data = "/Users/ale/H101r001/dat"
 const path_plat = "/Users/ale/automation/plat.txt"
 const path_results = "/Users/ale/Desktop/results_gevp"
 const path_rw = "/Users/ale/Desktop/data/ms1_dat"
@@ -26,18 +26,19 @@ const bdio_rec = [0, 2, 4, 6, 8, 10, 12, 14, 16 ]
 const sector = Dict("ll"=>true, "ls"=>false, "lh"=>true, "ss"=>false, "sh"=>true, "hh"=>true )
 const tau = 3
 const _t0 = 2
-const range_t_fit = [2,5]
-const rwf = true
-const compute_t0 = true
-const mass_shift = true 
+const range_t_fit = [3,3]
+const rwf = false
+const compute_t0 = false
+const mass_shift = false 
 #go to line 443 to change CL+chiral extrapolation fit model
 
 #========= DO NOT CHANGE AFTER THIS LINE ==========#
+##
 include("types.jl")
 include("tools.jl")
 include("const.jl")
 include("func_comb.jl")
-
+##
 #======== GET ENSEMBLE INFORMATION FROM DATABASE ==========#
 ensinfo = Vector{EnsInfo}(undef, length(ens_list))
 for i in 1:length(ens_list)
@@ -74,6 +75,8 @@ end
 
 ##
 
+
+##
 wpmm = Dict{Int64, Vector{Float64}}()
 wpmm[303] = [-1.0, 2.0,-1.0, -1.0]
 wpmm[300] = [-1.0, 1.5,-1.0, -1.0]
@@ -733,13 +736,6 @@ function plot_hist(cat::Vector{Cat} ; save::Bool=false)
     println(best)
     all_res_aux = vcat(getfield.(cat, :fit_res) .* hc... )
     all_res = [all_res_aux[i] /t0_ph[1] for i =1:length(all_res_aux)]
-    for i in 1:length(all_res)
-        if best -10 > all_res[i]
-            diff = best - all_res[i]  
-            rn=rand(60:130)
-            all_res[i] = all_res[i] + rn*diff/100
-        end
-    end
     println(length(all_res))
     uwerr.(all_res)
     uwerr(best)
@@ -771,15 +767,6 @@ function plot_cat(cat::Cat, ylab::LaTeXString ; cut::Float64=2.5, save::Bool=fal
     all_res_aux = getfield(cat, :fit_res) .* hc 
     all_res_aux = all_res_aux[aic_norm_index]
     all_res = [all_res_aux[i] /t0_ph[1] for i =1:length(all_res_aux)]
-    for i in 1:length(all_res)
-        if best -10 > all_res[i]
-            diff = best - all_res[i]  
-            println(diff)
-            rn=rand(60:130)
-            all_res[i] = all_res[i] + rn*diff/100
-            aic[i] -= minimum(aic)/(value(diff)*rn)*10
-        end
-    end
     uwerr.(all_res)
     uwerr(best)
     fig = figure()

src/tools.jl

@@ -51,7 +51,7 @@ function get_corr(ens::EnsInfo, g1::String="G5", g2::String="G5"; rw::Bool=false
         truncate_data!(aux, 739)
     end
     if ens.id =="H101"
-        truncate_data!(aux, [878, 999])
+        #truncate_data!(aux, [878, 999])
     end
     !rw ? obs = corr_obs.(aux, L=ens.L) : obs = corr_obs.(aux, L=ens.L, rw=read_rw(ens.id))
     return obs
@@ -146,6 +146,20 @@ function comp_mat(ensinfo::EnsInfo, tot_obs::Vector{juobs.Corr}, tau::Int64)
     end
     return mat_info
 end
+function comp_mat_lukas(ensinfo::EnsInfo, tot_obs::Vector{juobs.Corr}, dim::Int64)
+    mat_info = Vector{MatInfo}(undef, length(tot_obs))
+    for i in 1:length(tot_obs)
+        a11 = tot_obs[i].obs[1+1:end-6]
+        a22 = tot_obs[i].obs[3+1:end-4]
+        a33 = tot_obs[i].obs[5+1:end-2]
+        a12 = tot_obs[i].obs[2+1:end-5]
+        a13 = tot_obs[i].obs[3+1:end-4]
+        a23 = tot_obs[i].obs[4+1:end-3]
+        aux = juobs.get_matrix([a11, a22, a33], [a12, a13, a23])
+        mat_info[i] = MatInfo(ensinfo, aux, tot_obs[i].mu, tot_obs[i].y0)        
+    end
+    return mat_info
+end
 function gevp_mass(mat_obs::Vector{MatInfo}; t0::Int64=2, delta_t_in::Array{Int64}=[2,5], wpm::Union{Nothing, Dict{Int64,Vector{Float64}}}=nothing)
     y0s = getfield.(mat_obs, :y0)
     mu_list = getfield.(mat_obs, :mu) 
@@ -154,14 +168,15 @@ function gevp_mass(mat_obs::Vector{MatInfo}; t0::Int64=2, delta_t_in::Array{Int6
     for i in 1:length(mat_obs)
         mat = getfield(mat_obs[i], :mat_list)[y0s[i]+1:end-1] #matrices to diagonalise for a given sector [i]
         evals = getall_eigvals(mat, t0)
-        en = energies(evals) #ground and excited state energy
+        en = energies(evals, wpm=wpm) #ground and excited state energy
         aux_mass = Vector{uwreal}(undef,0)
         for t in delta_t_in[1]:delta_t_in[2]
             try
                 temp = fit_mass(en[1],t_in=t,  wpm=wpm)
                 push!(aux_mass, temp)
             catch
-                plat = select_plateau(mat_obs[i].ensinfo, mu_list)[end-2] .- y0s[i]
+                plat = select_plateau(mat_obs[i].ensinfo, mu_list)[end-2] .- y0s[i] 
+                println(plat)
                 temp = plat_av(en[1], plat)
                 uwerr(temp)
                 push!(aux_mass, temp)
@@ -173,15 +188,15 @@ function gevp_mass(mat_obs::Vector{MatInfo}; t0::Int64=2, delta_t_in::Array{Int6
     end
     return plat_en0
 end
-function fit_mass(en_i::Array{uwreal}; t_in::Int64=3, wpm::Union{Nothing, Dict{Int64,Vector{Float64}}}=nothing)
+function fit_mass(en_i::Array{uwreal}; t_in::Int64=3, t_fin::Int64=length(en_i)-1, wpm::Union{Nothing, Dict{String,Vector{Float64}}, Dict{Int64,Vector{Float64}}}=nothing)
     @.model(x,p) = p[1] + p[2]*exp(-x*p[3])
-    ydata = en_i[t_in:end-5]
+    ydata = en_i[t_in:t_fin]
     xdata = collect(1:length(ydata))
     par, _ = fit_routine(model, xdata, ydata, 3, wpm=wpm)
     println("\nPar[1] is ", par[1] )
     return par[1]
 end
-function gevp_dec(mat_obs::Vector{MatInfo}, mass::Vector{uwreal}; t0::Int64=2, wpm::Union{Nothing, Dict{Int64,Vector{Float64}}}=nothing, pl::Bool=true, n::Int64=1, pseudo::Bool=true)
+function gevp_dec(mat_obs::Vector{MatInfo}, mass::Vector{uwreal}; t0::Int64=2, wpm::Union{Nothing, Dict{Int64,Vector{Float64}}}=nothing, pl::Bool=true, n::Int64=1, pseudo::Bool=true, wilson::Bool=false)
     y0s = getfield.(mat_obs, :y0)
     mu_list = getfield.(mat_obs, :mu) 
     plat_dec0 = Vector{uwreal}(undef, length(mat_obs))
@@ -192,7 +207,7 @@ function gevp_dec(mat_obs::Vector{MatInfo}, mass::Vector{uwreal}; t0::Int64=2, w
         plat = select_plateau(mat_obs[i].ensinfo, mu_list)[end] .- y0s[i]
         try
             if pseudo
-                plat_dec0[i] = -dec(elem, plat, mass[i], mu_list[i], pl=pl)
+                plat_dec0[i] = -dec(elem, plat, mass[i], mu_list[i], pl=pl, wilson=wilson)
             else
                 plat_dec0[i] = -vec_dec(elem, plat, mass[i], pl=pl)
             end
@@ -204,7 +219,30 @@ function gevp_dec(mat_obs::Vector{MatInfo}, mass::Vector{uwreal}; t0::Int64=2, w
     end
     return plat_dec0
 end
-function getall_eig(a::Vector{Matrix{uwreal}}, t0; iter=30)
+#this second function accepts plateau
+function gevp_dec(mat_obs::Vector{MatInfo}, mass::Vector{uwreal}, plat::Vector{Vector{Int64}}; t0::Int64=2, wpm::Union{Nothing, Dict{Int64,Vector{Float64}}}=nothing, pl::Bool=true, n::Int64=1, pseudo::Bool=true, wilson::Bool=false)
+    y0s = getfield.(mat_obs, :y0)
+    mu_list = getfield.(mat_obs, :mu) 
+    plat_dec0 = Vector{uwreal}(undef, length(mat_obs))
+    for i in 1:length(mat_obs)
+        mat = getfield(mat_obs[i], :mat_list)[y0s[i]+2:end-1] #matrices to diagonalise for a given sector [i]
+        evecs = getall_eig(mat, t0)
+        elem = mat_elem(evecs, mat, mass[i], n) #ground and excited state energy
+        try
+            if pseudo
+                plat_dec0[i] = -dec(elem, plat[i], mass[i], mu_list[i], pl=pl, wilson=wilson)
+            else
+                plat_dec0[i] = -vec_dec(elem, plat[i], mass[i], pl=pl)
+            end
+        catch
+            println("Decay constant for the ensemble ", mat_obs[i].ensinfo.id, " failed in the sector ", mu_list[i], " pseudo sector?", pseudo)
+            println("The associated effective mass is ", mass[i])
+            plat_dec0[i] = uwreal(0)
+        end
+    end
+    return plat_dec0
+end
+function getall_eig(a::Vector{Matrix{uwreal}}, t0; iter=10)
     n = length(a)
     res = Vector{Matrix{uwreal}}(undef, n)
     [res[i] = uweigvecs(a[i], a[t0]) for i=1:n]
@@ -220,7 +258,7 @@ function mat_elem(evec::Vector{Matrix{uwreal}}, ct_mat::Vector{Matrix{uwreal}},
     end
     return res
 end
-function dec(mat_elem::Vector{uwreal}, plat::Array{Int64}, mass::uwreal, mu::Vector{Float64}; pl::Bool=true)
+function dec(mat_elem::Vector{uwreal}, plat::Array{Int64}, mass::uwreal, mu::Vector{Float64}; pl::Bool=true, wilson::Bool=false)
     aux = plat_av(mat_elem, plat)
     if pl
         uwerr(aux)
@@ -232,9 +270,15 @@ function dec(mat_elem::Vector{uwreal}, plat::Array{Int64}, mass::uwreal, mu::Vec
         figure()
         errorbar(x, y, yerr=dy, fmt="*", color="red")
         fill_between(plat[1]:plat[2], v-e, v+e, color="green")    
+        ylim(v - 10*e, v + 10*e)
         display(gcf())
+        close()
+    end
+    if !wilson
+        return sqrt(2/ (mass)^3) * sum(mu)* aux
+    else
+        return sqrt(2/mass) * aux 
     end
-    return sqrt(2/ (mass)^3) * sum(mu)* aux 
 end
 function vec_dec(mat_elem::Vector{uwreal}, plat::Array{Int64}, mass::uwreal; pl::Bool=true)
     aux = plat_av(mat_elem, plat)
@@ -246,6 +290,7 @@ function vec_dec(mat_elem::Vector{uwreal}, plat::Array{Int64}, mass::uwreal; pl:
         v = value(aux)
         e = err(aux)
         figure()
+        ylim(v - 10*e, v + 10*e)
         errorbar(x, y, yerr=dy, fmt="*", color="red")
         fill_between(plat[1]:plat[2], v-e, v+e, color="green")    
         display(gcf())