%% Train random forest for detecting noise
%
% Joel Karel
% Maastricht University
% 16-03-2017
data_dir = ['~/Documents/Research/CinCChallenge2017/CinC2017/training2017' filesep];
%% load labels training set

reffile = [data_dir '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);

%% generate parameters
x=[];
for cnt=1:size(RECORDS,1)
    [tm,ecg,fs,siginfo]=rdmat([data_dir RECORDS{cnt}]);
    x =[x;getNoiseFeatures(ecg)];
end
    

%% 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) = 1;
y(y==65) = 1;
y(y==126) = 0;

%% Exclude noisy signals since we will catch those earlier
Labels(y==1)='O';

%% 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);


%% 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
noiseTree = 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(noiseTree,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(noiseTree,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.
cmpctNoise = compact(noiseTree);

sz(1) = whos('noiseTree');
sz(2) = whos('cmpctNoise');
[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(cmpctNoise,[500:1000]);
sz(3) = whos('cmpctNoise');
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(cmpctNoise,val.x);
toc
tab = tabulate(val.y);
confusionmat(int8(val.Labels),int8(Yfit))


save noisemodel cmpctNoise
load model
save overallmodel cmpctNoise cmpctRus
