Merge pull request #9 from iitis/bMPS2

ready for merge

Project.toml

@@ -7,6 +7,7 @@ version = "0.0.1"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 GraphPlot = "a2cc645c-3eea-5389-862e-a155d0052231"
 LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d"

src/SpinGlassPEPS.jl

@@ -6,6 +6,8 @@ module SpinGlassPEPS
     using LightGraphs
     using MetaGraphs
     using CSV
+
+    using DocStringExtensions
     const product = Iterators.product
 
     include("base.jl")

src/base.jl

@@ -1,5 +1,5 @@
-export bond_dimension
-export _verify_bonds
+export bond_dimension, is_left_normalized, is_right_normalized
+export verify_bonds, verify_physical_dims, tensor, rank
 
 for (T, N) in ((:MPO, 4), (:MPS, 3))
     AT = Symbol(:Abstract, T)
@@ -37,15 +37,86 @@ Base.lastindex(a::AbstractMPSorMPO) = lastindex(a.tensors)
 Base.length(a::AbstractMPSorMPO) = length(a.tensors)
 Base.size(a::AbstractMPSorMPO) = (length(a.tensors), )
 
+LinearAlgebra.rank(ψ::MPS) = Tuple(size(A, 2) for A ∈ ψ)
 
-function MPS(vec::Vector{Vector{T}}) where  {T <: Number}
-    L = length(vec)
+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(Θ::AbstractArray{T}, Dcut::Int=typemax(Int), args...) where {T}
+    rank = ndims(Θ)
+    ψ = MPS(T, rank)
+
+    V = reshape(copy(conj(Θ)), (length(Θ), 1))
+
+    for i ∈ 1:rank
+        d = size(Θ, i)
+
+        # reshape
+        @cast M[(x, σ), y] |= V'[x, (σ, y)] (σ:d)
+
+        # decompose
+        U, Σ, V = svd(M, Dcut, args...)
+        V *= Diagonal(Σ)
+
+        # create MPS  
+        @cast A[x, σ, y] |= U[(x, σ), y] (σ:d)
+        ψ[i] = A
+    end
+    ψ
+end
+
+function _left_sweep_SVD(Θ::AbstractArray{T}, Dcut::Int=typemax(Int), args...) where {T}
+    rank = ndims(Θ)
+    ψ = MPS(T, rank)
+
+    U = reshape(copy(Θ), (length(Θ), 1))
+
+    for i ∈ rank:-1:1
+        d = size(Θ, i)
+
+        # reshape
+        @cast M[x, (σ, y)] |= U[(x, σ), y] (σ:d)
+
+        # decompose
+        U, Σ, V = svd(M, Dcut, args...)
+        U *= Diagonal(Σ)
+
+        # create MPS  
+        @cast B[x, σ, y] |= V'[x, (σ, y)] (σ:d)
+        ψ[i] = B
+    end
+    ψ
+end 
+
+function tensor(ψ::MPS, state::Union{Vector, NTuple})
+    C = I
+    for (A, σ) ∈ zip(ψ, state)
+        C *= A[:, idx(σ), :]
+    end
+    tr(C)
+end
+
+function tensor(ψ::MPS)
+    dims = rank(ψ)
+    Θ = Array{eltype(ψ)}(undef, dims)
+
+    for σ ∈ all_states(dims)
+        Θ[idx.(σ)...] = tensor(ψ, σ)
+    end 
+    Θ    
+end
+
+function MPS(states::Vector{Vector{T}}) where {T <: Number}
+    L = length(states)
     ψ = MPS(T, L)
     for i ∈ 1:L
-           A = reshape(vec[i], 1, :, 1)
-        ψ[i] = copy(A)
-    end    
-    return ψ
+        v = states[i]
+        ψ[i] = reshape(copy(v), (1, length(v), 1))
+    end
+    ψ
 end
 
 function MPO(ψ::MPS)
@@ -72,7 +143,7 @@ function MPS(O::MPO)
         @cast A[x, (σ, η), y] := W[x, σ, y, η]
         ψ[i] = A     
     end 
-    return ψ
+    ψ
 end  
 
 function Base.randn(::Type{MPS{T}}, L::Int, D::Int, d::Int) where {T}
@@ -82,20 +153,48 @@ function Base.randn(::Type{MPS{T}}, L::Int, D::Int, d::Int) where {T}
         ψ[i] = randn(T, D, d, D)
     end
     ψ[end] = randn(T, D, d, 1)
-    return ψ
+    ψ
 end
 
 function Base.randn(::Type{MPO{T}}, L::Int, D::Int, d::Int) where {T}
     ψ = randn(MPS{T}, L, D, d^2) 
     MPO(ψ)
+end  
+
+function is_left_normalized(ψ::MPS)
+    for i ∈ 1:length(ψ)
+        A = ψ[i]
+        DD = size(A, 3)
+
+        @tensor Id[x, y] := conj(A[α, σ, x]) * A[α, σ, y] order = (α, σ)
+        I(DD) ≈ Id ? () : return false
+    end  
+    true
+end
+
+function is_right_normalized(ϕ::MPS)   
+    for i ∈ 1:length(ϕ)
+        B = ϕ[i]
+        DD = size(B, 1)
+
+        @tensor Id[x, y] := B[x, σ, α] * conj(B[y, σ, α]) order = (α, σ)
+        I(DD) ≈ Id ? () : return false
+    end 
+    true
 end
 
 function _verify_square(ψ::AbstractMPS)
     arr = [size(A, 2) for A ∈ ψ]
     @assert isqrt.(arr) .^ 2 == arr "Incorrect MPS dimensions"
 end
 
-function _verify_bonds(ψ::AbstractMPS)
+function verify_physical_dims(ψ::AbstractMPS, dims::NTuple)
+    for i ∈ 1:length(ψ)
+        @assert size(ψ[i], 2) == dims[i] "Incorrect physical dim at site $(i)." 
+    end     
+end  
+
+function verify_bonds(ψ::AbstractMPS)
     L = length(ψ)
 
     @assert size(ψ[1], 1) == 1 "Incorrect size on the left boundary." 
@@ -104,19 +203,15 @@ function _verify_bonds(ψ::AbstractMPS)
     for i ∈ 1:L-1
         @assert size(ψ[i], 3) == size(ψ[i+1], 1) "Incorrect link between $i and $(i+1)." 
     end     
-end     
+end  
 
 function Base.show(::IO, ψ::AbstractMPS)
     L = length(ψ)
-    σ_list = [size(ψ[i], 2) for i ∈ 1:L] 
-    χ_list = [size(ψ[i][:, 1, :]) for i ∈ 1:L]
-
-    println("Matrix product state on $L sites:")
-    println("Physical dimensions: ")
-    _show_sizes(σ_list)
+    dims = [size(A) for A in ψ]
+
+    @info "Matrix product state on $L sites:" 
+    _show_sizes(dims)
     println("   ")
-    println("Bond dimensions:   ")
-    _show_sizes(χ_list)
 end
 
 function _show_sizes(dims::Vector, sep::String=" x ", Lcut::Int=8)
@@ -132,4 +227,4 @@ function _show_sizes(dims::Vector, sep::String=" x ", Lcut::Int=8)
         end
         println(dims[end])
     end
-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/compressions.jl

@@ -1,26 +1,4 @@
- export truncate!, canonise!, compress, compress!
-
-function LinearAlgebra.qr(M::AbstractMatrix, Dcut::Int) 
-    fact = pqrfact(M, rank=Dcut)
-    Q = fact[:Q]
-    R = fact[:R]
-    return _qr_fix!(Q, R)
-end     
-
-function rq(M::AbstractMatrix, Dcut::Int) 
-    fact = pqrfact(:c, conj.(M), rank=Dcut)
-    Q = fact[:Q]
-    R = fact[:R]
-    return _qr_fix!(Q, R)'
-end  
-
-function _qr_fix!(Q::AbstractMatrix, R::AbstractMatrix)
-    d = diag(R)
-    ph = d./abs.(d)
-    idim = size(R, 1)
-    q = Matrix(Q)[:, 1:idim]
-    return transpose(ph) .* q
-end
+export truncate!, canonise!, compress
 
 function canonise!(ψ::AbstractMPS)
     canonise!(ψ, :right)
@@ -47,7 +25,7 @@ function _right_sweep_SVD!(ψ::MPS, Dcut::Int=typemax(Int))
         @cast   M̃[(x, σ), y] |= M[x, σ, y]
 
         # decompose
-        U, Σ, V = psvd(M̃, rank=Dcut)
+        U, Σ, V = svd(M̃, Dcut)
 
         # create new    
         d = size(ψ[i], 2)
@@ -68,7 +46,7 @@ function _left_sweep_SVD!(ψ::MPS, Dcut::Int=typemax(Int))
         @cast   M̃[x, (σ, y)] |= M[x, σ, y]
 
         # decompose
-        U, Σ, V = psvd(M̃, rank=Dcut)
+        U, Σ, V = svd(M̃, Dcut)
 
         # create new 
         d = size(ψ[i], 2)
@@ -78,7 +56,7 @@ function _left_sweep_SVD!(ψ::MPS, Dcut::Int=typemax(Int))
 end
 
 
-function compress(ψ::AbstractMPS, Dcut::Int, tol::Number, max_sweeps::Int=4)
+function compress(ψ::AbstractMPS, Dcut::Int, tol::Number=1E-8, max_sweeps::Int=4)
 
     # Initial guess - truncated ψ 
     ϕ = copy(ψ)
@@ -91,30 +69,25 @@ function compress(ψ::AbstractMPS, Dcut::Int, tol::Number, max_sweeps::Int=4)
     overlap = 0 
     overlap_before = 1
 
-    println("Compressing down to: $Dcut") 
+    @info "Compressing down to" Dcut
 
     for sweep ∈ 1:max_sweeps            
         _left_sweep_var!!(ϕ, env, ψ, Dcut)
         overlap = _right_sweep_var!!(ϕ, env, ψ, Dcut)
 
         diff = abs(overlap_before - abs(overlap))
-        println("Convergence: ", diff)
+        @info "Convergence" diff
 
         if diff < tol
-            println("Finished in $sweep sweeps (of $max_sweeps).")
+            @info "Finished in $sweep sweeps of $(max_sweeps)."
             break
         else
             overlap_before = overlap
         end    
     end
-    return ϕ  
+    ϕ  
 end
 
-function compress!(ψ::MPS, Dcut::Int, tol::Number, max_sweeps::Int=4)
-    ϕ = compress(ψ, Dcut, tol, max_sweeps)
-    ψ = copy(ϕ)
-end 
-
 function _left_sweep_var!!(ϕ::MPS, env::Vector{<:AbstractMatrix}, ψ::MPS, Dcut::Int)
     S = eltype(ϕ)
 
@@ -173,6 +146,6 @@ function _right_sweep_var!!(ϕ::MPS, env::Vector{<:AbstractMatrix}, ψ::MPS, Dcu
         @tensor LL[x, y] := conj(A[β, σ, x]) * L[β, α] * B[α, σ, y] order = (α, β, σ)
         env[i+1] = LL
     end
-    return real(env[end][1])
+    real(env[end][1])
 end