Add files via upload

src/search_ad.jl

@@ -0,0 +1,150 @@
+#export MPS_from_gates, unique_neighbors, probable_states
+export MPS, unique_neighbors, probable_states
+export MPSControl
+
+struct MPSControl 
+    max_bond::Int
+    ϵ::Number
+    max_sweeps::Int
+end
+
+function unique_neighbors(ig::MetaGraph, i::Int)
+    nbrs = neighbors(ig::MetaGraph, i::Int)
+    filter(j -> j > i, nbrs)
+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 
+    ψ[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]]
+
+    if j == length(ψ)
+        @cast M̃[(x, a), σ, b] := C[x, σ] * δ[x, 1] * M[a, σ, b]  
+    else
+        @cast M̃[(x, a), σ, (y, b)] := C[x, σ] * δ[x, y] * M[a, σ, b]                     
+    end     
+    ψ[j] = M̃
+end
+
+function _apply_projector!(ψ::AbstractMPS, i::Int)
+    δ = I(2)
+    M = ψ[i]
+
+    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 
+    ψ[i] = M̃
+end
+
+function _apply_nothing!(ψ::AbstractMPS, i::Int)    
+    δ = I(2)       
+    M = ψ[i] 
+
+    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
+    ψ[i] = M̃    
+end
+
+function MPS(ig::MetaGraph, mps::MPSControl, gibbs::GibbsControl)
+    L = nv(ig)
+
+    # control for MPS
+    Dcut = mps.max_bond
+    tol = mps.ϵ
+    max_sweeps = mps.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)
+
+    for dβ ∈ schedule, i ∈ 1:L
+        _apply_bias!(ρ, ig, dβ, i) 
+
+        nbrs = unique_neighbors(ig, i)
+        if !isempty(nbrs)
+            for j ∈ nbrs 
+                _apply_exponent!(ρ, ig, dβ, i, j) 
+            end
+
+            _apply_projector!(ρ, i)
+
+            for l ∈ setdiff(1:L, union(i, nbrs)) 
+                _apply_nothing!(ρ, l) 
+            end
+        end
+
+        # reduce bond dimension
+        if bond_dimension(ρ) > Dcut
+            ρ = compress(ρ, Dcut, tol, max_sweeps) 
+        end
+    end
+    ρ
+end
+
+#function search_MPS(ψ::AbstractMPS, J::AdjacencyMatrix, N::Int, 
+#    MaxStates::Int, Cutoff::Int, β::Int, Iterations::Int,  ig::MetaGraph, 
+#    mps::MPSControl, gibbs::GibbsControl)
+#    """Search low energy states of classical Ising model using matrix product states"""
+
+#    for iii ∈ range(1, stop = Iterations)
+#        println(iii, '/', Iterations)
+#        ψ = MPS(ig, mps, gibbs)
+#    end
+#    prob, lprob, states = probable_states(ψ, N, Max_States)
+#    return prob, lprob, states
+#end
+
+function probable_states(ψ::Vector{Array{T, 3}}, N::Int64, Max_States::Int64) where T <: AbstractFloat
+    partial_conf = Array{Int64}(undef, 1, 0)
+    new_prob = zeros(0)
+    #for n in 1:N 
+        #new_conf = (p=(a,c)->c<2 ? a : [[x;y] for x=a for y=p(a,c-1)])(1:2,n)
+        new_conf = (p=(a,c)->c<2 ? a : [[x;y] for x=a for y=p(a,c-1)])(1:2,N)
+        #new_conf = reverse.(Iterators.product(fill(1:2,N)...))[:]
+        m=size(new_conf)[1]
+        for k in 1:m
+            #np = compute_probs(ψ[n], new_conf[k])
+            np = compute_probs(ψ, new_conf[k])
+            append!(new_prob, np[1])
+        end
+        idx = sortperm(new_prob; rev=true)
+        if length(idx) > Max_States
+            idx = idx[1:Max_States]
+        end
+        prob = new_prob[idx]
+        partial_conf = new_conf[idx]
+        println(prob)
+        println(partial_conf)
+
+    #end
+    sigma = partial_conf[1]
+    #println(sigma)
+    #E = energy(sigma, ig)
+    #println(E)
+    return prob, partial_conf
+end 

src/search_test.jl

@@ -0,0 +1,95 @@
+using CUDA
+using SpinGlassPEPS
+using LinearAlgebra
+using TensorOperations
+using MetaGraphs
+using LightGraphs
+using GraphPlot
+
+
+L = 3
+N = L^2
+instance = "./lattice_$L.txt"    
+ig = ising_graph(instance, N)
+Dcut = 16
+var_tol=1E-8
+max_sweeps = 4
+
+β = 1
+dβ = 0.25
+β_schedule = [dβ for _ ∈ 1:4]
+
+gibbs_param = GibbsControl(β, β_schedule)
+mps_param = MPSControl(Dcut, var_tol, max_sweeps) 
+
+#ρ = MPS(ig, mps_param, gibbs_param)
+#show(ρ)
+#println(size(ρ))
+
+using TensorOperations
+include("notation.jl")
+function make_qubo()
+    qubo = [(1,1) .5; (1,2) -0.5; (1,4) -1.5; (2,2) -1.; (2,3) -1.5; (2,5) -0.5; (3,3) 2.; (3,6) 1.5]
+    qubo = vcat(qubo, [(6,6) .05; (5,6) -0.25; (6,9) -0.52; (5,5) 0.75; (4,5) 0.5; (5,8) 0.5; (4,4) 0.; (4,7) -0.01])
+    qubo = vcat(qubo, [(7,7) 0.35; (7,8) 0.5; (8,8) -0.08; (8,9) -0.05; (9,9) 0.33])
+    [Qubo_el(qubo[i,1], qubo[i,2]) for i in 1:size(qubo, 1)]
+end
+
+qubo =  make_qubo()
+β = 2.
+include("mps_implementation.jl")
+mps = initialize_mps(L)
+
+grid = [1 2 3; 4 5 6; 7 8 9]
+ns = [Node_of_grid(i,grid) for i in 1:maximum(grid)]
+
+is, js = connections_for_mps(ns)
+all_is, all_js = cluster_conncetions(is,js)
+
+mps = construct_mps(qubo, β, 2, ns, all_is, all_js, 4, 0.)
+prob=compute_probs(mps, [1,2,1])
+#println(prob)
+
+a = [1, 1, 1]
+A = fill!(Vector{Vector{Int}}(undef, 1), a)
+new_conf = (p=(a,c)->c<2 ? a : [[x;y] for x=a for y=p(a,c-1)])(1:2,L)
+aa = new_conf[3][1]
+bb = new_conf[3][2]
+cc = new_conf[3][3]
+
+#prob=compute_probs(mps, [aa,bb,cc])
+#println(prob)
+
+
+
+ function probable_states(ψ::Vector{Array{T, 3}}, N::Int64, Max_States::Int64) where T <: AbstractFloat
+    partial_conf = Array{Int64}(undef, 1, 0)
+    new_prob = zeros(0)
+    #for n in 1:N 
+        #new_conf = (p=(a,c)->c<2 ? a : [[x;y] for x=a for y=p(a,c-1)])(1:2,n)
+        new_conf = (p=(a,c)->c<2 ? a : [[x;y] for x=a for y=p(a,c-1)])(1:2,N)
+        #new_conf = reverse.(Iterators.product(fill(1:2,N)...))[:]
+        m=size(new_conf)[1]
+        for k in 1:m
+            #np = compute_probs(ψ[n], new_conf[k])
+            np = compute_probs(ψ, new_conf[k])
+            append!(new_prob, np[1])
+        end
+        idx = sortperm(new_prob; rev=true)
+        if length(idx) > Max_States
+            idx = idx[1:Max_States]
+        end
+        prob = new_prob[idx]
+        partial_conf = new_conf[idx]
+        println(prob)
+        println(partial_conf)
+
+    #end
+    sigma = partial_conf[1]
+    #println(sigma)
+    #E = energy(sigma, ig)
+    #println(E)
+    return prob, partial_conf
+end 
+
+p = probable_states(mps, L, 5)