bug fixing and performance improvement

src/gevp.jl

@@ -5,9 +5,9 @@ using Revise, juobs, BDIO, DelimitedFiles, ADerrors, PyPlot, LaTeXStrings, LsqFi
 const path_data = "/Users/ale/Desktop/data"
 const path_plat = "/Users/ale/automation/plat.txt"
 const path_results = "/Users/ale/Desktop/results_gevp"
-const ensembles = [ "J303"]
+const ensembles = [ "J303","N300", "H400"]
 const sector = Dict("ll"=>false, "ls"=>false, "lh"=>true, "ss"=>false, "sh"=>true, "hh"=>true )
-const tau = 4
+const tau = 3
 const _t0 = 2
 const range_t_fit = [2,5]
 const rwf = false
@@ -65,12 +65,16 @@ println("\n Computing observables \n")
         @time begin
             ens = ensinfo[i]
             println("\n Ensemble: ", ens.id)
+            println("\nComputing spectrum...")
             en0_ps = gevp_mass(matinfo[i], t0=_t0, delta_t_in=range_t_fit, wpm=wpmm)
             en0_vec = gevp_mass(matinfo_vec[i], t0=_t0, delta_t_in=range_t_fit, wpm=wpmm)
             println( "\nground  ps mass: ", en0_ps)
             println( "\nground  vec mass: ", en0_vec)
-            dec0_ps = gevp_dec(matinfo[i], en0_ps, t0=_t0, delta_t_in=range_t_fit, wpm=wpmm)
-            dec0_vec = gevp_dec(matinfo_vec[i], en0_vec, t0=_t0, delta_t_in=range_t_fit, wpm=wpmm)            
+            println("\nComputing decay constants...")
+            dec0_ps = gevp_dec(matinfo[i], en0_ps, t0=_t0, n=1, wpm=wpmm, pl=false)
+            dec0_vec = gevp_dec(matinfo_vec[i], en0_vec, t0=_t0, n=1, wpm=wpmm, pl=false, pseudo=false)            
+            println( "\nground  ps decays: ", dec0_ps)
+            println( "\nground  vec decays: ", dec0_vec)
             ensobs[i] = EnsObs(ens, gen_mulist[i], en0_ps, dec0_ps, en0_vec, dec0_vec)
             if ens.deg
                 ensobs[i].m_ls = ensobs[i].m_ll
@@ -143,8 +147,8 @@ println("\n Matching Procedure \n")
             if ens.ensinfo.deg
                 ens.m_sh_match = ens.m_lh_match
                 ens.m_sh_vec_match = ens.m_lh_vec_match
-                #ens.f_sh_match = ens.f_lh_match
-                #ens.f_sh_vec_match = ens.f_lh_vec_match
+                ens.f_sh_match = ens.f_lh_match
+                ens.f_sh_vec_match = ens.f_lh_vec_match
             else
                 # m_sh
                 par, chi2exp = lin_fit(muh, ens.m_sh)     

src/tools.jl

@@ -164,16 +164,20 @@ function fit_mass(en_i::Array{uwreal}; t_in::Int64=3, wpm::Union{Nothing, Dict{I
     println("\nPar[1] is ", par[1] )
     return par[1]
 end
-function gevp_dec(mat_obs::Vector{MatInfo}, mass::Vector{uwreal}; t0::Int64=2, delta_t_in::Array{Int64}=[2,5], wpm::Union{Nothing, Dict{Int64,Vector{Float64}}}=nothing)
+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)
     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], 1) #ground and excited state energy
+        elem = mat_elem(evecs, mat, mass[i], n) #ground and excited state energy
         plat = select_plateau(mat_obs[i].ensinfo, mu_list)[end] .- y0s[i]
-        plat_dec0[i] = -dec(elem, plat, mass[i], mu_list[i])
+        if pseudo
+            plat_dec0[i] = -dec(elem, plat, mass[i], mu_list[i], pl=pl)
+        else
+            plat_dec0[i] = -vec_dec(elem, plat, mass[i], pl=pl)
+        end
     end
     return plat_dec0
 end
@@ -193,10 +197,38 @@ 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})
+function dec(mat_elem::Vector{uwreal}, plat::Array{Int64}, mass::uwreal, mu::Vector{Float64}; pl::Bool=true)
     aux = plat_av(mat_elem, plat)
+    if pl
+        uwerr(aux)
+        x = 1:length(mat_elem)
+        y = value.(mat_elem)
+        dy = err.(mat_elem)
+        v = value(aux)
+        e = err(aux)
+        figure()
+        errorbar(x, y, yerr=dy, fmt="*", color="red")
+        fill_between(plat[1]:plat[2], v-e, v+e, color="green")    
+        display(gcf())
+    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)
+    if pl
+        uwerr(aux)
+        x = 1:length(mat_elem)
+        y = value.(mat_elem)
+        dy = err.(mat_elem)
+        v = value(aux)
+        e = err(aux)
+        figure()
+        errorbar(x, y, yerr=dy, fmt="*", color="red")
+        fill_between(plat[1]:plat[2], v-e, v+e, color="green")    
+        display(gcf())
+    end
+    return aux / mass
+end
 function param_av(fit_masses::Vector{uwreal})
     w = 1 ./ err.(fit_masses).^2
     av = sum(w .* fit_masses) / sum(w)