diff --git a/.gitignore b/.gitignore index 0887050..2549c28 100644 --- a/.gitignore +++ b/.gitignore @@ -2,4 +2,4 @@ *.jl.cov *.jl.mem /docs/Manifest.toml -/docs/build/ +/docs/build/ \ No newline at end of file diff --git a/docs/src/index.md b/docs/src/index.md index d5cbd6d..c995c8b 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -22,6 +22,7 @@ Modules = [ ADAPT.Basics, ADAPT.Basics.Callbacks, ADAPT.Basics.Operators, + ADAPT.Basics.Pools, ] ``` @@ -30,6 +31,8 @@ Modules = [ Modules = [ ADAPT.OptimizationFreeADAPT, ADAPT.OverlapADAPT, + ADAPT.Degenerate_ADAPT, + ADAPT.Hamiltonians, ] ``` diff --git a/src/ADAPT.jl b/src/ADAPT.jl index 5922713..29b2653 100644 --- a/src/ADAPT.jl +++ b/src/ADAPT.jl @@ -76,6 +76,10 @@ module ADAPT # A suite of common operators, especially useful for constructing operator pools. include("base/Operators.jl") export Operators + + # A suite of common operator pools. + include("base/pools.jl") + export Pools end using .Basics export Ansatz @@ -83,6 +87,7 @@ module ADAPT export OptimOptimizer export Callbacks export Operators + export Pools module OptimizationFreeADAPT include("optimizationfree/OptimizationFree.jl") @@ -94,12 +99,21 @@ module ADAPT include("overlap/OverlapADAPT.jl") export Infidelity end + + module Degenerate_ADAPT + include("degenerateADAPT/DegenerateADAPT.jl") + export DEG_ADAPT + end + + module Hamiltonians + # A suite of common lattice Hamiltonians. + include("hamiltonians/latticemodels.jl") + end end #= TODO: -- Some pools. Grab from Diksha's code. - Lie rank calculation, in pauli-dedicated code. =# \ No newline at end of file diff --git a/src/base/__paulioperators.jl b/src/base/__paulioperators.jl index 1b2890c..e907e69 100644 --- a/src/base/__paulioperators.jl +++ b/src/base/__paulioperators.jl @@ -264,11 +264,29 @@ function measure_commutator( return expectation_value(commutator, Ψ) end +function Base.:≈(spv1::Vector{ScaledPauli{N}}, spv2::Vector{ScaledPauli{N}}) where {N} + if length(spv1) != length(spv2) return false end + for (sp1,sp2) in zip(spv1, spv2) + if !(sp1≈sp2) + return false + end + end + return true +end + +function otimes(sp::ScaledPauli{N}, p::Pauli{M}) where {N,M} + out_pauli = sp.pauli ⊗ p.pauli + out_coeff = sp.coeff * get_phase(p) + out = ScaledPauli{N+M}(out_coeff, out_pauli) + return out +end - - - - +function otimes(p::Pauli{M}, sp::ScaledPauli{N}) where {M,N} + out_pauli = p.pauli ⊗ sp.pauli + out_coeff = sp.coeff * get_phase(p) + out = ScaledPauli{N+M}(out_coeff, out_pauli) + return out +end #= TODO: Calculate DLA diff --git a/src/base/pools.jl b/src/base/pools.jl new file mode 100644 index 0000000..0080759 --- /dev/null +++ b/src/base/pools.jl @@ -0,0 +1,363 @@ +#= Set up various operator pools. =# +#= TO DO: do the extra functions in MyPauliOperators need to be added to the PauliOperators.jl package? =# +#= TO DO: tile operators by removing leading and trailing I's first. =# + +module Pools + + import PauliOperators + import PauliOperators: FixedPhasePauli, Pauli, ScaledPauli + import PauliOperators: ScaledPauliVector, PauliSum + + import ..MyPauliOperators: otimes, ≈ + import PauliOperators: ⊗, ≈ + + """ + qubitexcitation(n::Int, i::Int, k::Int) + qubitexcitation(n::Int, i::Int, j::Int, k::Int, l::Int) + + Qubit excitation operators as defined in Yordanov et al. 2021. + + Note that Yordanov's unitaries are defined as `exp(iθG)` rather than `exp(-iθG)`, + so variational parameters will be off by a sign. + + # Parameters + - `n`: total number of qubits + - `i,j,k,l`: qubit indices as defined in Yordanov's paper. + + # Returns + - `PauliOperators.ScaledPauliVector`: the qubit excitation operator + + Note that all Pauli terms in any single qubit excitation operator commute, + so the `ScaledPauliVector` representation is "safe". + + """ + function qubitexcitation(n::Int, i::Int, k::Int) + return (1/2) .* [ + ScaledPauli(Pauli(n; X=[i], Y=[k])), + -ScaledPauli(Pauli(n; X=[k], Y=[i])), + ] + end + + function qubitexcitation(n::Int, i::Int, j::Int, k::Int, l::Int) + return (1/8) .* [ + ScaledPauli(Pauli(n; X=[i,k,l], Y=[j])), + ScaledPauli(Pauli(n; X=[j,k,l], Y=[i])), + ScaledPauli(Pauli(n; X=[l], Y=[i,j,k])), + ScaledPauli(Pauli(n; X=[k], Y=[i,j,l])), + -ScaledPauli(Pauli(n; X=[i,j,l], Y=[k])), + -ScaledPauli(Pauli(n; X=[i,j,k], Y=[l])), + -ScaledPauli(Pauli(n; X=[j], Y=[i,k,l])), + -ScaledPauli(Pauli(n; X=[i], Y=[j,k,l])), + ] + end + + """ + qubitexcitationpool(n_system::Int) + + The number of singles excitations = (n 2), and the doubles = 3*(n 4). + + # Parameters + - `n_system`: Number of qubits in the system + + # Returns + - `pool`: the qubit-excitation-based pool as defined in Communications Physics 4, 1 (2021). + - `target_and_source`: Dict mapping each pool operator to the target and source orbitals involved in the excitation. + """ + function qubitexcitationpool(n_system::Int) + pool = ScaledPauliVector{n_system}[] + target_and_source = Dict{ScaledPauliVector{n_system}, Vector{Vector{Int64}}}() + + for i in 1:n_system + for j in i+1:n_system + # singles excitations + op = qubitexcitation(n_system, i, j) + push!(pool, op) + target_and_source[op] = [[i,j]] + + # doubles excitations + for k in j+1:n_system + for l in k+1:n_system + target_pair = [i,j]; source_pair = [k,l] + new_op = qubitexcitation(n_system, target_pair[1], target_pair[2], source_pair[1], source_pair[2]) + push!(pool, new_op) + target_and_source[new_op] = [target_pair,source_pair] + + target_pair = [i,k]; source_pair = [j,l] + new_op = qubitexcitation(n_system, target_pair[1], target_pair[2], source_pair[1], source_pair[2]) + push!(pool, new_op) + target_and_source[new_op] = [target_pair,source_pair] + + target_pair = [j,k]; source_pair = [i,l] + new_op = qubitexcitation(n_system, target_pair[1], target_pair[2], source_pair[1], source_pair[2]) + push!(pool, new_op) + target_and_source[new_op] = [target_pair,source_pair] + end + end + end + end + return pool, target_and_source + end + + + function qubitexcitation_complemented(n::Int, i::Int, k::Int) + return (1/2) .* [ + ScaledPauli(Pauli(n; X=[i,k])), + ScaledPauli(Pauli(n; Y=[i,k])), + ] + end + function qubitexcitation_complemented(n::Int, i::Int, j::Int, k::Int, l::Int) + return (1/8) .* [ + -ScaledPauli(Pauli(n; X=[i,j,k,l])), + ScaledPauli(Pauli(n; X=[i,j], Y=[k,l])), + -ScaledPauli(Pauli(n; X=[i,k], Y=[j,l])), + -ScaledPauli(Pauli(n; X=[i,l], Y=[j,k])), + -ScaledPauli(Pauli(n; X=[j,k], Y=[i,l])), + -ScaledPauli(Pauli(n; X=[j,l], Y=[i,k])), + ScaledPauli(Pauli(n; X=[k,l], Y=[i,j])), + -ScaledPauli(Pauli(n; Y=[i,j,k,l])), + ] + end + + """ + qubitexcitationpool_complemented(n_system::Int) + + Returns the complemented qubit excitation pool on n_system qubits, inspired from arXiv 2109.01318. + + # Parameters + - `n_system`: Number of qubits in the system + + # Returns + - `pool`: the complemented qubit-excitation-based pool. + - `target_and_source`: Dict mapping each pool operator to the target and source orbitals involved in the excitation. + """ + function qubitexcitationpool_complemented(n_system::Int) + pool = ScaledPauliVector{n_system}[] + target_and_source = Dict{ScaledPauliVector{n_system}, Vector{Vector{Int64}}}() + + for i in 1:n_system + for j in i+1:n_system + # singles excitations + op = qubitexcitation_complemented(n_system, i, j) + push!(pool, op) + target_and_source[op] = [[i,j]] + + # doubles excitations + for k in j+1:n_system + for l in k+1:n_system + target_pair = [i,j]; source_pair = [k,l] + new_op = qubitexcitation_complemented(n_system, target_pair[1], target_pair[2], source_pair[1], source_pair[2]) + push!(pool, new_op) + target_and_source[new_op] = [target_pair,source_pair] + + target_pair = [i,k]; source_pair = [j,l] + new_op = qubitexcitation_complemented(n_system, target_pair[1], target_pair[2], source_pair[1], source_pair[2]) + push!(pool, new_op) + target_and_source[new_op] = [target_pair,source_pair] + + target_pair = [j,k]; source_pair = [i,l] + new_op = qubitexcitation_complemented(n_system, target_pair[1], target_pair[2], source_pair[1], source_pair[2]) + push!(pool, new_op) + target_and_source[new_op] = [target_pair,source_pair] + end + end + end + end + return pool, target_and_source + end + + #= TODO: Fix definition of qubit ADAPT pool so duplicates don't appear. =# + """ + qubitadaptpool(n_system::Int) + + Returns the qubit ADAPT pool on n_system qubits as defined in PRX QUANTUM 2, 020310 (2021). + It is generated by taking each qubit-excitation-based operator and breaking it into individual Pauli terms. + + # Parameters + - `n_system`: Number of qubits in the system + + # Returns + - `pool`: the qubit-ADAPT pool. + """ + function qubitadaptpool(n_system::Int) + pool = ScaledPauliVector{n_system}[] + + for i in 1:n_system + for j in i+1:n_system + # singles operators + excitation_op = qubitexcitation(n_system, i, j) + for sp in excitation_op + sp_new = ScaledPauli(1.0+0.0im, sp.pauli) + push!(pool, [sp_new]) + end + + # doubles operators + for k in j+1:n_system + for l in k+1:n_system + target_pair = [i,j]; source_pair = [k,l] + new_op = qubitexcitation(n_system, target_pair[1], target_pair[2], source_pair[1], source_pair[2]) + for sp in new_op + sp_new = ScaledPauli(1.0+0.0im, sp.pauli) + push!(pool, [sp_new]) + end + + target_pair = [i,k]; source_pair = [j,l] + new_op = qubitexcitation(n_system, target_pair[1], target_pair[2], source_pair[1], source_pair[2]) + for sp in new_op + sp_new = ScaledPauli(1.0+0.0im, sp.pauli) + push!(pool, [sp_new]) + end + + target_pair = [j,k]; source_pair = [i,l] + new_op = qubitexcitation(n_system, target_pair[1], target_pair[2], source_pair[1], source_pair[2]) + for sp in new_op + sp_new = ScaledPauli(1.0+0.0im, sp.pauli) + push!(pool, [sp_new]) + end + end + end + end + end + return pool + end + + """ + fullpauli(n::Int) + + The pool of all (4^n) n-qubit Pauli operators. + + # Parameters + - `n`: Number of qubits in the system + + # Returns + - `pool`: the full pauli pool. + """ + function fullpauli(n::Int) + pool = ScaledPauliVector{n}[] + for plist in Iterators.product(ntuple(i->["I","X","Y","Z"],n)...) + pstr = join(plist) + pauli = [ScaledPauli(Pauli(pstr))] + push!(pool, pauli) + end + return pool + end + + """ + one_local_pool(n::Int64, axes=["I","X","Y","Z"]) + + Returns the one-local pool containing each one-local operator on n qubits. + + # Parameters + - `n`: Number of qubits in the system + + # Returns + - `pool`: the one-local pool. + """ + function one_local_pool(n::Int64, axes=["I","X","Y","Z"]) + pool = ScaledPauliVector(n) + for i in 1:n + "X" in axes && (push!(pool, ScaledPauli(Pauli(n; X=i)))) + "Y" in axes && (push!(pool, ScaledPauli(Pauli(n; Y=i)))) + "Z" in axes && (push!(pool, ScaledPauli(Pauli(n; Z=i)))) + end + return pool + end + + """ + two_local_pool(n::Int64, axes=["X","Y","Z"]) + + Returns the two-local pool containing each two-local operator on n qubits. + + # Parameters + - `n`: Number of qubits in the system + + # Returns + - `pool`: the two-local pool. + """ + function two_local_pool(n::Int64, axes=["X","Y","Z"]) + pool = ScaledPauliVector{n}[] + for pair in Iterators.product(ntuple(i->1:n, 2)...) + i,j = pair + if i < j + for pair2 in Iterators.product(ntuple(i->axes, 2)...) + a,b = pair2 + if a == "I" || b == "I" + continue + end + l = "I"^(i-1)*a*("I"^(j-i-1))*b*"I"^(n-j) + pauli = [ScaledPauli(Pauli(l))] + push!(pool, pauli) + end + end + end + return pool + end + + """ + oneandtwo_local_pool(n::Int64) + + Returns the union of the one-local and two-local pools on n qubits. + + # Parameters + - `n`: Number of qubits in the system + + # Returns + - `pool`: union of one-local and two-local pools. + """ + function oneandtwo_local_pool(n::Int64) + return vcat( + one_local_pool(n), + two_local_pool(n), + ) + end + + #= TODO: Where should this function live? =# + """ + tile_operators(L1::Int, L2::Int, chosen_operators::Vector{Vector{ScaledPauli{N}}}, PBCs) + + Constructs the tiled operators for a system of `L2` qubits, given a set of operators + defined for a smaller problem instance on `L1` qubits. + + # Parameters + - `L1`: number of qubits for small problem instance + - `L2`: number of qubits for large problem instance + - `chosen_operators`: list of operators for small problem instance + - `PBCs`: periodic boundary conditions + + # Returns + - `tiled_ops`: tiled operators as a Vector{Vector{ScaledPauli}} + """ + function tile_operators(L1::Int, L2::Int, chosen_operators::Vector{Vector{ScaledPauli{N}}}, PBCs) where {N} + @assert L2 >= L1 "L2 must be greater than L1." + if L2 == L1 + return chosen_operators + end + l = length(chosen_operators) + tiled_ops = ScaledPauliVector{L2}[] + permutations = PBCs ? (L2) : (L2-L1+1) + for i=1:l + spv = chosen_operators[i] + for j = 1:permutations + tiled_spv = ScaledPauli{L2}[] + id_left = Pauli("I"^(j-1)) + id_right = Pauli("I"^(L2-L1+1-j)) + for sp in spv + tiled_sp = (j-1)==0 ? sp : otimes(id_left,sp) + tiled_sp = (L2-L1+1-j)==0 ? tiled_sp : otimes(tiled_sp, id_right) + push!(tiled_spv, tiled_sp) + end + duplicate = false + for op in tiled_ops + if tiled_spv ≈ op + duplicate = true + break + end + end + if !duplicate + push!(tiled_ops, tiled_spv) + end + end + end + return tiled_ops + end +end + \ No newline at end of file diff --git a/src/degenerateADAPT/DegenerateADAPT.jl b/src/degenerateADAPT/DegenerateADAPT.jl new file mode 100644 index 0000000..08d199b --- /dev/null +++ b/src/degenerateADAPT/DegenerateADAPT.jl @@ -0,0 +1,79 @@ +import ..ADAPT + +""" + DegenerateADAPT + +Score pool operators by their initial gradients if they were to be appended to the ansatz. +Equivalently, score pool operators by the expectation value + of the commutator of the pool operator with the observable. +In the case where the largest scores (gradients) are degenerate between multiple pool operators, choose the +operator to append to the ansatz randomly. + +""" +struct DegenerateADAPT <: ADAPT.AdaptProtocol end +DEG_ADAPT = DegenerateADAPT() + +ADAPT.typeof_score(::DegenerateADAPT) = Float64 + +function ADAPT.adapt!( + ansatz::ADAPT.AbstractAnsatz, + trace::ADAPT.Trace, + adapt_type::DegenerateADAPT, + pool::ADAPT.GeneratorList, + observable::ADAPT.Observable, + reference::ADAPT.QuantumState, + callbacks::ADAPT.CallbackList, +) + # CALCULATE SCORES + scores = ADAPT.calculate_scores(ansatz, adapt_type, pool, observable, reference) + + # CHECK FOR CONVERGENCE + ε = eps(ADAPT.typeof_score(adapt_type)) + if all(score -> abs(score) < ε, scores) + ADAPT.set_converged!(ansatz, true) + return false + end + + # MAKE SELECTION + largest_score = maximum(scores); indices_of_largest_scores = findall(a->a==largest_score, scores) +# println("gradient degeneracy: ",length(indices_of_largest_scores)) +# println("operators with degenerate max. gradients: ",pool[indices_of_largest_scores]) + selected_index = rand(indices_of_largest_scores) + selected_score = scores[selected_index] + selected_generator = pool[selected_index] + selected_parameter = zero(ADAPT.typeof_parameter(ansatz)) + + # DEFER TO CALLBACKS + data = ADAPT.Data( + :scores => scores, + :selected_index => selected_index, + :selected_score => selected_score, + :selected_generator => selected_generator, + :selected_parameter => selected_parameter, + ) + + stop = false + for callback in callbacks + stop = stop || callback(data, ansatz, trace, adapt_type, pool, observable, reference) + # Note that, as soon as `stop` is true, subsequent callbacks are short-circuited. + end + (stop || ADAPT.is_converged(ansatz)) && return false + + # PERFORM ADAPTATION + push!(ansatz, selected_generator => selected_parameter) + ADAPT.set_optimized!(ansatz, false) + return true +end + +function ADAPT.calculate_score( + ansatz::ADAPT.AbstractAnsatz, + adapt_type::DegenerateADAPT, + generator::ADAPT.Generator, + observable::ADAPT.Observable, + reference::ADAPT.QuantumState, +) + L = length(ansatz) + candidate = deepcopy(ansatz) + push!(candidate, generator => zero(ADAPT.typeof_parameter(ansatz))) + return abs(ADAPT.partial(L+1, candidate, observable, reference)) +end \ No newline at end of file diff --git a/src/hamiltonians/latticemodels.jl b/src/hamiltonians/latticemodels.jl new file mode 100644 index 0000000..446620b --- /dev/null +++ b/src/hamiltonians/latticemodels.jl @@ -0,0 +1,100 @@ +import LinearAlgebra: I + +import PauliOperators +import PauliOperators: FixedPhasePauli, Pauli, ScaledPauli +import PauliOperators: ScaledPauliVector, PauliSum +import PauliOperators: clip!, jordan_wigner + +""" + xyz_model(L::Int, Jx::Float, Jy::Float, Jz::Float, PBCs::Bool) + +An XYZ Heisenberg Hamiltonian. + +# Parameters +- `L`: system size. +- `Jx`: coupling along X. +- `Jy`: coupling along Y. +- `Jz`: coupling along Z. +- `PBCs`: Periodic Boundary Conditions + +# Returns +- `PauliOperators.PauliSum`: the Hamiltonian + +""" +function xyz_model(L,Jx,Jy,Jz,PBCs) + H = PauliSum(L) + + for i=1:L-1 + term = PauliSum(L) + term += Jx*Pauli(L; X=[i,i+1]) + term += Jy*Pauli(L; Y=[i,i+1]) + term += Jz*Pauli(L; Z=[i,i+1]) + clip!(term) + sum!(H,term) + end + if PBCs + term = PauliSum(L) + term += Jx*Pauli(L; X=[L,1]) + term += Jy*Pauli(L; Y=[L,1]) + term += Jz*Pauli(L; Z=[L,1]) + clip!(term) + sum!(H,term) + end + + return H +end + +""" + hubbard_jw(graph::Array{T,2}, U, t) + +A Hubbard Hamiltonian in the Jordan-Wigner basis. + +Copied shamelessly from Diksha's ACSE repository. + +# Parameters +- `graph`: an adjacency matrix identifying couplings. Must be symmetric. +- `U`: Coulomb interaction for all sites +- `t`: hopping energy for all couplings + +# Returns +- `PauliOperators.PauliSum`: the Hamiltonian + +""" +function hubbard_hamiltonian(graph::Matrix{T}, U, t) where T + Ni, Nj = size(graph) + Ni == Nj || throw(DimensionMismatch) + Norb = Ni + N = 2*Norb + H = PauliSum(N) + for i in 1:Norb + ia = 2*i-1 + ib = 2*i + for j in i+1:Norb + ja = 2*j-1 + jb = 2*j + abs(graph[i,j]) > 1e-16 || continue + + tij = jordan_wigner(ia, N)*jordan_wigner(ja,N)' + tij += jordan_wigner(ib, N)*jordan_wigner(jb,N)' + tij += adjoint(tij) + clip!(tij) + sum!(H,t*tij) + end + + ni = jordan_wigner(ia, N)*jordan_wigner(ia,N)'*jordan_wigner(ib, N)*jordan_wigner(ib,N)' + sum!(H,U*ni) + end + return H +end + +""" + hubbard_hamiltonian(L::Int, U, t; pbc=false) + +Convenience constructor for a 1D nearest-neighbor Hubbard model with L sites. + +""" +function hubbard_hamiltonian(L::Int, U, t; pbc=false) + A = Matrix{Bool}(I, L, L) + pbc && (A[L,1] = A[1,L] = true) + return hubbard_hamiltonian(A, U, t) +end diff --git a/test/hubbard_qeb.jl b/test/hubbard_qeb.jl index 9793df7..e89c45a 100644 --- a/test/hubbard_qeb.jl +++ b/test/hubbard_qeb.jl @@ -9,19 +9,10 @@ u = 0.25 # BUILD OUT THE PROBLEM HAMILTONIAN: a periodic 1d Hubbard lattice at u=0.25 U = 4*u # Dimensionless parameter u ≡ U/4|t|, and we'll set units so |t|=1. -H = ADAPT.Operators.hubbard_hamiltonian(L, U, -1.0, pbc=true) +H = ADAPT.Hamiltonians.hubbard_hamiltonian(L, U, -1.0, pbc=true) # BUILD OUT THE QUBIT-EXCITATION POOL -pool = ScaledPauliVector{2L}[] -for i in 1:2L -for j in i+1:2L - push!(pool, ADAPT.Operators.qubitexcitation(2L, i, j)) - - for k in j+1:2L - for l in k+1:2L - push!(pool, ADAPT.Operators.qubitexcitation(2L, i, j, k, l)) - end; end -end; end +pool, target_and_source = ADAPT.Pools.qubitexcitationpool(2L) # CONSTRUCT A REFERENCE STATE neel = "0110"^(L >> 1) diff --git a/test/runtests.jl b/test/runtests.jl index 77c734f..5b3fac1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -136,6 +136,18 @@ end ) end + @testset "Degenerate-ADAPT" begin + ADAPT.validate( + ADAPT.Ansatz(Float64, pools[:ScaledPauliVector]), + ADAPT.Degenerate_ADAPT.DEG_ADAPT, + BFGS, + pools[:ScaledPauliVector], + observables[:PauliSum], + references[:Vector]; + label = "ScaledPauli[] Pool, Statevector", + ) + end + @testset "Overlap" begin overlap = let ψ = zeros(ComplexF64, 1<> 1); (L_small & 1 == 1) && (neel *= "0"); neel_index = parse(Int128, neel, base=2) + +# BUILD OUT THE POOL +if pooltype=="fullpauli" + pool = ADAPT.Pools.fullpauli(L_small) +elseif pooltype == "qubitexcitation" + pool, target_and_source = ADAPT.Pools.qubitexcitationpool(L_small) +elseif pooltype == "qubitadapt" + pool = ADAPT.Pools.qubitadaptpool(L_small) +end +println("pool size ",length(pool)); # println("pool: ",pool) + +# SELECT THE PROTOCOLS +adapt = ADAPT.Degenerate_ADAPT.DEG_ADAPT +vqe = ADAPT.OptimOptimizer(:BFGS; g_tol=1e-6) +gradient_cutoff = 1e-4; maxiters = 500 + +# SELECT THE CALLBACKS +callbacks = [ + ADAPT.Callbacks.Tracer(:energy, :selected_generator, :selected_index, :selected_score, :scores), +# ADAPT.Callbacks.Printer(:energy, :selected_generator, :selected_score), + ADAPT.Callbacks.ScoreStopper(gradient_cutoff), + ADAPT.Callbacks.ParameterStopper(maxiters), +] + +# EXACT DIAGONALIZATION +module Exact + import ..XXZHam, ..L_small + using LinearAlgebra + Hm = Matrix(XXZHam); E, U = eigen(Hm) # NOTE: Comment out after first run when debugging. + ψ0 = U[:,1] + E0 = E[1] +end +println("Exact gs energy: ",Exact.E0) + +# RUN MANY ADAPT-VQE TRIALS ON THE SMALL PROBLEM INSTANCE +chosen_operators = ScaledPauliVector{L_small}[] +trials = 100 +println("running trials on the small problem instance...") +for i=1:trials + # INITIALIZE THE REFERENCE STATE + ψ0 = zeros(ComplexF64,2^L_small); + ψ0[neel_index+1] = 1.0 # Neel + + # INITIALIZE THE ANSATZ AND TRACE + ansatz = ADAPT.Ansatz(Float64, pool) + trace = ADAPT.Trace() + + # RUN THE ALGORITHM + ADAPT.run!(ansatz, trace, adapt, vqe, pool, XXZHam, ψ0, callbacks) + selected_operators = trace[:selected_generator][1:end-1] +# energy_err = abs((Exact.E0-(trace[:energy][end-1]))/Exact.E0); print(" ",energy_err) + for selected_op in selected_operators + duplicate = false + for chosen_op in chosen_operators + if selected_op ≈ chosen_op + duplicate = true + break + end + end + if !duplicate + push!(chosen_operators, selected_op) + end + end +end +println("\nTrials on small problem instance complete.") +println(length(chosen_operators), " chosen operators: ",chosen_operators,"\n") + +# # FIND DLA OF OPERATORS +# t_0_Lie = time(); Lie_alg = Lie_algebra_elements(chosen_operators); t_f_Lie = time(); dt = (t_f_Lie - t_0_Lie)/60.0 +# # println("walltime to calculate Lie algebra for L = $L_tile was $dt minutes\n") +# println("DLA of chosen operators ",length(Lie_alg), " ",Lie_alg,"\n") + + +# RUN ADAPT-VQE ON THE LARGE PROBLEM INSTANCE + +# BUILD OUT THE PROBLEM HAMILTONIAN: an open XXZ model +XXZHam = ADAPT.Hamiltonians.xyz_model(L_large, Jxy, Jxy, Jz, PBCs) + +# EXACT DIAGONALIZATION +module Exact_large + import ..XXZHam, ..L_large + using LinearAlgebra + Hm = Matrix(XXZHam); E, U = eigen(Hm) # NOTE: Comment out after first run when debugging. + ψ0 = U[:,1] + E0 = E[1] +end +println("Exact gs energy: ",Exact_large.E0) + +# CONSTRUCT A REFERENCE STATE +neel = "01"^(L_large >> 1); (L_large & 1 == 1) && (neel *= "0"); neel_index = parse(Int128, neel, base=2) +ψ0 = zeros(ComplexF64,2^L_large); ψ0[neel_index+1] = 1.0 + +# BUILD OUT THE POOL +# Lie_algebra_elements(chosen_op_strings) # calculate the Lie algebra of the chosen operators +pool = ADAPT.Pools.tile_operators(L_small, L_large, chosen_operators, PBCs) + +# SELECT THE PROTOCOLS +adapt = ADAPT.VANILLA +vqe = ADAPT.OptimOptimizer(:BFGS; g_tol=1e-6) + +# INITIALIZE THE ANSATZ AND TRACE +ansatz = ADAPT.Ansatz(Float64, pool) +trace = ADAPT.Trace() + +# RUN THE ALGORITHM +ADAPT.run!(ansatz, trace, adapt, vqe, pool, XXZHam, ψ0, callbacks) + +# RESULTS +ψEND = ADAPT.evolve_state(ansatz, ψ0) +num_ADAPT_iters = length(ansatz) +E0 = trace[:energy][end-1]; rel_energy_err = abs((Exact_large.E0 - E0)/(Exact_large.E0)) +println("final ADAPT energy = ", E0) +println("ansatz length: ", length(ansatz)) +println("relative energy error = ",rel_energy_err) \ No newline at end of file