function [model timing] = run_bgm_2007_save_all(stim, comp)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameter set up
%
% Parameters were constant through all simulations 
%       
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tstep = 0.1;
model.p.NumSpeeds = 3;
% Level 1
  model.p.A1 = 0.001;
  model.p.B1 = 1;
  model.p.C1 = 2;
  model.p.D1 = 0.25;
  model.p.phi1 = 0.1;
  model.p.E1 = 10.225;    % 10, K3 1 , K4 1 gets good ogl results not yosem
  % Fxy is a 2D symmetric Gaussian with sigma = 1
  x = meshgrid(-3:3); xx = sqrt(x.^2 + (x').^2);
  model.p.Fxy = (model.p.E1/(2*pi))*exp(-xx.^2)';  
  clear  x xx;
  model.p.G1 = 0.031623;
% Level 2
  model.p.A2 = 10;
  model.p.B2 = 1;
  model.p.C2 = 2;
  model.p.D2 = 0.01;
  model.p.K2 = 20;
% Level 3
  model.p.A3 = 1;
  model.p.B3 = 1;
  model.p.C3 = 1;
  model.p.K3 = 2;
  model.p.A4 = 10;
  model.p.B4 = 1;
  model.p.C4 = 1;
  model.p.K4 = 2;
% Level 4
  model.p.A5 = 0.1;
  model.p.B5 = 1;
  model.p.C5 = 0.01;
% Level 5 & 6
  model.p.A6 = 0.5;   
  model.p.B6 = 1;
  model.p.C6 =0.5;
  model.p.D6 = 0.5;
  model.p.theta6 = 0.2;
  if     comp == 0, model.p.vdD = [0 0 0 0 0];
  elseif comp == 1, model.p.vdD = [0 0.5 1 1 10]; 
  elseif comp == 2, model.p.vdD = [0 0 0 0 5];
  elseif comp == 3, model.p.vdD = [0.25 0.25 1 0.25 10];
  else   error('no such competitive environment');
  end
  model.p.A7 = 0.5;
  model.p.B7 = 1;
  model.p.C7 = 4;
  model.p.D7 = 0.25;
  model.p.E7 = 0.25;
  model.p.G7 = 0.1;
  model.p.theta7 = 0.2;
  [model.p.w model.p.M model.p.N] = generate_w(size(stim)/2^(model.p.NumSpeeds-1), [0.5 0.625], 3);  % function below
  model.p.L = generate_L(2,2,3,0.25);  % function below

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Model initialization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ss=1:model.p.NumSpeeds
    s = ['s' num2str(ss)];
    if ss == 1, model.I.(s)(:,:,1) = stim(:,:,1);
    else model.I.(s)(:,:,1) = shrink(stim(:,:,1), 2^(ss-1));
    end
    model.I.(s)(:,:,2) =  1-model.I.(s)(:,:,1);
    model.a.(s) = zeros(size(model.I.(s),1), size(model.I.(s),2), 2);
    model.b.(s) = model.a.(s);
    model.x.(s) = model.a.(s); 
    model.y.(s) = model.a.(s);
    model.z.(s) = ones(size(model.I.(s),1), size(model.I.(s),2), 2);
    model.c.(s) = zeros(size(model.I.(s),1), size(model.I.(s),2), 2, 8);
    model.e.(s) = model.c.(s);
    model.f.(s) = zeros(size(model.I.(s),1), size(model.I.(s),2), 1, 8);
end
s = ['s' num2str(model.p.NumSpeeds)];
model.q = zeros(size(model.I.(s),1), size(model.I.(s),2), 1, 8);
model.Q = model.q;
model.m = model.q;
model.r = zeros(model.p.M, 1);
model.R = model.r;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Simulation Integration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    tic
for t = 1:size(stim,3)*(1/tstep)-1
    % update the stimulus frame if necessary
    tt = tstep*t;
    if tt == floor(tt)
        for ss=1:model.p.NumSpeeds
            s = ['s' num2str(ss)];
            if ss == 1, model.I.(s)(:,:,1) = stim(:,:,tt+1);
            else model.I.(s)(:,:,1) = shrink(stim(:,:,tt+1), 2^(ss-1));
            end
            model.I.(s)(:,:,2) =  1-model.I.(s)(:,:,1);
        end
    end
    % run the model   
    model = bgm_2007(model, tstep);
    % logging, reporting, scripts etc
   % show_res1(model);
     xs = size(model.I.s3);
     heading(:,t) = [2:3:xs(2) 2:3:xs(2)]*[model.R==max(model.R)];
     ocos(:,:,t) = model.b.s2(:,:,1)-model.b.s2(:,:,2);
     trans(:,:, t) = model.y.s2(:,:,1) - model.y.s2(:,:,2);
     dir(:,:,t, :) = model.f.s2(:,:,1, :);
     MT(:,:,t, :) = model.Q;
     MSTd(:,t) = model.R;
end
    model.heading = heading;
    model.ocos = ocos;
    model.trans = trans;
    model.dir = dir;
    model.MT = MT;
    model.MSTd = MSTd;
   % show_res(model, model.heading);
    %show_res1(model);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% get timing data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    timing = toc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% END CONTROL CODE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUB FUNCTION: shrink is used to sum over groups of cells to shrink the array.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function x = shrink(y, f)
   if f > 1, ff = 1/f; 
   else ff = f; f = 1/ff; end
   for t = 1:size(y,3)
       for i=1:floor(size(y,1)*ff)
           for j=1:floor(size(y,2)*ff)
              x(i,j,t,:) = sum(sum(y(i*f-f+1:i*f, j*f-f+1:j*f, t, :)))/f^2;
           end
       end
   end
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUB FUNCTION: generate_L, generates the long-range filter.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [L Nl err] = generate_L(l, sigx, sigy, k)
    wid=25;
    mid=13;
    th = 0.005;
    L = zeros(wid,wid, 8);
    rot45 = [cos(45*(pi/180)) -sin(45*(pi/180)); sin(45*(pi/180)) cos(45*(pi/180))];
    for i= 1:wid
        for j = 1:wid
            x = i-mid; y = j-mid;
            L(i,j,1) = (l/(2*pi*sigx*sigy))*exp(-k*((x/sigx)^2 + (y/sigy)^2));
            nxy = round(sqrt(2)*rot45*[x;y])+mid;
            L(max(min(nxy(1)-x ,wid),1),max(min(nxy(2) ,wid),1),8) = L(i,j,1);
        end
    end
    L(:,:,1) = (L(:,:,1)>th).*L(:,:,1);
    L(:,:,8) = (L(:,:,8)>th).*L(:,:,8);
    err = sum(sum(L(:,:,1)))-sum(sum(L(:,:,8)));
    L(:,:,5) = L(:,:,1); L(:,:,3) = rot90(L(:,:,1)); L(:,:,7) = L(:,:,3);
    L(:,:,4) = L(:,:,8); L(:,:,2) = fliplr(L(:,:,8)); L(:,:,6) = L(:,:,2);
    
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUB FUNCTION: generate_w, generates the template weight matrices
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [w M N] = generate_w(sz, vert_pos, horz_spac)
   sizx = round(2*max(sz));
   [U, V] = meshgrid(-sizx:sizx, -sizx:sizx); 
   % convert to 8 directions
   d = (4*atan2(V,U)/pi+1); d= d+(d<1)*8;
   tw(:,:,1,1) = double(d<2).*(2-d) + double(d > 8).*(d-8); 
   tw(:,:,1,2) = double(d<3&d>=2).*(3-d) + double(d>1&d<2).*(d-1);
   tw(:,:,1,3) = double(d<4&d>=3).*(4-d) + double(d>2&d<3).*(d-2);
   tw(:,:,1,4) = double(d<5&d>=4).*(5-d) + double(d>3&d<4).*(d-3);
   tw(:,:,1,5) = double(d<6&d>=5).*(6-d) + double(d>4&d<5).*(d-4);
   tw(:,:,1,6) = double(d<7&d>=6).*(7-d) + double(d>5&d<6).*(d-5);
   tw(:,:,1,7) = double(d<8&d>=7).*(8-d) + double(d>6&d<7).*(d-6);
   tw(:,:,1,8) = double(d>=8).*(9-d) + double(d>7&d<8).*(d-7);
   for d=1:8, tw(sizx+1, sizx+1,1,d) = 1; end
   counter=0;
   w=zeros(sz(1), sz(2), 1, 8);
   for vert=[round(vert_pos*sz(1))]
       for horz=round(horz_spac/2):horz_spac:sz(2)
           counter = counter + 1;
           w(:,:, counter,:) =  tw(sizx-vert:sizx+((size(w,1)-1) -vert), sizx-horz+2:sizx+((size(w,2)-1) - horz+2),1, :);
       end
   end
   M = counter;
   N = min(sum(sum(sum(w,4),2),1));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUB FUNCTION: show_res, some code to show results
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function show_res1(model)
figure
subplot(4,2,1); imagesc(model.a.s2(:,:,1)); set(gca, 'YDir', 'normal'); colormap('gray')
subplot(4,2,2); imagesc(model.a.s2(:,:,2));set(gca, 'YDir', 'normal'); colormap('gray')
subplot(4,2,3); imagesc(model.b.s2(:,:,1));set(gca, 'YDir', 'normal'); colormap('gray')
subplot(4,2,4); imagesc(model.b.s2(:,:,2));set(gca, 'YDir', 'normal'); colormap('gray')
subplot(4,2,5); imagesc(model.y.s2(:,:,1));set(gca, 'YDir', 'normal'); colormap('gray')
subplot(4,2,6); imagesc(model.y.s2(:,:,2));set(gca, 'YDir', 'normal'); colormap('gray')
subplot(4,2,7); vec_me(model.f.s2(:,:,1, :), 4, 0.1);

function show_res(model, mstd)
    figure
    xs = size(model.I.s3);
    subplot(2,2,1);  vec_me(model.Q, 4, 0.1);
    subplot(2,2,2);  vec_me8(model.Q, 4, 0.1);
    subplot(2,2,3);  plot([2:3:xs(2)], model.R(1:model.p.M/2), 'b', [2:3:xs(2)], model.R(model.p.M/2+1:end), 'c')
                     hold on
                     heading = [2:3:xs(2) 2:3:xs(2)]*(model.R == max(model.R))
                     text(heading, max(model.R), ['*' num2str(heading)]);
    subplot(2,2,4);  plot(mstd, 'b'); axis([0 140 0 xs(2)+1])
    
    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MAIN MODEL FUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
function o = bgm_2007(in, tstep)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  bgm_2007 - Primate motion model of Browning, Grossberg & Mingolla, 2007
%             notation conventions follow paper.  This code is released as
%             part of that publication and as such is subject to the
%             copyright restrictions defined by the publisher.
%
%  Parameters  in - the model state and model parameters from the previous
%                   timestep
%
%              out - the model and model parameters, iterated by 1 timestep
%
%  A simple Euler's method is used to integrate the model from the input
%  timestep to the output timestep.  The code attempts to be as generic as
%  possible, using parameters for control structures, however as published
%  there were for example 3 speeds, 8 directions etc, some modifications may be
%  necessary to deal with different structural values.
%
%  Date Created:    July 27, 2007
%  Author:          N Andrew Browning
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MODEL CONTROL
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Since s corresponds to the scale of the input it cannot be easily
% incorporated into a higher dimensional matrix, therefore a structure is
% used, indexes i and j are the first 2 dimensions of the matrix, on-off
% index (p) is the 3rd dimension of the matrix, direction index(d) is
% included as the 4th dimension of the matrix
%
% get the parameters from the input and reduce notation
p = in.p; o.p = p;
% run the simulation
for ss=1:p.NumSpeeds
     s = ['s' num2str(ss)];
     % The input needs to be scaled and the on-off channels need to be
     % segregated such that in.I.s holds the scaled inputs
     
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
     % Level 1, On-Center-Off-Surround Processing for Boundaries
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      I = in.I.(s); a = in.a.(s); b = in.b.(s); o.I.(s) = I;
      o.a.(s) = a + tstep*(-p.A1*a + (p.B1 - a).*I*p.C1 - (p.D1 + a).*sumFxyI(p.Fxy, I)   ); % 1.1, sumFxy is in functions section below
                tr = max(o.a.(s)-in.p.phi1, 0).^2;                                               
      o.b.(s) = tr ./ (p.G1^2 + tr);                                                       % 1.3
      clear I a tr;
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % Level 2, Transient Cells
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      x = in.x.(s); z = in.z.(s); y = in.y.(s);
      o.x.(s) = x + tstep*(p.A2*(-p.B2*x + (p.C2-x).*b));       % 2.1
      o.z.(s) = z + tstep*(p.D2*(1 - z - p.K2*x.*z));           % 2.2
      o.y.(s) = max(o.x.(s).*o.z.(s),0);                               % 2.3
      clear b x z;
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % Level 3, Directional Transient Cells 
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      c = in.c.(s); e = in.e.(s); 
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % It is possible to do this algorithmically as shown in this
      % commented code but it's faster to hard code each direction as
      % implemented below.
      %      for d=1:8
      %      if d < 5, D =d+4; else D = d-4; end
      %      theta = (d-1)*(2*pi/8); xd =round(cos(theta)); yd=round(sin(theta));
      %
      %      XY = [zeros((yd==-1), size(I,2)); zeros(size(I,1)-(yd~=0), xd==-1) c((yd==1)+1:end-(yd==-1),(xd==1)+1:end-(xd==-1),t,D) zeros(size(I,1)-(yd~=0), xd==1); zeros((yd==1), size(I,2))];
      %
      %      c(:,:,t+1, d) = max(c(:,:,t,d) + tstep*(Ac*(-Bc*c(:,:,t,d)+Cc*I(:,:,it) - Kc*XY )),0);
      %      e(:,:,t+1, d) = max(e(:,:,t,d) + tstep*(Ae*(-Be*e(:,:,t,d)+Ce*I(:,:,it) - Ke*XY )),0);
      %    end
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      nm = size(c);
      % Rightward direction
      XY = max( [c(:,2:end,:,5) zeros(nm(1), 1,2)],0);
      o.c.(s)(:,:,:,1) = c(:,:,:,1) + tstep*(p.A3*(-p.B3*c(:,:,:,1) + p.C3*y - p.K3*XY));  % 3.1
      o.e.(s)(:,:,:,1) = e(:,:,:,1) + tstep*(p.A4*(-p.B4*e(:,:,:,1) + p.C4*y - p.K4*XY));  % 3.2
      % UP-Rightward direction
      XY = max( [c(2:end,2:end,:,6) zeros(nm(1)-1, 1,2); zeros(1, nm(2),2)],0);
      o.c.(s)(:,:,:,2) = c(:,:,:,2) + tstep*(p.A3*(-p.B3*c(:,:,:,2) + p.C3*y - p.K3*XY));  % 3.1 
      o.e.(s)(:,:,:,2) = e(:,:,:,2) + tstep*(p.A4*(-p.B4*e(:,:,:,2) + p.C4*y - p.K4*XY));  % 3.2
      % Upward direction
      XY = max( [c(2:end,:,:,7); zeros(1, nm(2),2) ],0);
      o.c.(s)(:,:,:,3) = c(:,:,:,3) + tstep*(p.A3*(-p.B3*c(:,:,:,3) + p.C3*y - p.K3*XY));  % 3.1 
      o.e.(s)(:,:,:,3) = e(:,:,:,3) + tstep*(p.A4*(-p.B4*e(:,:,:,3) + p.C4*y - p.K4*XY));  % 3.2
      % Up-Leftward direction
      XY = max( [zeros(nm(1)-1, 1,2) c(2:end,1:end-1,:,8); zeros(1, nm(2),2) ],0);
      o.c.(s)(:,:,:,4) = c(:,:,:,4) + tstep*(p.A3*(-p.B3*c(:,:,:,4) + p.C3*y - p.K3*XY));  % 3.1 
      o.e.(s)(:,:,:,4) = e(:,:,:,4) + tstep*(p.A4*(-p.B4*e(:,:,:,4) + p.C4*y - p.K4*XY));  % 3.2
      % Leftward direction
      XY =  max( [zeros(nm(1), 1,2) c(:,1:end-1,:,1) ],0);
      o.c.(s)(:,:,:,5) = c(:,:,:,5) + tstep*(p.A3*(-p.B3*c(:,:,:,5) + p.C3*y - p.K3*XY));  % 3.1 
      o.e.(s)(:,:,:,5) = e(:,:,:,5) + tstep*(p.A4*(-p.B4*e(:,:,:,5) + p.C4*y - p.K4*XY));  % 3.2
      % Down-Leftward direction
      XY = max( [zeros(1, nm(2),2); zeros(nm(1)-1, 1,2) c(1:end-1,1:end-1,:,2) ],0);
      o.c.(s)(:,:,:,6) = c(:,:,:,6) + tstep*(p.A3*(-p.B3*c(:,:,:,6) + p.C3*y - p.K3*XY));  % 3.1 
      o.e.(s)(:,:,:,6) = e(:,:,:,6) + tstep*(p.A4*(-p.B4*e(:,:,:,6) + p.C4*y - p.K4*XY));  % 3.2
      % Downward direction
      XY = max( [zeros(1, nm(2),2); c(1:end-1,:,:,3) ],0);
      o.c.(s)(:,:,:,7) = c(:,:,:,7) + tstep*(p.A3*(-p.B3*c(:,:,:,7) + p.C3*y - p.K3*XY));  % 3.1 
      o.e.(s)(:,:,:,7) = e(:,:,:,7) + tstep*(p.A4*(-p.B4*e(:,:,:,7) + p.C4*y - p.K4*XY));  % 3.2
      % Downward-Rightward direction
      XY = max( [zeros(1, nm(2),2); c(1:end-1,2:end,:,4) zeros(nm(1)-1, 1,2) ],0); 
      o.c.(s)(:,:,:,8) = c(:,:,:,8) + tstep*(p.A3*(-p.B3*c(:,:,:,8) + p.C3*y - p.K3*XY));  % 3.1 
      o.e.(s)(:,:,:,8) = e(:,:,:,8) + tstep*(p.A4*(-p.B4*e(:,:,:,8) + p.C4*y - p.K4*XY));  % 3.2
      
      %o.c.(s) = max(o.c.(s), 0);
      %o.e.(s) = max(o.e.(s), 0);
      
      clear XY y c;
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % Level 4, Feature Enhancement via Directional Competition
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      f = in.f.(s); 
      sump = sum(max(e,0),3); sumDp = sum(sump ,4);
      for d = 1:8, sumDnedSump(:,:,:,d) = sumDp - sump(:,:,:,d); end
      o.f.(s) = f + tstep*(-p.A5*f + (p.B5-f).*sump - (p.C5+f).*sumDnedSump);              % 4.1
      clear e sump sumDp sumDnedSump f;
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % Levels 5 & 6 are not scale specific, resize and combine scales
      % this is a combination of eqn 5.1 and the sum over s in 5.2
      % shrink function is implemented in functions section below
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      n = 2^(p.NumSpeeds-ss);
      if ss == 1 
          o.m = (1/n)*combine(max(o.f.(s),0), n);  % first speed so initialize o.m
      elseif ss == p.NumSpeeds
          o.m = o.m + (1/n)*max(o.f.(s),0);      % last speed so no need to resize
      else
          o.m = o.m + (1/n)*combine(max(o.f.(s),0), n);  % resize and sum across scales
      end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% there is edge effect from the directional transients 1 pixel wide
% around the array, remove this
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
o.m([1 end], :, :, :) = 0;
o.m(:, [1 end], :, :) = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Level 5, MT+      Optic Flow
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
q=in.q; r = in.r; Q = in.Q; R = in.R;
qEx = sumLxy(p.L, in.m).*((p.C6/p.M)*sumRz(R,p.w)+1) + p.D6*Q;
o.q = q + tstep*(-p.A6*q + (p.B6-q).*qEx - q.*sumvD(p.vdD, Q));       % 5.2
o.Q = max(q-p.theta6, 0).^2;                                        % 5.3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Level 6, MSTd     Heading
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ri = 1:size(p.w,3)
    w = p.w(:,:,ri,:);
    rEx(ri,1) = w(1:end)*Q(1:end)';
end
o.r = r + tstep*(-p.A7*r + (p.B7-r).*((p.C7/p.N)*rEx + p.D7*R) - p.E7*r.*(sum(R)-R));  % 5.6
rminth = max(o.r-p.theta7, 0).^2;
o.R = rminth ./ (p.G7^2 + rminth);                                                  % 5.7

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% END MAIN MODEL FUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MAIN MODEL SUB FUNCTIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% sumFxy is used in the on-center off-surround networks to create a
% Guassian cell neighborhood
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = sumFxyI(F, I)
% 2D Gaussian Filter (F) across input (I), implemented as 2 1D filtering
% operations for efficiency.
for o = 1:2  % input has both on and off channels
    if length(F) == 1                           % scalar filter
        warning('SumFxyI: filter is scalar.');
        out(:,:,o) = F*I(:,:,o);
    elseif size(F, 1) == 1 | size(F, 2) == 1    % Allow F to be passed as a 1D Gaussian and then filter as 2D
        out(:,:,o) = filter2(F/sum(F), I(:,:,o), 'same');
        out(:,:,o) = filter2((F/sum(F))', out(:,:,o), 'same');
    else                                        % 2D Filter
        out(:,:,o) = filter2(F, I(:,:,o), 'same');
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% sumLxy is used for the long range filter in MT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = sumLxy(L, I)
    % can't split the diagonal filters without rotating the input which
    % introduces other issues so use the full 2D filter
    for d=1:8
       out(:,:,1,d)  = filter2(L(:,:,d), I(:,:,1,d), 'same');
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% sumRz is used for the feedback from MSTd to MT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%       
function out = sumRz(R,w)
    out(:,:,1,:) = zeros(size(w,1), size(w,2), 1, 8);
    for z = 1:length(R)
        out(:,:,1,:) = out(:,:,1,:) + R(z)*w(:,:,z,:);
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% sumvD is used for competititive interactions in MT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%       
function out = sumvD(v,I)
    for d = 1:8
        out(:,:,1,d) = v(1)*I(:,:,:,d) + v(2)*sum(I(:,:,:,dof([d+1 d-1])),4) + v(3)*sum(I(:,:,:,dof([d+2 d-2])),4) + v(4)*sum(I(:,:,:,dof([d+3 d-3])),4) + v(5)*I(:,:,:,dof([d+4]));
    end
function y = dof(x)
   y = mod(x,8);
   y = (y==0)*8+y;
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% combine is used to sum over groups of cells to shrink the array.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function x = combine(y, f)
   if f > 1, ff = 1/f; 
   else ff = f; f = 1/ff; end
   for t = 1:size(y,3)
       for i=1:floor(size(y,1)*ff)
           for j=1:floor(size(y,2)*ff)
              x(i,j,t,:) = sum(sum(y(i*f-f+1:i*f, j*f-f+1:j*f, t, :)));
           end
       end
   end