function spectratest(m0,n0)
global m n
randomize
m = m0;
n = n0;

%global myopt
%myopt = optimset('Diagnostics','off');
%myopt = optimset('Display','off');
%warning off
%warning('off');
pstr = mfilename
global pstr


global C
C=ones(m,1)/m
%C=rand(m,1); C=C/sum(C);

pdim = m*m;
xdim = pdim;
LB = zeros(xdim,1);
UB = LB*0;
UB(1:pdim) = 0*UB(1:pdim) + 1;

[Aeq, Beq] = normab(m,m);


global maxr
maxr = 0;

global P X B Xii Bii
Xii = [1:m^n]';
X=ind2coord(Xii,repmat([m],[1 n]));

while 1
  A = pnormdim(rand(m),1);
  P = markfillprob(A,C,n);
  % pick random subset B
  while 1
    BBii = randbigset;
    B = X(Bii,:);
    PB = sum(P(Bii));
    if PB > 1/2-1e-9
      break
    end
  end

  % get enlargement
  %t = randinseg(0,sqrt(n));
  %t = randinseg(0,1);
  %[At,Atii] = getAt(t);
  %t = minimize_t(Bii,t);
  global cBii dd ut
  [cBii,dd]= getcAdd;
  ut = unique(dd);

  
  x = [A(:)];
  y = fmincon(@myobj,x,[],[],Aeq,Beq,LB,UB);
  
  fprintf('%s: n=%d; m=%d; maxr = %1.9f \n',pstr,n,m,maxr);
  
end


return


function y = myobj(x)
global m n
global pstr
global C
pdim = m*m;
%fdim = m;
%xdim = pdim + fdim;
xdim = pdim;

A = makepos(reshape(x(1:pdim),[m m]));
P = markfillprob(A,C,n);
lam2 = lambda2(A);

global cBii dd ut
global X B Xii Bii
  
PB = sum(P(Bii));

mr = 0;
for ti = 1:length(ut)
  t = ut(ti);
  ai = find(dd>=t);
  Atii = cBii(ai);
  
  lhs = sum(P(Atii));
  rhs = exp(-n*(1-lam2)^2*t^2)/PB;
  r = lhs/rhs;
  mr = max(r,mr);
end

r = mr;
y = -r;
global maxr
if r > maxr
  global optA optB optBii
  optA = A;
  optB = optBii;
  maxr = r;
  fprintf('%s: n=%d; m=%d; maxr = %1.9f \n',pstr,n,m,maxr);
end


return




function Bii = randbigset
global m n P X A Xii Bii
Bii = [];
rp = randperm(m^n);
for i=1:length(rp)
  PA = sum(P(Bii));
  if PA >= 1/2
    break
  else
    Bii(end+1) = rp(i);
  end
end
return


function [cBii,dd]= getcAdd
global m n X B Xii Bii
cBii = setdiff(Xii,Bii);
dd = zeros(size(cBii));

for ai = 1:length(cBii)
  xi = cBii(ai);
  x = X(xi,:);
  d = hamset(x,B)/n;
  dd(ai) = d;
end
return


function d = hamset(x,A)
xA = repmat(x,[size(A,1),1]);
xA = (xA~=A);
d = min(sum(xA,2));
return



function [At,Atii] = getAt(t)
global m n X A Xii Bii
Atii0 = setdiff(Xii,Bii);
Atii = [];

for ai = 1:length(Atii0)
  xi = Atii0(ai);
  x = X(xi,:);
  [d,w] = talabrute(x,A,m,n);
  if d > t
    Atii(end+1) = xi;
  end
end
At = X(Atii,:);

return