added tests for inits

Author: MartinuzziFrancesco

src/ReservoirComputing.jl

@@ -40,7 +40,7 @@ export StandardStates, ExtendedStates, PaddedStates, PaddedExtendedStates
 export StandardRidge
 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
+       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,9 +282,7 @@ 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)
@@ -297,21 +295,24 @@ via the Chebyshev mapping. The first row is defined as:
     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).  
+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.  
+# 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
@@ -334,7 +335,6 @@ julia> input_matrix = chebyshev_mapping(10, 3)
  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.
@@ -343,25 +343,25 @@ julia> input_matrix = chebyshev_mapping(10, 3)
     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} 
+        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]))
+            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
 
 """
@@ -752,22 +752,21 @@ 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
+        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
 
 """
@@ -779,17 +778,19 @@ 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 
+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
@@ -812,8 +813,7 @@ julia> res_matrix = chaotic_init(8, 8)
    ⋅        -0.600945   ⋅          ⋅ 
    ⋅          ⋅        0.132667   2.21354
    ⋅        -2.60383    ⋅        -2.90391
- -0.578156    ⋅         ⋅          ⋅ 
-
+ -0.578156    ⋅         ⋅          ⋅
 ```
 
 [^xie2024]: Xie, Minzhi, Qianxue Wang, and Simin Yu.
@@ -822,9 +822,8 @@ julia> res_matrix = chaotic_init(8, 8)
     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}
-
+        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
@@ -833,10 +832,11 @@ function chaotic_init(rng::AbstractRNG, ::Type{T}, dims::Integer...;
     end
     d_estimate = log2(requested_order) / 2
     d_floor = max(floor(Int, d_estimate), 1)
-    d_ceil  = ceil(Int, d_estimate)
+    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
+    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
@@ -846,8 +846,10 @@ function chaotic_init(rng::AbstractRNG, ::Type{T}, dims::Integer...;
         """
     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)    
+    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

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