Merge 2762d34 into 6e6c40b

src/base.jl

@@ -146,6 +146,17 @@ function MPS(O::MPO)
     ψ
 end  
 
+function Base.randn(::Type{MPS{T}}, D::Int, rank::Union{Vector, NTuple}) where {T}
+    L = length(rank)
+    ψ = MPS(T, L)
+    ψ[1] = randn(T, 1, rank[1], D)
+    for i ∈ 2:(L-1)
+        ψ[i] = randn(T, D, rank[i], D)
+    end
+    ψ[end] = randn(T, D, rank[end], 1)
+    ψ
+end
+
 function Base.randn(::Type{MPS{T}}, L::Int, D::Int, d::Int) where {T}
     ψ = MPS(T, L)
     ψ[1] = randn(T, 1, d, D)

src/ising.jl

@@ -1,45 +1,71 @@
 export ising_graph, energy
-export gibbs_tensor
-export GibbsControl
-export brute_force
-export brute_force_lazy
-
-struct GibbsControl 
-    β::Number
-    β_schedule::Vector{<:Number}
-end
+export gibbs_tensor, brute_force
+export State
 
-function brute_force_lazy(ig::MetaGraph, k::Int=1)
-    L = nv(ig)
-    states = product(fill([-1, 1], L)...)
+const State = Union{Vector, NTuple}
+
+"""
+$(TYPEDSIGNATURES)
+
+Return the low energy spectrum
+
+# Details
+
+Calculates \$k\$ lowest energy states 
+together with the coresponding energies 
+of a classical Ising Hamiltonian
+"""
+
+function brute_force(ig::MetaGraph, k::Int=1)
+    states = all_states(get_prop(ig, :rank))
     energies = vec(energy.(states, Ref(ig)))
     perm = partialsortperm(energies, 1:k) 
     collect.(states)[perm], energies[perm]
 end    
 
-function brute_force(ig::MetaGraph, k::Int=1)
+_ising(σ::State) = 2 .* σ .- 1
+
+function _brute_force(ig::MetaGraph, k::Int=1)
     L = nv(ig)
-    states = ising.(digits.(0:2^L-1, base=2, pad=L))
+    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  
 
-function gibbs_tensor(ig::MetaGraph, opts::GibbsControl)
-    L = nv(ig)
-    β = opts.β
-    states = product(fill([-1, 1], L)...)
-    ρ = exp.(-β .* energy.(states, Ref(ig)))
+
+"""
+$(TYPEDSIGNATURES)
+
+Calculates Gibbs state of a classical Ising Hamiltonian
+
+# Details
+
+Calculates matrix elements (probabilities) of \$\\rho\$ 
+```math
+\$\\bra{\\σ}\\rho\\ket{\\sigma}\$
+```
+for all possible configurations \$\\σ\$.
+"""
+function gibbs_tensor(ig::MetaGraph)
+    β = get_prop(ig, :β)
+    rank = get_prop(ig, :rank)
+    ρ = exp.(-β .* energy.(all_states(rank), Ref(ig)))
     ρ ./ sum(ρ)
 end
 
 
 """
-Calculate the Ising energy as E = -sum_<i,j> s_i * J_ij * s_j - sum_j h_i * s_j.
+$(TYPEDSIGNATURES)
+
+Calculate the Ising energy 
+```math
+E = -\\sum_<i,j> s_i J_{ij} * s_j - \\sum_j h_i s_j.
+```
 """
-function energy(σ::Union{Vector, NTuple}, ig::MetaGraph)
+function energy(σ::State, ig::MetaGraph)
+    energy::Float64 = 0
 
-    energy = 0
     # quadratic
     for edge ∈ edges(ig)
         i, j = src(edge), dst(edge)         
@@ -56,19 +82,28 @@ function energy(σ::Union{Vector, NTuple}, ig::MetaGraph)
 end
 
 """
-Create a graph that represents the Ising Hamiltonian.
+$(TYPEDSIGNATURES)
+
+Create the Ising spin glass model.
+
+# Details
+
+Store extra information
 """
-function ising_graph(instance::String, L::Int, β::Number=1)
+function ising_graph(instance::Union{String, Dict}, L::Int, β::Number=1)
 
     # load the Ising instance
-    ising = CSV.File(instance, types=[Int, Int, Float64], comment = "#")
-    ig = MetaGraph(L, 0.0)
+    if typeof(instance) == String
+        ising = CSV.File(instance, types=[Int, Int, Float64], comment = "#")
+    else
+        ising = [ (i, j, J) for ((i, j), J) ∈ instance ] 
+    end
 
+    ig = MetaGraph(L, 0.0)
     set_prop!(ig, :description, "The Ising model.")
 
     # setup the model (J_ij, h_i)
-    for row ∈ ising 
-        i, j, v = row
+    for (i, j, v) ∈ ising 
         if i == j
             set_prop!(ig, i, :h, v) || error("Node $i missing!")
         else
@@ -84,18 +119,29 @@ function ising_graph(instance::String, L::Int, β::Number=1)
         end 
     end
 
-    # state and corresponding energy
+    # state (random by default) and corresponding energy
     state = 2(rand(L) .< 0.5) .- 1
 
     set_prop!(ig, :state, state)
     set_prop!(ig, :energy, energy(state, ig)) || error("Unable to calculate the Ising energy!")
 
-    # store extra information
+    # store extra information 
     set_prop!(ig, :β, β)
-
+    set_prop!(ig, :rank, fill(2, L))
+
     ig
 end
 
+"""
+$(TYPEDSIGNATURES)
+
+Calculate unique neighbors of node \$i\$
+
+# Details
+
+This is equivalent of taking the upper 
+diagonal of the adjacency matrix
+"""
 function unique_neighbors(ig::MetaGraph, i::Int)
     nbrs = neighbors(ig::MetaGraph, i::Int)
     filter(j -> j > i, nbrs)

src/search.jl

@@ -6,6 +6,7 @@ struct MPSControl
     max_bond::Int
     var_ϵ::Number
     max_sweeps::Int
+    β::Vector
 end
 
 # ρ needs to be in the right canonical form
@@ -72,57 +73,57 @@ function _apply_bias!(ψ::AbstractMPS, ig::MetaGraph, dβ::Number, i::Int)
     ψ[i] = M
 end
 
-function _apply_exponent!(ψ::AbstractMPS, ig::MetaGraph, dβ::Number, i::Int, j::Int)
+function _apply_exponent!(ψ::AbstractMPS, ig::MetaGraph, dβ::Number, i::Int, j::Int, last::Int)
     M = ψ[j]
-
-    d = size(M, 2)
-    basis = local_basis(d)
-
+    D = I(ψ, i)
+
     J = get_prop(ig, i, j, :J)  
-    C = [exp(0.5 * dβ * k * J * l) for k ∈ basis, l ∈ basis]
-    D = I(d)
+    C = [ exp(0.5 * dβ * k * J * l) for k ∈ local_basis(ψ, i), l ∈ local_basis(ψ, j) ]
 
-    if j == length(ψ)
-        @cast M̃[(x, a), σ, b] := C[σ, x] * M[a, σ, b]   
+    if j == last
+        @cast M̃[(x, a), σ, b] := C[x, σ] * M[a, σ, b]   
     else
-        @cast M̃[(x, a), σ, (y, b)] := C[σ, x] * D[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)
     M = ψ[i]
-    D = I(size(M, 2))
+    D = I(ψ, i)
 
     @cast M̃[a, σ, (y, b)] := D[σ, y] * M[a, σ, b]
     ψ[i] = M̃
 end
 
-function _apply_nothing!(ψ::AbstractMPS, i::Int) 
-    M = ψ[i] 
-    D = I(size(M, 2))
+function _apply_nothing!(ψ::AbstractMPS, l::Int, i::Int) 
+    M = ψ[l] 
+    D = I(ψ, i)
 
     @cast M̃[(x, a), σ, (y, b)] := D[x, y] * M[a, σ, b] 
-    ψ[i] = M̃    
+    ψ[l] = M̃    
 end
 
 _holes(nbrs::Vector) = setdiff(first(nbrs) : last(nbrs), nbrs)
 
-function MPS(ig::MetaGraph, mps::MPSControl, gibbs::GibbsControl)
+
+function MPS(ig::MetaGraph, control::MPSControl)
     L = nv(ig)
 
-    Dcut = mps.max_bond
-    tol = mps.var_ϵ
-    max_sweeps = mps.max_sweeps
+    Dcut = control.max_bond
+    tol = control.var_ϵ
+    max_sweeps = control.max_sweeps
+    schedule = control.β
     @info "Set control parameters for MPS" Dcut tol max_sweeps
 
-    β = gibbs.β
-    schedule = gibbs.β_schedule
+    β = get_prop(ig, :β)
+    rank = get_prop(ig, :rank)
 
     @assert β ≈ sum(schedule) "Incorrect β schedule."
 
     @info "Preparing Hadamard state as MPS"
-    ρ = HadamardMPS(L)
+    ρ = HadamardMPS(rank)
     is_right = true
 
     @info "Sweeping through β and σ" schedule
@@ -135,11 +136,11 @@ function MPS(ig::MetaGraph, mps::MPSControl, gibbs::GibbsControl)
             _apply_projector!(ρ, i)
 
             for j ∈ nbrs 
-                _apply_exponent!(ρ, ig, dβ, i, j) 
+                _apply_exponent!(ρ, ig, dβ, i, j, last(nbrs)) 
             end
 
             for l ∈ _holes(nbrs) 
-                _apply_nothing!(χ, l) 
+                _apply_nothing!(χ, l, i) 
             end
         end
 

src/utils.jl

@@ -2,12 +2,13 @@ export idx, ising, proj
 export HadamardMPS, rq
 export all_states, local_basis
 
-ising(σ::Union{Vector, NTuple}) = 2 .* σ .- 1
-
 idx(σ::Int) = (σ == -1) ? 1 : σ + 1
 _σ(idx::Int) = (idx == 1) ? -1 : idx - 1 
 
+LinearAlgebra.I(ψ::AbstractMPS, i::Int) = I(size(ψ[i], 2))
+
 local_basis(d::Int) = union(-1, 1:d-1)
+local_basis(ψ::AbstractMPS, i::Int) = local_basis(size(ψ[i], 2))
 
 function proj(state, dims::T) where {T <: Union{Vector, NTuple}}
     P = Matrix{Float64}[] 
@@ -25,9 +26,11 @@ function all_states(rank::T) where T <: Union{Vector, NTuple}
 end 
 
 function HadamardMPS(rank::T) where T <: Union{Vector, NTuple} 
-    MPS(fill(1, r) ./ sqrt(r) for r ∈ rank)
+    vec = [ fill(1, r) ./ sqrt(r) for r ∈ rank ]
+    MPS(vec)
 end
-HadamardMPS(L::Int) = MPS(fill([1., 1.] / sqrt(2), L))
+
+HadamardMPS(L::Int) = MPS(fill(2, L))
 
 function LinearAlgebra.qr(M::AbstractMatrix, Dcut::Int, args...) 
     fact = pqrfact(M, rank=Dcut, args...)

test/base.jl

@@ -23,6 +23,27 @@ T = ComplexF64
     @test ϕ ≈ ψ
 
     show(ψ)
+
+    dims = (3, 2, 5, 4) 
+    @info "Veryfing ψ of arbitrary rank" dims
+
+    ψ = randn(MPS{T}, D, dims)
+    @test verify_bonds(ψ) == nothing
+
+    @test ψ == ψ
+    @test ψ ≈ ψ
+
+    @test length(ψ) == length(dims)
+    @test size(ψ) == (length(dims), )
+    @test eltype(ψ) == ComplexF64
+    @test rank(ψ) == dims
+    @test bond_dimension(ψ) ≈ D
+
+    ϕ = copy(ψ) 
+    @test ϕ == ψ
+    @test ϕ ≈ ψ
+
+    show(ψ)
 end
 
 @testset "Random MPO" begin