function mutinfosketch
global m n
global doplots domask doopt myopt
myopt = optimset('MaxIter',1000);
doplots = 0;
domask = 0;
doopt = 1;
gatherstats(2,2)
return

function gatherstats(m0,n0)
global m n
global VALS BEES ma
global Ainf Binf n
global Xt mrat
global doplots domask doopt myopt
load GENSTR
VALS = []; BEES = [];
n = n0;
m = m0;
mrat = MRAT(m,n);
if ~domask
  [Ainf,Binf] = lipab(m,n);
end
F = zeros(m^n,1);

zz = F*0 + 1e-9;
oo = F*0+1;

    
P=rand(m^n,1);P=P/sum(P);
while 1
  %mi = maxpairmutinfo(P,m,n);
  %mi = genstrict(P,m,n);
  
  %t = round(randinseg(1,n));
  %dimx = repmat([m],[1 t]);
  %xtind = round(randinseg(1,m^n));
  %Xt = ind2coord(xtind,dimx);
  
  if doopt
    P=rand(m^n,1);P=P/sum(P);
    [x,fval] = fmincon(@myobj,P,[],[],oo',1,zz,oo,[],myopt);
    %P = pupdate(P);    
  else
    [Ef1,K1] = condexp(P,F,Xt         ,m,n);
    [Ef0,K0] = condexp(P,F,Xt(1:end-1),m,n);
    K = K1 - K0;
    
    if ~domask
      [phi,fv] = maxphi(K);
    else
      mask = optmask(K,m,n);
      phi = fillgmask(mask,m);
      fv = abs(phi*K);
    end
    
    if doplots
      hold off
      BEES(end+1) = mi;
      VALS(end+1) = fv;
      
      plot(BEES,VALS,'.b');
      
      hold on
      x = 0:.001:max(BEES);
      f = 0;
      for i=0:n-1
        f = f + x.^i;
      end
      %f = 1./(1-x);
      plot(x,f,'.r')
      drawnow
    end
  end
end
return


function y = myobj(P)
global m n
%global Xt
%global Ainf Binf
%global doplots domask

%m = 3;
%n = 3;
global mrat

P = abs(P); P = P/sum(P);

if 0
  F = zeros(m^n,1);
  [Ef1,K1] = condexp(P,F,Xt         ,m,n);
  [Ef0,K0] = condexp(P,F,Xt(1:end-1),m,n);
  K = K1 - K0;
  [phi,fv] = maxphi(K);
  %mi = maxpairmutinfo(P,m,n);
end
%h = genstrict(P,m,n);
h = neostrict1(P,m,n);
H = sum(h.^(0:n));

%[fv,Xt] = maxY(P,m,n);
[fv,Xt] = maxyhat(P,m,n);
%[fv,Xt] = maxYpsi(P,m,n);
  
y = -fv/H;

r = abs(y);
if r > mrat
  load GENSTR
  mrat = max(r,MRAT(m,n));
  MRAT(m,n) = mrat;
  global mP mX
  mP = P; mX = Xt;
  datum.P = P;
  datum.Xt = Xt;
  datum.h = h;
  DATA{m,n} = datum;
  %save GENSTR mP mX m n mrat h H
  save GENSTR MRAT DATA
end
fprintf('                            mrat = %1.9f\n',mrat);
return


function P = pupdate(P)
global m n
good = 0;
del = 1e-5;
P0 = P;
fv0 = strictrat(P);
while ((~good) | (del > 1e-9))
  DP = Ndiff('strictrat',1,{P});
  P1 = P + del*DP(:);
  P1 = pnormdim(abs(P1),1);
  fv1 = strictrat(P1);
  if fv1 > fv0
    good = 1;
  end
  del = del/2;
end
good
if good
  P = P1;
else
  P=rand(m^n,1);P=P/sum(P);
end
return


function r = strictrat(P)
global m n
global mrat

h = neostrict1(P,m,n);
H = sum(h.^(0:n));

%[fv,Xt] = maxY(P,m,n);
[fv,Xt] = maxyhat(P,m,n);
%[fv,Xt] = maxYpsi(P,m,n);
  
r = fv/H;

if r > mrat
  load GENSTR
  mrat = max(r,MRAT(m,n));
  MRAT(m,n) = mrat;
  global mP mX
  mP = P; mX = Xt;
  datum.P = P;
  datum.Xt = Xt;
  datum.h = h;
  DATA{m,n} = datum;
  %save GENSTR mP mX m n mrat h H
  save GENSTR MRAT DATA
end
fprintf('                            mrat = %1.9f\n',mrat);

return