clear all;
close all;
rand('state',2);

ntrials=50000;


xsize=128;
ysize=128;
freq = 2;
lambda = xsize/freq;
dotdens=0.25;
dotsize=round(0.05*lambda);
nbw = 2.5;
orbw = pi/180*30;
sx = sqrt(log(2)) * (lambda/2/pi) * ( 2.^nbw + 1) ./ ( 2.^nbw - 1) ;
sy = sqrt(log(2)) * (lambda/2/pi) * 2 / orbw ;


theta=0;
disparity=[-21:3:21];
halfdisparity=disparity/2
halfmaxdisparity=max(abs(disparity))/2
% THings for generating images
halfdot = fix(dotsize/2);
bigxsize = xsize+max(abs(disparity));
ndots = round(dotdens * (bigxsize+2*dotsize) * ysize / min(dotsize,xsize) / min(dotsize,ysize));

CORR=1;ANTI=2;
EVEN1=1;EVEN2=2;EVEN3=3;EVEN4=4;ODD1=5;ODD2=6;UNCORR=7;


x=[-xsize/2:xsize/2-1];
y=[-ysize/2:ysize/2-1];
[x2d,y2d]=meshgrid(x,y);


[tmp,nd]=size(disparity);

col={'b' 'r--'};
name={'even' 'odd'};

complex = zeros(ntrials,nd,2,7);

RF1 = 1/ (sx*sy*2*pi) .* exp(-0.5.*x2d.^2./sx^2) .* exp(-0.5.*y2d.^2./sy^2).* cos(2*pi/lambda.*x2d + pi/4) ;
RF2 = 1/ (sx*sy*2*pi) .* exp(-0.5.*x2d.^2./sx^2) .* exp(-0.5.*y2d.^2./sy^2) .* cos(2*pi/lambda.*x2d - pi/4);


Pos = inline('(x-2.8e-2).*(x>2.8e-2)');

figure('Units','centimeters','Position',[2 2 14.8 8],'PaperUnits','centimeters','paperposition',[0 0 14.8 8],'PaperPositionMode','auto')
    fs=10
ms=4;
    
cnt=0;
nonzero=0;
for trial=1:ntrials
	trial

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% Make a single big image to pick out:
	bigimage = zeros(ysize,bigxsize);
	for j=1:ndots
		% Pick position for dot:
		xpos = fix(rand(1)*(bigxsize+2*dotsize)) - halfdot;
		ypos = fix(rand(1)*ysize) - halfdot	;
		x1 = max(xpos,1);
		x2 = min(xpos + dotsize-1,bigxsize);
		y1 = max(ypos,1);
		y2 = min(ypos + dotsize-1,ysize);
		% Choose colour of dot (black or white)
		bigimage(y1:y2,x1:x2 ) = ((rand(1)>0.5) - 0.5 ) *2 .* ones(y2-y1+1,x2-x1+1);
	end
	% Make an unrealted image for uncorrelated case
	bigimage2 = zeros(ysize,bigxsize);
	for j=1:ndots
		% Pick position for dot:
		xpos = fix(rand(1)*(bigxsize+2*dotsize)) - halfdot;
		ypos = fix(rand(1)*ysize) - halfdot	;
		x1 = max(xpos,1);
		x2 = min(xpos + dotsize-1,bigxsize);
		y1 = max(ypos,1);
		y2 = min(ypos + dotsize-1,ysize);
		% Choose colour of dot (black or white)
		bigimage2(y1:y2,x1:x2 ) = ((rand(1)>0.5) - 0.5 ) *2 .* ones(y2-y1+1,x2-x1+1);
	end

	for j=1:nd
		% Now pick out images of the right size
		imL = bigimage(:,1+halfmaxdisparity-halfdisparity(j):xsize+halfmaxdisparity-halfdisparity(j));
		imR = bigimage(:,1+halfmaxdisparity+halfdisparity(j):xsize+halfmaxdisparity+halfdisparity(j));

		convL1 = sum(sum(RF1.*imL)); convL2 = sum(sum(RF2.*imL));
		convR1 = sum(sum(RF1.*imR)); convR2 = sum(sum(RF2.*imR));

		cnt=cnt+8;
		if Pos(convL1)>0; nonzero=nonzero+1;end
		if Pos(convL2)>0; nonzero=nonzero+1;end
		if Pos(convR1)>0; nonzero=nonzero+1;end
		if Pos(convR2)>0; nonzero=nonzero+1;end
		if Pos(-convL1)>0; nonzero=nonzero+1;end
		if Pos(-convL2)>0; nonzero=nonzero+1;end
		if Pos(-convR1)>0; nonzero=nonzero+1;end
		if Pos(-convR2)>0; nonzero=nonzero+1;end


		% Do odd/.even, correlated/anticorrelated:
		mono = Pos(convL1)^2+Pos(convR1)^2+Pos(-convL1)^2+Pos(-convR1)^2+Pos(convL2)^2+Pos(convR2)^2+Pos(-convL2)^2+Pos(-convR2)^2;
		halfcrossE1s = Pos(convL1)*Pos(convR1)+Pos(-convL1)*Pos(-convR1)+Pos(convL2)*Pos(convR2)+Pos(-convL2)*Pos(-convR2);
		halfcrossE1a = Pos(convL1)*Pos(-convR1)+Pos(-convL1)*Pos(convR1)+Pos(convL2)*Pos(-convR2)+Pos(-convL2)*Pos(convR2);
		complex(trial,j,CORR,EVEN1) = mono + 2*halfcrossE1s; %...
		%  					  ( Pos(convL1)+Pos(convR1) )^2 + ( Pos(-convL1)+Pos(-convR1) )^2 ...
		% 					+ ( Pos(convL2)+Pos(convR2) )^2 + ( Pos(-convL2)+Pos(-convR2) )^2 ;
		complex(trial,j,ANTI,EVEN1) = mono + 2*halfcrossE1a; %. ...
		%  					  ( Pos(convL1)+Pos(-convR1) )^2 + ( Pos(-convL1)+Pos(convR1) )^2 ...
		% 					+ ( Pos(convL2)+Pos(-convR2) )^2 + ( Pos(-convL2)+Pos(convR2) )^2 ;
		complex(trial,j,CORR,EVEN2) = mono - 2*halfcrossE1s; % ...
		%  					  ( Pos(convL1)-Pos(convR1) )^2 + ( Pos(-convL1)-Pos(-convR1) )^2 ...
		% 					+ ( Pos(convL2)-Pos(convR2) )^2 + ( Pos(-convL2)-Pos(-convR2) )^2 ;
		complex(trial,j,ANTI,EVEN2) = mono - 2*halfcrossE1a; % ...
		%  					  ( Pos(convL1)-Pos(-convR1) )^2 + ( Pos(-convL1)-Pos(convR1) )^2 ...
		% 					+ ( Pos(convL2)-Pos(-convR2) )^2 + ( Pos(-convL2)-Pos(convR2) )^2 ;
		halfcrossE3s = Pos(convL1)*Pos(convR2)+Pos(-convL1)*Pos(-convR2)+Pos(convL2)*Pos(convR1)+Pos(-convL2)*Pos(-convR1);
		halfcrossE3a = Pos(convL1)*Pos(-convR2)+Pos(-convL1)*Pos(convR2)+Pos(convL2)*Pos(-convR1)+Pos(-convL2)*Pos(convR1);
		complex(trial,j,CORR,EVEN3) = mono + 2*halfcrossE3s; % ...
		% 					  ( Pos(convL1)+Pos(convR2) )^2 + ( Pos(-convL1)+Pos(-convR2) )^2 ...
		% 					+ ( Pos(convL2)+Pos(convR1) )^2 + ( Pos(-convL2)+Pos(-convR1) )^2 ;
		complex(trial,j,ANTI,EVEN3) = mono + 2*halfcrossE3a; % ...
		%  					  ( Pos(convL1)+Pos(-convR2) )^2 + ( Pos(-convL1)+Pos(convR2) )^2 ...
		% 					+ ( Pos(convL2)+Pos(-convR1) )^2 + ( Pos(-convL2)+Pos(convR1) )^2 ;
		complex(trial,j,CORR,EVEN4) = mono - 2*halfcrossE3s; % ...
		%  					  ( Pos(convL1)-Pos(convR2) )^2 + ( Pos(-convL1)-Pos(-convR2) )^2 ...
		% 					+ ( Pos(convL2)-Pos(convR1) )^2 + ( Pos(-convL2)-Pos(-convR1) )^2 ;
		complex(trial,j,ANTI,EVEN4) = mono - 2*halfcrossE3a; % ...
		%				  ( Pos(convL1)-Pos(-convR2) )^2 + ( Pos(-convL1)-Pos(convR2) )^2 ...
		%				+ ( Pos(convL2)-Pos(-convR1) )^2 + ( Pos(-convL2)-Pos(convR1) )^2 ;
		halfcrossOs = -Pos(convL1)*Pos(convR2)-Pos(-convL1)*Pos(-convR2)+Pos(convL2)*Pos(convR1)+Pos(-convL2)*Pos(-convR1);
		halfcrossOa = -Pos(convL1)*Pos(-convR2)-Pos(-convL1)*Pos(convR2)+Pos(convL2)*Pos(-convR1)+Pos(-convL2)*Pos(convR1);
		complex(trial,j,CORR,ODD1)= mono + 2*halfcrossOs ; %...
		%                 	(Pos(convL1) - Pos(convR2))^2 + (Pos(-convL1) - Pos(-convR2))^2 ...
		%            	 + (Pos(convL2) + Pos(convR1))^2 + (Pos(-convL2) + Pos(-convR1))^2;
		complex(trial,j,ANTI,ODD1)= mono + 2*halfcrossOa ; %...
		%                	(Pos(convL1) - Pos(-convR2))^2 + (Pos(-convL1) - Pos(convR2))^2 ...
		%          	 + (Pos(convL2) + Pos(-convR1))^2 + (Pos(-convL2) + Pos(convR1))^2;
		complex(trial,j,CORR,ODD2)=  mono - 2*halfcrossOs ;%...
		%              	(Pos(convL1) + Pos(convR2))^2 + (Pos(-convL1) + Pos(-convR2))^2 ...
		%        	 + (Pos(convL2) - Pos(convR1))^2 + (Pos(-convL2) - Pos(-convR1))^2;
		complex(trial,j,ANTI,ODD2)= mono - 2*halfcrossOa ;% ...
		%            	(Pos(convL1) + Pos(-convR2))^2 + (Pos(-convL1) + Pos(convR2))^2 ...
		%      	 + (Pos(convL2) - Pos(-convR1))^2 + (Pos(-convL2) - Pos(convR1))^2;

		lmonoc(trial,j) = Pos(convL1)^2+Pos(-convL1)^2+Pos(convL2)^2+Pos(-convL2)^2;
		rmonoc(trial,j) = Pos(convR1)^2+Pos(-convR1)^2+Pos(convR2)^2+Pos(-convR2)^2;

		imR = bigimage2(:,1+halfmaxdisparity+halfdisparity(j):xsize+halfmaxdisparity+halfdisparity(j));
		convR1 = sum(sum(RF1.*imR)); convR2 = sum(sum(RF2.*imR));

		cnt = cnt+4;
		if Pos(convR1)>0; nonzero=nonzero+1; end
		if Pos(convR2)>0; nonzero=nonzero+1; end
		if Pos(-convR1)>0; nonzero=nonzero+1; end
		if Pos(-convR2)>0; nonzero=nonzero+1; end

		mono = Pos(convL1)^2+Pos(convR1)^2+Pos(-convL1)^2+Pos(-convR1)^2+Pos(convL2)^2+Pos(convR2)^2+Pos(-convL2)^2+Pos(-convR2)^2;
		halfcrossE1s = Pos(convL1)*Pos(convR1)+Pos(-convL1)*Pos(-convR1)+Pos(convL2)*Pos(convR2)+Pos(-convL2)*Pos(-convR2);
		halfcrossE1a = Pos(convL1)*Pos(-convR1)+Pos(-convL1)*Pos(convR1)+Pos(convL2)*Pos(-convR2)+Pos(-convL2)*Pos(convR2);
		complex(trial,j,UNCORR,EVEN1) = mono + 2*halfcrossE1s;
		complex(trial,j,UNCORR,EVEN1) = mono + 2*halfcrossE1a;
		complex(trial,j,UNCORR,EVEN2) = mono - 2*halfcrossE1s;
		complex(trial,j,UNCORR,EVEN2) = mono - 2*halfcrossE1a;
		halfcrossE3s = Pos(convL1)*Pos(convR2)+Pos(-convL1)*Pos(-convR2)+Pos(convL2)*Pos(convR1)+Pos(-convL2)*Pos(-convR1);
		halfcrossE3a = Pos(convL1)*Pos(-convR2)+Pos(-convL1)*Pos(convR2)+Pos(convL2)*Pos(-convR1)+Pos(-convL2)*Pos(convR1);
		complex(trial,j,UNCORR,EVEN3) = mono + 2*halfcrossE3s;
		complex(trial,j,UNCORR,EVEN3) = mono + 2*halfcrossE3a;
		complex(trial,j,UNCORR,EVEN4) = mono - 2*halfcrossE3s;
		complex(trial,j,UNCORR,EVEN4) = mono - 2*halfcrossE3a;
		halfcrossOs = -Pos(convL1)*Pos(convR2)-Pos(-convL1)*Pos(-convR2)+Pos(convL2)*Pos(convR1)+Pos(-convL2)*Pos(-convR1);
		halfcrossOa = -Pos(convL1)*Pos(-convR2)-Pos(-convL1)*Pos(convR2)+Pos(convL2)*Pos(-convR1)+Pos(-convL2)*Pos(convR1);
		complex(trial,j,UNCORR,ODD1)= mono + 2*halfcrossOs;
		complex(trial,j,UNCORR,ODD1)= mono + 2*halfcrossOa;
		complex(trial,j,UNCORR,ODD2)=  mono - 2*halfcrossOs;
		complex(trial,j,UNCORR,ODD2)= mono - 2*halfcrossOa;

	end % next disparity

	if fix(trial/500)  == trial/500
		Lmon = squeeze(mean(mean(lmonoc,2),1)) ;
		Rmon = squeeze(mean(mean(rmonoc,2),1)) ;
		unc = squeeze(mean(mean(complex(:,:,UNCORR,:),2),1)) ;

		disparity2 = [min(disparity) : 0.1 : max(disparity) ];
		subplot(2,3,1)
		curve = (squeeze(mean(complex(:,:,CORR,EVEN1),1))) ;
		mx=max(curve);
		amp1=max(curve)-min(curve);
		plot(disparity,curve/mx,'ko','MarkerFaceColor','k','markersize',ms);   hold on
		curve2 = spline(disparity,curve,disparity2);
		plot(disparity2,curve2/mx,'k')
		curve = (squeeze(mean(complex(:,:,ANTI,EVEN1),1))) ;
		amp2=max(curve)-min(curve);
		plot(disparity,curve/mx,'ko','markersize',ms);
		curve2 = spline(disparity,curve,disparity2);
		plot(disparity2,curve2/mx,'k:')
		ratio1 = amp2/amp1
		title('"Tuned excitatory"','Fontsize',fs)
		ylabel('firing rate','Fontsize',fs)
		axis tight
		plot([-73 73],[1 1]*unc(EVEN1)/mx,'k','clipping','off')
		text(77,unc(EVEN1)/mx,'U','fontsize',8)
		set(gca,'units','normalized','fontsize',8,'position',[0.1   0.58   0.23   0.34])

		subplot(2,3,4)
		curve = (squeeze(mean(complex(:,:,CORR,EVEN2),1))) ;
		mx=max(curve);
		amp1=max(curve)-min(curve);
		plot(disparity,curve/mx,'ko','MarkerFaceColor','k','markersize',ms);   hold on
		curve2 = spline(disparity,curve,disparity2);
		plot(disparity2,curve2/mx,'k')
		curve = (squeeze(mean(complex(:,:,ANTI,EVEN2),1))) ;
		amp2=max(curve)-min(curve);
		plot(disparity,curve/mx,'ko','markersize',ms);
		curve2 = spline(disparity,curve,disparity2);
		plot(disparity2,curve2/mx,'k:')
		ratio2 = amp2/amp1
		title('"Tuned inhibitory"','Fontsize',fs)
		ylabel('firing rate','Fontsize',fs)
		axis tight
		%plot([-73 73],[1 1]*Lmon/mx,'k','clipping','off')
		%plot([-73 73],[1 1]*Rmon/mx,'k','clipping','off')
		plot([-73 73],[1 1]*unc(EVEN2)/mx,'k','clipping','off')
		%text(77,Lmon/mx,'L','fontsize',8)
		%text(74,Rmon/mx,'R','fontsize',8)
		text(77,unc(EVEN2)/mx,'U','fontsize',8)
		set(gca,'units','normalized','fontsize',8,'position',[0.1   0.11   0.23   0.34])

		subplot(2,3,2)
		curve = (squeeze(mean(complex(:,:,CORR,ODD1),1))) ;
		mx=max(curve);
		amp1=max(curve)-min(curve);
		plot(disparity,curve/mx,'ko','MarkerFaceColor','k','markersize',ms);   hold on
		curve2 = spline(disparity,curve,disparity2);
		plot(disparity2,curve2/mx,'k')
		curve = (squeeze(mean(complex(:,:,ANTI,ODD1),1))) ;
		amp2=max(curve)-min(curve);
		plot(disparity,curve/mx,'ko','markersize',ms);
		curve2 = spline(disparity,curve,disparity2);
		plot(disparity2,curve2/mx,'k:')
		ratio3 = amp2/amp1
		title('"Near"','Fontsize',fs)
		axis tight
		plot([-73 73],[1 1]*unc(ODD1)/mx,'k','clipping','off')
		text(77,unc(ODD1)/mx,'U','fontsize',8)
		set(gca,'units','normalized','fontsize',8,'position',[0.41  0.58   0.23   0.34])

		subplot(2,3,5)
		curve = (squeeze(mean(complex(:,:,CORR,ODD2),1))) ;
		mx=max(curve);
		amp1=max(curve)-min(curve);
		plot(disparity,curve/mx,'ko','MarkerFaceColor','k','markersize',ms);   hold on
		curve2 = spline(disparity,curve,disparity2);
		plot(disparity2,curve2/mx,'k')
		curve = (squeeze(mean(complex(:,:,ANTI,ODD2),1))) ;
		amp2=max(curve)-min(curve);
		plot(disparity,curve/mx,'ko','markersize',ms);
		curve2 = spline(disparity,curve,disparity2);
		plot(disparity2,curve2/mx,'k:')
		ratio3 = amp2/amp1
		title('"Far"','Fontsize',fs)
		xlabel('disparity','Fontsize',fs)
		axis tight
		%plot([-73 73],[1 1]*Lmon/mx,'k','clipping','off')
		%plot([-73 73],[1 1]*Rmon/mx,'k','clipping','off')
		plot([-73 73],[1 1]*unc(ODD2)/mx,'k','clipping','off')
		%text(77,Lmon/mx,'L','fontsize',8)
		%text(60,Rmon/mx,'R','fontsize',8)
		text(77,unc(ODD2)/mx,'U','fontsize',8)
		set(gca,'units','normalized','fontsize',8,'position',[0.41  0.11   0.23   0.34])

		subplot(2,3,3)
		curve = (squeeze(mean(complex(:,:,CORR,EVEN3),1))) ;
		mx=max(curve);
		amp1=max(curve)-min(curve);
		plot(disparity,curve/mx,'ko','MarkerFaceColor','k','markersize',ms);   hold on
		curve2 = spline(disparity,curve,disparity2);
		plot(disparity2,curve2/mx,'k')
		curve = (squeeze(mean(complex(:,:,ANTI,EVEN3),1))) ;
		amp2=max(curve)-min(curve);
		plot(disparity,curve/mx,'ko','markersize',ms);
		curve2 = spline(disparity,curve,disparity2);
		plot(disparity2,curve2/mx,'k:')
		ratio1 = amp2/amp1
		title('"Notch"','Fontsize',fs)
		axis tight
		%plot([-73 73],[1 1]*Lmon/mx,'k','clipping','off')
		%plot([-73 73],[1 1]*Rmon/mx,'k','clipping','off')
		plot([-73 73],[1 1]*unc(EVEN3)/mx,'k','clipping','off')
		%text(77,Lmon/mx,'L','fontsize',8)
		%text(60,Rmon/mx,'R','fontsize',8)
		text(77,unc(EVEN3)/mx,'U','fontsize',8)
		set(gca,'units','normalized','fontsize',8,'position',[0.72  0.58   0.23   0.34])

		subplot(2,3,6)
		curve = (squeeze(mean(complex(:,:,CORR,EVEN4),1))) ;
		mx=max(curve);
		amp1=max(curve)-min(curve);
		plot(disparity,curve/mx,'ko','MarkerFaceColor','k','markersize',ms);   hold on
		curve2 = spline(disparity,curve,disparity2);
		plot(disparity2,curve2/mx,'k')
		curve = (squeeze(mean(complex(:,:,ANTI,EVEN4),1))) ;
		amp2=max(curve)-min(curve);
		plot(disparity,curve/mx,'ko','markersize',ms);
		curve2 = spline(disparity,curve,disparity2);
		plot(disparity2,curve2/mx,'k:')
		ratio2 = amp2/amp1
		title('"Notch"','Fontsize',fs)
		axis tight
		%plot([-73 73],[1 1]*Lmon/mx,'k','clipping','off')
		%plot([-73 73],[1 1]*Rmon/mx,'k','clipping','off')
		plot([-73 73],[1 1]*unc(EVEN4)/mx,'k','clipping','off')
		%text(77,Lmon/mx,'L','fontsize',8)
		%text(60,Rmon/mx,'R','fontsize',8)
		text(77,unc(EVEN4)/mx,'U','fontsize',8)
		set(gca,'units','normalized','fontsize',8,'position',[0.72   0.11   0.23   0.34])

		drawnow
	end


end % next trial

disp(['Convolutions were nonzero in ' num2str(nonzero/cnt*100) '% of cases'])
save SixExamples_Gabor_Thresh.mat



figure('Units','centimeters','Position',[2 2 8.8 8.8],'PaperUnits','centimeters','paperposition',[0 0 8.8 8.8],'PaperPositionMode','auto')
mx= max(max(RF1));
RF1 = RF1 / mx;
RF2 = RF2 / mx;
fs=10;
subplot(2,2,1)
imagesc(x,x,RF1)
ylabel('vertical position','Fontsize',fs)
axis square
subplot(2,2,3)
imagesc(x,x,RF2)
ylabel('vertical position','Fontsize',fs)
xlabel('horizontal position','Fontsize',fs);
axis square
subplot(2,2,2)
plot(x,RF1(64,:),'k');
set(gca,'xlim',[min(x) max(x)],'fontsize',fs)
%ylabel('RF function','Fontsize',fs)
subplot(2,2,4)
plot(x,RF2(64,:),'k');
%ylabel('RF function','Fontsize',fs)
xlabel('horizontal position','Fontsize',fs);
set(gca,'xlim',[min(x) max(x)],'fontsize',fs)
colormap(gray)


