Examples
In this section, we provide some examples where it is possible to apply RMA to analyze data.
Classifying data with a multi-layer perceptron and RMA
In this example, we demonstrate how to use a multi-layer perceptron with RMA to classify time series based on a parameter used to generate them. It is based on (Spezzatto et al., 2024), and we will use the package Flux.jl to perform the machine learning tasks.
Generating data
To exemplify this, we use the Lorenz system and define 5 classes, where each of them uses a different $\rho$ value to generate the time series. Our goal here is to classify the time series based on which $\rho$ value generated it. Thus, our classes are: ρ_cls = [26.0, 27.0, 28.0, 29.0, 30.0].
Each time series receives a different initial condition, which makes them different from each other. We also define a training set of 200 time series for each class, totaling 1,000 time series in the training set. Furthermore, we define an additional 50 time series for each class as our test set, for a total of 250 time series.
First, we can prepare some utilities to help generate our data:
using DynamicalSystemsBase, PredefinedDynamicalSystems
ρ_cls = [26.0, 27.0, 28.0, 29.0, 30.0]
num_samples_to_test = 50
num_samples_to_train = 200
train_timeseries = Vector{StateSpaceSet}(undef, length(ρ_cls) * num_samples_to_train)
test_timeseries = Vector{StateSpaceSet}(undef, length(ρ_cls) * num_samples_to_test)
train_labels = Vector{Float64}(undef, length(ρ_cls) * num_samples_to_train)
test_labels = Vector{Float64}(undef, length(ρ_cls) * num_samples_to_test)
function data_traj(ρ; u0 = rand(3), t = 250.0, Ttr = 1200.0, Δt = 0.2)
lorenz = PredefinedDynamicalSystems.lorenz(u0; ρ)
x, _ = trajectory(lorenz, t; Ttr, Δt)
return x
end
function generate_data!(labels, data, classes, num_elements_per_class)
c = 1
for i ∈ eachindex(labels)
labels[i] = classes[c]
data[i] = data_traj(classes[c])
if (i % num_elements_per_class == 0)
c += 1
end
end
endgenerate_data! (generic function with 1 method)Then, we can generate a training dataset:
generate_data!(train_labels, train_timeseries, ρ_cls, num_samples_to_train)
train_timeseries1000-element Vector{StateSpaceSet}:
3-dimensional StateSpaceSet{Float64} with 1251 points
3-dimensional StateSpaceSet{Float64} with 1251 points
3-dimensional StateSpaceSet{Float64} with 1251 points
3-dimensional StateSpaceSet{Float64} with 1251 points
3-dimensional StateSpaceSet{Float64} with 1251 points
3-dimensional StateSpaceSet{Float64} with 1251 points
3-dimensional StateSpaceSet{Float64} with 1251 points
3-dimensional StateSpaceSet{Float64} with 1251 points
3-dimensional StateSpaceSet{Float64} with 1251 points
3-dimensional StateSpaceSet{Float64} with 1251 points
⋮
3-dimensional StateSpaceSet{Float64} with 1251 points
3-dimensional StateSpaceSet{Float64} with 1251 points
3-dimensional StateSpaceSet{Float64} with 1251 points
3-dimensional StateSpaceSet{Float64} with 1251 points
3-dimensional StateSpaceSet{Float64} with 1251 points
3-dimensional StateSpaceSet{Float64} with 1251 points
3-dimensional StateSpaceSet{Float64} with 1251 points
3-dimensional StateSpaceSet{Float64} with 1251 points
3-dimensional StateSpaceSet{Float64} with 1251 pointsAnd a test dataset:
generate_data!(test_labels, test_timeseries, ρ_cls, num_samples_to_test)
test_timeseries250-element Vector{StateSpaceSet}:
3-dimensional StateSpaceSet{Float64} with 1251 points
3-dimensional StateSpaceSet{Float64} with 1251 points
3-dimensional StateSpaceSet{Float64} with 1251 points
3-dimensional StateSpaceSet{Float64} with 1251 points
3-dimensional StateSpaceSet{Float64} with 1251 points
3-dimensional StateSpaceSet{Float64} with 1251 points
3-dimensional StateSpaceSet{Float64} with 1251 points
3-dimensional StateSpaceSet{Float64} with 1251 points
3-dimensional StateSpaceSet{Float64} with 1251 points
3-dimensional StateSpaceSet{Float64} with 1251 points
⋮
3-dimensional StateSpaceSet{Float64} with 1251 points
3-dimensional StateSpaceSet{Float64} with 1251 points
3-dimensional StateSpaceSet{Float64} with 1251 points
3-dimensional StateSpaceSet{Float64} with 1251 points
3-dimensional StateSpaceSet{Float64} with 1251 points
3-dimensional StateSpaceSet{Float64} with 1251 points
3-dimensional StateSpaceSet{Float64} with 1251 points
3-dimensional StateSpaceSet{Float64} with 1251 points
3-dimensional StateSpaceSet{Float64} with 1251 pointsPreparing the input features
Our feature is the recurrence microstate distribution. Therefore, for each time series we must compute the probabilities and store them as a feature vector:
using RecurrenceMicrostatesAnalysis
N = 3
train_dataset = Matrix{Float64}(undef, 2^(N * N) + 2, length(train_labels))
test_dataset = Matrix{Float64}(undef, 2^(N * N) + 2, length(test_labels))
function get_probs!(dataset, timeseries, N)
for i ∈ eachindex(timeseries)
th, s = optimize(Threshold(), Shannon(), N, timeseries[i])
rmspace = RecurrenceMicrostates(th, N)
dist = probabilities(rmspace, timeseries[i])
dataset[1, i] = th
dataset[2, i] = s
dataset[3:end, i] = dist[1:end]
end
endget_probs! (generic function with 1 method)Note that we are also using the recurrence threshold and the recurrence entropy as features. This follows (Spezzatto et al., 2024), but it is not strictly necessary. If your input has a stable threshold range that maximizes the recurrence entropy, it is possible to use a mean threshold to compute all probability distributions and use only these distributions as features for the machine learning model.
Computing our feature vectors:
- Train
get_probs!(train_dataset, train_timeseries, N)
train_dataset514×1000 Matrix{Float64}:
16.7166 16.2533 16.1379 … 17.6861 17.7293
7.46743 7.42574 7.45215 7.60649 7.60631
0.00148716 0.00624351 0.00538455 0.00508968 0.00114101
0.0187305 0.0222818 0.0196921 0.0168331 0.0158331
0.00217946 0.00512814 0.00393585 0.00519224 0.00334611
0.0124229 0.013628 0.0138332 … 0.0126152 0.0121665
0.00726914 0.0135255 0.0119742 0.00915373 0.00676914
0.00571788 0.00680761 0.00715375 0.00742298 0.00682043
0.00701273 0.00871784 0.00784605 0.00716657 0.00578198
0.00251279 0.00279484 0.0030256 0.00146152 0.000858963
⋮ ⋱
0.00314099 0.00415379 0.00275637 … 0.0023333 0.00225638
0.00910245 0.00858963 0.00733324 0.00949988 0.00788451
0.00958962 0.00801272 0.00815374 0.00505122 0.00512814
8.97424e-5 0.000102563 8.97424e-5 0.000141024 8.97424e-5
0.0064871 0.00656402 0.00610249 0.00389739 0.00417943
0.000512814 0.000499994 0.000320509 … 0.000333329 0.00035897
0.00285894 0.00379482 0.00280766 0.00251279 0.00237176
0.0015769 0.00192305 0.00191023 0.00135896 0.00153844
0.0235382 0.0528968 0.0349611 0.0454481 0.0203972- Test
get_probs!(test_dataset, test_timeseries, N)
test_dataset514×250 Matrix{Float64}:
16.1356 17.1981 16.0439 … 17.8764 18.3128
7.42969 7.44596 7.45518 7.62026 7.64017
0.00925629 0.000217946 0.00639735 0.00293586 0.00133332
0.023051 0.0145511 0.0219484 0.0159998 0.0144229
0.00517942 0.00116665 0.0049743 0.00474353 0.00348713
0.0134229 0.011205 0.0131793 … 0.0111793 0.0104358
0.0137306 0.00608967 0.0132947 0.0079358 0.00629479
0.00507686 0.00603838 0.00611531 0.00767939 0.00687171
0.00843579 0.00512814 0.00776913 0.00693581 0.00624351
0.00344867 0.00191023 0.00299996 0.00143588 0.000961526
⋮ ⋱
0.00364098 0.0029615 0.0028974 … 0.00285894 0.00315381
0.00782041 0.0100255 0.00660248 0.00880758 0.00991013
0.00751272 0.0108588 0.00810246 0.00533326 0.00623069
0.000141024 0.000166665 6.41017e-5 0.000333329 0.000333329
0.00566659 0.00803836 0.00576916 0.0044743 0.00519224
0.000615377 0.000807682 0.000641017 … 0.000461533 0.000551275
0.0037179 0.00312817 0.00323073 0.00292304 0.00260253
0.00224356 0.00187177 0.00253843 0.00189741 0.00191023
0.0483327 0.0273714 0.0376662 0.0462045 0.0480763Defining the neural network model
Here, we use Flux.jl to define our model. As mentioned before, we are using a multi-layer perceptron. However, it is also possible to use other approaches, e.g., random forests (Silveira et al., 2025).
using Flux
model = Chain(
Dense(2^(N * N) + 2 => 512, identity),
Dense(512 => 256, selu),
Dense(256 => 64, selu),
Dense(64 => length(ρ_cls)),
softmax
)
model = f64(model)Chain(
Dense(514 => 512), # 263_680 parameters
Dense(512 => 256, selu), # 131_328 parameters
Dense(256 => 64, selu), # 16_448 parameters
Dense(64 => 5), # 325 parameters
NNlib.softmax,
) # Total: 8 arrays, 411_781 parameters, 3.142 MiB.Training the MLP
Now it is just standard machine learning procedure...
train_labels = Flux.onehotbatch(train_labels, ρ_cls)
test_labels = Flux.onehotbatch(test_labels, ρ_cls)
loader = Flux.DataLoader((train_dataset, train_labels), batchsize = 32, shuffle = true)
opt = Flux.setup(Flux.Adam(0.005), model)
for epc ∈ 1:50
for (x, y) ∈ loader
_, grads = Flux.withgradient(model) do m
y_hat = m(x)
Flux.crossentropy(y_hat, y)
end
Flux.update!(opt, model, grads[1])
end
endModel evaluation
Finally, let's check our accuracy 🙂
using LinearAlgebra
function get_quantifiers(predict, trusty, classes)
conf = zeros(Int, length(classes), length(classes))
sz = size(predict, 2)
for i in 1:sz
mx_prd = findmax(predict[:, i])
mx_trt = findmax(trusty[:, i])
conf[mx_prd[2], mx_trt[2]] += 1
end
return tr(conf) / (sum(conf) + eps())
end
accuracy = get_quantifiers(model(test_dataset), test_labels, ρ_cls)
accuracy * 100100.0