clear
close all
load noduli_per_master


% Caricati: n e p (negativi e positivi)

Nn = size(n, 1);
Np = size(p, 1);

% LE FEATURE SONO:
% volume - raggio mm - sfericita' - intensita' media (colonne 1 2 3 4)

% Istogrammi

% VOLUME
xmin = min([min(n(:,1)) ; min(p(:,1))]); % il minimo valore di volume, considerando n e p
xmax = max([max(n(:,1)) ; max(p(:,1))]); % il massimo....
centers = linspace(xmin, xmax, 11);
hn = hist(n(:,1), centers);
hp = hist(p(:,1), centers);
subplot(1, 3, 1), plot(centers, hn / Nn, 'b-o', centers, hp / Np, 'r-o')
title('VOLUME')

% RAGGIO
xmin = min([min(n(:,2)) ; min(p(:,2))]);
xmax = max([max(n(:,2)) ; max(p(:,2))]);
centers = linspace(xmin, xmax, 11);
hn = hist(n(:,2), centers);
hp = hist(p(:,2), centers);
subplot(1, 3, 2), plot(centers, hn / Nn, 'b-o', centers, hp / Np, 'r-o')
title('RAGGIO')

% SFERICITA'
xmin = min([min(n(:,3)) ; min(p(:,3))]);
xmax = max([max(n(:,3)) ; max(p(:,3))]);
centers = linspace(xmin, xmax, 11);
hn = hist(n(:,3), centers);
hp = hist(p(:,3), centers);
subplot(1, 3, 3), plot(centers, hn / Nn, 'b-o', centers, hp / Np, 'r-o')
title('SFERICITA')

%%%%%%%%%%%%%%%%%%%%

% - SCATTER PLOT - SCATTER PLOT - SCATTER PLOT - SCATTER PLOT - SCATTER PLOT

figure

%  volume vs raggio
subplot(1, 3, 1), plot(n(:,1), n(:,2), 'b.', p(:,1), p(:,2), 'r.'), xlabel('VOLUME'), ylabel('RAGGIO')
%  volume vs sfericita'
subplot(1, 3, 2), plot(n(:,1), n(:,3), 'b.', p(:,1), p(:,3), 'r.'), xlabel('VOLUME'), ylabel('SFERICITA')
%  raggio vs sfericita'
subplot(1, 3, 3), plot(n(:,2), n(:,3), 'b.', p(:,2), p(:,3), 'r.'), xlabel('RAGGIO'), ylabel('SFERICITA')

% SCELGO LE FEATURES: volume - raggio mm - sfericita'  (colonne 1 2 3)
% (nel file c'e' anche l'intensita' media)

% NOTARE CHE E' STATO APPLICATO GIA' UN FILTRO SUL VOLUME, MINIMO E MASSIMO
% (tra 30 e 1700 circa) in base ai valori ricavati dai positivi

% Ora il dataset va diviso in due parti, una usata per l'allenamento del
% classificatore, e una per il test

% Prima mescoliamo i gruppi

nIdx = randperm(Nn);
n = n(nIdx, :);

pIdx  = randperm(Np);
p  = p(pIdx, :);

% calcoliamo il numero di positivi diviso 2, e ugualmente per i negativi

Npd2  = floor(Np/2);
Nnd2  = floor(Nn/2);

% Ora definiamo il set di allenamento e quello di test (hold out):

nTr = n(1:Nnd2,:);       % feature negativi di training
pTr = p(1:Npd2,:);       % feature positivi di training
nTe = n(Nnd2+1:end,:);   % feature negativi di test
pTe = p(Npd2+1:end,:);   % feature positivi di test

gTr = [zeros(Nnd2,1)    ; ones(Npd2,1)];      % gruppi (etichette o target) di training
gTe = [zeros(Nn-Nnd2,1) ; ones(Np-Npd2,1)];   % gruppi (etichette o target) di test

% metto insieme i dati

x       = [nTr;pTr];      % feature di train (matrice N x 4)
xtest   = [nTe;pTe];      % feature di test  (idem)
lab     = gTr;            % label/target di train
labtest = gTe;            % label/target di test

% Classificazione

if 0  % 1 = classificatore a discriminante (lineare)
    
    y = classify(xtest, x, lab, 'linear');
    
else  % rete neurale, devo trasporre!
    x = x';          % adesso ogni riga contiene una feature per tutti gli oggetti
    lab = lab';
    xtest = xtest';  % adesso ogni riga contiene una feature per tutti gli oggetti
    
    net=patternnet(5,'trainrp');   % 5 neuroni nascosti
    %net=patternnet([],'trainrp');   % nessun  neurone nascosto, rete lineare
    
    net.trainParam.epochs = 1000;  % numero massimo di epoche di allenamento (criterio di arresto)
    net.trainParam.goal = 1e-10;
    net.trainParam.min_grad = 1e-50;

    net.divideParam.trainRatio = 0.9;
    net.divideParam.valRatio   = 0.1;
    net.divideParam.testRatio  = 0;

    net = train(net, x, lab);  % training
    y = net(xtest);   % applico la rete allenata ai dati di test
    y = y';        % costretto a trasporre per tornare alla situazione di vettore colonna 
    lab = lab';
    
    [xROC, yROC, tROC, AUC, optrocpt] = perfcurve(labtest, y, 1);  % calcolo della curva roc
    figure, plot(xROC,yROC);
    disp(['AUC = ' num2str(AUC)])
    
    % Youden's index
    % https://en.wikipedia.org/wiki/Youden's_J_statistic
    
    youd = yROC - xROC;
    [maxYoud, maxYoudIdx] = max(youd);
    tROCYoud = tROC(maxYoudIdx);
    xROCYoud = xROC(maxYoudIdx);
    yROCYoud = yROC(maxYoudIdx);
    hold on
    plot(xROCYoud, yROCYoud, 'ro')
    hold off
    
    y = y > tROCYoud;
end

% Ora confrontiamo i gruppi noti a priori labtest con le classi generate y

%  gruppo classe
%  1      1     TP
%  1      0     FN
%  0      1     FP
%  0      0     TN

TP = sum(    labtest  .*     y)
FN = sum(    labtest  .* not(y))
FP = sum(not(labtest) .*     y)
TN = sum(not(labtest) .* not(y))

sensitivity = TP / (TP + FN);
specificity = TN / (TN + FP);

accuracy = (TP + TN) / (TP + FN + FP + TN);

disp(['sensitivity = ' num2str(sensitivity)])
disp(['specificity = ' num2str(specificity)])
disp(['accuracy = ' num2str(accuracy)])
