function gensketchpad(mhid0,mobs0,n0)
global doplots domask mobs mhid n m
randomize
mobs = mobs0;
mhid = mhid0;
m = min(mobs,mhid);
n = n0;
%domask = 0;
%sketchmn
%sketchpsi
%sketchdif
%stricheck
%hmmcheck
%latentcheck
latentopti
%etaijconv -- it's not
%Htt1
%varcheck
%mccheck
%trexplore
%chainxplore
%cliquetest
%chaincheck
%reordxplore
return


function reordxplore
global n m P T minpi maxpi
global mrat mT mpsi H0
global longpsi 
longpsi = 0;
mrat = 0;
T = chaingraph(n);
while 1
  % get a random graph
  %k = round(randinseg(2,n));
  %T = randkgraph(n,2,k);
  % get random potentials
  if longpsi
    [nodesi,nodesj]=find(T);
    ii = find(nodesi<nodesj);
    nodesi = nodesi(ii);
    nodesj = nodesj(ii);
    ne = length(nodesi);
    psi = rand(m,m,ne); 
  else
    psi = rand(m); 
  end
  psi = psi/min(psi(:));
  P = treefillprob(m,n,T,psi);
  %[hh,Ht,H0] = gethhn(P,m,n);  
  % see which ordering of the nodes gives greatest stricture
  %[H1,pi1] = maxhrand(P,m,n,100,H0);
  [minH,minpi,maxH,maxpi] = minmaxhrand(P,m,n,100);  
  x = psi(:);
  [x,fval] = fmincon(@reordopt,x,[] ,[]  ,[] ,[]  ,0*x+1,[]);
end
return

function y = reordopt(x)
global n m P T minpi maxpi
global longpsi
x = abs(x);
if longpsi
  [nodesi,nodesj]=find(T);
  ii = find(nodesi<nodesj);
  nodesi = nodesi(ii);
  nodesj = nodesj(ii);
  ne = length(nodesi);
  psi = reshape(x,[m m ne]);
else
  psi = reshape(x,[m m]);
end
P = treefillprob(m,n,T,psi);
Pmin = reorderP(P,m,n,minpi);
[hh,Ht,Hmin] = gethhn(Pmin,m,n);
Pmax = reorderP(P,m,n,maxpi);
[hh,Ht,Hmax] = gethhn(Pmax,m,n);
r = Hmax/Hmin;

y = -r;

global mrat mT mpsi mpi0 mpi1
if r > mrat
  mrat = r;
  mT = T;
  mpsi = psi;
  mpi0 = minpi;
  mpi1 = maxpi;
end
fprintf('ordxplore n=%d   MRAT = %1.9f\n',n,mrat);
return



function chaincheck
global n m P T psi
mrat = 0;
T = chaingraph(n);
while 1
  [nodesi,nodesj]=find(T);
  ii = find(nodesi<nodesj);
  nodesi = nodesi(ii);
  nodesj = nodesj(ii);
  ne = length(nodesi);
  psi = rand(m,m,ne)
  for e=1:ne
    pse = rand(m,m);
    pse = pse/min(pse(:));
    psi(:,:,e) = pse;
  end
  P = treefillprob(m,n,T,psi);  
  for i=1:n-1
    %for x=1:m
    for y1=1:m
      for y2=y1+1:m
        %checkixy(i,x,y1,y2)
        checkiy(i,y1,y2)
      end
    end
    %end
  end
end
return


function p = PSI(x,y,i)
global psi
% for technical reasons, we take x and y to be strings
% and PSI is psi(x(end),y(1))
if (isempty(x) | isempty(y))
  p = 1;
else
  p = psi(x(end),y(1),i);
end
return

function a = axy(i,x,y)
global n m P T psi
a = 0;
for vi=1:m^(i-1)
  v = ind2coord(vi,dimn(i-1));
  for zi=1:m^(n-i-1)
    z = ind2coord(zi,dimn(n-i-1));
    %a = a + pi(v,1,i-2)*PSI(v(1:i-1),y,i-1)*           PSI(x,z,i+1)*pi(z,i+2,n);    
    a =  a + pi(v,1,i-2)*PSI(v(1:i-1),y,i-1)*PSI(y,x,i)*PSI(x,z,i+1)*pi(z,i+2,n);
  end
end
return

function b = bxy(i,x,y)
global n m P T psi
b = 0;
for vi=1:m^(i-1)
  v = ind2coord(vi,dimn(i-1));
  for zi=1:m^(n-i-1)
    z = ind2coord(zi,dimn(n-i-1));
    for x1=1:m
      if x1~=x
        %b = b + pi(v,1,i-2)*PSI(v(1:i-1),y,i-1)*            PSI(x1,z,i+1)*pi(z,i+2,n);
        b =  b + pi(v,1,i-2)*PSI(v(1:i-1),y,i-1)*PSI(y,x1,i)*PSI(x1,z,i+1)*pi(z,i+2,n);      
      end
    end
  end
end
return

function checkiy(i,y,y1)
global n m P T psi
[R,rr] = psi2r(psi);
R = rr(i);

[g0,A,C,dd] = is_markov(P,m,n);
A = A(:,:,i);

s = 0;
for x=1:m
  p0 = A(x,y);
  p1 = A(x,y1);
  if p0>p1
    a0 = axy(i,x,y);
    b0 = bxy(i,x,y);
    c0 = b0/a0;
    %UB = a0*R/(a0*R+b0);
    %UB = R/(R+c0);
    UB = 1/(1+c0);
    %good = (p0<=UB)
    good = abs(p0-UB)<=1e-9
    if ~good
      'ub fails'
      keyboard
    end
    a1 = axy(i,x,y1);
    b1 = bxy(i,x,y1);
    c1 = b1/a1;
    %LB = a1/(a1+b1*R);
    %LB = 1/(1+c1*R);
    LB = 1/(1+c1);
    %good = (p1>=LB)
    good = abs(p1-LB)<=1e-9
    if ~good
      'lb fails'
      keyboard
    end
    
    d = (UB-LB);
    %d0 = R/(R+c0) - 1/(1+c0*R);
    %d1 = R/(R+c1) - 1/(1+c1*R);
    
    if 0
      gmin = (min(d0,d1) >= d);
      gmax = (max(d0,d1) >= d);
      good = gmin & gmax;
      if ~good
        'gminmax fails'
        keyboard
      end
      
      dd = abs(d0-d1)
      if dd > 1e-9
        'dd too big'
        keyboard
      end

    end
    
    %s = s + max(d0,d1);
    s = s + d;
  end
end
%s = s/2;
th = sum(abs(A(:,y)-A(:,y1)))/2;
%good = (th <= s)
good = abs(th-s)<=1e-9
if ~good
  'bound fails'
  keyboard
end
B = (R-1)/(R+1);
good = (s <= B)
if ~good
  'too crude fails'
  keyboard
end


return



function checkixy0(i,x,y,y1)
global n m P T psi
[R,rr] = psi2r(psi);
%R = rr(i);

num1 = 0;
a = 0;
for vi=1:m^(i-1)
  v = ind2coord(vi,dimn(i-1));
  for zi=1:m^(n-i-1)
    z = ind2coord(zi,dimn(n-i-1));
    %num1 = num1 + getPX(P,m,n,[v y x z],1:n);
    %num1 = num1 + pi([v y x z],1,n);
    %num1 = num1 + pi(v,1,i-2)*PSI(v(1:i-1),y,i-1)*PSI(y,x,i)*PSI(x,z,i+1)*pi(z,i+2,n);
    %fprintf('[v y x z] = [%s] \n',num2str([v y x z]));
    % replace psi(y,x) with R
    num1 = num1 + pi(v,1,i-2)*PSI(v(1:i-1),y,i-1)* R *PSI(x,z,i+1)*pi(z,i+2,n);    
    a = a + pi(v,1,i-2)*PSI(v(1:i-1),y,i-1)*PSI(x,z,i+1)*pi(z,i+2,n);    
  end
end
den1 = 0;
b = 0;
for vi=1:m^(i-1)
  v = ind2coord(vi,dimn(i-1));
  for zi=1:m^(n-i-1)
    z = ind2coord(zi,dimn(n-i-1));
    for x1=1:m
      %den1 = den1 + getPX(P,m,n,[v y x1 z],1:n);
      %den1 = den1 + pi([v y x1 z],1,n);
      %fprintf('[v y x1 z] = [%s] \n',num2str([v y x1 z]));      
      %den1 = den1 + pi(v,1,i-2)*PSI(v(1:i-1),y,i-1)*PSI(y,x1,i)*PSI(x1,z,i+1)*pi(z,i+2,n);
      % replace with R+m-1
      den1 = den1 + pi(v,1,i-2)*PSI(v(1:i-1),y,i-1)* (R+1) *PSI(x1,z,i+1)*pi(z,i+2,n);
      if x1~=x
        b = b + pi(v,1,i-2)*PSI(v(1:i-1),y,i-1)*PSI(x1,z,i+1)*pi(z,i+2,n);
      end
    end
  end
end

num2 = 0;
for vi=1:m^(i-1)
  v = ind2coord(vi,dimn(i-1));
  for zi=1:m^(n-i-1)
    z = ind2coord(zi,dimn(n-i-1));
    %num2 = num2 + pi([v y1 x z],1,n);
    num2 = num2 + pi(v,1,i-2)*PSI(v(1:i-1),y1,i-1)*PSI(y1,x,i)*PSI(x,z,i+1)*pi(z,i+2,n);    
  end
end
den2 = 0;
for vi=1:m^(i-1)
  v = ind2coord(vi,dimn(i-1));
  for zi=1:m^(n-i-1)
    z = ind2coord(zi,dimn(n-i-1));
    for x1=1:m
      %den2 = den2 + pi([v y1 x1 z],1,n);
      den2 = den2 + pi(v,1,i-2)*PSI(v(1:i-1),y1,i-1)*PSI(y1,x1,i)*PSI(x1,z,i+1)*pi(z,i+2,n);
      %now make den2 bigger, so num2/den2 is smaller
      %den2 = den2 + pi(v,1,i-2)*PSI(v(1:i-1),y1,i-1)*R*PSI(x1,z,i+1)*pi(z,i+2,n);    
    end
  end
end

%p1 = num1/den1
th1 = num1/den1 - num2/den2;
[g0,A,C,dd] = is_markov(P,m,n);
p1 = A(x,y,i);
if 0
  [a b]
  UB = a*R / (a*R + b);
  %(p1 <= R/(R+(m-1)/R))
  good = (p1 <= UB)
  if ~good
    'noop1'
    keyboard
  end
end
%UB = R/(R+(m-1)/R);
UB = R/(R+(m-1));
good = (p1 <= UB)
if ~good
  'noopU'
  keyboard
end


p2 = A(x,y1,i);
%LB = 1/(1+(m-1)*R);
LB = 1/(1+(m-1)*R);
good = (p2 >= LB)
if ~good
  'noopL'
  keyboard
end

HB = (UB-LB)*m/2;
good = (HB <= (R-1)/(R+1))
if ~good
  'noopR'
  keyboard
end



if 0
  th2 = A(x,y,i) - A(x,y1,i);
  %good = abs(p1-p2)<1e-9
  %good = (p1>=p2)
  %good = abs(th1-th2)<1e-9
  %good = (th1>th2)
  if ~good
    'puncher'
    keyboard
  end
end

if 0
  B = 1-exp(.5*(1-R));
  th1 = abs(th1)
  good = (th1 <= B)
  th1/B
  if ~good
    'too crude'
    keyboard
  end
end

return

function d = dimn(n)
global m
d = repmat([m],[1 n]);
return

function p = pi(u,i,j)
global n m T psi
p = 1;
for k=i:j-1
  p = p * psi(u(k-i+1),u(k-i+2),k);
end
return

function chainxplore
global n m P T pi
global mrat mT mpsi
global longpsi
longpsi = 1;
mrat = 0;
T = chaingraph(n);
while 1
  %pi = randperm(n);
  if longpsi
    [nodesi,nodesj]=find(T);
    ii = find(nodesi<nodesj);
    nodesi = nodesi(ii);
    nodesj = nodesj(ii);
    ne = length(nodesi);
    psi = rand(m,m,ne); 
  else
    psi = rand(m); 
  end
  psi = psi/min(psi(:));
  x = psi(:);
  [x,fval] = fmincon(@chainopt,x,[] ,[]  ,[] ,[]  ,0*x+1,[]);
end
return

function y = chainopt(x)
global n m P T pi
global longpsi
x = abs(x);
if longpsi
  [nodesi,nodesj]=find(T);
  ii = find(nodesi<nodesj);
  nodesi = nodesi(ii);
  nodesj = nodesj(ii);
  ne = length(nodesi);
  psi = reshape(x,[m m ne]);
else
  psi = reshape(x,[m m]);
end
P = treefillprob(m,n,T,psi);
%P = reorderP(P,m,n,pi);
[R,rr] = psi2r(psi);
[b,A,C,dd] = is_markov(P,m,n); 
[h,hh] = stricthh(A);
%th = strict(A);
%th = strict(A(:,:,2));
%r = log(1-th)./(1-R);

%B = (R+1)/2;
BB = (rr-1)./(rr+1);
%[hh,Ht,H] = gethhn(P,m,n);
%r = H/B;
r = max(hh./rr);

y = -r;

global mrat mT mpsi mpi
if r > mrat
  mrat = r;
  mT = T;
  mpsi = psi;
  mpi = pi;
end
fprintf('chainxplore n=%d   MRAT = %1.9f\n',n,mrat);

return




function trexplore
global n m P T pi
global mrat mT mpsi
global longpsi
longpsi = 0;
mrat = 0;
%global PI
%PI = generate_perms(n);
%T = gridgraph(n);
%T = bipart(n);
%T = stargraph(n);
%T = chaingraph(n);
while 1
  %pi = randperm(n);
  %pi = scanord(n);  
  T = randtree(n);
  %k = round(randinseg(2,n));
  %T = randkgraph(n,2,k);
  %T = makecycle(n);
  if longpsi
    [nodesi,nodesj]=find(T);
    ii = find(nodesi<nodesj);
    nodesi = nodesi(ii);
    nodesj = nodesj(ii);    
    ne = length(nodesi);
    psi = rand(m,m,ne); 
  else
    psi = rand(m); 
  end
  psi = psi/min(psi(:));
  x = psi(:);
  [x,fval] = fmincon(@treeopt,x,[] ,[]  ,[] ,[]  ,0*x+1,[]);
end
return

function y = treeopt(x)
global n m P T pi
global longpsi
x = abs(x);
if longpsi
  [nodesi,nodesj]=find(T);
  ii = find(nodesi<nodesj);
  nodesi = nodesi(ii);
  nodesj = nodesj(ii);  
  ne = length(nodesi);
  psi = reshape(x,[m m ne]);
else
  psi = reshape(x,[m m]);
end
P = treefillprob(m,n,T,psi);
%[hh,Ht,H] = gethhn(P,m,n);
%pi = scanord(n);
%P = reorderP(P,m,n,pi);
maxd = max(sum(T));
[hh,Ht,H] = gethhn(P,m,n);
%hval = max(hh(:));
%[H,pi] = maxhord(P,m,n);
%[H,pi] = maxhrand(P,m,n,20,0);


R = psi2r(psi);

%r = H/R/maxd;
%B = ((maxd+1)*R-1)/(n-1);
%r = hval / B;
%r = hval*(n-1)/R^2;
%r = (H-1)/(R-1);
%B = (R+1)/2;
%B = (R+1);
%r = H / B / maxd;
%r = H / B;
%HH = [H0 H1];
%r = max(HH)/min(HH);
%r = H1/H0;
B = (R-1)/(R+1);

%r = maxhi/B/(maxd-1);

if 1
  maxhr = 0;
  for i=1:n-1
    %for j=i+1:n
    %  Bij = B^(j-i)*maxd;
    %maxhr = max(maxhr,hh(i,j)/Bij);
    %end
    maxhr = max(maxhr,hh(i,i+1));
  end
  %r = maxhr;
end


r = maxhr / B / (maxd-1);

%r = max(hh(:))/(R-1);
%r = 0;
%for i=1:n-1
%  mh = max(hh(i,:));
%  ri = mh * (n-i) / (R-1+1e-7);
%  r = max(r,ri)
%end


y = -r;

global mrat mT mpsi mpi
if r > mrat
  mrat = r;
  mT = T;
  mpsi = psi;
  mpi = pi;
end
fprintf('trexplore n=%d   k=%d      MRAT = %1.9f\n',n,maxd,mrat);

return




function hh = getTH(A)
[n,n,T] = size(A);
hh = zeros(1,T);
for t=1:T
  h = 0;
  for i=1:n
    for j=i+1:n
      x = A(:,i,t);
      y = A(:,j,t);
      d = .5*sum(abs(x-y));
      if d>h
        h=d;
      end
    end
  end
  hh(t) = h;
end
return


function varcheck
global n m phi Ainf Binf D mind
[Ainf,Binf] = lipab(m,n);
P = rand(m^n,1); P = P/sum(P);
zz = P*0;
oo = zz+1;
mask = repmat({[1]},[1 n]);
phi = fillgmask(mask,m);
global mrat
mrat = 0;
mind = 1;
D = pairwisedist(m,n);
while 1
  P = rand(m^n,1); P = P/sum(P);
  %[x,fval] = fmincon(@varobj,P,[],[],oo',1,zz+1e-9,(1/2)*oo);
  [x,fval] = fmincon(@varobj,P,[],[],oo',1,zz+1e-9,(1/2)*oo,@mynonl);  
end
return

function [C,Ceq] = mynonl(P)
global n m phi Ainf Binf D mind
Ceq = -1;
C = mind - P'*D*P;
return

function y = varobj(P)
global n m mrat mP mphi 
%pdim = m^n;
%fdim = m^n;
%xdim = pdim + fdim;
%P = x(1:pdim);
%P = makepos(P);
%phi = x(pdim+1:end);

[V,phi] = maxvar(P);
%E = phi * P;
%V = (phi-E).^2 * P;

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

r = H/V;

y = -r;

if r > mrat
  mrat = r;
  mP = P;
  mphi = phi;
end
fprintf('varobj n=%d                            MRAT = %1.9f\n',n,mrat);
return



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

global mrat
mrat = 0;

while 1
  P = rand(m^n,1); P = P/sum(P);
  [x,fval] = fmincon(@htt1obj,P,[],[],oo',1,zz+1e-9,oo);
end
return

function y = htt1obj(P)
global n m mrat mP
[hh,Ht,H] = gethhn(P,m,n);
%r = Ht(2)/Ht(1);
r = hh(1,3)/hh(1,2);

y = -r;

if r > mrat
  mrat = r;
  mP = P;
end
fprintf('MRAT = %1.9f\n',mrat);
return




function latentcheck
global n mhid mobs Po Ph Poh
while 1
  %  cnt = cnt+1;
  
  Ph  = pnormdim(rand(mhid^n,     1),1);
  Poh = pnormdim(rand(mobs^n,mhid^n),1);
  
  if 0
    % give P(h) very low stricture
    Ph = Ph + 1e5; Ph = Ph/sum(Ph);
    % give P(o|h) high stricture
    [a,b] = size(Poh);
    c = min(a,b);
    Poh(1:c,1:c) = Poh(1,c) + diag(diag(Poh(1:c,1:c)))+1e5;
    Poh = pnormdim(Poh,1);
  elseif 0
    % give P(h) high stricture
    ii=randperm(length(Ph));
    ii=ii(1:length(ii)/2);
    Ph(ii)=Ph(ii)+1e5;
    Ph = Ph/sum(Ph);
    % give P(o|h) low stricture
    Poh = Poh + 1e6;
    Poh = pnormdim(Poh,1);
  elseif 0
    % give P(h) high stricture
    ii=randperm(length(Ph));
    ii=ii(1:length(ii)/2);
    Ph(ii)=Ph(ii)+1e5;
    Ph = Ph/sum(Ph);
    % give P(o|h) high stricture
    [a,b] = size(Poh);
    c = min(a,b);
    Poh(1:c,1:c) = Poh(1,c) + diag(diag(Poh(1:c,1:c)))+1e5;
    Poh = pnormdim(Poh,1);
  elseif 0
    % give P(h) low stricture
    Ph = Ph + 1; Ph = Ph/sum(Ph);
    % give P(o|h) low stricture
    Poh = Poh + 1;
    Poh = pnormdim(Poh,1);
  end
  Po  = Poh * Ph;
  
    
  %Poh = pnormdim(rand(mobs,mhid,n),1);
  %Po = latentfillprob(Ph,Poh,mhid,mobs,n);
  
  
  [hh0,Hto,Ho] = gethhn(Po,mobs,n);
  [hh1,Hth,Hh] = gethhn(Ph,mhid,n);
  
  HH = condH(Poh,mobs,n);
  %Hoh = HH*Ph;
  %Hoh = max(HH);
  %good = (Ho < Hh * Hoh)
  %good = (Ho<=Hh)
  %good = all(Hto<=Hth)
  %good = all(all( hh0<=hh1)) % ok, so that's not true
  
  fprintf(' Ho = %1.5f   Hh = %1.5f    Hoh = %1.5f \n',Ho,Hh,Hoh);
  
  if Ho > max(Hh,Hoh)
    'hot dog'
    keyboard
  end
  
  if 0
    for i=1:n
      nui = min(hh1(i,i+1:end));
      for j=i+1:n
        muij = hh0(i,j);
        good = (muij < nui);
        if ~good
          'munu fails'
          keyboard
        end      
      end
    end
 end
  if 0
    'bitin off too much there pardner'
    keyboard
  end
end
return





function hmmcheck
global A B C n mhid mobs Pm Ph
mrat = 0;
while 1
  %  cnt = cnt+1;
  
  A = pnormdim(rand(mhid,mhid,n-1),1);
  C = pnormdim(rand(mhid,1   ,1  ),1);
  B = pnormdim(rand(mobs,mhid,n  ),1);

  Pm = markovfillprob(A,C,n);
  Ph = hmmfillprob1(A,B,C);
  
  [hh0,Ht,H] = gethhn(Ph,mobs,n);
  [hh1,Ht,H] = gethhn(Pm,mhid,n);

  good = all(all( hh0<=hh1));
  rat = max(max(hh0./(hh1+realmin)));
  mrat = max(mrat,rat);
  fprintf('mrat = %1.9f \n',mrat);
  if ~good
    'no good'
    keyboard
  end
end

return


function [hh,Ht,H] = gethhn1(P,m,n)
hh = zeros(n);
for i=1:n
  for j=i+1:n
    h = etaijmin(P,m,n,i,j);
    hh(i,j) = h;
  end
end
Ht = sum(hh,2)';
H = max(Ht);
return


function h = etaij(P,m,n,i,j)
h = 0;
dims = repmat([m],[1 i]);
for xi = 1:m^i
  X = ind2coord(xi,dims);
  hX = etaijX(P,m,n,i,j,X);
  h0 = byhand0(i,j,X);
  g = (abs(hX-h0)<1e-9);
  if ~g
    'h0'
  end  
  h1 = byhand1(i,j,X);
  g = (h0 <= h1);
  if ~g
    'h1'
  end  
  h2 = byhand2(i,j,X);  
  g = (abs(h1-h2)<1e-9);
  if ~g
    'h2'
  end  
  h3 = byhand3(i,j,X);  
  g = (abs(h2-h3)<1e-9);
  if ~g
    'h3'
  end  
  if 0
    
  end
  h = max(h,hX);
end
return



% yes
function h = byhand3(i,j,y)
global n mhid mobs Po Ph Poh
dimz = repmat([mhid],[1 n]);
h = 0;
NUM = 0;
DEN = 0;
for zi = 1:mhid^n
  nuz = Ph(zi);

  if 0
    z = ind2coord(zi,dimz);
    y1z = 1;
    for k=1:i-1
      y1z = y1z * Poh(y(k),z(k),k);
    end
    yz = y1z * Poh(y(i),z(i),i);
  end
  
  yz  = getPX(Poh(:,zi),mobs,n,[y       ],[1:i  ]);
  y1z = getPX(Poh(:,zi),mobs,n,[y(1:i-1)],[1:i-1]);
    
  NUM = NUM + (y1z-yz) * nuz;
  DEN = DEN + yz * nuz;
end
h = h + abs(NUM/DEN);
h = h/2;
return


function h = byhand2(i,j,y)
global n mhid mobs Po Ph Poh
dimz = repmat([mhid],[1 n]);
dimx = repmat([mobs],[1 n-j+1]);
h = 0;
for xi = 1:mobs^(n-j+1)
  x = ind2coord(xi,dimx);
  NUM = 0;
  DEN = 0;
  
  for zi = 1:mhid^n
    nuz = Ph(zi);

    if 0
      z = ind2coord(zi,dimz);
      y1z = 1;
      for k=1:i-1
        y1z = y1z * Poh(y(k),z(k),k);
      end
      yz = y1z * Poh(y(i),z(i),i);
      xyz = yz;
      xy1z = y1z;
      for l=j:n
        xyz = xyz * Poh(x(l-j+1),z(l),l);
        xy1z = xy1z * Poh(x(l-j+1),z(l),l);
      end
    end

    xy1z = getPX(Poh(:,zi),mobs,n,[y(1:i-1) x],[1:i-1, j:n]);
    xyz  = getPX(Poh(:,zi),mobs,n,[y        x],[1:i  , j:n]);
    yz  = getPX(Poh(:,zi),mobs,n,[y       ],[1:i  ]);
    y1z = getPX(Poh(:,zi),mobs,n,[y(1:i-1)],[1:i-1]);
    
    NUM = NUM + (xy1z-xyz) * nuz;
    DEN = DEN + yz * nuz;

  end
  h = h + abs(NUM/DEN);
end
h = h/2;
return


% this also
function h = byhand1(i,j,y)
global n mhid mobs Po Ph Poh
dimz = repmat([mhid],[1 n]);
dimx = repmat([mobs],[1 n-j+1]);
h = 0;
for xi = 1:mobs^(n-j+1)
  x = ind2coord(xi,dimx);
  NUM1 = 0;
  DEN1 = 0;
  NUM2 = 0;
  DEN2 = 0;
  
  for zi = 1:mhid^n
    nuz = Ph(zi);
    
    if 0
      z = ind2coord(zi,dimz);
      y1z = 1;
      for k=1:i-1
        y1z = y1z * Poh(y(k),z(k),k);
      end
      yz = y1z * Poh(y(i),z(i),i);
      xyz = yz;
      xy1z = y1z;
      for l=j:n
        xyz = xyz * Poh(x(l-j+1),z(l),l);
        xy1z = xy1z * Poh(x(l-j+1),z(l),l);
      end
    end
    
    xy1z = getPX(Poh(:,zi),mobs,n,[y(1:i-1) x],[1:i-1, j:n]);
    xyz  = getPX(Poh(:,zi),mobs,n,[y        x],[1:i  , j:n]);
    yz  = getPX(Poh(:,zi),mobs,n,[y       ],[1:i  ]);
    y1z = getPX(Poh(:,zi),mobs,n,[y(1:i-1)],[1:i-1]);
    
    NUM1 = NUM1 + xyz * nuz;
    DEN1 = DEN1 + yz * nuz;

    NUM2 = NUM2 + xy1z * nuz;
    DEN2 = DEN2 + y1z * nuz;
  end
  h = h + abs(NUM1/DEN1 - NUM2/DEN1);
  %h = h + pl(NUM1/DEN1 - NUM2/DEN1);
end
h = h/2;
return



% this works!!
function h = byhand0(i,j,y)
global n mhid mobs Po Ph Poh
dimz = repmat([mhid],[1 n]);
dimx = repmat([mobs],[1 n-j+1]);
h = 0;
for xi = 1:mobs^(n-j+1)
  x = ind2coord(xi,dimx);
  NUM1 = 0;
  DEN1 = 0;
  NUM2 = 0;
  DEN2 = 0;
  
  for zi = 1:mhid^n
    %z = ind2coord(zi,dimz);
    nuz = Ph(zi);
    
    %y1z = 1;
    %for k=1:i-1
    %  y1z = y1z * Poh(y(k),z(k),k);
    %end
    %yz = y1z * Poh(y(i),z(i),i);
    xy1z = getPX(Poh(:,zi),mobs,n,[y(1:i-1) x],[1:i-1, j:n]);
    
    %xyz = yz;
    %xy1z = y1z;
    %for l=j:n
    %  xyz = xyz * Poh(x(l-j+1),z(l),l);
    %  xy1z = xy1z * Poh(x(l-j+1),z(l),l);
    %end
    xyz = getPX(Poh(:,zi),mobs,n,[y x],[1:i, j:n]);
    
    yz  = getPX(Poh(:,zi),mobs,n,[y       ],[1:i  ]);
    y1z = getPX(Poh(:,zi),mobs,n,[y(1:i-1)],[1:i-1]);
    
    NUM1 = NUM1 + xyz * nuz;
    DEN1 = DEN1 + yz * nuz;

    NUM2 = NUM2 + xy1z * nuz;
    DEN2 = DEN2 + y1z * nuz;
  end
  h = h + abs(NUM1/DEN1 - NUM2/DEN2);
  %h = h + pl(NUM1/DEN1 - NUM2/DEN2);
end
h = h/2;
return

function stricheck
global m n mrat mP

m = 2; n=4;

while 1
  P=rand(m^n,1);  P=P/sum(P);
  [Y,K,TV] = condexphat(P,0*P,1,m,n,2);
  [KK,ss] = getkk(K,m,n);
  %[h,hh] = neostrict1(P,m,n);
  %hh = hh*0 + h;
  hh = spcstrict(P,m,n);
  %ch = cumprod(hh);
  %gg = (ch>ss(2:end));
  gg = (abs(ss(2:end)-hh)<1e-10);
  good = all(gg)
  if ~good
    'no sir'
    keyboard
  end
end

return

function sketchpsi
global m n mrat mP
global doplots domask
x = clock;
x = 1/x(end)*1e7;
rand('seed',x)

%m=2
%n=5

mrat = 0;

pdim = m^n;
xdim = pdim;

LB = zeros(xdim,1)+1e-10;
UB = 0*LB;
UB(1:pdim) = UB(1:pdim)+1;

Aeqp = ones(1,pdim);
Beqp = 1;

while 1
  P=rand(m^n,1); P = P/sum(P);
  P = abs(P); P = P/sum(P);
  [P,fv] = fmincon(@myobjpsi,P,[],[],Aeqp,Beqp,LB,UB);
  if abs(fv) > mrat
    mP = P
  end
  fprintf('m=%d n=%d       mrat = %1.9f\n',m,n,mrat);
end      
return


function y = myobjpsi(P)
global m n

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

[Y,K,TV] = condexphat(P,0*P,1,m,n,2);
%[Ef1,K1] = condexp(P,phi,[1],m,n);
%[Ef0,K0] = condexp(P,phi,[ ],m,n);
%K = K1 - K0;
%Y = Ef1 - Ef0;
Y = psinorm(K,m,n);

%[h,hh] = neostrict1(P,m,n); % this used to be the one!!!
[h,hh] = neostrict1(P,m,n);
%[h,hh] = genstrict1(P,m,n);
%H = sum(h.^(0:n-1));
H = hh2H(hh);

y = -abs(Y)/H;

global mrat mP
if abs(y) > mrat
  mrat = abs(y);
  mP = P;
  %updatedat(P,m,n,y)
end
%fprintf('                            mrat = %1.9f\n',mrat);
fprintf('m=%d n=%d       mrat = %1.9f\n',m,n,mrat);
return



function sketchmn
global m n mrat mP mphi
global Ainf Binf
global doplots domask
x = clock;
x = 1/x(end)*1e7;
rand('seed',x)

m=2
n=3


load GENSTR
mrat = MRAT(m,n);


pdim = m^n;
fdim = m^n;
xdim = pdim + fdim;

LB = zeros(xdim,1)+1e-10;
UB = 0*LB;
UB(1:pdim) = UB(1:pdim)+1;
UB(pdim+1:end) = UB(pdim+1:end)+n;


[Ainf,Binf] = lipab(m,n);

Aeqp = ones(1,pdim);
Beqp = 1;

Ain = assembleblock({zeros(0,pdim),Ainf});
Aeq = assembleblock({Aeqp,zeros(0,fdim)});

global mP mphi
while 1

  P=rand(m^n,1); P = P/sum(P);
  mask = randmask(m,n);
  phi = fillgmask(mask,m);
  
  x = [P(:); phi(:)];
  
  [x,fv] = fmincon(@myobjmn,x,Ain,Binf,Aeq,Beqp,LB,UB);
  P = x(1:pdim);
  P = abs(P); P = P/sum(P);
  phi = x(pdim+1:end);
  if abs(fv) > mrat
    mP = P
    mphi = phi
    updatedat(P,m,n,fv)
  end
  fprintf('                            mrat = %1.9f\n',mrat);
end      
    

return


function updatedat(P,m,n,fv)
load GENSTR
mrat = MRAT(m,n);
if abs(fv) > mrat
  mrat = abs(fv);
  MRAT(m,n) = mrat;
  DATA{m,n} = P;
  save GENSTR MRAT DATA
end
return


function y = myobjmn(x)
global m n
pdim = m^n;
fdim = m^n;
xdim = pdim + fdim;

P = x(1:pdim);
P = abs(P); P = P/sum(P);
phi = x(pdim+1:end);

[Y,K,TV] = condexphat(P,phi,1,m,n,2);
%[Ef1,K1] = condexp(P,phi,[1],m,n);
%[Ef0,K0] = condexp(P,phi,[ ],m,n);
%K = K1 - K0;
%Y = Ef1 - Ef0;

[h,hh] = neostrict1(P,m,n);
%[h,hh] = genstrict1(P,m,n);
%H = sum(h.^(0:n-1));
H = hh2H(hh);

y = -abs(Y)/H;

global mrat mP mphi
if abs(y) > mrat
  mrat = abs(y);
  mP = P
  mphi = phi
  updatedat(P,m,n,y)
end
%fprintf('                            mrat = %1.9f\n',mrat);


return



function sketch2
global m n
global Ainf Binf
global doplots domask

m = 2; n = 2;

LB = zeros(8,1);
UB = zeros(8,1); UB(1:4) = UB(1:4)+1; UB(5:end)=UB(5:end)+2;

[Ainf,Binf] = lipab(m,n);
Aeqp = ones(1,4);
Beqp = 1;

Ain = assembleblock({zeros(0,4),Ainf});
Aeq = assembleblock({Aeqp,zeros(0,4)});

global mrat mP mphi
mrat = 0;

while 1

  P=rand(m^n,1); P = P/sum(P);
  mask = randmask(m,n);
  phi = fillgmask(mask,m);
  
  x = [P(:); phi(:)];
  
  [x,fv] = fmincon(@myobj2,x,Ain,Binf,Aeq,Beqp,LB,UB);
  P = x(1:4);
  phi = x(5:end);
  if abs(fv) > mrat
    mrat = abs(fv);
    mP = P
    mphi = phi
  end
  fprintf('                            mrat = %1.9f\n',mrat);
end      
    

return

function y = myobj2(x)
P = x(1:4);
phi = x(5:end);

[Yhat,K,TV] = condexphat(P,phi,1,2,2,2);

P0 = getcondP(P,2,2,2,1 ,1 );
P1 = getcondP(P,2,2,2,[],[]);  
h = .5*sum(abs(P0-P1));

y = -abs(Yhat)/(1+h);

return

function sketch1
global m n
global Ainf Binf
global doplots domask
global Xt
F = zeros(m^n,1);
zz = F*0;
oo = zz+1;


if ~domask
  [Ainf,Binf] = lipab(m,n);
end

global mrat

while 1
  P = rand(m^n,1); P = P/sum(P);
  t = round(randinseg(1,n));
  dimx = repmat([m],[1 t]);
  xtind = round(randinseg(1,m^n));
  Xt = ind2coord(xtind,dimx)
  
  K = randn(size(P));
  [x,fval] = fmincon(@myobj,P,[],[],oo',1,zz,oo);
  %[x,fval] = fminunc(@myobj1,K);
  r = abs(fval);
  P = abs(x);
  K = x;
  if r > mrat
    mrat = r;    
    save GENDATA P r Xt K
  end
end
return
[Yhat,K,TV] = condexphat(P,F,Xt,m,n,1);
[phi,fv1] = maxphi(K);
fv1
%fv2 = sum(abs(K))
fv2 = TV/2;
%good = fv1 < fv2 + 1e-7

return
%[fv,xhat,K] = boundYxhat(P,phi,Xt,m,n)

[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

yy = [];
for xhat = 1:m
  yy(xhat) = condexphat(P,phi,Xt,m,n,xhat);
end
fv1 = max(abs(yy));

fv
fv1
return



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

global mrat

F = zeros(m^n,1);


[Yhat,K,TV] = condexphat(P,F,Xt,m,n,1);
if 0
  [phi,fv1] = maxphi(K);
  %fv1
  fv2 = sum(abs(K));
  %fv2 = TV;
  %good = fv1 < fv2 + 1e-7
  y = -fv1/fv2;
end

%Npsi = psinorm(K);
[phi,Nphi] = maxphi(K);
Nmas = masknorm(K);
%L1 = L1norm(K);

%y = -Npsi/L1;
%y = -Npsi/TV;
y = -Nphi/Nmas;

r = abs(y);
if r > mrat
  mrat = r;
  save GENDATA P r Xt
end
fprintf('                            mrat = %1.9f\n',mrat);
return



function y = psinorm(K,m,n)
y = max(Psi(K,m,n),Psi(-K,m,n));
return

function y = masknorm(K)
global m n
phi1 = fillgmask(optmask( K,m,n),m);
phi2 = fillgmask(optmask(-K,m,n),m);
y = max(abs(phi1*K),abs(phi2*K));
return


function y = L1norm(K)
y = sum(abs(K));
return




function y = myobj1(K)
global m n

global mrat

Npsi = psinorm(K);
L1 = L1norm(K);

y = -Npsi/L1;

r = abs(y);
if r > mrat
  mrat = r;
  save GENDATA K r
end
fprintf('                            mrat = %1.9f\n',mrat);
return


function sketchdif
global m n mrat mP
global doplots domask
x = clock;
x = 1/x(end)*1e7;
rand('seed',x)

%m=2
%n=3
mrat = 0;
pdim = m^n;
xdim = pdim;

LB = zeros(xdim,1)+1e-10;
UB = 0*LB;
UB(1:pdim) = UB(1:pdim)+1;

Aeqp = ones(1,pdim);
Beqp = 1;
while 1
  P=rand(m^n,1); P = P/sum(P);
  [P,fv] = fmincon(@myobjdif,P,[],[],Aeqp,Beqp,LB,UB);
  P = abs(P); P = P/sum(P);
  if abs(fv) > mrat
    mP = P
  end
  fprintf('                            mrat = %1.9f\n',mrat);
end      
return

function y = myobjdif(P)
global m n
P = abs(P); P = P/sum(P);

%P1 = getmargX(P,m,n,[1:n-1]); % marginalize out last component
%[Y,K,TV] = condexphat(P1,0*P1,1,m,n-1,2);

%[Ef1,K1] = condexp(P,phi,[1],m,n);
%[Ef0,K0] = condexp(P,phi,[ ],m,n);
%K = K1 - K0;
%Y = Ef1 - Ef0;
%Y1 = psinorm(K,m,n-1);

[Y,K,TV] = condexphat(P,0*P,1,m,n,2);
Y2 = psinorm(K,m,n);

[h,hh] = neostrict1(P,m,n);
%[h,hh] = genstrict1(P,m,n);
%H = sum(h.^(0:n-1));
%H = hh2H(hh);

delh = prod(hh);
delh = delh + 1e-5;


y = -((Y2-Y1))/delh;

global mrat mP
if abs(y) > mrat
  mrat = abs(y);
  mP = P;
  %updatedat(P,m,n,y)
end
fprintf('                            mrat = %1.9f\n',mrat);
return



function latentoptii
global n mhid mobs Po Ph pstr
hdim = mhid^n;
ohdim = mobs*mhid;
xdim = hdim + ohdim;

LB = zeros(xdim,1);
LB = LB+1e-9;
UB = ones(xdim,1);

Aeqh = ones(1,hdim);
Beqh = 1;
Aeqo = ones(0,ohdim);
rind = 0;
for h = 1:mhid
  rind = rind+1;
  for o = 1:mobs
    cind = coord2ind([o h],[mobs,mhid]);
    Aeqo(rind,cind) = 1;
  end
end
Beqo = ones(size(Aeqo,1),1);

Aeq = assembleblock({Aeqh,Aeqo});
Beq = [Beqh; Beqo];

global mrat
mrat = 0;
while 1
  %  cnt = cnt+1;
  
  Ph  = pnormdim(rand(mhid^n,     1),1);
  Poh = pnormdim(rand(mobs,mhid),1);
  
  x = [Ph(:); Poh(:)];
  [x,fv] = fmincon(@mylatentii,x,[],[],Aeq,Beq,LB,UB);
end
return

function y = mylatentii(x)
global n mhid mobs Po Ph mrat mPh mPoh pstr
hdim = mhid^n;
ohdim = mobs*mhid;
xdim = hdim + ohdim;

Ph = reshape(x(1:hdim),[hdim,1]);
Poh = reshape(x(hdim+1:end),[mobs mhid]);
Ph = makepos(Ph);
Poh = makepos(Poh);

tha = strict(Poh);

Po = latentiifillprob(Ph,Poh,mhid,mobs,n);

[hh0,Hto,Ho] = gethhn(Po,mobs,n);
[hh1,Hth,Hh] = gethhn(Ph,mhid,n);

r = 0;
%pstr = 'htratio';
%pstr = 'hhratio';
%pstr = 'hh13ratio';
pstr = 'mainineq';
%r = max(Hto(1:end-1)./Hth(1:end-1)/(1+tha));
%r = Ho/Hh/(1+tha);
if 1
  %for i=1:n
  for i=1
    %r = max(r,Hto(i)/(1+tha)/Hth(i));
    for j=i+1:n  
    %for j=3
      %r = max(r,hh0(i,j)/tha/(hh1(i,j)+realmin));
      eta_o = hh0(i,j);
      eta_h = hh1(i,j);
      %rij = eta_o / (tha + eta_h + tha*eta_h);
      %rij = eta_o / ( (1+tha)*eta_h);
      if max(eta_o,eta_h) > 1e-10
        %rij = eta_o / ( (1+tha)*eta_h);
        rij = eta_o / eta_h;
        r = max(r,rij);
      end
    end
  end
end
%r = Ho / Hh / tha;
%pstr = 'with tha';
%r = Ho / Hh;
%pstr = 'no tha';

y = -r;


if r > mrat
  mrat = r;
  mPh = Ph;
  mPoh = Poh;
  save GENPROCDAT0 mhid mobs n mPh mPoh mrat Ho Hh tha
end
fprintf('n=%d %s: Ho = %1.5f   Hh = %1.5f        MRAT = %1.9f\n',n,pstr,Ho,Hh,mrat);

% do diagnostics
if 0
  for i=1:n
    nui = min(hh1(i,i+1:end));
    for j=i+1:n
      muij = hh0(i,j);
      good = (muij < nui);
      if ~good
        'munu fails'
        keyboard
      end      
    end
  end
  good = all(all( hh0<=hh1*tha));
  if ~good
    'hh fails'
    keyboard
  end      
end

return



function latentopti
global n mhid mobs Po Ph
hdim = mhid^n;
ohdim = mobs*mhid*n;
xdim = hdim + ohdim;

LB = zeros(xdim,1);
LB = LB+1e-9;
UB = ones(xdim,1);

Aeqh = ones(1,hdim);
Beqh = 1;
Aeqo = ones(0,ohdim);
rind = 0;
for h = 1:mhid
  for i = 1:n
    rind = rind+1;
    for o = 1:mobs
      cind = coord2ind([o h i],[mobs,mhid,n]);
      Aeqo(rind,cind) = 1;
    end
  end
end
Beqo = ones(size(Aeqo,1),1);

Aeq = assembleblock({Aeqh,Aeqo});
Beq = [Beqh; Beqo];

global mrat
mrat = 0;
while 1
  %  cnt = cnt+1;
  
  Ph  = pnormdim(rand(mhid^n,     1),1);
  Poh = pnormdim(rand(mobs,mhid,n),1);
  
  x = [Ph(:); Poh(:)];
  [x,fv] = fmincon(@mylatenti,x,[],[],Aeq,Beq,LB,UB);
end
return

function y = mylatenti(x)
global n mhid mobs Po Ph mrat mPh mPoh
hdim = mhid^n;
ohdim = mobs*mhid*n;
xdim = hdim + ohdim;

Ph = reshape(x(1:hdim),[hdim,1]);
Poh = reshape(x(hdim+1:end),[mobs mhid n]);
Ph = makepos(Ph);
Poh = makepos(Poh);

Po = latentfillprob(Ph,Poh,mhid,mobs,n);

[hh0,Hto,Ho] = gethhn(Po,mobs,n);
[hh1,Hth,Hh] = gethhn(Ph,mhid,n);

tha = strict(Poh);

r = 0;

pstr = 'latenti mainineq 28.7';
%r = Ho / Hh / (1+tha);
% "true"
%r = Ho / Hh;


%pstr = 'latenti htratio';
%r = max(Hto(1:end-1)./Hth(1:end-1)/(1+tha));
% true
% FALSE!! load HTCOUNTEREX

%pstr = 'latenti hhratio';
%r = max(max(hh0./(hh1+realmin)))/(1+tha);
% false!

if 1
  pstr = 'latenti hhii2';
  for i=1:n
    for j=i+2:n
      eta_o = hh0(i,j);
      eta_h = hh1(i,j);
      if max(eta_o,eta_h) > 1e-9
        %rij = eta_o / ( (1+tha)*eta_h);
        rij = eta_o / eta_h;
        r = max(r,rij);
      end
    end
  end
end

y = -r;

if r > mrat
  mrat = r;
  mPh = Ph;
  mPoh = Poh;
  save GENPROCDAT0 mhid mobs n mPh mPoh mrat Ho Hh
end
fprintf('n=%d %s: Ho = %1.5f   Hh = %1.5f        MRAT = %1.9f\n',n,pstr,Ho,Hh,mrat);

if 0
  for t=2:n
    hth = htxx(Ph,mhid,n,t);
    hto = htxx(Po,mobs,n,t);
    %good = all(hto <= (1+tha)*hth);
    good = (max(hto) <= (1+tha)*min(hth));
    if ~good
      'crude fails'
      keyboard
    end
  end
end
return



function latentopt
global n mhid mobs Po Ph

hdim = mhid^n;
ohdim = mobs^n*mhid^n;
xdim = hdim + ohdim;

LB = zeros(xdim,1);
LB = LB+1e-9;
UB = ones(xdim,1);

Aeqh = ones(1,hdim);
Beqh = 1;
Aeqo = ones(0,ohdim);
oo = [1:mobs^n]';
for h = 1:hdim
  ooh = [oo repmat(h,[mobs^n,1])];
  cc = coord2ind(ooh,[mobs^n,mhid^n]);
  Aeqo(h,cc) = ones(1,mobs^n);
end
Beqo = ones(hdim,1);

Aeq = assembleblock({Aeqh,Aeqo});
Beq = [Beqh; Beqo];

global mrat
mrat = 0;
while 1
  %  cnt = cnt+1;
  
  Ph  = pnormdim(rand(mhid^n,     1),1);
  Poh = pnormdim(rand(mobs^n,mhid^n),1);
  
  x = [Ph(:); Poh(:)];
  [x,fv] = fmincon(@mylatent,x,[],[],Aeq,Beq,LB,UB);
end
return


function y = mylatent(x)
global n mhid mobs Po Ph
hdim = mhid^n;
ohdim = mobs^n*mhid^n;
xdim = hdim + ohdim;

Ph = reshape(x(1:hdim),[hdim,1]);
Poh = reshape(x(hdim+1:end),[mobs^n,mhid^n]);

Po  = Poh * Ph;

[hh0,Hto,Ho] = gethhn(Po,mobs,n);
[hh1,Hth,Hh] = gethhn(Ph,mhid,n);

HH = condH(Poh,mobs,n);
%Hoh = HH*Ph;
Hoh = max(HH);

%r = Ho / max(Hh,Hoh); % NO
%r = Ho / (Hh*Hoh);    % NO
r = Ho / (Hh+Hoh);

y = -r;

global mrat mPh mPoh
if r > mrat
  mrat = r;
  mPh = Ph;
  mPoh = Poh;
  save GENPROCDAT0 mhid mobs n mPh mPoh mrat Ho Hh Hoh
end
fprintf('Ho = %1.5f   Hh = %1.5f    Hoh = %1.5f       MRAT = %1.9f\n',Ho,Hh,Hoh,mrat);
return


function h = strict(A)
[n,n,T] = size(A);
h = 0;
for t=1:T
  for i=1:n
    for j=i+1:n
      x = A(:,i,t);
      y = A(:,j,t);
      d = .5*sum(abs(x-y));
      if d>h
        h=d;
      end
    end
  end
end
return




function [v,phi] = maxvar(P)
global m n Ainf Binf
myopt = optimset('Diagnostics','off');
myopt = optimset('Display','off');

LB = 0*P;
UB = LB + n;

V = diag(P) - P*P';
% F' * V * F = Var[F]

[x,fval] = quadprog(-V,zeros(size(P)),Ainf,Binf,[],[],LB,UB,[],myopt);
%[x,fval] = quadprog(-V,zeros(size(P)),[],[],[],[],LB,UB,[],myopt);
v = -fval*2;
phi = x';
return


function cliquetest
global n m P G
global mrat mG mpsi
P = rand(m^n,1); P = P/sum(P);
zz = P*0+1e-9;
oo = zz+1;
mrat = 0;
%global PI
%PI = generate_perms(n);
G = 1-eye(n); % K_n
while 1
  %psi = rand(m^n,1); 
  %psi = psi/min(psi(:));
  %x = psi(:);
  P=rand(m^n,1);  P=P/sum(P);
  %[x,fval] = fmincon(@cliqopt,P,[] ,[]  ,[] ,0*P+1e-9,0*P+1,[]);
  [x,fval] = fmincon(@cliqopt,P,[],[],oo',1,zz+1e-9,oo);    
end
return

function y = cliqopt(x)
global n m P G
global longpsi
%psi = abs(x);
%P = psi / sum(psi);
%maxd = max(sum(G));
P = makepos(x);

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

h = max(hh(:));

R = max(P(:))/min(P(:));

%B = ((maxd+1)*R-1)/(n-1);
%r = h / B;
r = h/(R-1);

y = -r;
global mrat mG mpsi mP
if r > mrat
  mrat = r;
  mG = G;
  %mpsi = psi;
  mP = P;
end
fprintf('cliquetest n=%d      MRAT = %1.9f\n',n,mrat);

return