function [M,ind,n_swaps] = swapping(m)
%SWAPPING  Confluent permutation by swapping adjacent elements.
%          [ISWAP,IND,N_SWAPS] = SWAPPING(M) takes a vector M containing
%          the integers 1:k (some repeated if K < LENGTH(M))
%          and constructs a swapping scheme that produces
%          a confluent permutation, with elements ordered by ascending
%          average position. The confluent permutation is obtained by using
%          the LAPACK routine ZTREX to move m(ISWAP(i,2)) to m(ISWAP(i,1))
%          by swapping adjacent elements, for i = 1:SIZE(M,1).
%          The cell array vector IND defines the resulting block form:
%          IND{i} contains the indices of the i'th block in the permuted form.
%          N_SWAPS is the total number of swaps required.

mmax = max(m); M = []; ind = {}; h = zeros(1,mmax);
g = zeros(1,mmax);

for i = 1:mmax
    p = find(m==i);
    h(i) = length(p);
    g(i) = sum(p)/h(i);
end

[x,y] = sort(g);
mdone = 1;

for i = y
    if any(m(mdone:mdone+h(i)-1) ~= i)
        f = find(m==i); g = mdone:mdone+h(i)-1;
        ff = f(f~=g); gg = g(f~=g);

      % Create vector v = mdone:f(end) with all elements of f deleted.
        v = mdone-1 + find(m(mdone:f(end)) ~= i);

      %  v = zeros(1,f(end)-g(1)+1);
      %  v(f-g(1)+1) = 1; v = g(1)-1 + find(v==0);

        M(end+1:end+length(gg),:) = [gg' ff'];

        m(g(end)+1:f(end)) = m(v);
        m(g) = i*ones(1,h(i));
        ind = cat(2,ind,{mdone:mdone+h(i)-1}); mdone = mdone + h(i);
    else
        ind = cat(2,ind,{mdone:mdone+h(i)-1}); mdone = mdone + h(i);
    end
end

n_swaps = sum(abs(diff(M')));