function [X,D] = ordercols(M,S)
% X = ordercols(M,S)
%      X returns an indexing into the columns of M that 
%      tries to order them to put similar columns next to each 
%      other.  Greedy, heuristic, home-baked algorithm.
%      Optional S forces start point
% 2004-05-09 dpwe@ee.columbia.edu

[nr,nc] = size(M);

X = zeros(1,nc);

% Form similarity between all columns as normalized xcorr (rho)
DD = (M'*M)./sqrt(repmat(sum(M.^2),nc,1).*repmat(sum(M.^2)',1,nc));
%DD = M'*M;

DD = zeros(nc, nc);
for i = 1:nc;
  DD(i,:) = 1./max(.001,sum((M - repmat(M(:,i),1,nc)).^2));
end



mind = min(min(DD))-1;
DD = DD .* (1-eye(nc));
DD = DD + eye(nc)*mind;
D = DD;

%% Find the most similar pair (off leading diag)
%bc = find(max(DD) == max(max(DD)));
%br = find(DD(:,bc) == max(max(DD)));
%bc = bc(1);
%br = br(1);
%% Start with the one with the less-good next option
%if max(DD(:,bc)*([1:nc]~=br)) > max(DD(:,br)*([1:nc]~=bc))
%  X(1) = br;
%else
%  X(1) = bc;
%end

if nargin < 2

  % Instead, start from all possible initial columns...
XX = zeros(nc,nc);
XX(:,1) = [1:nc]';

% .. and record total similarity for each row (to choose best at end)
ts = zeros(1,nc);

for rr = 1:nc

  X = XX(rr,:);
  totscore = 0;
  
  % Remove first column choice
  DD = D;
  DD(:,X(1)) = mind;

  for i = 2:nc
    % Now just choose greedily
    bc = find( DD(X(i-1),:) == max(DD(X(i-1),:)) );
    X(i) = bc(1);
    %  .. and remove cols as we go along
    %  DD(X(i),:) = mind;
    totscore = totscore + DD(X(i-1),X(i));
    DD(:,X(i)) = mind;
  end

  XX(rr,:) = X;
  ts(rr) = totscore;

end

  br = find(ts == max(ts));
  br = br(1);
  X = XX(br, :);
else
  X = zeros(1,nc);
  X(1) = S;
  totscore = 0;
  DD = D;
  DD(:,X(1)) = mind;
  for i = 2:nc
    % Now just choose greedily
    bc = find( DD(X(i-1),:) == max(DD(X(i-1),:)) );
    X(i) = bc(1);
    %  .. and remove cols as we go along
    %  DD(X(i),:) = mind;
    totscore = totscore + DD(X(i-1),X(i));
    DD(:,X(i)) = mind;
  end
end