V1_ReverseCorrelation/V1_reverse_correl.m

function [w,rcImage,rcFreq,rcOrient,rcBlank,rcJointOrient,fcStim]=V1_reverse_correl()
sigmaV1=4;
sigmaLGN=1;
plotResults=1;
plotNodes=[1:5,37];
maxTime=10;

%imageType='sparse';%'hartley';%'dense';%'gratings';%
imageType='RHS2003';%'BR2002';%'Felsen02';%'CDL2005';%'RHS1997';%
imageBorder=13; %number of blank pixels surrounding image
angle_correl=0;
freq_correl=0;
ave_activity=0;
norm_correl=0;

%define v1 prediction node receptive fields
w=filter_definitions_V1_simple(sigmaV1);

%CREATE THE IMAGE DATA SET
disp('generating image dataset')
switch imageType
 
 case 'dense'     %The need to store all images makes this method impractical!
  numTrials=1;    %number of times each image is presented
  imageSize=32;   %size, in pixels, of each image
  gridSize=16;
  numImages=20000;
  presentationTime=3; %number of iterations for which each image is presented each time
  Iraw=images_dense_noise(imageSize,imageBorder,gridSize,numImages);
 
 case 'msequence' %The need to store all images makes this method impractical!
				  %Cleverer code would avoid this. 
				  %to calculate correlation, at a time lag of zero only, use
				  %V1_reverse_correl_m_seq.m
  numTrials=1;
  presentationTime=3;
  gridSize=16;
  Iraw=images_m_sequence(imageSize,imageBorder,gridSize);
 
 case 'sparse'
  numTrials=15; %number of times each image is presented
  imageSize=32; %size, in pixels, of each image
  gridSize=16;
  barLen=1
  % imageSize=48 %size, in pixels, of each image
  % gridSize=12
  presentationTime=3; %number of iterations for which each image is presented each time
  Iraw=images_sparse_noise(imageSize,imageBorder,gridSize,barLen);
  [a,b,numImages]=size(Iraw);
  ave_activity=1;
 
 case 'hartley'
  numTrials=15; %number of times each image is presented
  imageSize=32; %size, in pixels, of each image
  presentationTime=3; %number of iterations for which each image is presented each time
  Iraw=images_hartley_subspace(imageSize,imageBorder,floor(imageSize/4),1); 

 case 'RHS1997'
  numTrials=50; %number of times each image is presented
  imageSize=48; %size, in pixels, of each image (approx 8*grating wavelength)
  presentationTime=2; %number of iterations for which each image is presented each time
  angle_correl=1;
  angleStep=10;
  angle_range=[0-angleStep/2,180-angleStep/2];  
  num_bins=18;
  k=0;
  for phaseVal=0:90:270;
	for angleVal=0:angleStep:180-angleStep;
	  k=k+1;
	  sf(k)=0.17;
	  angle(k)=angleVal;
	  phase(k)=phaseVal;
	  Iraw(:,:,k)=image_circular_grating(imageSize,imageBorder,1./sf(k),angle(k),phase(k),1);
	end
  end
 
 case 'CDL2005'
  numTrials=10; %number of times each image is presented
  %imageSize=15 %inner diameter of annular grating generating zero response
  imageSize=60 %four times 
  imageSize=7
  presentationTime=2; %number of iterations for which each image is presented each time
  angle_correl=1;
  angleStep=15;
  angle_range=[0-angleStep/2,180-angleStep/2];  
  num_bins=12;
  norm_correl=1;
  k=0;
  for numSets=1:10
	for phaseVal=-135:45:180;
	  for angleVal=0:angleStep:180-angleStep;
		k=k+1;
		sf(k)=0.17;
		angle(k)=angleVal;
		phase(k)=phaseVal;
		Iraw(:,:,k)=image_circular_grating(imageSize,imageBorder,1./sf(k),angle(k),phase(k),1);
	  end
	end
  end
  
 case 'RHS2003'
  numTrials=15 %number of times each image is presented
  imageSize=13 %size, in pixels, of each image (approx =peak of SF)
  %imageSize=45 %size, in pixels, of each image (approx =3.5*peak of SF)
  presentationTime=2; %number of iterations for which each image is presented each time
  angle_correl=1;
  norm_correl=1;
  angleStep=10;
  angle_range=[0-angleStep/2,180-angleStep/2];  
  num_bins=18;
  k=0;
  for numSets=1:10
	for phaseVal=0:45:315;
	  for angleVal=0:angleStep:180-angleStep;
		k=k+1;
		sf(k)=0.17;
		angle(k)=angleVal;
		phase(k)=phaseVal;
		Iraw(:,:,k)=image_circular_grating(imageSize,imageBorder,1./sf(k),angle(k),phase(k),1);
	  end
	end
	for blank=1:8
	  k=k+1;
	  sf(k)=NaN;
	  angle(k)=NaN;
	  phase(k)=NaN;
	  Iraw(:,:,k)=zeros(imageSize+2*imageBorder)+0.5;
	end
  end

 case 'Felsen02'
  numTrials=25; %number of times each image is presented
  imageSize=45; %size, in pixels, of each image
  presentationTime=3; %number of iterations for which each image is presented each time
  angle_correl=1;
  angleStep=15;
  angle_range=[0-angleStep/2,180-angleStep/2];  
  num_bins=12;
  k=0;
  for numSets=1:40
	for phaseVal=0:90:270;
	  for angleVal=0:angleStep:180-angleStep;
		k=k+1;
		sf(k)=0.17;
		angle(k)=angleVal;
		phase(k)=phaseVal;
		Iraw(:,:,k)=image_circular_grating(imageSize,imageBorder,1./sf(k),angle(k),phase(k),1);
	  end
	end
  end
 
 case 'BR2002'
  numTrials=15 %number of times each image is presented
  imageSize=29 %size, in pixels, of each image 
  presentationTime=3; %number of iterations for which each image is presented each time
  freq_correl=1;
  norm_correl=1;
  octaves=5;
  freq_range=((2.^[1/octaves:1/octaves:octaves])/2^(octaves+1));
  num_bins=length(freq_range);
  k=0;
  for numSets=1:10
	for phaseVal=0:45:315;
	  for sfVal=freq_range
		k=k+1;
		sf(k)=sfVal;;
		angle(k)=90;
		phase(k)=phaseVal;
		Iraw(:,:,k)=image_circular_grating(imageSize,imageBorder,1./sf(k),angle(k),phase(k),1);
	  end
	end
	for blank=1:length(freq_range);
	  k=k+1;
	  sf(k)=NaN;
	  angle(k)=NaN;
	  phase(k)=NaN;
	  Iraw(:,:,k)=zeros(imageSize+2*imageBorder)+0.5;
	end
  end
end
Iraw=single(Iraw);
disp(imageType);
[a,b,numImages]=size(Iraw)
Isequence=zeros(a,b,numImages,'single');
binned_sf_seq=zeros(1,numImages*presentationTime,'single');
binned_angle_seq=zeros(1,numImages*presentationTime,'single');
if angle_correl
  [binned_angle,bin_labels_angle]=binned_param(angle,angle_range,num_bins);
end
if freq_correl
  [binned_sf]=binned_param(log(sf),log(freq_range),num_bins);
  bin_labels_sf=freq_range;
end
%PERFORM REVERSE CORRELATION EXPERIMENT
rcImage=[];
rcFreq=[];
rcOrient=[];
rcBlank=[];
rcJointOrient=[];
fcStim=[];
y=[];
for trial=1:numTrials
  trial
  disp('generating image sequence')
  %compile images into a sequence where each image is shown for in random order
  %for the specified number of iterations
  presentationOrder=randperm(numImages);
  t=0;
  for k=1:numImages
	for i=1:presentationTime
	  t=t+1;
	  Isequence(:,:,t)=Iraw(:,:,presentationOrder(k));
	  if angle_correl
		binned_angle_seq(t)=binned_angle(presentationOrder(k));
	  end
	  if freq_correl
		binned_sf_seq(t)=binned_sf(presentationOrder(k));
	  end
	end
  end
  [a,b,iterations]=size(Isequence);

  %pre-process image sequence to generate input to V1 model
  X=preprocess_V1_input(Isequence,sigmaLGN);
  
  disp('recording responses')
  %Present images to V1 model
  [y,r,e,ytrace]=dim_activation_conv(w,X,iterations,y);
  numSimple=length(ytrace);
  %calculate the response of complex cells by taking the maximum over the
  %outputs of the 4 simple cells with same orientation
  for complex=1:8
	ytrace{numSimple+complex}=calc_resp_complex(ytrace,1,1,complex+[0:8:31]);
  end

  disp('calculating reverse correlations for trial');
  %spatial reverse correlation 
  rcImage_this_trial=calc_reverse_correl_image(Isequence,ytrace,maxTime);
  %take mean of correlation across all trials
  rcImage=mean_over_trials(trial,rcImage,rcImage_this_trial);

  if angle_correl
	%orientation and spatial frequence reverse correlation
	[rcOrient_this_trial,rcBlank_this_trial]=...
		calc_reverse_correl_param(binned_angle_seq,num_bins,ytrace,maxTime,norm_correl);
	rcJointOrient_this_trial=...
		calc_reverse_correl_joint_param(binned_angle_seq,num_bins,ytrace,maxTime);
	
	%take mean of correlation across all trials
	rcOrient=mean_over_trials(trial,rcOrient,rcOrient_this_trial);
	rcBlank=mean_over_trials(trial,rcBlank,rcBlank_this_trial);
	rcJointOrient=mean_over_trials(trial,rcJointOrient,rcJointOrient_this_trial);
  end
  if freq_correl
	%orientation and spatial frequence reverse correlation
	[rcFreq_this_trial,rcBlank_this_trial]=...
		calc_reverse_correl_param(binned_sf_seq,num_bins,ytrace,maxTime,norm_correl);

	%take mean of correlation across all trials
	rcFreq=mean_over_trials(trial,rcFreq,rcFreq_this_trial);
	rcBlank=mean_over_trials(trial,rcBlank,rcBlank_this_trial);
  end
  
  if ave_activity
	%stimulus identity forward correlation 
	fcStim_this_trial=...
		calc_average_activity(presentationOrder,presentationTime,ytrace,maxTime);
	%take mean of correlation across all trials
	fcStim=mean_over_trials(trial,fcStim,fcStim_this_trial);
  end


  if plotResults || trial==numTrials
	disp('plotting reverse correlations')
	figure(1), clf
	plot_reverse_correl_image(w,rcImage,imageBorder,plotNodes);
	if angle_correl
	  figure(2), clf
	  plot_reverse_correl_orient(w,rcOrient,bin_labels_angle,plotNodes);
	  %figure(3), clf
	  %plot_reverse_correl_joint_orient(w,rcJointOrient,bin_labels_angle,plotNodes);
	end
	if freq_correl
	  figure(2), clf
	  plot_reverse_correl_freq(w,rcFreq,bin_labels_sf,plotNodes);
	end
  end
end



function rcImage=calc_reverse_correl_image(Isequence,ytrace,maxTime)
[a,b,iterations]=size(Isequence);
Isequence=Isequence-0.5; %subtrace mean luminance from images
%Isequence=abs(Isequence); %for complex cell RF mapping
for node=1:length(ytrace)
  %for each node with an RF centred at the middle of the image 
  resp(1:iterations)=ytrace{node};
  for tOffset=-maxTime:maxTime
	%calculate correlation between response and input image at fixed offset times
	i=max(1-tOffset,1)-1;
	for t=max(1+tOffset,1):min(iterations+tOffset,iterations)
	  i=i+1;
	  %correlation between response at time t and image at time i
	  rcImage_comp(:,:,i)=resp(t).*Isequence(:,:,i); 
	end
	%mean correlation over all images
	rcImage{node}(:,:,tOffset+maxTime+1)=mean(rcImage_comp,3);
  end
end



function [correl,blank_correl]=calc_reverse_correl_param(binned_param_seq,num_bins,ytrace,maxTime,norm)
if nargin<5, norm=0; end
[iterations]=length(binned_param_seq);
numNodes=length(ytrace);
numTimes=2*maxTime+1;
normOccur=zeros(num_bins,numTimes);
normOccurBlank=zeros(numTimes);
%normResp=zeros(numNodes,numTimes);

for node=1:numNodes
  %for each node with an RF centred at the middle of the image 
  resp(1:iterations)=ytrace{node};
  correl{node}(1:num_bins,1:numTimes)=0;
  blank_correl{node}(1:numTimes)=0;
  for tOffset=-maxTime:maxTime
	%calculate correlation between response and the parameter at previous times
	time=tOffset+maxTime+1;
	i=max(1-tOffset,1)-1;
	for t=max(1+tOffset,1):min(iterations+tOffset,iterations)
	  i=i+1;
	  bin=binned_param_seq(i);
	  %cummulative correlation between response at time t and parameter at time i
	  if isnan(bin)
		blank_correl{node}(time)=blank_correl{node}(time)+resp(t);
		if node==1, normOccurBlank(time)=normOccurBlank(time)+1;, end
	  elseif bin>0
		correl{node}(bin,time)=correl{node}(bin,time)+resp(t);
		if node==1,	normOccur(bin,time)=normOccur(bin,time)+1; end
	  end
	  %normResp(node,time)=normResp(node,time)+resp(t);
	end
  end
end
%normalise by total response from the node - can do this during post processing
%for node=1:numNodes
%  for tOffset=-maxTime:maxTime
%	time=tOffset+maxTime+1;
%	correl{node}(:,time)=correl{node}(:,time)./normResp(node,time);
%  end
%end
if norm
  %normalise by total number of times each parameter was shown
  for node=1:numNodes
	for tOffset=-maxTime:maxTime
	  time=tOffset+maxTime+1;
	  correl{node}(:,time)=correl{node}(:,time)./normOccur(:,time);
	  blank_correl{node}(time)=blank_correl{node}(time)./normOccurBlank(time);
	 % disp(['node=',int2str(node),' time=',int2str(time),...
	%		' normBlank=',int2str(normOccurBlank(time)),...
	%		' sum(normOccur)=',int2str(sum(normOccur(:,time)')),...
	%		' blankCorrel=',num2str(blank_correl{node}(time))]);
	end
  end
end



function jcorrel=calc_reverse_correl_joint_param(binned_param_seq,num_bins,ytrace,maxTime)
[iterations]=length(binned_param_seq);
numNodes=length(ytrace);
numTimes=2*maxTime+1;
normOccur=zeros(num_bins,num_bins,numTimes)+1e-9;

for node=1:numNodes
  %for each node with an RF centred at the middle of the image 
  resp(1:iterations)=ytrace{node};
  %resp=max(0,resp-mean(resp));
  jcorrel{node}(1:num_bins,1:num_bins,1:numTimes)=0;
  for tOffset=-maxTime:maxTime
	%calculate correlation between response and the parameter at previous times
	time=tOffset+maxTime+1;
	i=max(1-tOffset,1)-1;
	for t=max(1+tOffset,1):min(iterations+tOffset-1,iterations)
	  i=i+1;
	  param1=binned_param_seq(i);
	  param2=binned_param_seq(i+1);
	  %cummulative correlation between response at time t and parameter pair at
      %time i and i+1
	  if param1>0 & param2>0 %%& param1~=param2
		jcorrel{node}(param1,param2,time)=...
			jcorrel{node}(param1,param2,time)+resp(t);
		if node==1, 
		  normOccur(param1,param2,time)=normOccur(param1,param2,time)+1; 
		end
	  end
	end
  end
end
%normalise by total number of times each parameter pair was shown
for node=1:numNodes
  for tOffset=-maxTime:maxTime
	time=tOffset+maxTime+1;
	jcorrel{node}(:,:,time)=jcorrel{node}(:,:,time)./normOccur(:,:,time);
  end
end



function fcStim=calc_average_activity(presentationOrder,presentationTime,ytrace,maxTime)
numNodes=length(ytrace);
iterations=length(ytrace{1});
%yMUA=mean(cat(1,ytrace{:}));%the mean activity of all nodes across time

for node=1:numNodes
  fcStim{node}=zeros(max(presentationOrder),1+2*maxTime);

  for tOffset=1+maxTime:iterations-maxTime
	fcStim{node}(presentationOrder(ceil(tOffset/presentationTime)),:)=...
		ytrace{node}(tOffset-maxTime:tOffset+maxTime);
  end
end




function meanData=mean_over_trials(trial,meanData,newData)
fac=1/trial;
for node=1:length(newData)
  if trial==1
	meanData{node}=fac.*newData{node};
  else
	meanData{node}=fac.*newData{node}+(1-fac).*meanData{node};
  end  
end