% not numerically stable
function fv = hmmsum1(Xt,phi)
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)

fv = 0;

for sind = 1:mhid^n
  s = ind2coord(sind,dimS);

  xsums = zeros(mobs^n,1);
  xprob = zeros(mobs^n,1);
  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);
    ft  = phi(tind);
    ft1 = phi(tind1);
  
    pxs = Pcondxs1(x,s);
    psXt = Pjoinsx1(s,Xt) / pXt;
    psXt1 = Pjoinsx1(s,Xt1) / pXt1;
    
    %xsums(xind) = pxs * (psXt*ft - psXt1*ft1);
    xsums(xind) = (psXt*ft - psXt1*ft1);
    xprob(xind) = pxs;
  end
  %fv = fv + max(xsums); % nope -- too crude
  fv = fv + xsums'*xprob; 
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


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