A simple and flexible code for Reservoir Computing architectures like Echo State Networks.
This toolbox works for Python 3 (and should be compatible for Python2). We just updated it from Python 2 to 3, so tell us if you have any issue with it.
Run and analyse these two files to see how to make timeseries prediction with Echo State Networks:
-
simple_example_MackeyGlass.py (using the ESN class)
python simple_example_MackeyGlass.py
-
minimalESN_MackeyGlass.py (without the ESN class)
python minimalESN_MackeyGlass.py
You can generate and train a reservoir to predict the MackeyGlass timeseries in a few steps:
-
Define the number of dimension of input, recurrent and outputs layers:
n_inputs = 1 input_bias = True # add a constant input to 1 n_outputs = 1 N = 300 # number of recurrent units
-
Define the random input Win and recurrent W matrices (quick method with generation tools):
import mat_gen W = mat_gen.generate_internal_weights(N=N, spectral_radius=1.0, proba=1.0, Wstd=1.0) # Normal distribution with mean 0 and standard deviation 0 Win = mat_gen.generate_input_weights(nbr_neuron=N, dim_input=n_inputs, input_scaling=1.0, proba=1.0, input_bias=input_bias)
-
(instead of previous step) Define yourself the random input Win and recurrent W matrices (customize method):
import numpy as np # Generating matrices Win and W W = np.random.rand(N,N) - 0.5 if input_bias: Win = np.random.rand(N,dim_inp+1) - 0.5 else: Win = np.random.rand(N,dim_inp) - 0.5 # Apply mask to make matrices sparse proba_non_zero_connec_W = 0.2 # set the probability of non-zero connections mask = np.random.rand(N,N) # create a mask with Uniform[0;1] distribution W[mask > proba_non_zero_connec_W] = 0 # apply mask on W: set to zero some connections given by the mask mask = np.random.rand(N,Win.shape[1]) # Do the same for input matrix Win[mask > proba_non_zero_connec_Win] = 0 # Scaling matrices Win and W input_scaling = 1.0 # Define the scaling of the input matrix Win Win = Win * input_scaling # Apply scaling spectral_radius = 1.0 # Define the scaling of the recurrent matrix W print 'Computing spectral radius ...', original_spectral_radius = np.max(np.abs(np.linalg.eigvals(W))) W = W * (spectral_radius / original_spectral_radius) # Rescale W to reach the requested spectral radius
-
Define the Echo State Network (ESN):
import ESN reservoir = ESN.ESN(lr=leak_rate, W=W, Win=Win, input_bias=input_bias, ridge=regularization_coef, Wfb=None, fbfunc=None)
-
Define your input/output training and testing data:
In this step, we load the dataset to perform the prediction of the chaotic MackeyGlass timeseries, and we split the data into the different subsets.
data = np.loadtxt('MackeyGlass_t17.txt') train_in = data[None,0:trainLen].T # input data (TRAINING PHASE) train_out = data[None,0+1:trainLen+1].T # output to be predicted (TRAINING PHASE) test_in = data[None,trainLen:trainLen+testLen].T # input data (TESTING PHASE) test_out = data[None,trainLen+1:trainLen+testLen+1].T # output to be predicted (TESTING PHASE)
-
Train the ESN:
Be careful to give lists for the input and output (i.e. teachers) training data. Here we are training with only one timeseries, but you actually can provide a list of timeseries segments to train from.
wash_nr_time_step defines the initial warming-up period: corresponding reservoir states are discarded for training.
internal_trained = reservoir.train(inputs=[train_in,], teachers=[train_out,], wash_nr_time_step=100)
-
Test the ESN (i.e. predict the next value in the timeseries):
output_pred, internal_pred = reservoir.run(inputs=[test_in,], reset_state=False)
-
Compute the error made on test data:
print("\nRoot Mean Squared error:") print(np.sqrt(np.mean((output_pred[0] - test_out)**2))/testLen)
-
Plot the internal states of the ESN and the outputs for test data.
If the training was sucessful, predicted output and real curves should overlap:
import matplotlib.pyplot as plt ## Plot internal states of the ESN plt.figure() plt.plot( internal_trained[0][:200,:12]) plt.ylim([-1.1,1.1]) plt.title('Activations $\mathbf{x}(n)$ from Reservoir Neurons ID 0 to 11 for 200 time steps') ## Plot the predicted output with the real one. plt.figure(figsize=(12,4)) plt.plot(test_out, color='0.75', lw=1.0) plt.plot(output_pred[0], color='0.00', lw=1.5) plt.title("Ouput predictions against real timeseries") plt.legend(["real timeseries", "output predictions"]) # plt.ylim(-1.1,1.1) plt.show()
If you want to have more information on all the steps and more option (for example, have a reservoir with output feedback), please have a look at simple_example_MackeyGlass.py.