function [result] = SubsetDev(head, tail, count)
%--------------------------------------------------------------------------
% for every subset of components having the given number of elements,
% calculate the square deviation of the scalp potential resulting from that
% subset from the grand average scalp potential, summed over the whole
% scalp for all samples
% arguments:
% head:  the components already chosen; initially []
% tail:  the components available to choose from; intially [1:comps]
% count: the number of elements in each subset
%--- BCT 05.02.15 ---------------------------------------------------------

  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

  result = [];
  if count > 1
    for i = tail
      tail = tail(tail ~= i);
      result = [result SubsetDev([head i], tail, count - 1)];
    end
  else
    for i = tail
      tuple = [head i];
      result = [result ...
        [tuple sum(sum((A(:,tuple) * sm(tuple,:) - xm).^2))]' ...
      ];
    end
  end
end
