function myrade3(n0,lb)
randomize
global n
n = n0; % dimensionality of R^n

global P0 P1
P0 = ones(2^n,1); P0 = P0/sum(P0);
zz = P0*0 + lb; oo = P0*0+1;


  %Aeq = is_stat(2,n);
  Aeq = is_stat_sym(n);
  Aeq = [Aeq;oo'];
  Beq = zeros(size(Aeq,1),1); Beq(end)=1;
  
  Aeq = assembleblock({Aeq,Aeq});
  Beq = [Beq; Beq];
  
zz = [zz;zz];
oo = [oo;oo];
  
global maxr
maxr = 0;
while 1
  P0 = randbinsym(n);
  P1 = randbinsym(n);
  
  x = [P0(:);P1(:)];
  
  y = fmincon(@radrat,x,[],[],Aeq,Beq,zz,oo);
  
  fprintf('n=%d;  maxr = %1.9f \n',n,maxr);
end

return



function y = radrat(PP)
global n maxr
P0 = makepos(PP(1:2^n));
P1 = makepos(PP(2^n+1:end));


[hh,Ht,H0] = gethhn(P0,2,n);
[hh,Ht,H1] = gethhn(P1,2,n);

v1 = Psi(P0-P1,2,n);
%v2 = Ht(1)+1;
v2 = H0*H1;
%v2 = 1+sum(max(hh));
r = v1/v2;

y = -r;

if r > maxr
  global optP
  optP = P1;
  maxr = r;
  fprintf('n=%d;  maxr = %1.9f \n',n,maxr);
end


return

function R = radave(A,P)
global SIGMA
[n,m] = size(A);
R = P'*max(abs(SIGMA * A),[],2);
%R = P'*max((SIGMA * A),[],2); % w/o abs it's zero!!!
return



function B = radbound(A)
[n,N] = size(A);
ma = max(sqrt(sum(A.^2)));
B = ma * sqrt(2*log(N))/n;
return

function w = A2w(A)
%w = max(sum(abs(A)));
w = max(abs(A),[],2);
return
[n,N] = size(A);
w = zeros(n,1);
A = abs(A);
for i=1:N
    for j=i+1:N
        a = A(:,i);
        b = A(:,j);
        %d = abs(a-b);
        d = a+b;
        w = max(w,d);
    end
end
return