% [Perm, U, V, n] = partition (Sx, nx, Sy, ny, tol)
% Finds U and V such that  Perm = U' * Sx' * Sy * V 
% = [I 0 0 0;
%    0 0 I 0;
%    0 I 0 0;
%    0 0 0 I];
% where the number of elements in each block of rows have
% dimensions given by n
%
% n(3)+n(4) = ns-nx(1) and n(2)+n(4) = ns-ny(1)
% ns = size (Sx,1) = size (Sy,1)

function [Perm, U, V, nn] = partition (Sx, nx, Sy, ny, tol)

if nargin < 2,
  disp ('usage:  [Perm, U, V, n] = partition (Sx, nx, Sy, ny, tol)');
  return;
end

if nargin < 5, tol = 1e-12; end

[ns, nss] = size (Sx);
if ns ~= nss, error ('Sx must be square'); end
if [ns, nss] ~= size (Sy), error ('Sx and Sy do not match'); end
if sum (nx) ~= ns, error ('Sx and nx do not match'); end
if sum (ny) ~= ns, error ('Sy and ny do not match'); end

RR0 = Sx' * Sy;
RR = RR0;
U = eye(ns);
V = eye(ns);
nn = zeros(length(nx),length(ny));

for k = 1:length(nx),
  for i = 1:length(ny),
    dx = sum(nn(k,1:i-1));
    dy = sum(nn(1:k-1,i));
    indx = 1+sum(nx(1:k-1))+dx:ns;
    indy = 1+sum(ny(1:i-1))+dy:ns;
    [RR(indx,indy), U_, V_, ni] = ...
       partition1 (RR(indx,indy), nx(k)-dx, ny(i)-dy, tol);
    U(:,indx) = U(:,indx) * U_;
    V(:,indy) = V(:,indy) * V_;
    nn(k,i) = ni;
  end
end

Perm = double (RR > 0.5);

if any (sum(Perm) ~= 1) || any (sum(Perm') ~= 1),
  error ('internal: Perm has bad structure');
end
if (norm (U' * RR0 * V - Perm) > tol*sqrt(ns)),
  warning ('High tolerances in partitioning');
  norm (U' * RR0 * V - Perm)
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [RR, U, V, n] = partition1 (RR0, nx, ny, tol)
[u, s, v] = svd (RR0(1:nx,1:ny));
n = min (size(s));
n = sum (diag (s(1:n,1:n)) > tol);
Si = diag (1./sqrt(diag(s(1:n,1:n))));

[nsx, nsy] = size (RR0);
U1 = blkdiag (u*blkdiag(Si,eye(nx-n)), eye(nsx-nx));
V1 = blkdiag (v*blkdiag(Si,eye(ny-n)), eye(nsy-ny));
RR = U1' * RR0 * V1;

U2 = [eye(n) -RR(n+1:nsx,1:n)'; zeros(nsx-n,n) eye(nsx-n)];
V2 = [eye(n) -RR(1:n,n+1:nsy); zeros(nsy-n,n) eye(nsy-n)];
RR = U2' * RR * V2;
U = U1 * U2;
V = V1 * V2;
