function strongstatsketch(m0,n0,w0)
global m n w
m = m0;
n = n0;
w = w0;

global COMBS
COMBS = generate_combs(n,w);

P = rand(m^n,1); P = P/sum(P);
zz = P*0 + 1e-12; oo = P*0+1;

if 1 % require stationarity
  Aeq = is_stat(P,m,n);
  Aeq = [Aeq;oo'];
  Beq = zeros(size(Aeq,1),1); Beq(end)=1;
else
  Aeq = oo';
  Beq = 1;
end


global maxH rmin
maxH = 0;
rmin = 1e6;

while 1
  %P = rand(m^n,1); 
  P = ones(m^n,1); P = P + randn(m^n,1)/m^n/1000;
  P = P/sum(P);
  
  %y = fmincon(@myobj0,P,[],[],oo',1,zz,oo,@strong_stat_constr);
  %fprintf(' maxH = %1.9f \n',maxH);

  y = fmincon(@slidw,P,[],[],Aeq,Beq,zz,oo);
  %y = fmincon(@windep,P,[],[],Aeq,Beq,zz,oo);
  fprintf(' rmin = %1.9f \n',rmin);
end

return


function y = windep(P)
% all w-tuples
global m n w
global rmin optP

P = makepos(P);

global COMBS
pH = 1;
for c = 1:size(COMBS,1)
  %for i=1:n
  %  for j=i+1:n
  %Pk = getmargX(P,m,n,[i j]);
  ii = find(COMBS(c,:));
  Pk = getmargX(P,m,n,ii);
  [hh,Ht,H] = gethhn(Pk,m,w);
  pH = pH * H;
  %  end
  %end
end

[hh,Ht,H] = gethhn(P,m,n);
y = pH / H;

if y < rmin
  optP = P;
  rmin = y;
  fprintf(' rmin = %1.9f \n',rmin);
end
return


function y = slidw(P)
% sliding window
global m n w
global rmin optP

P = makepos(P);

pH = 1;
for k=1:n-w+1
  Pk = getmargX(P,m,n,[k:k+w-1]);
  [hh,Ht,H] = gethhn(Pk,m,w);
  pH = pH * H;
end

[hh,Ht,H] = gethhn(P,m,n);
y = pH / H;

if y < rmin
  optP = P;
  rmin = y;
  fprintf(' rmin = %1.9f \n',rmin);
end
return



function y = myobj0(X)
global m n
global maxH maxP
X = makepos(X);
[hhx,Hx,H] = gethhn(X,m,n);
if H > maxH
  maxP = X;
  maxH = H;
  fprintf(' maxH = %1.9f \n',maxH);
end
y = -H;
return