Add a Sobol sequence wrapper that includes the initial point (0, 0, ...)

C.f. JuliaMath/Sobol.jl#31

Author: krivenko

src/inchworm.jl

@@ -2,12 +2,12 @@ module inchworm
 
 using Printf
 
-import Sobol: SobolSeq
-
 import Keldysh; kd = Keldysh
 
 import QInchworm; teval = QInchworm.topology_eval
 
+using  QInchworm.utility: BetterSobolSeq, next!
+
 import QInchworm.configuration: Expansion,
                                 Configurations,
                                 operator,
@@ -81,7 +81,7 @@ function inchworm_step(expansion::Expansion,
 
             teval.update_inch_times!(od.configurations, t_i, t_w, t_f)
 
-            seq = SobolSeq(2 * od.order)
+            seq = BetterSobolSeq(2 * od.order)
             N = 0
             order_contrib = deepcopy(zero_sector_block_matrix)
             order_contrib_prev = deepcopy(zero_sector_block_matrix)
@@ -156,7 +156,7 @@ function inchworm_step_bare(expansion::Expansion,
             teval.update_inch_times!(od.configurations, t_i, t_i, t_f)
 
             d = 2 * od.order
-            seq = SobolSeq(d)
+            seq = BetterSobolSeq(d)
             N = 0
             order_contrib = deepcopy(zero_sector_block_matrix)
             order_contrib_prev = deepcopy(zero_sector_block_matrix)

src/qmc_integrate.jl

@@ -1,9 +1,8 @@
 module qmc_integrate
 
-import Sobol: SobolSeq, next!
 import Keldysh; kd = Keldysh
 
-import QInchworm.utility: get_ref
+using QInchworm.utility: get_ref, BetterSobolSeq, next!
 
 #
 # I use the notations introduced in https://arxiv.org/pdf/2002.12372.pdf.
@@ -170,7 +169,7 @@ function qmc_time_ordered_integral(f,
                                    t_i::kd.BranchPoint,
                                    t_f::kd.BranchPoint;
                                    init = zero(contour_function_return_type(f)),
-                                   seq = SobolSeq(d),
+                                   seq = BetterSobolSeq(d),
                                    τ::Real,
                                    N::Int)
     @assert kd.heaviside(c, t_f, t_i)
@@ -225,7 +224,7 @@ function qmc_time_ordered_integral_n_samples(
     t_i::kd.BranchPoint,
     t_f::kd.BranchPoint;
     init = zero(contour_function_return_type(f)),
-    seq = SobolSeq(d),
+    seq = BetterSobolSeq(d),
     τ::Real,
     N_samples::Int)
     @assert kd.heaviside(c, t_f, t_i)
@@ -276,7 +275,7 @@ function qmc_time_ordered_integral_sort(f,
                                         t_i::kd.BranchPoint,
                                         t_f::kd.BranchPoint;
                                         init = zero(contour_function_return_type(f)),
-                                        seq = SobolSeq(d),
+                                        seq = BetterSobolSeq(d),
                                         N::Int)
     @assert kd.heaviside(c, t_f, t_i)
 
@@ -321,7 +320,7 @@ function qmc_time_ordered_integral_root(f,
                                         t_i::kd.BranchPoint,
                                         t_f::kd.BranchPoint;
                                         init = zero(contour_function_return_type(f)),
-                                        seq = SobolSeq(d),
+                                        seq = BetterSobolSeq(d),
                                         N::Int)
     @assert kd.heaviside(c, t_f, t_i)
 
@@ -394,7 +393,7 @@ function qmc_inchworm_integral_root(f,
                                     t_w::kd.BranchPoint,
                                     t_f::kd.BranchPoint;
                                     init = zero(contour_function_return_type(f)),
-                                    seq = SobolSeq(d_bold + d_bare),
+                                    seq = BetterSobolSeq(d_bold + d_bare),
                                     N::Int)
     @assert kd.heaviside(c, t_w, t_i)
     @assert kd.heaviside(c, t_f, t_w)

src/utility.jl

@@ -9,6 +9,8 @@ using Interpolations: BoundaryCondition,
                       WoodburyMatrices,
                       lut!
 
+import Sobol: AbstractSobolSeq, SobolSeq, next!, ndims
+
 """
     Inverse of get_point(c::AbstractContour, ref)
 
@@ -92,4 +94,25 @@ function (spline::IncrementalSpline{T})(z) where {T<:Number}
     @inbounds spline.data[i] * (1-δx^2) + spline.data[i + 1] * (1-(1-δx)^2) + c3 * (0.25-(δx-0.5)^2)
 end
 
+"""
+    Sobol sequence including the initial point (0, 0, ...)
+
+    C.f. https://github.com/stevengj/Sobol.jl/issues/31
+"""
+mutable struct BetterSobolSeq{N} <: AbstractSobolSeq{N}
+    seq::SobolSeq{N}
+    init_pt_returned::Bool
+
+    BetterSobolSeq(N::Int) = new{N}(SobolSeq(N), false)
+end
+
+function next!(bseq::BetterSobolSeq)
+    if bseq.init_pt_returned
+        next!(bseq.seq)
+    else
+        bseq.init_pt_returned = true
+        zeros(Float64, ndims(bseq.seq))
+    end
+end
+
 end

test/Project.toml

@@ -1,4 +1,5 @@
 [deps]
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 Keldysh = "50f2bc7e-5fd7-11e9-13c7-85fb88b4f34f"
 KeldyshED = "675b6b9c-7c2f-11e9-3bf3-dfd4d61640f7"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

test/nca_equil.jl

@@ -46,7 +46,7 @@ end
     tmax = 1.
 
     # -- Quasi Monte Carlo
-    N = 10000
+    N = 2^13
 
     # -- Exact Diagonalization
 

test/qmc_integrate.jl

@@ -21,7 +21,7 @@ import QInchworm.qmc_integrate: qmc_time_ordered_integral,
         τ = 5.0
 
         @testset "d = 1, constant integrand" begin
-            let d = 1, c = contour, f = t -> 1.0, N = 200000
+            let d = 1, c = contour, f = t -> 1.0, N = 2^17
                 val, N_samples = qmc_time_ordered_integral(f, d, c, c(0.1), c(0.5); τ=τ, N=N)
                 @show N_samples / N
                 @test isapprox(val, -0.4, rtol=1e-4)
@@ -41,7 +41,7 @@ import QInchworm.qmc_integrate: qmc_time_ordered_integral,
         end
 
         @testset "d = 1, linear in t integrand" begin
-            let d = 1, c = contour, f = t -> t[1].val, N = 200000
+            let d = 1, c = contour, f = t -> t[1].val, N = 2^17
                 val, N_samples = qmc_time_ordered_integral(f, d, c, c(0.1), c(0.5), τ=τ, N=N)
                 @show N_samples / N
                 @test isapprox(val, -0.28, rtol=1e-4)
@@ -61,7 +61,7 @@ import QInchworm.qmc_integrate: qmc_time_ordered_integral,
         end
 
         @testset "d = 2, constant integrand" begin
-            let d = 2, c = contour, f = t -> (@assert kd.heaviside(c, t[1], t[2]); 1.0), N = 500000
+            let d = 2, c = contour, f = t -> (@assert kd.heaviside(c, t[1], t[2]); 1.0), N = 2^18
                 val, N_samples = qmc_time_ordered_integral(f, d, c, c(0.1), c(0.5), τ=τ, N=N)
                 @show N_samples / N
                 @test isapprox(val, 0.08, rtol=5e-3)
@@ -81,7 +81,7 @@ import QInchworm.qmc_integrate: qmc_time_ordered_integral,
         end
 
         @testset "d = 2, bilinear in t1, t2 integrand" begin
-            let d = 2, c = contour, f = t -> (@assert kd.heaviside(c, t[1], t[2]); t[1].val * t[2].val), N = 1000000
+            let d = 2, c = contour, f = t -> (@assert kd.heaviside(c, t[1], t[2]); t[1].val * t[2].val), N = 2^19
                 val, N_samples = qmc_time_ordered_integral(f, d, c, c(0.1), c(0.5), τ=τ, N=N)
                 @show N_samples / N
                 @test isapprox(val, 0.0392, rtol=5e-3)
@@ -106,7 +106,7 @@ import QInchworm.qmc_integrate: qmc_time_ordered_integral,
         τ = 5.0
 
         @testset "d = 1, constant integrand" begin
-            let d = 1, c = contour, f = t -> 1.0, N_samples = 50000
+            let d = 1, c = contour, f = t -> 1.0, N_samples = 2^16
                 val, N = qmc_time_ordered_integral_n_samples(f, d, c, c(0.1), c(0.5); τ=τ, N_samples=N_samples)
                 @show N_samples / N
                 @test isapprox(val, -0.4, rtol=1e-4)
@@ -126,7 +126,7 @@ import QInchworm.qmc_integrate: qmc_time_ordered_integral,
         end
 
         @testset "d = 1, linear in t integrand" begin
-            let d = 1, c = contour, f = t -> t[1].val, N_samples = 50000
+            let d = 1, c = contour, f = t -> t[1].val, N_samples = 2^15
                 val, N = qmc_time_ordered_integral_n_samples(f, d, c, c(0.1), c(0.5), τ=τ, N_samples=N_samples)
                 @show N_samples / N
                 @test isapprox(val, -0.28, rtol=1e-4)
@@ -146,7 +146,7 @@ import QInchworm.qmc_integrate: qmc_time_ordered_integral,
         end
 
         @testset "d = 2, constant integrand" begin
-            let d = 2, c = contour, f = t -> (@assert kd.heaviside(c, t[1], t[2]); 1.0), N_samples = 50000
+            let d = 2, c = contour, f = t -> (@assert kd.heaviside(c, t[1], t[2]); 1.0), N_samples = 2^15
                 val, N = qmc_time_ordered_integral_n_samples(f, d, c, c(0.1), c(0.5), τ=τ, N_samples=N_samples)
                 @show N_samples / N
                 @test isapprox(val, 0.08, rtol=5e-3)
@@ -166,7 +166,7 @@ import QInchworm.qmc_integrate: qmc_time_ordered_integral,
         end
 
         @testset "d = 2, bilinear in t1, t2 integrand" begin
-            let d = 2, c = contour, f = t -> (@assert kd.heaviside(c, t[1], t[2]); t[1].val * t[2].val), N_samples = 50000
+            let d = 2, c = contour, f = t -> (@assert kd.heaviside(c, t[1], t[2]); t[1].val * t[2].val), N_samples = 2^15
                 val, N = qmc_time_ordered_integral_n_samples(f, d, c, c(0.1), c(0.5), τ=τ, N_samples=N_samples)
                 @show N_samples / N
                 @test isapprox(val, 0.0392, rtol=5e-3)
@@ -189,7 +189,7 @@ import QInchworm.qmc_integrate: qmc_time_ordered_integral,
     @testset "qmc_time_ordered_integral_sort()" begin
 
         @testset "d = 1, constant integrand" begin
-            let d = 1, c = contour, f = t -> 1.0, N = 200000
+            let d = 1, c = contour, f = t -> 1.0, N = 2^17
                 val = qmc_time_ordered_integral_sort(f, d, c, c(0.1), c(0.5); N=N)
                 @test isapprox(val, -0.4, rtol=1e-4)
                 val = qmc_time_ordered_integral_sort(f, d, c, c(0.1), c(2.0); N=N)
@@ -204,7 +204,7 @@ import QInchworm.qmc_integrate: qmc_time_ordered_integral,
         end
 
         @testset "d = 1, linear in t integrand" begin
-            let d = 1, c = contour, f = t -> t[1].val, N = 200000
+            let d = 1, c = contour, f = t -> t[1].val, N = 2^18
                 val = qmc_time_ordered_integral_sort(f, d, c, c(0.1), c(0.5), N=N)
                 @test isapprox(val, -0.28, rtol=1e-5)
                 val = qmc_time_ordered_integral_sort(f, d, c, c(0.1), c(2.0), N=N)
@@ -219,7 +219,7 @@ import QInchworm.qmc_integrate: qmc_time_ordered_integral,
         end
 
         @testset "d = 2, constant integrand" begin
-            let d = 2, c = contour, f = t -> (@assert kd.heaviside(c, t[1], t[2]); 1.0), N = 500000
+            let d = 2, c = contour, f = t -> (@assert kd.heaviside(c, t[1], t[2]); 1.0), N = 2^18
                 val = qmc_time_ordered_integral_sort(f, d, c, c(0.1), c(0.5), N=N)
                 @test isapprox(val, 0.08, rtol=1e-4)
                 val = qmc_time_ordered_integral_sort(f, d, c, c(0.1), c(2.0), N=N)
@@ -234,7 +234,7 @@ import QInchworm.qmc_integrate: qmc_time_ordered_integral,
         end
 
         @testset "d = 2, bilinear in t1, t2 integrand" begin
-            let d = 2, c = contour, f = t -> (@assert kd.heaviside(c, t[1], t[2]); t[1].val * t[2].val), N = 1000000
+            let d = 2, c = contour, f = t -> (@assert kd.heaviside(c, t[1], t[2]); t[1].val * t[2].val), N = 2^19
                 val = qmc_time_ordered_integral_sort(f, d, c, c(0.1), c(0.5), N=N)
                 @test isapprox(val, 0.0392, rtol=1e-4)
                 val = qmc_time_ordered_integral_sort(f, d, c, c(0.1), c(2.0), N=N)
@@ -252,7 +252,7 @@ import QInchworm.qmc_integrate: qmc_time_ordered_integral,
     @testset "qmc_time_ordered_integral_root()" begin
 
         @testset "d = 1, constant integrand" begin
-            let d = 1, c = contour, f = t -> 1.0, N = 200000
+            let d = 1, c = contour, f = t -> 1.0, N = 2^17
                 val = qmc_time_ordered_integral_root(f, d, c, c(0.1), c(0.5); N=N)
                 @test isapprox(val, -0.4, rtol=1e-4)
                 val = qmc_time_ordered_integral_root(f, d, c, c(0.1), c(2.0); N=N)
@@ -267,7 +267,7 @@ import QInchworm.qmc_integrate: qmc_time_ordered_integral,
         end
 
         @testset "d = 1, linear in t integrand" begin
-            let d = 1, c = contour, f = t -> t[1].val, N = 200000
+            let d = 1, c = contour, f = t -> t[1].val, N = 2^18
                 val = qmc_time_ordered_integral_root(f, d, c, c(0.1), c(0.5), N=N)
                 @test isapprox(val, -0.28, rtol=1e-5)
                 val = qmc_time_ordered_integral_root(f, d, c, c(0.1), c(2.0), N=N)
@@ -282,7 +282,7 @@ import QInchworm.qmc_integrate: qmc_time_ordered_integral,
         end
 
         @testset "d = 2, constant integrand" begin
-            let d = 2, c = contour, f = t -> (@assert kd.heaviside(c, t[1], t[2]); 1.0), N = 500000
+            let d = 2, c = contour, f = t -> (@assert kd.heaviside(c, t[1], t[2]); 1.0), N = 2^18
                 val = qmc_time_ordered_integral_root(f, d, c, c(0.1), c(0.5), N=N)
                 @test isapprox(val, 0.08, rtol=1e-4)
                 val = qmc_time_ordered_integral_root(f, d, c, c(0.1), c(2.0), N=N)
@@ -297,7 +297,7 @@ import QInchworm.qmc_integrate: qmc_time_ordered_integral,
         end
 
         @testset "d = 2, bilinear in t1, t2 integrand" begin
-            let d = 2, c = contour, f = t -> (@assert kd.heaviside(c, t[1], t[2]); t[1].val * t[2].val), N = 1000000
+            let d = 2, c = contour, f = t -> (@assert kd.heaviside(c, t[1], t[2]); t[1].val * t[2].val), N = 2^19
                 val = qmc_time_ordered_integral_root(f, d, c, c(0.1), c(0.5), N=N)
                 @test isapprox(val, 0.0392, rtol=1e-4)
                 val = qmc_time_ordered_integral_root(f, d, c, c(0.1), c(2.0), N=N)
@@ -315,7 +315,7 @@ import QInchworm.qmc_integrate: qmc_time_ordered_integral,
     @testset "qmc_inchworm_integral_root()" begin
 
         @testset "d_bold = 0, d_bare = 3, constant integrand" begin
-            let d_bold = 0, d_bare = 3, c = contour, t_i = c(1.1), t_w = c(5.0), t_f = c(5.5), N = 1000000
+            let d_bold = 0, d_bare = 3, c = contour, t_i = c(1.1), t_w = c(5.0), t_f = c(5.5), N = 2^19
                 function f(t)
                     @assert kd.heaviside(c, t_f, t[1])
                     @assert kd.heaviside(c, t[1], t[2])
@@ -329,7 +329,7 @@ import QInchworm.qmc_integrate: qmc_time_ordered_integral,
         end
 
         @testset "d_bold = 0, d_bare = 3, multilinear integrand" begin
-            let d_bold = 0, d_bare = 3, c = contour, t_i = c(1.1), t_w = c(5.0), t_f = c(5.5), N = 1000000
+            let d_bold = 0, d_bare = 3, c = contour, t_i = c(1.1), t_w = c(5.0), t_f = c(5.5), N = 2^19
                 function f(t)
                     @assert kd.heaviside(c, t_f, t[1])
                     @assert kd.heaviside(c, t[1], t[2])
@@ -343,7 +343,7 @@ import QInchworm.qmc_integrate: qmc_time_ordered_integral,
         end
 
         @testset "d_bold = 3, d_bare = 2, constant integrand" begin
-            let d_bold = 3, d_bare = 2, c = contour, t_i = c(1.1), t_w = c(5.0), t_f = c(5.5), N = 2000000
+            let d_bold = 3, d_bare = 2, c = contour, t_i = c(1.1), t_w = c(5.0), t_f = c(5.5), N = 2^20
                 function f(t)
                     @assert kd.heaviside(c, t_f, t[1])
                     @assert kd.heaviside(c, t[1], t[2])
@@ -360,7 +360,7 @@ import QInchworm.qmc_integrate: qmc_time_ordered_integral,
         end
 
         @testset "d_bold = 3, d_bare = 2, multilinear integrand" begin
-            let d_bold = 3, d_bare = 2, c = contour, t_i = c(1.1), t_w = c(5.0), t_f = c(5.5), N = 2000000
+            let d_bold = 3, d_bare = 2, c = contour, t_i = c(1.1), t_w = c(5.0), t_f = c(5.5), N = 2^20
                 function f(t)
                     @assert kd.heaviside(c, t_f, t[1])
                     @assert kd.heaviside(c, t[1], t[2])