function gensketchpad(mhid0,mobs0,n0)
global doplots domask mobs mhid n m
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
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 etaijconv
global n mhid mobs Po Ph Poh

while 1
  Ph   = pnormdim(rand(mhid^n,     1),1);
  Poh0 = pnormdim(rand(mobs,mhid),1);
  %Poh1 = pnormdim(rand(mobs,mhid),1);
  %Poh2 = pnormdim(rand(mobs,mhid),1);
  
  for i=1:n
    dims = repmat([mobs],[1 i]);
    for j=i+1:n
      rho1 = pnormdim(rand(mobs,1),1);
      rho2 = pnormdim(rand(mobs,1),1);
      Poh1 = Poh0; Poh1(:,h) = rho1;
      Poh2 = Poh0; Poh2(:,h) = rho2;
      a = rand;
      b = 1-a;
      Pohc = a*Poh1 + b*Poh2;
      Po1 = latentiifillprob(Ph,Poh1,mhid,mobs,n);
      Po2 = latentiifillprob(Ph,Poh2,mhid,mobs,n);
      Poc = latentiifillprob(Ph,Pohc,mhid,mobs,n);
      
      for xi = 1:mobs^i
        X = ind2coord(xi,dims);
        h1 = etaijX(Po1,mobs,n,i,j,X);
        h2 = etaijX(Po2,mobs,n,i,j,X);    
        hc = etaijX(Poc,mobs,n,i,j,X);
        g = (hc < a*h1 + b*h2);
        if ~g
          'not convex'
          keyboard
        end  
      end
    end
  end
end
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
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] = gethhn1(Ph,mobs,n)
  [hh1,Ht,H] = gethhn(Pm,mhid,n)

  good = all(all( hh0<=hh1)) 
  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 = etaij(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 0
  %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);
      r = max(r,rij);
    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 GENPROCDAT 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);

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

pstr = 'latenti htratio';
r = max(Hto(1:end-1)./Hth(1:end-1)/(1+tha));



y = -r;

if r > mrat
  mrat = r;
  mPh = Ph;
  mPoh = Poh;
  save GENPROCDAT 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);
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 GENPROCDAT 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

%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------


%--------------------------------------------------------------
function r = F(x)
global nstate
%global lenseq
if ~x(end), x=x(1:end-1); end;
r = length(find(x~=1));
return
%r = sum(x);
%return
global nstate
%global lenseq
%if ~x(end), x=x(1:end-1); end;
x=x+1;
lenseq = length(x);
global FN
FF = FN{lenseq};
ind = coord2ind(x,(repmat(nstate,1,lenseq)));
r = FF(ind);
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------


%--------------------------------------------------------------
function ind = coord2ind(coords,dims)
% dims is a ROW vector of dimensions
% coords is a matrix of ROW coordinates vectors
% inds is a list of indices
% this is a recent (16/02/2002) change

% compare with matlab's sub2ind routine

if (length(dims)<size(coords,2))
  dims = repmat(dims,[1,size(coords,2)]);
end

cD = cumprod([1,dims(1:end-1)]);
ind = (coords-1) * cD' + 1;

return


function ind = coord2ind_old(coords,dims)
% NOTE: coords and dims are both ROW vectors

cD = cumprod([1,dims(1:end-1)]);
ind = (coords-1) * cD' + 1;

% REALLY OLD VERSION, to illustrate the idea
%ind = 0;
%coords = coords - 1;
%f = 1;
%for ij = 1:length(dims)
%  ind = ind + coords(ij)*f;
%  f = f*dims(ij);
%end
%ind = ind+1;
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------




%--------------------------------------------------------------
function p = P(x,A,C)
x=x+1;
if isempty(x)
  p = 1;
else
  p = C(x(1));
end
for j=2:length(x)
  p = p * A(x(j),x(j-1)); % markov case
end
return

%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function S = Psi(K,m,n)
S = 0;
for t=1:n
  S = S + sum(K.*(K>0));
  K = getK1(K,m,n-t+1);
end
%S = S+pl(K);
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------


%--------------------------------------------------------------
function K1 = getK1(K,m,n)
dims = repmat([m],[1 n]);
dims1 = repmat([m],[1 n-1]);

K1 = zeros(m^(n-1),1);

for xind2n = 1:m^(n-1)
  x2n = ind2coord(xind2n,dims1);
  for x1 = 1:m
    xind = coord2ind([x1 x2n],dims);
    K1(xind2n) = K1(xind2n) + K(xind);
  end
end
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------


%--------------------------------------------------------------
function c = ind2coord(ind,dim)
% dims is a ROW vector of dimensions
% ind is a column vector of indices
% this is a recent (16/02/2002) change

% compare with matlab's ind2sub routine

ind = ind-1;
c = zeros(size(ind,1),size(dim,2));
for ij = 1:length(dim)
  c(:,ij) = mod(ind,dim(ij));
  ind = (ind - c(:,ij))/dim(ij);
end
c = round(c) + 1;
return

function c = ind2coord_old(ind,dim)
ind = ind-1;
c = zeros(size(dim));
for ij = 1:length(dim)
  c(ij) = mod(ind,dim(ij));
  ind = (ind - c(ij))/dim(ij);
end
c = round(c) + 1;
return



%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------





%--------------------------------------------------------------
function A = assembleblock(AA)
% takes a cell array of matrices
% and assembles them into a block matrix

nb = length(AA);
mm = zeros(nb,1);
nn = zeros(nb,1);
for i=1:nb
  mm(i) = size(AA{i},1);
  nn(i) = size(AA{i},2);  
end

A = zeros(sum(mm),sum(nn));
row = 0;
col = 0;
for i=1:nb
  A(row+1:row+mm(i),col+1:col+nn(i)) = AA{i};
  row = row + mm(i);
  col = col + nn(i);
end
return

%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function HH = condH(Poh,mobs,n)
nz = size(Poh,2);
HH = zeros(1,nz);
for zi = 1:nz
  [hh,Ht,H] = gethhn(Poh(:,zi),mobs,n);
  HH(zi) = H;
end
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------


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


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------


%--------------------------------------------------------------
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);
  h = max(h,hX);
end
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------


%--------------------------------------------------------------
function h = etaijX(P,m,n,i,j,X)
P0 = getcondP(P,m,n,[j:n],[1:i  ],[X(1:i  )]);
P1 = getcondP(P,m,n,[j:n],[1:i-1],[X(1:i-1)]);
h = .5*sum(abs(P0-P1));
return


%function Pjj = getcondP(P,m,n,jj,ii,y)
% Pjj(x) = P{ X[jj]=x | X[ii]=y }


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------


%--------------------------------------------------------------
function Pjj = getcondP(P,m,n,jj,ii,y)
% Pjj(x) = P{ X[jj]=x | X[ii]=y }
if ~isempty(intersect(ii,jj))
  error('ii and jj must not intersect');
end

ljj = length(jj);
dimi = repmat([m],[1 ljj]);
Pjj = zeros(m^ljj,1);

for xi = 1:m^ljj
  x = ind2coord(xi,dimi);
  Pjj(xi) = getPX(P,m,n,[y x],[ii jj]);
end
Pjj = Pjj / (sum(Pjj)+realmin);
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------


%--------------------------------------------------------------
function [p,inds] = getPX(P,m,n,X,ii)
% p = P(X(ii) = x)
% works for all distrubutions P on m^n
if length(ii)==0
  p = 1;
  return
end

dims = repmat([m],[1 n]);
v = repmat({1:m},[1 n]);
for j=1:length(ii)
  v{ii(j)} = X(j);
end
xx = list_all_ind_g(v);
inds = coord2ind(xx,dims);

p = sum(P(inds));
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------


%--------------------------------------------------------------
% a generalized "list_all_ind"
% takes a cell array v of vectors of the same object type (though possibly different lengths)
% and returns all the permutations of those objects
% NOTE: v must be a cell array of ROW vectors

function list = list_all_ind_g(v)
list = [];
if ~isempty(v)
  d1 = v{1};

  L1 = list_all_ind_g(v(2:end));
  [m,n] = size(L1);
  m = max(m,1);
  D1 = repmat(d1,[m,1]); D1 = D1(:);
  L2 = repmat(L1,[length(d1),1]);
  list = [D1,L2];

end
return


%  EXAMPLE:
%
% v = 
% 
%     'abc'    'de'    'fg'
%
% >> ind = list_all_ind_g(v)
% 
% ind =
%
% adf
% adg
% aef
% aeg
% bdf
% bdg
% bef
% beg
% cdf
% cdg
% cef
% ceg

%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------









%--------------------------------------------------------------
function [Ef,K] = condexp(P,F,Xt,m,n)
% computes the conditional expectation
% E[f(X)|Xt]
% a numerically stable routine

t = length(Xt);
dimt = repmat([m],[1 n-t]);
dims = repmat([m],[1 n]);
%Snt =  list_all_ind(dimt);
%Sn  =  list_all_ind(dims);

pXt = getPXt(Xt,P,m,n);
Ef = 0;
K = zeros(size(P));
for xint = 1:m^(n-t)
  xtn = ind2coord(xint,dimt);
  x = [Xt xtn];
  xins = coord2ind(x,dims);
  f = F(xins);
  p = P(xins)/(pXt+realmin);
  K(xins) = p;
  Ef = Ef + f*p;
end
return



function p = getPXt(Xt,P,m,n)
t = length(Xt);
if t==0
  p = 1;
  return
end
dimt = repmat([m],[1 n-t]);
dims = repmat([m],[1 n]);
p = 0;
for xint = 1:m^(n-t)
  xtn = ind2coord(xint,dimt);
  x = [Xt xtn];
  xins = coord2ind(x,dims);
  p = p + P(xins);
end
return

% BAD, saved only as an illustration of the 
% indexation discrepancy between list_all_ind
% and ind2coord
function p = getPXt00(Xt,Sn,P,m,n)
t = length(Xt);
if t==0
  p = 1;
  return
end
inds = 1:m^n;
for i=1:t
  indi = find(Sn(:,i)==Xt(i));
  inds = intersect(inds,indi);
end
p = sum(P(inds));
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function [Yhat,K,TV] = condexphat(P,F,Xt0,m,n,xhat)
% computes Yhat with \hat x_t in place of x_t -- an upper bound

t = length(Xt0);
dimt = repmat([m],[1 n-t]);
dims = repmat([m],[1 n]);

tt0 = 1:t;
tt1 = 1:t-1;
Xt1 = Xt0(tt1);

pXt0 = getPXt(Xt0,P,m,n) + realmin;
pXt1 = getPXt(Xt1,P,m,n) + realmin;

Yhat = 0;
K = zeros(m^n,1);
TV = 0;
for xint = 1:m^(n-t)
  xtn = ind2coord(xint,dimt); % xtn = \seq{x}{t+1}{n}
  % first term, simple:
  x = [Xt0 xtn];
  xins = coord2ind(x,dims);
  join = P(xins);
  cond = join / pXt0;
  Yhat = Yhat + cond * F(xins);
  K(xins) = K(xins) + cond;
  c1 = cond;
  % second term, trickier:
  join = 0;
  for xt=1:m
    x = [Xt1 xt xtn];
    xins = coord2ind(x,dims);
    join = join + P(xins);
  end
  cond = join / pXt1;
  x = [Xt1 xhat xtn];
  xins = coord2ind(x,dims);
  Yhat = Yhat - cond * F(xins);
  K(xins) = K(xins) - cond;  
  c2 = cond;
  
  TV = TV + abs(c1-c2);
end
TV = TV/2;
return

function p = getPXt(Xt,P,m,n)
t = length(Xt);
if t==0
  p = 1;
  return
end
dimt = repmat([m],[1 n-t]);
dims = repmat([m],[1 n]);
p = 0;
for xint = 1:m^(n-t)
  xtn = ind2coord(xint,dimt);
  x = [Xt xtn];
  xins = coord2ind(x,dims);
  p = p + P(xins);
end
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function phi = fillgmask(mask,m)
n = length(mask);
phi = zeros(1,m^n);
dims = repmat([m],[1 n]);
for xind = 1:m^n
  x = ind2coord(xind,dims);
  phi(xind) = applymask(mask,x);
end
return


function r = applymask(mask,x)
n = length(x);
r = 0;
for i=1:n
  r = r + ismember(x(i),mask{i});
end
return

%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function [KK,ss] = getkk(K,m,n)
KK = {};
ss = [];
for t=1:n
  KK{t} = K;
  ss(t) = sum(pl(K));
  K = getK1(K,m,n-t+1);
end
return



dims = repmat([m],[1 n]);
dims1 = repmat([m],[1 n-1]);

K1 = zeros(m^(n-1),1);

for xind2n = 1:m^(n-1)
  x2n = ind2coord(xind2n,dims1);
  for x1 = 1:m
    xind = coord2ind([x1 x2n],dims);
    K1(xind2n) = K1(xind2n) + K(xind);
  end
end
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------


%--------------------------------------------------------------
function y=pl(x)
y = max(0,x);
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------




%--------------------------------------------------------------
function H = hh2H(hh)
H = hh2Hi(hh);
return
HH = [];
n = length(hh);
%for i=2:2^n
%  setstr = dec2bin(i-1,n);
%  ii = find(setstr=='1');
for i=1:n
  ii = i:n;
  HH(end+1) = hh2Hi(hh(ii));
end
H = max(HH);
return

function H = hh2Hi(hh)
H = 1;
for i=1:length(hh);
  x = 1;
  for j=1:i
    x = x * hh(j);
  end
  H = H+x;
end
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function P = hmmfillprob1(A0,B0,C0)
global A B C n mhid mobs
A=A0; B=B0; C=C0;
[mobs,mhid,n] = size(B);


dimX = repmat([mobs],[1 n]);
dimS = repmat([mhid],[1 n]);

P = zeros(mobs^n,1);

for xind = 1:mobs^(n)
  x = ind2coord(xind,dimX);
  for sind = 1:mhid^n
    s = ind2coord(sind,dimS);
    P(xind) = P(xind) + Pjoinsx1(s,x);
  end
end

return

function p = Pcondxs1(x,s)
% p = P(x|s)
global A B C n mhid mobs
p = 1;
for i=1:length(x)
  p = p * B(x(i),s(i),i);
end
return


function p = Pjoinsx1(s,x)
% p = P(s,x)
global A B C n mhid mobs
p = C(s(1));
for i=2:length(s)
  p = p * A(s(i),s(i-1),i-1);
  %p = p * A(s(i),s(i-1));
end
p = p * Pcondxs1(x,s);
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function P = latentfillprob(Ph,Poh,mhid,mobs,n)


dimX = repmat([mobs],[1 n]);
dimS = repmat([mhid],[1 n]);

P = zeros(mobs^n,1);

for xind = 1:mobs^(n)
  x = ind2coord(xind,dimX);
  for sind = 1:mhid^n
    s = ind2coord(sind,dimS);
    P(xind) = P(xind) + Ph(sind)*Pcond(Poh,x,s);
  end
end

return

function p = Pcond(Poh,x,s)
% p = P(x|s)
p = 1;
for i=1:length(x)
  p = p * Poh(x(i),s(i),i);
end
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function P = latentiifillprob(Ph,Poh,mhid,mobs,n)
dimX = repmat([mobs],[1 n]);
dimS = repmat([mhid],[1 n]);
P = zeros(mobs^n,1);
for xind = 1:mobs^(n)
  x = ind2coord(xind,dimX);
  for sind = 1:mhid^n
    s = ind2coord(sind,dimS);
    P(xind) = P(xind) + Ph(sind)*Pcond(Poh,x,s);
  end
end
return

function p = Pcond(Poh,x,s)
% p = P(x|s)
p = 1;
for i=1:length(x)
  p = p * Poh(x(i),s(i));
end
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function [Ainf,Binf] = lipab(m,n)
d = 1;
mn = m^n;
L = mn;
Ainf = zeros(L,mn);
Binf = zeros(L,1);
MN = repmat([m],[1 n]);
rind = 0;
for ind1 = 1:mn
  for ind2 = ind1+1:mn
    x = ind2coord(ind1,MN);
    y = ind2coord(ind2,MN);
    if length(find(x~=y))==1  % Hamming dist == 1
      for s = [-1,1]
        rind = rind+1;
        Ainf(rind,ind1) =  s;
        Ainf(rind,ind2) = -s;
        Binf(rind,1) = d;
      end
    end
  end
end
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function P = makepos(P)
P = pnormdim(pl(P)+1e-10,1);
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------


%--------------------------------------------------------------
function A = pnormdim(B,d)
% takes a tensor B of nonnegative numbers
% normalizes B along the dimension d
% such that sum(B,d) = ones(...)


DIM = size(B);

%repm = repmat([1],[1,length(DIM)]);
repm = ones(size(DIM));
repm(d) = DIM(d);

sB = sum(B,d);
nB = (sB <= 0);
sB = sB + nB;  % to prevent divide-by-zero errors
% this is a terrible solution for large sparse matrices!
% use spnormdim

%A=B./repmat(sum(B,d)+realmin,repm);
A=B./repmat(sB,repm);

% the action of this routine if A and B are rank-2 tensors:

%if d == 1
%   A=B./repmat(sum(B,1),[size(B,1),1]);
%else
%   A=B./repmat(sum(B,2),[1,size(B,2)]);
%end

return


function cut
% important! if there is a "row"/"column"/n-dimensional "whatever" of all zeros,
% this routine keeps them as zeros after the normalization!
% here, we use matlab's eps as the "effective zero"
% this is a dirty trick to be avoided in general
% but then again, normalization is extremely unstable when it comes to trying

B2 = B + eps/10;
B1 = B2.*(B2>eps);

A=B1./repmat(sum(B2,d),repm);
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------




%--------------------------------------------------------------
function P = markovfillprob(A0,C0,n0)
global A C n m
A=A0; C=C0; n=n0;
m = size(A,1);
if size(A,3)==1
  A = repmat(A,[1 1 n-1]);
end

dimS = repmat([m],[1 n]);

P = zeros(m^n,1);
for sind = 1:m^n
  s = ind2coord(sind,dimS);
  P(sind) = Ps(s);
end

return

function p = Ps(s)
% p = P(s,x)
global A C n
p = C(s(1));
for i=2:length(s)
  p = p * A(s(i),s(i-1),i-1);
end
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function [phi,fv] = maxphi(K)
global Ainf Binf n
if isempty(Ainf)
    error('Ainf not initialized')
end
[xp,fvalp] = linprog( K,Ainf,Binf,[],[],0*K,0*K+n);
[xn,fvaln] = linprog(-K,Ainf,Binf,[],[],0*K,0*K+n);
if abs(fvalp)>abs(fvaln)
  fv = fvalp;
  x = xp;
else
  fv = fvaln;
  x = xn;
end
phi = x';
fv = abs(phi*K);
return



function phi = maxphi0(K);
global Ainl Binl Ain1 Bin1 GG iif m n mylinopt

[x,fval] = linprog(-K,Ainl,Binl,[],[],0*K,0*K+n,0,mylinopt);
phi = round(x');
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function [h,hh] = neostrict1(P,m,n)

%P = P + 1e-12;
%P = P/sum(P);
h = 0;
hh = [];
for i=2:n
  hh(i-1) = stricti(P,m,n,i);
  %h = max(h,stricti(P,m,n,i));
end
h = max(hh);
return

function h = stricti(P,m,n,i)
%ii = setdiff(1:n,i);
ii = 1:i-1;
lii = length(ii);
PP = [];

dims = repmat([m],[1 lii]);
for xi = 1:m^lii
  X = ind2coord(xi,dims);
  [pX,Xinds] = getPX(P,m,n,X,ii);
  Pi = P(Xinds)/(pX+realmin);
  PP = [PP Pi];
end
h = strict(PP);
return


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



%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function mask = optmask(K,m,n)
%mask = {};
mask = repmat({[]},[1 n]);
for i=1:n
  mu = getmu(K,m,n);
  %phim = fillgmask(mu,m);
  mask{i} = mu;
  K = getK1(K,m,n);
  n = n-1;
end
return


function mu = getmu(K,m,n)
maxmu = [];
maxv = 0;
mask0 = repmat({[]},[1 n]);

for i=1:2^m
  setstr = dec2bin(i-1,m);
  ii1 = find(setstr=='1');
  mask = mask0;
  mask{1} = ii1;
  phim = fillgmask(mask,m);
  fv = phim*K;
  if fv > maxv
    maxv = fv;
    maxmu = ii1;
  end
end
mu = maxmu;
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function D = pairwisedist(m,n)
dims = repmat([m],[1 n]);
N = m^n;
D = zeros(N);
for ii = 1:N
  x = ind2coord(ii,dims);
  for jj = 1:N
    y = ind2coord(jj,dims);
    D(ii,jj) = length(find(x~=y));
  end
end
return

%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function r = randinseg(a,b)
r = rand*(b-a) + a;
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function mask = randmask(m,n)
mask = cell(1,n);
for i=1:n
  mm = find(round(rand(1,m)));
  % avoid triv. masks
  if isempty(mm)
    k = randperm(m);
    mm = k(1);
  end
  if length(mm)==m
    k = randperm(m);
    mm = setdiff(mm,k(1));
  end
  mask{i} = mm;
end
return

%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function hh = spcstrict(P,m,n)
% conditioned on X1=0

hh = [];
for i=2:n
  hh(end+1) = stricti(P,m,n,i);
end
return

function h = stricti(P,m,n,i)
h = 0;

P0 = getcondP(P,m,n,[i:n],[1],[1]);
P1 = getcondP(P,m,n,[i:n],[ ],[ ]);
h = .5*sum(abs(P0-P1));

return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------