CellsortICA

 [ica_sig, ica_filters, ica_A, numiter] = CellsortICA(mixedsig, mixedfilters, PCuse, mu, nIC, ica_A_guess, termtol, maxrounds)

CELLSORT
 Perform ICA with a standard set of parameters, including skewness as the
 objective function

 Inputs:
   mixedsig - N x T matrix of N temporal signal mixtures sampled at T
   points.
   mixedfilters - N x X x Y array of N spatial signal mixtures sampled at
   X x Y spatial points.
   CovEvals - eigenvalues of the covariance matrix
   PCuse - vector of indices of the components to be included. If empty,
   use all the components
   mu - parameter (between 0 and 1) specifying weight of temporal
   information in spatio-temporal ICA
   nIC - number of ICs to derive
   termtol - termination tolerance; fractional change in output at which
   to end iteration of the fixed point algorithm.
   maxrounds - maximum number of rounds of iterations

 Outputs:
     ica_sig - nIC x T matrix of ICA temporal signals
     ica_filters - nIC x X x Y array of ICA spatial filters
     ica_A - nIC x N orthogonal unmixing matrix to convert the input to output signals
     numiter - number of rounds of iteration before termination

 Routine is based on the fastICA package (Hugo Gävert, Jarmo Hurri, Jaakko Särelä, and Aapo
 Hyvärinen, http://www.cis.hut.fi/projects/ica/fastica)

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

0001 function [ica_sig, ica_filters, ica_A, numiter] = CellsortICA(mixedsig, ...
0002     mixedfilters, CovEvals, PCuse, mu, nIC, ica_A_guess, termtol, maxrounds)
0003 % [ica_sig, ica_filters, ica_A, numiter] = CellsortICA(mixedsig, mixedfilters, PCuse, mu, nIC, ica_A_guess, termtol, maxrounds)
0004 %
0005 %CELLSORT
0006 % Perform ICA with a standard set of parameters, including skewness as the
0007 % objective function
0008 %
0009 % Inputs:
0010 %   mixedsig - N x T matrix of N temporal signal mixtures sampled at T
0011 %   points.
0012 %   mixedfilters - N x X x Y array of N spatial signal mixtures sampled at
0013 %   X x Y spatial points.
0014 %   CovEvals - eigenvalues of the covariance matrix
0015 %   PCuse - vector of indices of the components to be included. If empty,
0016 %   use all the components
0017 %   mu - parameter (between 0 and 1) specifying weight of temporal
0018 %   information in spatio-temporal ICA
0019 %   nIC - number of ICs to derive
0020 %   termtol - termination tolerance; fractional change in output at which
0021 %   to end iteration of the fixed point algorithm.
0022 %   maxrounds - maximum number of rounds of iterations
0023 %
0024 % Outputs:
0025 %     ica_sig - nIC x T matrix of ICA temporal signals
0026 %     ica_filters - nIC x X x Y array of ICA spatial filters
0027 %     ica_A - nIC x N orthogonal unmixing matrix to convert the input to output signals
0028 %     numiter - number of rounds of iteration before termination
0029 %
0030 % Routine is based on the fastICA package (Hugo Gävert, Jarmo Hurri, Jaakko Särelä, and Aapo
0031 % Hyvärinen, http://www.cis.hut.fi/projects/ica/fastica)
0032 %
0033 % Eran Mukamel, Axel Nimmerjahn and Mark Schnitzer, 2009
0034 % Email: eran@post.harvard.edu, mschnitz@stanford.edu
0035 
0036 fprintf('-------------- CellsortICA %s -------------- \n', date)
0037 
0038 if (nargin<4) || isempty(PCuse)
0039     PCuse = [1:size(mixedsig,1)];
0040 end
0041 if (nargin<6) || isempty(nIC)
0042     nIC = length(PCuse);
0043 end
0044 if (nargin<7) || isempty(ica_A_guess)
0045     ica_A_guess = randn(length(PCuse), nIC);
0046 end
0047 if (nargin<8) || isempty(termtol)
0048     termtol = 1e-6;
0049 end
0050 if (nargin<9) || isempty(maxrounds)
0051     maxrounds = 100;
0052 end
0053 if isempty(mu)||(mu>1)||(mu<0)
0054     error('Spatio-temporal parameter, mu, must be between 0 and 1.')
0055 end
0056 
0057 % Check that ica_A_guess is the right size
0058 if size(ica_A_guess,1)~= length(PCuse) || size(ica_A_guess,2)~=nIC
0059     error('Initial guess for ica_A is the wrong size.')
0060 end
0061 if nIC>length(PCuse)
0062     error('Cannot estimate more ICs than the number of PCs.')
0063 end
0064     
0065 [pixw,pixh] = size(mixedfilters(:,:,1));
0066 npix = pixw*pixh;
0067 
0068 % Select PCs
0069 if mu > 0 || ~isempty(mixedsig)
0070     mixedsig = mixedsig(PCuse,:);
0071 end
0072 if mu < 1 || ~isempty(mixedfilters)
0073     mixedfilters = reshape(mixedfilters(:,:,PCuse),npix,length(PCuse));
0074 end
0075 CovEvals = CovEvals(PCuse);
0076 
0077 % Center the data by removing the mean of each PC
0078 mixedmean = mean(mixedsig,2);
0079 mixedsig = mixedsig - mixedmean * ones(1, size(mixedsig,2));
0080 
0081 % Create concatenated data for spatio-temporal ICA
0082 nx = size(mixedfilters,1);
0083 nt = size(mixedsig,2);
0084 if mu == 1
0085     % Pure temporal ICA
0086     sig_use = mixedsig;
0087 elseif mu == 0
0088     % Pure spatial ICA
0089     sig_use = mixedfilters';
0090 else
0091     % Spatial-temporal ICA
0092     sig_use = [(1-mu)*mixedfilters', mu*mixedsig];
0093     sig_use = sig_use / sqrt(1-2*mu+2*mu^2); % This normalization ensures that, if both mixedfilters and mixedsig have unit covariance, then so will sig_use
0094 end
0095 
0096 % Perform ICA
0097 [ica_A, numiter] = fpica_standardica(sig_use, nIC, ica_A_guess, termtol, maxrounds);
0098 
0099 % Sort ICs according to skewness of the temporal component
0100 ica_W = ica_A';
0101 
0102 ica_sig = ica_W * mixedsig;
0103 ica_filters = reshape((mixedfilters*diag(CovEvals.^(-1/2))*ica_A)', nIC, nx);  % This is the matrix of the generators of the ICs
0104 ica_filters = ica_filters / npix^2;
0105 
0106 icskew = skewness(ica_sig');
0107 [icskew, ICord] = sort(icskew, 'descend');
0108 ica_A = ica_A(:,ICord);
0109 ica_sig = ica_sig(ICord,:);
0110 ica_filters = ica_filters(ICord,:);
0111 ica_filters = reshape(ica_filters, nIC, pixw, pixh);
0112 
0113 % Note that with these definitions of ica_filters and ica_sig, we can decompose
0114 % the sphered and original movie data matrices as:
0115 %     mov_sphere ~ mixedfilters * mixedsig = ica_filters * ica_sig = (mixedfilters*ica_A') * (ica_A*mixedsig),
0116 %     mov ~ mixedfilters * pca_D * mixedsig.
0117 % This gives:
0118 %     ica_filters = mixedfilters * ica_A' = mov * mixedsig' * inv(diag(pca_D.^(1/2)) * ica_A'
0119 %     ica_sig = ica_A * mixedsig = ica_A * inv(diag(pca_D.^(1/2))) * mixedfilters' * mov
0120 
0121     function [B, iternum] = fpica_standardica(X, nIC, ica_A_guess, termtol, maxrounds)
0122         
0123         numSamples = size(X,2);
0124         
0125         B = ica_A_guess;
0126         BOld = zeros(size(B));
0127         
0128         iternum = 0;
0129         minAbsCos = 0;
0130         
0131         errvec = zeros(maxrounds,1);
0132         while (iternum < maxrounds) && ((1 - minAbsCos)>termtol)
0133             iternum = iternum + 1;
0134             
0135             % Symmetric orthogonalization.
0136             B = (X * ((X' * B) .^ 2)) / numSamples;
0137             B = B * real(inv(B' * B)^(1/2));
0138             
0139             % Test for termination condition.
0140             minAbsCos = min(abs(diag(B' * BOld)));
0141             
0142             BOld = B;
0143             errvec(iternum) = (1 - minAbsCos);
0144         end
0145         
0146         if iternum<maxrounds
0147             fprintf('Convergence in %d rounds.\n', iternum)
0148         else
0149             fprintf('Failed to converge; terminating after %d rounds, current change in estimate %3.3g.\n', ...
0150                 iternum, 1-minAbsCos)
0151         end
0152     end
0153 end
0154