Add PreBinner

README.md

@@ -84,6 +84,12 @@ g_mean = mean(ep, f)
 Δg = std_error(ep, grad_f)
 ```
 
+### PreBinner
+
+`PreBinner(binner, K)` is a wrapper for other binners. It acts as an initial averaging stage, calculating averages of K pushed values before forwarding them to wrapper binner. This removes are certain number of low binning levels (With the unbinned statistics still present in PreBinner) reducing overall memory requirements. 
+
+`PreBinner` follows the common interface with a few adjustments. As examplified above the constructor takes another binner and there is no method for directly processing a time series. The interpretation of binning levels/bin sizes is slightly adjusted, so that 1 refers to unbinned statistics of the PreBinner and 2..N+1 to the levels/bin sizes 1..N in the wrapper binner. 
+
 ### Interface
 
 You can construct each binner with or without initial data:

src/BinningAnalysis.jl

@@ -9,7 +9,7 @@ include("generic.jl")
 export std_error, tau
 
 
-# Logarithmic binning
+# Logarithmic binning (keeping log(N) data points)
 include("log/accumulators.jl")
 include("log/binning.jl")
 include("log/statistics.jl")
@@ -19,27 +19,37 @@ export all_vars, all_varNs, all_taus, all_std_errors, all_autocorrelations, all_
 export convergence, has_converged
 
 
-# "Full" binning
+# Full binning (keeping all data points)
 include("full/binning.jl")
 include("full/statistics.jl")
 export FullBinner
 
 
-# Jackknife resampling
+# Jackknife resampling (estimate function result variance via resampling)
 include("Jackknife.jl")
 @reexport using .Jackknife
 
 
+# log binning with covariance matrix to estimate variance of function result
 include("ErrorPropagation/binning.jl")
 include("ErrorPropagation/statistics.jl")
 export ErrorPropagator
 export means, vars, varNs, taus, std_errors, covmat, all_means
 
 
+# temporal binning with increasing bin size for estimating convergence
 include("incrementing/IncrementBinner.jl")
 export IncrementBinner
 export indices
 
+
+# (expensive) autocorrelation estimation via progression of correlation
 include("direct_tau.jl")
 
+
+# Wrapper for Binners which accumulates N-element means before passing to binners 
+# (compression of low binning levels)
+include("PreBinner.jl")
+export PreBinner
+
 end # module

src/PreBinner.jl

@@ -0,0 +1,131 @@
+mutable struct PreBinner{T, BT, AT} <: AbstractBinner{T}
+    binner::BT
+    buffer::T
+    accumulator::AT
+    n::Int # 0..N-1, push and reset on N
+    N::Int
+end
+
+function PreBinner(B::AbstractBinner{T}, N::Int) where {T}
+    N > 0 || error("Must average at least 1 value.")
+
+    if T <: Number
+        buffer = zero(T)
+    else
+        dummy = mean(B)
+        buffer = zeros(eltype(dummy), size(dummy))
+    end
+
+    return PreBinner(B, buffer, Variance{T}(deepcopy(buffer)), 0, N)
+end
+
+################################################################################
+### Util
+################################################################################
+
+function Base.count(B::PreBinner{T, <: LogBinner}, lvl::Integer = 1) where T
+    return B.N * count(B.binner, max(1, lvl-1)) + (lvl == 1) * B.n
+end
+Base.length(B::PreBinner) = B.N * length(B.binner) + B.n
+Base.ndims(B::PreBinner) = ndims(B.binner)
+Base.isempty(B::PreBinner) = length(B) == 0
+
+function Base.isapprox(a::PreBinner{T, BT}, b::PreBinner{T, BT}; kwargs...) where {T, BT}
+    ((a.n == b.n) && (a.N == b.N)) || return false
+    ((a.n == 0) || isapprox(a.accumulator, b.accumulator; kwargs...)) || return false
+    return isapprox(a.binner, b.binner; kwargs...)
+end
+
+function Base.:(==)(a::PreBinner{T, BT}, b::PreBinner{T, BT}) where {T, BT}
+    ((a.n == b.n) && (a.N == b.N)) || return false
+    ((a.n == 0) || (a.accumulator == b.accumulator)) || return false
+    return a.binner == b.binner
+end
+
+Base.:(!=)(a::PreBinner, b::PreBinner) = !(a == b)
+
+
+# short version (shows up in arrays etc.)
+function Base.show(io::IO, B::PreBinner{T, BT}) where {T, BT}
+    print(io, "PreBinner{$(T), $(BT.name.name)}()")
+end
+# verbose version (shows up in the REPL)
+function Base.show(io::IO, ::MIME"text/plain", B::PreBinner)
+    print(io, "PreBinned ($(B.N)) ")
+    _print_header(io, B.binner)
+
+    n = length(B.binner)
+    println(io)
+    print(io, "| Count: ", B.N, " × ", n, " + ", B.n, " (", length(B), ")")
+    if n > 0 && ndims(B) == 0
+        print(io, "\n| Mean: ", round.(mean(B), digits=5))
+        print(io, "\n| StdError: ", round.(std_error(B), digits=5))
+    end
+    nothing
+end
+
+
+function Base.empty!(B::PreBinner)
+    empty!(B.binner)
+    B.n = 0
+    return B
+end
+
+nlevels(B::PreBinner) = 1 + nlevels(B.binner)
+capacity(B::PreBinner) = B.N * capacity(B.binner) + B.N - 1
+
+################################################################################
+### Pushing
+################################################################################
+
+function Base.push!(B::PreBinner{T,N}, value::S) where {N, T, S}
+    ndims(T) == ndims(S) || throw(DimensionMismatch("Expected $(ndims(T)) dimensions but got $(ndims(S))."))
+
+    # 0/N -> 1 write
+    # N-1 -> N add, div, push
+    # else add
+    n = B.n % B.N
+    _push!(B.accumulator, value)
+    if n == 0
+        _set!(B, value)
+        B.n = 0
+    elseif n == B.N - 1
+        _push!(B, value)
+    else
+        _add!(B, value)
+    end
+    B.n += 1
+    return nothing
+end
+
+_set!(B::PreBinner, value::Number) = B.buffer = value
+_set!(B::PreBinner, value::AbstractArray) = B.buffer .= value
+_add!(B::PreBinner, value::Number) = B.buffer += value
+_add!(B::PreBinner, value::AbstractArray) = B.buffer .+= value
+function _push!(B::PreBinner, value::Number)
+    B.buffer += value
+    push!(B.binner, B.buffer / B.N)
+end
+function _push!(B::PreBinner, value::AbstractArray)
+    B.buffer .+= value
+    B.buffer ./= B.N
+    push!(B.binner, B.buffer)
+end
+
+################################################################################
+### Statistics
+################################################################################
+
+function mean(B::PreBinner, lvl = 1)
+    lvl == 1 ? mean(B.accumulator) : mean(B.binner, max(1, lvl-1))
+end
+
+_reliable_level(B::PreBinner) = 1 + _reliable_level(B.binner)
+
+function var(B::PreBinner, lvl =  _reliable_level(B))
+    lvl == 1 ? var(B.accumulator) : var(B.binner, max(1, lvl-1))
+end
+
+function varN(B::PreBinner, lvl = _reliable_level(B))
+    lvl == 1 ? varN(B.accumulator) : varN(B.binner, max(1, lvl-1))
+end

src/full/binning.jl

@@ -6,7 +6,7 @@ end
 FullBinner(::Type{T} = Float64) where T = FullBinner(Vector{T}(undef, 0))
 
 
-@forward FullBinner.x (Base.length, Base.size, Base.lastindex, Base.ndims, Base.iterate,
+@forward FullBinner.x (Base.length, Base.lastindex,  Base.iterate,
                         Base.getindex, Base.setindex!, Base.view, Base.axes,
                         Base.resize!, Base.sizehint!, Base.empty!, Base.isempty)
 
@@ -16,7 +16,8 @@ function Base.isapprox(a::FullBinner, b::FullBinner; kwargs...)
     (length(a) == length(b)) && (isempty(a) || isapprox(a.x, b.x; kwargs...))
 end
 
-
+Base.ndims(::FullBinner{T}) where T = ndims(T)
+Base.size(B::FullBinner) = size(first(B.x))
 # Cosmetics
 
 function _print_header(io::IO, B::FullBinner{T,A}) where {T, A}

test/PreBinner.jl

@@ -0,0 +1,88 @@
+@testset "Cosmetics (show, print, etc.)" begin
+    B = PreBinner(LogBinner(), 10)
+
+    # empty binner
+    io = IOBuffer()
+    println(io, B) # compact
+    show(io, MIME"text/plain"(), B) # full
+
+    l = String(take!(io))
+    @test l == "PreBinner{Float64, LogBinner}()\nPreBinned (10) LogBinner{Float64,32}\n| Count: 10 × 0 + 0 (0)"
+    @test length(readlines(io)) == 0
+
+
+    # filled binner
+    StableRNGs.seed!(rng, 123)
+    append!(B, rand(rng, 1000))
+    show(io, MIME"text/plain"(), B)
+
+    l = String(take!(io))
+    @test l == "PreBinned (10) LogBinner{Float64,32}\n| Count: 10 × 100 + 10 (1010)\n| Mean: 0.49387\n| StdError: 0.00946"
+    @test length(readlines(io)) == 0
+    close(io);
+end
+
+for wrapped in (LogBinner, FullBinner)
+    @testset "Wrapped $wrapped" begin
+        B = PreBinner(wrapped(), 10)
+        B2 = PreBinner(wrapped(), 10)
+
+        @test B.binner isa wrapped
+        @test length(B) == 0
+        if wrapped == LogBinner
+            @test count(B) == 0
+            @test BinningAnalysis.nlevels(B) == 1 + BinningAnalysis.nlevels(B.binner)
+            @test BinningAnalysis.capacity(B) == B.N * BinningAnalysis.capacity(B.binner) + B.N - 1
+        end
+        @test ndims(B) == 0
+        @test isempty(B)
+        @test eltype(B) == Float64
+        @test B == B2
+        @test B ≈ B2
+        @test !(B != B2)
+
+        @test B.N == 10
+        @test B.n == 0
+
+
+        StableRNGs.seed!(rng, 123)
+        xs = rand(rng, 789)
+        append!(B, xs)
+
+        @test !(B == B2)
+        @test !(B ≈ B2)
+        @test B != B2
+        @test B.n == 9
+        @test length(B) == 789
+        @test length(B.binner) == 78
+
+        @test mean(B) ≈ mean(xs) # Float errors
+        @test std_error(B) == std_error(B.binner)
+        @test var(B) == var(B.binner)
+        @test BinningAnalysis._reliable_level(B) == 1 + BinningAnalysis._reliable_level(B.binner)
+
+
+        empty!(B)
+
+        @test B == B2
+        @test B.n == 0
+        @test length(B) == 0
+        @test length(B.binner) == 0
+    end
+end
+
+@testset "Array elements" begin
+    B = PreBinner(LogBinner(zeros(3, 3)), 7)
+    StableRNGs.seed!(rng, 123)
+    xs = [rand(rng, 3, 3) for _ in 1:100]
+    append!(B, xs)
+
+    @test B.n == 100 % 7
+    @test length(B) == 100
+    @test length(B.binner) == div(100, 7)
+
+    @test mean(B) ≈ mean(xs) # Float errors
+    @test std_error(B) == std_error(B.binner)
+    @test var(B) == var(B.binner)
+    @test BinningAnalysis._reliable_level(B) == 1 + BinningAnalysis._reliable_level(B.binner)
+end

test/fullbinning.jl

@@ -4,9 +4,9 @@
         F2 = FullBinner()
         @test typeof(F) <: BinningAnalysis.AbstractBinner{Float64}
         @test eltype(F) == Float64
-        @test ndims(F) == 1
+        @test ndims(F) == 0
         @test length(F) == 0
-        @test size(F) == (0,)
+        @test_throws BoundsError size(F)
         @test lastindex(F) == 0
         @test axes(F) == (Base.OneTo(0),)
         @test isempty(F)

test/runtests.jl

@@ -10,12 +10,10 @@ const rng = StableRNG(123)
         include("logbinning.jl")
     end
 
-
     @testset "Full Binning" begin
         include("fullbinning.jl")
     end
 
-
     @testset "Error Propagator" begin
         include("errorpropagator.jl")
     end
@@ -26,6 +24,10 @@ const rng = StableRNG(123)
 
     include("systematic_correlation.jl")
 
+    @testset "PreBinner" begin
+        include("PreBinner.jl")
+    end
+
     @testset "Generic functions" begin
         x = [0.0561823, 0.846613, 0.813439, 0.357134, 0.157445, 0.103298, 0.948842, 0.629425, 0.290206, 0.00695332, 0.869828, 0.949165, 0.897995, 0.916239, 0.457564, 0.349827, 0.398683, 0.264218, 0.72754, 0.934315, 0.666448, 0.134813, 0.364933, 0.829088, 0.256443, 0.595029, 0.172097, 0.241686, 0.489935, 0.239663, 0.391291, 0.00751015, 0.138935, 0.0569876, 0.571786, 0.694996, 0.798602, 0.923308, 0.73978, 0.414774, 0.835145, 0.731303, 0.271647, 0.707796, 0.00348624, 0.0905812, 0.316176, 0.921054, 0.131037, 0.599667, 0.805071, 0.440813, 0.086516, 0.363658, 0.476161, 0.931257, 0.28974, 0.78717, 0.60822, 0.144024, 0.214432, 0.061922, 0.626495, 0.512072, 0.758078, 0.840485, 0.242576, 0.147441, 0.599222, 0.993569, 0.0365044, 0.0983033, 0.713144, 0.422394, 0.480044, 0.968745, 0.518475, 0.431319, 0.4432, 0.526007, 0.612975, 0.468387, 0.262145, 0.888011, 0.105744, 0.325821, 0.769525, 0.289073, 0.336083, 0.443037, 0.489698, 0.141654, 0.915284, 0.319068, 0.341001, 0.704346, 0.0794996, 0.0412352, 0.70016, 0.0195158]
         @test std_error(x) ≈ std_error(x; method=:log)