% set up the problem

global A;  % [chans x comps] xm = A * sm; the mixing matrix
global xm; % [chans x samps] grand average scalp potentials
global sm; % [comps x samps] independent components

comps = 3; chans = 5; samps = 7; epochs = 9;
% s = rand(comps, samps, epochs);
A = [
  .4 .3 .3
  .8 .1 .1
  .1 .8 .1
  .1 .1 .8
  .3 .3 .4
];
for i = 1:epochs; x(:,:,i) = A * s(:,:,i); end;
xm = mean(x, 3); % grand average scalp potentials
sm = mean(s, 3); % grand average sources

% solve it

result = SubsetDev([], [1:comps], 2);
[val ndx] = sort(result(3,:));
result = result(:,ndx);

% [chans samps] = size(xm); comps = size(sm, 1);
% result = [];
% for k = 1:comps
%   for j = k+1:comps
%     if j ~= k
%       xs = A(:,j)*sm(j,:) + A(:,k)*sm(k,:);
%       c = [];
%       for i = 1:chans
% %         if sum(xs(i,:)) ~= 0
% %           c(i) = xcov(xs(i,:), xm(i,:), 0, 'coeff');
% %         else c(i) = 0; end;
%          c(i) = sum((xs(i,:) - xm(i,:)).^2);
%       end;
%       result = [result [j k sum(c)]'];
%     end;
%   end;
% end;
% [val ndx] = sort(result(3,:));
% fwd = result(:,ndx);
% rev = result(:,fliplr(ndx));
