function phimixcheck
global m n i j y w
m = 2; n=4;
global maxr maxP maxi maxj maxy maxw
maxr = 0;

if 0
  y=maxy
  w=maxw
  i=maxi
  j=maxj
  P=maxP  
end

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

while 1

  P = rand(m^n,1); P = P/sum(P);
  for i=1:n-1
    for j=i+1:n

      for yi = 1:m^(i)
        y = ind2coord(yi,repmat(m,1,i));
        for w = 1:m
          [x,fval] = fmincon(@myobj,P,[],[],oo',1,zz+1e-5,oo);
          fprintf('maxr = %1.9f \n',maxr);  
        end
      end
      
    end
  end

end

return


function phimixcheck2
global m n i j y w
m = 2; n=4;
global maxr maxP maxA
maxr = 0;

T = n-1;

% box constraints
LB = zeros(m^2*T,1)+1e-10;
UB = ones(m^2*T,1);
% normalization constraints
Aeq = zeros(m*T,m^2*T);
for i=1:m*T
  for j=1:m
    Aeq(i,m*(i-1)+j) = 1;
  end
end
Beq = ones(m*T,1);

while 1

  A = pnormdim(rand(m,m,T),1);
  [x,fval] = fmincon(@myobj1,A(:),[],[],Aeq,Beq,LB,UB);
  fprintf('maxr = %1.9f \n',maxr);  
  
end

return


function val = myobj1(x)
global m n i j y w
global maxr maxP maxA

T = n-1;

A = reshape(x,[m m T]);
C = pnormdim(ones(m,1),1);

P = markfillprob(A,C,n);

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

% to show phi-mixing & not eta-mixing
% need PHI(k) small and H big
r = (H-1) / PHI(end);

if r > maxr
  maxr = r;
  %maxP = P;
  maxA = A;
  fprintf('maxr = %1.9f \n',maxr);
end
val = -r;

return







function phimixcheck0
global m n
m = 2; n=4;
global maxr maxP maxi maxj maxy maxw
maxr = 0;

if 0
  y=maxy
  w=maxw
  i=maxi
  j=maxj
  P=maxP  
end

while 1

  P = rand(m^n,1); P = P/sum(P);
  Eta = 0;
  Phi = 0;
  for i=1:n-1
    for j=i+1:n
      for yi = 1:m^(i)
        y = ind2coord(yi,repmat(m,1,i));
        for w = 1:m
          P0 = getcondP(P,m,n,[j:n],[1:i],[y(1:i  )  ]);
          P1 = getcondP(P,m,n,[j:n],[1:i],[y(1:i-1) w]);
          eta = .5*sum(abs(P0-P1));
          eta = eta + eps;
          Eta = max(Eta,eta);
        end
        P2 = getmargX(P,m,n,j:n);
        phi0 = .5*sum(abs(P0-P2));
        phi1 = .5*sum(abs(P1-P2));
        phi = max(phi0,phi1);
        phi = phi + eps;
        Phi = max(Phi,phi);
      end
      
      r = Eta / Phi;
      
      if r > maxr
        maxr = r;
        maxP = P;
        maxi = i;
        maxj = j;
        maxy = y;
        maxw = w;
        fprintf('maxr = %1.9f \n',maxr);
      end
    
    end
  end
  
end

return



function val = myobj(P)
global m n i j y w
global maxr maxP maxi maxj maxy maxw
Eta = 0;
Phi = 0;

%for yi = 1:m^(i)
%  y = ind2coord(yi,repmat(m,1,i));
%  for w = 1:m
    
    P0 = getcondP(P,m,n,[j:n],[1:i],[y(1:i  )  ]);
    P1 = getcondP(P,m,n,[j:n],[1:i],[y(1:i-1) w]);
    eta = .5*sum(abs(P0-P1));
    eta = eta + eps;
    Eta = max(Eta,eta);
    
    P2 = getmargX(P,m,n,j:n);
    %phi = .5*sum(abs(P0-P2));
    phi0 = .5*sum(abs(P0-P2));
    phi1 = .5*sum(abs(P1-P2));
    %phi = max(phi0,phi1);
    phi = phi0+phi1;
    phi = phi + eps;
    Phi = max(Phi,phi);
    
%  end
%end

r = Eta / Phi;
%r = Phi / Eta;

if r > maxr
  maxr = r;
  maxP = P;
  maxi = i;
  maxj = j;
  maxy = y;
  maxw = w;
  fprintf('maxr = %1.9f \n',maxr);
end
val = -r;

return