Merge pull request #5 from williamjsdavis/add-turbulence-moments

Add turbulence moments

src/OXBR_turbulence.jl

@@ -0,0 +1,192 @@
+# (Online) Histogram/Kernel Based Regression (OKBR) for turbulence moments
+#NOTES:
+# - No single variants for online methods
+# - No offline kernel methods
+
+
+## Offline algorithms
+
+## Multi step algorithm, 1D
+
+""" Histogram Based Regression moments, Algorithm A """
+function HBR_moments_turb_A(X, tau1_samples, tau2, edge_vector)
+    nτ = length(tau1_samples)
+    nx = length(edge_vector) - 1
+    M1 = zeros(nτ,nx)
+    M2 = zeros(nτ,nx)
+
+    incr2 = X[(tau2+1):end] - X[1:end-tau2]
+    bin_index = map(incr2) do y
+        in_range(edge_vector,y) ? find_bin(edge_vector,y) : 0
+    end 
+    for (i,tau1) in enumerate(tau1_samples)
+        incr1 = X[(tau1+1):(end-(tau2-tau1))] - X[1:(end-tau2)]
+
+        for j in 1:(length(edge_vector)-1)
+            iX = findall(bin_index[1:end] .== j)
+            ΔX = incr1[iX] .- incr2[iX]
+            M1[i,j] = mean(ΔX)
+            residuals = ΔX
+            M2[i,j] = mean(residuals.^2)
+        end
+    end
+    return M1, M2
+end
+
+""" Histogram Based Regression moments, Algorithm C """
+function HBR_moments_turb_C(X, tau1_samples, tau2, edge_vector)
+    nτ = length(tau1_samples)
+    nx = length(edge_vector) - 1
+    nX = length(X)
+    N = zeros(nτ,nx)
+    M1 = zeros(nτ,nx)
+    M2 = zeros(nτ,nx)
+
+    for (i_left, X_left) in enumerate(X[1:end-tau2])
+        incr2 = X[i_left+tau2] - X_left
+        if in_range(edge_vector, incr2)
+            k = find_bin(edge_vector, incr2)
+            for (j,tau1) in enumerate(tau1_samples)
+                ii = i_left + j
+                if ii <= nX
+                    incr1 = X[i_left+tau1] - X_left
+                    ΔX = incr1 - incr2
+                    setindex!(N, N[j,k] + 1, j, k)
+                    setindex!(
+                        M1,
+                        update_mean(M1[j, k], ΔX, N[j, k]),
+                        j, k
+                    )
+                    setindex!(
+                        M2,
+                        update_ss(M2[j, k], ΔX, N[j, k]),
+                        j, k
+                    )
+                end
+            end
+        end
+    end
+
+    return M1, M2
+end
+
+## Online algorithms
+
+## Online Hisogram algrithms
+
+## Multi step algorithm, 1D
+
+mutable struct OHBRu_turb_multiple{T<:AbstractRange}
+    edges::T
+    tau_i::UnitRange{Int}
+    N::Array{Int64,2}
+    M1::Array{Float64,2}
+    M2::Array{Float64,2}
+    mem::Array{Float64,1}
+end
+function OHBRu_turb_multiple(x_range, tau_i)
+    Nx = length(x_range) - 1
+    τ_len = length(tau_i)
+    mem = zeros(Float64, τ_len)
+    mem .= NaN
+    OHBRu_turb_multiple(
+        x_range,
+        tau_i,
+        zeros(Int, τ_len, Nx),
+        zeros(Float64, τ_len, Nx),
+        zeros(Float64, τ_len, Nx),
+        mem
+    )
+end
+
+function add_data!(ohbr::OHBRu_turb_multiple, X_new)
+    X_left = ohbr.mem[end]
+    incr2 = X_new - X_left
+    if in_range(ohbr.edges, incr2)
+        j_bin = find_bin(ohbr.edges, incr2) 
+        for (i_tau, X_right) in enumerate(view(ohbr.mem,1:length(ohbr.mem)-1))
+            incr1 = X_right - X_left
+            ΔX = incr1 - incr2
+            setindex!(ohbr.N, ohbr.N[i_tau,j_bin] + 1, i_tau, j_bin)
+            setindex!(
+                ohbr.M1,
+                update_mean(ohbr.M1[i_tau,j_bin], ΔX, ohbr.N[i_tau,j_bin]),
+                i_tau, j_bin
+            )
+            setindex!(
+                ohbr.M2,
+                update_ss(ohbr.M2[i_tau,j_bin], ΔX, ohbr.N[i_tau,j_bin]),
+                i_tau, j_bin
+            )
+        end
+    end
+    update_mem!(ohbr.mem, X_new)
+    return nothing
+end
+
+## Online Kernel algrithms
+
+## Multi step algorithm, 1D
+
+mutable struct OKBRu_turb_multiple{T<:AbstractRange}
+    x_eval_points::T
+    tau_i::UnitRange{Int}
+    w::Array{Float64,2}
+    M1::Array{Float64,2}
+    M2::Array{Float64,2}
+    mem::Array{Float64,1}
+    w_mem::Array{Float64,1}
+    kernel::Kernel
+    hinv::Float64
+end
+function OKBRu_turb_multiple(x_eval_points, tau_i, kernel, h::Float64)
+    Nx = length(x_eval_points)
+    τ_len = length(tau_i)
+    mem = zeros(Float64, τ_len)
+    mem .= NaN
+    hinv = inv(h)
+    OKBRu_turb_multiple(
+        x_eval_points,
+        tau_i,
+        zeros(Float64, τ_len, Nx),
+        zeros(Float64, τ_len, Nx),
+        zeros(Float64, τ_len, Nx),
+        mem,
+        zeros(Float64, τ_len),
+        kernel,
+        hinv
+    )
+end
+
+function add_data!(okbr::OKBRu_turb_multiple, X_new)
+    X_left = okbr.mem[end]
+    if isnan(X_left)
+        update_mem!(okbr.mem, X_new)
+        return nothing
+    end
+
+    incr2 = X_new - X_left
+    for (j_ind, x_eval) in enumerate(okbr.x_eval_points)
+        K_weight = apply_kernel(x_eval - incr2, okbr.kernel, okbr.hinv)
+        if K_weight > 0.0
+            for (i_tau, X_right) in enumerate(view(okbr.mem,1:length(okbr.mem)-1))
+                incr1 = X_right - X_left
+                ΔX = incr1 - incr2
+                w_old = okbr.w[i_tau,j_ind]
+                setindex!(okbr.w, w_old + K_weight, i_tau, j_ind)
+                setindex!(
+                    okbr.M1,
+                    update_wmean(okbr.M1[i_tau, j_ind], w_old, ΔX, K_weight),
+                    i_tau, j_ind
+                )
+                setindex!(
+                    okbr.M2,
+                    update_wss(okbr.M2[i_tau, j_ind], w_old, ΔX, K_weight),
+                    i_tau, j_ind
+                )
+            end
+        end
+    end
+    update_mem!(okbr.mem, X_new)
+    return nothing
+end

src/OnlineMoments.jl

@@ -37,6 +37,10 @@ export OKBRu_mod_single, OKBRu_mod_multiple
 export OHBR_welford_single
 export M2
 
+# Turbulence moments
+export HBR_moments_turb_A, HBR_moments_turb_C
+export OHBRu_turb_multiple, OKBRu_turb_multiple
+
 # Autocorrelation
 export offline_autocorr, online_autocorr
 export AutoCov
@@ -60,6 +64,7 @@ include("OKBR.jl")
 include("OKBR_uncorrected.jl")
 include("XKBR_mod.jl")
 include("XKBR_mod_uncorrected.jl")
+include("OXBR_turbulence.jl")
 include("autocorr.jl")
 
 end # module

test/runtests.jl

@@ -116,6 +116,23 @@ include("./test_OKBR_uncorrected.jl")
 
 include("./compare_online.jl")
 
+## Testing turbulence moments
+tau1_samples = [4,3,2]
+N_tau1 = length(tau1_samples)
+tau2 = 6
+
+## Histogram methods
+M1_tu_est_A, M2_tu_est_A = HBR_moments_turb_A(X_small, tau1_samples, tau2, x_edges)
+M1_tu_est_C, M2_tu_est_C = HBR_moments_turb_C(X_small, tau1_samples, tau2, x_edges)
+
+## Online Histogram methods
+
+include("./test_OHBR_turbulence.jl")
+
+## Online Kernel methods
+
+include("./test_OKBR_turbulence.jl")
+
 ## Testing autocorrelation
 
 # Offline method

test/test_HBR_turbulence.jl

@@ -0,0 +1,18 @@
+# Testing original algorithms for turbulence calculations
+
+## Histogram Based Regression, 1D (normal methods)
+
+@testset "HBR moments, turbulence" begin
+    @testset "Size" begin
+        @test size(M1_tu_est_A) == (N_tau1, N_x)
+        @test size(M2_tu_est_A) == (N_tau1, N_x)
+        @test size(M1_tu_est_C) == (N_tau1, N_x)
+        @test size(M2_tu_est_C) == (N_tau1, N_x)
+    end
+
+    @testset "Values" begin
+        # Algorithms A and C give almost the same results
+        @test all(M1_tu_est_A .≈ M1_tu_est_C)
+        @test all(M2_tu_est_A .≈ M2_tu_est_C)
+    end
+end

test/test_OHBR_turbulence.jl

@@ -0,0 +1,36 @@
+# Testing streaming algorithms for turbulence calculations
+
+## Online Histogram Based Regression, 1D (online methods)
+
+@testset "OHBR (multiple, turbulence)" begin
+    X_stream = stream_data(X_small)
+
+    # Time-difference sampling
+    tau2_range = 1:tau2
+    N_tau2 = length(tau2_range)
+
+    ohbr_u_multiple = OHBRu_turb_multiple(x_edges, tau2_range)
+    @testset "Structs" begin
+        @test ohbr_u_multiple.edges == x_edges
+        @test ohbr_u_multiple.tau_i == tau_i_range
+        @test size(ohbr_u_multiple.N) == (N_tau2, N_x)
+        @test size(ohbr_u_multiple.M1) == (N_tau2, N_x)
+        @test size(ohbr_u_multiple.M2) == (N_tau2, N_x)
+        @test size(ohbr_u_multiple.mem) == (N_tau2,)
+    end
+
+    for _ in 1:N_data
+        add_data!(ohbr_u_multiple, X_stream())
+    end
+
+    # Comparing with offline methods
+    tau_i = tau2 .- tau1_samples
+
+    @testset "Moments" begin
+        # This streaming algorithm should be almost identical to algorithm A/C
+        @test all(ohbr_u_multiple.M1[tau_i,:] .≈ M1_tu_est_A)
+        @test all(ohbr_u_multiple.M2[tau_i,:] .≈ M2_tu_est_A)
+        @test all(ohbr_u_multiple.M1[tau_i,:] .≈ M1_tu_est_C)
+        @test all(ohbr_u_multiple.M2[tau_i,:] .≈ M2_tu_est_C)
+    end
+end

test/test_OKBR_turbulence.jl

@@ -0,0 +1,35 @@
+# Testing streaming algorithms for turbulence calculations
+
+## Online Kernel Based Regression, 1D (online methods), uncorrected
+
+@testset "OKBR (multiple, turbulence)" begin
+    X_stream = stream_data(X_small)
+
+    # Time-difference sampling
+    tau2_range = 1:tau2
+    N_tau2 = length(tau2_range)
+
+    okbr_u_multiple = OKBRu_turb_multiple(x_centers, tau2_range, kernel_boxcar, h)
+    @testset "Structs" begin
+        @test okbr_u_multiple.x_eval_points == x_centers
+        @test size(okbr_u_multiple.w) == (N_tau, N_x)
+        @test size(okbr_u_multiple.M1) == (N_tau, N_x)
+        @test size(okbr_u_multiple.M2) == (N_tau, N_x)
+        @test size(okbr_u_multiple.mem) == (N_tau,)
+    end
+
+    for _ in 1:N_data
+        add_data!(okbr_u_multiple, X_stream())
+    end
+
+    # Comparing with offline methods
+    tau_i = tau2 .- tau1_samples
+
+    @testset "Moments" begin
+        # This streaming algorithm and HBR algorithm A/C should be almost the same
+        @test all(okbr_u_multiple.M1[tau_i,:] .≈ M1_tu_est_A)
+        @test all(okbr_u_multiple.M2[tau_i,:] .≈ M2_tu_est_A)
+        @test all(okbr_u_multiple.M1[tau_i,:] .≈ M1_tu_est_C)
+        @test all(okbr_u_multiple.M2[tau_i,:] .≈ M2_tu_est_C)
+    end
+end