function [B,G] = beliefs(E,H,S,M,A,C,T,P)
% BELIEFS Calculate node and edge beliefs
%
%   [B,G] = BELIEFS(E,H,S,M) where E is an NxN adjacency matrix and M
%   is an NxN cell containing the Lx1 messages. B is an LxN
% matrix of resulting node beliefs, and G is the NxN cell
% of LxL edge beliefs.
%
%   [B,G] = 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,G] = BELIEFS(E,H,S,M,A,C) with C==1 will not exponentiate
%   beliefs, but return them in logspace.
%
%   [B,G] = 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,G] = BELIEFS(E,H,S,M,A,C,T,P) where P contains the message products

% Jerod Weinman
% jerod@acm.org
    
% $Id: beliefs.m 1 2009-05-07 13:41:00Z weinman $
 
% 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 = size(E,1); % Number of nodes
  L = size(H,1); % Number of labels

  if nargin<5
	A = 0;
  end;

  if nargin<6
	C = 0;
  end;
  
  if nargin<7
	T = [];
  end;

  if nargin<8
	% Product of all messages into each node
	msgProd = zeros(L,N);
	
	mp = zeros(L,1);
  
	for n=1:N
	  
	  if isempty(T)
		nbr = neighbors(E,n);
	  else
		nbr = T{n};
	  end;
	  
	  mp(:) = 0;
	  
	  for k=nbr(:)'
		mp = mp + M{k,n};
	  end;
	  
	  msgProd(:,n) = mp;
	  
	end;
  else
	msgProd = P;
  end;
  
	
  B = msgProd + H; % LxN

  mx = max(B,[],1); % 1xN
  
  if A
	% MAX PRODUCT
	[ll,nn] = find(mx(ones(1,L),:)==B);
	[nns,nni] = unique(nn);
	B = ll(nni)';
  else
	% SUM PRODUCT

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

	Z = mx + log ( sum( exp(B-mx(ones(1,L),:)), 1));
	
	B = B - Z(ones(1,L),:);
	
	if ~C
	  B = exp(B);
	  
	  % doubly-normalize
	  Z = sum(B);
	  B = B./Z(ones(1,L),:);
	end;
	
  end;
  
  if A || nargout==1
	return;
  end;
  
  [N0,N1] = find(E);
  

  G = cell(N,N);
  
  for e=1:length(N0)
	
	n0 = N0(e);
	n1 = N1(e);
	
	
	% Product of all messages into n0 (sans n1)
	prod0 = msgProd(:,n0) - M{n1,n0};   % 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{n0,n1};
	
	prod1 = prod1'; % 1xL
	
	% repmat( prod1',[L 1])
	prod1 = prod1(ones(1,L),:); % LxL -- n1 changes in col
	
	phi0 = H(:,n0*ones(1,L));  % LxL -- n0 changes in row
	phi1 = H(:,n1*ones(1,L))'; % LxL -- n1 changes in col
	
	prod = phi0 + phi1 + S{n0,n1} + 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
	
	G{n0,n1} = be;
	
  end;
