% This program will train LSTM with x y z axis raw data from 
% vehicle starting faults
% 8/9/2024
% Coded in MATLAB R2024b
% Place file under MATLAB folder NIS/Audiofaults

%Remove previous setting
close all 
clear 

%Location for input fault vibration files
dataFolder = "./NIS/Axisfaults";

[files,flabels] = getmyfileslog(dataFolder);
[tsize,z] = size(files);

%Get data from logs
for idx=1:1:tsize
    tdata = readtable(files(idx),'FileType','text');
    data(idx,1) = {table2array(tdata)};
    labels(idx,1) = flabels(idx);
end

[asize,z]=size(data);
numChannels = size(data{1},2);

%Show vibration data in tiled time series graph
figure
tiledlayout(4,4)
for i = 1:1:asize
    nexttile
    stackedplot(data{i},DisplayLabels="Channel "+string(1:numChannels))
    xlabel("Time Step")
    title("Class: " + string(labels(i)))
end

classNames = categories(labels)

%Set training and validation partition as 90 10
numObservations = numel(data);
[idxTrain,idxTest] = trainingPartitions(numObservations,[0.9 0.1]);
XTrain = data(idxTrain);
TTrain = labels(idxTrain);

XTest = data(idxTest);
TTest = labels(idxTest);

numObservationsXTrain = numel(XTrain);
for i=1:numObservationsXTrain
    sequence = XTrain{i};
    sequenceLengths(i) = size(sequence,1);
end

[sequenceLengths,idx] = sort(sequenceLengths);
dataTrain = XTrain(idx);
textTrain = TTrain(idx);

%Show vibration fault data sequence length
figure
bar(sequenceLengths)
xlabel("Sequence")
ylabel("Length")
title("Sorted Data")

%Set training environment values
numHiddenUnits = 400;
numClasses = 5;

layers = [
    sequenceInputLayer(numChannels)
    bilstmLayer(numHiddenUnits,OutputMode="last")
    fullyConnectedLayer(numClasses)
    softmaxLayer]

options = trainingOptions("adam", ...
    MaxEpochs=100, ...
    InitialLearnRate=0.0002,...
    GradientThreshold=1, ...
    Shuffle="never", ...
    Plots="training-progress", ...
    Metrics="accuracy", ...
    Verbose=false);

%Train vibration tool
VibFN = trainnet(dataTrain,textTrain,layers,"crossentropy",options);

%Save trained vibration tool
save VibFN.mat

%Test classification accuracy with singular files
%Edit and adjust as needed
tdata = readtable('axisCoilFault1.log','FileType','text');
data1 = table2array(tdata); 
tdata = readtable('axisCoilFault2.log','FileType','text');
data2 = table2array(tdata); 
tdata = readtable('axisCoilFault3.log','FileType','text');
data3 = table2array(tdata); 
XTest={data1,data2,data3};

TTest = categorical(["BadCoil","BadCoil","BadCoil"]);
numObservationsTest = numel(XTest);
for i=1:numObservationsTest
    sequence = XTest{i};
    sequenceLengthsTest(i) = size(sequence,1);
end

[sequenceLengthsTest,idx] = sort(sequenceLengthsTest);
XTest = XTest(idx);

%Show data sequence length
figure
bar(sequenceLengthsTest)
xlabel("Sequence")
ylabel("Length")
title("Sorted Data")

acc = testnet(VibFN,XTest,TTest,"accuracy")

scores = minibatchpredict(VibFN,XTest);
YTest = scores2label(scores,classNames);

%Show confusion chart matrix with prediction accuracy
figure
confusionchart(TTest,YTest)