linalg implemented in juobs

juobs_linalg added, juobs deps extended (LinearAlgebra)

Project.toml

@@ -7,5 +7,6 @@ version = "0.1.0"
 ADerrors = "5e92007d-7bf1-471c-8ceb-4591b8b567a9"
 BDIO = "375f315e-f2c4-11e9-2ef9-134f02f79e27"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

src/juobs.jl

@@ -8,6 +8,7 @@ include("juobs_tools.jl")
 include("juobs_obs.jl")
 
 export read_mesons, read_ms1
+export get_matrix, uwgevp_tot, energies
 export corr_obs, plat_av, lin_fit, x_lin_fit, y_lin_fit
 export meff, dec_const_pcvc 
 

src/juobs_linalg.jl

@@ -11,10 +11,8 @@
 #
 #
 ##################################################################
-## Setting up packages
-using PyPlot, ADerrors, LsqFit, Distributions, Statistics, DelimitedFiles, ForwardDiff, LinearAlgebra
 
-"""
+@doc raw"""
 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 
@@ -24,7 +22,7 @@ at fixed time i=1..T.  It takes as input:
     Each correlator is an vector of uwreal variables of dimension T.
 
 Example:
-
+```@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
 
@@ -51,16 +49,21 @@ 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} )
+function get_matrix(corr_diag::Vector{Vector{uwreal}}, corr_upper::Vector{Vector{uwreal}})
     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.
+    if d != (n*(n-1) / 2)
+        println("Error get_matrix")
+        return nothing
+    end
+    res = Vector{Matrix}(undef, time) # array with all the matrices at each time step.
+    aux = Matrix{uwreal}(undef, n, n)
     for t in 1:time  
-        aux = Matrix{uwreal}(undef, n, n)
         count=0
-        testing = corr_upper[1][1]
+        #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]
@@ -70,35 +73,35 @@ function get_matrix(corr_diag::Vector{Array}, corr_upper::Vector{Array} )
             aux[i,i] = corr_diag[i][t]
         end
         aux[n,n] = corr_diag[n][t]
-        push!(res, aux)
+        res[t] = aux
     end
     return res
 end
-
-"""
+get_matrix(corr_diag::Vector{Corr}, corr_upper::Vector{Corr}) = get_matrix(getfield.(corr_diag, :obs), getfield.(corr_upper, :obs))
+@doc raw"""
 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})
+function energies(evals::Vector{Vector{uwreal}})
     time = length(evals) 
     n = length(evals[1])
-    eff_en = Vector{Array}()
+    eff_en = Vector{Vector}(undef, n)
+    aux_en = Vector{uwreal}(undef, time-1)
     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)))
+            aux_en[t] = log(sqrt(ratio * ratio))
         end
-        uwerr.(aux_en)
-        push!(eff_en, aux_en)
+        #uwerr.(aux_en)
+        eff_en[i] = 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
@@ -112,23 +115,16 @@ The columns of such matrix are the eigenvectors associated with eigenvalues eval
 
 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
+function uwgevp_tot(mat_list::Vector{Array{uwreal, 2}}, tnot::Int64; iter::Int64 = 30, evec::Bool = false)
+    if !evec
+        aux_evals = uwgevp.(mat_list, c_tnot=mat_list[tnot], iter = iter)
+        evals = get_diag.(aux_evals)
         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
+        aux_evals, aux_evec =  uwgevp.(mat_list, c_tnot=mat_list[tnot], iter = iter, evec= true)
+        evals = get_diag.(aux_evals)
+        evecs = aux_evec
+        return evals, evecs
     end
 end
 
@@ -142,14 +138,10 @@ C(t)V = C(t_0)D V, where V is the eigenvector matrix and D is the diagonal matri
 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)
+function uwgevp(c_t:: Matrix{uwreal}; c_tnot::Matrix{uwreal}, iter::Int64 = 30, evec::Bool = 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
+    return uwevp(c, iter = iter, evec = evec)
 end
 
 """
@@ -161,7 +153,7 @@ The keyword iter, set by default to 30, selects the number of iterations of the
 """
 function uwevp(a::Matrix{uwreal}; iter = 30, evec = false)
     n = size(a,1)
-    if evec == true
+    if evec
         u = idty(n)
         for k in 1:iter
             q_k, r_k = qr(a)
@@ -181,10 +173,7 @@ 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
+    res = [a[k, k] for k = 1:n]
     return res
 end
 
@@ -255,7 +244,7 @@ function make_householder(a::Vector{uwreal})
     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 / uwdot(v,v)) * uwdot(v, v)
     #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
@@ -276,7 +265,7 @@ function make_householder(a::Vector{Float64})
     println(typeof(v))
     y = similar(H)
     println(typeof(y))
-    H = H .-  (2.0 / dot(v,v)) * mul!( y, reshape(v, :,1), reshape(v, 1 ,:))
+    H = H .-  (2.0 / dot(v,v)) * mul!(y, v, transpose(v))
     return H
 end
 
@@ -349,12 +338,9 @@ 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 = zeros(n)
     e[k] = 1.0
-    return e
+    return [uwreal(e[i]) for i=1:n]
 end
 
 """
@@ -366,30 +352,19 @@ Depending on the arguments, uwdot returns the vector vector or the matrix matrix
 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
+uwdot(a::Vector{uwreal}, b::Vector{uwreal}) = sum(a .* b)
 
 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
+    c = Matrix{uwreal}(undef, n, p)
+    if m != k
+        println("Error uwdot")
+        return nothing
+    end
+    for i = 1:n
+        for j = 1:p
+            c[i, j] = sum(a[i, :] .* b[:, j])
         end
     end
     return c
@@ -409,9 +384,9 @@ end
 """
 Return the transpose of a matrix of uwreals
 """
-function Base.transpose(a::Matrix{uwreal})
+function transpose(a::Matrix{uwreal})
     res = similar(a)
-    n = size(a, 1 )
+    n = size(a, 1)
     for i in 1:n-1
         for j in i+1:n
             temp = a[i,j]
@@ -429,9 +404,10 @@ uwcopysign(a::uwreal, b::uwreal)
 It implements  the copysign function for uwreal variables
 """
 function uwcopysign(a::uwreal, b::uwreal)
+    c = deepcopy(a)
     aux =  copysign(value(a), value(b))
-    a.mean = aux
-    return a
+    c.mean = aux
+    return c
 end
 
 
@@ -456,7 +432,7 @@ end
 
 
 # operator overloading
-
+#=
 """
 Base.:*(x::uwreal, y::Matrix{uwreal})
 
@@ -475,3 +451,4 @@ function Base.:*(x::uwreal, y::Matrix{uwreal})
     end 
     return res   
 end
+=#
\ No newline at end of file