Daniel@0: function CPD = update_ess(CPD, fmarginal, evidence, ns, cnodes, hidden_bitv)
Daniel@0: % UPDATE_ESS Update the Expected Sufficient Statistics of a softmax node
Daniel@0: % function CPD = update_ess(CPD, fmarginal, evidence, ns, cnodes, hidden_bitv)
Daniel@0: %
Daniel@0: % fmarginal = overall posterior distribution of self and its parents
Daniel@0: % fmarginal(i1,i2...,ik,s)=prob(Pa1=i1,...,Pak=ik, self=s| X)
Daniel@0: % 
Daniel@0: % => 1) prob(self|Pa1,...,Pak)=fmarginal/prob(Pa1,...,Pak) with prob(Pa1,...,Pak)=sum{s,fmarginal}
Daniel@0: %       [self estimation -> CPD.self_vals]
Daniel@0: % 	  2) prob(Pa1,...,Pak) [WIRLS weights -> CPD.eso_weights]
Daniel@0: %
Daniel@0: % Hidden_bitv is ignored
Daniel@0: 
Daniel@0: % Written by Pierpaolo Brutti
Daniel@0: 
Daniel@0: if ~adjustable_CPD(CPD), return; end
Daniel@0: 
Daniel@0: domain     = fmarginal.domain;                              
Daniel@0: self       = domain(end);          
Daniel@0: ps         = domain(1:end-1);                                     
Daniel@0: cnodes     = domain(CPD.cpndx);
Daniel@0: cps        = myintersect(domain, cnodes);                     
Daniel@0: dps        = mysetdiff(ps, cps);                            
Daniel@0: dn_use     = dps;
Daniel@0: if isempty(evidence{self}) dn_use = [dn_use self]; end % if self is hidden we must consider its dimension  
Daniel@0: dps_as_cps = domain(CPD.dps_as_cps.ndx);
Daniel@0: odom       = domain(~isemptycell(evidence(domain))); 
Daniel@0: 
Daniel@0: ns = zeros(1, max(domain));
Daniel@0: ns(domain) = CPD.sizes;     % CPD.sizes = bnet.node_sizes([ps self]);
Daniel@0: ens = ns;                   % effective node sizes
Daniel@0: ens(odom) = 1;              
Daniel@0: dpsize = prod(ns(dps));
Daniel@0: 
Daniel@0: % Extract the params compatible with the observations (if any) on the discrete parents (if any)
Daniel@0: dops = myintersect(dps, odom);
Daniel@0: dpvals = cat(1, evidence{dops});
Daniel@0: 
Daniel@0: subs = ind2subv(ens(dn_use), 1:prod(ens(dn_use)));
Daniel@0: dpmap = find_equiv_posns(dops, dn_use);
Daniel@0: if ~isempty(dpmap), subs(:,dpmap) = subs(:,dpmap)+repmat(dpvals(:)',[size(subs,1) 1])-1; end
Daniel@0: supportedQs = subv2ind(ns(dn_use), subs); subs=subs(1:prod(ens(dps)),1:length(dps));
Daniel@0: Qarity = prod(ns(dn_use));
Daniel@0: if isempty(dn_use), Qarity = 1; end   
Daniel@0: 
Daniel@0: fullm.T              = zeros(Qarity, 1);
Daniel@0: fullm.T(supportedQs) = fmarginal.T(:);
Daniel@0: rs_dim = CPD.sizes;    rs_dim(CPD.cpndx) = 1;           %
Daniel@0: if ~isempty(evidence{self}), rs_dim(end)=1; end         % reshaping the marginal
Daniel@0: fullm.T              = reshape(fullm.T, rs_dim);        %
Daniel@0: 
Daniel@0: % --------------------------------------------------------------------------------UPDATE--
Daniel@0: 
Daniel@0: CPD.nsamples = CPD.nsamples + 1;
Daniel@0: 
Daniel@0: % 1) observations vector -> CPD.parents_vals ---------------------------------------------
Daniel@0: cpvals = cat(1, evidence{cps});
Daniel@0: 
Daniel@0: if ~isempty(dps_as_cps),   % ...get in the dp_as_cp parents... 
Daniel@0:     separator          = CPD.dps_as_cps.separator;
Daniel@0:     dp_as_cpmap        = find_equiv_posns(dps_as_cps, dps);       
Daniel@0:     for i=1:dpsize,
Daniel@0:         dp_as_cpvals=zeros(1,sum(ns(dps_as_cps)));
Daniel@0:         possible_vals = ind2subv(ns(dps),i);
Daniel@0:         ll=find(ismember(subs(:,dp_as_cpmap), possible_vals(dp_as_cpmap), 'rows')==1);   
Daniel@0:         if ~isempty(ll),
Daniel@0:             where_one = separator + possible_vals(dp_as_cpmap);
Daniel@0:             dp_as_cpvals(where_one)=1;                            
Daniel@0:         end
Daniel@0:         CPD.parent_vals(CPD.nsamples,:,i) = [dp_as_cpvals(:); cpvals(:)]';
Daniel@0:     end
Daniel@0: else
Daniel@0:     CPD.parent_vals(CPD.nsamples,:) = cpvals(:)';
Daniel@0: end
Daniel@0: 
Daniel@0: % 2) weights vector -> CPD.eso_weights ----------------------------------------------------
Daniel@0: if isempty(evidence{self}),             % self is hidden
Daniel@0:     pesi=reshape(sum(fullm.T, length(rs_dim)),[dpsize,1]);
Daniel@0: else
Daniel@0:     pesi=reshape(fullm.T,[dpsize,1]);
Daniel@0: end
Daniel@0: assert(approxeq(sum(pesi),1));          % check
Daniel@0: 
Daniel@0: % 3) estimate (if R is hidden) or recover (if R is obs) self'value-------------------------
Daniel@0: if isempty(evidence{self})                                  % P(self|Pa1,...,Pak)=fmarginal/prob(Pa1,...,Pak)
Daniel@0:     r=reshape(mk_stochastic(fullm.T), [dpsize ns(self)]);   % matrix size: prod{j,ns(Paj)} x ns(self)      
Daniel@0: else
Daniel@0:     r = zeros(dpsize,ns(self));
Daniel@0:     for i=1:dpsize, if pesi(i)~=0, r(i,evidence{self}) = 1; end; end
Daniel@0: end
Daniel@0: for i=1:dpsize, if pesi(i)~=0, assert(approxeq(sum(r(i,:)),1)); end; end     % check
Daniel@0: 
Daniel@0: % 4) save the previous values --------------------------------------------------------------
Daniel@0: for i=1:dpsize
Daniel@0:     CPD.eso_weights(CPD.nsamples,:,i)=pesi(i);
Daniel@0:     CPD.self_vals(CPD.nsamples,:,i) = r(i,:); 
Daniel@0: end