From 2b0a8fd34358dc108d91c6c4e56b77bfb138dfb5 Mon Sep 17 00:00:00 2001 From: odunbar Date: Wed, 18 Jun 2025 11:44:24 -0700 Subject: [PATCH 01/25] New utilities to transform data --- src/Utilities.jl | 372 +++++++++++++++++++++++++++++++++++++ test/Utilities/runtests.jl | 35 ++++ 2 files changed, 407 insertions(+) diff --git a/src/Utilities.jl b/src/Utilities.jl index 1d3e34a20..b3c14543d 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -90,4 +90,376 @@ function zscore2orig(Z::AbstractMatrix{FT}, mean::AbstractVector{FT}, std::Abstr return X end +# Data processing tooling: +# 1) Define a struct subset to the abstract types +# 2) Define a constructor and getters +# 3) Define methods: +# - initialize() +# - process_data() +# - reverse_process_data() + +abstract type PairedDataContainerProcessor end # tools that operate on inputs and outputs +abstract type DataContainerProcessor end # tools that operate on only inputs or outputs + +abstract type UnivariateAffineScaling end +abstract type QuartileScaling <: UnivariateAffineScaling end +abstract type MinMaxScaling <: UnivariateAffineScaling end +abstract type ZScoreScaling <: UnivariateAffineScaling end + +# DCProcessors +""" +$(TYPEDEF) + +The AffineScaler{T} will create an encoding of the data_container via affine transformations: +- quartile_scale() : creates `QuartileScaling`, encoding with the transform (x - Q2(x))/(Q3(x) - Q1(x)) +- minmax_scale() : creates `MinMaxScaling`, encoding with the transform (x - min(x))/(max(x) - min(x)) +- zscore_scale() : creates `ZScoreScaling`, encoding with the (univariate) transform (x-μ)/σ +""" +struct AffineScaler{T, VV <: AbstractVector} <: DataContainerProcessor + "storage for the shift applied to data" + shift::VV + "storage for the scaling" + scale::VV +end + +quartile_scale() = AffineScaler(QuartileScaling) +quartile_scale(T::Type) = AffineScaler(QuartileScaling, T) +minmax_scale() = AffineScaler(MinMaxScaling) +minmax_scale(T::Type) = AffineScaler(MinMaxScaling, T) +zscore_scale() = AffineScaler(ZScoreScaling) +zscore_scale(T::Type) = AffineScaler(ZScoreScaling, T) + +AffineScaler(::Type{UAS}) where {UAS <: UnivariateAffineScaling} = + AffineScaler{UAS, Vector{Float64}}(Float64[], Float64[]) +AffineScaler(::Type{UAS}, T::Type) where {UAS <: UnivariateAffineScaling} = AffineScaler{UAS, Vector{T}}(T[], T[]) + +get_type(as::AffineScaler{T, VV}) where {T, VV} = T +get_shift(as::AffineScaler) = as.shift +get_scale(as::AffineScaler) = as.scale + +function Base.show(io::IO, as::AffineScaler) + out = "AffineScaler: $(get_type(as))" + print(io, out) +end + + + +function initialize_processor!(as::AffineScaler, data::MM) where {MM <: AbstractMatrix} + if length(get_shift(as)) == 0 + T = get_type(as) + initialize_processor!(as, data, T) + end +end + +function initialize_processor!( + as::AffineScaler, + data::MM, + T::Type{QS}, +) where {MM <: AbstractMatrix, QS <: QuartileScaling} + quartiles_vec = [quantile(dd, [0.25, 0.5, 0.75]) for dd in eachrow(data)] + quartiles_mat = reduce(hcat, quartiles_vec) # 3 rows: Q1, Q2, and Q3 + append!(get_shift(as), quartiles_mat[2, :]) + append!(get_scale(as), (quartiles_mat[3, :] - quartiles_mat[1, :])) +end + +function initialize_processor!( + as::AffineScaler, + data::MM, + T::Type{MMS}, +) where {MM <: AbstractMatrix, MMS <: MinMaxScaling} + minmax_vec = [[minimum(dd), maximum(dd)] for dd in eachrow(data)] + minmax_mat = reduce(hcat, minmax_vec) # 3 rows: Q1, Q2, and Q3 + append!(get_shift(as), minmax_mat[1, :]) + append!(get_scale(as), (minmax_mat[2, :] - minmax_mat[1, :])) +end + +function initialize_processor!( + as::AffineScaler, + data::MM, + T::Type{ZSS}, +) where {MM <: AbstractMatrix, ZSS <: ZScoreScaling} + stat_vec = [[mean(dd), std(dd)] for dd in eachrow(data)] + stat_mat = reduce(hcat, stat_vec) # 3 rows: Q1, Q2, and Q3 + append!(get_shift(as), stat_mat[1, :]) + append!(get_scale(as), stat_mat[2, :]) +end + + +function encode_data(as::AffineScaler, data::MM) where {MM <: AbstractMatrix} + out = deepcopy(data) + for i in 1:size(out, 1) + out[i, :] .-= get_shift(as)[i] + out[i, :] /= get_scale(as)[i] + end + return out +end + +function decode_data(as::AffineScaler, data::MM) where {MM <: AbstractMatrix} + out = deepcopy(data) + for i in 1:size(out, 1) + out[i, :] *= get_scale(as)[i] + out[i, :] .+= get_shift(as)[i] + end + return out +end + +""" +$(TYPEDEF) + +Standardizes the data to a multivariate N(0,I) distribution via `(xcov)^{-1/2}*(x-x_mean)`, along with rank reduction if data is low rank, or if the user provides a rank +""" +struct Standardizer{VV1 <: AbstractVector, VV2 <: AbstractVector, VV3 <: AbstractVector} <: DataContainerProcessor + "user-provided rank for additional truncation of the space. used if rank < rank(data_cov)" + rank::Int + "storage for the data mean" + data_mean::VV1 + "storage for wide-matrix representing `data_cov^{-1/2}`" + data_reduction_mat::VV2 + "storage for tall-matrix representing `data_cov^{1/2}`" + data_inflation_mat::VV3 +end + +standardize() = Standardizer(typemax(Int), Float64[], Any[], Any[]) +standardize(T::Type) = Standardizer(typemax(Int), T[], Any[], Any[]) +standardize(T1::Type, T2::Type, T3::Type) = Standardizer(typemax(Int), T1[], T2[], T3[]) +standardize(rk::Int) = Standardizer(rk, Float64[], Any[], Any[]) +standardize(rk::Int, T::Type) = Standardizer(rk, T[], Any[], Any[]) +standardize(rk::Int, T1::Type, T2::Type, T3::Type) = Standardizer(rk, T1[], T2[], T3[]) + +get_data_mean(ss::Standardizer) = ss.data_mean +get_data_reduction_mat(ss::Standardizer) = ss.data_reduction_mat +get_data_inflation_mat(ss::Standardizer) = ss.data_inflation_mat +get_rank(ss::Standardizer) = ss.rank + +function Base.show(io::IO, ss::Standardizer) + if get_rank(ss) < typemax(Int) + out = "Standardizer, user-rank=$(get_rank(ss))" + else + out = "Standardizer" + end + print(io, out) +end + +function initialize_processor!(ss::Standardizer, data::MM) where {MM <: AbstractMatrix} + if length(get_data_mean(ss)) == 0 + append!(get_data_mean(ss), mean(data, dims = 2)) + end + if length(get_data_reduction_mat(ss)) == 0 + # can use tsvd here so we only compute up to rk + data_cov = cov(data, dims = 2) + svdc = svd(data_cov) + rk = min(rank(data_cov), max(get_rank(ss), 1)) + # Want to use the nonsquare singular vector + if size(svdc.U, 1) == size(svdc.U, 2) + mat = svdc.Vt + else + mat = svdc.U' + end + + reduction_mat = Diagonal(1 ./ sqrt.(svdc.S[1:rk])) * mat[1:rk, :] + inflation_mat = mat[1:rk, :]' * Diagonal(sqrt.(svdc.S[1:rk])) + push!(get_data_reduction_mat(ss), reduction_mat) + push!(get_data_inflation_mat(ss), inflation_mat) + end +end + +function encode_data(ss::Standardizer, data::MM) where {MM <: AbstractMatrix} + data_mean = get_data_mean(ss) + reduction_mat = get_data_reduction_mat(ss)[1] + return reduction_mat * (data .- data_mean) +end + +function decode_data(ss::Standardizer, data::MM) where {MM <: AbstractMatrix} + data_mean = get_data_mean(ss) + data_mean = get_data_mean(ss) + inflation_mat = get_data_inflation_mat(ss)[1] + return inflation_mat * data .+ data_mean +end + +""" +$(TYPEDEF) + +Decorrelate the data via an SVD decomposition, with optional truncation to space of singular vectors corresponding to largest singular values. + +# Fields +$(TYPEDFIELDS) +""" +struct Decorrelater{MM1, MM2, FT} <: DataContainerProcessor + "the decorrelation matrix - provided by the user" + decorrelation_matrix::MM1 + "the inverse of the decorrelation matrix" + inv_decorrelation_matrix::MM2 + "the fraction of variance to be retained after truncating singular values (1 implies no truncation)" + retained_variance::FT +end + +# ...struct VariationalAutoEncoder <: DataContainerProcessor end + +# PDCProcessors +struct DataInformedReducer <: PairedDataContainerProcessor end + +struct LikelihoodInformedReducer <: PairedDataContainerProcessor end + +struct CanonicalCorrelationReducer <: PairedDataContainerProcessor end + + +#### + + +# generic function to build and encode data +function initialize_and_encode_data!(dcp::DCP, data::MM) where {DCP <: DataContainerProcessor, MM <: AbstractMatrix} + initialize_processor!(dcp, data) + return encode_data(dcp, data) +end + +function decode_data(dcp::DCP, data::MM) where {DCP <: DataContainerProcessor, MM <: AbstractMatrix} + return decode_data(dcp, data) +end + +initialize_and_encode_data!( + dcp::DCP, + data::MM, + apply_to::AS, +) where {DCP <: DataContainerProcessor, MM <: AbstractMatrix, AS <: AbstractString} = + initialize_and_encode_data!(dcp, data) + +decode_data( + dcp::DCP, + data::MM, + apply_to::AS, +) where {DCP <: DataContainerProcessor, MM <: AbstractMatrix, AS <: AbstractString} = decode_data(dcp, data) + +function initialize_and_encode_data!( + dcp::PDCP, + data, + apply_to::AS, +) where {PDCP <: PairedDataContainerProcessor, AS <: AbstractString} + input_data, output_data = data + initialize_processor!(dcp, input_data, output_data, apply_to) + return encode_data(dcp, input_data, output_data, apply_to) +end + +function decode_data(dcp::PDCP, data, apply_to::AS) where {PDCP <: PairedDataContainerProcessor, AS <: AbstractString} + input_data, output_data = data + return decode_data(dcp, input_data, output_data, apply_to) +end + +""" +$TYPEDSIGNATURES + +Create a flatter encoder schedule for the +from the user's proposed schedule of the form: +```julia +enc_schedule = [ + (DataProcessor1(...), "in"), + (DataProcessor2(...), "out"), + (PairedDataProcessor3(...), "in"), + (DataProcessor4(...), "in_and_out"), +] +``` +This function creates the encoder scheduler that is also machine readable +```julia +enc_schedule = [ + (DataProcessor1(...), x -> get_inputs(x), "in"), + (DataProcessor2(...), x -> get_outputs(x), "out"), + (DataProcessor2(...), x -> get_outputs(x), "out"), + (PairedDataProcessor3(...), x -> (get_outputs(x), get_outputs(x)), "in"), + (DataProcessor4(...), x -> get_inputs(x), "in"), + (DataProcessor4(...), x -> get_outputs(x), "out"), +] +``` +and the decoder schedule is a copy of the encoder schedule reversed (and processors copied) +""" +function create_encoder_schedule(schedule_in::VV) where {VV <: AbstractVector} + + encoder_schedule = [] + for (processor, apply_to) in schedule_in + # converts the string into the extraction of data + if isa(processor, DataContainerProcessor) + if apply_to == "in" + func = x -> get_inputs(x) + push!(encoder_schedule, (processor, func, apply_to)) + + elseif apply_to == "out" + func = x -> get_outputs(x) + push!(encoder_schedule, (processor, func, apply_to)) + + elseif apply_to == "in_and_out" + func1 = x -> get_inputs(x) + func2 = x -> get_outputs(x) + push!(encoder_schedule, (processor, func1, "in")) + push!(encoder_schedule, (deepcopy(processor), func2, "out")) + else + @warn( + "Expected schedule keywords ∈ {\"in\",\"out\",\"in_and_out\"}. Received $(str), ignoring processor $(processor)..." + ) + end + + # extract all the data (needed for the paired processor, but still pass through what you apply to) + elseif isa(processor, PairedDataContainerProcessor) + if apply_to ∈ ["in", "out"] + func = x -> (get_inputs(x), get_outputs(x)) + push!(encoder_schedule, (processor, func, apply_to)) + elseif apply_to == "in_and_out" + func = x -> (get_inputs(x), get_outputs(x)) + push!(encoder_schedule, (processor, func, "in")) + push!(encoder_schedule, (deepcopy(processor), func, "out")) + else + @warn( + "Expected schedule keywords ∈ {\"in\",\"out\",\"in_and_out\"}. Received $(str), ignoring processor $(processor)..." + ) + end + end + end + + return encoder_schedule +end + +function encode_with_schedule( + encoder_schedule::VV, + io_pairs::PDC, +) where {VV <: AbstractVector, PDC <: PairedDataContainer} + + processed_io_pairs = deepcopy(io_pairs) + + # apply_to is the string "in", "out" etc. + for (processor, extract_data, apply_to) in encoder_schedule + @info "Encoding data: $(apply_to) with $(processor)" + processed = initialize_and_encode_data!(processor, extract_data(processed_io_pairs), apply_to) + + if apply_to == "in" + processed_io_pairs = PairedDataContainer(processed, get_outputs(processed_io_pairs)) + elseif apply_to == "out" + processed_io_pairs = PairedDataContainer(get_inputs(processed_io_pairs), processed) + end + end + + return processed_io_pairs +end + +function decode_with_schedule( + encoder_schedule::VV, + io_pairs::PDC, +) where {VV <: AbstractVector, PDC <: PairedDataContainer} + + processed_io_pairs = deepcopy(io_pairs) + + # apply_to is the string "in", "out" etc. + for idx in reverse(eachindex(encoder_schedule)) + (processor, extract_data, apply_to) = encoder_schedule[idx] + @info "Decoding data: $(apply_to) with $(processor)" + processed = decode_data(processor, extract_data(processed_io_pairs), apply_to) + + if apply_to == "in" + processed_io_pairs = PairedDataContainer(processed, get_outputs(processed_io_pairs)) + elseif apply_to == "out" + processed_io_pairs = PairedDataContainer(get_inputs(processed_io_pairs), processed) + end + end + + return processed_io_pairs +end + + + end # module diff --git a/test/Utilities/runtests.jl b/test/Utilities/runtests.jl index 7fe46f522..c60c23a5b 100644 --- a/test/Utilities/runtests.jl +++ b/test/Utilities/runtests.jl @@ -69,3 +69,38 @@ using CalibrateEmulateSample.DataContainers @test minimum(eigvals(pdmat2)) >= (1 - 1e-4) * tol end + + +@testset "Data Preprocessing" begin + + # get some data as IO pairs + in_dim = 10 + out_dim = 50 + samples = 100 + + in_data = 20 * randn(in_dim, samples) + out_data = 3 * randn(in_dim, samples) .- 10 + + io_pairs = PairedDataContainer(in_data, out_data) + + # Create a lossless encoding schedule + schedule_builder = [ + (zscore_scale(), "in_and_out"), # + (quartile_scale(), "in"), + (standardize(), "out"), + (minmax_scale(), "in_and_out"), + ] + + # make schedule more parsable + encoder_schedule = create_encoder_schedule(schedule_builder) + + # encode the data using the schedule + encoded_io_pairs = encode_with_schedule(encoder_schedule, io_pairs) + + # decode the data using the schedule + decoded_io_pairs = decode_with_schedule(encoder_schedule, encoded_io_pairs) + + tol = 1e-12 + @test all(isapprox.(get_inputs(io_pairs), get_inputs(decoded_io_pairs), atol = tol)) + @test all(isapprox.(get_outputs(io_pairs), get_outputs(decoded_io_pairs), atol = tol)) +end From e7c6b8ac11616ce5ccf5673c55a625cfa511f22f Mon Sep 17 00:00:00 2001 From: odunbar Date: Wed, 18 Jun 2025 14:12:33 -0700 Subject: [PATCH 02/25] add exports --- src/Utilities.jl | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/src/Utilities.jl b/src/Utilities.jl index b3c14543d..6364a86a1 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -1,5 +1,7 @@ module Utilities + + using DocStringExtensions using LinearAlgebra using Statistics @@ -12,6 +14,17 @@ using ..DataContainers export get_training_points export orig2zscore export zscore2orig + +export PairedDataContainerProcessor,DataContainerProcessor +export UnivariateAffineScaling, QuartileScaling, MinMaxScaling, ZScoreScaling, Standardizer + +export quartile_scale, minmax_scale, zscore_scale, standardize +export get_type, get_shift, get_scale, get_data_mean, get_data_reduction_mat, get_data_inflation_mat, get_rank +export create_encoder_schedule, encode_with_schedule, decode_with_schedule +export initialize_processor!, initialize_and_encode_data!, encode_data, decode_data + + + """ $(DocStringExtensions.TYPEDSIGNATURES) @@ -203,6 +216,17 @@ function decode_data(as::AffineScaler, data::MM) where {MM <: AbstractMatrix} return out end +function encode_matrix(as::AffineScaler, data::MM) where {MM <: AbstractMatrix} + out = deepcopy(data) + for i in 1:size(out, 1) + out[i, :] .-= get_shift(as)[i] + out[i, :] /= get_scale(as)[i] + end + return out +end + + + """ $(TYPEDEF) From d638860f5eccb55e5d8fc5a2808754876a7efd08 Mon Sep 17 00:00:00 2001 From: odunbar Date: Wed, 18 Jun 2025 15:18:59 -0700 Subject: [PATCH 03/25] added obs-noise-cov encoding/decoding alongside data processing tools --- src/Utilities.jl | 136 +++++++++++++++++++++++-------------- test/Utilities/runtests.jl | 11 ++- 2 files changed, 94 insertions(+), 53 deletions(-) diff --git a/src/Utilities.jl b/src/Utilities.jl index 6364a86a1..7c620ae8c 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -15,13 +15,14 @@ export get_training_points export orig2zscore export zscore2orig -export PairedDataContainerProcessor,DataContainerProcessor -export UnivariateAffineScaling, QuartileScaling, MinMaxScaling, ZScoreScaling, Standardizer +export PairedDataContainerProcessor, DataContainerProcessor +export UnivariateAffineScaling, QuartileScaling, MinMaxScaling, ZScoreScaling, Standardizer export quartile_scale, minmax_scale, zscore_scale, standardize export get_type, get_shift, get_scale, get_data_mean, get_data_reduction_mat, get_data_inflation_mat, get_rank export create_encoder_schedule, encode_with_schedule, decode_with_schedule -export initialize_processor!, initialize_and_encode_data!, encode_data, decode_data +export initialize_processor!, + initialize_and_encode_data!, encode_data, decode_data, encode_obs_noise_cov, decode_obs_noise_cov @@ -104,12 +105,6 @@ function zscore2orig(Z::AbstractMatrix{FT}, mean::AbstractVector{FT}, std::Abstr end # Data processing tooling: -# 1) Define a struct subset to the abstract types -# 2) Define a constructor and getters -# 3) Define methods: -# - initialize() -# - process_data() -# - reverse_process_data() abstract type PairedDataContainerProcessor end # tools that operate on inputs and outputs abstract type DataContainerProcessor end # tools that operate on only inputs or outputs @@ -119,7 +114,7 @@ abstract type QuartileScaling <: UnivariateAffineScaling end abstract type MinMaxScaling <: UnivariateAffineScaling end abstract type ZScoreScaling <: UnivariateAffineScaling end -# DCProcessors +# Processors """ $(TYPEDEF) @@ -155,15 +150,6 @@ function Base.show(io::IO, as::AffineScaler) print(io, out) end - - -function initialize_processor!(as::AffineScaler, data::MM) where {MM <: AbstractMatrix} - if length(get_shift(as)) == 0 - T = get_type(as) - initialize_processor!(as, data, T) - end -end - function initialize_processor!( as::AffineScaler, data::MM, @@ -197,6 +183,12 @@ function initialize_processor!( append!(get_scale(as), stat_mat[2, :]) end +function initialize_processor!(as::AffineScaler, data::MM) where {MM <: AbstractMatrix} + if length(get_shift(as)) == 0 + T = get_type(as) + initialize_processor!(as, data, T) + end +end function encode_data(as::AffineScaler, data::MM) where {MM <: AbstractMatrix} out = deepcopy(data) @@ -216,13 +208,17 @@ function decode_data(as::AffineScaler, data::MM) where {MM <: AbstractMatrix} return out end -function encode_matrix(as::AffineScaler, data::MM) where {MM <: AbstractMatrix} - out = deepcopy(data) - for i in 1:size(out, 1) - out[i, :] .-= get_shift(as)[i] - out[i, :] /= get_scale(as)[i] - end - return out +initialize_processor!(as::AffineScaler, data::MM, obs_noise_cov) where {MM <: AbstractMatrix} = + initialize_processor!(as, data) +encode_data(as::AffineScaler, data::MM, obs_noise_cov) where {MM <: AbstractMatrix} = encode_data(as, data) +decode_data(as::AffineScaler, data::MM, obs_noise_cov) where {MM <: AbstractMatrix} = decode_data(as, data) + +function encode_obs_noise_cov(as::AffineScaler, obs_noise_cov::MM) where {MM <: AbstractMatrix} + return Diagonal(1 ./ get_scale(as) .^ 2) * obs_noise_cov +end + +function decode_obs_noise_cov(as::AffineScaler, enc_obs_noise_cov::MM) where {MM <: AbstractMatrix} + return Diagonal(get_scale(as) .^ 2) * enc_obs_noise_cov end @@ -294,16 +290,31 @@ function encode_data(ss::Standardizer, data::MM) where {MM <: AbstractMatrix} end function decode_data(ss::Standardizer, data::MM) where {MM <: AbstractMatrix} - data_mean = get_data_mean(ss) data_mean = get_data_mean(ss) inflation_mat = get_data_inflation_mat(ss)[1] return inflation_mat * data .+ data_mean end +initialize_processor!(ss::Standardizer, data::MM, obs_noise_cov) where {MM <: AbstractMatrix} = + initialize_processor!(ss, data) +encode_data(ss::Standardizer, data::MM, obs_noise_cov) where {MM <: AbstractMatrix} = encode_data(ss, data) +decode_data(ss::Standardizer, data::MM, obs_noise_cov) where {MM <: AbstractMatrix} = decode_data(ss, data) + +function encode_obs_noise_cov(ss::Standardizer, obs_noise_cov::MM) where {MM <: AbstractMatrix} + reduction_mat = get_data_reduction_mat(ss)[1] + return reduction_mat * obs_noise_cov * reduction_mat' +end + +function decode_obs_noise_cov(ss::Standardizer, enc_obs_noise_cov::MM) where {MM <: AbstractMatrix} + inflation_mat = get_data_inflation_mat(ss)[1] + return inflation_mat * enc_obs_noise_cov * inflation_mat' +end + + """ $(TYPEDEF) -Decorrelate the data via an SVD decomposition, with optional truncation to space of singular vectors corresponding to largest singular values. +Decorrelate the data via an SVD decomposition using the observational noise `obs_noise_cov`, with optional truncation of singular vectors corresponding to largest singular values. # Fields $(TYPEDFIELDS) @@ -317,6 +328,8 @@ struct Decorrelater{MM1, MM2, FT} <: DataContainerProcessor retained_variance::FT end + + # ...struct VariationalAutoEncoder <: DataContainerProcessor end # PDCProcessors @@ -331,41 +344,58 @@ struct CanonicalCorrelationReducer <: PairedDataContainerProcessor end # generic function to build and encode data -function initialize_and_encode_data!(dcp::DCP, data::MM) where {DCP <: DataContainerProcessor, MM <: AbstractMatrix} - initialize_processor!(dcp, data) - return encode_data(dcp, data) -end - -function decode_data(dcp::DCP, data::MM) where {DCP <: DataContainerProcessor, MM <: AbstractMatrix} - return decode_data(dcp, data) +function initialize_and_encode_data!( + dcp::DCP, + data::MM, + obs_noise_cov::USorM, +) where {DCP <: DataContainerProcessor, USorM <: Union{UniformScaling, AbstractMatrix}, MM <: AbstractMatrix} + initialize_processor!(dcp, data, obs_noise_cov) + return encode_data(dcp, data, obs_noise_cov) end initialize_and_encode_data!( dcp::DCP, data::MM, + obs_noise_cov::USorM, apply_to::AS, -) where {DCP <: DataContainerProcessor, MM <: AbstractMatrix, AS <: AbstractString} = - initialize_and_encode_data!(dcp, data) +) where { + DCP <: DataContainerProcessor, + MM <: AbstractMatrix, + USorM <: Union{UniformScaling, AbstractMatrix}, + AS <: AbstractString, +} = initialize_and_encode_data!(dcp, data, obs_noise_cov) decode_data( dcp::DCP, data::MM, + obs_noise_cov::USorM, apply_to::AS, -) where {DCP <: DataContainerProcessor, MM <: AbstractMatrix, AS <: AbstractString} = decode_data(dcp, data) +) where { + DCP <: DataContainerProcessor, + MM <: AbstractMatrix, + USorM <: Union{UniformScaling, AbstractMatrix}, + AS <: AbstractString, +} = decode_data(dcp, data, obs_noise_cov) function initialize_and_encode_data!( dcp::PDCP, data, + obs_noise_cov::USorM, apply_to::AS, -) where {PDCP <: PairedDataContainerProcessor, AS <: AbstractString} +) where {PDCP <: PairedDataContainerProcessor, USorM <: Union{UniformScaling, AbstractMatrix}, AS <: AbstractString} input_data, output_data = data - initialize_processor!(dcp, input_data, output_data, apply_to) - return encode_data(dcp, input_data, output_data, apply_to) + initialize_processor!(dcp, input_data, output_data, obs_noise_cov, apply_to) + return encode_data(dcp, input_data, output_data, obs_noise_cov, apply_to) end -function decode_data(dcp::PDCP, data, apply_to::AS) where {PDCP <: PairedDataContainerProcessor, AS <: AbstractString} +function decode_data( + dcp::PDCP, + data, + obs_noise_cov::USorM, + apply_to::AS, +) where {PDCP <: PairedDataContainerProcessor, USorM <: Union{UniformScaling, AbstractMatrix}, AS <: AbstractString} input_data, output_data = data - return decode_data(dcp, input_data, output_data, apply_to) + return decode_data(dcp, input_data, output_data, obs_noise_cov, apply_to) end """ @@ -442,46 +472,52 @@ end function encode_with_schedule( encoder_schedule::VV, io_pairs::PDC, -) where {VV <: AbstractVector, PDC <: PairedDataContainer} + obs_noise_cov_in::USorM, +) where {VV <: AbstractVector, PDC <: PairedDataContainer, USorM <: Union{UniformScaling, AbstractMatrix}} processed_io_pairs = deepcopy(io_pairs) + processed_obs_noise_cov = deepcopy(obs_noise_cov_in) # apply_to is the string "in", "out" etc. for (processor, extract_data, apply_to) in encoder_schedule @info "Encoding data: $(apply_to) with $(processor)" - processed = initialize_and_encode_data!(processor, extract_data(processed_io_pairs), apply_to) - + processed = + initialize_and_encode_data!(processor, extract_data(processed_io_pairs), processed_obs_noise_cov, apply_to) if apply_to == "in" processed_io_pairs = PairedDataContainer(processed, get_outputs(processed_io_pairs)) elseif apply_to == "out" + processed_obs_noise_cov = encode_obs_noise_cov(processor, processed_obs_noise_cov) processed_io_pairs = PairedDataContainer(get_inputs(processed_io_pairs), processed) end end - return processed_io_pairs + return processed_io_pairs, processed_obs_noise_cov end function decode_with_schedule( encoder_schedule::VV, io_pairs::PDC, -) where {VV <: AbstractVector, PDC <: PairedDataContainer} + obs_noise_cov_in::USorM, +) where {VV <: AbstractVector, USorM <: Union{UniformScaling, AbstractMatrix}, PDC <: PairedDataContainer} processed_io_pairs = deepcopy(io_pairs) + processed_obs_noise_cov = deepcopy(obs_noise_cov_in) # apply_to is the string "in", "out" etc. for idx in reverse(eachindex(encoder_schedule)) (processor, extract_data, apply_to) = encoder_schedule[idx] @info "Decoding data: $(apply_to) with $(processor)" - processed = decode_data(processor, extract_data(processed_io_pairs), apply_to) + processed = decode_data(processor, extract_data(processed_io_pairs), processed_obs_noise_cov, apply_to) if apply_to == "in" processed_io_pairs = PairedDataContainer(processed, get_outputs(processed_io_pairs)) elseif apply_to == "out" + processed_obs_noise_cov = decode_obs_noise_cov(processor, processed_obs_noise_cov) processed_io_pairs = PairedDataContainer(get_inputs(processed_io_pairs), processed) end end - return processed_io_pairs + return processed_io_pairs, processed_obs_noise_cov end diff --git a/test/Utilities/runtests.jl b/test/Utilities/runtests.jl index c60c23a5b..7caeab393 100644 --- a/test/Utilities/runtests.jl +++ b/test/Utilities/runtests.jl @@ -1,6 +1,7 @@ using Test using Random using Statistics +using Distributions using LinearAlgebra using CalibrateEmulateSample.Utilities @@ -79,7 +80,8 @@ end samples = 100 in_data = 20 * randn(in_dim, samples) - out_data = 3 * randn(in_dim, samples) .- 10 + obs_noise_cov = [max(5.0 - abs(i - j), 0.0) for i in 1:out_dim, j in 1:out_dim] # [5 4 3 2 1 0 0 ...] off diagonal + out_data = rand(MvNormal(-10 * ones(out_dim), obs_noise_cov), samples) .- 10 io_pairs = PairedDataContainer(in_data, out_data) @@ -95,12 +97,15 @@ end encoder_schedule = create_encoder_schedule(schedule_builder) # encode the data using the schedule - encoded_io_pairs = encode_with_schedule(encoder_schedule, io_pairs) + encoded_io_pairs, encoded_obs_noise_cov = encode_with_schedule(encoder_schedule, io_pairs, obs_noise_cov) # decode the data using the schedule - decoded_io_pairs = decode_with_schedule(encoder_schedule, encoded_io_pairs) + decoded_io_pairs, decoded_obs_noise_cov = + decode_with_schedule(encoder_schedule, encoded_io_pairs, encoded_obs_noise_cov) tol = 1e-12 @test all(isapprox.(get_inputs(io_pairs), get_inputs(decoded_io_pairs), atol = tol)) @test all(isapprox.(get_outputs(io_pairs), get_outputs(decoded_io_pairs), atol = tol)) + @test all(isapprox.(obs_noise_cov, decoded_obs_noise_cov, atol = tol)) + end From 3a559b06fc3d8576406ff6c43e2ef637e5f3db08 Mon Sep 17 00:00:00 2001 From: odunbar Date: Wed, 18 Jun 2025 16:16:20 -0700 Subject: [PATCH 04/25] prior cov and obs cov as structure matrices --- src/Utilities.jl | 87 +++++++++++++++++++++++++------------- test/Utilities/runtests.jl | 18 +++++--- 2 files changed, 68 insertions(+), 37 deletions(-) diff --git a/src/Utilities.jl b/src/Utilities.jl index 7c620ae8c..8b8189a22 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -208,17 +208,17 @@ function decode_data(as::AffineScaler, data::MM) where {MM <: AbstractMatrix} return out end -initialize_processor!(as::AffineScaler, data::MM, obs_noise_cov) where {MM <: AbstractMatrix} = +initialize_processor!(as::AffineScaler, data::MM, structure_matrix) where {MM <: AbstractMatrix} = initialize_processor!(as, data) -encode_data(as::AffineScaler, data::MM, obs_noise_cov) where {MM <: AbstractMatrix} = encode_data(as, data) -decode_data(as::AffineScaler, data::MM, obs_noise_cov) where {MM <: AbstractMatrix} = decode_data(as, data) +encode_data(as::AffineScaler, data::MM, structure_matrix) where {MM <: AbstractMatrix} = encode_data(as, data) +decode_data(as::AffineScaler, data::MM, structure_matrix) where {MM <: AbstractMatrix} = decode_data(as, data) -function encode_obs_noise_cov(as::AffineScaler, obs_noise_cov::MM) where {MM <: AbstractMatrix} - return Diagonal(1 ./ get_scale(as) .^ 2) * obs_noise_cov +function encode_structure_matrix(as::AffineScaler, structure_matrix::MM) where {MM <: AbstractMatrix} + return Diagonal(1 ./ get_scale(as) .^ 2) * structure_matrix end -function decode_obs_noise_cov(as::AffineScaler, enc_obs_noise_cov::MM) where {MM <: AbstractMatrix} - return Diagonal(get_scale(as) .^ 2) * enc_obs_noise_cov +function decode_structure_matrix(as::AffineScaler, enc_structure_matrix::MM) where {MM <: AbstractMatrix} + return Diagonal(get_scale(as) .^ 2) * enc_structure_matrix end @@ -295,35 +295,35 @@ function decode_data(ss::Standardizer, data::MM) where {MM <: AbstractMatrix} return inflation_mat * data .+ data_mean end -initialize_processor!(ss::Standardizer, data::MM, obs_noise_cov) where {MM <: AbstractMatrix} = +initialize_processor!(ss::Standardizer, data::MM, structure_matrix) where {MM <: AbstractMatrix} = initialize_processor!(ss, data) -encode_data(ss::Standardizer, data::MM, obs_noise_cov) where {MM <: AbstractMatrix} = encode_data(ss, data) -decode_data(ss::Standardizer, data::MM, obs_noise_cov) where {MM <: AbstractMatrix} = decode_data(ss, data) +encode_data(ss::Standardizer, data::MM, structure_matrix) where {MM <: AbstractMatrix} = encode_data(ss, data) +decode_data(ss::Standardizer, data::MM, structure_matrix) where {MM <: AbstractMatrix} = decode_data(ss, data) -function encode_obs_noise_cov(ss::Standardizer, obs_noise_cov::MM) where {MM <: AbstractMatrix} +function encode_structure_matrix(ss::Standardizer, structure_matrix::MM) where {MM <: AbstractMatrix} reduction_mat = get_data_reduction_mat(ss)[1] - return reduction_mat * obs_noise_cov * reduction_mat' + return reduction_mat * structure_matrix * reduction_mat' end -function decode_obs_noise_cov(ss::Standardizer, enc_obs_noise_cov::MM) where {MM <: AbstractMatrix} +function decode_structure_matrix(ss::Standardizer, enc_structure_matrix::MM) where {MM <: AbstractMatrix} inflation_mat = get_data_inflation_mat(ss)[1] - return inflation_mat * enc_obs_noise_cov * inflation_mat' + return inflation_mat * enc_structure_matrix * inflation_mat' end """ $(TYPEDEF) -Decorrelate the data via an SVD decomposition using the observational noise `obs_noise_cov`, with optional truncation of singular vectors corresponding to largest singular values. +Decorrelate the data via an SVD decomposition using a structure_matrix (`prior_cov` for inputs, `obs_noise_cov` for outputs), with optional truncation of singular vectors corresponding to largest singular values. # Fields $(TYPEDFIELDS) """ struct Decorrelater{MM1, MM2, FT} <: DataContainerProcessor - "the decorrelation matrix - provided by the user" - decorrelation_matrix::MM1 - "the inverse of the decorrelation matrix" - inv_decorrelation_matrix::MM2 + "the structure matrix - provided by the user" + structure_matrix::MM1 + "the inverse of the structure matrix" + inv_structure_matrix::MM2 "the fraction of variance to be retained after truncating singular values (1 implies no truncation)" retained_variance::FT end @@ -472,54 +472,81 @@ end function encode_with_schedule( encoder_schedule::VV, io_pairs::PDC, - obs_noise_cov_in::USorM, -) where {VV <: AbstractVector, PDC <: PairedDataContainer, USorM <: Union{UniformScaling, AbstractMatrix}} + prior_cov_in::USorM1, + obs_noise_cov_in::USorM2, +) where { + VV <: AbstractVector, + PDC <: PairedDataContainer, + USorM1 <: Union{UniformScaling, AbstractMatrix}, + USorM2 <: Union{UniformScaling, AbstractMatrix}, +} processed_io_pairs = deepcopy(io_pairs) + processed_prior_cov = deepcopy(prior_cov_in) processed_obs_noise_cov = deepcopy(obs_noise_cov_in) # apply_to is the string "in", "out" etc. for (processor, extract_data, apply_to) in encoder_schedule @info "Encoding data: $(apply_to) with $(processor)" - processed = - initialize_and_encode_data!(processor, extract_data(processed_io_pairs), processed_obs_noise_cov, apply_to) if apply_to == "in" + structure_matrix = processed_prior_cov + elseif apply_to == "out" + structure_matrix = processed_obs_noise_cov + end + processed = initialize_and_encode_data!(processor, extract_data(processed_io_pairs), structure_matrix, apply_to) + if apply_to == "in" + processed_prior_cov = encode_structure_matrix(processor, structure_matrix) processed_io_pairs = PairedDataContainer(processed, get_outputs(processed_io_pairs)) elseif apply_to == "out" - processed_obs_noise_cov = encode_obs_noise_cov(processor, processed_obs_noise_cov) + processed_obs_noise_cov = encode_structure_matrix(processor, structure_matrix) processed_io_pairs = PairedDataContainer(get_inputs(processed_io_pairs), processed) end end - return processed_io_pairs, processed_obs_noise_cov + return processed_io_pairs, processed_prior_cov, processed_obs_noise_cov end function decode_with_schedule( encoder_schedule::VV, io_pairs::PDC, - obs_noise_cov_in::USorM, -) where {VV <: AbstractVector, USorM <: Union{UniformScaling, AbstractMatrix}, PDC <: PairedDataContainer} + prior_cov_in::USorM1, + obs_noise_cov_in::USorM2, +) where { + VV <: AbstractVector, + USorM1 <: Union{UniformScaling, AbstractMatrix}, + USorM2 <: Union{UniformScaling, AbstractMatrix}, + PDC <: PairedDataContainer, +} processed_io_pairs = deepcopy(io_pairs) + processed_prior_cov = deepcopy(prior_cov_in) processed_obs_noise_cov = deepcopy(obs_noise_cov_in) # apply_to is the string "in", "out" etc. for idx in reverse(eachindex(encoder_schedule)) (processor, extract_data, apply_to) = encoder_schedule[idx] + if apply_to == "in" + structure_matrix = processed_prior_cov + elseif apply_to == "out" + structure_matrix = processed_obs_noise_cov + end + @info "Decoding data: $(apply_to) with $(processor)" - processed = decode_data(processor, extract_data(processed_io_pairs), processed_obs_noise_cov, apply_to) + processed = decode_data(processor, extract_data(processed_io_pairs), structure_matrix, apply_to) if apply_to == "in" + processed_prior_cov = decode_structure_matrix(processor, structure_matrix) processed_io_pairs = PairedDataContainer(processed, get_outputs(processed_io_pairs)) elseif apply_to == "out" - processed_obs_noise_cov = decode_obs_noise_cov(processor, processed_obs_noise_cov) + processed_obs_noise_cov = decode_structure_matrix(processor, structure_matrix) processed_io_pairs = PairedDataContainer(get_inputs(processed_io_pairs), processed) end end - return processed_io_pairs, processed_obs_noise_cov + return processed_io_pairs, processed_prior_cov, processed_obs_noise_cov end + end # module diff --git a/test/Utilities/runtests.jl b/test/Utilities/runtests.jl index 7caeab393..ae5689398 100644 --- a/test/Utilities/runtests.jl +++ b/test/Utilities/runtests.jl @@ -79,9 +79,11 @@ end out_dim = 50 samples = 100 - in_data = 20 * randn(in_dim, samples) + x = randn(in_dim, in_dim) + prior_cov = x * x' + in_data = rand(MvNormal(zeros(in_dim), prior_cov), samples) obs_noise_cov = [max(5.0 - abs(i - j), 0.0) for i in 1:out_dim, j in 1:out_dim] # [5 4 3 2 1 0 0 ...] off diagonal - out_data = rand(MvNormal(-10 * ones(out_dim), obs_noise_cov), samples) .- 10 + out_data = rand(MvNormal(-10 * ones(out_dim), obs_noise_cov), samples) io_pairs = PairedDataContainer(in_data, out_data) @@ -89,23 +91,25 @@ end schedule_builder = [ (zscore_scale(), "in_and_out"), # (quartile_scale(), "in"), - (standardize(), "out"), - (minmax_scale(), "in_and_out"), + (standardize(), "in_and_out"), + (minmax_scale(), "out"), ] # make schedule more parsable encoder_schedule = create_encoder_schedule(schedule_builder) # encode the data using the schedule - encoded_io_pairs, encoded_obs_noise_cov = encode_with_schedule(encoder_schedule, io_pairs, obs_noise_cov) + (encoded_io_pairs, encoded_prior_cov, encoded_obs_noise_cov) = + encode_with_schedule(encoder_schedule, io_pairs, prior_cov, obs_noise_cov) # decode the data using the schedule - decoded_io_pairs, decoded_obs_noise_cov = - decode_with_schedule(encoder_schedule, encoded_io_pairs, encoded_obs_noise_cov) + (decoded_io_pairs, decoded_prior_cov, decoded_obs_noise_cov) = + decode_with_schedule(encoder_schedule, encoded_io_pairs, encoded_prior_cov, encoded_obs_noise_cov) tol = 1e-12 @test all(isapprox.(get_inputs(io_pairs), get_inputs(decoded_io_pairs), atol = tol)) @test all(isapprox.(get_outputs(io_pairs), get_outputs(decoded_io_pairs), atol = tol)) + @test all(isapprox.(prior_cov, decoded_prior_cov, atol = tol)) @test all(isapprox.(obs_noise_cov, decoded_obs_noise_cov, atol = tol)) end From 6160ec4f7fb9f207f56410dae6fdffd684eb27b4 Mon Sep 17 00:00:00 2001 From: odunbar Date: Thu, 19 Jun 2025 14:33:08 -0700 Subject: [PATCH 05/25] added Decorrelation, and a bunch of unit tests --- src/Utilities.jl | 156 +++++++++++++++++++++++++++++-------- test/Utilities/runtests.jl | 127 +++++++++++++++++++++++++++++- 2 files changed, 251 insertions(+), 32 deletions(-) diff --git a/src/Utilities.jl b/src/Utilities.jl index 8b8189a22..c39364b3f 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -18,7 +18,7 @@ export zscore2orig export PairedDataContainerProcessor, DataContainerProcessor export UnivariateAffineScaling, QuartileScaling, MinMaxScaling, ZScoreScaling, Standardizer -export quartile_scale, minmax_scale, zscore_scale, standardize +export quartile_scale, minmax_scale, zscore_scale, standardize, decorrelate export get_type, get_shift, get_scale, get_data_mean, get_data_reduction_mat, get_data_inflation_mat, get_rank export create_encoder_schedule, encode_with_schedule, decode_with_schedule export initialize_processor!, @@ -167,7 +167,7 @@ function initialize_processor!( T::Type{MMS}, ) where {MM <: AbstractMatrix, MMS <: MinMaxScaling} minmax_vec = [[minimum(dd), maximum(dd)] for dd in eachrow(data)] - minmax_mat = reduce(hcat, minmax_vec) # 3 rows: Q1, Q2, and Q3 + minmax_mat = reduce(hcat, minmax_vec) # 2 rows: min max append!(get_shift(as), minmax_mat[1, :]) append!(get_scale(as), (minmax_mat[2, :] - minmax_mat[1, :])) end @@ -178,7 +178,7 @@ function initialize_processor!( T::Type{ZSS}, ) where {MM <: AbstractMatrix, ZSS <: ZScoreScaling} stat_vec = [[mean(dd), std(dd)] for dd in eachrow(data)] - stat_mat = reduce(hcat, stat_vec) # 3 rows: Q1, Q2, and Q3 + stat_mat = reduce(hcat, stat_vec) # 2 rows: mean, std append!(get_shift(as), stat_mat[1, :]) append!(get_scale(as), stat_mat[2, :]) end @@ -226,7 +226,9 @@ end """ $(TYPEDEF) -Standardizes the data to a multivariate N(0,I) distribution via `(xcov)^{-1/2}*(x-x_mean)`, along with rank reduction if data is low rank, or if the user provides a rank +Standardizes the data to a multivariate N(0,I) distribution via `(xcov)^{-1/2}*(x-x_mean)`, along with rank reduction if data is low rank, or if the user provides a rank. This guarantees that the data samples will have sample mean `0` and covariance `I` after processing. + +Similar to the `Decorrelater`, but there the structure matrices will directly be used to perform decorrelation, and will become `I` after processing. In that case, the data samples may not have covariance `I` after processing. """ struct Standardizer{VV1 <: AbstractVector, VV2 <: AbstractVector, VV3 <: AbstractVector} <: DataContainerProcessor "user-provided rank for additional truncation of the space. used if rank < rank(data_cov)" @@ -234,9 +236,9 @@ struct Standardizer{VV1 <: AbstractVector, VV2 <: AbstractVector, VV3 <: Abstrac "storage for the data mean" data_mean::VV1 "storage for wide-matrix representing `data_cov^{-1/2}`" - data_reduction_mat::VV2 + encoder_mat::VV2 "storage for tall-matrix representing `data_cov^{1/2}`" - data_inflation_mat::VV3 + decoder_mat::VV3 end standardize() = Standardizer(typemax(Int), Float64[], Any[], Any[]) @@ -247,13 +249,13 @@ standardize(rk::Int, T::Type) = Standardizer(rk, T[], Any[], Any[]) standardize(rk::Int, T1::Type, T2::Type, T3::Type) = Standardizer(rk, T1[], T2[], T3[]) get_data_mean(ss::Standardizer) = ss.data_mean -get_data_reduction_mat(ss::Standardizer) = ss.data_reduction_mat -get_data_inflation_mat(ss::Standardizer) = ss.data_inflation_mat +get_encoder_mat(ss::Standardizer) = ss.encoder_mat +get_decoder_mat(ss::Standardizer) = ss.decoder_mat get_rank(ss::Standardizer) = ss.rank function Base.show(io::IO, ss::Standardizer) if get_rank(ss) < typemax(Int) - out = "Standardizer, user-rank=$(get_rank(ss))" + out = "Standardizer: user-rank=$(get_rank(ss))" else out = "Standardizer" end @@ -264,7 +266,7 @@ function initialize_processor!(ss::Standardizer, data::MM) where {MM <: Abstract if length(get_data_mean(ss)) == 0 append!(get_data_mean(ss), mean(data, dims = 2)) end - if length(get_data_reduction_mat(ss)) == 0 + if length(get_encoder_mat(ss)) == 0 # can use tsvd here so we only compute up to rk data_cov = cov(data, dims = 2) svdc = svd(data_cov) @@ -276,23 +278,23 @@ function initialize_processor!(ss::Standardizer, data::MM) where {MM <: Abstract mat = svdc.U' end - reduction_mat = Diagonal(1 ./ sqrt.(svdc.S[1:rk])) * mat[1:rk, :] - inflation_mat = mat[1:rk, :]' * Diagonal(sqrt.(svdc.S[1:rk])) - push!(get_data_reduction_mat(ss), reduction_mat) - push!(get_data_inflation_mat(ss), inflation_mat) + encoder_mat = Diagonal(1 ./ sqrt.(svdc.S[1:rk])) * mat[1:rk, :] + decoder_mat = mat[1:rk, :]' * Diagonal(sqrt.(svdc.S[1:rk])) + push!(get_encoder_mat(ss), encoder_mat) + push!(get_decoder_mat(ss), decoder_mat) end end function encode_data(ss::Standardizer, data::MM) where {MM <: AbstractMatrix} data_mean = get_data_mean(ss) - reduction_mat = get_data_reduction_mat(ss)[1] - return reduction_mat * (data .- data_mean) + encoder_mat = get_encoder_mat(ss)[1] + return encoder_mat * (data .- data_mean) end function decode_data(ss::Standardizer, data::MM) where {MM <: AbstractMatrix} data_mean = get_data_mean(ss) - inflation_mat = get_data_inflation_mat(ss)[1] - return inflation_mat * data .+ data_mean + decoder_mat = get_decoder_mat(ss)[1] + return decoder_mat * data .+ data_mean end initialize_processor!(ss::Standardizer, data::MM, structure_matrix) where {MM <: AbstractMatrix} = @@ -301,34 +303,126 @@ encode_data(ss::Standardizer, data::MM, structure_matrix) where {MM <: AbstractM decode_data(ss::Standardizer, data::MM, structure_matrix) where {MM <: AbstractMatrix} = decode_data(ss, data) function encode_structure_matrix(ss::Standardizer, structure_matrix::MM) where {MM <: AbstractMatrix} - reduction_mat = get_data_reduction_mat(ss)[1] - return reduction_mat * structure_matrix * reduction_mat' + encoder_mat = get_encoder_mat(ss)[1] + return encoder_mat * structure_matrix * encoder_mat' end function decode_structure_matrix(ss::Standardizer, enc_structure_matrix::MM) where {MM <: AbstractMatrix} - inflation_mat = get_data_inflation_mat(ss)[1] - return inflation_mat * enc_structure_matrix * inflation_mat' + decoder_mat = get_decoder_mat(ss)[1] + return decoder_mat * enc_structure_matrix * decoder_mat' end """ $(TYPEDEF) -Decorrelate the data via an SVD decomposition using a structure_matrix (`prior_cov` for inputs, `obs_noise_cov` for outputs), with optional truncation of singular vectors corresponding to largest singular values. +Decorrelate the data via an SVD decomposition using a structure_matrix (`prior_cov` for inputs, `obs_noise_cov` for outputs), with optional truncation of singular vectors corresponding to largest singular values. The structure matrix will also become exactly `I` after processing. + +Similar to the `Standardizer`, except that the `Standardizer` uses the estimated covariance of the data to decorrelate in-place of the structure matrix. There, the data samples will have sample mean `0` and covariance `I` after processing. + +In the `Decorrelater` the user can do a combined approach where one uses cov(data) + structure matrix for decorrelation with the `add_estimated_data_cov=true` # Fields $(TYPEDFIELDS) """ -struct Decorrelater{MM1, MM2, FT} <: DataContainerProcessor - "the structure matrix - provided by the user" - structure_matrix::MM1 - "the inverse of the structure matrix" - inv_structure_matrix::MM2 +struct Decorrelater{VV1, VV2, VV3, FT} <: DataContainerProcessor + "storage for the data mean" + data_mean::VV1 + "the (possibly-encoded) structure matrix - provided by the user" + encoder_mat::VV2 + "the inverse of the (possibly-encoded) structure matrix" + decoder_mat::VV3 "the fraction of variance to be retained after truncating singular values (1 implies no truncation)" - retained_variance::FT + retain_var::FT + "If true, add the estimated cov(data) to the structure matrix for to modify the decorrelation" + add_estimated_cov::Bool +end + +decorrelate(; retain_var::FT = Float64(1.0), add_estimated_cov = false) where {FT} = + Decorrelater([], [], [], min(max(retain_var, FT(0)), FT(1)), add_estimated_cov) + +get_data_mean(dd::Decorrelater) = dd.data_mean +get_encoder_mat(dd::Decorrelater) = dd.encoder_mat +get_decoder_mat(dd::Decorrelater) = dd.decoder_mat +get_retain_var(dd::Decorrelater) = dd.retain_var +get_add_estimated_cov(dd::Decorrelater) = dd.add_estimated_cov + +function Base.show(io::IO, dd::Decorrelater) + out = "Decorrelater" + if get_retain_var(dd) < 1.0 + out *= ": retain_var=$(get_retain_var(dd)) " + end + if get_add_estimated_cov(dd) + out *= ": add_estimated_cov=$(get_add_estimated_cov(dd)) " + end + print(io, out) end +function initialize_processor!( + dd::Decorrelater, + data::MM, + structure_matrix::USorM, +) where {MM <: AbstractMatrix, USorM <: Union{UniformScaling, AbstractMatrix}} + if length(get_data_mean(dd)) == 0 + append!(get_data_mean(dd), mean(data, dims = 2)) + end + + if length(get_encoder_mat(dd)) == 0 + # Can do tsvd here for large matrices + if get_add_estimated_cov(dd) + svdA = svd(structure_matrix + cov(data .- mean(data, dims = 2), dims = 2)) # with data covariance? + else + svdA = svd(structure_matrix) + end + ret_var = get_retain_var(dd) + if ret_var < 1.0 + sv_cumsum = cumsum(svdA.S .^ 2) / sum(svdA.S .^ 2) # variance contributions are (sing_val)^2 + trunc = minimum(findall(x -> (x > ret_var), sv_cumsum)) + @info " truncating at $(trunc)/$(length(sv_cumsum)) retaining $(100.0*sv_cumsum[trunc])% of the variance of the structure matrix" + else + trunc = rank(structure_matrix) + end + sqrt_inv_sv = Diagonal(1.0 ./ sqrt.(svdA.S[1:trunc])) + sqrt_sv = Diagonal(sqrt.(svdA.S[1:trunc])) + + if size(svdA.U, 1) == size(svdA.U, 2) + mat = svdA.Vt + else + mat = svdA.U + end + + encoder_mat = sqrt_inv_sv * mat[1:trunc, :] + decoder_mat = mat[1:trunc, :]' * sqrt_sv + + push!(get_encoder_mat(dd), encoder_mat) + push!(get_decoder_mat(dd), decoder_mat) + end +end + +function encode_data(dd::Decorrelater, data::MM) where {MM <: AbstractMatrix} + data_mean = get_data_mean(dd) + encoder_mat = get_encoder_mat(dd)[1] + return encoder_mat * (data .- data_mean) +end + +function decode_data(dd::Decorrelater, data::MM) where {MM <: AbstractMatrix} + data_mean = get_data_mean(dd) + decoder_mat = get_decoder_mat(dd)[1] + return decoder_mat * data .+ data_mean +end +encode_data(dd::Decorrelater, data::MM, structure_matrix) where {MM <: AbstractMatrix} = encode_data(dd, data) +decode_data(dd::Decorrelater, data::MM, structure_matrix) where {MM <: AbstractMatrix} = decode_data(dd, data) + +function encode_structure_matrix(dd::Decorrelater, structure_matrix::MM) where {MM <: AbstractMatrix} + encoder_mat = get_encoder_mat(dd)[1] + return encoder_mat * structure_matrix * encoder_mat' +end + +function decode_structure_matrix(dd::Decorrelater, enc_structure_matrix::MM) where {MM <: AbstractMatrix} + decoder_mat = get_decoder_mat(dd)[1] + return decoder_mat * enc_structure_matrix * decoder_mat' +end # ...struct VariationalAutoEncoder <: DataContainerProcessor end @@ -487,7 +581,7 @@ function encode_with_schedule( # apply_to is the string "in", "out" etc. for (processor, extract_data, apply_to) in encoder_schedule - @info "Encoding data: $(apply_to) with $(processor)" + @info "Encoding data: \"$(apply_to)\" with $(processor)" if apply_to == "in" structure_matrix = processed_prior_cov elseif apply_to == "out" @@ -531,7 +625,7 @@ function decode_with_schedule( structure_matrix = processed_obs_noise_cov end - @info "Decoding data: $(apply_to) with $(processor)" + @info "Decoding data: \"$(apply_to)\" with $(processor)" processed = decode_data(processor, extract_data(processed_io_pairs), structure_matrix, apply_to) if apply_to == "in" diff --git a/test/Utilities/runtests.jl b/test/Utilities/runtests.jl index ae5689398..be0a47cef 100644 --- a/test/Utilities/runtests.jl +++ b/test/Utilities/runtests.jl @@ -74,6 +74,10 @@ end @testset "Data Preprocessing" begin + using LinearAlgebra, Distributions, Statistics + using CalibrateEmulateSample.DataContainers + using CalibrateEmulateSample.Utilities + # get some data as IO pairs in_dim = 10 out_dim = 50 @@ -86,13 +90,134 @@ end out_data = rand(MvNormal(-10 * ones(out_dim), obs_noise_cov), samples) io_pairs = PairedDataContainer(in_data, out_data) + test_names = ["zscore", "quartile", "minmax"] + + # Test the encodings + univariate_tests = + [[(zscore_scale(), "in_and_out")], [(quartile_scale(), "in_and_out")], [(minmax_scale(), "in_and_out")]] + + tol = 1e-12 + + for (name, sch) in zip(test_names, univariate_tests) + encoder_schedule = create_encoder_schedule(sch) + (encoded_io_pairs, encoded_prior_cov, encoded_obs_noise_cov) = + encode_with_schedule(encoder_schedule, io_pairs, prior_cov, obs_noise_cov) + + (decoded_io_pairs, decoded_prior_cov, decoded_obs_noise_cov) = + decode_with_schedule(encoder_schedule, encoded_io_pairs, encoded_prior_cov, encoded_obs_noise_cov) + for (enc_dat, dec_dat, test_dat, enc_covv, dec_covv, test_covv, dim) in zip( + (get_inputs(encoded_io_pairs), get_outputs(encoded_io_pairs)), + (get_inputs(decoded_io_pairs), get_outputs(decoded_io_pairs)), + (get_inputs(io_pairs), get_outputs(io_pairs)), + (encoded_prior_cov, encoded_obs_noise_cov), + (decoded_prior_cov, decoded_obs_noise_cov), + (prior_cov, obs_noise_cov), + (in_dim, out_dim), + ) + + if name == "zscore" + stat_vec = [[mean(dd), std(dd)] for dd in eachrow(enc_dat)] + stat_mat = reduce(hcat, stat_vec) + @test all(isapprox.(stat_mat[1, :], zeros(dim), atol = tol)) + @test all(isapprox.(stat_mat[2, :], ones(dim), atol = tol)) + + test_vec = [[mean(dd), std(dd)] for dd in eachrow(test_dat)] + test_mat = reduce(hcat, test_vec) + @test isapprox(norm(enc_covv - Diagonal(1 ./ test_mat[2, :] .^ 2) * test_covv), 0.0, atol = tol * dim^2) + elseif name == "quartile" + quartiles_vec = [quantile(dd, [0.25, 0.5, 0.75]) for dd in eachrow(enc_dat)] + quartiles_mat = reduce(hcat, quartiles_vec) # 3 rows: Q1, Q2, and Q3 + @test all(isapprox.(quartiles_mat[2, :], zeros(dim), atol = tol)) + @test all(isapprox.(quartiles_mat[3, :] - quartiles_mat[1, :], ones(dim), atol = tol)) + test_vec = [quantile(dd, [0.25, 0.5, 0.75]) for dd in eachrow(test_dat)] + test_mat = reduce(hcat, test_vec) + @test isapprox( + norm(enc_covv - Diagonal(1 ./ (test_mat[3, :] - test_mat[1, :]) .^ 2) * test_covv), + 0.0, + atol = tol * dim^2, + ) + elseif name == "minmax" + minmax_vec = [[minimum(dd), maximum(dd)] for dd in eachrow(enc_dat)] + minmax_mat = reduce(hcat, minmax_vec) # 2 rows: min, max + @test all(isapprox.(minmax_mat[1, :], zeros(dim), atol = tol)) + @test all(isapprox.(minmax_mat[2, :], ones(dim), atol = tol)) + test_vec = [[minimum(dd), maximum(dd)] for dd in eachrow(test_dat)] + test_mat = reduce(hcat, test_vec) + @test isapprox( + norm(enc_covv - Diagonal(1 ./ (test_mat[2, :] - test_mat[1, :]) .^ 2) * test_covv), + 0.0, + atol = tol * dim^2, + ) + end + + @test isapprox(norm(dec_dat - test_dat), 0.0, atol = tol * dim) + @test isapprox(norm(dec_covv - test_covv), 0.0, atol = tol * dim^2) + + end + + end + + + # multivariate lossless encodings + test_names = ["standardize", "decorrelate"] + multivariate_lossless_tests = [[(standardize(), "in_and_out")], [(decorrelate(), "in_and_out")]] + + for (name, sch) in zip(test_names, multivariate_lossless_tests) + encoder_schedule = create_encoder_schedule(sch) + (encoded_io_pairs, encoded_prior_cov, encoded_obs_noise_cov) = + encode_with_schedule(encoder_schedule, io_pairs, prior_cov, obs_noise_cov) + + (decoded_io_pairs, decoded_prior_cov, decoded_obs_noise_cov) = + decode_with_schedule(encoder_schedule, encoded_io_pairs, encoded_prior_cov, encoded_obs_noise_cov) + for (enc_dat, dec_dat, test_dat, enc_covv, dec_covv, test_covv, dim) in zip( + (get_inputs(encoded_io_pairs), get_outputs(encoded_io_pairs)), + (get_inputs(decoded_io_pairs), get_outputs(decoded_io_pairs)), + (get_inputs(io_pairs), get_outputs(io_pairs)), + (encoded_prior_cov, encoded_obs_noise_cov), + (decoded_prior_cov, decoded_obs_noise_cov), + (prior_cov, obs_noise_cov), + (in_dim, out_dim), + ) + + if name == "standardize" + pop_mean = mean(enc_dat, dims = 2) + pop_cov = cov(enc_dat, dims = 2) + dimm = size(pop_cov, 1) + @test all(isapprox.(pop_mean, zeros(dimm), atol = tol)) + @test isapprox(norm(pop_cov - I), 0.0, atol = tol * dimm^2) # expect very accurate + @test isapprox(norm(enc_covv - I), 0.0, atol = 0.2 * dimm^2) # expect poorly accurate + + elseif name == "decorrelate" + pop_mean = mean(enc_dat, dims = 2) + pop_cov = cov(enc_dat, dims = 2) + dimm = size(pop_cov, 1) + @test all(isapprox.(pop_mean, zeros(dimm), atol = tol)) + @test isapprox(norm(pop_cov - I), 0.0, atol = 0.2 * dimm^2) # expect poorly accurate + @test isapprox(norm(enc_covv - I), 0.0, atol = tol * dimm^2) # expect very accurate + end + + @test isapprox(norm(dec_dat - test_dat), 0.0, atol = tol * dim) + @test isapprox(norm(dec_covv - test_covv), 0.0, atol = tol * dim^2) + + end + end + + multivariate_lossy_tests = [ + [(standardize(5), "in_and_out")], + [(decorrelate(retain_var = 0.95), "in_and_out")], + [(decorrelate(add_estimated_cov = true), "in_and_out")], + ] + + # LOSSY TESTS HERE + - # Create a lossless encoding schedule + # combine a few lossless encoding schedules schedule_builder = [ (zscore_scale(), "in_and_out"), # (quartile_scale(), "in"), (standardize(), "in_and_out"), (minmax_scale(), "out"), + (decorrelate(), "in_and_out"), ] # make schedule more parsable From 9e0aa25e4f727bdfc2bee0b3a225e1c57606fed6 Mon Sep 17 00:00:00 2001 From: odunbar Date: Thu, 19 Jun 2025 15:52:18 -0700 Subject: [PATCH 06/25] completed unit tests for encoder-decoders --- test/Utilities/runtests.jl | 100 +++++++++++++++++++++++++++++++++---- 1 file changed, 91 insertions(+), 9 deletions(-) diff --git a/test/Utilities/runtests.jl b/test/Utilities/runtests.jl index be0a47cef..d33a83e49 100644 --- a/test/Utilities/runtests.jl +++ b/test/Utilities/runtests.jl @@ -90,15 +90,33 @@ end out_data = rand(MvNormal(-10 * ones(out_dim), obs_noise_cov), samples) io_pairs = PairedDataContainer(in_data, out_data) - test_names = ["zscore", "quartile", "minmax"] + test_names = [ + "zscore", + "quartile", + "minmax", + "standardize", + "decorrelate", + "decorrelate-estimate-cov", + "standardize-truncate-to-5", + "decorrelate-retain-0.95-var", + ] - # Test the encodings - univariate_tests = - [[(zscore_scale(), "in_and_out")], [(quartile_scale(), "in_and_out")], [(minmax_scale(), "in_and_out")]] + # Test encodings-decodings individually + univariate_tests = [ + [(zscore_scale(), "in_and_out")], + [(quartile_scale(), "in_and_out")], + [(minmax_scale(), "in_and_out")], + [(standardize(), "in_and_out")], + [(decorrelate(), "in_and_out")], + [(decorrelate(add_estimated_cov = true), "in_and_out")], + [(standardize(5), "in_and_out")], + [(decorrelate(retain_var = 0.95), "in_and_out")], + ] + lossless = [fill(true, 6); fill(false, 2)] # are these lossy approximations tol = 1e-12 - for (name, sch) in zip(test_names, univariate_tests) + for (name, sch, ll_flag) in zip(test_names, univariate_tests, lossless) encoder_schedule = create_encoder_schedule(sch) (encoded_io_pairs, encoded_prior_cov, encoded_obs_noise_cov) = encode_with_schedule(encoder_schedule, io_pairs, prior_cov, obs_noise_cov) @@ -114,7 +132,7 @@ end (prior_cov, obs_noise_cov), (in_dim, out_dim), ) - + # univariate "rescaling" tests if name == "zscore" stat_vec = [[mean(dd), std(dd)] for dd in eachrow(enc_dat)] stat_mat = reduce(hcat, stat_vec) @@ -150,8 +168,47 @@ end ) end - @test isapprox(norm(dec_dat - test_dat), 0.0, atol = tol * dim) - @test isapprox(norm(dec_covv - test_covv), 0.0, atol = tol * dim^2) + # Multivariate lossless tests + pop_mean = mean(enc_dat, dims = 2) + pop_cov = cov(enc_dat, dims = 2) + dimm = size(pop_cov, 1) + if name == "standardize" + @test all(isapprox.(pop_mean, zeros(dimm), atol = tol)) + @test isapprox(norm(pop_cov - I), 0.0, atol = tol * dimm^2) # expect very accurate + @test isapprox(norm(enc_covv - I), 0.0, atol = 0.2 * dimm^2) # expect poorly accurate + + elseif name == "decorrelate" + @test all(isapprox.(pop_mean, zeros(dimm), atol = tol)) + @test isapprox(norm(pop_cov - I), 0.0, atol = 0.2 * dimm^2) # expect poorly accurate + @test isapprox(norm(enc_covv - I), 0.0, atol = tol * dimm^2) # expect very accurate + elseif name == "decorrelate-estimate-cov" + @test all(isapprox.(pop_mean, zeros(dimm), atol = tol)) + @test isapprox(norm(pop_cov - I), 0.0, atol = 0.2 * dimm^2) # expect poorly accurate + @test isapprox(norm(enc_covv - I), 0.0, atol = 0.2 * dimm^2) # expect poorly accurate + end + + # Multivariate lossy dim-reduction tests + if name == "standardize-truncate-to-5" + @test dimm == 5 + @test all(isapprox.(pop_mean, zeros(dimm), atol = tol)) + @test isapprox(norm(pop_cov - I), 0.0, atol = tol * dimm^2) # expect very accurate + @test isapprox(norm(enc_covv - I), 0.0, atol = 0.2 * dimm^2) # expect poorly accurate + + elseif name == "decorrelate-retain-0.95-var" + svdc = svd(test_covv) + var_cumsum = cumsum(svdc.S .^ 2) ./ sum(svdc.S .^ 2) + @test var_cumsum[dimm] > 0.95 + @test var_cumsum[dimm - 1] < 0.95 + @test all(isapprox.(pop_mean, zeros(dimm), atol = tol)) + @test isapprox(norm(pop_cov - I), 0.0, atol = 0.2 * dimm^2) # expect poorly accurate + @test isapprox(norm(enc_covv - I), 0.0, atol = tol * dimm^2) # expect very accurate + end + + # test decode approximation of lossless options + if ll_flag + @test isapprox(norm(dec_dat - test_dat), 0.0, atol = tol * dim) + @test isapprox(norm(dec_covv - test_covv), 0.0, atol = tol * dim^2) + end end @@ -160,7 +217,7 @@ end # multivariate lossless encodings test_names = ["standardize", "decorrelate"] - multivariate_lossless_tests = [[(standardize(), "in_and_out")], [(decorrelate(), "in_and_out")]] + multivariate_lossless_tests = [] for (name, sch) in zip(test_names, multivariate_lossless_tests) encoder_schedule = create_encoder_schedule(sch) @@ -202,12 +259,37 @@ end end end + test_names = ["standardize-truncate-to-5", "decorrelate-retain-partial-var", "decorrelate-estimate-cov"] + multivariate_lossy_tests = [ [(standardize(5), "in_and_out")], [(decorrelate(retain_var = 0.95), "in_and_out")], [(decorrelate(add_estimated_cov = true), "in_and_out")], ] + for (name, sch) in zip(test_names, multivariate_lossless_tests) + encoder_schedule = create_encoder_schedule(sch) + (encoded_io_pairs, encoded_prior_cov, encoded_obs_noise_cov) = + encode_with_schedule(encoder_schedule, io_pairs, prior_cov, obs_noise_cov) + + (decoded_io_pairs, decoded_prior_cov, decoded_obs_noise_cov) = + decode_with_schedule(encoder_schedule, encoded_io_pairs, encoded_prior_cov, encoded_obs_noise_cov) + for (enc_dat, dec_dat, test_dat, enc_covv, dec_covv, test_covv, dim) in zip( + (get_inputs(encoded_io_pairs), get_outputs(encoded_io_pairs)), + (get_inputs(decoded_io_pairs), get_outputs(decoded_io_pairs)), + (get_inputs(io_pairs), get_outputs(io_pairs)), + (encoded_prior_cov, encoded_obs_noise_cov), + (decoded_prior_cov, decoded_obs_noise_cov), + (prior_cov, obs_noise_cov), + (in_dim, out_dim), + ) + + + end + end + + + # LOSSY TESTS HERE From 2441e943c8c47335c43e34c83a8bab74e55dbee2 Mon Sep 17 00:00:00 2001 From: odunbar Date: Thu, 19 Jun 2025 15:55:08 -0700 Subject: [PATCH 07/25] add getters --- src/Utilities.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Utilities.jl b/src/Utilities.jl index c39364b3f..09d3d9d31 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -19,7 +19,7 @@ export PairedDataContainerProcessor, DataContainerProcessor export UnivariateAffineScaling, QuartileScaling, MinMaxScaling, ZScoreScaling, Standardizer export quartile_scale, minmax_scale, zscore_scale, standardize, decorrelate -export get_type, get_shift, get_scale, get_data_mean, get_data_reduction_mat, get_data_inflation_mat, get_rank +export get_type, get_shift, get_scale, get_data_mean, get_data_encoder_mat, get_data_decoder_mat, get_rank, get_retain_var, get_add_estimated_cov export create_encoder_schedule, encode_with_schedule, decode_with_schedule export initialize_processor!, initialize_and_encode_data!, encode_data, decode_data, encode_obs_noise_cov, decode_obs_noise_cov @@ -371,7 +371,7 @@ function initialize_processor!( # Can do tsvd here for large matrices if get_add_estimated_cov(dd) - svdA = svd(structure_matrix + cov(data .- mean(data, dims = 2), dims = 2)) # with data covariance? + svdA = svd(structure_matrix + cov(data, dims = 2)) # with data covariance? else svdA = svd(structure_matrix) end From dda82ffada3633836c5d77aeb6f29bfb5195b6d1 Mon Sep 17 00:00:00 2001 From: odunbar Date: Thu, 19 Jun 2025 16:25:31 -0700 Subject: [PATCH 08/25] full unit tests --- src/Utilities.jl | 20 ++++++++--------- test/Utilities/runtests.jl | 44 +++++++++++++++++++++++++++++++++----- 2 files changed, 49 insertions(+), 15 deletions(-) diff --git a/src/Utilities.jl b/src/Utilities.jl index 09d3d9d31..e5beda369 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -16,10 +16,18 @@ export orig2zscore export zscore2orig export PairedDataContainerProcessor, DataContainerProcessor -export UnivariateAffineScaling, QuartileScaling, MinMaxScaling, ZScoreScaling, Standardizer +export UnivariateAffineScaling, AffineScaler, QuartileScaling, MinMaxScaling, ZScoreScaling, Standardizer, Decorrelater export quartile_scale, minmax_scale, zscore_scale, standardize, decorrelate -export get_type, get_shift, get_scale, get_data_mean, get_data_encoder_mat, get_data_decoder_mat, get_rank, get_retain_var, get_add_estimated_cov +export get_type, + get_shift, + get_scale, + get_data_mean, + get_encoder_mat, + get_decoder_mat, + get_rank, + get_retain_var, + get_add_estimated_cov export create_encoder_schedule, encode_with_schedule, decode_with_schedule export initialize_processor!, initialize_and_encode_data!, encode_data, decode_data, encode_obs_noise_cov, decode_obs_noise_cov @@ -131,15 +139,11 @@ struct AffineScaler{T, VV <: AbstractVector} <: DataContainerProcessor end quartile_scale() = AffineScaler(QuartileScaling) -quartile_scale(T::Type) = AffineScaler(QuartileScaling, T) minmax_scale() = AffineScaler(MinMaxScaling) -minmax_scale(T::Type) = AffineScaler(MinMaxScaling, T) zscore_scale() = AffineScaler(ZScoreScaling) -zscore_scale(T::Type) = AffineScaler(ZScoreScaling, T) AffineScaler(::Type{UAS}) where {UAS <: UnivariateAffineScaling} = AffineScaler{UAS, Vector{Float64}}(Float64[], Float64[]) -AffineScaler(::Type{UAS}, T::Type) where {UAS <: UnivariateAffineScaling} = AffineScaler{UAS, Vector{T}}(T[], T[]) get_type(as::AffineScaler{T, VV}) where {T, VV} = T get_shift(as::AffineScaler) = as.shift @@ -242,11 +246,7 @@ struct Standardizer{VV1 <: AbstractVector, VV2 <: AbstractVector, VV3 <: Abstrac end standardize() = Standardizer(typemax(Int), Float64[], Any[], Any[]) -standardize(T::Type) = Standardizer(typemax(Int), T[], Any[], Any[]) -standardize(T1::Type, T2::Type, T3::Type) = Standardizer(typemax(Int), T1[], T2[], T3[]) standardize(rk::Int) = Standardizer(rk, Float64[], Any[], Any[]) -standardize(rk::Int, T::Type) = Standardizer(rk, T[], Any[], Any[]) -standardize(rk::Int, T1::Type, T2::Type, T3::Type) = Standardizer(rk, T1[], T2[], T3[]) get_data_mean(ss::Standardizer) = ss.data_mean get_encoder_mat(ss::Standardizer) = ss.encoder_mat diff --git a/test/Utilities/runtests.jl b/test/Utilities/runtests.jl index d33a83e49..3ad5aa8ca 100644 --- a/test/Utilities/runtests.jl +++ b/test/Utilities/runtests.jl @@ -74,11 +74,42 @@ end @testset "Data Preprocessing" begin - using LinearAlgebra, Distributions, Statistics - using CalibrateEmulateSample.DataContainers - using CalibrateEmulateSample.Utilities - - # get some data as IO pairs + # quick build tests and test getters + zs = zscore_scale() + mm = minmax_scale() + qq = quartile_scale() + QQ = AffineScaler{QuartileScaling, Vector{Int}}([1], [2]) + @test isa(zs, AffineScaler) + @test get_type(zs) == ZScoreScaling + @test isa(mm, AffineScaler) + @test get_type(mm) == MinMaxScaling + @test isa(qq, AffineScaler) + @test get_type(qq) == QuartileScaling + @test get_shift(QQ) == [1] + @test get_scale(QQ) == [2] + + ss = standardize() + @test isa(ss, Standardizer) + @test get_rank(ss) == typemax(Int) + ss2 = standardize(10) + @test get_rank(ss2) == 10 + SS = Standardizer(1, [1], [2], [3]) + @test get_data_mean(SS) == [1] + @test get_encoder_mat(SS) == [2] + @test get_decoder_mat(SS) == [3] + + dd = decorrelate() + @test get_retain_var(dd) == 1.0 + @test get_add_estimated_cov(dd) == false + dd2 = decorrelate(retain_var = 0.7, add_estimated_cov = true) + @test get_retain_var(dd2) == 0.7 + @test get_add_estimated_cov(dd2) == true + DD = Decorrelater([1], [2], [3], 1.0, false) + @test get_data_mean(DD) == [1] + @test get_encoder_mat(DD) == [2] + @test get_decoder_mat(DD) == [3] + + # get some data as IO pairs for functional tests in_dim = 10 out_dim = 50 samples = 100 @@ -101,6 +132,8 @@ end "decorrelate-retain-0.95-var", ] + + # Test encodings-decodings individually univariate_tests = [ [(zscore_scale(), "in_and_out")], @@ -114,6 +147,7 @@ end ] lossless = [fill(true, 6); fill(false, 2)] # are these lossy approximations + # functional test pipeline tol = 1e-12 for (name, sch, ll_flag) in zip(test_names, univariate_tests, lossless) From 322106d2201d5b7b0d4e744aab9331c3e14959c0 Mon Sep 17 00:00:00 2001 From: odunbar Date: Thu, 19 Jun 2025 16:57:12 -0700 Subject: [PATCH 09/25] new encode/decode methods once initialized --- src/Utilities.jl | 144 +++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 119 insertions(+), 25 deletions(-) diff --git a/src/Utilities.jl b/src/Utilities.jl index e5beda369..bf110b6a8 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -214,8 +214,6 @@ end initialize_processor!(as::AffineScaler, data::MM, structure_matrix) where {MM <: AbstractMatrix} = initialize_processor!(as, data) -encode_data(as::AffineScaler, data::MM, structure_matrix) where {MM <: AbstractMatrix} = encode_data(as, data) -decode_data(as::AffineScaler, data::MM, structure_matrix) where {MM <: AbstractMatrix} = decode_data(as, data) function encode_structure_matrix(as::AffineScaler, structure_matrix::MM) where {MM <: AbstractMatrix} return Diagonal(1 ./ get_scale(as) .^ 2) * structure_matrix @@ -299,8 +297,6 @@ end initialize_processor!(ss::Standardizer, data::MM, structure_matrix) where {MM <: AbstractMatrix} = initialize_processor!(ss, data) -encode_data(ss::Standardizer, data::MM, structure_matrix) where {MM <: AbstractMatrix} = encode_data(ss, data) -decode_data(ss::Standardizer, data::MM, structure_matrix) where {MM <: AbstractMatrix} = decode_data(ss, data) function encode_structure_matrix(ss::Standardizer, structure_matrix::MM) where {MM <: AbstractMatrix} encoder_mat = get_encoder_mat(ss)[1] @@ -411,8 +407,6 @@ function decode_data(dd::Decorrelater, data::MM) where {MM <: AbstractMatrix} decoder_mat = get_decoder_mat(dd)[1] return decoder_mat * data .+ data_mean end -encode_data(dd::Decorrelater, data::MM, structure_matrix) where {MM <: AbstractMatrix} = encode_data(dd, data) -decode_data(dd::Decorrelater, data::MM, structure_matrix) where {MM <: AbstractMatrix} = decode_data(dd, data) function encode_structure_matrix(dd::Decorrelater, structure_matrix::MM) where {MM <: AbstractMatrix} encoder_mat = get_encoder_mat(dd)[1] @@ -444,7 +438,7 @@ function initialize_and_encode_data!( obs_noise_cov::USorM, ) where {DCP <: DataContainerProcessor, USorM <: Union{UniformScaling, AbstractMatrix}, MM <: AbstractMatrix} initialize_processor!(dcp, data, obs_noise_cov) - return encode_data(dcp, data, obs_noise_cov) + return encode_data(dcp, data) end initialize_and_encode_data!( @@ -462,14 +456,12 @@ initialize_and_encode_data!( decode_data( dcp::DCP, data::MM, - obs_noise_cov::USorM, apply_to::AS, ) where { DCP <: DataContainerProcessor, MM <: AbstractMatrix, - USorM <: Union{UniformScaling, AbstractMatrix}, AS <: AbstractString, -} = decode_data(dcp, data, obs_noise_cov) +} = decode_data(dcp, data) function initialize_and_encode_data!( dcp::PDCP, @@ -479,17 +471,16 @@ function initialize_and_encode_data!( ) where {PDCP <: PairedDataContainerProcessor, USorM <: Union{UniformScaling, AbstractMatrix}, AS <: AbstractString} input_data, output_data = data initialize_processor!(dcp, input_data, output_data, obs_noise_cov, apply_to) - return encode_data(dcp, input_data, output_data, obs_noise_cov, apply_to) + return encode_data(dcp, input_data, output_data, apply_to) end function decode_data( dcp::PDCP, data, - obs_noise_cov::USorM, apply_to::AS, -) where {PDCP <: PairedDataContainerProcessor, USorM <: Union{UniformScaling, AbstractMatrix}, AS <: AbstractString} +) where {PDCP <: PairedDataContainerProcessor, AS <: AbstractString} input_data, output_data = data - return decode_data(dcp, input_data, output_data, obs_noise_cov, apply_to) + return decode_data(dcp, input_data, output_data, apply_to) end """ @@ -563,6 +554,8 @@ function create_encoder_schedule(schedule_in::VV) where {VV <: AbstractVector} return encoder_schedule end +# Functions to encode/decode with uninitialized schedule (require structure matrices as input) + function encode_with_schedule( encoder_schedule::VV, io_pairs::PDC, @@ -581,7 +574,7 @@ function encode_with_schedule( # apply_to is the string "in", "out" etc. for (processor, extract_data, apply_to) in encoder_schedule - @info "Encoding data: \"$(apply_to)\" with $(processor)" + @info "Initialize encoding of data: \"$(apply_to)\" with $(processor)" if apply_to == "in" structure_matrix = processed_prior_cov elseif apply_to == "out" @@ -600,6 +593,61 @@ function encode_with_schedule( return processed_io_pairs, processed_prior_cov, processed_obs_noise_cov end +# Functions to encode/decode with initialized schedule + +function encode_with_schedule( + encoder_schedule::VV, + data_container::DC, + in_or_out::AS, +) where { + VV <: AbstractVector, + DC <: DataContainer, + AS <: AbstractString, +} + + if !(in_or_out ∈ ["in","out"]) + throw(ArgumentError("`in_or_out` must be either \"in\" (data is an input) or \"out\" (data is an output). Received $(in_or_out)")) + end + processed_container = deepcopy(data_container) + + # apply_to is the string "in", "out" etc. + for (processor, extract_data, apply_to) in encoder_schedule + if apply_to == in_or_out + @info "Encoding data: \"$(apply_to)\" with $(processor)" + processed = encode_data(processor, extract_data(processed_io_pairs), apply_to) + processed_container = DataContainer(processed) + + end + end + + return processed_container +end + +function encode_with_schedule( + encoder_schedule::VV, + structure_matrix::USorM, + in_or_out::AS, +) where { + VV <: AbstractVector, + USorM <: Union{UniformScaling, AbstractMatrix}, + AS <: AbstractString, +} + + if !(in_or_out ∈ ["in","out"]) + throw(ArgumentError("`in_or_out` must be either \"in\" (data is an input) or \"out\" (data is an output). Received $(in_or_out)")) + end + processed_structure_matrix = deepcopy(structure_matrix) + + # apply_to is the string "in", "out" etc. + for (processor, extract_data, apply_to) in encoder_schedule + if apply_to == in_or_out + processed_structure_matrix = encode_structure_matrix(processor, processed_structure_matrix) + end + end + + return processed_structure_matrix +end + function decode_with_schedule( encoder_schedule::VV, io_pairs::PDC, @@ -619,20 +667,13 @@ function decode_with_schedule( # apply_to is the string "in", "out" etc. for idx in reverse(eachindex(encoder_schedule)) (processor, extract_data, apply_to) = encoder_schedule[idx] - if apply_to == "in" - structure_matrix = processed_prior_cov - elseif apply_to == "out" - structure_matrix = processed_obs_noise_cov - end - - @info "Decoding data: \"$(apply_to)\" with $(processor)" - processed = decode_data(processor, extract_data(processed_io_pairs), structure_matrix, apply_to) + processed = decode_data(processor, extract_data(processed_io_pairs), apply_to) if apply_to == "in" - processed_prior_cov = decode_structure_matrix(processor, structure_matrix) + processed_prior_cov = decode_structure_matrix(processor, processed_prior_cov) processed_io_pairs = PairedDataContainer(processed, get_outputs(processed_io_pairs)) elseif apply_to == "out" - processed_obs_noise_cov = decode_structure_matrix(processor, structure_matrix) + processed_obs_noise_cov = decode_structure_matrix(processor, processed_obs_noise_cov) processed_io_pairs = PairedDataContainer(get_inputs(processed_io_pairs), processed) end end @@ -640,6 +681,59 @@ function decode_with_schedule( return processed_io_pairs, processed_prior_cov, processed_obs_noise_cov end +function decode_with_schedule( + encoder_schedule::VV, + data_container::DC, + in_or_out::AS, +) where { + VV <: AbstractVector, + DC <: DataContainer, + AS <: AbstractString, +} + + if !(in_or_out ∈ ["in","out"]) + throw(ArgumentError("`in_or_out` must be either \"in\" (data is an input) or \"out\" (data is an output). Received $(in_or_out)")) + end + processed_container = deepcopy(data_container) + + # apply_to is the string "in", "out" etc. + for idx in reverse(eachindex(encoder_schedule)) + (processor, extract_data, apply_to) = encoder_schedule[idx] + if apply_to == in_or_out + processed = decode_data(processor, extract_data(processed_io_pairs), apply_to) + processed_container = DataContainer(processed) + + end + end + + return processed_container +end + +function decode_with_schedule( + encoder_schedule::VV, + structure_matrix::USorM, + in_or_out::AS, +) where { + VV <: AbstractVector, + USorM <: Union{UniformScaling, AbstractMatrix}, + AS <: AbstractString, +} + + if !(in_or_out ∈ ["in","out"]) + throw(ArgumentError("`in_or_out` must be either \"in\" (data is an input) or \"out\" (data is an output). Received $(in_or_out)")) + end + processed_structure_matrix = deepcopy(structure_matrix) + + # apply_to is the string "in", "out" etc. + for idx in reverse(eachindex(encoder_schedule)) + (processor, extract_data, apply_to) = encoder_schedule[idx] + if apply_to == in_or_out + processed_structure_matrix = decode_structure_matrix(processor, processed_structure_matrix) + end + end + + return processed_structure_matrix +end From 270e575f5b95cd4a62315206322f82bc498c0eef Mon Sep 17 00:00:00 2001 From: odunbar Date: Thu, 19 Jun 2025 17:01:01 -0700 Subject: [PATCH 10/25] remove duplicate tests --- test/Utilities/runtests.jl | 81 ++------------------------------------ 1 file changed, 4 insertions(+), 77 deletions(-) diff --git a/test/Utilities/runtests.jl b/test/Utilities/runtests.jl index 3ad5aa8ca..22b56ae60 100644 --- a/test/Utilities/runtests.jl +++ b/test/Utilities/runtests.jl @@ -249,83 +249,6 @@ end end - # multivariate lossless encodings - test_names = ["standardize", "decorrelate"] - multivariate_lossless_tests = [] - - for (name, sch) in zip(test_names, multivariate_lossless_tests) - encoder_schedule = create_encoder_schedule(sch) - (encoded_io_pairs, encoded_prior_cov, encoded_obs_noise_cov) = - encode_with_schedule(encoder_schedule, io_pairs, prior_cov, obs_noise_cov) - - (decoded_io_pairs, decoded_prior_cov, decoded_obs_noise_cov) = - decode_with_schedule(encoder_schedule, encoded_io_pairs, encoded_prior_cov, encoded_obs_noise_cov) - for (enc_dat, dec_dat, test_dat, enc_covv, dec_covv, test_covv, dim) in zip( - (get_inputs(encoded_io_pairs), get_outputs(encoded_io_pairs)), - (get_inputs(decoded_io_pairs), get_outputs(decoded_io_pairs)), - (get_inputs(io_pairs), get_outputs(io_pairs)), - (encoded_prior_cov, encoded_obs_noise_cov), - (decoded_prior_cov, decoded_obs_noise_cov), - (prior_cov, obs_noise_cov), - (in_dim, out_dim), - ) - - if name == "standardize" - pop_mean = mean(enc_dat, dims = 2) - pop_cov = cov(enc_dat, dims = 2) - dimm = size(pop_cov, 1) - @test all(isapprox.(pop_mean, zeros(dimm), atol = tol)) - @test isapprox(norm(pop_cov - I), 0.0, atol = tol * dimm^2) # expect very accurate - @test isapprox(norm(enc_covv - I), 0.0, atol = 0.2 * dimm^2) # expect poorly accurate - - elseif name == "decorrelate" - pop_mean = mean(enc_dat, dims = 2) - pop_cov = cov(enc_dat, dims = 2) - dimm = size(pop_cov, 1) - @test all(isapprox.(pop_mean, zeros(dimm), atol = tol)) - @test isapprox(norm(pop_cov - I), 0.0, atol = 0.2 * dimm^2) # expect poorly accurate - @test isapprox(norm(enc_covv - I), 0.0, atol = tol * dimm^2) # expect very accurate - end - - @test isapprox(norm(dec_dat - test_dat), 0.0, atol = tol * dim) - @test isapprox(norm(dec_covv - test_covv), 0.0, atol = tol * dim^2) - - end - end - - test_names = ["standardize-truncate-to-5", "decorrelate-retain-partial-var", "decorrelate-estimate-cov"] - - multivariate_lossy_tests = [ - [(standardize(5), "in_and_out")], - [(decorrelate(retain_var = 0.95), "in_and_out")], - [(decorrelate(add_estimated_cov = true), "in_and_out")], - ] - - for (name, sch) in zip(test_names, multivariate_lossless_tests) - encoder_schedule = create_encoder_schedule(sch) - (encoded_io_pairs, encoded_prior_cov, encoded_obs_noise_cov) = - encode_with_schedule(encoder_schedule, io_pairs, prior_cov, obs_noise_cov) - - (decoded_io_pairs, decoded_prior_cov, decoded_obs_noise_cov) = - decode_with_schedule(encoder_schedule, encoded_io_pairs, encoded_prior_cov, encoded_obs_noise_cov) - for (enc_dat, dec_dat, test_dat, enc_covv, dec_covv, test_covv, dim) in zip( - (get_inputs(encoded_io_pairs), get_outputs(encoded_io_pairs)), - (get_inputs(decoded_io_pairs), get_outputs(decoded_io_pairs)), - (get_inputs(io_pairs), get_outputs(io_pairs)), - (encoded_prior_cov, encoded_obs_noise_cov), - (decoded_prior_cov, decoded_obs_noise_cov), - (prior_cov, obs_noise_cov), - (in_dim, out_dim), - ) - - - end - end - - - - # LOSSY TESTS HERE - # combine a few lossless encoding schedules schedule_builder = [ @@ -353,4 +276,8 @@ end @test all(isapprox.(prior_cov, decoded_prior_cov, atol = tol)) @test all(isapprox.(obs_noise_cov, decoded_obs_noise_cov, atol = tol)) + + + + end From 0b709765e778873813becca433cdb6afc9ba8e85 Mon Sep 17 00:00:00 2001 From: odunbar Date: Fri, 20 Jun 2025 10:54:06 -0700 Subject: [PATCH 11/25] rng, better tols, more coverage --- test/Utilities/runtests.jl | 49 +++++++++++++++++++++----------------- 1 file changed, 27 insertions(+), 22 deletions(-) diff --git a/test/Utilities/runtests.jl b/test/Utilities/runtests.jl index 22b56ae60..48efdcb16 100644 --- a/test/Utilities/runtests.jl +++ b/test/Utilities/runtests.jl @@ -74,6 +74,9 @@ end @testset "Data Preprocessing" begin + # Seed for pseudo-random number generator + rng = Random.MersenneTwister(4154) + # quick build tests and test getters zs = zscore_scale() mm = minmax_scale() @@ -110,15 +113,16 @@ end @test get_decoder_mat(DD) == [3] # get some data as IO pairs for functional tests + in_dim = 10 out_dim = 50 - samples = 100 + samples = 40 # for full test coverage have samples in_dim < samples < out_dim - x = randn(in_dim, in_dim) + x = randn(rng, in_dim, in_dim) prior_cov = x * x' - in_data = rand(MvNormal(zeros(in_dim), prior_cov), samples) + in_data = rand(rng, MvNormal(zeros(in_dim), prior_cov), samples) obs_noise_cov = [max(5.0 - abs(i - j), 0.0) for i in 1:out_dim, j in 1:out_dim] # [5 4 3 2 1 0 0 ...] off diagonal - out_data = rand(MvNormal(-10 * ones(out_dim), obs_noise_cov), samples) + out_data = rand(rng, MvNormal(-10 * ones(out_dim), obs_noise_cov), samples) io_pairs = PairedDataContainer(in_data, out_data) test_names = [ @@ -206,19 +210,19 @@ end pop_mean = mean(enc_dat, dims = 2) pop_cov = cov(enc_dat, dims = 2) dimm = size(pop_cov, 1) + big_tol = 0.1 if name == "standardize" @test all(isapprox.(pop_mean, zeros(dimm), atol = tol)) @test isapprox(norm(pop_cov - I), 0.0, atol = tol * dimm^2) # expect very accurate - @test isapprox(norm(enc_covv - I), 0.0, atol = 0.2 * dimm^2) # expect poorly accurate - + @test isapprox(norm(enc_covv - I), 0.0, atol = big_tol * dimm^2) # expect poorly accurate elseif name == "decorrelate" @test all(isapprox.(pop_mean, zeros(dimm), atol = tol)) - @test isapprox(norm(pop_cov - I), 0.0, atol = 0.2 * dimm^2) # expect poorly accurate + @test isapprox(norm(pop_cov - I), 0.0, atol = big_tol * dimm^2) # expect poorly accurate @test isapprox(norm(enc_covv - I), 0.0, atol = tol * dimm^2) # expect very accurate elseif name == "decorrelate-estimate-cov" @test all(isapprox.(pop_mean, zeros(dimm), atol = tol)) - @test isapprox(norm(pop_cov - I), 0.0, atol = 0.2 * dimm^2) # expect poorly accurate - @test isapprox(norm(enc_covv - I), 0.0, atol = 0.2 * dimm^2) # expect poorly accurate + @test isapprox(norm(pop_cov - I), 0.0, atol = big_tol * dimm^2) # expect poorly accurate + @test isapprox(norm(enc_covv - I), 0.0, atol = big_tol * dimm^2) # expect poorly accurate end # Multivariate lossy dim-reduction tests @@ -226,31 +230,36 @@ end @test dimm == 5 @test all(isapprox.(pop_mean, zeros(dimm), atol = tol)) @test isapprox(norm(pop_cov - I), 0.0, atol = tol * dimm^2) # expect very accurate - @test isapprox(norm(enc_covv - I), 0.0, atol = 0.2 * dimm^2) # expect poorly accurate - + @test isapprox(norm(enc_covv - I), 0.0, atol = big_tol * dimm^2) # expect poorly accurate elseif name == "decorrelate-retain-0.95-var" svdc = svd(test_covv) var_cumsum = cumsum(svdc.S .^ 2) ./ sum(svdc.S .^ 2) @test var_cumsum[dimm] > 0.95 @test var_cumsum[dimm - 1] < 0.95 @test all(isapprox.(pop_mean, zeros(dimm), atol = tol)) - @test isapprox(norm(pop_cov - I), 0.0, atol = 0.2 * dimm^2) # expect poorly accurate + @test isapprox(norm(pop_cov - I), 0.0, atol = big_tol * dimm^2) # expect poorly accurate @test isapprox(norm(enc_covv - I), 0.0, atol = tol * dimm^2) # expect very accurate end # test decode approximation of lossless options if ll_flag - @test isapprox(norm(dec_dat - test_dat), 0.0, atol = tol * dim) - @test isapprox(norm(dec_covv - test_covv), 0.0, atol = tol * dim^2) + # when dimm < dim, loss can occur in some tests + tol1 = (name == "decorrelate" && dimm < dim) ? big_tol : tol + tol2 = (name == "standardize" && dimm < dim) ? big_tol : tol + @test isapprox(norm(dec_dat - test_dat), 0.0, atol = tol1 * dim * samples) + @test isapprox(norm(dec_covv - test_covv), 0.0, atol = tol2 * dim^2) end end end + # combine a few lossless encoding schedules (lossless requires samples>dims) + samples = 150 # for full test coverage have samples in_dim < samples < out_dim + in_data = rand(MvNormal(zeros(in_dim), prior_cov), samples) + out_data = rand(MvNormal(-10 * ones(out_dim), obs_noise_cov), samples) + io_pairs = PairedDataContainer(in_data, out_data) - - # combine a few lossless encoding schedules schedule_builder = [ (zscore_scale(), "in_and_out"), # (quartile_scale(), "in"), @@ -273,11 +282,7 @@ end tol = 1e-12 @test all(isapprox.(get_inputs(io_pairs), get_inputs(decoded_io_pairs), atol = tol)) @test all(isapprox.(get_outputs(io_pairs), get_outputs(decoded_io_pairs), atol = tol)) - @test all(isapprox.(prior_cov, decoded_prior_cov, atol = tol)) - @test all(isapprox.(obs_noise_cov, decoded_obs_noise_cov, atol = tol)) - + @test isapprox(norm(prior_cov - decoded_prior_cov), 0.0, atol = tol*in_dim^2) + @test isapprox(norm(obs_noise_cov - decoded_obs_noise_cov), 0.0, atol = tol*out_dim^2) - - - end From 6d97354682fffb5b3fa805dc9dea3f26c2a2284a Mon Sep 17 00:00:00 2001 From: odunbar Date: Fri, 20 Jun 2025 11:33:46 -0700 Subject: [PATCH 12/25] added new encode/decode methods for data or matrix only --- src/Utilities.jl | 103 +++++++++++++++++-------------------- test/Utilities/runtests.jl | 38 ++++++++++++-- 2 files changed, 81 insertions(+), 60 deletions(-) diff --git a/src/Utilities.jl b/src/Utilities.jl index bf110b6a8..f5d2f367c 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -30,7 +30,7 @@ export get_type, get_add_estimated_cov export create_encoder_schedule, encode_with_schedule, decode_with_schedule export initialize_processor!, - initialize_and_encode_data!, encode_data, decode_data, encode_obs_noise_cov, decode_obs_noise_cov + initialize_and_encode_data!, encode_data, decode_data, encode_structure_matrix, decode_structure_matrix @@ -435,50 +435,42 @@ struct CanonicalCorrelationReducer <: PairedDataContainerProcessor end function initialize_and_encode_data!( dcp::DCP, data::MM, - obs_noise_cov::USorM, + structure_mat::USorM, ) where {DCP <: DataContainerProcessor, USorM <: Union{UniformScaling, AbstractMatrix}, MM <: AbstractMatrix} - initialize_processor!(dcp, data, obs_noise_cov) + initialize_processor!(dcp, data, structure_mat) return encode_data(dcp, data) end initialize_and_encode_data!( dcp::DCP, data::MM, - obs_noise_cov::USorM, + structure_mat::USorM, apply_to::AS, ) where { DCP <: DataContainerProcessor, MM <: AbstractMatrix, USorM <: Union{UniformScaling, AbstractMatrix}, AS <: AbstractString, -} = initialize_and_encode_data!(dcp, data, obs_noise_cov) +} = initialize_and_encode_data!(dcp, data, structure_mat) decode_data( dcp::DCP, data::MM, apply_to::AS, -) where { - DCP <: DataContainerProcessor, - MM <: AbstractMatrix, - AS <: AbstractString, -} = decode_data(dcp, data) +) where {DCP <: DataContainerProcessor, MM <: AbstractMatrix, AS <: AbstractString} = decode_data(dcp, data) function initialize_and_encode_data!( dcp::PDCP, data, - obs_noise_cov::USorM, + structure_mat::USorM, apply_to::AS, ) where {PDCP <: PairedDataContainerProcessor, USorM <: Union{UniformScaling, AbstractMatrix}, AS <: AbstractString} input_data, output_data = data - initialize_processor!(dcp, input_data, output_data, obs_noise_cov, apply_to) + initialize_processor!(dcp, input_data, output_data, structure_mat, apply_to) return encode_data(dcp, input_data, output_data, apply_to) end -function decode_data( - dcp::PDCP, - data, - apply_to::AS, -) where {PDCP <: PairedDataContainerProcessor, AS <: AbstractString} +function decode_data(dcp::PDCP, data, apply_to::AS) where {PDCP <: PairedDataContainerProcessor, AS <: AbstractString} input_data, output_data = data return decode_data(dcp, input_data, output_data, apply_to) end @@ -599,24 +591,23 @@ function encode_with_schedule( encoder_schedule::VV, data_container::DC, in_or_out::AS, -) where { - VV <: AbstractVector, - DC <: DataContainer, - AS <: AbstractString, -} - - if !(in_or_out ∈ ["in","out"]) - throw(ArgumentError("`in_or_out` must be either \"in\" (data is an input) or \"out\" (data is an output). Received $(in_or_out)")) +) where {VV <: AbstractVector, DC <: DataContainer, AS <: AbstractString} + + if !(in_or_out ∈ ["in", "out"]) + throw( + ArgumentError( + "`in_or_out` must be either \"in\" (data is an input) or \"out\" (data is an output). Received $(in_or_out)", + ), + ) end processed_container = deepcopy(data_container) # apply_to is the string "in", "out" etc. for (processor, extract_data, apply_to) in encoder_schedule if apply_to == in_or_out - @info "Encoding data: \"$(apply_to)\" with $(processor)" - processed = encode_data(processor, extract_data(processed_io_pairs), apply_to) + processed = encode_data(processor, get_data(processed_container)) processed_container = DataContainer(processed) - + end end @@ -627,21 +618,21 @@ function encode_with_schedule( encoder_schedule::VV, structure_matrix::USorM, in_or_out::AS, -) where { - VV <: AbstractVector, - USorM <: Union{UniformScaling, AbstractMatrix}, - AS <: AbstractString, -} - - if !(in_or_out ∈ ["in","out"]) - throw(ArgumentError("`in_or_out` must be either \"in\" (data is an input) or \"out\" (data is an output). Received $(in_or_out)")) +) where {VV <: AbstractVector, USorM <: Union{UniformScaling, AbstractMatrix}, AS <: AbstractString} + + if !(in_or_out ∈ ["in", "out"]) + throw( + ArgumentError( + "`in_or_out` must be either \"in\" (data is an input) or \"out\" (data is an output). Received $(in_or_out)", + ), + ) end processed_structure_matrix = deepcopy(structure_matrix) # apply_to is the string "in", "out" etc. for (processor, extract_data, apply_to) in encoder_schedule if apply_to == in_or_out - processed_structure_matrix = encode_structure_matrix(processor, processed_structure_matrix) + processed_structure_matrix = encode_structure_matrix(processor, processed_structure_matrix) end end @@ -685,14 +676,14 @@ function decode_with_schedule( encoder_schedule::VV, data_container::DC, in_or_out::AS, -) where { - VV <: AbstractVector, - DC <: DataContainer, - AS <: AbstractString, -} - - if !(in_or_out ∈ ["in","out"]) - throw(ArgumentError("`in_or_out` must be either \"in\" (data is an input) or \"out\" (data is an output). Received $(in_or_out)")) +) where {VV <: AbstractVector, DC <: DataContainer, AS <: AbstractString} + + if !(in_or_out ∈ ["in", "out"]) + throw( + ArgumentError( + "`in_or_out` must be either \"in\" (data is an input) or \"out\" (data is an output). Received $(in_or_out)", + ), + ) end processed_container = deepcopy(data_container) @@ -700,9 +691,9 @@ function decode_with_schedule( for idx in reverse(eachindex(encoder_schedule)) (processor, extract_data, apply_to) = encoder_schedule[idx] if apply_to == in_or_out - processed = decode_data(processor, extract_data(processed_io_pairs), apply_to) + processed = decode_data(processor, get_data(processed_container)) processed_container = DataContainer(processed) - + end end @@ -713,14 +704,14 @@ function decode_with_schedule( encoder_schedule::VV, structure_matrix::USorM, in_or_out::AS, -) where { - VV <: AbstractVector, - USorM <: Union{UniformScaling, AbstractMatrix}, - AS <: AbstractString, -} - - if !(in_or_out ∈ ["in","out"]) - throw(ArgumentError("`in_or_out` must be either \"in\" (data is an input) or \"out\" (data is an output). Received $(in_or_out)")) +) where {VV <: AbstractVector, USorM <: Union{UniformScaling, AbstractMatrix}, AS <: AbstractString} + + if !(in_or_out ∈ ["in", "out"]) + throw( + ArgumentError( + "`in_or_out` must be either \"in\" (data is an input) or \"out\" (data is an output). Received $(in_or_out)", + ), + ) end processed_structure_matrix = deepcopy(structure_matrix) @@ -728,7 +719,7 @@ function decode_with_schedule( for idx in reverse(eachindex(encoder_schedule)) (processor, extract_data, apply_to) = encoder_schedule[idx] if apply_to == in_or_out - processed_structure_matrix = decode_structure_matrix(processor, processed_structure_matrix) + processed_structure_matrix = decode_structure_matrix(processor, processed_structure_matrix) end end diff --git a/test/Utilities/runtests.jl b/test/Utilities/runtests.jl index 48efdcb16..207c3b676 100644 --- a/test/Utilities/runtests.jl +++ b/test/Utilities/runtests.jl @@ -74,7 +74,7 @@ end @testset "Data Preprocessing" begin - # Seed for pseudo-random number generator + # Seed for pseudo-random number generator rng = Random.MersenneTwister(4154) # quick build tests and test getters @@ -244,7 +244,7 @@ end # test decode approximation of lossless options if ll_flag # when dimm < dim, loss can occur in some tests - tol1 = (name == "decorrelate" && dimm < dim) ? big_tol : tol + tol1 = (name == "decorrelate" && dimm < dim) ? big_tol : tol tol2 = (name == "standardize" && dimm < dim) ? big_tol : tol @test isapprox(norm(dec_dat - test_dat), 0.0, atol = tol1 * dim * samples) @test isapprox(norm(dec_covv - test_covv), 0.0, atol = tol2 * dim^2) @@ -282,7 +282,37 @@ end tol = 1e-12 @test all(isapprox.(get_inputs(io_pairs), get_inputs(decoded_io_pairs), atol = tol)) @test all(isapprox.(get_outputs(io_pairs), get_outputs(decoded_io_pairs), atol = tol)) - @test isapprox(norm(prior_cov - decoded_prior_cov), 0.0, atol = tol*in_dim^2) - @test isapprox(norm(obs_noise_cov - decoded_obs_noise_cov), 0.0, atol = tol*out_dim^2) + @test isapprox(norm(prior_cov - decoded_prior_cov), 0.0, atol = tol * in_dim^2) + @test isapprox(norm(obs_noise_cov - decoded_obs_noise_cov), 0.0, atol = tol * out_dim^2) + + # enc/dec just data + samples = 1 # try one sample + new_in_data = rand(rng, MvNormal(zeros(in_dim), prior_cov), samples) + new_out_data = rand(rng, MvNormal(-10 * ones(out_dim), obs_noise_cov), samples) + id = DataContainer(new_in_data) + od = DataContainer(new_out_data) + + encoded_id = encode_with_schedule(encoder_schedule, id, "in") + encoded_od = encode_with_schedule(encoder_schedule, od, "out") + decoded_id = decode_with_schedule(encoder_schedule, encoded_id, "in") + decoded_od = decode_with_schedule(encoder_schedule, encoded_od, "out") + # tests in latent space? + @test isapprox(norm(get_data(decoded_id) - get_data(id)), 0.0, atol = tol * in_dim * samples) + @test isapprox(norm(get_data(decoded_od) - get_data(od)), 0.0, atol = tol * out_dim * samples) + @test_throws ArgumentError encode_with_schedule(encoder_schedule, id, "bad") + @test_throws ArgumentError decode_with_schedule(encoder_schedule, id, "bad") + + # enc/dec just mats + encoded_pc = encode_with_schedule(encoder_schedule, prior_cov, "in") + encoded_oc = encode_with_schedule(encoder_schedule, obs_noise_cov, "out") + decoded_pc = decode_with_schedule(encoder_schedule, encoded_pc, "in") + decoded_oc = decode_with_schedule(encoder_schedule, encoded_oc, "out") + @test isapprox(norm(encoded_pc - encoded_prior_cov), 0.0, atol = tol * in_dim^2) + @test isapprox(norm(encoded_oc - encoded_obs_noise_cov), 0.0, atol = tol * out_dim^2) + @test isapprox(norm(decoded_pc - decoded_prior_cov), 0.0, atol = tol * in_dim^2) + @test isapprox(norm(decoded_oc - decoded_obs_noise_cov), 0.0, atol = tol * out_dim^2) + @test_throws ArgumentError encode_with_schedule(encoder_schedule, prior_cov, "bad") + @test_throws ArgumentError decode_with_schedule(encoder_schedule, encoded_pc, "bad") + end From 6db49f99d6f7ac9d128d650d70ec97f5c008c003 Mon Sep 17 00:00:00 2001 From: odunbar Date: Fri, 20 Jun 2025 18:57:48 -0700 Subject: [PATCH 13/25] added CCA --- src/Utilities.jl | 203 ++++++++++++++++++++++++++++++++----- test/Utilities/runtests.jl | 39 +++++-- 2 files changed, 207 insertions(+), 35 deletions(-) diff --git a/src/Utilities.jl b/src/Utilities.jl index f5d2f367c..39aed2746 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -17,7 +17,6 @@ export zscore2orig export PairedDataContainerProcessor, DataContainerProcessor export UnivariateAffineScaling, AffineScaler, QuartileScaling, MinMaxScaling, ZScoreScaling, Standardizer, Decorrelater - export quartile_scale, minmax_scale, zscore_scale, standardize, decorrelate export get_type, get_shift, @@ -32,6 +31,10 @@ export create_encoder_schedule, encode_with_schedule, decode_with_schedule export initialize_processor!, initialize_and_encode_data!, encode_data, decode_data, encode_structure_matrix, decode_structure_matrix +export CanonicalCorrelation +export canonical_correlation +export get_apply_to + """ @@ -266,18 +269,13 @@ function initialize_processor!(ss::Standardizer, data::MM) where {MM <: Abstract end if length(get_encoder_mat(ss)) == 0 # can use tsvd here so we only compute up to rk - data_cov = cov(data, dims = 2) - svdc = svd(data_cov) - rk = min(rank(data_cov), max(get_rank(ss), 1)) - # Want to use the nonsquare singular vector - if size(svdc.U, 1) == size(svdc.U, 2) - mat = svdc.Vt - else - mat = svdc.U' - end + cov_data = cov(data, dims = 2) + svdc = svd(cov_data) + rk = min(rank(cov_data), max(get_rank(ss), 1)) - encoder_mat = Diagonal(1 ./ sqrt.(svdc.S[1:rk])) * mat[1:rk, :] - decoder_mat = mat[1:rk, :]' * Diagonal(sqrt.(svdc.S[1:rk])) + # rescale with sqrt(size(data,2)) to normalize + encoder_mat = Diagonal(1 ./ sqrt.(svdc.S[1:rk])) * svdc.Vt[1:rk, :] + decoder_mat = svdc.Vt[1:rk, :]' * Diagonal(sqrt.(svdc.S[1:rk])) push!(get_encoder_mat(ss), encoder_mat) push!(get_decoder_mat(ss), decoder_mat) end @@ -382,14 +380,9 @@ function initialize_processor!( sqrt_inv_sv = Diagonal(1.0 ./ sqrt.(svdA.S[1:trunc])) sqrt_sv = Diagonal(sqrt.(svdA.S[1:trunc])) - if size(svdA.U, 1) == size(svdA.U, 2) - mat = svdA.Vt - else - mat = svdA.U - end - - encoder_mat = sqrt_inv_sv * mat[1:trunc, :] - decoder_mat = mat[1:trunc, :]' * sqrt_sv + # as we have svd of cov-matrix we can use U or Vt + encoder_mat = sqrt_inv_sv * svdA.Vt[1:trunc, :] + decoder_mat = svdA.Vt[1:trunc, :]' * sqrt_sv push!(get_encoder_mat(dd), encoder_mat) push!(get_decoder_mat(dd), decoder_mat) @@ -421,17 +414,166 @@ end # ...struct VariationalAutoEncoder <: DataContainerProcessor end # PDCProcessors -struct DataInformedReducer <: PairedDataContainerProcessor end +# struct InverseProblemInformed <: PairedDataContainerProcessor end +# struct LikelihoodInformed <: PairedDataContainerProcessor end + + +""" +$(TYPEDEF) + +Uses both input and output data to learn a subspace of maximal correlation between inputs and outputs. The subspace for a pair (X,Y) will be of size minimum(rank(X),rank(Y)), computed using SVD-based method +e.g. See e.g., https://numerical.recipes/whp/notes/CanonCorrBySVD.pdf + +# Fields +$(TYPEDFIELDS) +""" +struct CanonicalCorrelation{VV1, VV2, VV3, FT, VV4} <: PairedDataContainerProcessor + "storage for the input or output data mean" + data_mean::VV1 + "the encoding matrix of input or output canonical correlations" + encoder_mat::VV2 + "the decoding matrix of input or output canonical correlations" + decoder_mat::VV3 + "the fraction of variance to be retained after truncating singular values (1 implies no truncation)" + retain_var::FT + "Stores whether this is an input or output encoder (vector with string \"in\" or \"out\")" + apply_to::VV4 +end + +canonical_correlation(; retain_var = Float64(1.0)) = + CanonicalCorrelation(Any[], Any[], Any[], retain_var, AbstractString[]) + +get_data_mean(cc::CanonicalCorrelation) = cc.data_mean +get_encoder_mat(cc::CanonicalCorrelation) = cc.encoder_mat +get_decoder_mat(cc::CanonicalCorrelation) = cc.decoder_mat +get_retain_var(cc::CanonicalCorrelation) = cc.retain_var +get_apply_to(cc::CanonicalCorrelation) = cc.apply_to + +function Base.show(io::IO, cc::CanonicalCorrelation) + if get_retain_var(cc) < 1.0 + out = "CanonicalCorrelation: retain_var=$(get_retain_var(cc))" + else + out = "CanonicalCorrelation" + end + print(io, out) +end + +function initialize_processor!( + cc::CanonicalCorrelation, + in_data::MM, + out_data::MM, + apply_to::AS, +) where {MM <: AbstractMatrix, AS <: AbstractString} + + if length(get_apply_to(cc)) == 0 + push!(get_apply_to(cc), apply_to) + end + + if length(get_data_mean(cc)) == 0 + if apply_to == "in" + append!(get_data_mean(cc), mean(in_data, dims = 2)) + elseif apply_to == "out" + append!(get_data_mean(cc), mean(out_data, dims = 2)) + end + end + + if length(get_encoder_mat(cc)) == 0 + + if size(in_data, 2) < size(in_data, 1) || size(out_data, 2) < size(out_data, 1) + throw( + ArgumentError( + "CanonicalCorrelation implementation not defined for # data samples < dimensions, please obtain more samples, or perform prior dimension reduction approaches until this is satisfied", + ), + ) + end + + # Individually decompose in and out + # Want to use the nonsquare singular vector + svdi = svd(in_data .- mean(in_data, dims = 2)) + svdo = svd(out_data .- mean(out_data, dims = 2)) + # determine regime + + # ensure we non-square sv (in_mat = (in_dim x n_samples), out_mat = (out_dim x n_samples)) + in_mat_sq, in_mat_nonsq = (size(svdi.U, 1) == size(svdi.U, 2)) ? (svdi.U, svdi.Vt) : (svdi.Vt, svdi.U) + out_mat_sq, out_mat_nonsq = (size(svdo.U, 1) == size(svdo.U, 2)) ? (svdo.U, svdo.Vt) : (svdo.Vt, svdo.U) + + svdio = svd(out_mat_nonsq * in_mat_nonsq') + + # retain variance + ret_var = get_retain_var(cc) + if ret_var < 1.0 + sv_cumsum = cumsum(svdio.S .^ 2) / sum(svdio.S .^ 2) # variance contributions are (sing_val)^2 + trunc = minimum(findall(x -> (x > ret_var), sv_cumsum)) + @info " truncating at $(trunc)/$(length(sv_cumsum)) retaining $(100.0*sv_cumsum[trunc])% of the variance in the joint space" + else + trunc = min(rank(in_data), rank(out_data)) + end + + if apply_to == "in" + in_dim = size(in_data, 1) + svdio_mat = (size(svdio.U, 1) == in_dim) ? svdio.U : svdio.V + # mat' * Sx⁻¹ * Uxt + encoder_mat = svdio_mat[:, 1:trunc]' * Diagonal(1 ./ svdi.S) * in_mat_sq' + decoder_mat = in_mat_sq * Diagonal(svdi.S) * svdio_mat[:, 1:trunc] + else + apply_to == "out" + out_dim = size(out_data, 1) + svdio_mat = (size(svdio.U, 1) == out_dim) ? svdio.U : svdio.V + # Vt * Sy⁻¹ * Uyt + encoder_mat = svdio_mat[:, 1:trunc]' * Diagonal(1 ./ svdo.S) * out_mat_sq' + decoder_mat = out_mat_sq * Diagonal(svdo.S) * svdio_mat[:, 1:trunc] + end + + push!(get_encoder_mat(cc), encoder_mat) + push!(get_decoder_mat(cc), decoder_mat) + + # Note: To check CCA: + # u = in_encoder * (in - mean(in)) + # v = out_encoder * (out - mean(out)) + # u * u' = v * v' = I, + # v * u' = u * v' = Diagonal(svdio.S[1:trunc]) + end +end + +initialize_processor!( + cc::CanonicalCorrelation, + in_data::MM, + out_data::MM, + structure_matrix, + apply_to::AS, +) where {MM <: AbstractMatrix, AS <: AbstractString} = initialize_processor!(cc, in_data, out_data, apply_to) + + +function encode_data(cc::CanonicalCorrelation, data::MM) where {MM <: AbstractMatrix} + data_mean = get_data_mean(cc) + encoder_mat = get_encoder_mat(cc)[1] + return encoder_mat * (data .- data_mean) +end + +function decode_data(cc::CanonicalCorrelation, data::MM) where {MM <: AbstractMatrix} + data_mean = get_data_mean(cc) + decoder_mat = get_decoder_mat(cc)[1] + return decoder_mat * data .+ data_mean +end + +function encode_structure_matrix(cc::CanonicalCorrelation, structure_matrix::MM) where {MM <: AbstractMatrix} + encoder_mat = get_encoder_mat(cc)[1] + return encoder_mat * structure_matrix * encoder_mat' +end + +function decode_structure_matrix(cc::CanonicalCorrelation, enc_structure_matrix::MM) where {MM <: AbstractMatrix} + decoder_mat = get_decoder_mat(cc)[1] + return decoder_mat * enc_structure_matrix * decoder_mat' +end + -struct LikelihoodInformedReducer <: PairedDataContainerProcessor end -struct CanonicalCorrelationReducer <: PairedDataContainerProcessor end #### -# generic function to build and encode data +# generic functions to initialize, encode, and decode data. With Data processors or PairedData processors function initialize_and_encode_data!( dcp::DCP, data::MM, @@ -467,12 +609,20 @@ function initialize_and_encode_data!( ) where {PDCP <: PairedDataContainerProcessor, USorM <: Union{UniformScaling, AbstractMatrix}, AS <: AbstractString} input_data, output_data = data initialize_processor!(dcp, input_data, output_data, structure_mat, apply_to) - return encode_data(dcp, input_data, output_data, apply_to) + if apply_to == "in" + return encode_data(dcp, input_data) + elseif apply_to == "out" + return encode_data(dcp, output_data) + end end function decode_data(dcp::PDCP, data, apply_to::AS) where {PDCP <: PairedDataContainerProcessor, AS <: AbstractString} input_data, output_data = data - return decode_data(dcp, input_data, output_data, apply_to) + if apply_to == "in" + return decode_data(dcp, input_data) + elseif apply_to == "out" + return decode_data(dcp, output_data) + end end """ @@ -559,7 +709,6 @@ function encode_with_schedule( USorM1 <: Union{UniformScaling, AbstractMatrix}, USorM2 <: Union{UniformScaling, AbstractMatrix}, } - processed_io_pairs = deepcopy(io_pairs) processed_prior_cov = deepcopy(prior_cov_in) processed_obs_noise_cov = deepcopy(obs_noise_cov_in) diff --git a/test/Utilities/runtests.jl b/test/Utilities/runtests.jl index 207c3b676..298e2aed8 100644 --- a/test/Utilities/runtests.jl +++ b/test/Utilities/runtests.jl @@ -116,7 +116,7 @@ end in_dim = 10 out_dim = 50 - samples = 40 # for full test coverage have samples in_dim < samples < out_dim + samples = 120 # for full test coverage have samples in_dim < samples < out_dim x = randn(rng, in_dim, in_dim) prior_cov = x * x' @@ -125,36 +125,36 @@ end out_data = rand(rng, MvNormal(-10 * ones(out_dim), obs_noise_cov), samples) io_pairs = PairedDataContainer(in_data, out_data) - test_names = [ + test_names = [ # order as in tests "zscore", "quartile", "minmax", "standardize", "decorrelate", "decorrelate-estimate-cov", + "canonical-correlation", "standardize-truncate-to-5", "decorrelate-retain-0.95-var", ] - - # Test encodings-decodings individually - univariate_tests = [ + schedules = [ [(zscore_scale(), "in_and_out")], [(quartile_scale(), "in_and_out")], [(minmax_scale(), "in_and_out")], [(standardize(), "in_and_out")], [(decorrelate(), "in_and_out")], [(decorrelate(add_estimated_cov = true), "in_and_out")], + [(canonical_correlation(), "in_and_out")], [(standardize(5), "in_and_out")], [(decorrelate(retain_var = 0.95), "in_and_out")], ] - lossless = [fill(true, 6); fill(false, 2)] # are these lossy approximations + lossless = [fill(true, 6); fill(false, 3)] # are these lossy approximations? # functional test pipeline tol = 1e-12 - for (name, sch, ll_flag) in zip(test_names, univariate_tests, lossless) + for (name, sch, ll_flag) in zip(test_names, schedules, lossless) encoder_schedule = create_encoder_schedule(sch) (encoded_io_pairs, encoded_prior_cov, encoded_obs_noise_cov) = encode_with_schedule(encoder_schedule, io_pairs, prior_cov, obs_noise_cov) @@ -214,7 +214,7 @@ end if name == "standardize" @test all(isapprox.(pop_mean, zeros(dimm), atol = tol)) @test isapprox(norm(pop_cov - I), 0.0, atol = tol * dimm^2) # expect very accurate - @test isapprox(norm(enc_covv - I), 0.0, atol = big_tol * dimm^2) # expect poorly accurate + @test isapprox(norm(enc_covv - I), 0.0, atol = big_tol * dimm^2) # expect poorly accurate, particularly if dimm < dim elseif name == "decorrelate" @test all(isapprox.(pop_mean, zeros(dimm), atol = tol)) @test isapprox(norm(pop_cov - I), 0.0, atol = big_tol * dimm^2) # expect poorly accurate @@ -241,6 +241,26 @@ end @test isapprox(norm(enc_covv - I), 0.0, atol = tol * dimm^2) # expect very accurate end + # Paired data processor reduction: + if name == "canonical-correlation" + @test dimm == min(rank(get_inputs(io_pairs)), rank(get_outputs(io_pairs))) + @test isapprox(norm(enc_dat * enc_dat' - I), 0.0, atol = tol * dimm^2) # test in or out orthogonality + + # check cross-orthogonality is diagonal (nb this test will be duplicate) + enc_in = get_inputs(encoded_io_pairs) + enc_out = get_outputs(encoded_io_pairs) + @test isapprox(norm(enc_in * enc_out' - Diagonal(diag(enc_in * enc_out'))), 0.0, atol = tol * dimm^2) # test cross - orthogonality + + @test isapprox(norm(enc_out * enc_in' - Diagonal(diag(enc_out * enc_in'))), 0.0, atol = tol * dimm^2) # test cross - orthogonality + + # decoder test is lossless only for smaller dimension + if dim == min(in_dim, out_dim) + @test isapprox(norm(dec_dat - test_dat), 0.0, atol = tol * dim * samples) + @test isapprox(norm(dec_covv - test_covv), 0.0, atol = tol * dim^2) + end + + end + # test decode approximation of lossless options if ll_flag # when dimm < dim, loss can occur in some tests @@ -248,6 +268,7 @@ end tol2 = (name == "standardize" && dimm < dim) ? big_tol : tol @test isapprox(norm(dec_dat - test_dat), 0.0, atol = tol1 * dim * samples) @test isapprox(norm(dec_covv - test_covv), 0.0, atol = tol2 * dim^2) + end end @@ -315,4 +336,6 @@ end @test_throws ArgumentError decode_with_schedule(encoder_schedule, encoded_pc, "bad") + + end From 1c3cbda320463b1be35709d55692eb1b31e09485 Mon Sep 17 00:00:00 2001 From: odunbar Date: Tue, 24 Jun 2025 10:24:38 -0700 Subject: [PATCH 14/25] API and codecov --- src/Utilities.jl | 7 +++++-- test/Utilities/runtests.jl | 26 ++++++++++++++++---------- 2 files changed, 21 insertions(+), 12 deletions(-) diff --git a/src/Utilities.jl b/src/Utilities.jl index 39aed2746..d540e3fd9 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -672,7 +672,7 @@ function create_encoder_schedule(schedule_in::VV) where {VV <: AbstractVector} push!(encoder_schedule, (deepcopy(processor), func2, "out")) else @warn( - "Expected schedule keywords ∈ {\"in\",\"out\",\"in_and_out\"}. Received $(str), ignoring processor $(processor)..." + "Expected schedule keywords ∈ {\"in\",\"out\",\"in_and_out\"}. Received $(apply_to), ignoring processor $(processor)..." ) end @@ -687,7 +687,7 @@ function create_encoder_schedule(schedule_in::VV) where {VV <: AbstractVector} push!(encoder_schedule, (deepcopy(processor), func, "out")) else @warn( - "Expected schedule keywords ∈ {\"in\",\"out\",\"in_and_out\"}. Received $(str), ignoring processor $(processor)..." + "Expected schedule keywords ∈ {\"in\",\"out\",\"in_and_out\"}. Received $(apply_to), ignoring processor $(processor)..." ) end end @@ -696,6 +696,9 @@ function create_encoder_schedule(schedule_in::VV) where {VV <: AbstractVector} return encoder_schedule end +create_encoder_schedule(schedule_in::TT) where {TT <: Tuple} = create_encoder_schedule([schedule_in]) + + # Functions to encode/decode with uninitialized schedule (require structure matrices as input) function encode_with_schedule( diff --git a/test/Utilities/runtests.jl b/test/Utilities/runtests.jl index 298e2aed8..419aefb8c 100644 --- a/test/Utilities/runtests.jl +++ b/test/Utilities/runtests.jl @@ -135,21 +135,23 @@ end "canonical-correlation", "standardize-truncate-to-5", "decorrelate-retain-0.95-var", + "canonical-correlation-0.95-var", ] # Test encodings-decodings individually schedules = [ - [(zscore_scale(), "in_and_out")], - [(quartile_scale(), "in_and_out")], - [(minmax_scale(), "in_and_out")], - [(standardize(), "in_and_out")], - [(decorrelate(), "in_and_out")], - [(decorrelate(add_estimated_cov = true), "in_and_out")], - [(canonical_correlation(), "in_and_out")], - [(standardize(5), "in_and_out")], - [(decorrelate(retain_var = 0.95), "in_and_out")], + (zscore_scale(), "in_and_out"), + (quartile_scale(), "in_and_out"), + (minmax_scale(), "in_and_out"), + (standardize(), "in_and_out"), + (decorrelate(), "in_and_out"), + (decorrelate(add_estimated_cov = true), "in_and_out"), + (canonical_correlation(), "in_and_out"), + (standardize(5), "in_and_out"), + (decorrelate(retain_var = 0.95), "in_and_out"), + (canonical_correlation(retain_var = 0.95), "in_and_out"), ] - lossless = [fill(true, 6); fill(false, 3)] # are these lossy approximations? + lossless = [fill(true, 6); fill(false, 4)] # are these lossy approximations? # functional test pipeline tol = 1e-12 @@ -287,9 +289,13 @@ end (standardize(), "in_and_out"), (minmax_scale(), "out"), (decorrelate(), "in_and_out"), + (canonical_correlation(), "in"), ] # make schedule more parsable + + @test_logs (:warn,) create_encoder_schedule((canonical_correlation(), "bad")) + @test_logs (:warn,) create_encoder_schedule((zscore_scale(), "bad")) encoder_schedule = create_encoder_schedule(schedule_builder) # encode the data using the schedule From a8eb1180ddc55aefc66f5f3713c8807ae8da3653 Mon Sep 17 00:00:00 2001 From: odunbar Date: Tue, 24 Jun 2025 12:09:56 -0700 Subject: [PATCH 15/25] added many docstrings --- src/Utilities.jl | 325 ++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 319 insertions(+), 6 deletions(-) diff --git a/src/Utilities.jl b/src/Utilities.jl index d540e3fd9..b502de4af 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -129,10 +129,14 @@ abstract type ZScoreScaling <: UnivariateAffineScaling end """ $(TYPEDEF) -The AffineScaler{T} will create an encoding of the data_container via affine transformations: -- quartile_scale() : creates `QuartileScaling`, encoding with the transform (x - Q2(x))/(Q3(x) - Q1(x)) -- minmax_scale() : creates `MinMaxScaling`, encoding with the transform (x - min(x))/(max(x) - min(x)) -- zscore_scale() : creates `ZScoreScaling`, encoding with the (univariate) transform (x-μ)/σ +The AffineScaler{T} will create an encoding of the data_container via affine transformations. + +Different methods `T` will build different transformations: +- [`quartile_scale`](@ref) : creates `QuartileScaling`, +- [`minmax_scale`](@ref) : creates `MinMaxScaling` +- [`zscore_scale`](@ref) : creates `ZScoreScaling` + +and are accessed with [`get_type`](@ref) """ struct AffineScaler{T, VV <: AbstractVector} <: DataContainerProcessor "storage for the shift applied to data" @@ -141,15 +145,54 @@ struct AffineScaler{T, VV <: AbstractVector} <: DataContainerProcessor scale::VV end +""" +$(TYPEDSIGNATURES) + +Constructs `AffineScaler{QuartileScaling}` processor. +As part of an encoder schedule, it will apply the transform ``\\frac{x - Q2(x)}{Q3(x) - Q1(x)}`` to each data dimension. +Also known as "robust scaling" +""" quartile_scale() = AffineScaler(QuartileScaling) + +""" +$(TYPEDSIGNATURES) + +Constructs `AffineScaler{MinMaxScaling}` processor. +As part of an encoder schedule, this will apply the transform ``\\frac{x - \\min(x)}{\\max(x) - \\min(x)}`` to each data dimension. +""" minmax_scale() = AffineScaler(MinMaxScaling) + +""" +$(TYPEDSIGNATURES) + +Constructs `AffineScaler{ZScoreScaling}` processor. +As part of an encoder schedule, this will apply the transform ``\\frac{x-\\mu}{\\sigma}``, (where ``x\\sim N(\\mu,\\sigma)``), to each data dimension. +For multivariate standardization, see `Standardizer` +""" zscore_scale() = AffineScaler(ZScoreScaling) AffineScaler(::Type{UAS}) where {UAS <: UnivariateAffineScaling} = AffineScaler{UAS, Vector{Float64}}(Float64[], Float64[]) +""" +$(TYPEDSIGNATURES) + +Gets the UnivariateAffineScaling type `T` +""" get_type(as::AffineScaler{T, VV}) where {T, VV} = T + +""" +$(TYPEDSIGNATURES) + +Gets the `shift` field of the `AffineScaler` +""" get_shift(as::AffineScaler) = as.shift + +""" +$(TYPEDSIGNATURES) + +Gets the `scale` field of the `AffineScaler` +""" get_scale(as::AffineScaler) = as.scale function Base.show(io::IO, as::AffineScaler) @@ -197,6 +240,11 @@ function initialize_processor!(as::AffineScaler, data::MM) where {MM <: Abstract end end +""" +$(TYPEDSIGNATURES) + +Apply the `AffineScaler` encoder, on a columns-are-data matrix +""" function encode_data(as::AffineScaler, data::MM) where {MM <: AbstractMatrix} out = deepcopy(data) for i in 1:size(out, 1) @@ -206,6 +254,11 @@ function encode_data(as::AffineScaler, data::MM) where {MM <: AbstractMatrix} return out end +""" +$(TYPEDSIGNATURES) + +Apply the `AffineScaler` decoder, on a columns-are-data matrix +""" function decode_data(as::AffineScaler, data::MM) where {MM <: AbstractMatrix} out = deepcopy(data) for i in 1:size(out, 1) @@ -215,13 +268,29 @@ function decode_data(as::AffineScaler, data::MM) where {MM <: AbstractMatrix} return out end +""" +$(TYPEDSIGNATURES) + +Computes and populates the `shift` and `scale` fields for the `AffineScaler` +""" initialize_processor!(as::AffineScaler, data::MM, structure_matrix) where {MM <: AbstractMatrix} = initialize_processor!(as, data) + +""" +$(TYPEDSIGNATURES) + +Apply the `AffineScaler` encoder to a provided structure matrix +""" function encode_structure_matrix(as::AffineScaler, structure_matrix::MM) where {MM <: AbstractMatrix} return Diagonal(1 ./ get_scale(as) .^ 2) * structure_matrix end +""" +$(TYPEDSIGNATURES) + +Apply the `AffineScaler` decoder to a provided structure matrix +""" function decode_structure_matrix(as::AffineScaler, enc_structure_matrix::MM) where {MM <: AbstractMatrix} return Diagonal(get_scale(as) .^ 2) * enc_structure_matrix end @@ -233,7 +302,12 @@ $(TYPEDEF) Standardizes the data to a multivariate N(0,I) distribution via `(xcov)^{-1/2}*(x-x_mean)`, along with rank reduction if data is low rank, or if the user provides a rank. This guarantees that the data samples will have sample mean `0` and covariance `I` after processing. -Similar to the `Decorrelater`, but there the structure matrices will directly be used to perform decorrelation, and will become `I` after processing. In that case, the data samples may not have covariance `I` after processing. +Similar to the [`Decorrelater`](@ref), but there the structure matrices will directly be used to perform decorrelation, and will become `I` after processing. In that case, the data samples may not have covariance `I` after processing. + +Preferred construction is using the [`standardize`](@ref) method. + +# Fields +$(TYPEDFIELDS) """ struct Standardizer{VV1 <: AbstractVector, VV2 <: AbstractVector, VV3 <: AbstractVector} <: DataContainerProcessor "user-provided rank for additional truncation of the space. used if rank < rank(data_cov)" @@ -246,12 +320,46 @@ struct Standardizer{VV1 <: AbstractVector, VV2 <: AbstractVector, VV3 <: Abstrac decoder_mat::VV3 end +""" +$(TYPEDSIGNATURES) + +Constructs the [`Standardizer`](@ref) struct with no truncation of rank +""" standardize() = Standardizer(typemax(Int), Float64[], Any[], Any[]) +""" +$(TYPEDSIGNATURES) + +Constructs the `Standardizer` struct with a provided truncation to rank `rk`. +""" standardize(rk::Int) = Standardizer(rk, Float64[], Any[], Any[]) + +""" +$(TYPEDSIGNATURES) + +returns the `data_mean` field of the `Standardizer`. +""" get_data_mean(ss::Standardizer) = ss.data_mean + +""" +$(TYPEDSIGNATURES) + +returns the `encoder_mat` field of the `Standardizer`. +""" get_encoder_mat(ss::Standardizer) = ss.encoder_mat + +""" +$(TYPEDSIGNATURES) + +returns the `decoder_mat` field of the `Standardizer`. +""" get_decoder_mat(ss::Standardizer) = ss.decoder_mat + +""" +$(TYPEDSIGNATURES) + +returns the `rank` field of the `Standardizer`. +""" get_rank(ss::Standardizer) = ss.rank function Base.show(io::IO, ss::Standardizer) @@ -281,26 +389,52 @@ function initialize_processor!(ss::Standardizer, data::MM) where {MM <: Abstract end end +""" +$(TYPEDSIGNATURES) + +Apply the `Standardizer` encoder, on a columns-are-data matrix +""" function encode_data(ss::Standardizer, data::MM) where {MM <: AbstractMatrix} data_mean = get_data_mean(ss) encoder_mat = get_encoder_mat(ss)[1] return encoder_mat * (data .- data_mean) end + +""" +$(TYPEDSIGNATURES) + +Apply the `Standardizer` decoder, on a columns-are-data matrix +""" function decode_data(ss::Standardizer, data::MM) where {MM <: AbstractMatrix} data_mean = get_data_mean(ss) decoder_mat = get_decoder_mat(ss)[1] return decoder_mat * data .+ data_mean end +""" +$(TYPEDSIGNATURES) + +Computes and populates the `data_mean` and `encoder_mat` and `decoder_mat` fields for the `Standardizer` +""" initialize_processor!(ss::Standardizer, data::MM, structure_matrix) where {MM <: AbstractMatrix} = initialize_processor!(ss, data) +""" +$(TYPEDSIGNATURES) + +Apply the `Standardizer` encoder to a provided structure matrix +""" function encode_structure_matrix(ss::Standardizer, structure_matrix::MM) where {MM <: AbstractMatrix} encoder_mat = get_encoder_mat(ss)[1] return encoder_mat * structure_matrix * encoder_mat' end +""" +$(TYPEDSIGNATURES) + +Apply the `Standardizer` decoder to a provided structure matrix +""" function decode_structure_matrix(ss::Standardizer, enc_structure_matrix::MM) where {MM <: AbstractMatrix} decoder_mat = get_decoder_mat(ss)[1] return decoder_mat * enc_structure_matrix * decoder_mat' @@ -312,10 +446,12 @@ $(TYPEDEF) Decorrelate the data via an SVD decomposition using a structure_matrix (`prior_cov` for inputs, `obs_noise_cov` for outputs), with optional truncation of singular vectors corresponding to largest singular values. The structure matrix will also become exactly `I` after processing. -Similar to the `Standardizer`, except that the `Standardizer` uses the estimated covariance of the data to decorrelate in-place of the structure matrix. There, the data samples will have sample mean `0` and covariance `I` after processing. +Similar to the [`Standardizer`](@ref), except that the `Standardizer` uses the estimated covariance of the data to decorrelate in-place of the structure matrix. There, the data samples will have sample mean `0` and covariance `I` after processing. In the `Decorrelater` the user can do a combined approach where one uses cov(data) + structure matrix for decorrelation with the `add_estimated_data_cov=true` +Preferred construction is with the [`decorrelate`](@ref) method + # Fields $(TYPEDFIELDS) """ @@ -332,13 +468,49 @@ struct Decorrelater{VV1, VV2, VV3, FT} <: DataContainerProcessor add_estimated_cov::Bool end +""" +$(TYPEDSIGNATURES) + +Constructs the `Decorrelater` struct. Users can add optional keyword arguments: +- `retain_var`[=1.0]: to project onto the leading singular vectors such that `retain_var` variance is retained +- `add_estimated_cov`[=false]: to add the estimated covariance to the structure matrix, and use this summation to create the new subspace. (false uses just the structure matrix) +""" decorrelate(; retain_var::FT = Float64(1.0), add_estimated_cov = false) where {FT} = Decorrelater([], [], [], min(max(retain_var, FT(0)), FT(1)), add_estimated_cov) +""" +$(TYPEDSIGNATURES) + +returns the `data_mean` field of the `Decorrelater`. +""" get_data_mean(dd::Decorrelater) = dd.data_mean + +""" +$(TYPEDSIGNATURES) + +returns the `encoder_mat` field of the `Decorrelater`. +""" get_encoder_mat(dd::Decorrelater) = dd.encoder_mat + +""" +$(TYPEDSIGNATURES) + +returns the `decoder_mat` field of the `Decorrelater`. +""" get_decoder_mat(dd::Decorrelater) = dd.decoder_mat + +""" +$(TYPEDSIGNATURES) + +returns the `retain_var` field of the `Decorrelater`. +""" get_retain_var(dd::Decorrelater) = dd.retain_var + +""" +$(TYPEDSIGNATURES) + +returns the `add_estimated_cov` field of the `Decorrelater`. +""" get_add_estimated_cov(dd::Decorrelater) = dd.add_estimated_cov function Base.show(io::IO, dd::Decorrelater) @@ -352,6 +524,11 @@ function Base.show(io::IO, dd::Decorrelater) print(io, out) end +""" +$(TYPEDSIGNATURES) + +Computes and populates the `data_mean` and `encoder_mat` and `decoder_mat` fields for the `Decorrelater` +""" function initialize_processor!( dd::Decorrelater, data::MM, @@ -389,23 +566,44 @@ function initialize_processor!( end end + +""" +$(TYPEDSIGNATURES) + +Apply the `Decorrelater` encoder, on a columns-are-data matrix +""" function encode_data(dd::Decorrelater, data::MM) where {MM <: AbstractMatrix} data_mean = get_data_mean(dd) encoder_mat = get_encoder_mat(dd)[1] return encoder_mat * (data .- data_mean) end +""" +$(TYPEDSIGNATURES) + +Apply the `Decorrelater` decoder, on a columns-are-data matrix +""" function decode_data(dd::Decorrelater, data::MM) where {MM <: AbstractMatrix} data_mean = get_data_mean(dd) decoder_mat = get_decoder_mat(dd)[1] return decoder_mat * data .+ data_mean end +""" +$(TYPEDSIGNATURES) + +Apply the `Decorrelater` encoder to a provided structure matrix +""" function encode_structure_matrix(dd::Decorrelater, structure_matrix::MM) where {MM <: AbstractMatrix} encoder_mat = get_encoder_mat(dd)[1] return encoder_mat * structure_matrix * encoder_mat' end +""" +$(TYPEDSIGNATURES) + +Apply the `Decorrelater` decoder to a provided structure matrix +""" function decode_structure_matrix(dd::Decorrelater, enc_structure_matrix::MM) where {MM <: AbstractMatrix} decoder_mat = get_decoder_mat(dd)[1] return decoder_mat * enc_structure_matrix * decoder_mat' @@ -424,6 +622,8 @@ $(TYPEDEF) Uses both input and output data to learn a subspace of maximal correlation between inputs and outputs. The subspace for a pair (X,Y) will be of size minimum(rank(X),rank(Y)), computed using SVD-based method e.g. See e.g., https://numerical.recipes/whp/notes/CanonCorrBySVD.pdf +Preferred construction is with the [`canonical_correlation`](@ref) method + # Fields $(TYPEDFIELDS) """ @@ -440,13 +640,48 @@ struct CanonicalCorrelation{VV1, VV2, VV3, FT, VV4} <: PairedDataContainerProces apply_to::VV4 end +""" +$(TYPEDSIGNATURES) + +Constructs the `CanonicalCorrelation` struct. Can optionally provide the keyword +- `retain_var`[=1.0]: to project onto the leading singular vectors (of the input-output product) such that `retain_var` variance is retained. +""" canonical_correlation(; retain_var = Float64(1.0)) = CanonicalCorrelation(Any[], Any[], Any[], retain_var, AbstractString[]) +""" +$(TYPEDSIGNATURES) + +returns the `data_mean` field of the `CanonicalCorrelation`. +""" get_data_mean(cc::CanonicalCorrelation) = cc.data_mean + +""" +$(TYPEDSIGNATURES) + +returns the `encoder_mat` field of the `CanonicalCorrelation`. +""" get_encoder_mat(cc::CanonicalCorrelation) = cc.encoder_mat + +""" +$(TYPEDSIGNATURES) + +returns the `decoder_mat` field of the `CanonicalCorrelation`. +""" get_decoder_mat(cc::CanonicalCorrelation) = cc.decoder_mat + +""" +$(TYPEDSIGNATURES) + +returns the `retain_var` field of the `CanonicalCorrelation`. +""" get_retain_var(cc::CanonicalCorrelation) = cc.retain_var + +""" +$(TYPEDSIGNATURES) + +returns the `apply_to` field of the `CanonicalCorrelation`. +""" get_apply_to(cc::CanonicalCorrelation) = cc.apply_to function Base.show(io::IO, cc::CanonicalCorrelation) @@ -458,6 +693,7 @@ function Base.show(io::IO, cc::CanonicalCorrelation) print(io, out) end + function initialize_processor!( cc::CanonicalCorrelation, in_data::MM, @@ -535,6 +771,11 @@ function initialize_processor!( end end +""" +$(TYPEDSIGNATURES) + +Computes and populates the `data_mean`, `encoder_mat`, `decoder_mat` and `apply_to` fields for the `CanonicalCorrelation` +""" initialize_processor!( cc::CanonicalCorrelation, in_data::MM, @@ -544,23 +785,43 @@ initialize_processor!( ) where {MM <: AbstractMatrix, AS <: AbstractString} = initialize_processor!(cc, in_data, out_data, apply_to) +""" +$(TYPEDSIGNATURES) + +Apply the `CanonicalCorrelation` encoder, on a columns-are-data matrix +""" function encode_data(cc::CanonicalCorrelation, data::MM) where {MM <: AbstractMatrix} data_mean = get_data_mean(cc) encoder_mat = get_encoder_mat(cc)[1] return encoder_mat * (data .- data_mean) end +""" +$(TYPEDSIGNATURES) + +Apply the `CanonicalCorrelation` decoder, on a columns-are-data matrix +""" function decode_data(cc::CanonicalCorrelation, data::MM) where {MM <: AbstractMatrix} data_mean = get_data_mean(cc) decoder_mat = get_decoder_mat(cc)[1] return decoder_mat * data .+ data_mean end +""" +$(TYPEDSIGNATURES) + +Apply the `CanonicalCorrelation` encoder to a provided structure matrix +""" function encode_structure_matrix(cc::CanonicalCorrelation, structure_matrix::MM) where {MM <: AbstractMatrix} encoder_mat = get_encoder_mat(cc)[1] return encoder_mat * structure_matrix * encoder_mat' end +""" +$(TYPEDSIGNATURES) + +Apply the `CanonicalCorrelation` decoder to a provided structure matrix +""" function decode_structure_matrix(cc::CanonicalCorrelation, enc_structure_matrix::MM) where {MM <: AbstractMatrix} decoder_mat = get_decoder_mat(cc)[1] return decoder_mat * enc_structure_matrix * decoder_mat' @@ -583,6 +844,11 @@ function initialize_and_encode_data!( return encode_data(dcp, data) end +""" +$(TYPEDSIGNATURES) + +Initializes the `DataContainerProcessor` encoder (often requires data, and structure matrices), then encodes the provided columns-are-data matrix +""" initialize_and_encode_data!( dcp::DCP, data::MM, @@ -595,12 +861,23 @@ initialize_and_encode_data!( AS <: AbstractString, } = initialize_and_encode_data!(dcp, data, structure_mat) + +""" +$(TYPEDSIGNATURES) + +decodes the columns-are-data matrix with the processor. +""" decode_data( dcp::DCP, data::MM, apply_to::AS, ) where {DCP <: DataContainerProcessor, MM <: AbstractMatrix, AS <: AbstractString} = decode_data(dcp, data) +""" +$(TYPEDSIGNATURES) + +Initializes the `PairedDataContainerProcesser` encoder (often requires input & output data, and structure matrices), then encodes either the input or output data (pair of columns-are-data matrices) based on `apply_to`. +""" function initialize_and_encode_data!( dcp::PDCP, data, @@ -616,6 +893,11 @@ function initialize_and_encode_data!( end end +""" +$(TYPEDSIGNATURES) + +decodes the input or output dat (pair of columns-are-data matrices) a with the processor, based on `apply_to`. +""" function decode_data(dcp::PDCP, data, apply_to::AS) where {PDCP <: PairedDataContainerProcessor, AS <: AbstractString} input_data, output_data = data if apply_to == "in" @@ -696,11 +978,18 @@ function create_encoder_schedule(schedule_in::VV) where {VV <: AbstractVector} return encoder_schedule end +""" +Create size-1 encoder schedule with a tuple of `(DataProcessor1(...), apply_to)` with `apply_to = "in"`, `"out"` or `"in_and_out"`. +""" create_encoder_schedule(schedule_in::TT) where {TT <: Tuple} = create_encoder_schedule([schedule_in]) # Functions to encode/decode with uninitialized schedule (require structure matrices as input) +""" +$TYPEDSIGNATURES +Takes in the created encoder schedule (See [`create_encoder_schedule`](@ref)), and initializes it, and encodes the paired data container, and structure matrices with it. +""" function encode_with_schedule( encoder_schedule::VV, io_pairs::PDC, @@ -738,7 +1027,11 @@ function encode_with_schedule( end # Functions to encode/decode with initialized schedule +""" +$TYPEDSIGNATURES +Takes in an already initialized encoder schedule, and encodes a `DataContainer`, the `in_or_out` string indicates if the data is input `"in"` or output `"out"` data (and thus encoded differently) +""" function encode_with_schedule( encoder_schedule::VV, data_container::DC, @@ -766,6 +1059,11 @@ function encode_with_schedule( return processed_container end +""" +$TYPEDSIGNATURES + +Takes in an already initialized encoder schedule, and encodes a structure matrix, the `in_or_out` string indicates if the structure matrix is for input `"in"` or output `"out"` space (and thus encoded differently) +""" function encode_with_schedule( encoder_schedule::VV, structure_matrix::USorM, @@ -791,6 +1089,11 @@ function encode_with_schedule( return processed_structure_matrix end +""" +$TYPEDSIGNATURES + +Takes in an already initialized encoder schedule, and decodes a `DataContainer`, and structure matrices with it, the `in_or_out` string indicates if the data is input `"in"` or output `"out"` data (and thus decoded differently) +""" function decode_with_schedule( encoder_schedule::VV, io_pairs::PDC, @@ -824,6 +1127,11 @@ function decode_with_schedule( return processed_io_pairs, processed_prior_cov, processed_obs_noise_cov end +""" +$TYPEDSIGNATURES + +Takes in an already initialized encoder schedule, and decodes a `DataContainer`, the `in_or_out` string indicates if the data is input `"in"` or output `"out"` data (and thus decoded differently) +""" function decode_with_schedule( encoder_schedule::VV, data_container::DC, @@ -852,6 +1160,11 @@ function decode_with_schedule( return processed_container end +""" +$TYPEDSIGNATURES + +Takes in an already initialized encoder schedule, and decodes a structure matrix, the `in_or_out` string indicates if the structure matrix is for input `"in"` or output `"out"` space (and thus decoded differently) +""" function decode_with_schedule( encoder_schedule::VV, structure_matrix::USorM, From ef4348662a271ec25fcba95496d407987047cc83 Mon Sep 17 00:00:00 2001 From: odunbar Date: Tue, 24 Jun 2025 14:10:35 -0700 Subject: [PATCH 16/25] removed standardization -> now only decorrelation --- src/Utilities.jl | 84 ++++++++++++++++++++++++++------------ test/Utilities/runtests.jl | 69 +++++++++++++------------------ 2 files changed, 86 insertions(+), 67 deletions(-) diff --git a/src/Utilities.jl b/src/Utilities.jl index b502de4af..23bd3602e 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -16,17 +16,17 @@ export orig2zscore export zscore2orig export PairedDataContainerProcessor, DataContainerProcessor -export UnivariateAffineScaling, AffineScaler, QuartileScaling, MinMaxScaling, ZScoreScaling, Standardizer, Decorrelater -export quartile_scale, minmax_scale, zscore_scale, standardize, decorrelate +export UnivariateAffineScaling, AffineScaler, QuartileScaling, MinMaxScaling, ZScoreScaling +export quartile_scale, minmax_scale, zscore_scale +export Decorrelater, decorrelate_sample_cov, decorrelate_structure_mat, decorrelate export get_type, get_shift, get_scale, get_data_mean, get_encoder_mat, get_decoder_mat, - get_rank, get_retain_var, - get_add_estimated_cov + get_decorrelate_with export create_encoder_schedule, encode_with_schedule, decode_with_schedule export initialize_processor!, initialize_and_encode_data!, encode_data, decode_data, encode_structure_matrix, decode_structure_matrix @@ -296,7 +296,7 @@ function decode_structure_matrix(as::AffineScaler, enc_structure_matrix::MM) whe end - +#= """ $(TYPEDEF) @@ -439,33 +439,41 @@ function decode_structure_matrix(ss::Standardizer, enc_structure_matrix::MM) whe decoder_mat = get_decoder_mat(ss)[1] return decoder_mat * enc_structure_matrix * decoder_mat' end - +=# """ $(TYPEDEF) -Decorrelate the data via an SVD decomposition using a structure_matrix (`prior_cov` for inputs, `obs_noise_cov` for outputs), with optional truncation of singular vectors corresponding to largest singular values. The structure matrix will also become exactly `I` after processing. +Decorrelate the data via tkaing an SVD decomposition and projecting onto the singular-vectors. -Similar to the [`Standardizer`](@ref), except that the `Standardizer` uses the estimated covariance of the data to decorrelate in-place of the structure matrix. There, the data samples will have sample mean `0` and covariance `I` after processing. +Preferred construction is with the methods +- [`decorrelate_structure_mat`](@ref) +- [`decorrelate_sample_cov`](@ref) +- [`decorrelate`](@ref) -In the `Decorrelater` the user can do a combined approach where one uses cov(data) + structure matrix for decorrelation with the `add_estimated_data_cov=true` +For `decorrelate_structure_mat`: +The SVD is taken over a structure matrix (e.g., `prior_cov` for inputs, `obs_noise_cov` for outputs). The structure matrix will become exactly `I` after processing. -Preferred construction is with the [`decorrelate`](@ref) method +For `decorrelate_sample_cov`: +The SVD is taken over the estimated covariance of the data. The data samples will have a `Normal(0,I)` distribution after processing. + +For `decorrelate(;decorrelate_with="combined")` (default): +The SVD is taken to be the sum of structure matrix and estimated covariance. This may be more robust to ill-specification of structure matrix, or poor estimation of the sample covariance. # Fields $(TYPEDFIELDS) """ -struct Decorrelater{VV1, VV2, VV3, FT} <: DataContainerProcessor +struct Decorrelater{VV1, VV2, VV3, FT, AS <: AbstractString} <: DataContainerProcessor "storage for the data mean" data_mean::VV1 - "the (possibly-encoded) structure matrix - provided by the user" + "the matrix used to perform encoding" encoder_mat::VV2 - "the inverse of the (possibly-encoded) structure matrix" + "the inverse of the the matrix used to perform encoding" decoder_mat::VV3 "the fraction of variance to be retained after truncating singular values (1 implies no truncation)" retain_var::FT - "If true, add the estimated cov(data) to the structure matrix for to modify the decorrelation" - add_estimated_cov::Bool + "Switch to choose what form of matrix to use to decorrelate the data" + decorrelate_with::AS end """ @@ -473,10 +481,31 @@ $(TYPEDSIGNATURES) Constructs the `Decorrelater` struct. Users can add optional keyword arguments: - `retain_var`[=1.0]: to project onto the leading singular vectors such that `retain_var` variance is retained -- `add_estimated_cov`[=false]: to add the estimated covariance to the structure matrix, and use this summation to create the new subspace. (false uses just the structure matrix) +- `decorrelate_with` [="structure_matrix"]: from which matrix do we provide subspace directions, options are + - "structure_mat", see [`decorrelate_structure_mat`] + - "sample_cov", see [`decorrelate_sample_cov`] + - "combined", sums the "sample_cov" and "structure_mat" matrices +""" +decorrelate(; retain_var::FT = Float64(1.0), decorrelate_with="combined") where {FT} = + Decorrelater([], [], [], min(max(retain_var, FT(0)), FT(1)), decorrelate_with) + +""" +$(TYPEDSIGNATURES) + +Constructs the `Decorrelater` struct, setting decorrelate_with = "sample_cov". Encoding data with this will ensure that the distribution of data samples after encoding will be `Normal(0,I)`. One can additionally add keywords: +- `retain_var`[=1.0]: to project onto the leading singular vectors such that `retain_var` variance is retained +""" +decorrelate_sample_cov(; retain_var::FT = Float64(1.0)) where {FT} = + Decorrelater([], [], [], min(max(retain_var, FT(0)), FT(1)), "sample_cov") + +""" +$(TYPEDSIGNATURES) + +Constructs the `Decorrelater` struct, setting decorrelate_with = "structure_mat". This encoding will transform a provided structure matrix into `I`. One can additionally add keywords: +- `retain_var`[=1.0]: to project onto the leading singular vectors such that `retain_var` variance is retained """ -decorrelate(; retain_var::FT = Float64(1.0), add_estimated_cov = false) where {FT} = - Decorrelater([], [], [], min(max(retain_var, FT(0)), FT(1)), add_estimated_cov) +decorrelate_structure_mat(; retain_var::FT = Float64(1.0)) where {FT} = + Decorrelater([], [], [], min(max(retain_var, FT(0)), FT(1)), "structure_mat") """ $(TYPEDSIGNATURES) @@ -509,18 +538,16 @@ get_retain_var(dd::Decorrelater) = dd.retain_var """ $(TYPEDSIGNATURES) -returns the `add_estimated_cov` field of the `Decorrelater`. +returns the `decorrelate_with` field of the `Decorrelater`. """ -get_add_estimated_cov(dd::Decorrelater) = dd.add_estimated_cov +get_decorrelate_with(dd::Decorrelater) = dd.decorrelate_with function Base.show(io::IO, dd::Decorrelater) out = "Decorrelater" if get_retain_var(dd) < 1.0 out *= ": retain_var=$(get_retain_var(dd)) " end - if get_add_estimated_cov(dd) - out *= ": add_estimated_cov=$(get_add_estimated_cov(dd)) " - end + out *= ": decorrelate_with=$(get_decorrelate_with(dd)) " print(io, out) end @@ -541,10 +568,15 @@ function initialize_processor!( if length(get_encoder_mat(dd)) == 0 # Can do tsvd here for large matrices - if get_add_estimated_cov(dd) - svdA = svd(structure_matrix + cov(data, dims = 2)) # with data covariance? + decorrelate_with = get_decorrelate_with(dd) + if decorrelate_with == "structure_mat" + svdA = svd(structure_matrix) + elseif decorrelate_with == "sample_cov" + svdA = svd(cov(data, dims = 2)) + elseif decorrelate_with == "combined" + svdA = svd(structure_matrix + cov(data, dims = 2)) else - svdA = svd(structure_matrix) + throw(ArgumentError("Keyword `decorrelate_with` must be taken from [\"sample_cov\", \"structure_mat\", \"combined\"]. Received $(decorrelate_with)")) end ret_var = get_retain_var(dd) if ret_var < 1.0 diff --git a/test/Utilities/runtests.jl b/test/Utilities/runtests.jl index 419aefb8c..f2a6337f5 100644 --- a/test/Utilities/runtests.jl +++ b/test/Utilities/runtests.jl @@ -91,27 +91,19 @@ end @test get_shift(QQ) == [1] @test get_scale(QQ) == [2] - ss = standardize() - @test isa(ss, Standardizer) - @test get_rank(ss) == typemax(Int) - ss2 = standardize(10) - @test get_rank(ss2) == 10 - SS = Standardizer(1, [1], [2], [3]) - @test get_data_mean(SS) == [1] - @test get_encoder_mat(SS) == [2] - @test get_decoder_mat(SS) == [3] - dd = decorrelate() @test get_retain_var(dd) == 1.0 - @test get_add_estimated_cov(dd) == false - dd2 = decorrelate(retain_var = 0.7, add_estimated_cov = true) + @test get_decorrelate_with(dd) == "combined" + dd2 = decorrelate_sample_cov(retain_var = 0.7) @test get_retain_var(dd2) == 0.7 - @test get_add_estimated_cov(dd2) == true - DD = Decorrelater([1], [2], [3], 1.0, false) + @test get_decorrelate_with(dd2) == "sample_cov" + dd3 = decorrelate_structure_mat(retain_var = 0.7) + @test get_retain_var(dd3) == 0.7 + @test get_decorrelate_with(dd3) == "structure_mat" + DD = Decorrelater([1], [2], [3], 1.0, "test") @test get_data_mean(DD) == [1] @test get_encoder_mat(DD) == [2] @test get_decoder_mat(DD) == [3] - # get some data as IO pairs for functional tests in_dim = 10 @@ -129,12 +121,11 @@ end "zscore", "quartile", "minmax", - "standardize", - "decorrelate", - "decorrelate-estimate-cov", + "decorrelate-sample-cov", + "decorrelate-structure-mat", + "decorrelate-combined", "canonical-correlation", - "standardize-truncate-to-5", - "decorrelate-retain-0.95-var", + "decorrelate-structure-mat-retain-0.95-var", "canonical-correlation-0.95-var", ] @@ -143,12 +134,11 @@ end (zscore_scale(), "in_and_out"), (quartile_scale(), "in_and_out"), (minmax_scale(), "in_and_out"), - (standardize(), "in_and_out"), - (decorrelate(), "in_and_out"), - (decorrelate(add_estimated_cov = true), "in_and_out"), + (decorrelate_sample_cov(), "in_and_out"), + (decorrelate_structure_mat(), "in_and_out"), + (decorrelate(), "in_and_out"), # combined (canonical_correlation(), "in_and_out"), - (standardize(5), "in_and_out"), - (decorrelate(retain_var = 0.95), "in_and_out"), + (decorrelate_structure_mat(retain_var = 0.95), "in_and_out"), (canonical_correlation(retain_var = 0.95), "in_and_out"), ] lossless = [fill(true, 6); fill(false, 4)] # are these lossy approximations? @@ -213,27 +203,24 @@ end pop_cov = cov(enc_dat, dims = 2) dimm = size(pop_cov, 1) big_tol = 0.1 - if name == "standardize" - @test all(isapprox.(pop_mean, zeros(dimm), atol = tol)) - @test isapprox(norm(pop_cov - I), 0.0, atol = tol * dimm^2) # expect very accurate - @test isapprox(norm(enc_covv - I), 0.0, atol = big_tol * dimm^2) # expect poorly accurate, particularly if dimm < dim - elseif name == "decorrelate" + if name == "decorrelate-structure-mat" @test all(isapprox.(pop_mean, zeros(dimm), atol = tol)) @test isapprox(norm(pop_cov - I), 0.0, atol = big_tol * dimm^2) # expect poorly accurate @test isapprox(norm(enc_covv - I), 0.0, atol = tol * dimm^2) # expect very accurate - elseif name == "decorrelate-estimate-cov" + + elseif name == "decorrelate-sample-cov" + @test all(isapprox.(pop_mean, zeros(dimm), atol = tol)) + @test isapprox(norm(pop_cov - I), 0.0, atol = tol * dimm^2) # expect very accurate + @test isapprox(norm(enc_covv - I), 0.0, atol = big_tol * dimm^2) # expect poorly accurate, particularly if dimm < dim + + elseif name == "decorrelate-combined" @test all(isapprox.(pop_mean, zeros(dimm), atol = tol)) @test isapprox(norm(pop_cov - I), 0.0, atol = big_tol * dimm^2) # expect poorly accurate @test isapprox(norm(enc_covv - I), 0.0, atol = big_tol * dimm^2) # expect poorly accurate end # Multivariate lossy dim-reduction tests - if name == "standardize-truncate-to-5" - @test dimm == 5 - @test all(isapprox.(pop_mean, zeros(dimm), atol = tol)) - @test isapprox(norm(pop_cov - I), 0.0, atol = tol * dimm^2) # expect very accurate - @test isapprox(norm(enc_covv - I), 0.0, atol = big_tol * dimm^2) # expect poorly accurate - elseif name == "decorrelate-retain-0.95-var" + if name == "decorrelate-structure-mat-retain-0.95-var" svdc = svd(test_covv) var_cumsum = cumsum(svdc.S .^ 2) ./ sum(svdc.S .^ 2) @test var_cumsum[dimm] > 0.95 @@ -266,8 +253,8 @@ end # test decode approximation of lossless options if ll_flag # when dimm < dim, loss can occur in some tests - tol1 = (name == "decorrelate" && dimm < dim) ? big_tol : tol - tol2 = (name == "standardize" && dimm < dim) ? big_tol : tol + tol1 = (name == "decorrelate-structure-mat" && dimm < dim) ? big_tol : tol + tol2 = (name == "decorrelate-sample-cov" && dimm < dim) ? big_tol : tol @test isapprox(norm(dec_dat - test_dat), 0.0, atol = tol1 * dim * samples) @test isapprox(norm(dec_covv - test_covv), 0.0, atol = tol2 * dim^2) @@ -286,9 +273,9 @@ end schedule_builder = [ (zscore_scale(), "in_and_out"), # (quartile_scale(), "in"), - (standardize(), "in_and_out"), + (decorrelate_sample_cov(), "in_and_out"), (minmax_scale(), "out"), - (decorrelate(), "in_and_out"), + (decorrelate_structure_mat(), "in_and_out"), (canonical_correlation(), "in"), ] From c795600be463bbcc65feb7b2eae568cc1d35f41c Mon Sep 17 00:00:00 2001 From: odunbar Date: Tue, 24 Jun 2025 14:11:36 -0700 Subject: [PATCH 17/25] removed commented code --- src/Utilities.jl | 147 +---------------------------------------------- 1 file changed, 1 insertion(+), 146 deletions(-) diff --git a/src/Utilities.jl b/src/Utilities.jl index 23bd3602e..8583579df 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -167,7 +167,7 @@ $(TYPEDSIGNATURES) Constructs `AffineScaler{ZScoreScaling}` processor. As part of an encoder schedule, this will apply the transform ``\\frac{x-\\mu}{\\sigma}``, (where ``x\\sim N(\\mu,\\sigma)``), to each data dimension. -For multivariate standardization, see `Standardizer` +For multivariate standardization, see [`Decorrelater`](@ref) """ zscore_scale() = AffineScaler(ZScoreScaling) @@ -296,151 +296,6 @@ function decode_structure_matrix(as::AffineScaler, enc_structure_matrix::MM) whe end -#= -""" -$(TYPEDEF) - -Standardizes the data to a multivariate N(0,I) distribution via `(xcov)^{-1/2}*(x-x_mean)`, along with rank reduction if data is low rank, or if the user provides a rank. This guarantees that the data samples will have sample mean `0` and covariance `I` after processing. - -Similar to the [`Decorrelater`](@ref), but there the structure matrices will directly be used to perform decorrelation, and will become `I` after processing. In that case, the data samples may not have covariance `I` after processing. - -Preferred construction is using the [`standardize`](@ref) method. - -# Fields -$(TYPEDFIELDS) -""" -struct Standardizer{VV1 <: AbstractVector, VV2 <: AbstractVector, VV3 <: AbstractVector} <: DataContainerProcessor - "user-provided rank for additional truncation of the space. used if rank < rank(data_cov)" - rank::Int - "storage for the data mean" - data_mean::VV1 - "storage for wide-matrix representing `data_cov^{-1/2}`" - encoder_mat::VV2 - "storage for tall-matrix representing `data_cov^{1/2}`" - decoder_mat::VV3 -end - -""" -$(TYPEDSIGNATURES) - -Constructs the [`Standardizer`](@ref) struct with no truncation of rank -""" -standardize() = Standardizer(typemax(Int), Float64[], Any[], Any[]) -""" -$(TYPEDSIGNATURES) - -Constructs the `Standardizer` struct with a provided truncation to rank `rk`. -""" -standardize(rk::Int) = Standardizer(rk, Float64[], Any[], Any[]) - - -""" -$(TYPEDSIGNATURES) - -returns the `data_mean` field of the `Standardizer`. -""" -get_data_mean(ss::Standardizer) = ss.data_mean - -""" -$(TYPEDSIGNATURES) - -returns the `encoder_mat` field of the `Standardizer`. -""" -get_encoder_mat(ss::Standardizer) = ss.encoder_mat - -""" -$(TYPEDSIGNATURES) - -returns the `decoder_mat` field of the `Standardizer`. -""" -get_decoder_mat(ss::Standardizer) = ss.decoder_mat - -""" -$(TYPEDSIGNATURES) - -returns the `rank` field of the `Standardizer`. -""" -get_rank(ss::Standardizer) = ss.rank - -function Base.show(io::IO, ss::Standardizer) - if get_rank(ss) < typemax(Int) - out = "Standardizer: user-rank=$(get_rank(ss))" - else - out = "Standardizer" - end - print(io, out) -end - -function initialize_processor!(ss::Standardizer, data::MM) where {MM <: AbstractMatrix} - if length(get_data_mean(ss)) == 0 - append!(get_data_mean(ss), mean(data, dims = 2)) - end - if length(get_encoder_mat(ss)) == 0 - # can use tsvd here so we only compute up to rk - cov_data = cov(data, dims = 2) - svdc = svd(cov_data) - rk = min(rank(cov_data), max(get_rank(ss), 1)) - - # rescale with sqrt(size(data,2)) to normalize - encoder_mat = Diagonal(1 ./ sqrt.(svdc.S[1:rk])) * svdc.Vt[1:rk, :] - decoder_mat = svdc.Vt[1:rk, :]' * Diagonal(sqrt.(svdc.S[1:rk])) - push!(get_encoder_mat(ss), encoder_mat) - push!(get_decoder_mat(ss), decoder_mat) - end -end - -""" -$(TYPEDSIGNATURES) - -Apply the `Standardizer` encoder, on a columns-are-data matrix -""" -function encode_data(ss::Standardizer, data::MM) where {MM <: AbstractMatrix} - data_mean = get_data_mean(ss) - encoder_mat = get_encoder_mat(ss)[1] - return encoder_mat * (data .- data_mean) -end - - -""" -$(TYPEDSIGNATURES) - -Apply the `Standardizer` decoder, on a columns-are-data matrix -""" -function decode_data(ss::Standardizer, data::MM) where {MM <: AbstractMatrix} - data_mean = get_data_mean(ss) - decoder_mat = get_decoder_mat(ss)[1] - return decoder_mat * data .+ data_mean -end - -""" -$(TYPEDSIGNATURES) - -Computes and populates the `data_mean` and `encoder_mat` and `decoder_mat` fields for the `Standardizer` -""" -initialize_processor!(ss::Standardizer, data::MM, structure_matrix) where {MM <: AbstractMatrix} = - initialize_processor!(ss, data) - -""" -$(TYPEDSIGNATURES) - -Apply the `Standardizer` encoder to a provided structure matrix -""" -function encode_structure_matrix(ss::Standardizer, structure_matrix::MM) where {MM <: AbstractMatrix} - encoder_mat = get_encoder_mat(ss)[1] - return encoder_mat * structure_matrix * encoder_mat' -end - -""" -$(TYPEDSIGNATURES) - -Apply the `Standardizer` decoder to a provided structure matrix -""" -function decode_structure_matrix(ss::Standardizer, enc_structure_matrix::MM) where {MM <: AbstractMatrix} - decoder_mat = get_decoder_mat(ss)[1] - return decoder_mat * enc_structure_matrix * decoder_mat' -end -=# - """ $(TYPEDEF) From 9d6e5738b2e8a6abb0bf3e29721e53f654637692 Mon Sep 17 00:00:00 2001 From: odunbar Date: Tue, 24 Jun 2025 14:17:04 -0700 Subject: [PATCH 18/25] prior->input, obs_noise_cov -> output structure mats --- src/Utilities.jl | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/src/Utilities.jl b/src/Utilities.jl index 8583579df..67dc91271 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -880,8 +880,8 @@ Takes in the created encoder schedule (See [`create_encoder_schedule`](@ref)), a function encode_with_schedule( encoder_schedule::VV, io_pairs::PDC, - prior_cov_in::USorM1, - obs_noise_cov_in::USorM2, + input_structure_mat::USorM1, + output_structure_mat::USorM2, ) where { VV <: AbstractVector, PDC <: PairedDataContainer, @@ -889,28 +889,28 @@ function encode_with_schedule( USorM2 <: Union{UniformScaling, AbstractMatrix}, } processed_io_pairs = deepcopy(io_pairs) - processed_prior_cov = deepcopy(prior_cov_in) - processed_obs_noise_cov = deepcopy(obs_noise_cov_in) + processed_input_structure_mat = deepcopy(input_structure_mat) + processed_output_structure_mat = deepcopy(output_structure_mat) # apply_to is the string "in", "out" etc. for (processor, extract_data, apply_to) in encoder_schedule @info "Initialize encoding of data: \"$(apply_to)\" with $(processor)" if apply_to == "in" - structure_matrix = processed_prior_cov + structure_matrix = processed_input_structure_mat elseif apply_to == "out" - structure_matrix = processed_obs_noise_cov + structure_matrix = processed_output_structure_mat end processed = initialize_and_encode_data!(processor, extract_data(processed_io_pairs), structure_matrix, apply_to) if apply_to == "in" - processed_prior_cov = encode_structure_matrix(processor, structure_matrix) + processed_input_structure_mat = encode_structure_matrix(processor, structure_matrix) processed_io_pairs = PairedDataContainer(processed, get_outputs(processed_io_pairs)) elseif apply_to == "out" - processed_obs_noise_cov = encode_structure_matrix(processor, structure_matrix) + processed_output_structure_mat = encode_structure_matrix(processor, structure_matrix) processed_io_pairs = PairedDataContainer(get_inputs(processed_io_pairs), processed) end end - return processed_io_pairs, processed_prior_cov, processed_obs_noise_cov + return processed_io_pairs, processed_input_structure_mat, processed_output_structure_mat end # Functions to encode/decode with initialized schedule @@ -984,8 +984,8 @@ Takes in an already initialized encoder schedule, and decodes a `DataContainer`, function decode_with_schedule( encoder_schedule::VV, io_pairs::PDC, - prior_cov_in::USorM1, - obs_noise_cov_in::USorM2, + input_structure_mat::USorM1, + output_structure_mat::USorM2, ) where { VV <: AbstractVector, USorM1 <: Union{UniformScaling, AbstractMatrix}, @@ -994,8 +994,8 @@ function decode_with_schedule( } processed_io_pairs = deepcopy(io_pairs) - processed_prior_cov = deepcopy(prior_cov_in) - processed_obs_noise_cov = deepcopy(obs_noise_cov_in) + processed_input_structure_mat = deepcopy(input_structure_mat) + processed_output_structure_mat = deepcopy(output_structure_mat) # apply_to is the string "in", "out" etc. for idx in reverse(eachindex(encoder_schedule)) @@ -1003,15 +1003,15 @@ function decode_with_schedule( processed = decode_data(processor, extract_data(processed_io_pairs), apply_to) if apply_to == "in" - processed_prior_cov = decode_structure_matrix(processor, processed_prior_cov) + processed_input_structure_mat = decode_structure_matrix(processor, processed_input_structure_mat) processed_io_pairs = PairedDataContainer(processed, get_outputs(processed_io_pairs)) elseif apply_to == "out" - processed_obs_noise_cov = decode_structure_matrix(processor, processed_obs_noise_cov) + processed_output_structure_mat = decode_structure_matrix(processor, processed_output_structure_mat) processed_io_pairs = PairedDataContainer(get_inputs(processed_io_pairs), processed) end end - return processed_io_pairs, processed_prior_cov, processed_obs_noise_cov + return processed_io_pairs, processed_input_structure_mat, processed_output_structure_mat end """ From 4ce4af5a0605e6eae024a421f846a9200160211b Mon Sep 17 00:00:00 2001 From: odunbar Date: Tue, 24 Jun 2025 14:23:29 -0700 Subject: [PATCH 19/25] docstring typos --- src/Utilities.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/Utilities.jl b/src/Utilities.jl index 67dc91271..755e60bd7 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -335,11 +335,11 @@ end $(TYPEDSIGNATURES) Constructs the `Decorrelater` struct. Users can add optional keyword arguments: -- `retain_var`[=1.0]: to project onto the leading singular vectors such that `retain_var` variance is retained -- `decorrelate_with` [="structure_matrix"]: from which matrix do we provide subspace directions, options are - - "structure_mat", see [`decorrelate_structure_mat`] - - "sample_cov", see [`decorrelate_sample_cov`] - - "combined", sums the "sample_cov" and "structure_mat" matrices +- `retain_var`[=`1.0`]: to project onto the leading singular vectors such that `retain_var` variance is retained +- `decorrelate_with` [=`"structure_matrix"`]: from which matrix do we provide subspace directions, options are + - `"structure_mat"`, see [`decorrelate_structure_mat`](@ref) + - `"sample_cov"`, see [`decorrelate_sample_cov`](@ref) + - `"combined"`, sums the `"sample_cov"` and `"structure_mat"` matrices """ decorrelate(; retain_var::FT = Float64(1.0), decorrelate_with="combined") where {FT} = Decorrelater([], [], [], min(max(retain_var, FT(0)), FT(1)), decorrelate_with) @@ -348,7 +348,7 @@ decorrelate(; retain_var::FT = Float64(1.0), decorrelate_with="combined") where $(TYPEDSIGNATURES) Constructs the `Decorrelater` struct, setting decorrelate_with = "sample_cov". Encoding data with this will ensure that the distribution of data samples after encoding will be `Normal(0,I)`. One can additionally add keywords: -- `retain_var`[=1.0]: to project onto the leading singular vectors such that `retain_var` variance is retained +- `retain_var`[=`1.0`]: to project onto the leading singular vectors such that `retain_var` variance is retained """ decorrelate_sample_cov(; retain_var::FT = Float64(1.0)) where {FT} = Decorrelater([], [], [], min(max(retain_var, FT(0)), FT(1)), "sample_cov") @@ -357,7 +357,7 @@ decorrelate_sample_cov(; retain_var::FT = Float64(1.0)) where {FT} = $(TYPEDSIGNATURES) Constructs the `Decorrelater` struct, setting decorrelate_with = "structure_mat". This encoding will transform a provided structure matrix into `I`. One can additionally add keywords: -- `retain_var`[=1.0]: to project onto the leading singular vectors such that `retain_var` variance is retained +- `retain_var`[=`1.0`]: to project onto the leading singular vectors such that `retain_var` variance is retained """ decorrelate_structure_mat(; retain_var::FT = Float64(1.0)) where {FT} = Decorrelater([], [], [], min(max(retain_var, FT(0)), FT(1)), "structure_mat") From 41e6ecf63069684f3d539496561ac77634a3d627 Mon Sep 17 00:00:00 2001 From: odunbar Date: Tue, 24 Jun 2025 15:30:35 -0700 Subject: [PATCH 20/25] format --- src/Utilities.jl | 22 ++++++++++------------ test/Utilities/runtests.jl | 4 ++-- 2 files changed, 12 insertions(+), 14 deletions(-) diff --git a/src/Utilities.jl b/src/Utilities.jl index 755e60bd7..007ea6b39 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -20,13 +20,7 @@ export UnivariateAffineScaling, AffineScaler, QuartileScaling, MinMaxScaling, ZS export quartile_scale, minmax_scale, zscore_scale export Decorrelater, decorrelate_sample_cov, decorrelate_structure_mat, decorrelate export get_type, - get_shift, - get_scale, - get_data_mean, - get_encoder_mat, - get_decoder_mat, - get_retain_var, - get_decorrelate_with + get_shift, get_scale, get_data_mean, get_encoder_mat, get_decoder_mat, get_retain_var, get_decorrelate_with export create_encoder_schedule, encode_with_schedule, decode_with_schedule export initialize_processor!, initialize_and_encode_data!, encode_data, decode_data, encode_structure_matrix, decode_structure_matrix @@ -341,8 +335,8 @@ Constructs the `Decorrelater` struct. Users can add optional keyword arguments: - `"sample_cov"`, see [`decorrelate_sample_cov`](@ref) - `"combined"`, sums the `"sample_cov"` and `"structure_mat"` matrices """ -decorrelate(; retain_var::FT = Float64(1.0), decorrelate_with="combined") where {FT} = - Decorrelater([], [], [], min(max(retain_var, FT(0)), FT(1)), decorrelate_with) +decorrelate(; retain_var::FT = Float64(1.0), decorrelate_with = "combined") where {FT} = + Decorrelater([], [], [], min(max(retain_var, FT(0)), FT(1)), decorrelate_with) """ $(TYPEDSIGNATURES) @@ -425,13 +419,17 @@ function initialize_processor!( # Can do tsvd here for large matrices decorrelate_with = get_decorrelate_with(dd) if decorrelate_with == "structure_mat" - svdA = svd(structure_matrix) + svdA = svd(structure_matrix) elseif decorrelate_with == "sample_cov" - svdA = svd(cov(data, dims = 2)) + svdA = svd(cov(data, dims = 2)) elseif decorrelate_with == "combined" svdA = svd(structure_matrix + cov(data, dims = 2)) else - throw(ArgumentError("Keyword `decorrelate_with` must be taken from [\"sample_cov\", \"structure_mat\", \"combined\"]. Received $(decorrelate_with)")) + throw( + ArgumentError( + "Keyword `decorrelate_with` must be taken from [\"sample_cov\", \"structure_mat\", \"combined\"]. Received $(decorrelate_with)", + ), + ) end ret_var = get_retain_var(dd) if ret_var < 1.0 diff --git a/test/Utilities/runtests.jl b/test/Utilities/runtests.jl index f2a6337f5..6e771afc2 100644 --- a/test/Utilities/runtests.jl +++ b/test/Utilities/runtests.jl @@ -208,12 +208,12 @@ end @test isapprox(norm(pop_cov - I), 0.0, atol = big_tol * dimm^2) # expect poorly accurate @test isapprox(norm(enc_covv - I), 0.0, atol = tol * dimm^2) # expect very accurate - elseif name == "decorrelate-sample-cov" + elseif name == "decorrelate-sample-cov" @test all(isapprox.(pop_mean, zeros(dimm), atol = tol)) @test isapprox(norm(pop_cov - I), 0.0, atol = tol * dimm^2) # expect very accurate @test isapprox(norm(enc_covv - I), 0.0, atol = big_tol * dimm^2) # expect poorly accurate, particularly if dimm < dim - elseif name == "decorrelate-combined" + elseif name == "decorrelate-combined" @test all(isapprox.(pop_mean, zeros(dimm), atol = tol)) @test isapprox(norm(pop_cov - I), 0.0, atol = big_tol * dimm^2) # expect poorly accurate @test isapprox(norm(enc_covv - I), 0.0, atol = big_tol * dimm^2) # expect poorly accurate From bba4cfbfd91ec9bd8ef3f7b632d00792f75c37d6 Mon Sep 17 00:00:00 2001 From: odunbar Date: Wed, 25 Jun 2025 10:14:57 -0700 Subject: [PATCH 21/25] symmetrise the affine transformation --- src/Utilities.jl | 8 ++++---- test/Utilities/runtests.jl | 6 +++--- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/Utilities.jl b/src/Utilities.jl index 007ea6b39..631840d1b 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -173,7 +173,7 @@ $(TYPEDSIGNATURES) Gets the UnivariateAffineScaling type `T` """ -get_type(as::AffineScaler{T, VV}) where {T, VV} = T +get_type(as::AffineScaler{T}) where {T} = T """ $(TYPEDSIGNATURES) @@ -277,7 +277,7 @@ $(TYPEDSIGNATURES) Apply the `AffineScaler` encoder to a provided structure matrix """ function encode_structure_matrix(as::AffineScaler, structure_matrix::MM) where {MM <: AbstractMatrix} - return Diagonal(1 ./ get_scale(as) .^ 2) * structure_matrix + return Diagonal(1 ./ get_scale(as)) * structure_matrix * Diagonal(1 ./ get_scale(as)) end """ @@ -286,14 +286,14 @@ $(TYPEDSIGNATURES) Apply the `AffineScaler` decoder to a provided structure matrix """ function decode_structure_matrix(as::AffineScaler, enc_structure_matrix::MM) where {MM <: AbstractMatrix} - return Diagonal(get_scale(as) .^ 2) * enc_structure_matrix + return Diagonal(get_scale(as)) * enc_structure_matrix * Diagonal(get_scale(as)) end """ $(TYPEDEF) -Decorrelate the data via tkaing an SVD decomposition and projecting onto the singular-vectors. +Decorrelate the data via taking an SVD decomposition and projecting onto the singular-vectors. Preferred construction is with the methods - [`decorrelate_structure_mat`](@ref) diff --git a/test/Utilities/runtests.jl b/test/Utilities/runtests.jl index 6e771afc2..c48ebc4ef 100644 --- a/test/Utilities/runtests.jl +++ b/test/Utilities/runtests.jl @@ -171,7 +171,7 @@ end test_vec = [[mean(dd), std(dd)] for dd in eachrow(test_dat)] test_mat = reduce(hcat, test_vec) - @test isapprox(norm(enc_covv - Diagonal(1 ./ test_mat[2, :] .^ 2) * test_covv), 0.0, atol = tol * dim^2) + @test isapprox(norm(enc_covv - Diagonal(1 ./ test_mat[2, :]) * test_covv * Diagonal(1 ./ test_mat[2, :] )), 0.0, atol = tol * dim^2) elseif name == "quartile" quartiles_vec = [quantile(dd, [0.25, 0.5, 0.75]) for dd in eachrow(enc_dat)] quartiles_mat = reduce(hcat, quartiles_vec) # 3 rows: Q1, Q2, and Q3 @@ -180,7 +180,7 @@ end test_vec = [quantile(dd, [0.25, 0.5, 0.75]) for dd in eachrow(test_dat)] test_mat = reduce(hcat, test_vec) @test isapprox( - norm(enc_covv - Diagonal(1 ./ (test_mat[3, :] - test_mat[1, :]) .^ 2) * test_covv), + norm(enc_covv - Diagonal(1 ./ (test_mat[3, :] - test_mat[1, :])) * test_covv * Diagonal(1 ./ (test_mat[3, :] - test_mat[1, :]))), 0.0, atol = tol * dim^2, ) @@ -192,7 +192,7 @@ end test_vec = [[minimum(dd), maximum(dd)] for dd in eachrow(test_dat)] test_mat = reduce(hcat, test_vec) @test isapprox( - norm(enc_covv - Diagonal(1 ./ (test_mat[2, :] - test_mat[1, :]) .^ 2) * test_covv), + norm(enc_covv - Diagonal(1 ./ (test_mat[2, :] - test_mat[1, :])) * test_covv * Diagonal(1 ./ (test_mat[2, :] - test_mat[1, :]))), 0.0, atol = tol * dim^2, ) From 2db4e383b3b707ae23e9edac3060c80f61114f01 Mon Sep 17 00:00:00 2001 From: odunbar Date: Wed, 25 Jun 2025 10:32:32 -0700 Subject: [PATCH 22/25] rank varies with choice of matrix. --- src/Utilities.jl | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/src/Utilities.jl b/src/Utilities.jl index 631840d1b..22c07650e 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -304,7 +304,7 @@ For `decorrelate_structure_mat`: The SVD is taken over a structure matrix (e.g., `prior_cov` for inputs, `obs_noise_cov` for outputs). The structure matrix will become exactly `I` after processing. For `decorrelate_sample_cov`: -The SVD is taken over the estimated covariance of the data. The data samples will have a `Normal(0,I)` distribution after processing. +The SVD is taken over the estimated covariance of the data. The data samples will have a `Normal(0,I)` distribution after processing, For `decorrelate(;decorrelate_with="combined")` (default): The SVD is taken to be the sum of structure matrix and estimated covariance. This may be more robust to ill-specification of structure matrix, or poor estimation of the sample covariance. @@ -420,10 +420,15 @@ function initialize_processor!( decorrelate_with = get_decorrelate_with(dd) if decorrelate_with == "structure_mat" svdA = svd(structure_matrix) + rk = rank(structure_matrix) elseif decorrelate_with == "sample_cov" - svdA = svd(cov(data, dims = 2)) + cd = cov(data, dims = 2) + svdA = svd(cd) + rk = rank(cd) elseif decorrelate_with == "combined" - svdA = svd(structure_matrix + cov(data, dims = 2)) + spluscd = structure_matrix + cov(data, dims = 2) + svdA = svd(spluscd) + rk = rank(spluscd) else throw( ArgumentError( @@ -437,7 +442,7 @@ function initialize_processor!( trunc = minimum(findall(x -> (x > ret_var), sv_cumsum)) @info " truncating at $(trunc)/$(length(sv_cumsum)) retaining $(100.0*sv_cumsum[trunc])% of the variance of the structure matrix" else - trunc = rank(structure_matrix) + trunc = rk end sqrt_inv_sv = Diagonal(1.0 ./ sqrt.(svdA.S[1:trunc])) sqrt_sv = Diagonal(sqrt.(svdA.S[1:trunc])) From 939ecb8e0419bf5e3385abda850a8f8210cdb8c0 Mon Sep 17 00:00:00 2001 From: odunbar Date: Wed, 25 Jun 2025 11:00:35 -0700 Subject: [PATCH 23/25] neaten code based on comments --- src/Utilities.jl | 31 ++++++++++++++++--------------- test/Utilities/runtests.jl | 6 +++--- 2 files changed, 19 insertions(+), 18 deletions(-) diff --git a/src/Utilities.jl b/src/Utilities.jl index 22c07650e..be1653f11 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -330,13 +330,13 @@ $(TYPEDSIGNATURES) Constructs the `Decorrelater` struct. Users can add optional keyword arguments: - `retain_var`[=`1.0`]: to project onto the leading singular vectors such that `retain_var` variance is retained -- `decorrelate_with` [=`"structure_matrix"`]: from which matrix do we provide subspace directions, options are +- `decorrelate_with` [=`"combined"`]: from which matrix do we provide subspace directions, options are - `"structure_mat"`, see [`decorrelate_structure_mat`](@ref) - `"sample_cov"`, see [`decorrelate_sample_cov`](@ref) - `"combined"`, sums the `"sample_cov"` and `"structure_mat"` matrices """ decorrelate(; retain_var::FT = Float64(1.0), decorrelate_with = "combined") where {FT} = - Decorrelater([], [], [], min(max(retain_var, FT(0)), FT(1)), decorrelate_with) + Decorrelater([], [], [], clamp(retain_var, FT(0), FT(1)), decorrelate_with) """ $(TYPEDSIGNATURES) @@ -345,7 +345,7 @@ Constructs the `Decorrelater` struct, setting decorrelate_with = "sample_cov". E - `retain_var`[=`1.0`]: to project onto the leading singular vectors such that `retain_var` variance is retained """ decorrelate_sample_cov(; retain_var::FT = Float64(1.0)) where {FT} = - Decorrelater([], [], [], min(max(retain_var, FT(0)), FT(1)), "sample_cov") + Decorrelater([], [], [], clamp(retain_var, FT(0), FT(1)), "sample_cov") """ $(TYPEDSIGNATURES) @@ -354,7 +354,7 @@ Constructs the `Decorrelater` struct, setting decorrelate_with = "structure_mat" - `retain_var`[=`1.0`]: to project onto the leading singular vectors such that `retain_var` variance is retained """ decorrelate_structure_mat(; retain_var::FT = Float64(1.0)) where {FT} = - Decorrelater([], [], [], min(max(retain_var, FT(0)), FT(1)), "structure_mat") + Decorrelater([], [], [], clamp(retain_var, FT(0), FT(1)), "structure_mat") """ $(TYPEDSIGNATURES) @@ -393,10 +393,10 @@ get_decorrelate_with(dd::Decorrelater) = dd.decorrelate_with function Base.show(io::IO, dd::Decorrelater) out = "Decorrelater" + out *= ": decorrelate_with=$(get_decorrelate_with(dd))" if get_retain_var(dd) < 1.0 - out *= ": retain_var=$(get_retain_var(dd)) " + out *= ", retain_var=$(get_retain_var(dd))" end - out *= ": decorrelate_with=$(get_decorrelate_with(dd)) " print(io, out) end @@ -536,8 +536,8 @@ $(TYPEDSIGNATURES) Constructs the `CanonicalCorrelation` struct. Can optionally provide the keyword - `retain_var`[=1.0]: to project onto the leading singular vectors (of the input-output product) such that `retain_var` variance is retained. """ -canonical_correlation(; retain_var = Float64(1.0)) = - CanonicalCorrelation(Any[], Any[], Any[], retain_var, AbstractString[]) +canonical_correlation(; retain_var::FT = Float64(1.0)) where {FT} = + CanonicalCorrelation(Any[], Any[], Any[], clamp(retain_var, FT(0), FT(1)) , AbstractString[]) """ $(TYPEDSIGNATURES) @@ -575,10 +575,13 @@ returns the `apply_to` field of the `CanonicalCorrelation`. get_apply_to(cc::CanonicalCorrelation) = cc.apply_to function Base.show(io::IO, cc::CanonicalCorrelation) + + out = "CanonicalCorrelation:" + if length(get_apply_to(cc))>0 + out *=" apply_to=$(get_apply_to(cc)[1])" + end if get_retain_var(cc) < 1.0 - out = "CanonicalCorrelation: retain_var=$(get_retain_var(cc))" - else - out = "CanonicalCorrelation" + out *= " retain_var=$(get_retain_var(cc))" end print(io, out) end @@ -614,12 +617,10 @@ function initialize_processor!( end # Individually decompose in and out - # Want to use the nonsquare singular vector svdi = svd(in_data .- mean(in_data, dims = 2)) svdo = svd(out_data .- mean(out_data, dims = 2)) - # determine regime - # ensure we non-square sv (in_mat = (in_dim x n_samples), out_mat = (out_dim x n_samples)) + # ensure correct shaping (in_mat = (in_dim x n_samples), out_mat = (out_dim x n_samples)) in_mat_sq, in_mat_nonsq = (size(svdi.U, 1) == size(svdi.U, 2)) ? (svdi.U, svdi.Vt) : (svdi.Vt, svdi.U) out_mat_sq, out_mat_nonsq = (size(svdo.U, 1) == size(svdo.U, 2)) ? (svdo.U, svdo.Vt) : (svdo.Vt, svdo.U) @@ -928,7 +929,7 @@ function encode_with_schedule( in_or_out::AS, ) where {VV <: AbstractVector, DC <: DataContainer, AS <: AbstractString} - if !(in_or_out ∈ ["in", "out"]) + if in_or_out ∉ ["in", "out"] throw( ArgumentError( "`in_or_out` must be either \"in\" (data is an input) or \"out\" (data is an output). Received $(in_or_out)", diff --git a/test/Utilities/runtests.jl b/test/Utilities/runtests.jl index c48ebc4ef..c83ecc278 100644 --- a/test/Utilities/runtests.jl +++ b/test/Utilities/runtests.jl @@ -108,7 +108,7 @@ end in_dim = 10 out_dim = 50 - samples = 120 # for full test coverage have samples in_dim < samples < out_dim + samples = 120 x = randn(rng, in_dim, in_dim) prior_cov = x * x' @@ -117,7 +117,7 @@ end out_data = rand(rng, MvNormal(-10 * ones(out_dim), obs_noise_cov), samples) io_pairs = PairedDataContainer(in_data, out_data) - test_names = [ # order as in tests + test_names = [ # order as in schedules below "zscore", "quartile", "minmax", @@ -271,7 +271,7 @@ end io_pairs = PairedDataContainer(in_data, out_data) schedule_builder = [ - (zscore_scale(), "in_and_out"), # + (zscore_scale(), "in_and_out"), (quartile_scale(), "in"), (decorrelate_sample_cov(), "in_and_out"), (minmax_scale(), "out"), From 71c2a6993af531a59270b8f830b8403952fe6465 Mon Sep 17 00:00:00 2001 From: odunbar Date: Wed, 25 Jun 2025 11:45:38 -0700 Subject: [PATCH 24/25] consistent error messages and more intrusive testing for them --- src/Utilities.jl | 81 ++++++++++++++++++++++---------------- test/Utilities/runtests.jl | 33 ++++++++++++++-- 2 files changed, 76 insertions(+), 38 deletions(-) diff --git a/src/Utilities.jl b/src/Utilities.jl index be1653f11..d8ce2498e 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -411,7 +411,7 @@ function initialize_processor!( structure_matrix::USorM, ) where {MM <: AbstractMatrix, USorM <: Union{UniformScaling, AbstractMatrix}} if length(get_data_mean(dd)) == 0 - append!(get_data_mean(dd), mean(data, dims = 2)) + push!(get_data_mean(dd), vec(mean(data, dims = 2))) end if length(get_encoder_mat(dd)) == 0 @@ -463,7 +463,7 @@ $(TYPEDSIGNATURES) Apply the `Decorrelater` encoder, on a columns-are-data matrix """ function encode_data(dd::Decorrelater, data::MM) where {MM <: AbstractMatrix} - data_mean = get_data_mean(dd) + data_mean = get_data_mean(dd)[1] encoder_mat = get_encoder_mat(dd)[1] return encoder_mat * (data .- data_mean) end @@ -474,7 +474,7 @@ $(TYPEDSIGNATURES) Apply the `Decorrelater` decoder, on a columns-are-data matrix """ function decode_data(dd::Decorrelater, data::MM) where {MM <: AbstractMatrix} - data_mean = get_data_mean(dd) + data_mean = get_data_mean(dd)[1] decoder_mat = get_decoder_mat(dd)[1] return decoder_mat * data .+ data_mean end @@ -537,7 +537,7 @@ Constructs the `CanonicalCorrelation` struct. Can optionally provide the keyword - `retain_var`[=1.0]: to project onto the leading singular vectors (of the input-output product) such that `retain_var` variance is retained. """ canonical_correlation(; retain_var::FT = Float64(1.0)) where {FT} = - CanonicalCorrelation(Any[], Any[], Any[], clamp(retain_var, FT(0), FT(1)) , AbstractString[]) + CanonicalCorrelation(Any[], Any[], Any[], clamp(retain_var, FT(0), FT(1)), AbstractString[]) """ $(TYPEDSIGNATURES) @@ -577,8 +577,8 @@ get_apply_to(cc::CanonicalCorrelation) = cc.apply_to function Base.show(io::IO, cc::CanonicalCorrelation) out = "CanonicalCorrelation:" - if length(get_apply_to(cc))>0 - out *=" apply_to=$(get_apply_to(cc)[1])" + if length(get_apply_to(cc)) > 0 + out *= " apply_to=$(get_apply_to(cc)[1])" end if get_retain_var(cc) < 1.0 out *= " retain_var=$(get_retain_var(cc))" @@ -594,15 +594,19 @@ function initialize_processor!( apply_to::AS, ) where {MM <: AbstractMatrix, AS <: AbstractString} + if apply_to ∉ ["in", "out"] + bad_apply_to(apply_to) + end + if length(get_apply_to(cc)) == 0 push!(get_apply_to(cc), apply_to) end if length(get_data_mean(cc)) == 0 if apply_to == "in" - append!(get_data_mean(cc), mean(in_data, dims = 2)) + push!(get_data_mean(cc), vec(mean(in_data, dims = 2))) elseif apply_to == "out" - append!(get_data_mean(cc), mean(out_data, dims = 2)) + push!(get_data_mean(cc), vec(mean(out_data, dims = 2))) end end @@ -642,8 +646,7 @@ function initialize_processor!( # mat' * Sx⁻¹ * Uxt encoder_mat = svdio_mat[:, 1:trunc]' * Diagonal(1 ./ svdi.S) * in_mat_sq' decoder_mat = in_mat_sq * Diagonal(svdi.S) * svdio_mat[:, 1:trunc] - else - apply_to == "out" + elseif apply_to == "out" out_dim = size(out_data, 1) svdio_mat = (size(svdio.U, 1) == out_dim) ? svdio.U : svdio.V # Vt * Sy⁻¹ * Uyt @@ -682,7 +685,7 @@ $(TYPEDSIGNATURES) Apply the `CanonicalCorrelation` encoder, on a columns-are-data matrix """ function encode_data(cc::CanonicalCorrelation, data::MM) where {MM <: AbstractMatrix} - data_mean = get_data_mean(cc) + data_mean = get_data_mean(cc)[1] encoder_mat = get_encoder_mat(cc)[1] return encoder_mat * (data .- data_mean) end @@ -693,7 +696,7 @@ $(TYPEDSIGNATURES) Apply the `CanonicalCorrelation` decoder, on a columns-are-data matrix """ function decode_data(cc::CanonicalCorrelation, data::MM) where {MM <: AbstractMatrix} - data_mean = get_data_mean(cc) + data_mean = get_data_mean(cc)[1] decoder_mat = get_decoder_mat(cc)[1] return decoder_mat * data .+ data_mean end @@ -777,10 +780,13 @@ function initialize_and_encode_data!( ) where {PDCP <: PairedDataContainerProcessor, USorM <: Union{UniformScaling, AbstractMatrix}, AS <: AbstractString} input_data, output_data = data initialize_processor!(dcp, input_data, output_data, structure_mat, apply_to) + if apply_to == "in" return encode_data(dcp, input_data) elseif apply_to == "out" return encode_data(dcp, output_data) + else + bad_apply_to(apply_to) end end @@ -795,6 +801,8 @@ function decode_data(dcp::PDCP, data, apply_to::AS) where {PDCP <: PairedDataCon return decode_data(dcp, input_data) elseif apply_to == "out" return decode_data(dcp, output_data) + else + bad_apply_to(apply_to) end end @@ -875,6 +883,14 @@ Create size-1 encoder schedule with a tuple of `(DataProcessor1(...), apply_to)` create_encoder_schedule(schedule_in::TT) where {TT <: Tuple} = create_encoder_schedule([schedule_in]) +function bad_apply_to(apply_to::AS) where {AS <: AbstractString} + throw( + ArgumentError( + "processer can only be applied to inputs (\"in\") or outputs (\"out\"). received $(apply_to). \n Please use `create_encoder_schedule` prior to encoding/decoding to ensure correct schedule format", + ), + ) +end + # Functions to encode/decode with uninitialized schedule (require structure matrices as input) """ $TYPEDSIGNATURES @@ -903,6 +919,8 @@ function encode_with_schedule( structure_matrix = processed_input_structure_mat elseif apply_to == "out" structure_matrix = processed_output_structure_mat + else + bad_apply_to(apply_to) end processed = initialize_and_encode_data!(processor, extract_data(processed_io_pairs), structure_matrix, apply_to) if apply_to == "in" @@ -930,11 +948,7 @@ function encode_with_schedule( ) where {VV <: AbstractVector, DC <: DataContainer, AS <: AbstractString} if in_or_out ∉ ["in", "out"] - throw( - ArgumentError( - "`in_or_out` must be either \"in\" (data is an input) or \"out\" (data is an output). Received $(in_or_out)", - ), - ) + bad_in_or_out(in_or_out) end processed_container = deepcopy(data_container) @@ -950,6 +964,15 @@ function encode_with_schedule( return processed_container end +function bad_in_or_out(in_or_out::AS) where {AS <: AbstractString} + throw( + ArgumentError( + "`in_or_out` must be either \"in\" (data is an input) or \"out\" (data is an output). Received $(in_or_out)", + ), + ) +end + + """ $TYPEDSIGNATURES @@ -962,11 +985,7 @@ function encode_with_schedule( ) where {VV <: AbstractVector, USorM <: Union{UniformScaling, AbstractMatrix}, AS <: AbstractString} if !(in_or_out ∈ ["in", "out"]) - throw( - ArgumentError( - "`in_or_out` must be either \"in\" (data is an input) or \"out\" (data is an output). Received $(in_or_out)", - ), - ) + bad_in_or_out(in_or_out) end processed_structure_matrix = deepcopy(structure_matrix) @@ -1012,6 +1031,8 @@ function decode_with_schedule( elseif apply_to == "out" processed_output_structure_mat = decode_structure_matrix(processor, processed_output_structure_mat) processed_io_pairs = PairedDataContainer(get_inputs(processed_io_pairs), processed) + else + bad_apply_to(apply_to) end end @@ -1029,12 +1050,8 @@ function decode_with_schedule( in_or_out::AS, ) where {VV <: AbstractVector, DC <: DataContainer, AS <: AbstractString} - if !(in_or_out ∈ ["in", "out"]) - throw( - ArgumentError( - "`in_or_out` must be either \"in\" (data is an input) or \"out\" (data is an output). Received $(in_or_out)", - ), - ) + if in_or_out ∉ ["in", "out"] + bad_in_or_out(in_or_out) end processed_container = deepcopy(data_container) @@ -1062,12 +1079,8 @@ function decode_with_schedule( in_or_out::AS, ) where {VV <: AbstractVector, USorM <: Union{UniformScaling, AbstractMatrix}, AS <: AbstractString} - if !(in_or_out ∈ ["in", "out"]) - throw( - ArgumentError( - "`in_or_out` must be either \"in\" (data is an input) or \"out\" (data is an output). Received $(in_or_out)", - ), - ) + if in_or_out ∉ ["in", "out"] + bad_in_or_out(in_or_out) end processed_structure_matrix = deepcopy(structure_matrix) diff --git a/test/Utilities/runtests.jl b/test/Utilities/runtests.jl index c83ecc278..490f77722 100644 --- a/test/Utilities/runtests.jl +++ b/test/Utilities/runtests.jl @@ -141,6 +141,7 @@ end (decorrelate_structure_mat(retain_var = 0.95), "in_and_out"), (canonical_correlation(retain_var = 0.95), "in_and_out"), ] + lossless = [fill(true, 6); fill(false, 4)] # are these lossy approximations? # functional test pipeline @@ -171,7 +172,11 @@ end test_vec = [[mean(dd), std(dd)] for dd in eachrow(test_dat)] test_mat = reduce(hcat, test_vec) - @test isapprox(norm(enc_covv - Diagonal(1 ./ test_mat[2, :]) * test_covv * Diagonal(1 ./ test_mat[2, :] )), 0.0, atol = tol * dim^2) + @test isapprox( + norm(enc_covv - Diagonal(1 ./ test_mat[2, :]) * test_covv * Diagonal(1 ./ test_mat[2, :])), + 0.0, + atol = tol * dim^2, + ) elseif name == "quartile" quartiles_vec = [quantile(dd, [0.25, 0.5, 0.75]) for dd in eachrow(enc_dat)] quartiles_mat = reduce(hcat, quartiles_vec) # 3 rows: Q1, Q2, and Q3 @@ -180,7 +185,12 @@ end test_vec = [quantile(dd, [0.25, 0.5, 0.75]) for dd in eachrow(test_dat)] test_mat = reduce(hcat, test_vec) @test isapprox( - norm(enc_covv - Diagonal(1 ./ (test_mat[3, :] - test_mat[1, :])) * test_covv * Diagonal(1 ./ (test_mat[3, :] - test_mat[1, :]))), + norm( + enc_covv - + Diagonal(1 ./ (test_mat[3, :] - test_mat[1, :])) * + test_covv * + Diagonal(1 ./ (test_mat[3, :] - test_mat[1, :])), + ), 0.0, atol = tol * dim^2, ) @@ -192,7 +202,12 @@ end test_vec = [[minimum(dd), maximum(dd)] for dd in eachrow(test_dat)] test_mat = reduce(hcat, test_vec) @test isapprox( - norm(enc_covv - Diagonal(1 ./ (test_mat[2, :] - test_mat[1, :])) * test_covv * Diagonal(1 ./ (test_mat[2, :] - test_mat[1, :]))), + norm( + enc_covv - + Diagonal(1 ./ (test_mat[2, :] - test_mat[1, :])) * + test_covv * + Diagonal(1 ./ (test_mat[2, :] - test_mat[1, :])), + ), 0.0, atol = tol * dim^2, ) @@ -264,6 +279,8 @@ end end + + # combine a few lossless encoding schedules (lossless requires samples>dims) samples = 150 # for full test coverage have samples in_dim < samples < out_dim in_data = rand(MvNormal(zeros(in_dim), prior_cov), samples) @@ -271,7 +288,7 @@ end io_pairs = PairedDataContainer(in_data, out_data) schedule_builder = [ - (zscore_scale(), "in_and_out"), + (zscore_scale(), "in_and_out"), (quartile_scale(), "in"), (decorrelate_sample_cov(), "in_and_out"), (minmax_scale(), "out"), @@ -283,6 +300,14 @@ end @test_logs (:warn,) create_encoder_schedule((canonical_correlation(), "bad")) @test_logs (:warn,) create_encoder_schedule((zscore_scale(), "bad")) + func = x -> (get_inputs(x), get_outputs(x)) + bad_encoder_schedule = [(canonical_correlation(), func, "bad")] + @test_throws ArgumentError encode_with_schedule(bad_encoder_schedule, io_pairs, prior_cov, obs_noise_cov) + @test_throws ArgumentError decode_with_schedule(bad_encoder_schedule, io_pairs, prior_cov, obs_noise_cov) + @test_throws ArgumentError initialize_and_encode_data!(canonical_correlation(), func(io_pairs), prior_cov, "bad") + @test_throws ArgumentError decode_data(canonical_correlation(), func(io_pairs), "bad") + + encoder_schedule = create_encoder_schedule(schedule_builder) # encode the data using the schedule From 93caaaafb63d2993fbb947737d91199884bf62a0 Mon Sep 17 00:00:00 2001 From: odunbar Date: Wed, 25 Jun 2025 13:07:30 -0700 Subject: [PATCH 25/25] affine -> elementwise --- src/Utilities.jl | 96 +++++++++++++++++++------------------- test/Utilities/runtests.jl | 8 ++-- 2 files changed, 52 insertions(+), 52 deletions(-) diff --git a/src/Utilities.jl b/src/Utilities.jl index d8ce2498e..7a85a811d 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -16,7 +16,7 @@ export orig2zscore export zscore2orig export PairedDataContainerProcessor, DataContainerProcessor -export UnivariateAffineScaling, AffineScaler, QuartileScaling, MinMaxScaling, ZScoreScaling +export UnivariateAffineScaling, ElementwiseScaler, QuartileScaling, MinMaxScaling, ZScoreScaling export quartile_scale, minmax_scale, zscore_scale export Decorrelater, decorrelate_sample_cov, decorrelate_structure_mat, decorrelate export get_type, @@ -123,7 +123,7 @@ abstract type ZScoreScaling <: UnivariateAffineScaling end """ $(TYPEDEF) -The AffineScaler{T} will create an encoding of the data_container via affine transformations. +The ElementwiseScaler{T} will create an encoding of the data_container via elementwise affine transformations. Different methods `T` will build different transformations: - [`quartile_scale`](@ref) : creates `QuartileScaling`, @@ -132,7 +132,7 @@ Different methods `T` will build different transformations: and are accessed with [`get_type`](@ref) """ -struct AffineScaler{T, VV <: AbstractVector} <: DataContainerProcessor +struct ElementwiseScaler{T, VV <: AbstractVector} <: DataContainerProcessor "storage for the shift applied to data" shift::VV "storage for the scaling" @@ -142,108 +142,108 @@ end """ $(TYPEDSIGNATURES) -Constructs `AffineScaler{QuartileScaling}` processor. +Constructs `ElementwiseScaler{QuartileScaling}` processor. As part of an encoder schedule, it will apply the transform ``\\frac{x - Q2(x)}{Q3(x) - Q1(x)}`` to each data dimension. Also known as "robust scaling" """ -quartile_scale() = AffineScaler(QuartileScaling) +quartile_scale() = ElementwiseScaler(QuartileScaling) """ $(TYPEDSIGNATURES) -Constructs `AffineScaler{MinMaxScaling}` processor. +Constructs `ElementwiseScaler{MinMaxScaling}` processor. As part of an encoder schedule, this will apply the transform ``\\frac{x - \\min(x)}{\\max(x) - \\min(x)}`` to each data dimension. """ -minmax_scale() = AffineScaler(MinMaxScaling) +minmax_scale() = ElementwiseScaler(MinMaxScaling) """ $(TYPEDSIGNATURES) -Constructs `AffineScaler{ZScoreScaling}` processor. +Constructs `ElementwiseScaler{ZScoreScaling}` processor. As part of an encoder schedule, this will apply the transform ``\\frac{x-\\mu}{\\sigma}``, (where ``x\\sim N(\\mu,\\sigma)``), to each data dimension. For multivariate standardization, see [`Decorrelater`](@ref) """ -zscore_scale() = AffineScaler(ZScoreScaling) +zscore_scale() = ElementwiseScaler(ZScoreScaling) -AffineScaler(::Type{UAS}) where {UAS <: UnivariateAffineScaling} = - AffineScaler{UAS, Vector{Float64}}(Float64[], Float64[]) +ElementwiseScaler(::Type{UAS}) where {UAS <: UnivariateAffineScaling} = + ElementwiseScaler{UAS, Vector{Float64}}(Float64[], Float64[]) """ $(TYPEDSIGNATURES) Gets the UnivariateAffineScaling type `T` """ -get_type(as::AffineScaler{T}) where {T} = T +get_type(es::ElementwiseScaler{T}) where {T} = T """ $(TYPEDSIGNATURES) -Gets the `shift` field of the `AffineScaler` +Gets the `shift` field of the `ElementwiseScaler` """ -get_shift(as::AffineScaler) = as.shift +get_shift(es::ElementwiseScaler) = es.shift """ $(TYPEDSIGNATURES) -Gets the `scale` field of the `AffineScaler` +Gets the `scale` field of the `ElementwiseScaler` """ -get_scale(as::AffineScaler) = as.scale +get_scale(es::ElementwiseScaler) = es.scale -function Base.show(io::IO, as::AffineScaler) - out = "AffineScaler: $(get_type(as))" +function Base.show(io::IO, es::ElementwiseScaler) + out = "ElementwiseScaler: $(get_type(es))" print(io, out) end function initialize_processor!( - as::AffineScaler, + es::ElementwiseScaler, data::MM, T::Type{QS}, ) where {MM <: AbstractMatrix, QS <: QuartileScaling} quartiles_vec = [quantile(dd, [0.25, 0.5, 0.75]) for dd in eachrow(data)] quartiles_mat = reduce(hcat, quartiles_vec) # 3 rows: Q1, Q2, and Q3 - append!(get_shift(as), quartiles_mat[2, :]) - append!(get_scale(as), (quartiles_mat[3, :] - quartiles_mat[1, :])) + append!(get_shift(es), quartiles_mat[2, :]) + append!(get_scale(es), (quartiles_mat[3, :] - quartiles_mat[1, :])) end function initialize_processor!( - as::AffineScaler, + es::ElementwiseScaler, data::MM, T::Type{MMS}, ) where {MM <: AbstractMatrix, MMS <: MinMaxScaling} minmax_vec = [[minimum(dd), maximum(dd)] for dd in eachrow(data)] minmax_mat = reduce(hcat, minmax_vec) # 2 rows: min max - append!(get_shift(as), minmax_mat[1, :]) - append!(get_scale(as), (minmax_mat[2, :] - minmax_mat[1, :])) + append!(get_shift(es), minmax_mat[1, :]) + append!(get_scale(es), (minmax_mat[2, :] - minmax_mat[1, :])) end function initialize_processor!( - as::AffineScaler, + es::ElementwiseScaler, data::MM, T::Type{ZSS}, ) where {MM <: AbstractMatrix, ZSS <: ZScoreScaling} stat_vec = [[mean(dd), std(dd)] for dd in eachrow(data)] stat_mat = reduce(hcat, stat_vec) # 2 rows: mean, std - append!(get_shift(as), stat_mat[1, :]) - append!(get_scale(as), stat_mat[2, :]) + append!(get_shift(es), stat_mat[1, :]) + append!(get_scale(es), stat_mat[2, :]) end -function initialize_processor!(as::AffineScaler, data::MM) where {MM <: AbstractMatrix} - if length(get_shift(as)) == 0 - T = get_type(as) - initialize_processor!(as, data, T) +function initialize_processor!(es::ElementwiseScaler, data::MM) where {MM <: AbstractMatrix} + if length(get_shift(es)) == 0 + T = get_type(es) + initialize_processor!(es, data, T) end end """ $(TYPEDSIGNATURES) -Apply the `AffineScaler` encoder, on a columns-are-data matrix +Apply the `ElementwiseScaler` encoder, on a columns-are-data matrix """ -function encode_data(as::AffineScaler, data::MM) where {MM <: AbstractMatrix} +function encode_data(es::ElementwiseScaler, data::MM) where {MM <: AbstractMatrix} out = deepcopy(data) for i in 1:size(out, 1) - out[i, :] .-= get_shift(as)[i] - out[i, :] /= get_scale(as)[i] + out[i, :] .-= get_shift(es)[i] + out[i, :] /= get_scale(es)[i] end return out end @@ -251,13 +251,13 @@ end """ $(TYPEDSIGNATURES) -Apply the `AffineScaler` decoder, on a columns-are-data matrix +Apply the `ElementwiseScaler` decoder, on a columns-are-data matrix """ -function decode_data(as::AffineScaler, data::MM) where {MM <: AbstractMatrix} +function decode_data(es::ElementwiseScaler, data::MM) where {MM <: AbstractMatrix} out = deepcopy(data) for i in 1:size(out, 1) - out[i, :] *= get_scale(as)[i] - out[i, :] .+= get_shift(as)[i] + out[i, :] *= get_scale(es)[i] + out[i, :] .+= get_shift(es)[i] end return out end @@ -265,28 +265,28 @@ end """ $(TYPEDSIGNATURES) -Computes and populates the `shift` and `scale` fields for the `AffineScaler` +Computes and populates the `shift` and `scale` fields for the `ElementwiseScaler` """ -initialize_processor!(as::AffineScaler, data::MM, structure_matrix) where {MM <: AbstractMatrix} = - initialize_processor!(as, data) +initialize_processor!(es::ElementwiseScaler, data::MM, structure_matrix) where {MM <: AbstractMatrix} = + initialize_processor!(es, data) """ $(TYPEDSIGNATURES) -Apply the `AffineScaler` encoder to a provided structure matrix +Apply the `ElementwiseScaler` encoder to a provided structure matrix """ -function encode_structure_matrix(as::AffineScaler, structure_matrix::MM) where {MM <: AbstractMatrix} - return Diagonal(1 ./ get_scale(as)) * structure_matrix * Diagonal(1 ./ get_scale(as)) +function encode_structure_matrix(es::ElementwiseScaler, structure_matrix::MM) where {MM <: AbstractMatrix} + return Diagonal(1 ./ get_scale(es)) * structure_matrix * Diagonal(1 ./ get_scale(es)) end """ $(TYPEDSIGNATURES) -Apply the `AffineScaler` decoder to a provided structure matrix +Apply the `ElementwiseScaler` decoder to a provided structure matrix """ -function decode_structure_matrix(as::AffineScaler, enc_structure_matrix::MM) where {MM <: AbstractMatrix} - return Diagonal(get_scale(as)) * enc_structure_matrix * Diagonal(get_scale(as)) +function decode_structure_matrix(es::ElementwiseScaler, enc_structure_matrix::MM) where {MM <: AbstractMatrix} + return Diagonal(get_scale(es)) * enc_structure_matrix * Diagonal(get_scale(es)) end diff --git a/test/Utilities/runtests.jl b/test/Utilities/runtests.jl index 490f77722..687532b22 100644 --- a/test/Utilities/runtests.jl +++ b/test/Utilities/runtests.jl @@ -81,12 +81,12 @@ end zs = zscore_scale() mm = minmax_scale() qq = quartile_scale() - QQ = AffineScaler{QuartileScaling, Vector{Int}}([1], [2]) - @test isa(zs, AffineScaler) + QQ = ElementwiseScaler{QuartileScaling, Vector{Int}}([1], [2]) + @test isa(zs, ElementwiseScaler) @test get_type(zs) == ZScoreScaling - @test isa(mm, AffineScaler) + @test isa(mm, ElementwiseScaler) @test get_type(mm) == MinMaxScaling - @test isa(qq, AffineScaler) + @test isa(qq, ElementwiseScaler) @test get_type(qq) == QuartileScaling @test get_shift(QQ) == [1] @test get_scale(QQ) == [2]