Block average fixes

CHANGELOG.md

@@ -12,6 +12,7 @@ MolSimToolkit.jl Changelog
 
 Version 1.32.1-DEV
 --------------
+- ![FIX][badge-fix] all BlockAverage plots are now printed with correct time units, and the use of fractional time units (fractional delays between samples) is properly handled.
 
 Version 1.32.0
 --------------

docs/src/block_averages.md

@@ -105,7 +105,7 @@ These functions support `Unitful.jl` quantities both for the response vector and
 ```@example block_averages
 using Unitful
 xu = x .* u"Å" # assign units to the values of x
-b = block_average(xu; dt=1u"ns") # provide the time-step, in ns
+b = block_average(xu; dt=0.0025u"ns") # provide the time delay between data points
 plot(b)
 ```
 

ext/BlockAverages.jl

@@ -72,11 +72,13 @@ function plot(
     xscale=:identity,
     title=""
 )
+    tu = data.dt / oneunit(data.dt)
     l = @layout [ a{0.2h} ; b c ; d e]
     p = plot(layout=l)
     plot!(subplot=1, 
+        collect(1:length(data.x)) * data.dt,
         data.x,
-        xlabel="step",
+        xlabel="time",
         ylabel="value",
         label="",
         color=:black,
@@ -91,9 +93,9 @@ function plot(
         subplot=2,
     )
     plot!(
-        data.blocksize, data.xmean_maxerr,
+        data.dt * data.blocksize, data.xmean_maxerr,
         ylabel="worst block value",
-        xlabel=L"\textrm{block~size}",
+        xlabel="block size",
         label=nothing,
         linewidth=2,
         marker=:circle,
@@ -102,14 +104,14 @@ function plot(
         subplot=2
     )
     annotate!(
-        maximum(data.blocksize) - 0.1 * maximum(data.blocksize),
+        tu * maximum(data.blocksize) - 0.1 * tu * maximum(data.blocksize),
         (max(data.xmean_maxerr[end], maximum(data.xmean_maxerr)) - 0.1 * (maximum(data.xmean_maxerr) - minimum(data.xmean_maxerr))) / oneunit(data.xmean),
         text("mean = $(_round(data.xmean, digits=2))", "Computer Modern", 12, :right),
         subplot=2,
     )
-    plot!(data.blocksize, data.xmean_stderr,
+    plot!(data.dt * data.blocksize, data.xmean_stderr,
         ylabel=L"SD / \sqrt{N_{blocks}}",
-        xlabel=L"\textrm{block~size}",
+        xlabel="block size",
         label=nothing,
         linewidth=2,
         marker=:circle,
@@ -119,8 +121,8 @@ function plot(
     )
     # Auto correlation function
     plot!(
-        data.lags * oneunit(data.tau),
-        data.autocor,
+        data.lags * data.dt,
+        data.autocor * data.dt,
         ylabel=L"c(\Delta t)",
         xlabel=L"\Delta t",
         label=nothing,
@@ -129,20 +131,20 @@ function plot(
         subplot=4
     )
     t95 = 1.96 / sqrt(length(data.x))
-    hline!([t95], subplot=4, ls=:dash, label="", color=:grey)
-    exp_fit = exp.(-inv((data.tau / oneunit(data.tau))) .* data.lags) * oneunit(data.xmean)
+    hline!([t95 * tu], subplot=4, ls=:dash, label="", color=:grey)
+    exp_fit = exp.(-inv((data.tau/oneunit(data.tau))) .* tu .* data.lags ) * oneunit(data.xmean)
     plot!(
-        data.lags * oneunit(data.tau),
-        exp_fit,
+        data.lags * data.dt,
+        exp_fit * data.dt,
         label=nothing,
         linewidth=2,
         color=:black,
         alpha=0.5,
         subplot=4,
     )
     annotate!(
-        (data.lags[end] - 0.2 * data.lags[end]),
-        (0.8 * max(maximum(data.autocor), maximum(exp_fit))) / oneunit(data.xmean),
+        (tu * data.lags[end] - 0.2 * tu * data.lags[end]),
+        (0.8 * tu * max(maximum(data.autocor), maximum(exp_fit))) / oneunit(data.xmean),
         text("τ = $(_round(data.tau; digits=2))", "Computer Modern", 12, :right),
         subplot=4,
     )
@@ -153,8 +155,8 @@ function plot(
         fontfamily="Computer Modern",
         xlims=xlims,
         ylims=ylims,
-        leftmargin=0.1cm,
-        rightmargin=0.1cm,
+        leftmargin=0.3cm,
+        rightmargin=0.3cm,
     )
     plot!(subplot=5,
         xticks=nothing,
@@ -169,11 +171,11 @@ function plot(
     i95 = findfirst(i -> (data.autocor[i] / oneunit(data.autocor[i])) <= t95, eachindex(data.lags))
     isnothing(i95) && (i95 = length(data.lags))
     i95 -= 1
-    plot!((1,1), subplot=5, lc=:white, label="\n"*latexstring("\\Delta t (0.95) = $(data.lags[i95] * oneunit(data.tau))"))
-    plot!((1,1), subplot=5, lc=:white, label=latexstring("\\textrm{Integrated-}\\tau  = $(_round(data.tau_int; digits=4))"))
-    plot!((1,1), subplot=5, lc=:white, label=latexstring("N = $(length(data.x))"))
-    plot!((1,1), subplot=5, lc=:white, label=latexstring("N_{eff}  = $(round(Int, data.n_effective))"))
-    plot!((1,1), subplot=5, lc=:white, label=latexstring("SEM(N_{eff}) = $(_round(data.xmean_stderr_neff; digits=4))"))
+    plot!((1,1), subplot=5, lc=:white, label="\n"*latexstring("\\textrm{\\Delta t (0.95) = $(_round(data.lags[i95] * data.tau; digits=4))}"))
+    plot!((1,1), subplot=5, lc=:white, label=latexstring("\\textrm{Integrated-\\tau = $(_round(data.tau_int; digits=4))}"))
+    plot!((1,1), subplot=5, lc=:white, label=latexstring("\\textrm{N = $(length(data.x))}"))
+    plot!((1,1), subplot=5, lc=:white, label=latexstring("\\textrm{N_{eff}  = $(round(Int, data.n_effective))}"))
+    plot!((1,1), subplot=5, lc=:white, label=latexstring("\\textrm{SEM(N_{eff}) = $(_round(data.xmean_stderr_neff; digits=4))}"))
     return p
 end
 

src/BlockAverages.jl

@@ -50,6 +50,7 @@ struct BlockAverageData{T,DT}
     tau_int::DT
     n_effective::Float64
     xmean_stderr_neff::T
+    dt::DT
 end
 
 function _print_block_sizes(blocksize)
@@ -154,7 +155,7 @@ julia> x = BlockAverages.test_data(10^6); # example data generator
 
 julia> b = block_average(x, lags=0:100:10^5)
 -------------------------------------------------------------------
-BlockAverageData{Float64}
+BlockAverageData
 -------------------------------------------------------------------
 Estimated value (mean by default) = -0.13673023261452855
 Length of data series: 1000000
@@ -246,7 +247,7 @@ function block_average(
     tau_int = 1 + 2 * sum(auto_cor[i] for i in 2:i95-1; init=0.0)
     n_eff = n / tau_int
     xmean_stderr_neff = std(x) / sqrt(n_eff)
-    
+
     return BlockAverageData{TINPUT,DT}(
         x_input,
         xmean * oneunit(TINPUT),
@@ -255,10 +256,11 @@ function block_average(
         xmean_stderr * oneunit(TINPUT),
         lags,
         auto_cor * oneunit(TINPUT),
-        tau * oneunit(dt),
-        tau_int * oneunit(dt),
+        tau * dt,
+        tau_int * dt,
         n_eff,
         xmean_stderr_neff * oneunit(TINPUT),
+        dt,
     )
 
 end
@@ -427,8 +429,9 @@ end # module BlockAverage
     @test_throws "block_size not" block_distribution(sin.(range(0.0, 10.0; length=10)); block_size=-1)
 
     # Test output with units and the definition of dt
-    x = x .* 1u"cm"
-    @test parse_show(block_average(x; dt=1u"s"); repl=Dict("MolSimToolkit." => "", "BlockAverages." => "")) ≈
+    xu = x .* 1u"cm"
+    b1 = block_average(xu; dt=1u"s")
+    @test parse_show(b1; repl=Dict("MolSimToolkit." => "", "BlockAverages." => "")) ≈
           """
           -------------------------------------------------------------------
           BlockAverageData
@@ -453,7 +456,7 @@ end # module BlockAverage
           -------------------------------------------------------------------
           """ float_match = (a, b) -> isapprox(a, b; rtol=0.1)
 
-    @test parse_show(block_distribution(x; block_size=2); repl=Dict("MolSimToolkit." => "", "BlockAverages." => "")) ≈
+    @test parse_show(block_distribution(xu; block_size=2); repl=Dict("MolSimToolkit." => "", "BlockAverages." => "")) ≈
           """
           -------------------------------------------------------------------
           BlockDistribution
@@ -466,4 +469,20 @@ end # module BlockAverage
           -------------------------------------------------------------------
           """
 
+    # Test output with fractional unit of dt
+    xu = x .* 1u"cm"
+    b2 = block_average(xu; dt=0.25u"s")
+    @test b1.autocor ≈ b2.autocor
+    @test b1.n_effective ≈ b2.n_effective
+    @test b1.xmean ≈ b2.xmean
+    @test b1.tau ≈ 4 * b2.tau rtol=0.05
+    @test b1.tau_int ≈ 4 * b2.tau_int rtol=0.05
+
+    b3 = block_average(xu; dt=2u"s")
+    @test b1.autocor ≈ b2.autocor
+    @test b1.n_effective ≈ b2.n_effective
+    @test b1.xmean ≈ b2.xmean
+    @test b3.tau ≈ 8 * b2.tau rtol=0.05
+    @test b3.tau_int ≈ 8 * b2.tau_int rtol=0.05
+
 end