function gensketchpad1(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 d = dimn(n)
global m
d = repmat([m],[1 n]);
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 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 = max(Hto(1:end-1)./Hth(1:end-1));
%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);
      if max(eta_o,eta_h) > 1e-6
        %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 02.8';
%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!
%r = max(Hto(1:end-1)./Hth(1:end-1));

if 0
  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 GENPROCDAT2 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 [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


function [h,xim] = etaij(P,m,n,i,j)
global A B C mobs mhid hhat

h = 0;
xim = 0;
dims = repmat([m],[1 i-1]);
for xi = 1:m^(i-1)
  X = ind2coord(xi,dims);
  hX = etaijX0(P,m,n,i,j,X);
  if hX > h
    h = hX;
    xim = xi;
  end
end
return



function h = etaijX0(P,m,n,i,j,X)
h = 0;
for y1 = 1:m
  for y2 = y1+1:m
    P0 = getcondP(P,m,n,[j:n],[1:i],[X y1]);
    P1 = getcondP(P,m,n,[j:n],[1:i],[X y2]);
    h12 = .5*sum(abs(P0-P1));
    h = max(h,h12);
  end
end
return