function radebasic(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 1 % 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 = radbasic(x)
global F n
A = reshape(x,[n F]);

global P0 maxr

r = radave(A,P0);
B = radbound(A);
r = r/B;


y = -r;

if r > maxr
  global optA
  optA = A;
  maxr = r;
  fprintf('m=%d n=%d;  maxr = %1.9f \n',m,n,maxr);
end


return


function y = radrat(x)
%function y = radrat(P)
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 A P0 maxr
global P0 maxr

r0 = radave(A,P0);
r1 = radave(A,P);

%B = radbound(A);

[hh,Ht,H] = gethhn(P,2,n);

r = r1/r0;
%r = r0/r1;  % unbounded
%r = r1/B;
%r = r0/B;
r = r/H;
y = -r;

%global Asyn
%s = sum(abs(Asyn * P));

%y = -r/(eps+s);
%y = -r + 1e7*s;

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)/n;
R = P'*max((SIGMA * A),[],2)/n; % 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