www.pudn.com > ConstrainedEM.zip > mpe_local_net.m

function [mpe, ll] = mpe_local_net(engine,pots)
% origin is calc_mpe_bucket of BNT

% engine is can object containing the mnet and the elimination order.
% pots is the evidence to put in the nodes.


max_over=engine.elim_order;
mnet=engine.mnet;
n = length(mnet);
k=size(pots,2);

% enter the evidence
B = cell(1,n); 
for i=1:n
  B{i} = dpot( i, k, pots(i,:) );
end

order = max_over; % no attempt to optimize this


% enter the constraints potentials

[x y]=find(triu(engine.mnet)==1);
for i=1:length(y)
  b = bucket_num([ x(i) y(i) ] , order);
  tmpot=dpot([ x(i) y(i) ],[ k k ],ones(k,k)-eye(k));
  B{b} = multiply_pots(B{b}, tmpot );
end


% Do backward phase
max_over = max_over(length(max_over):-1:1); % reverse
for i=max_over(1:end-1)        
  % max-ing over variable i which occurs in bucket j
  j = bucket_num(i, order);
  rest = mysetdiff(domain_pot(B{j}), i);
  %temp = marginalize_pot_max(B{j}, rest);
  temp = marginalize_pot(B{j}, rest, 1);
  b = bucket_num(domain_pot(temp), order);
          sprintf('maxing over bucket %d (var %d), putting result into bucket %d\n', j, i, b)
  sB=struct(B{b});  % violate object privacy
  if ~isempty(sB.domain)
    B{b} = multiply_pots(B{b}, temp);
  else
    B{b} = temp;
  end
end
result = B{1};
marginal = pot_to_marginal(result);
[prob, mpe] = max(marginal.T);

% handle impossible cases
if ~(prob>0)
  mpe = [];    
  ll = -inf;
  %warning('evidence has zero probability')
  return
end

ll = log(prob);

% Do forward phase    
for ii=2:n
  marginal = pot_to_marginal(B{ii});
  mpeidx = [];
  for jj=order(1:length(mpe))
    assert(ismember(jj, marginal.domain)) %%% bug
    temp = find_equiv_posns(jj, marginal.domain);
    mpeidx = [mpeidx, temp] ;
    if isempty(temp)
      mpeidx = [mpeidx, Inf] ;
    end
  end
  [mpeidxsorted sortedtompe] = sort(mpeidx) ;
  
  % maximize the matrix obtained from assigning values from previous buckets.
  % this is done by building a string and using eval.
  
  kk=1;
  sargs = '(';
  for jj=1:length(marginal.domain)
    if (jj~=1)
      sargs = [sargs, ','];
    end
    if (mpeidxsorted(kk)==jj)
      sargs = [sargs, num2str(mpe(sortedtompe(kk)))];
      if (kk