Merge branch 'bMPS2' of https://github.com/iitis/SpinGlassPEPS.jl int…

…o bMPS2

src/base.jl

@@ -1,6 +1,5 @@
 export bond_dimension, is_left_normalized, is_right_normalized
 export verify_bonds, verify_physical_dims, tensor, rank
-using Base
 
 for (T, N) in ((:MPO, 4), (:MPS, 3))
     AT = Symbol(:Abstract, T)
@@ -40,17 +39,17 @@ Base.size(a::AbstractMPSorMPO) = (length(a.tensors), )
 
 LinearAlgebra.rank(ψ::MPS) = Tuple(size(A, 2) for A ∈ ψ)
 
-MPS(A::Array) = MPS(A, :right)
-MPS(A::Array, s::Symbol, args...) = MPS(A, Val(s), typemax(Int), args...)
-MPS(A::Array, s::Symbol, Dcut::Int, args...) = MPS(A, Val(s), Dcut, args...)
-MPS(A::Array, ::Val{:right}, Dcut::Int, args...) = _left_sweep_SVD(A, Dcut, args...)
-MPS(A::Array, ::Val{:left}, Dcut::Int, args...) = _right_sweep_SVD(A, Dcut, args...)
+MPS(A::AbstractArray) = MPS(A, :right)
+MPS(A::AbstractArray, s::Symbol, args...) = MPS(A, Val(s), typemax(Int), args...)
+MPS(A::AbstractArray, s::Symbol, Dcut::Int, args...) = MPS(A, Val(s), Dcut, args...)
+MPS(A::AbstractArray, ::Val{:right}, Dcut::Int, args...) = _left_sweep_SVD(A, Dcut, args...)
+MPS(A::AbstractArray, ::Val{:left}, Dcut::Int, args...) = _right_sweep_SVD(A, Dcut, args...)
 
-function _right_sweep_SVD(Θ::Array{T}, Dcut::Int=typemax(Int), args...) where {T}
+function _right_sweep_SVD(Θ::AbstractArray{T}, Dcut::Int=typemax(Int), args...) where {T}
     rank = ndims(Θ)
     ψ = MPS(T, rank)
 
-    V = Reshape(copy(conj(Θ)), (length(Θ), 1))
+    V = reshape(copy(conj(Θ)), (length(Θ), 1))
 
     for i ∈ 1:rank
         d = size(Θ, i)
@@ -69,11 +68,11 @@ function _right_sweep_SVD(Θ::Array{T}, Dcut::Int=typemax(Int), args...) where {
     ψ
 end
 
-function _left_sweep_SVD(Θ::Array{T}, Dcut::Int=typemax(Int), args...) where {T}
+function _left_sweep_SVD(Θ::AbstractArray{T}, Dcut::Int=typemax(Int), args...) where {T}
     rank = ndims(Θ)
     ψ = MPS(T, rank)
 
-    U = Reshape(copy(Θ), (length(Θ), 1))
+    U = reshape(copy(Θ), (length(Θ), 1))
 
     for i ∈ rank:-1:1
         d = size(Θ, i)
@@ -93,8 +92,11 @@ function _left_sweep_SVD(Θ::Array{T}, Dcut::Int=typemax(Int), args...) where {T
 end 
 
 function tensor(ψ::MPS, state::Union{Vector, NTuple})
-    ρ = [A[:, idx(σ), :] for (A, σ) ∈ zip(ψ, state)]
-    tr(*(ρ...))
+    C = I
+    for  (A, σ) ∈ zip(ψ, state)
+        C *= A[:, idx(σ), :]
+    end
+    tr(C)
 end
 
 function tensor(ψ::MPS)
@@ -112,7 +114,7 @@ function MPS(states::Vector{Vector{T}}) where {T <: Number}
     ψ = MPS(T, L)
     for i ∈ 1:L
         v = states[i]
-        ψ[i] = Reshape(copy(v), (1, length(v), 1))
+        ψ[i] = reshape(copy(v), (1, length(v), 1))
     end
     ψ
 end

src/compression.jl

@@ -1,5 +1,5 @@
 
-function make_left_canonical(t1::Array{T, 3}, t2::Array{T, 3}) where T <: AbstractFloat
+function make_left_canonical(t1::AbstractArray{T, 3}, t2::AbstractArray{T, 3}) where T <: AbstractFloat
     s = size(t1)
 
     p1 = [1,3,2]
@@ -21,7 +21,7 @@ end
 
 
 
-function make_right_canonical(t1::Array{T, 3}, t2::Array{T, 3}) where T <: AbstractFloat
+function make_right_canonical(t1::AbstractArray{T, 3}, t2::AbstractArray{T, 3}) where T <: AbstractFloat
 
     s = size(t2)
 
@@ -91,7 +91,7 @@ function right_canonical_approx(mps::Vector{Array{T, 3}}, χ::Int) where T <: Ab
 end
 
 
-function QR_make_right_canonical(t2::Array{T, 3}) where T <: AbstractFloat
+function QR_make_right_canonical(t2::AbstractArray{T, 3}) where T <: AbstractFloat
 
     s = size(t2)
     p = [2,3,1]
@@ -109,7 +109,7 @@ function QR_make_right_canonical(t2::Array{T, 3}) where T <: AbstractFloat
     Bnew, R
 end
 
-function QR_make_left_canonical(t2::Array{T, 3}) where T <: AbstractFloat
+function QR_make_left_canonical(t2::AbstractArray{T, 3}) where T <: AbstractFloat
 
     s = size(t2)
     p = [1,3,2]
@@ -127,12 +127,12 @@ function QR_make_left_canonical(t2::Array{T, 3}) where T <: AbstractFloat
 end
 
 """
-    function R_update(U::Array{T, 3}, U_exact::Array{T, 3}, R::Array{T, 2}) where T <: AbstractFloat
+    function R_update(U::AbstractArray{T, 3}, U_exact::AbstractArray{T, 3}, R::AbstractArray{T, 2}) where T <: AbstractFloat
 
 update the right enviroment in the approximation,
 Return matrix, the updated enviroment
 """
-function R_update(U::Array{T, 3}, U_exact::Array{T, 3}, R::Array{T, 2}) where T <: AbstractFloat
+function R_update(U::AbstractArray{T, 3}, U_exact::AbstractArray{T, 3}, R::AbstractArray{T, 2}) where T <: AbstractFloat
     C = zeros(T, size(U, 1), size(U_exact, 1))
     @tensor begin
         #v concers contracting modes of size 1 in C
@@ -142,12 +142,12 @@ function R_update(U::Array{T, 3}, U_exact::Array{T, 3}, R::Array{T, 2}) where T
 end
 
 """
-    L_update(U::Array{T, 3}, U_exact::Array{T, 3}, R::Array{T, 2}) where T <: AbstractFloat
+    L_update(U::AbstractArray{T, 3}, U_exact::AbstractArray{T, 3}, R::AbstractArray{T, 2}) where T <: AbstractFloat
 
 update the left enviroment in the approximation,
 Return matrix, the updated enviroment
 """
-function L_update(U::Array{T, 3}, U_exact::Array{T, 3}, R::Array{T, 2}) where T <: AbstractFloat
+function L_update(U::AbstractArray{T, 3}, U_exact::AbstractArray{T, 3}, R::AbstractArray{T, 2}) where T <: AbstractFloat
     C = zeros(T, size(U, 2), size(U_exact, 2))
     @tensor begin
         C[a,b] = U[x,a,v]*U_exact[y,b,v]*R[x,y]

src/contractions.jl

@@ -10,7 +10,7 @@ export left_env, right_env, dot!
 #
 
 function LinearAlgebra.dot(ψ::MPS, state::Union{Vector, NTuple}) 
-    C = ones(eltype(ψ), 1, 1)
+    C = I
 
     for (M, σ) ∈ zip(ψ, state)        
         i = idx(σ)
@@ -28,7 +28,7 @@ function LinearAlgebra.dot(ϕ::MPS, ψ::MPS)
         M̃ = conj(ϕ[i])
         @tensor C[x, y] := M̃[β, σ, x] * C[β, α] * M[α, σ, y] order = (α, β, σ) 
     end
-    C[1]
+    tr(C)
 end
 
 function left_env(ϕ::MPS, ψ::MPS) 
@@ -80,7 +80,7 @@ function LinearAlgebra.dot(ϕ::MPS, O::Union{Vector, NTuple}, ψ::MPS) #where T
         Mat = O[i]
         @tensor C[x, y] := M̃[β, σ, x] * Mat[σ, η] * C[β, α] * M[α, η, y] order = (α, η, β, σ)
     end
-    C[1]
+    tr(C)
 end
 
 function LinearAlgebra.dot(O::AbstractMPO, ψ::T) where {T <: AbstractMPS}

src/search.jl

@@ -75,9 +75,9 @@ function _apply_exponent!(ψ::AbstractMPS, ig::MetaGraph, dβ::Number, i::Int, j
     J = get_prop(ig, i, j, :J) 
     C = [exp(0.5 * dβ * k * J * l) for k ∈ [-1, 1], l ∈ [-1, 1]]
 
-    D = I(2)
-    P = [[1.] [0.]]
-    if j == length(ψ) δ = P' else δ = D end
+
+    δ = j == length(ψ) ? P' : D
+    j == length(ψ) ? δ = P' : δ = D
 
     @cast M̃[(x, a), σ, (y, b)] := C[σ, x] * δ[x, y] * M[a, σ, b]                      
     ψ[j] = M̃
@@ -86,19 +86,15 @@ end
 function _apply_projector!(ψ::AbstractMPS, i::Int)
     M = ψ[i]
 
-    D = I(2)
-    P = [[1.] [0.]]
-    if i == 1 δ = P else δ = D end
+    δ = i == 1 ? P : D
 
     @cast M̃[(x, a), σ, (y, b)] := D[σ, y] * δ[x, y] * M[a, σ, b]
     ψ[i] = M̃
 end
 
 function _apply_nothing!(ψ::AbstractMPS, i::Int) 
     M = ψ[i] 
-
-    D = I(2)
-    P = [[1.] [0.]]     
+
     if i == 1  δ = P  elseif  i == length(ψ)  δ = P'  else  δ = D end  
 
     @cast M̃[(x, a), σ, (y, b)] := δ[x, y] * M[a, σ, b] 

src/utils.jl

@@ -1,6 +1,6 @@
 export idx, ising, proj
 export HadamardMPS, rq
-export Reshape, all_states
+export all_states
 
 ising(σ::Union{Vector, NTuple}) = 2 .* σ .- 1
 
@@ -24,13 +24,6 @@ end
 
 HadamardMPS(L::Int) = MPS(fill([1., 1.] / sqrt(2), L))
 
-function Reshape(A::AbstractArray{T}, dims::Tuple) where {T <: Number}
-    ord = reverse(1:length(dims))
-
-    A = reshape(A, reverse(dims))
-    permutedims(A, ord)
-end 
-
 function LinearAlgebra.qr(M::AbstractMatrix, Dcut::Int, args...) 
     fact = pqrfact(M, rank=Dcut, args...)
     Q = fact[:Q]

test/base.jl

@@ -43,7 +43,7 @@ end
 @testset "Reshaping (row-wise)" begin
     vec = Vector(1:6)
 
-    A = Reshape(vec, (2, 3))
+    A = reshape_row(vec, (2, 3))
     B = [1 2 3; 4 5 6]
 
     @test A ≈ B
@@ -55,7 +55,7 @@ end
     vec = kron(states...)
 
     ψ = tensor(MPS(states))
-    ϕ = Reshape(vec, dims)
+    ϕ = reshape_row(vec, dims)
 
     @test ψ ≈ ϕ
 end 

test/runtests.jl

@@ -5,6 +5,13 @@ using TensorOperations
 
 using Test
 
+function reshape_row(A::AbstractArray{T}, dims::Tuple) where {T <: Number}
+    ord = reverse(1:length(dims))
+
+    A = reshape(A, reverse(dims))
+    permutedims(A, ord)
+end 
+
 my_tests = []
 if CUDA.functional() && CUDA.has_cutensor() && false
     push!(my_tests,