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

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);
  
  x = [A(:);P(:)];
  
  y = fmincon(@radrat,x,[],[],Aeq,Beq,LB,UB);
  
  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
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

lhs = maxdev(P,A,m,n);

global P0
rhs = maxdev_sym3(P,P0,A,n);

r = lhs / rhs;

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

y = -r;

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


return



function E = maxdev_sym3(P,Psig,calF,n)
% computes max_{f\in F} (1/n) sum_i sig_i(f(X_i)-f(Y_i))
% i.e., second symmetrization
% where f: X -> R
% each column of F is such an f

%f = [1:m]';
%Psig = stat_sym_sig(P,m,n,f);
%Psig = ones(2^n,1); Psig = Psig/sum(Psig);
[m,F] = size(calF);

E = 0;
for xi=1:m^n
  x = i2c(xi,m,n);
  FX = getFX(calF,x);
  for si=1:2^n
    sig = 2*(i2c(si,2,n)-1)-1;
    SIG = repmat(sig,[F 1]);
    Z = max(abs(sum(SIG.*(FX),2)/n));
    E = E + Z * Psig(si) * P(xi);
  end
end

return


function E = maxdev(P,F,m,n)
% computes max_{f\in F} Pf - P_n f
% where f: X -> R
% each column of F is such an f

P1 = getmargX(P,m,n,1); % process is stationary

E = 0;
for xi=1:m^n
  x = i2c(xi,m,n);
  PF = P1' * F;
  
  Fx = getFx(F,x);
  
  Z = max(abs(PF-Fx'));
  E = E + Z * P(xi);
  
end

return

function Fx = getFx(calF,x)
[m,F] = size(calF);
n = length(x);
FX = zeros(F,n);
for f=1:F
  for i=1:n
    FX(f,i) = calF(x(i),f);
  end
end
Fx = mean(FX,2);
return



function FX = getFX(calF,x)
[m,F] = size(calF);
n = length(x);
FX = zeros(F,n);
for f=1:F
  for i=1:n
    FX(f,i) = calF(x(i),f);
  end
end
return


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