From 5c176e59c01fbd7911cfd9fd81b6a13112df2ed1 Mon Sep 17 00:00:00 2001 From: Arne Bouillon Date: Mon, 7 Jul 2025 18:15:55 -0700 Subject: [PATCH 1/5] Reduce number of methods; give data processors separate files --- src/Utilities.jl | 841 +++---------------------- src/Utilities/canonical_correlation.jl | 218 +++++++ src/Utilities/decorrelator.jl | 216 +++++++ src/Utilities/elementwise_scaler.jl | 179 ++++++ test/Utilities/runtests.jl | 3 +- 5 files changed, 689 insertions(+), 768 deletions(-) create mode 100644 src/Utilities/canonical_correlation.jl create mode 100644 src/Utilities/decorrelator.jl create mode 100644 src/Utilities/elementwise_scaler.jl diff --git a/src/Utilities.jl b/src/Utilities.jl index e222fdb4b..1272b21ae 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -12,22 +12,9 @@ EnsembleKalmanProcess = EnsembleKalmanProcesses.EnsembleKalmanProcess using ..DataContainers export get_training_points -export orig2zscore -export zscore2orig export PairedDataContainerProcessor, DataContainerProcessor -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, - 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, 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 +export create_encoder_schedule, initialize_and_encode_with_schedule!, encode_with_schedule, decode_with_schedule @@ -67,15 +54,12 @@ function get_training_points( return training_points end -# Data processing tooling: -abstract type PairedDataContainerProcessor end # tools that operate on inputs and outputs -abstract type DataContainerProcessor end # tools that operate on only inputs or outputs +# Data processing tooling: -abstract type UnivariateAffineScaling end -abstract type QuartileScaling <: UnivariateAffineScaling end -abstract type MinMaxScaling <: UnivariateAffineScaling end -abstract type ZScoreScaling <: UnivariateAffineScaling end +abstract type DataProcessor end +abstract type PairedDataContainerProcessor <: DataProcessor end # tools that operate on inputs and outputs +abstract type DataContainerProcessor <: DataProcessor end # tools that operate on only inputs or outputs # define how to have equality Base.:(==)(a::DCP, b::DCP) where {DCP <: DataContainerProcessor} = @@ -83,693 +67,48 @@ Base.:(==)(a::DCP, b::DCP) where {DCP <: DataContainerProcessor} = Base.:(==)(a::PDCP, b::PDCP) where {PDCP <: PairedDataContainerProcessor} = all(getfield(a, f) == getfield(b, f) for f in fieldnames(PDCP)) -# Processors -""" -$(TYPEDEF) - -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`, -- [`minmax_scale`](@ref) : creates `MinMaxScaling` -- [`zscore_scale`](@ref) : creates `ZScoreScaling` - -and are accessed with [`get_type`](@ref) -""" -struct ElementwiseScaler{T, VV <: AbstractVector} <: DataContainerProcessor - "storage for the shift applied to data" - shift::VV - "storage for the scaling" - scale::VV -end - -""" -$(TYPEDSIGNATURES) - -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() = ElementwiseScaler(QuartileScaling) - -""" -$(TYPEDSIGNATURES) - -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() = ElementwiseScaler(MinMaxScaling) - -""" -$(TYPEDSIGNATURES) - -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() = ElementwiseScaler(ZScoreScaling) - -ElementwiseScaler(::Type{UAS}) where {UAS <: UnivariateAffineScaling} = - ElementwiseScaler{UAS, Vector{Float64}}(Float64[], Float64[]) - -""" -$(TYPEDSIGNATURES) - -Gets the UnivariateAffineScaling type `T` -""" -get_type(es::ElementwiseScaler{T}) where {T} = T - -""" -$(TYPEDSIGNATURES) - -Gets the `shift` field of the `ElementwiseScaler` -""" -get_shift(es::ElementwiseScaler) = es.shift - -""" -$(TYPEDSIGNATURES) - -Gets the `scale` field of the `ElementwiseScaler` -""" -get_scale(es::ElementwiseScaler) = es.scale - -function Base.show(io::IO, es::ElementwiseScaler) - out = "ElementwiseScaler: $(get_type(es))" - print(io, out) -end - -function initialize_processor!( - 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(es), quartiles_mat[2, :]) - append!(get_scale(es), (quartiles_mat[3, :] - quartiles_mat[1, :])) -end - -function initialize_processor!( - 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(es), minmax_mat[1, :]) - append!(get_scale(es), (minmax_mat[2, :] - minmax_mat[1, :])) -end - -function initialize_processor!( - 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(es), stat_mat[1, :]) - append!(get_scale(es), stat_mat[2, :]) -end - -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 `ElementwiseScaler` encoder, on a columns-are-data matrix -""" -function encode_data(es::ElementwiseScaler, data::MM) where {MM <: AbstractMatrix} - out = deepcopy(data) - for i in 1:size(out, 1) - out[i, :] .-= get_shift(es)[i] - out[i, :] /= get_scale(es)[i] - end - return out -end - -""" -$(TYPEDSIGNATURES) - -Apply the `ElementwiseScaler` decoder, on a columns-are-data matrix -""" -function decode_data(es::ElementwiseScaler, data::MM) where {MM <: AbstractMatrix} - out = deepcopy(data) - for i in 1:size(out, 1) - out[i, :] *= get_scale(es)[i] - out[i, :] .+= get_shift(es)[i] - end - return out -end - -""" -$(TYPEDSIGNATURES) - -Computes and populates the `shift` and `scale` fields for the `ElementwiseScaler` -""" -initialize_processor!(es::ElementwiseScaler, data::MM, structure_matrix) where {MM <: AbstractMatrix} = - initialize_processor!(es, data) - - -""" -$(TYPEDSIGNATURES) - -Apply the `ElementwiseScaler` encoder to a provided structure matrix -""" -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 `ElementwiseScaler` decoder to a provided structure matrix -""" -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 - - -""" -$(TYPEDEF) - -Decorrelate the data via taking an SVD decomposition and projecting onto the singular-vectors. - -Preferred construction is with the methods -- [`decorrelate_structure_mat`](@ref) -- [`decorrelate_sample_cov`](@ref) -- [`decorrelate`](@ref) - -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, - -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, AS <: AbstractString} <: DataContainerProcessor - "storage for the data mean" - data_mean::VV1 - "the matrix used to perform encoding" - encoder_mat::VV2 - "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 - "Switch to choose what form of matrix to use to decorrelate the data" - decorrelate_with::AS -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` [=`"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([], [], [], clamp(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([], [], [], clamp(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_structure_mat(; retain_var::FT = Float64(1.0)) where {FT} = - Decorrelater([], [], [], clamp(retain_var, FT(0), FT(1)), "structure_mat") - -""" -$(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 `decorrelate_with` field of the `Decorrelater`. -""" -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))" - end - 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, - structure_matrix::USorM, -) where {MM <: AbstractMatrix, USorM <: Union{UniformScaling, AbstractMatrix}} - if length(get_data_mean(dd)) == 0 - push!(get_data_mean(dd), vec(mean(data, dims = 2))) - end - - if length(get_encoder_mat(dd)) == 0 - - # Can do tsvd here for large matrices - decorrelate_with = get_decorrelate_with(dd) - if decorrelate_with == "structure_mat" - svdA = svd(structure_matrix) - rk = rank(structure_matrix) - elseif decorrelate_with == "sample_cov" - cd = cov(data, dims = 2) - svdA = svd(cd) - rk = rank(cd) - elseif decorrelate_with == "combined" - spluscd = structure_matrix + cov(data, dims = 2) - svdA = svd(spluscd) - rk = rank(spluscd) - else - 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 - sv_cumsum = cumsum(svdA.S) / sum(svdA.S) # variance contributions are (sing_val) for these matrices - trunc_val = minimum(findall(x -> (x > ret_var), sv_cumsum)) - @info " truncating at $(trunc_val)/$(length(sv_cumsum)) retaining $(100.0*sv_cumsum[trunc_val])% of the variance of the structure matrix" - else - trunc_val = rk - if rk < size(data, 1) - @info " truncating at $(trunc_val)/$(size(data,1)), as low-rank data detected" - end - end - - sqrt_inv_sv = Diagonal(1.0 ./ sqrt.(svdA.S[1:trunc_val])) - sqrt_sv = Diagonal(sqrt.(svdA.S[1:trunc_val])) - # as we have svd of cov-matrix we can use U or Vt - encoder_mat = sqrt_inv_sv * svdA.Vt[1:trunc_val, :] - decoder_mat = svdA.Vt[1:trunc_val, :]' * sqrt_sv - - push!(get_encoder_mat(dd), encoder_mat) - push!(get_decoder_mat(dd), decoder_mat) - 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)[1] - 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)[1] - 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' -end - -# ...struct VariationalAutoEncoder <: DataContainerProcessor end - -# PDCProcessors -# 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 - -Preferred construction is with the [`canonical_correlation`](@ref) method - -# 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 - -""" -$(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::FT = Float64(1.0)) where {FT} = - CanonicalCorrelation(Any[], Any[], Any[], clamp(retain_var, FT(0), FT(1)), 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) - - 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 *= " retain_var=$(get_retain_var(cc))" - 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 apply_to ∉ ["in", "out"] +function _encode_data(proc::P, data, apply_to::AS) where {P <: DataProcessor, AS <: AbstractString} + input_data, output_data = data + if apply_to == "in" + return encode_data(proc, input_data) + elseif apply_to == "out" + return encode_data(proc, output_data) + else 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" - push!(get_data_mean(cc), vec(mean(in_data, dims = 2))) - elseif apply_to == "out" - push!(get_data_mean(cc), vec(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 - svdi = svd(in_data .- mean(in_data, dims = 2)) - svdo = svd(out_data .- mean(out_data, dims = 2)) - - # 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) - - svdio = svd(in_mat_nonsq * out_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)^2for these matrices - trunc_val = minimum(findall(x -> (x > ret_var), sv_cumsum)) - @info " truncating at $(trunc_val)/$(length(sv_cumsum)) retaining $(100.0*sv_cumsum[trunc_val])% of the variance in the joint space" - else - trunc_val = min(rank(in_data), rank(out_data)) - end - in_dim = size(in_data, 1) - in_svdio_mat, out_svdio_mat = (size(svdio.U, 1) == in_dim) ? (svdio.U, svdio.V) : (svdio.V, svdio.U') - if apply_to == "in" - # mat' * Sx⁻¹ * Uxt - encoder_mat = in_svdio_mat[:, 1:trunc_val]' * Diagonal(1 ./ svdi.S) * in_mat_sq' - decoder_mat = in_mat_sq * Diagonal(svdi.S) * in_svdio_mat[:, 1:trunc_val] - - elseif apply_to == "out" - out_dim = size(out_data, 1) - # Vt * Sy⁻¹ * Uyt - encoder_mat = out_svdio_mat[:, 1:trunc_val]' * Diagonal(1 ./ svdo.S) * out_mat_sq' - decoder_mat = out_mat_sq * Diagonal(svdo.S) * out_svdio_mat[:, 1:trunc_val] - end - - push!(get_encoder_mat(cc), encoder_mat) - push!(get_decoder_mat(cc), decoder_mat) - - # Note: To check CCA: - # u = in_encoder * (in_data .- mean(in_data, dims=2)) - # v = out_encoder * (out_data .- mean(out_data, dims=2)) - # u * u' = v * v' = I, - # v * u' = u * v' = Diagonal(svdio.S[1:trunc_val]) - 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, - out_data::MM, - structure_matrix, - apply_to::AS, -) 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)[1] - 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)[1] - 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' -end - - - - - -#### - - -# generic functions to initialize, encode, and decode data. With Data processors or PairedData processors -function initialize_and_encode_data!( - dcp::DCP, - data::MM, - structure_mat::USorM, -) where {DCP <: DataContainerProcessor, USorM <: Union{UniformScaling, AbstractMatrix}, MM <: AbstractMatrix} - initialize_processor!(dcp, data, structure_mat) - 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, - structure_mat::USorM, - apply_to::AS, -) where { - DCP <: DataContainerProcessor, - MM <: AbstractMatrix, - USorM <: Union{UniformScaling, AbstractMatrix}, - 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, - structure_mat::USorM, - apply_to::AS, -) where {PDCP <: PairedDataContainerProcessor, USorM <: Union{UniformScaling, AbstractMatrix}, AS <: AbstractString} +function _decode_data(proc::P, data, apply_to::AS) where {P <: DataProcessor, 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) + return decode_data(proc, input_data) elseif apply_to == "out" - return encode_data(dcp, output_data) + return decode_data(proc, output_data) else bad_apply_to(apply_to) 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} +function _initialize_and_encode_data!(proc::P, data, structure_mats, apply_to::AS) where {P <: DataProcessor, AS <: AbstractString} input_data, output_data = data - if apply_to == "in" - return decode_data(dcp, input_data) - elseif apply_to == "out" - return decode_data(dcp, output_data) + input_structure_mat, output_structure_mat = structure_mats + + if P isa PairedDataContainerProcessor + initialize_processor!(proc, input_data, output_data, apply_to == "in" ? input_structure_mat : output_structure_mat, apply_to) else - bad_apply_to(apply_to) + if apply_to == "in" + initialize_processor!(proc, input_data, input_structure_mat) + else + initialize_processor!(proc, output_data, output_structure_mat) + end end + + _encode_data(proc, data, apply_to) end """ @@ -799,44 +138,18 @@ enc_schedule = [ 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 $(apply_to), 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 $(apply_to), ignoring processor $(processor)..." - ) - end + if apply_to ∈ ["in", "out"] + push!(encoder_schedule, (processor, apply_to)) + elseif apply_to == "in_and_out" + push!(encoder_schedule, (processor, "in")) + push!(encoder_schedule, (deepcopy(processor), "out")) + else + @warn( + "Expected schedule keywords ∈ {\"in\",\"out\",\"in_and_out\"}. Received $(apply_to), ignoring processor $(processor)..." + ) end end @@ -848,21 +161,13 @@ 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 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!( +function initialize_and_encode_with_schedule!( encoder_schedule::VV, io_pairs::PDC, input_structure_mat::USorM1, @@ -878,21 +183,16 @@ function encode_with_schedule!( 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 + for (processor, apply_to) in encoder_schedule @info "Initialize encoding of data: \"$(apply_to)\" with $(processor)" + + processed = _initialize_and_encode_data!(processor, processed_io_pairs, (processed_input_structure_mat, processed_output_structure_mat), apply_to) + if apply_to == "in" - 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" - processed_input_structure_mat = encode_structure_matrix(processor, structure_matrix) + processed_input_structure_mat = encode_structure_matrix(processor, processed_input_structure_mat) processed_io_pairs = PairedDataContainer(processed, get_outputs(processed_io_pairs)) elseif apply_to == "out" - processed_output_structure_mat = encode_structure_matrix(processor, structure_matrix) + processed_output_structure_mat = encode_structure_matrix(processor, processed_output_structure_mat) processed_io_pairs = PairedDataContainer(get_inputs(processed_io_pairs), processed) end end @@ -900,6 +200,7 @@ function encode_with_schedule!( return processed_io_pairs, processed_input_structure_mat, processed_output_structure_mat end + # Functions to encode/decode with initialized schedule """ $TYPEDSIGNATURES @@ -911,33 +212,22 @@ function encode_with_schedule( data_container::DC, in_or_out::AS, ) where {VV <: AbstractVector, DC <: DataContainer, AS <: AbstractString} - if in_or_out ∉ ["in", "out"] bad_in_or_out(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 + for (processor, apply_to) in encoder_schedule if apply_to == in_or_out processed = encode_data(processor, get_data(processed_container)) processed_container = DataContainer(processed) - end end 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 @@ -948,14 +238,13 @@ function encode_with_schedule( structure_matrix::USorM, in_or_out::AS, ) where {VV <: AbstractVector, USorM <: Union{UniformScaling, AbstractMatrix}, AS <: AbstractString} - if !(in_or_out ∈ ["in", "out"]) bad_in_or_out(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 + for (processor, apply_to) in encoder_schedule if apply_to == in_or_out processed_structure_matrix = encode_structure_matrix(processor, processed_structure_matrix) end @@ -980,15 +269,14 @@ function decode_with_schedule( USorM2 <: Union{UniformScaling, AbstractMatrix}, PDC <: PairedDataContainer, } - processed_io_pairs = deepcopy(io_pairs) 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)) - (processor, extract_data, apply_to) = encoder_schedule[idx] - processed = decode_data(processor, extract_data(processed_io_pairs), apply_to) + (processor, apply_to) = encoder_schedule[idx] + processed = _decode_data(processor, processed_io_pairs, apply_to) if apply_to == "in" processed_input_structure_mat = decode_structure_matrix(processor, processed_input_structure_mat) @@ -1014,7 +302,6 @@ function decode_with_schedule( data_container::DC, in_or_out::AS, ) where {VV <: AbstractVector, DC <: DataContainer, AS <: AbstractString} - if in_or_out ∉ ["in", "out"] bad_in_or_out(in_or_out) end @@ -1022,11 +309,10 @@ 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] + (processor, apply_to) = encoder_schedule[idx] if apply_to == in_or_out processed = decode_data(processor, get_data(processed_container)) processed_container = DataContainer(processed) - end end @@ -1043,7 +329,6 @@ function decode_with_schedule( structure_matrix::USorM, in_or_out::AS, ) where {VV <: AbstractVector, USorM <: Union{UniformScaling, AbstractMatrix}, AS <: AbstractString} - if in_or_out ∉ ["in", "out"] bad_in_or_out(in_or_out) end @@ -1051,7 +336,7 @@ 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] + (processor, apply_to) = encoder_schedule[idx] if apply_to == in_or_out processed_structure_matrix = decode_structure_matrix(processor, processed_structure_matrix) end @@ -1061,5 +346,29 @@ function decode_with_schedule( end +# Errors + +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 + +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 + + +# Processors + +include("Utilities/canonical_correlation.jl") +include("Utilities/decorrelator.jl") +include("Utilities/elementwise_scaler.jl") end # module diff --git a/src/Utilities/canonical_correlation.jl b/src/Utilities/canonical_correlation.jl new file mode 100644 index 000000000..5ee1426c0 --- /dev/null +++ b/src/Utilities/canonical_correlation.jl @@ -0,0 +1,218 @@ +# included in Utilities.jl + +export CanonicalCorrelation, canonical_correlation, get_apply_to +export get_shift, get_scale, get_data_mean, get_encoder_mat, get_decoder_mat, get_retain_var, get_decorrelate_with + +""" +$(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) +""" +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 + +""" +$(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::FT = Float64(1.0)) where {FT} = + CanonicalCorrelation(Any[], Any[], Any[], clamp(retain_var, FT(0), FT(1)), 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) + + 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 *= " retain_var=$(get_retain_var(cc))" + 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 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" + push!(get_data_mean(cc), vec(mean(in_data, dims = 2))) + elseif apply_to == "out" + push!(get_data_mean(cc), vec(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 + svdi = svd(in_data .- mean(in_data, dims = 2)) + svdo = svd(out_data .- mean(out_data, dims = 2)) + + # 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) + + svdio = svd(in_mat_nonsq * out_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)^2for these matrices + trunc_val = minimum(findall(x -> (x > ret_var), sv_cumsum)) + @info " truncating at $(trunc_val)/$(length(sv_cumsum)) retaining $(100.0*sv_cumsum[trunc_val])% of the variance in the joint space" + else + trunc_val = min(rank(in_data), rank(out_data)) + end + in_dim = size(in_data, 1) + in_svdio_mat, out_svdio_mat = (size(svdio.U, 1) == in_dim) ? (svdio.U, svdio.V) : (svdio.V, svdio.U') + if apply_to == "in" + # mat' * Sx⁻¹ * Uxt + encoder_mat = in_svdio_mat[:, 1:trunc_val]' * Diagonal(1 ./ svdi.S) * in_mat_sq' + decoder_mat = in_mat_sq * Diagonal(svdi.S) * in_svdio_mat[:, 1:trunc_val] + + elseif apply_to == "out" + out_dim = size(out_data, 1) + # Vt * Sy⁻¹ * Uyt + encoder_mat = out_svdio_mat[:, 1:trunc_val]' * Diagonal(1 ./ svdo.S) * out_mat_sq' + decoder_mat = out_mat_sq * Diagonal(svdo.S) * out_svdio_mat[:, 1:trunc_val] + end + + push!(get_encoder_mat(cc), encoder_mat) + push!(get_decoder_mat(cc), decoder_mat) + + # Note: To check CCA: + # u = in_encoder * (in_data .- mean(in_data, dims=2)) + # v = out_encoder * (out_data .- mean(out_data, dims=2)) + # u * u' = v * v' = I, + # v * u' = u * v' = Diagonal(svdio.S[1:trunc_val]) + 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, + out_data::MM, + structure_matrix, + apply_to::AS, +) 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)[1] + 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)[1] + 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' +end diff --git a/src/Utilities/decorrelator.jl b/src/Utilities/decorrelator.jl new file mode 100644 index 000000000..69906a9b6 --- /dev/null +++ b/src/Utilities/decorrelator.jl @@ -0,0 +1,216 @@ +# included in Utilities.jl + +export Decorrelator, decorrelate_sample_cov, decorrelate_structure_mat, decorrelate +export get_data_mean, get_encoder_mat, get_decoder_mat, get_retain_var, get_decorrelate_with + +""" +$(TYPEDEF) + +Decorrelate the data via taking an SVD decomposition and projecting onto the singular-vectors. + +Preferred construction is with the methods +- [`decorrelate_structure_mat`](@ref) +- [`decorrelate_sample_cov`](@ref) +- [`decorrelate`](@ref) + +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, + +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 Decorrelator{VV1, VV2, VV3, FT, AS <: AbstractString} <: DataContainerProcessor + "storage for the data mean" + data_mean::VV1 + "the matrix used to perform encoding" + encoder_mat::VV2 + "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 + "Switch to choose what form of matrix to use to decorrelate the data" + decorrelate_with::AS +end + +""" +$(TYPEDSIGNATURES) + +Constructs the `Decorrelator` 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` [=`"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} = + Decorrelator([], [], [], clamp(retain_var, FT(0), FT(1)), decorrelate_with) + +""" +$(TYPEDSIGNATURES) + +Constructs the `Decorrelator` 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} = + Decorrelator([], [], [], clamp(retain_var, FT(0), FT(1)), "sample_cov") + +""" +$(TYPEDSIGNATURES) + +Constructs the `Decorrelator` 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_structure_mat(; retain_var::FT = Float64(1.0)) where {FT} = + Decorrelator([], [], [], clamp(retain_var, FT(0), FT(1)), "structure_mat") + +""" +$(TYPEDSIGNATURES) + +returns the `data_mean` field of the `Decorrelator`. +""" +get_data_mean(dd::Decorrelator) = dd.data_mean + +""" +$(TYPEDSIGNATURES) + +returns the `encoder_mat` field of the `Decorrelator`. +""" +get_encoder_mat(dd::Decorrelator) = dd.encoder_mat + +""" +$(TYPEDSIGNATURES) + +returns the `decoder_mat` field of the `Decorrelator`. +""" +get_decoder_mat(dd::Decorrelator) = dd.decoder_mat + +""" +$(TYPEDSIGNATURES) + +returns the `retain_var` field of the `Decorrelator`. +""" +get_retain_var(dd::Decorrelator) = dd.retain_var + +""" +$(TYPEDSIGNATURES) + +returns the `decorrelate_with` field of the `Decorrelator`. +""" +get_decorrelate_with(dd::Decorrelator) = dd.decorrelate_with + +function Base.show(io::IO, dd::Decorrelator) + out = "Decorrelator" + out *= ": decorrelate_with=$(get_decorrelate_with(dd))" + if get_retain_var(dd) < 1.0 + out *= ", retain_var=$(get_retain_var(dd))" + end + print(io, out) +end + +""" +$(TYPEDSIGNATURES) + +Computes and populates the `data_mean` and `encoder_mat` and `decoder_mat` fields for the `Decorrelator` +""" +function initialize_processor!( + dd::Decorrelator, + data::MM, + structure_matrix::USorM, +) where {MM <: AbstractMatrix, USorM <: Union{UniformScaling, AbstractMatrix}} + if length(get_data_mean(dd)) == 0 + push!(get_data_mean(dd), vec(mean(data, dims = 2))) + end + + if length(get_encoder_mat(dd)) == 0 + + # Can do tsvd here for large matrices + decorrelate_with = get_decorrelate_with(dd) + if decorrelate_with == "structure_mat" + svdA = svd(structure_matrix) + rk = rank(structure_matrix) + elseif decorrelate_with == "sample_cov" + cd = cov(data, dims = 2) + svdA = svd(cd) + rk = rank(cd) + elseif decorrelate_with == "combined" + spluscd = structure_matrix + cov(data, dims = 2) + svdA = svd(spluscd) + rk = rank(spluscd) + else + 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 + sv_cumsum = cumsum(svdA.S) / sum(svdA.S) # variance contributions are (sing_val) for these matrices + trunc_val = minimum(findall(x -> (x > ret_var), sv_cumsum)) + @info " truncating at $(trunc_val)/$(length(sv_cumsum)) retaining $(100.0*sv_cumsum[trunc_val])% of the variance of the structure matrix" + else + trunc_val = rk + if rk < size(data, 1) + @info " truncating at $(trunc_val)/$(size(data,1)), as low-rank data detected" + end + end + + sqrt_inv_sv = Diagonal(1.0 ./ sqrt.(svdA.S[1:trunc_val])) + sqrt_sv = Diagonal(sqrt.(svdA.S[1:trunc_val])) + # as we have svd of cov-matrix we can use U or Vt + encoder_mat = sqrt_inv_sv * svdA.Vt[1:trunc_val, :] + decoder_mat = svdA.Vt[1:trunc_val, :]' * sqrt_sv + + push!(get_encoder_mat(dd), encoder_mat) + push!(get_decoder_mat(dd), decoder_mat) + end +end + + +""" +$(TYPEDSIGNATURES) + +Apply the `Decorrelator` encoder, on a columns-are-data matrix +""" +function encode_data(dd::Decorrelator, data::MM) where {MM <: AbstractMatrix} + data_mean = get_data_mean(dd)[1] + encoder_mat = get_encoder_mat(dd)[1] + return encoder_mat * (data .- data_mean) +end + +""" +$(TYPEDSIGNATURES) + +Apply the `Decorrelator` decoder, on a columns-are-data matrix +""" +function decode_data(dd::Decorrelator, data::MM) where {MM <: AbstractMatrix} + data_mean = get_data_mean(dd)[1] + decoder_mat = get_decoder_mat(dd)[1] + return decoder_mat * data .+ data_mean +end + +""" +$(TYPEDSIGNATURES) + +Apply the `Decorrelator` encoder to a provided structure matrix +""" +function encode_structure_matrix(dd::Decorrelator, structure_matrix::MM) where {MM <: AbstractMatrix} + encoder_mat = get_encoder_mat(dd)[1] + return encoder_mat * structure_matrix * encoder_mat' +end + +""" +$(TYPEDSIGNATURES) + +Apply the `Decorrelator` decoder to a provided structure matrix +""" +function decode_structure_matrix(dd::Decorrelator, enc_structure_matrix::MM) where {MM <: AbstractMatrix} + decoder_mat = get_decoder_mat(dd)[1] + return decoder_mat * enc_structure_matrix * decoder_mat' +end diff --git a/src/Utilities/elementwise_scaler.jl b/src/Utilities/elementwise_scaler.jl new file mode 100644 index 000000000..c8b84454b --- /dev/null +++ b/src/Utilities/elementwise_scaler.jl @@ -0,0 +1,179 @@ +# included in Utilities.jl + +export UnivariateAffineScaling, ElementwiseScaler, QuartileScaling, MinMaxScaling, ZScoreScaling +export quartile_scale, minmax_scale, zscore_scale +export get_type, get_shift, get_scale + +""" +$(TYPEDEF) + +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`, +- [`minmax_scale`](@ref) : creates `MinMaxScaling` +- [`zscore_scale`](@ref) : creates `ZScoreScaling` + +and are accessed with [`get_type`](@ref) +""" +struct ElementwiseScaler{T, VV <: AbstractVector} <: DataContainerProcessor + "storage for the shift applied to data" + shift::VV + "storage for the scaling" + scale::VV +end + +abstract type UnivariateAffineScaling end +abstract type QuartileScaling <: UnivariateAffineScaling end +abstract type MinMaxScaling <: UnivariateAffineScaling end +abstract type ZScoreScaling <: UnivariateAffineScaling end + +""" +$(TYPEDSIGNATURES) + +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() = ElementwiseScaler(QuartileScaling) + +""" +$(TYPEDSIGNATURES) + +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() = ElementwiseScaler(MinMaxScaling) + +""" +$(TYPEDSIGNATURES) + +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 [`Decorrelator`](@ref) +""" +zscore_scale() = ElementwiseScaler(ZScoreScaling) + +ElementwiseScaler(::Type{UAS}) where {UAS <: UnivariateAffineScaling} = + ElementwiseScaler{UAS, Vector{Float64}}(Float64[], Float64[]) + +""" +$(TYPEDSIGNATURES) + +Gets the UnivariateAffineScaling type `T` +""" +get_type(es::ElementwiseScaler{T}) where {T} = T + +""" +$(TYPEDSIGNATURES) + +Gets the `shift` field of the `ElementwiseScaler` +""" +get_shift(es::ElementwiseScaler) = es.shift + +""" +$(TYPEDSIGNATURES) + +Gets the `scale` field of the `ElementwiseScaler` +""" +get_scale(es::ElementwiseScaler) = es.scale + +function Base.show(io::IO, es::ElementwiseScaler) + out = "ElementwiseScaler: $(get_type(es))" + print(io, out) +end + +function initialize_processor!( + 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(es), quartiles_mat[2, :]) + append!(get_scale(es), (quartiles_mat[3, :] - quartiles_mat[1, :])) +end + +function initialize_processor!( + 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(es), minmax_mat[1, :]) + append!(get_scale(es), (minmax_mat[2, :] - minmax_mat[1, :])) +end + +function initialize_processor!( + 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(es), stat_mat[1, :]) + append!(get_scale(es), stat_mat[2, :]) +end + +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 `ElementwiseScaler` encoder, on a columns-are-data matrix +""" +function encode_data(es::ElementwiseScaler, data::MM) where {MM <: AbstractMatrix} + out = deepcopy(data) + for i in 1:size(out, 1) + out[i, :] .-= get_shift(es)[i] + out[i, :] /= get_scale(es)[i] + end + return out +end + +""" +$(TYPEDSIGNATURES) + +Apply the `ElementwiseScaler` decoder, on a columns-are-data matrix +""" +function decode_data(es::ElementwiseScaler, data::MM) where {MM <: AbstractMatrix} + out = deepcopy(data) + for i in 1:size(out, 1) + out[i, :] *= get_scale(es)[i] + out[i, :] .+= get_shift(es)[i] + end + return out +end + +""" +$(TYPEDSIGNATURES) + +Computes and populates the `shift` and `scale` fields for the `ElementwiseScaler` +""" +initialize_processor!(es::ElementwiseScaler, data::MM, structure_matrix) where {MM <: AbstractMatrix} = + initialize_processor!(es, data) + + +""" +$(TYPEDSIGNATURES) + +Apply the `ElementwiseScaler` encoder to a provided structure matrix +""" +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 `ElementwiseScaler` decoder to a provided structure matrix +""" +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 9e560ac6b..aebcc63b6 100644 --- a/test/Utilities/runtests.jl +++ b/test/Utilities/runtests.jl @@ -72,7 +72,7 @@ end 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") + DD = Decorrelator([1], [2], [3], 1.0, "test") @test get_data_mean(DD) == [1] @test get_encoder_mat(DD) == [2] @test get_decoder_mat(DD) == [3] @@ -297,7 +297,6 @@ end 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") From 2110fea9ff42ee59a9d7b6c67093b2022b8fd9bb Mon Sep 17 00:00:00 2001 From: Arne Bouillon Date: Mon, 7 Jul 2025 18:19:00 -0700 Subject: [PATCH 2/5] Pass both structure matrices to PDCPs --- src/Utilities.jl | 12 +++++++----- src/Utilities/canonical_correlation.jl | 3 ++- 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/src/Utilities.jl b/src/Utilities.jl index 1272b21ae..7156af1cb 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -95,16 +95,18 @@ function _decode_data(proc::P, data, apply_to::AS) where {P <: DataProcessor, AS end function _initialize_and_encode_data!(proc::P, data, structure_mats, apply_to::AS) where {P <: DataProcessor, AS <: AbstractString} - input_data, output_data = data - input_structure_mat, output_structure_mat = structure_mats - if P isa PairedDataContainerProcessor - initialize_processor!(proc, input_data, output_data, apply_to == "in" ? input_structure_mat : output_structure_mat, apply_to) + initialize_processor!(proc, data..., structure_mats..., apply_to) else + input_data, output_data = data + input_structure_mat, output_structure_mat = structure_mats + if apply_to == "in" initialize_processor!(proc, input_data, input_structure_mat) - else + elseif apply_to == "out" initialize_processor!(proc, output_data, output_structure_mat) + else + bad_apply_to(apply_to) end end diff --git a/src/Utilities/canonical_correlation.jl b/src/Utilities/canonical_correlation.jl index 5ee1426c0..de08b4db2 100644 --- a/src/Utilities/canonical_correlation.jl +++ b/src/Utilities/canonical_correlation.jl @@ -170,7 +170,8 @@ initialize_processor!( cc::CanonicalCorrelation, in_data::MM, out_data::MM, - structure_matrix, + input_structure_matrix, + output_structure_matrix, apply_to::AS, ) where {MM <: AbstractMatrix, AS <: AbstractString} = initialize_processor!(cc, in_data, out_data, apply_to) From 01112fea757624da0018740f39ea10ec41160974 Mon Sep 17 00:00:00 2001 From: Arne Bouillon Date: Mon, 7 Jul 2025 18:31:33 -0700 Subject: [PATCH 3/5] Format code --- src/Utilities.jl | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/src/Utilities.jl b/src/Utilities.jl index 7156af1cb..53a23cfc7 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -94,7 +94,12 @@ function _decode_data(proc::P, data, apply_to::AS) where {P <: DataProcessor, AS end end -function _initialize_and_encode_data!(proc::P, data, structure_mats, apply_to::AS) where {P <: DataProcessor, AS <: AbstractString} +function _initialize_and_encode_data!( + proc::P, + data, + structure_mats, + apply_to::AS, +) where {P <: DataProcessor, AS <: AbstractString} if P isa PairedDataContainerProcessor initialize_processor!(proc, data..., structure_mats..., apply_to) else @@ -188,7 +193,12 @@ function initialize_and_encode_with_schedule!( for (processor, apply_to) in encoder_schedule @info "Initialize encoding of data: \"$(apply_to)\" with $(processor)" - processed = _initialize_and_encode_data!(processor, processed_io_pairs, (processed_input_structure_mat, processed_output_structure_mat), apply_to) + processed = _initialize_and_encode_data!( + processor, + processed_io_pairs, + (processed_input_structure_mat, processed_output_structure_mat), + apply_to, + ) if apply_to == "in" processed_input_structure_mat = encode_structure_matrix(processor, processed_input_structure_mat) From 90e26e30b44c250e58c60657ba8c7c39846c5660 Mon Sep 17 00:00:00 2001 From: Arne Bouillon Date: Wed, 9 Jul 2025 11:33:59 -0700 Subject: [PATCH 4/5] Fix bugs --- src/Utilities.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Utilities.jl b/src/Utilities.jl index 53a23cfc7..a6cfc1657 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -72,7 +72,7 @@ Base.:(==)(a::PDCP, b::PDCP) where {PDCP <: PairedDataContainerProcessor} = function _encode_data(proc::P, data, apply_to::AS) where {P <: DataProcessor, AS <: AbstractString} - input_data, output_data = data + input_data, output_data = get_data(data) if apply_to == "in" return encode_data(proc, input_data) elseif apply_to == "out" @@ -83,7 +83,7 @@ function _encode_data(proc::P, data, apply_to::AS) where {P <: DataProcessor, AS end function _decode_data(proc::P, data, apply_to::AS) where {P <: DataProcessor, AS <: AbstractString} - input_data, output_data = data + input_data, output_data = get_data(data) if apply_to == "in" return decode_data(proc, input_data) @@ -100,10 +100,10 @@ function _initialize_and_encode_data!( structure_mats, apply_to::AS, ) where {P <: DataProcessor, AS <: AbstractString} - if P isa PairedDataContainerProcessor + if proc isa PairedDataContainerProcessor initialize_processor!(proc, data..., structure_mats..., apply_to) else - input_data, output_data = data + input_data, output_data = get_data(data) input_structure_mat, output_structure_mat = structure_mats if apply_to == "in" From 8d33f8fd31812c0060974b4ec5f9e1f6153225cf Mon Sep 17 00:00:00 2001 From: Arne Bouillon Date: Wed, 9 Jul 2025 11:34:26 -0700 Subject: [PATCH 5/5] Adapt tests to new functions --- src/Emulator.jl | 7 ++++++- src/Utilities.jl | 13 ++++++++++--- test/Emulator/runtests.jl | 12 ++++++------ test/Utilities/runtests.jl | 15 +++++++++------ 4 files changed, 31 insertions(+), 16 deletions(-) diff --git a/src/Emulator.jl b/src/Emulator.jl index 652f984f8..288ed472f 100644 --- a/src/Emulator.jl +++ b/src/Emulator.jl @@ -174,7 +174,12 @@ function Emulator( enc_schedule = create_encoder_schedule(encoder_schedule) (encoded_io_pairs, encoded_input_structure_mat, encoded_output_structure_mat) = - encode_with_schedule!(enc_schedule, input_output_pairs, input_structure_mat, output_structure_mat) + initialize_and_encode_with_schedule!( + enc_schedule, + input_output_pairs, + input_structure_mat, + output_structure_mat, + ) # build the machine learning tool in the encoded space build_models!( diff --git a/src/Utilities.jl b/src/Utilities.jl index a6cfc1657..f8c1ae154 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -14,7 +14,14 @@ using ..DataContainers export get_training_points export PairedDataContainerProcessor, DataContainerProcessor -export create_encoder_schedule, initialize_and_encode_with_schedule!, encode_with_schedule, decode_with_schedule +export create_encoder_schedule, + initialize_and_encode_with_schedule!, + encode_with_schedule, + decode_with_schedule, + encode_data, + encode_structure_matrix, + decode_data, + decode_structure_matrix @@ -101,7 +108,7 @@ function _initialize_and_encode_data!( apply_to::AS, ) where {P <: DataProcessor, AS <: AbstractString} if proc isa PairedDataContainerProcessor - initialize_processor!(proc, data..., structure_mats..., apply_to) + initialize_processor!(proc, get_data(data)..., structure_mats..., apply_to) else input_data, output_data = get_data(data) input_structure_mat, output_structure_mat = structure_mats @@ -250,7 +257,7 @@ function encode_with_schedule( structure_matrix::USorM, in_or_out::AS, ) where {VV <: AbstractVector, USorM <: Union{UniformScaling, AbstractMatrix}, AS <: AbstractString} - if !(in_or_out ∈ ["in", "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/Emulator/runtests.jl b/test/Emulator/runtests.jl index 5c97c0033..5fae44e86 100644 --- a/test/Emulator/runtests.jl +++ b/test/Emulator/runtests.jl @@ -44,11 +44,11 @@ struct MLTester <: Emulators.MachineLearningTool end @test get_io_pairs(em) == io_pairs default_encoder = (decorrelate_sample_cov(), "in_and_out") # for these inputs this is the default enc_sch = create_encoder_schedule(default_encoder) - enc_io_pairs, enc_I_in, enc_I_out = encode_with_schedule!(enc_sch, io_pairs, 1.0 * I(p), 1.0 * I(d)) + enc_io_pairs, enc_I_in, enc_I_out = initialize_and_encode_with_schedule!(enc_sch, io_pairs, 1.0 * I(p), 1.0 * I(d)) @test get_encoder_schedule(em)[1][1] == enc_sch[1][1] # inputs: proc - @test get_encoder_schedule(em)[1][3] == enc_sch[1][3] # inputs: apply_to + @test get_encoder_schedule(em)[1][2] == enc_sch[1][2] # inputs: apply_to @test get_encoder_schedule(em)[2][1] == enc_sch[2][1] # outputs... - @test get_encoder_schedule(em)[2][3] == enc_sch[2][3] + @test get_encoder_schedule(em)[2][2] == enc_sch[2][2] @test get_data(get_encoded_io_pairs(em)) == get_data(enc_io_pairs) @test get_data(get_encoded_io_pairs(em)) == get_data(enc_io_pairs) @@ -91,7 +91,7 @@ struct MLTester <: Emulators.MachineLearningTool end enc_sch1 = create_encoder_schedule([(decorrelate_sample_cov(), "in"), (decorrelate_structure_mat(), "out")]) enc_sch2 = create_encoder_schedule([(decorrelate_structure_mat(), "in"), (decorrelate_structure_mat(), "out")]) enc_sch3 = create_encoder_schedule([(decorrelate_structure_mat(), "in"), (decorrelate_sample_cov(), "out")]) - (_, _, _) = encode_with_schedule!( + (_, _, _) = initialize_and_encode_with_schedule!( enc_sch1, io_pairs, 1.0 * I(p), @@ -99,7 +99,7 @@ struct MLTester <: Emulators.MachineLearningTool end ) @test get_encoder_schedule(em1) == enc_sch1 - (_, _, _) = encode_with_schedule!( + (_, _, _) = initialize_and_encode_with_schedule!( enc_sch2, io_pairs, 4.0 * I(p), @@ -107,7 +107,7 @@ struct MLTester <: Emulators.MachineLearningTool end ) @test get_encoder_schedule(em2) == enc_sch2 - (_, _, _) = encode_with_schedule!(enc_sch3, io_pairs, Γ, I(d)) + (_, _, _) = initialize_and_encode_with_schedule!(enc_sch3, io_pairs, Γ, I(d)) @test get_encoder_schedule(em3) == enc_sch3 end diff --git a/test/Utilities/runtests.jl b/test/Utilities/runtests.jl index aebcc63b6..1d5250197 100644 --- a/test/Utilities/runtests.jl +++ b/test/Utilities/runtests.jl @@ -143,7 +143,7 @@ end 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) + initialize_and_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) @@ -293,18 +293,21 @@ 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) + bad_encoder_schedule = [(canonical_correlation(), "bad")] + @test_throws ArgumentError initialize_and_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 decode_data(canonical_correlation(), func(io_pairs), "bad") encoder_schedule = create_encoder_schedule(schedule_builder) # encode the data using the schedule (encoded_io_pairs, encoded_prior_cov, encoded_obs_noise_cov) = - encode_with_schedule!(encoder_schedule, io_pairs, prior_cov, obs_noise_cov) + initialize_and_encode_with_schedule!(encoder_schedule, io_pairs, prior_cov, obs_noise_cov) # decode the data using the schedule (decoded_io_pairs, decoded_prior_cov, decoded_obs_noise_cov) =