%% Train random forest for distinguishing rhythms
%
% Joel Karel
% Maastricht University
% 16-03-2017
%% load labels training set

reffile = ['~/Documents/Research/CinCChallenge2017/CinC2017/training2017'  filesep 'REFERENCE.csv'];
% reffile = ['~/Documents/Research/CinCChallenge2017/CinC2017/training2017' filesep 'REFERENCE.csv'];

fid = fopen(reffile, 'r');
if(fid ~= -1)
    Ref = textscan(fid,'%s %s','Delimiter',',');
else
    error(['Could not open ' reffile ' for scoring. Exiting...'])
end
fclose(fid);

RECORDS = Ref{1};
target  = Ref{2};
N       = length(RECORDS);


% find indices of different groups in REFERENCE file
idx1 = ismember(Ref{2},'N');
idxN=find(idx1);
idx2 = ismember(Ref{2},'A');
idxA=find(idx2);
idx3 = ismember(Ref{2},'O');
idxO=find(idx3);
idx4 = ismember(Ref{2},'~');
idxnoise=find(idx4);

%% load parameters

%load('~/Documents/Research/CinCChallenge2017/workinprogress/newresults.mat') 
load newparameters
x=newparameters(:,1:46);
%remove saen1 and saen2
x(:,18:19)=[];

Parameters = {'P_LF'; 'f_LF'; 'P_HF'; 'f_HF'; 'meanRR'; 'SDRR'; 'rMSSDRR'; 'pNN10'; 'pNN20'; ...
    'pNN30'; 'pNN40'; 'pNN50'; 'pNN60'; 'pNN70'; 'pNN80';'pNN90'; 'pNN100'; ...
    'V0'; 'V1'; 'LV2'; 'UV2'; 'numbercomponents95percQRST'; 'SimQRST1'; 'SimQRST2'; ...
    'SimQRST3'; 'SimQRST4'; 'SimQRST5'; 'numbercomponents95percP'; 'SimP1'; 'SimP2'; ...
    'SimP3'; 'SimP4'; 'SimP5'; 'SCN'; 'DF'; 'fWP'; 'fWPMAW'; 'rhe';  'SampEn'; 'SampEnMAW'; ... 
    'OI'; 'SE'; 'FWA'; 'differenceMediaRRmeanRR'};


%% convert labels to numerical values
Labels = char(Ref{2});
y = double(Labels);
numel(find(y == 78));
% 78 = N
numel(find(y == 79));
% 79 = O
numel(find(y == 65));
% 65 = A
numel(find(y == 126));
% 126 = ~

y(y==78) = 1;
y(y==79) = 3;
y(y==65) = 2;
y(y==126) = 4;

%% Exclude noisy signals since we will catch those earlier
x(y==4,:)=[];
Labels(y==4)=[];
y(y==4)=[];

%% Split in training and validation set - 75% vs 25%
indx=rand(size(x,1),1)<0.75;
train.x=x(indx,:);
train.y=y(indx);
train.Labels=Labels(indx);
val.x=x(~indx,:);
val.y=y(~indx);
val.Labels=Labels(~indx);

%% The actual model building
rng(1); % For reproducibility

%c.ClassNames = ['A','N','O'];
%c.ClassificationCosts = [0 5 3; 3 0 1; 3 1 0];
C=[0 5 3; 3 0 1; 3 1 0];

Mdl = TreeBagger(50,train.x,train.Labels,'OOBPrediction','On',...
    'Method','classification','Cost',C,'ClassNames',{'A','N','O'})

figure;
oobErrorBaggedEnsemble = oobError(Mdl);
plot(oobErrorBaggedEnsemble)
xlabel 'Number of grown trees';
ylabel 'Out-of-bag classification error';

% https://nl.mathworks.com/help/stats/ensemble-methods.html


%% Learn ensembles with RUSBoost to handle class imbalance
% Note that the results vary from training to training

% see how much each class is represented
tabulate(train.Labels)
tabulate(val.Labels)
tabulate(Labels)

% fit the ensemble with deep trees
t = templateTree('MinLeafSize',5);
tic
rusTree = fitensemble(train.x,train.Labels,'RUSBoost',1000,t,...
    'LearnRate',0.1,'nprint',100);
toc

% Plot the classification error against the number of members in the ensemble
figure;
tic
plot(loss(rusTree,val.x,val.Labels,'mode','cumulative'));
toc
grid on;
xlabel('Number of trees');
ylabel('Test classification error');

% Examine the confusion matrix for each class as a percentage of the true class.
tic
Yfit = predict(rusTree,val.x);
toc
tab = tabulate(val.y);
bsxfun(@rdivide,confusionmat(int8(val.Labels),int8(Yfit)),tab(:,2))*100

% The ensemble is large. Remove the data using the compact method.
cmpctRus = compact(rusTree);

sz(1) = whos('rusTree');
sz(2) = whos('cmpctRus');
[sz(1).bytes sz(2).bytes]

% Remove half the trees from cmpctRus. This action is likely to have minimal effect on the predictive performance, based on the observation that 400 out of 1000 trees give nearly optimal accuracy.
cmpctRus = removeLearners(cmpctRus,[500:1000]);
sz(3) = whos('cmpctRus');
sz(3).bytes

% The reduced compact ensemble takes about a quarter the memory of the full ensemble. Its overall loss rate is under 24%:
tic
Yfit = predict(cmpctRus,val.x);
toc
tab = tabulate(val.y);
confusionmat(int8(val.Labels),int8(Yfit))


%save model cmpctRus