function d = wasser1(P,Q,m,n)
mn = m^n;
Aeqc = colsum(mn,mn);
Aeqr = rowsum(mn,mn);
%Aeq = assembleblock({Aeqc,Aeqr});
Aeq = [Aeqc;Aeqr;ones(1,mn^2)];
Beq = [P;Q;1];

D = zeros(mn);
for xi=1:mn
  x = i2c(xi,m,n);
  for yi=xi+1:mn
    y = i2c(yi,m,n);
    D(xi,yi) = rho(x,y);
  end
end
D = D + D';
%D = D/n;
D = D(:);
LB = D*0;
UB = LB + 1;

[U,d] = linprog(D,[],[],Aeq,Beq,LB,UB);
U = reshape(U,[mn mn])
return


function Aeq = colsum(m,n)
rind = 0;
Aeq = zeros(0,m*n);
for c=1:n
  rind = rind+1;
  for r=1:m
    ind = coord2ind([r c],[m n]);
    Aeq(rind,ind) = 1;
  end
end
return


function Aeq = rowsum(m,n)
rind = 0;
Aeq = zeros(0,m*n);
for r=1:n
  rind = rind+1;
  for c=1:m
    ind = coord2ind([r c],[m n]);
    Aeq(rind,ind) = 1;
  end
end
return


function c=i2c(i,m,n)
dims = repmat(m,[1 n]);
c = ind2coord(i,dims);
return