From 8d6a10c126afe33284cf77958a717acbd31486e4 Mon Sep 17 00:00:00 2001 From: Anna Knoerr Date: Mon, 30 May 2022 23:31:01 +0200 Subject: [PATCH 01/21] Test file --- Test/test.jl | 1 + 1 file changed, 1 insertion(+) create mode 100644 Test/test.jl diff --git a/Test/test.jl b/Test/test.jl new file mode 100644 index 000000000..d6459e005 --- /dev/null +++ b/Test/test.jl @@ -0,0 +1 @@ +xxx From b6deabf6f67a048fdeaf568dcfecc1b45fb90491 Mon Sep 17 00:00:00 2001 From: Anna Knoerr Date: Mon, 5 Sep 2022 09:56:53 +0200 Subject: [PATCH 02/21] updated dependencies for testing jackknife function --- lib/BloqadeQMC/Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/lib/BloqadeQMC/Project.toml b/lib/BloqadeQMC/Project.toml index 1161a316f..640e91988 100644 --- a/lib/BloqadeQMC/Project.toml +++ b/lib/BloqadeQMC/Project.toml @@ -10,6 +10,7 @@ Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" +Jackknife = "133c4774-c510-5f25-b650-129312f89f69" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7" OnlineStats = "a15396b6-48d5-5d58-9928-6d29437db91e" From 95c4e9dbb323933d4789b94fdb82cb94d98a353f Mon Sep 17 00:00:00 2001 From: Anna Knoerr Date: Mon, 5 Sep 2022 09:58:39 +0200 Subject: [PATCH 03/21] added jackknife unit tests --- lib/BloqadeQMC/test/error.jl | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) create mode 100644 lib/BloqadeQMC/test/error.jl diff --git a/lib/BloqadeQMC/test/error.jl b/lib/BloqadeQMC/test/error.jl new file mode 100644 index 000000000..5def0ef2c --- /dev/null +++ b/lib/BloqadeQMC/test/error.jl @@ -0,0 +1,23 @@ +using Test +using Statistics: mean, std +using Measurements: value, uncertainty +# using BloqadeQMC: jackknife +using Jackknife: estimate, variance, bias + +include("../src/error.jl") + +@testset "simple numerical tests for jackknife" begin + x = [3,2,11,7,6] + + # Resampling should not change estimate when f = mean. + m = mean(x) + @test value(jackknife(mean, x)) ≈ m + + # Resampling should give a better estimate for e.g. f = inv. We check this using the Julia library Jackknife + g(x) = 1/mean(x) # The estimate function from Jackknife takes needs g to be explicity defined as a function of the mean + @test value(jackknife(inv, x)) ≈ estimate(g, x) + @test uncertainty(jackknife(inv, x)) ≈ sqrt(variance(g, x)) + + # Test jackknife when dispatched on f = identity + @test std(x) > uncertainty(jackknife(x)) +end From d09626828931184c244522737871ca31ae4e45de Mon Sep 17 00:00:00 2001 From: Anna Knoerr Date: Mon, 5 Sep 2022 10:14:39 +0200 Subject: [PATCH 04/21] updated dependencies to include Jackknife --- lib/BloqadeQMC/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/BloqadeQMC/Project.toml b/lib/BloqadeQMC/Project.toml index 640e91988..1d1671b35 100644 --- a/lib/BloqadeQMC/Project.toml +++ b/lib/BloqadeQMC/Project.toml @@ -10,7 +10,6 @@ Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" -Jackknife = "133c4774-c510-5f25-b650-129312f89f69" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7" OnlineStats = "a15396b6-48d5-5d58-9928-6d29437db91e" @@ -21,6 +20,7 @@ RandomNumbers = "e6cf234a-135c-5ec9-84dd-332b85af5143" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +Jackknife = "133c4774-c510-5f25-b650-129312f89f69" [compat] BinningAnalysis = "0.5" From e72c4f16812658b28ae721b3696a5a836f555928 Mon Sep 17 00:00:00 2001 From: Anna Knoerr Date: Mon, 5 Sep 2022 10:15:05 +0200 Subject: [PATCH 05/21] unit tests for jackknife resampling --- lib/BloqadeQMC/test/error.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/lib/BloqadeQMC/test/error.jl b/lib/BloqadeQMC/test/error.jl index 5def0ef2c..2d29128e5 100644 --- a/lib/BloqadeQMC/test/error.jl +++ b/lib/BloqadeQMC/test/error.jl @@ -1,11 +1,9 @@ -using Test +using Test, BloqadeQMC using Statistics: mean, std using Measurements: value, uncertainty -# using BloqadeQMC: jackknife +using BloqadeQMC: jackknife using Jackknife: estimate, variance, bias -include("../src/error.jl") - @testset "simple numerical tests for jackknife" begin x = [3,2,11,7,6] From 12261cbf562487be7f36a5c4a9b57d564d82587a Mon Sep 17 00:00:00 2001 From: Phillip Weinberg Date: Tue, 6 Sep 2022 09:55:00 -0400 Subject: [PATCH 06/21] Delete test.jl removing extra file not needed. --- Test/test.jl | 1 - 1 file changed, 1 deletion(-) delete mode 100644 Test/test.jl diff --git a/Test/test.jl b/Test/test.jl deleted file mode 100644 index d6459e005..000000000 --- a/Test/test.jl +++ /dev/null @@ -1 +0,0 @@ -xxx From 034f8f42c962481884b349cbe6491e43ca54c157 Mon Sep 17 00:00:00 2001 From: Anna Knoerr Date: Tue, 6 Sep 2022 19:38:15 +0200 Subject: [PATCH 07/21] updated comments defining the TFIM ops to match the code --- lib/BloqadeQMC/src/ising/TFIM.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/BloqadeQMC/src/ising/TFIM.jl b/lib/BloqadeQMC/src/ising/TFIM.jl index 465bdafb3..1293d23f8 100644 --- a/lib/BloqadeQMC/src/ising/TFIM.jl +++ b/lib/BloqadeQMC/src/ising/TFIM.jl @@ -24,8 +24,8 @@ const ISING_OP_SIZE = 5 ############################################################################### # TFIM ops: -# (1,2,1,0,i) is an off-diagonal site operator h*sigma^x_i -# (1,1,1,0,i) is a diagonal site operator h +# (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) @inline getoperatorlocality(::Type{<:AbstractIsing}, op::NTuple{ISING_OP_SIZE, Int}) = @inbounds op[1] From 7de44c7532b39218a496c1eaaaa220c68bc231cb Mon Sep 17 00:00:00 2001 From: Anna Knoerr Date: Mon, 19 Sep 2022 16:43:43 +0200 Subject: [PATCH 08/21] more jackknife tests --- lib/BloqadeQMC/test/error.jl | 29 +++++++++++++++++++++++------ 1 file changed, 23 insertions(+), 6 deletions(-) diff --git a/lib/BloqadeQMC/test/error.jl b/lib/BloqadeQMC/test/error.jl index 2d29128e5..541de8ee6 100644 --- a/lib/BloqadeQMC/test/error.jl +++ b/lib/BloqadeQMC/test/error.jl @@ -2,7 +2,10 @@ using Test, BloqadeQMC using Statistics: mean, std using Measurements: value, uncertainty using BloqadeQMC: jackknife -using Jackknife: estimate, variance, bias +using Random +using RandomNumbers + +rng = Xorshifts.Xoroshiro128Plus(1234) @testset "simple numerical tests for jackknife" begin x = [3,2,11,7,6] @@ -12,10 +15,24 @@ using Jackknife: estimate, variance, bias @test value(jackknife(mean, x)) ≈ m # Resampling should give a better estimate for e.g. f = inv. We check this using the Julia library Jackknife - g(x) = 1/mean(x) # The estimate function from Jackknife takes needs g to be explicity defined as a function of the mean - @test value(jackknife(inv, x)) ≈ estimate(g, x) - @test uncertainty(jackknife(inv, x)) ≈ sqrt(variance(g, x)) + @test value(jackknife(inv, x)) ≈ 0.158110765 + @test uncertainty(jackknife(inv, x)) ≈ 0.052468798 +end + +@testset "functional tests for jackknife" begin + + mean_val = 5.0 + x = fill(mean_val, 100) + randn!(rng, zeros(100)) + j_error = uncertainty(jackknife(x)) + #check jackknifed std < original sample std + @test j_error < std(x) + + add(x,y,z) = x+y+z + num_measurements_list = [10,100,1000,10000,100000] + #create a set of normally random values about a mean value for x,y,z for different measuremensts + xs_list = [[fill(mean_val, num_measurements_list[n]) + randn!(rng, zeros(num_measurements_list[n])) for v in 1:3] for n in 1:length(num_measurements_list)] + jackknife_errors = [uncertainty(jackknife(add, xs_list[n][1], xs_list[n][2], xs_list[n][3])) for n in 1:length(num_measurements_list)] + #test that error decreases as the number of observations gets larger + @test issorted(jackknife_errors, lt=Base.isgreater) - # Test jackknife when dispatched on f = identity - @test std(x) > uncertainty(jackknife(x)) end From 1828614f17f3ee87d372233bd5c24eb1d427b225 Mon Sep 17 00:00:00 2001 From: Anna Knoerr Date: Tue, 27 Sep 2022 21:15:54 +0200 Subject: [PATCH 09/21] removed Jackknife dependency --- lib/BloqadeQMC/Project.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/lib/BloqadeQMC/Project.toml b/lib/BloqadeQMC/Project.toml index 1d1671b35..1161a316f 100644 --- a/lib/BloqadeQMC/Project.toml +++ b/lib/BloqadeQMC/Project.toml @@ -20,7 +20,6 @@ RandomNumbers = "e6cf234a-135c-5ec9-84dd-332b85af5143" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -Jackknife = "133c4774-c510-5f25-b650-129312f89f69" [compat] BinningAnalysis = "0.5" From 7c8c851c066ec92914c915c11b50d3a62b41e89e Mon Sep 17 00:00:00 2001 From: Anna Knoerr Date: Thu, 29 Sep 2022 01:36:10 +0200 Subject: [PATCH 10/21] t-test with binning for single QMC run --- lib/BloqadeQMC/test/ltfim.jl | 86 +++++++++++++++++++++++++++++++----- 1 file changed, 75 insertions(+), 11 deletions(-) diff --git a/lib/BloqadeQMC/test/ltfim.jl b/lib/BloqadeQMC/test/ltfim.jl index 07fddd3f1..35ffae412 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,86 @@ 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 / τ_abs + 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 / τ_abs2 + 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 / τ_energy + println(N_eff_E) + energy_binned = measurement(mean(BE), std_error(BE)*sqrt(N_eff_E)) # converting standard error to std dev + println() + + + + 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) + + 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 From 20a14478a45c1275ea471601e3e2db8678465edb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Anna=20Kn=C3=B6rr?= <61286790+Atomyka@users.noreply.github.com> Date: Thu, 29 Sep 2022 01:59:00 +0200 Subject: [PATCH 11/21] Delete TFIM.jl --- lib/BloqadeQMC/src/ising/TFIM.jl | 172 ------------------------------- 1 file changed, 172 deletions(-) delete mode 100644 lib/BloqadeQMC/src/ising/TFIM.jl diff --git a/lib/BloqadeQMC/src/ising/TFIM.jl b/lib/BloqadeQMC/src/ising/TFIM.jl deleted file mode 100644 index 1293d23f8..000000000 --- a/lib/BloqadeQMC/src/ising/TFIM.jl +++ /dev/null @@ -1,172 +0,0 @@ -abstract type AbstractIsing{O <: AbstractOperatorSampler} <: Hamiltonian{2,O} end -abstract type AbstractTFIM{O <: AbstractOperatorSampler} <: AbstractIsing{O} end - -# H = -J ∑_{⟨ij⟩} σ_i σ_j - h ∑_i σ^x_i -struct NearestNeighbourTFIM{O} <: AbstractTFIM{O} - op_sampler::O - J::Float64 - hx::Float64 - Ns::Int - energy_shift::Float64 -end - - -# H = ∑_{ij} J_ij σ_i σ_j - ∑_i h_i σ^x_i -struct TFIM{O,M <: UpperTriangular{Float64},V <: AbstractVector{Float64}} <: AbstractTFIM{O} - op_sampler::O - J::M - hx::V - Ns::Int - energy_shift::Float64 -end - -const ISING_OP_SIZE = 5 -############################################################################### - -# 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) -@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)) - -############################################################################### - -function make_prob_vector(J::UpperTriangular{T}, hx::AbstractVector{T}) where T - @assert length(hx) == size(J, 1) == size(J, 2) - - ops = Vector{NTuple{ISING_OP_SIZE, Int}}(undef, 0) - p = Vector{T}(undef, 0) - energy_shift = zero(T) - - for i in eachindex(hx) - if !iszero(hx[i]) - push!(ops, makediagonalsiteop(AbstractTFIM, i)) - push!(p, hx[i]) - energy_shift += hx[i] - end - end - - # only take J_ij terms from upper triangle - for j in axes(J, 2), i in axes(J, 1) - if i < j # i != j: we don't want self-interactions - if J[i, j] != 0 - push!(p, 2*abs(J[i, j])) - push!(ops, (2, 1, length(p), i, j)) - energy_shift += abs(J[i, j]) - end - end - end - - return ops, p, energy_shift -end - -function make_uniform_tfim(bond_spins::Vector{NTuple{2,Int}}, Ns::Int, J::T, hx::T) where T - hx_ = hx * ones(T, Ns) - J_ = zeros(T, Ns, Ns) - for (i, j) in bond_spins - i, j = (i <= j) ? (i, j) : (j, i) - J_[i, j] = -J - end - - return UpperTriangular(triu!(J_, 1)), hx_ -end - -############################################################################### - -function TFIM(J::UpperTriangular{Float64}, hx::AbstractVector{Float64}) - @assert length(hx) == size(J, 1) == size(J, 2) - - ops, p, energy_shift = make_prob_vector(J, hx) - Ns = length(hx) - op_sampler = OperatorSampler(ops, p) - - return TFIM{typeof(op_sampler), typeof(J), typeof(hx)}( - op_sampler, J, hx, Ns, energy_shift - ) -end - - -function TFIM(bond_spin, Ns::Int, Nb::Int, hx::Float64, J::Float64) - J_, hx_ = make_uniform_tfim(bond_spin, Ns, J, hx) - return TFIM(J_, hx_) -end - -abstract type HXField; end -struct ConstantHX <: HXField; end -struct VaryingHX <: HXField; end - -hxfield(::AbstractIsing) = VaryingHX() -hxfield(::NearestNeighbourTFIM) = ConstantHX() - -total_hx(::ConstantHX, H::AbstractIsing) = nspins(H) * H.hx -total_hx(::VaryingHX, H::AbstractIsing) = sum(H.hx) -total_hx(H::AbstractIsing) = total_hx(hxfield(H), H) - -# returns true if H.J[site1, site2] is negative -Base.@propagate_inbounds isferromagnetic(H::TFIM, (site1, site2)::NTuple{2, Int}) = signbit(H.J[site1, site2]) -haslongitudinalfield(::AbstractTFIM) = false - - -############################################################################### - - -function energy(::Type{<:BinaryGroundState}, H::AbstractIsing, ns::Vector{<:Real}; resampler::Function=jackknife) - hx = total_hx(H) - - if !iszero(hx) - E = -hx * resampler(inv, ns) - else - E = measurement(zero(H.energy_shift)) - end - - return H.energy_shift + E -end - -function energy(::Type{<:BinaryGroundState}, H::AbstractIsing, ns_mean::Real) - hx = total_hx(H) - - if !iszero(hx) - E = -hx / ns_mean - else - E = zero(H.energy_shift) - end - - return H.energy_shift + E -end From 417093745b644952064e867ef952686d4e406a64 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Anna=20Kn=C3=B6rr?= <61286790+Atomyka@users.noreply.github.com> Date: Thu, 29 Sep 2022 01:59:42 +0200 Subject: [PATCH 12/21] Delete error.jl --- lib/BloqadeQMC/test/error.jl | 69 ------------------------------------ 1 file changed, 69 deletions(-) delete mode 100644 lib/BloqadeQMC/test/error.jl diff --git a/lib/BloqadeQMC/test/error.jl b/lib/BloqadeQMC/test/error.jl deleted file mode 100644 index 834e686cb..000000000 --- a/lib/BloqadeQMC/test/error.jl +++ /dev/null @@ -1,69 +0,0 @@ -using Test, BloqadeQMC -using Statistics: mean, std -using Measurements: value, uncertainty -using BloqadeQMC: jackknife, bootstrap -using Random -using RandomNumbers - -rng = Xorshifts.Xoroshiro128Plus(1234) - -@testset "simple numerical tests for jackknife" begin - x = [3,2,11,7,6] - - # Resampling should not change estimate when f = mean. - m = mean(x) - @test value(jackknife(mean, x)) ≈ m - - # Resampling should give a better estimate for e.g. f = inv. We check this using the Julia library Jackknife - @test value(jackknife(inv, x)) ≈ 0.158110765 - @test uncertainty(jackknife(inv, x)) ≈ 0.052468798 -end - -@testset "functional tests for jackknife" begin - - mean_val = 5.0 - x = fill(mean_val, 100) + randn!(rng, zeros(100)) - j_error = uncertainty(jackknife(x)) - #check jackknifed std < original sample std - @test j_error < std(x) - - add(x,y,z) = x+y+z - num_measurements_list = [10,100,1000,10000,100000] - #create a set of normally random values about a mean value for x,y,z for different measuremensts - xs_list = [[fill(mean_val, num_measurements_list[n]) + randn!(rng, zeros(num_measurements_list[n])) for v in 1:3] for n in 1:length(num_measurements_list)] - jackknife_errors = [uncertainty(jackknife(add, xs_list[n][1], xs_list[n][2], xs_list[n][3])) for n in 1:length(num_measurements_list)] - #test that error decreases as the number of observations gets larger - @test issorted(jackknife_errors, lt=Base.isgreater) - -end - -@testset "test for bootstrap error throwing" begin - cube_sum(x,y) = x^3+y^3 - x = [1,2,3] - y = [4,5] - @test_throws ErrorException bootstrap(cube_sum, x, y) - - x = [1,2,3] - y = [4,5,6] - z = [7,8,9] - @test_throws MethodError bootstrap(cube_sum, x, y, z) - - cube_each(x,y) = [x^3, y^3] - @test_throws ErrorException bootstrap(cube_each, x, y) -end - -@testset "functional bootstrap test" begin - mean_val = 5.0 - x = fill(mean_val, 100) + randn!(rng, zeros(100)) - bs_error = uncertainty(bootstrap(x)) - #check bootstraped std < original sample std - @test bs_error < std(x) - - add(x,y,z) = x+y+z - num_measurements_list = [10,100,1000,10000,100000] - #create a set of normally random values about a mean value for x,y,z for different measuremets - xs_list = [[fill(mean_val, num_measurements_list[n]) + randn!(rng, zeros(num_measurements_list[n])) for v in 1:3] for n in 1:length(num_measurements_list)] - bootstrap_errors = [uncertainty(bootstrap(add, xs_list[n][1], xs_list[n][2], xs_list[n][3])) for n in 1:length(num_measurements_list)] - #test that error decreases as the number of observations gets larger - @test issorted(bootstrap_errors, lt=Base.isgreater) -end \ No newline at end of file From f376328c0671d6ab0273a0bd6fe4f6ab9598f5c3 Mon Sep 17 00:00:00 2001 From: Anna Knoerr Date: Thu, 29 Sep 2022 16:18:09 +0200 Subject: [PATCH 13/21] fixed bug in energy binning --- lib/BloqadeQMC/test/ltfim.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/BloqadeQMC/test/ltfim.jl b/lib/BloqadeQMC/test/ltfim.jl index 35ffae412..6c4f60eef 100644 --- a/lib/BloqadeQMC/test/ltfim.jl +++ b/lib/BloqadeQMC/test/ltfim.jl @@ -3,7 +3,7 @@ using Statistics using Random using RandomNumbers using Measurements -using Plots: histogram, display, savefig +# using Plots: histogram, display, savefig using BinningAnalysis using BloqadeQMC: jackknife @@ -144,7 +144,7 @@ THRESHOLD = 2.576 # 99% Two-sided CI of the t-distribution with infinite dofs # Binning analysis for observables based on ns energy(x) = -x / beta + H.energy_shift - BE = LogBinner(energy.(ns) /= nspins(H)) + BE = LogBinner(energy.(ns) / nspins(H)) println("Correlation time for energy_binned samples") τ_energy = tau(BE) println(τ_energy) From 0d9148d6f2f935829d2f0b904331f6c2d8f33556 Mon Sep 17 00:00:00 2001 From: Anna Knoerr Date: Thu, 29 Sep 2022 16:18:56 +0200 Subject: [PATCH 14/21] included ltfim.jl and error.jl in runtests.jl --- lib/BloqadeQMC/test/runtests.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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") From e9c73cc479876498d1ad831390bf4fec696aeb45 Mon Sep 17 00:00:00 2001 From: Anna Knoerr Date: Thu, 29 Sep 2022 17:14:11 +0200 Subject: [PATCH 15/21] readded error.jl --- lib/BloqadeQMC/test/error.jl | 69 ++++++++++++++++++++++++++++++++++++ 1 file changed, 69 insertions(+) create mode 100644 lib/BloqadeQMC/test/error.jl diff --git a/lib/BloqadeQMC/test/error.jl b/lib/BloqadeQMC/test/error.jl new file mode 100644 index 000000000..834e686cb --- /dev/null +++ b/lib/BloqadeQMC/test/error.jl @@ -0,0 +1,69 @@ +using Test, BloqadeQMC +using Statistics: mean, std +using Measurements: value, uncertainty +using BloqadeQMC: jackknife, bootstrap +using Random +using RandomNumbers + +rng = Xorshifts.Xoroshiro128Plus(1234) + +@testset "simple numerical tests for jackknife" begin + x = [3,2,11,7,6] + + # Resampling should not change estimate when f = mean. + m = mean(x) + @test value(jackknife(mean, x)) ≈ m + + # Resampling should give a better estimate for e.g. f = inv. We check this using the Julia library Jackknife + @test value(jackknife(inv, x)) ≈ 0.158110765 + @test uncertainty(jackknife(inv, x)) ≈ 0.052468798 +end + +@testset "functional tests for jackknife" begin + + mean_val = 5.0 + x = fill(mean_val, 100) + randn!(rng, zeros(100)) + j_error = uncertainty(jackknife(x)) + #check jackknifed std < original sample std + @test j_error < std(x) + + add(x,y,z) = x+y+z + num_measurements_list = [10,100,1000,10000,100000] + #create a set of normally random values about a mean value for x,y,z for different measuremensts + xs_list = [[fill(mean_val, num_measurements_list[n]) + randn!(rng, zeros(num_measurements_list[n])) for v in 1:3] for n in 1:length(num_measurements_list)] + jackknife_errors = [uncertainty(jackknife(add, xs_list[n][1], xs_list[n][2], xs_list[n][3])) for n in 1:length(num_measurements_list)] + #test that error decreases as the number of observations gets larger + @test issorted(jackknife_errors, lt=Base.isgreater) + +end + +@testset "test for bootstrap error throwing" begin + cube_sum(x,y) = x^3+y^3 + x = [1,2,3] + y = [4,5] + @test_throws ErrorException bootstrap(cube_sum, x, y) + + x = [1,2,3] + y = [4,5,6] + z = [7,8,9] + @test_throws MethodError bootstrap(cube_sum, x, y, z) + + cube_each(x,y) = [x^3, y^3] + @test_throws ErrorException bootstrap(cube_each, x, y) +end + +@testset "functional bootstrap test" begin + mean_val = 5.0 + x = fill(mean_val, 100) + randn!(rng, zeros(100)) + bs_error = uncertainty(bootstrap(x)) + #check bootstraped std < original sample std + @test bs_error < std(x) + + add(x,y,z) = x+y+z + num_measurements_list = [10,100,1000,10000,100000] + #create a set of normally random values about a mean value for x,y,z for different measuremets + xs_list = [[fill(mean_val, num_measurements_list[n]) + randn!(rng, zeros(num_measurements_list[n])) for v in 1:3] for n in 1:length(num_measurements_list)] + bootstrap_errors = [uncertainty(bootstrap(add, xs_list[n][1], xs_list[n][2], xs_list[n][3])) for n in 1:length(num_measurements_list)] + #test that error decreases as the number of observations gets larger + @test issorted(bootstrap_errors, lt=Base.isgreater) +end \ No newline at end of file From f5a3aebe7006f20c29a4b3b09711ba09af949125 Mon Sep 17 00:00:00 2001 From: Anna Knoerr Date: Thu, 29 Sep 2022 17:16:37 +0200 Subject: [PATCH 16/21] readded TFIM.jl --- lib/BloqadeQMC/src/ising/TFIM.jl | 172 +++++++++++++++++++++++++++++++ 1 file changed, 172 insertions(+) create mode 100644 lib/BloqadeQMC/src/ising/TFIM.jl diff --git a/lib/BloqadeQMC/src/ising/TFIM.jl b/lib/BloqadeQMC/src/ising/TFIM.jl new file mode 100644 index 000000000..1293d23f8 --- /dev/null +++ b/lib/BloqadeQMC/src/ising/TFIM.jl @@ -0,0 +1,172 @@ +abstract type AbstractIsing{O <: AbstractOperatorSampler} <: Hamiltonian{2,O} end +abstract type AbstractTFIM{O <: AbstractOperatorSampler} <: AbstractIsing{O} end + +# H = -J ∑_{⟨ij⟩} σ_i σ_j - h ∑_i σ^x_i +struct NearestNeighbourTFIM{O} <: AbstractTFIM{O} + op_sampler::O + J::Float64 + hx::Float64 + Ns::Int + energy_shift::Float64 +end + + +# H = ∑_{ij} J_ij σ_i σ_j - ∑_i h_i σ^x_i +struct TFIM{O,M <: UpperTriangular{Float64},V <: AbstractVector{Float64}} <: AbstractTFIM{O} + op_sampler::O + J::M + hx::V + Ns::Int + energy_shift::Float64 +end + +const ISING_OP_SIZE = 5 +############################################################################### + +# 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) +@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)) + +############################################################################### + +function make_prob_vector(J::UpperTriangular{T}, hx::AbstractVector{T}) where T + @assert length(hx) == size(J, 1) == size(J, 2) + + ops = Vector{NTuple{ISING_OP_SIZE, Int}}(undef, 0) + p = Vector{T}(undef, 0) + energy_shift = zero(T) + + for i in eachindex(hx) + if !iszero(hx[i]) + push!(ops, makediagonalsiteop(AbstractTFIM, i)) + push!(p, hx[i]) + energy_shift += hx[i] + end + end + + # only take J_ij terms from upper triangle + for j in axes(J, 2), i in axes(J, 1) + if i < j # i != j: we don't want self-interactions + if J[i, j] != 0 + push!(p, 2*abs(J[i, j])) + push!(ops, (2, 1, length(p), i, j)) + energy_shift += abs(J[i, j]) + end + end + end + + return ops, p, energy_shift +end + +function make_uniform_tfim(bond_spins::Vector{NTuple{2,Int}}, Ns::Int, J::T, hx::T) where T + hx_ = hx * ones(T, Ns) + J_ = zeros(T, Ns, Ns) + for (i, j) in bond_spins + i, j = (i <= j) ? (i, j) : (j, i) + J_[i, j] = -J + end + + return UpperTriangular(triu!(J_, 1)), hx_ +end + +############################################################################### + +function TFIM(J::UpperTriangular{Float64}, hx::AbstractVector{Float64}) + @assert length(hx) == size(J, 1) == size(J, 2) + + ops, p, energy_shift = make_prob_vector(J, hx) + Ns = length(hx) + op_sampler = OperatorSampler(ops, p) + + return TFIM{typeof(op_sampler), typeof(J), typeof(hx)}( + op_sampler, J, hx, Ns, energy_shift + ) +end + + +function TFIM(bond_spin, Ns::Int, Nb::Int, hx::Float64, J::Float64) + J_, hx_ = make_uniform_tfim(bond_spin, Ns, J, hx) + return TFIM(J_, hx_) +end + +abstract type HXField; end +struct ConstantHX <: HXField; end +struct VaryingHX <: HXField; end + +hxfield(::AbstractIsing) = VaryingHX() +hxfield(::NearestNeighbourTFIM) = ConstantHX() + +total_hx(::ConstantHX, H::AbstractIsing) = nspins(H) * H.hx +total_hx(::VaryingHX, H::AbstractIsing) = sum(H.hx) +total_hx(H::AbstractIsing) = total_hx(hxfield(H), H) + +# returns true if H.J[site1, site2] is negative +Base.@propagate_inbounds isferromagnetic(H::TFIM, (site1, site2)::NTuple{2, Int}) = signbit(H.J[site1, site2]) +haslongitudinalfield(::AbstractTFIM) = false + + +############################################################################### + + +function energy(::Type{<:BinaryGroundState}, H::AbstractIsing, ns::Vector{<:Real}; resampler::Function=jackknife) + hx = total_hx(H) + + if !iszero(hx) + E = -hx * resampler(inv, ns) + else + E = measurement(zero(H.energy_shift)) + end + + return H.energy_shift + E +end + +function energy(::Type{<:BinaryGroundState}, H::AbstractIsing, ns_mean::Real) + hx = total_hx(H) + + if !iszero(hx) + E = -hx / ns_mean + else + E = zero(H.energy_shift) + end + + return H.energy_shift + E +end From 226020b1f04aa856118cb1205ceb5beb1b299a42 Mon Sep 17 00:00:00 2001 From: Anna Knoerr Date: Fri, 30 Sep 2022 00:05:24 +0200 Subject: [PATCH 17/21] modified N_eff calculation --- lib/BloqadeQMC/test/ltfim.jl | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/lib/BloqadeQMC/test/ltfim.jl b/lib/BloqadeQMC/test/ltfim.jl index 6c4f60eef..46b99f540 100644 --- a/lib/BloqadeQMC/test/ltfim.jl +++ b/lib/BloqadeQMC/test/ltfim.jl @@ -127,7 +127,7 @@ THRESHOLD = 2.576 # 99% Two-sided CI of the t-distribution with infinite dofs τ_abs = tau(B) println(τ_abs) println("N_eff") - N_eff = MCS / τ_abs + 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() @@ -137,7 +137,7 @@ THRESHOLD = 2.576 # 99% Two-sided CI of the t-distribution with infinite dofs τ_abs2 = tau(B) println(τ_abs2) println("N_eff2") - N_eff2 = MCS / τ_abs2 + 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() @@ -149,13 +149,20 @@ THRESHOLD = 2.576 # 99% Two-sided CI of the t-distribution with infinite dofs τ_energy = tau(BE) println(τ_energy) println("N_eff_E") - N_eff_E = MCS / τ_energy + 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") @@ -176,7 +183,7 @@ THRESHOLD = 2.576 # 99% Two-sided CI of the t-distribution with infinite dofs println(energy) println("energy_binned") println(energy_binned) - + println() heat_capacity = jackknife(ns .^ 2, ns) do nsqr, n From dd67c24f9b8ac32f52a69c7e2d2679114528bcc3 Mon Sep 17 00:00:00 2001 From: Anna Knoerr Date: Sat, 1 Oct 2022 03:46:45 +0200 Subject: [PATCH 18/21] restructured updated rules within new updates folder --- lib/BloqadeQMC/src/updates/cluster.jl | 333 ++++++++++++++++++++++ lib/BloqadeQMC/src/updates/diagonal.jl | 135 +++++++++ lib/BloqadeQMC/src/updates/line.jl | 47 +++ lib/BloqadeQMC/src/updates/mc_step.jl | 26 ++ lib/BloqadeQMC/src/updates/multibranch.jl | 47 +++ 5 files changed, 588 insertions(+) create mode 100644 lib/BloqadeQMC/src/updates/cluster.jl create mode 100644 lib/BloqadeQMC/src/updates/diagonal.jl create mode 100644 lib/BloqadeQMC/src/updates/line.jl create mode 100644 lib/BloqadeQMC/src/updates/mc_step.jl create mode 100644 lib/BloqadeQMC/src/updates/multibranch.jl diff --git a/lib/BloqadeQMC/src/updates/cluster.jl b/lib/BloqadeQMC/src/updates/cluster.jl new file mode 100644 index 000000000..628f2ef17 --- /dev/null +++ b/lib/BloqadeQMC/src/updates/cluster.jl @@ -0,0 +1,333 @@ +# 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!() + +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 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..2e2caee56 --- /dev/null +++ b/lib/BloqadeQMC/src/updates/mc_step.jl @@ -0,0 +1,26 @@ +# Todo: Merge the following two mc_step functions + +# mc_step from groundstate.jl + +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 From 82daf9d52b96c5bc7781c89717a093ff5369404d Mon Sep 17 00:00:00 2001 From: Anna Knoerr Date: Mon, 3 Oct 2022 19:24:42 +0200 Subject: [PATCH 19/21] modified include structure to test new update files --- lib/BloqadeQMC/src/BloqadeQMC.jl | 1 + lib/BloqadeQMC/src/updates/updates.jl | 5 +++++ 2 files changed, 6 insertions(+) create mode 100644 lib/BloqadeQMC/src/updates/updates.jl 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/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 From 6fb0f5f80ffa42fada02fcba050e5167cad4b7da Mon Sep 17 00:00:00 2001 From: Anna Knoerr Date: Mon, 3 Oct 2022 19:28:25 +0200 Subject: [PATCH 20/21] copied over missing cluster_update!() function --- lib/BloqadeQMC/src/updates/cluster.jl | 12 +++++++++++- lib/BloqadeQMC/src/updates/mc_step.jl | 10 ++++++++++ 2 files changed, 21 insertions(+), 1 deletion(-) diff --git a/lib/BloqadeQMC/src/updates/cluster.jl b/lib/BloqadeQMC/src/updates/cluster.jl index 628f2ef17..8f0dedfcc 100644 --- a/lib/BloqadeQMC/src/updates/cluster.jl +++ b/lib/BloqadeQMC/src/updates/cluster.jl @@ -235,7 +235,7 @@ trialstate_weight_change(qmc_state::BinaryThermalState, lsize::Int, Ns::Int, i:: ############################################################################# -# changed the function name below from cluster_update!() to generic_cluster_update!() +# 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) @@ -331,3 +331,13 @@ function generic_cluster_update!(rng::AbstractRNG, update_kernel!::Function, acc 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/mc_step.jl b/lib/BloqadeQMC/src/updates/mc_step.jl index 2e2caee56..40f3a3bcd 100644 --- a/lib/BloqadeQMC/src/updates/mc_step.jl +++ b/lib/BloqadeQMC/src/updates/mc_step.jl @@ -2,6 +2,16 @@ # 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...) From d664212e1da72cd73c9b751b0fd4c7990c9c91f1 Mon Sep 17 00:00:00 2001 From: Anna Knoerr Date: Tue, 4 Oct 2022 20:54:49 +0200 Subject: [PATCH 21/21] file suggestions for grouping evth related to QMC datastructures --- lib/BloqadeQMC/src/datastructures/Rydberg.jl | 155 +++++++++++++++++ lib/BloqadeQMC/src/datastructures/alias.jl | 90 ++++++++++ .../src/datastructures/hamiltonian.jl | 52 ++++++ .../src/datastructures/improved_op_sampler.jl | 75 +++++++++ .../src/datastructures/operators.jl | 73 ++++++++ .../src/datastructures/probabilityvector.jl | 20 +++ .../src/datastructures/qmc_state.jl | 157 ++++++++++++++++++ 7 files changed, 622 insertions(+) create mode 100644 lib/BloqadeQMC/src/datastructures/Rydberg.jl create mode 100644 lib/BloqadeQMC/src/datastructures/alias.jl create mode 100644 lib/BloqadeQMC/src/datastructures/hamiltonian.jl create mode 100644 lib/BloqadeQMC/src/datastructures/improved_op_sampler.jl create mode 100644 lib/BloqadeQMC/src/datastructures/operators.jl create mode 100644 lib/BloqadeQMC/src/datastructures/probabilityvector.jl create mode 100644 lib/BloqadeQMC/src/datastructures/qmc_state.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