function [block,LP,LP_best]=chmmdecode(X,Y,T,chmm);% function [block,LP,LP_best]=chmmdecode(X,Y,T,chmm)%% Viterbi and single-state decoding for chmm% X,Y       observations% T         length of each sequence (N must evenly divide by T, default T=N)% chmm       chmm data structure%% block().q_star    maximum posterior state sequence % block().gamma     the posterior: p(q_t=i given X)% block().delta     proby of each previous state: see eq 33a Rabiner (1989)% block().psi       most likely pre-cursor state: see eq 33b Rabiner (1989)% LP                log posterior of model% LP_best           log posterior of best sequencep=length(X(1,:));N=length(X(:,1));tiny=exp(-700);if (rem(N,T)~=0)  disp('Error: Data matrix length must be multiple of sequence length T');  return;end;N=N/T;K=chmm.K;P=chmm.P;Pi=chmm.Pi;alpha=zeros(T,K);beta=zeros(T,K);gamma=zeros(T,K);% Initialise Viterbi bitsdelta=zeros(T,K);psi=zeros(T,K);Gamma=[];Gammasum=zeros(1,K);Xi=zeros(T-1,K*K);likv=zeros(1,N);for n=1:N    Bone = obslike(X,T,n,chmm.hmmone);  Btwo = obslike(Y,T,n,chmm.hmmtwo);	  % Augmenting  for i=1:T,    Btemp=Bone(i,:)'*ones(1,chmm.hmmtwo.K);     Bcompone(i,:)=reshape(Btemp,1,K);    Btemp=ones(1,chmm.hmmone.K)'*Btwo(i,:);     Bcomptwo(i,:)=reshape(Btemp,1,K);  end;  B=[Bcompone.* Bcomptwo];  scale=zeros(T,1);  % Scaling for delta  dscale=zeros(T,1);  alpha(1,:)=Pi(:)'.*B(1,:);  scale(1)=sum(alpha(1,:));   alpha(1,:)=alpha(1,:)/(scale(1)+tiny);  % For viterbi decoding  delta(1,:) = alpha(1,:);    % Eq. 32(a) Rabiner (1989)                              % Eq. 32(b) Psi already zero  for i=2:T    alpha(i,:)=(alpha(i-1,:)*P).*B(i,:);     scale(i)=sum(alpha(i,:));    alpha(i,:)=alpha(i,:)/(scale(i)+tiny);    for k=1:K,      v=delta(i-1,:).*P(:,k)';      mv=max(v);      delta(i,k)=mv*B(i,k);  % Eq 33a Rabiner (1989)      if length(find(v==mv)) > 1        % no unique maximum - so pick one at random	tmp1=find(v==mv);	tmp2=rand(length(tmp1),1);	[tmp3,tmp4]=max(tmp2);	psi(i,k)=tmp4;      else      	psi(i,k)=find(v==mv);  % ARGMAX; Eq 33b Rabiner (1989)      end    end;      % SCALING FOR DELTA ????    dscale(i)=sum(delta(i,:));    delta(i,:)=delta(i,:)/(dscale(i)+tiny);      end;  % Get beta values for single state decoding  beta(T,:)=ones(1,K)/scale(T);  for i=T-1:-1:1    beta(i,:)=(beta(i+1,:).*B(i+1,:))*(P')/scale(i);   end;      % Get gamma values for single state decoding  gamma=(alpha.*beta);   gamma=rdiv(gamma,rsum(gamma));  gammasum=sum(gamma);      xi=zeros(T-1,K*K);  for i=1:T-1    t=P.*( alpha(i,:)' * (beta(i+1,:).*B(i+1,:)));    xi(i,:)=t(:)'/sum(t(:));  end;      likv(n)=sum(log(scale+(scale==0)*tiny));  lik_best(n)=sum(log(dscale+(dscale==0)*tiny));  Gamma=[Gamma; gamma];  Gammasum=Gammasum+gammasum;  block(n).delta=delta;  block(n).gamma=gamma;  block(n).psi=psi;end;% Backtracking for Viterbi decodingfor n=1:N  block(n).q_star (T) = find(block(n).delta(T,:)==max(block(n).delta(T,:))); % Eq 34b Rabiner  for i=T-1:-1:1,    block(n).q_star(i) = block(n).psi(i+1,block(n).q_star(i+1));  endendLP=sum(likv);LP_best=sum(lik_best);