function M = messages(F,S,R,T)
%   MESSAGES Calculates messages from the factor to variables
%
% M = MESSAGES(F,S) where F is the factor, S are the variable
% to factor messages for the graph the factor belongs to.
%
% M = MESSAGES(F,S,R) where R=1 indicates that the max-product (rather
% than sum-product) algorithm should be used.
%
% M = MESSAGES(F,S,R,T) where T contain the pre-calculated
% Cartesian outer product indices for multiplying messages and the
% factor. 
 
% Jerod Weinman
% jerod@acm.org
    
% $Id: messages.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<3
	R = 0;
  end;
  
  % Arguments / Node names
  A = F.arguments;
  
  % 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;
  
  % If a two-variable factor, handle specially (faster? should
  % be. but not)
  if 0%length(A)==2
	
	for a=1:2
	  if a==1
		P = F.table + S{A(1)}(:,ones(1,L(2)));
	  else
		P = (F.table + S{A(2)}(:,ones(1,L(1)))')'; 
	  end;
	  
	  mx = max(P,[],1);
	
	  if R
		% MAX PRODUCT
		mm = mx;
	  else
		% intermediate table (max factor removed, exponentiated)
		mt = exp(P-mx(ones(1,L(a)),:));
		mm = mx + log(sum(mt,1)); % take marginal
	  mm = reshape(mm,[L(a) 1]); % vectorize appropriately
	end;
	
	mx = max(mm); % normalize
	M{a} = mm - mx - log(sum(exp(mm-mx)));
	
	end;
	return;
  end;
	
  if nargin<4
	T = combindices(F);
  end;
  
  % Take (log) product of combinatorial combinations
  P = F.table;
  
  for a=1:length(A)
    Sa{a} = S{A(a)}(T{a});
	P = P + Sa{a};
  end;
  
  % Marginalize each dimension

  for a=1:length(A)
	
	% Product of factor with all messages except for the current variable
	%mm = P - S{A(a)}(T{a});
    mm = P - Sa{a};
    
	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{a} = mm - mx - log(sum(exp(mm-mx)));
	
  end;
	
  
