function psixy(m0)
global m mrat
%m0 = 2;
m = m0;
mrat = 0;

xdim = 2*m;
R = 1;
%r = 1e-10;
r = .01;
LB = zeros(xdim,1)+r;
UB = zeros(xdim,1)+R;

%A = zeros(0,xdim);
A = assembleblock({ones(1,m),ones(1,m)});
B = ones(2,1);

%LB = LB+.1;
%UB = UB-.1;

while 1
%for i=1
  %psi = rand(m);
  %psi = rand(2,m); psi = psi/min(psi(:));
  a = rand(m,1); a = a/sum(a);
  b = rand(m,1); b = b/sum(b);
  psi = [a,b];  
  
  
  if 0
    psi1 = psi;
    psi = rand(2,m);  psi = psi/min(psi(:));
    psi2 = psi;
    lam = rand;
    good = tha(lam*psi1 + (1-lam)*psi2) <= lam*tha(psi1) + (1-lam)*tha(psi2)
    if ~good
      'not conv'
      keyboard
    end
    a = rand(m,1); a = a/sum(a);
    b = rand(m,1); b = b/sum(b);
    R = max([a;b]); r = min([a;b]); 
    th = .5*sum(abs(a-b));
    %B = (R-r)/(R+r);
    B = (R-r);
    good = (th<=B);
    rat = th/B;
    if ~good
      keyboard
    end
    mrat = max(mrat,rat);
    fprintf('rat = %1.5f     mrat = %1.5f \n',rat,mrat);
  end
  [x,fval] = fmincon(@tha,psi(:),[],[],A,B,LB,UB);
  ab = reshape(x,[m 2]);
  a = ab(:,1)
  b = ab(:,2)
  keyboard
  %x = psi(:);
  %[x,fval] = fmincon(@psiopt,x,[] ,[]  ,[] ,[]  ,0*x+1,0*x+1e3);  
end

return


function y = psiopt(x)
global m
%psi = reshape(x,[m m]);
psi = reshape(x,[2 m]);
R = max(psi(:))/min(psi(:));

%x = 1;
y = 1;
y1 = 2;

th = 0;
for x=1:m
  p1 = psi(y,x) / sum(psi(y,:));
  p2 = psi(y1,x) / sum(psi(y1,:));
  th = th + abs(p1-p2)/2;
end
%th = abs(p1-p2);
%B = 1-exp(.5*(1-R));
B = (R-1)/(R+1);
%r = th / B;
r = th/B;

y = -r;

global mrat mpsi
if r > mrat
  mrat = r;
  mpsi = psi;
end
fprintf('psixy   MRAT = %1.9f\n',mrat);

return


function y = tha(ab)
global m
ab = reshape(ab,[m 2]);

a = ab(:,1);
b = ab(:,2);
th = .5*sum(abs(a-b));

R = max(ab(:));
r = min(ab(:));

B=(R-r)/(R+r);

%y = -th/B;
y = -th;
return