diff --git a/docs/make.jl b/docs/make.jl index dcc31176..8dd919c0 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -36,6 +36,7 @@ makedocs( ], "Approximators" => [ "Approximators API" => "api/approximators.md", + "Selectors API" => "api/selectors.md", "ApproximatorErrors API" => "api/approximator_errors.md", ], ], diff --git a/docs/src/api/selectors.md b/docs/src/api/selectors.md new file mode 100644 index 00000000..81458dd6 --- /dev/null +++ b/docs/src/api/selectors.md @@ -0,0 +1,33 @@ +# Selectors +```@contents +Pages = ["selectors.md"] +``` + +## Abstract Types +```@docs +Selector + +SelectorRecipe +``` + +## Selector Structures +```@docs +RLinearAlgebra.LUPP + +LUPPRecipe + +``` + +## Exported Functions +```@docs +complete_selector + +update_selector! + +select_indices! +``` + +## Internal Functions +```@docs + +``` diff --git a/src/Approximators.jl b/src/Approximators.jl index 7b8b61e9..441df78d 100644 --- a/src/Approximators.jl +++ b/src/Approximators.jl @@ -56,6 +56,7 @@ randomized rangefinder. """ abstract type RangeApproximatorRecipe <: ApproximatorRecipe end + ################################### # Docstring Components ################################### @@ -72,6 +73,7 @@ approx_arg_list = Dict{Symbol,String}( :A => "`A::AbstractMatrix`, a target matrix for approximation.", :compressor_recipe => "`S::CompressorRecipe`, a fully initialized realization for a compression method for a specific matrix or collection of matrices and vectors.", + ) approx_output_list = Dict{Symbol,String}( @@ -82,8 +84,10 @@ approx_output_list = Dict{Symbol,String}( approx_method_description = Dict{Symbol,String}( :complete_approximator => "A function that generates an `ApproximatorRecipe` given arguments.", + :update_approximator => "A function that updates the `ApproximatorRecipe` in place given the arguments.", + :rapproximate => "A function that computes a low-rank approximation of the matrix `A` using the information in the provided `Approximator` data structure.", :complete_approximator_error => "A function that generates an `ApproximatorErrorRecipe` @@ -161,6 +165,7 @@ function complete_approximator(approximator::Approximator, A::AbstractMatrix) ) end + ################################### # rapproximate Interface ################################### @@ -430,3 +435,8 @@ end include("Approximators/RangeApproximators/rangefinder.jl") include("Approximators/RangeApproximators/randsvd.jl") include("Approximators/RangeApproximators/helpers/power_its.jl") + +############################################ +# Include the selector files +############################################ +include("Approximators/Selectors.jl") diff --git a/src/Approximators/Selectors.jl b/src/Approximators/Selectors.jl new file mode 100644 index 00000000..7e656b3c --- /dev/null +++ b/src/Approximators/Selectors.jl @@ -0,0 +1,147 @@ +################################### +# Abstract Types +################################### +""" + Selector + +An abstract type containing user controlled parameters for a technique that selects indices +from a matrix. +""" +abstract type Selector end + +""" + SelectorRecipe + +An abstract type containing user controlled parameters and preallocated memory for a +technique that selects indices from a matrix. +""" +abstract type SelectorRecipe end + +select_arg_list = Dict{Symbol,String}( + :selector => "`selector::Selector`, a data structure containing the user-defined + parameters associated with a particular selection method.", + :selector_recipe => "`selector::SelectorRecipe`, a fully initialized realization + for a selector method for a particular matrix.", + :idx => "`idx::vector`, a vector where selected indices will be placed.", + :n_idx => "`n_idx::Int64`, the number of indices to be selected.", + :start_idx => "`start_idx::Int64`. the starting location in `idx` where the indices + will be placed.", + :A => "`A::AbstractMatrix`, a target matrix for approximation.", +) +select_output_list = Dict{Symbol,String}( + :selector_recipe => "A `SelectorRecipe` object.", +) +select_method_description = Dict{Symbol,String}( + :update_selector => "A function that updates the `SelectorRecipe` in place given the + arguments.", + :complete_selector => "A function that generates a `SelectorRecipe` given + arguments.", + :select_indices => "A function that selects indices from a matrix `A` using a specific + `SelectorRecipe`. It updates the vector `idx` in place with `n_idx` new indices starting + at index `start_idx`." +) +################################## +# Complete Selector Interface +################################## +""" + complete_selector(selector::Selector, A::AbstractMatrix) + +$(select_method_description[:complete_selector]) + +# Arguments +- $(select_arg_list[:selector]) +- $(select_arg_list[:A]) + +# Outputs +- $(select_output_list[:selector_recipe]) +""" +function complete_selector(selector::Selector, A::AbstractMatrix) + return throw( + ArgumentError( + "No method `complete_selector` exists for selector of type\ + $(typeof(selector)) and matrix of type $(typeof(A))." + ) + ) +end + +################################## +# update_selector! +################################## +""" + update_selector!(selector::SelectorRecipe) + +$(select_method_description[:update_selector]) + +# Arguments +- $(select_arg_list[:selector_recipe]) + +# Outputs +- $(select_output_list[:selector_recipe]) +""" +function update_selector!(selector::SelectorRecipe) + return throw( + ArgumentError( + "No method `update_selector!` exists for selector of type\ + $(typeof(selector))." + ) + ) +end + +""" + update_selector!(selector::SelectorRecipe, A::AbstractMatrix) + +$(select_method_description[:update_selector]) + +# Arguments +- $(select_arg_list[:selector_recipe]) +- $(select_arg_list[:A]) + +# Outputs +- $(select_output_list[:selector_recipe]) +""" +function update_selector!(selector::SelectorRecipe, A::AbstractMatrix) + return update_selector!(selector) +end + +#################################### +# select_indices! +#################################### +""" + select_indices!( + idx::AbstractVector, + selector::SelectorRecipe, + A::AbstractMatrix, + n_idx::Int64, + start_idx::Int64 + ) + +$(select_method_description[:select_indices]) + +# Arguments +- $(select_arg_list[:idx]) +- $(select_arg_list[:selector_recipe]) +- $(select_arg_list[:A]) +- $(select_arg_list[:n_idx]) +- $(select_arg_list[:start_idx]) + +# Outputs +- Returns `nothing` +""" +function select_indices!( + idx::AbstractVector, + selector::SelectorRecipe, + A::AbstractMatrix, + n_idx::Int64, + start_idx::Int64 +) + return throw( + ArgumentError( + "No method `select_indices` exists for selector of type $(typeof(selector)),\ + matrix of type $(typeof(A)), idx of type $(typeof(idx)), n_idx of type \ + $(typeof(n_idx)), and start_idx of type $(typeof(start_idx))." + ) + ) +end + +# Include the selector files +include("Selectors/lupp.jl") diff --git a/src/Approximators/Selectors/lupp.jl b/src/Approximators/Selectors/lupp.jl new file mode 100644 index 00000000..1b8cd789 --- /dev/null +++ b/src/Approximators/Selectors/lupp.jl @@ -0,0 +1,110 @@ +""" + LUPP <: Selector + +A `Selector` that implements LU with partial pivoting for selecting column indices from a +matrix. + +# Fields +- `compressor::Compressor`, the compression technique that will be applied to the matrix, + before selecting indices. + +# Constructor + LUPP(;compressor = Identity()) + +## Keywords +- `compressor::Compressor`, the compression technique that will be applied to the matrix, + before selecting indices. Defaults to the `Identity` compressor. + +## Returns +- A `LUPP` object. + +!!! note "Implementation Note" + LU with partial pivoting is classically implemented to select rows of a matrix. Here we + apply LU with partial pivoting to the transpose of the inputted matrix to select + columns. +""" +mutable struct LUPP <: Selector + compressor::Compressor +end + +function LUPP(;compressor = Identity()) + return LUPP(compressor) +end + +""" + LUPPRecipe <: SelectorRecipe + +A `SelectorRecipe` that contains all the necessary preallocations for selecting column +indices from a matrix using LU with partial pivoting. + + +# Fields +- `compressor::CompressorRecipe`, the compression technique that will applied to the matrix, + before selecting indices. +- `SA::AbstractMatrix`, a buffer matrix for storing the sketched matrix. + +!!! note "Implementation Note" + LU with partial pivoting is classically implemented to select rows of a matrix. Here we + apply LU with partial pivoting to the transpose of the inputted matrix to select + columns. +""" +mutable struct LUPPRecipe <: SelectorRecipe + compressor::CompressorRecipe + SA::AbstractMatrix +end + +function complete_selector(ingredients::LUPP, A::AbstractMatrix) + compressor = complete_compressor(ingredients.compressor, A) + n_rows, n_cols = size(compressor) + SA = Matrix{eltype(A)}(undef, n_rows, n_cols) + return LUPPRecipe(compressor, SA) +end + +function update_selector!(selector::LUPPRecipe) + update_compressor!(selector.compressor) + return nothing +end + +function select_indices!( + idx::AbstractVector, + selector::LUPPRecipe, + A::AbstractMatrix, + n_idx::Int64, + start_idx::Int64 +) + # you cannot select more column indices than there are columns in the matrix + if n_idx > size(A, 2) + throw( + DimensionMismatch( + "`n_idx` cannot be larger than the number of columns in `A`." + ) + ) + end + + # start_idx + n_idx must be less than the length of the idx vector + if start_idx + n_idx - 1 > size(idx, 1) + throw( + DimensionMismatch( + "`start_idx` + `n_idx` - 1 cannot be larger than the lenght of `idx`." + ) + ) + end + + # you cannot select more indices than the compression dimension because that is when + # LUPP will stop selecting new pivots because the LU factorization will have been formed + if n_idx > selector.compressor.n_rows + throw( + DimensionMismatch( + "Must select fewer indices then the `compression_dim`." + ) + ) + + end + + mul!(selector.SA, selector.compressor, A) + # because LUPP selects rows and selectors select columns we need to pivot on A' + p = lu!(selector.SA').p + # store n_idx indices in the appropriate part of the idx + idx[start_idx:start_idx + n_idx - 1] = p[1:n_idx] + return nothing +end diff --git a/src/RLinearAlgebra.jl b/src/RLinearAlgebra.jl index b93cfb51..4addd7c8 100644 --- a/src/RLinearAlgebra.jl +++ b/src/RLinearAlgebra.jl @@ -1,7 +1,8 @@ module RLinearAlgebra import Base.:* import Base: transpose, adjoint -import LinearAlgebra: Adjoint, axpby!, dot, I, ldiv!, lmul!, lq!, lq, LQ, mul!, norm, qr!, svd +import LinearAlgebra: Adjoint, axpby!, dot, I, ldiv!, lmul!, lq!, lq, LQ, lu! +import LinearAlgebra: mul!, norm, qr!, svd import StatsBase: sample, sample!, ProbabilityWeights, wsample! import Random: bitrand, rand!, randn! import SparseArrays: SparseMatrixCSC, sprandn, sparse @@ -62,4 +63,8 @@ export FullResidual, FullResidualRecipe export ApproximatorError, ApproximatorErrorRecipe export complete_approximator_error, compute_approximator_error, compute_approximator_error! +# Export Selector types and functions +export Selector, SelectorRecipe +export LUPP, LUPPRecipe +export complete_selector, update_selector!, select_indices! end #module diff --git a/test/Approximators/Selectors/lupp.jl b/test/Approximators/Selectors/lupp.jl new file mode 100644 index 00000000..9cb0a310 --- /dev/null +++ b/test/Approximators/Selectors/lupp.jl @@ -0,0 +1,185 @@ +module LUPP_tests +using Test +using ..FieldTest +using ..ApproxTol +using RLinearAlgebra +import LinearAlgebra: mul! + +@testset "LUPP Tests" begin + @testset "LUPP" begin + @test supertype(LUPP) == Selector + + # test fieldnames and types + @test fieldnames(LUPP) == (:compressor,) + @test fieldtypes(LUPP) == (Compressor,) + + # Test Constructor + let + sel = LUPP() + @test typeof(sel) == LUPP + @test typeof(sel.compressor) == Identity + end + + let + sel = LUPP(compressor = SparseSign()) + @test typeof(sel) == LUPP + @test typeof(sel.compressor) == SparseSign + end + + end + + @testset "LUPPRecipe" begin + @test supertype(LUPPRecipe) == SelectorRecipe + # test fieldnames and types + @test fieldnames(LUPPRecipe) == (:compressor, :SA) + @test fieldtypes(LUPPRecipe) == (CompressorRecipe, AbstractMatrix) + end + + @testset "LUPP: Complete Selector" begin + # test with identity compressor + let n_rows = 2, + n_cols = 2, + A = zeros(n_rows, n_cols), + sel = LUPP() + + sel_rec = complete_selector(sel, A) + @test typeof(sel_rec) == LUPPRecipe + @test typeof(sel_rec.compressor) == IdentityRecipe{Left} + @test typeof(sel_rec.SA) <: AbstractMatrix + @test size(sel_rec.SA) == (n_rows, n_cols) + end + + # test with gaussian compressor + let n_rows = 3, + n_cols = 3, + A = zeros(n_rows, n_cols), + comp_dim = 2, + sel = LUPP(compressor = Gaussian(compression_dim = comp_dim)) + + sel_rec = complete_selector(sel, A) + @test typeof(sel_rec) == LUPPRecipe + @test typeof(sel_rec.compressor) == GaussianRecipe{Left} + @test typeof(sel_rec.SA) <: AbstractMatrix + @test size(sel_rec.SA) == (comp_dim, n_cols) + end + + end + + @testset "LUPP: Update Selector" begin + let n_rows = 2, + n_cols = 2, + A = zeros(n_rows, n_cols), + sel_rec = complete_selector(LUPP(compressor = Gaussian()), A) + G_old = deepcopy(sel_rec.compressor.op) + update_selector!(sel_rec) + @test typeof(sel_rec) == LUPPRecipe + # check that Gaussian Matrix actually changes + @test sel_rec.compressor.op != G_old + end + + end + + @testset "LUPP: Select Indices" begin + A = [0.0 10.0 0.0; + 0.0 0.0 20.0; + 30.0 0.0 0.0] + # test the error checking + # start with checking n_idx not being larger than number of columns + let A = A, + idx = zeros(Int64, 3), + start_idx = 1, + n_idx = 4 + + @test_throws DimensionMismatch select_indices!( + idx, + complete_selector(LUPP(), A), + A, + n_idx, + start_idx + ) + end + + # check that n_idx will not go over index vector + let A = A, + idx = zeros(Int64, 3), + start_idx = 3, + n_idx = 2 + + @test_throws DimensionMismatch select_indices!( + idx, + complete_selector(LUPP(), A), + A, + n_idx, + start_idx + ) + end + + # check that n_idx is not larger than the compression_dim + let A = A, + idx = zeros(Int64, 3), + start_idx = 1, + n_idx = 3 + + @test_throws DimensionMismatch select_indices!( + idx, + complete_selector( + LUPP(compressor = Gaussian(compression_dim = 2)), + A + ), + A, + n_idx, + start_idx + ) + end + + # Test with identity compressor + # notice that in this selection it has nothing to do with column norm + # test selecting only one index + let A = deepcopy(A), + idx = zeros(Int64, 3), + start_idx = 2, + n_idx = 1, + sel_rec = complete_selector(LUPP(), A) + + select_indices!(idx, sel_rec, A, n_idx, start_idx) + @test idx[2] == 2 + end + + # test selecting three indices + let A = deepcopy(A), + idx = zeros(Int64, 3), + start_idx = 1, + n_idx = 3, + sel_rec = complete_selector(LUPP(), A) + + select_indices!(idx, sel_rec, A, n_idx, start_idx) + @test idx == [2; 3; 1] + end + + # test with gaussian compressor + # test selecting only one index + let A = deepcopy(A), + idx = zeros(Int64, 3), + start_idx = 2, + n_idx = 1, + sel_rec = complete_selector(LUPP(compressor = Gaussian()), A) + + select_indices!(idx, sel_rec, A, n_idx, start_idx) + @test idx[2] == 1 + end + + # test selecting two indices + let A = deepcopy(A), + idx = zeros(Int64, 3), + start_idx = 1, + n_idx = 2, + sel_rec = complete_selector(LUPP(compressor = Gaussian()), A) + + select_indices!(idx, sel_rec, A, n_idx, start_idx) + @test idx == [1; 3; 0] + end + end + +end + +end diff --git a/test/Approximators/selector_abstract_types.jl b/test/Approximators/selector_abstract_types.jl new file mode 100644 index 00000000..6c4dd9e6 --- /dev/null +++ b/test/Approximators/selector_abstract_types.jl @@ -0,0 +1,101 @@ +module selector_abstract_types +using Test, RLinearAlgebra +import RLinearAlgebra: complete_selector, update_selector!, select_indices! + + +############################# +# Initial Testing Parameters +############################# +struct TestSelector <: Selector end +mutable struct TestSelectorRecipe <: SelectorRecipe + code::Int64 +end + +######################################## +# Tests for Selector Abstract Types +######################################## +@testset "Selector Abstract Types" begin + @test isdefined(Main, :Selector) + @test isdefined(Main, :SelectorRecipe) +end + +@testset "Selector Completion Errors" begin + A = ones(2,2) + + @test_throws ArgumentError complete_selector(TestSelector(), A) +end + +############################# +# Update Testing Parameters +############################# +complete_selector(selector::TestSelector, A::AbstractMatrix) = + TestSelectorRecipe(1) + +@testset "Selector Completion" begin + A = ones(2, 2) + + select = complete_selector(TestSelector(), A) + @test select isa TestSelectorRecipe + @test select.code == 1 +end + +############################### +# Test Updating selector errors +############################### + +@testset "Selector Update Errors" begin + A = ones(2,2) + + select = complete_selector(TestSelector(), A) + @test_throws ArgumentError update_selector!(select, A) + @test_throws ArgumentError update_selector!(select) + +end + +############################### +# Test update_selector +############################### + + +@testset "Selector Updating" begin + A = ones(2, 2) + function update_selector!(selector::TestSelectorRecipe, A::AbstractMatrix) + selector.code = 2 + end + + function update_selector!(selector::TestSelectorRecipe) + selector.code = 2 + end + + # test 2 arg selector + let A = A + select = complete_selector(TestSelector(), A) + update_selector!(select, A) + @test select isa TestSelectorRecipe + @test select.code == 2 + end + + # test 1 arg selector + let A = A + select = complete_selector(TestSelector(), A) + update_selector!(select) + @test select isa TestSelectorRecipe + @test select.code == 2 + end + +end + +################################# +# Test selection errors +################################# +@testset "Select Indices" begin + A = ones(2,2) + select = complete_selector(TestSelector(), A) + idx = zeros(Int64, 2) + n_idx = 2 + start_idx = 1 + + @test_throws ArgumentError select_indices!(idx, select, A, n_idx, start_idx) +end + +end diff --git a/test/runtests.jl b/test/runtests.jl index 31698e5e..e9c8017d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -16,6 +16,8 @@ directs = "Compressors/", "Approximators/", "Approximators/RangeApproximator/", + "Approximators/Selectors/", + "Solvers/", "Solvers/helpers/", "Solvers/ErrorMethods/", "Solvers/Loggers/",