wolffd@0: function [mpe, ll] = calc_mpe_bucket(bnet, new_evidence, max_over)
wolffd@0: %
wolffd@0: % PURPOSE:
wolffd@0: %       CALC_MPE Computes the most probable explanation to the network nodes
wolffd@0: %       given the evidence.
wolffd@0: %       
wolffd@0: %       [mpe, ll] = calc_mpe(engine, new_evidence, max_over)
wolffd@0: %
wolffd@0: % INPUT:
wolffd@0: %       bnet  - the bayesian network
wolffd@0: %       new_evidence - optional, if specified - evidence to be incorporated [cell(1,n)]
wolffd@0: %       max_over - optional, if specified determines the variable elimination order [1:n]
wolffd@0: %
wolffd@0: % OUTPUT:
wolffd@0: %       mpe - the MPE assignmet for the net variables (or [] if no satisfying assignment)
wolffd@0: %       ll - log assignment probability.
wolffd@0: %
wolffd@0: % Notes:
wolffd@0: % 1. Adapted from '@var_elim_inf_engine\marginal_nodes' for MPE by Ron Zohar, 8/7/01
wolffd@0: % 2. Only discrete potentials are supported at this time.
wolffd@0: % 3. Complexity: O(nw*) where n is the number of nodes and w* is the induced tree width.
wolffd@0: % 4. Implementation based on:
wolffd@0: %  - R. Dechter, "Bucket Elimination: A Unifying Framework for Probabilistic Inference", 
wolffd@0: %                 UA1 96, pp. 211-219.
wolffd@0: 
wolffd@0: 
wolffd@0: ns = bnet.node_sizes;
wolffd@0: n = length(bnet.dag);
wolffd@0: evidence = cell(1,n);
wolffd@0: if (nargin<2)
wolffd@0:     new_evidence = evidence;
wolffd@0: end
wolffd@0: 
wolffd@0: onodes = find(~isemptycell(new_evidence));  % observed nodes
wolffd@0: hnodes = find(isemptycell(new_evidence));  % hidden nodes
wolffd@0: pot_type = determine_pot_type(bnet, onodes);
wolffd@0: 
wolffd@0: if pot_type ~= 'd'
wolffd@0:   error('only disrete potentials supported at this time')    
wolffd@0: end
wolffd@0: 
wolffd@0: for i=1:n
wolffd@0:   fam = family(bnet.dag, i);
wolffd@0:   CPT{i} = convert_to_pot(bnet.CPD{bnet.equiv_class(i)}, pot_type, fam(:), evidence);        
wolffd@0: end 
wolffd@0: 
wolffd@0: % handle observed nodes: set impossible cases' probability to zero
wolffd@0: % rather than prun matrix (this makes backtracking easier)
wolffd@0: 
wolffd@0: for ii=onodes
wolffd@0:   lIdx = 1:ns(ii);
wolffd@0:   lIdx = setdiff(lIdx, new_evidence{ii});
wolffd@0:   
wolffd@0:   sCPT=struct(CPT{ii});  % violate object privacy
wolffd@0:   
wolffd@0:   sargs = '';
wolffd@0:   for jj=1:(length(sCPT.domain)-1)
wolffd@0:     sargs = [sargs, ':,']; 
wolffd@0:   end        
wolffd@0:   for jj=lIdx
wolffd@0:     eval(['sCPT.T(', sargs, num2str(jj), ')=0;']);
wolffd@0:   end
wolffd@0:   CPT{ii}=dpot(sCPT.domain, sCPT.sizes, sCPT.T);        
wolffd@0: end
wolffd@0: 
wolffd@0: B = cell(1,n); 
wolffd@0: for b=1:n
wolffd@0:   B{b} = mk_initial_pot(pot_type, [], [], [], []);
wolffd@0: end
wolffd@0: 
wolffd@0: if (nargin<3)
wolffd@0:   max_over = (1:n);
wolffd@0: end   
wolffd@0: order = max_over; % no attempt to optimize this
wolffd@0: 
wolffd@0: 
wolffd@0: % Initialize the buckets with the CPDs assigned to them
wolffd@0: for i=1:n
wolffd@0:   b = bucket_num(domain_pot(CPT{i}), order);
wolffd@0:   B{b} = multiply_pots(B{b}, CPT{i});
wolffd@0: end
wolffd@0: 
wolffd@0: % Do backward phase
wolffd@0: max_over = max_over(length(max_over):-1:1); % reverse
wolffd@0: for i=max_over(1:end-1)        
wolffd@0:   % max-ing over variable i which occurs in bucket j
wolffd@0:   j = bucket_num(i, order);
wolffd@0:   rest = mysetdiff(domain_pot(B{j}), i);
wolffd@0:   %temp = marginalize_pot_max(B{j}, rest);
wolffd@0:   temp = marginalize_pot(B{j}, rest, 1);
wolffd@0:   b = bucket_num(domain_pot(temp), order);
wolffd@0:   %        fprintf('maxing over bucket %d (var %d), putting result into bucket %d\n', j, i, b);
wolffd@0:   sB=struct(B{b});  % violate object privacy
wolffd@0:   if ~isempty(sB.domain)
wolffd@0:     B{b} = multiply_pots(B{b}, temp);
wolffd@0:   else
wolffd@0:     B{b} = temp;
wolffd@0:   end
wolffd@0: end
wolffd@0: result = B{1};
wolffd@0: marginal = pot_to_marginal(result);
wolffd@0: [prob, mpe] = max(marginal.T);
wolffd@0: 
wolffd@0: % handle impossible cases
wolffd@0: if ~(prob>0)
wolffd@0:   mpe = [];    
wolffd@0:   ll = -inf;
wolffd@0:   %warning('evidence has zero probability')
wolffd@0:   return
wolffd@0: end
wolffd@0: 
wolffd@0: ll = log(prob);
wolffd@0: 
wolffd@0: % Do forward phase    
wolffd@0: for ii=2:n
wolffd@0:   marginal = pot_to_marginal(B{ii});
wolffd@0:   mpeidx = [];
wolffd@0:   for jj=order(1:length(mpe))
wolffd@0:     assert(ismember(jj, marginal.domain)) %%% bug
wolffd@0:     temp = find_equiv_posns(jj, marginal.domain);
wolffd@0:     mpeidx = [mpeidx, temp] ;
wolffd@0:     if isempty(temp)
wolffd@0:       mpeidx = [mpeidx, Inf] ;
wolffd@0:     end
wolffd@0:   end
wolffd@0:   [mpeidxsorted sortedtompe] = sort(mpeidx) ;
wolffd@0:   
wolffd@0:   % maximize the matrix obtained from assigning values from previous buckets.
wolffd@0:   % this is done by building a string and using eval.
wolffd@0:   
wolffd@0:   kk=1;
wolffd@0:   sargs = '(';
wolffd@0:   for jj=1:length(marginal.domain)
wolffd@0:     if (jj~=1)
wolffd@0:       sargs = [sargs, ','];
wolffd@0:     end
wolffd@0:     if (mpeidxsorted(kk)==jj)
wolffd@0:       sargs = [sargs, num2str(mpe(sortedtompe(kk)))];
wolffd@0:       if (kk<length(mpe))
wolffd@0: 	kk = kk+1 ;
wolffd@0:       end
wolffd@0:     else
wolffd@0:       sargs = [sargs, ':'];
wolffd@0:     end
wolffd@0:   end
wolffd@0:   sargs = [sargs, ')'] ;   
wolffd@0:   eval(['[val, loc] = max(marginal.T', sargs, ');'])        
wolffd@0:   mpe = [mpe loc];
wolffd@0: end     
wolffd@0: [I,J] = sort(order);
wolffd@0: mpe = mpe(J);
wolffd@0: 
wolffd@0: 
wolffd@0: 
wolffd@0: %%%%%%%%%
wolffd@0: 
wolffd@0: function b = bucket_num(domain, order)
wolffd@0: 
wolffd@0: b = max(find_equiv_posns(domain, order));
wolffd@0: