including gevp code and linear algebra routines

src/uwLinAlg.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
+#
+#
+##################################################################
+## Setting up packages
+using PyPlot, ADerrors, LsqFit, Distributions, Statistics, DelimitedFiles, ForwardDiff, LinearAlgebra
+
+"""
+get_matrix(corr_diag::Vector{Array}, corr_upper::Vector{Array} )
+
+This method  returns an array of dim T where each element is a  symmetrix matrix  of dimension n of uwreal correlators 
+at fixed time i=1..T.  It takes as input:
+    corr_diag  :  vector of dimension n of correlators liying on the diagonal 
+    corr_upper :  vector of correlators liying on the upper diagonal.
+    Each correlator is an vector of uwreal variables of dimension T.
+
+Example:
+
+for i in 1:n
+a[i,i]  = Vector{uwreal} # vector of uwreal variables of dimension T. They will constitute the diagonal elements of the matrices
+
+for i in 1:n-1 
+    for j in i+1:n
+        a[i,j]  = Vector{uwreal} # vector of uwreal variables of dimension T. They will constitute the upper diagonal elements of the matrices. A matrix
+        of dimension n*n has n(n-1)/2 upper diagonal elements.
+
+Assume n=4
+
+diagonal = Vector{Array}()
+push!(diagonal, a[1,1],a[2,2],a[3,3],a[4,4])
+upsize = Vector{Array}()
+push!(upsize, a[1,2], a[1,3], a[1,4], a[2,3], a[2,4], a[3,4])
+
+array_of_matrices = get_matrix(diagonal, upsize)
+
+Julia> T-element Array{Array,1}
+
+size(array_of_matrices)
+
+Julia> (T,)
+
+array_of_matrices[t] # t in 1:T
+
+Julia> 4*4 Array{uwreal,2}
+"""
+function get_matrix(corr_diag::Vector{Array}, corr_upper::Vector{Array} )
+    time = length(corr_diag[1]) # total time slices 
+    n = length(corr_diag) # matrix dimension
+    d = length(corr_upper) # upper elements
+    res = Vector{Array}() # array with all the matrices at each time step.
+    for t in 1:time  
+        aux = Matrix{uwreal}(undef, n, n)
+        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]
+                aux[j,i] = corr_upper[d-count][t]
+                count +=1
+            end
+            aux[i,i] = corr_diag[i][t]
+        end
+        aux[n,n] = corr_diag[n][t]
+        push!(res, aux)
+    end
+    return res
+end
+
+"""
+energies(evals::Vector{Array})
+
+Given a vector where each entry evals[t] is a uwreal array of eigenvalues, this method computes the effective energies of the first N states, where N=dim(evals[t]).
+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::Vector{Array})
+    time = length(evals) 
+    n = length(evals[1])
+    eff_en = Vector{Array}()
+    for i in 1:n
+        aux_en = Vector{uwreal}()
+        for t in 1:time-1
+            ratio = evals[t][i] / evals[t+1][i]
+            push!(aux_en, log(ratio*uwcopysign(uwreal(1.0),ratio)))
+        end
+        uwerr.(aux_en)
+        push!(eff_en, aux_en)
+    end
+    return eff_en
+end
+
+"""
+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].
+
+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].
+
+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].
+
+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{Array}, tnot::Int64; iter = 30, evec = false)
+    n = length(mat_list)
+    evals = Vector{Array}()
+    if evec == false
+        for i in 1:n
+           aux_evals =  uwgevp(mat_list[i], mat_list[tnot], iter = iter)
+           push!(evals, get_diag(aux_evals))
+        end
+        return evals
+    else
+        evecs = Vector{Array}()
+        for i in 1:n
+            aux_evals, aux_evec=  uwgevp(mat_list[i], mat_list[tnot], iter = iter, evec = true)
+            push!(evals, get_diag(aux_evals))
+            push!(evecs, aux_evec)
+         end
+         return evals, evecs
+    end
+end
+
+
+"""
+uwgevp(c_t:: Matrix{uwreal}, c_tnot::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.
+
+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::Matrix{uwreal}; iter = 30, evec = false)
+    c_tnot_inv = invert(c_tnot)
+    c = uwdot(c_tnot_inv, c_t) #matrix to diagonalize
+    if evec == false
+        return uwevp(c, iter = iter)
+    else
+        return uwevp(c, iter = iter, evec = true)
+    end
+end
+
+"""
+uwevp(a::Matrix{uwreal}; iter = 30, evec = false)
+
+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)
+    n = size(a,1)
+    if evec == true
+        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
+        for k in 1:iter
+            q_k, r_k = qr(a)
+            a = uwdot(r_k, q_k)
+        end
+        return a
+    end
+end
+
+
+function get_diag(a::Matrix{uwreal})
+    n = size(a,1)
+    res = Vector{uwreal}()
+    for i in 1:n
+        push!(res, a[i,i])
+    end
+    return res
+end
+
+
+"""
+Given an input matrix of uwreals, this function returns the inverse of that matrix using QR decomposition and
+backward substitution.
+"""
+function invert(a::Matrix{uwreal})
+    #q,r = qr(hess_red(a))
+    q,r = qr(a)
+    r_inv = upper_inversion(r)
+    return uwdot(r_inv, transpose(q))
+end
+
+
+"""
+qr(A::Matrix{uwreal})
+
+qr(A::Matrix{Float64})
+
+This methods returns the QR decomposition of a matrix A of either uwreal or Float64 variables
+"""
+function qr(A::Matrix{uwreal})
+    m,n = size(A)
+    Q = idty(m)
+    for i in 1:(n-(m==n))
+        H = idty(m)
+        H[i:m,i:m] = make_householder(A[i:m,i])
+        Q = uwdot(Q,H)
+        A = uwdot(H,A)
+    end
+    for i in 2:n
+        for j in 1:i-1
+            A[i,j] = 0.0
+        end
+    end
+    return Q,A
+end
+function qr(A::Matrix{Float64})
+    m,n = size(A)
+    Q = Matrix{Float64}(I,m,m)
+    for i in 1:(n-(m==n))
+        H =  Matrix{Float64}(I,m,m)
+        H[i:m,i:m] = make_householder(A[i:m,i])
+        X = similar(Q)
+        Y = similar(A)
+        Q = mul!(X,Q,H)
+        A = mul!(Y,H,A)
+    end
+    return Q,A
+end
+
+"""
+make_householder(a::Vector{uwreal})
+
+make_householder(a::Vector{Float64})
+
+This method returns the householder matrix used in qr decomposition to get an upper triangular matrix
+
+It takes as input either a vector of uwreal or a vector of Float64
+"""
+function make_householder(a::Vector{uwreal})
+    n = length(a)
+    v = Vector{uwreal}(undef, n)
+    for i in 1:n
+        v[i] = a[i]  / (a[1] + uwcopysign(uwnorm(a), a[1]))
+    end
+    v[1] = uwreal(1.0) 
+    H = idty(n)
+    H = H -  (2.0 / uwdot(v,v)) * uwdot(reshape(v, :,1), reshape(v, 1 ,:))
+    #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
+    return H
+end
+
+function make_householder(a::Vector{Float64})
+    n = length(a)
+    v = a / (a[1] + copysign(norm(a), a[1])) 
+    v[1] = 1.0
+    H = Matrix{Float64}(I, n, n)
+    println(typeof(v))
+    y = similar(H)
+    println(typeof(y))
+    H = H .-  (2.0 / dot(v,v)) * mul!( y, reshape(v, :,1), reshape(v, 1 ,:))
+    return H
+end
+
+"""
+hess_red(a::Matrix{uwreal})
+
+Given a matrix of uwreal variables, this function returns the Hessenberg form of the matrix.
+"""
+function hess_red(a::Matrix{uwreal})
+    n = size(a, 1)
+    res = Matrix{uwreal}(undef, n, n)
+    res = a
+    for j in 1:n-2
+        mi = idty(n)
+        mi_inv = idty(n)
+        for k in j+2:n
+            mi[k,j+1] = res[k,j] / res[j+1, j]
+            mi_inv[k,j+1] = -mi[k,j+1]
+        end
+        res = uwdot(mi_inv, uwdot(res, mi) )
+    end
+    return res
+end
+
+"""
+upper_inversion(u::Matrix{uwreal})
+
+This method uses backward propagation to compute the inverse of an upper triangular matrix u. Note that the inverse of u must be upper triangular as well.
+Tests executed so far always confirmed this property.
+"""
+function upper_inversion(u::Matrix{uwreal})
+    n = size(u,1)
+    u_inv = Matrix{uwreal}(undef, n, n)
+    for i in 1:n
+        e_i = canonical_basis(n, i)
+        x = backward_sub(u, e_i)
+        u_inv[:,i] = x
+    end
+    return u_inv
+end
+"""
+backward_sub(u::Matrix{uwreal}, y::Vector{uwreal})
+
+This method operates backward substitution to solve a liner system Ux = y where U is an upper triangular uwreal matrix 
+and y is a uwreal vector. It returns a vector of uwreal x that might be used to inverte the upper triangular matrix U.
+"""
+function backward_sub(u::Matrix{uwreal}, y::Vector{uwreal})
+    n = length(y)
+    x =  Vector{uwreal}(undef, n)
+    x[n] = y[n] / u[n,n]
+    for i in range(n-1,1, step=-1)
+        temp = 0.0
+        for j in i+1:n
+            temp += u[i,j] * x[j]
+        end
+        x[i] = (y[i] - temp)  /u[i,i]
+        if x[i] !=0
+            #uwerr(x[i])
+        end
+    end
+    
+    return x
+end
+
+"""
+canonical_basis(n::Int64, k::Int64)
+
+Given the dimension n and the position k, canonical_basis returns a vector of uwreal with all entries set to zero
+except for the k-th entry which is set to 1.0.
+This method is used as input vector in backward_sub to compute the inverse of an upper triangular matrix.
+"""
+function canonical_basis(n::Int64, k::Int64)
+    e = Vector{uwreal}(undef,n)
+    for i in 1:n
+        e[i] = 0.0
+    end
+    e[k] = 1.0
+    return e
+end
+
+"""
+uwdot(a::Vector{uwreal}, b::Vector{uwreal})
+
+uwdot(a::Matrix{uwreal}, b::Matrix{uwreal})
+
+Depending on the arguments, uwdot returns the vector vector or the matrix matrix product 
+of uwreal data type.
+
+"""
+function uwdot(a::Vector{uwreal}, b::Vector{uwreal})
+    n = length(a)  
+    c = uwreal(0) #Vector{uwreal}(undef, n)
+    for i in 1:n
+       c += a[i] * b[i] 
+    end
+    #uwerr(c)
+    return c
+end
+
+function uwdot(a::Matrix{uwreal}, b::Matrix{uwreal})
+    n,m = size(a)
+    k,p = size(b)
+    c = Matrix{uwreal}(undef, n,p)
+    for i in 1:n
+        for j in 1:p
+            temp = 0.0 # uwreal()
+            for l in 1:m
+                temp += a[i,l] * b[l,j]
+            end
+            c[i,j] = temp
+            if c[i,j] !=0 && c[i,j]!=1
+                #uwerr(c[i,j])
+            end
+        end
+    end
+    return c
+end
+    
+"""
+uwnorm(a::Vector{uwreal})
+
+It returns the norm of a vector of uwreal
+"""
+function uwnorm(a::Vector{uwreal})
+    norm = sqrt(uwdot(a,a))
+    uwerr(norm)
+    return norm
+end
+
+"""
+Return the transpose of a matrix of uwreals
+"""
+function Base.transpose(a::Matrix{uwreal})
+    res = similar(a)
+    n = size(a, 1 )
+    for i in 1:n-1
+        for j in i+1:n
+            temp = a[i,j]
+            res[i,j] = a[j,i]
+            res[j,i] = temp
+        end
+        res[i,i] = a[i,i]
+    end
+    res[n,n] = a[n,n]
+    return res
+end
+"""
+uwcopysign(a::uwreal, b::uwreal)
+
+It implements  the copysign function for uwreal variables
+"""
+function uwcopysign(a::uwreal, b::uwreal)
+    aux =  copysign(value(a), value(b))
+    a.mean = aux
+    return a
+end
+
+
+"""
+idty(n::Int64)
+
+It returns an indetity matrix of uwreal variables given the size n 
+"""
+function idty(n::Int64)
+    res = Matrix{uwreal}(undef,n,n)
+    for i in 1:n
+        for j in 1:n
+            if i==j
+                res[i,j]=1.0
+            else
+                res[i,j]=0.0
+            end
+        end
+    end
+    return res
+end
+
+
+# operator overloading
+
+"""
+Base.:*(x::uwreal, y::Matrix{uwreal})
+
+Multiplication operator overloading between uwreal variable and matrix of uwreal 
+"""
+function Base.:*(x::uwreal, y::Matrix{uwreal})
+    n,m = size(y)
+    res = Matrix{uwreal}(undef, n,m)
+    for i in 1:n
+        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 
+    return res   
+end
\ No newline at end of file