%% Sequence-to-Sequence LSTM Sequence Length Diagrams
%% Download Data
% Download and unzip the Turbofan Engine Degradation Simulation Data Set from https://ti.arc.nasa.gov/tech/dash/groups/pcoe/prognostic-data-repository/.

filename = "CMAPSSData.zip";
dataFolder = "data";
unzip(filename,dataFolder)

%
%% Prepare Training Data
% Load the data using the function preprocessDataTrain listed at the end of this example. The function prepareDataTrain extracts the data from filenamePredictors and returns the cell arrays XTrain and YTrain which contain the training predictor and response sequences.

filenamePredictors = fullfile(dataFolder,"train_FD001.txt");
[XTrain,YTrain] = prepareDataTrain(filenamePredictors);

% Remove Features with Constant Values. Find the rows of the data which have the same minimum and maximum values and remove them.

m = min([XTrain{:}],[],2);
M = max([XTrain{:}],[],2);
idxConstant = M == m;

for i = 1:numel(XTrain)
    XTrain{i}(idxConstant,:) = [];
end

% Normalize Training Predictors. Normalize the training predictors to have zero mean and unit variance. To calculate the mean and standard deviation over all observations, concatenate the sequence data horizontally.

mu = mean([XTrain{:}],2);
sig = std([XTrain{:}],0,2);

for i = 1:numel(XTrain)
    XTrain{i} = (XTrain{i} - mu) ./ sig;
end

% Clip Responses. To learn more from the sequence data when the engines are close to failing. Clip the responses at the threshold 150. This makes the network treat instances with higher remaining useful life values as equal.

thr = 150;
for i = 1:numel(YTrain)
    YTrain{i}(YTrain{i} > thr) = thr;
end

% Prepare Data for Padding. To minimize the amount of padding added to the mini-batches, sort the training data by sequence length.

sequenceLengthsUnsorted = cellfun(@(x) size(x,2),XTrain);

[sequenceLengths,idx] = sort(sequenceLengthsUnsorted,'descend');
miniBatchSize = 20;
dataColor = [0 0.4470 0.7410];
paddingColor = [0.9290 0.6940 0.1250];
numSequences = size(XTrain,1);

%

figure
plotRUL(XTrain{1},YTrain{1},thr)

%


figure
i=0;
while (i+1)*miniBatchSize <= numel(sequenceLengths)
    
    idx = i*miniBatchSize+1:(i+1)*miniBatchSize;
    miniBatchSequenceLengths = sequenceLengths(idx)';
    maxLength = max(miniBatchSequenceLengths);
    
    data = [ ...
        miniBatchSequenceLengths;
        maxLength-miniBatchSequenceLengths]';
    
    h = bar(idx,data,1,'stacked','LineStyle','none');
    hold on
    
    h(1).FaceColor = dataColor;
    h(2).FaceColor = paddingColor;
    i=i+1;
end

for i = miniBatchSize:miniBatchSize:numel(sequenceLengths)-1
    h(3) = plot([i i]+.5,ylim,'r--');
end

title("Sorted Data")
xlabel("Sequence")
ylabel("Length")
legend(h,["Sequence Data","Padding","Mini-Batch"+newline+"Boundaries"], ...
    'Location','northeast')

%

h = figure;
h.Position(3) = 2*h.Position(3);
subplot(1,2,2)
i=0;
while (i+1)*miniBatchSize <= numel(sequenceLengths)
    
    idx = i*miniBatchSize+1:(i+1)*miniBatchSize;
    miniBatchSequenceLengths = sequenceLengths(idx)';
    maxLength = max(miniBatchSequenceLengths);
    
    data = [ ...
        miniBatchSequenceLengths;
        maxLength-miniBatchSequenceLengths]';
    
    h = bar(idx,data,1,'stacked','LineStyle','none');
    hold on
    
    h(1).FaceColor = dataColor;
    h(2).FaceColor = paddingColor;
    i=i+1;
end

for i = miniBatchSize:miniBatchSize:numel(sequenceLengths)-1
    h(3) = plot([i i]+.5,ylim,'r--');
end

title("Sorted Data")
xlabel("Sequence")
ylabel("Length")

legend(h,["Sequence Data","Padding","Mini-Batch"+newline+"Boundaries"], ...
    'Location','southeast')

subplot(1,2,1)
i=0;
while (i+1)*miniBatchSize <= numel(sequenceLengthsUnsorted)
    
    idx = i*miniBatchSize+1:(i+1)*miniBatchSize;
    miniBatchSequenceLengths = sequenceLengthsUnsorted(idx)';
    maxLength = max(miniBatchSequenceLengths);
    
    data = [ ...
        miniBatchSequenceLengths;
        maxLength-miniBatchSequenceLengths]';
    
    h = bar(idx,data,1,'stacked','LineStyle','none');
    hold on
    
    h(1).FaceColor = dataColor;
    h(2).FaceColor = paddingColor;
    i=i+1;
end

for i = miniBatchSize:miniBatchSize:numel(sequenceLengthsUnsorted)-1
    h(3) = plot([i i]+.5,ylim,'r--');
end

title("Unsorted Data")
xlabel("Sequence")
ylabel("Length")

%
%% Example Functions
% The function prepareDataTrain extracts the data from filenamePredictors and returns the cell arrays XTrain and YTrain which contain the training predictor and response sequences.

function [XTrain,YTrain] = prepareDataTrain(filenamePredictors)

dataTrain = dlmread(filenamePredictors);

numObservations = max(dataTrain(:,1));

XTrain = cell(numObservations,1);
YTrain = cell(numObservations,1);
for i = 1:numObservations
    idx = dataTrain(:,1) == i;
    
    X = dataTrain(idx,3:end)';
    XTrain{i} = X;
    
    timeSteps = dataTrain(idx,2)';
    Y = fliplr(timeSteps);
    YTrain{i} = Y;
end

end

% The function prepareDataTest extracts the data from filenamePredictors and filenameResponses and returns the cell arrays XTest and YTest which contain the test predictor and response sequences.

function [XTest,YTest] = prepareDataTest(filenamePredictors, filenameResponses)

XTest = prepareDataTrain(filenamePredictors);

RULTest = dlmread(filenameResponses);

numObservations = numel(RULTest);

YTest = cell(numObservations,1);
for i = 1:numObservations
    X = XTest{i};
    sequenceLength = size(X,2);
    
    rul = RULTest(i);
    YTest{i} = rul+sequenceLength-1:-1:rul;
end

end

% The function plotRUL plots the sensor data X and the remaining useful life Y, and overlays the healthy region given by the threshold thr.

function plotRUL(X,Y,thr)

figure
subplot(2,1,1)

% Plot sensor data.
h1 = plot(X','Color',[0.65 0.65 0.65]);

% Plot failure.
hold on
sequenceLength = numel(Y);
h2 = plot([sequenceLength sequenceLength],ylim,'r','LineWidth',2);

% Plot healthy patch.
yy = ylim;
v = sequenceLength - thr + 1;
x = [0 v v 0];
y = [yy(1) yy(1) yy(2) yy(2)];
h3 = patch(x,y,'green', ...
    'LineStyle','none', ...
    'FaceAlpha',0.25);
hold off

title("Sensor Data")
legend([h1(1) h2 h3],["Sensor Data" "Failure" "Heathy"],'Location','southeast')

% Plot RUL.
subplot(2,1,2)
plot(Y)

% Plot failure.
yMax = thr + 25;
hold on
plot([sequenceLength sequenceLength],[0 yMax],'r','LineWidth',2);

% Plot healthy patch.
x = [0 v v 0];
y = [0 0 yMax yMax];
patch(x,y,'green', ...
    'LineStyle','none', ...
    'FaceAlpha',0.25);

hold off

ylim([0 yMax])
xlabel("Time Step")
ylabel("RUL")
title("Remaining Useful Life")

end

% Copyright 2012 The MathWorks, Inc.