diff --git a/Project.toml b/Project.toml index 773a172f..b629b063 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ReservoirComputing" uuid = "7c2d2b1e-3dd4-11ea-355a-8f6a8116e294" authors = ["Francesco Martinuzzi"] -version = "0.10.8" +version = "0.10.9" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" diff --git a/docs/src/api/inits.md b/docs/src/api/inits.md index dd423518..bd52797e 100644 --- a/docs/src/api/inits.md +++ b/docs/src/api/inits.md @@ -7,6 +7,7 @@ weighted_init informed_init minimal_init + chebyshev_mapping ``` ## Reservoirs @@ -18,4 +19,5 @@ cycle_jumps simple_cycle pseudo_svd + chaotic_init ``` diff --git a/src/ReservoirComputing.jl b/src/ReservoirComputing.jl index 46fd362b..61c505b8 100644 --- a/src/ReservoirComputing.jl +++ b/src/ReservoirComputing.jl @@ -38,8 +38,9 @@ include("reca/reca_input_encodings.jl") export NLADefault, NLAT1, NLAT2, NLAT3 export StandardStates, ExtendedStates, PaddedStates, PaddedExtendedStates export StandardRidge -export scaled_rand, weighted_init, informed_init, minimal_init -export rand_sparse, delay_line, delay_line_backward, cycle_jumps, simple_cycle, pseudo_svd +export scaled_rand, weighted_init, informed_init, minimal_init, chebyshev_mapping +export rand_sparse, delay_line, delay_line_backward, cycle_jumps, + simple_cycle, pseudo_svd, chaotic_init export RNN, MRNN, GRU, GRUParams, FullyGated, Minimal export train export ESN, HybridESN, KnowledgeModel, DeepESN diff --git a/src/esn/esn_inits.jl b/src/esn/esn_inits.jl index 79805381..f9ad66c2 100644 --- a/src/esn/esn_inits.jl +++ b/src/esn/esn_inits.jl @@ -282,6 +282,86 @@ function _create_irrational(irrational::Irrational, start::Int, res_size::Int, return T.(input_matrix) end +""" + chebyshev_mapping([rng], [T], dims...; + amplitude=one(T), sine_divisor=one(T), + chebyshev_parameter=one(T), return_sparse=true) + +Generate a Chebyshev-mapped matrix [^xie2024]. The first row is initialized +using a sine function and subsequent rows are iteratively generated +via the Chebyshev mapping. The first row is defined as: + +```math + w(1, j) = amplitude * sin(j * π / (sine_divisor * n_cols)) +``` + +for j = 1, 2, …, n_cols (with n_cols typically equal to K+1, where K is the number of input layer neurons). +Subsequent rows are generated by applying the mapping: + +```math + w(i+1, j) = cos(chebyshev_parameter * acos(w(i, j))) +``` + +# Arguments + + - `rng`: Random number generator. Default is `Utils.default_rng()` + from WeightInitializers. + - `T`: Type of the elements in the reservoir matrix. + Default is `Float32`. + - `dims`: Dimensions of the matrix. Should follow `res_size x in_size`. Here, res_size is assumed to be K+1. + +# Keyword arguments + + - `amplitude`: Scaling factor used to initialize the first row. + This parameter adjusts the amplitude of the sine function. Default value is one. + - `sine_divisor`: Divisor applied in the sine function's phase. Default value is one. + - `chebyshev_parameter`: Control parameter for the Chebyshev mapping in + subsequent rows. This parameter influences the distribution of the + matrix elements. Default is one. + - `return_sparse`: If `true`, the function returns the matrix as a sparse matrix. Default is `false`. + +# Examples + +```jldoctest +julia> input_matrix = chebyshev_mapping(10, 3) +10×3 Matrix{Float32}: + 0.866025 0.866025 1.22465f-16 + 0.866025 0.866025 -4.37114f-8 + 0.866025 0.866025 -4.37114f-8 + 0.866025 0.866025 -4.37114f-8 + 0.866025 0.866025 -4.37114f-8 + 0.866025 0.866025 -4.37114f-8 + 0.866025 0.866025 -4.37114f-8 + 0.866025 0.866025 -4.37114f-8 + 0.866025 0.866025 -4.37114f-8 + 0.866025 0.866025 -4.37114f-8 +``` + +[^xie2024]: Xie, Minzhi, Qianxue Wang, and Simin Yu. + "Time Series Prediction of ESN Based on Chebyshev Mapping and Strongly + Connected Topology." + Neural Processing Letters 56.1 (2024): 30. +""" +function chebyshev_mapping(rng::AbstractRNG, ::Type{T}, dims::Integer...; + amplitude::AbstractFloat=one(T), sine_divisor::AbstractFloat=one(T), + chebyshev_parameter::AbstractFloat=one(T), + return_sparse::Bool=false) where {T <: Number} + input_matrix = DeviceAgnostic.zeros(rng, T, dims...) + n_rows, n_cols = dims[1], dims[2] + + for idx_cols in 1:n_cols + input_matrix[1, idx_cols] = amplitude * sin(idx_cols * pi / (sine_divisor * n_cols)) + end + for idx_rows in 2:n_rows + for idx_cols in 1:n_cols + input_matrix[idx_rows, idx_cols] = cos(chebyshev_parameter * acos(input_matrix[ + idx_rows - 1, idx_cols])) + end + end + + return return_sparse ? sparse(input_matrix) : input_matrix +end + ### reservoirs """ @@ -672,11 +752,118 @@ function get_sparsity(M, dim) return size(M[M .!= 0], 1) / (dim * dim - size(M[M .!= 0], 1)) #nonzero/zero elements end +function digital_chaotic_adjacency(rng::AbstractRNG, bit_precision::Integer; + extra_edge_probability::AbstractFloat=0.1) + matrix_order = 2^(2 * bit_precision) + adjacency_matrix = zeros(Int, matrix_order, matrix_order) + for row_index in 1:(matrix_order - 1) + adjacency_matrix[row_index, row_index + 1] = 1 + end + adjacency_matrix[matrix_order, 1] = 1 + for row_index in 1:matrix_order, column_index in 1:matrix_order + if row_index != column_index && rand(rng) < extra_edge_probability + adjacency_matrix[row_index, column_index] = 1 + end + end + + return adjacency_matrix +end + +""" + chaotic_init([rng], [T], dims...; + extra_edge_probability=T(0.1), spectral_radius=one(T), + return_sparse_matrix=true) + +Construct a chaotic reservoir matrix using a digital chaotic system [^xie2024]. + +The matrix topology is derived from a strongly connected adjacency +matrix based on a digital chaotic system operating at finite precision. +If the requested matrix order does not exactly match a valid order the +closest valid order is used. + +# Arguments + + - `rng`: Random number generator. Default is `Utils.default_rng()` + from WeightInitializers. + - `T`: Type of the elements in the reservoir matrix. + Default is `Float32`. + - `dims`: Dimensions of the reservoir matrix. + +# Keyword arguments + + - `extra_edge_probability`: Probability of adding extra random edges in + the adjacency matrix to enhance connectivity. Default is 0.1. + - `desired_spectral_radius`: The target spectral radius for the + reservoir matrix. Default is one. + - `return_sparse_matrix`: If `true`, the function returns the + reservoir matrix as a sparse matrix. Default is `true`. + +# Examples + +```jldoctest +julia> res_matrix = chaotic_init(8, 8) +┌ Warning: +│ +│ Adjusting reservoir matrix order: +│ from 8 (requested) to 4 +│ based on computed bit precision = 1. +│ +└ @ ReservoirComputing ~/.julia/dev/ReservoirComputing/src/esn/esn_inits.jl:805 +4×4 SparseArrays.SparseMatrixCSC{Float32, Int64} with 6 stored entries: + ⋅ -0.600945 ⋅ ⋅ + ⋅ ⋅ 0.132667 2.21354 + ⋅ -2.60383 ⋅ -2.90391 + -0.578156 ⋅ ⋅ ⋅ +``` + +[^xie2024]: Xie, Minzhi, Qianxue Wang, and Simin Yu. + "Time Series Prediction of ESN Based on Chebyshev Mapping and Strongly + Connected Topology." + Neural Processing Letters 56.1 (2024): 30. +""" +function chaotic_init(rng::AbstractRNG, ::Type{T}, dims::Integer...; + extra_edge_probability::AbstractFloat=T(0.1), spectral_radius::AbstractFloat=one(T), + return_sparse_matrix::Bool=true) where {T <: Number} + requested_order = first(dims) + if length(dims) > 1 && dims[2] != requested_order + @warn """\n + Using dims[1] = $requested_order for the chaotic reservoir matrix order.\n + """ + end + d_estimate = log2(requested_order) / 2 + d_floor = max(floor(Int, d_estimate), 1) + d_ceil = ceil(Int, d_estimate) + candidate_order_floor = 2^(2 * d_floor) + candidate_order_ceil = 2^(2 * d_ceil) + chosen_bit_precision = abs(candidate_order_floor - requested_order) <= + abs(candidate_order_ceil - requested_order) ? d_floor : d_ceil + actual_matrix_order = 2^(2 * chosen_bit_precision) + if actual_matrix_order != requested_order + @warn """\n + Adjusting reservoir matrix order: + from $requested_order (requested) to $actual_matrix_order + based on computed bit precision = $chosen_bit_precision. \n + """ + end + + random_weight_matrix = T(2) * rand(rng, T, actual_matrix_order, actual_matrix_order) .- + T(1) + adjacency_matrix = digital_chaotic_adjacency( + rng, chosen_bit_precision; extra_edge_probability=extra_edge_probability) + reservoir_matrix = random_weight_matrix .* adjacency_matrix + current_spectral_radius = maximum(abs, eigvals(reservoir_matrix)) + if current_spectral_radius != 0 + reservoir_matrix .*= spectral_radius / current_spectral_radius + end + + return return_sparse_matrix ? sparse(reservoir_matrix) : reservoir_matrix +end + ### fallbacks #fallbacks for initializers #eventually to remove once migrated to WeightInitializers.jl for initializer in (:rand_sparse, :delay_line, :delay_line_backward, :cycle_jumps, - :simple_cycle, :pseudo_svd, - :scaled_rand, :weighted_init, :informed_init, :minimal_init) + :simple_cycle, :pseudo_svd, :chaotic_init, + :scaled_rand, :weighted_init, :informed_init, :minimal_init, :chebyshev_mapping) @eval begin function ($initializer)(dims::Integer...; kwargs...) return $initializer(Utils.default_rng(), Float32, dims...; kwargs...) diff --git a/test/esn/test_inits.jl b/test/esn/test_inits.jl index 3a337f42..944a4613 100644 --- a/test/esn/test_inits.jl +++ b/test/esn/test_inits.jl @@ -1,7 +1,7 @@ using ReservoirComputing, LinearAlgebra, Random, SparseArrays -const res_size = 30 -const in_size = 3 +const res_size = 16 +const in_size = 4 const radius = 1.0 const sparsity = 0.1 const weight = 0.2 @@ -24,13 +24,15 @@ reservoir_inits = [ delay_line_backward, cycle_jumps, simple_cycle, - pseudo_svd + pseudo_svd, + chaotic_init ] input_inits = [ scaled_rand, weighted_init, minimal_init, - minimal_init(; sampling_type=:irrational) + minimal_init(; sampling_type=:irrational), + chebyshev_mapping ] @testset "Reservoir Initializers" begin