function M = spmessages(F,S,R,T)
%   SPMESSAGES Calculates messages sparsely from the factor to variables
%
% M = SPMESSAGES(F,S) where F is the factor, S are the variable
% to factor messages for the graph the factor belongs to.
%
% M = SPMESSAGES(F,S,R) where R=1 indicates that the max-product (rather
% than sum-product) algorithm should be used.
%
% M = SPMESSAGES(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: spmessages.m 33 2011-07-26 20:28:07Z 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;

  % Find finite indices
  I = cellfun(@(ss) find(~isinf(ss)),S(A),'UniformOutput',0)';

  % Calculate message to each variable

  for a=1:length(A)

    % Include all indices of the message for this variable
    I2 = I;
    I2{a} = 1:L(a);
	  
    % Indices for all other valid combindations
    T = combindices(F,I2);
	  
	  
    T2 = cellfun(@mkvector,T,'UniformOutput',0); % shape into vectors
		
    ndx = sub2ind(size(F.table),T2{:}); % and get linear indices
    
    % Take (log) product of combinatorial combinations
    P = reshape(F.table(ndx),cellfun('length',I2));

    % Product of factor with all messages except for the current variable
    for b=1:length(A)

      if b==a
        continue;
      end;
      
      % Matlab is stupid. This reshape is unecessary but for it's
      % vector indicing's tendency to give you back something the
      % same orientation (row/col) as the vector, rather than the
      % size of the index, which is the behavior when the index is
      % anything but a vector.
      
      P = P + reshape(S{A(b)}(T{b}),size(T{b}));
      
    end;
    

    
    % Marginalize other dimensions
    for b=1:length(A)
      
      if b==a
        continue;
      end;
      
      % max along the dimension
      mx = max(P,[],b);
      
      if R
        % MAX PRODUCT
        P = 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) = length(I{b});
        
        % intermediate table
        mt = exp(P - repmat(mx,rd));
        
        % step 3: marginalize dimension,
        P = mx + log(sum(mt,b));
      end;
    end;
    
    P = reshape(P,L(a),1);
    
    % Normalize message
    mx = max(P);
    
    M{a} = P - mx - log(sum(exp(P-mx)));
    
  end;
	
  
