function [network, Xtraining, Ytraining, Xtesting, Ytesting]= Bidirectional_LSTM_network(input_variables, opts)
nt= length(input_variables);

% divide data in Test and Train parts.
trained_timepoints=opts.trained_timepoints;

train_size= trained_timepoints/nt;                 
delays= opts.delays;
% % % % % % % % % % % % Normalize data
mu= mean(input_variables,'omitnan'); sigma= std(input_variables,'omitnan');
data= (input_variables'-mu)./sigma;

% % % Divide input data into training and testing dataset % % % % % % % % % 
numTrainPoints = floor(train_size*numel(data));
Xtraining   = data(1:numTrainPoints);
Xtesting   = data(numTrainPoints+1:end);

% % % Delays on the time series % % % % % % % % % % % % % % % % % % % % % % 
ntrain = length(Xtraining);
max_delay = max(delays);
range = max_delay+1:ntrain;
X= []; %tX=[];
for id = delays
    X=[X; Xtraining(:,range-id)];
%     tX=[tX; chain(:,range-id)];
end
Y = Xtraining(:,range);
Xtraining = X;
Ytraining = Y;
% 
% figure, plot(Xtraining,'r'); hold on; plot(Ytraining,'b');
% return;
ntest = length(Xtesting);
range = max_delay+1:ntest;
X= []; %tX=[];
for id = delays; 
    X=[X; Xtesting(:,range-id)];
%     tX=[tX; chain(:,range-id)];
end
Y = Xtesting(:,range);
Xtesting = X;
Ytesting = Y;
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %     

% Bi-LSTM parameters   
layers = [ ...
    sequenceInputLayer(size(delays,2))
    bilstmLayer(opts.h1,'OutputMode','sequence')
    bilstmLayer(opts.h2,'OutputMode','last')  
    fullyConnectedLayer(1)
    regressionLayer];

opts = trainingOptions('adam', ...
    'MaxEpochs',opts.maxepochs, ...
    'MiniBatchSize',opts.minibatchsize, ...
    'GradientThreshold',Inf, ...
    'L2Regularization', opts.l2reg, ...
    'InitialLearnRate',0.001, ...
    'LearnRateSchedule','piecewise', ...
    'LearnRateDropPeriod',125, ...
    'LearnRateDropFactor',opts.learndropout, ...
    'Shuffle','never', ...     
    'Verbose',0, ...
    'ExecutionEnvironment','cpu', ...
    'Plots','none');

%     'SequenceLength', 'longest', ...
%     'SequencePaddingValue', 20, ...
% Train BiLSTM network 
% for i=1:length(Xtraining)
%     Xtrain{i,1}= Xtraining(:,i); 
%     Ytrain(i,1) = Ytraining(:,i);
% end
% whos Xtrain
% whos Ytrain
for i=1:size(Xtraining,2)
    Xtrain{i,1} = Xtraining(:,i);  
    Ytrain(i,1) = Ytraining(:,i);  
end
% for i=1:size(Xtesting,2)
%     Xtest{i,1} = Xtesting(:,i);  
%     Ytest(i,1) = Ytesting(:,i); 
% end
network = trainNetwork(Xtrain,Ytrain,layers,opts);

return;
