function radesketch(m0,n0,F0)
global m n F
m = m0;
n = n0;
F = F0;

global pstr
pstr = 'radesketch'



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

radinit(n);

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

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

% require stationarity
Aeq = is_stat(m,n);
%%Py = pnormdim(ones(m,1),1);
%%[Aeq,Beq] = given_stat(m,n,Py);
Aeq = assembleblock({zeros(0,adim),Aeq});
Aeq = [Aeq;Aeq0];
%Aeq = Aeq0; %don't require stationarity
Beq = zeros(size(Aeq,1),1); Beq(end)=1;


global A maxr
maxr = 0;
while 1
  A = randn(m,F);
  %Py = pnormdim(rand(m,1),1);
  %Py = pnormdim(ones(m,1),1);
  %P = randgivenstat(m,n,Py);
  %P = pnormdim(rand(m^n,1),1);
  P = randbinsym(n);
  
  x = [A(:);P(:)];
  
  y = fmincon(@radrat,x,[],[],Aeq,Beq,LB,UB);
  
  fprintf('%s: n=%d m=%d F=%d; rat = %s/%s;  maxr = %1.9f \n',pstr,n,m,F,lstr,rstr,maxr);
  %fprintf('m=%d n=%d  F=%d;  maxr = %1.9f \n',m,n,F,maxr);
end

return


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

P = makepos(P);

global maxr

%r0 = maxdev(P,A,m,n);
%r = r0/r1;
%r1 = maxdev_sym1(P,A,m,n);
%r2 = maxdev_sym2(P,A,m,n);
%r = r1/r2;

%B = radbound(A);
global P0
%B = radave(A,P0);
%B = radaveF(A,P);

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

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

r = r2 / r1 / H;

y = -r;

if r > maxr
  global optP optA
  optP = P;
  optA = A;
  maxr = r;
  fprintf('%s: n=%d m=%d F=%d; rat = %s/%s;  maxr = %1.9f \n',pstr,n,m,F,lstr,rstr,maxr);
  %fprintf('m=%d n=%d  F=%d;  maxr = %1.9f \n',m,n,F,maxr);
end


return



function B = radbound(A)
[n,N] = size(A);
ma = max(sqrt(sum(A.^2)));
B = ma * sqrt(2*log(N))/n;
return

function [Aeq,Beq] = given_stat(m,n,Py)
mn = m^n;

Aeq = zeros(0,mn);
Beq = zeros(0,1);
rind = 0;
for i=1:n
  % enforce constraint: marg(X_i) = Py
  for y = 1:m
    rind = rind + 1;
    Beq(rind,1) = Py(y);
    for xi = 1:m^n
      x = i2c(xi,m,n);
      if x(i) == y
        Aeq(rind,xi) = 1;
      end
    end
  end
end

return


function c=i2c(i,m,n)
dims = repmat(m,[1 n]);
c = ind2coord(i,dims);
return


function Aeq = is_stat_sym(n)
Aeq = zeros(0,2^n);
rind = 0;
for i=1:n
  rind = rind + 1;
  for xi = 1:2^n
    x = i2c(xi,n);
    
    Aeq(rind,xi) = x(i);
  end
end
return