ad/mps (#30)

* splitting MPS function into MPS and _compress_MPS

* Renaming _compress_MPS to _apply_layer_of_gates, adding multiply_purifications

* MPS2 does not work

* some changes in MPS2 - working

* add empty MPS into multiply_purifications

* start splitting test files, corrections in MPS2

* correct holes in _apply_layer_of_gates, correct MPS2, add some tests

* unification of notation

* changes after pulling request

* unification of MPS and MPS2, change type from String to Symbol, add some tests

* some changes in tests

* Update Project.toml

* Update Project.toml

* remove Infiltrator

* corrections in full_spectrum

* add more tests

src/spectrum.jl

@@ -2,6 +2,7 @@ export unique_neighbors
 export full_spectrum, brute_force
 export MPSControl
 export solve, solve_new
+export MPS2
 
 struct Spectrum
     energies::Array{<:Number}
@@ -12,7 +13,8 @@ struct MPSControl
     max_bond::Int
     var_ϵ::Number
     max_sweeps::Int
-    β::Vector
+    β::Vector 
+    dβ::Vector
 end
 
 # ρ needs to be ∈ the right canonical form
@@ -170,30 +172,32 @@ function _apply_nothing!(ψ::AbstractMPS, l::Int, i::Int)
     ψ[l] = M̃
 end
 
-_holes(nbrs::Vector, i::Int) = setdiff(i + 1 : last(nbrs), nbrs)
 
-function MPS(ig::MetaGraph, control::MPSControl, β::Float64=1.0)
-    L = nv(ig)
-
-    Dcut = control.max_bond
-    tol = control.var_ϵ
-    max_sweeps = control.max_sweeps
-    schedule = control.β
-    @info "Set control parameters for MPS" Dcut tol max_sweeps
+function multiply_purifications(χ::AbstractMPS, ϕ::AbstractMPS, L::Int)
+    T = eltype(χ)
+    ψ = MPS(T, L)
 
-    rank = get_prop(ig, :rank)
+    for i ∈ 1:L 
+        A1 = χ[i]
+        A2 = ϕ[i]
+
+        @cast B[(l, x), σ, (r, y)] := A1[l, σ, r] * A2[x, σ, y]
+        ψ[i] = B
+    end
+    ψ
 
-    @assert β ≈ sum(schedule) "Incorrect β schedule."
+end
 
-    @info "Preparing Hadamard state as MPS"
-    ρ = HadamardMPS(rank)
-    is_right = true
+_holes(l::Int, nbrs::Vector) = setdiff(l+1 : last(nbrs), nbrs)
 
-    @info "Sweeping through β and σ" schedule
-    for dβ ∈ schedule, i ∈ 1:L
-        _apply_bias!(ρ, ig, dβ, i)
+function _apply_layer_of_gates(ig::MetaGraph, ρ::AbstractMPS, control::MPSControl, dβ::Number)
+    L = nv(ig)
+    Dcut = control.max_bond
+    tol = control.var_ϵ
+    max_sweeps = control.max_sweeps
+    for i ∈ 1:L
+        _apply_bias!(ρ, ig, dβ, i) 
         is_right = false
-
         nbrs = unique_neighbors(ig, i)
         if !isempty(nbrs)
             _apply_projector!(ρ, i)
@@ -202,8 +206,8 @@ function MPS(ig::MetaGraph, control::MPSControl, β::Float64=1.0)
                 _apply_exponent!(ρ, ig, dβ, i, j, last(nbrs))
             end
 
-            for l ∈ _holes(nbrs, i)
-                _apply_nothing!(ρ, l, i) 
+            for l ∈ _holes(i, nbrs) 
+                _apply_nothing!(ρ, l, i)  
             end
         end
 
@@ -212,15 +216,83 @@ function MPS(ig::MetaGraph, control::MPSControl, β::Float64=1.0)
             ρ = compress(ρ, Dcut, tol, max_sweeps)
             is_right = true
         end
-    end
-
+        
+    end 
     if !is_right
         canonise!(ρ, :right)
         is_right = true
     end
     ρ
 end
 
+function MPS(ig::MetaGraph, control::MPSControl)
+
+    Dcut = control.max_bond
+    tol = control.var_ϵ
+    max_sweeps = control.max_sweeps
+    schedule = control.β
+    @info "Set control parameters for MPS" Dcut tol max_sweeps
+
+    β = get_prop(ig, :β)
+    rank = get_prop(ig, :rank)
+
+    @assert β ≈ sum(schedule) "Incorrect β schedule."
+
+    @info "Preparing Hadamard state as MPS"
+    ρ = HadamardMPS(rank)
+    is_right = true
+    @info "Sweeping through β and σ" schedule
+    for dβ ∈ schedule
+        ρ = _apply_layer_of_gates(ig, ρ, control, dβ)
+    end
+    ρ
+end
+
+function MPS(ig::MetaGraph, control::MPSControl, type::Symbol) 
+    L = nv(ig)
+    Dcut = control.max_bond
+    tol = control.var_ϵ
+    max_sweeps = control.max_sweeps
+    @info "Set control parameters for MPS" Dcut tol max_sweeps
+    dβ = get_prop(ig, :dβ)
+    β = get_prop(ig, :β)
+    rank = get_prop(ig, :rank)
+
+    @info "Preparing Hadamard state as MPS"
+    ρ = HadamardMPS(rank)
+    is_right = true
+    @info "Sweeping through β and σ" dβ
+
+    if type == :log
+        k = ceil(log2(β/dβ))
+        dβmax = β/(2^k)
+        ρ = _apply_layer_of_gates(ig, ρ, control, dβmax)
+        for j ∈ 1:k
+            ρ = multiply_purifications(ρ, ρ, L)
+            if bond_dimension(ρ) > Dcut
+                @info "Compresing MPS" bond_dimension(ρ), Dcut
+                ρ = compress(ρ, Dcut, tol, max_sweeps) 
+                is_right = true
+            end
+        end
+        ρ
+    elseif type == :lin
+        k = β/dβ
+        dβmax = β/k
+        ρ = _apply_layer_of_gates(ig, ρ, control, dβmax)
+        ρ0 = copy(ρ)
+        for j ∈ 1:k
+            ρ = multiply_purifications(ρ, ρ0, L)
+            if bond_dimension(ρ) > Dcut
+                @info "Compresing MPS" bond_dimension(ρ), Dcut
+                ρ = compress(ρ, Dcut, tol, max_sweeps) 
+                is_right = true
+            end
+        end
+    end
+    ρ
+
+end
 
 """
 $(TYPEDSIGNATURES)
@@ -255,4 +327,4 @@ function full_spectrum(cl::Cluster; num_states::Int=prod(cl.rank))
     σ = collect.(all_states(cl.rank))
     energies = energy.(σ, Ref(cl))
     Spectrum(energies[1:num_states], σ[1:num_states])   
-end
+end

test/runtests.jl

@@ -37,6 +37,7 @@ push!(my_tests,
     "compressions.jl",
     "ising.jl",
     "indexing.jl",
+    "searchMPS.jl",
     "spectrum.jl",
     "graph.jl",
     "PEPS.jl",

test/searchMPS.jl

@@ -0,0 +1,112 @@
+using MetaGraphs
+using LightGraphs
+using GraphPlot
+
+L = 4
+N = L^2
+
+instance = "$(@__DIR__)/instances/$(N)_001.txt"  
+
+ig = ising_graph(instance, N)
+set_prop!(ig, :β, 1.) #rand(Float64))
+r = fill(2, N)
+set_prop!(ig, :rank, r)
+set_prop!(ig, :dβ, 0.001)
+
+ϵ = 1E-5
+D = 16
+var_ϵ = 1E-8
+sweeps = 40
+type1 = :log
+type2 = :lin
+β = [get_prop(ig, :β)]
+dβ = [get_prop(ig, :dβ)] 
+control = MPSControl(D, var_ϵ, sweeps, β, dβ) 
+
+states = all_states(get_prop(ig, :rank))
+ϱ = gibbs_tensor(ig)
+@test sum(ϱ) ≈ 1
+
+@testset "Verifying gate operations 2" begin
+    rank = get_prop(ig, :rank)
+    @info "Testing MPS"
+
+    rψ1 = MPS(ig, control, type1)
+    rψ2 = MPS(ig, control, type2)
+    rψ3 = MPS(ig, control)
+    overlap12 = dot(rψ1, rψ2)
+    @test abs(1 - overlap12) < ϵ
+    overlap13 = dot(rψ1, rψ3)
+    @test abs(1 - overlap13) < ϵ
+    overlap23 = dot(rψ2, rψ3)
+    @test abs(1 - overlap23) < ϵ
+
+    for max_states ∈ [1, N, 2*N, N^2]
+        @info "Testing spectrum_new"
+        states_new1, prob_new1, pCut_new1 = solve_new(rψ1, max_states)   
+        states_new2, prob_new2, pCut_new2 = solve_new(rψ2, max_states)
+        states_new3, prob_new3, pCut_new3 = solve_new(rψ3, max_states)         
+
+        for st ∈ states_new1
+            @test st ∈ states_new2
+        end
+        eng_new1 = zeros(length(prob_new1))
+        eng_new2 = zeros(length(prob_new2))
+        eng_new3 = zeros(length(prob_new3))
+        for (j, p) ∈ enumerate(prob_new1)
+            σ = states_new1[j, :]
+            eng_new1[j] = energy(σ, ig)
+        end
+        for (j, p) ∈ enumerate(prob_new2)
+            σ = states_new2[j, :]
+            eng_new2[j] = energy(σ, ig)
+        end
+        for (j, p) ∈ enumerate(prob_new3)
+            σ = states_new3[j, :]
+            eng_new3[j] = energy(σ, ig)
+        end
+
+        perm1 = partialsortperm(eng_new1, 1:max_states)
+        eng_new1 = eng_new1[perm1]
+        states_new1 = states_new1[perm1, :]
+        prob_new1 = prob_new1[perm1]
+        state1 = states_new1[1, :]
+        @info "Testing MPS2 - logarithmic"
+        @info "The largest discarded probability" pCut_new1
+        @info "State with the lowest energy" state1
+        @info "Probability of the state with the lowest energy" prob_new1[1]
+        @info "The lowest energy" eng_new1[1]
+
+        perm2 = partialsortperm(eng_new2, 1:max_states)
+        eng_new2 = eng_new2[perm2]
+        states_new2 = states_new2[perm2, :]
+        prob_new2 = prob_new2[perm2]
+        state2 = states_new2[1, :]
+        @info "Testing MPS2 - linear"
+        @info "The largest discarded probability" pCut_new2
+        @info "State with the lowest energy" state2
+        @info "Probability of the state with the lowest energy" prob_new2[1]
+        @info "The lowest energy" eng_new2[1]
+
+        perm3 = partialsortperm(eng_new3, 1:max_states)
+        eng_new3 = eng_new3[perm3]
+        states_new3 = states_new3[perm3, :]
+        prob_new3 = prob_new3[perm3]
+        state3 = states_new3[1, :]
+        @info "Testing MPS"
+        @info "The largest discarded probability" pCut_new3
+        @info "State with the lowest energy" state3
+        @info "Probability of the state with the lowest energy" prob_new3[1]
+        @info "The lowest energy" eng_new3[1]
+
+        if state1 == state2
+            @test eng_new1[1] == eng_new2[1]
+        elseif state2 == state3
+            @test eng_new2[1] == eng_new3[1]
+        elseif state1 == state3
+            @test eng_new1[1] == eng_new3[1]
+        end
+
+
+    end
+end

test/spectrum.jl

@@ -2,30 +2,33 @@ using MetaGraphs
 using LightGraphs
 using GraphPlot
 
-L = 4
+L = 2
 N = L^2
 
 instance = "$(@__DIR__)/instances/$(N)_001.txt"  
 
 ig = ising_graph(instance, N)
-β = 5.
+set_prop!(ig, :β, 1.) #rand(Float64))
 #r = [3, 2, 5, 4]
 r = fill(2, N)
 set_prop!(ig, :rank, r)
+set_prop!(ig, :dβ, 0.01)
 
 sgn = -1.
-ϵ = 1E-7
+ϵ = 1E-8
 D = prod(r) + 1
 var_ϵ = 1E-8
 sweeps = 4
-schedule = [β]
-control = MPSControl(D, var_ϵ, sweeps, schedule) 
+schedule = [get_prop(ig, :β)]
+dβ = [get_prop(ig, :dβ)]
+control = MPSControl(D, var_ϵ, sweeps, schedule, dβ) 
 
 states = all_states(get_prop(ig, :rank))
-ϱ = gibbs_tensor(ig, β)
+ϱ = gibbs_tensor(ig)
 @test sum(ϱ) ≈ 1
 
 @testset "Verifying gate operations" begin
+    β = get_prop(ig, :β)
     rank = get_prop(ig, :rank)
 
     χ = HadamardMPS(rank)
@@ -56,7 +59,7 @@ states = all_states(get_prop(ig, :rank))
                 end
             end
 
-            for l ∈ SpinGlassPEPS._holes(nbrs, i) 
+            for l ∈ SpinGlassPEPS._holes(i, nbrs)
                 SpinGlassPEPS._apply_nothing!(χ, l, i) 
             end
         end
@@ -75,6 +78,7 @@ end
 
     @testset "Exact Gibbs pure state (MPS)" begin
         L = nv(ig)
+        β = get_prop(ig, :β)
         rank = get_prop(ig, :rank)
 
         @info "Generating Gibbs state - |ρ>" L rank β ϵ
@@ -117,7 +121,7 @@ end
 
         @info "Verifying MPS from gates"
 
-        Gψ = MPS(ig, control, β) 
+        Gψ = MPS(ig, control) 
 
         @test_nowarn is_right_normalized(Gψ)
         @test bond_dimension(Gψ) > 1
@@ -136,16 +140,14 @@ end
             @test ϱ[idx.(σ)...] ≈ p
         end 
 
-        for max_states ∈ [N, 2*N, N^2]
+        for max_states ∈ [1, N, 2*N, N^2]
 
             @info "Verifying low energy spectrum" max_states
             @info "Testing spectrum"
             states, prob, pCut = solve(rψ, max_states)
             sp = brute_force(ig, num_states = max_states)
 
             @info "The largest discarded probability" pCut
-            @test maximum(prob) > pCut
-
 
             for (j, (p, e)) ∈ enumerate(zip(prob, sp.energies))
                 σ = states[:, j]
@@ -174,5 +176,6 @@ end
             @info "The lowest energy" eng_new[1]
 
         end
-    end    
+
+    end
 end