function samson
global m n mrat
mrat = 0;
m=2; n=5;
xdim = m^n;
Aeq = ones(1,xdim);
Beq = ones(1,1);
LB = zeros(xdim,1);
LB = LB+1e-9;
UB = ones(xdim,1);

while 1
  P = pnormdim(rand(m^n,1),1);
  x = P;
  [x,fv] = fmincon(@myopt,x,[],[],Aeq,Beq,LB,UB);
  fprintf('mrat = %1.9f \n',mrat);
end


return

function y = myopt(x)
global m n mrat
P = x;
[hh0,Ht0,H0] = gethhn(P,m,n);
[hh1,Ht1,H1] = gethhn1(P,m,n);

r = 0;
for i=1:n-1
  for j=i+1:n
    r = max(r,hh0(i,j)/hh1(i,j));
  end
end
if r > 1.2
  keyboard
end
mrat = max(r,mrat);
fprintf('MRAT = %1.9f\n',mrat);
y = -r;
return

%----
mrat = 0;
while 1
  P = pnormdim(rand(m^n,1),1);
  
  [hh0,Ht0,H0] = gethhn(P,m,n);
  [hh1,Ht1,H1] = gethhn1(P,m,n);
  
  %good = all(all(hh0<=hh1))
  %if ~good
  %  keyboard
  %end
  r = 0;
  for i=1:n-1
    for j=i+1:n
      r = max(r,hh0(i,j)/hh1(i,j));
    end
  end
  mrat = max(r,mrat);
  fprintf('MRAT = %1.9f\n',mrat);
end
return

function [hh,Ht,H] = gethhn1(P,m,n)
hh = zeros(n);
for i=1:n
  for j=i+1:n
    h = etaij(P,m,n,i,j);
    hh(i,j) = h;
  end
end
Ht = 1+sum(hh,2)';
H = max(Ht);
return


function [h,xim] = etaij(P,m,n,i,j)
h = 0;
xim = 0;
dims = repmat([m],[1 i-1]);
for xi = 1:m^(i-1)
  X = ind2coord(xi,dims);
  hX = etaijX(P,m,n,i,j,X);
  %hX = mxetaijX(P,m,n,i,j,X);  % might as well do it fast!
  %h = max(h,hX);
  if hX > h
    h = hX;
    xim = xi;
  end
end
return

function h = etaijX(P,m,n,i,j,X)
h = 0;
for y1 = 1:m
  for y2 = y1+1:m
    P0 = getcondP(P,m,n,[j:n],[1:i],[X y1]);
    P1 = getcondP(P,m,n,[j:n],[1:i],[X y2]);
    h12 = .5*sum(abs(P0-P1));
    h = max(h,h12);
  end
end
return