allen-classification/g_binary_task.m at main · nadjastojanovic/allen-classification · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
close all;
clear all;

% writing program output to file
log_filename = 'logs/g_F3_595263152_GCamp6f.txt';
diary(log_filename);
diary on;

filename = "data/F3_595263152_GCamp6f.mat";

fprintf('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n');
fprintf('Processing file: %s\n', filename);
load(filename);

num_neurons = size(data, 1) - 3;

% filter to keep only data with temporal frequency of 2 or 8 Hz
data = data(:, data(end-2, :) == 2 | data(end-2, :) == 8);

intervals = unique(data(end-1, :));
features = zeros(num_neurons, length(intervals));

% create the features table
col = 1;
for i = intervals
    for n = 1:num_neurons+1 % +1 so the last row will be 2 or 8 (y)
        features(n, col) = mean(data(n, data(end-1, :) == i));
    end
    col = col + 1;
end

% assign labels: 0 for 2 Hz and 1 for 8 Hz
features(end, features(end, :) == 2) = 0;
features(end, features(end, :) == 8) = 1;

% separate features into X and y
X = features(1:end-1, :)';
y = features(end, :)';

% standardize the data
X = zscore(X); % <- need to take a further look at this

% cross validation before RFE
k = 5;
cv = cvpartition(y, 'KFold', k);

before_accuracies = zeros(1, k);
for i = 1:k
    trainIdx = cv.training(i);
    testIdx = cv.test(i);

    % partition into train and test sets
    XTrain = X(trainIdx, :);
    yTrain = y(trainIdx);
    XTest = X(testIdx, :);
    yTest = y(testIdx);

    % fit to an LDA model
    lda_model = fitcdiscr(XTrain, yTrain);
    yPred = predict(lda_model, XTest);

    before_accuracies(i) = sum(yPred == yTest) / numel(yTest);
end

before_accuracy = mean(before_accuracies) * 100;
fprintf('Cross validation accuracy before RFE: %.2f%%\n', before_accuracy);

% recursive feature elimination
opts = statset("Display", "off"); % "iter" or "off"
tf = sequentialfs(@myfun, X, y, "Direction", "backward", "Options", opts);
X = X(:, tf);

features_kept = sum(tf);
fprintf('Number of features kept in RFE: %d/%d (%d%)\n', features_kept, num_neurons, (features_kept/num_neurons)*100);

% cross validation after RFE
cv = cvpartition(y, 'KFold', k);

after_accuracies = zeros(1, 5);
for i = 1:k
    trainIdx = cv.training(i);
    testIdx = cv.test(i);

    % partition into train and test sets
    XTrain = X(trainIdx, :);
    yTrain = y(trainIdx);
    XTest = X(testIdx, :);
    yTest = y(testIdx);

    % fit to an LDA model
    lda_model = fitcdiscr(XTrain, yTrain);
    yPred = predict(lda_model, XTest);

    after_accuracies(i) = sum(yPred == yTest) / numel(yTest);
end

after_accuracy = mean(after_accuracies) * 100;
improvement = after_accuracy - before_accuracy;
fprintf('Accuracy improvement: %.2f%%\n', improvement);
fprintf('Cross validation accuracy after RFE: %.2f%%\n', after_accuracy);

fprintf('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n');

function criterion = myfun(XTrain, yTrain, XTest, yTest)
    mdl = fitcdiscr(XTrain, yTrain);
    yPred = predict(mdl, XTest);
    criterion = sum(yTest ~= yPred);  % count misclassified samples
end

diary off