function hmmsketchpad1
global A B C n mhid mobs phi Ainf Binf
global mrat
%mhid = 2;
%mobs = 2;
%n = 7;
hmmsketchpad01
return



function hmmsketchpad0
global A B C n mhid mobs phi Ainf Binf
global mrat
myopt = optimset('MaxIter',1000);

mhid = 2;
mobs = 2;
n = 7;

adim = mhid^2;
%bdim = mhid*mobs;
bdim = mhid*mobs*n;
xdim = adim+bdim;

% box constraints
%LB = zeros(adim,1)+1e-10;
%UB = ones(adim,1);
LB = zeros(xdim,1)+1e-10;
UB = ones(xdim,1);


% normalization constraints
Aeqa = zeros(mhid,adim);
for i=1:mhid
  for j=1:mhid
    Aeqa(i,mhid*(i-1)+j) = 1;
  end
end
Beqa = ones(mhid,1);

Aeqb = zeros(mhid,bdim);
rind = 0;
for c=1:mhid
  for i=1:n
    rind = rind+1;
    for r=1:mobs
      ind = coord2ind([r c i],[mobs mhid n]);
      Aeqb(rind,ind) = 1;
    end
  end
end
Beqb = ones(size(Aeqb,1),1);


mrat = 0;

while 1

  h = rand;
  [Aina,Bina] = constrict(mhid,h);
  
  Ain = assembleblock({Aina,zeros(0,bdim)});
  Bin = [Bina];
  Aeq = assembleblock({Aeqa,Aeqb});
  Beq = [Beqa;Beqb];

  A = pnormdim(rand(mhid,mhid),1);
  %A = pnormdim(A+100,1);
  %C = pnormdim(rand(mhid,1   ),1);
  C = pnormdim(ones(mhid,1   ),1);
  %C = C*0; C(2) = 1;
  %B = pnormdim(rand(mobs,mhid),1);
  B = pnormdim(rand(mobs,mhid,n),1);
  %B = eye(mhid);

  x = [A(:); B(:)];
  [x,fval] = fmincon(@myobjeta,x,Ain,Bin,Aeq,Beq,LB,UB,[],myopt);
  fprintf('HMM mrat = %1.9f \n',mrat);

  A = reshape(x(     1:adim     ),[mhid mhid]);
  B = reshape(x(adim+1:adim+bdim),[mobs mhid n]);

end
  
return


function y = myobjeta(x)
global A B C n mhid mobs phi Ainf Binf
global mrat
adim = mhid^2;
%bdim = mhid*mobs;
bdim = mhid*mobs*n;
xdim = adim+bdim;

A = reshape(x(     1:adim     ),[mhid mhid]);
B = reshape(x(adim+1:adim+bdim),[mobs mhid n]);

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

tol = 1e-10;
hh0 = max(hh0,tol);
hh1 = max(hh1,tol);

%good = all(all( hh0<=hh1));
%y = max(max(hh0./hh1));
r = 0;
for i=1:n-1
  for j=i+1:n
    r = max(r,hh0(i,j)/hh1(i,j));
  end
end
y = -r;
if r > mrat
  mrat = r;
  global mA mB mC
  mA=A;mB=B;mC=C;
  fprintf('HMM mrat = %1.9f \n',mrat);
  save HMMDAT mA mB mC mrat
end

return



function hmmsketchpad01
global A B C n mhid mobs phi Ainf Binf
global mrat
myopt = optimset('MaxIter',1000);

load GENPROCDAT01
B = Poh;
[mhid,mobs,n] = size(B);
m=mhid;
T = n-1;

adim = m^2*T;
xdim = adim;

% box constraints
LB = zeros(xdim,1)+1e-10;
UB = ones(xdim,1);

% normalization constraints
Aeq = zeros(m*T,xdim);
for i=1:m*T
  for j=1:m
    Aeq(i,m*(i-1)+j) = 1;
  end
end
Beq = ones(m*T,1);

mrat = 0;

while 1

  A = pnormdim(rand(mhid,mhid,T),1);
  %A = pnormdim(A+100,1);
  %C = pnormdim(rand(mhid,1   ),1);
  C = pnormdim(ones(mhid,1   ),1);
  %C = C*0; C(2) = 1;
  %B = pnormdim(rand(mobs,mhid),1);
  %B = pnormdim(rand(mobs,mhid,n),1);
  %B = eye(mhid);

  x = A(:);
  [x,fval] = fmincon(@myobjeta01,x,[],[],Aeq,Beq,LB,UB,[],myopt);
  fprintf('HMM01 mrat = %1.9f \n',mrat);

  A = reshape(x(     1:adim     ),[mhid mhid T]);
end
  
return


function y = myobjeta01(x)
global A B C n mhid mobs phi Ainf Binf
global mrat
T = n-1;
m=mhid;
adim = m^2*T;
xdim = adim;

A = reshape(x(     1:adim     ),[mhid mhid T]);

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

tol = 1e-10;
hh0 = max(hh0,tol);
hh1 = max(hh1,tol);

%good = all(all( hh0<=hh1));
%y = max(max(hh0./hh1));
r = 0;
for i=1:n-1
  for j=i+1:n
    r = max(r,hh0(i,j)/hh1(i,j));
  end
end
y = -r;
if r > mrat
  mrat = r;
  global mA mB mC
  mA=A;mC=C;
  fprintf('HMM01 mrat = %1.9f \n',mrat);
  save HMMDAT01 mA mC mrat
end

return