function martonsketchpad
global m n
m = 5;
n = 1;

while 1
  p = rand(m^n,1); p = p/sum(p);
  q = rand(m^n,1); q = q/sum(q);

  d1 = dL1(p,q)
  %d2 = tv(p,q)
  %d3 = dbar(p,q)
  d4 = dbar1(p,q)
  tryfail(abs(d1-d4),0,'bad');
  % good! established that dL1 = tv
  
  %tryfail(d3,d1,'bad');
  
end
return


function d = dbar1(p,q)
% this is a linear program
%m = length(p);
global m n
MN  = repmat([m],[1 n ]);
udim = m^n;
LB = zeros(udim,1);
UB =  ones(udim,1);

rind = 0;
Ain = zeros(rind,udim);
Bin = zeros(rind,1);

Ain = eye(udim);
Bin = min([p q]')';

c = ones(udim,1);

%[u,s] = linprog(-c,Ain,Bin,[],[],LB,UB);  
%s = -s;
u = Bin;
s = sum(u);

d = 1-s;

x = p-q;

ip = find(x>0);
iq = find(x<0);

tryfail(abs(1/2*sum(x(ip))-sum(u(ip))),0,'bad1');
tryfail(abs(1/2*sum(x(iq))-sum(u(iq))),0,'bad2');


return



function d = dbar(p,q)
% this is a linear program
%m = length(p);
global m n
MN  = repmat([m],[1 n ]);
udim = m^(2*n);
LB = zeros(udim,1);
UB =  ones(udim,1);

rind = 0;
Aeq = zeros(rind,udim);
Beq = zeros(rind,1);

rind=rind+1;
Aeq(rind,:) = ones(1,udim);
Beq(rind,1) = 1;
% it's a joint probability

for xind=1:m^n
  rind = rind+1;
  %x = ind2coord(xind,MN);
  for yind=1:m^n
    %y = ind2coord(yind,MN);
    ind = coord2ind([xind yind],[m^n m^n]);
    Aeq(rind,ind) = 1;
  end
  Beq(rind,1) = p(xind);
end

for yind=1:m^n
  rind = rind+1;
  %y = ind2coord(yind,MN);
  for xind=1:m^n
    %x = ind2coord(xind,MN);
    ind = coord2ind([xind yind],[m^n m^n]);
    Aeq(rind,ind) = 1;
  end
  Beq(rind,1) = q(yind);
end

c = zeros(udim,1);
for xind=1:m^n
  x = ind2coord(xind,MN);
  for yind=1:m^n
    y = ind2coord(yind,MN);
    ind = coord2ind([xind yind],[m^n m^n]);
    c(ind) = rho(x,y);
  end
end

[u,d] = linprog(c,[],[],Aeq,Beq,LB,UB);  
u = reshape(u,[m^n m^n]);
tryfail(abs(1-sum(diag(u))-d),0,'hmm');
return


function r=rho(x,y)
%r = (x~=y);
%r = ~all(x==y);
%return
l = length(x);
r = sum(x~=y)/l;
return


function d = dL1(p,q)
d = 1/2*sum(abs(p-q));
return

function d = tv(p,q)
% total variation norm
%p=p(:)';q=q(:)';
n = length(p);
d = 0;
for i=1:2^n
  setstr = dec2bin(i-1,n);
  B = find(setstr=='1');
  d = max(d,(sum(p(B)-q(B))));
end
return




function tryfail(lhs,rhs,failstr)
good = (lhs <= rhs+1e-7);
if ~good
  failstr
  keyboard
end
return