wolffd@0: function [gamma, loglik, marginals, marginalsT] = bk_ff_fb(prior, transmat, obslik, filter_only, hnodes, ns)
wolffd@0: % BK_FF_FB Fully factored Boyen-Koller version of forwards-backwards
wolffd@0: % [gamma, loglik, marginals, marginalsT] = bk_ff_hmm(prior, transmat, obslik, filter_only, hnodes, ns)
wolffd@0: 
wolffd@0: ss  = length(ns);
wolffd@0: S = length(prior);
wolffd@0: T = size(obslik, 2);
wolffd@0: marginals = cell(ss,T);
wolffd@0: marginalsT = cell(ss,T);
wolffd@0: scale = zeros(1,T);
wolffd@0: alpha = zeros(S, T);
wolffd@0: 
wolffd@0: transmat2 = transmat';
wolffd@0: for t=1:T
wolffd@0:   if t==1
wolffd@0:     [alpha(:,t), scale(t)] = normalise(prior(:) .* obslik(:,t));
wolffd@0:   else
wolffd@0:     [alpha(:,t), scale(t)] = normalise((transmat2 * alpha(:,t-1)) .* obslik(:,t));
wolffd@0:   end
wolffd@0:   [marginals(:,t), marginalsT(:,t)] = project_joint_onto_marginals(alpha(:,t), hnodes, ns);
wolffd@0:   alpha(:,t) = combine_marginals_into_joint(marginalsT(:,t), hnodes, ns);
wolffd@0:   %fprintf('alpha t=%d\n', t);
wolffd@0:   %celldisp(marginals(1:8,t))
wolffd@0: end
wolffd@0: loglik = sum(log(scale));
wolffd@0: 
wolffd@0: if filter_only
wolffd@0:   gamma = alpha;
wolffd@0:   return;
wolffd@0: end
wolffd@0: 
wolffd@0: beta = zeros(S,T);
wolffd@0: gamma = zeros(S,T);
wolffd@0: t = T;
wolffd@0: beta(:,t) = ones(S,1);
wolffd@0: gamma(:,t) = normalise(alpha(:,t) .* beta(:,t));
wolffd@0: [marginals(:,t), marginalsT(:,t)] = project_joint_onto_marginals(gamma(:,t), hnodes, ns);
wolffd@0: 
wolffd@0: for t=T-1:-1:1
wolffd@0:   b = beta(:,t+1) .* obslik(:,t+1); 
wolffd@0:   beta(:,t) = normalise((transmat * b));
wolffd@0:   [junk, tempT] = project_joint_onto_marginals(beta(:,t), hnodes, ns);
wolffd@0:   beta(:,t) = combine_marginals_into_joint(tempT, hnodes, ns);
wolffd@0:   %gamma(:,t) = normalise(alpha(:,t) .* beta(:,t));
wolffd@0:   %[marginals(:,t), marginalsT(:,t)] = project_joint_onto_marginals(gamma(:,t), hnodes, ns);
wolffd@0: end
wolffd@0: 
wolffd@0: gamma2 = zeros(S,T);
wolffd@0: for t=T-1:-1:1
wolffd@0:   b = beta(:,t+1) .* obslik(:,t+1); 
wolffd@0:   xi(:,:,t) = normalise((transmat .* (alpha(:,t) * b')));      
wolffd@0:   if t==T-1
wolffd@0:     gamma2(:,T) = sum(xi(:,:,T-1), 1)';
wolffd@0:   end
wolffd@0:   gamma2(:,t) = sum(xi(:,:,t), 2);
wolffd@0:   [marginals(:,t), marginalsT(:,t)] = project_joint_onto_marginals(gamma2(:,t), hnodes, ns);
wolffd@0: end
wolffd@0: 
wolffd@0: