function TrainClassifier(feature_file)
% This function extracts features for each record present  in a folder
%
%  Input:
%       - feature_file:         file containing table with extracted
%                               features in different records
%       
% --
% ECG classification from single-lead segments using Deep Convolutional Neural 
% Networks and Feature-Based Approaches - December 2017
% 
% Released under the GNU General Public License
%
% Copyright (C) 2017  Fernando Andreotti, Oliver Carr
% University of Oxford, Insitute of Biomedical Engineering, CIBIM Lab - Oxford 2017
% fernando.andreotti@eng.ox.ac.uk
%
% 
% For more information visit: https://github.com/fernandoandreotti/cinc-challenge2017
% 
% Referencing this work
%
% Andreotti, F., Carr, O., Pimentel, M.A.F., Mahdi, A., & De Vos, M. (2017). 
% Comparing Feature Based Classifiers and Convolutional Neural Networks to Detect 
% Arrhythmia from Short Segments of ECG. In Computing in Cardiology. Rennes (France).
%
% Last updated : December 2017
% 
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
% 
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
% 
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

load(feature_file)
NFEAT=size(allfeats,2);
NFEAT=NFEAT-2;

% Get summary statistics on the distribution of the features in each
% signal, using following:
% - median
% - inter quartile range
% - range
% - min value
% - max value
% - 25% perctile
% - 50% perctile
% - 75% percentile
% - Real coefficients of Hilbert transform 
% - Absolute values of Hilbert transform
% - Skewness
% - Kurtosis
% 
feat = zeros(max(allfeats.rec_number),16*NFEAT);
for i=1:max(allfeats.rec_number)
    fprintf('Processing record %d .. \n',i)
    ind=find(table2array(allfeats(:,1))==i);
    feat(i,1:NFEAT)=nanmean(table2array(allfeats(ind,3:end)));
    feat(i,1*NFEAT+1:2*NFEAT)=nanstd(table2array(allfeats(ind,3:end)));
    if length(ind)>2
        PCAn=pca(table2array(allfeats(ind,3:end)));
        feat(i,2*NFEAT+1:3*NFEAT)=PCAn(:,1);
        feat(i,3*NFEAT+1:4*NFEAT)=PCAn(:,2);
    else
        feat(i,2*NFEAT+1:3*NFEAT)=NaN;
        feat(i,3*NFEAT+1:4*NFEAT)=NaN;
    end
    feat(i,4*NFEAT+1:5*NFEAT)=nanmedian(table2array(allfeats(ind,3:end)));
    feat(i,5*NFEAT+1:6*NFEAT)=iqr(table2array(allfeats(ind,3:end)));
    feat(i,6*NFEAT+1:7*NFEAT)=range(table2array(allfeats(ind,3:end)));
    feat(i,7*NFEAT+1:8*NFEAT)=min(table2array(allfeats(ind,3:end)));
    feat(i,8*NFEAT+1:9*NFEAT)=max(table2array(allfeats(ind,3:end)));
    feat(i,9*NFEAT+1:10*NFEAT)=prctile(table2array(allfeats(ind,3:end)),25);
    feat(i,10*NFEAT+1:11*NFEAT)=prctile(table2array(allfeats(ind,3:end)),50);
    feat(i,11*NFEAT+1:12*NFEAT)=prctile(table2array(allfeats(ind,3:end)),75);
    HIL=hilbert(table2array(allfeats(ind,3:end)));
    feat(i,12*NFEAT+1:13*NFEAT)=real(HIL(1,:));
    feat(i,13*NFEAT+1:14*NFEAT)=abs(HIL(1,:));
    feat(i,14*NFEAT+1:15*NFEAT)=skewness(table2array(allfeats(ind,3:end)));
    feat(i,15*NFEAT+1:16*NFEAT)=kurtosis(table2array(allfeats(ind,3:end))); 
end

In = feat;
Ntrain = size(In,1);
In(isnan(In)) = 0;
% Standardizing input
In = In - mean(In);
In = In./std(In);

labels = {'A' 'N' 'O' '~'};
Out = reference_tab{:,2};
Outbi = cell2mat(cellfun(@(x) strcmp(x,labels),Out,'UniformOutput',0));
Outde = bi2de(Outbi);
Outde(Outde == 4) = 3;
Outde(Outde == 8) = 4;
clear Out
rng(1); % For reproducibility
%% Perform cross-validation
%== Subset sampling
k = 5;
cv = cvpartition(Outde,'kfold',k);
confusion = zeros(4,4,k);
F1save = zeros(k,4);
F1_best = 0;
for i=1:k
    fprintf('Cross-validation loop %d \n',i)
    trainidx = find(training(cv,i));
    trainidx = trainidx(randperm(length(trainidx)));
    testidx  = find(test(cv,i));
    %% Bagged trees (oversampled)
    ens = fitensemble(In(trainidx,:),Outde(trainidx),'Bag',50,'Tree','type','classification');
    [~,probTree] = predict(ens,In(testidx,:));
    
    %% Neural networks
    net = patternnet(10);
    net = train(net,In(trainidx,:)',Outbi(trainidx,:)');            
    probNN = net(In(testidx,:)')';    
    
    %% Combining methods
    C = cat(3,probTree,probNN);
    C = mean(C,3);
    estimate = zeros(size(C,1),1);
    for r = 1:size(C,1)
        [~,estimate(r)] = max(C(r,:));
    end
    confmat = confusionmat(Outde(testidx),estimate);
    confusion(:,:,i) = confmat;    
    F1 = zeros(1,4);
    for j = 1:4
        F1(j)=2*confmat(j,j)/(sum(confmat(j,:))+sum(confmat(:,j)));
        fprintf('F1 measure for %s rhythm: %1.4f \n',labels{j},F1(j))
    end
    F1save(i,:) = F1;
    
    if F1 > F1_best
        F1_best = F1;
        ensTree_best = compact(ens);
        nnet_best = net;
    end
end
%% Producing statistics
confusion = sum(confusion,3);
F1 = zeros(1,4);
for i = 1:4
    F1(i)=2*confusion(i,i)/(sum(confusion(i,:))+sum(confusion(:,i)));
    fprintf('F1 measure for %s rhythm: %1.4f \n',labels{i},F1(i))
end
fprintf('Final F1 measure:  %1.4f\n',mean(F1))

%% Save output
save('results_allfeat.mat','F1save','F1_best')
save('ensTree.mat','ensTree_best')
save('nNets.mat','nnet_best')