%This is a program that plots pooled data from all fish
close all
clear all
%% Calculate the average firing rate of cells with more than 10 spikes in each fish
filenames1=dir('*_Stage1_output.mat');
filenames2=dir('*_Stage2_output.mat');
for i=1:length(filenames1)
    filenames1(i).name
    load(filenames1(i).name)
    load(filenames2(i).name)
    fr_list=[];
    avefr{i}
    for j=1:length(avefr)
        fr_list=[fr_list,avefr{j}(nuraw_count{j}>10)];
    end
    mean_fr{i}=mean(fr_list);
    sd_fr{i}=std(fr_list);
end
% histogram(fr_list,10,'Normalization','Probability')
% xlabel('Firing rate (Hz)');
% ylabel('Probability');
%% Plot the spike/ swim direction relation
Poptheta_rat=cell(5,6);
for i=1:length(filenames1)
    load(filenames1(i).name);
    for j=1:length(Theta_crct)
        if ~isempty(Theta_crct{j})
            NPOPtheta{i,j}=histcounts(Theta_interp_POP{j},linspace(0,180,5));%number of spikes at each engle
            Nalltheta{i,j}=histcounts(Theta_crct{j},linspace(0,180,5))/15;%time spent at each angle.
            Poptheta_rat{i,j}=NPOPtheta{i,j}./Nalltheta{i,j};%spike rate at each angle

        else
            Poptheta_rat{i,j}=NaN(1,4);
        end
    end
    tnum=length(cell2mat(Poptheta_rat(i,:)))/4;
    group=[];
    for k=1:4
        group=[group,k*ones(1,tnum)];
    end
    [P(i),A,stats]=kruskalwallis(reshape(cell2mat(Poptheta_rat(i,:)'),[1,4*tnum]),group);
    figure
    multcompare(stats)
end
%%
% for i=1:length(filenames1)
%     for j=1:length(Theta_crct)
%         if isempty(Poptheta_rat{i,j})
%             Poptheta_rat{i,j}=NaN(1,4);
%         end
%     end
%     Poptheta_rat_ave{i}=nanmean(cell2mat(Poptheta_rat(i,:)'));
%     Poptheta_rat_std{i}=nanstd(cell2mat(Poptheta_rat(i,:)'));
%     figure
%     bar(linspace(0,180,5),Poptheta_rat_ave{i})
%     hold on
%     errorbar(linspace(0,180,5),Poptheta_rat_ave{i},Poptheta_rat_std{i})
%     xlim([-10,190])
%     ylim([0,0.2])
% end
%% look at each cell individually
close all

for i=1:length(filenames1)%fish number
    m=1;
    load(filenames1(i).name);
    for j=1:length(Theta_crct)%trial number
        [a,b]=size(uname);
        for w=1:a
            if isempty(uname{w,j})
                uname{w,j}='NO';
            end
        end

        if ~isempty(Theta_crct{j})
            Nalltheta{i,j}=histcounts(Theta_crct{j},linspace(0,180,4))/15;%time spent at each angle for that trial
            for k=1:length(nuraw_count{j})%cell number, 
                if nuraw_count{j}(k)>10
                    Ntheta{i,j,k}=histcounts(Theta_interp{k,j},linspace(0,180,4));%number of spikes for each cell at a particular angle
                   
                    theta_rat{i,j,k}=Ntheta{i,j,k}./Nalltheta{i,j};%spike rate for each cell at a given angle
%                     figure
%                     bar(linspace(0,180,5),theta_rat{i,j,k})
%                     hold on
%                     title(strcat(uname{k,j},'trial',num2str(j)))
%                     xlim([-10,190])
%                     ylim([0,0.2])
                else
                    theta_rat{i,j,k}=NaN(1,3);
                end
                %%save the name of all units present    
                unitnames(m,:)=uname{k,j};
                unitcounts(m)=nuraw_count{j}(k);
                m=m+1;

            end

        else
            theta_rat{i,j,1}=NaN(1,3);
            Ntheta{i,j,k}=NaN(1,3);

        end
        
        %pool across trials for each unit
        
        
        %close all

    end
    %for each fish and each unit calculate a direction preference
    goodcells=unique(unitnames(unitcounts>10,2));
    temp=zeros(1,6);
    

    for w=1:length(goodcells)
        cellname=strcat('U',goodcells(w));
        index=strmatch(cellname,uname);
        [a,b]=ind2sub(size(uname),index);
        tempmat=zeros(1,3);
        totcnt=zeros(1,3);
        for z=1:length(a)
            if ~isnan(sum(cell2mat(Ntheta(i,b(z),a(z)))))&&~isempty(cell2mat(Ntheta(i,b(z),a(z))))
                tempmat=tempmat+cell2mat(Ntheta(i,b(z),a(z)));%sum all the counts for each bin for that cell across trials
                totcnt=totcnt+cell2mat(Nalltheta(i,b(z)));
            end
        end
        mean_theta(i,w,:)=tempmat./totcnt;%normalized by the amount of time the animal spent swimming at that angle
        norm_mean_theta(i,w,:)=mean_theta(i,w,:)/max(mean_theta(i,w,:));%normalize to the largest value so that the max is 1
        %caluclate mean and standard deviation across trials for a given cell, no SD is defined if the cell appears  only on one trial        
        figure
        bar([0,90,180],squeeze(norm_mean_theta(i,w,:)))
%         hold on
%         errorbar([0,90,180],squeeze(norm_mean_theta(i,w,:)),squeeze(std_theta(i,w,:)/max(std_theta(i,w,:))),'.')
%         xlim([-10,190])
%         ylim([0,0.2])
        xlabel('heading angle (degrees)')
        ylabel('spikes per second')
        title(strcat('fish',filenames1(i).name(1:8),cellname))
%         %saveas(gcf,strcat(filenames(i).name(1:8),uname(w),'DirTune.pdf'))
%         tnum=length(cell2mat(tempmat'))/3;
% 
%         group=[];
%         for k=1:3
%             group=[group,k*ones(1,tnum)];
%         end
%         [P2{i,w},A,stats]=kruskalwallis(reshape(cell2mat(tempmat),[1,tnum*3]),group);
%         figure
%         multcompare(stats)
        %for each cell , note the bin with largest amplitude and also %save
        %the p value
        [maxang, maxangindex]=max(mean_theta(i,w,:));
        MM=moment(norm_mean_theta(i,w,:),2);
        Theta_stats{i,w}=[maxangindex,length(a),MM];
        
    end
    close all
    clear unitnames;
    clear unitcounts;
    clear m;
end
%%calculate the percentage of direction tuned cells
%Theta_stats first column the direction at which largest number of spikes
%were found.1=forward,2=turn,3=backward
thetastats=cell2mat((reshape(Theta_stats,[35,1])));
bkspkfrac=sum(thetastats(:,1)==3)/length(thetastats(:,1))
bkspk=thetastats(find(thetastats(:,1)==3),3);
bkspkmeMM=mean(bkspk)
bkspksdMM=std(bkspk)
frspkfrac=sum(thetastats(:,1)==1)/length(thetastats(:,1))
frspk=thetastats(find(thetastats(:,1)==1),3);
frspkmeMM=mean(frspk)
frspksdMM=std(frspk)
mdspkfrac=sum(thetastats(:,1)==2)/length(thetastats(:,1))
mdspk=thetastats(find(thetastats(:,1)==2),3);
mdspkmeMM=mean(mdspk)
mdspksdMM=std(mdspk)
%save('SwimDirSpkStats.mat','bkspkfrac','bkspkmeMM','bkspksdMM','frspkfrac','frspkmeMM', 'frspksdMM', 'mdspkfrac','mdspkmeMM','mdspksdMM')
%%estimate weather the peakiness of the distribution is generally more for
%%the cells that prefer backward swims or turns or forward
[p,a,stats]=kruskalwallis(thetastats(:,3),thetastats(:,1))
% multcompare(stats)
figure
title(strcat('backward swim 2nd moment',num2str(100*bkspkfrac),'%'))
hold on
histogram(bkspk,10,'Normalization','Probability')
%saveas(gcf,'spkatbackward.pdf')

figure
title(strcat('forward swim 2nd moment',num2str(100*frspkfrac),'%'))
hold on
histogram(frspk,10,'Normalization','Probability')
%saveas(gcf,'spkatforward.pdf')

figure
title(strcat('turn swim 2nd moment',num2str(100*mdspkfrac),'%'))
hold on
histogram(mdspk,10,'Normalization','Probability')
%saveas(gcf,'spkatturn.pdf')

%%
%plot the pooled histogram of the distance from the max spike triggered
%landmark

filenames4=dir('*_Stage9_output.mat');
SpkTrigLMPool=[];
for i=1:length(filenames4)
    load(filenames4(i).name);
    SpkTrigLMPool=[SpkTrigLMPool,dist_cm(:)'];
end
histogram(cell2mat(SpkTrigLMPool(:)),10,'Normalization','Probability')
%saveas(gcf,'SpikeTriggeredLMPoolDisttoCent.pdf')