% this agrees with the generic hmmphicoef derived from
% condexp -- that's good!

function K = hmmphicoef0(A0,B0,C0,n0,Xt)
global A B C n mhid mobs
A = A0; B = B0; C = C0; n = n0;
mhid = size(A,1);
mobs = size(B,1);

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

K = zeros(mobs^n,1);
t = length(Xt);
Xt1 = Xt(1:t-1);
pXt = forwprob(A,B,C,Xt);
pXt1 = forwprob(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 = Pcondxs(x,s);
    psXt = Pjoinsx(s,Xt) / pXt;
    psXt1 = Pjoinsx(s,Xt1) / pXt1;    
    
    K(tind)  = K(tind)  + pxs * psXt;
    K(tind1) = K(tind1) - pxs * psXt1;
  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 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
b = zeros(N,T);         % the beta 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