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

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


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

global maxr
maxr = 0;

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

while 1
  P = unifP(m,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);
  
  y = fmincon(@myobj,P,[],[],0*P'+1,1,0*P+1e-9,0*P+1);
  
  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;

P = makepos(x);
[hh,Ht,H] = gethhn(P,m,n);


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(-t^2/4/H)/PB;
  
  r = lhs/rhs;
  mr = max(r,mr);
end

r = mr;
y = -r;
global maxr
if r > maxr
  global optP optB optBii
  optP = P;
  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 = talabrute(x,B,m,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