diff --git a/Project.toml b/Project.toml index 571edf55..c0d3851a 100644 --- a/Project.toml +++ b/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" diff --git a/src/SpinGlassPEPS.jl b/src/SpinGlassPEPS.jl index 5473ea85..fe841cac 100644 --- a/src/SpinGlassPEPS.jl +++ b/src/SpinGlassPEPS.jl @@ -6,6 +6,8 @@ module SpinGlassPEPS using LightGraphs using MetaGraphs using CSV + + using DocStringExtensions const product = Iterators.product include("base.jl") diff --git a/src/base.jl b/src/base.jl index 9f24a624..703b4143 100644 --- a/src/base.jl +++ b/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,12 +153,34 @@ 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) @@ -95,7 +188,13 @@ function _verify_square(ψ::AbstractMPS) @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 \ No newline at end of file diff --git a/src/compression.jl b/src/compression.jl index dc937dee..78fc9303 100644 --- a/src/compression.jl +++ b/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] diff --git a/src/compressions.jl b/src/compressions.jl index 8f572718..bd358340 100644 --- a/src/compressions.jl +++ b/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 diff --git a/src/contractions.jl b/src/contractions.jl index 0faa3171..bdb0dd6c 100644 --- a/src/contractions.jl +++ b/src/contractions.jl @@ -1,4 +1,3 @@ -using Base export left_env, right_env, dot! # --------------------------- Conventions ------------------------ @@ -10,6 +9,16 @@ export left_env, right_env, dot! # --------------------------------------------------------------- # +function LinearAlgebra.dot(ψ::MPS, state::Union{Vector, NTuple}) + C = I + + for (M, σ) ∈ zip(ψ, state) + i = idx(σ) + C = M[:, i, :]' * (C * M[:, i, :]) + end + tr(C) +end + function LinearAlgebra.dot(ϕ::MPS, ψ::MPS) T = promote_type(eltype(ψ), eltype(ϕ)) C = ones(T, 1, 1) @@ -19,7 +28,7 @@ function LinearAlgebra.dot(ϕ::MPS, ψ::MPS) M̃ = conj(ϕ[i]) @tensor C[x, y] := M̃[β, σ, x] * C[β, α] * M[α, σ, y] order = (α, β, σ) end - return C[1] + tr(C) end function left_env(ϕ::MPS, ψ::MPS) @@ -37,7 +46,7 @@ function left_env(ϕ::MPS, ψ::MPS) @tensor C[x, y] := M̃[β, σ, x] * C[β, α] * M[α, σ, y] order = (α, β, σ) L[i+1] = C end - return L + L end # NOT tested yet @@ -56,12 +65,31 @@ function right_env(ϕ::MPS, ψ::MPS) @tensor D[x, y] := M[x, σ, α] * D[α, β] * M̃[y, σ, β] order = (β, α, σ) R[i] = D end - return R + R end +""" +$(TYPEDSIGNATURES) + +Calculates the norm of an MPS \$\\ket{\\phi}\$ +""" LinearAlgebra.norm(ψ::AbstractMPS) = sqrt(abs(dot(ψ, ψ))) -function LinearAlgebra.dot(ϕ::MPS, O::Vector{T}, ψ::MPS) where {T <: AbstractMatrix} + +""" +$(TYPEDSIGNATURES) + +Calculates \$\\bra{\\phi} O \\ket{\\psi}\$ + +# Details + +Calculates the matrix element of \$O\$ +```math +\\bra{\\phi} O \\ket{\\psi} +``` +in one pass, utlizing `TensorOperations`. +""" +function LinearAlgebra.dot(ϕ::MPS, O::Union{Vector, NTuple}, ψ::MPS) #where T <: AbstractMatrix S = promote_type(eltype(ψ), eltype(ϕ), eltype(O[1])) C = ones(S, 1, 1) @@ -70,8 +98,8 @@ function LinearAlgebra.dot(ϕ::MPS, O::Vector{T}, ψ::MPS) where {T <: AbstractM M̃ = conj.(ϕ[i]) Mat = O[i] @tensor C[x, y] := M̃[β, σ, x] * Mat[σ, η] * C[β, α] * M[α, η, y] order = (α, η, β, σ) -end - return C[1] + end + tr(C) end function LinearAlgebra.dot(O::AbstractMPO, ψ::T) where {T <: AbstractMPS} @@ -83,19 +111,19 @@ function LinearAlgebra.dot(O::AbstractMPO, ψ::T) where {T <: AbstractMPS} W = O[i] M = ψ[i] - @reduce N[(x, a), (y, b), σ] := sum(η) W[x, σ, y, η] * M[a, η, b] + @reduce N[(x, a), σ, (y, b)] := sum(η) W[x, σ, y, η] * M[a, η, b] ϕ[i] = N end - return ϕ + ϕ end -function dot!(O::AbstractMPO, ψ::AbstractMPS) +function dot!(ψ::AbstractMPS, O::AbstractMPO) L = length(ψ) for i in 1:L W = O[i] M = ψ[i] - @reduce N[(x, a), (y, b), σ] := sum(η) W[x, σ, y, η] * M[a, η, b] + @reduce N[(x, a), σ, (y, b)] := sum(η) W[x, σ, y, η] * M[a, η, b] ψ[i] = N end end @@ -114,7 +142,7 @@ function LinearAlgebra.dot(O1::AbstractMPO, O2::AbstractMPO) @reduce V[(x, a), σ, (y, b), η] := sum(γ) W1[x, σ, y, γ] * W2[a, γ, b, η] O[i] = V end - return O + O end Base.:(*)(O1::AbstractMPO, O2::AbstractMPO) = dot(O1, O2) diff --git a/src/cuda/base.jl b/src/cuda/base.jl index 2b77c2dd..61b165d5 100644 --- a/src/cuda/base.jl +++ b/src/cuda/base.jl @@ -63,4 +63,26 @@ function CuMPS(O::CuMPO) ψ[i] = A end return ψ -end \ No newline at end of file +end + +function is_left_normalized(ψ::CuMPS) + for i ∈ 1:length(ψ) + A = ψ[i] + DD = size(A, 3) + + @cutensor Id[x, y] := conj(A[α, σ, x]) * A[α, σ, y] order = (α, σ) + I(DD) ≈ Id ? () : return false + end + true +end + +function is_right_normalized(ϕ::CuMPS) + for i ∈ 1:length(ϕ) + B = ϕ[i] + DD = size(B, 1) + + @cutensor Id[x, y] := B[x, σ, α] * conj(B[y, σ, α]) order = (α, σ) + I(DD) ≈ Id ? () : return false + end + true +end \ No newline at end of file diff --git a/src/ising.jl b/src/ising.jl index cdc5f454..b5e78806 100644 --- a/src/ising.jl +++ b/src/ising.jl @@ -1,22 +1,36 @@ export ising_graph, energy export gibbs_tensor export GibbsControl +export brute_force +export brute_force_lazy struct GibbsControl β::Number β_schedule::Vector{<:Number} end -function gibbs_tensor(ig::MetaGraph, opts::GibbsControl) +function brute_force_lazy(ig::MetaGraph, k::Int=1) L = nv(ig) - β = opts.β + states = product(fill([-1, 1], L)...) + energies = vec(energy.(states, Ref(ig))) + perm = partialsortperm(energies, 1:k) + collect.(states)[perm], energies[perm] +end - all_states = product(fill([-1, 1], L)...) - rank = fill(2, L) +function brute_force(ig::MetaGraph, k::Int=1) + L = nv(ig) + states = ising.(digits.(0:2^L-1, base=2, pad=L)) + energies = energy.(states, Ref(ig)) + perm = partialsortperm(energies, 1:k) + states[perm], energies[perm] +end - r = exp.(-β * energy.(all_states, Ref(ig))) - ρ = reshape(r, rank...) - ρ / sum(ρ) +function gibbs_tensor(ig::MetaGraph, opts::GibbsControl) + L = nv(ig) + β = opts.β + states = product(fill([-1, 1], L)...) + ρ = exp.(-β .* energy.(states, Ref(ig))) + ρ ./ sum(ρ) end @@ -35,8 +49,8 @@ function energy(σ::Union{Vector, NTuple}, ig::MetaGraph) # linear for i ∈ vertices(ig) - h = get_prop(ig, i, :h) - energy += h * σ[i] + h = get_prop(ig, i, :h) + energy += h * σ[i] end -energy end @@ -47,7 +61,7 @@ Create a graph that represents the Ising Hamiltonian. function ising_graph(instance::String, L::Int, β::Number=1) # load the Ising instance - ising = CSV.File(instance, types=[Int, Int, Float64]) + ising = CSV.File(instance, types=[Int, Int, Float64], comment = "#") ig = MetaGraph(L, 0.0) set_prop!(ig, :description, "The Ising model.") @@ -63,6 +77,13 @@ function ising_graph(instance::String, L::Int, β::Number=1) end end + # by default h should be zero + for i ∈ 1:nv(ig) + if !has_prop(ig, i, :h) + set_prop!(ig, i, :h, 0.) || error("Cannot set bias at node $(i).") + end + end + # state and corresponding energy state = 2(rand(L) .< 0.5) .- 1 @@ -73,4 +94,9 @@ function ising_graph(instance::String, L::Int, β::Number=1) set_prop!(ig, :β, β) ig +end + +function unique_neighbors(ig::MetaGraph, i::Int) + nbrs = neighbors(ig::MetaGraph, i::Int) + filter(j -> j > i, nbrs) end \ No newline at end of file diff --git a/src/search.jl b/src/search.jl index 617a0da2..0eebbf12 100644 --- a/src/search.jl +++ b/src/search.jl @@ -1,106 +1,159 @@ export MPS_from_gates, unique_neighbors export MPSControl +export spectrum struct MPSControl max_bond::Int - ϵ::Number + var_ϵ::Number max_sweeps::Int end -function unique_neighbors(ig::MetaGraph, i::Int) - nbrs = neighbors(ig::MetaGraph, i::Int) - filter(j -> j > i, nbrs) +# ρ needs to be in the right canonical form +function spectrum(ψ::MPS, keep::Int) + @assert keep > 0 "Number of states has to be > 0" + T = eltype(ψ) + + keep_extra = keep + pCut = prob = 0. + k = 1 + + if keep < prod(rank(ψ)) + keep_extra += 1 + end + + states = fill([], 1, k) + left_env = ones(T, 1, 1, k) + + for (i, M) ∈ enumerate(ψ) + _, d, β = size(M) + + pdo = zeros(T, k, d) + LL = zeros(T, β, β, k, d) + config = zeros(Int, i, k, d) + + for j ∈ 1:k + L = left_env[:, :, j] + + for σ ∈ local_basis(d) + m = idx(σ) + + LL[:, :, j, m] = M[:, m, :]' * (L * M[:, m, :]) + pdo[j, m] = tr(LL[:, :, j, m]) + config[:, j, m] = vcat(states[:, j]..., σ) + end + end + + perm = collect(1: k * d) + k = min(k * d, keep_extra) + + if k >= keep_extra + partialsortperm!(perm, vec(pdo), 1:k, rev=true) + prob = vec(pdo)[perm] + pCut < last(prob) ? pCut = last(prob) : () + end + + @cast A[α, β, (l, d)] |= LL[α, β, l, d] + left_env = A[:, :, perm] + + @cast B[α, (l, d)] |= config[α, l, d] + states = B[:, perm] + end + states[:, 1:keep], prob[1:keep], pCut end function _apply_bias!(ψ::AbstractMPS, ig::MetaGraph, dβ::Number, i::Int) M = ψ[i] - if has_prop(ig, i, :h) - h = get_prop(ig, i, :h) - v = [exp(-0.5 * dβ * h * σ) for σ ∈ [-1, 1]] - @cast M[x, σ, y] = M[x, σ, y] * v[σ] - end + d = size(M, 2) + + h = get_prop(ig, i, :h) + v = [exp(0.5 * dβ * h * σ) for σ ∈ local_basis(d)] + + @cast M[x, σ, y] = M[x, σ, y] * v[σ] ψ[i] = M end function _apply_exponent!(ψ::AbstractMPS, ig::MetaGraph, dβ::Number, i::Int, j::Int) - δ = I(2) M = ψ[j] - J = get_prop(ig, i, j, :J) - C = [exp(-0.5 * dβ * k * J * l) for k ∈ [-1, 1], l ∈ [-1, 1]] + d = size(M, 2) + basis = local_basis(d) + + J = get_prop(ig, i, j, :J) + C = [exp(0.5 * dβ * k * J * l) for k ∈ basis, l ∈ basis] + D = I(d) if j == length(ψ) - @cast M̃[(x, a), σ, b] := C[x, σ] * δ[x, 1] * M[a, σ, b] + @cast M̃[(x, a), σ, b] := C[σ, x] * M[a, σ, b] else - @cast M̃[(x, a), σ, (y, b)] := C[x, σ] * δ[x, y] * M[a, σ, b] - end + @cast M̃[(x, a), σ, (y, b)] := C[σ, x] * D[x, y] * M[a, σ, b] + end ψ[j] = M̃ end function _apply_projector!(ψ::AbstractMPS, i::Int) - δ = I(2) M = ψ[i] + D = I(size(M, 2)) - if i == 1 - @cast M̃[a, σ, (y, b)] := δ[σ, y] * δ[1, y] * M[a, σ, b] - else - @cast M̃[(x, a), σ, (y, b)] := δ[σ, y] * δ[x, y] * M[a, σ, b] - end + @cast M̃[a, σ, (y, b)] := D[σ, y] * M[a, σ, b] ψ[i] = M̃ end -function _apply_nothing!(ψ::AbstractMPS, i::Int) - δ = I(2) +function _apply_nothing!(ψ::AbstractMPS, i::Int) M = ψ[i] + D = I(size(M, 2)) - if i == 1 - @cast M̃[a, σ, (y, b)] := δ[1, y] * M[a, σ, b] - elseif i == length(ψ) - @cast M̃[(x, a), σ, b] := δ[x, 1] * M[a, σ, b] - else - @cast M̃[(x, a), σ, (y, b)] := δ[x, y] * M[a, σ, b] - end + @cast M̃[(x, a), σ, (y, b)] := D[x, y] * M[a, σ, b] ψ[i] = M̃ end +_holes(nbrs::Vector) = setdiff(first(nbrs) : last(nbrs), nbrs) + function MPS(ig::MetaGraph, mps::MPSControl, gibbs::GibbsControl) L = nv(ig) - # control for MPS Dcut = mps.max_bond - tol = mps.ϵ + tol = mps.var_ϵ max_sweeps = mps.max_sweeps + @info "Set control parameters for MPS" Dcut tol max_sweeps - # control for Gibbs state β = gibbs.β schedule = gibbs.β_schedule @assert β ≈ sum(schedule) "Incorrect β schedule." - # prepare ~ Hadamard state as MPS - prod_state = fill([1., 1.], nv(ig)) - ρ = MPS(prod_state) + @info "Preparing Hadamard state as MPS" + ρ = HadamardMPS(L) + is_right = true + @info "Sweeping through β and σ" schedule for dβ ∈ schedule, i ∈ 1:L _apply_bias!(ρ, ig, dβ, i) + is_right = false nbrs = unique_neighbors(ig, i) if !isempty(nbrs) + _apply_projector!(ρ, i) + for j ∈ nbrs _apply_exponent!(ρ, ig, dβ, i, j) end - _apply_projector!(ρ, i) - - for l ∈ setdiff(1:L, union(i, nbrs)) - _apply_nothing!(ρ, l) + for l ∈ _holes(nbrs) + _apply_nothing!(χ, l) end end - # reduce bond dimension if bond_dimension(ρ) > Dcut + @info "Compresing MPS" bond_dimension(ρ), Dcut ρ = compress(ρ, Dcut, tol, max_sweeps) + is_right = true end end + + if !is_right + canonise!(ρ, :right) + is_right = true + end ρ -end \ No newline at end of file +end + diff --git a/src/tests/notation_tests.jl b/src/tests/notation_tests.jl index 6e14849e..5769a641 100644 --- a/src/tests/notation_tests.jl +++ b/src/tests/notation_tests.jl @@ -200,7 +200,6 @@ end M = rand(10,10) m = matrix2qubo_vec(M) - end @testset "test notation" begin diff --git a/src/utils.jl b/src/utils.jl index 1c67a92a..5ade3d28 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,3 +1,59 @@ -# export newdim +export idx, ising, proj +export HadamardMPS, rq +export all_states, local_basis -# newdim(::Type{T}, dims) where {T<:AbstractArray} = T.name.wrapper{eltype(T), dims} +ising(σ::Union{Vector, NTuple}) = 2 .* σ .- 1 + +idx(σ::Int) = (σ == -1) ? 1 : σ + 1 +_σ(idx::Int) = (idx == 1) ? -1 : idx - 1 + +local_basis(d::Int) = union(-1, 1:d-1) + +function proj(state, dims::T) where {T <: Union{Vector, NTuple}} + P = Matrix{Float64}[] + for (σ, r) ∈ zip(state, dims) + v = zeros(r) + v[idx(σ)...] = 1. + push!(P, v * v') + end + P +end + +function all_states(rank::T) where T <: Union{Vector, NTuple} + basis = [local_basis(r) for r ∈ rank] + product(basis...) +end + +function HadamardMPS(rank::T) where T <: Union{Vector, NTuple} + MPS(fill(1, r) ./ sqrt(r) for r ∈ rank) +end +HadamardMPS(L::Int) = MPS(fill([1., 1.] / sqrt(2), L)) + +function LinearAlgebra.qr(M::AbstractMatrix, Dcut::Int, args...) + fact = pqrfact(M, rank=Dcut, args...) + Q = fact[:Q] + R = fact[:R] + return _qr_fix(Q, R) +end + +function rq(M::AbstractMatrix, Dcut::Int, args...) + fact = pqrfact(:c, conj.(M), rank=Dcut, args...) + 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 + +function LinearAlgebra.svd(A::AbstractMatrix, Dcut::Int, args...) + U, Σ, V = psvd(A, rank=Dcut, args...) + d = diag(U) + ph = d ./ abs.(d) + return U * Diagonal(ph), Σ, V * Diagonal(ph) +end \ No newline at end of file diff --git a/test/base.jl b/test/base.jl index 067875b6..e178e843 100644 --- a/test/base.jl +++ b/test/base.jl @@ -7,7 +7,7 @@ T = ComplexF64 @testset "Random MPS" begin ψ = randn(MPS{T}, sites, D, d) - @test _verify_bonds(ψ) == nothing + @test verify_bonds(ψ) == nothing @test ψ == ψ @test ψ ≈ ψ @@ -15,6 +15,8 @@ T = ComplexF64 @test length(ψ) == sites @test size(ψ) == (sites, ) @test eltype(ψ) == ComplexF64 + @test rank(ψ) == Tuple(fill(d, 1:sites)) + @test bond_dimension(ψ) ≈ D ϕ = copy(ψ) @test ϕ == ψ @@ -23,25 +25,6 @@ T = ComplexF64 show(ψ) end -@testset "Retrieving bond dimension" begin - S = Float64 - size = 5 - D = 6 - d = 2 - ψ = randn(MPS{S}, sites, D, d) - - @test bond_dimension(ψ) ≈ D -end - -@testset "MPS from Product state" begin - L = 3 - - prod_state = [ [1, 1.] / sqrt(2) for _ ∈ 1:L] - ϕ = MPS(prod_state) - - show(ϕ) -end - @testset "Random MPO" begin O = randn(MPO{T}, sites, D, d) @@ -57,4 +40,81 @@ end @test P ≈ O end +@testset "Reshaping (row-wise)" begin + vec = Vector(1:6) + + A = reshape_row(vec, (2, 3)) + B = [1 2 3; 4 5 6] + + @test A ≈ B +end + +@testset "Basic vector to tensor reshaping" begin + dims = (2, 3, 4, 5) + states = [randn(T, d) for d ∈ dims] + vec = kron(states...) + + ψ = tensor(MPS(states)) + ϕ = reshape_row(vec, dims) + + @test ψ ≈ ϕ +end + +@testset "MPS from tensor" begin + ϵ = 1E-14 + + dims = (2,3,4,3,5) + sites = length(dims) + A = randn(T, dims) + + @test sqrt(sum(abs.(A) .^ 2)) ≈ norm(A) + + @test ndims(A) == sites + @test size(A) == dims + + ψ = MPS(A, :right) + + @test norm(ψ) ≈ 1 + @test_nowarn verify_bonds(ψ) + @test_nowarn verify_physical_dims(ψ, dims) + @test is_right_normalized(ψ) + show(ψ) + + AA = tensor(ψ) + + @test rank(ψ) == size(AA) + @test norm(AA) ≈ 1 + @test size(AA) == size(A) + + vA = vec(A) + nA = norm(vA) + @test abs(1 - abs(dot(vec(AA), vA ./ nA))) < ϵ + #@test AA ≈ A ./ norm(A) # this is true "module phase" + + B = randn(T, dims...) + ϕ = MPS(B, :left) + + @test norm(ϕ) ≈ 1 + @test_nowarn verify_bonds(ϕ) + @test_nowarn verify_physical_dims(ϕ, dims) + @test is_left_normalized(ϕ) + show(ϕ) + + BB = tensor(ϕ) + + @test rank(ϕ) == size(BB) + @test norm(BB) ≈ 1 + @test sqrt(sum(abs.(B) .^ 2)) ≈ norm(B) + + vB = vec(B) + nB = norm(vB) + @test abs(1 - abs(dot(vec(BB), vB ./ nB))) < ϵ + #@test BB ≈ B ./ norm(B) # this is true "module phase" + + χ = MPS(A, :left) + + @test norm(χ) ≈ 1 + @test abs(1 - abs(dot(ψ, χ))) < ϵ +end + end \ No newline at end of file diff --git a/test/compressions.jl b/test/compressions.jl index 492320f2..abe7b89f 100644 --- a/test/compressions.jl +++ b/test/compressions.jl @@ -18,34 +18,14 @@ T = Float64 @testset "Canonisation (left)" begin canonise!(ψ, :left) show(ψ) - - is_left_normalized = true - for i ∈ 1:length(ψ) - A = ψ[i] - DD = size(A, 3) - - @tensor Id[x, y] := conj(A[α, σ, x]) * A[α, σ, y] order = (α, σ) - is_left_normalized *= I(DD) ≈ Id - end - - @test is_left_normalized + @test is_left_normalized(ψ) @test dot(ψ, ψ) ≈ 1 end @testset "Canonisation (right)" begin canonise!(ϕ, :right) show(ϕ) - - is_right_normalized = true - for i ∈ 1:length(ϕ) - B = ϕ[i] - DD = size(B, 1) - - @tensor Id[x, y] := B[x, σ, α] * conj(B[y, σ, α]) order = (α, σ) - is_right_normalized *= I(DD) ≈ Id - end - - @test is_right_normalized + @test is_right_normalized(ϕ) @test dot(ϕ, ϕ) ≈ 1 end @@ -71,6 +51,21 @@ end @test dot(ψ, ψ) ≈ 1 end +@testset "" begin + ϵ = 1E-14 + ψ = randn(MPS{T}, sites, D, d) + + l = copy(ψ) + r = copy(l) + canonise!(l, :left) + canonise!(r, :right) + + @test dot(l, l) ≈ 1 + @test dot(r, r) ≈ 1 + + @test abs(1 - abs(dot(l, r))) < ϵ +end + @testset "Variational compression" begin Dcut = 5 tol = 1E-4 @@ -88,45 +83,12 @@ end println("(Φ, Φ) = ", dot(Φ, Φ)) overlap = dot(Ψ, Φ) - dist1 = 2 - 2 * real(overlap) - dist2 = norm(Ψ)^2 + norm(Φ)^2 - 2 * real(overlap) + dist1 = 2 - 2 * abs(overlap) + dist2 = norm(Ψ)^2 + norm(Φ)^2 - 2 * abs(overlap) @test abs(dist1 - dist2) < 1e-14 println("(Φ, Ψ) = ", overlap) println("dist(Φ, Ψ)^2 = ", dist2) end - -@testset "Compare with Krzysiek's implementation" begin - Dcut = 2 - tol = 1E-4 - max_sweeps = 5 - - canonise!(Φ, :right) - @test dot(Φ, Φ) ≈ 1 - - Ψ = compress(Φ, Dcut, tol, max_sweeps) - - Φ_trunc = copy(Φ) - truncate!(Φ_trunc, :right, Dcut) - - permuted_mps = map(x->permutedims(x, (1,3,2)), Φ.tensors) - # tensors = compress_mps_itterativelly(, Φ_trunc.tensors, Dcut, tol) - tensors = compress_iter(permuted_mps, Dcut, tol) - tensors = map(x->permutedims(x, (1,3,2)), tensors) - ξ = MPS(tensors) - # ξ.tensors = tensors - - @test dot(ξ, ξ) ≈ 1 - - println("(ξ, ξ) = ", dot(ξ, ξ)) - - overlap = dot(Ψ, ξ) - dist1 = 2 - 2 * real(overlap) - dist2 = norm(Ψ)^2 + norm(ξ)^2 - 2 * real(overlap) - - @test abs(dist1 - dist2) < 1e-14 - - println("Krzysiek wins - flawless victory, fatality.") -end end \ No newline at end of file diff --git a/test/contractions.jl b/test/contractions.jl index d69df337..a2f20bc8 100644 --- a/test/contractions.jl +++ b/test/contractions.jl @@ -18,6 +18,12 @@ Id = [I(d) for _ ∈ 1:length(ϕ)] @test dot(ψ, Id, ϕ) ≈ dot(ψ, ϕ) @test norm(ψ) ≈ sqrt(abs(dot(ψ, ψ))) + + ψ[end] *= 1/norm(ψ) + @test dot(ψ, ψ) ≈ 1 + + ϕ[1] *= 1/norm(ϕ) + @test dot(ϕ, ϕ) ≈ 1 end @testset "left environment" begin diff --git a/test/instances/128_001.txt b/test/instances/128_001.txt new file mode 100644 index 00000000..6e6454bd --- /dev/null +++ b/test/instances/128_001.txt @@ -0,0 +1,481 @@ +# +1 1 0.200000 +2 2 0.146667 +3 3 -0.173333 +4 4 -0.200000 +5 5 0.013333 +6 6 0.040000 +7 7 0.013333 +8 8 0.120000 +9 9 -0.093333 +10 10 -0.013333 +11 11 0.120000 +12 12 0.173333 +13 13 0.093333 +14 14 0.040000 +15 15 -0.146667 +16 16 -0.146667 +17 17 -0.146667 +18 18 -0.013333 +19 19 0.173333 +20 20 -0.066667 +21 21 0.066667 +22 22 -0.120000 +23 23 -0.173333 +24 24 0.120000 +25 25 -0.146667 +26 26 -0.093333 +27 27 0.146667 +28 28 -0.173333 +29 29 -0.120000 +30 30 -0.146667 +31 31 0.120000 +32 32 -0.120000 +33 33 -0.200000 +34 34 -0.173333 +35 35 -0.066667 +36 36 -0.120000 +37 37 -0.040000 +38 38 0.173333 +39 39 0.040000 +40 40 0.120000 +41 41 -0.066667 +42 42 0.120000 +43 43 0.066667 +44 44 -0.173333 +45 45 0.120000 +46 46 0.120000 +47 47 -0.013333 +48 48 -0.120000 +49 49 -0.120000 +50 50 0.146667 +51 51 -0.066667 +52 52 0.066667 +53 53 0.146667 +54 54 0.120000 +55 55 -0.173333 +56 56 -0.146667 +57 57 -0.066667 +58 58 0.173333 +59 59 -0.173333 +60 60 -0.200000 +61 61 0.093333 +62 62 0.200000 +63 63 -0.200000 +64 64 -0.040000 +65 65 0.146667 +66 66 -0.040000 +67 67 0.120000 +68 68 -0.146667 +69 69 -0.120000 +70 70 -0.093333 +71 71 -0.013333 +72 72 -0.013333 +73 73 0.200000 +74 74 -0.173333 +75 75 -0.173333 +76 76 -0.200000 +77 77 0.066667 +78 78 0.200000 +79 79 -0.146667 +80 80 -0.173333 +81 81 -0.146667 +82 82 0.120000 +83 83 0.120000 +84 84 0.173333 +85 85 -0.013333 +86 86 -0.146667 +87 87 0.200000 +88 88 0.013333 +89 89 -0.173333 +90 90 0.200000 +91 91 0.173333 +92 92 -0.040000 +93 93 -0.120000 +94 94 -0.120000 +95 95 0.120000 +96 96 0.066667 +97 97 -0.066667 +98 98 -0.173333 +99 99 0.200000 +100 100 -0.066667 +101 101 -0.013333 +102 102 0.066667 +103 103 0.013333 +104 104 0.200000 +105 105 -0.120000 +106 106 0.120000 +107 107 -0.093333 +108 108 -0.013333 +109 109 -0.093333 +110 110 -0.173333 +111 111 -0.120000 +112 112 0.093333 +113 113 -0.120000 +114 114 -0.040000 +115 115 0.120000 +116 116 0.093333 +117 117 -0.146667 +118 118 -0.066667 +119 119 -0.066667 +120 120 -0.040000 +121 121 -0.200000 +122 122 -0.120000 +123 123 -0.093333 +124 124 -0.013333 +125 125 0.066667 +126 126 -0.200000 +127 127 0.173333 +128 128 -0.146667 +1 5 0.200000 +1 6 -0.333333 +1 7 0.866667 +1 8 -1.000000 +1 33 -0.200000 +2 5 -0.466667 +2 6 1.000000 +2 7 0.600000 +2 8 -0.066667 +2 34 4.333333 +3 5 0.466667 +3 6 -0.600000 +3 7 -0.866667 +3 8 -0.466667 +3 35 -0.066667 +4 5 0.600000 +4 6 -1.000000 +4 7 0.066667 +4 8 -3.666667 +4 36 -0.466667 +5 13 0.333333 +6 14 -0.733333 +7 15 1.000000 +8 16 -0.333333 +9 13 -0.600000 +9 14 0.866667 +9 15 -0.600000 +9 16 -0.466667 +9 41 -0.200000 +10 13 0.466667 +10 14 -0.600000 +10 15 0.733333 +10 16 -0.866667 +10 42 -1.000000 +11 13 0.200000 +11 14 0.333333 +11 15 1.000000 +11 16 -0.866667 +11 43 0.333333 +12 13 -0.866667 +12 14 0.066667 +12 15 0.466667 +12 16 0.600000 +12 44 -0.333333 +13 21 -0.466667 +14 22 0.466667 +15 23 -0.200000 +16 24 0.733333 +17 21 3.000000 +17 22 0.200000 +17 23 0.066667 +17 24 -0.200000 +17 49 -1.000000 +18 21 2.333333 +18 22 -0.866667 +18 23 2.333333 +18 24 1.000000 +18 50 0.866667 +19 21 -1.666667 +19 22 -0.866667 +19 23 -0.466667 +19 24 -1.000000 +19 51 -0.466667 +20 21 -0.333333 +20 22 -0.066667 +20 23 -0.866667 +20 24 -0.600000 +20 52 -1.000000 +21 29 -0.066667 +22 30 -0.066667 +23 31 0.466667 +24 32 1.000000 +25 29 0.200000 +25 30 0.733333 +25 31 -1.000000 +25 32 -0.066667 +25 57 0.466667 +26 29 -0.466667 +26 30 -0.733333 +26 31 0.733333 +26 32 -0.866667 +26 58 0.600000 +27 29 -0.066667 +27 30 -0.333333 +27 31 -1.000000 +27 32 -0.600000 +27 59 0.333333 +28 29 0.733333 +28 30 -0.066667 +28 31 -0.600000 +28 32 -0.733333 +28 60 0.333333 +33 37 0.733333 +33 38 4.333333 +33 39 0.600000 +33 40 -0.200000 +33 65 0.866667 +34 37 0.866667 +34 38 -1.000000 +34 39 -0.200000 +34 40 0.733333 +34 66 0.333333 +35 37 -0.333333 +35 38 -0.333333 +35 39 0.600000 +35 40 0.733333 +35 67 -0.600000 +36 37 0.466667 +36 38 0.733333 +36 39 0.200000 +36 40 -1.000000 +36 68 -0.333333 +37 45 0.333333 +38 46 -0.733333 +39 47 5.000000 +40 48 0.333333 +41 45 0.600000 +41 46 -0.333333 +41 47 -0.466667 +41 48 -4.333333 +41 73 -0.066667 +42 45 -0.600000 +42 46 -0.200000 +42 47 0.733333 +42 48 -0.466667 +42 74 0.866667 +43 45 -0.733333 +43 46 0.333333 +43 47 -2.333333 +43 48 0.200000 +43 75 0.333333 +44 45 0.866667 +44 46 -1.000000 +44 47 -0.866667 +44 48 -1.000000 +44 76 0.200000 +45 53 0.866667 +46 54 -0.200000 +47 55 0.333333 +48 56 -0.866667 +49 53 -0.466667 +49 54 -0.733333 +49 55 -0.066667 +49 56 -0.066667 +49 81 -0.600000 +50 53 -1.000000 +50 54 -1.000000 +50 55 -0.066667 +50 56 -0.600000 +50 82 3.000000 +51 53 1.000000 +51 54 0.200000 +51 55 0.200000 +51 56 0.733333 +51 83 -0.600000 +52 53 -0.733333 +52 54 -0.733333 +52 55 0.066667 +52 56 -0.733333 +52 84 0.733333 +53 61 0.733333 +54 62 0.733333 +55 63 -0.066667 +56 64 -2.333333 +57 61 5.000000 +57 62 0.333333 +57 63 -0.866667 +57 64 -0.466667 +57 89 1.000000 +58 61 -0.200000 +58 62 -0.200000 +58 63 0.066667 +58 64 0.466667 +58 90 -0.733333 +59 61 -0.733333 +59 62 -0.200000 +59 63 0.333333 +59 64 -0.733333 +59 91 -0.200000 +60 61 -0.333333 +60 62 -0.600000 +60 63 -0.200000 +60 64 5.000000 +60 92 0.866667 +65 69 -3.000000 +65 70 1.000000 +65 71 1.000000 +65 72 0.733333 +65 97 0.466667 +66 69 1.000000 +66 70 -0.333333 +66 71 -0.066667 +66 72 0.066667 +66 98 0.733333 +67 69 -0.733333 +67 70 0.600000 +67 71 -3.666667 +67 72 4.333333 +67 99 0.600000 +68 69 -0.866667 +68 70 0.466667 +68 71 -0.333333 +68 72 -0.600000 +68 100 0.200000 +69 77 -0.600000 +70 78 -0.600000 +71 79 0.866667 +72 80 -0.200000 +73 77 -0.200000 +73 78 -1.666667 +73 79 0.333333 +73 80 -0.466667 +73 105 0.733333 +74 77 -0.066667 +74 78 0.066667 +74 79 0.200000 +74 80 -0.200000 +74 106 0.466667 +75 77 0.866667 +75 78 0.333333 +75 79 -0.066667 +75 80 3.000000 +75 107 -1.000000 +76 77 1.666667 +76 78 1.000000 +76 79 -0.600000 +76 80 0.333333 +76 108 1.000000 +77 85 -0.866667 +78 86 0.866667 +79 87 0.466667 +80 88 0.733333 +81 85 -1.000000 +81 86 1.000000 +81 87 -0.466667 +81 88 3.666667 +81 113 -4.333333 +82 85 0.200000 +82 86 -0.200000 +82 87 0.866667 +82 88 -0.733333 +82 114 0.333333 +83 85 0.733333 +83 86 -0.200000 +83 87 0.066667 +83 88 0.333333 +83 115 -0.200000 +84 85 0.200000 +84 86 0.066667 +84 87 -0.866667 +84 88 -0.733333 +84 116 -0.333333 +85 93 0.333333 +86 94 -1.000000 +87 95 -0.866667 +88 96 -0.200000 +89 93 -0.066667 +89 94 0.600000 +89 95 -0.866667 +89 96 1.000000 +89 121 0.333333 +90 93 -0.466667 +90 94 0.733333 +90 95 1.000000 +90 96 -0.466667 +90 122 0.866667 +91 93 0.733333 +91 94 0.200000 +91 95 0.333333 +91 96 -0.733333 +91 123 -0.600000 +92 93 -0.600000 +92 94 -0.466667 +92 95 0.200000 +92 96 -0.733333 +92 124 0.600000 +97 101 -0.866667 +97 102 -0.333333 +97 103 -0.466667 +97 104 0.733333 +98 101 -2.333333 +98 102 -3.000000 +98 103 0.866667 +98 104 0.733333 +99 101 1.000000 +99 102 -0.866667 +99 103 0.733333 +99 104 -0.733333 +100 101 0.866667 +100 102 -0.066667 +100 103 -1.000000 +100 104 -0.200000 +101 109 0.466667 +102 110 0.600000 +103 111 0.066667 +104 112 0.600000 +105 109 -3.666667 +105 110 1.000000 +105 111 -1.000000 +105 112 0.600000 +106 109 0.466667 +106 110 0.733333 +106 111 0.600000 +106 112 0.866667 +107 109 3.666667 +107 110 -0.466667 +107 111 0.600000 +107 112 0.333333 +108 109 0.866667 +108 110 0.066667 +108 111 0.733333 +108 112 -0.733333 +109 117 0.333333 +110 118 -0.600000 +111 119 0.866667 +112 120 0.333333 +113 117 5.000000 +113 118 -0.866667 +113 119 -1.666667 +113 120 -0.600000 +114 117 1.000000 +114 118 0.333333 +114 119 0.333333 +114 120 -0.600000 +115 117 -0.466667 +115 118 -0.466667 +115 119 0.600000 +115 120 0.066667 +116 117 0.600000 +116 118 0.733333 +116 119 -0.600000 +116 120 0.333333 +117 125 -1.000000 +118 126 0.466667 +119 127 -0.466667 +120 128 0.466667 +121 125 -0.200000 +121 126 0.600000 +121 127 0.066667 +121 128 4.333333 +122 125 0.066667 +122 126 -0.600000 +122 127 1.000000 +122 128 0.333333 +123 125 0.466667 +123 126 -0.733333 +123 127 -0.066667 +123 128 0.600000 +124 125 -0.333333 +124 126 -0.466667 +124 127 -0.200000 +124 128 -0.866667 \ No newline at end of file diff --git a/test/instances/16_001.txt b/test/instances/16_001.txt new file mode 100644 index 00000000..35b36c65 --- /dev/null +++ b/test/instances/16_001.txt @@ -0,0 +1,36 @@ +# 4x4 +1 1 0.062500 +3 3 0.062500 +6 6 -0.062500 +7 7 0.062500 +8 8 0.031250 +9 9 -0.031250 +10 10 -0.031250 +12 12 -0.062500 +13 13 0.062500 +14 14 0.062500 +15 15 0.062500 +16 16 -0.031250 +1 2 0.187500 +1 5 -0.156250 +2 6 0.937500 +3 4 0.187500 +3 7 -0.031250 +4 8 -0.156250 +5 6 -1.000000 +5 9 0.156250 +6 7 0.187500 +6 10 0.250000 +7 8 0.093750 +7 11 0.156250 +8 12 0.187500 +9 10 0.093750 +9 13 -0.187500 +10 11 0.156250 +10 14 0.250000 +11 12 -0.156250 +11 15 -0.218750 +12 16 0.062500 +13 14 0.156250 +14 15 0.218750 +15 16 0.250000 \ No newline at end of file diff --git a/test/instances/4_001.txt b/test/instances/4_001.txt new file mode 100644 index 00000000..6efa7550 --- /dev/null +++ b/test/instances/4_001.txt @@ -0,0 +1,13 @@ +# bias +1 1 0.1 +2 2 0.2 +3 3 0.3 +4 4 0.4 +# square +1 2 0.7 +1 3 0.4 +2 4 0.3 +3 4 0.5 +# diag +1 4 0.6 +2 3 0.4 diff --git a/test/instances/4_002.txt b/test/instances/4_002.txt new file mode 100644 index 00000000..b8f2e15e --- /dev/null +++ b/test/instances/4_002.txt @@ -0,0 +1,6 @@ +# square +1 2 0.7 +1 3 0.0 +2 4 0.0 +3 4 0.0 + diff --git a/test/lattice_3.txt b/test/instances/9_001.txt similarity index 95% rename from test/lattice_3.txt rename to test/instances/9_001.txt index 52761369..784a75bc 100644 --- a/test/lattice_3.txt +++ b/test/instances/9_001.txt @@ -1,4 +1,4 @@ -i j v +# i j v 1 2 0.1 2 3 0.3 4 5 0.7 diff --git a/test/instances/Egs.txt b/test/instances/Egs.txt new file mode 100644 index 00000000..3af03f42 --- /dev/null +++ b/test/instances/Egs.txt @@ -0,0 +1,3 @@ +# ins E_gs +16_001 -4.6875 +128_001 -210.9333333 diff --git a/test/ising.jl b/test/ising.jl index 7ade12e5..4416a9e3 100644 --- a/test/ising.jl +++ b/test/ising.jl @@ -1,13 +1,13 @@ using MetaGraphs using LightGraphs using GraphPlot -using Base @testset "Ising" begin - L = 3 - N = L^2 - instance = "./lattice_$L.txt" + L = 4 + N = L^2 + instance = "$(@__DIR__)/instances/$(N)_001.txt" + ig = ising_graph(instance, N) E = get_prop(ig, :energy) @@ -40,4 +40,42 @@ using Base @test B+B' == A gplot(ig, nodelabel=1:N) + + @testset "Naive brute force" begin + k = 2^N + + states, energies = brute_force(ig, k) + + display(states[1:5]) + println(" ") + display(energies[1:5]) + println(" ") + + @test energies ≈ energy.(states, Ref(ig)) + + states_lazy, energies_lazy = brute_force_lazy(ig, k) + + @test energies_lazy ≈ energies + @test states_lazy == states + + if k == 2^N + + β = rand(Float64) + opts = GibbsControl(β, [β]) + + ρ = gibbs_tensor(ig, opts) + @test size(ρ) == Tuple(fill(2, N)) + + r = exp.(-β .* energies) + R = r ./ sum(r) + + @test sum(R) ≈ 1 + @test sum(ρ) ≈ 1 + + @test maximum(R) ≈ maximum(ρ) + @test minimum(R) ≈ minimum(ρ) + + @test [ρ[idx.(σ)...] for σ ∈ states] ≈ R + end + end end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 91584bf0..d4fa65a4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,8 +5,15 @@ 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() +if CUDA.functional() && CUDA.has_cutensor() && false push!(my_tests, "cuda/base.jl", "cuda/contractions.jl", @@ -21,7 +28,7 @@ push!(my_tests, "ising.jl", "search.jl" ) - + for my_test in my_tests include(my_test) end diff --git a/test/search.jl b/test/search.jl index 9fe56b8c..a393df59 100644 --- a/test/search.jl +++ b/test/search.jl @@ -2,36 +2,158 @@ using MetaGraphs using LightGraphs using GraphPlot -L = 3 +L = 2 N = L^2 -instance = "./lattice_$L.txt" + +instance = "$(@__DIR__)/instances/$(N)_001.txt" ig = ising_graph(instance, N) -@testset "MPS from gates" begin +ϵ = 1E-14 - Dcut = 16 - var_tol=1E-8 - max_sweeps = 4 +Dcut = N^2 + 1 +var_tol = 1E-8 +max_sweeps = 4 - β = 1 - dβ = 0.25 - β_schedule = [dβ for _ ∈ 1:4] +β = 1 +dβ = 0.25 +β_schedule = [β] #[dβ for _ ∈ 1:4] - gibbs_param = GibbsControl(β, β_schedule) - mps_param = MPSControl(Dcut, var_tol, max_sweeps) +gibbs_param = GibbsControl(β, β_schedule) +mps_param = MPSControl(Dcut, var_tol, max_sweeps) - @testset "Generate ρ" begin - ρ = MPS(ig, mps_param, gibbs_param) +ϱ = gibbs_tensor(ig, gibbs_param) +@test sum(ϱ) ≈ 1 - show(ρ) - @test_nowarn _verify_bonds(ρ) - canonise!(ρ) - @test dot(ρ, ρ) ≈ 1 - end +@testset "Verifying gate operations" begin + χ = HadamardMPS(N) + + d = 2 + rank = fill(d, N) + states = all_states(rank) + T = ones(rank...) ./ d^N + + @test sum(T) ≈ 1 + + for i ∈ 1:N + SpinGlassPEPS._apply_bias!(χ, ig, β, i) + + h = get_prop(ig, i, :h) + for σ ∈ states + T[idx.(σ)...] *= exp(β * σ[i] * h) + end + + nbrs = unique_neighbors(ig, i) + + if !isempty(nbrs) + + SpinGlassPEPS._apply_projector!(χ, i) + for j ∈ nbrs + SpinGlassPEPS._apply_exponent!(χ, ig, β, i, j) - @testset "Generate Gibbs state exactly" begin - r = gibbs_tensor(ig, gibbs_param) - @test sum(r) ≈ 1 + J = get_prop(ig, i, j, :J) + for σ ∈ states + T[idx.(σ)...] *= exp(β * σ[i] * J * σ[j]) + end + end + + for l ∈ i+1:maximum(nbrs) + if l ∉ nbrs + SpinGlassPEPS._apply_nothing!(χ, l) + end + end + end + + show(χ) + verify_bonds(χ) + + @test abs(dot(χ, χ) - sum(T)) < ϵ end + + @test T ./ sum(T) ≈ ϱ +end + +@testset "MPS from gates" begin + + @testset "Exact Gibbs pure state (MPS)" begin + d = 2 + L = nv(ig) + dims = fill(d, L) + + @info "Generating Gibbs state - |ρ>" d L dims β ϵ + + states = all_states(dims) + ψ = ones(dims...) + + for σ ∈ states + for i ∈ 1:L + h = get_prop(ig, i, :h) + + nbrs = unique_neighbors(ig, i) + ψ[idx.(σ)...] *= exp(0.5 * β * h * σ[i]) + + for j ∈ nbrs + J = get_prop(ig, i, j, :J) + ψ[idx.(σ)...] *= exp(0.5 * β * σ[i] * J * σ[j]) + end + end + end + + ρ = abs.(ψ) .^ 2 + @test ρ / sum(ρ) ≈ ϱ + + @info "Generating MPS from |ρ>" + rψ = MPS(ψ) + @test dot(rψ, rψ) ≈ 1 + @test_nowarn is_right_normalized(rψ) + + lψ = MPS(ψ, :left) + @test dot(lψ, lψ) ≈ 1 + + vlψ = vec(tensor(lψ)) + vrψ = vec(tensor(rψ)) + + vψ = vec(ψ) + vψ /= norm(vψ) + + @test abs(1 - abs(dot(vlψ, vrψ))) < ϵ + @test abs(1 - abs(dot(vlψ, vψ))) < ϵ + + @info "Verifying MPS from gates" + + Gψ = MPS(ig, mps_param, gibbs_param) + + @test_nowarn is_right_normalized(Gψ) + @test bond_dimension(Gψ) > 1 + @test dot(Gψ, Gψ) ≈ 1 + @test_nowarn verify_bonds(Gψ) + + @test abs(1 - abs(dot(Gψ, rψ))) < ϵ + + @info "Verifying probabilities" L β + + for σ ∈ states + p = dot(rψ, σ) + r = dot(rψ, proj(σ, dims), rψ) + + @test p ≈ r + @test ϱ[idx.(σ)...] ≈ p + end + + for max_states ∈ [1, N, 2*N, N^2] + @info "Verifying low energy spectrum" max_states + + states, prob, pCut = spectrum(rψ, max_states) + states_bf, energies = brute_force(ig, max_states) + + @info "The largest discarded probability" pCut + @test maximum(prob) > pCut + + for (j, (p, e)) ∈ enumerate(zip(prob, energies)) + σ = states[:, j] + @test ϱ[idx.(σ)...] ≈ p + @test abs(energy(σ, ig) - e) < ϵ + end + end + end end \ No newline at end of file