% load data
data = load('MackeyGlass_t17.txt');
initLen = 100;%前100项用来初始化储备池
trainLen = 2000;
testLen = 2000;
% plot MackeyGlass: 1:1000
figure(10);
plot(data(1:1000));
title('A sample of data');
% generate ESN
inSize = 1;
outSize = 1;
resSize = 1000;
a = 0.93;
rand( 'seed' ,40);%把随机数保存了一样,下次用的时候再拿出来。
Win = (rand(resSize,inSize+1)-0.5) .* 1;
W = rand(resSize,resSize)-0.5;
%normalizing and setting spectral radius
disp 'Computing spectral radius'
opt.disp = 0;
rhoW = abs(eigs(W,1,'lm',opt));
disp 'done'
W = W .* (1.25 / rhoW);
% collect state matrix
X = zeros(1+inSize+resSize,trainLen-initLen);
% target matrix
Yt = data(initLen+2:trainLen+1)';
% collect X by run the reservoir
x = zeros(resSize,1);
for t = 1:trainLen
u = data(t);
x = (1-a)*x + a*(eye(1000) ./ (1 + exp(-Win*[1;u] - W*x )));
%x = (1-a)*x + a*tanh( Win*[1;u] + W*x);
if t > initLen
X(:,t-initLen) = [1;u;x];
end
end
% train output
reg = 1e-8;
Wout = ((X*X' + reg*eye(1+resSize+inSize)) \ (X*Yt'))';
% by test data collect Y
Y = zeros(outSize,testLen);
u = data(trainLen+1);
for t = 1:testLen
x = (1-a)*x + a*(eye(1000) / (1 + exp(-Win*[1;u] - W*x )));
%x = (1-a)*x + a*tanh( Win*[1;u] + W*x );%1./(1 + exp())
y = Wout*[1;u;x];
Y(:,t) = y;
u = y;
end
% compute error
errorLen = 500;
mse = sum((data(trainLen+2:trainLen+errorLen+1)'-Y(1,1:errorLen)).^2)./errorLen;
disp( ['MSE = ' ,num2str( mse )]);
% plot signals
figure(1);
plot( data(trainLen+2:trainLen+testLen+1),'color', [0,0.75,0] );
hold on;
plot(Y','b');
hold off;
axis tight;
title('From n=0 target and generate signals');
legend('target signals','free-running predicted signals');
figure(2);
plot(X(1:20,1:200)');
title('reservoir x(n)');
figure(3);
bar(Wout');
title('output weight Wout');