CellsortPlotPCspectrum

CellsortPlotPCspectrum(fn, CovEvals, pcuse)

 Plot the principal component (PC) spectrum and compare with the
 corresponding random-matrix noise floor

 Inputs:
   fn - movie file name. Must be in TIFF format.
   CovEvals - eigenvalues of the covariance matrix
   pcuse - [optional] - indices of PCs included in dimensionally reduced
   data set

 Eran Mukamel, Axel Nimmerjahn and Mark Schnitzer, 2009
 Email: eran@post.harvard.edu, mschnitz@stanford.edu

0001 function CellsortPlotPCspectrum(fn, CovEvals, pcuse)
0002 % function CellsortPlotPCspectrum(fn, CovEvals, pcuse)
0003 %
0004 % Plot the principal component (PC) spectrum and compare with the
0005 % corresponding random-matrix noise floor
0006 %
0007 % Inputs:
0008 %   fn - movie file name. Must be in TIFF format.
0009 %   CovEvals - eigenvalues of the covariance matrix
0010 %   pcuse - [optional] - indices of PCs included in dimensionally reduced
0011 %   data set
0012 %
0013 % Eran Mukamel, Axel Nimmerjahn and Mark Schnitzer, 2009
0014 % Email: eran@post.harvard.edu, mschnitz@stanford.edu
0015 
0016 if nargin<3
0017     pcuse = [];
0018 end
0019 
0020 [pixw,pixh] = size(imread(fn,1));
0021 npix = pixw*pixh;
0022 nt = tiff_frames(fn);
0023 
0024 % Random matrix prediction (Sengupta & Mitra)
0025 p1 = npix; % Number of pixels
0026 q1 = nt; % Number of time frames
0027 q = max(p1,q1);
0028 p = min(p1,q1);
0029 sigma = 1;
0030 lmax = sigma*sqrt(p+q + 2*sqrt(p*q));
0031 lmin = sigma*sqrt(p+q - 2*sqrt(p*q));
0032 lambda = [lmin: (lmax-lmin)/100.0123423421: lmax];
0033 rho = (1./(pi*lambda*(sigma^2))).*sqrt((lmax^2-lambda.^2).*(lambda.^2-lmin^2));
0034 rho(isnan(rho)) = 0;
0035 rhocdf = cumsum(rho)/sum(rho);
0036 noiseigs = interp1(rhocdf, lambda, [p:-1:1]'/p, 'linear', 'extrap').^2 ;
0037 
0038 % Normalize the PC spectrum
0039 normrank = min(nt-1,length(CovEvals));
0040 pca_norm = CovEvals*noiseigs(normrank) / (CovEvals(normrank)*noiseigs(1));
0041 
0042 clf
0043 plot(pca_norm, 'o-', 'Color', [1,1,1]*0.3, 'MarkerFaceColor', [1,1,1]*0.3, 'LineWidth',2)
0044 hold on
0045 plot(noiseigs / noiseigs(1), 'b-', 'LineWidth',2)
0046 plot(2*noiseigs / noiseigs(1), 'b--', 'LineWidth',2)
0047 if ~isempty(pcuse)
0048     plot(pcuse, pca_norm(pcuse), 'rs', 'LineWidth',2)
0049 end
0050 hold off
0051 formataxes
0052 set(gca,'XScale','log','YScale','log', 'Color','none')
0053 xlabel('PC rank')
0054 ylabel('Normalized variance')
0055 axis tight
0056 if isempty(pcuse)
0057     legend('Data variance','Noise floor','2 x Noise floor')
0058 else
0059     legend('Data variance','Noise floor','2 x Noise floor','Retained PCs')
0060 end
0061 
0062 fntitle = fn;
0063 fntitle(fn=='_') = ' ';
0064 title(fntitle)
0065 
0066 function formataxes
0067 
0068 set(gca,'FontSize',12,'FontWeight','bold','FontName','Helvetica','LineWidth',2,'TickLength',[1,1]*.02,'tickdir','out')
0069 set(gcf,'Color','w','PaperPositionMode','auto')
0070 
0071 function j = tiff_frames(fn)
0072 %
0073 % n = tiff_frames(filename)
0074 %
0075 % Returns the number of slices in a TIFF stack.
0076 %
0077 %
0078 
0079 status = 1; j=0;
0080 jstep = 10^3;
0081 while status
0082     try
0083         j=j+jstep;
0084         imread(fn,j);
0085     catch
0086         if jstep>1
0087             j=j-jstep;
0088             jstep = jstep/10;
0089         else
0090             j=j-1;
0091             status = 0;
0092         end
0093     end
0094 end