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

% Excitable System Model
% Dynamic pattern formation and collisions in networks of excitable elements 
% 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 = 200;
% number of time steps   
   NT = 100;
% Izhikevich parameters;   
%   a = 0.01; b = 0.10; c = -65; dd = 2; vE = 1;
% Threshold value   
   vTH = 1.0;
% Coupling Strength function
   CE = 0.4;
   CI = 0.1;
   dE = 14;
   dI = 42;
   dM = 15;
   d0 = 5.4;
   WE = 35;   % 25.28  % 50
   WI = -70;     % -70  -90
   Iext = 0.03;  % 0.03
   tau = 20;
   
% =======================================================================   
% Initial states for active elements
% =======================================================================
    tMax = 100;                % sumulation time  [ms]
    t = linspace(0,tMax,NT);  % time  [ms] 
    dt = t(2) - t(1);         % time step  [ms];
    K = dt/tau;               % constant
    
    NA = 10;                  % number of initial activated neurons 
    %   load('Sact50v1','Srows','Scols');
    
    v = zeros(N,N,NT);         % membrane potential
    
%   v(48,50,1) = 30;
     v(100,100,1) = 1;
%    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;


  a = 20; b = 180;
  Srows = round(a + (b-a) .* rand(NA,1));
  Scols = round(a + (b-a) .* rand(NA,1));
%  
%  % load('Sact10v2','Srows','Scols');
%   
  for c = 1 : NA
    v(Srows(c),Scols(c),1) = 1;      
  end
  %save('Sact10v2','Srows','Scols');

   
% =======================================================================
% 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);
   d = sqrt(X.^2 + Y.^2);
  
   indexE = find(d < d0);
   indexI = find(d >= d0);
   
% Coupling Strengths
   W = CE .* exp(-d.^2 ./ dE) - CI .* exp(-d.^2 ./ dI);
   W(d>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(d > dM) = 0;
   
  
% ======================================================================
% TIME EVOLUTION OF STATES
% ======================================================================

% time step dt = 1;
  
 for tS = 2 : NT-5
    
 for m1 = 1 : N
 for n1 = 1 : N
           
   flagS = 2;
   if v(m1,n1,tS-1) >= 1;    flagS = 1; end;
   if  v(m1,n1,tS-1) < 0;    flagS = 3; end;    
   
   switch flagS
   
   case 1
     % 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 v(m2,n2,tS-1) >= 1; A(cR,cC) = 1; end;
         end    % n2
       end      % m2
     end         % cC
     end         % cR
           
     Stot = sum(sum(A .* Wn));
     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


%%
% figure(2)
%  for tS = 1 : NT
%     contourf(xG,yG,v(:,:,tS)) ;
%     %shading interp;
%     shading flat;
%     colorbar;
%     title(num2str(tS));
%     pause(0.001);
%  end


%%
% flag: 0 animated gif NOT saved or 1 saved /  file name
   f_gif = 1;
   ag_name = 'ag_LIF004.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]);

for tS = 1 : NT-5
    hold off
    
    z = find(v(:,:,tS) == 1);
    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);
   
   
  
    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)

% vN = zeros(NT,1);
% for cc = 1 : NT
%     vN(cc) = v(50,51,cc);
% end
    yP = zeros(NT,1);
    xP = t; 
    for tS = 1 :NT
     yP(tS) = v(100,100,tS);
    end
    plot(xP, yP);
    
    hold on
    for tS = 1 :NT
     yP(tS) = v(100,101,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
