function M = message(F,P,S,N,R)
% MESSAGES Calculates message from the factor to a variable
%
%  M = MESSAGES(F,P,S,N,R) where F is the factor, P is ?

    
% Jerod Weinman
% jerod@acm.org
    
% $Id: message.m 1 2009-05-07 13:41:00Z weinman $
 
% Jerod Weinman
% jerod@acm.org
% Copyright 2008, 2009, 2011
 
%    This file is part of MatLibFactorGraph.
%
%    MatLibFactorGraph 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.
%
%    MatLibFactorGraph 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 MatLibFactorGraph.  If not, see <http://www.gnu.org/licenses/>.

  if nargin<4
	R = 0;
  end;
  
  % Arguments / Node names
  A = F.arguments;
  
  if ~any(A==N)
	error('Node %d is not a factor argument',N);
  end;
  
	
  % Label dimensions for argument nodes
  L = size(F.table);
  
  % If a singleton factor, return normalized vector
  if length(A)==1
	mm = F.table;
	mx = max(mm);
	M = mm - mx - log(sum(exp(mm-mx)));
	return;
  end;
  
  % Linear indices spanning each variable of the factor
  for a=1:length(A)
	I{a} = 1:L(a);
  end;
  
  % Combinatorial subscripts for Cartesian product
  T = cellndgrid(I);

	
  P = P - S{N}(T{A==N});
  
  
  % Marginalize each dimension
  for b=1:length(A)
	  
	if b==a
	  continue;
	end;
	  
	% max along the dimension
	mx = max(mm,[],b);
	
	if R
	  % MAX PRODUCT
	  mm = mx;
	else
	  % SUM PRODUCT
	  
	  % Marginalize dimension b.
	  
	  % step 1: find maximum along dimension. (done above)
	  
	  % step 2: for accuracy, remove largest factor before exponentiating
	  
	  % replication dimension: we have to replicate the maximum
	  % along the dimension over which we took it
	  rd = ones(1,length(A));
	  rd(b) = L(b);
	  
	  % intermediate table
	  mt = exp(mm - repmat(mx,rd));
	  
	  % step 3: marginalize dimension,
	  mm = mx + log(sum(mt,b));
	end;
  end;
  
  mm = reshape(mm,L(a),1);
  
  % Normalize message
  mx = max(mm);
  
  M = mm - mx - log(sum(exp(mm-mx)));
  

	
  
