function myrade(F0,n0)
global F n
F = F0; % size of |A| = F
n = n0; % dimensionality of R^n

radinit(n);
global P0
zz = P0*0 + 1e-12; oo = P0*0+1;

pdim = 2^n;
adim = F*n;
xdim = pdim + adim;
LB = zeros(xdim,1);
LB(1:adim) = LB(1:adim) - 1;
LB(adim+1:end) = zz;
UB = LB*0;
UB(1:adim) = UB(1:adim) + 1;
UB(adim+1:end) = oo;

Aeq0 = zeros(1,xdim);
Aeq0(adim+1:end) = oo';

if 0 % require sym. stationarity
  %Aeq = is_stat(2,n);
  Aeq = is_stat_sym(n);
  Aeq = assembleblock({zeros(0,adim),Aeq});
  Aeq = [Aeq;Aeq0];
  Beq = zeros(size(Aeq,1),1); Beq(end)=1;
else
  %Aeq = oo';
  Aeq = zeros(1,xdim);
  Aeq(adim+1:end) = oo';
  Beq = 1;
end


global A maxr
maxr = 0;
while 1
  A = randn(n,F);
  P = rand(2^n,1);
  P = P/sum(P);
  %%P = ones(2^n,1); P = P + randn(2^n,1)/2^n/10000;
  %P = randbinsym(n);
  
  x = [A(:);P(:)];
  
  y = fmincon(@radrat,x,[],[],Aeq,Beq,LB,UB);
  %%y = fmincon(@radrat,x,Ain,Bin,Aeq,Beq,LB,UB);
  
  %y = fmincon(@radbasic,A(:),[],[],[],[],0*A(:)-N,0*A(:)+N);
  
  fprintf('F=%d n=%d;  maxr = %1.9f \n',F,n,maxr);
end

return



function y = radrat(x)
global F n
pdim = 2^n;
adim = F*n;
xdim = pdim + adim;
A = reshape(x(1:adim),[n F]);
P = x(adim+1:end);
P = makepos(P);

global maxr
r = radave(A,P);
w = A2w(A);
B = sum(w);
r = r/B;

y = -r;


if r > maxr
  global optP optA
  optP = P;
  optA = A;
  maxr = r;
  fprintf('F=%d n=%d;  maxr = %1.9f \n',F,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