wolffd@0: function CPD = update_ess(CPD, fmarginal, evidence, ns, cnodes, hidden_bitv)
wolffd@0: % UPDATE_ESS Update the Expected Sufficient Statistics of a CPD (MLP)
wolffd@0: % CPD = update_ess(CPD, family_marginal, evidence, node_sizes, cnodes, hidden_bitv)
wolffd@0: %
wolffd@0: % fmarginal = overall posterior distribution of self and its parents
wolffd@0: % fmarginal(i1,i2...,ik,s)=prob(Pa1=i1,...,Pak=ik, self=s| X)
wolffd@0: % 
wolffd@0: % => 1) prob(self|Pa1,...,Pak)=fmarginal/prob(Pa1,...,Pak) with prob(Pa1,...,Pak)=sum{s,fmarginal}
wolffd@0: %       [self estimation -> CPD.self_vals]
wolffd@0: % 	  2) prob(Pa1,...,Pak) [SCG weights -> CPD.eso_weights]
wolffd@0: %
wolffd@0: % Hidden_bitv is ignored
wolffd@0: 
wolffd@0: % Written by Pierpaolo Brutti
wolffd@0: 
wolffd@0: if ~adjustable_CPD(CPD), return; end
wolffd@0: 
wolffd@0: dom = fmarginal.domain;                              
wolffd@0: cdom = myintersect(dom, cnodes);                     
wolffd@0: assert(~any(isemptycell(evidence(cdom))));           
wolffd@0: ns(cdom)=1;
wolffd@0: 
wolffd@0: self = dom(end);                                  
wolffd@0: ps=dom(1:end-1);                                     
wolffd@0: dpdom=mysetdiff(ps,cdom);                            
wolffd@0: 
wolffd@0: dnodes = mysetdiff(1:length(ns), cnodes);            
wolffd@0: 
wolffd@0: ddom = myintersect(ps, dnodes);                      %
wolffd@0: if isempty(evidence{self}),                          % if self is hidden in what follow we must 
wolffd@0:     ddom = myintersect(dom, dnodes);                 % consider its dimension
wolffd@0: end                                                  % 
wolffd@0: 
wolffd@0: odom = dom(~isemptycell(evidence(dom)));    
wolffd@0: hdom = dom(isemptycell(evidence(dom)));              % hidden parents in domain
wolffd@0:  
wolffd@0: dobs = myintersect(ddom, odom);             
wolffd@0: dvals = cat(1, evidence{dobs});             
wolffd@0: ens = ns;                                            % effective node sizes              
wolffd@0: ens(dobs) = 1;                              
wolffd@0:                                             
wolffd@0: dpsz=prod(ns(dpdom));
wolffd@0: S=prod(ens(ddom));
wolffd@0: subs = ind2subv(ens(ddom), 1:S);
wolffd@0: mask = find_equiv_posns(dobs, ddom);
wolffd@0: for i=1:length(mask),
wolffd@0:     subs(:,mask(i)) = dvals(i);
wolffd@0: end
wolffd@0: supportedQs = subv2ind(ns(ddom), subs);
wolffd@0: 
wolffd@0: Qarity = prod(ns(ddom));
wolffd@0: if isempty(ddom),                      
wolffd@0:   Qarity = 1;                         
wolffd@0: end                                
wolffd@0: fullm.T = zeros(Qarity, 1);
wolffd@0: fullm.T(supportedQs) = fmarginal.T(:);
wolffd@0: 
wolffd@0: % For dynamic (recurrent) net-------------------------------------------------------------
wolffd@0: % ----------------------------------------------------------------------------------------
wolffd@0: high=size(evidence,1);                                  % slice height
wolffd@0: ss_ns=ns(1:high);                                       % single slice nodes sizes
wolffd@0: pos=self;                                               %
wolffd@0: slice_num=0;                                            %
wolffd@0: while pos>high,                                         % 
wolffd@0:     slice_num=slice_num+1;                              % find active slice
wolffd@0:     pos=pos-high;                                       % pos=self posistion into a single slice
wolffd@0: end                                                     %
wolffd@0: 
wolffd@0: last_dim=pos-1;                                         % 
wolffd@0: if isempty(evidence{self}),                             % 
wolffd@0:     last_dim=pos;                                       %
wolffd@0: end                                                     % last_dim=last reshaping dimension      
wolffd@0: reg=dom-slice_num*high;
wolffd@0: dex=myintersect(reg(find(reg>=0)), [1:last_dim]);       %           
wolffd@0: rs_dim=ss_ns(dex);                                      % reshaping dimensions
wolffd@0: 
wolffd@0: if slice_num>0,
wolffd@0:     act_slice=[]; past_ancest=[];                       %
wolffd@0:     act_slice=slice_num*high+[1:high];                  % recover the active slice nodes
wolffd@0:     % past_ancest=mysetdiff(ddom, act_slice);
wolffd@0:     past_ancest=mysetdiff(ps, act_slice);               % recover ancestors contained into past slices
wolffd@0:     app=ns(past_ancest);
wolffd@0:     rs_dim=[app(:)' rs_dim(:)'];                        %
wolffd@0: end                                                     %
wolffd@0: if length(rs_dim)==1, rs_dim=[1 rs_dim]; end            %
wolffd@0: if size(rs_dim,1)~=1, rs_dim=rs_dim';    end            %
wolffd@0: 
wolffd@0: fullm.T=reshape(fullm.T, rs_dim);                       % reshaping the marginal
wolffd@0: 
wolffd@0: % ----------------------------------------------------------------------------------------
wolffd@0: % ----------------------------------------------------------------------------------------
wolffd@0: 
wolffd@0: % X = cts parent, R = discrete self
wolffd@0: 
wolffd@0: % 1) observations vector -> CPD.parents_vals -------------------------------------------------
wolffd@0: x = cat(1, evidence{cdom});
wolffd@0: 
wolffd@0: % 2) weights vector -> CPD.eso_weights -------------------------------------------------------
wolffd@0: if isempty(evidence{self}) % R is hidden
wolffd@0:     sum_over=length(rs_dim);
wolffd@0:     app=sum(fullm.T, sum_over);    
wolffd@0:     pesi=reshape(app,[dpsz,1]);
wolffd@0:     clear app;
wolffd@0: else
wolffd@0:     pesi=reshape(fullm.T,[dpsz,1]);
wolffd@0: end
wolffd@0: 
wolffd@0: assert(approxeq(sum(pesi),1));
wolffd@0: 
wolffd@0: % 3) estimate (if R is hidden) or recover (if R is obs) self'value----------------------------
wolffd@0: if isempty(evidence{self})              % R is hidden    
wolffd@0:     app=mk_stochastic(fullm.T);         % P(self|Pa1,...,Pak)=fmarginal/prob(Pa1,...,Pak)
wolffd@0:     app=reshape(app,[dpsz ns(self)]);   % matrix size: prod{j,ns(Paj)} x ns(self)      
wolffd@0:     r=app;
wolffd@0:     clear app;
wolffd@0: else
wolffd@0:     r = zeros(dpsz,ns(self));
wolffd@0:     for i=1:dpsz
wolffd@0:         if pesi(i)~=0, r(i,evidence{self}) = 1; end
wolffd@0:     end
wolffd@0: end
wolffd@0: for i=1:dpsz
wolffd@0:     if pesi(i) ~=0, assert(approxeq(sum(r(i,:)),1)); end
wolffd@0: end
wolffd@0: 
wolffd@0: CPD.nsamples = CPD.nsamples + 1;            
wolffd@0: CPD.parent_vals(CPD.nsamples,:) = x(:)';
wolffd@0: for i=1:dpsz
wolffd@0:     CPD.eso_weights(CPD.nsamples,:,i)=pesi(i);
wolffd@0:     CPD.self_vals(CPD.nsamples,:,i) = r(i,:); 
wolffd@0: end