function MH_grid(N_grid,Ni,swn,range)
% Mexican-Hat grid of Izhikevich neurons
%
% Script by Filip Piekniewski 2009
%
% N_grid - size of the grid
% Ni - number of inhibitory neurons
% swn - number of random connections
% range - time of simulation (ms)
% Sample use: 
% MH_grid(50,1,0,1000)
%
%

Ne=N_grid*N_grid;
h=0;
MH=0.3;

K=min(floor((N_grid)/2),10);
D=zeros(K+1,K+1);
 for i=-K:K
    for j=-K:K
        d=0.5*sqrt(i^2+j^2);
        D(K+1+i,K+1+j)=50*MH*(exp(-MH*d^2) - 0.5* exp(-0.5*MH*d^2));
    end;
 end;   

h=(((2*K)^2)+2*Ni)*Ne+Ni*(2*Ne)+Ni*(Ni-1)+swn;
I=zeros(h,1);
J=zeros(h,1);
V=zeros(h,1);
h=1; 
 % Computing sparse matrix elements
 for n1=1:Ne % for every excitory neuron    
          
          x=floor((n1-1) ./ N_grid)+1;
          y=mod(n1-1, N_grid)+1;
          
          for i=-K:K-1
              for j=-K:K-1
                  x1=mod(x+i-1,N_grid)+1;
                  y1=mod(y+j-1,N_grid)+1;
                  I(h)=n1;
                  J(h)=(x1-1)*N_grid+y1;
                  V(h)=D(K+1+i,K+1+j);
                  h=h+1;           
              end; 
          end;
          for n2=(Ne+1):(Ne+Ni)
              I(h)=n1;
              J(h)=n2;
              V(h)=10 * (1/(Ne+Ni));
              h=h+1;
              I(h)=n2;
              J(h)=n1;
              V(h)=10 * (1/(Ne+Ni));
              h=h+1;
              
          end;    
      end;
      for n1=(Ne+1):(Ne+Ni)
          for n2=1:Ne
              I(h)=n1;
              J(h)=n2;
              V(h)=-20;
              h=h+1;
              I(h)=n2;
              J(h)=n1;
              V(h)=-20;
              h=h+1;
          end;
          for n2=(Ne+1):(Ne+Ni)
              if n1~=n2
                   I(h)=n1;
                   J(h)=n2;
                   V(h)=1;
                   h=h+1;
              end;
          end;    
      end;    
      for n=1:swn % small world network stuff
          n1=floor(Ne*rand(1));
          n2=floor(Ne*rand(1));
          if n1==0 n1=1; end;
          if n2==0 n2=1; end;
          if (n2~=n1)
                  I(h)=n1;
                  J(h)=n2;
                  V(h)=-10*rand(1);
                  h=h+1; 
                  I(h)=n1;
                  J(h)=n2;
                  V(h)=V(h-1);
                  h=h+1;
          end;   
      end;
     S=sparse(I,J,V);  
     % The weight matrix (sparse) is created.
     
     
cont=1;     
last_fired=zeros(N_grid+1,N_grid+1);
last_fired_aux=ones(N_grid+1,N_grid+1);
stimuli=zeros(N_grid+1,N_grid+1);


blob_count=[];
firing_count=[];
blob_list_x=[];
blob_list_y=[];
blob_list_timespan=[];
blob_list_status=[];
blob_trajectory_x=[];
blob_trajectory_y=[];
blob_trajectory_l=[];
blob_distance=[];
fired_x=[];
fired_y=[];
blob_spikes=[];

distances=[];
timespans=[];
spikes=[];


Ni=1;
re=rand(Ne,1); % Generate Ne random values from uniform dist
ri=rand(Ni,1); % Generate Ni random values from uniform dist
a=[0.02*ones(Ne,1);    0.02+0.08*ri];
b=[0.2*ones(Ne,1);     0.25-0.05*ri];
c=[-65*ones(Ne,1);      -65*ones(Ni,1)]; % reseting voltage for the fired neurons
d=[2*ones(Ne,1);          2*ones(Ni,1)]; % reseting u for the fired neurons
v=-65*ones(Ne+Ni,1);  % Initial values of v
u=b.*v;               % Initial values of u

       for i=1:N_grid
         for j=1:N_grid % initialization of voltage (modify if you like)
        %  v(transform2(i,j,N_grid))=150; u(transform2(i,j,N_grid))=u(transform2(i,j,N_grid))+2; 
         end;    
       end; 

               
firings=[];                            % spike timings
firings2=[ones(N_grid)*ones(N_grid)];  % spike patterns
voltage=[N_grid,N_grid];
voltage=[ones(N_grid)*ones(N_grid)]; 
uu=[N_grid,N_grid+1];
uu=[ones(N_grid)*ones(N_grid)]; 
 
ff=figure;

for t=1:range          % 
   blobs=zeros(N_grid,N_grid);
   if mod(t,100)==0 
       disp(['t=',num2str(t)]); 
   end; 
   I=[5*randn(Ne,1);2*ones(Ni,1)]; % thalamic input 
   for nneuron=1:Ne
       stimuli(getX(nneuron,N_grid),getY(nneuron,N_grid))=I(nneuron);
   end;    
   fired=find(v>=30); % indices of spikes
   firing_count(t)=size(fired,1);
   if ~isempty(fired)  % is not empty fired     
      firings=[firings; t+0*fired, fired];     
      I=I+sum(S(:,fired),2);
      v(fired)=c(fired);  
      u(fired)=u(fired)+d(fired);
      fired_x=[];
      fired_y=[];
      for fneuron=1:size(fired);
          neuron_x=getX(fired(fneuron),N_grid);
          neuron_y=getY(fired(fneuron),N_grid);
          fired_x(end+1)=neuron_x;
          fired_y(end+1)=neuron_y;
          last_fired(neuron_x,neuron_y)=t;
          [m,ind]=min((blob_list_x-neuron_x).^2+(blob_list_y-neuron_y).^2);
          if sqrt(m)<5
              blob_spikes(ind)=blob_spikes(ind)+1;
          end;
      end;    
      
   end;
   v=v+0.5*(0.04*v.^2+5*v+140-u+I); %step 0.5ms
   v=v+0.5*(0.04*v.^2+5*v+140-u+I); %for stability
   u=u+a.*(b.*v-u);  
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   
         for xx=1:Ne
          if v(xx)>10 
              voltage(getX(xx,N_grid),getY(xx,N_grid))=10;
          else
              voltage(getX(xx,N_grid),getY(xx,N_grid))=v(xx);
          end
          uu(getX(xx,N_grid),getY(xx,N_grid))=u(xx);
          if v(xx)>=-30 
              blobs(getX(xx,N_grid),getY(xx,N_grid))=1;
          end;    
      end;
      blob_centers_x=[];
      blob_centers_y=[];
      blobs1=bwmorph(blobs,'clean');
      blob_list_status=0*blob_list_x;
      [blobs,n]=bwlabel(blobs1);
      for i=1:n
          [rr, cc] = find(blobs==i);
          b_x=mean(rr);
          b_y=mean(cc);
          
          blob_centers_x(end+1)=b_x;
          blob_centers_y(end+1)=b_y;
          
          [m,ind]=min((blob_list_x-b_x).^2+(blob_list_y-b_y).^2);
          if sqrt(m)<3 
              
              blob_trajectory_l(ind)=blob_trajectory_l(ind)+1;
              blob_trajectory_x(ind,blob_trajectory_l(ind))=b_x;
              blob_trajectory_y(ind,blob_trajectory_l(ind))=b_y;
              
              blob_distance(ind)=blob_distance(ind)+sqrt((blob_list_x(ind)-b_x)^2+(blob_list_y(ind)-b_y)^2);
              
              blob_list_x(ind)=b_x;
              blob_list_y(ind)=b_y;
              blob_list_timespan(ind)=blob_list_timespan(ind)+1;
              blob_list_status(ind)=2; %supported;
          else % New blob is born
                
              blob_list_x(end+1)=b_x;
              blob_list_y(end+1)=b_y;
              blob_list_timespan(end+1)=1;
              blob_list_status(end+1)=1; % created;
              blob_trajectory_l(end+1)=1;
              blob_distance(end+1)=0;
              blob_trajectory_x(end+1,blob_trajectory_l(end))=b_x;
              blob_trajectory_y(end+1,blob_trajectory_l(end))=b_y;
              blob_spikes(end+1)=0;

          end;    
      end;    
      blob_count(t)=n;
      surviving=find(blob_list_status>0);
      dying=find(blob_list_status==0);
      if ~isempty(dying)
           for dead=1:size(dying,2) 
                 timespans(end+1)=blob_list_timespan(dying(dead)); 
                 distances(end+1)=blob_distance(dying(dead)); 
                 spikes(end+1)=blob_spikes(dying(dead))/blob_list_timespan(dying(dead));
                 
           end;
      end;
      blob_list_x=blob_list_x(surviving);
      blob_list_y=blob_list_y(surviving);
      blob_list_timespan=blob_list_timespan(surviving);
      blob_trajectory_l=blob_trajectory_l(surviving);
      blob_spikes=blob_spikes(surviving);
      blob_trajectory_x=blob_trajectory_x(surviving,:);
      blob_trajectory_y=blob_trajectory_y(surviving,:);
      subplot(2,2,1);
      if (mod(t,10)==0) % draw every 10 steps
        if (cont==1)
          contourf(voltage)
        else     
          surf(voltage)
          f = findobj('Type','surface');
          set(f,'FaceLighting','phong');
          shading interp
          az=-90;el=90;
          view(az,el);
        end;
        title('V (voltage)')
        colormap(jet)
        axis([0,N_grid,0,N_grid]);
        axis square;
        subplot(2,2,2);
        plot(fired_y,fired_x,'.');
        hold on;
        plot(blob_centers_y,blob_centers_x,'ro','MarkerSize',10);
      
        for i=1:size(blob_trajectory_l,2)
          if (blob_trajectory_l(i)>1)
            plot(blob_trajectory_y(i,1:blob_trajectory_l(i)),blob_trajectory_x(i,1:blob_trajectory_l(i)),'k','LineWidth',0.5);
          end;
        end;    
      
        hold off;
        title('Firings and blob centers')
        axis([0,N_grid,0,N_grid],'square');      
        subplot(2,2,3);      
        if (cont==1)
          contourf(uu);
        else
          surf(uu)
          f = findobj('Type','surface');
          set(f,'FaceLighting','phong');
          set(gca,'Zlim',[-90,50])
          shading interp
          az=-90;el=90;
          view(az,el);  
        end;
        title('U (recovery)')
        colormap(jet)
        axis([0,N_grid,0,N_grid]);
        axis square;
           
        subplot(2,2,4);
        pcolor(blobs);
        axis square;
        drawnow;
        title('Blobs')
        % Uncoment to print animation frames 
        %print(ff, '-r150', '-dpng', ['figure_oct', num2str(t/5),'.png']);
      end;
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;
print(ff, '-r300', '-depsc', ['domain', num2str(N_grid)]);
saveas(ff,['domain',num2str(N_grid),'.fig']);

ff1=figure;
plot(firings(:,1),firings(:,2),'.');
title('Spike raster');
xlabel 'Time (ms)';

print(ff1, '-r300', '-depsc', ['spike_raster', num2str(N_grid)]);
saveas(ff1,['spike_raster',num2str(N_grid),'.fig']);

ff2=figure;
hist(timespans,40);
title('Blob timespans (all)');
print(ff2, '-r300', '-depsc', ['Timespans_all', num2str(N_grid)]);
disp(['Median blob timespans (all) :',num2str(median(timespans))]);
saveas(ff1,['blob_timespans_all',num2str(N_grid),'.fig']);

ff5=figure;
timespans1=timespans(timespans>2);
hist(timespans1,40);
title('Blob timespans');
print(ff5, '-r300', '-depsc', ['Timespans', num2str(N_grid)]);
disp(['Median blob timespans (>2) :',num2str(median(timespans))]);
saveas(ff5,['blob_timespans',num2str(N_grid),'.fig']);

ff4=figure;
hist(distances,40);
title('Blob distances (all)');
print(ff4, '-r300', '-depsc', ['Distances_all', num2str(N_grid)]);
saveas(ff4,['blob_distances_all',num2str(N_grid),'.fig']);
disp(['Median distances travelled (all) :',num2str(median(distances))]);



ff6=figure;
distances=distances(distances>1);
hist(distances,40);
title('Blob distances');
print(ff6, '-r300', '-depsc', ['Distances', num2str(N_grid)]);
saveas(ff6,['blob_timespans',num2str(N_grid),'.fig']);
disp(['Median distances travelled (>1) :',num2str(median(distances))]);

ff7=figure;
hist(spikes,40);
title('Average spikes/blob/ms (all)');
print(ff7, '-r300', '-depsc', ['Spikes_all', num2str(N_grid)]);
saveas(ff7,['blob_spikes_all',num2str(N_grid),'.fig']);

ff8=figure;
spikes1=spikes(timespans>10);
hist(spikes1,40);
title('Spikes for blobs surviving 10 steps');
print(ff8, '-r300', '-depsc', ['Spikes_10', num2str(N_grid)]);
saveas(ff8,['blob_spikes_10_',num2str(N_grid),'.fig']);

ff3=figure;
plot(2:range,blob_count(2:range));
hold on
plot(2:range,firing_count(2:range),'r');
title('Firing activity and blob count');
xlabel 'Time (ms)';
print(ff3, '-r300', '-depsc', ['Firing_and_blobs', num2str(N_grid)]);
saveas(ff1,['Firing_and_blobs',num2str(N_grid),'.fig']);     
     
     
     
     
     
     
function A = getX(nr, size)
 
   A=floor((nr-1) ./ size)+1;

function A = getY(nr, size)
 
   A=mod(nr-1, size)+1;
