%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));
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');
    
    load(strcat('landmark_coordinates_',filename{tn}(1:8),'.mat'));
    
    %plot the landmarks
    figure(1)
    plot(tank_circle(1,:),(1200-tank_circle(2,:)),'k')
    tank_circle2(1,:)=tank_circle(1,:)-90*cos(ang);
    tank_circle2(2,:)=tank_circle(2,:)-90*sin(ang);

    hold on
    plot(hc(:,1),(1200-hc(:,2)))
    plot(stp(:,1),(1200-stp(:,2)))
    plot(ltp(:,1),(1200-ltp(:,2)))
    plot(sqp(:,1),(1200-sqp(:,2)))
    plot(sc(1,:),(1200-sc(2,:)))
    plot(lc(1,:),(1200-lc(2,:)))
    plot(filter(:,1),(1200-filter(:,2)))
    
    %% plot a 10 cm area around the landmark
    
    %look at a 75 pixel (10 cm) distance around the home
    hc2(1,:)=hc(1,:)-75;
    hc2(2,2)=hc(2,2)-75;
    hc2(2,1)=hc(2,1)+75;
    hc2(3,1)=hc2(2,1);
    hc2(3,2)=hc(3,2)+75;
    hc2(4,1)=hc2(1,1);
    hc2(4,2)=hc2(3,2);
    hc2(5,:)=hc2(1,:);
    plot(hc2(:,1),(1200-hc2(:,2)))
    hc2grid=polygrid(hc2(:,1),1200-hc2(:,2),(1/80));
    %plot(hc2grid(:,1),hc2grid(:,2),'*')
    
    %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)))
    stp2grid=polygrid(stp2(:,1),1200-stp2(:,2),(1/80));
    %plot(stp2grid(:,1),stp2grid(:,2),'*')
    
    %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)))
    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)))
    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,:))
    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,:))
    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)))
    filter2grid=polygrid(filter2(:,1),1200-filter2(:,2),(1/80));
    %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,1:a)=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));
                    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)=NaN;
                    fr_in_ltp(tn,k,j)=NaN;
                end
            end
        else
            fr_in_ltp(tn,k,1:a)=NaN;
        end
        
    end
    figure(2)
    for j=1:a
        [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
        %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));
                    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)=NaN;
                    fr_in_stp(tn,k,j)=NaN;
                end
            end
        else
            fr_in_stp(tn,k,1:a)=NaN;
        end
        
    end
    figure(2)
    for j=1:a
        [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
        %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)=NaN;
                    fr_in_sqp(tn,k,j)=NaN;
                end
            end
        else
            fr_in_sqp(tn,k,1:a)=NaN;
        end
        
    end
    figure(2)
    for j=1:a
        [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
        %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)=NaN;
                    fr_in_sc(tn,k,j)=NaN;
                end
            end
        else
            fr_in_sc(tn,k,1:a)=NaN;
        end
        
    end
    figure(2)
    for j=1:a
        [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
        %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)=NaN;
                    fr_in_lc(tn,k,j)=NaN;
                end
            end
        else
            fr_in_lc(tn,k,1:a)=NaN;
        end
        
    end
    figure(2)
    for j=1:a
        [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
        %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)=NaN;
                    fr_in_filter(tn,k,j)=NaN;
                end
            end
        else
            fr_in_filter(tn,k,1:a)=NaN;
        end
        
    end
    figure(2)
    for j=1:a
        [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
        %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 tn==1
        index_out=find(in_home==0 & in_filter==0 & out_tank==1  );
    elseif tn<6
        index_out=find(in_home==0 & in_filter==0 & in_ltp==0 & in_stp==0 & in_sqp==0 & in_sc==0 & in_lc==0& out_tank==1);
    else
        index_out=find(in_filter==0 & in_ltp==0 & in_stp==0 & in_sqp==0 & in_sc==0 & in_lc==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 tn==1
                    index_spk_out=find(in_spk_home==0 & in_spk_filter==0 & out_spk_tank==1);
                    plot(x_spike_temp(index_spk_out),y_spike_temp(index_spk_out),'bo')
                
                elseif tn<6
                   index_spk_out=find(in_spk_home==0 & 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);
                    plot(x_spike_temp(index_spk_out),y_spike_temp(index_spk_out),'bo')

                else
                    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);
                    plot(x_spike_temp(index_spk_out),y_spike_temp(index_spk_out),'bo')

                end
                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)=NaN;
                fr_out_all(tn,k,j)=NaN;
            end
        end
    else
        fr_out_all(tn,k,1:a)=NaN;
    end
    
end

%% 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));
    in_filter=inpolygon(x_temp,y_temp,filter2(:,1),1200-filter2(:,2));
   out_tank=inpolygon(x_temp,y_temp,tank_circle2(1,:),1200-tank_circle2(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 tn==1
        index_in=find(in_home== 1 | in_filter==1 | out_tank==0 )
    elseif tn<6
         index_in=find(in_home== 1 | in_filter==1 | in_ltp==1 | in_stp==1 | in_sqp==1 | in_sc==1 | in_lc==1| out_tank==0 )
    else
        index_in=find(in_filter==1 | in_ltp==1 | in_stp==1 | in_sqp==1 | in_sc==1 | in_lc==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 tn==1
                    index_spk_in=find(in_spk_home== 1 | in_spk_filter==1| out_spk_tank==0  )
                elseif tn<6
                    index_spk_in=find(in_spk_home== 1 | 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 )
                else
                    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 )
                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
%%plot landmark configs
close all
figure(2)
for tn=1:length(filename)
    scale1=1;
    scale2=1;
    a=0;
    
    subplot(1,length(filename),tn)

    plot(filter(:,1)/scale1+a,(1200-filter(:,2))/(scale2)+a,'k')
    hold on
    plot(tank_circle(1,:)/scale1+a,(1200-tank_circle(2,:))/(scale2)+a,'k')
    if filename{tn}(11)=='2'
        plot(hc(:,1)/scale1+a,(1200-hc(:,2))/(scale2)+a)
    elseif filename{tn}(11)~='7'
        plot(hc(:,1)/scale1+a,(1200-hc(:,2))/(scale2)+a)
        plot(stp(:,1)/scale1+a,(1200-stp(:,2))/(scale2)+a);
        plot(ltp(:,1)/scale1+a,(1200-ltp(:,2))/(scale2)+a);
        plot(sqp(:,1)/scale1+a,(1200-sqp(:,2))/(scale2)+a);
        plot(sc(1,:)/scale1+a,(1200-sc(2,:))/(scale2)+a);
        plot(lc(1,:)/scale1+a,(1200-lc(2,:))/(scale2)+a);
    else
        %plot(hc(:,1)/scale1,(1200-hc(:,2))/(scale2),':')
        plot(stp(:,1)/scale1+a,(1200-stp(:,2))/(scale2)+a);
        plot(ltp(:,1)/scale1+a,(1200-ltp(:,2))/(scale2)+a);
        plot(sqp(:,1)/scale1+a,(1200-sqp(:,2))/(scale2)+a);
        plot(sc(1,:)/scale1+a,(1200-sc(2,:))/(scale2)+a);
        plot(lc(1,:)/scale1+a,(1200-lc(2,:))/(scale2)+a);
    end
    axis square
    axis off
end
figure(2)
%saveas(gcf,'landmarkconfigs.pdf');
%%
group=[];

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

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

for i=1:12
    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:6
    
    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 for U1,U2,U4, U6
U1_index=[1,1,1,1,5,1];
U2_index=[2,2,2,2,1,NaN];
U4_index=[NaN,4,3,4,4,NaN];
U6_index=[NaN,5,5,3,2,2];

group_cons=[ones(1,10),2*ones(1,10)];
group_within=[ones(1,5),2*ones(1,5)];

%for U2 around the house
%1.home
fr_hc_U1=[];
for i=1:6
    if ~isnan(U1_index(i))
        fr_hc_U1=[fr_hc_U1,fr_in_hc(i,:,U1_index(i))];
    else
        fr_hc_U1=[fr_hc_U1,NaN(1,10)];
    end
  
end
figure(1);subplot(4,2,1);
boxplot(fr_hc_U1,group,'ColorGroup',group_color,'Positions',axis_pos,'Colors',box_color,'BoxStyle','filled','Labels',axis_label,'LabelOrientation','horizontal')
ylabel('Firing rate (Hz)')
% title('U2')
%xlabel('Trial number')
box off
clear frs
for i=1:5
    if ~isnan(U1_index(i)) && ~isnan(U1_index(i+1))
        frs=[fr_in_hc(i,:,U1_index(i)),fr_in_hc(i+1,:,U1_index(i+1))];
        p_h_U1(i)=kruskalwallis(frs,group_cons)
        clear frs
    else
        p_h_U1(i)=NaN;
    end
end
clear frs_w
%withih trial stats
for i=1:6
    if ~isnan(U1_index(i))
        frs_w=[fr_in_hc(i,1:5,U1_index(i)),fr_in_hc(i,6:10,U1_index(i))];
        p_h_U1_w(i)=kruskalwallis(frs_w,group_within)
        clear frs_w
    else
        p_h_U1_w(i)=NaN;
    end
end

fr_hc_U2=[];
for i=1:6
    if ~isnan(U2_index(i))
        fr_hc_U2=[fr_hc_U2,fr_in_hc(i,:,U2_index(i))];
    else
        fr_hc_U2=[fr_hc_U2,NaN(1,10)];
    end
  
end
figure(1);subplot(4,2,3);
boxplot(fr_hc_U2,group,'ColorGroup',group_color,'Positions',axis_pos,'Colors',box_color,'BoxStyle','filled','Labels',axis_label,'LabelOrientation','horizontal')
ylabel('Firing rate (Hz)')
% title('U2')
%xlabel('Trial number')
box off
for i=1:5
    if ~isnan(U2_index(i)) && ~isnan(U2_index(i+1))
        frs=[fr_in_hc(i,:,U2_index(i)),fr_in_hc(i+1,:,U2_index(i+1))];
        p_h_U2(i)=kruskalwallis(frs,group_cons)
        clear frs
    else
        p_h_U2(i)=NaN;
    end
end
clear frs_w
%withih trial stats
for i=1:6
    if ~isnan(U2_index(i))
        frs_w=[fr_in_hc(i,1:5,U2_index(i)),fr_in_hc(i,6:10,U2_index(i))];
        p_h_U2_w(i)=kruskalwallis(frs_w,group_within)
        clear frs_w
    else
        p_h_U2_w(i)=NaN;
    end
end

fr_hc_U4=[];
for i=1:6
    if ~isnan(U4_index(i))
        fr_hc_U4=[fr_hc_U4,fr_in_hc(i,:,U4_index(i))];
    else
        fr_hc_U4=[fr_hc_U4,NaN(1,10)];
    end
  
end
figure(1);subplot(4,2,5);
boxplot(fr_hc_U4,group,'ColorGroup',group_color,'Positions',axis_pos,'Colors',box_color,'BoxStyle','filled','Labels',axis_label,'LabelOrientation','horizontal')
ylabel('Firing rate (Hz)')
% title('U2')
%xlabel('Trial number')
box off
for i=1:5
    if ~isnan(U4_index(i)) && ~isnan(U4_index(i+1))
        frs=[fr_in_hc(i,:,U4_index(i)),fr_in_hc(i+1,:,U4_index(i+1))];
        p_h_U4(i)=kruskalwallis(frs,group_cons)
        clear frs
    else
        p_h_U4(i)=NaN;
    end
end
clear frs_w
%withih trial stats
for i=1:6
    if ~isnan(U4_index(i))
        frs_w=[fr_in_hc(i,1:5,U4_index(i)),fr_in_hc(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

fr_hc_U6=[];
for i=1:6
    if ~isnan(U6_index(i))
        fr_hc_U6=[fr_hc_U6,fr_in_hc(i,:,U6_index(i))];
    else
        fr_hc_U6=[fr_hc_U6,NaN(1,10)];
    end
  
end
figure(1);subplot(4,2,7);
boxplot(fr_hc_U6,group,'ColorGroup',group_color,'Positions',axis_pos,'Colors',box_color,'BoxStyle','filled','Labels',axis_label,'LabelOrientation','horizontal')
ylabel('Firing rate (Hz)')
% title('U2')
%xlabel('Trial number')
box off
for i=1:5
    if ~isnan(U6_index(i)) && ~isnan(U6_index(i+1))
        frs=[fr_in_hc(i,:,U6_index(i)),fr_in_hc(i+1,:,U6_index(i+1))];
        p_h_U6(i)=kruskalwallis(frs,group_cons)
        clear frs
    else
        p_h_U6(i)=NaN;
    end
end
clear frs_w
%withih trial stats
for i=1:6
    if ~isnan(U6_index(i))
        frs_w=[fr_in_hc(i,1:5,U6_index(i)),fr_in_hc(i,6:10,U6_index(i))];
        p_h_U6_w(i)=kruskalwallis(frs_w,group_within)
        clear frs_w
    else
        p_h_U6_w(i)=NaN;
    end
end


% all firing rates
fr_all_U1=[];
for i=1:6
    if ~isnan(U1_index(i))
        if 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(1,10)];
        end
    else
        fr_all_U1=[fr_all_U1,NaN(1,10)];
    end
end
figure(1);subplot(4,2,2);
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)')
% title('U2')
%xlabel('Trial number')
box off
for i=1:5
    if ~isnan(U1_index(i)) && ~isnan(U1_index(i+1))
        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
    else
        p_all_U1(i)=NaN;
    end
end
clear frs_w
%withih trial stats
for i=1:6
    if ~isnan(U1_index(i))
        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
    else
        p_all_U1_w(i)=NaN;
    end
end

fr_all_U2=[];
for i=1:6
    if ~isnan(U2_index(i))
        if sum(~isnan(fr_all_spks(i,1:5,U2_index(i))))>0 && sum(~isnan(fr_all_spks(i,6:10,U2_index(i))))>0 %the fish should have traversed at least once in the two halves
                fr_all_U2=[fr_all_U2,fr_all_spks(i,:,U2_index(i))];
        else
                fr_all_U2=[fr_all_U2,NaN(1,10)];
        end
    else
        fr_all_U2=[fr_all_U2,NaN(1,10)];
    end
end
figure(1);subplot(4,2,4);
boxplot(fr_all_U2,group,'ColorGroup',group_color,'Positions',axis_pos,'Colors',box_color,'BoxStyle','filled','Labels',axis_label,'LabelOrientation','horizontal')
ylabel('Firing rate (Hz)')
% title('U2')
%xlabel('Trial number')
box off
clear frs
for i=1:5
    if ~isnan(U2_index(i)) && ~isnan(U2_index(i+1))
        frs=[fr_all_spks(i,:,U2_index(i)),fr_all_spks(i+1,:,U2_index(i+1))];
        p_all_U2(i)=kruskalwallis(frs,group_cons)
        clear frs
    else
        p_all_U2(i)=NaN;
    end
end
clear frs_w
%withih trial stats
for i=1:6
    if ~isnan(U2_index(i))
        frs_w=[fr_all_spks(i,1:5,U2_index(i)),fr_all_spks(i,6:10,U2_index(i))];
        p_all_U2_w(i)=kruskalwallis(frs_w,group_within)
        clear frs_w
    else
        p_all_U2_w(i)=NaN;
    end
end

fr_all_U4=[];
for i=1:6
    if ~isnan(U4_index(i))
        if 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
    else
        fr_all_U4=[fr_all_U4,NaN(1,10)];
    end
end
figure(1);subplot(4,2,6);
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)')
% title('U2')
%xlabel('Trial number')
box off
for i=1:5
    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
clear frs_w
%withih trial stats
for i=1:6
    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

fr_all_U6=[];
for i=1:6
    if ~isnan(U6_index(i))
        if sum(~isnan(fr_all_spks(i,1:5,U6_index(i))))>0 && sum(~isnan(fr_all_spks(i,6:10,U6_index(i))))>0 %the fish should have traversed at least once in the two halves
                fr_all_U6=[fr_all_U6,fr_all_spks(i,:,U6_index(i))];
        else
                fr_all_U6=[fr_all_U6,NaN(1,10)];
        end
    else
        fr_all_U6=[fr_all_U6,NaN(1,10)];
    end
end
figure(1);subplot(4,2,8);
boxplot(fr_all_U6,group,'ColorGroup',group_color,'Positions',axis_pos,'Colors',box_color,'BoxStyle','filled','Labels',axis_label,'LabelOrientation','horizontal')
ylabel('Firing rate (Hz)')
% title('U2')
%xlabel('Trial number')
box off
clear frs
for i=1:5
    if ~isnan(U6_index(i)) && ~isnan(U6_index(i+1))
        frs=[fr_all_spks(i,:,U6_index(i)),fr_all_spks(i+1,:,U6_index(i+1))];
        p_all_U6(i)=kruskalwallis(frs,group_cons)
        clear frs
    else
        p_all_U6(i)=NaN;
    end
end
clear frs_w
%withih trial stats
for i=1:6
    if ~isnan(U6_index(i))
        frs_w=[fr_all_spks(i,1:5,U6_index(i)),fr_all_spks(i,6:10,U6_index(i))];
        p_all_U6_w(i)=kruskalwallis(frs_w,group_within)
        clear frs_w
    else
        p_all_U6_w(i)=NaN;
    end
end

figure(1)
%saveas(gcf,strcat(exp_date,'L-frhome-R-frall.pdf'))
%%
%% plot comparing only across trials
group=[];

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

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

for i=1:12
    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:6
    
    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 for U1,U2,U4, U6
U1_index=[1,1,1,1,5,1];
U2_index=[2,2,2,2,1,NaN];
U4_index=[NaN,4,3,4,4,NaN];
U6_index=[NaN,5,5,3,2,2];

group_cons=[ones(1,10),2*ones(1,10)];
group_within=[ones(1,5),2*ones(1,5)];

%for U2 around the house
%1.home
fr_hc_U1=[];
for i=1:6
    if ~isnan(U1_index(i))
        fr_hc_U1=[fr_hc_U1,fr_in_hc(i,:,U1_index(i))];
    else
        fr_hc_U1=[fr_hc_U1,NaN(1,10)];
    end
  
end
figure(1);subplot(4,1,1);
boxplot(fr_hc_U1,group2)
ylabel('Firing rate (Hz)')
box off
clear frs
%withih trial stats
hold on
fr_hc_U2=[];
for i=1:6
    if ~isnan(U2_index(i))
        fr_hc_U2=[fr_hc_U2,fr_in_hc(i,:,U2_index(i))];
    else
        fr_hc_U2=[fr_hc_U2,NaN(1,10)];
    end
  
end
figure(1);subplot(4,1,2);
boxplot(fr_hc_U2,group2)
ylabel('Firing rate (Hz)')
% title('U2')
%xlabel('Trial number')
box off
fr_hc_U4=[];
for i=1:6
    if ~isnan(U4_index(i))
        fr_hc_U4=[fr_hc_U4,fr_in_hc(i,:,U4_index(i))];
    else
        fr_hc_U4=[fr_hc_U4,NaN(1,10)];
    end
  
end
figure(1);subplot(4,1,3);
boxplot(fr_hc_U4,group2)
ylabel('Firing rate (Hz)')
% title('U2')
%xlabel('Trial number')
box off

fr_hc_U6=[];
for i=1:6
    if ~isnan(U6_index(i))
        fr_hc_U6=[fr_hc_U6,fr_in_hc(i,:,U6_index(i))];
    else
        fr_hc_U6=[fr_hc_U6,NaN(1,10)];
    end
  
end
figure(1);subplot(4,1,4);
boxplot(fr_hc_U6,group2)
ylabel('Firing rate (Hz)')
% title('U2')
%xlabel('Trial number')

%%
% all firing rates
fr_all_U1=[];
for i=1:6
    if ~isnan(U1_index(i))
        if 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(1,10)];
        end
    else
        fr_all_U1=[fr_all_U1,NaN(1,10)];
    end
end
figure(1);subplot(4,2,2);
boxplot(fr_all_U1,group2,'ColorGroup',group_color,'Positions',axis_pos,'Colors',box_color,'BoxStyle','filled','Labels',axis_label,'LabelOrientation','horizontal')
ylabel('Firing rate (Hz)')
% title('U2')
%xlabel('Trial number')
box off
for i=1:5
    if ~isnan(U1_index(i)) && ~isnan(U1_index(i+1))
        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
    else
        p_all_U1(i)=NaN;
    end
end
clear frs_w
%withih trial stats
for i=1:6
    if ~isnan(U1_index(i))
        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
    else
        p_all_U1_w(i)=NaN;
    end
end

fr_all_U2=[];
for i=1:6
    if ~isnan(U2_index(i))
        if sum(~isnan(fr_all_spks(i,1:5,U2_index(i))))>0 && sum(~isnan(fr_all_spks(i,6:10,U2_index(i))))>0 %the fish should have traversed at least once in the two halves
                fr_all_U2=[fr_all_U2,fr_all_spks(i,:,U2_index(i))];
        else
                fr_all_U2=[fr_all_U2,NaN(1,10)];
        end
    else
        fr_all_U2=[fr_all_U2,NaN(1,10)];
    end
end
figure(1);subplot(4,2,4);
boxplot(fr_all_U2,group2,'ColorGroup',group_color,'Positions',axis_pos,'Colors',box_color,'BoxStyle','filled','Labels',axis_label,'LabelOrientation','horizontal')
ylabel('Firing rate (Hz)')
% title('U2')
%xlabel('Trial number')
box off
clear frs
for i=1:5
    if ~isnan(U2_index(i)) && ~isnan(U2_index(i+1))
        frs=[fr_all_spks(i,:,U2_index(i)),fr_all_spks(i+1,:,U2_index(i+1))];
        p_all_U2(i)=kruskalwallis(frs,group_cons)
        clear frs
    else
        p_all_U2(i)=NaN;
    end
end
clear frs_w
%withih trial stats
for i=1:6
    if ~isnan(U2_index(i))
        frs_w=[fr_all_spks(i,1:5,U2_index(i)),fr_all_spks(i,6:10,U2_index(i))];
        p_all_U2_w(i)=kruskalwallis(frs_w,group_within)
        clear frs_w
    else
        p_all_U2_w(i)=NaN;
    end
end

fr_all_U4=[];
for i=1:6
    if ~isnan(U4_index(i))
        if 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
    else
        fr_all_U4=[fr_all_U4,NaN(1,10)];
    end
end
figure(1);subplot(4,2,6);
boxplot(fr_all_U4,group2,'ColorGroup',group_color,'Positions',axis_pos,'Colors',box_color,'BoxStyle','filled','Labels',axis_label,'LabelOrientation','horizontal')
ylabel('Firing rate (Hz)')
% title('U2')
%xlabel('Trial number')
box off
for i=1:5
    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
clear frs_w
%withih trial stats
for i=1:6
    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

fr_all_U6=[];
for i=1:6
    if ~isnan(U6_index(i))
        if sum(~isnan(fr_all_spks(i,1:5,U6_index(i))))>0 && sum(~isnan(fr_all_spks(i,6:10,U6_index(i))))>0 %the fish should have traversed at least once in the two halves
                fr_all_U6=[fr_all_U6,fr_all_spks(i,:,U6_index(i))];
        else
                fr_all_U6=[fr_all_U6,NaN(1,10)];
        end
    else
        fr_all_U6=[fr_all_U6,NaN(1,10)];
    end
end
figure(1);subplot(4,2,8);
boxplot(fr_all_U6,group2,'ColorGroup',group_color,'Positions',axis_pos,'Colors',box_color,'BoxStyle','filled','Labels',axis_label,'LabelOrientation','horizontal')
ylabel('Firing rate (Hz)')
% title('U2')
%xlabel('Trial number')
box off
clear frs
for i=1:5
    if ~isnan(U6_index(i)) && ~isnan(U6_index(i+1))
        frs=[fr_all_spks(i,:,U6_index(i)),fr_all_spks(i+1,:,U6_index(i+1))];
        p_all_U6(i)=kruskalwallis(frs,group_cons)
        clear frs
    else
        p_all_U6(i)=NaN;
    end
end
clear frs_w
%withih trial stats
for i=1:6
    if ~isnan(U6_index(i))
        frs_w=[fr_all_spks(i,1:5,U6_index(i)),fr_all_spks(i,6:10,U6_index(i))];
        p_all_U6_w(i)=kruskalwallis(frs_w,group_within)
        clear frs_w
    else
        p_all_U6_w(i)=NaN;
    end
end

figure(2)

%%

%do the same thing for other 
%2. Large triangle

%3. Small triangle

%4. Small circle

%5. Large circle

%6. small square

%7. filter

% 8. outside all landmarks

%% compare the firing rate inside and outside LMs

group3=[];
for i=1:6
    group3=[group3,[(2*i-1)*ones(1,10),(2*i)*ones(1,10)]];
end

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


fr_in_out_all_U1=[];
for i=1:6
    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(4,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:6
    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_U2=[];
for i=1:6
    if ~isnan(U2_index(i))
        fr_in_out_all_U2=[fr_in_out_all_U2,fr_in_all(i,:,U2_index(i)),fr_out_all(i,:,U2_index(i))];
    else
        fr_in_out_all_U2=[fr_in_out_all_U2,NaN(1,20)];
    end
  
end
figure(1);subplot(4,1,2);
boxplot(fr_in_out_all_U2,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:6
    if ~isnan(U2_index(i))
        frs=[fr_in_all(i,:,U2_index(i)),fr_out_all(i,:,U2_index(i))];
        p_in_out_U2(i)=kruskalwallis(frs,group_cons)
        clear frs
    else
        p_in_out_U2(i)=NaN;
    end
end

fr_in_out_all_U4=[];
for i=1:6
    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(4,1,3);
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)')
box off
clear frs
for i=1:6
    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_U6=[];
for i=1:6
    if ~isnan(U6_index(i))
        fr_in_out_all_U6=[fr_in_out_all_U6,fr_in_all(i,:,U6_index(i)),fr_out_all(i,:,U6_index(i))];
    else
        fr_in_out_all_U6=[fr_in_out_all_U6,NaN(1,20)];
    end
  
end
figure(1);subplot(4,1,4);
boxplot(fr_in_out_all_U6,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:6
    if ~isnan(U6_index(i))
        frs=[fr_in_all(i,:,U6_index(i)),fr_out_all(i,:,U6_index(i))];
        p_in_out_U6(i)=kruskalwallis(frs,group_cons)
        clear frs
    else
        p_in_out_U6(i)=NaN;
    end
end

figure(1)
%saveas(gcf,strcat(exp_date,'in_out_fr.pdf'))

%save('C4_stationarytest_vars.mat')