% extraction des data pour une waveform analysis puis k-mean methods pour
% clusteriser les differents de neurons


clear all; clc; close all;
load('RearlyTraining.mat')


WFlength=800; %total duration of the waevform 800microseconds
bins=25;%bins=25microsecond
time=(bins:bins:WFlength);
runAnalysis=1;


for i=1:size(R.WF,1)
    if sum(R.WF(i,:))==0
        R.WFanalysis(i,2)=0;
        R.WFanalysis(i,3)=0;
    else
        [maxValue, MaxIndex] = max(R.WF(i,:));
        [minValue, MinIndex] = min(R.WF(i,:));
        
        PeakValleyW=abs(time(MaxIndex)-time(MinIndex));
        PeakValleyM=abs(maxValue-minValue);
        R.WFanalysis(i,2)=PeakValleyM;
        
        idx = find(diff(R.WF(i,:) >= minValue/2));
        x2 = time(idx) + (minValue/2 - R.WF(i,idx)) .* (time(idx+1) - time(idx)) ./ (R.WF(i,idx+1) - R.WF(i,idx));
        if length(x2)==1
            R.WFanalysis(i,3)=NaN;
        else
            ValleyFWHM=x2(2)-x2(1);
            R.WFanalysis(i,3)=ValleyFWHM;
        end
    end
end

selection=R.WFanalysis(:,2)~=0;
X=R.WFanalysis(selection,:);

idx(1:sum(selection),1)=1;
idx(X(:,1)>20,1)=2; % threshold 20Hz to separate MSN from FSI
idx(idx(:,1)==2 & X(:,3)>150,1)=0; 
idx(X(:,1)<20 & X(:,1)>12.5 & idx(:,1)==1,1)=0; % intermediate frequencies = not classified


for Y=1:size(idx,1)
    if idx(Y,1)==1
        color(Y,:)=[255/255 0/255 0/255];%MSN
        sizeplot(Y,1)=50;
    elseif idx(Y,1)==2
        color(Y,:)=[0 255/255 0/255];%FSI
        sizeplot(Y,1)=50;
    elseif idx(Y,1)==3
        color(Y,:)=[0/255 0/255 255/255]; %TAN
        sizeplot(Y,1)=50;
    elseif idx(Y,1)==0
        color(Y,:)=[0/255 0/255 0/255];%unclassified
        sizeplot(Y,1)=5;
    end
end


figure;
subplot(3,4,[1 2 3 5 6 7 9 10 11])
scatter(X(:,1),X(:,3),sizeplot,color)
hold on
title(strcat('all neurons clustering kmean'))
xlabel('FR') % x-axis label
ylabel('half-width')% y-axis label
hold off

Celltype(1:length(selection),1)=R.SESSION;
Celltype(1:length(selection),2)=0;
Celltype(selection,2)=idx(:,1);

save('Celltype_earlyTraining.mat','Celltype')