%U5bzgNN2005EquationNumsForPosting.m - 
%This version includes Equation numbers matching those used in Hasselmo and
%Eichenbaum, 2005 and can create all but one of the figures from that paper.
%In addition, it can save a movie of spatial alternation behavior with a brown oval rat with a tail.
%Graphics of the virtual rat are plotted at intervals k=10 (k can be changed).  
%I also inserted a wSpltHipGr=1 option (with Splitter Hippocampal Activity Graphis) to select whether I
%want to save and show the hippocampal activity which when on slows things down a
%bit for large simulations.
%I ran a simulation of 7000 steps and got 102 correct trials and nError=0.

%FIGURES for Hasselmo and Eichenbaum (2005) Neural Networks are summarized
%in the section titled CHOOSE FIGURES TO SHOW.
%This summary addresses all simulations for that article except Figure 10, which was made with:
%CurvesaapTwoWayJustPlotsUsedForFIGURE6fewerColoredBACKUP.m - 
%From folder AHmatlabMATLABscriptsThetaTheoryPaperSep2004neurNet2005

%To get rid of previous windows, use "close all" - comment it out
%if you want to compare windows from successive simulations.
close all
clear all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%IMPORTANT - CHOOSE FIGURES TO SHOW - 
%Creates figs used in Hasselmo && Eichenbaum (2005) Neural Networks 18:1172-1190 
%showFastSpatAlt demonstrates a simple simulation of spatial alternation at low spatial resolution.
%This option creates some figures used in HE2005: Fig 3B uses Vaca1 (lower right). 
%Fig. 6A uses StateActionVal, Fig. 6B uses Q from Hip retrieval.
showFastSpatAlt=0;
%showSplitters plots splitter cells with higher resolution in both x and y direction
%This creates figs. used in HE2005 Fig. 8B1 left splitter, right splitter (middle)
%Fig. 8B2 uses left and right non-splitter. 
showSplitters=1; 
%showPrecess plots theta phase precession with in Skaggs plots (two cycles)
%with separate plots for first track run and later track runs.
%This creates HE2005 Fig. 9A2 (with Final Single Neuron Precession Plot (mid left).
%For sparse plot in Fig. 9A2 set Tsteps=85.
%Creates Fig. 9B1 and B2 in windows Color Final Single and Color First Pass.
%For denser color plots can set Tsteps=145 so get more passes.
showPrecess=0;  %Uses randomAfferent =0.6 and gives higher resolution in y direction.
%oldGraphics shows activity in many regions over time in old graphics
%format (used in Fig. 5 of HE2005).
oldGraphics=0;

%SPEED UP BY PLOTTING RAT LESS FREQUENTLY.
if showSplitters==1;
    k=10; %Interval for plotting rat. Very fast with k=20
    num=1; %Start number for interval which increments.
else
    k=1;
    num=1;
end
nCorrect=0;
nError=0;

%Set details of spatial alternation.
showSpatalt=0; showSpatalt3by5=0; showSpatalt5by5=0; showSpatalt3by7=0;
if showFastSpatAlt==1 || oldGraphics==1
showSpatalt=1; %Select size below.
showSpatalt3by5=1; showSpatalt5by5=0; showSpatalt3by7=0;
end

wGr=0;  %wGr=withGraphics - gives old style graphics.
if oldGraphics==1; wGr=1; end

%DEFINE DURATION OF TIME OF SIMULATION
%Tsteps is the number of step of movement.
if showSpatalt3by5==1  
Tsteps=150; forceperiod=40; 
%Forceperiod is the shaping time that behavior is "forced" in particular pathways.
elseif showSpatalt5by5==1 
    Tsteps=140; forceperiod=100;
elseif showPrecess==1
Tsteps=85; forceperiod=0;
elseif showSplitters==1
Tsteps=2000; forceperiod=800;
end

if oldGraphics==1; Tsteps=60; forceperiod=40; end

%DEFINE TASK BEING PERFORMED
openfield=0;
%Circular track used for plotting theta phase precession
circulartrack=0;
if showPrecess==1
    circulartrack=1;
end
%Spatial Alternation - NOTE must NOT have openfield=1 with this.
spatalt=0;
%Higher resolution with longer stem.
tallspatalt=0;
if showSplitters==1 || showSpatalt==1
     tallspatalt=1;
end

%DEFINE STARTING SIZE OF ENVIRONMENT.
%The first option can create spatial alternation of any size.
if tallspatalt==1 %IMPORTANT: ysz must be ODD number.
    xsz=3; ysz=5; %ysz must be ODD number only
    if showSpatalt5by5==1
        xsz=5; ysz=5; end
    if showSpatalt3by7==1
        xsz=3; ysz=7; end

    %Set scale of location relative to state.
      lxsz=xsz*1;
      lysz=ysz*1;
 
    if showSplitters==1 %Automatically make it 5 by 5.
      xsz=5; ysz=5; lxsz=xsz*4; lysz=ysz*4; end
end

if spatalt==1 
xsz=3; ysz=5; end

if openfield==1
    xsz=3; ysz=4; lxsz=xsz*4; lysz=ysz*4; end

if circulartrack==1
    xsz=3; ysz=7;
    lxsz=xsz*1;  %Have same size for circular track!
    lysz=ysz*4; end

%DEFINE STARTING BARRIERS
%This section creates the maze environment with size xsz by ysz
%and then creates barriers with the Sp(sx,sy)=1 and draws thebarriers.
Sp=zeros(xsz,ysz);
Lp=zeros(lxsz,lysz);
%Create a separate array for representing barriers with B(sx,sy)=1.
%This was necessary to ensure that indexes could be tested with sx+1 etc.
B=zeros(xsz+1,ysz+1);
LB=zeros(lxsz+1,lysz+1);

%Set up barriers for different tasks within the overall space.
if spatalt==1; Sp(2,2)=1; B(2,2)=1; Sp(2,4)=1; B(2,4)=1; end

if tallspatalt==1 %Define starting state walls.
    for ct=2:xsz-1; 
    for ct2=2:(ysz-1)/2; Sp(ct,ct2)=1; B(ct,ct2)=1; end
    for ct2=((ysz-1)/2)+2:(ysz-1); Sp(ct,ct2)=1; B(ct,ct2)=1; end; end; end

if tallspatalt==1 %Define starting location walls.
    for ct=lxsz/xsz+1:lxsz-lxsz/xsz; 
    for ct2=lysz/ysz+1:((lysz)/2)-(lysz/ysz)/2; Lp(ct,ct2)=1; LB(ct,ct2)=1; end
    for ct2=((lysz)/2)+(lysz/ysz)/2+1:(lysz-lysz/ysz); Lp(ct,ct2)=1; LB(ct,ct2)=1; end; end; end

%Circular track starting walls
if circulartrack==1; for ct=2:xsz-1; for ct2=2:ysz-1; Sp(ct,ct2)=1; B(ct,ct2)=1; end; end; end

if circulartrack==1;
    for ct=(lxsz/xsz)+1:lxsz-(lxsz/xsz);
        for ct2=(lysz/ysz)+1:lysz-(lysz/ysz);
            Lp(ct,ct2)=1; LB(ct,ct2)=1; end; end; end


%DEFINE START LOCATION AND REWARD LOCATION.
%Set the starting location of the virtual rat.
if openfield==1 || circulartrack==1
    stlx=lxsz/xsz+1; stly=1; stsx=2; stsy=1;
%Set the goal location in the field.
    glx=2; gly=lysz; gsx=2; gsy=ysz; end

if spatalt==1; stsx=1; stsy=3; gsx=3; glft=1; glrt=5; wntlft=0; wntrt=0; glg=1; grg=1; end

if tallspatalt==1
    stsx=1; stsy=(ysz-1)/2+1; stlx=1; stly=(lysz/ysz)*(ysz-1)/2+1;
    gsx=xsz; glft=1; glrt=ysz; wntlft=0; wntrt=0; glg=1; grg=1; end
    
%The following prevents a reward being given for the wrong location with
%readyLev=-2 which is the start value in spatial alternation.
if openfield==1  || circulartrack==1; readyLev=1; end
if spatalt==1 || tallspatalt==1; readyLev=-1; end

%The following define the appearance of reward (cheese - yellow circle) in
%movie.
 kpRew=0;
 kpMx=10;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%SET STARTING LOCATION AND STATE
%The current location/state of the "rat" is defined by sx and sy (set to starting location).
sx=stsx; sy=stsy; lx=stlx; ly=stly;
%The previous location/state of the "rat" is defined by osx and osy.
olx=stlx; oly=stly; osx=stsx; osy=stsy; oosx=osx; oosy=osy;

%DEFINE SIZE OF STATE VECTORS ON BASIS OF ENVIRONMENT AND ACTIONS.
%Number of output/action states No is 4 for up, down, left, right.
    No=4;

%Hip is the memory of previous output/action at a particular state, so has
%the same size as the number of output states.
    Hip=xsz*ysz+1;

%Total number "n" of states and "nloc" number of locations. 
    n=xsz*ysz+1;
    nloc=lxsz*lysz+1;

%DEFINE VALUE AND ACTION-VALUE VECTORS
%Network uses action values to compute decisions.
%V=value function (value of each state - number of states=xsz*ysz)
%Q=action-value function (value of each action at each state = xsz*ysz*No.
V=zeros(1,n);  %Value function for each state - total states=xsz*ysz
dV=zeros(1,n);
Q=zeros(1,n*No);  %Action value function - total number=xsz*ysz*No
Qhip=Q; Output=3; oldOutput=1; Rew=0; rewAvail=0; alpha=0.01; memalpha=0.01; 
oldgate=zeros(1,n);
gate=zeros(1,n); %Gate for indexing retrieved memory sx,sy
    
%SELECTION OF GRAPHICS AND CURRENT LOCATION PLOT
%IMPORTANT - wGr is now regulated by SuppFig1==1 at start of script.
strengCurrBehav=0.0; 
wSpltGr=0; 
if showPrecess==1 || showSplitters==1 || showSpatalt==1
    strengCurrBehav=0.5; %This indicates the strength of the gray scale plotting of current
%behavioral state in the plot of retrieval activity in Vaca1.
    wSpltGr=1; %Turns on precession and splitter graphics for the cell chosen 
%with coordinates pcsx, pcsy
    wSpltHipGr=0; %Regulates saving and presentation of hippocampal activity plots at end.
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%
%DEFINE PARAMETERS OF HIPPOCAMPAL FUNCTION

%DETERMINE USE OF DENTATE GYRUS - CURRENT CONTEXT OR MATCHING CONTEXT.
%This parameter determines the use of DG.  For UseDG=0, it just sets adg=step.
%For UseDG=1 it uses DG to grab a more local temporal context to do episodic retrieval.
UseDG=0;

%INCLUDE AFFERENT IN aecIII
%This includes the afferent input in the computation of spreading activity
%in aecIII.  This strongly drives ca1 due to convergence with current
%context from ca3, unless I set UseDG=1.
includeAfferent=1;  %If set to 1 get current context activity in splitters and precession.  
%If zero don't get current context in splitters (but also lose flt part of precession).
%IMPORTANT - If want to avoid current sensory input driving splitters, use includeAfferent=0.

%DEFINE RETRIEVAL SPREAD.
%Number of step of retrieval for computation of Vaca1 at each location.
rtSprMx=3;  %NOTE spatalt 3 by 5 FAILS if rtSprMx= or > 5, only works <=4.
%Random afferent input timing for precession plots (0.6)
randAfferent=0.0;     %Set to 0.0 if want deterministic input causing uniform forward spread.
if showPrecess==1; randAfferent=0.6; end   %Set to 0.6 for Fig 4, 
%Background activity in CA3 must be smaller than mu^fullpath = mu^lysz*2=0.04^36   
CA3background=0; %
if showPrecess==1; CA3background=1.0e-120; rtSprMx=4; end %Only use this background for showing the Mehta effect.          
%tao=number of steps of spreading computed in entorhinal cortex (usually
%activity will fall below threshold before this limit is reached).
%tao=tao in the NN2005 equations.
tao=3*rtSprMx; %tao=12
%Add a flt period for encoding with thetaEC=eta and thetaCA3=mu.
flt=1;

%DEFINE PARAMETERS OF HIPPOCAMPAL FUNCTION
        eta=0.04; %Regulates forward spread in ECIII.
        epsilon=0.0001;
        threshEC=eta-epsilon;
        mu=0.01;  %Regulates context spread in CA3.
        inhibCA1=0.25; %Reduced from 0.25 - 0.15 makes massive errors.
        
%DEFINE HIPPOCAMPAL VECTORS AND MATRICES 
    maxt=(rtSprMx+flt)*tao;
    a=zeros(1,n);
    olda=zeros(1,n);
    aecII=zeros(1,n);
    aecIII=zeros(1,n);  %Specific activity in aecIII at a certain phase of theta.
    aecIIIfull=zeros(1,n); %Activity in aecIII after full spread.
    adg=zeros(1,Tsteps);
    aca3=zeros(1,n);
    aca1=zeros(1,n);
    WecIII=zeros(n,n);
    Wdg=zeros(Tsteps,n);
    Wca3=zeros(n,Tsteps);
    VaecIII=zeros(n,Tsteps*maxt);
    Vaca1=zeros(n,(Tsteps+1)*maxt);
    if wGr~=1; Vaca3=zeros(n,(Tsteps+1)*maxt); end
    inhibSumCA1=0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%INITIAL DEFINITION OF GRAPHICS/FIGURES OF NETWORK FUNCTION.

%CREATE TICK MARKS.
for mmt=1:n; yticks(1,mmt)=mmt*3+0.5; end
for mmt2=1:n; xticks(1,mmt2)=mmt2*3+0.5;end

%STATE AND LOCATION PLOTS   
%Plot state as perceived by agent (lower resolution sx, sy coordinates)
xL=64*(1-0.4*Sp); figure('Position',[20 560 140 140]);
figh=gcf; axh=gca; colormap(gray);
set(gca,'Drawmode','fast'); image(xL)
set(figh,'Name','State (sx,sy)','numbertitle','off');

%Plot of exact location of agent in lx, ly coordinates.
xL=64*(1-0.4*Lp); figure('Position',[170 560 140 140]);
figH2=gcf; axH2=gca; colormap(gray);
set(gca,'Drawmode','fast'); set(figH2,'Name','Location (lx,ly)','numbertitle','off');

%ACTIVITY IN ECIII AND CA1
if wSpltHipGr==1 %This conditional selects hippocampus activity graphics for end.
%Graphics for ECIII activity (far right side of screen)
xVaecIII=64*(1-VaecIII);
figure('Position',[950 640 400 200]);
figVaecIII=gcf; axVaecIII=gca; colormap(gray);
set(gca,'Drawmode','fast'); image(xVaecIII)
set(figVaecIII,'Name','VaecIII','numbertitle','off');

%Graphics for CA1 spiking activity (far right side of screen)
xVaca1=64*(1-Vaca1);
figure('Position',[950 340 400 200]);
figVaca1=gcf; axVaca1=gca; colormap(gray);
set(gca,'Drawmode','fast'); image(xVaca1)
set(figVaca1,'Name','Vaca1','numbertitle','off');

%Graphics for subthreshold vector aca3=aecII
xVaca3=64*(1-aca3);
figure('Position',[950 40 400 200]);
figaca3=gcf; axVaca3=gca; colormap(gray);
set(gca,'Drawmode','fast'); image(xVaca3)
set(figaca3,'Name','aca3','numbertitle','off');
end %End conditional of wSpltHipGr==1

%VALUE FUNCTION GRAPHICS
%Graphics of value function (lower left side of screen).
xV=Sp; xV=64*(1-xV);
figure('Position',[20 40 140 120]);
figk=gcf; axk=gca; colormap(gray);
set(gca,'Drawmode','fast'); image(xV)
set(figk,'Name','StateVal(V) or Wib','numbertitle','off');

%Graphics for action value functions (middle left side of screen).
xQ=zeros(xsz*3,ysz*3); xQ=64*(1-xQ);
figure('Position',[20 280 140 120]);
figj=gcf; axj=gca; colormap(gray);
set(gca,'Drawmode','fast'); image(xQ)
set(figj,'Name','StateActionVal(Q) or Wo','numbertitle','off');
  set(axj,'XTick',xticks, 'YTick',yticks);
  set(axj,'XGrid','on', 'YGrid','on','ZGrid','off','GridLineStyle','-');

%Graphics for memory action value function (middle left side of screen)
xQ2=zeros(xsz*3,ysz*3);
xQ2=64*(1-xQ2);
figure('Position',[170 280 140 120]);
figg=gcf; axg=gca; colormap(gray);
set(gca,'Drawmode','fast'); image(xQ2)
set(figg,'Name','Q from Hip retrieval','numbertitle','off');
  set(axg,'XTick',xticks, 'YTick',yticks);
  set(axg,'XGrid','on', 'YGrid','on','ZGrid','off','GridLineStyle','-');

%PRECESSION AND SPLITTER CELL GRAPHICS
if wSpltGr==1;
    
    Pgate=zeros(1,n); %Separate gate for data for precessPlace) 
    plchp=zeros(1,nloc); %Place cell plot datqa
    plchpl=zeros(1,nloc); %Splitter cell plot data
    plchpr=zeros(1,nloc);
    plchpln=zeros(1,nloc);
    plchprn=zeros(1,nloc);
    prcplt=zeros(maxt+1,nloc); %Precession plot data

%Graphics for theta phase precession at end
xQ3Bend=zeros(2*maxt,nloc);
xQ3Bend=64*(1-xQ3Bend);
figure('Position',[320 720 300 100]);
figibend=gcf; axibend=gca; colormap(gray);
set(gca,'Drawmode','fast'); image(xQ3Bend)
set(figibend,'Name','Final Single neuron Precession Plot','numbertitle','off');

%Graphics for theta phase precession at first pass
xQ3Bfirst=zeros(2*maxt,nloc);
xQ3Bfirst=64*(1-xQ3Bfirst);
figure('Position',[320 520 300 100]);
figibf=gcf; axibf=gca; colormap(gray);
set(gca,'Drawmode','fast'); image(xQ3Bfirst)
set(figibf,'Name','first pass (non-)precession','numbertitle','off');

%Colored graphics for theta phase precession at end
xQ3BendB=zeros(2*maxt,nloc);
xQ3BendB=64*(1-xQ3BendB);
figure('Position',[320 320 300 100]);
figibendb=gcf; axibendb=gca; colormap(gray);
set(gca,'Drawmode','fast'); image(xQ3BendB)
set(figibendb,'Name','Color final Single neuron Precession Plot','numbertitle','off');

%Colored graphics for theta phase precession at first pass
xQ3BfirstB=zeros(2*maxt,nloc);
xQ3BfirstB=64*(1-xQ3BfirstB);
figure('Position',[320 120 300 100]);
figibfb=gcf; axibfb=gca; colormap(gray);
set(gca,'Drawmode','fast'); image(xQ3BfirstB)
set(figibfb,'Name','Color first pass (non-)precession','numbertitle','off');

%Graphics for placecell
xQ3E=64*(1-Lp);
figure('Position',[630 600 140 100]);
figie=gcf; axie=gca; colormap(gray);
set(gca,'Drawmode','fast'); image(xQ3E)
set(figie,'Name','place cell','numbertitle','off');

%Graphics for splitter left placecell
xQ3F=64*(1-Lp);
figure('Position',[630 400 140 100]);
figif=gcf; axif=gca; colormap(gray);
set(gca,'Drawmode','fast'); image(xQ3F)
set(figif,'Name','left splitter','numbertitle','off');

%Graphics for right splitter place cell
xQ3G=64*(1-Lp);
figure('Position',[780 400 140 100]);
figig=gcf; axig=gca; colormap(gray);
set(gca,'Drawmode','fast'); image(xQ3G)
set(figig,'Name','right splitter','numbertitle','off');

%Graphics for non-splitter left placecell
xQ3Fn=64*(1-Lp);
figure('Position',[630 200 140 100]);
figifn=gcf; axifn=gca; colormap(gray);
set(gca,'Drawmode','fast'); image(xQ3Fn)
set(figifn,'Name','left non-splitter cell','numbertitle','off');

%Graphics for right non-splitter place cell
xQ3Gn=64*(1-Lp);
figure('Position',[780 200 140 100]);
figign=gcf; axign=gca; colormap(gray);
set(gca,'Drawmode','fast'); image(xQ3Gn)
set(figign,'Name','right non-splitter cell','numbertitle','off');
end %END PRECESSION and SPLITTER graphics - wSpltGr conditional.

%PLOTS FOR OLD VERSION OF GRAPHICS (wGr=1)
if   wGr==1
%Define graphics for hippocampal output plots of aecII, aecIII by space.
%Graphics for aecII
xQ3D=zeros(1,n);
xQ3D=64*(1-0.4*Sp);
figure('Position',[170 560 160 160]);
figid=gcf; axid=gca; colormap(gray);
set(gca,'Drawmode','fast'); image(xQ3D)
set(figid,'Name','aecII by location (state)','numbertitle','off');

%Graphics for aecIII
xQ3H=zeros(1,n);
xQ3H=64*(1-0.4*Sp);
figure('Position',[340 560 160 160]);
figih=gcf; axih=gca; colormap(gray);
set(gca,'Drawmode','fast'); image(xQ3H)
set(figih,'Name','aecIII by location (state)','numbertitle','off');

%Graphics for aca1
xQ3K=zeros(1,n);
xQ3K=64*(1-0.4*Sp);
figure('Position',[510 560 160 160]);
figik=gcf; axik=gca; colormap(gray);
set(gca,'Drawmode','fast'); image(xQ3K)
set(figik,'Name','aca1 by location (state)','numbertitle','off');

%Graphics for WecIII
xW=64*(1-WecIII);
figure('Position',[680 560 160 160]);
figW=gcf; axW=gca; colormap(gray);
set(gca,'Drawmode','fast'); image(xW)
set(figW,'Name','WecIII','numbertitle','off');

%Graphics for Wca3
xWca3=64*(1-Wca3);
figure('Position',[850 560 160 160]);
figWca3=gcf; axWca3=gca; colormap(gray);
set(gca,'Drawmode','fast'); image(xWca3)
set(figWca3,'Name','Wca3','numbertitle','off');

%Graphics for vector of aecIII
xVaecIII=64*(1-VaecIII);
figure('Position',[170 220 160 200]);
figVaecIII=gcf; axVaecIII=gca; colormap(gray);
set(gca,'Drawmode','fast'); image(xVaecIII)
set(figVaecIII,'Name','VaecIII','numbertitle','off');

%Graphics for subthreshold vector aca1
xaca1=64*(1-aca1);
figure('Position',[340 220 160 200]);
figaca1=gcf; axaca1=gca; colormap(gray);
set(gca,'Drawmode','fast'); image(xaca1)
set(figaca1,'Name','aca1','numbertitle','off');

%Graphics for vector of spiking activity in ca1= Vaca1
xVaca1=64*(1-Vaca1);
figure('Position',[510 220 160 200]);
figVaca1=gcf; axVaca1=gca; colormap(gray);
set(gca,'Drawmode','fast'); image(xVaca1)
set(figVaca1,'Name','Vaca1','numbertitle','off');

%Graphics for subthreshold vector aca3=aecII
xVaca3=64*(1-aca3);
figure('Position',[680 220 160 200]);
figaca3=gcf; axVaca3=gca; colormap(gray);
set(gca,'Drawmode','fast'); image(xVaca3)
set(figaca3,'Name','aca3','numbertitle','off');

%Graphics for adg activity
xadg=64*(1-adg);
figure('Position',[850 220 160 200]);
figadg=gcf; axadg=gca; colormap(gray);
set(gca,'Drawmode','fast'); image(xadg)
set(figadg,'Name','adg','numbertitle','off');

end %END PLOTS FOR OLD GRAPHICS wGr==1 conditional


%This part of the script runs the "rat" through the environment.
%START MAIN PHASE/LOOP  -- This is the Main loop.
for step = 1:Tsteps;

%LOGIC FOR COMPUTING REWARD PRESENTATION.
%Compute finding of reward if state=goal location.
if (spatalt==1 || tallspatalt==1)
   
if sx==gsx && sy==glft && grg==1 && (readyLev==1 || readyLev==-1);
    rewAvail=1;goalfound=1;grg=0;glg=1;readyLev=0; Correct=1; nCorrect=nCorrect+1
    step
elseif sx==gsx && sy==glrt && glg==1 && (readyLev==1 || readyLev==-1); 
    rewAvail=1;goalfound=1;grg=1;glg=0;readyLev=0; Correct=1; nCorrect=nCorrect+1
    step
elseif sx==gsx && sy==glft && glg==1 && readyLev==1
    readyLev=0; Error=1; nError=nError+1
    step
elseif sx==gsx && sy==glrt && grg==1 && readyLev==1
    readyLev=0; Error=1; nError=nError+1
    step
elseif sx==stsx && sy==stsy && readyLev==0
    readyLev=1;
end
end %End spatalt,tallspatalt conditional

if(openfield==1 && sx==gsx && sy==gsy && (readyLev==1));
    rewAvail=1;goalfound=1;lx=stlx;ly=stly;sx=stsx;sy=stsy;
end

%SETTING BLOCKS TO GUIDE BEHAVIOR.
%Put block so rat goes counterclockwise, so that same place cell
%coordinates for spat alt splitter cell can be used in circulartrack to get
%precession going the correct direction to match Skaggs plots.
if circulartrack==1
    if sx==stsx && sy==stsy && step<2*(lxsz/xsz)*(lysz/ysz)
       Sp(1,1)=1; B(1,1)=1;
       Lp(lxsz/xsz,1:lysz/ysz)=1;
       LB(lxsz/xsz,1:lysz/ysz)=1;
   elseif step > 2
       Sp(1,1)=0; B(1,1)=0;
       Lp(lxsz/xsz,1:lysz/ysz)=0;
       LB(lxsz/xsz,1:lysz/ysz)=0;
   end
end

%SET BLOCKS TO GUIDE MOVEMENT DOWN STEM IN SPATALT
if (spatalt==1 && step < forceperiod) || (spatalt==1);
   if sx==1 && sy==3 
       Sp(1,2)=1; B(1,2)=1; Sp(1,4)=1; B(1,4)=1;
   else
       Sp(1,2)=0;B(1,2)=0;Sp(1,4)=0;B(1,4)=0;
   end
   %SET BLOCKS TO GUIDE ALTERNATION
   if sx==3 && sy==2
       wntlft=1;wntrt=0;Sp(3,4)=0;B(3,4)=0;
   elseif sx==3 && sy==4
       wntlft=0;wntrt=1;Sp(3,2)=0;B(3,2)=0;
   end
   if wntlft==1 && sx==3 && sy==3
       Sp(3,2)=1; B(3,2)=1;
   elseif wntrt==1 && sx==3 && sy==3
       Sp(3,4)=1; B(3,4)=1;
   end
end

if (spatalt==1 && step >= forceperiod)
    Sp(1:xsz,1:ysz)=0;
    B(1:xsz,1:ysz)=0;
    Sp(2,2)=1;B(2,2)=1;
    Sp(2,4)=1;B(2,4)=1;
    if sx==1 && sy==3 
       Sp(1,2)=1;B(1,2)=1;
       Sp(1,4)=1;B(1,4)=1;
   end
end

%SET BLOCKS TO GUIDE MOVEMENT DOWN STEM IN SPATALT
if (tallspatalt==1) 
   if sx==stsx && sy==stsy
       Sp(1,(ysz-1)/2)=1;
       B(1,(ysz-1)/2)=1;
       Lp(1:lxsz/xsz,(lysz/ysz)*((ysz-1)/2))=1;
       LB(1:lxsz/xsz,(lysz/ysz)*(ysz-1)/2)=1;
       Sp(1,(ysz-1)/2+2)=1;
       B(1,(ysz-1)/2+2)=1;
       Lp(1:lxsz/xsz,(lysz/ysz)*((ysz-1)/2+1)+1)=1;
       LB(1:lxsz/xsz,(lysz/ysz)*((ysz-1)/2+1)+1)=1;
   else
       Sp(1,(ysz-1)/2)=0;
       B(1,(ysz-1)/2)=0;
       Lp(1:lxsz/xsz,(lysz/ysz)*((ysz-1)/2))=0;
       LB(1:lxsz/xsz,(lysz/ysz)*(ysz-1)/2)=0;
       Sp(1,(ysz-1)/2+2)=0;
       B(1,(ysz-1)/2+2)=0;
       Lp(1:lxsz/xsz,(lysz/ysz)*((ysz-1)/2+1)+1)=0;
       LB(1:lxsz/xsz,(lysz/ysz)*((ysz-1)/2+1)+1)=0;
   end
end
   %SET BLOCKS TO GUIDE ALTERNATION
if (tallspatalt==1 && step < forceperiod) 
    if sx==xsz && sy==(ysz-1)/2
       wntlft=1;wntrt=0;
       Sp(xsz,(ysz-1)/2+2)=0;B(xsz,(ysz-1)/2+2)=0;
       Lp(lxsz-(lxsz/xsz)+1:lxsz,(lysz/ysz)*((ysz-1)/2))=0;
       LB(lxsz-(lxsz/xsz)+1:lxsz,(lysz/ysz)*(ysz-1)/2)=0;
   elseif sx==xsz && sy==(ysz-1)/2+2
       wntlft=0; wntrt=1;
       Sp(xsz,(ysz-1)/2)=0;B(xsz,(ysz-1)/2)=0;
       Lp(lxsz-(lxsz/xsz)+1:lxsz,(lysz/ysz)*((ysz-1)/2+1)+1)=0;
       LB(lxsz-(lxsz/xsz)+1:lxsz,(lysz/ysz)*((ysz-1)/2+1)+1)=0;
   end
   if wntlft==1 && sx==xsz && sy==(ysz-1)/2
       Sp(xsz,(ysz-1)/2)=1;B(xsz,(ysz-1)/2)=1;
       Lp(lxsz-lxsz/xsz+1:lxsz,(lysz/ysz)*((ysz-1)/2+1)+1)=0;
       LB(lxsz-lxsz/xsz+1:lxsz,(lysz/ysz)*((ysz-1)/2+1)+1)=0;
       Lp(lxsz-lxsz/xsz+1:lxsz,(lysz/ysz)*((ysz-1)/2))=1;
       LB(lxsz-lxsz/xsz+1:lxsz,(lysz/ysz)*(ysz-1)/2)=1;
   elseif wntrt==1 && sx==xsz && sy==(ysz-1)/2+2
       Sp(xsz,(ysz-1)/2+2)=1;B(xsz,(ysz-1)/2+2)=1;
       Lp(lxsz-lxsz/xsz+1:lxsz,(lysz/ysz)*((ysz-1)/2))=0;
       LB(lxsz-lxsz/xsz+1:lxsz,(lysz/ysz)*(ysz-1)/2)=0;
       Lp(lxsz-lxsz/xsz+1:lxsz,(lysz/ysz)*((ysz-1)/2+1)+1)=1;
       LB(lxsz-lxsz/xsz+1:lxsz,(lysz/ysz)*((ysz-1)/2+1)+1)=1;
   end
end

if (tallspatalt==1 && step >= forceperiod)   
       Sp(xsz,(ysz-1)/2)=0;
       Lp(lxsz-lxsz/xsz+1:lxsz,(lysz/ysz)*((ysz-1)/2))=0;
       LB(lxsz-lxsz/xsz+1:lxsz,(lysz/ysz)*(ysz-1)/2)=0;
       Lp(lxsz-lxsz/xsz+1:lxsz,(lysz/ysz)*((ysz-1)/2+1)+1)=0;
       LB(lxsz-lxsz/xsz+1:lxsz,(lysz/ysz)*((ysz-1)/2+1)+1)=0;
    B(xsz,(ysz-1)/2)=0;
    Sp(xsz,(ysz-1)/2+2)=0;
    B(xsz,(ysz-1)/2+2)=0;
end

if (tallspatalt==1 && step >= forceperiod)
   if sx==1 && sy==stsy
       Sp(1,(ysz-1)/2)=1;B(1,(ysz-1)/2)=1;
       Sp(1,(ysz-1)/2+2)=1;B(1,(ysz-1)/2+2)=1;
   elseif (sx==1 && sy==1) || (sx==1 && sy==ysz)
       Sp(1,(ysz-1)/2)=0;B(1,(ysz-1)/2)=0;
       Sp(1,(ysz-1)/2+2)=0;B(1,(ysz-1)/2+2)=0;
   end
end

if openfield==1;
if (olx==glx && oly==gly); lx=stlx; ly=stly; olx=stlx; oly=stly; end
end

%SHOW LOCATION OF AGENT IN ENVIRONMENT.
%These graphics show the moving location (state) of the agent (rat) in the
%environment.
 xL=64*(1-0.4*Sp); %Insert state environment.
 xL(sx,sy)=0; %Zero appears black - shows agent in perceived state (sx,sy).

  if step==k*num; %step==1 || step>forceperiod; %Only show rat after forceperiod ends.

 xLp=64*(1-0.4*Lp); %Insert location environment.
 
 if showSplitters~=1;
    xLp(lx,ly)=0; %Zero appears black - shows location (lx,ly).
 end %End black rectangle for nonsplitters.
 axes(axH2); 
 image(xLp); 

 if showSplitters==1 ; %Show virtual rat for splitter task.
     
 if Output == 1 %These lines draw the brown rat (Facecolor) with tail (line)
 rectangle('position',[ly-0.8,lx-0.8,1.2,1.2+1.8],'curvature',[1,1],'Facecolor',[0.5,0.3,0.2],'EraseMode','none')
 line([ly ly-1.5],[lx+2.2 lx+3.7],'Color',[0.5,0.3,0.2],'EraseMode','none');
elseif Output == 2
 rectangle('position',[ly-0.8,lx-0.8,1.2,1.2+1.8],'curvature',[1,1],'Facecolor',[0.5,0.3,0.2],'EraseMode','none')
 line([ly ly-1.5],[lx lx-1.5],'Color',[0.5,0.3,0.2],'EraseMode','none');
elseif Output == 3
 rectangle('position',[ly-0.8,lx-0.8,1.2+1.8,1.2],'curvature',[1,1],'Facecolor',[0.5,0.3,0.2],'EraseMode','none')
 line([ly+2.2 ly+3.7],[lx lx-1.5],'Color',[0.5,0.3,0.2],'EraseMode','none');
elseif Output == 4
 rectangle('position',[ly-0.8,lx-0.8,1.2+1.8,1.2],'curvature',[1,1],'Facecolor',[0.5,0.3,0.2],'EraseMode','none')
 line([ly ly-1.5],[lx lx+1.5],'Color',[0.5,0.3,0.2],'EraseMode','none');
end

 %Draws cheese in this section.
 if Rew==1 || (kpRew>=1 && kpRew<=ceil(kpMx/k))
   rectangle('position',[ly-0.8, lx-0.8, 0.1+3*(kpMx-kpRew)/kpMx,0.1+3*(kpMx-kpRew)/kpMx],'curvature',[0.8,0.6],'Facecolor','y');
     kpRew=kpRew+1;
 elseif kpRew>ceil(kpMx/k); kpRew=0;
 end
 
 end %End of showSplitters==1 conditional.
 
 drawnow
 
  elseif step~=k*num && Rew==1;
      kpRew=kpRew+1;
 end %End step=k*num conditional.
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ENCODE AND RETRIEVE EPISODIC MEMORY IN HIPPOCAMPUS CIRCUIT
%Start Episodic memory CONDITIONAL loop.
   %ESSENTIAL COMPONENTS
   %Update "a" based on location sx, sy and put previous "a" into "olda".
   olda=a;

   %PUT PERCEIVED STATE INTO THE ACTIVITY VECTOR (based on actual location lx,ly) 
   a((osx-1)*ysz+(osy-1)+1)=0;
   a((sx-1)*ysz+(sy-1)+1)=1;
   a=double(a);
   
%UPDATE the "memory" of previous retrieval oldgate from the previous
%behavioral step before cycling through rtspr.  Will eventually
%probably want to buffer the whole retrieved sequence for its usefulness.
      oldgate=gate;

%EQUATION 1 - ENCODE FORWARD ASSOC IN WecIII - provides forward spread based on prior associations.
if a*olda'==0  %This conditional was added for case where multiple step lx,ly inside single state sx,sy.
    %This is a dot product which is zero only when the new activity "a" does not overlap the previous activity "olda"
    %When used, this conditional prevents self-connections (bs*bsT).  Makes activity less messy.
     WecIII=(a)'*(olda) + WecIII;  %Equation 1 - NN2005. Store forward associations on WecIII.

    %Maximum limit on growth of WecIII - prevents explosive spread.
    for wn=1:n; for wn2=1:n; if WecIII(wn,wn2) > 1; WecIII(wn,wn2) = 1; end; end; end  

end %End if a*olda'==0 conditional.

%EQUATION 3 - ENCODING PROCESSES.
   %TEMPORAL CONTEXT IN aecII - aecII provides decaying representation of activity - temporal context
   %(Use a(1:n-1) to leave out Rew which causes memory problems (theta off)
   aecII(1:n)=a(1:n) + mu*aecII(1:n); %Equation 3 in NN2005.

%EQUATION 4 - Wca3 (and Wdg create input to DG for each point of aecII activity)
   Wdg(step, 1:n)=aecII;
   %Add CA3 - create connection between units in CA3 activated by DG and the context input from ECII.
   Wca3(1:n, step)=aecII' + CA3background*ones(1,n)'; %Implements equation 4 in NN2005.
   
%COMPUTE DENTATE GYRUS ACTIVITY DUE TO Wdg*acII.
   %Activation of episodic representations in DG by current context
          if UseDG==1  %This retrieves context in DG then projects forward for temporal context.
              adg=(Wdg*aecII')';
              if step > 3; for ss=0:3; adg(step - ss)=0; end; end; %Zero out most recent association.
              [c,i]=max(adg); adg=zeros(1,Tsteps);
              %Following step chooses DG context which is the length of the retrieval cycle forward from the current input.
              if step > rtSprMx && i < step-(rtSprMx+flt)
                 adg(i+rtSprMx)=1.0; end
          elseif UseDG==0  %This version just uses current context, setting everything to zero and setting adg(step)=1
              adg=zeros(1,Tsteps); 
              if step > 3*(lysz/ysz*lxsz/xsz)  %This guards against cases where agents stays in start state sx,sy for two full traversals of the state.
                %Note that having this start from step-lysz/ysz was necessary for use of "a" in equation.
                adg(step-3*(lysz/ysz*lxsz/xsz))=1.0; end
          end
   
%THETA CYCLE - UPDATE HIPPOCAMPAL ACTIVITY 
for time = 1:maxt;
    
%EQUATION 2 - COMPUTE FORWARD SPREAD IN ECIII.
              if rand > randAfferent  %If randAfferent=0, then aff=1 the whole time.
                   aff=1; else; aff=0;  end
 
           if time > flt*tao 
               if includeAfferent==1
               aecIIIfull=(((eta^(1/((time-flt*tao)/tao))*WecIII))*aecIII')' + aff*a; %Equation 2.
               elseif includeAfferent==0
               aecIIIfull=(((eta^(1/((time-flt*tao)/tao))*WecIII))*(aecIII + aff*a)')'; end;   
           else; aecIIIfull=(eta^tao*WecIII*aecIII')' + aff*a; end
   
            for wn=1:n %Inclusion of threshold PsyEC in Equation 2.
               if aecIIIfull(wn) > threshEC; 
                   aecIII(wn)=aecIIIfull(wn)-threshEC; 
               else; aecIII(wn)=0; end; end

%EQUATION 5 - COMPUTE CA3 ACTIVITY.
          %Add oscillations of aca3 using the theta oscillations
          if time >= flt*tao
               thetaCA3=mu^((time-flt*tao)/tao); %Theta in Equation 5 (flt*tao=phi in equations).
          else; thetaCA3=mu^tao; end      
          %Multiply CA3 theta times retrieval from DG.
          aca3=thetaCA3*(Wca3*adg')'; %Equation 5 in NN2005.
          
          
%EQUATION 6 - COMPUTE CA1 ACTIVITY - Multiplication in CA1 
          oldaca1=aca1;
          aca1=aecIII.*aca3; %Subthreshold part of Equation 6 in NN2005.
          %COMPUTE INHIBITION - Set threshold on basis of multiplication in CA1 
          inhibSumCA1=0;
          inhibSumCA1=inhibCA1*sum(aca1);
          %COMPUTE CA1 OUTPUT - Compute thresholded output of CA1. 
          %Hold activity Vaca1 as a row vector for use in Learning below.
          Pgate=zeros(1,n);
          Pgate=(aca1 > inhibSumCA1); %Thresholded Equation 6.

          if time==maxt; %rtSprMx*tao-1   %Select end of retrieved sequence for TD learning.
              gate=zeros(1,n);
              gate=(aca1 > inhibSumCA1);
          end
          
 %CREATE VECTORS FOR OLDGRAPHICS
          %Create vectors for graphics and performance measure
          if wGr==1 
          VaecIII(1:n,(step-1)*maxt+time)=5*aecIII'; 
          Vadg(1:Tsteps,(step-1)*maxt+time)=adg';
          Vaca3(1:n,(step-1)*maxt+time)=aca3';
          subaca1(1:n,(step-1)*maxt+time)=aca1'; 
          end
          
 %CREATE VECTORS TO SHOW ECIII and CA1 ACTIVITIES AT END.
           if wSpltHipGr==1
           VaecIII(1:n,(step-1)*maxt+time)=5*aecIII'; 
           Vaca1(1:n,(step-1)*maxt+time)=(aca1' > inhibSumCA1) + strengCurrBehav*(a'); 
           Vaca3(1:n,(step-1)*maxt+time)=aca3';
           end
          
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%SECTION FOR PLACE CELL, SPLITTER AND PRECESSION GRAPHICS.
if wSpltGr==1 %Start wSpltGr conditional.

%SELECTION OF INDEX FOR SPECIFIC PLACE CELLS TO BE PLOTTED.
if spatalt==1
    %Splitter cell in start of left arm, sx, sy coordinates.
    pcsx=3; pcsy=2; %place sx, sy
    %Non-splitter cell - exactly at base of stem - choice point.
    nspltcsx=3; nspltcsy=3;
    xnsplt=ysz*(nspltcsx-1)+nspltcsy;
elseif circulartrack==1
    pcsx=3; %place sx
    pcsy=ysz-2; %place sy
elseif tallspatalt==1 
    %Splitter cell - just to right of bottom of stem.
    pcsx=xsz;
    pcsy=(ysz-1)/2+2;
    %Non-splitter cell - exactly at base of stem - choice point.
    nspltcsx=xsz;
    nspltcsy=(ysz-1)/2+1;
    xnsplt=ysz*(nspltcsx-1)+nspltcsy;
end

%SET INDICES FOR PLOTTING PLACE CELL AND SPLITTER CELL RESPONSES.
xplace=ysz*(pcsx-1)+pcsy; %xplace is index for place cell being plotted.
xs=ysz*(sx-1)+sy; %Index for current state of rat used when plotting cell activity.
xl=lysz*(lx-1)+ly; %Index for current position within states
sT=(xsz+ysz)/Tsteps; %Scale spike averaging for rate plots to the size of environment.
sTe=(xsz+ysz)/(Tsteps*tao); %Scale to full size.
LTe=(lxsz+lysz) /(Tsteps*tao); 
    
%PLACE CELL GRAPHICS - Plots how much a selected cell (xplace pcsx, pcsy) is
%active when rat in different locations.
%PLOT ONE PLACE
if step > 2*ysz;
     plchp(xl)=1*Pgate(xplace) + plchp(xl);
     xQ3E(lx,ly)=64*(1-2*(LTe)*plchp(xl)); % Can't add because pushes to white+ xQ3E(sx,sy);
end

%PLOT FIRST PASS (<2*lysz) PRECESSION 
if step < 2*lysz && Pgate(xplace) > 0;  %For tallspatalt, need to make > 4*lysz to avoid first cycle no precession.
    %For first pass in linear track for Mehta effect use step < 2*lysz
    prcplt(maxt+1-time, xl)=Pgate(xplace) + prcplt(maxt+1-time,xl);
      if 0.5*maxt+1-time > 0;
      xQ3Bfirst(0.5*maxt+1-time,xl)= 64*(1-10*(sT)*prcplt(maxt+1-time,xl));  %Not scaled by: (sT)*
      else
      xQ3Bfirst(2.5*maxt+1-time,xl)= 64*(1-10*(sT)*prcplt(maxt+1-time,xl));  %Not scaled by: (sT)*
      end
    xQ3Bfirst(1.5*(maxt)+1-time,xl)= 64*(1-10*(sT)*prcplt(maxt+1-time,xl));  %Not scaled by: (sT)*
%PLOT SECOND PASS (>2*lysz) PRECESSION AT END
elseif step > 2*lysz && Pgate(xplace) > 0;  
    prcplt(maxt+1-time, xl)=Pgate(xplace) + prcplt(maxt+1-time,xl);
      if 0.5*maxt+1-time > 0;
      xQ3Bend(0.5*maxt+1-time,xl)= 64*(1-10*(sT)*prcplt(maxt+1-time,xl));  %Not scaled by: (sT)*
      else
      xQ3Bend(2.5*maxt+1-time,xl)= 64*(1-10*(sT)*prcplt(maxt+1-time,xl));  %Not scaled by: (sT)*
      end
    xQ3Bend(1.5*(maxt)+1-time,xl)= 64*(1-10*(sT)*prcplt(maxt+1-time,xl));  %Not scaled by: (sT)*
end


%PLOT EXAMPLES OF SPLITTER AND NONSPLITTER CELLS
if (spatalt==1 || tallspatalt==1)
   if step > 2*ysz && grg==1; %Data to plot right splitter and nonsplitter
   plchpr(xl)=1*Pgate(xplace) + plchpr(xl);
   xQ3F(lx,ly)=64*(1-2*(LTe)*plchpr(xl));
   plchprn(xl)=1*Pgate(xnsplt) + plchprn(xl);
   xQ3Fn(lx,ly)=64*(1-2*(LTe)*plchprn(xl));
   elseif step > 2*ysz && glg==1; %Data to plot left splitter and nonsplitter
   plchpl(xl)=1*Pgate(xplace) + plchpl(xl);
   xQ3G(lx,ly)=64*(1-2*(LTe)*plchpl(xl));
   plchpln(xl)=1*Pgate(xnsplt) + plchpln(xl);
   xQ3Gn(lx,ly)=64*(1-2*(LTe)*plchpln(xl)); %
   end
%End spatalt conditional.
end

%AT END PLOT PLACE CELL, SPLITTER CELL AND NONSPLITTER GRAPHICS
%These graphics show the activity of a single neuron in region CA1
%cumulative for each trial based on the location of the virtual rat.
   if (step==Tsteps && time==maxt) %Already started wSpltGr conditional above.
   %End plot of precession of indexed cell based on data from the first trial.
   axes(axibf); image(xQ3Bfirst); drawnow
   
   %End plot of precession of indexed cell on later trials.
   axes(axibend); image(xQ3Bend); drawnow
   
   %End plot of place cell in top middle.
   axes(axie); image(xQ3E); drawnow
   
   %End plots left and right splitter in middle.
   axes(axif); image(xQ3F); drawnow
   axes(axig); image(xQ3G); drawnow
   
   %End plots left and right nonsplitter.
   axes(axifn); image(xQ3Fn); drawnow
   axes(axign); image(xQ3Gn); drawnow
   
   %PLOT COLORED AND GAUSSIAN SMOOTHED PLOTS OF PRECESSION
gaussian=1; Tmax=2*maxt;
if gaussian==1
  for tt=1:Tmax; for ss=1:lxsz*lysz
    if xQ3Bfirst(tt,ss)~=64 && tt>8 && ss > 8 && tt<(Tmax-8) && ss<((lxsz*lysz)-8)
      for ttt=tt-8:tt+8; for sss=ss-8:ss+8; if ttt~=tt || sss~=ss
        xQ3BfirstB(ttt,sss)= xQ3BfirstB(ttt,sss)-5*(exp(-(ttt-tt)^2/5)*exp(-(sss-ss)^2/5));
end; end; end; end; end; end; end;
xQ3BfirstB=(64-xQ3BfirstB);
          
gaussian=1; Tmax=2*maxt;
if gaussian==1
  for tt=1:Tmax; for ss=1:lxsz*lysz
    if xQ3Bend(tt,ss)~=64 && tt>8 && ss > 8 && tt<(Tmax-8) && ss<((lxsz*lysz)-8)
      for ttt=tt-8:tt+8; for sss=ss-8:ss+8; if ttt~=tt || sss~=ss
        xQ3BendB(ttt,sss)= xQ3BendB(ttt,sss)-5*(exp(-(ttt-tt)^2/5)*exp(-(sss-ss)^2/5));
end; end; end; end; end; end; end;
xQ3BendB=(64-xQ3BendB);

%End plot of precession with colors and Gaussian smoothing.
   axes(axibfb); colormap('default'); image(xQ3BfirstB); drawnow
   axes(axibendb); colormap('default'); image(xQ3BendB); drawnow
   
   end %End conditional that graphics only appear at step=Tsteps.
  
end %End wSpltGr conditional.   
  
end %End THETA CYCLE - end of time for time=1:maxt loop



%CREATE FINAL CA1 AND ECIII ACTIVITY PLOTS (Vaca1 and VaecIII activity versus time).
if step==Tsteps && wSpltHipGr==1
   axes(axVaecIII); xVaecIII=64*(1-VaecIII);
   image(xVaecIII); drawnow
   axes(axVaca1); xVaca1=64*(1-Vaca1); 
   image(xVaca1); drawnow
   axes(axVaca3); xVaca3=64*(1-Vaca3);
   image(xVaca3); drawnow
end

%SHOW OLD GRAPHICS FORMAT.  VIEW ALL ACTIVITY AS SIMULATION PROGRESSES
if wGr==1  
%These plots show activity of all neurons according to the location being encoded
%by each neuron.  (Not place cell maps, as they do not show the firing of
%single neurons relative to the location of the virtual rat when they fire.

   %Show the activity in ECII. Create the graphical file for aecII by location.
   for msy=1:ysz; for msx=1:xsz; xQ3D(msx,msy)=64*(1-1000*aecII(ysz*(msx-1)+msy)); end; end
   axes(axid); image(xQ3D); drawnow
   %Show the connectivity of WecIII
   xW=64*(1-WecIII); axes(axW); image(xW); drawnow
   %Show the connectivity of Wca3
   xWca3=64*(1-Wca3); axes(axWca3); image(xWca3); drawnow
   %Show the retrieval vector in aecIII
   xVaecIII=64*(1-0.2*VaecIII); axes(axVaecIII); image(xVaecIII); drawnow
   %Show the retrieval spiking vector in ca1=Vaca1
   xaca1=64*(1 - 1E+17*subaca1); axes(axaca1); image(xaca1); drawnow
   %Show the retrieval spiking vector in ca1=Vaca1
   xVaca1=64*(1-Vaca1); axes(axVaca1); image(xVaca1); drawnow
   %Show the retrieval subthreshold vector in ca3 = Vaca3=aecII
   xVaca3=64*(1-1000*Vaca3); axes(axVaca3); image(xVaca3); drawnow
   %Show the retrieval activity in DG
   xadg=64*(1-Vadg); axes(axadg); image(xadg); drawnow
%End wGr==1 conditional
end

%LOGIC FOR REWARD DELIVERY (INPUT TO LAST ELEMENT OF VECTOR "a")
   %This line states that if reward was not available OR correct output not
   %produced, then NO reward
if rewAvail==0; Rew=0; rewAvail=0;
%This line states that if reward was available and correct output produced, then reward is received.
elseif rewAvail==1; Rew=1; rewAvail=0; goalfound=0; 
end  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%CREATE VALUE FUNCTION FOR BEHAVIORAL STATE ELEMENTS.
%Define index for old state-action io and
io=(osx-1)*ysz*No +(osy-1)*No;
is=(osx-1)*ysz + osy;
%Spiking - change representation to match size of a - number of spiking
%neurons.
dV(is)=V((sx-1)*ysz+sy)+Rew-V(is);
V(is)=V(is)+alpha*dV(is);

%UPDATE THE ACTION VALUE FUNCTION FOR BEHAVIORAL STATE.
%Use TD error to update current Q value (updates action Output with state).
Q(io+Output)=Q(io+Output)+alpha*dV(is);

%UPDATE ACTION-VALUE FUNCTION FOR MEMORY STATES.
%Use the state index of the one retrieved memory state from
%the END of the retrieval sequence.  Give this index to Qhip to be updated using dV
%computed for the prior position.
for hh=1:Hip
  if oldgate(hh)==1;  
  Qhip((hh-1)*No+Output)=Qhip((hh-1)*No+Output)+memalpha*oldgate(hh)*dV(is); end
end

%GRAPH THE VALUE FUNCTION.
%These graphics show the value function V of each state (not the Q)
xV(osx,osy)=64*(1.0-5*V(is));
if (step==Tsteps)
  axes(axk); image(xV); drawnow; end

%GRAPH THE ACTION VALUE FUNCTION.
%Update xQ, the graphical representation of Q.
if Output==1; xQ(osx*3-2,osy*3-1)=64*(1-10*Q(io+Output));
elseif Output==2; xQ(osx*3,osy*3-1)=64*(1-10*Q(io+Output));
elseif Output==3; xQ(osx*3-1,osy*3-2)=64*(1-10*Q(io+Output));
elseif Output==4; xQ(osx*3-1,osy*3)=64*(1-10*Q(io+Output));
%End creating the Q action-value function output
end

%This command plots the action-value function Q at end
if step==Tsteps;
   axes(axj); image(xQ)
     set(axj,'XTick',xticks, 'YTick',yticks);
     set(axj,'XGrid','on', 'YGrid','on','ZGrid','off','GridLineStyle','-');
   drawnow
end

%This section plots Qhip in the format of the location grid with
%up,down,left,right segments for each location in state-action-value format.
for msx=1:xsz
   for msy=1:ysz
      im=(msx-1)*ysz*No+(msy-1)*No;
   xQ2(msx*3-2,msy*3-1)=64*(1-10*Qhip(im + 1));
   xQ2(msx*3,msy*3-1)=64*(1-10*Qhip(im + 2));
   xQ2(msx*3-1,msy*3-2)=64*(1-10*Qhip(im + 3));
   xQ2(msx*3-1,msy*3)=64*(1-10*Qhip(im + 4));
   end
end

%Plots the action-value function Q for the memory Hip retrieval at end.
if step==Tsteps
   axes(axg);
   image(xQ2 )
      set(axg,'XTick',xticks, 'YTick',yticks);
      set(axg,'XGrid','on', 'YGrid','on','ZGrid','off','GridLineStyle','-');
   drawnow
end

if step==k*num;
    num=num+1;
end

if step==Tsteps;
    nCorrect
    nError
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%COMPUTATION OF OUTPUT/ACTION.    
    z=sprand(1,1,1);
    rando=sprand(1,4,1);  %Creates a 4 element random vector with values between 0 and 1.
    
    memgate=zeros(1,No);
    for hhh=1:Hip
        if gate(hhh)==1; memgate=Qhip(((hhh-1)*No + 1):((hhh-1)*No + 4));
        end
    end

    o=0.0001*rando + Q((sx-1)*ysz*No+(sy-1)*No+1:(sx-1)*ysz*No+(sy-1)*No+No) + 10*memgate;

       
%Make sure that the output/action "o" won't push the rat into a barrier or
%off the edge of the field.  Use of oosx, oosy keeps it progressing forward
    if lx<=1 || LB(lx-1,ly)==1 || (sx-1==oosx && sy==oosy) || (lx-1==olx && ly==oly)
        o(1,1)=-100; end
    if LB(lx+1,ly)==1 || lx>=lxsz  || (sx+1==oosx && sy==oosy) || (lx+1==olx && ly==oly) 
        o(1,2)=-100; end
    if ly<=1  || LB(lx,ly-1)==1 || (sy-1==oosy && sx==oosx) || (ly-1==oly && lx==olx)
        o(1,3)=-100; end
    if LB(lx,ly+1)==1 || ly>=lysz  || (sy+1==oosy && sx==oosx) || (ly+1==oly && lx==olx)
        o(1,4)=-100; end
    
%CHOOSE THE MAXIMUM OUTPUT/ACTION IN THE OUTPUT VECTOR.
%The following line finds the max of the output vector for selection of next movement. 
   [C,I]=max(o);
 
%Before altering the sx and sy values on basis of output, save the osx and osy values for computing TDlearning.
   olx=lx; oly=ly; osx=sx; osy=sy;
   oldOutput=Output;
   
    %For the output "I" indicating component of output vector "o" which is
    %maximum, update the state of the network.
    if o(1,1)==-100 && o(1,2)==-100 && o(1,3)==-100 && o(1,4)==-100  
       %This line prevents situations where all options are removed and agent gets stuck.
        lx=olx; ly=oly; o=zeros(1,No);
   elseif I==1 
       lx=lx-1; o=zeros(1,No); o(1,1)=1; Output=1;
   elseif I==2 
       lx=lx+1; o=zeros(1,No); o(1,2)=1; Output=2;
   elseif I==3 
       ly=ly-1; o=zeros(1,No); o(1,3)=1; Output=3;
   elseif I==4 
       ly=ly+1; o=zeros(1,No); o(1,4)=1; Output=4;
    end

%UPDATE SX,SY FROM LX,LY (state cue sx,sy from location lx,ly)
for st=1:xsz; if st > (lx*xsz)/lxsz-0.1 && st <= lx*xsz/lxsz+0.9; sx=st; end; end
for st=1:ysz; if st > ly*ysz/lysz-0.1 && st <= ly*ysz/lysz+0.9; sy=st; end; end

%If the value of osx or osy was updated, then change oosx and oosy
if osx~=sx || osy~=sy; oosx=osx; oosy=osy; end

%End the whole MAIN (ENCODING) LOOP.
end

%END OF SCRIPT