Home > NoiseTools > nt_pca_big_old.m

# nt_pca_big_old

## PURPOSE

[topcs,pwr,y]=nt_pca_big(x,N,aggregate,normflag) - approximate PCA of big data

## SYNOPSIS

function [topcs,pwr,y]=nt_pca_big(x,N,aggregate,normflag)

## DESCRIPTION

    [topcs,pwr,y]=nt_pca_big(x,N,aggregate,normflag) - approximate PCA of big data

      r: matrix to project to PCs
      pwr: power per PC
      y: PCs

      x: data
      N: desired number of PCs
      aggregate: 'local' [default], 'global', or 'globmean'
      normflag: normalize at each level [default: don't]

     Data are processed recursively in chunks of 2*N columns.  Each chunk is
     submitted to PCA and either N PCs are retained ('local') or 1 PC is
     retained ('global').  If 'globmean' the data are averaged over chunks
     of nchans/N channels, then PCA

     NoiseTools

## CROSS-REFERENCE INFORMATION

This function calls:

- nt_fold y=fold(x,epochsize) - fold 2D to 3D
- nt_greetings nt_greetings - display message the first time the toolbox is used
- nt_imagescc imagescc - plot image with symmetric scaling
- nt_normcol [y,norm]=nt_normcol(x,w) - normalize each column so its weighted msq is 1
- nt_pca0 [topcs,pwr,y]=nt_pca0(x,shifts,nkeep,threshold,w) - time-shift pca
- nt_pca_big [topcs,pwr,y]=nt_pca_big(x,N,normflag) - approximate PCA of big data
- nt_unfold y=nt_fold(x) - unfold 3D to 2D

This function is called by:

## SOURCE CODE

    0001 function [topcs,pwr,y]=nt_pca_big(x,N,aggregate,normflag)
    0002 %[topcs,pwr,y]=nt_pca_big(x,N,aggregate,normflag) - approximate PCA of big data
    0003 %
    0004 %  r: matrix to project to PCs
    0005 %  pwr: power per PC
    0006 %  y: PCs
    0007 %
    0008 %  x: data
    0009 %  N: desired number of PCs
    0010 %  aggregate: 'local' [default], 'global', or 'globmean'
    0011 %  normflag: normalize at each level [default: don't]
    0012 %
    0013 % Data are processed recursively in chunks of 2*N columns.  Each chunk is
    0014 % submitted to PCA and either N PCs are retained ('local') or 1 PC is
    0015 % retained ('global').  If 'globmean' the data are averaged over chunks
    0016 % of nchans/N channels, then PCA
    0017 %
    0018 % NoiseTools
    0019 nt_greetings();
    0020
    0021 if nargin<1; error('!'); end
    0022 if nargin<2||isempty(N);
    0023     N=100;
    0024     disp('defaulting to N=30');
    0025 end
    0026 if nargin<3||isempty(aggregate); aggregate='local'; end
    0027 if nargin<4||isempty(normflag); normflag=0; end
    0028
    0029 [nsamples,nchans,ntrials]=size(x);
    0030 x=nt_unfold(x);
    0031
    0032 if normflag;
    0033     tonorm=1./sqrt(mean(x.^2));
    0034     tonorm(find(isnan(tonorm)))=0;
    0035     x=bsxfun(@times,x,tonorm);
    0036 end
    0037 if ~ (strcmp(aggregate,'local') || strcmp(aggregate,'global') || strcmp(aggregate,'globmean'))
    0038     error('!');
    0039 end
    0040
    0041
    0042
    0043 if nchans<=2*N;
    0044     [topcs,pwr]=nt_pca0(x);
    0045 else
    0046     if strcmp(aggregate,'local')
    0047         topcs=zeros(nchans,N*2);
    0048         a=round(nchans/2);
    0049         [topcs(1:a,1:N)]=nt_pca_big(x(:,1:a),N,aggregate,normflag);
    0050         [topcs(a+1:end,N+1:2*N)]=nt_pca_big(x(:,a+1:end),N,aggregate,normflag);
    0051         [topcs2,pwr,z]=nt_pca_big(x*topcs,N,aggregate,normflag);
    0052         topcs=topcs*topcs2;
    0053     elseif strcmp(aggregate,'global')
    0054         nChunks=ceil(nchans/(2*N));
    0055         chunkSize=ceil(nchans/nChunks);
    0056         nKeep=ceil(N/nChunks);
    0057         topcs=zeros(nchans,nChunks*nKeep);
    0058         for iChunk=1:nChunks
    0059             a=(iChunk-1)*chunkSize;
    0060             idx=(a+1):min(a+chunkSize,nchans);
    0061             tp=nt_pca_big(x(:,idx),N,aggregate,normflag);
    0062             topcs(idx,(iChunk-1)*nKeep+(1:nKeep))=tp(:,1:nKeep);
    0063         end
    0064         %figure(1); clf; nt_imagescc(topcs);
    0065         [topcs2,pwr,z]=nt_pca_big(x*topcs,N,aggregate,normflag);
    0066         topcs=topcs*topcs2;
    0067     else % globmean
    0068         chunkSize=ceil(nchans/(2*N));
    0069         topcs=zeros(nchans,2*N);
    0070         for iChunk=1:2*N
    0071             a=(iChunk-1)*chunkSize;
    0072             idx=(a+1):min(a+chunkSize,nchans);
    0073             topcs(idx,iChunk)=1/numel(idx);
    0074             %disp(idx)
    0075         end
    0076         %figure(1); clf; nt_imagescc(topcs);
    0077         [topcs2,pwr,z]=nt_pca_big(x*topcs,N,aggregate,normflag);
    0078         topcs=topcs*topcs2;
    0079     end
    0080
    0081 end
    0082 if size(topcs,2)>=N;
    0083     topcs=topcs(:,1:N); % truncate to save space
    0084     pwr=pwr(1:N);
    0085 else
    0086     topcs(:,end:N)=0; % should rarely happen
    0087     pwr(end:N)=0;
    0088 end
    0089
    0090 %disp([size(topcs)])
    0091
    0092 if normflag; topcs=bsxfun(@times,topcs,tonorm'); end
    0093
    0094 if nargout>1; y=x*topcs; y=nt_fold(y,nsamples);end
    0095
    0096 if ~nargout % don't return, just plot something
    0097     y=nt_unfold(y);
    0098     semilogy(mean(y.^2), '.-');
    0099     xlabel('pc'); ylabel('power');
    0100     figure(2); clf
    0101     plot(y); xlabel('sample');
    0102     figure(3); clf
    0103     nt_imagescc(nt_normcol(x'*y)');
    0104     ylabel('pc'); xlabel('channel');
    0105     clear topcs y
    0106 end
    0107
    0108
    0109
    0110
    0111
    0112
    0113 %% test code
    0114 if 0
    0115     % basic
    0116     x=randn(10000,3000);
    0117     [topcs,pwr,y]=nt_pca_big(x,10,'local');
    0118     figure(1); clf; plot(pwr);
    0119 end
    0120 if 0
    0121     x=sin(2*pi*3*(1:1000)'/1000)*randn(1,1000);
    0122     x=2*x + randn(size(x));
    0123      [topcs,pwr,y]=nt_pca_big(x);
    0124     figure(1); plot(pwr);
    0125     figure(2); subplot 121; plot(x); subplot 122; plot(x*topcs);
    0126 end
    0127 if 0
    0128     nchans=100000;
    0129     nnoise=900;
    0130     x=repmat(sin(2*pi*3*(1:1000)'/1000), 1, nchans);
    0131     x= 0.06*nt_normcol(x) + nt_normcol(randn(size(x,1),nnoise)*randn(nnoise,nchans));
    0132     [topcs,pwr,y]=nt_pca_big(x,100,'global');
    0133     figure(1); plot(pwr);
    0134     figure(2); subplot 121; plot(x(:,1:100)); subplot 122; plot(y(:,1));
    0135 end
    0136 if 0
    0137     nchans=100000;
    0138     nnoise=900;
    0139     x=repmat(sin(2*pi*3*(1:1000)'/1000), 1, nchans);
    0140     x= 0.01*nt_normcol(x) + nt_normcol(randn(size(x,1),nnoise)*randn(nnoise,nchans));
    0141     [topcs,pwr,y]=nt_pca_big(x,100,'globmean');
    0142     figure(1); plot(pwr);
    0143     figure(2); subplot 121; plot(x(:,1:100)); subplot 122; plot(y(:,1));
    0144 end
    0145