% COMPSORT - reorder ICA components, first largest to smallest by the size
% of their maximum variance in the single-component projections,
% then (if specified) the nlargest component projections are
% reordered by the (within-epoch) time point at which they reach
% their max variance.
%
% Usage:
% >> [windex,maxvar,maxframe,maxepoch,maxmap] ...
% = compsort(data,weights,sphere,datamean, ...
% frames,nlargest,chanlist);
% Inputs:
% data = (chans,frames*epochs) the input data set decomposed by RUNICA
% weights = ica weight matrix returned by RUNICA
% sphere = sphering matrix returned by RUNICA
%
% Optional:
% datamean = means removed from each row*epoch in RUNICA
% (Note: 0 -> input data means are distributed among components
% : 1 -> input data means are removed from the components (default))
% frames = frames per epoch in data (0 -> all)
% nlargest = number of largest ICA components to order by latency to max var
% (other returned in reverse order of variance) (0 -> none)
% chanlist = list of channel numbers to sort on (0 -> all)
%
% Outputs:
% windex = permuted order of rows in output weights (1-chans)
% maxvar = maximum variance of projected components (in perm. order)
% maxframe = frame of maximum variance (1-frames) (in perm. order)
% maxepoch = epoch number of frame of maxvar (1-nepochs) (in perm. order)
% maxmap = projected scalp map at max (in perm. order)
%
% Authors: Scott Makeig, SCCN/INC/UCSD, La Jolla, 1996
%
% See also: RUNICA
% Copyright (C) 1996 Scott Makeig, SCCN/INC/UCSD, scott@sccn.ucsd.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.
% 11-30-96 as compsort.m
% 12-19-96 moved def of epochs below default frames def -sm
% 02-18-97 added chanlist -sm
% 03-11-97 added default chanlist=0 -> all chans, fixed datamean and VAR -sm
% 03-17-97 made nlargest default -> nlargest=0 -sm
% 04-03-97 fixed problems and removed permweights, permcomps outputs -sm
% 04-04-97 shortened name to COMPSORT -sm
% 06-05-97 corrected variance computation -sm
% 06-07-97 changed order of args to conform to runica -sm
% 06-10-97 fixed recent bug in maxvar order -sm
% 07-23-97 made datamean==0 distribute means among components -sm
% 08-04-97 added datamean=1 option -sm
% 01-25-02 reformated help & license, added links -ad
function [windex,maxvar,maxframe,maxepoch,maxmap] = compsort(data,weights,sphere,datamean,frames,nlargest,chanlist)
if nargin<3,
fprintf('compsort(): needs at least three arguments.\n\n');
help compsort
return
end
[chans,framestot] = size(data);
if framestot==0,
fprintf('Gcompsort(): cannot process an empty data array.\n');
return
end
[srows,scols] = size(sphere);
[wrows,wcols] = size(weights);
if nargin<7,
chanlist = [1:chans];
end
if chanlist==0,
chanlist = [1:chans];
end
if length(chanlist)~=chans && wrows<chans, % if non-square weight matrix
fprintf('compsort(): chanlist not allowed with non-square weights.\n');
return
end
if size(chanlist,1)>1
chanlist = chanlist'; % make a row vector
end
if size(chanlist,1)>1
fprintf('compsort(): chanlist must be a vector.\n');
return
end
if nargin<6,
nlargest=0;
end
if nargin<5,
frames = 0;
end
if nargin<4,
datamean = 1;
end
if frames ==0,
frames = framestot;
end
epochs = framestot/frames;
% activations = (wrows,wcols)x(srows,scols)x(chans,framestot)
if chans ~= scols || srows ~= wcols,
fprintf('compsort(): input data dimensions do not match.\n');
fprintf(...
' i.e., Either chans %d ~= sphere cols %d or sphere rows %d ~= weights cols %d\n',...
chans,scols,srows,wcols);
return
end
if wrows ~= chans && nlargest ~= 0 && nlargest ~= wrows,
fprintf(...
'compsort(): cannot project components back to order by size - nchans ~= ncomponents.\n');
return
end
if floor(framestot/frames)*frames ~= framestot
fprintf(...
'compsort(): input data frames does not divide data length.\n');
return
end
if nlargest > wrows,
fprintf(...
'compsort(): there are only %d rows in the weight matrix.\n',wrows);
return
end
if epochs ~= floor(epochs),
fprintf(...
'compsort(): input frames does not subdivide data length.\n');
return
end
%
%%%%%%%%%%%%%%%%%%%% Reorder weight matrix %%%%%%%%%%%%%%%%%%%%%
%
if datamean == 1,
data = data - mean(data')'*ones(1,framestot); % remove channel means
elseif datamean~=0, % remove given means
if size(datamean,2) ~= epochs || size(datamean,1) ~= chans,
fprintf('compsort(): datamean must be 0, 1, or (chans,epochs)\n');
return
end
for e=1:epochs
data(:,(e-1)*frames+1:e*frames) = data(:,(e-1)*frames+1:e*frames)...
- datamean(:,e)*ones(1,frames);
end
end % compute mean data matrix inherited from runica()
comps = weights*sphere*data; % Note: distributes means if datamean==0
maxvar = zeros(wrows,1); % size of the projections
maxframe = zeros(wrows,1); % frame of the abs(max) projection
maxepoch = zeros(wrows,1); % epoch of the abs(max) projection
maxmap = zeros(wrows,chans); % leave 0s unless weights is square
if chans==wrows,
icainv = inv(weights*sphere);
fprintf('Computing projected variance for all %d components:\n',wrows);
for s=1:wrows
fprintf('%d ',s); % construct single-component data matrix
% project to scalp
compproj = icainv(:,s)*comps(s,:);
compv = zeros(frames*epochs,1);
compv = (sum(compproj(chanlist,:).*compproj(chanlist,:)))' ...
/(length(chanlist)-1);
[m,mi] = max(compv); % find max variance
maxvar(s) = m;
maxframe(s)=rem(mi,frames); % count from beginning of each epoch!
maxepoch(s)=floor(mi/frames)+1;% record epoch number of max
if maxframe(s)==0, % if max var is in last frame . . .
maxframe(s) = frames;
maxepoch(s) = maxepoch(s)-1;
end
maxmap(:,s) = compproj(:,mi); % record scalp projection at max(var(proj))
end
else % weight matrix is non-square, sort components by latency only
fprintf('compsort() - non-square weights - finding max latencies.\n');
for s=1:wrows
compv = comps(s).*comps(s,:)'; % get variance analogues at each time point
[m,mi] = max(abs(compv)'); % find abs(max)'s
maxvar(s) = m;
maxframe(s)=rem(mi,frames); % count from beginning of each epoch!
maxepoch(s)=floor(mi/frames)+1;% record epoch number of max
if maxframe(s)==0, % if max var is in last frame . . .
maxframe(s) = frames;
maxepoch(s) = maxepoch(s)-1;
end
end
end
fprintf('\n');
[maxvar windex] = sort(maxvar');% sort components by relative size
windex = windex(wrows:-1:1)'; % reverse order
maxvar = maxvar(wrows:-1:1)'; % make each returned vector a column vector
% Note: maxvar reordered by SORT above
maxframe = maxframe(windex); % reorder maxframes
maxepoch = maxepoch(windex); % reorder maxepoch
maxmap = maxmap(:,windex); % reorder maxmap columns
if nlargest>0,
fprintf('Ordering largest %d components by frame of abs(max): ',nlargest);
[m mfi] = sort(maxframe(1:nlargest)); % sort by frame order
windex(1:nlargest) = windex(mfi);
maxframe(1:nlargest) = m;
maxvar(1:nlargest) = maxvar(mfi);
maxepoch(1:nlargest) = maxepoch(mfi);
maxmap(:,1:nlargest) = maxmap(:,mfi); % reorder largest maxmap columns
fprintf('\n');
end
Datasets
Models
