From d04a9b7aaaf7d6be8bc0e0676e226b2f70f9c418 Mon Sep 17 00:00:00 2001 From: annamariadziubyna <73058800+annamariadziubyna@users.noreply.github.com> Date: Tue, 17 Nov 2020 21:36:46 +0100 Subject: [PATCH] Add files via upload --- src/search_ad.jl | 150 +++++++++++++++++++++++++++++++++++++++++++++ src/search_test.jl | 95 ++++++++++++++++++++++++++++ 2 files changed, 245 insertions(+) create mode 100644 src/search_ad.jl create mode 100644 src/search_test.jl diff --git a/src/search_ad.jl b/src/search_ad.jl new file mode 100644 index 00000000..0defdf3f --- /dev/null +++ b/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 \ No newline at end of file diff --git a/src/search_test.jl b/src/search_test.jl new file mode 100644 index 00000000..616dd6b0 --- /dev/null +++ b/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) \ No newline at end of file