V1_ReverseCorrelation/plot_reverse_correl_orient_logprob.m

function Apeak=plot_reverse_correl_orient_logprob(w,rcOrient,rcBlank,bin_labels,plotNodes,plotTimes,plotA,plotR)
colPlus=1; 
if nargin<5
  numNodes=length(rcOrient);
else
  numNodes=length(plotNodes);
end
[poo,numTimes]=size(rcOrient{1});
maxNodes=size(w,1);
maxTime=floor(numTimes/2);
if nargin<6
  plotTimes=[-maxTime:maxTime];
else
  numTimes=length(plotTimes);
end
if nargin<7, plotA=0; end
if nargin<8, plotR=0; end


numBins=size(rcOrient{1},1);
%binsToLabel=[1:floor(numBins/2):numBins+1];
binsToLabel=[floor(numBins/4),1+ceil(numBins/2),2+ceil(numBins*0.75)];

meanA=0;
meanRmax=0;
meanRmin=0;
meanRorth=0;
clf
nodeNum=0;
for node=plotNodes
  plotAnode=0;
  plotRnode=0;
  nodeNum=nodeNum+1;
  %for each node show the the reverse correlation at each time lag
  toff=0;
  maxData=0;
  for time=plotTimes
	toff=toff+1;

	data=rcOrient{node}(:,time+maxTime+1);
	%blank=rcBlank{node}(time+maxTime+1);
	%blank=mean(rcOrient{node}(:,time+maxTime+1));
	tblank=plotTimes(1)+maxTime+1;
	blank=rcOrient{node}(:,tblank);
	
	%calculate R values (the "tuning curve") and parameters describing this curve,
    %as used by Ringach et al
	R(:,toff)=log10(data./blank);
	[Rmax(toff),oriPref]=max(R(:,toff));
	Rmin(toff)=min(R(:,toff));
	oriOrth=mod(oriPref-numBins/2-1,numBins)+1;
	Rorth(toff)=R(oriOrth,toff);
	A(toff)=Rmax(toff)-Rmin(toff);
	
	if ~plotA & ~plotR
	  %plot tuning curve a current time step
	  maxsubplot(numNodes+1,colPlus+numTimes+2,(colPlus+numTimes+2)*(nodeNum)+toff+colPlus);
	  plot(R(:,toff),'LineWidth',2), axis([1,numBins+1,-0.15,0.5])
	  hold on, plot([0,numBins+1],[0,0],'k')
	  %plot([0,numBins],[1,1].*rcBlank{node}(time+maxTime+1),'g')
	  %plot(rcOrient{node}(:,time+maxTime+1),'r');
	  drawnow;
	  if nodeNum==1 & ~isempty(bin_labels)
		set(gca,'XTick',binsToLabel,'XTickLabel',num2str(bin_labels(binsToLabel)'));
	  else
		set(gca,'XTick',[]);
	  end
	  if toff==1
		set(gca,'YTick',[-0.25:0.25:1]);
		ylabel('OT(\tau,\theta)','FontSize',12)
	  else
		set(gca,'YTick',[]);
	  end
	  if nodeNum==1 & toff==ceil(numTimes/2)
		xlabel('Orientation (degrees)','FontSize',12)
	  end
	  set(gca,'FontSize',10);
	end
  end
  %plot tuning curve parameters through time
  if ~plotA & ~plotR
	maxsubplot(numNodes+1,colPlus+numTimes+2,(colPlus+numTimes+2)*(nodeNum)+toff+colPlus+1);
	plotAnode=1;
  end
  if plotA
	maxsubplot(numNodes+1,colPlus+2,(colPlus+2)*(nodeNum)+colPlus+1);
	plotAnode=1;
  end
  if plotAnode
	plot(A,'r','LineWidth',3), hold on, plot(Rmax,'k'); 
	plot(Rmin,'b--','LineWidth',3), plot(Rorth,'g','LineWidth',3);
	plot([0,numTimes+2],[0,0],'k')
	ax=axis; axis([1,numTimes,ax(3),ax(4)])
	if nodeNum==1, 
	  toLabel=unique([min(plotTimes),0,max(plotTimes)]);
	  set(gca,'XTick',unique(1+toLabel-min(plotTimes)),'XTickLabel',int2str(toLabel'))
	  xlabel('Time (iterations)','FontSize',12)
	else
	  set(gca,'XTick',[]);	
	end
	ylabel('OT_{min}(\tau), OT_{orth}(\tau), A(\tau)','FontSize',12)
	set(gca,'FontSize',10);
  end
  
  %find tuning curves at R(dev) and R(dec)
  [Amax,tmax]=max(A);
  indThres=find(A>0.5*Amax);
  tdev=indThres(1); %tdev=1st point with value >=half peak value
  diff_indThres=diff(indThres);
  if max(diff_indThres)>1 %trace has double peak: tdec is half way down 1st peak
	tdec=find(diff_indThres>1);
	tdec=indThres(tdec(1));
  else %trace has one peak: tdec=last point with value >=half peak value
	tdec=indThres(length(indThres));
  end

  %linearly interpolate between tdev and preceeding point
  distA=abs(A(tdev)-0.5*Amax);
  distB=abs(A(max(1,tdev-1))-0.5*Amax);
  Rtdev=(distA*R(:,max(1,tdev-1))+distB*R(:,tdev))./(distA+distB);
  %if distA>distB, Rtdev=R(:,tdev-1); end
  
  %linearly interpolate between tdec and following point
  distA=abs(A(tdec)-0.5*Amax);
  distB=abs(A(min(numTimes,tdec+1))-0.5*Amax);
  Rtdec=(distA*R(:,min(numTimes,tdec+1))+distB*R(:,tdec))./(distA+distB);
  %if distA>distB, Rtdec=R(:,tdec+1); end
  
  if ~plotA & ~plotR
	maxsubplot(numNodes+1,colPlus+numTimes+2,(colPlus+numTimes+2)*(nodeNum)+toff+colPlus+2);
	plotRnode=1;
  end
  if plotR
	maxsubplot(numNodes+1,colPlus+2,(colPlus+2)*(nodeNum)+colPlus+2);
	plotRnode=1;
  end
  if plotRnode
	plot(Rtdev,'r','LineWidth',3), hold on, 
	plot(R(:,tmax),'k--','LineWidth',3), plot(Rtdec,'b','LineWidth',3), 
	plot([0,numBins+2],[0,0],'k')
	ax=axis; axis([1,numBins+1,ax(3),ax(4)])
	if nodeNum==1 & ~isempty(bin_labels)
	  set(gca,'XTick',binsToLabel,'XTickLabel',num2str(bin_labels(binsToLabel)'));
	  xlabel('Orientation (degrees)','FontSize',25)
	else
	  set(gca,'XTick',[]);	
	end
 	ylabel('OT(\tau,\theta)','FontSize',25)
	set(gca,'FontSize',22);
  end
  
  %align temporal parameters to peak of amplitude
  [Apeak(node),tpeak]=max(A);
  reqdpeak=2-min(plotTimes);
  peakdiff=tpeak-reqdpeak;
  if peakdiff<0
	disp(['ERROR: poorly written code, make redqpeak the time of the earliest' ...
		  ' peak in the data set'])
  elseif peakdiff>0
	A=A(peakdiff+1:numTimes); A(length(A):numTimes)=A(length(A));
	Rmax=Rmax(peakdiff+1:numTimes); Rmax(length(Rmax):numTimes)=Rmax(length(Rmax));
	Rmin=Rmin(peakdiff+1:numTimes); Rmin(length(Rmin):numTimes)=Rmin(length(Rmin));
	Rorth=Rorth(peakdiff+1:numTimes); Rorth(length(Rorth):numTimes)=Rorth(length(Rorth));
  end
  [poo,tpeak]=max(A); 
  if tpeak~=reqdpeak, disp('ERROR: failed to align curves'), end

  %record average, across nodes, of tuning curve parameters through time 
  meanA=meanA+A;
  meanRmax=meanRmax+Rmax;
  meanRmin=meanRmin+Rmin;
  meanRorth=meanRorth+Rorth;
  
  if node<=maxNodes & ~isempty(w)
	%show corresponding synaptic weigths
	RF=w{node,1}-w{node,2};
	lim=max(1e-9,0.9.*max(max(abs(RF))));
	maxsubplot(numNodes+1,2+numTimes,(numTimes+2)*(nodeNum)+1);
	imagesc(RF,[-lim,lim]);
	axis('equal','tight'), set(gca,'XTick',[],'YTick',[]); drawnow;
  end
end
set(gcf,'PaperPosition',21.5.*[0 0 1.25*max(1,plotA+plotR) 0.96*numNodes]);

%plot mean tuning curve parameters through time 
meanA=meanA./length(plotNodes)
meanRmax=meanRmax./length(plotNodes)
meanRmin=meanRmin./length(plotNodes);
meanRorth=meanRorth./length(plotNodes);
figure
plot(meanA,'r','LineWidth',4), hold on, plot(meanRmax,'m-.','LineWidth',2); 
plot(meanRmin,'b--','LineWidth',4), plot(meanRorth,'g','LineWidth',4);
plot([0,numTimes+2],[0,0],'k')
ax=axis; axis([1,numTimes,ax(3),ax(4)])
[Rmax,tpeak]=max(meanA);
set(gca,'FontSize',22);
toLabel=unique([min(plotTimes),0:abs(min(plotTimes)):max(plotTimes)]);
set(gca,'XTick',unique(1+toLabel-min(plotTimes)),'XTickLabel',int2str(toLabel'))
xlabel('Time (iterations)','FontSize',25)
ylabel('OT_{min}(\tau), OT_{orth}(\tau), OT_{max}(\tau), A(\tau)','FontSize',22)
set(gcf,'PaperPosition',[0 0 12.9 12.8*0.99]);