function B = beliefsf(G,M,A,C,P)
% BELIEFSF Calculate node and edge beliefs
%
%   B = BELIEFS(E,G,M) where E is an edge list, G is the model
% containing the set of node and edge energies (log compatibilities)
% M is an |E|xLx2 array containing the edge messages.
  
%   B = BELIEFS(E,H,S,M,A) where A==0 indicates M contains the
%   usual sum-product messages and marginals will be returned. A==1
%   indicates M contains max-product messages and a maximal
%   configuration will be returned in B.
%
%   B = BELIEFS(E,H,S,M,A,C) with C==1 will not exponentiate
%   beliefs, but return them in logspace.
%
%   B = BELIEFS(E,H,S,M,A,C,T) where T is a struct of length N
%   containing the list of neighbors for each node.
%
%   B = BELIEFS(E,H,S,M,A,C,T,P) where P contains the message products

% Jerod Weinman
% jerod@acm.org
    
% $Id%
 
% Jerod Weinman
% jerod@acm.org
% Copyright 2004, 2008, 2009
 
%    This file is part of Matlab GridCRF.
%
%    Matlab GridCRF is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    Matlab GridCRF is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with Matlab GridCRF.  If not, see <http://www.gnu.org/licenses/>.
  
  N = G.numNodes;
  L = G.numLabels;

  if nargin<2
    M = [];
  end;
  
  if nargin<3
	A = 0;
  end;

  if nargin<4
	C = 0;
  end;

  E0 = G.adjmat(:,1);
  E1 = G.adjmat(:,2);
  
  if nargin<5
    
	% Product of all messages into each node
	msgProd = zeros(L,N);

    if nargin>3 && ~isempty(M)
      msgProd = calcMsgProd(E0,E1,M,N);
    end;
  else
	msgProd = P;
  end;
  
  B = G;
	
  Bn = msgProd' + G.node; % NxL

  mx = max(Bn,[],2); % Nx1
  
  if A
	% MAX PRODUCT
    Bn = argmax(Bn,2);
  else
	% SUM PRODUCT

	% To normalize, find the sum in log space and subtract it.

	Z = mx + log ( sum( exp(Bn-mx(:,ones(1,L))), 2));
	
	Bn = Bn - Z(:,ones(1,L));
	
	if ~C
	  Bn = exp(Bn);
	  
	  % doubly-normalize
	  Z = sum(Bn,2);
	  Bn = Bn./Z(:,ones(1,L));
	end;
	
  end;
  
  if A
    B = Bn;
	return;
  end;
  
  B.node = Bn;
  
  B.edge = zeros(size(G.edge));
  
  if ~isempty(M)
    for e=1:length(E0)
      
      n0 = E0(e);
      n1 = E1(e);
	
      % Product of all messages into n0 (sans n1)
      prod0 = msgProd(:,n0) - M(:,e,2);   % Lx1
      
      % repmat( prod0,[1 L])
      prod0 = prod0(:,ones(1,L));    % LxL -- n0 changes in row
      
      % Product of all messages into n1 (sans n0)
      prod1 = msgProd(:,n1) - M(:,e,1);
      
      prod1 = prod1'; % 1xL
      
      % repmat( prod1',[L 1])
      prod1 = prod1(ones(1,L),:); % LxL -- n1 changes in col
      
      phi0 = G.node(n0*ones(1,L),:)';  % LxL -- n0 changes in row
      phi1 = G.node(n1*ones(1,L),:); % LxL -- n1 changes in col
  
      psi = G.edge(e,:);
      
      psi = reshape(psi,[L L]);
      
      prod = phi0 + phi1 + psi + prod0 + prod1;
      
      mx = max(max(prod));
      Z = mx + log(sumv(exp(prod-mx)));
      
      be = prod - Z;
      
      if ~C
        be = exp(be);
        
        % doubly normalize
        be = be./sumv(be);
      end
      
      B.edge(e,:) = be(:);
      
    end;
  else
    for e=1:length(E0)
      
      n0 = E0(e);
      n1 = E1(e);
	
      
      phi0 = G.node(:,n0*ones(1,L));  % LxL -- n0 changes in row
      phi1 = G.node(:,n1*ones(1,L))'; % LxL -- n1 changes in col
      
      prod = phi0 + phi1;
      
      mx = max(max(prod));
      Z = mx + log(sumv(exp(prod-mx)));
      
      be = prod - Z;
      
      if ~C
        be = exp(be);
        
        % doubly normalize
        be = be./sumv(be);
      end
      
      B.edge(e,:) = be(:);
      
    end;
  end;
