%ConteScriptHasselmo2008EJN.m - This script generates the figures presented in a 2008 
%European Journal of Neuroscience paper.
%CITATION: Hasselmo, M.E. (2008) Temporally structured replay of neural activity in a 
%model of entorhinal cortex, hippocampus and postsubiculum. Eur. J. Neurosci. 28:1301-1315.

%Script written by
%Michael E. Hasselmo, Boston University, copyright 2007, 2008, 2009, 2010.

%This script does a single run as in the earlier figures of that paper (but can be reset
%to run multiple runs of 10 as in later figures.  Note that the connectivity from grid
%cells to place cells is reset randomly, so the pattern is different on each run.
%Activity during encoding (waking) is plotted in graphics on left.  
%Activity during Retrieval (REM) is plotted on the right.  
%Each row has different information:
%1. Location and speed versus time are plotted in the top plots.
%Location looks like a sine wave because it is one dimension of location in the circle.
%Speed is a flat line with downward deflections for slowing at reward.
%2. Place cell spiking activity vs time is in the middle plots (Encoding and Retrieval),
%3. head direction cell activity vs. time are in the bottom plots, 
%4. head direction vectors vs. space are in the smaller square plots next to 
%head direction vs. time.
%5. The plot on the upper right shows the perceived x,y position during encoding (gray line)
%and during retrieval (black line).  This shows where the rat "thinks" it is located
%based on the inverse of the grid cell oscillation phase.  
%Can increase to 10 examples by setting NumPtestsTot=10 instead of 1.

%Note that there is a "while" function that makes it continue running tests until
%the internal retrieval (right) matches the number of cycles of running on the track during
%encoding.  

%The pattern of grid cell activity vs space can be seen by setting plotGrids=1.
%Note that there are 75 grid cells plotted during creation, encoding and retrieval,
%so this slows the script down.
%The pattern of place cell activity vs spcace can be seen by setting plotPlaceCells=1.
%This plots an array of numPlotPlace rows and numPlotPlace columns for cells during encoding and
%retrieval (not during creation).
%Setting plotAllPlace=1 plots the sum of all 400 place cells at once.
%Note that the plots of place cell spiking versus time plot the number of
%cells determined by numPlaceToPlot=10.
 
close all
clear all

tic

%TIME
fullT=1200; %START.  FULLT.  
kplot=60; k=0; oldt=1; %Plotting
kplace=1; kp=kplace; %Interval of place cell formation. Set kp to kplace so that first place will be same?

%HEAD DIRECTION CELL PARAMETERS
dt=0.02;
h=[0; pi/3; 2*pi/3; pi; 4*pi/3; 5*pi/3];
HHst=[cos(h) sin(h)];

hd=zeros(length(h),fullT); %Head direction cell activity at each time.
HDnum=length(h);

%RUN REPEATED TESTS WITH DIFFERENT PARAMETERS
numPsizesTot=1; %Number of different sizes of numP or other variable. Get memory errors when numP>36.
numPtestsTot=1; % 3%Number of individual tests at a single variable value.
numCircRun=zeros(numPsizesTot,numPtestsTot+1);
numPlaceToPlot=10; %This sets number of place cells that will be plotted in place cell plots.

%FIGURE FOR PLOTTING TRAJECTORIES
figSize=numPsizesTot*numPtestsTot;
figure('name','x y','Numbertitle','off','Position',[910 700 160 100+70*figSize]); ax14=gca;
axes(ax14);

rc=[1 0 0; 1 1 0; 1 0 1; 1 1 1; 0 1 1; 0 1 0.5; 0 1 0; 1 1 0; 1 0.5 0; 1 0 0];
colorSpikes=0;
figure('name','x,y,speed vs time','Numbertitle','off','Position',[50 700 300 100+70]); ax111=gca;
figure('name','x,y,speed vs time retrieval','Numbertitle','off','Position',[600 700 300 100+70]); ax112=gca;
figure('name','Place Enc v time','Numbertitle','off','Position',[50 400 300 100+70*figSize]); ax16b=gca;
figure('name','Place Ret v time','Numbertitle','off','Position',[600 400 300 100+70*figSize]); ax21b=gca;
figure('name','HD encoding vs time','Numbertitle','off','Position',[50 100 300 100+70]); ax18b=gca;
figure('name','HD retrieval vs time','Numbertitle','off','Position',[600 100 300 100+70*figSize]); ax19b=gca;
figure('name','HD vectors Enc','Numbertitle','off','Position',[360 100 160 100+70]); ax18=gca;
figure('name','HD vectors Ret','Numbertitle','off','Position',[910 100 160 100+70*figSize]); ax19=gca;

strengHD=1.0;

for numPsizes=1:numPsizesTot;
    numP=20+numPsizes*2; %Started with 4+numPsizes*2 to give 6,8,10,12,14,16....
    %Get good function with numP=18+numPsizes*2 (e.g. starting point of numP=20!!
    
    for numPtests=1:numPtestsTot;
        while numCircRun(numPsizes,numPtests+1)<floor(strengHD*4)-1;
            numCircRun(numPsizes,numPtests+1)=0;
            placeSpike2=zeros(numP^2,fullT);
            placeSpike3=zeros(numP^2,fullT);
            plotrphi=zeros(length(h)/2,fullT);
            hd3=zeros(length(h),fullT);

%CREATE TRAJECTORIES - 
speed=5; revs=0.5; mom=0.99; %Rat movement
xv=zeros(2,fullT); 
dxv=zeros(2,fullT); %Position vectors.

datainput=5; %datainput=1 - Hafting, datainput=4  = circular, 5=circular with stopping
 
minT=1;

if datainput==1;
    minT=8; %1208 
load Hafting_Fig2c_Trial2.mat
dt=0.02; sumspeed=0;
xv(1,1:fullT+1)=pos_x(minT:minT+fullT,1); %Take position from data file pos_x at step t+minT.
xv(2,1:fullT+1)=pos_y(minT:minT+fullT,1); %Take position from data file pos_y at step t
dxv(1,1:fullT+1)=(pos_x(minT:minT+fullT,1)-pos_x(minT-1:minT+fullT-1,1))/dt;
dxv(2,1:fullT+1)=(pos_y(minT:minT+fullT,1)-pos_y(minT-1:minT+fullT-1,1))/dt; %Velocity vector over time.
        for t=1:fullT
            speed=sqrt(dxv(1,t+1)^2+dxv(2,t+1)^2);
            sumspeed=sumspeed+speed;
        end
        sumspeed*dt %Print total distance covered.
        fullT*dt     %Print total time run.
        avespeed=sumspeed/fullT %Print average speed.
        xmax=100; xmin=-100;  ymax=100; ymin=-100;
%elseif datainput=0;
%    minT=2; %For random trajectory, must set minT=2 so that t-1 is at index 1.
elseif datainput==4;  %Circle without stopping points.
    dxv=[0.35; 0.35];
    dt=0.02; sumspeed=0;
    xmax=60; xmin=-60; ymax=60; ymin=-60; %Environment size for circular track.
    radius=47.5;
    for t=1:fullT
        if t==1 xv=[radius; 0]; phi(:,1)=(HHst*xv(:,1)>0).*(HHst*xv(:,1)); 
        end; 
        cirangle=atan2(xv(2,t),xv(1,t));
        dxv(1,t+1)=(1.0+0.0*rand)*sin(cirangle-1/96)/dt; %/dt; %1/4 larger inward vs. 1/ 64 larger outward.  Use 1/32 to get larger inward movement?
        dxv(2,t+1)=-(1.0+0.0*rand)*cos(cirangle-1/96)/dt; %/dt; %For 20*sin -1/5 is best.  for 10*sin
        xv(:,t+1) = xv(:,t) + dt*dxv(:,t+1) ;
        speed=sqrt(dxv(1,t+1)^2+dxv(2,t+1)^2);
        sumspeed=sumspeed+speed;
    end
    sumspeed*dt  %Print total distance covered.
    fullT*dt  %Print total time run.
avespeed=sumspeed/fullT  %Print average speed.
pos_x=xv(1,:);
pos_y=xv(2,:);
elseif datainput==5; %Circle with stopping points.
    dxv=[0; -0.35];
    reward=zeros(1,fullT);
    speed=zeros(1,fullT);
    dt=0.02; sumspeed=0;
    xmax=60; xmin=-60; ymax=60; ymin=-60; %Environment size for circular track.
    radius=47.5;
    stopangle=0; 
    ccstop=1;
    for t=1:fullT
        if t==1 xv=[radius; 0]; phi(:,1)=(HHst*xv(:,1)>0).*(HHst*xv(:,1)); 
        end; 
        cirangle=atan2(xv(2,t),xv(1,t));
        if t>1
        cirangleOLD=atan2(xv(2,t-1),xv(1,t-1));
        dangle=cirangleOLD-cirangle;
        if abs(dangle)<pi/2
        stopangle=stopangle+dangle;
        end
        end
        if xv(2,t)<12 && xv(2,t)>2 && xv(1,t)>20
        ccstop=ccstop+1;
            dxv(1,t+1)=(0.5+0.0*rand)*sin(cirangle-1/96)/dt; %This causes Slower movement
            dxv(2,t+1)=-(0.5+0.0*rand)*cos(cirangle-1/96)/dt;
            xv(:,t+1) = xv(:,t) + dt*dxv(:,t+1) ; %This updates with slow movement
            [ccstop t+1 xv(1,t) xv(2,t)]
            if ccstop==30;
                reward(1,t)=1;
                stopangle=0;
                ccstop=1;
            end
        else
        dxv(1,t+1)=(1.0+0.0*rand)*sin(cirangle-1/96)/dt; %/dt; %1/4 larger inward vs. 1/ 64 larger outward.  Use 1/32 to get larger inward movement?
        dxv(2,t+1)=-(1.0+0.0*rand)*cos(cirangle-1/96)/dt; %/dt; %For 20*sin -1/5 is best.  for 10*sin
        xv(:,t+1) = xv(:,t) + dt*dxv(:,t+1) ;
        end
        speed(1,t+1)=sqrt(dxv(1,t+1)^2+dxv(2,t+1)^2);
        sumspeed=sumspeed+speed(1,t+1);
    end
    sumspeed*dt  %Print total distance covered.
    fullT*dt  %Print total time run.
avespeed=sumspeed/fullT  %Print average speed.
pos_x=xv(1,:);
pos_y=xv(2,:);
end

axes(ax111);
tt=[1:fullT+1];
plot(tt,pos_x,tt,speed);
xlim([1 fullT]);
ylim([-55 55]);


%SOME GRAPHICS PARAMETERS
plotH=1; %Sets whether head direction vectors are plotted at end.
gy=0.7; %Sets trajectory color to gray.
figsize=200;
plotGrids=0; %Plot grid cell firing patterns (slow).
plotPlaceCells=0; %Plot a subset of place cells (slow).
numPlotPlace=5; %Number of rows and columns of place cells to plot when plotPlaceCells=1;
plotAllPlace=1; %Plot the sum of all place cell firing.
plotEnc=1;
plotRet=0;


%******INITIAL CELL CREATION************************
%CREATE GRID CELLS FOR INITIAL PLACE CELL CREATION
phi=zeros(HDnum,fullT); %Vector for phase of dendrites (in radians).
phis=zeros(HDnum/2,fullT);  %Same soma phase in all parts of vector, but easier to compute with this.

T=zeros(1,fullT);  T(1,1)=1;
iphi=[0; 0; 0]; %Vector of initial phases of dendrites (in radians) %iphi=[pi -0.5*pi -0.5*pi];
B=0.00385; %Scaling factor for effect of head direction on membrane oscillation frequency.
%strengHD=0; %Strength of head direction input during retrieval.
numF=3; %Number of different frequencies (usually 3 in most simulations)
baseF=0; %Baseline frequency for range. (usually 0 in most simulations using cos(phi-phis)
multFF=1; %Addition for going through frequencies.
numPhi=5; %Number of different starting phases in each dimension (default=5 for most simulations)
%numP=20; %Number of rows (and columns) of place cells.
thresh=0.3; %Threshold for grid cells.
pOff=0.2; %Grid cell phase offset.
stdlim=2.7; %2.7 best for circle %Determines size of accepted place fields.
xspike=zeros(2,numF,numPhi^2,fullT); %Location of grid cellspiking
tspike=zeros(numF*numPhi^2,fullT); %Spiking of grid cells during creation.
tspike2=zeros(1,numF*numPhi^2,fullT); %Spiking of grid cells during encoding.
timep(1,numF*numPhi^2).spike=[];
allPlotPlaceX=zeros(1,1); %fullT
allPlotPlaceY=zeros(1,1);
sumPlace=zeros(1,fullT);

if plotGrids==1
figure('name','diff grids','Numbertitle','off','Position',[250 50 numPhi*150 numPhi*150]); ax15=gca;
end

%CREATION OF GRID CELLS AND PLACE CELLS.

useDXV=1;
%Cycle through frequencies and phases for creating grid cells.
for ff=1:numF; %GRID CELL INITIAL CREATION.
  f=baseF+multFF*ff;  
  for pp=1:numPhi;
    for p2=1:numPhi;
    ix=[xv(1,1)+(pp+ff+pOff)*(sqrt(3)/(f*B*(numPhi))); xv(2,1)+(p2+ff+pOff)*(sqrt(3)/(f*B*(numPhi)))]; %f*G=H=sqrt(3)/2*B. G=sqrt(3)/2*f*B
    hx=((HHst*ix)>0).*(HHst*(ix));
    LH=length(h)/2; MH=length(h);
       if useDXV==1; phis(:,1)=2*pi*f*dt; phi(:,1)=2*pi*f*(dt+B*hx); end;

         for t=2:fullT;
           if useDXV==1; phis(:,t)=phis(:,t-1)+2*pi*f*dt; phi(:,t)=phi(:,t-1)+2*pi*f*dt*(1+B*(HHst*dxv(:,t)>0).*(HHst*dxv(:,t))); 
           elseif useDXV==0; phis(:,t)=2*pi*f*t; phi(:,t)=2*pi*f*(t+B*HHst*(xv(:,t)-ix));
           end; %Creation of phase representation.
         end

LH=length(h)/2; MH=length(h);
g=prod(cos((phi(1:LH,:)-phi(LH+1:MH,:))));

spiket=g>thresh; %Spiking of grid cell.

tspike((ff-1)*numPhi^2+(pp-1)*numPhi+p2,1:fullT)=spiket; %Store spiking of all grid cells.

if plotGrids==1 %OPTIONAL GRID CELL FIGURE
axes(ax15); 
axes(ax15);
plotSpikeX=xv(1,1:fullT).*spiket;
plotSpikeY=xv(2,1:fullT).*spiket;
subplot(numPhi,numPhi,(p2-1)*numPhi+pp,'replace'); ax15=gca;
plot(plotSpikeX,plotSpikeY,'r.');
xlim([xmin-1 xmax+1])
ylim([ymin-1 ymax+1])
set(gca,'XTickLabel',{' '})
set(gca,'YTickLabel',{' '})
drawnow
end %End if plotGrids==1

toc
        end %End for p2=1:numPhi
    end %End for pp=1:numPhi
end %End for f=1:numF

%INITIAL PLACE CELLS CREATION - place cells.
%numP=6; Number set with numP above.
placeSpike=zeros(numP^2,fullT);

WHP=zeros(HDnum,numP^2); %Creates matrix from place cells to HD cells.
WPG=zeros(numP^2,numF*(numPhi^2)); %WPG matrix will store associations between grid cells and place cells.
numGtoP=3;
nP=1;
nnnn=1;
while nP<(numP^2+1)
    nnnn=nnnn+1;
    if nnnn>50000; break; end;
    gp=1;
    pPhi=ones(gp,1);
for ppp=1:numF*numPhi^2;  %SELECT THREE GRID CELLS AT RANDOM WITHOUT OVERLAP.
    pRand=rand>0.9; %Needs to be 0.9 or larger.  0.8 made it VERY slow - too many.
    if pRand==1; pPhi(gp,1)=ppp; gp=gp+1; end %Sequentially activates elements.
    if gp==numGtoP+1; break; end; %ppp=1:numF*numPhi^2; end;
end

pspike=prod(tspike(pPhi,1:fullT));  %Multiplies together random set of grid cells.
plotPlaceX=xv(1,1:fullT).*pspike;
plotPlaceY=xv(2,1:fullT).*pspike;

if length(find(pspike))>3 && gp==numGtoP+1; %Only take place cells with at least three spikes from three grid inputs gp.

    stX=std(plotPlaceX(1,find(pspike)));
    stY=std(plotPlaceY(1,find(pspike)));
    
    if stX<stdlim && stY<stdlim && stX>0.1 && stY>0.1; %Selects only place cells with three grid inputs and >3 spikes
        
placeSpike(nP,:)=pspike;

%Keep track of times of spikes of selected place cells.
sumPlace=sumPlace+pspike;

WPG(nP,pPhi)=1; %For a successful place cell, adds connectivity to WPG for input from grids.
    nP=nP+1;
    
    end %End std loop
end %End find(pspike) loop
end %End nP loop.

nnnn
toc


if plotAllPlace==1;
    %Always plot the coverage of all place cells
figure('name','all place cells','Numbertitle','off','Position',[1000 400 200 200]); ax17=gca;
set(gcf,'DefaultAxesColorOrder',[gy gy gy])
axes(ax17);
plot(xv(1,1:fullT).*(sumPlace>0.5),xv(2,1:fullT).*(sumPlace>0.5),'r.'); %Plots the sum of all place cell firing to see coverage.
drawnow
end

%********ENCODING OF TRAJECTORY WITH PLACE, HD AND GRID************************
eeeT=fullT; %100; %fullT;
tspike2=zeros(numF*numPhi^2,eeeT);
placeSpike2=zeros(numP^2,eeeT);

%Grid cell phase arrays
ephis=zeros(length(h)/2,numF*numPhi^2); %,eeeT);
ephi=zeros(length(h),numF*numPhi^2); %,eeeT);
ephiOLD=zeros(length(h),numF*numPhi^2);
ephisOLD=zeros(length(h)/2,numF*numPhi^2);

useDXV=1;
%Cycle through frequencies and phases for creating grid cells.
for ff=1:numF; %GRID CELLS DURING ENCODING (USING SAME AS DURING CREATION.
    f=baseF+multFF*ff;  
    for pp=1:numPhi;
        for p2=1:numPhi;
        ix=[xv(1,1)+(pp+ff+pOff)*(sqrt(3)/(f*B*(numPhi))); xv(2,1)+(p2+ff+pOff)*(sqrt(3)/(f*B*(numPhi)))]; %f*G=H=sqrt(3)/2*B. G=sqrt(3)/2*f*B
        hx=((HHst*ix)>0).*(HHst*(ix));
        LH=length(h)/2; MH=length(h);
        if useDXV==1; ephis(:,1)=2*pi*f*dt; ephi(:,1)=2*pi*f*(dt+B*hx); end;

for t=2:eeeT;
if useDXV==1; ephis(:,t)=ephis(:,t-1)+2*pi*f*dt; ephi(:,t)=ephi(:,t-1)+2*pi*f*dt*(1+B*(HHst*dxv(:,t)>0).*(HHst*dxv(:,t))); 
elseif useDXV==0; ephis(:,t)=2*pi*f*t; ephi(:,t)=2*pi*f*(t+B*HHst*(xv(:,t)-ix));
end; %Creation of phase representation.
end

LH=length(h)/2; MH=length(h);
g=prod(cos((ephi(1:LH,:)-ephi(LH+1:MH,:))));

spiket=g>thresh; %Spiking of grid cell.

tspike2((ff-1)*numPhi^2+(pp-1)*numPhi+p2,1:fullT)=spiket; %Store spiking of all grid cells.

        end %End for p2=1:numPhi
    end %End for pp=1:numPhi
end %End for f=1:numF


for t=1:eeeT;

%PLACE CELLS DURING ENCODING
placeSpike2(:,t)=((WPG*tspike2(:,t))./numGtoP)>0.9; %Spiking of place cells during encoding.

speed(1,t)=sqrt(dxv(1,t)^2+dxv(2,t)^2);

%HD CELLS DURING ENCODING
hd(:,t)=(HHst*dxv(:,t)>0).*(HHst*dxv(:,t)); %Add the head direction cell activity during encoding.

hd(:,t);

    plothd(:,t)=hd(1:LH,t)-hd(LH+1:MH,t);

norm=sum(placeSpike2(:,t));
if norm~=0;
    for xx=1:numP^2;
        if placeSpike2(xx,t)==1;
WHP(:,xx)=(WHP(:,xx)+hd(:,t).*placeSpike2(xx,t))/2; %Create connections from place to head direction.
%WHP(:,xx)=(WHP(:,xx)+hd(:,t).*placeSpike2(xx,t)./norm)/2; %Create connections from place to head direction.
        end  %This procedure takes elements with presynaptic input and averages the existing and new weight.
    end
end

end %End for t=1:eeeT


if plotGrids==1
figure('name','grid during encoding','Numbertitle','off','Position',[250 50 numPhi*150 numPhi*150]); ax15b=gca;
axes(ax15b); %Grid cell figure

%Plot Grid cells during encoding.
for ff=1:numF;
    f=baseF+multFF*ff;
    for pp=1:numPhi;
        for p2=1:numPhi;
plotSpikeX=xv(1,1:eeeT).*tspike2((ff-1)*numPhi^2+(pp-1)*numPhi+p2,1:eeeT);
plotSpikeY=xv(2,1:eeeT).*tspike2((ff-1)*numPhi^2+(pp-1)*numPhi+p2,1:eeeT);
axes(ax15b);
%subplot(numF*numPhi,numPhi,(f-1)*numPhi^2+(pp-1)*numPhi+p2); ax15b=gca;
subplot(numPhi,numPhi,(pp-1)*numPhi+p2,'replace'); ax15b=gca;
plot(plotSpikeX,plotSpikeY,'r.');
xlim([xmin-1 xmax+1])
ylim([ymin-1 ymax+1])
set(gca,'XTickLabel',{' '})
set(gca,'YTickLabel',{' '})
drawnow
       end %End for p2=1:numPhi
    end %End for pp=1:numPhi
end %End for f=1:numF
end %End if plotGrids==1

if plotPlaceCells==1
figure('name','place cells Encoding','Numbertitle','off','Position',[800 50 numPlotPlace*70 numPlotPlace*70]); ax16bbb=gca;

%Plot Place Cells during encoding.
for pp=1:numPlotPlace^2;
plotPlaceX=xv(1,1:eeeT).*placeSpike2(pp,:);
plotPlaceY=xv(2,1:eeeT).*placeSpike2(pp,:);
axes(ax16bbb);
subplot(numPlotPlace,numPlotPlace,pp,'replace'); ax16bbb=gca;
plot(plotPlaceX,plotPlaceY,'r.');
xlim([xmin-1 xmax+1])
ylim([ymin-1 ymax+1])
set(gca,'XTickLabel',{' '})
set(gca,'YTickLabel',{' '})
end %End for pp
end %End if plotPlaceCells==1


%PLOT PLACE CELLS VS TIME DURING ENCODING spike times with colors.
axes(ax16b);
subplot(numPsizesTot*numPtestsTot,1,(numPsizes-1)*numPtestsTot+numPtests,'replace'); ax16b=gca;

numSING=numPlaceToPlot;
for sing=1:numSING;
    fps=find(placeSpike2(sing,:));
    for mm=1:length(fps);
        if colorSpikes==1;
            plot([fps(mm) fps(mm)],[sing sing-0.5],'Color',[rc(sing,:)],'LineWidth',1); hold on;
        elseif colorSpikes==0;
            plot([fps(mm) fps(mm)],[sing sing-0.5],'Color','k','LineWidth',1); hold on;
        end
    end
end
if colorSpikes==1;
    set(gca,'Color','k');
elseif colorSpikes==0
    set(gca,'Color','w');
end
xlim([1 eeeT]);
ylim([0 numSING+0.5]);
drawnow

if plotEnc==1;

%PLOT HD cells during encoding
axes(ax18);

for m=1:5:eeeT-1;   %Plot the head direction cell vectors. %If want to plot those during encoding start m=1:
     denom=cos(h(1,1))*sin(h(2,1))-sin(h(1,1))*cos(h(2,1));
     hdx=((plothd(1,m+1))*sin(h(2,1))-(plothd(2,m+1))*sin(h(1,1))); %./(denom*wd*B);
     hdy=((plothd(1,m+1))*-cos(h(2,1))+(plothd(2,m+1))*cos(h(1,1))); %./(denom*wd*B);
    
      AA=[xv(1,m) xv(1,m)+hdx];  
      BB=[xv(2,m) xv(2,m)+hdy];
      
      CC=[xv(1,m) xv(1,m)+dxv(1,m+1)*dt];
      DD=[xv(2,m) xv(2,m)+dxv(2,m+1)*dt];
      if sum(AA)~=0 && sum(BB)~=0;
      plot(AA,BB,'k-',CC,DD,'r-','LineWidth',1); %AA,BB,'k-',
      hold on
      end
end

%PLOT HD CELLS VERSUS TIME DURING ENCODING
axes(ax18b)

spikehd=zeros(length(h),length(hd));
for ht=1:length(hd); 
    maxhd=mean(max(hd));
    spikehd(:,ht)=rand(length(h),1)>(maxhd*ones(length(h),1)-hd(:,ht))/maxhd; 
end
numSING=length(h);
for sing=1:numSING;
    fps=find(spikehd(sing,:));
    for mm=1:length(fps);
%plot([fps(mm) fps(mm)],[sing sing-0.5],'Color',[rc(sing,:)],'LineWidth',1); hold on;
        if colorSpikes==1;
            plot([fps(mm) fps(mm)],[sing sing-0.5],'Color',[rc(sing,:)],'LineWidth',1); hold on;
        elseif colorSpikes==0;
            plot([fps(mm) fps(mm)],[sing sing-0.5],'Color','k','LineWidth',1); hold on;
        end
    end
end
if colorSpikes==1;
    set(gca,'Color','k');
elseif colorSpikes==0
    set(gca,'Color','w');
end
xlim([1 eeeT]);
ylim([0 numSING+0.5]);
drawnow

end %End plotEnc==1

toc

%********RETRIEVAL OF TRAJECTORY WITH PLACE, HD AND GRID************************
%Add RETRIEVAL period.

rrrT=eeeT; %eeeT;
hd3=zeros(length(h),rrrT);
hd3(:,1:5)=hd(:,5:9); %Set first step of retrieval hd to same as encoding hd.
tspike3=zeros(numF*numPhi^2,rrrT); 
tspike3(:,1:5)=tspike2(:,5:9); %Set initial grid cell spiking to be same as during encoding
%Grid cell spiking initial conditions don't matter. Can set to tspike2(:,21:22) and still function.  But does initial phase set by ix matter?
placeSpike3=zeros(numP^2,rrrT);
placeSpike3(:,1:5)=placeSpike2(:,5:9); %set initial place cell to be same as during encoding.

rphis=zeros(length(h)/2,numF*numPhi^2); 
rphi=zeros(length(h),numF*numPhi^2); 
plotrphi=zeros(length(h)/2,rrrT);
plotiG=zeros(length(h),1);

rphiOLD=zeros(length(h),numF*numPhi^2);
rphisOLD=zeros(length(h)/2,numF*numPhi^2);

if useDXV==1; %Set initial conditions for grid cells.
for ff=1:numF; 
    f=baseF+multFF*ff;
    for pp=1:numPhi;
        for p2=1:numPhi;
            cG=(ff-1)*numPhi^2+(pp-1)*numPhi+p2; %Number of current grid cell.
        ix=[xv(1,1)+(pp+ff+pOff)*ceil(sqrt(3)/(f*B*(numPhi))); xv(2,1)+(p2+ff+pOff)*ceil(sqrt(3)/(f*B*(numPhi)))]; %f*G=H=sqrt(3)/2*B. G=sqrt(3)/2*f*B
           iG=((HHst*ix)>0).*(HHst*ix);
           for nG=1:length(h);
             rphis(nG,cG)=2*pi*f*dt; rphi(nG,cG)=2*pi*f*(dt+B*iG(nG)); 
           end
           if cG==1
               plotiG=(HHst*(ix-xv(:,1))>0).*(HHst*(ix-xv(:,1)));
           end
           end %End for p2=1:numPhi
    end %End for pp=1:numPhi
end %End for f=1:numF
end %End if useDXV.

%RETRIEVAL ACTIVITY OF GRID CELLS, PLACE CELLS AND HD CELLS
for t=2:rrrT;
%GRID CELLS
for ff=1:numF; 
    f=baseF+multFF*ff;
    for pp=1:numPhi;
        for p2=1:numPhi;
            cG=(ff-1)*numPhi^2+(pp-1)*numPhi+p2; %Number of current grid cell.
        ix=[xv(1,1)+(pp+ff+pOff)*ceil(sqrt(3)/(f*B*(numPhi))); xv(2,1)+(p2+ff+pOff)*ceil(sqrt(3)/(f*B*(numPhi)))]; %f*G=H=sqrt(3)/2*B. G=sqrt(3)/2*f*B
    
if useDXV==1; 
    rphiOLD=rphi;
    rphisOLD=rphis; 
    rphis(:,cG)=rphisOLD(:,cG)+2*pi*f*dt; 
    rphi(:,cG)=rphiOLD(:,cG)+2*pi*f*dt*(1+B*strengHD*hd3(:,t));  %Use head direction vector hd3 retrieved from WHP.

elseif useDXV==0; rphis(:,cG)=2*pi*f*t; rphi(:,cG)=2*pi*f*(t+B*hd3(:,t)); 
end

LH=length(h)/2; MH=length(h);
tspike3(cG,t)=prod(cos(rphi(1:LH,cG)-rphi(LH+1:MH,cG)))>thresh; 

if cG==1;
    plotrphi(:,t)=rphi(1:LH,cG)-rphi(LH+1:MH,cG);
    plothd3(:,t)=hd3(1:LH,t)-hd3(LH+1:MH,t);
end
       end %End for p2=1:numPhi
    end %End for pp=1:numPhi
end %End for f=1:numF

%PLACE CELLS DURING RETRIEVAL
placeSpike3(:,t)=((WPG*tspike3(:,t))./numGtoP)>0.9; %Create place cell activity during encoding.
%HD CELLS
if max(placeSpike3(:,t))==0
    hd3(:,t+1)=hd3(:,t);
else 
    norm=sum(placeSpike3(:,t));
hd3(:,t+1)=WHP*placeSpike3(:,t)/norm;
end

end %End for t=1:rrrT

%PLOT INVERSE TRANSFORM TO SEE WHERE RAT "THINKS" IT IS LOCATED.
axes(ax14); %Go back to x y plot and plot inverse.
%subplot(numPsizesTot, numPtestsTot, (numPsizes-1)*numPtestsTot+numPtests); ax14=gca; %This does rows of tests.
subplot(numPsizesTot*numPtestsTot,1,(numPsizes-1)*numPtestsTot+numPtests,'replace'); ax14=gca; %This does one column of all.

%PLOT ACTUAL TRAJECTORY IN GRAY
set(gcf,'DefaultAxesColorOrder',[gy gy gy])

plot(pos_x(minT:minT+fullT), pos_y(minT:minT+fullT));
hold on

%PLOT RETRIEVED TRAJECTORY IN BLACK
set(gcf,'DefaultAxesColorOrder',[0 0 0])

 wd=(baseF+multFF*1)*2*pi;
 
 %Calculate inverse transform of grid cell phase (plotrphi(1,:) and plotrphi(2,:)
    denom=cos(h(1,1))*sin(h(2,1))-sin(h(1,1))*cos(h(2,1));
    xph=((plotrphi(1,:)-B*wd*plotiG(1,1))*sin(h(2,1))-(plotrphi(2,:)-B*wd*plotiG(2,1))*sin(h(1,1)))./(denom*wd*B);
    yph=((plotrphi(1,:)-B*wd*plotiG(1,1))*-cos(h(2,1))+(plotrphi(2,:)-B*wd*plotiG(2,1))*cos(h(1,1)))./(denom*wd*B);
    plot(xph(1,2:rrrT),yph(1,2:rrrT),'k','LineWidth',2); %xlim([xmin+30 xmax+30]); ylim([ymin-60 ymax-60]); 
xlim([xmin-1 xmax+1])
ylim([ymin-1 ymax+1])
drawnow

toc

%PLOT GRID CELLS DURING RETRIEVAL.
if plotGrids==1
figure('name','grid Retrieval','Numbertitle','off','Position',[250 50 numPhi*150 numF*150]); ax20=gca;

for ff=1:numF;
    f=baseF+multFF*ff;
    for pp=1:numPhi;
        for p2=1:numPhi;
plotSpikeX=xph(1,1:rrrT).*tspike3((ff-1)*numPhi^2+(pp-1)*numPhi+p2,1:rrrT);
plotSpikeY=yph(1,1:rrrT).*tspike3((ff-1)*numPhi^2+(pp-1)*numPhi+p2,1:rrrT);
axes(ax20);
subplot(numF,numPhi,(ff-1)*numPhi+pp,'replace'); ax20=gca;
plot(plotSpikeX,plotSpikeY,'r.');
xlim([xmin-1 xmax+1])
ylim([ymin-1 ymax+1])
drawnow
       end %End for p2=1:numPhi
    end %End for pp=1:numPhi
end %End for f=1:numF
end %End if plotGrids==1

%Plot Place Cells during retrieval.
if plotPlaceCells==1
figure('name','place cells Retrieval','Numbertitle','off','Position',[800 50 numPlotPlace*70 numPlotPlace*70]); ax21=gca;

for pp=1:numPlotPlace^2;
plotPlaceX=xph(1,1:rrrT).*placeSpike3(pp,:);
plotPlaceY=yph(1,1:rrrT).*placeSpike3(pp,:);
axes(ax21);
subplot(numPlotPlace,numPlotPlace,pp,'replace'); ax21=gca;
plot(plotPlaceX,plotPlaceY,'r.');
xlim([xmin-1 xmax+1])
ylim([ymin-1 ymax+1])
set(gca,'XTickLabel',{' '})
set(gca,'YTickLabel',{' '})
end %End for pp
end %End if plotPlaceCells==1

%PLOT PLACE VS TIME DURING RETRIEVAL.
axes(ax21b);
subplot(numPsizesTot*numPtestsTot,1,(numPsizes-1)*numPtestsTot+numPtests,'replace'); ax21b=gca;

numSING=numPlaceToPlot;
for sing=1:numSING;
    %for mm=1:length(timep(sing).spike);
    fps=find(placeSpike3(sing,:));
    for mm=1:length(fps);
        if colorSpikes==1;
            plot([fps(mm) fps(mm)],[sing sing-0.5],'Color',[rc(sing,:)],'LineWidth',1); hold on;
        elseif colorSpikes==0;
            plot([fps(mm) fps(mm)],[sing sing-0.5],'Color','k','LineWidth',1); hold on;
        end
        %plot([fps(mm) fps(mm)],[sing sing-0.5],'Color',[rc(sing,:)],'LineWidth',1); hold on;
        %plot([fps(mm) fps(mm)],[sing sing-0.5],'Color','k','LineWidth',1); hold on;
    end
end
if colorSpikes==1;
    set(gca,'Color','k');
elseif colorSpikes==0
    set(gca,'Color','w');
end
xlim([1 rrrT]);
ylim([0 numSING+0.5]);
drawnow

%PLOT HD CELLS VS SPACE during retrieval
%figure('name','HD cells Retrieval','Numbertitle','off','Position',[750 400 300 300]); ax19=gca;
axes(ax19);
subplot(numPsizesTot*numPtestsTot,1,(numPsizes-1)*numPtestsTot+numPtests,'replace'); ax19=gca; %This does one column of all.

plot(xph(1,2:rrrT),yph(1,2:rrrT),'r-','LineWidth',1); %xlim([xmin+30 xmax+30]); ylim([ymin-60 ymax-60]); 
hold on

for m=2:5:rrrT-1;   %Plot the head direction cell vectors. %If want to plot those during encoding start m=1:
     denom=cos(h(1,1))*sin(h(2,1))-sin(h(1,1))*cos(h(2,1));
     hdx=((plothd3(1,m+1))*sin(h(2,1))-(plothd3(2,m+1))*sin(h(1,1))); %./(denom*wd*B);
     hdy=((plothd3(1,m+1))*-cos(h(2,1))+(plothd3(2,m+1))*cos(h(1,1))); %./(denom*wd*B);
    
      %hdx=hd3(1,m+1);
      %hdy=-cot(h(2,1))*hd3(1,m+1) + hd3(2,m+1)/sin(h(2,1));  %Translates back out of head direction coordinates.
      AA=[xph(1,m) xph(1,m)+hdx];  
      BB=[yph(1,m) yph(1,m)+hdy];   

      if sum(AA)~=0 && sum(BB)~=0;
      plot(AA,BB,'k-','LineWidth',1); %AA,BB,'k-',
      xlim([xmin-10 xmax+10])
      ylim([ymin-10 ymax+10])
      hold on
      end
end
drawnow;

%PLOT HD CELLS VS TIME DURING RETRIEVAL
axes(ax19b); %Plot HD retrieval versus time.
subplot(numPsizesTot*numPtestsTot,1,(numPsizes-1)*numPtestsTot+numPtests,'replace'); ax19b=gca;

spikehd3=zeros(length(h),length(hd3));
for ht=1:5:length(hd3); 
    maxhd3=mean(max(hd3)); %Create a spiking version of the HD activity in hd3.
    spikehd3(:,ht)=rand(length(h),1)>(maxhd3*ones(length(h),1)-hd3(:,ht))/maxhd3; 
    %Generate a spike when rand is greater than 1-hd3.
end
numSING=length(h);
for sing=1:numSING;
    fps=find(spikehd3(sing,:));
    for mm=1:length(fps);
%plot([fps(mm) fps(mm)],[sing sing-0.5],'Color','k','LineWidth',1); hold on;
        if colorSpikes==1;
            plot([fps(mm) fps(mm)],[sing sing-0.5],'Color',[rc(sing,:)],'LineWidth',1); hold on;
        elseif colorSpikes==0;
            plot([fps(mm) fps(mm)],[sing sing-0.5],'Color','k','LineWidth',1); hold on;
        end
    end
end
if colorSpikes==1;
    set(gca,'Color','k');
elseif colorSpikes==0
    set(gca,'Color','w');
end
xlim([1 rrrT]);
ylim([0 numSING+0.5]);
drawnow

%PLOT X POS AND RUNNING SPEED FROM INVERSE TRANSFORM OF PHASE DURING
%RETRIEVAL
rspeed=zeros(1,rrrT);
for rt=3:rrrT;
    rspeed(1,rt)=sqrt((xph(rt)-xph(rt-1))^2 + (yph(rt)-yph(rt-1))^2)/dt;
end
axes(ax112);
tt=[1:rrrT];
plot(tt,xph(1,1:rrrT),tt,rspeed);
xlim([1 rrrT]);
ylim([-55 55]);

    halfrun=0
    numCircRun(numPsizes,1)=numP;
    for xyt=1:rrrT-1
        if yph(xyt)<-5 && yph(xyt+1)>-5 && xph(xyt)<-10 && halfrun==0;
            numCircRun(numPsizes,numPtests+1)=numCircRun(numPsizes,numPtests+1)+0.5;
            halfrun=1;
        end
        if yph(xyt)>5 && yph(xyt+1)<5 && xph(xyt)>10 && halfrun==1;
            numCircRun(numPsizes,numPtests+1)=numCircRun(numPsizes,numPtests+1)+0.5;
            halfrun=0;
        end 
    end
    
    strengHD
    numCircRun
    
        end %End while numCircRun<2.0
    end %End numPsizes
end %End numPtests
    
fid=fopen('numCircRunAA','a');
fprintf(fid, '%g\t',numCircRun);
fprintf(fid, '\n');

fclose(fid);

%END OF SCRIPT