From 37980b22018bc109a01ae747962fac6aba85d6cb Mon Sep 17 00:00:00 2001
From: David Widmann <david.widmann@it.uu.se>
Date: Thu, 17 Jun 2021 02:27:05 +0200
Subject: [PATCH] Reduce allocations

---
 src/mcse.jl | 42 ++++++++++++++++++++++++++----------------
 1 file changed, 26 insertions(+), 16 deletions(-)

diff --git a/src/mcse.jl b/src/mcse.jl
index 3db9852..e222d22 100644
--- a/src/mcse.jl
+++ b/src/mcse.jl
@@ -21,35 +21,45 @@ function mcse_bm(x::AbstractVector{<:Real}; size::Int=floor(Int, sqrt(length(x))
     n = length(x)
     m = min(div(n, 2), size)
     m == size || @warn "batch size was reduced to $m"
-    mbar = [Statistics.mean(@view(x[(i + 1):(i + m)])) for i in 0:m:(n - m)]
-    return StatsBase.sem(mbar)
+    mcse = StatsBase.sem(Statistics.mean(@view(x[(i + 1):(i + m)])) for i in 0:m:(n - m))
+    return mcse
 end
 
 function mcse_imse(x::AbstractVector{<:Real})
     n = length(x)
-    m = div(n - 2, 2)
-    x_ = map(Float64, x)
-    ghat = StatsBase.autocov(x_, [0, 1])
+    lags = [0, 1]
+    ghat = StatsBase.autocov(x, lags)
     Ghat = sum(ghat)
-    value = -ghat[1] + 2 * Ghat
-    for i in 1:m
-        Ghat = min(Ghat, sum(StatsBase.autocov(x_, [2 * i, 2 * i + 1])))
+    @inbounds value = Ghat + ghat[2]
+    @inbounds for i in 2:2:(n - 2)
+        lags[1] = i
+        lags[2] = i + 1
+        StatsBase.autocov!(ghat, x, lags)
+        Ghat = min(Ghat, sum(ghat))
         Ghat > 0 || break
         value += 2 * Ghat
     end
-    return sqrt(value / n)
+
+    mcse = sqrt(value / n)
+
+    return mcse
 end
 
 function mcse_ipse(x::AbstractVector{<:Real})
     n = length(x)
-    m = div(n - 2, 2)
-    x_ = map(Float64, x)
-    ghat = StatsBase.autocov(x_, [0, 1])
-    value = ghat[1] + 2 * ghat[2]
-    for i in 1:m
-        Ghat = sum(StatsBase.autocov(x_, [2 * i, 2 * i + 1]))
+    lags = [0, 1]
+    ghat = StatsBase.autocov(x, lags)
+    @inbounds value = ghat[1] + 2 * ghat[2]
+    @inbounds for i in 2:2:(n - 2)
+        lags[1] = i
+        lags[2] = i + 1
+        StatsBase.autocov!(ghat, x, lags)
+        Ghat = sum(ghat)
         Ghat > 0 || break
         value += 2 * Ghat
     end
-    return sqrt(value / n)
+
+    mcse = sqrt(value / n)
+
+    return mcse
 end