function hmmsketchpad
hmmsketchpad0
return


function [phi,fv] = maxphi(K)
global A B C n mhid mobs phi mask mrat Ainf Binf
[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 = round(x');
fv = abs(phi*K);
return

function hmmsketchpad2
global A B C n mhid mobs phi mask mrat Ainf Binf

mhid = 3;
mobs = 2;
n = 4;


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

%cnt = 0;
while 1
  %  cnt = cnt+1;
  
  A = pnormdim(rand(mhid,mhid,n-1),1);
  %A = pnormdim(rand(mhid,mhid),1); A = repmat(A,[1 1 n-1]);
  C = pnormdim(rand(mhid,1   ,1  ),1);
  B = pnormdim(rand(mobs,mhid,n  ),1);
  B0 = B;
  B1 = pnormdim(rand(mobs,mhid,n  ),1);
  %B = pnormdim(rand(mobs,mhid,1  ),1); B = repmat(B,[1 1 n]);
  
  mask = randmask(mobs,n);
  phi = fillgmask(mask,mobs);
  
  K = hmmphicoef0([1]);
  %[phi,fv0] = maxphi(K);
  %mask = optmask(K,mobs,n);
  %phi = fillgmask(mask,mobs);
  fv0 = phi*K;
  fv0 = abs(fv0)

  %K1 = hmmphicoef1([1]);
  %fv1 = abs(phi*K1);

  fv1 = hmmsum1([1],phi);
  fv1 = abs(fv1)

  h = strict(A);
  Hn = sum(h.^(0:n-1));
  good = fv1 < Hn

  if 0
    %fv0 = hmmsum1([1],phi)
    B = optB
    
    fv1 = hmmsum1([1],phi)
    
    good = abs(fv0) < abs(fv1);  
    
    %h = strict(A);
    %fv0
    %fv1
    %Hn = sum(h.^(0:n-1))
  end  

  if 0
    i = 1;
    %i = round(randinseg(1,n));
    %i = round(randinseg(2,n));
    s = round(randinseg(1,mhid));
    
    lam = rand;
    B = B0;
    fv0 = abs(hmmsum1([1],phi));
    
    %B = B1;
    B = B0;
    B(:,s,i) = B1(:,s,i);
    fv1 = abs(hmmsum1([1],phi));
    
    %B = lam*B0 + (1-lam)*B1;
    B = B0;
    B(:,s,i) = lam*B0(:,s,i) + (1-lam)*B1(:,s,i);
    fvlam = abs(hmmsum1([1],phi));
    
    good = fvlam < lam*fv0 + (1-lam)*fv1 + 1e-9
    %if 2*Hn < fv1
  end


  if ~good
    keyboard
  end
end

return


function oB = optB
global A B C n mhid mobs phi mask mrat Ainf Binf
B0 = B;
oB = B0;
for i=2:n
  for s=1:mhid
    %oB = B0;
    fv = [];
    for x=1:mobs
      oB(:,s,i) = zeros(mobs,1);
      oB(x,s,i) = 1;
      B = oB;
      fv(end+1) = abs(hmmsum1([1],phi));
    end
    [a,x] = max(fv);
    oB(:,s,i) = zeros(mobs,1);
    oB(x,s,i) = 1;
  end
end
B = B0;
return


function K = hmmphicoef0(Xt)
global A B C n mhid mobs
dimX = repmat([mobs],[1 n]);
dimS = repmat([mhid],[1 n]);
K = zeros(mobs^n,1);
t = length(Xt);
Xt1 = Xt(1:t-1);
pXt = forwprob1(A,B,C,Xt);
pXt1 = forwprob1(A,B,C,Xt1);
for xind = 1:mobs^n
  x = ind2coord(xind,dimX);
  xt  = x; xt(1:t) = Xt;
  xt1 = x; xt1(1:t-1) = Xt1;
  tind  = coord2ind(xt, dimX);  
  tind1 = coord2ind(xt1,dimX);
  for sind = 1:mhid^n
    s = ind2coord(sind,dimS);
    pxs = Pcondxs1(x,s);
    psXt = Pjoinsx1(s,Xt) / pXt;
    psXt1 = Pjoinsx1(s,Xt1) / pXt1;    
    K(tind)  = K(tind)  + pxs * psXt;
    K(tind1) = K(tind1) - pxs * psXt1;
  end
end
return


function K = hmmphicoef1(Xt)
global A B C n mhid mobs
dimX = repmat([mobs],[1 n]);
dimS = repmat([mhid],[1 n]);
K = zeros(mobs^n,1);
t = length(Xt);
Xt1 = Xt(1:t-1);
pXt = forwprob1(A,B,C,Xt);
pXt1 = forwprob1(A,B,C,Xt1);
for xind = 1:mobs^n
  x = ind2coord(xind,dimX);
  xt  = x; xt(1:t) = Xt;
  xt1 = x; xt1(1:t-1) = Xt1;
  tind  = coord2ind(xt, dimX);  
  tind1 = coord2ind(xt1,dimX);
  for sind = 1:mhid^n
    s = ind2coord(sind,dimS);
    %pxs = Pcondxs1(x,s);
    pxs = 1;
    psXt = Pjoinsx1(s,Xt) / pXt;
    psXt1 = Pjoinsx1(s,Xt1) / pXt1;    
    K(tind)  = K(tind)  + pxs * psXt;
    K(tind1) = K(tind1) - pxs * psXt1;
  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);
end
p = p * Pcondxs1(x,s);
return


function p = forwprob1(A,B,C,u)
if isempty(u)
  p = 1;
  return
end
T = size(u,2);          % the length of the observation string
N = size(A,1);          % the number of nodes in the graph
a = zeros(N,T);         % the alpha matrix
a(:,1) = B(u(1),:,1)'.*C;
for t = 2:T,
  a(:,t) = (A(:,:,t-1)*a(:,t-1)).*B(u(t),:,t-1)';
end
p = sum(a(:,end));
return


function hmmsketchpad1
global A B C n mhid mobs phi
global mask mrat

mhid = 3;
mobs = 3;
n = 4;

cnt = 0;
while 1
  cnt = cnt+1;

  A = pnormdim(rand(mhid,mhid),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 = eye(mhid);


  mask = randmask(mobs,n);
  %mask = repmat({[1]},[1 n]);
  phi = fillgmask(mask,mobs);
  
  P = fillprobs;
  Ef1 = condexp(P,phi,[1],mobs,n);
  Ef0 = condexp(P,phi,[ ],mobs,n);
  Yhmm = abs(Ef1-Ef0);
  
  
  P = markfillprob(A,C,n);
  Ef1 = condexp(P,phi,[1],mobs,n);
  Ef0 = condexp(P,phi,[ ],mobs,n);
  Ymrk = abs(Ef1-Ef0);

  good = Yhmm < 2*Ymrk + 1e-9;
  if ~good
    'nope'
    keyboard
  end
  fprintf('cnt = %d, so far so good \n',cnt);
end

return

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



mhid = 3;
mobs = 2;
n = 3;

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

% box constraints
%LB = zeros(adim,1)+1e-10;
%UB = ones(adim,1);
LB = zeros(xdim,1)+1e-10;
%LB(adim+1:end) = LB(adim+1:end) - 1e6;  % no box constraints on phi
UB = ones(xdim,1);
UB(adim+bdim+1:end) = UB(adim+bdim+1:end) + n;  % no box constraints on phi


% 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);

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

mrat = 0;

while 1

  %mask = randmask(mobs,n);
  mask = repmat({[1]},[1 n]);
  phi = fillgmask(mask,mobs);
  
  h = rand;
  [Aina,Bina] = constrict(mhid,h);
  
  %Ain = assembleblock({Aina,zeros(0,bdim),Ainf});
  %Bin = [Bina;Binf];
  %Aeq = assembleblock({Aeqa,Aeqb,zeros(0,fdim)});
  %Beq = [Beqa;Beqb];
  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);

  %if ~exist('firsttime','var')
  %  firsttime = 0
  %  load HMMFAIL
  %  A = A(:,:,1)
  %end
  
  %x = [A(:); B(:); phi(:)];
  x = [A(:); B(:)];
  %[x,fval] = fmincon(@myobj,x,Ain,Bin,Aeq,Beq,LB,UB,[],myopt);
  [x,fval] = fmincon(@myobjeta,x,Ain,Bin,Aeq,Beq,LB,UB,[],myopt);
  %[x,fval] = fmincon(@objhmmark,x,Ain,Bin,Aeq,Beq,LB,UB);
  %[x,fval] = fmincon(@mysum,x,Ain,Bin,Aeq,Beq,LB,UB,[],myopt);  


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

end
  
return





function y = myobj(x)
global A B C n mhid mobs phi Ainf Binf
global mask

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

A = reshape(x(     1:adim     ),[mhid mhid]);
A = repmat(A,[1 1 n-1]);
%B = reshape(x(adim+1:adim+bdim),[mobs mhid]);
B = reshape(x(adim+1:adim+bdim),[mobs mhid n]);
%phi = x(adim+bdim+1:end);

%y = -abs(E([1])-E([]));

P = hmmfillprob1(A,B,C);
%P = fillprobs;
[Ef1,K1] = condexp(P,phi,[1],mobs,n);
[Ef0,K0] = condexp(P,phi,[ ],mobs,n);

if 0
  K = K1 - K0;
  [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 = round(x');
  y = -abs(phi*K);

  h = strict(A);
  Hn = sum(h.^(0:n-1));
  global mrat
  r = abs(y)/Hn;
  if r>mrat
    mrat = r;
  end
  fprintf('                   maxr = %1.5f \n',mrat);
else
  y = -abs(Ef1-Ef0);
  
  %global H0
  %h = strict(A);
  %Hn = sum(h.^(0:n-1));
  %fv = boundYxhat(P,phi,[1],mobs,n);
  %if fv > H0 + 1e-9
  %  'too crude'
  %  keyboard
  %end
  
  %y = -abs(hmmsum1([1],phi));
end

return


function y = mysum(x)
global A B C n mhid mobs phi Ainf Binf
global mask

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

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

A = repmat(A,[1 1 n-1]);
B = repmat(B,[1 1 n]);
  
y = -abs(hmmsum1([1],phi));
return



function y = objhmmark(x)
global A B C n mhid mobs phi
global mask

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

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

  P = fillprobs;
  Ef1 = condexp(P,phi,[1],mobs,n);
  Ef0 = condexp(P,phi,[ ],mobs,n);
  Yhmm = abs(Ef1-Ef0);
  
  
  P = markfillprob(A,C,n);
  Ef1 = condexp(P,phi,[1],mobs,n);
  Ef0 = condexp(P,phi,[ ],mobs,n);
  Ymrk = abs(Ef1-Ef0);



y = -Yhmm/Ymrk;
return




%obsolete
function Ef = E(Xt)
% E[f(X)|Xt]
global A B C n mhid mobs

t = length(Xt);
dimX = repmat([mobs],[1 n-t]);
dimS = repmat([mhid],[1 n]);

Ef = 0;

pXt = forwprob(A,B,C,Xt);
for xind = 1:mobs^(n-t)
  x = ind2coord(xind,dimX);
  %xt  = x; xt(1:t) = Xt;
  xt = [Xt x];
  f = F(xt);
  for sind = 1:mhid^n
    s = ind2coord(sind,dimS);
    s1t = s(1:t);
    stn = s(t+1:n);
    SgnXt = Pjoinsx(s,Xt) / pXt;
    XgnS  = Pcondxs(x,stn);
    Ef = Ef + SgnXt * XgnS * f;
    %Ef = Ef + Pcondxs(x,s) * Pjoinsx(s,Xt) / pXt * f;
  end
end
return


function P = fillprobs
global A B C n mhid mobs

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) + Pjoinsx(s,x);
  end
end

return

function p = Pcondxs(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));
end
return


function p = Pjoinsx(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));
end
p = p * Pcondxs(x,s);
return



function r = F(x)
r = length(find(x~=1));
return
global mobs phi
lenseq = length(x);
ind = coord2ind(x,(repmat(mobs,1,lenseq)));
r = phi(ind);
return
global mask
%r = 1;
%return
r = applymask(mask,x);
return

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




function p = forwprob(A,B,C,u)

if isempty(u)
  p = 1;
  return
end

T = size(u,2);          % the length of the observation string
N = size(A,1);          % the number of nodes in the graph
a = zeros(N,T);         % the alpha matrix

a(:,1) = B(u(1),:)'.*C;
for t = 2:T,
  a(:,t) = (A*a(:,t-1)).*B(u(t),:)';
end

p = sum(a(:,end));

return




function [Aina,Bina] = constrict(m,h,T)
% stricture constraints

if ~exist('T','var')
  T = 1;
end

SIGNS = 2*list_all_ind(repmat([2],[1 m]))-3;
Aina = zeros(size(SIGNS,1)*m*(m-1)/2*T,m^2*T);
rind = 0;

for t=1:T
  for i1=1:m
    for i2=i1+1:m
      for s = 1:size(SIGNS,1)
        rind = rind+1;
        ss = SIGNS(s,:);
        for j=1:m
          ind1 = coord2ind([j i1 t],[m m T]);
          ind2 = coord2ind([j i2 t],[m m T]);
          Aina(rind,ind1) =  ss(j);
          Aina(rind,ind2) = -ss(j);
        end
      end
    end
  end
end
Bina = ones(size(Aina,1),1)*h;
Bina = Bina * 2;
return


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