adding chaotic and chebyshev inits

Author: MartinuzziFrancesco

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

src/esn/esn_inits.jl

@@ -282,6 +282,68 @@ 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 `true`.
+
+[^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 = true) 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 +734,59 @@ 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
+
+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 "Using dims[1] = $requested_order for the chaotic reservoir matrix order."
+    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 "Chaotic reservoir matrix order adjusted from $requested_order to $actual_matrix_order based on computed bit precision = $chosen_bit_precision."
+    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...)