add multi-threading for unbiasedEstimates

demo/lgDemoEstimateComparison.jl

@@ -0,0 +1,85 @@
+import RData
+using RNGPool
+using CoupledConditionalSMC
+
+# data generated with A = 0.95 (not 0.9)
+dataJLS = RData.load("demo/ar1data.RData")
+ys = vec(dataJLS["observations"])
+
+include("lgModel.jl")
+
+# they start at time 0 with observations from time 1,2,...
+# so we have an initial distribution of N(0,1+0.9^2)
+theta = LGTheta(0.9, 1.0, 1.0, 1.0, 0.0, 1.81)
+model = makeLGModel(theta, ys)
+lM = LinearGaussian.makelM(theta)
+ko = kalman(theta, ys)
+
+import Statistics: mean, std
+
+setRNGs(12345)
+
+function runDemo(model, N, n, m, AT, AS, BS)
+  println("LG Model with N = ", N, ", n = ", n)
+
+  function h(path::Vector{Float64Particle})
+    return path[1].x
+  end
+
+  println("\n-----\nTrue value = ", ko.smoothingMeans[1])
+
+  if AT
+    println("Ancestral Tracing:")
+    times, values = CoupledConditionalSMC.unbiasedEstimates(model, h, N, n, 1,
+      m, true, true)
+    println(mean(times), ", ", std(times))
+    println(mean(values), ", ", std(values))
+    println()
+  end
+
+  if AS
+    println("Ancestor Sampling:")
+    times, values = CoupledConditionalSMC.unbiasedEstimates(model, lM, h, N, n, 1,
+      m, :AS, true, true)
+    println(mean(times), ", ", std(times))
+    println(mean(values), ", ", std(values))
+    println()
+  end
+
+  if BS
+    println("Backward Sampling:")
+    times, values  = CoupledConditionalSMC.unbiasedEstimates(model, lM, h, N, n, 1,
+      m, :BS, true, true)
+    println(mean(times), ", ", std(times))
+    println(mean(values), ", ", std(values))
+    println()
+  end
+
+end
+
+# N = 128 T = 50 17.84 (17.13) 7.73 (5.11)
+# N = 256 T = 100 13.16 (11.09) 7.59 (5.05)
+# N = 512 T = 200 12.52 (10.64) 6.77 (3.85)
+# N = 1024 T = 400 12.74 (10.96) 6.77 (3.47)
+# N = 2048 T = 800 13.58 (9.56) 6.34 (2.95)
+
+m = 1000
+
+runDemo(model, 64, 50, m, true, true, true)
+runDemo(model, 128, 50, m, true, true, true)
+
+runDemo(model, 128, 100, m, true, true, true)
+runDemo(model, 256, 100, m, true, true, true)
+
+runDemo(model, 256, 200, m, true, true, true)
+runDemo(model, 512, 200, m, true, true, true)
+
+runDemo(model, 512, 400, m, true, true, true)
+runDemo(model, 1024, 400, m, true, true, true)
+
+runDemo(model, 64, 800, m, false, true, true)
+runDemo(model, 128, 800, m, false, true, true)
+runDemo(model, 256, 800, m, false, true, true)
+runDemo(model, 512, 800, m, false, true, true)
+runDemo(model, 1024, 800, m, true, true, true)
+runDemo(model, 2048, 800, m, true, true, true)

demo/lgDemoEstimates.jl

@@ -7,62 +7,60 @@ include("lgModel.jl")
 function runDemo(n::Int64, N::Int64, b::Int64, m::Int64, maxit::Int64)
   model, lM, ko = setupLGModel(n, 1.0)
 
-  ccsmcio = CCSMCIO{model.particle, model.pScratch}(N, model.maxn)
-
   function h(path::Vector{Float64Particle})
     return path[1].x
   end
 
   println("\n-----\nTrue value = ", ko.smoothingMeans[1])
 
   println("\nBackward sampling, dependent initialization")
-  resultsBS = CoupledConditionalSMC.unbiasedEstimates(model, lM, h, ccsmcio,
-    b, m, :BS, false, maxit)
+  resultsBS = CoupledConditionalSMC.unbiasedEstimates(model, lM, h, N, n,
+    b, m, :BS, false, true, maxit)
   println("mean of estimates = ", mean(resultsBS[2]))
   println("mean no. iterations = ", mean(resultsBS[1]))
   println("estimated variance of estimator = ", var(resultsBS[2]))
   println("should be approx. standard normal: ",
     (mean(resultsBS[2]) - ko.smoothingMeans[1])/sqrt(var(resultsBS[2])/m))
 
   println("\nBackward sampling, independent initialization")
-  resultsBS = CoupledConditionalSMC.unbiasedEstimates(model, lM, h, ccsmcio,
-    b, m, :BS, true, maxit)
+  resultsBS = CoupledConditionalSMC.unbiasedEstimates(model, lM, h, N, n,
+    b, m, :BS, true, true, maxit)
   println("mean of estimates = ", mean(resultsBS[2]))
   println("mean no. iterations = ", mean(resultsBS[1]))
   println("estimated variance of estimator = ", var(resultsBS[2]))
   println("should be approx. standard normal: ",
     (mean(resultsBS[2]) - ko.smoothingMeans[1])/sqrt(var(resultsBS[2])/m))
 
   println("\nAncestral tracing, dependent initialization")
-  results = CoupledConditionalSMC.unbiasedEstimates(model, h, ccsmcio, b, m,
-    false, maxit)
+  results = CoupledConditionalSMC.unbiasedEstimates(model, h, N, n, b, m,
+    false, true, maxit)
   println("mean of estimates = ", mean(results[2]))
   println("mean no. iterations = ", mean(results[1]))
   println("estimated variance of estimator = ", var(results[2]))
   println("should be approx. standard normal: ",
     (mean(results[2]) - ko.smoothingMeans[1])/sqrt(var(results[2])/m))
 
   println("\nAncestral tracing, independent initialization")
-  results = CoupledConditionalSMC.unbiasedEstimates(model, h, ccsmcio, b, m,
-    true, maxit)
+  results = CoupledConditionalSMC.unbiasedEstimates(model, h, N, n, b, m,
+    true, true, maxit)
   println("mean of estimates = ", mean(results[2]))
   println("mean no. iterations = ", mean(results[1]))
   println("estimated variance of estimator = ", var(results[2]))
   println("should be approx. standard normal: ",
     (mean(results[2]) - ko.smoothingMeans[1])/sqrt(var(results[2])/m))
 
   println("\nAncestor sampling, dependent initialization")
-  results = CoupledConditionalSMC.unbiasedEstimates(model, lM, h, ccsmcio, b, m,
-    :AS, false, maxit)
+  results = CoupledConditionalSMC.unbiasedEstimates(model, lM, h, N, n, b, m,
+    :AS, false, true, maxit)
   println("mean of estimates = ", mean(results[2]))
   println("mean no. iterations = ", mean(results[1]))
   println("estimated variance of estimator = ", var(results[2]))
   println("should be approx. standard normal: ",
     (mean(results[2]) - ko.smoothingMeans[1])/sqrt(var(results[2])/m))
 
   println("\nAncestor sampling, independent initialization")
-  results = CoupledConditionalSMC.unbiasedEstimates(model, lM, h, ccsmcio, b, m,
-    :AS, true, maxit)
+  results = CoupledConditionalSMC.unbiasedEstimates(model, lM, h, N, n, b, m,
+    :AS, true, true, maxit)
   println("mean of estimates = ", mean(results[2]))
   println("mean no. iterations = ", mean(results[1]))
   println("estimated variance of estimator = ", var(results[2]))

demo/lgDemoJLSComparison.jl

@@ -63,7 +63,7 @@ end
 # N = 1024 T = 400 12.74 (10.96) 6.77 (3.47)
 # N = 2048 T = 800 13.58 (9.56) 6.34 (2.95)
 
-m = 128
+m = 1000
 
 runDemo(model, 64, 50, m, true, true, true)
 runDemo(model, 128, 50, m, true, true, true)

demo/simpleDemoVisualize.jl

@@ -31,14 +31,16 @@ visualizeCCSMC(simpleModel, 128, true, true, 100)
 visualizeCCSMC(simpleModel, 2^14, true, false, 100)
 visualizeCCSMC(simpleModel, 2^14, true, true, 100)
 
-## CCAPF performs well for N ∈ {2,4,8} but not the rest
+## CCAPF performs well for N ∈ {2,4} but struggles for moderate values
 visualizeCCSMC(simpleModel, simplelM, 2, :AS, true)
 visualizeCCSMC(simpleModel, simplelM, 2, :AS, true, true)
 visualizeCCSMC(simpleModel, simplelM, 4, :AS, true)
 visualizeCCSMC(simpleModel, simplelM, 4, :AS, true, true)
 visualizeCCSMC(simpleModel, simplelM, 8, :AS, true)
 visualizeCCSMC(simpleModel, simplelM, 8, :AS, true, true)
 visualizeCCSMC(simpleModel, simplelM, 16, :AS, true, true, 1000)
-visualizeCCSMC(simpleModel, simplelM, 32, :AS, true, true, 1000)
 visualizeCCSMC(simpleModel, simplelM, 64, :AS, true, true, 1000)
-visualizeCCSMC(simpleModel, simplelM, 128, :AS, true, true, 1000)
+visualizeCCSMC(simpleModel, simplelM, 256, :AS, true, true, 100)
+visualizeCCSMC(simpleModel, simplelM, 1024, :AS, true, true, 100)
+visualizeCCSMC(simpleModel, simplelM, 4096, :AS, true, true, 100)
+visualizeCCSMC(simpleModel, simplelM, 2^14, :AS, true, true, 100)

src/unbiased.jl

@@ -37,32 +37,42 @@ function unbiasedEstimate(model::SMCModel, h::F, ccsmcio::CCSMCIO{Particle},
 end
 
 function unbiasedEstimates(model::SMCModel, lM::F1, h::F2,
-  ccsmcio::CCSMCIO{Particle}, b::Int64, m::Int64, algorithm::Symbol = :BS,
+  N::Int64, n::Int64, b::Int64, m::Int64, algorithm::Symbol = :BS,
   independentInitialization::Bool = false, rngCouple::Bool = false,
   maxit::Int64 = typemax(Int64)) where {F1<:Function, F2<:Function, Particle}
 
+  nt = Threads.nthreads()
+  ccsmcios::Vector{CCSMCIO{model.particle, model.pScratch}} =
+    Vector{CCSMCIO{model.particle, model.pScratch}}(undef, nt)
+  Threads.@threads for i in 1:nt
+    ccsmcios[i] = CCSMCIO{model.particle, model.pScratch}(N, n)
+  end
+
   # just to get the type of v
-  initializeCCSMC(model, ccsmcio, independentInitialization)
-  v = h(ccsmcio.ref1)
+  initializeCCSMC(model, ccsmcios[1], independentInitialization)
+  v = h(ccsmcios[1].ref1)
   T = typeof(v)
 
   iterations = Vector{Int64}(undef, m)
   values = Vector{T}(undef, m)
 
-  @showprogress 10 for i in 1:m
-    v = unbiasedEstimate(model, lM, h, ccsmcio, b, algorithm,
+  p = Progress(div(m, nt), 10)
+  Threads.@threads for i in 1:m
+    tid = Threads.threadid()
+    v = unbiasedEstimate(model, lM, h, ccsmcios[tid], b, algorithm,
       independentInitialization, rngCouple, maxit)
     iterations[i] = v[1]
     values[i] = v[2]
+    tid == 1 && update!(p, i)
   end
 
   return iterations, values
 end
 
-function unbiasedEstimates(model::SMCModel, h::F, ccsmcio::CCSMCIO{Particle},
+function unbiasedEstimates(model::SMCModel, h::F, N::Int64, n::Int64,
   b::Int64, m::Int64, independentInitialization::Bool = false,
   rngCouple::Bool = false, maxit::Int64 = typemax(Int64)) where
   {F<:Function, Particle}
-  return unbiasedEstimates(model, error, h, ccsmcio, b, m, :AT,
+  return unbiasedEstimates(model, error, h, N, n, b, m, :AT,
     independentInitialization, rngCouple, maxit)
 end

test/estimate_test.jl

@@ -31,17 +31,17 @@ function testEstimate(n::Int64, N::Int64, b::Int64, m::Int64, maxit::Int64)
 
   trueValue = ko.smoothingMeans[1]
 
-  resultsBS = CoupledConditionalSMC.unbiasedEstimates(model, lM, h, ccsmcio, b, m,
+  resultsBS = CoupledConditionalSMC.unbiasedEstimates(model, lM, h, N, n, b, m,
     :BS)
   @test mean(resultsBS[2]) ≈ trueValue atol=0.1
   Zval = (mean(resultsBS[2]) - trueValue)/sqrt(var(resultsBS[2])/m)
   @test abs(Zval) < 3
 
-  results = CoupledConditionalSMC.unbiasedEstimates(model, h, ccsmcio, b, m)
+  results = CoupledConditionalSMC.unbiasedEstimates(model, h, N, n, b, m)
   Zval = (mean(results[2]) - trueValue)/sqrt(var(results[2])/m)
   @test abs(Zval) < 3
 
-  results = CoupledConditionalSMC.unbiasedEstimates(model, lM, h, ccsmcio, b, m,
+  results = CoupledConditionalSMC.unbiasedEstimates(model, lM, h, N, n, b, m,
     :AS)
   Zval = (mean(results[2]) - trueValue)/sqrt(var(results[2])/m)
   @test abs(Zval) < 3