From aa2efb6c5e2aa98e160afe265aef2fff7b8fa74f Mon Sep 17 00:00:00 2001 From: williamjsdavis Date: Sat, 1 Jul 2023 02:21:41 -0700 Subject: [PATCH 1/2] Added HBR + OHBR turbulence methods --- src/OXBR_turbulence.jl | 192 +++++++++++++++++++++++++++++++++++ src/OnlineMoments.jl | 5 + test/runtests.jl | 13 +++ test/test_HBR_turbulence.jl | 18 ++++ test/test_OHBR_turbulence.jl | 35 +++++++ 5 files changed, 263 insertions(+) create mode 100644 src/OXBR_turbulence.jl create mode 100644 test/test_HBR_turbulence.jl create mode 100644 test/test_OHBR_turbulence.jl diff --git a/src/OXBR_turbulence.jl b/src/OXBR_turbulence.jl new file mode 100644 index 0000000..a674da6 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/src/OnlineMoments.jl b/src/OnlineMoments.jl index a7a3212..8913b59 100644 --- a/src/OnlineMoments.jl +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index dc923cf..bb5fae4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -116,6 +116,19 @@ 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") + ## Testing autocorrelation # Offline method diff --git a/test/test_HBR_turbulence.jl b/test/test_HBR_turbulence.jl new file mode 100644 index 0000000..f1b3dc7 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/test/test_OHBR_turbulence.jl b/test/test_OHBR_turbulence.jl new file mode 100644 index 0000000..726d8dc --- /dev/null +++ b/test/test_OHBR_turbulence.jl @@ -0,0 +1,35 @@ +# Testing streaming algorithms for turbulence calculations + +## Online Histogram Based Regression, 1D (online methods) + +@testset "OHBR (multiple, turbulence)" begin + X_stream = stream_data(X_small) + + 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 From 6a9ea124cfc8444960c244ab11d990c959c143d0 Mon Sep 17 00:00:00 2001 From: williamjsdavis Date: Sat, 1 Jul 2023 02:28:54 -0700 Subject: [PATCH 2/2] Added OKBR turbulence methods --- test/runtests.jl | 4 ++++ test/test_OHBR_turbulence.jl | 1 + test/test_OKBR_turbulence.jl | 35 +++++++++++++++++++++++++++++++++++ 3 files changed, 40 insertions(+) create mode 100644 test/test_OKBR_turbulence.jl diff --git a/test/runtests.jl b/test/runtests.jl index bb5fae4..9cb2652 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -129,6 +129,10 @@ M1_tu_est_C, M2_tu_est_C = HBR_moments_turb_C(X_small, tau1_samples, tau2, x_edg include("./test_OHBR_turbulence.jl") +## Online Kernel methods + +include("./test_OKBR_turbulence.jl") + ## Testing autocorrelation # Offline method diff --git a/test/test_OHBR_turbulence.jl b/test/test_OHBR_turbulence.jl index 726d8dc..197e36a 100644 --- a/test/test_OHBR_turbulence.jl +++ b/test/test_OHBR_turbulence.jl @@ -5,6 +5,7 @@ @testset "OHBR (multiple, turbulence)" begin X_stream = stream_data(X_small) + # Time-difference sampling tau2_range = 1:tau2 N_tau2 = length(tau2_range) diff --git a/test/test_OKBR_turbulence.jl b/test/test_OKBR_turbulence.jl new file mode 100644 index 0000000..ec12b48 --- /dev/null +++ b/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