Merge branch 'rearrange-linalg'

src/juobs.jl

@@ -9,7 +9,7 @@ include("juobs_tools.jl")
 include("juobs_obs.jl")
 
 export read_mesons, read_ms1, read_ms, read_md, truncate_data!
-export get_matrix, uwgevp_tot, energies, uwdot
+export get_matrix, energies, uwdot, uweigvals, uweigvecs, uweigen, invert, getall_eigvals, getall_eigvecs
 export corr_obs, md_sea, plat_av, lin_fit, x_lin_fit, y_lin_fit, fit_routine
 export meff, dec_const_pcvc, comp_t0
 

src/juobs_linalg.jl

-##################################################################
-# Alessandro Conigli wrote this script.
-# It implements multiple linear algebra utilites optimised to work 
-# with  uwreal data types.
-# The major features so fare implemented are:
-#   - matrix matrix and vector vector multiplication
-#   - norm of uwreal  vectors
-#   - copysign function for uwreal variables 
-#   - creation of identity matrix of uwreal
-#   - element wise multiplication between uwreal variables and matrix of uwreal
-#
-#
-##################################################################
+#=
+using ForwardDiff, LinearAlgebra
+for op in (:eigvals, :eigvecs)
+    @eval function LinearAlgebra.$op(a::Matrix{uwreal})
+
+        function fvec(x::Matrix)
+            return LinearAlgebra.$op(x)
+        end
+        return fvec(value.(a))
+        #return uwreal(LinearAlgebra.$op(a.mean), a.prop, ForwardDiff.derivative($op, a.mean)*a.der)
+    end
+end
+=#
 
 @doc raw"""
 get_matrix(corr_diag::Vector{Array}, corr_upper::Vector{Array} )
@@ -62,7 +62,6 @@ function get_matrix(corr_diag::Vector{Vector{uwreal}}, corr_upper::Vector{Vector
     aux = Matrix{uwreal}(undef, n, n)
     for t in 1:time  
         count=0
-        #testing = corr_upper[1][1]
         for i in range(n-1,1, step=-1)
             for j in range(n, i+1, step=-1)
                 aux[i,j] = corr_upper[d-count][t]
@@ -84,103 +83,161 @@ Given a vector where each entry evals[t] is a uwreal array of eigenvalues, this
 The index t here runs from 1:T=lenght(evals), while the index i stands for the number of energy levels computed: i = length(evals[t])
 It returns a vector array eff_en where each entry  eff_en[t] contains the first N states energies as uwreal objects 
 """
-function energies(evals::Union{Vector{Vector{uwreal}},Array{Array{uwreal}} })
+function energies(evals::Union{Vector{Vector{uwreal}},Array{Array{uwreal}} }; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)     
     time = length(evals) 
     n = length(evals[1])
     eff_en = Array{Array{uwreal}}(undef, n)
     aux_en = Vector{uwreal}(undef, time-1)
     for i in 1:n
-        #aux_en = Vector{uwreal}(undef, time-1)
         for t in 1:time-1
             ratio = evals[t][i] / evals[t+1][i]
             aux_en[t] = 0.5*log(ratio * ratio)
         end
-        uwerr.(aux_en)
+        isnothing(wpm) ? uwerr.(aux_en) : uwerr.(aux_en, wpm)
         eff_en[i] = copy(aux_en)
     end
     return eff_en
 end
 
 @doc raw"""
-uwgevp_tot(mat_list::Vector{Array}, tnot::Int64; iter = 30, evec = false)
-
-Given a list of uwreal matrices and an integer value tnot, this method solves the generalised eigenvalue problem
-for all the matrices in the list, using the same auxiliary matrix mat_list[tnot].
+getall_eigvals(a::Vector{Matrix}, t0; iter=30 )
+
+This function solves a GEVP problem, returning the eigenvalues, for a list of matrices, taking as generalised matrix the one 
+at index t0, i.e:
+        C(t_i)v_i = λ_i C(t_0) v_i, with i=1:lenght(a)
+It takes as input:
+    - a::Vector{Matrix} : a vector of matrices
+    - t0::Int64 : idex value at which the fixed matrix is taken
+    - iter=30 : the number of iterations of the qr algorithm used to extract the eigenvalues 
+
+It returns:
+    - res = Vector{Vector{uwreal}}
+    where res[i] are the generalised eigenvalues of the i-th matrix of the input array. 
+
+Examples:
+'''@example
+mat_array = get_matrix(diag, upper_diag)
+evals = getall_eigvals(mat_array, 5)
+'''
+"""
+function getall_eigvals(a::Vector{Matrix}, t0::Int64; iter=30 )
+    n = length(a)
+    res = Vector{Vector{uwreal}}(undef, n)
+    [res[i] = uweigvals(a[i], a[t0]) for i=1:n]
+    return res
+end
 
-By default, it returns a vector evals where each entry evals[i] is a vector of uwreal containing the 
-generalised eigenvalues  of the matrix mat_list[i].
+@doc raw"""
+uweigvals(a::Matrix{uwreal}; iter = 30)
+uweigvals(a::Matrix{uwreal}, b::Matrix{uwreal};  iter = 30)
 
-If evec is set to true, the methods returns also an array evecs where each entry evecs[i] is a matrix of uwreal. 
-The columns of such matrix are the eigenvectors associated with eigenvalues evals[i].
+This function computes the eigenvalues of a matrix of uwreal objects.
+If a second matrix b is given as input, it returns the generalised eigenvalues instead.
+It takes as input:
+    - a::Matrix{uwreal} : a matrix of uwreal
+    - b::Matrix{uwreal} : a matrix of uwreal, optional
+It returns:
+    - res = Vector{uwreal}: a vector where each elements is an eigenvalue 
 
-The keyword iter, set by default to 30, selects the number of iterations of the qr algorithm before stopping.
 """
-function uwgevp_tot(mat_list::Vector{Matrix}, tnot::Int64; iter::Int64 = 30, evec::Bool = false)
-    n = length(mat_list)
-    evals = Array{Array{uwreal}}(undef, 0 )
-    if !evec
-        println(tnot)
-        c_inv = invert(mat_list[tnot])
-        for i =1:n
-            aux_evals = uwgevp(mat_list[i], c_inv, iter = iter)
-            push!(evals , get_diag(aux_evals))
+function uweigvals(a::Matrix{uwreal}; iter = 30)
+    n = size(a,1)
+    for k in 1:iter
+        q_k, r_k = qr(a)
+        a = uwdot(r_k, q_k)
     end
-        return evals
-    else
-        evecs = Array{Matrix}(undef, 0)
-        for i = 2:n-4
-            c_inv = invert(mat_list[i])
-            aux_evals, aux_evec =  uwgevp(mat_list[i+3], c_inv,  evec=evec, iter=iter)
-            push!(evals,  get_diag(aux_evals))
-            ptilda = norm_evec(aux_evec, mat_list[i])
-            m = uwdot(mat_list[i], ptilda)
-            push!(evecs, m)
-            println( "\n t is ", i+3, " tnot is ", i)            
-            #if  2*tnot < i && i<n-5
-            #    tnot+=5
-            #end
-        end
-        return evals, evecs
+    return get_diag(a)
+end
+function uweigvals(a::Matrix{uwreal}, b::Matrix{uwreal};  iter = 30)
+    c = uwdot(invert(b),a)
+    for k in 1:iter
+        q_k, r_k = qr(c)
+        c = uwdot(r_k, q_k)
     end
+    return get_diag(c)
 end
-"""
-uwgevp(c_t:: Matrix{uwreal}, c_tnot_inv::Matrix{uwreal}; iter = 30, evec = false)
 
-Given two matrices of uwreal C(t) and C(t0), this method solves the generalized eigenvalue problem
-C(t)V = C(t_0)D V, where V is the eigenvector matrix and D is the diagonal matrix of eigenvalues.
+@doc raw"""
+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:
+'''@example
+mat_array = get_matrix(diag, upper_diag)
+evecs = getall_eigvecs(mat_array, 5)
+'''
+"""
+function getall_eigvecs(a::Vector{Matrix}, delta_t; iter=30)
+    n = length(a)-delta_t
+    res = Vector{Matrix{uwreal}}(undef, n)
+    [res[i] = uweigvecs(a[i+delta_t], a[i]) for i=1:n]
+    return res 
 
-If evec = true also the eigenvectors are returned as columns of the orthogonal matrix u.
-The keyword iter, set by default to 30, selects the number of iterations of the qr algorithm before stopping.
-"""
-function uwgevp(c_t:: Matrix{uwreal}, c_tnot_inv::Matrix{uwreal}; iter::Int64 = 30, evec::Bool = false)
-    c = uwdot(c_tnot_inv, c_t) #matrix to diagonalize
-    return uwevp(c, iter = iter, evec = evec)
 end
 
-"""
-uwevp(a::Matrix{uwreal}; iter = 30, evec = false)
+@doc raw"""
+uweigvecs(a::Matrix{uwreal}; iter = 30)
+uweigvecs(a::Matrix{uwreal}, b::Matrix{uwreal};  iter = 30)
+
+This function computes the eigenvectors of a matrix of uwreal objects.
+If a second matrix b is given as input, it returns the generalised eigenvectors instead.
+It takes as input:
+    - a::Matrix{uwreal} : a matrix of uwreal
+    - b::Matrix{uwreal} : a matrix of uwreal, optional
+It returns:
+    - res = Matrix{uwreal}: a matrix where each column is an eigenvector 
 
-Return the eigenvalues of the uwreal matrix a.
-If evec = true also the eigenvectors are returned as columns of the orthogonal matrix u.
-The keyword iter, set by default to 30, selects the number of iterations of the qr algorithm before stopping.
 """
-function uwevp(a::Matrix{uwreal}; iter = 30, evec = false)
+function uweigvecs(a::Matrix{uwreal}; iter = 30)
     n = size(a,1)
-    if evec
     u = idty(n)
     for k in 1:iter
         q_k, r_k = qr(a)
         a = uwdot(r_k, q_k)
         u = uwdot(u, q_k)
     end
-        return a, u
-    else
+    return  u
+end
+function uweigvecs(a::Matrix{uwreal}, b::Matrix{uwreal}; iter = 30)
+    c  = uwdot(invert(b), a)
+    n = size(c,1)
+    u = idty(n)
     for k in 1:iter
-            q_k, r_k = qr(a)
-            a = uwdot(r_k, q_k)
-        end
-        return a
+        q_k, r_k = qr(c)
+        c = uwdot(r_k, q_k)
+        u = uwdot(u, q_k)
     end
+    return  u
+end
+
+@doc raw"""
+uweigen(a::Matrix{uwreal}; iter = 30)
+uweigen(a::Matrix{uwreal}, b::Matrix{uwreal};  iter = 30)
+
+This function computes the eigenvalues and the eigenvectors  of a matrix of uwreal objects.
+If a second matrix b is given as input, it returns the generalised eigenvalues and eigenvectors instead.
+It takes as input:
+    - a::Matrix{uwreal} : a matrix of uwreal
+    - b::Matrix{uwreal} : a matrix of uwreal, optional
+It returns:
+    - evals = Vector{uwreal}: a vector where each elements is an eigenvalue 
+    -evecs  = Matrix{uwreal}: a matrix where the i-th column is the eigenvector of the i-th eigenvalue
+"""
+function uweigen(a::Matrix{uwreal}; iter = 30)
+    return uweigvals(a, iter=iter), uweigvecs(a, iter=iter)
+end
+function uweigen(a::Matrix{uwreal}, b::Matrix{uwreal}; iter = 30)
+    return uweigvals(a, b, iter=iter), uweigvecs(a, b, iter=iter)
 end
 
 function norm_evec(evec::Matrix{uwreal}, ctnot::Matrix{uwreal})
@@ -190,9 +247,6 @@ function norm_evec(evec::Matrix{uwreal}, ctnot::Matrix{uwreal})
         aux_norm = (uwdot(evec[:,i], uwdot(ctnot, evec[:,i])).^2).^0.25
         res_evec[:,i] = 1 ./aux_norm .* evec[:,i]
     end
-    #println("res_evec is ", res_evec)
-    #println("evec is ", evec)
-    #println(res_evec)
     return res_evec
 end
 function get_diag(a::Matrix{uwreal})
@@ -200,6 +254,14 @@ function get_diag(a::Matrix{uwreal})
     res = [a[k, k] for k = 1:n]
     return res
 end
+function get_diag(v::Vector{uwreal})
+    n = length(v)
+    res = idty(n)
+    for i =1:n
+        res[i,i] = v[i]
+    end
+    return res
+end
 
 
 """
@@ -269,12 +331,10 @@ function make_householder(a::Vector{uwreal})
     v[1] = uwreal(1.0) 
     H = idty(n)
     H = H -  (2.0 / uwdot(v,v)) * uwdot(reshape(v, :,1), transpose(v)) #this last term is a vector product  
-    #H = H -  (2.0 / dot(value.(v),value.(v))) * uwdot(reshape(v, :,1), reshape(v, 1 ,:))
     for i in 1:n
         for j in 1:n
-            #println(H[i,j].ids)
             if H[i,j]!= 0 &&  H[i,j]!=1 
-                #uwerr(H[i,j])
+
             end
         end
     end
@@ -347,7 +407,6 @@ function backward_sub(u::Matrix{uwreal}, y::Vector{uwreal})
         end
         x[i] = (y[i] - temp)  /u[i,i]
         if x[i] !=0
-            #uwerr(x[i])
         end
     end
     
@@ -402,9 +461,9 @@ uwnorm(a::Vector{uwreal})
 
 It returns the norm of a vector of uwreal
 """
-function uwnorm(a::Vector{uwreal})
+function uwnorm(a::Vector{uwreal}; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
     norm = sqrt(uwdot(a,a))
-    uwerr(norm)
+    isnothing(wpm) ? uwerr(norm) : uwerr(norm, wpm)
     return norm 
 end
 
@@ -476,7 +535,6 @@ function Base.:*(x::uwreal, y::Matrix{uwreal})
         for j in 1:m
             res[i,j] = x * y[i,j]
             if res[i,j] !=0 && res[i,j] !=1
-                #uwerr(res[i,j])
             end
         end        
     end 

src/juobs_reader.jl

@@ -73,7 +73,7 @@ end
 read_mesons(path::String, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{Int64, Nothing}=nothing)
 read_mesons(path::Vector{String}, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{Int64, Nothing}=nothing)
 
-This faction read a mesons dat file at a given path and returns a vector of CData structures for different masses and Dirac structures.
+This function read a mesons dat file at a given path and returns a vector of CData structures for different masses and Dirac structures.
 Dirac structures g1 and/or g2 can be passed as string arguments in order to filter correaltors.
 ADerrors id can be specified as argument. If is not specified, the id is fixed according to the ensemble name (example: "H400"-> id = 400)
 

src/juobs_tools.jl

@@ -247,6 +247,14 @@ fit_routine(model::Function, xdata::Array{uwreal}, ydata::Array{uwreal}, param::
 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 Γ(t) is summed up a point where its error δΓ(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.
+
 '''@example
 @. 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]