
%% PC Antonio

clc
clear


dirData = 'Y:\OptoMECLEC\';
animal  = {'OML18','OML18','OML18','OML18','OML19','OML19'};
sessions  = {'day1','day2','day4','day5','day2','day3'};
folder='Metrics';

ses=1;
dirses = [dirData animal{ses} '\'  sessions{ses}];
cd(dirses);
filename=[animal{ses} sessions{ses}]
load([filename '.mat'])


dir = dirses;
dirsave = 'PCsAnt\';
load([sessions{ses} '_Pos.mat']);
%%
% speed thrs
for i = 1:length(Pos.ON)
    if Pos.ON(i,3) < 0.2
       Pos.ON(i,2) = NaN;
    end
end
for i = 1:length(Pos.Off)
    if Pos.Off(i,3) < 0.2
       Pos.Off(i,2) = NaN;
    end
end   
 figure; plot(Pos.ON(:,1),Pos.ON(:,2),'.b');hold on;  plot(Pos.Off(:,1),Pos.Off(:,2),'.r');hold on;
posTrials{1} = Pos.ON(:,1:2);  
posTrials{2} = Pos.Off(:,1:2);  

% get PCs
for unit = 1:length(spikes.times)
    for t = 1:2
    [curve{unit}{t},stats{unit}{t}] = FiringCurve(posTrials{t},spikes.times{unit},'nBins',100,'threshold',0.2,'minSize',10,'minPeak',1,'smooth',3); 
    end
end

%2- Refine PFs 
secondaryPF = 0.6; %secondary PF should have peak FR hihger than secondayPF*peakFR of primary PF
trackCm = {'0','50','100','150','200'}; 
for i=1:length(spikes.times)
    figure('visible','off');
    subplot(1,2,1);plot(curve{i}{1}.x,curve{i}{1}.rate);title('ON');ylabel('FR(Hz)');xlabel('track (cm)');hold on;
    set(gca,'XTick',[0 0.25 0.5 0.75 1],'XTickLabel',trackCm);
    if stats{i}{1}.peak ~= 0
       for k=1:length(stats{i}{1}.peak) % for each PF of this cell
           if stats{i}{1}.peak(k)>stats{i}{1}.peak(1)*secondaryPF % destroy small secondary PFs
               peakLoc(k)=curve{i}{1}.x(stats{i}{1}.x(k));
               scatter(peakLoc(k),stats{i}{1}.peak(k),'MarkerFaceColor','r','LineWidth',1);
               fieldIni(k)=curve{i}{1}.x(stats{i}{1}.fieldX(k,1));
               fieldFin(k)=curve{i}{1}.x(stats{i}{1}.fieldX(k,2));
               x1=fieldIni(k);y1=get(gca,'ylim');hold on; plot([x1 x1],y1,'-.k','LineWidth',1);
               x2=fieldFin(k);y1=get(gca,'ylim');hold on; plot([x2 x2],y1,'-.k','LineWidth',1);
           else
               stats{i}{1}.x(k) = NaN;      stats{i}{1}.y(k) = NaN; 
               stats{i}{1}.field(k) = 0;    stats{i}{1}.size(k) = 0; 
               stats{i}{1}.peak(k) = 0;     stats{i}{1}.mean(k) = 0;  
               stats{i}{1}.fieldX(k,1:2) = [NaN,NaN]; stats{i}{1}.fieldY(k,1:2) = [NaN,NaN];                 
           end
       end
       for k=1:length(stats{i}{1}.peak)% destroy too large PF (CHANGE for recalculate!!)
           if stats{i}{1}.peak(k)~=0 && ... 
                   curve{i}{1}.x(stats{i}{1}.fieldX(k,2))-curve{i}{1}.x(stats{i}{1}.fieldX(k,1)) > 0.5;
               stats{i}{1}.x(k) = NaN;      stats{i}{1}.y(k) = NaN; 
               stats{i}{1}.field(k) = 0;    stats{i}{1}.size(k) = 0; 
               stats{i}{1}.peak(k) = 0;     stats{i}{1}.mean(k) = 0;  
               stats{i}{1}.fieldX(k,1:2) = [NaN,NaN]; stats{i}{1}.fieldY(k,1:2) = [NaN,NaN];   
           end     
       end
       clear temp peakLoc fieldIni fieldFin x1 x2 y1 y2 k
    end
    subplot(1,2,2);plot(curve{i}{1}.x,curve{i}{2}.rate);title('OFF');ylabel('FR(Hz)');xlabel('track (cm)');hold on;
    set(gca,'XTick',[0 0.25 0.5 0.75 1],'XTickLabel',trackCm);
    if stats{i}{2}.peak ~= 0
       for k=1:length(stats{i}{2}.peak) 
           if stats{i}{2}.peak(k)>stats{i}{2}.peak(1)*secondaryPF
               peakLoc(k)=curve{i}{2}.x(stats{i}{2}.x(k));
               scatter(peakLoc(k),stats{i}{2}.peak(k),'MarkerFaceColor','r','LineWidth',1);
               fieldIni(k)=curve{i}{2}.x(stats{i}{2}.fieldX(k,1));
               fieldFin(k)=curve{i}{2}.x(stats{i}{2}.fieldX(k,2));
               x1=fieldIni(k);y1=get(gca,'ylim');hold on; plot([x1 x1],y1,'-.k','LineWidth',1);
               x2=fieldFin(k);y1=get(gca,'ylim');hold on; plot([x2 x2],y1,'-.k','LineWidth',1);   
           else % destroy too small secondary fields
               stats{i}{2}.x(k) = NaN;      stats{i}{2}.y(k) = NaN; 
               stats{i}{2}.field(k) = 0;    stats{i}{2}.size(k) = 0; 
               stats{i}{2}.peak(k) = 0;     stats{i}{2}.mean(k) = 0;  
               stats{i}{2}.fieldX(k,1:2) = [NaN,NaN]; stats{i}{2}.fieldY(k,1:2) = [NaN,NaN]; 
           end  
       end
       for k=1:length(stats{i}{2}.peak) 
           if stats{i}{2}.peak(k)~=0 && ... % destroy too large PF (CHANGE for recalculate!!)
               curve{i}{2}.x(stats{i}{2}.fieldX(k,2))-curve{i}{2}.x(stats{i}{2}.fieldX(k,1)) > 0.5;
               stats{i}{2}.x(k) = NaN;      stats{i}{2}.y(k) = NaN; 
               stats{i}{2}.field(k) = 0;    stats{i}{2}.size(k) = 0; 
               stats{i}{2}.peak(k) = 0;     stats{i}{2}.mean(k) = 0;  
               stats{i}{2}.fieldX(k,1:2) = [NaN,NaN]; stats{i}{2}.fieldY(k,1:2) = [NaN,NaN];   
           end    
       end
       clear temp peakLoc fieldIni fieldFin x1 x2 y1 y2 k
    end
    suptitle(['cell' num2str(i)]);
    saveas(gcf,[dirsave 'cell' num2str(i)],'png'); close all;
end
save([dirsave 'PCs.mat'],'curve','stats'); 
 
% boundaries of each PF
for i=1:length(spikes.times)
    for j=1:2
        for k=1:length(stats{i}{j}.peak) 
            if stats{i}{j}.peak(k) ~= 0
                boundaries{i}{j}(k,1)= curve{i}{j}.x(stats{i}{j}.fieldX(k,1));
                boundaries{i}{j}(k,2)= curve{i}{j}.x(stats{i}{j}.fieldX(k,2));
            else
                boundaries{i}{j}(k,1)= NaN;
                boundaries{i}{j}(k,2)= NaN;
            end
        end   
    end
end

% count nPFs
count = zeros(2,1); % one cell per each PF 
for i=1:length(spikes.times)
   for j=1:2
          for k=1:length(stats{i}{j}.x) %size(boundaries{i}{j},1)
              if ~isnan (stats{i}{j}.x(k))%(boundaries{i}{j}(k,1)) % for each PF not removed
                count(j)=count(j)+1;
              end
          end
   end
end
totalPFs = count;

% thetaP from Farnaz -- better keep only maze part
thetaP(:,1) = Theta.Phase.timestamps;
thetaP(:,2) = Theta.Phase.phase;
save('thetaP.mat','thetaP');

countON = 0; % per each PF individually
for i=1:length(spikes.times)
   for j=1
          for k=1:length(stats{i}{j}.x) %size(boundaries{i}{j},1)
              if ~isnan (stats{i}{j}.x(k))%(boundaries{i}{j}(k,1)) % for each PF not removed
                countON=countON+1;
                [dataPPon{countON},statsPPon{countON}] = PhasePrecession(posTrials{j},spikes.times{i},thetaP,'boundaries',boundaries{i}{j}(k,:));
              end
          end
   end
end
countOFF = 0;
for i=1:length(spikes.times)
   for j=2
          for k=1:length(stats{i}{j}.x) %size(boundaries{i}{j},1)
              if ~isnan (stats{i}{j}.x(k))%(boundaries{i}{j}(k,1)) % for each PF not removed
                countOFF=countOFF+1;
                [dataPPoff{countOFF},statsPPoff{countOFF}] = PhasePrecession(posTrials{j},spikes.times{i},thetaP,'boundaries',boundaries{i}{j}(k,:));
              end
          end
   end   
end
save([dirsave 'PhasePrec.mat'],'dataPPon','statsPPon','dataPPoff','statsPPoff');

% PLOT phase precession 
for i=1:totalPFs(1)
   figure('visible','on');
   PlotPhasePrecession(dataPPon{i},statsPPon{i});
   suptitle(['cell ' num2str(i)]);
   saveas(gcf,[dirsave 'ON' 'cell' num2str(i)],'png'); 
   close all;       
end
for i=1:totalPFs(2)
   figure('visible','on');
   PlotPhasePrecession(dataPPoff{i},statsPPoff{i});
   suptitle(['cell ' num2str(i)]);
   saveas(gcf,[dirsave 'OFF' 'cell' num2str(i)],'png'); 
   close all;       
end

% Clasify PF in preccesing and no preccesing 
for i=1:totalPFs(1)
    if statsPPon{i}.p < 0.05 && statsPPon{i}.slope < 0
       precPFon(i,1) = 1; 
    else 
       precPFon(i,1) = 0;
    end
end
for i=1:totalPFs(2)
    if statsPPoff{i}.p < 0.05 && statsPPoff{i}.slope < 0
       precPFoff(i,1) = 1; 
    else 
       precPFoff(i,1) = 0;
    end
end
fracPP(1) = sum(precPFon)/totalPFs(1);
fracPP(2) = sum(precPFoff)/totalPFs(2);
figure;
subplot(1,2,1);bar(totalPFs);title('nPFs');set(gca,'XTickLabel',{'ON','OFF'});
subplot(1,2,2);bar(fracPP);title('% preces PFs');set(gca,'XTickLabel',{'ON','OFF'});

save([dirsave 'nPFs.mat'],'totalPFs','fracPP');

%% pool
dir = 'A:\OptoMECLEC\';
sessions = {'OML18\day1\';'OML18\day2\';'OML18\day4\';'OML19\day2\';'OML19\day3\'};
ppON = []; ppOFF = [];
for ses = 1:length(sessions)
    load([dir sessions{ses} 'PCsAnt\PhasePrec']);
    ppON = cat(2,ppON,statsPPon);
    ppOFF = cat(2,ppOFF,statsPPoff);
    clear statsPPoff statsPPon
end

% Clasify PF in preccesing and no preccesing 
for i=1:length(ppON)
    if ppON{i}.p < 0.05 && ppON{i}.slope < 0
       precPFon(i,1) = 1; 
    else 
       precPFon(i,1) = 0;
    end
end
for i=1:length(ppOFF)
    if ppOFF{i}.p < 0.05 && ppOFF{i}.slope < 0
       precPFoff(i,1) = 1; 
    else 
       precPFoff(i,1) = 0;
    end
end
fracPP(1) = sum(precPFon)/length(ppON);
fracPP(2) = sum(precPFoff)/length(ppOFF);
figure;
subplot(1,2,1);bar([length(ppON) length(ppOFF)]);title('nPFs');set(gca,'XTickLabel',{'ON','OFF'});
subplot(1,2,2);bar(fracPP);title('% preces PFs');set(gca,'XTickLabel',{'ON','OFF'});

% pool stats
slope = cell(2,2); intercept = cell(2,1); r2 = cell(2,1);
for i = 1:length(ppON)
    if ppON{i}.p < 0.05 && ppON{i}.slope < 0
    slope{1} = cat(1,slope{1},ppON{i}.slope);
    intercept{1} = cat(1,intercept{1},ppON{i}.intercept);
    r2{1} = cat(1,r2{1},ppON{i}.r2);
    end
end
for i = 1:length(ppOFF)
    if ppOFF{i}.p < 0.05 && ppOFF{i}.slope < 0
    slope{2} = cat(1,slope{2},ppOFF{i}.slope);
    intercept{2} = cat(1,intercept{2},ppOFF{i}.intercept);
    r2{2} = cat(1,r2{2},ppOFF{i}.r2);
    end
end

figure;
for i = 1:2
subplot(1,2,i);boxplot(r2{i});ylim([0 0.3]);
end
figure;
for i = 1:2
subplot(1,2,i);boxplot(slope{i});ylim([-30 0]);
end
pS = ranksum(slope{1},slope{2})
figure;
for i = 1:2
subplot(1,2,i);boxplot(rad2deg(mod(intercept{i},2*pi)));ylim([0 350]);
end




