diff --git a/lib/BloqadeQMC/src/BloqadeQMC.jl b/lib/BloqadeQMC/src/BloqadeQMC.jl index af1e3269d..b55739db3 100644 --- a/lib/BloqadeQMC/src/BloqadeQMC.jl +++ b/lib/BloqadeQMC/src/BloqadeQMC.jl @@ -59,6 +59,7 @@ include("hamiltonian.jl") include("operatorsamplers/improved_op_sampler.jl") include("diagnostics/Diagnostics.jl") include("ising/Ising.jl") +include("updates/updates.jl") include("measurements.jl") include("error.jl") diff --git a/lib/BloqadeQMC/src/datastructures/Rydberg.jl b/lib/BloqadeQMC/src/datastructures/Rydberg.jl new file mode 100644 index 000000000..47837321c --- /dev/null +++ b/lib/BloqadeQMC/src/datastructures/Rydberg.jl @@ -0,0 +1,155 @@ +using Base.Iterators + +# changed type tree s.t. AbstractRydberg is a direct subtype of Hamiltonian +abstract type AbstractRydberg{O <: AbstractOperatorSampler} <: Hamiltonian{2,O} end + +struct Rydberg{O,M <: UpperTriangular{Float64},UΩ <: AbstractVector{Float64}, Uδ <: AbstractVector{Float64}, L <: Lattice} <: AbstractRydberg{O} + op_sampler::O + V::M + Ω::UΩ + δ::Uδ + lattice::L + energy_shift::Float64 +end + +nspins(H::Rydberg) = nspins(H.lattice) + + +function make_prob_vector(H::Type{<:AbstractRydberg}, V::UpperTriangular{T}, Ω::AbstractVector{T}, δ::AbstractVector{T}; epsilon=0.0) where T + @assert length(Ω) == length(δ) == size(V, 1) == size(V, 2) + @assert (0.0 <= epsilon <= 1.0) "epsilon must be in the range [0, 1]!" + + ops = Vector{NTuple{ISING_OP_SIZE, Int}}() + p = Vector{T}() + energy_shift = zero(T) + + for i in eachindex(Ω) + if !iszero(Ω[i]) + push!(ops, makediagonalsiteop(AbstractLTFIM, i)) + push!(p, Ω[i] / 2) + energy_shift += Ω[i] / 2 + end + end + + Ns = length(Ω) + bond_spins = Set{NTuple{2,Int}}() + coordination_numbers = zeros(Int, Ns) + for j in axes(V, 2), i in axes(V, 1) + if i < j && !iszero(V[i, j]) + push!(bond_spins, (i, j)) + coordination_numbers[i] += 1 + coordination_numbers[j] += 1 + end + end + + n = diagonaloperator(H) + I = Diagonal(LinearAlgebra.I, 2) + + # TODO: add fictitious bonds if there's a z-field on an "unbonded" site + for (site1, site2) in bond_spins + # by this point we can assume site1 < site2 + δb1 = δ[site1] / coordination_numbers[site1] + δb2 = δ[site2] / coordination_numbers[site2] + local_H = V[site1, site2]*kron(n, n) - δb1*kron(n, I) - δb2*kron(I, n) + + p_spins = -diag(local_H) + C = abs(min(0, minimum(p_spins))) + epsilon*abs(minimum(p_spins[2:end])) + #dont use the zero matrix element for the epsilon shift + p_spins .+= C + energy_shift += C + + for (t, p_t) in enumerate(p_spins) + push!(p, p_t) + push!(ops, (2, t, length(p), site1, site2)) + end + end + + return ops, p, energy_shift +end + + +############################################################################### + +# function Rydberg(dims::NTuple{D, Int}, R_b, Ω, δ; pbc=true, trunc::Int=0, epsilon::Float64=0) where D +# if D == 1 +# lat = Chain(dims[1], 1.0, pbc) +# elseif D == 2 +# lat = Rectangle(dims[1], dims[2], 1.0, 1.0, pbc) +# else +# error("Unsupported number of dimensions. 1- and 2-dimensional lattices are supported.") +# end +# return Rydberg(lat, R_b, Ω, δ; trunc=trunc, epsilon=epsilon) +# end + + +# We only want one Rydberg() interface, integrated into BloqadeExpr. +# Probably, we'll want one function generate_interaction_matrix() that creates V and then feeds that matrix into Rydberg(). +# Check with Phil how that function is coming along. + +function Rydberg(lat::Lattice, R_b::Float64, Ω::Float64, δ::Float64; trunc::Int=0, epsilon=0) + Ns = nspins(lat) + V = zeros(Float64, Ns, Ns) + + if trunc > 0 + _dist = sort!(collect(Set(lat.distance_matrix))) + uniq_dist = Vector{Float64}(undef, 0) + for i in eachindex(_dist) + if length(uniq_dist) == 0 + push!(uniq_dist, _dist[i]) + elseif !(last(uniq_dist) ≈ _dist[i]) + push!(uniq_dist, _dist[i]) + end + end + smallest_k = sort!(uniq_dist)[1:(trunc+1)] + dist = copy(lat.distance_matrix) + for i in eachindex(dist) + if dist[i] > last(smallest_k) && !(dist[i] ≈ last(smallest_k)) + dist[i] = zero(dist[i]) + end + end + elseif lat isa Rectangle && all(lat.PBC) + V = zeros(Ns, Ns) + K = 3 # figure out an efficient way to set this dynamically + + dist = zeros(Ns, Ns) + for v2 in -K:K, v1 in -K:K + dV = zeros(Ns, Ns) + for x2 in axes(dV, 2), x1 in axes(dV, 1) + i1, j1 = divrem(x1 - 1, lat.n1) + i2, j2 = divrem(x2 - 1, lat.n1) + r = [i2 - i1 + v1*lat.n1, j2 - j1 + v2*lat.n2] + dV[x1, x2] += Ω * (R_b/norm(r, 2))^6 + end + # @show v2, v1, maximum(abs, dV) + + V += dV + end + + V = (V + V') / 2 # should already be symmetric but just in case + + return Rydberg(UpperTriangular(triu!(V, 1)), Ω*ones(Ns), δ*ones(Ns), lat; epsilon=epsilon) + else + dist = lat.distance_matrix + end + + @inbounds for i in 1:(Ns-1) + for j in (i+1):Ns + # a zero entry in distance_matrix means there should be no bond + V[i, j] = dist[i, j] != 0.0 ? Ω * (R_b / dist[i, j])^6 : 0.0 + end + end + V = UpperTriangular(triu!(V, 1)) + + return Rydberg(V, Ω*ones(Ns), δ*ones(Ns), lat; epsilon=epsilon) +end + +function Rydberg(V::AbstractMatrix{T}, Ω::AbstractVector{T}, δ::AbstractVector{T}, lattice::Lattice; epsilon=zero(T)) where T + ops, p, energy_shift = make_prob_vector(AbstractRydberg, V, Ω, δ, epsilon=epsilon) + op_sampler = ImprovedOperatorSampler(AbstractLTFIM, ops, p) + return Rydberg{typeof(op_sampler), typeof(V), typeof(Ω), typeof(δ), typeof(lattice)}(op_sampler, V, Ω, δ, lattice, energy_shift) +end + +# Check whether the two functions below are needed in updates. + +total_hx(H::Rydberg)::Float64 = sum(H.Ω) / 2 +haslongitudinalfield(H::AbstractRydberg) = !iszero(H.δ) diff --git a/lib/BloqadeQMC/src/datastructures/alias.jl b/lib/BloqadeQMC/src/datastructures/alias.jl new file mode 100644 index 000000000..33551a487 --- /dev/null +++ b/lib/BloqadeQMC/src/datastructures/alias.jl @@ -0,0 +1,90 @@ +############################################################################### +# Walker's Alias Method: draw samples in O(1) time, initialize vector in O(n) +# caveat: unsure if it's possible to build an efficient update scheme +# for these types of probability vectors + +struct ProbabilityAlias{T} <: AbstractProbabilityVector{T} + normalization::T + + cutoffs::Vector{T} + alias::Vector{Int} + + # initialize using Vose's algorithm + function ProbabilityAlias{T}(weights::Vector{T}) where {T <: AbstractFloat} + if length(weights) == 0 + throw(ArgumentError("weight vector must have non-zero length!")) + end + if any(x -> x < zero(T), weights) + throw(ArgumentError("weights must be non-negative!")) + end + + weights_ = copy(weights) + normalization = sum(weights) + N = length(weights) + avg = normalization / N + + underfull = Stack{Int}() + overfull = Stack{Int}() + + for (i, w) in enumerate(weights_) + if w >= avg + push!(overfull, i) + else + push!(underfull, i) + end + end + + cutoffs = zeros(float(T), N) + alias = zeros(Int, N) + + while !isempty(underfull) && !isempty(overfull) + less = pop!(underfull) + more = pop!(overfull) + + cutoffs[less] = weights_[less] * N + alias[less] = more + + weights_[more] += weights_[less] + weights_[more] -= avg + + if weights_[more] >= avg + push!(overfull, more) + else + push!(underfull, more) + end + end + + while !isempty(underfull) + cutoffs[pop!(underfull)] = normalization + end + + while !isempty(overfull) + cutoffs[pop!(overfull)] = normalization + end + + cutoffs /= normalization + + new{T}(sum(weights), cutoffs, alias) + end +end + +ProbabilityAlias(p::Vector{T}) where T = ProbabilityAlias{T}(p) +@inline length(pvec::ProbabilityAlias) = length(pvec.cutoffs) +@inline normalization(pvec::ProbabilityAlias) = pvec.normalization + +function show(io::IO, p::ProbabilityAlias{T}) where T + r = repr(p.normalization; context=IOContext(io, :limit=>true)) + print(io, "ProbabilityAlias{$T}($r)") +end + +# function Base.rand(rng::AbstractRNG, pvec::ProbabilityAlias{T}) where T +# u, i::Int = modf(muladd(length(pvec), rand(rng), 1.0)) +# return @inbounds (u < pvec.cutoffs[i]) ? i : pvec.alias[i] +# end + +# xorshift prngs seem to be noticably faster if you just sample from them twice +function Base.rand(rng::AbstractRNG, pvec::ProbabilityAlias{T}) where T + u = rand(rng) + i = rand(rng, 1:length(pvec)) + return @inbounds (u < pvec.cutoffs[i]) ? i : pvec.alias[i] +end diff --git a/lib/BloqadeQMC/src/datastructures/hamiltonian.jl b/lib/BloqadeQMC/src/datastructures/hamiltonian.jl new file mode 100644 index 000000000..edd48eedf --- /dev/null +++ b/lib/BloqadeQMC/src/datastructures/hamiltonian.jl @@ -0,0 +1,52 @@ +# To befirst merged with Rydberg (the only specific Hamiltonian we need) and eventually integrated with BloqadeExpr interfaces +abstract type Hamiltonian{D,O<:AbstractOperatorSampler} end + +localdim(::Hamiltonian{D}) where {D} = D + +zero(H::Hamiltonian{2}) = zeros(Bool, nspins(H)) +zero(H::Hamiltonian) = zeros(Int, nspins(H)) +one(H::Hamiltonian{2}) = ones(Bool, nspins(H)) +one(H::Hamiltonian) = ones(Int, nspins(H)) + +# nspins(H::Hamiltonian) = H.Ns # Rydberg has its own version +nbonds(H::Hamiltonian) = H.Nb # Is this needed for Rydberg? + +@inline isdiagonal(H) = op -> isdiagonal(H, op) +@inline isidentity(H) = op -> isidentity(H, op) +@inline issiteoperator(H) = op -> issiteoperator(H, op) +@inline isbondoperator(H) = op -> isbondoperator(H, op) +@inline getsite(H) = op -> getsite(H, op) +@inline getbondsites(H) = op -> getbondsites(H, op) + +@inline diag_update_normalization(H::Hamiltonian) = normalization(H.op_sampler) + +############################################################################### + +# These energy functions are only used in Rydberg ED not in LTFIM QMC + +energy(::Type{<:BinaryThermalState}, H::Hamiltonian, β::Float64, n::Int) = H.energy_shift - (n / β) +function energy(::Type{<:BinaryThermalState}, H::Hamiltonian, β::Float64, ns::Vector{T}) where {T <: Real} + E = -mean_and_stderr(ns) / β + return H.energy_shift + E +end +energy(qmc_state::BinaryQMCState, args...; kwargs...) = energy(typeof(qmc_state), args...; kwargs...) + +energy_density(S::Type{<:BinaryQMCState}, H::Hamiltonian, args...; kwargs...) = + energy(S, H, args...; kwargs...) / nspins(H) + +energy_density(qmc_state::BinaryQMCState, args...; kwargs...) = energy_density(typeof(qmc_state), args...; kwargs...) + +############################################################################### + +# Move this stuff to state files? + +function BinaryGroundState(H::Hamiltonian{2,O}, M::Int, trialstate::Union{Nothing, AbstractTrialState}=nothing) where {K, O <: AbstractOperatorSampler{K}} + z = zero(H) # Here we are initializing state into all zeros + BinaryGroundState(z, init_op_list(2*M, Val{K}()), trialstate)::BinaryGroundState{K, typeof(z)} +end + + +function BinaryThermalState(H::Hamiltonian{2,O}, cutoff::Int) where {K, O <: AbstractOperatorSampler{K}} # cutoff is essentially 2M but here we tune it via an algorithm + z = zero(H) + BinaryThermalState(z, init_op_list(cutoff, Val{K}()))::BinaryThermalState{K, typeof(z)} +end diff --git a/lib/BloqadeQMC/src/datastructures/improved_op_sampler.jl b/lib/BloqadeQMC/src/datastructures/improved_op_sampler.jl new file mode 100644 index 000000000..2b7b589f1 --- /dev/null +++ b/lib/BloqadeQMC/src/datastructures/improved_op_sampler.jl @@ -0,0 +1,75 @@ +# Same file as previous improved_operator_sampler.jl just with Abstract type copied over from operatorsampler.jl + +# Do we still need AbstractOperatorSampler at all? Probably we don't want a structure called "Improvedxxx" once we don't have the structure called "xxx" anymore :) + +abstract type AbstractOperatorSampler{K, T, P <: AbstractProbabilityVector{T}} end +firstindex(::AbstractOperatorSampler) = 1 +lastindex(os::AbstractOperatorSampler) = length(os) +@inline normalization(os::AbstractOperatorSampler) = normalization(os.pvec) + +abstract type AbstractImprovedOperatorSampler{K, T, P} <: AbstractOperatorSampler{K, T, P} end + +struct ImprovedOperatorSampler{K, T, P} <: AbstractImprovedOperatorSampler{K, T, P} + operators::Vector{NTuple{K, Int}} + pvec::P + op_log_weights::Vector{T} +end + +# only supports the LTFIM/Rydberg cases for now +function ImprovedOperatorSampler(H::Type{<:Hamiltonian{2, <:AbstractOperatorSampler}}, operators::Vector{NTuple{K, Int}}, p::Vector{T}) where {T <: AbstractFloat, K} + @assert length(operators) == length(p) "Given vectors must have the same length!" + + op_log_weights = log.(p) + + max_mel_ops = Vector{NTuple{K, Int}}() + p_modified = Vector{T}() + + # fill with all the site operators first + for (i, op) in enumerate(operators) + if issiteoperator(H, op) && isdiagonal(H, op) + push!(max_mel_ops, op) + push!(p_modified, @inbounds p[i]) + end + end + idx = findall(isbondoperator(H), operators) + ops = operators[idx] + p_mod = p[idx] + + perm = sortperm(ops, by=getbondsites(H)) + ops = ops[perm] + p_mod = p_mod[perm] + + op_groups = Dict{NTuple{2, Int}, Vector{NTuple{K, Int}}}() + p_groups = Dict{NTuple{2, Int}, Vector{T}}() + + while !isempty(ops) + op = pop!(ops) + site1, site2 = getbondsites(H, op) + p_ = pop!(p_mod) + + if haskey(op_groups, (site1, site2)) + push!(op_groups[(site1, site2)], op) + push!(p_groups[(site1, site2)], p_) + else + op_groups[(site1, site2)] = [op] + p_groups[(site1, site2)] = [p_] + end + end + + for (k, p_gr) in p_groups + i = argmax(p_gr) + push!(max_mel_ops, op_groups[k][i]) + push!(p_modified, p_gr[i]) + end + + pvec = probability_vector(p_modified) + return ImprovedOperatorSampler{K, T, typeof(pvec)}(max_mel_ops, pvec, op_log_weights) +end + +@inline rand(rng::AbstractRNG, os::ImprovedOperatorSampler) = @inbounds os.operators[rand(rng, os.pvec)] + +Base.@propagate_inbounds getlogweight(os::ImprovedOperatorSampler, w::Int) = os.op_log_weights[w] + +@inline length(os::ImprovedOperatorSampler) = length(os.operators) + +############################################################################## diff --git a/lib/BloqadeQMC/src/datastructures/operators.jl b/lib/BloqadeQMC/src/datastructures/operators.jl new file mode 100644 index 000000000..95d3c2e90 --- /dev/null +++ b/lib/BloqadeQMC/src/datastructures/operators.jl @@ -0,0 +1,73 @@ +############################################################################### + +# TFIM ops: +# (1,-2,i,0,i) is an off-diagonal site operator h*sigma^x_i +# (1,1,i,0,i) is a diagonal site operator h +# (0,0,0,0,0) is the identity operator I - NOT USED IN THE PROJECTOR CASE +# (2,1,w,i,j) is a diagonal bond operator J(sigma^z_i sigma^z_j) + + +############################################################################### + +# LTFIM ops: +# (-2,i,i) is an off-diagonal site operator h(sigma^+_i + sigma^-_i) +# (-1,i,i) is a diagonal site operator h +# (0,0,0) is the identity operator I - NOT USED IN THE PROJECTOR CASE +# (t,i,j) is a diagonal bond operator J(sigma^z_i sigma^z_j) + hzb(sigma^z_i + sigma^z_j) +# t denotes the spin config at sites i,j, subtract 1 then convert to binary +# t = 1 -> 00 -> down-down +# t = 2 -> 01 -> down-up +# t = 3 -> 10 -> up-down +# t = 4 -> 11 -> up-up +# spin_config(t) = divrem(t - 1, 2) + +@inline getbondtype(::AbstractLTFIM, s1::Bool, s2::Bool) = (s1<<1 | s2) + 1 +@inline spin_config(::AbstractLTFIM, t::Int)::NTuple{2,Int} = divrem(t - 1, 2) +@inline spin_config(H::AbstractLTFIM, op::NTuple{ISING_OP_SIZE, Int}) = @inbounds spin_config(H, getoperatortype(H, op)) + +############################################################################### + +@inline getoperatorlocality(::Type{<:AbstractIsing}, op::NTuple{ISING_OP_SIZE, Int}) = @inbounds op[1] +@inline getoperatorlocality(H::AbstractIsing, op::NTuple{ISING_OP_SIZE, Int}) = getoperatorlocality(typeof(H), op) +@inline getoperatortype(::Type{<:AbstractIsing}, op::NTuple{ISING_OP_SIZE, Int}) = @inbounds op[2] +@inline getoperatortype(H::AbstractIsing, op::NTuple{ISING_OP_SIZE, Int}) = getoperatortype(typeof(H), op) +@inline getweightindex(::Type{<:AbstractIsing}, op::NTuple{ISING_OP_SIZE, Int}) = @inbounds op[3] +@inline getweightindex(H::AbstractIsing, op::NTuple{ISING_OP_SIZE, Int}) = getweightindex(typeof(H), op) +@inline getsite(::Type{<:AbstractIsing}, op::NTuple{ISING_OP_SIZE, Int}) = @inbounds op[end] +@inline getsite(H::AbstractIsing, op::NTuple{ISING_OP_SIZE, Int}) = getsite(typeof(H), op) +@inline getbondsites(::Type{<:AbstractIsing}, op::NTuple{ISING_OP_SIZE, Int}) = @inbounds (op[end-1], op[end]) +@inline getbondsites(H::AbstractIsing, op::NTuple{ISING_OP_SIZE, Int}) = getbondsites(typeof(H), op) + +@inline convertoperatortype(H::Type{<:AbstractIsing}, op::NTuple{ISING_OP_SIZE, Int}, new_t::Int) = + @inbounds (op[1], new_t, getweightindex(H, op) - getoperatortype(H, op) + new_t, op[4], op[5]) +@inline convertoperatortype(H::AbstractIsing, op::NTuple{ISING_OP_SIZE, Int}, new_t::Int) = convertoperatortype(typeof(H), op, new_t) + +@inline isidentity(H::Type{<:AbstractIsing}, op::NTuple{ISING_OP_SIZE, Int}) = (getoperatorlocality(H, op) == 0) +@inline isidentity(H::AbstractIsing, op::NTuple{ISING_OP_SIZE, Int}) = isidentity(typeof(H), op) +@inline issiteoperator(H::Type{<:AbstractIsing}, op::NTuple{ISING_OP_SIZE, Int}) = (getoperatorlocality(H, op) == 1) +@inline issiteoperator(H::AbstractIsing, op::NTuple{ISING_OP_SIZE, Int}) = issiteoperator(typeof(H), op) +@inline isbondoperator(H::Type{<:AbstractIsing}, op::NTuple{ISING_OP_SIZE, Int}) = (getoperatorlocality(H, op) == 2) +@inline isbondoperator(H::AbstractIsing, op::NTuple{ISING_OP_SIZE, Int}) = isbondoperator(typeof(H), op) +@inline isdiagonal(H::Type{<:AbstractIsing}, op::NTuple{ISING_OP_SIZE, Int}) = (getoperatortype(H, op) != -2) +@inline isdiagonal(H::AbstractIsing, op::NTuple{ISING_OP_SIZE, Int}) = isdiagonal(typeof(H), op) + +@inline makeidentity(::Type{<:AbstractIsing})::NTuple{ISING_OP_SIZE, Int} = (0, 0, 0, 0, 0) +@inline makediagonalsiteop(::Type{<:AbstractIsing}, i::Int)::NTuple{ISING_OP_SIZE, Int} = (1, 1, i, 0, i) +@inline makeoffdiagonalsiteop(::Type{<:AbstractIsing}, i::Int)::NTuple{ISING_OP_SIZE, Int} = (1, -2, i, 0, i) +@inline makeidentity(H::AbstractIsing) = makeidentity(typeof(H)) +@inline makediagonalsiteop(H::AbstractIsing, i::Int) = makediagonalsiteop(typeof(H), i) +@inline makeoffdiagonalsiteop(H::AbstractIsing, i::Int) = makeoffdiagonalsiteop(typeof(H), i) + +@inline getbondtype(::AbstractTFIM, s1::Bool, s2::Bool) = 1 + +@inline diagonaloperator(::Type{<:AbstractIsing}) = Diagonal([-1, 1]) +@inline diagonaloperator(H::AbstractIsing) = diagonaloperator(typeof(H)) + +@inline getlogweight(H::AbstractIsing, op::NTuple{ISING_OP_SIZE, Int}) = @inbounds getlogweight(H.op_sampler, getweightindex(H, op)) + +############################################################################### + +# from previous Rydberg.jl: + +@inline diagonaloperator(::Type{<:AbstractRydberg}) = Diagonal([0, 1]) +@inline diagonaloperator(H::AbstractRydberg) = diagonaloperator(typeof(H)) diff --git a/lib/BloqadeQMC/src/datastructures/probabilityvector.jl b/lib/BloqadeQMC/src/datastructures/probabilityvector.jl new file mode 100644 index 000000000..c060dcd99 --- /dev/null +++ b/lib/BloqadeQMC/src/datastructures/probabilityvector.jl @@ -0,0 +1,20 @@ +# Understand: What is special about this vector? Some sort of efficiency? Was there something about push! / stack / memory? + +abstract type AbstractProbabilityVector{T <: Real} <: AbstractVector{T} end +Base.@propagate_inbounds getindex(pvec::AbstractProbabilityVector, i) = getweight(pvec, i) +Base.@propagate_inbounds setindex!(pvec::AbstractProbabilityVector, w, i) = setweight!(pvec, w, i) +Base.@propagate_inbounds setprobability!(pvec::AbstractProbabilityVector{T}, p::T, i::Int) where T = setweight!(pvec, p*normalization(pvec), i) +Base.@propagate_inbounds getprobability(pvec::AbstractProbabilityVector, i) = getweight(pvec, i) / normalization(pvec) + +firstindex(::AbstractProbabilityVector) = 1 +lastindex(pvec::AbstractProbabilityVector) = length(pvec) +size(pvec::AbstractProbabilityVector) = (length(pvec),) + +############################################################################### + +include("alias.jl") + +############################################################################### + +# Get the "recommended" probability vector type +probability_vector(p::Vector{T}) where T = ProbabilityAlias(p) diff --git a/lib/BloqadeQMC/src/datastructures/qmc_state.jl b/lib/BloqadeQMC/src/datastructures/qmc_state.jl new file mode 100644 index 000000000..c9b0175b6 --- /dev/null +++ b/lib/BloqadeQMC/src/datastructures/qmc_state.jl @@ -0,0 +1,157 @@ +function init_op_list(length, K=Val{3}()) # I'm guessing the 3 is the default value since the LTFIM operators are defined by 3 indices? + operator_list = [ntuple(_->0, K) for _ in 1:length] + return operator_list +end + + +abstract type AbstractStateType end +struct Thermal <: AbstractStateType; end +struct Ground <: AbstractStateType; end + +abstract type AbstractQMCState{S<:AbstractStateType,T,K} end + +const AbstractGroundState = AbstractQMCState{Ground} +const AbstractThermalState = AbstractQMCState{Thermal} + +struct QMCState{S,T,K,V <: AbstractVector{T},P <: Union{Nothing, AbstractTrialState{Float64, T}}} <: AbstractQMCState{S,T,K} + left_config::V + right_config::V + propagated_config::V + + operator_list::Vector{NTuple{K,Int}} + + linked_list::Vector{Int} + leg_types::V + associates::Vector{Int} + leg_sites::Vector{Int} + op_indices::Vector{Int} + + in_cluster::Vector{Int} + cstack::PushVector{Int, Vector{Int}} + current_cluster::PushVector{Int, Vector{Int}} + + first::Vector{Int} + last::Union{Vector{Int}, Nothing} + + trialstate::P # This is where AbstractTrialState gets used in case of GroundState simulation + + function QMCState{S, T, K, V, P}( + left_config::V, right_config::V, propagated_config::V, + operator_list, + link_list, leg_types, associates, leg_sites, op_indices, + in_cluster, cstack, current_cluster, + first, last, trialstate + ) where {S, K, T, V, P} + + if S isa Type{<:Thermal} + @assert last isa Vector{Int} + @assert P == Nothing + else + @assert last === nothing + @assert P <: AbstractTrialState + end + + new{S, T, K, V, P}( + left_config, right_config, propagated_config, + operator_list, + link_list, leg_types, associates, leg_sites, op_indices, + in_cluster, cstack, current_cluster, + first, last, trialstate + ) + end + + function QMCState{S, T, K, V}( + left_config::V, right_config::V, operator_list::Vector{NTuple{K,Int}}, + trialstate::Union{Nothing, AbstractTrialState{Float64, T}}=nothing + ) where {S, T, K, V <: AbstractVector{T}} + @assert left_config !== right_config "left_config and right_config can't be the same array!" + @assert size(left_config) === size(right_config) "left_config and right_config must have the same size!" + + if S isa Type{<:Ground} + len = 2*length(left_config) + 4*length(operator_list) + if trialstate === nothing + trialstate = PlusState{Float64, T}() + end + else + len = 4*length(operator_list) + trialstate = nothing + end + link_list = zeros(Int, len) + leg_types = similar(left_config, T, len) + associates = zeros(Int, len) + leg_sites = zeros(Int, len) + op_indices = zeros(Int, len) + + in_cluster = zeros(Int, len) + cstack = PushVector{Int}(nextpow(2, length(left_config))) + current_cluster = PushVector{Int}(nextpow(2, length(left_config))) + + first = zeros(Int, length(left_config)) + last = (S isa Type{<:Thermal}) ? copy(first) : nothing + args = [ + left_config, right_config, copy(left_config), + operator_list, + link_list, leg_types, associates, leg_sites, op_indices, + in_cluster, cstack, current_cluster, + first, last, trialstate + ] + + QMCState{S, T, K, V, typeof(trialstate)}(args...) + end +end + +QMCState{S, T, K, V}(left_config::V, operator_list, trialstate::Union{Nothing, AbstractTrialState{Float64, T}}=nothing) where {S, T, K, V} = + QMCState{S, T, K, V}(left_config, copy(left_config), operator_list, trialstate) + +QMCState{S, T, K}(left_config::V, right_config::V, operator_list, trialstate::Union{Nothing, AbstractTrialState{Float64, T}}=nothing) where {S, T, K, V} = + QMCState{S, T, K, V}(left_config, right_config, operator_list, trialstate) +QMCState{S, T, K}(left_config::V, operator_list, trialstate::Union{Nothing, AbstractTrialState{Float64, T}}=nothing) where {S, T, K, V} = + QMCState{S, T, K, V}(left_config, operator_list, trialstate) + +QMCState{S, T}(left_config::V, right_config::V, operator_list::Vector{NTuple{K,Int}}, trialstate::Union{Nothing, AbstractTrialState{Float64, T}}=nothing) where {S, T, K, V} = + QMCState{S, T, K}(left_config, right_config, operator_list, trialstate) +QMCState{S, T}(left_config, operator_list::Vector{NTuple{K,Int}}, trialstate::Union{Nothing, AbstractTrialState{Float64, T}}=nothing) where {S, T, K} = + QMCState{S, T, K}(left_config, operator_list, trialstate) + +QMCState{S}(left_config::V, right_config::V, operator_list, trialstate::Union{Nothing, AbstractTrialState}=nothing) where {S, V} = + QMCState{S, eltype(left_config)}(left_config, right_config, operator_list, trialstate) +QMCState{S}(left_config, operator_list, trialstate::Union{Nothing, AbstractTrialState}=nothing) where S = + QMCState{S, eltype(left_config)}(left_config, operator_list, trialstate) + + +const GroundState{T,K,V} = QMCState{Ground,T,K,V} +const ThermalState{T,K,V} = QMCState{Thermal,T,K,V} + +const BinaryQMCState{K,V <: AbstractVector{Bool}} = QMCState{S, Bool, K, V} where {S <: AbstractStateType} +const BinaryGroundState{K,V <: AbstractVector{Bool}} = GroundState{Bool, K, V} +const BinaryThermalState{K,V <: AbstractVector{Bool}} = ThermalState{Bool, K, V} + + +function convert(::Type{QMCState{S′, T, K, V}}, state::QMCState{S, T, K, V}) where {S, S′,T, K, V} + if S′ == S + return state + elseif S′ isa Type{<:Thermal} + len = 4*length(state.operator_list) + last = copy(state.first) + trialstate = nothing + else + # make the operator list length even by adding one identity operator + if isodd(length(state.operator_list)) + push!(state.operator_list, ntuple(_ -> 0, K)) + end + len = 2*length(state.left_config) + 4*length(state.operator_list) + last = nothing + trialstate = PlusState() + end + + resize!(state.linked_list, len) + resize!(state.leg_types, len) + resize!(state.associates, len) + resize!(state.leg_sites, len) + resize!(state.op_indices, len) + resize!(state.in_cluster, len) + + args = [getfield(state, field) for field in fieldnames(typeof(state))] + + return QMCState{S′, T, K, V}(args[1:end-2]..., last, trialstate) +end diff --git a/lib/BloqadeQMC/src/updates/cluster.jl b/lib/BloqadeQMC/src/updates/cluster.jl new file mode 100644 index 000000000..8f0dedfcc --- /dev/null +++ b/lib/BloqadeQMC/src/updates/cluster.jl @@ -0,0 +1,343 @@ +# Main parts here: huge link_list_update!() and cluster_update!() functions. +# They don´t seem to call each other -> understand when each is used & how they interact. + +function link_list_update!(::AbstractRNG, qmc_state::BinaryQMCState, H::AbstractIsing, ::Diagnostics) # fairly standard, builds necessary data structure for both line and multibranch + Ns = nspins(H) + spin_left = qmc_state.left_config + + # retrieve linked list data structures + LinkList = qmc_state.linked_list # needed for cluster update + LegType = qmc_state.leg_types + + # A diagonal bond operator has non trivial associates for cluster building + Associates = qmc_state.associates + + op_indices = qmc_state.op_indices + leg_sites = qmc_state.leg_sites + + if qmc_state isa BinaryGroundState + First = qmc_state.first + + # The first N elements of the linked list are the spins of the LHS basis state + @inbounds for i in 1:Ns + LegType[i] = spin_left[i] + Associates[i] = 0 + First[i] = i + if H isa AbstractLTFIM + op_indices[i] = 0 + leg_sites[i] = 0 + end + end + idx = Ns + else + First = fill!(qmc_state.first, 0) #initialize the First list + Last = fill!(qmc_state.last, 0) #initialize the Last list + idx = 0 + end + + spin_prop = copyto!(qmc_state.propagated_config, spin_left) # the propagated spin state + + # Now, add the 2M operators to the linked list. Each has either 2 or 4 legs + @inbounds for (n, op) in enumerate(qmc_state.operator_list) + if issiteoperator(H, op) + site = getsite(H, op) + # lower or left leg + idx += 1 + F = First[site] + LinkList[idx] = F + if qmc_state isa BinaryGroundState || !iszero(F) + LinkList[F] = idx # completes backwards link + else + Last[site] = idx + end + + LegType[idx] = spin_prop[site] + Associates[idx] = 0 + if H isa AbstractLTFIM + op_indices[idx] = n + leg_sites[idx] = 1 + end + + if !isdiagonal(H, op) # off-diagonal site operator + spin_prop[site] ⊻= 1 # spinflip + end + + # upper or right leg + idx += 1 + First[site] = idx + LegType[idx] = spin_prop[site] + Associates[idx] = 0 + if H isa AbstractLTFIM + op_indices[idx] = n + leg_sites[idx] = 1 + end + elseif qmc_state isa BinaryGroundState || isbondoperator(H, op) # diagonal bond operator + site1, site2 = bond = getbondsites(H, op) + spins = spin_prop[site1], spin_prop[site2] + num_sites = 2 # length(bond) + num_legs = 4 # 2*num_sites + + # vertex leg indices + # thinking of imaginary time as increasing as we move upward, + # these indices refer to the + # lower left, lower right, upper left, upper right legs respectively + # v1, v2, v3, v4 = idx + 1, idx + 2, idx + 3, idx + 4 + + @simd for i in 1:num_sites + v = idx + i + st = bond[i] + F = First[st] + LinkList[v] = F + if qmc_state isa BinaryGroundState || !iszero(F) + LinkList[F] = v # completes backwards link + else + Last[st] = v + end + First[st] = v + num_sites + end + + @simd for i in 1:num_legs + v = idx + i + m = mod(i, 1:num_sites) + LegType[v] = spins[m] + Associates[v] = v + 1 + if H isa AbstractLTFIM + op_indices[v] = n + leg_sites[v] = m + end + end + Associates[idx + num_legs] = idx + 1 + idx += num_legs + end + end + + if qmc_state isa BinaryGroundState + # The last N elements of the linked list are the final spin state + @inbounds for i in 1:Ns + idx += 1 + F = First[i] + LinkList[idx] = F + LinkList[F] = idx + LegType[idx] = spin_prop[i] + Associates[idx] = 0 + if H isa AbstractLTFIM + op_indices[idx] = 0 + leg_sites[idx] = 0 + end + end + else + #Periodic boundary conditions for finite-beta + @inbounds for i in 1:Ns + F = First[i] + if !iszero(F) #This might be encountered at high temperatures + L = Last[i] + LinkList[F] = L + LinkList[L] = F + end + end + end + # @debug statements are not run unless debug logging is enabled + @debug("Link List basis state propagation status: $(spin_prop == qmc_state.right_config)", + spin_prop, + qmc_state.right_config) + + return idx +end +link_list_update!(qmc_state, H, d::Diagnostics) = + link_list_update!(Random.GLOBAL_RNG, qmc_state, H, d) + + +############################################################################# + + +@inline function _map_back_operator_list!(ocount::Int, qmc_state::BinaryQMCState, H::AbstractIsing, d::Diagnostics) + operator_list = qmc_state.operator_list + LegType = qmc_state.leg_types + + # if we build an array that maps n to ocount, this loop will become + # very easy to parallelize + @inbounds for (n, op) in enumerate(operator_list) + if isbondoperator(H, op) + if H isa AbstractLTFIM + s1, s2 = LegType[ocount], LegType[ocount+1] + t = getbondtype(H, s1, s2) + operator_list[n] = convertoperatortype(H, op, t) + fit!(d.tmatrix, op, operator_list[n]) + end + ocount += 4 + elseif issiteoperator(H, op) + if LegType[ocount] == LegType[ocount+1] # diagonal + operator_list[n] = makediagonalsiteop(H, getsite(H, op)) + else # off-diagonal + operator_list[n] = makeoffdiagonalsiteop(H, getsite(H, op)) + end + fit!(d.tmatrix, op, operator_list[n]) + ocount += 2 + end + end +end + + +@inline function _map_back_basis_states!(rng::AbstractRNG, lsize::Int, qmc_state::BinaryGroundState, H::AbstractIsing) + Ns = nspins(H) + spin_left, spin_right = qmc_state.left_config, qmc_state.right_config + LegType = qmc_state.leg_types + + @inbounds for i in 1:Ns + spin_left[i] = LegType[i] # left basis state + spin_right[i] = LegType[lsize-Ns+i] # right basis state + end + return Ns + 1 # next one is leg Ns + 1 +end + +@inline function _map_back_basis_states!(rng::AbstractRNG, lsize::Int, qmc_state::BinaryThermalState, H::AbstractIsing) + Ns = nspins(H) + spin_left, spin_right = qmc_state.left_config, qmc_state.right_config + LegType = qmc_state.leg_types + + First = qmc_state.first + Last = qmc_state.last + @inbounds for i in 1:Ns + F = First[i] + if !iszero(F) + spin_left[i] = LegType[Last[i]] # left basis state + spin_right[i] = LegType[F] # right basis state + else + #randomly flip spins not connected to operators + spin_left[i] = spin_right[i] = rand(rng, Bool) + end + end + return 1 # first leg +end + + +function trialstate_weight_change(qmc_state::BinaryGroundState, lsize::Int, Ns::Int, i::Int) # This we can get rid of because it's based on TrialState, only returns logweights -> before cutting this, just let it return 0 + if !(qmc_state.trialstate isa AbstractProductState) + if i <= Ns + push!(qmc_state.trialstate.left_flips, i) + elseif i > (lsize - Ns) + push!(qmc_state.trialstate.right_flips, i - lsize + Ns) + end + return 0.0 + else + if i <= Ns + # return logweightchange(qmc_state.trialstate, qmc_state.left_config[i]) + return 0.0 + elseif i > (lsize - Ns) + # return logweightchange(qmc_state.trialstate, qmc_state.right_config[i]) + return 0.0 + else + return 0.0 + end + end +end +trialstate_weight_change(qmc_state::BinaryThermalState, lsize::Int, Ns::Int, i::Int) = 0.0 + +############################################################################# + +# changed the function name below from cluster_update!() to generic_cluster_update!() in order to distinguish it from the other cluster_update function defined at the very bottom + +function generic_cluster_update!(rng::AbstractRNG, update_kernel!::Function, acceptance::Function, lsize::Int, qmc_state::BinaryQMCState, H::AbstractIsing, d::Diagnostics) + Ns = nspins(H) + operator_list = qmc_state.operator_list + + LinkList = qmc_state.linked_list + LegType = qmc_state.leg_types + Associates = qmc_state.associates + op_indices = qmc_state.op_indices + + in_cluster = fill!(qmc_state.in_cluster, 0) + cstack = qmc_state.cstack # This is the stack of vertices in a cluster + current_cluster = qmc_state.current_cluster + runstats = d.runstats + + ccount = 0 # cluster number counter + + num_accept = 0 + num_reject = 0 + + @inbounds for i in 1:lsize + # Add a new leg onto the cluster + if iszero(in_cluster[i]) && iszero(Associates[i]) + ccount += 1 + push!(cstack, i) + in_cluster[i] = ccount + if qmc_state isa BinaryGroundState && !(qmc_state.trialstate isa AbstractProductState) + left_flips = empty!(qmc_state.trialstate.left_flips) + right_flips = empty!(qmc_state.trialstate.right_flips) + end + + empty!(current_cluster) + push!(current_cluster, i) + lnA = 0.0 + if qmc_state isa BinaryGroundState + lnA += trialstate_weight_change(qmc_state, lsize, Ns, i) + end + + while !isempty(cstack) + leg = LinkList[pop!(cstack)] + + if iszero(in_cluster[leg]) + in_cluster[leg] = ccount # add the new leg and flip it + push!(current_cluster, leg) + if qmc_state isa BinaryGroundState + lnA += trialstate_weight_change(qmc_state, lsize, Ns, i) + end + a = Associates[leg] + + a == 0 && continue + # from this point on, we know we're on a bond op + op = operator_list[op_indices[leg]] + w = getweightindex(H, op) - getoperatortype(H, op) + preflip_bond_type, postflip_bond_type = update_kernel!(qmc_state, H, ccount, leg, a) + lnA += ( + getlogweight(H.op_sampler, w + postflip_bond_type) + - getlogweight(H.op_sampler, w + preflip_bond_type) + ) + end + end + + if qmc_state isa BinaryGroundState && !(qmc_state.trialstate isa AbstractProductState) + lnA += logweightchange(qmc_state.trialstate, qmc_state.left_config, left_flips) + lnA += logweightchange(qmc_state.trialstate, qmc_state.right_config, right_flips) + end + + A = acceptance(H, lnA) + fit!(runstats, :cluster_update_accept, A) + + if rand(rng) < A + @inbounds for j in current_cluster + LegType[j] ⊻= 1 # spinflip + end + fit!(runstats, :accepted_cluster_sizes, length(current_cluster)) + num_accept += 1 + else + fit!(runstats, :rejected_cluster_sizes, length(current_cluster)) + num_reject += 1 + end + fit!(runstats, :cluster_sizes, length(current_cluster)) + end + end + + if !(runstats isa NoStats) + fit!(runstats, :accepted_cluster_count, num_accept+1) + fit!(runstats, :rejected_cluster_count, num_reject+1) + fit!(runstats, :cluster_count, ccount+1) + end + + # map back basis states and operator list + ocount = _map_back_basis_states!(rng, lsize, qmc_state, H) + _map_back_operator_list!(ocount, qmc_state, H, d) + + return lsize +end + +# might put this cluster_update function in a more visible position somewhere + +function cluster_update!(rng, qmc_state, H::Hamiltonian, d::Diagnostics; p::Float64=0.0, kw...) + if rand(rng) < p + multibranch_update!(rng, qmc_state, H, d) + else + line_update!(rng, qmc_state, H, d) + end +end diff --git a/lib/BloqadeQMC/src/updates/diagonal.jl b/lib/BloqadeQMC/src/updates/diagonal.jl new file mode 100644 index 000000000..342e7fc7e --- /dev/null +++ b/lib/BloqadeQMC/src/updates/diagonal.jl @@ -0,0 +1,135 @@ +# deleted alignment_check and insert_diagonal_operator for TFIM case + +function insert_diagonal_operator!(rng::AbstractRNG, qmc_state::BinaryQMCState{K, V}, H::AbstractLTFIM{<:AbstractImprovedOperatorSampler{K, T}}, spin_prop::V, n::Int) where {K, T, V} +# Gets called in full_diagonal_update() below + op = rand(rng, H.op_sampler) + + @inbounds if issiteoperator(H, op) + qmc_state.operator_list[n] = op + return op, zero(T) + else + t = getoperatortype(H, op) + lw1 = getlogweight(H, op) + site1, site2 = getbondsites(H, op) + real_t = getbondtype(H, spin_prop[site1], spin_prop[site2]) + if t == real_t + qmc_state.operator_list[n] = op + return op, lw1 + end + + op2 = convertoperatortype(H, op, real_t) + lw2 = getlogweight(H, op2) + + if rand(rng) < exp(lw2 - lw1) + qmc_state.operator_list[n] = op2 + return op2, lw2 + else + return nothing, zero(lw2) + end + end +end + +# Diagonal update from groundstate.jl + +function full_diagonal_update!(rng::AbstractRNG, qmc_state::BinaryGroundState, H::AbstractIsing, d::Diagnostics) +# gets called in mc_step() in mc_step.jl + spin_prop = copyto!(qmc_state.propagated_config, qmc_state.left_config) # the propagated spin state + runstats = d.runstats + + if !(runstats isa NoStats) # Technically, this runstats stuff could be cut, but it is useful for checking MC on developer side + failures = 0 + count = 0 + end + + for (n, op) in enumerate(qmc_state.operator_list) + if !isdiagonal(H, op) + @inbounds spin_prop[getsite(H, op)] ⊻= 1 # spinflip + else + op = nothing + if !(runstats isa NoStats); i = -1; end + while op === nothing # here we are attempting step 3 until we achieve a successful insertion + op, _ = insert_diagonal_operator!(rng, qmc_state, H, spin_prop, n) + if !(runstats isa NoStats); i += 1; end + end + if !(runstats isa NoStats); failures += i; count += 1; end + end + end + + # @debug statements are not run unless debug logging is enabled + @debug("Diagonal Update basis state propagation status: $(spin_prop == qmc_state.right_config)", + spin_prop, + qmc_state.right_config) + + if !(runstats isa NoStats) + fit!(runstats, :diag_update_fails, failures/count) + end +end + +######################################################################################################################################### + +# Diagonal updates (full_diagonal_update_beta and resize_op_list from mixedstate.jl) + +function resize_op_list!(qmc_state::BinaryThermalState{K}, H::AbstractIsing, new_size::Int) where {K} + operator_list = filter!(!isidentity(H), qmc_state.operator_list) + len = length(operator_list) + + if len < new_size + tail = init_op_list(new_size - len, K) + append!(operator_list, tail) + end + + len = 4*length(operator_list) + # these are going to be overwritten by the cluster update which will be + # called right after the diagonal update that called this function + resize!(qmc_state.linked_list, len) + resize!(qmc_state.leg_types, len) + resize!(qmc_state.associates, len) + resize!(qmc_state.leg_sites, len) + resize!(qmc_state.op_indices, len) + resize!(qmc_state.in_cluster, len) +end + + +function full_diagonal_update_beta!(rng::AbstractRNG, qmc_state::BinaryThermalState, H::AbstractIsing, beta::Real; eq::Bool = false) # eq = true must be set during eq phase for tuning M (resizing) +# gets called in mc_step() in mc_step.jl + # As we can see: thermal case does not call runstats + P_norm = beta * diag_update_normalization(H) + + num_ids = count(isidentity(H), qmc_state.operator_list) + + spin_prop = copyto!(qmc_state.propagated_config, qmc_state.left_config) # the propagated spin state + + @inbounds for (n, op) in enumerate(qmc_state.operator_list) + if !isdiagonal(H, op) # Step 4 + spin_prop[getsite(H, op)] ⊻= 1 # spinflip + elseif !isidentity(H, op) + if rand(rng)*P_norm < (num_ids + 1) # implements condition (16) + qmc_state.operator_list[n] = makeidentity(H) # This corresponds to removing the operator in Step 1 + num_ids += 1 + end + else + if rand(rng)*num_ids < P_norm # implements condition (17) + op, _ = insert_diagonal_operator!(rng, qmc_state, H, spin_prop, n) # This corresponds to inserting an operator in Step 2. + # Sampling step 3 is implemented in insert_diagonal_operator!() + if op !== nothing + num_ids -= 1 + end + end + end + end + + # DEBUG + # if spin_prop != qmc_state.right_config # check the spin propagation for error + # error("Basis state propagation error in diagonal update!") + # end + + total_list_size = length(qmc_state.operator_list) + num_ops = total_list_size - num_ids + + if eq && 1.2*num_ops > total_list_size + resize_op_list!(qmc_state, H, round(Int, 1.5*num_ops, RoundUp)) # Only happens in equilibration phase + end + + return num_ops +end + diff --git a/lib/BloqadeQMC/src/updates/line.jl b/lib/BloqadeQMC/src/updates/line.jl new file mode 100644 index 000000000..5b0463cd3 --- /dev/null +++ b/lib/BloqadeQMC/src/updates/line.jl @@ -0,0 +1,47 @@ +# TODO: there might be a slight performance degradation with this +# functional approach, might wanna try to inject this function directly +# into the body of the cluster update + +@inline function line_kernel!(qmc_state::BinaryQMCState, H::AbstractIsing, ccount::Int, leg::Int, a::Int) + Ns = nspins(H) + LegType, Associates, leg_sites = qmc_state.leg_types, qmc_state.associates, qmc_state.leg_sites + in_cluster, cstack, current_cluster = qmc_state.in_cluster, qmc_state.cstack, qmc_state.current_cluster + + @inbounds ll, la = LegType[leg], LegType[a] + @inbounds sl, sa = leg_sites[leg], leg_sites[a] + # TODO: check if this is inputting the spins in the correct order + if sl > sa + preflip_bond_type = getbondtype(H, ll, la) + postflip_bond_type = getbondtype(H, !ll, la) + else + preflip_bond_type = getbondtype(H, la, ll) + postflip_bond_type = getbondtype(H, la, !ll) + end + # ^using circshift for this is too slow, gonna have to manually expand out + # a bunch of ifs for k-local (might need to do some metaprogramming!) + # using SVectors and sortperm should work (it's super fast!) and is easier + # also, should short-circuit if order doesn't matter, i.e. + # the longitudinal field is uniform + coordination numbers are all equal + + # now add the straight-through leg to the cluster + @inbounds straight_thru = Associates[a] + @inbounds in_cluster[straight_thru] = ccount + push!(cstack, straight_thru) + push!(current_cluster, straight_thru) + + return preflip_bond_type, postflip_bond_type +end + +line_acceptance(H::AbstractIsing, lnA::T) where {T <: Real} = exp(min(lnA, zero(T))) + +function line_cluster_update!(rng::AbstractRNG, lsize::Int, qmc_state::BinaryQMCState, H::AbstractIsing, d::Diagnostics) + return generic_cluster_update!(rng, line_kernel!, line_acceptance, lsize, qmc_state, H, d) # call cluster_update with line kernel +end + +############################################################################# + +function line_update!(rng::AbstractRNG, qmc_state::BinaryQMCState, H::AbstractIsing, d::Diagnostics) + lsize = link_list_update!(rng, qmc_state, H, d) + return line_cluster_update!(rng, lsize, qmc_state, H, d) +end + diff --git a/lib/BloqadeQMC/src/updates/mc_step.jl b/lib/BloqadeQMC/src/updates/mc_step.jl new file mode 100644 index 000000000..40f3a3bcd --- /dev/null +++ b/lib/BloqadeQMC/src/updates/mc_step.jl @@ -0,0 +1,36 @@ +# Todo: Merge the following two mc_step functions + +# mc_step from groundstate.jl + +# I am keeping this function name cluster_update!() to keep it separate from the generic_cluster_update() function in cluster.jl. Previously, the two functions had the same name. + +function cluster_update!(rng, qmc_state, H::Hamiltonian, d::Diagnostics; p::Float64=0.0, kw...) + if rand(rng) < p + multibranch_update!(rng, qmc_state, H, d) + else + line_update!(rng, qmc_state, H, d) + end +end + +function mc_step!(f::Function, rng::AbstractRNG, qmc_state::BinaryGroundState, H::Hamiltonian, d::Diagnostics; kw...) + full_diagonal_update!(rng, qmc_state, H, d) + lsize = cluster_update!(rng, qmc_state, H, d; kw...) + f(lsize, qmc_state, H) +end +# mc_step!(f::Function, qmc_state::BinaryGroundState, H::Hamiltonian, d::Diagnostics; kw...) = mc_step!(f, Random.GLOBAL_RNG, qmc_state, H, d; kw...) +mc_step!(rng::AbstractRNG, qmc_state::BinaryGroundState, H::Hamiltonian, d::Diagnostics; kw...) = mc_step!((args...) -> nothing, rng, qmc_state, H, d; kw...) +# mc_step!(qmc_state::BinaryGroundState, H::Hamiltonian, d::Diagnostics; kw...) = mc_step!(Random.GLOBAL_RNG, qmc_state, H, d; kw...) + + +# mc_step from mixedstate.jl + +function mc_step_beta!(f::Function, rng::AbstractRNG, qmc_state::BinaryThermalState, H::AbstractIsing, beta::Real, d::Diagnostics; eq::Bool = false, kw...) + num_ops = full_diagonal_update_beta!(rng, qmc_state, H, beta; eq=eq) + lsize = cluster_update!(rng, qmc_state, H, d; kw...) + f(lsize, qmc_state, H) + return num_ops +end + +# mc_step_beta!(f::Function, qmc_state, H, beta, d::Diagnostics; eq = false, kw...) = mc_step_beta!(f, Random.GLOBAL_RNG, qmc_state, H, beta, d; eq = eq, kw...) +mc_step_beta!(rng::AbstractRNG, qmc_state, H, beta, d::Diagnostics; eq = false, kw...) = mc_step_beta!((args...) -> nothing, rng, qmc_state, H, beta, d; eq = eq, kw...) +# mc_step_beta!(qmc_state, H, beta, d::Diagnostics; eq = false, kw...) = mc_step_beta!(Random.GLOBAL_RNG, qmc_state, H, beta, d; eq = eq, kw...) \ No newline at end of file diff --git a/lib/BloqadeQMC/src/updates/multibranch.jl b/lib/BloqadeQMC/src/updates/multibranch.jl new file mode 100644 index 000000000..8f6e95d3c --- /dev/null +++ b/lib/BloqadeQMC/src/updates/multibranch.jl @@ -0,0 +1,47 @@ +@inline function multibranch_kernel!(qmc_state::BinaryQMCState, H::AbstractIsing, ccount::Int, leg::Int, a::Int) + Ns = nspins(H) + LegType, Associates, leg_sites = qmc_state.leg_types, qmc_state.associates, qmc_state.leg_sites + in_cluster, cstack, current_cluster = qmc_state.in_cluster, qmc_state.cstack, qmc_state.current_cluster + + @inbounds ll, la = LegType[leg], LegType[a] + @inbounds sl, sa = leg_sites[leg], leg_sites[a] + # TODO: check if this is inputting the spins in the correct order + if sl > sa + preflip_bond_type = getbondtype(H, ll, la) + postflip_bond_type = getbondtype(H, !ll, !la) + else + preflip_bond_type = getbondtype(H, la, ll) + postflip_bond_type = getbondtype(H, !la, !ll) + end + + # now check all associates and add them to the cluster + @inbounds while a != 0 && iszero(in_cluster[a]) + push!(cstack, a) + in_cluster[a] = ccount + push!(current_cluster, a) + a = Associates[a] + end + + return preflip_bond_type, postflip_bond_type +end + +multibranch_acceptance(H::AbstractIsing, lnA::T) where {T <: Real} = + haslongitudinalfield(H) ? exp(min(lnA, zero(lnA))) : T(0.5) +multibranch_acceptance(H::AbstractRydberg, lnA::T) where {T <: Real} = exp(min(lnA, zero(T))) +# in the TFIM case, acceptance rate is exactly 1 so we set it to 1/2 to ensure ergodicity +multibranch_acceptance(H::AbstractTFIM, lnA::T) where {T <: Real} = T(0.5) + + +function multibranch_cluster_update!(rng::AbstractRNG, lsize::Int, qmc_state::BinaryQMCState, H::AbstractIsing, d::Diagnostics) + return generic_cluster_update!(rng, multibranch_kernel!, multibranch_acceptance, lsize, qmc_state, H, d) # call cluster_update with the proper kernel +end + +# Here, I cut the simplified implementation for TFIM + +############################################################################# + + +function multibranch_update!(rng::AbstractRNG, qmc_state::BinaryQMCState, H::AbstractIsing, d::Diagnostics) + lsize = link_list_update!(rng, qmc_state, H, d) + return multibranch_cluster_update!(rng, lsize, qmc_state, H, d) +end diff --git a/lib/BloqadeQMC/src/updates/updates.jl b/lib/BloqadeQMC/src/updates/updates.jl new file mode 100644 index 000000000..639506ee3 --- /dev/null +++ b/lib/BloqadeQMC/src/updates/updates.jl @@ -0,0 +1,5 @@ +include("cluster.jl") +include("diagonal.jl") +include("line.jl") +include("mc_step.jl") +include("multibranch.jl") \ No newline at end of file diff --git a/lib/BloqadeQMC/test/ltfim.jl b/lib/BloqadeQMC/test/ltfim.jl index 07fddd3f1..46b99f540 100644 --- a/lib/BloqadeQMC/test/ltfim.jl +++ b/lib/BloqadeQMC/test/ltfim.jl @@ -3,7 +3,9 @@ using Statistics using Random using RandomNumbers using Measurements - +# using Plots: histogram, display, savefig +using BinningAnalysis +using BloqadeQMC: jackknife expected_values = Dict{Tuple{Bool, Int, Float64}, Dict{String, Float64}}( (false, 10, Inf64) => Dict("M" => 0.9300787137837009, @@ -94,14 +96,16 @@ expected_values = Dict{Tuple{Bool, Int, Float64}, Dict{String, Float64}}( # @test stdscore(energy, expected_vals["H"]) < THRESHOLD # end +THRESHOLD = 2.576 # 99% Two-sided CI of the t-distribution with infinite dofs -@testset "1D LTFIM Thermal State $N-sites, PBC=$PBC, β=1.0" for PBC in [true, false], N in [5, 10] +@testset "1D LTFIM Thermal State $N-sites, PBC=$PBC, β=1.0" for PBC in [true], N in [5] # rng = Xorshifts.Xoroshiro128Plus(1234) - rng = MersenneTwister(4321) + rng = MersenneTwister(1432) - H = LTFIM((N,), 1.0, 1.0, 1.0, PBC) + H = LTFIM((N,), 1.0, 0.0, 0.0, PBC) th = BinaryThermalState(H, 1000) beta = 1.0 + d = Diagnostics() MCS = 1_000_000 EQ_MCS = 100_000 @@ -109,26 +113,93 @@ expected_values = Dict{Tuple{Bool, Int, Float64}, Dict{String, Float64}}( mags = zeros(MCS) ns = zeros(MCS) - [mc_step_beta!(rng, th, H, beta) for i in 1:EQ_MCS] + [mc_step_beta!(rng, th, H, beta, d) for i in 1:EQ_MCS] for i in 1:MCS # Monte Carlo Steps - ns[i] = mc_step_beta!(rng, th, H, beta) do lsize, th, H - mags[i] = magnetization(sample(H, th)) + ns[i] = mc_step_beta!(rng, th, H, beta, d) do lsize, th, H + mags[i] = magnetization(sample(H, th, 1)) end end + + # Binning Analysis for observables based on mags + B = LogBinner(abs.(mags)) + println("Correlation time for abs_mag_binned samples") + τ_abs = tau(B) + println(τ_abs) + println("N_eff") + N_eff = MCS / (2 * τ_abs + 1) + println(N_eff) + abs_mag_binned = measurement(mean(B), std_error(B)*sqrt(N_eff)) # converting standard error to std dev + println() + + B2 = LogBinner(abs2.(mags)) + println("Correlation time for mag_sqr_binned samples") + τ_abs2 = tau(B) + println(τ_abs2) + println("N_eff2") + N_eff2 = MCS / (2 * τ_abs2 + 1) + println(N_eff2) + mag_sqr_binned = measurement(mean(B2), std_error(B2)*sqrt(N_eff2)) # converting standard error to std dev + println() + + # Binning analysis for observables based on ns + energy(x) = -x / beta + H.energy_shift + BE = LogBinner(energy.(ns) / nspins(H)) + println("Correlation time for energy_binned samples") + τ_energy = tau(BE) + println(τ_energy) + println("N_eff_E") + N_eff_E = MCS / (2 * τ_energy + 1) + println(N_eff_E) + energy_binned = measurement(mean(BE), std_error(BE)*sqrt(N_eff_E)) # converting standard error to std dev + println() + # Mean looks good, but I'm worried about the std_dev... Huge increase!!! + + # println("Number of binning levels") + # println(length(BE.accumulators)) + # println(has_converged(BE)) + # println("Number of entries in last binning level") + # println(BE.accumulators[20]) + # println() + + # Unbinned calculations + abs_mag = mean_and_stderr(abs, mags) + println("abs_mag") + println(abs_mag) + println("abs_mag_binned") + println(abs_mag_binned) + println() + mag_sqr = mean_and_stderr(abs2, mags) + println("mag_abs") + println(mag_sqr) + println("mag_sqr_binned") + println(mag_sqr_binned) energy = mean_and_stderr(x -> -x/beta, ns) + H.energy_shift energy /= nspins(H) + println("energy") + println(energy) + println("energy_binned") + println(energy_binned) + println() + heat_capacity = jackknife(ns .^ 2, ns) do nsqr, n nsqr - n^2 - n end expected_vals = expected_values[(PBC, N, 1.0)] - @test abs(stdscore(abs_mag, expected_vals["|M|"])) < THRESHOLD - @test abs(stdscore(mag_sqr, expected_vals["M^2"])) < THRESHOLD - @test abs(stdscore(energy, expected_vals["H"])) < THRESHOLD - @test abs(stdscore(heat_capacity, expected_vals["C"])) < THRESHOLD + # @test abs(stdscore(abs_mag, expected_vals["|M|"])) < THRESHOLD + @test abs(stdscore(abs_mag_binned, expected_vals["|M|"])) < THRESHOLD + # @test abs(stdscore(mag_sqr, expected_vals["M^2"])) < THRESHOLD + @test abs(stdscore(mag_sqr_binned, expected_vals["M^2"])) < THRESHOLD + # @test abs(stdscore(energy, expected_vals["H"])) < THRESHOLD + @test abs(stdscore(energy_binned, expected_vals["H"])) < THRESHOLD + # @test abs(stdscore(heat_capacity, expected_vals["C"])) < THRESHOLD + + # histJlogPBC = histogram(mags, bins = :scott, weights = repeat(1:5000, outer = 200), yscale=:log10) + # savefig(histJlogPBC,"histJlogPBC.png") + # display(hist) end diff --git a/lib/BloqadeQMC/test/runtests.jl b/lib/BloqadeQMC/test/runtests.jl index f1fcb408b..a572ddf83 100644 --- a/lib/BloqadeQMC/test/runtests.jl +++ b/lib/BloqadeQMC/test/runtests.jl @@ -24,6 +24,6 @@ H = @inferred LTFIM((10,), 1.0, 1.0, 1.0) THRESHOLD = 2.576 # 99% Two-sided CI of the t-distribution with infinite dofs - -# include("ltfim.jl") +include("error.jl") +include("ltfim.jl") # include("tfim.jl")