function jforest5();

clear all;


%Each cell has different connectivity
%connectivities updated regularly based on firing correlations


%PARAMETER SPECIFICATION

N=40;  % forest is NxN
it=600; % # iterations
%v=round(it/5); %update display every v frames
v=1;

t=8; %max number turns spent in refractory period
d=.6; %decay constant: the amount by which old activity decays each turn, equal to exp(-1/tau)
th=1; %threshold for a cell to become active
c=0; %rectification level.  any activation below c becomes equal to c
a=1; %maximum magnitude of the random activation
	 %specify the distribution below, within the main loop
	 
histrecall=80; %Number of displayed past histograms (all are saved and disp at end)
histhist=zeros(it,1);
	 
h=50; %update weights every h timesteps

uact=.5; %maximum increase in weight for positively correlated firing
Updateact(1)=0.95;  %decay multiplier for weights that contributed to no correlated firing
%update instructions are stored in the matrix Updateact.  It is displayed in figure 14
ua=h; Updateact(2:ua+1)=(uact)*(-1*exp(-.8*((1:ua)))+ones(1,ua))+1;
%ua is the maximum number of positive correlations that any single weight can have

uinh=.5; %maximum decrease in weight (make more negative) for successful inhibition
Updateinh(1)=1;  %decay multiplier for weights that did not successfully inhibit
%update instructions are stored in the matrix Updateinh.  It is displayed in figure 15
ui=2*h; Updateinh(2:ui+1)=(uinh)*(-1*exp(-.1*((1:ui)))+ones(1,ui))+1;
%ui is the maximum number of successful inhibition points that any single weight can earn
%(these points are defined below)


%PERIODIC CONDITIONS

p=0; %periodic firing patterns come up every p epochs  (set to 0 for never)
pend=0;
%they take the form of this matrix:
	%A big ol' X
%P=eye(N)+fliplr(eye(N));
	%Spiral conditions
%P=zeros(N);P(round(N/2),1:round(2*N/3))=ones(1,round(2*N/3));
	%???
%P(round(N/2)+3,1:round(N/3))=-10000*ones(1,round(N/3));
	%A line up and down the middle
P=zeros(N);P(1:N,round(N/2))=ones(N,1);

Pend=zeros(N);Pend(round(N/4):round(N/4)+1,round(N/4):round(N/4)+1)=ones(2);



%MATRIX INITIALIZATIONS
B=zeros(N);  %burning matrix: 1 = firing right now; 0 = not firing right now
A=zeros(N);  %activation matrix: nonnegative real numbers
        %represents the pre-firing membrane potential (for trees, temperature)
T=th*i*ones(N);  %thresholds matrix: positive real numbers
R=zeros(N);	%random activations matrix: this gets filled in (randomly)
			%and added to A each epoch


%INITIAL CONDITIONS

%T(round(N/2)+3,1:round(N/2))=T(round(N/2)+3,1:round(N/2))+102400*ones(1,round(N/2));



%DISPLAY SETUP
	
%set all the colormaps and positions
%the network display matrix
figure(3)
colormap(gray)
%shows which cells are firing, color-coded by cluster
figure(7)
colormap(cool)
%Cell (25,15)'s firing state (either on or off)
figure(1)
colormap(hot)
%firing correlations for cell (25,15)
figure(2)
colormap(gray)
%connectivity matrix for cell (25,15)
figure(5)
colormap(gray)
%current histogram of cluster sizes
figure(6)
%shows most recent 80 histograms, line by line.
figure(8)
colormap(pink)
%shows ALL histograms, line by line.  Displayed only at the very end of all epochs
figure(9)
colormap(pink)

figure(10)
colormap(hot)
figure(11)
colormap(hot)
figure(12)
colormap(gray)
figure(13)
colormap(gray)

%shows how to update weights for positive correlations
figure(14)
plot(0:ua,Updateact)
%shows how to update weights for successful inhibitions (defined below)
figure(15)
plot(0:ui,Updateinh)

%make border for individual connectivity display (used in DRAW IT section of MAIN LOOP)
disppadded=50*ones(2*c+3);


pause(4) %show what we've drawn so far



%SET UP THE CONNECTIVITY MATRICES
tic

%All cells connect to their neighbors in different ways
%Cell cc,dd has the connectivity matrix C(:,:,cc,dd)

%define the connectivity at each cell
%chose from among the following patterns:
hc=0;	%	hard circle
nb=0;	%	nearset neighbor
dog=1;	%	DOG

%parameters of the connectivity (note: some connectivity types don't use every parameter)
c=5; %radius of the connectivity matrix
ri=2; %radius of inner activation
ro=3; %radius of outer inhibition
ni=.4; %noise on inner radius (for DOG, = max proportion to deviate from smooth gaussian)
no=.1; %noise on outer radius (for DOG, = max proportion to deviate from smooth gaussian)

si=4; %strength of inner activation (positive for excit, negative for inhib)
so=-1; %strength of outer activation

%premade parameter sets
%DOG: strong act,strong inh
%dog=1; c=5; ri=1; ro=3; ni=.4; no=.1; si=4; so=-0.1;


figure(17)
colormap(gray)
X=[-1*c:1*c]/1;
plot(X,(si*exp(-((X/ri).^2))+so*exp(-((X/ro).^2))).*abs(sign(X)),X,ones(2*c+1),X,zeros(2*c+1))
pause(0.01)

for cc=1:N;for dd=1:N; %each cell has diff initial connectivity
%for cc=1:1;for dd=1:1;  %each cell has same initial connectivity (also set 0 to 1 on line 135)
	
	%hard circle of excit, hard circle of inhib (noisy boundary, but all weights = 1 or -1)
	if hc;
		%Defines C, the connectivity matrix.
		for ii=1:(2*c+1); for jj=1:(2*c+1);
			C(ii,jj,cc,dd)=2*(((ii-c-1)^2+(jj-c-1)^2)^0.5<ri+ni*rand(1))-(((ii-c-1)^2+(jj-c-1)^2)^0.5<ro+no*rand(1));
		end; end;
		C(:,:,cc,dd)=0.5*C(:,:,cc,dd);
	end;
	
	%nearest neighbor connectivity
	if nb;
		%C=[0.7 1 .7;1 0 1;0.7 1 0.7];
		C(:,:,cc,dd)=[0 1 0;1 0 1;0 1 0];
	end;
	
	%Difference of Gaussians
	if dog;
		for ii=1:(2*c+1); for jj=1:(2*c+1);
			r=((ii-c-1)^2+(jj-c-1)^2)^0.5;
			C(ii,jj,cc,dd)=si*exp(-(r/ri).^2)*(1-ni+ni*rand)+(so*exp(-(r/ro).^2))*(1-no+no*rand);
			%C(ii,jj,cc,dd)=(si*exp(-(r/ri).^2)+so*exp(-(r/ro).^2));
		end; end;
	end;
	
	%display each of the connectivity  matrices (time consuming!)
	if 0;
		disp=zeros(2*c+3);  %disp = C(:,:,:,:) plus a border whose color corresponds to 0.
							%it's useful for getting a sense of the magnitudes
							%excitation and inhibition
		for ii=2:(2*c+2); for jj=2:(2*c+2);
			disp(ii,jj)=C(ii-1,jj-1,cc,dd);
		end; end;
		figure(3)
		colormap(gray)
		image(50+100./(max(so,si)).*disp)
		pause(0.1)
	end;

end; end;  %for cc=1:N;dd=1:N;


%make all other connectivities like the first
if 0;
	for cc=2:N;for dd=2:N;
		for ii=1:(2*c+1); for jj=1:(2*c+1);
			r=((ii-c-1)^2+(jj-c-1)^2)^0.5;
			C(ii,jj,cc,dd)=C(ii,jj,1,1);
			%C(ii,jj,cc)=(si*exp(-(r/ri).^2)+so*exp(-(r/ro).^2));
		end; end;
	end;end;
end;


%these matrices store firing correlations and successful inhibitions
for cc=1:N; for dd=1:N;
	Hact(:,:,cc,dd)=zeros(2*c+1); %Hact stores positive correlations (I fire, you fire)
	Hinh(:,:,cc,dd)=zeros(2*c+1); %Hinh stores negative (I fire, you don't)
end;end;


%DISPLAY THE CONNECTIVITY MATRICES FOR SOME CELLS

if 0;  %(time consuming)
	for cc=1:rem(N,2*c+1); for dd=1:rem(N,2*c+1);
		tempdisp=zeros(N);
		for jj=0:floor(N/(2*c+1));for kk=0:floor(N/(2*c+1));
			temp=zeros(N);
			temp(cc+(2*c+1)*jj,dd+(2*c+1)*kk)=1;
			tempdisp=tempdisp+conv2(temp,C(:,:,cc+(2*c+1)*jj,dd+(2*c+1)*kk),'same');
		end;end;
		figure(3)
		image(25+50*tempdisp);
		pause(1)
	end;end;
end;

%print out how long it took to set up the matrices
MatrixSetupTime=1
toc



%MAIN LOOP
%tic
for ii=1:it; %count off iterations (epochs)

		
	%DECIDE WHO'S GOING TO FIRE THIS TIME
	
	%first remember who fired last time, so we can later find the correlations
	Bprev=B;
	
	%take present activation level divided by threshold, i.e. A./abs(T)
	%The following albegra is equivalent to a heaviside function centered at 1.
	At=0.5*(sign(A).*A+A);
	B=ceil(0.5.*(  ones(N)  +  sign( (At./abs(T))-ones(N) )  ));
	
	%toss in the periodic firing pattern if it's time
	if (p>0)*(ii<pend); if rem(ii,p)==1; B=B+P; end; end;  %add periodic input, if any
	if (ii==pend); %pend is the iteration on which a one-time event occurs, ending period input
		B=B+Pend;  %Pend is the activity pattern to add in at this time
	end;
	

	
	%CALCULATE NEW ACTIVATION LEVELS
	
	%remember activity from last time
	Aprev=A;
	
	%make new random inputs  -- to change, comment/uncomment "R=..."
	%R=a*rand(N);							%uniform distribution
	%R=zeros(N);							%no random input at all!
	k=50; R=a*(exp(-k))*exp(k*rand(N)); 	%k changes the convexity of the
					%exponential curve (but not maximum height)
					%1: almost no convexity;  4: medium;  20: very convex!
					%(I usually leave it at 50)
					
	%account for the contribution of firing cells
	cnv=zeros(N);  %cnv will be the total contribution of cells that fired
	for cc=1:N; for dd=1:N;  %go through every cell in the NxN grid...
		if B(cc,dd)>0;    %... but only worry about the cells that fired!
			%add its activation to other cells, according to the connectivity matric.
			here=zeros(N);  here(cc,dd)=1;
			cnv=cnv+conv2(here,C(:,:,cc,dd),'same');
		end;
	end; end;
	
	%decay the old activation, add the incoming contributions of cells that fired,
	%and throw in the random activations
	A=d.*A + cnv + R;
	
	
	
	%UPDATE THRESHOLDS AND ACTIVATION
	
	%update thresholds to place those who just fired in the refractory state
	T=T+100*(2.^(round(t.*rand(N)))).*B;
	%the real part is shot into the stratosphere (to a random height)
	
	%remove one layer of delay
	T=T-0.5.*real(T)-0.5.*(rem(real(T),100));
	%the real part descends by one half each epoch
	
	%zero the activation of anyone who is in the refractory period
	A=A.*(ones(N)-sign(abs(real(T))));
	%sodium won't go in if potassium's already out!
	
	%rectify
	A=0.5*(sign(A).*A+A);
	%excess negative activity cannot be tolerated!  This algebra makes any
	%cells below l equal to l.  Currently only supported for l = 0.
	
	
	
	%LEARN, BEAST!
	
	%make some padded matrices, to make the computations easier
	Bpadded=zeros(N+2*c); Bpadded(c+1:N+c,c+1:N+c)=B(1:N,1:N); %just B plus a border	
	Apadded=zeros(N+2*c); Apadded(c+1:N+c,c+1:N+c)=A(1:N,1:N); %just A plus a border
	Appadded=zeros(N+2*c); Appadded(c+1:N+c,c+1:N+c)=Aprev(1:N,1:N); %just Aprev plus a border	
	Tpadded=zeros(N+2*c); Tpadded(c+1:N+c,c+1:N+c)=sign(real(T(1:N,1:N)));
	%a border around this NxN matrix: 1 = refractory, 0 = not refractory
	
	%update correlations (every epoch)
	for cc=1:N; for dd=1:N;  %go through every cell in the NxN grid...
		if Bprev(cc,dd)>0;  %...and if it fired last time
			
			%(these calculations are performed on the whole connectivity area at once
			%to save processing time)
			
			%tally one additional positive correlation
			Hact(:,:,cc,dd)=Hact(:,:,cc,dd)+Bpadded(cc:cc+2*c,dd:dd+2*c);
			
			%for successful inhibitions, do something more complicated
			%a cell in the connectivity radius is "successfully inhibitted" if:
			%	1) its activation decreased by more than just the decay multiplier
			%	2) that decrease in activation was NOT due to firing
			%		(recall that firing zeros the activation level)
			%among "successfully inhibitted" cells, distinguish several magnitudes of
			%success based on the proportion by which activation decreased:
			%(d is the proportion by which the cell automatically decays each epoch)
			%	above  90% of d: do not change the successful inhibition tally
			%	90% to 71% of d: add 1 to the successful inhibition tally
			%	70% to 21% of d: add 2
			%	20% to  1% of d: add 4
			Hinh(:,:,cc,dd)=Hinh(:,:,cc,dd);
			
			%these tallies, of positive correlations and successful inhibitions,
			%are periocially used to update the weights in the
			%connectivity matrix.  (period = h)
			
		end;
	end; end;
	
	%update weights (every h epochs)
	if rem(ii,h)==0;
		for cc=1:N; for dd=1:N; %go through every cell
			%for 
			C(:,:,cc,dd)=Updateact(round(Hact(:,:,cc,dd)+ones(2*c+1))).*C(:,:,cc,dd);
			
			%clear the record of positive and negative correlations
			Hact(:,:,cc,dd)=zeros(2*c+1);
			Hinh(:,:,cc,dd)=zeros(2*c+1);
		end;end;
	end;

	
		
	%MAKE HISTOGRAM OF CLUSTER SIZES
	
	individualdisplay=1; %shows each cluster in diff color
	if (individualdisplay); counter=1; disp=zeros(N); end;
	done=ones(N); %(1 = not done, 0 = done)
	%conn=[0 1 0;1 1 1;0 1 0]; %nearest neighbor
	conn=[1 1 1;1 1 1;1 1 1];  %including diagonals
	maxclustersize=1;
	clear histo; histo(1)=0;
	
	for cc=1:N; for dd=1:N; %cycle through forest
		if B(cc,dd).*done(cc,dd)==1;
			%burn down the cluster
			burner=zeros(N); burner(cc,dd)=1;
			burnerold=ones(N);
			while sum(sum(burner ~= burnerold));
				burnerold=burner;
				burner=ceil(B.*conv2(burner,conn,'same')/10);
			end; %burner loop
		
			%fix what's left
			done=done-burner;
			if (individualdisplay); disp=disp+counter*burner; counter=counter+1; end;
			if(0)
				figure(7)
				image(100*disp/(max(max(disp))+.1))
				pause(0.1)
			end;
			clustersize=sum(sum(burner));
			if clustersize>maxclustersize;  %if it's the first time
				histo(clustersize)=1;
			else;  %if it's not the first time
				histo(clustersize)=histo(clustersize)+1;
			end;
			maxclustersize=max(maxclustersize,clustersize);
		end;  %if this square is on fire
	end;end;   %cycle through forest
	histhist(ii,1:maxclustersize)=histo;

	
	
	%DRAW IT
	
	if rem(ii,v)==0;
		tic
		%figure(1)
		%image(100*B);  %Display burning matrix
		%figure(2)
		%image(20*imag(T)+80*real(T)/(100*2^(t/2)));  %The pretty thresholds matrix
		figure(4)
		colormap(gray)
		image(50+50*10*cnv)
		figure(3)
		%image(30+20*A+50*B-.35*(20*imag(T)+80*real(T)/(100*2^(t/2))));
		image(30+20*A+50*B-100*real(T)/(100*2^(t/2)));
		figure(5)
		disppadded(2:2*c+2,2:2*c+2)=50+50*C(:,:,25,15);
		image(disppadded)
		figure(1)
		image(50*B(25,15))
		%figure(5)
		%image(100*R);  %The random activations (boring!)
		figure(2)
		image(50+10*Hact(:,:,25,15)-5*Hinh(:,:,25,15))
		%Hinh(:,:,25,15)
		
		%Plot histogram of cluster sizes
		figure(7)
		image(100*disp/(max(max(disp))+.1))
		pause(0.1)
		
		figure(6)
		if maxclustersize<40;
			histo(40)=0;
			plot(1:40,histo);
		else;
			plot(1:maxclustersize,histo)
		end;
	
		figure(8)
		image(8*histhist(max(1,ii-histrecall):max(histrecall,ii),:))
		
		%for kk=0:N-1;
		%	grammer(1+kk*N:N*kk+N)=disp(1:N,kk+1);
		%end;
		%figure(8)
		%hist(grammer)
		
		figure(10)
		image(35+30*Apadded(25:25+2*c,15:15+2*c))
		figure(11)
		image(35+30*Appadded(25:25+2*c,15:15+2*c))
		figure(12)
		quot=Apadded(25:25+2*c,15:15+2*c)./(Appadded(25:25+2*c,15:15+2*c)+0.00000000000001*ones(2*c+1));
		image(50+20*quot)
		figure(13)
		total=2*ones(2*c+1)-0.5*sign(quot-.9*d*ones(2*c+1))-0.5*sign(quot-.7*d*ones(2*c+1))-1*sign(quot-.2*d*ones(2*c+1));
		image(100-20*total-100*Tpadded(25:25+2*c,15:15+2*c))

		
		pause(0.0000001)
		DrawingTimeForIteration=ii
		toc
	end; %DRAW IT

end;

%toc


%PUT UP ALL THE HISTOGRAMS
figure(16)
colormap(copper)
image(10*histhist)


%DISPLAY THE CONNECTIVITY MATRICES AT SOME POINTS

if 0;
	tempcount=zeros(N);
	for cc=1:rem(N,2*c+1)+c; for dd=1:rem(N,2*c+1)+c; %cycle lattice definitions
		tempdisp=zeros(N+2*c);
		tempdisp(cc,dd)=1;tempdisp(cc+floor(N/(2*c+1))*(2*c+1),dd+floor(N/(2*c+1))*(2*c+1))=1;
		for jj=0:floor((N-cc)/(2*c+1));for kk=0:floor((N-dd)/(2*c+1)); %cycle points of this lattice
			temp=zeros(N+2*c);
			tempcount(cc+(2*c+1)*jj,dd+(2*c+1)*kk)=1+tempcount(cc+(2*c+1)*jj,dd+(2*c+1)*kk);
			temp(cc+(2*c+1)*jj+c,dd+(2*c+1)*kk+c)=1;
			tempdisp=tempdisp+conv2(temp,C(:,:,cc+(2*c+1)*jj,dd+(2*c+1)*kk),'same');
		end;end;
		figure(3)
		image(25+50*tempdisp(c+1:N+c,c+1:N+c));
		pause
		pause(0.01)
	end;end;
end;

