% d:ap16/ns/scripts/ns_IZH_N012.m

% Izhikevich Quadratic Neuron spiking Model
% [2D] lattice of neurons with Mexicn Hat Function Coupling 
% Change code for different initial states at time step 1

% Ian Cooper
% 01 nov 2016

clear all
close all
clc
tic


% ========================================================================
% MODEL PARAMETERS
% ========================================================================
% save file flagSave == 1 otherwise flagSave == 0
   flagSave = 0;   
   fileName = 'N200T200P30M10v2.mat';
   
% Number of N x N lattice elements
   N = 100;
% Number of time steps   
   NT = 200;
% Izhikevich Model parameters    
    C = 50; vr = -60; vt = -40; k = 0.7; 
    a = 0.01; b = 5; c = -55; d = 150; 
    vPeak = 35;
    Sext = 20;
% Coupling Strength function
   CE = 0.4;
   CI = 0.1;
   dE = 14;
   dI = 42;
   dM = 15;
   d0 = 5.4;
  
   WE = 8000;   
   WI = -8000;  
   
   
% =======================================================================   
% Initial states for active elements
% =======================================================================
    dt = 1;                % time step   [ms]
    t = 0 : dt: (NT-1)*dt;  % time  [ms] 
   
  %  K = dt/tau;               % constant
    
    NA = 500;                  % number of initial activated neurons 
    %   load('Sact50v1','Srows','Scols');
    
    v = c .* ones(N,N,NT);         % membrane potential
   
    u = zeros(N,N,NT);         % membrane potential
%   v(48,50,1) = 30;
     v(50,50,1) = vPeak;
%    u(25, 25, 1) = 1;
%    u(10, 30, 1) = 1;
%    u(30, 13, 1) = 1;
%    u(38, 48, 1) = 1;
%    u(88, 88, 1) = 1;
%    u(70,  8, 1)  = 1;
%    u(60, 78, 1) = 1;


    aR = 20; bR = 80;
    Srows = round(aR + (bR-aR) .* rand(NA,1));
    Scols = round(aR + (bR-aR) .* rand(NA,1));
%  
%  % load('Sact10v2','Srows','Scols');
%   
  for cc = 1 : NA
    v(Srows(cc),Scols(cc),1) = vPeak;      
  end
  %save('Sact10v2','Srows','Scols');

    vP = v;
    u = 0 .* v;
    
% =======================================================================
% CALCULATIONS
% =======================================================================

% Lattice sites
   x = 1: N; y = 1: N;
   [xG,yG] = meshgrid(x,y); 
   
 % Distances
   xd = -dM:dM;    yd = xd;  NW = length(xd);
   [X,Y] = meshgrid(xd,yd);
   dXY = sqrt(X.^2 + Y.^2);
  
   indexE = find(dXY < d0);
   indexI = find(dXY >= d0);
   
% Coupling Strengths
   W = CE .* exp(-dXY.^2 ./ dE) - CI .* exp(-dXY.^2 ./ dI);
   W(dXY>dM) = 0;
   
% Normalised Coupling Strengths
   Wn = zeros(NW,NW);
   WnE = sum(sum(W(indexE)));
   WnI = sum(sum(W(indexI)));
   Wn(indexE) = WE .* W(indexE) ./ WnE;
   Wn(indexI) = WI .* W(indexI) ./ WnI;
   Wn(dXY > dM) = 0;
   
  
% ======================================================================
% TIME EVOLUTION OF STATES
% ======================================================================

% time step dt = 1;
  
 for tS = 2 : NT
    
 for m1 = 1 : N
 for n1 = 1 : N
           
%    flagS = 1;
%    if v(m1,n1,tS-1) >= vPeak; flagS = 2; end;
%   
   
   % if  v(m1,n1,tS-1) < 0;    flagS = 3; end;    
   
%    switch flagS
%    
%    case 1
%     vT = v(m1,n1,tS-1); uT = u(m1,n1,tS-1);
%     vT = vT + (dt/2) * (k*(vT - vr)*(vT - vt) - uT ) / C;
%     vT = vT + (dt/2) * (k*(vT - vr)*(vT - vt) - uT ) / C;
%     uT = uT + dt * a*(b*(vT - vr) - uT);  
       
%      % 5  
%      tSR =  round(1.01 + (10 - 1) * rand(1,1));
%     % tSR = 10;
%      v(m1,n1,tS:tS+tSR) = -1;
           
%    case 2
     A = zeros(NW,NW);   % activity matrix
     for cR = 1 : NW
     for cC = 1 : NW
       m2 = m1 - (dM+1) + cR;
       n2 = n1 - (dM+1) + cC;
       if m2 > 0 && m2 < (N+1)
         if n2 > 0 && n2 < (N+1)   
         if vP(m2,n2,tS-1) >= vPeak; A(cR,cC) = 1; end;
         end    % n2
       end      % m2
     end         % cC
     end         % cR
           
    Stot = sum(sum(A .* Wn)) + Sext;
    
    vT = v(m1,n1,tS-1); uT = u(m1,n1,tS-1);
    vT = vT + (dt/2) * (k*(vT - vr)*(vT - vt) - uT + Stot) / C;
    vT = vT + (dt/2) * (k*(vT - vr)*(vT - vt) - uT + Stot) / C;
    uT = uT + dt * a*(b*(vT - vr) - uT);
  
%    end   % switch 
   
    vP(m1,n1,tS) = vT;
    
   if vT >= vPeak    % a spike is fired!
       vP(m1,n1,tS) = vPeak; % padding the spike amplitude
       vT = c; % membrane voltage reset
       uT = uT + d; % recovery variable update
   end; % if
 
     
   v(m1,n1,tS) = vT; u(m1,n1,tS) = uT;
   
  
   
    %vT = v(m1,n1,tS-1);
    % vT = vT - K * vT + Iext + Stot; 
     
%      if vT > 1; vT = 1; end;
%      v(m1,n1,tS) = vT;
%     % if vT < -1; disp(vT); pause; end;
%      
%    case 3
%         
%    end   % switch
%   
%   if m1 == 50 && n1 == 49 
%         fprintf('tS = %2.0f \n',tS);
%         fprintf('flagS = %2.0f \n',flagS);
%         fprintf('flagT = %2.0f \n',flagT);
%         fprintf('1 =  %2.1f \n',v(m1,n1,tS-1));
%         fprintf('2 =  %2.1f \n',v(m1,n1,tS));
%         disp('  ');
%   end 
  
%   if v(m1,n1,tS) > 1; v(m1,n1,tS) = 1; end;
   
 end  % m1
 end  % n1 

end   % tS
%   for t = 1 : NT  
%   v(v(:,:,t)>=30) = 30;
%   end  
%if flagSave == 1;
%   save(fileName);
%end

%%
% GRAPHICS   =============================================================


% flag: 0 animated gif NOT saved or 1 saved /  file name
   f_gif = 0;
   ag_name = 'ag_IZHN.gif';
% Delay in seconds before displaying the next image  
   delay = 0.1;  
% Frame counter start
   nt = 1; 

%%  
figure(1)
 set(gcf,'units','normalized','position',[0.01 0.2 0.23 0.32]);
 ag_name = 'ag_IZH_N1.gif';  nt = 1;
for tS = 1 : NT   
    vP(1,1,tS) = -200;
    vP(N,N,tS) = vPeak;
   contourf(xG,yG,vP(:,:,tS));
   shading interp
   colorbar;
   tm1 = 't  =  ';  tm3 = '  ms';
   tm2 = num2str(tS*dt,'%2.1f\n');
   tm = [tm1 tm2 tm3];
   title(tm);
   pause(0.01) 
   
        if f_gif > 0 
          frame = getframe(1);
          im = frame2im(frame);
          [imind,cm] = rgb2ind(im,256);
     %  On the first loop, create the file. In subsequent loops, append.
          if nt == 1
            imwrite(imind,cm,ag_name,'gif','DelayTime',delay,'loopcount',inf);
          else
           imwrite(imind,cm,ag_name,'gif','DelayTime',delay,'writemode','append');
          end
       nt = nt+1;
       end
   
   
end
%%   
figure(5)
set(gcf,'units','normalized','position',[0.3 0.2 0.23 0.32]);
ag_name = 'ag_IZH_N5.gif';  nt = 1;
for tS = 1 : NT
    hold off
    
    z = find(vP(:,:,tS) == 35);
    xP = xG(z); yP = yG(z);
    
    hPlot = plot(xP, yP, 'o');
    set(hPlot,'markerFaceColor',[1 0 0],'markerEdgeColor',[1 0 0], ...
        'markerSize',3);
%     hold on
%     z = find(v(:,:,tS) == -1);
%     xP = xG(z); yP = yG(z);
%     
%     hPlot = plot(xP, yP, 'o');
%     set(hPlot,'markerFaceColor',[0 0 0],'markerEdgeColor',[0 0 0], ...
%         'markerSize',3);
   
    axis([0 N 0 N]);
    tm1 = num2str(t(tS),'%2.1f \n');
    tm2 = '   ms';
    tm = [tm1 tm2];
    title(tm);
    pause(0.1);
    
       if f_gif > 0 
          frame = getframe(5);
          im = frame2im(frame);
          [imind,cm] = rgb2ind(im,256);
     %  On the first loop, create the file. In subsequent loops, append.
          if nt == 1
            imwrite(imind,cm,ag_name,'gif','DelayTime',delay,'loopcount',inf);
          else
           imwrite(imind,cm,ag_name,'gif','DelayTime',delay,'writemode','append');
          end
       nt = nt+1;
      end
end


%     box on 
%     
%      if f_gif > 0 
%        frame = getframe(1);
%        im = frame2im(frame);
%        [imind,cm] = rgb2ind(im,256);
%      %  On the first loop, create the file. In subsequent loops, append.
%        if nt == 1
%           imwrite(imind,cm,ag_name,'gif','DelayTime',delay,'loopcount',inf);
%        else
%          imwrite(imind,cm,ag_name,'gif','DelayTime',delay,'writemode','append');
%        end
%        nt = nt+1;
%      end
%     
%     z = find(v(:,:,tS) <0);
%     xP = xG(z); yP = yG(z);
%     hPlot = plot(xP, yP, 'o');
%     set(hPlot,'markerFaceColor',[1 1 1],'markerEdgeColor',[1 1 1], ...
%     'markerSize',3); 
% end
%     
% 
% %%

%% -----------------------------------------------------------------------
figure(3)
   set(gcf,'units','normalized','position',[0.6 0.2 0.23 0.32]);
%  vN = zeros(NT,1);
%  for cc = 1 : NT
%      vN(cc) = v(51,51,cc);
%  end
    yP = zeros(NT,1);
    xP = t; 
    for tS = 1 :NT
     yP(tS) = vP(51,51,tS);
    end
    plot(xP, yP);
    
%     hold on
%     for tS = 1 :NT
%      yP(tS) = v(50,51,tS);
%     end
%      plot(xP, yP);

% %      for tS = 1 :NT
% %      yP(tS) = v(100,99,tS);
% %     end
% %      plot(xP, yP);
% %   legend('100','101','99');
% %    

%%

% END ====================================================================
   
toc
