%ben-ishai net
%version with weight matrix

%close all;

clear all;
n  = 100;
dth= pi/n;

% parameters
lam0 = 5; % >= 0
lam1 = 2.01; % >= 0
c = 2.5;
a = 50;
eps = 0.001; % 0 ... 0.5

% create input
h=zeros(1,n);
h2=h;
phi0 = pi/4;
dphi= pi/2;

h	= a*c*(1-eps+eps*cos(2*((1:n)*dth+phi0)) );
h2	= a*c*(1-eps+eps*cos(2*((1:n)*dth+phi0+dphi)) ); % h2 is shifted wrt h

% create weight matrix
hmat = zeros(n,n);
th_ar = (-pi/2:pi/n:pi/2-pi/n);
for i=1:n
    thp=(i-1)*dth-pi/2; 
    hmat(i,:)= (-lam0+lam1*cos(2*(th_ar-thp)))/n;
end

% the rates
v=ones(n,1);

nrepmax = 10000;
err_eps = 0.0003;
thr= 25;

peakloc_vs_t=[];
peakampl_vs_t=[];


for icase=1:2
	
	if (icase==2)
		h=h2;
	endif	

	irep =0;
	err=1;
	% iterate dynamics until err is small, or max £ reps is reached
	while (err > err_eps & irep < nrepmax)

	  	in = hmat*v + h';
	   	in = max(in-thr,0); % rectifty

	   	vold = v;
	   	v = v+0.2*(-v+in);     % integrate

	   	err = mean(abs(vold-v));%/mean(v);
	   	irep = irep+1;
	 	%  plot (v);  hold on;
		[m, im]= max(v);
   		peakloc_vs_t =  [peakloc_vs_t im];
		peakampl_vs_t = [peakampl_vs_t m];
	end
	if (irep == nrepmax)
    	'nrepmax reached !'
	end

	plot(v);
	hold on;
end

figure(2);
hold off
plot (peakloc_vs_t)

figure(3);
hold off
plot (peakampl_vs_t)


