add Gibbs_tensor

src/ising.jl

@@ -1,4 +1,19 @@
 export ising_graph, energy
+export Gibbs_tensor
+
+function Gibbs_tensor(ig::MetaGraph, opts::Gibbs_control)
+    L = nv(ig)
+    β = opts.β
+
+    all_states = Iterators.product([[-1, 1] for _ ∈ 1:L]...)
+    rank = [2 for i ∈ 1:L]
+
+    r = [exp(-β * energy(ig, collect(σ))) for σ ∈ all_states]
+    ρ = reshape(r, rank...)
+
+    Z = sum(ρ)
+    return ρ / Z
+end
 
 function energy(ig::MetaGraph, σ::Vector{<:Number})
     """

test/search.jl

@@ -15,16 +15,23 @@ ig = ising_graph(instance, N)
 
     β = 1
     dβ = 0.25
-    β_schedule = [dβ, dβ, dβ, dβ]
+    β_schedule = [dβ for _ ∈ 1:4]
 
-    mps = MPS_control(Dcut, var_tol, max_sweeps) 
-    gibbs = Gibbs_control(β, β_schedule)
+    gibbs_param = Gibbs_control(β, β_schedule)
+    mps_param = MPS_control(Dcut, var_tol, max_sweeps) 
 
-    ρ = MPS_from_gates(ig, mps, gibbs) 
+    @testset "Generate ρ" begin
+        ρ = MPS_from_gates(ig, mps_param, gibbs_param) 
 
-    show(ρ)
-    @test _verify_bonds(ρ)
+        show(ρ)
+        @test _verify_bonds(ρ)
 
-    canonise!(ρ)
-    @test dot(ρ, ρ) ≈ 1
+        canonise!(ρ)
+        @test dot(ρ, ρ) ≈ 1
+    end
+
+    @testset "Generate Gibbs state exactly" begin
+        r = Gibbs_tensor(ig, gibbs_param)
+        @test sum(r) ≈ 1
+    end
 end