wolffd@0: function CPD = update_ess2(CPD, fmarginal, evidence, ns, cnodes, hidden_bitv)
wolffd@0: % UPDATE_ESS Update the Expected Sufficient Statistics of a hhmm Q node.
wolffd@0: % function CPD = update_ess(CPD, fmarginal, evidence, ns, cnodes, idden_bitv)
wolffd@0: 
wolffd@0: % Figure out the node numbers associated with each parent
wolffd@0: dom = fmarginal.domain;
wolffd@0: self = dom(end); % by assumption
wolffd@0: old_self = dom(CPD.old_self_ndx);
wolffd@0: Fself = dom(CPD.Fself_ndx);
wolffd@0: Fbelow = dom(CPD.Fbelow_ndx);
wolffd@0: Qps = dom(CPD.Qps_ndx);
wolffd@0: 
wolffd@0: Qsz = CPD.Qsz;
wolffd@0: Qpsz = CPD.Qpsz;
wolffd@0: 
wolffd@0: 
wolffd@0: fmarg = add_ev_to_dmarginal(fmarginal, evidence, ns);
wolffd@0: 
wolffd@0: 
wolffd@0: 
wolffd@0: % hor_counts(old_self, Qps, self),
wolffd@0: % fmarginal(old_self, Fbelow, Fself, Qps, self)
wolffd@0: % hor_counts(i,k,j) = fmarginal(i,2,1,k,j) % below has finished, self has not
wolffd@0: % ver_counts(i,k,j) = fmarginal(i,2,2,k,j) % below has finished, and so has self (reset)
wolffd@0: % Since any of i,j,k may be observed, we write
wolffd@0: % hor_counts(counts_ndx{:}) = fmarginal(fmarg_ndx{:})
wolffd@0: % where e.g., counts_ndx = {1, ':', 2} if Qps is hidden but we observe old_self=1, self=2.
wolffd@0: % To create this counts_ndx, we write counts_ndx = mk_multi_ndx(3, obs_dim, obs_val)
wolffd@0: % where counts_obs_dim = [1 3], counts_obs_val = [1 2] specifies the values of dimensions 1 and 3.
wolffd@0: 
wolffd@0: counts_obs_dim = [];
wolffd@0: fmarg_obs_dim = [];
wolffd@0: obs_val = []; 
wolffd@0: if hidden_bitv(self)
wolffd@0:   effQsz = Qsz;
wolffd@0: else
wolffd@0:   effQsz = 1;
wolffd@0:   counts_obs_dim = [counts_obs_dim 3];
wolffd@0:   fmarg_obs_dim = [fmarg_obs_dim 5];
wolffd@0:   obs_val = [obs_val evidence{self}];
wolffd@0: end
wolffd@0:   
wolffd@0: % e.g., D=4, d=3, Qps = all Qs above, so dom = [Q3(t-1) F4(t-1) F3(t-1) Q1(t) Q2(t) Q3(t)].
wolffd@0: % so self = Q3(t), old_self = Q3(t-1), CPD.Qps = [1 2], Qps = [Q1(t) Q2(t)]
wolffd@0: dom = fmarginal.domain;
wolffd@0: self = dom(end);
wolffd@0: old_self = dom(1);
wolffd@0: Qps = dom(length(dom)-length(CPD.Qps):end-1);
wolffd@0: 
wolffd@0: Qsz = CPD.Qsizes(CPD.d);
wolffd@0: Qpsz = prod(CPD.Qsizes(CPD.Qps));
wolffd@0: 
wolffd@0: % If some of the Q nodes are observed (which happens during supervised training)
wolffd@0: % the counts will only be non-zero in positions
wolffd@0: % consistent with the evidence. We put the computed marginal responsibilities
wolffd@0: % into the appropriate slots of the big counts array.
wolffd@0: % (Recall that observed discrete nodes only have a single effective value.)
wolffd@0: % (A more general, but much slower, way is to call add_evidence_to_dmarginal.)
wolffd@0: % We assume the F nodes are never observed.
wolffd@0: 
wolffd@0: obs_self = ~hidden_bitv(self);
wolffd@0: obs_Qps = (~isempty(Qps)) & (~any(hidden_bitv(Qps))); % we assume that all or none of the Q parents are observed
wolffd@0: 
wolffd@0: if obs_self
wolffd@0:   self_val = evidence{self};
wolffd@0:   oldself_val = evidence{old_self};
wolffd@0: end
wolffd@0: 
wolffd@0: if obs_Qps
wolffd@0:   Qps_val = subv2ind(Qpsz, cat(1, evidence{Qps}));
wolffd@0:   if Qps_val == 0
wolffd@0:     keyboard
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: if CPD.d==1 % no Qps from above
wolffd@0:   if ~CPD.F1toQ1 % no F from self
wolffd@0:     % marg(Q1(t-1), F2(t-1), Q1(t))                            
wolffd@0:     % F2(t-1) P(Q1(t)=j | Q1(t-1)=i)
wolffd@0:     % 1       delta(i,j)
wolffd@0:     % 2       transprob(i,j)
wolffd@0:     if obs_self
wolffd@0:       hor_counts = zeros(Qsz, Qsz);
wolffd@0:       hor_counts(oldself_val, self_val) = fmarginal.T(2);
wolffd@0:     else
wolffd@0:       marg = reshape(fmarginal.T, [Qsz 2 Qsz]);
wolffd@0:       hor_counts = squeeze(marg(:,2,:));
wolffd@0:     end
wolffd@0:   else
wolffd@0:     % marg(Q1(t-1), F2(t-1), F1(t-1), Q1(t))                            
wolffd@0:     % F2(t-1) F1(t-1)  P(Qd(t)=j| Qd(t-1)=i)
wolffd@0:     % ------------------------------------------------------
wolffd@0:     % 1        1         delta(i,j)
wolffd@0:     % 2        1         transprob(i,j)
wolffd@0:     % 1        2         impossible
wolffd@0:     % 2        2         startprob(j)
wolffd@0:     if obs_self
wolffd@0:       marg = myreshape(fmarginal.T, [1 2 2 1]);
wolffd@0:       hor_counts = zeros(Qsz, Qsz);
wolffd@0:       hor_counts(oldself_val, self_val) = marg(1,2,1,1);
wolffd@0:       ver_counts = zeros(Qsz, 1);
wolffd@0:       %ver_counts(self_val) = marg(1,2,2,1);
wolffd@0:       ver_counts(self_val) = marg(1,2,2,1) + marg(1,1,2,1);
wolffd@0:     else
wolffd@0:       marg = reshape(fmarginal.T, [Qsz 2 2 Qsz]);
wolffd@0:       hor_counts = squeeze(marg(:,2,1,:));
wolffd@0:       %ver_counts = squeeze(sum(marg(:,2,2,:),1)); % sum over i
wolffd@0:       ver_counts = squeeze(sum(marg(:,2,2,:),1)) + squeeze(sum(marg(:,1,2,:),1)); % sum i,b
wolffd@0:     end
wolffd@0:   end % F1toQ1
wolffd@0: else % d ~= 1
wolffd@0:   if CPD.d < CPD.D % general case
wolffd@0:     % marg(Qd(t-1), Fd+1(t-1), Fd(t-1), Qps(t), Qd(t))                            
wolffd@0:     % Fd+1(t-1) Fd(t-1)  P(Qd(t)=j| Qd(t-1)=i, Qps(t)=k)
wolffd@0:     % ------------------------------------------------------
wolffd@0:     % 1        1         delta(i,j)
wolffd@0:     % 2        1         transprob(i,k,j)
wolffd@0:     % 1        2         impossible
wolffd@0:     % 2        2         startprob(k,j)
wolffd@0:     if obs_Qps & obs_self
wolffd@0:       marg = myreshape(fmarginal.T, [1 2 2 1 1]);
wolffd@0:       k = 1;
wolffd@0:       hor_counts = zeros(Qsz, Qpsz, Qsz);
wolffd@0:       hor_counts(oldself_val, Qps_val, self_val) = marg(1, 2,1, k,1);
wolffd@0:       ver_counts = zeros(Qpsz, Qsz);
wolffd@0:       %ver_counts(Qps_val, self_val) = marg(1, 2,2, k,1);
wolffd@0:       ver_counts(Qps_val, self_val) = marg(1, 2,2, k,1) + marg(1, 1,2, k,1);
wolffd@0:     elseif obs_Qps & ~obs_self
wolffd@0:       marg = myreshape(fmarginal.T, [Qsz 2 2 1 Qsz]);
wolffd@0:       k = 1;
wolffd@0:       hor_counts = zeros(Qsz, Qpsz, Qsz);
wolffd@0:       hor_counts(:, Qps_val, :) = marg(:, 2,1, k,:);
wolffd@0:       ver_counts = zeros(Qpsz, Qsz);
wolffd@0:       %ver_counts(Qps_val, :) = sum(marg(:, 2,2, k,:), 1);
wolffd@0:       ver_counts(Qps_val, :) = sum(marg(:, 2,2, k,:), 1) + sum(marg(:, 1,2, k,:), 1);
wolffd@0:     elseif ~obs_Qps & obs_self
wolffd@0:       error('not yet implemented')
wolffd@0:     else % everything is hidden
wolffd@0:       marg = reshape(fmarginal.T, [Qsz 2 2 Qpsz Qsz]);
wolffd@0:       hor_counts = squeeze(marg(:,2,1,:,:)); % i,k,j
wolffd@0:       %ver_counts = squeeze(sum(marg(:,2,2,:,:),1)); % sum over i
wolffd@0:       ver_counts = squeeze(sum(marg(:,2,2,:,:),1)) + squeeze(sum(marg(:,1,2,:,:),1)); % sum over i,b
wolffd@0:     end
wolffd@0:   else % d == D, so no F from below
wolffd@0:     % marg(QD(t-1), FD(t-1), Qps(t), QD(t))                            
wolffd@0:     % FD(t-1) P(QD(t)=j | QD(t-1)=i, Qps(t)=k)
wolffd@0:     % 1      transprob(i,k,j) 
wolffd@0:     % 2      startprob(k,j)
wolffd@0:     if obs_Qps & obs_self
wolffd@0:       marg = myreshape(fmarginal.T, [1 2 1 1]);
wolffd@0:       k = 1;
wolffd@0:       hor_counts = zeros(Qsz, Qpsz, Qsz);
wolffd@0:       hor_counts(oldself_val, Qps_val, self_val) = marg(1, 1, k,1);
wolffd@0:       ver_counts = zeros(Qpsz, Qsz);
wolffd@0:       ver_counts(Qps_val, self_val) = marg(1, 2, k,1);
wolffd@0:     elseif obs_Qps & ~obs_self
wolffd@0:       marg = myreshape(fmarginal.T, [Qsz 2 1 Qsz]);
wolffd@0:       k = 1;
wolffd@0:       hor_counts = zeros(Qsz, Qpsz, Qsz);
wolffd@0:       hor_counts(:, Qps_val, :) = marg(:, 1, k,:);
wolffd@0:       ver_counts = zeros(Qpsz, Qsz);
wolffd@0:       ver_counts(Qps_val, :) = sum(marg(:, 2, k, :), 1);
wolffd@0:     elseif ~obs_Qps & obs_self
wolffd@0:       error('not yet implemented')
wolffd@0:     else % everything is hidden
wolffd@0:       marg = reshape(fmarginal.T, [Qsz 2 Qpsz Qsz]);
wolffd@0:       hor_counts = squeeze(marg(:,1,:,:));
wolffd@0:       ver_counts = squeeze(sum(marg(:,2,:,:),1)); % sum over i
wolffd@0:     end
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: CPD.sub_CPD_trans = update_ess_simple(CPD.sub_CPD_trans, hor_counts);
wolffd@0: 
wolffd@0: if ~isempty(CPD.sub_CPD_start)
wolffd@0:   CPD.sub_CPD_start = update_ess_simple(CPD.sub_CPD_start, ver_counts);
wolffd@0: end
wolffd@0: