% random_forest_classifier.m

% ******************************************************************

%  Script to classify the Greenland valley dataset using random forest

% *******************************************************************

clear
close all

%% Assemble training dataset

trainingdata = readtable('TrainingDataValleys.csv','VariableNamingRule','preserve');

% Extract morphometric variables and labels
X = table2array(trainingdata(:,4:7));
label = table2array(trainingdata(:,8));

Nx = size(X,1);

% 1 = glacial; 2 = fluvial
Y = zeros(Nx,1);
glac_idx = strcmp(label,'Glacial');
fluv_idx = strcmp(label,'Fluvial');
Y(glac_idx,1) = 1;
Y(fluv_idx,1) = 2;

% Shuffle the order of the dataset
a = randperm(Nx)';
for i=1:numel(X(1,:))
    X(:,i) = X(a,i);
end
Y = Y(a);

%% Extract full unlabelled valley dataset from Greenland

fulldata = readtable('GreenlandValleys.csv','VariableNamingRule','preserve');

% Extract morphometrics and examine cross-correlation
A = table2array(fulldata(:,10:13));
CC = corrcoef(A);

%% Train model and estimate misclassification error using cross validation

% Set number of trees in forest and number of folds for cross validation
NumTrees = 1000;
folds = 10;

% Train random forest model
Mdl = TreeBagger(NumTrees,X,Y,'SampleWithReplacement','on',...
    'OOBPrediction','on','OOBPredictorImportance','on','Method','classification',...
    'PredictorNames',{'D','W','V','C'});

% K-fold validation
cvp = cvpartition(Y,'KFold',folds);
cvError = crossval('mcr',X,Y,'Predfun',@classtb,'Partition',cvp);

% Out-of-bag classification error
oobce = oobError(Mdl);
oobce_final = oobce(end);

%% Apply trained model to full dataset

% Classification of full data
[Yfit, scores] = predict(Mdl,A);

Yfit = str2double(Yfit);    % 1 = glacial, 2 = fluvial
glac_score = scores(:,1);   % 1 = glacial in every tree, 0 = fluvial in every tree
fluv_score = scores(:,2);   % 0 = glacial in every tree, 1 = fluvial in every tree

%% Number of valleys of different classes

N = numel(Yfit);
Ng = sum(Yfit(:) == 1);
Nf = sum(Yfit(:) == 2);

%% Crossval function

function Yfit = classtb(XTrain,YTrain,XTest)

NumTrees = 100;  % Note this value can be different to the main code to speed up results
Mdl = TreeBagger(NumTrees,XTrain,YTrain,'SampleWithReplacement','on',...
    'OOBPrediction','on','OOBPredictorImportance','on','Method','classification',...
    'PredictorNames',{'D','W','V','C'});

[Yfit, ~] = predict(Mdl,XTest);
Yfit = str2double(Yfit);

end