%AYHasselmoJCogNeurPFCforPosting.m - Script for posting on the Hasselmo lab web site.
%Script of goal-directed behavior guided by neocortical circuit model.
%%y Michael E. Hasselmo, Department of Psychology Center for Memory and
%Brain, Boston University, 2 Cummington St., Boston, MA 02215.
%hasselmo@bu.edu, http://people.bu.edu/hasselmo/

%This MATLAB script demonstrates examples of goal-directed behavior by a
%virtual rat (represented by a black square) in an open field environment.
%In this example, the virtual rat is started on the left (sv=2, sh=1).  
%The agent must find a goal/reward location on the right side of the field
%sv=2, sh=svsize.  Each time it finds the goal, it is reset to the starting location.
%This script is a simplified replication of the script used for the
%publication on prefrontal cortical function presented in:
%Hasselmo, M.E. (2005) A model of prefrontal cortical mechanisms for goal
%directed behavior.  J. Cogn. Neurosci.
%This script shows same properties of learning tasks, and the same general action-values
%as those in the older scripts used for the J. Cogn. Neurosci. paper.  Note
%that I've scaled the action value image for the neocortex version by 0.3,
%so what previously showed up as black now has gray scale values.
%The script is divided up into separate functions for the
%Neocortex/prefrontal cortex model (Neocortex=1) and the TDlearning example
%(Neocortex=0).  They share the formation of environment, plotting of
%location and action values, and the computation of output.
%These functions have not necessarily been implemented in an optimal way
%within MATLAB.  I appreciate receiving revised versions of the script
%which incorporate more efficient techniques for implementation.

%TDlearning section (Neocortex==0).
%This section is from TDlearningExample.m - 
%It uses temporal difference learning (TDlearning) and the actor-critic
%process to build a value function (V) and an action value function (Q) with each transition, 
%so that it initially learns to value of actions at locations adjacent to the reward, and
%then gradually learns to follow the optimal pathway (usually) which is
%straight across the middle from left to right.  Greater exploration of
%the environment can be obtained by setting RandomExploreThreshold to lower values
%which will cause it to generate random outputs.
%For details of reinforcement learning see:
%Sutton, R.S. and Barto, A.G. (1998) Reinforcement Learning. MIT Press,
%Cambridge, MA.
%This book will give most of the background of this example implementation,
%which is basically just an implementation of reinforcement learning in a
%gridworld environment and does not contain the additional features of my neurobiological
%implementation of neocortical function or hippocampal episodic retrieval.

%***************************
%Clear the screen.
close all
clear

%***************************
%CHOICES CHOICES CHOICES CHOICES CHOICES
%SET-UP CHOICES.

%CHOOSE NEOCORTEX (Neocortex=1) VS. REINFORCEMENT LEARNING (Neocortex=0).
Neocortex=1;

%CHOOSE AMOUNT OF GRAPHICS TO SHOW
%withgraphics=1 runs slower, but shows value function and action-value function 
%along with the location of the agent during learning.  
%withgraphics=0 shows only location of the agent.
withgraphics=1; 

%DEFINE TASK
openfield=1;

%DEFINE DURATION OF TIME OF SIMULATION
%Tsteps is the number of time steps of movement.
Tsteps=300;


%*******************
%ENVIRONMENT ENVIRONMENT ENVIRONMENT ENVIRONMENT ENVIRONMENT
%DEFINE STRUCTURE OF ENVIRONMENT.
%Sets the size of the space being explored.  State coordinates are labeled
%with sv=vertical state coordinate, sh=horizontal state coordinate
svsize=3;
shsize=3;
%This section creates the maze environment with size svsize by shsize
%and then creates barriers with the Sp(sv,sh)=1 and draws the barriers.
Sp=zeros(svsize,shsize);
%Create a separate array for representing barriers with B(sv,sh)=1.
%This was necessary to ensure that indexes could be tested with sv+1 etc.
B=zeros(svsize+1,shsize+1);

%Use this format to set up barriers within the open field space.
%Sp(2,2)=1;
%B(2,2)=1;
%Sp(2,4)=1;
%B(2,4)=1;

%DEFINE START LOCATION AND REWARD LOCATION.
%Set the starting location of the virtual rat.
    startsv=2;
    startsh=1;
%Set the goal location in the field.
    goalsv=2;
    goalsh=shsize;
    
%If svsize==1, then linear track requires different goal location.
if svsize==1;  
    startsv=1;
    startsh=1;
    goalsv=1;
    goalsh=shsize;
end

%PUT AGENT INTO START LOCATION.
%The current location/state of the "rat" appears as a black square
%defined by sv and sh (set to starting location).
%sv=vertical position in image, sh=horizontal position in image.
sv=startsv;
sh=startsh;
%The previous location/state of the "rat" is defined by oldsv and oldsh.
oldsv=startsv;
oldsh=startsh;


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

%Total number "n" of states and actions = number of place representations plus number of
%outputs (for proprioceptive feedback) plus reward sensor (one unit).
n=svsize*shsize + No + 1;

%Number of connections between columns is n*n
N=n*n;



%****************
%GRAPHICS GRAPHICS GRAPHICS GRAPHICS GRAPHICS GRAPHICS
%INITIAL DEFINITION OF GRAPHICS/FIGURES OF NETWORK FUNCTION.
%The following is used to plot location.
xL=zeros(svsize,shsize);
%The following is used to plot the value function (V).
xC=zeros(svsize,shsize);
%The following will be used to plot the action value (Q) function.
xQ=zeros(svsize*3,shsize*3);

%Graphics showing location of agent in environment.
xL=64*(1-0.4*Sp);
figure('Position',[360 20 300 300], 'doublebuffer','on');
figh=gcf;
axh=gca;
colormap(gray);
image(xL, 'EraseMode','none')
set(figh,'Name','Location (state)','numbertitle','off');

%Graphics showing value function of each state.
xC=Sp;
xC=64*(1-xC);
figure('Position',[40 20 300 300], 'doublebuffer','on');
figk=gcf;
axk=gca;
colormap(gray);
set(gca,'Drawmode','fast')
image(xC, 'EraseMode','none')
set(figk,'Name','StateVal(V) or Wib','numbertitle','off');
   
%Graphics showing action-value function Q plotted by state/location
xQ=64*(1-xQ);
figure('Position',[360 420 300 300], 'doublebuffer','on');
figj=gcf;
axj=gca;
colormap(gray);
set(gca,'Drawmode','fast')
image(xQ, 'EraseMode','none')
set(figj,'Name','Plot of State-Action-Value (Q)','numbertitle','off');
   set(axj,'XTick',[3.5,6.5,9.5,12.5], 'YTick',[3.5,6.5,9.5,12.5]);
   set(axj,'XGrid','on', 'YGrid','on','ZGrid','off','GridLineStyle','-');

 
%***************

%SET NEOCORTEX VECTORS AND MATRICES.
%SET NEOCORTEX ACTIVITY VECTORS.
if Neocortex==1
%Set up the Neocortex learning vectors for the minicolumns are set up with length N=n*n.
%a is sensory input, go is output for goal spread, gi is input for goal spread.  
%co is output for current state, ci is input fro current state.
a=zeros(1,N);
olda=zeros(1,N);
gi=zeros(1,N);
go=zeros(1,N);
ci=zeros(1,N);
ci2=zeros(1,N);
co=zeros(1,N);
goR=zeros(1,N);  
%Both TD learn and Neocortex use output vector o.
o=zeros(1,No);
%The network will transition between action and state conditions.
%When action=1 the network input "a" contains the proprioceptive feedback but
%not state information.  When action=0 and state=1, the input vector "a"
%contains information about the sensory world.
action=0;
%Must have state=1 at start so that the first location will appearin "a".
state=1;
%Start Reward with default of zero.
Reward=0;

%SET NEOCORTEX GRAPHICS VECTORS FOR ACTIVITY.
%These vectors allow viewing of vectors over all time steps at end of
%simulation.
cA=zeros(1,N);
coldA=zeros(1,N);
cGO=zeros(1,N);
cGI=zeros(1,N);
cCO=zeros(1,N);
cQ=zeros(1,N);
cO=zeros(1,No);

%SET NEOCORTEX WEIGHT MATRICES.  
%The matrix of connections have size N*N.  Note that Wif are the n*n squares on the
%diagonal of the matrix, and the Wf are the other connections off that set of squares.
Wg=zeros(N,N);
Wig=zeros(N,N);
Wc=zeros(N,N);
Wic=zeros(N,N);
%The matrix of connections to the output units has size 2*N.
Wo=zeros(No,N);

%Create connectivity for converting single element of a or goR into a large scale vector
%which multiplies the single active element in the input A matrix into a block of n active units.  
Wa=zeros(n,N);
for m=0:n-1
    Wa(m+1,m*n+1:m*n+n)=1;
end

%INITIAL WEIGHT MATRIX STRENGTHS
%Create the selective Wg and Wc connectivity which is selective from unit i in minicolumn j to unit j in minicolumn i
for m2=1:n;
%m2 is the index of the minicolumn 
for m=1:n;
    %m2 indexes row minicolumn m2, in this minicolumn neuron m+1 receives from neuron m2+1 of minicolumn m
    %the row i is minicolumn m2 (which must be multiplied by the number of neurons n) + individual neuron m+1 
    %which receives from minicolumn m (multiplied by n) individual neuron m2+1
        Wg((m2-1)*n+m,(m-1)*n+m2)=0.5;
        Wc((m2-1)*n+m,(m-1)*n+m2)=0.5;
    end
end

%for m2=1:N;  %Does NOT work with all to all connectivity.
%    for m=1:N;
%        Wg(m,m2)=0.5;
%        Wc(m,m2)=0.5;
%        Wig(m,m2)=0.5;
%end
%end

%Create the Wig connectivity which is only internal to each column
for m2=0:n-1;
%m2 is the index of the minicolumn 
    %m2 indexes row minicolumn m2, and everything from m2*n to m2*n+n+1 should be filled in both row and column
    Wig(m2*n+1:m2*n+n,m2*n+1:m2*n+n)=0.5;
end

%End neocortex Neocortex==1 conditional.
end


%***************

%SET TD LEARNING VECTORS (LOOK-UP TABLES). 
%V=value function (value of each state - number of states=svsize*shsize)
%Q=action-value function (value of each action at each state = svsize*shsize*No.
if Neocortex==0;
V=zeros(1,n);
dV=zeros(1,n);
Q=zeros(1,n*No);
Reward=0;
o=zeros(1,No);
Output=3;
%Learning rate: alpha sets the magnitude of change with each computation of error.  Will
%learn much faster with alpha=0.2 than with smaller alpha (e.g. =0.01).
alpha=0.2;

%End Neocortex==0 conditional.
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN
%This part of the script runs the agent ("virtual rat") through the environment.
%START MAIN PHASE/LOOP  -- This is the Main loop.
for steps = 1:Tsteps

%LOGIC FOR COMPUTING REWARD PRESENTATION.

%Compute finding of reward if state=goal location.
if(sv==goalsv & sh==goalsh)
    Reward=1;
end
    
if openfield==1;
if (oldsv==goalsv & oldsh==goalsh);
    sv=startsv;
    sh=startsh;
    oldsv=startsv;
    oldsh=startsh;
    Reward=0;   %Commenting this out can allow local loops to form! 
    if Neocortex==1;
        a=zeros(1,N);
        o=zeros(1,No); 
        q=zeros(1,N);  
        %gi=zeros(1,N);
        %go=zeros(1,N);
        %ci=zeros(1,N);
        %co=zeros(1,N);
    end
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);
 xL(sv,sh)=0;
 %This command makes the new image on each step
 axes(axh);
 image(xL,'EraseMode','none')
drawnow



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%NEOCORTEX NEOCORTEX NEOCORTEX NEOCORTEX
%Neocortex CONDITIONAL
if Neocortex==1;

%COMPUTE STATE INPUT BASED ON CURRENT STATE
%This section updates the input vector "a" based on the current location if
%the state vector is 1 (and action=0).
if state==1;
    a((sv-1)*shsize*n+(sh-1)*n+1:(sv-1)*shsize*n+(sh-1)*n+n)=1;
    action=0; 
%End of state==1 conditional.
end

%COMPUTE THE INPUT "a" BASED ON PROPRIOCEPTION OF OUTPUT(ACTION)
%The following two lines provide input in activity y which corresponds to the sensory perception of
%the output produced by the network.  
if action==1;
a((svsize*shsize)*n+1:(svsize*shsize)*n+n)=o(1);
a((svsize*shsize+1)*n+1:(svsize*shsize+1)*n+n)=o(2);
a((svsize*shsize+2)*n+1:(svsize*shsize+2)*n+n)=o(3);
a((svsize*shsize+3)*n+1:(svsize*shsize+3)*n+n)=o(4);
state=1;
end

%COMPUTE INPUT "a" BASED ON REWARD (LAST COLUMN)
if(Reward==0 & action==0);
   a((svsize*shsize+No)*n+1:(svsize*shsize+No)*n+n)=0;
   Reward=0;
end
if(Reward==1 & action==0);
   a((svsize*shsize+No)*n+1:(svsize*shsize+No)*n+n)=1;
   Reward=0;
end
 
%COMPUTATIONS FOR NEOCORTEX (TDlean==0) OF ACTIVITY FOR ACTION-VALUES

%IMPORTANT: Update output back population "go" on basis of input "a"
go=(0.6*a + 0.6*go  >= 0.7);  %This is Equation E1b version.
%go=a;  %This equation would be used to replace E1b if using Equation E1 in paper.

%ENCODING LOOP -- Update back population "gi" with previous input
%"olda" and backward subthreshold activity Wg*go.
gi=(0.6*olda + 0.6*(Wg*go')' >= 0.7);  %Equation E2.
gi=double(gi);

%ENCODING LOOP -- Compute the change in backward connections Wg (from go to gi)
dWg=gi'*go;  %Equation E3.
for m2=0:n-1;
    for m=0:n-1;
        Wg(m2*n+m+1,m*n+m2+1)=Wg(m2*n+m+1,m*n+m2+1)+dWg(m2*n+m+1,m*n+m2+1);
        if Wg(m2*n+m+1,m*n+m2+1)>1.0
            Wg(m2*n+m+1,m*n+m2+1)=1.0;
        end
    end
end

%ENCODING LOOP - Recompute the activity of output backward/associative
%population go
go=go+ci2;  %Equation E4.

%ENCODING LOOP -- Compute the change in INTERNAL backward connections Wig
for m3=0:n-1;
    for m2=1:n;
      for m=1:n;
       digW(m3*n+m2,m3*n+m)=go(1,m3*n+m2)*gi(1,m3*n+m);  %Equation E5.
       Wig(m3*n+m2,m3*n+m)=Wig(m3*n+m2,m3*n+m)+ digW(m3*n+m2,m3*n+m);
       if Wig(m3*n+m2,m3*n+m)>1.0;
           Wig(m3*n+m2,m3*n+m)=1.0;
       end
     end
   end
%End the m3 for loop.
end

%ENCODING LOOP - Compute the activity of output forward/planning population
co=olda;  %Equation E6.

%ENCODING LOOP - Update forward population ci
ci=(0.6*a + 0.6*co*Wc  >= 0.7);    %Equation E7
ci=double(ci);

%ENCODING LOOP -- Compute the change in forward connections Wc
dcW=ci'*co;   %Equation E8.

for m2=0:n-1;
    for m=0:n-1;
       Wc(m2*n+m+1,m*n+m2+1)=Wc(m2*n+m+1,m*n+m2+1)+dcW(m2*n+m+1,m*n+m2+1);
   end
end

co=gi;  %Equation E9.

%ENCODING LOOP -- Compute the change in output connections Wo
dWo=o'*co;   %Equation E10.
for m2=1:No;
    for m=1:N;
        Wo(m2,m)=Wo(m2,m)+ 0.3*dWo(m2,m); 
        if Wo(m2,m) >= 1.0;
            Wo(m2,m)=1.0;
        end
    end
end

%GRAPHICS VECTOR UPDATE.
if withgraphics==1
cA=vertcat(cA,a);
coldA=vertcat(coldA,olda);
cGO=vertcat(cGO,go);
cGI=vertcat(cGI,gi);
cCO=vertcat(cCO,co);
cO=vertcat(cO,o);
end

%ENCODING LOOP - Create a buffer for holding previous value of "a" which 
%does not get zeroed between encoding and retrieval.
olda=zeros(1,N);
olda=a;

%ENCODING LOOP - Update the buffer "ci2" which does NOT get zeroed between
%encoding and retrieval.
ci2=ci;

%ZERO OUT VECTORS
%The zeroing between encoding and retrieval may not be necessary - may only
%be necessary to zero after retrieval.
a=zeros(1,N);
o=zeros(1,No);
%ci=zeros(1,N);
%co=zeros(1,N);
%gi=zeros(1,N);
%go=zeros(1,N);

%RETRIEVAL LOOP
%Rsteps are the number of steps for retrieval between each goal input.  
Rsteps=2;
Gsteps=svsize+shsize;
%Increase this by increasing Gsteps so goal activity is sweeping
%backward continuously.

%GSTEPS LOOP
for goal=1:Gsteps;

%Set the goal vector goR before going through the retrieval loops.
goR((svsize*shsize+No)*n+1:(svsize*shsize+No)*n+n)=1;

%RSTEPS LOOP
for ret=1:Rsteps;
 
gi=((Wg*(goR + go >= 0.7)')' >= 0.7);  %Equation R1.
gi=double(gi);
%Calculate local feedforward inhibition arising from the portion of x in
%INDIVIDUAL minicolumns.
inhsumlgi=zeros(1,n);
for inh=0:n-1;
    inhsumlgi(inh+1)=sum(gi(inh*n+1:inh*n+n));
end
inhsumgi=inhsumlgi*Wa;

%RETRIEVAL LOOP - Recompute population "go"
go=(goR + (Wig*gi')' - 0.4*inhsumgi  >= 0.6);  %Equation R2 with goR added (CORRECTION).

%RETRIEVAL LOOP - Update population co where backward spread from the goal
%through populations "gi" and "go" interacts with current location ("olda").
co=(0.6*olda + 0.6*gi  >= 0.7);  %Equation R3.

%End the RSTEPS loop.
end
%End the GSTEPS loop.
end

cCO=vertcat(cCO,co);

%VIEW "ACTION VALUES" FOR NEOCORTEX
%This command shows the output function oq=Wo*co due to backward retrieval
%spread and Wo on each step (Wo*co=Wo*(a+gi)=Wo*(a+Wg*(Wig*Wg)^n*goR)
if action==0 & withgraphics==1;  %This conditional keeps the graphics from turning off on action cycles.
 oq=(Wo*co(1,1:N)')' ;     
   xQ(sv*3-2,sh*3-1)=64*(1-oq(1));  
   xQ(sv*3,sh*3-1)=64*(1-oq(2));
   xQ(sv*3-1,sh*3-2)=64*(1-oq(3));
   xQ(sv*3-1,sh*3)=64*(1-oq(4));
   %Note that I scaled it to 0.3 so large changes will appear as gray-scale.

   axes(axj);
   image(xQ,'EraseMode','none')
   set(axj,'XTick',[3.5,6.5,9.5,12.5], 'YTick',[3.5,6.5,9.5,12.5]);
   set(axj,'XGrid','on', 'YGrid','on','ZGrid','off','GridLineStyle','-');
   drawnow
%End of xQ output viewing (end action==0 conditional).
end

%End Neocortex==1 conditional.
end



%%%%%%%%%%%%%%%%%%%%%%%%%%

%TDLEARN TDLEARN TDLEARN TDLEARN TDLEARN
%COMPUTATIONS FOR TDLEARNING OF VALUE AND ACTION-VALUES
%TDlearning (Neocortex==0) CONDITIONAL
if Neocortex==0;
%Create value function for state elements.
dV((oldsv-1)*shsize+oldsh)=V((sv-1)*shsize+sh)+Reward-V((oldsv-1)*shsize+oldsh);
V((oldsv-1)*shsize+oldsh)=V((oldsv-1)*shsize+oldsh)+alpha*dV((oldsv-1)*shsize+oldsh);

%These graphics show the value function V of each state (not the Q)
xC(oldsv,oldsh)=64*(1-2*V((oldsv-1)*shsize+oldsh));
if withgraphics==1 | steps==Tsteps
    %This command makes the new image of the value function on each step
    axes(axk);
    image(xC,'EraseMode','none')
    drawnow
    %end the withgraphics conditional
end

%Update the graphics file for plotting the action-value function Q.
if Output==1
   xQ(oldsv*3-2,oldsh*3-1)=64*(1-Q((oldsv-1)*shsize*No+(oldsh-1)*No+Output));
elseif Output==2
   xQ(oldsv*3,oldsh*3-1)=64*(1-Q((oldsv-1)*shsize*No+(oldsh-1)*No+Output));
elseif Output==3
   xQ(oldsv*3-1,oldsh*3-2)=64*(1-Q((oldsv-1)*shsize*No+(oldsh-1)*No+Output));
elseif Output==4
   xQ(oldsv*3-1,oldsh*3)=64*(1-Q((oldsv-1)*shsize*No+(oldsh-1)*No+Output));
%End creating the Q action-value function output
end

%This command makes the new image of the action-value function Q on each step
if withgraphics==1 | steps==Tsteps
    axes(axj);
   image(xQ,'EraseMode','none')
   set(axj,'XTick',[3.5,6.5,9.5,12.5], 'YTick',[3.5,6.5,9.5,12.5]);
   set(axj,'XGrid','on', 'YGrid','on','ZGrid','off','GridLineStyle','-');
   drawnow
   %End the withgraphics conditional.
end
   
%Now implement the actor critic architecture by taking the TD error and
%using it to update the current Q value (based on the index of latest
%Output determined in the output section).  Note that the structure of the
%Q vector is basically that each element of the state vector has No
%possible outputs, so this vector involves going to the correct state
%location and then to the index for the output that was previously generated. 

Q((oldsv-1)*shsize*No+(oldsh-1)*No+Output)=Q((oldsv-1)*shsize*No+(oldsh-1)*No+Output)+alpha*dV((oldsv-1)*shsize+oldsh);

%End the Neocortex==0 conditional segment.
end
%END TDLEARN



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%OUTPUT OUTPUT OUTPUT OUTPUT OUTPUT OUTPUT
%COMPUTATION OF OUTPUT/ACTION. 

if Neocortex==0 | (Neocortex==1 & action==0);  
    %Conditional for Neocortex so only computes output on state cycles.
    
%Create a five element vector "rando" with random numbers between 0 and 1.
    z=sprand(1,1,1);
    z2=sprand(1,1,1);
    z3=sprand(1,1,1);
    z4=sprand(1,1,1);
    z5=sprand(1,1,1);
    
    rando=[z z2 z3 z4];
   
%This function sets how much random exploration there is even when there is
%prior learning.  If set it to 1.0 it will still explore initially, but
%learning will take over completely when available.
randomExploreThreshold=0.99;

if z5 > randomExploreThreshold;
       o=rando;
   elseif Neocortex==1;    
       o=(0.001*(rando)  + (Wo*co(1,1:N)')'  );   
   elseif Neocortex==0;
       o=0.001*rando + Q((oldsv-1)*shsize*No+(oldsh-1)*No+1:(oldsv-1)*shsize*No+(oldsh-1)*No+No);
   end
       
%Make sure that the output/action "o" won't push the rat into a barrier or
%off the edge of the field.
    if sv<=1 | B(sv-1,sh)==1 %| (sv-1==oldsv & sh==oldsh)
        o(1,1)=-100;
    end
    if sv>=svsize | B(sv+1,sh)==1 %| (sv+1==oldsv & sh==oldsh) 
        o(1,2)=-100;
    end
    if sh<=1 | B(sv,sh-1)==1 %| (sh-1==oldsh & sv==oldsv)
        o(1,3)=-100;
    end
    if sh>=shsize | B(sv,sh+1)==1 %| (sh+1==oldsh & sv==oldsv)
        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.  This section finds 
   [C,I]=max(o);
  
    %Before altering the sv and sh values on basis of output, save the
    %oldsv and oldsh values for computing TDlearning.
    oldsv=sv;
    oldsh=sh;
    
    %The output "I" was produced by the "max" function above (representing
    %the element of the vector which was the largest).
    %For the output "I" indicating component of output vector "o" which is
    %maximum, update the state of the network.
    
    %IF first element is largest, go Down.
       if I==1 
       sv=sv-1;
       sh=sh;
       o=zeros(1,No);
       o(1,1)=1;
       Output=1;
   %If second element is largest go Up.
   elseif I==2 
       sv=sv+1;
       sh=sh;
       o=zeros(1,No);
       o(1,2)=1;
       Output=2;
   %If third element is largest go Left.
   elseif I==3 
       sh=sh-1;
       sv=sv;
        o=zeros(1,No);
       o(1,3)=1;
       Output=3;
   %If fourth element is largest to Right.
   elseif I==4 
       sh=sh+1;
       sv=sv;
       o=zeros(1,No);
       o(1,4)=1;
       Output=4;
   end

%NEOCORTEX - Update action flag.
    %Update action flag so know an action generated and
    %that the "a" vector should be updated by proprioception and not by
    %location information.
    action=1;
    state=0;

%End action==0 conditional i.e. if Neocortex==0 | (Neocortex==1 & action==0)
end

%End the whole MAIN (ENCODING) LOOP (for steps=1:Tsteps).
end


%SHOW END GRAPHICS SUMMARIZING ACTIVITY DURING SIMULATION.
if withgraphics==1 && Neocortex==1
 
cGO=vertcat(cGO,Wg);
%Show the activity of pop gi with Wig at end.
xB=64*(1-cGO);
figure('Position',[40 0 400 400]);
h=gcf;
set(h,'Name','"go" vector with Wg at end');
set(h,'numbertitle','off');
colormap(gray);
image(xB')     

cGI=vertcat(cGI,Wig);
%Show the activity of pop gi with Wig at end.
xB=64*(1-cGI);
figure('Position',[80 0 400 400]);
h=gcf;
set(h,'Name','"gi" vector with Wig at end');
set(h,'numbertitle','off');
colormap(gray);
image(xB')     

cCO=vertcat(cCO,Wo);
%Show the activity of pop co with Wo at end.
xB=64*(1-cCO);
figure('Position',[120 0 400 400]);
h=gcf;
set(h,'Name','"co" vector with Wo at end');
set(h,'numbertitle','off');
colormap(gray);
image(xB')     

%Show the activity of pop "a"
xB=64*(1-cA);
figure('Position',[200 0 400 400]);
h=gcf;
set(h,'Name','"a" vector all steps');
set(h,'numbertitle','off');
colormap(gray);
image(xB')    
    
%Show the activity of pop olda
xB=64*(1-coldA);
figure('Position',[240 0 400 400]);
h=gcf;
set(h,'Name','Olda vector all steps');
set(h,'numbertitle','off');
colormap(gray);
image(xB')


%End final withgraphics==1 conditional.
end