%This is a program to find out if the spiking pattern of given units
%changes over time in a 10 cm radius around landmarks.
clearvars
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;
ang=0:0.01:2*pi;
fs=500;%sampling rate for the position of the fish (interpolated)
time_spent_stp2=zeros(1,length(filename));
spkcnt_stp2=zeros(size(uname));
fr_stp2=zeros(size(uname));
%plot all spikes
[a,b]=size(nutimes);
U1_index=[NaN,1,1,1,1,1,1];
U4_index=[NaN,NaN,2,2,2,2,2];
U5_index=[NaN,NaN,3,3,3,3,3];

for tn=1:length(filename)
    X_C{tn}=interp1(xvalues{tn},hd_crct{tn}(:,2),time{tn},'spline');
    Y_C{tn}=interp1(xvalues{tn},hd_crct{tn}(:,1),time{tn},'spline');
    
    fig=figure(3);
    fig.Renderer='painters'
    for i=1:4
        
        subplot(4,length(filename),tn+(i-1)*length(filename))
        if filename{tn}(11)=='1' || filename{tn}(11)=='3' || filename{tn}(11)=='4' || filename{tn}(11)=='5'
            load(strcat('landmark_coordinates_',filename{tn}(1:8),'_t3.mat'));
            plot(tank_circle(1,:),(1200-tank_circle(2,:)),'Color',[0.5,0.5,0.5],'LineWidth',0.5)
            hold on
            %plot(hc(:,1),(1200-hc(:,2)))
            plot(stp(:,1),(1200-stp(:,2)),'Color',[0.5,0.5,0.5],'LineWidth',0.5)
            plot(ltp(:,1),(1200-ltp(:,2)),'Color',[0.5,0.5,0.5],'LineWidth',0.5)
            plot(sqp(:,1),(1200-sqp(:,2)),'Color',[0.5,0.5,0.5],'LineWidth',0.5)
            plot(sc(1,:),(1200-sc(2,:)),'Color',[0.5,0.5,0.5],'LineWidth',0.5)
            plot(lc(1,:),(1200-lc(2,:)),'Color',[0.5,0.5,0.5],'LineWidth',0.5)
            plot(filter(:,1),(1200-filter(:,2)),'Color',[0.5,0.5,0.5],'LineWidth',0.5)
            
        else
            load(strcat('landmark_coordinates_',filename{tn},'.mat'));
            plot(tank_circle(1,:),(1200-tank_circle(2,:)),'Color',[0.5,0.5,0.5],'LineWidth',0.5)
            hold on
            %plot(hc(:,1),(1200-hc(:,2)))
            plot(stp(:,1),(1200-stp(:,2)),'Color',[0.5,0.5,0.5],'LineWidth',0.5)
            plot(ltp(:,1),(1200-ltp(:,2)),'Color',[0.5,0.5,0.5],'LineWidth',0.5)
            plot(sqp(:,1),(1200-sqp(:,2)),'Color',[0.5,0.5,0.5],'LineWidth',0.5)
            plot(sc(1,:),(1200-sc(2,:)),'Color',[0.5,0.5,0.5],'LineWidth',0.5)
            plot(lc(1,:),(1200-lc(2,:)),'Color',[0.5,0.5,0.5],'LineWidth',0.5)
            plot(filter(:,1),(1200-filter(:,2)),'Color',[0.5,0.5,0.5],'LineWidth',0.5)
            
        end
        axis equal
        axis off
        
        

    end
    if ~isnan(U1_index(tn))
        subplot(4,length(filename),tn+length(filename))
        x_spike_p1=Place_int_unit{tn,U1_index(tn),1}(nutimes{U1_index(tn),tn}>=0 & nutimes{U1_index(tn),tn}<(time{tn}(end)/2));
        y_spike_p1=Place_int_unit{tn,U1_index(tn),2}(nutimes{U1_index(tn),tn}>=0 & nutimes{U1_index(tn),tn}<(time{tn}(end)/2));
        plot(x_spike_p1,y_spike_p1,'r.','MarkerSize',0.5)
        hold on
        x_spike_p2=Place_int_unit{tn,U1_index(tn),1}(nutimes{U1_index(tn),tn}>=(time{tn}(end)/2) & nutimes{U1_index(tn),tn}<(time{tn}(end)));
        y_spike_p2=Place_int_unit{tn,U1_index(tn),2}(nutimes{U1_index(tn),tn}>=(time{tn}(end)/2) & nutimes{U1_index(tn),tn}<(time{tn}(end)));
        plot(x_spike_p2,y_spike_p2,'k.','MarkerSize',0.5)
        hold on
    else
           axis equal
    axis off
    
    end
    axis equal
    axis off
    
    
    if ~isnan(U4_index(tn))
        subplot(4,length(filename),tn+2*length(filename))
        x_spike_p1=Place_int_unit{tn,U4_index(tn),1}(nutimes{U4_index(tn),tn}>=0 & nutimes{U4_index(tn),tn}<(time{tn}(end)/2));
        y_spike_p1=Place_int_unit{tn,U4_index(tn),2}(nutimes{U4_index(tn),tn}>=0 & nutimes{U4_index(tn),tn}<(time{tn}(end)/2));
        plot(x_spike_p1,y_spike_p1,'r.','MarkerSize',0.5)
        hold on
        x_spike_p2=Place_int_unit{tn,U4_index(tn),1}(nutimes{U4_index(tn),tn}>=(time{tn}(end)/2) & nutimes{U4_index(tn),tn}<(time{tn}(end)));
        y_spike_p2=Place_int_unit{tn,U4_index(tn),2}(nutimes{U4_index(tn),tn}>=(time{tn}(end)/2) & nutimes{U4_index(tn),tn}<(time{tn}(end)));
        plot(x_spike_p2,y_spike_p2,'k.','MarkerSize',0.5)
        hold on
    else
           axis equal
    axis off
    
    end
    
    axis equal
    axis off
    
    
    if ~isnan(U5_index(tn))
        subplot(4,length(filename),tn+3*length(filename))
        x_spike_p1=Place_int_unit{tn,U5_index(tn),1}(nutimes{U5_index(tn),tn}>=0 & nutimes{U5_index(tn),tn}<(time{tn}(end)/2));
        y_spike_p1=Place_int_unit{tn,U5_index(tn),2}(nutimes{U5_index(tn),tn}>=0 & nutimes{U5_index(tn),tn}<(time{tn}(end)/2));
        plot(x_spike_p1,y_spike_p1,'r.','MarkerSize',0.5)
        hold on
        x_spike_p2=Place_int_unit{tn,U5_index(tn),1}(nutimes{U5_index(tn),tn}>=(time{tn}(end)/2) & nutimes{U5_index(tn),tn}<(time{tn}(end)));
        y_spike_p2=Place_int_unit{tn,U5_index(tn),2}(nutimes{U5_index(tn),tn}>=(time{tn}(end)/2) & nutimes{U5_index(tn),tn}<(time{tn}(end)));
        plot(x_spike_p2,y_spike_p2,'k.','MarkerSize',0.5)
        hold on
    else
           axis equal
    axis off
    
    end
    
    axis equal
    axis off
    
    
    %%plot trajectory bellow
    
    figure(3)
    subplot(4,length(filename),tn)
    
    if filename{tn}(11)=='1' || filename{tn}(11)=='3' || filename{tn}(11)=='4' || filename{tn}(11)=='5'
        load(strcat('landmark_coordinates_',filename{tn}(1:8),'_t3.mat'));
        plot(tank_circle(1,:),(1200-tank_circle(2,:)),'Color',[0.5,0.5,0.5],'LineWidth',0.5)
        hold on
        %plot(hc(:,1),(1200-hc(:,2)))
        plot(stp(:,1),(1200-stp(:,2)),'Color',[0.5,0.5,0.5],'LineWidth',0.5)
        plot(ltp(:,1),(1200-ltp(:,2)),'Color',[0.5,0.5,0.5],'LineWidth',0.5)
        plot(sqp(:,1),(1200-sqp(:,2)),'Color',[0.5,0.5,0.5],'LineWidth',0.5)
        plot(sc(1,:),(1200-sc(2,:)),'Color',[0.5,0.5,0.5],'LineWidth',0.5)
        plot(lc(1,:),(1200-lc(2,:)),'Color',[0.5,0.5,0.5],'LineWidth',0.5)
        plot(filter(:,1),(1200-filter(:,2)),'Color',[0.5,0.5,0.5],'LineWidth',0.5)
        
    else
        load(strcat('landmark_coordinates_',filename{tn},'.mat'));
        plot(tank_circle(1,:),(1200-tank_circle(2,:)),'Color',[0.5,0.5,0.5],'LineWidth',0.5)
        hold on
        %plot(hc(:,1),(1200-hc(:,2)))
        plot(stp(:,1),(1200-stp(:,2)),'Color',[0.5,0.5,0.5],'LineWidth',0.5)
        plot(ltp(:,1),(1200-ltp(:,2)),'Color',[0.5,0.5,0.5],'LineWidth',0.5)
        plot(sqp(:,1),(1200-sqp(:,2)),'Color',[0.5,0.5,0.5],'LineWidth',0.5)
        plot(sc(1,:),(1200-sc(2,:)),'Color',[0.5,0.5,0.5],'LineWidth',0.5)
        plot(lc(1,:),(1200-lc(2,:)),'Color',[0.5,0.5,0.5],'LineWidth',0.5)
        plot(filter(:,1),(1200-filter(:,2)),'Color',[0.5,0.5,0.5],'LineWidth',0.5)
        
    end
    
    
    x_p1=Y_C{tn}(time{tn}>=0 & time{tn}<(time{tn}(end)/2));
    y_p1=X_C{tn}(time{tn}>=0 & time{tn}<(time{tn}(end)/2));
    
    plot(x_p1,y_p1,'r','Linewidth',0.1)
    hold on
    x_p2=Y_C{tn}(time{tn}>=(time{tn}(end)/2) & time{tn}<(time{tn}(end)));
    y_p2=X_C{tn}(time{tn}>=(time{tn}(end)/2) & time{tn}<(time{tn}(end)));
    plot(x_p2,y_p2,'k','Linewidth',0.1)
    
    axis equal
    axis off
    
    
end
   saveas(fig,strcat(exp_date,'AllspikesTrajScatter.pdf'))


for tn=1:length(filename)
    X_C{tn}=interp1(xvalues{tn},hd_crct{tn}(:,2),time{tn},'spline');
    Y_C{tn}=interp1(xvalues{tn},hd_crct{tn}(:,1),time{tn},'spline');

    if filename{tn}(11)=='1' || filename{tn}(11)=='3' || filename{tn}(11)=='4' || filename{tn}(11)=='5'
        load(strcat('landmark_coordinates_',filename{tn}(1:8),'_t3.mat'));
    
    %plot the landmarks
    figure(1)
    plot(tank_circle(1,:),(1200-tank_circle(2,:)),'b','LineWidth',0.5)
    hold on
    tank_circle2(1,:)=tank_circle(1,:)-90*cos(ang);
    tank_circle2(2,:)=tank_circle(2,:)-90*sin(ang);

    %plot(hc(:,1),(1200-hc(:,2)))
    plot(stp(:,1),(1200-stp(:,2)),'b','LineWidth',0.5)
    plot(ltp(:,1),(1200-ltp(:,2)),'b','LineWidth',0.5)
    plot(sqp(:,1),(1200-sqp(:,2)),'b','LineWidth',0.5)
    plot(sc(1,:),(1200-sc(2,:)),'b','LineWidth',0.5)
    plot(lc(1,:),(1200-lc(2,:)),'b','LineWidth',0.5)
    plot(filter(:,1),(1200-filter(:,2)),'b','LineWidth',0.5)
       
    %10 cm area around small triangle
    stp2(1,1)=stp(1,1)-(75*sind(30));
    stp2(1,2)=stp(1,2)+(75*cosd(30));
    
    stp2(3,1)=stp(3,1)+75;
    stp2(3,2)=stp(3,2);
    
    stp2(2,1)=stp(2,1)-(75*sind(30));
    stp2(2,2)=stp(2,2)-(75*cosd(30));
    
    stp2(4,:)=stp2(1,:);
    
    plot(stp2(:,1),(1200-stp2(:,2)),'r')
    stp2grid=polygrid(stp2(:,1),1200-stp2(:,2),(1/80));
    
    %same around the large triangle
    
    ltp2(1,1)=ltp(1,1)-(75*sind(30));
    ltp2(1,2)=ltp(1,2)+(75*cosd(30));
    
    ltp2(3,1)=ltp(3,1)+75;
    ltp2(3,2)=ltp(3,2);
    
    ltp2(2,1)=ltp(2,1)-(75*sind(30));
    ltp2(2,2)=ltp(2,2)-(75*cosd(30));
    
    ltp2(4,:)=ltp2(1,:);
    
    plot(ltp2(:,1),(1200-ltp2(:,2)),'r')
    ltp2grid=polygrid(ltp2(:,1),1200-ltp2(:,2),(1/80));
    % plot(ltp2grid(:,1),ltp2grid(:,2),'*')
    
    %10 cm area around the small square
    sqp2(1,:)=sqp(1,:)-75;
    sqp2(2,2)=sqp(2,2)-75;
    sqp2(2,1)=sqp(2,1)+75;
    sqp2(3,1)=sqp2(2,1);
    sqp2(3,2)=sqp(3,2)+75;
    sqp2(4,1)=sqp2(1,1);
    sqp2(4,2)=sqp2(3,2);
    sqp2(5,:)=sqp2(1,:);
    plot(sqp2(:,1),(1200-sqp2(:,2)),'r')
    sqp2grid=polygrid(sqp2(:,1),1200-sqp2(:,2),(1/80));
    %plot(sqp2grid(:,1),sqp2grid(:,2),'*')
    
    %10 cm radius around the circles
    sc2(1,:)=sc(1,:)+75*cos(ang);
    sc2(2,:)=sc(2,:)+75*sin(ang);
    plot(sc2(1,:),1200-sc2(2,:),'r')
    sc2grid=polygrid(sc2(1,:),1200-sc2(2,:),(1/80));
    %plot(sc2grid(:,1),sc2grid(:,2),'*')
    
    lc2(1,:)=lc(1,:)+75*cos(ang);
    lc2(2,:)=lc(2,:)+75*sin(ang);
    plot(lc2(1,:),1200-lc2(2,:),'r')
    lc2grid=polygrid(lc2(1,:),1200-lc2(2,:),(1/80));
    %plot(lc2grid(:,1),lc2grid(:,2),'*')
    
    %10 cm radius around the filter staying within the tank
    filter2(1,1)=filter(1,1)-75;
    filter2(1,2)=filter(1,2)-75;
    filter2(2,1)=filter(2,1)+75;
    filter2(2,2)=filter(2,2)-75;
    filter2(3,1)=filter(3,1)+75;
    filter2(3,2)=filter(3,2);
    filter2(4,1)=filter(4,1)-75;
    filter2(4,2)=filter(4,2);
    filter2(5,:)=filter2(1,:);
    
    plot(filter2(:,1),(1200-filter2(:,2)),'r')
    filter2grid=polygrid(filter2(:,1),1200-filter2(:,2),(1/80));
   else
        load(strcat('landmark_coordinates_',filename{tn},'.mat'));
        %plot the landmarks
        figure(1)
        plot(tank_circle(1,:),(1200-tank_circle(2,:)),'k')
        hold on
        %plot(hc(:,1),(1200-hc(:,2)))
        plot(stp(:,1),(1200-stp(:,2)),'k')
        plot(ltp(:,1),(1200-ltp(:,2)),'k')
        plot(sqp(:,1),(1200-sqp(:,2)),'k')
        plot(sc(1,:),(1200-sc(2,:)),'k')
        plot(lc(1,:),(1200-lc(2,:)),'k')
        plot(filter(:,1),(1200-filter(:,2)),'k')

        %10 cm area around small triangle
        stp2(1,1)=stp(1,1)-(75*sind(30));
        stp2(1,2)=stp(1,2)-(75*cosd(30));

        stp2(3,1)=stp(3,1);
        stp2(3,2)=stp(3,2)+75;

        stp2(2,1)=stp(2,1)+(75*sind(30));
        stp2(2,2)=stp(2,2)-(75*cosd(30));

        stp2(4,:)=stp2(1,:);

        plot(stp2(:,1),(1200-stp2(:,2)),'r')
        stp2grid=polygrid(stp2(:,1),1200-stp2(:,2),(1/80));

        %same around the large triangle

        ltp2(1,1)=ltp(1,1)-(75*sind(30));
        ltp2(1,2)=ltp(1,2)+(75*cosd(30));

        ltp2(3,1)=ltp(3,1)+75;
        ltp2(3,2)=ltp(3,2);

        ltp2(2,1)=ltp(2,1)-(75*sind(30));
        ltp2(2,2)=ltp(2,2)-(75*cosd(30));

        ltp2(4,:)=ltp2(1,:);

        plot(ltp2(:,1),(1200-ltp2(:,2)),'r')
        ltp2grid=polygrid(ltp2(:,1),1200-ltp2(:,2),(1/80));
        % plot(ltp2grid(:,1),ltp2grid(:,2),'*')

        %10 cm area around the small square
        sqp2(1,:)=sqp(1,:)-75;
        sqp2(2,2)=sqp(2,2)-75;
        sqp2(2,1)=sqp(2,1)+75;
        sqp2(3,1)=sqp2(2,1);
        sqp2(3,2)=sqp(3,2)+75;
        sqp2(4,1)=sqp2(1,1);
        sqp2(4,2)=sqp2(3,2);
        sqp2(5,:)=sqp2(1,:);
        plot(sqp2(:,1),(1200-sqp2(:,2)),'r')
        sqp2grid=polygrid(sqp2(:,1),1200-sqp2(:,2),(1/80));
        %plot(sqp2grid(:,1),sqp2grid(:,2),'*')

        %10 cm radius around the circles
        sc2(1,:)=sc(1,:)+75*cos(ang);
        sc2(2,:)=sc(2,:)+75*sin(ang);
        plot(sc2(1,:),1200-sc2(2,:),'r')
        sc2grid=polygrid(sc2(1,:),1200-sc2(2,:),(1/80));
        %plot(sc2grid(:,1),sc2grid(:,2),'*')

        lc2(1,:)=lc(1,:)+75*cos(ang);
        lc2(2,:)=lc(2,:)+75*sin(ang);
        plot(lc2(1,:),1200-lc2(2,:),'r')
        lc2grid=polygrid(lc2(1,:),1200-lc2(2,:),(1/80));
        %plot(lc2grid(:,1),lc2grid(:,2),'*')

        %10 cm radius around the filter staying within the tank
        filter2(1,1)=filter(1,1)-75;
        filter2(1,2)=filter(1,2)-75;
        filter2(2,1)=filter(2,1)+75;
        filter2(2,2)=filter(2,2)-75;
        filter2(3,1)=filter(3,1)+75;
        filter2(3,2)=filter(3,2);
        filter2(4,1)=filter(4,1)-75;
        filter2(4,2)=filter(4,2);
        filter2(5,:)=filter2(1,:);

        plot(filter2(:,1),(1200-filter2(:,2)),'r')
        filter2grid=polygrid(filter2(:,1),1200-filter2(:,2),(1/80));
   
    end

    %plot(filter2grid(:,1),filter2grid(:,2),'*')
    %for each unit, divide the trial int
%     plot(Place_int_unit{tn,1,1},Place_int_unit{tn,1,2},'ko')
%     plot(Y_C{tn},X_C{tn},'.')
    
    %divide the trial into 10 min periods and calculate the number of
    %visits to the area around each land mark
    %different landmarks are present in different trials.
    timebins=time{tn}(1):(time{tn}(end)-time{tn}(1))/10:time{tn}(end);
    %calculate time spent in the first and second half of the trial
%     for k=1:length(timebins)-1
%         x_temp=Y_C{tn}(time{tn}>=timebins(k) & time{tn}<timebins(k+1));
%         y_temp=X_C{tn}(time{tn}>=timebins(k) & time{tn}<timebins(k+1));
%         in=inpolygon(x_temp,y_temp,hc2(:,1),1200-hc2(:,2));
%         plot(x_temp(in==1),y_temp(in==1),'k.');
%         %calculate the time spent in the roi
%         time_in_hc(tn,k)=sum(inpolygon(x_temp,y_temp,hc2(:,1),1200-hc2(:,2)))/fs;%time spent within the roi during this part of the trial
%         [a,b]=size(nutimes);
%         tn
%         if time_in_hc(tn,k)>0
%             for j=1:a
%                 if length(nutimes{j,tn})>=1
%                     x_spike_temp=Place_int_unit{tn,j,1}(nutimes{j,tn}>=timebins(k) & nutimes{j,tn}<timebins(k+1));
%                     y_spike_temp=Place_int_unit{tn,j,2}(nutimes{j,tn}>=timebins(k) & nutimes{j,tn}<timebins(k+1));
%                     in=inpolygon(x_spike_temp,y_spike_temp,hc2(:,1),1200-hc2(:,2));
%                     plot(x_spike_temp(in==1),y_spike_temp(in==1),'r*')
%                     spkcnt_in_hc(tn,k,j)=sum(inpolygon(x_spike_temp,y_spike_temp,hc2(:,1),1200-hc2(:,2)))
%                     fr_in_hc(tn,k,j)=spkcnt_in_hc(tn,k,j)/time_in_hc(tn,k);
%                 else
%                     spkcnt_in_hc(tn,k,j)=NaN;
%                     fr_in_hc(tn,k,j)=NaN;
% 
%                 end
%             end
%         else
%             
%             fr_in_hc(tn,k,:)=NaN;
%         end
%         
%     end
%     figure(2)
%     title('home')
%     for j=1:a
%         [m,ind]=max((fr_in_hc(tn,:,j)));
%         max_ind{tn,j}=ind;
%         if  tn==2 
%             plot([0,(fr_in_hc(tn,2,j)-fr_in_hc(tn,1,j))/fr_in_hc(tn,1,j)],'r*-')
%             hold on
%         elseif tn==1 || tn==3 || tn==4
%             plot([0,(fr_in_hc(tn,2,j)-fr_in_hc(tn,1,j))/fr_in_hc(tn,1,j)],'b*-')
%         else
%             plot([0,(fr_in_hc(tn,2,j)-fr_in_hc(tn,1,j))/fr_in_hc(tn,1,j)],'m*-')
%         end
%         %calculate the slope which is after minus before
%         fr_diff_hc(tn,j)=fr_in_hc(tn,2,j)-fr_in_hc(tn,1,j);
%     end
%     
    %% do the same thing for the large triangle
    figure(1)
    for k=1:length(timebins)-1
        x_temp=Y_C{tn}(time{tn}>=timebins(k) & time{tn}<timebins(k+1));
        y_temp=X_C{tn}(time{tn}>=timebins(k) & time{tn}<timebins(k+1));
        in=inpolygon(x_temp,y_temp,ltp2(:,1),1200-ltp2(:,2));
        plot(x_temp(in==1),y_temp(in==1),'k.');
        %calculate the time spent in the roi
        time_in_ltp(tn,k)=sum(inpolygon(x_temp,y_temp,ltp2(:,1),1200-ltp2(:,2)))/fs;%time spent within the roi during this part of the trial
        [a,b]=size(nutimes);
        tn
        if time_in_ltp(tn,k)>0
            for j=1:a
                if length(nutimes{j,tn})>=1
                    x_spike_temp=Place_int_unit{tn,j,1}(nutimes{j,tn}>=timebins(k) & nutimes{j,tn}<timebins(k+1));
                    y_spike_temp=Place_int_unit{tn,j,2}(nutimes{j,tn}>=timebins(k) & nutimes{j,tn}<timebins(k+1));
                    if ~isempty(x_spike_temp)
                        in=inpolygon(x_spike_temp,y_spike_temp,ltp2(:,1),1200-ltp2(:,2));
                        plot(x_spike_temp(in==1),y_spike_temp(in==1),'r*')
                        spkcnt_in_ltp(tn,k,j)=sum(inpolygon(x_spike_temp,y_spike_temp,ltp2(:,1),1200-ltp2(:,2)))
                        fr_in_ltp(tn,k,j)=spkcnt_in_ltp(tn,k,j)/time_in_ltp(tn,k);
                    else
                        spkcnt_in_ltp(tn,k,j)=0;
                        fr_in_ltp(tn,k,j)=spkcnt_in_ltp(tn,k,j)/time_in_ltp(tn,k);

                    end
                else
                    spkcnt_in_ltp(tn,k,j)=0;
                    fr_in_ltp(tn,k,j)=spkcnt_in_ltp(tn,k,j)/time_in_ltp(tn,k);
                end
            end
        else
            fr_in_ltp(tn,k,1:a)=NaN;
        end
        
    end
    figure(2)
    for j=1:a
        if sum(~isnan(fr_in_ltp(tn,:,j)))>0
            [m,ind]=max((fr_in_ltp(tn,:,j)));
            max_ind{tn,j}=ind;
            if tn==1 %not present never seen
                %plot([0,(fr_in_ltp(tn,2,j)-fr_in_ltp(tn,1,j))/fr_in_ltp(tn,1,j)],'k*-')
                hold on
            elseif tn==2
                plot([0,(fr_in_ltp(tn,2,j)-fr_in_ltp(tn,1,j))/fr_in_ltp(tn,1,j)],'r*-')
            elseif tn==3 || tn==4
                plot([0,(fr_in_ltp(tn,2,j)-fr_in_ltp(tn,1,j))/fr_in_ltp(tn,1,j)],'b*-')
            elseif tn==5
                plot([0,(fr_in_ltp(tn,2,j)-fr_in_ltp(tn,1,j))/fr_in_ltp(tn,1,j)],'m*-')
            end
        end
        %calculate the slope which is after minus before
    end
    
    %% do the same for small triangle
figure(1)
    for k=1:length(timebins)-1
        x_temp=Y_C{tn}(time{tn}>=timebins(k) & time{tn}<timebins(k+1));
        y_temp=X_C{tn}(time{tn}>=timebins(k) & time{tn}<timebins(k+1));
        in=inpolygon(x_temp,y_temp,stp2(:,1),1200-stp2(:,2));
        plot(x_temp(in==1),y_temp(in==1),'k.');
        %calculate the time spent in the roi
        time_in_stp(tn,k)=sum(inpolygon(x_temp,y_temp,stp2(:,1),1200-stp2(:,2)))/fs;%time spent within the roi during this part of the trial
        [a,b]=size(nutimes);
        tn
        if time_in_stp(tn,k)>0
            for j=1:a
                if length(nutimes{j,tn})>=1
                    
                    x_spike_temp=Place_int_unit{tn,j,1}(nutimes{j,tn}>=timebins(k) & nutimes{j,tn}<timebins(k+1));
                    y_spike_temp=Place_int_unit{tn,j,2}(nutimes{j,tn}>=timebins(k) & nutimes{j,tn}<timebins(k+1));
                    if ~isempty(x_spike_temp)
                        in=inpolygon(x_spike_temp,y_spike_temp,stp2(:,1),1200-stp2(:,2));
                        plot(x_spike_temp(in==1),y_spike_temp(in==1),'r*')
                        spkcnt_in_stp(tn,k,j)=sum(inpolygon(x_spike_temp,y_spike_temp,stp2(:,1),1200-stp2(:,2)))
                        fr_in_stp(tn,k,j)=spkcnt_in_stp(tn,k,j)/time_in_stp(tn,k);
                    else
                        spkcnt_in_stp(tn,k,j)=0;
                        fr_in_stp(tn,k,j)=spkcnt_in_stp(tn,k,j)/time_in_stp(tn,k);
                    end
                else
                    spkcnt_in_stp(tn,k,j)=0;
                    fr_in_stp(tn,k,j)=spkcnt_in_stp(tn,k,j)/time_in_stp(tn,k);
                end
            end
        else
            fr_in_stp(tn,k,1:a)=NaN;
        end
        
    end
    figure(2)
    for j=1:a
        if sum(~isnan(fr_in_stp(tn,:,j)))>0
            [m,ind]=max((fr_in_stp(tn,:,j)));
            max_ind{tn,j}=ind;
            if tn==1 %not present never seen
                %plot([0,(fr_in_stp(tn,2,j)-fr_in_stp(tn,1,j))/fr_in_stp(tn,1,j)],'k*-')
                hold on
            elseif tn==2
                plot([0,(fr_in_stp(tn,2,j)-fr_in_stp(tn,1,j))/fr_in_stp(tn,1,j)],'r*-')
            elseif tn==3 || tn==4
                plot([0,(fr_in_stp(tn,2,j)-fr_in_stp(tn,1,j))/fr_in_stp(tn,1,j)],'b*-')
            elseif tn==5
                plot([0,(fr_in_stp(tn,2,j)-fr_in_stp(tn,1,j))/fr_in_stp(tn,1,j)],'m*-')
            end
        end
        %calculate the slope which is after minus before
    end
        %% do the same for small square
figure(1)
    for k=1:length(timebins)-1
        x_temp=Y_C{tn}(time{tn}>=timebins(k) & time{tn}<timebins(k+1));
        y_temp=X_C{tn}(time{tn}>=timebins(k) & time{tn}<timebins(k+1));
        in=inpolygon(x_temp,y_temp,sqp2(:,1),1200-sqp2(:,2));
        plot(x_temp(in==1),y_temp(in==1),'k.');
        %calculate the time spent in the roi
        time_in_sqp(tn,k)=sum(inpolygon(x_temp,y_temp,sqp2(:,1),1200-sqp2(:,2)))/fs;%time spent within the roi during this part of the trial
        [a,b]=size(nutimes);
        tn
        if time_in_sqp(tn,k)>0
            for j=1:a
                if length(nutimes{j,tn})>=1
                    x_spike_temp=Place_int_unit{tn,j,1}(nutimes{j,tn}>=timebins(k) & nutimes{j,tn}<timebins(k+1));
                    y_spike_temp=Place_int_unit{tn,j,2}(nutimes{j,tn}>=timebins(k) & nutimes{j,tn}<timebins(k+1));
                    in=inpolygon(x_spike_temp,y_spike_temp,sqp2(:,1),1200-sqp2(:,2));
                    plot(x_spike_temp(in==1),y_spike_temp(in==1),'r*')
                    spkcnt_in_sqp(tn,k,j)=sum(inpolygon(x_spike_temp,y_spike_temp,sqp2(:,1),1200-sqp2(:,2)))
                    fr_in_sqp(tn,k,j)=spkcnt_in_sqp(tn,k,j)/time_in_sqp(tn,k);
                else
                    spkcnt_in_sqp(tn,k,j)=0;
                    fr_in_sqp(tn,k,j)=spkcnt_in_sqp(tn,k,j)/time_in_sqp(tn,k);
                end
            end
        else
            fr_in_sqp(tn,k,1:a)=NaN;
        end
        
    end
    figure(2)
    for j=1:a
        if sum(~isnan(fr_in_sqp(tn,:,j)))>0
            [m,ind]=max((fr_in_sqp(tn,:,j)));
            max_ind{tn,j}=ind;
            if tn==1 %not present never seen
                %plot([0,(fr_in_sqp(tn,2,j)-fr_in_sqp(tn,1,j))/fr_in_sqp(tn,1,j)],'k*-')
                hold on
            elseif tn==2
                plot([0,(fr_in_sqp(tn,2,j)-fr_in_sqp(tn,1,j))/fr_in_sqp(tn,1,j)],'r*-')
            elseif tn==3 || tn==4
                plot([0,(fr_in_sqp(tn,2,j)-fr_in_sqp(tn,1,j))/fr_in_sqp(tn,1,j)],'b*-')
            elseif tn==5
                plot([0,(fr_in_sqp(tn,2,j)-fr_in_sqp(tn,1,j))/fr_in_sqp(tn,1,j)],'m*-')
            end
        end
        %calculate the slope which is after minus before
    end

%%
    %% do the same for small circle
figure(1)
    for k=1:length(timebins)-1
        x_temp=Y_C{tn}(time{tn}>=timebins(k) & time{tn}<timebins(k+1));
        y_temp=X_C{tn}(time{tn}>=timebins(k) & time{tn}<timebins(k+1));
        in=inpolygon(x_temp,y_temp,sc2(1,:),1200-sc2(2,:));
        plot(x_temp(in==1),y_temp(in==1),'k.');
        %calculate the time spent in the roi
        time_in_sc(tn,k)=sum(inpolygon(x_temp,y_temp,sc2(1,:),1200-sc2(2,:)))/fs;%time spent within the roi during this part of the trial
        [a,b]=size(nutimes);
        tn
        if time_in_sc(tn,k)>0
            for j=1:a
                if length(nutimes{j,tn})>=1
                    x_spike_temp=Place_int_unit{tn,j,1}(nutimes{j,tn}>=timebins(k) & nutimes{j,tn}<timebins(k+1));
                    y_spike_temp=Place_int_unit{tn,j,2}(nutimes{j,tn}>=timebins(k) & nutimes{j,tn}<timebins(k+1));
                    in=inpolygon(x_spike_temp,y_spike_temp,sc2(1,:),1200-sc2(2,:));
                    plot(x_spike_temp(in==1),y_spike_temp(in==1),'r*')
                    spkcnt_in_sc(tn,k,j)=sum(inpolygon(x_spike_temp,y_spike_temp,sc2(1,:),1200-sc2(2,:)))
                    fr_in_sc(tn,k,j)=spkcnt_in_sc(tn,k,j)/time_in_sc(tn,k);
                else
                    spkcnt_in_sc(tn,k,j)=0;
                    fr_in_sc(tn,k,j)=spkcnt_in_sc(tn,k,j)/time_in_sc(tn,k);
                end
            end
        else
            fr_in_sc(tn,k,1:a)=NaN;
        end
        
    end
    figure(2)
    for j=1:a
        if sum(~isnan(fr_in_sc(tn,:,j)))>0
            [m,ind]=max((fr_in_sc(tn,:,j)));
            max_ind{tn,j}=ind;
            if tn==1 %not present never seen
                %plot([0,(fr_in_sc(tn,2,j)-fr_in_sc(tn,1,j))/fr_in_sc(tn,1,j)],'k*-')
                hold on
            elseif tn==2
                plot([0,(fr_in_sc(tn,2,j)-fr_in_sc(tn,1,j))/fr_in_sc(tn,1,j)],'r*-')
            elseif tn==3 || tn==4
                plot([0,(fr_in_sc(tn,2,j)-fr_in_sc(tn,1,j))/fr_in_sc(tn,1,j)],'b*-')
            elseif tn==5
                plot([0,(fr_in_sc(tn,2,j)-fr_in_sc(tn,1,j))/fr_in_sc(tn,1,j)],'m*-')
            end
        end
        %calculate the slope which is after minus before
    end
%% do the same for the large circle
figure(1)
    for k=1:length(timebins)-1
        x_temp=Y_C{tn}(time{tn}>=timebins(k) & time{tn}<timebins(k+1));
        y_temp=X_C{tn}(time{tn}>=timebins(k) & time{tn}<timebins(k+1));
        in=inpolygon(x_temp,y_temp,lc2(1,:),1200-lc2(2,:));
        plot(x_temp(in==1),y_temp(in==1),'k.');
        %calculate the time spent in the roi
        time_in_lc(tn,k)=sum(inpolygon(x_temp,y_temp,lc2(1,:),1200-lc2(2,:)))/fs;%time spent within the roi during this part of the trial
        [a,b]=size(nutimes);
        tn
        if time_in_lc(tn,k)>0
            for j=1:a
                if length(nutimes{j,tn})>=1
                    x_spike_temp=Place_int_unit{tn,j,1}(nutimes{j,tn}>=timebins(k) & nutimes{j,tn}<timebins(k+1));
                    y_spike_temp=Place_int_unit{tn,j,2}(nutimes{j,tn}>=timebins(k) & nutimes{j,tn}<timebins(k+1));
                    in=inpolygon(x_spike_temp,y_spike_temp,lc2(1,:),1200-lc2(2,:));
                    plot(x_spike_temp(in==1),y_spike_temp(in==1),'r*')
                    spkcnt_in_lc(tn,k,j)=sum(inpolygon(x_spike_temp,y_spike_temp,lc2(1,:),1200-lc2(2,:)))
                    fr_in_lc(tn,k,j)=spkcnt_in_lc(tn,k,j)/time_in_lc(tn,k);
                else
                    spkcnt_in_lc(tn,k,j)=0;
                    fr_in_lc(tn,k,j)=spkcnt_in_lc(tn,k,j)/time_in_lc(tn,k);
                end
            end
        else
            fr_in_lc(tn,k,1:a)=NaN;
        end
        
    end
    figure(2)
    for j=1:a
        if sum(~isnan(fr_in_lc(tn,:,j)))>0
            [m,ind]=max((fr_in_lc(tn,:,j)));
            max_ind{tn,j}=ind;
            if tn==1 %not present never seen
                %plot([0,(fr_in_lc(tn,2,j)-fr_in_lc(tn,1,j))/fr_in_lc(tn,1,j)],'k*-')
                hold on
            elseif tn==2
                plot([0,(fr_in_lc(tn,2,j)-fr_in_lc(tn,1,j))/fr_in_lc(tn,1,j)],'r*-')
            elseif tn==3 || tn==4
                plot([0,(fr_in_lc(tn,2,j)-fr_in_lc(tn,1,j))/fr_in_lc(tn,1,j)],'b*-')
            elseif tn==5
                plot([0,(fr_in_lc(tn,2,j)-fr_in_lc(tn,1,j))/fr_in_lc(tn,1,j)],'m*-')
            end
        end
        %calculate the slope which is after minus before
    end
%% do the same for filter
figure(1)
    for k=1:length(timebins)-1
        x_temp=Y_C{tn}(time{tn}>=timebins(k) & time{tn}<timebins(k+1));
        y_temp=X_C{tn}(time{tn}>=timebins(k) & time{tn}<timebins(k+1));
        in=inpolygon(x_temp,y_temp,filter2(:,1),1200-filter2(:,2));
        plot(x_temp(in==1),y_temp(in==1),'k.');
        %calculate the time spent in the roi
        time_in_filter(tn,k)=sum(inpolygon(x_temp,y_temp,filter2(:,1),1200-filter2(:,2)))/fs;%time spent within the roi during this part of the trial
        [a,b]=size(nutimes);
        tn
        if time_in_filter(tn,k)>0
            for j=1:a
                if length(nutimes{j,tn})>=1
                    x_spike_temp=Place_int_unit{tn,j,1}(nutimes{j,tn}>=timebins(k) & nutimes{j,tn}<timebins(k+1));
                    y_spike_temp=Place_int_unit{tn,j,2}(nutimes{j,tn}>=timebins(k) & nutimes{j,tn}<timebins(k+1));
                    in=inpolygon(x_spike_temp,y_spike_temp,filter2(:,1),1200-filter2(:,2));
                    plot(x_spike_temp(in==1),y_spike_temp(in==1),'r*')
                    spkcnt_in_filter(tn,k,j)=sum(inpolygon(x_spike_temp,y_spike_temp,filter2(:,1),1200-filter2(:,2)))
                    fr_in_filter(tn,k,j)=spkcnt_in_filter(tn,k,j)/time_in_filter(tn,k);
                else
                    spkcnt_in_filter(tn,k,j)=0;
                    fr_in_filter(tn,k,j)=spkcnt_in_filter(tn,k,j)/time_in_filter(tn,k);
                end
            end
        else
            fr_in_filter(tn,k,1:a)=NaN;
        end
        
    end
    figure(2)
    for j=1:a
        if sum(~isnan(fr_in_filter(tn,:,j)))>0
            [m,ind]=max((fr_in_filter(tn,:,j)));
            max_ind{tn,j}=ind;
            if tn==1 %not present never seen
                plot([0,(fr_in_filter(tn,2,j)-fr_in_filter(tn,1,j))/fr_in_filter(tn,1,j)],'k*-')
                hold on
            elseif tn==2
                plot([0,(fr_in_filter(tn,2,j)-fr_in_filter(tn,1,j))/fr_in_filter(tn,1,j)],'r*-')
            elseif tn==3 || tn==4
                plot([0,(fr_in_filter(tn,2,j)-fr_in_filter(tn,1,j))/fr_in_filter(tn,1,j)],'b*-')
            elseif tn==5
                plot([0,(fr_in_filter(tn,2,j)-fr_in_filter(tn,1,j))/fr_in_filter(tn,1,j)],'m*-')
            end
        end
        %calculate the slope which is after minus before
    end


%% now do this for spikes that are outside of all polygons
figure(1)
for k=1:length(timebins)-1
    x_temp=Y_C{tn}(time{tn}>=timebins(k) & time{tn}<timebins(k+1));
    y_temp=X_C{tn}(time{tn}>=timebins(k) & time{tn}<timebins(k+1));
%     in_home=inpolygon(x_temp,y_temp,hc2(:,1),1200-hc2(:,2));
    out_tank=inpolygon(x_temp,y_temp,tank_circle2(1,:),1200-tank_circle2(2,:));
    
    in_filter=inpolygon(x_temp,y_temp,filter2(:,1),1200-filter2(:,2));
    in_ltp=inpolygon(x_temp,y_temp,ltp2(:,1),1200-ltp2(:,2));
    in_stp=inpolygon(x_temp,y_temp,stp2(:,1),1200-stp2(:,2));
    in_sqp=inpolygon(x_temp,y_temp,sqp2(:,1),1200-sqp2(:,2));
    in_sc=inpolygon(x_temp,y_temp,sc2(1,:),1200-sc2(2,:));
    in_lc=inpolygon(x_temp,y_temp,lc2(1,:),1200-lc2(2,:));

    if filename{tn}(11)=='1' || filename{tn}(11)=='3' %not including t1 anyway
        index_out=find(in_filter==0 & in_ltp==0 & in_stp==0 & in_sqp==0 & in_sc==0 & in_lc==0 & out_tank==1);

    elseif filename{tn}(11)=='4'
        
        index_out=find(in_filter==0 & in_stp==0 & in_sqp==0 & in_sc==0 & in_lc==0& out_tank==1);
        
    elseif filename{tn}(11)== '5' || filename{tn}(11)== '6' || filename{tn}(11)== '7'
        
        index_out=find(in_filter==0 & in_ltp==0 & in_stp==0 & in_sqp==0 & in_lc==0& out_tank==1);
        
    else
        index_out=find(in_filter==0 & in_ltp==0 & in_stp==0& out_tank==1);
       
    end
    plot(x_temp(index_out),y_temp(index_out),'y.');
    %calculate the time spent in the roi
    time_out_all(tn,k)=length(index_out)/fs;%time spent within the roi during this part of the trial
    [a,b]=size(nutimes);
    tn
    if time_out_all(tn,k)>0
        for j=1:a
            if length(nutimes{j,tn})>=1
                x_spike_temp=Place_int_unit{tn,j,1}(nutimes{j,tn}>=timebins(k) & nutimes{j,tn}<timebins(k+1));
                y_spike_temp=Place_int_unit{tn,j,2}(nutimes{j,tn}>=timebins(k) & nutimes{j,tn}<timebins(k+1));
                %in_spk_home=inpolygon(x_spike_temp,y_spike_temp,hc2(:,1),1200-hc2(:,2));
                out_spk_tank=inpolygon(x_spike_temp,y_spike_temp,tank_circle2(1,:),1200-tank_circle2(2,:));
                in_spk_filter=inpolygon(x_spike_temp,y_spike_temp,filter2(:,1),1200-filter2(:,2));
                in_spk_ltp=inpolygon(x_spike_temp,y_spike_temp,ltp2(:,1),1200-ltp2(:,2));
                in_spk_stp=inpolygon(x_spike_temp,y_spike_temp,stp2(:,1),1200-stp2(:,2));
                in_spk_sqp=inpolygon(x_spike_temp,y_spike_temp,sqp2(:,1),1200-sqp2(:,2));
                in_spk_sc=inpolygon(x_spike_temp,y_spike_temp,sc2(1,:),1200-sc2(2,:));
                in_spk_lc=inpolygon(x_spike_temp,y_spike_temp,lc2(1,:),1200-lc2(2,:));

                if filename{tn}(11)=='1' || filename{tn}(11)=='3'
                    index_spk_out=find(in_spk_filter==0 & in_spk_ltp==0 & in_spk_stp==0 & in_spk_sqp==0 & in_spk_sc==0 & in_spk_lc==0 & out_spk_tank==1);
                    
                elseif filename{tn}(11)=='4'
                    
                    index_spk_out=find(in_spk_filter==0 & in_spk_stp==0 & in_spk_sqp==0 & in_spk_sc==0 & in_spk_lc==0 & out_spk_tank==1);
                    
                elseif filename{tn}(11)== '5' || filename{tn}(11)== '6' || filename{tn}(11)== '7'
                    
                    index_spk_out=find(in_spk_filter==0 & in_spk_ltp==0 & in_spk_stp==0 & in_spk_sqp==0 & in_spk_lc==0 & out_spk_tank==1);
                    
                else
                    index_spk_out=find(in_spk_filter==0 & in_spk_ltp==0 & in_spk_stp==0 & out_spk_tank==1);
                    
                end
                plot(x_spike_temp(index_spk_out),y_spike_temp(index_spk_out),'bo')
                spkcnt_out_all(tn,k,j)=length(index_spk_out);
                fr_out_all(tn,k,j)=spkcnt_out_all(tn,k,j)/time_out_all(tn,k);
            else
                spkcnt_out_all(tn,k,j)=0;
                fr_out_all(tn,k,j)=spkcnt_out_all(tn,k,j)/time_out_all(tn,k);
            end
        end
    else
        fr_out_all(tn,k,1:a)=NaN;
    end
    
end%% now do this for spikes that are inside of all polygons
%% now do this for spikes that are inside of all polygons
for k=1:length(timebins)-1
    x_temp=Y_C{tn}(time{tn}>=timebins(k) & time{tn}<timebins(k+1));
    y_temp=X_C{tn}(time{tn}>=timebins(k) & time{tn}<timebins(k+1));
    %in_home=inpolygon(x_temp,y_temp,hc2(:,1),1200-hc2(:,2));
    out_tank=inpolygon(x_temp,y_temp,tank_circle2(1,:),1200-tank_circle2(2,:));
    in_filter=inpolygon(x_temp,y_temp,filter2(:,1),1200-filter2(:,2));
    in_ltp=inpolygon(x_temp,y_temp,ltp2(:,1),1200-ltp2(:,2));
    in_stp=inpolygon(x_temp,y_temp,stp2(:,1),1200-stp2(:,2));
    in_sqp=inpolygon(x_temp,y_temp,sqp2(:,1),1200-sqp2(:,2));
    in_sc=inpolygon(x_temp,y_temp,sc2(1,:),1200-sc2(2,:));
    in_lc=inpolygon(x_temp,y_temp,lc2(1,:),1200-lc2(2,:));
        
    if filename{tn}(11)=='1' || filename{tn}(11)=='3'
        index_in=find(in_filter==1 | in_ltp==1 | in_stp==1 | in_sqp==1 | in_sc==1 | in_lc==1 | out_tank==0);

    elseif filename{tn}(11)=='4'
        
        index_in=find(in_filter==1 | in_stp==1 | in_sqp==1 | in_sc==1 | in_lc==1| out_tank==0);
        
    elseif filename{tn}(11)== '5' || filename{tn}(11)== '6' || filename{tn}(11)== '7'
        
        index_in=find(in_filter==1 | in_ltp==1 | in_stp==1 | in_sqp==1 | in_lc==1| out_tank==0);
        
    else
        index_in=find(in_filter==1 | in_ltp==1 | in_stp==1| out_tank==0);
       
    end
    %plot(x_temp(index_in),y_temp(index_in),'y.');
    %calculate the time spent in the roi
    time_in_all(tn,k)=length(index_in)/fs;%time spent within the roi during this part of the trial
    [a,b]=size(nutimes);
    tn
    if time_in_all(tn,k)>0
        for j=1:a
            if length(nutimes{j,tn})>=1
                x_spike_temp=Place_int_unit{tn,j,1}(nutimes{j,tn}>=timebins(k) & nutimes{j,tn}<timebins(k+1));
                y_spike_temp=Place_int_unit{tn,j,2}(nutimes{j,tn}>=timebins(k) & nutimes{j,tn}<timebins(k+1));
                %in_spk_home=inpolygon(x_spike_temp,y_spike_temp,hc2(:,1),1200-hc2(:,2));
                 out_spk_tank=inpolygon(x_spike_temp,y_spike_temp,tank_circle2(1,:),1200-tank_circle2(2,:));

                in_spk_filter=inpolygon(x_spike_temp,y_spike_temp,filter2(:,1),1200-filter2(:,2));
                in_spk_ltp=inpolygon(x_spike_temp,y_spike_temp,ltp2(:,1),1200-ltp2(:,2));
                in_spk_stp=inpolygon(x_spike_temp,y_spike_temp,stp2(:,1),1200-stp2(:,2));
                in_spk_sqp=inpolygon(x_spike_temp,y_spike_temp,sqp2(:,1),1200-sqp2(:,2));
                in_spk_sc=inpolygon(x_spike_temp,y_spike_temp,sc2(1,:),1200-sc2(2,:));
                in_spk_lc=inpolygon(x_spike_temp,y_spike_temp,lc2(1,:),1200-lc2(2,:));
                
                if filename{tn}(11)=='1' || filename{tn}(11)=='3'
                    index_spk_in=find(in_spk_filter==1 | in_spk_ltp==1 | in_spk_stp==1 | in_spk_sqp==1 | in_spk_sc==1 | in_spk_lc==1 | out_spk_tank==0);
                    
                elseif filename{tn}(11)=='4'
                    
                    index_spk_in=find(in_spk_filter==1 | in_spk_stp==1 | in_spk_sqp==1 | in_spk_sc==1 | in_spk_lc==1| out_spk_tank==0);
                    
                elseif filename{tn}(11)== '5' || filename{tn}(11)== '6' || filename{tn}(11)== '7'
                    
                    index_spk_in=find(in_spk_filter==1 | in_spk_ltp==1 | in_spk_stp==1 | in_spk_sqp==1 | in_spk_lc==1| out_spk_tank==0);
                    
                else
                    index_spk_in=find(in_spk_filter==1 | in_spk_ltp==1 | in_spk_stp==1| out_spk_tank==0);
                    
                end

                plot(x_spike_temp(index_spk_in),y_spike_temp(index_spk_in),'ro')
                spkcnt_in_all(tn,k,j)=length(index_spk_in);
                fr_in_all(tn,k,j)=spkcnt_in_all(tn,k,j)/time_in_all(tn,k);
            else
                spkcnt_in_all(tn,k,j)=NaN;
                fr_out_all(tn,k,j)=NaN;
            end
        end
    else
        fr_in_all(tn,k,1:a)=NaN;
    end
    
end


%%calculate the firing rate trend over time regardless of location
for k=1:length(timebins)-1
    for j=1:a
        if length(nutimes{j,tn})>=1
            x_spike_temp=Place_int_unit{tn,j,1}(nutimes{j,tn}>=timebins(k) & nutimes{j,tn}<timebins(k+1));
            all_spks(tn,k,j)=length(x_spike_temp);
            fr_all_spks(tn,k,j)=all_spks(tn,k,j)/mean(diff(timebins));
        else
            
            all_spks(tn,k,j)=NaN;
            fr_all_spks(tn,k,j)=NaN;
        end
    end
end

end
%for each cell, construct a vector with the following groups
%%
close all
figure(2)
for tn=1:length(filename)
    scale1=1;
    scale2=1;
    a=0;
    
    if filename{tn}(11)=='1' || filename{tn}(11)=='3' || filename{tn}(11)=='4' || filename{tn}(11)=='5'
        load(strcat('landmark_coordinates_',filename{tn}(1:8),'_t3.mat'));
    else
        
        load(strcat('landmark_coordinates_',filename{tn},'.mat'));
    end
    subplot(1,7,tn)
    hold on
    plot((filter(:,1))/scale1,(1200-filter(:,2))/scale2)
    plot((tank_circle(1,:))/scale1,(1200-tank_circle(2,:))/scale2,'k')
    if  filename{tn}(11)=='1'
        plot((hc(:,1))/scale1,(1200-hc(:,2))/scale2)
        
    elseif filename{tn}(11)=='3'
        plot((stp(:,1))/scale1,(1200-stp(:,2))/scale2);
        plot((ltp(:,1))/scale1,(1200-ltp(:,2))/scale2);
        plot((sqp(:,1))/scale1,(1200-sqp(:,2))/scale2);
        plot((sc(1,:))/scale1,(1200-sc(2,:))/scale2);
        plot((lc(1,:))/scale1,(1200-lc(2,:))/scale2);
    elseif filename{tn}(11)=='4'
        
        plot(stp(:,1)/scale1,(1200-stp(:,2))/scale2);
        plot(sqp(:,1)/scale1,(1200-sqp(:,2))/scale2);
        plot(sc(1,:)/scale1,(1200-sc(2,:))/scale2);
        plot(lc(1,:)/scale1,(1200-lc(2,:))/scale2);
        
    elseif filename{tn}(11)== '5'
        plot(stp(:,1)/scale1,(1200-stp(:,2))/scale2);
        plot(ltp(:,1)/scale1,(1200-ltp(:,2))/scale2);
        plot(sqp(:,1)/scale1,(1200-sqp(:,2))/scale2);
        plot(lc(1,:)/scale1,(1200-lc(2,:))/scale2);
    elseif filename{tn}(11)== '6' || filename{tn}(11)== '7'
        plot(stp(:,1)/scale1,(1200-stp(:,2))/scale2);
        plot(ltp(:,1)/scale1,(1200-ltp(:,2))/scale2);
        plot(sqp(:,1)/scale1,(1200-sqp(:,2))/scale2);
        plot(lc(1,:)/scale1,(1200-lc(2,:))/scale2);
    else
        plot(stp(:,1)/scale1,(1200-stp(:,2))/scale2);
        plot(ltp(:,1)/scale1,(1200-ltp(:,2))/scale2);
       
    end
    axis off
    axis square
    
end
figure(2)
saveas(gcf,'landmarkconfigs.pdf');

%%
group=[];
for i=1:7
    group=[group,[(2*i-1)*ones(1,5),(2*i)*ones(1,5)]];
end
group2=[];
for i=1:7
    group2=[group2,i*ones(1,10)];
end

group_color=[];
for i=1:7
    group_color=[group_color,[ones(1,5),2*ones(1,5)]];
end
for i=1:14
    if mod(i,2)==1
        group_separation(i)=1;
    else
        group_separation(i)=2;
    end
end
        

for i=1:14
    if mod(i,2)==1
        axis_pos(i)=i;

        box_color(i)='r';
    else
        axis_pos(i)=(i-1)+0.5;
        box_color(i)='k';
    end
end

for i=1:7
    
    axis_label{2*i-1}=pad(strcat('T',num2str(i)),5,'left');
    axis_label{2*i}='';
    %axis_label{2*i}=strcat(num2str(i),'_p2');

end

% do the analysis around small square which is stable and changes in the
% last trial
% 6. small square


fr_sqp_U1=[];
for i=1:7
    if ~isnan(U1_index(i)) && (sum(~isnan(fr_in_sqp(i,1:5,U1_index(i))))>0 || sum(~isnan(fr_in_sqp(i,6:10,U1_index(i))))>0) %the fish should have traversed at least once in the two halves
            fr_sqp_U1=[fr_sqp_U1,fr_in_sqp(i,:,U1_index(i))];
    else
            fr_sqp_U1=[fr_sqp_U1,NaN(1,10)];
    end
end

figure(1);subplot(3,2,1)
boxplot(fr_sqp_U1,group,'ColorGroup',group_color,'Positions',axis_pos,'Colors',box_color,'BoxStyle','filled','Labels',axis_label,'LabelOrientation','horizontal')
ylabel('Firing rate (Hz)')
box off;%xlim([2 13.75]);%xlim([2 13.75])
%calculate significance across consecutive trials
group_cons=[ones(1,10),2*ones(1,10)];
%across trial stats
clear frs
for i=2:6
    
    frs=[fr_in_sqp(i,:,U1_index(i)),fr_in_sqp(i+1,:,U1_index(i+1))];
    p_h_U1(i)=kruskalwallis(frs,group_cons)
    clear frs
end
%withih trial stats
group_within=[ones(1,5),2*ones(1,5)];
for i=2:7
    frs_w=[fr_in_sqp(i,1:5,U1_index(i)),fr_in_sqp(i,6:10,U1_index(i))];
    p_h_U1_w(i)=kruskalwallis(frs_w,group_within)
    clear frs_w
end

figure(1);subplot(3,2,3)
fr_sqp_U4=[];
for i=1:7
    
    if ~isnan(U4_index(i)) && (sum(~isnan(fr_in_sqp(i,1:5,U4_index(i))))>0 || sum(~isnan(fr_in_sqp(i,6:10,U4_index(i))))>0) %the fish should have traversed at least once in the two halves
            fr_sqp_U4=[fr_sqp_U4,fr_in_sqp(i,:,U4_index(i))];
    else
            fr_sqp_U4=[fr_sqp_U4,NaN(1,10)];
    end
end
boxplot(fr_sqp_U4,group,'ColorGroup',group_color,'Positions',axis_pos,'Colors',box_color,'BoxStyle','filled','Labels',axis_label,'LabelOrientation','horizontal')
ylabel('Firing rate (Hz)')
box off;%xlim([2 13.75]);%xlim([2 13.75])

group_cons=[ones(1,10),2*ones(1,10)];
%across trial stats
for i=1:6
    if ~isnan(U4_index(i))&& ~isnan(U4_index(i+1))
        frs=[fr_in_sqp(i,:,U4_index(i)),fr_in_sqp(i+1,:,U4_index(i+1))];
        p_h_U4(i)=kruskalwallis(frs,group_cons)
        clear frs
    else
        p_h_U4(i)=NaN;
    end
end
%withih trial stats
group_within=[ones(1,5),2*ones(1,5)];
for i=1:7
    if ~isnan(U4_index(i)) 
        frs_w=[fr_in_sqp(i,1:5,U4_index(i)),fr_in_sqp(i,6:10,U4_index(i))];
        p_h_U4_w(i)=kruskalwallis(frs_w,group_within)
        clear frs_w
    else
        p_h_U4_w(i)=NaN;
    end
end

figure(1);subplot(3,2,5)

fr_sqp_U5=[];
for i=1:7
    if ~isnan(U5_index(i)) && (sum(~isnan(fr_in_sqp(i,1:5,U5_index(i))))>0 || sum(~isnan(fr_in_sqp(i,6:10,U5_index(i))))>0) %the fish should have traversed at least once in the two halves
            fr_sqp_U5=[fr_sqp_U5,fr_in_sqp(i,:,U5_index(i))];
    else
            fr_sqp_U5=[fr_sqp_U5,NaN(1,10)];
    end
end

boxplot(fr_sqp_U5,group,'ColorGroup',group_color,'Positions',axis_pos,'Colors',box_color,'BoxStyle','filled','Labels',axis_label,'LabelOrientation','horizontal')
ylabel('Firing rate (Hz)')
box off;%xlim([2 13.75])
%calculate significance across consecutive trials
group_cons=[ones(1,10),2*ones(1,10)];
%across trial stats
for i=1:6
    if ~isnan(U5_index(i)) && ~isnan(U5_index(i+1))
        frs=[fr_in_sqp(i,:,U5_index(i)),fr_in_sqp(i+1,:,U5_index(i+1))];
        p_h_U5(i)=kruskalwallis(frs,group_cons)
        clear frs
    else
         p_h_U5(i)=NaN
    end
end
%withih trial stats
group_within=[ones(1,5),2*ones(1,5)];
for i=1:7
    if ~isnan(U5_index(i)) 
        frs_w=[fr_in_sqp(i,1:5,U5_index(i)),fr_in_sqp(i,6:10,U5_index(i))];
        p_h_U5_w(i)=kruskalwallis(frs_w,group_within)
        clear frs_w
    else
        p_h_U5_w(i)=NaN;
    end
end






%all firing rates

%%do the stationarity test on the time series within each trial regardless
%%of position
%calculate the firing rate in the same time bins
figure(1);subplot(3,2,2)

fr_all_U1=[];
for i=1:7
    if ~isnan(U1_index(i)) && (sum(~isnan(fr_all_spks(i,1:5,U1_index(i))))>0 || sum(~isnan(fr_all_spks(i,6:10,U1_index(i))))>0) %the fish should have traversed at least once in the two halves
            fr_all_U1=[fr_all_U1,fr_all_spks(i,:,U1_index(i))];
    else
            fr_all_U1=[fr_all_U1,NaN*ones(1,10)];
    end
end
boxplot(fr_all_U1,group,'ColorGroup',group_color,'Positions',axis_pos,'Colors',box_color,'BoxStyle','filled','Labels',axis_label,'LabelOrientation','horizontal')
ylabel('Firing rate (Hz)')
%xlabel('Trial number')
box off;%xlim([2 13.75])

group_cons=[ones(1,10),2*ones(1,10)];
%across trial stats
for i=2:6
    frs=[fr_all_spks(i,:,U1_index(i)),fr_all_spks(i+1,:,U1_index(i+1))];
    p_all_U1(i)=kruskalwallis(frs,group_cons)
    clear frs
end
%withih trial stats
group_within=[ones(1,5),2*ones(1,5)];
for i=2:7
    frs_w=[fr_all_spks(i,1:5,U1_index(i)),fr_all_spks(i,6:10,U1_index(i))];
    p_all_U1_w(i)=kruskalwallis(frs_w,group_within)
    clear frs_w
end

figure(1);subplot(3,2,4)
fr_all_U4=[];
for i=1:7
    if ~isnan(U4_index(i)) && (sum(~isnan(fr_all_spks(i,1:5,U4_index(i))))>0 || sum(~isnan(fr_all_spks(i,6:10,U4_index(i))))>0) %the fish should have traversed at least once in the two halves
            fr_all_U4=[fr_all_U4,fr_all_spks(i,:,U4_index(i))];
    else
            fr_all_U4=[fr_all_U4,NaN(1,10)];
    end
end
boxplot(fr_all_U4,group,'ColorGroup',group_color,'Positions',axis_pos,'Colors',box_color,'BoxStyle','filled','Labels',axis_label,'LabelOrientation','horizontal')
ylabel('Firing rate (Hz)')
%xlabel('Trial number')
box off;%xlim([2 13.75])

group_cons=[ones(1,10),2*ones(1,10)];
%across trial stats
for i=1:6
    if ~isnan(U4_index(i)) && ~isnan(U4_index(i+1))
        frs=[fr_all_spks(i,:,U4_index(i)),fr_all_spks(i+1,:,U4_index(i+1))];
        p_all_U4(i)=kruskalwallis(frs,group_cons)
        clear frs
    else
        p_all_U4(i)=NaN;
    end
end
%withih trial stats
group_within=[ones(1,5),2*ones(1,5)];
for i=1:7
    if ~isnan(U4_index(i))
        frs_w=[fr_all_spks(i,1:5,U4_index(i)),fr_all_spks(i,6:10,U4_index(i))];
        p_all_U4_w(i)=kruskalwallis(frs_w,group_within)
        clear frs_w
    else
        p_all_U4_w(i)=NaN;
    end
end


figure(1);subplot(3,2,6)
fr_all_U5=[];
for i=1:7

    if ~isnan(U5_index(i)) && (sum(~isnan(fr_all_spks(i,1:5,U5_index(i))))>0 || sum(~isnan(fr_all_spks(i,6:10,U5_index(i))))>0) %the fish should have traversed at least once in the two halves
            fr_all_U5=[fr_all_U5,fr_all_spks(i,:,U5_index(i))];
    else
            fr_all_U5=[fr_all_U5,NaN(1,10)];
    end
end
boxplot(fr_all_U5,group,'ColorGroup',group_color,'Positions',axis_pos,'Colors',box_color,'BoxStyle','filled','Labels',axis_label,'LabelOrientation','horizontal')
ylabel('Firing rate (Hz)')
%xlabel('Trial number')
box off;%xlim([2 13.75])

group_cons=[ones(1,10),2*ones(1,10)];
%across trial stats
for i=1:6
    if ~isnan(U5_index(i)) && ~isnan(U5_index(i+1))
        frs=[fr_all_spks(i,:,U5_index(i)),fr_all_spks(i+1,:,U5_index(i+1))];
        p_all_U5(i)=kruskalwallis(frs,group_cons)
        clear frs
    else
        p_all_U5(i)=NaN;
    end
end
%withih trial stats
group_within=[ones(1,5),2*ones(1,5)];
for i=1:7
    if ~isnan(U5_index(i))
    frs_w=[fr_all_spks(i,1:5,U5_index(i)),fr_all_spks(i,6:10,U5_index(i))];
    p_all_U5_w(i)=kruskalwallis(frs_w,group_within)
    clear frs_w
    else
       p_all_U5_w(i)=NaN;
    end
end




figure(1)
saveas(gcf,strcat(exp_date,'L-frhome-R-frall.pdf'))

%% compare the firing rate inside and outside LMs
close all
group3=[];
for i=1:7
    group3=[group3,[(2*i-1)*ones(1,10),(2*i)*ones(1,10)]];
end

group3_color=[];
for i=1:7
    group3_color=[group3_color,[ones(1,10),2*ones(1,10)]];
end


fr_in_out_all_U1=[];
for i=1:7
    if ~isnan(U1_index(i))
        fr_in_out_all_U1=[fr_in_out_all_U1,fr_in_all(i,:,U1_index(i)),fr_out_all(i,:,U1_index(i))];
    else
        fr_in_out_all_U1=[fr_in_out_all_U1,NaN(1,20)];
    end
  
end
figure(1);subplot(3,1,1);
boxplot(fr_in_out_all_U1,group3,'ColorGroup',group3_color,'Positions',axis_pos,'Colors',box_color,'BoxStyle','filled','Labels',axis_label,'LabelOrientation','horizontal')
ylabel('Firing rate (Hz)')
box off
clear frs
for i=1:7
    if ~isnan(U1_index(i))
        frs=[fr_in_all(i,:,U1_index(i)),fr_out_all(i,:,U1_index(i))];
        p_in_out_U1(i)=kruskalwallis(frs,group_cons)
        clear frs
    else
        p_in_out_U1(i)=NaN;
    end
end




fr_in_out_all_U4=[];
for i=1:7
    if ~isnan(U4_index(i))
        fr_in_out_all_U4=[fr_in_out_all_U4,fr_in_all(i,:,U4_index(i)),fr_out_all(i,:,U4_index(i))];
    else
        fr_in_out_all_U4=[fr_in_out_all_U4,NaN(1,20)];
    end
  
end
figure(1);subplot(3,1,2);
boxplot(fr_in_out_all_U4,group3,'ColorGroup',group3_color,'Positions',axis_pos,'Colors',box_color,'BoxStyle','filled','Labels',axis_label,'LabelOrientation','horizontal')
ylabel('Firing rate (Hz)')
ylim([-0.05,0.2])
box off
clear frs
for i=1:7
    if ~isnan(U4_index(i))
        frs=[fr_in_all(i,:,U4_index(i)),fr_out_all(i,:,U4_index(i))];
        p_in_out_U4(i)=kruskalwallis(frs,group_cons)
        clear frs
    else
        p_in_out_U4(i)=NaN;
    end
end

fr_in_out_all_U5=[];
for i=1:7
    if ~isnan(U5_index(i))
        fr_in_out_all_U5=[fr_in_out_all_U5,fr_in_all(i,:,U5_index(i)),fr_out_all(i,:,U5_index(i))];
    else
        fr_in_out_all_U5=[fr_in_out_all_U5,NaN(1,20)];
    end
  
end
figure(1);subplot(3,1,3);
boxplot(fr_in_out_all_U5,group3,'ColorGroup',group3_color,'Positions',axis_pos,'Colors',box_color,'BoxStyle','filled','Labels',axis_label,'LabelOrientation','horizontal')
ylabel('Firing rate (Hz)')
box off
clear frs
for i=1:7
    if ~isnan(U5_index(i))
        frs=[fr_in_all(i,:,U5_index(i)),fr_out_all(i,:,U5_index(i))];
        p_in_out_U5(i)=kruskalwallis(frs,group_cons)
        clear frs
    else
        p_in_out_U5(i)=NaN;
    end
end

figure(1)
saveas(gcf,strcat(exp_date,'in_out_fr.pdf'))
%% do the same for spike count
close all
spkcnt_in_out_all_U1=[];
for i=1:7
    if ~isnan(U1_index(i))
        spkcnt_in_out_all_U1=[spkcnt_in_out_all_U1,spkcnt_in_all(i,:,U1_index(i)),spkcnt_out_all(i,:,U1_index(i))];
    else
        spkcnt_in_out_all_U1=[spkcnt_in_out_all_U1,NaN(1,20)];
    end
  
end

figure(2);subplot(3,1,1);
boxplot(spkcnt_in_out_all_U1,group3,'ColorGroup',group3_color,'Positions',axis_pos,'Colors',box_color,'BoxStyle','filled','Labels',axis_label,'LabelOrientation','horizontal')
ylabel('Spike count')
box off
clear frs
for i=1:7
    if ~isnan(U1_index(i))
        frs=[spkcnt_in_all(i,:,U1_index(i)),spkcnt_out_all(i,:,U1_index(i))];
        p_in_out_U1(i)=kruskalwallis(frs,group_cons)
        clear frs
    else
        p_in_out_U1(i)=NaN;
    end
end




spkcnt_out_all_U4=[];
for i=1:7
    if ~isnan(U4_index(i))
        spkcnt_out_all_U4=[spkcnt_out_all_U4,spkcnt_in_all(i,:,U4_index(i)),spkcnt_out_all(i,:,U4_index(i))];
    else
        spkcnt_out_all_U4=[spkcnt_out_all_U4,NaN(1,20)];
    end
  
end
figure(2);subplot(3,1,2);
boxplot(spkcnt_out_all_U4,group3,'ColorGroup',group3_color,'Positions',axis_pos,'Colors',box_color,'BoxStyle','filled','Labels',axis_label,'LabelOrientation','horizontal')
ylabel('Spike count')
%ylim([-0.05,0.2])
box off
clear frs
for i=1:7
    if ~isnan(U4_index(i))
        frs=[spkcnt_in_all(i,:,U4_index(i)),spkcnt_out_all(i,:,U4_index(i))];
        p_in_out_U4(i)=kruskalwallis(frs,group_cons)
        clear frs
    else
        p_in_out_U4(i)=NaN;
    end
end

spkcnt_out_all_U5=[];
for i=1:7
    if ~isnan(U5_index(i))
        spkcnt_out_all_U5=[spkcnt_out_all_U5,spkcnt_in_all(i,:,U5_index(i)),spkcnt_out_all(i,:,U5_index(i))];
    else
        spkcnt_out_all_U5=[spkcnt_out_all_U5,NaN(1,20)];
    end
  
end
figure(2);subplot(3,1,3);
boxplot(spkcnt_out_all_U5,group3,'ColorGroup',group3_color,'Positions',axis_pos,'Colors',box_color,'BoxStyle','filled','Labels',axis_label,'LabelOrientation','horizontal')
ylabel('Spike count')
box off
clear frs
for i=1:7
    if ~isnan(U5_index(i))
        frs=[spkcnt_in_all(i,:,U5_index(i)),spkcnt_out_all(i,:,U5_index(i))];
        p_in_out_U5(i)=kruskalwallis(frs,group_cons)
        clear frs
    else
        p_in_out_U5(i)=NaN;
    end
end

figure(2)
saveas(gcf,strcat(exp_date,'in_out_spkcnt.pdf'))
%%
close all


All_in_U1=[]
for i=1:7
    if ~isnan(U1_index(i))
        All_in_U1=[All_in_U1, fr_in_all(i,:,U1_index(i))];
    end
end

All_out_U1=[]
for i=1:7
    if ~isnan(U1_index(i))
        All_out_U1=[All_out_U1, fr_out_all(i,:,U1_index(i))];
    end
end
group4=[ones(size(All_in_U1)),2*ones(size(All_out_U1))];
kruskalwallis([All_in_U1,All_out_U1],group4)
figure(2)
saveas(gcf,strcat(exp_date,'_U1_in_out_fr_alltrials.pdf'))
%pool all in trials

%%

%% calculate the fr near other landmarks

fr_ltp_U1=[];
for i=1:7
    if i>1 || sum(~isnan(fr_in_ltp(i,1:5,1)))>0 || sum(~isnan(fr_in_ltp(i,6:10,1)))>0 %the fish should have traversed at least once in the two halves
            fr_ltp_U1=[fr_ltp_U1,fr_in_ltp(i,:,1)];
    else
            fr_ltp_U1=[fr_ltp_U1,NaN(1,10)];
    end
end
[p_U1_ltp,anovatab,stats_U1_ltp]=kruskalwallis(fr_ltp_U1,group);
multcompare(stats_U1_ltp)

fr_ltp_U5=[];
for i=1:7
    if i==1
        index=1;
    else
        index=3;
    end
    if i~=2 || sum(~isnan(fr_in_ltp(i,1:5,index)))>0 || sum(~isnan(fr_in_ltp(i,6:10,index)))>0 %the fish should have traversed at least once in the two halves
            fr_ltp_U5=[fr_ltp_U5,fr_in_ltp(i,:,index)];
    else
            fr_ltp_U5=[fr_ltp_U5,NaN(1,10)];
    end
end
[p_U5_ltp,anovatab,stats_U5_ltp]=kruskalwallis(fr_ltp_U5,group);
multcompare(stats_U5_ltp)

fr_ltp_U4=[];
for i=1:7
    
    if i~=2 || sum(~isnan(fr_in_ltp(i,1:5,2)))>0 || sum(~isnan(fr_in_ltp(i,6:10,2)))>0 %the fish should have traversed at least once in the two halves
            fr_ltp_U4=[fr_ltp_U4,fr_in_ltp(i,:,2)];
    else
            fr_ltp_U4=[fr_ltp_U4,NaN(1,10)];
    end
end
[p_U4_ltp,anovatab,stats_U4_ltp]=kruskalwallis(fr_ltp_U4,group);
multcompare(stats_U4_ltp)

%% 3. Small triangle


fr_stp_U1=[];
for i=1:7
    if i>1 || sum(~isnan(fr_in_stp(i,1:5,1)))>0 || sum(~isnan(fr_in_stp(i,6:10,1)))>0 %the fish should have traversed at least once in the two halves
            fr_stp_U1=[fr_stp_U1,fr_in_stp(i,:,1)];
    else
            fr_stp_U1=[fr_stp_U1,NaN(1,10)];
    end
end
[p_U1_stp,anovatab,stats_U1_stp]=kruskalwallis(fr_stp_U1,group);
multcompare(stats_U1_stp)

fr_stp_U5=[];
for i=1:7
    if i==1
        index=1;
    else
        index=3;
    end
    if i~=2 || sum(~isnan(fr_in_stp(i,1:5,index)))>0 || sum(~isnan(fr_in_stp(i,6:10,index)))>0 %the fish should have traversed at least once in the two halves
            fr_stp_U5=[fr_stp_U5,fr_in_stp(i,:,index)];
    else
            fr_stp_U5=[fr_stp_U5,NaN(1,10)];
    end
end
[p_U5_stp,anovatab,stats_U5_stp]=kruskalwallis(fr_stp_U5,group);
multcompare(stats_U5_stp)

fr_stp_U4=[];
for i=1:7
    
    if i~=2 || sum(~isnan(fr_in_stp(i,1:5,2)))>0 || sum(~isnan(fr_in_stp(i,6:10,2)))>0 %the fish should have traversed at least once in the two halves
            fr_stp_U4=[fr_stp_U4,fr_in_stp(i,:,2)];
    else
            fr_stp_U4=[fr_stp_U4,NaN(1,10)];
    end
end
[p_U4_stp,anovatab,stats_U4_stp]=kruskalwallis(fr_stp_U4,group);
multcompare(stats_U4_stp)

%% 3. Small circle

fr_sc_U1=[];
for i=1:7
    if i>1 || sum(~isnan(fr_in_sc(i,1:5,1)))>0 || sum(~isnan(fr_in_sc(i,6:10,1)))>0 %the fish should have traversed at least once in the two halves
            fr_sc_U1=[fr_sc_U1,fr_in_sc(i,:,1)];
    else
            fr_sc_U1=[fr_sc_U1,NaN(1,10)];
    end
end
[p_U1_sc,anovatab,stats_U1_sc]=kruskalwallis(fr_sc_U1,group);
multcompare(stats_U1_sc)

fr_sc_U5=[];
for i=1:7
    if i==1
        index=1;
    else
        index=3;
    end
    if i~=2 || sum(~isnan(fr_in_sc(i,1:5,index)))>0 || sum(~isnan(fr_in_sc(i,6:10,index)))>0 %the fish should have traversed at least once in the two halves
            fr_sc_U5=[fr_sc_U5,fr_in_sc(i,:,index)];
    else
            fr_sc_U5=[fr_sc_U5,NaN(1,10)];
    end
end
[p_U5_sc,anovatab,stats_U5_sc]=kruskalwallis(fr_sc_U5,group);
multcompare(stats_U5_sc)

fr_sc_U4=[];
for i=1:7
    
    if i~=2 || sum(~isnan(fr_in_sc(i,1:5,2)))>0 || sum(~isnan(fr_in_sc(i,6:10,2)))>0 %the fish should have traversed at least once in the two halves
            fr_sc_U4=[fr_sc_U4,fr_in_sc(i,:,2)];
    else
            fr_sc_U4=[fr_sc_U4,NaN(1,10)];
    end
end
[p_U4_sc,anovatab,stats_U4_sc]=kruskalwallis(fr_sc_U4,group);
multcompare(stats_U4_sc)

%% 5. Large circle

fr_lc_U1=[];
for i=1:7
    if i>1 || sum(~isnan(fr_in_lc(i,1:5,1)))>0 || sum(~isnan(fr_in_lc(i,6:10,1)))>0 %the fish should have traversed at least once in the two halves
            fr_lc_U1=[fr_lc_U1,fr_in_lc(i,:,1)];
    else
            fr_lc_U1=[fr_lc_U1,NaN(1,10)];
    end
end
[p_U1_lc,anovatab,stats_U1_lc]=kruskalwallis(fr_lc_U1,group);
multcompare(stats_U1_lc)

fr_lc_U5=[];
for i=1:7
    if i==1
        index=1;
    else
        index=3;
    end
    if i~=2 || sum(~isnan(fr_in_lc(i,1:5,index)))>0 || sum(~isnan(fr_in_lc(i,6:10,index)))>0 %the fish should have traversed at least once in the two halves
            fr_lc_U5=[fr_lc_U5,fr_in_lc(i,:,index)];
    else
            fr_lc_U5=[fr_lc_U5,NaN(1,10)];
    end
end
[p_U5_lc,anovatab,stats_U5_lc]=kruskalwallis(fr_lc_U5,group);
multcompare(stats_U5_lc)

fr_lc_U4=[];
for i=1:7
    
    if i~=2 || sum(~isnan(fr_in_lc(i,1:5,2)))>0 || sum(~isnan(fr_in_lc(i,6:10,2)))>0 %the fish should have traversed at least once in the two halves
            fr_lc_U4=[fr_lc_U4,fr_in_lc(i,:,2)];
    else
            fr_lc_U4=[fr_lc_U4,NaN(1,10)];
    end
end
[p_U4_lc,anovatab,stats_U4_lc]=kruskalwallis(fr_lc_U4,group);
multcompare(stats_U4_lc)



%% 7. filter


fr_filter_U1=[];
for i=1:7
    if i>1 || sum(~isnan(fr_in_filter(i,1:5,1)))>0 || sum(~isnan(fr_in_filter(i,6:10,1)))>0 %the fish should have traversed at least once in the two halves
            fr_filter_U1=[fr_filter_U1,fr_in_filter(i,:,1)];
    else
            fr_filter_U1=[fr_filter_U1,NaN(1,10)];
    end
end
[p_U1_filter,anovatab,stats_U1_filter]=kruskalwallis(fr_filter_U1,group);
multcompare(stats_U1_filter)

fr_filter_U5=[];
for i=1:7
    if i==1
        index=1;
    else
        index=3;
    end
    if i~=2 || sum(~isnan(fr_in_filter(i,1:5,index)))>0 || sum(~isnan(fr_in_filter(i,6:10,index)))>0 %the fish should have traversed at least once in the two halves
            fr_filter_U5=[fr_filter_U5,fr_in_filter(i,:,index)];
    else
            fr_filter_U5=[fr_filter_U5,NaN(1,10)];
    end
end
[p_U5_filter,anovatab,stats_U5_filter]=kruskalwallis(fr_filter_U5,group);
multcompare(stats_U5_filter)

fr_filter_U4=[];
for i=1:7
    
    if i~=2 || sum(~isnan(fr_in_filter(i,1:5,2)))>0 || sum(~isnan(fr_in_filter(i,6:10,2)))>0 %the fish should have traversed at least once in the two halves
            fr_filter_U4=[fr_filter_U4,fr_in_filter(i,:,2)];
    else
            fr_filter_U4=[fr_filter_U4,NaN(1,10)];
    end
end
[p_U4_filter,anovatab,stats_U4_filter]=kruskalwallis(fr_filter_U4,group);
multcompare(stats_U4_filter)






