function Saccades = netmovie(Type);
% Saccades = netmovie(Type) .... makes plot of activity in network over time
%  Type - 0 for sequential trial, 1 for parallel trial, 
%  returns Saccades - Sx2 matrix where S is number of saccades made, and
%                     for each is given 2 values, target choice and latency

N = 64;   % 2-D array activity
T = 1000;  % number of frames in movie
H = 5;   % number of arrays including input array
         % H==1 is the input
         % H==2 is the select all array
         % H==3 is the select two array
         % H==4 is the select one array, or build-up
         % H==5 is the burst outputs
         
Saccades = [];         
Alpha = 2.0;   % scaling param for remapping
Bias = -2.8;   % membrane threshold
Noise = 0.15;  % multiplicative noise
VisDelay = 60; % delay for retinal input to reach FEF, ms
Lambda = (1.-exp(-1./6.));  % membrane constant
Param = netpar;  % load weight parameters

%******** set experimental condition
if (Type==0)   % sequential presentation, 60 ms flash of each target
  WaitPeriod = 500;    % activate stop signal to simulate memory delay  
  SlowBuildUp = 0.0;  
  Targets = [ [50,110,(0.1*N),1.0];...
              [150,210,(0.3*N),1.0];...
              [250,310,(0.7*N),1.0]];
end
if (Type==1)  % parallel search trial, targets remain illuminated
  WaitPeriod = 0;
  SlowBuildUp = 2.0;   % 0.0 fast sac onset, 3.0 slower onset 
  Targets = [ [50,1000,(0.1*N),1.3];...
              [50,1000,(0.3*N),1.0];...
              [50,1000,(0.7*N),0.7]];
end

%********* set-up palette *********
map = zeros(100,3);
for i = 1:100
   map(i,:) = ((100-i)/100.);
end
colormap(map);

%******** weights are loaded ********
FWeights = cell(H,H);
for i = 1:H
   for j = 1:H
      it = (((i-1)*H)+j);
      FWeights{i,j} = fft(gaussian_weight_1d(N,Param(it,:)));
      if ((i==2)&(j==2))
         FOdd22 = Alpha * fft(dgaussian_weight_1d(N,Param(it,:)));
      end
      if ((i==3)&(j==3))
         FOdd33 = Alpha * fft(dgaussian_weight_1d(N,Param(it,:)));
      end
   end
end

%********* negative bias at fixation ***********************
NegFixBias22 = -2. * gaussian_1d(N,((N/2)+0.5),2.0);
NegFixBias33 = -1. * gaussian_1d(N,((N/2)+0.5),2.0);

%********* set-up activity ******
Activity = cell(1,H); 
FActivity = cell(1,H);
Preactivity = cell(1,H);
U = cell(1,H);
for j = 1:H
   Activity{j} = zeros(1,N);
   Preactivity{j} = zeros(1,N);
   U{j} = zeros(1,N);
   U{j} = Bias;
end
%********* fixation activity variables ********
fixact = [];
UFix = 0.0;
Fix = 0.5;

%*********** arrays for movies *******
Record = cell(T,H);
Acton = zeros(T,N);
M = moviein(1);

%*********** loop through step by step *****
making_saccade = 0;
for i = 1:T
   %************** remapping and saccades *********
   if ((max(Activity{H}) > 0.5)&(making_saccade==0))  % burst is on-going
      making_saccade = 1;
      tmpact = Activity{H};
      sacdir = sign( sum(tmpact(1,1:(N/2))) - sum(tmpact(1,((N/2)+1):N)) );
      FWeights{2,2} = FWeights{2,2} + (sacdir * FOdd22 );
      FWeights{3,3} = FWeights{3,3} + (sacdir * FOdd33 );
      sacpeak = find(Activity{H}==max(Activity{H}));
      %************* record the saccade: target and latency
      latency = i;
      target = 0;
      mindist = N;
      sacpeak = find(Activity{H}==max(Activity{H}));
      for z = 1:size(Targets,1)
         dist = abs( sacpeak - Targets(z,3) );
         if (dist>(N/2))        %check for toroid distance ...
            dist = N - dist;
         end
         if (dist<mindist)
            mindist = dist;
            target = z;
         end
      end
      Saccades = [Saccades ; [target,latency] ];
      %****************************************************
   end
   if ((max(Activity{H}) < 0.5)&(making_saccade==1))  % turn off
      making_saccade = 0;
      FWeights{2,2} = FWeights{2,2} - (sacdir * FOdd22 );
      FWeights{3,3} = FWeights{3,3} - (sacdir * FOdd33 );
   end
   %************** compute stop input *************
   XStop = 0.0;
   if (max(Activity{H}) > 0.5)
      %track shifting activity in select all array
      act = Activity{2};
      maxo = 0.0;
      for z = max(1,(sacpeak-4)):min((sacpeak+4),N)
         if (act(1,z)>maxo)
            maxo = act(1,z);
            maxi = z;
         end
      end
      shift = maxi-sacpeak;
      sacpeak = maxi;
      %************ update visual inputs with shift ****
      for z = 1:size(Targets,1)
         Targets(z,3) = Targets(z,3)+shift;
         if (Targets(z,3)>N)
            Targets(z,3) = Targets(z,3) - N;
         end
         if (Targets(z,3)<1)
            Targets(z,3) = Targets(z,3) + N;
         end
      end
      %********** update the stop input based on % saccade done ****
      D = -sacdir * ( sacpeak - ((N/2)+0.5) );
      if (D<8)
         XStop = 3.1*(8-D);   %turn off signal risesly sharply near fixation
      else
         XStop = 0.0;
      end
   end
   if (i<WaitPeriod)
      XStop = 4.0;
   end
   %*************** fix activation ****************
   uf = XStop + (2.5 * Fix) ...
        - (6.0 * mean( Activity{H})) ...   %inhib from burst
        - (18.0 * mean(Activity{H-1}));    %inhib from build-up
   UFix = UFix * (1.-Lambda);
   UFix = UFix + (Lambda * uf);
   Fix = activation_function(UFix);
   if (Noise)
      Fix = Fix * random('Normal',1.0,Noise);
   end
   fixact = [fixact Fix];
   %*************** compute visual input **********  
   Activity{1} = input_function(N,i,Targets,VisDelay);
   %*************** transfer act to fourier *******
   for j = 1:H
      FActivity{j} = fft(Activity{j});
   end
   %*************** compute convolutions **********
   for j = 2:H
      Preactivity{j} = 0.0;
      for k = 1:H
         Preactivity{j} = Preactivity{j} + ...
                FWeights{j,k} .* FActivity{k};
      end
      Preactivity{j} = real( ifft(Preactivity{j}) );
      Preactivity{j} = Preactivity{j} + Bias;
      %********** include neg bias at fixation in select all + two
      if (j==2)
         Preactivity{j} = Preactivity{j} + NegFixBias22;
      end
      if (j==3)
         Preactivity{j} = Preactivity{j} + NegFixBias33;
      end
      %********** include interaction from fixation to burst+build-up
      if (j==(H-1))
         Preactivity{j} = Preactivity{j} - (4. * Fix);  %build-up inhib from fix
         Preactivity{j} = Preactivity{j} - SlowBuildUp; % from 0-3, slows sac onset   
      end
      if (j==H)
         Preactivity{j} = Preactivity{j} - (12. * Fix); %burst inhib from fix
      end
      %*********************************************************
      U{j} = U{j} * (1.-Lambda);
      U{j} = U{j} + (Lambda * Preactivity{j} );
      act = activation_function(U{j});
      if (Noise)
        nu = random('Normal',1.0,Noise,1,N);
        Activity{j} = nu .* act;
     else
        Activity{j} = act;
     end
   end
   %************** record activity at time step ******
   for j=1:H
      Record{i,j} = Activity{j};
   end
end   

%*************** plot the results *********************
for j = 1:H+1
   subplot(H+1,1,(H+2-j));
   if (j<(H+1))
     for i = 1:T
      Acton(i,:) = Record{i,j};
     end   
     Acton = 50. * Acton;
     image(Acton');
     if (j==5)
         ylabel('Burst');
     end
     if (j==4) 
         ylabel('Build-Up');
     end
     if (j==3)
         ylabel('TakeTwo');
     end
     if (j==2)
         ylabel('TakeAll');
     end
     if (j==1)
         xlabel('Time (milliseconds)');
         ylabel('Inputs');
     end
     
  else
     plot(1:T,fixact,'k-');
     ylabel('Fixation');
     title('See Figure 7 of Manuscript For Details');
  end
end