%This is a program that calculates the spatial information content for each
%cell and determines whether or not it is statistically significant
%comapring it with randomly shuffled spikes.
clear all
close all
exp_date=input('please type in the experiment date : ','s');
load(strcat(exp_date,'_Stage1_output.mat'));
load(strcat(exp_date,'_Stage2_output.mat'));
y=0:80:1600;
x=0:80:1200;
resy=1600/80;
resx=1200/80;
%% Calculate the information content of each cell
for tn=1:length(filename)
    ndata=length(nuraw_count{tn});
    for w=1:ndata
        temp=0;
        for m=1:resx
            for n=1:resy
                if avefrbin{tn}(m,n,w)>0 
                    temp=temp+occup{tn}(m,n,w)*(avefrbin{tn}(m,n,w)/avefr{tn}(w))*log2(avefrbin{tn}(m,n,w)/avefr{tn}(w));
                end
            end
        end
        Infcont{tn}(w)=temp;
        if tn~=1 
        % calculate skewness of the heading angle
        %plot polar histogram of angles at the time of spike
        [a1,e]=histcounts(Theta_interp{w,tn},linspace(0,180,4),'Normalization','count');
        [a2,e]=histcounts(Theta_crct{tn},linspace(0,180,4),'Normalization','count');
        FR_Theta_norm{w,tn}=a1./(a2/15);
        skew_theta{tn}(w)=skewness(FR_Theta_norm{w,tn});

        histogram('BinEdges',e,'BinCounts',FR_Theta_norm{w,tn})
        title(strcat(filename{tn},'_',uname{w,tn}))
        %saveas(gcf,strcat('ThetaHist',filename{tn},'_',uname{w,tn},'.pdf'))
        bin_cent=e(1:end-1)+mean(diff(e))/2;
        all_vect=[bin_cent;FR_Theta_norm{w,tn}];
        [X,Y]=pol2cart(all_vect(1,:),all_vect(2,:));
        Z=[X',Y'];
        Norm_VectAve=sum(Z)/sum(FR_Theta_norm{w,tn});
        [TH,Rvec_l{w,tn}]=cart2pol(Norm_VectAve(1),Norm_VectAve(2));
        Ave_TH{tn}(w)=180*TH/pi;
        end
    %%do the same calculation for shuffled spikes
    [nurawshuffcnt,nurawshufftmp]=hist(nutimes{w,tn},time2{tn});
    %shuffle the location of the spikes

    for i=1:1000
         rng('shuffle') 
%         vt = rand(size(time2{tn}));
%         spikes = (avefr{tn}(w)*0.002) > vt;
%         poissspktimes=time2{tn}(spikes==1);
%         poissspkcnt(i)=length(poissspktimes);
%         shuffspktimes=poissspktimes;
        
        shift_time=30:10:time2{tn}(end)-30;
        rnd_shift_time=shift_time(randi(length(shift_time)));
        %duration of the recording -10s
        nbins_shift=rnd_shift_time*fs;
        shuff=circshift(nurawshuffcnt,nbins_shift);
        shuffspktimes=nurawshufftmp(shuff==1);
        if tn~=1
            Theta_interp_shuff=interp1(xvalues{tn},Theta_crct{tn},shuffspktimes,'spline');
            [a1,e]=histcounts(Theta_interp_shuff,linspace(0,180,4),'Normalization','count');
            [a2,e]=histcounts(Theta_crct{tn},linspace(0,180,4),'Normalization','count');
            FR_Theta_norm_shuff=a1./(a2/15);
            skew_theta_shuff{tn}(w,i)=skewness(FR_Theta_norm_shuff);      
        end
        for k=1:2
         Place_int_unit_shuff{tn,i,w,k}=interp1(xvalues{tn},hd_crct{tn}(:,k),shuffspktimes,'spline');
        end            
    end
%      if w==1
%         [poprawshuffcnt,poprawshufftmp]=hist(poptimes{tn}, time2{tn});
%         for i=1:1000
%             rng('shuffle');
%             shuffpop=poprawshuffcnt(randperm(length(poprawshuffcnt),length(poprawshuffcnt)));
%             shuffspktimespop=poprawshufftmp(shuffpop==1);
%             for k=1:2
%                 Place_int_pop_shuff{tn,i,k}=interp1(xvalues{tn},hd_crct{tn}(:,k),shuffspktimespop,'spline');
%             end            
%         end
%      end
    clear nurawshuffcnt
    clear nurawshufftmp

    end
    spkcnt_shuff{tn}=zeros(resx,resy,ndata,1000);
    spkcntpop_shuff{tn}=zeros(resx,resy,1000);
    avefrbin_shuff{tn}=NaN(resx,resy,ndata,1000);
    avefrbinpop{tn}=NaN(resx,resy,1000);

   %calculate the suffled firing rate per bins     
    for z=1:1000
        for w=1:ndata
            for k=1:length(Place_int_unit_shuff{tn,z,w,1})
                for i=1:resx
                    for j=1:resy
                        if (Place_int_unit_shuff{tn,z,w,2}(k)>=x(i)&&Place_int_unit_shuff{tn,z,w,2}(k)<x(i+1)&& Place_int_unit_shuff{tn,z,w,1}(k)>=y(j)&&Place_int_unit_shuff{tn,z,w,1}(k)<y(j+1))
                            spkcnt_shuff{tn}(i,j,w,z)=spkcnt_shuff{tn}(i,j,w,z)+1;
                        end

                    end
                end
            end
            %remove isolated bins 
            for ii=1:resx
                for jj=1:resy
                    if spkcnt_shuff{tn}(ii,jj,w)>0
                        if ii>1 && jj>1 && ii<15
                            temp=nansum([spkcnt_shuff{tn}(ii-1,jj-1,w,z),spkcnt_shuff{tn}(ii,jj-1,w,z),spkcnt_shuff{tn}(ii+1,jj-1,w,z),spkcnt_shuff{tn}(ii-1,jj,w,z),spkcnt_shuff{tn}(ii+1,jj,w,z),spkcnt_shuff{tn}(ii-1,jj+1,w,z),spkcnt_shuff{tn}(ii,jj+1,w,z),spkcnt_shuff{tn}(ii+1,jj+1,w,z)]);
                        elseif ii==1 && jj>1
                            temp=nansum([spkcnt_shuff{tn}(ii,jj-1,w,z),spkcnt_shuff{tn}(ii+1,jj-1,w,z),spkcnt_shuff{tn}(ii+1,jj,w,z),spkcnt_shuff{tn}(ii,jj+1,w,z),spkcnt_shuff{tn}(ii+1,jj+1,w,z)]);
                        elseif ii>1 && jj==1
                            temp=nansum([spkcnt_shuff{tn}(ii-1,jj,w,z),spkcnt_shuff{tn}(ii+1,jj,w,z),spkcnt_shuff{tn}(ii-1,jj+1,w,z),spkcnt_shuff{tn}(ii,jj+1,w,z),spkcnt_shuff{tn}(ii+1,jj+1,w,z)]);
                        elseif ii==15
                            temp=nansum([spkcnt_shuff{tn}(ii-1,jj-1,w,z),spkcnt_shuff{tn}(ii,jj-1,w,z),spkcnt_shuff{tn}(ii-1,jj,w,z),spkcnt_shuff{tn}(ii-1,jj+1,w,z),spkcnt_shuff{tn}(ii,jj+1,w,z)]);
                        end
                        if temp<=0
                           spkcnt_shuff{tn}(ii,jj,w,z)=NaN;
                        end
                    end
                end
            end
            %do the same for population
%             if w==1
%                 for k=1:length(Place_int_pop_shuff{tn,z,1})
%                     for i=1:resx
%                         for j=1:resy
%                             if (Place_int_pop_shuff{tn,z,2}(k)>=x(i)&&Place_int_pop_shuff{tn,z,2}(k)<x(i+1)&& Place_int_pop_shuff{tn,z,1}(k)>=y(j)&&Place_int_pop_shuff{tn,z,1}(k)<y(j+1))
%                                 spkcntpop_shuff{tn}(i,j,z)=spkcntpop_shuff{tn}(i,j,z)+1;
%                             end
% 
%                         end
%                     end
%                 end
%             end
        for m=1:resx
            for n=1:resy
                if  isnan(visitd{tn}(m,n)) || visitd{tn}(m,n)<5 
                    spkcntpop_shuff{tn}(m,n,z)=NaN;
                    spkcnt_shuff{tn}(m,n,w,z)=NaN;
                    
                end
            end
        end

        end
    end
    
    %calculate the avefrbin for each cell and each random shuffle
    for z=1:1000
        for w=1:ndata
            timespent_unit_shuff=timespent{tn}(:,:).*(~isnan(spkcnt_shuff{tn}(:,:,w,z)));
            total_time_shuff=nansum(nansum(timespent_unit_shuff));
            occup_shuff{tn}(:,:,w,z)=timespent_unit_shuff/total_time_shuff;
            avefr_shuff{tn}(w,z)=nansum(nansum(spkcnt_shuff{tn}(:,:,w,z)))/total_time_shuff;
            avefrbin_shuff{tn}(:,:,w,z)=spkcnt_shuff{tn}(:,:,w,z)./ timespent_unit_shuff;
            DiffMat{tn}(:,:,w,z)=(avefrbin_shuff{tn}(:,:,w,z)/max(max(avefrbin_shuff{tn}(:,:,w,z))))-(avefrbin{tn}(:,:,w)/max(max(avefrbin{tn}(:,:,w))));
            DiffMat{tn}(:,:,w,z)=DiffMat{tn}(:,:,w,z)/(resx*resy-sum(sum(isnan(DiffMat{tn}(:,:,w,z)))));
            sum_DiffMat{tn}(w,z)=nansum(nansum(DiffMat{tn}(:,:,w,z)));
            
        end
    end
    for w=1:ndata
        Stats_shuff_DiffMat{tn}(w)=signrank(sum_DiffMat{tn}(w,:));
    end

    %% calculate the difference between the shuffled image and the actual image
    
    
    %%calculate information content for the shuffled data and compare with
    %%the actual data
    temp=0;
    Infcont_shuff{tn}=zeros(ndata,1000);

    for i=1:ndata
        for z=1:1000
            for m=1:resx
                for n=1:resy
                    if avefrbin_shuff{tn}(m,n,i,z)>0 
                        temp=temp+occup_shuff{tn}(m,n,i,z)*(avefrbin_shuff{tn}(m,n,i,z)/avefr_shuff{tn}(i,z))*log2(avefrbin_shuff{tn}(m,n,i,z)/avefr_shuff{tn}(i,z));
                    end
                end
            end
            Infcont_shuff{tn}(i,z)=temp;
            temp=0;
        end
    end
    if tn~=1
    %find which cells show significant skewness
        skewshuffmean{tn}=mean(skew_theta_shuff{tn},2);
        skewshuffstd{tn}=std(skew_theta_shuff{tn},[],2);
    end
    
%find which cell has significant place information:
    shuffmean{tn}=mean(Infcont_shuff{tn},2);
    shuffstd{tn}=std(Infcont_shuff{tn},[],2);
    for i=1:ndata
        [H{tn}(i),P{tn}(i)]=ztest(Infcont{tn}(i),shuffmean{tn}(i),shuffstd{tn}(i),'alpha',0.05);
        if tn~=1
         [Hsk{tn}(i),Psk{tn}(i)]=ztest(skew_theta{tn}(i),skewshuffmean{tn}(i),skewshuffstd{tn}(i),'alpha',0.05);
        end
    end
    %do the same for population
%     shuffmeanpop{tn}=mean(Infcontpop_shuff{tn});
%     shuffstdpop{tn}=std(Infcontpop_shuff{tn});

%     [Hpop{tn},Ppop{tn}]=ztest(Infocontpop(tn),shuffmeanpop{tn},shuffstdpop{tn},'alpha',0.05);
%     placecellpop{tn}=Hpop{tn}>0
end
for tn=1:length(filename)
    placecell{tn}=H{tn}>0
    skewedcell{tn}=Hsk{tn}>0
end

%save(strcat(exp_date,'_Stage4_output.mat'),'placecell','P','Infcont','Infcont_shuff','Ave_TH','skew_theta','skewedcell','Psk','Stats_shuff_DiffMat');