function M = sumprodmessages(F,S)
% SUMPRODMESSAGES Calculate sum-product messages from a factor
% 
% M = SUMPRODMESSAGES(F,S) were F is the wordfactor and S contains
% the variable to factor messages. M is the factor to variable messages.

  
% Jerod Weinman
% jerod@acm.org
  
% $Id: sumprodmessages.m 50 2011-07-29 16:01:00Z weinman $
 
% Jerod Weinman
% jerod@acm.org
% Copyright 2006, 2008, 2009, 2011
 
%    This file is part of VIDI.
%
%    VIDI 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.
%
%    VIDI 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 VIDI.  If not, see <http://www.gnu.org/licenses/>.
  
% Arguments / Node names
  A = chararguments(F);
  
  
  % Marginalize each dimension - sum over each character message

  Mc = zeros(1,length(A));

  for a=1:length(A)

    % We need to log-sum-exp character node to factor messages
    
    mx = max(S{A(a)});

    Mc(a) = mx + log(sum(exp(S{A(a)}-mx)));
    
  end;

  % (log) product of all message sums 
  Mw0 = sum(Mc);

  % product (of all message sums) except the node at a given index.
  Mc = Mw0-Mc;
  
  % Add word node to factor message for non-word case
  Mc = S{F.wordnode}(1) + Mc;
  
  % Calculate word message products
  Mw = zeros(size(F.words));
  
  if 1
    Sc = cell2mat(S(A)');

    [ndxW,ndxC] = ndgrid(1:numwords(F),1:wordlen(F));
    ndxW = ndxW(:);
    ndxC = ndxC(:);
    
    indWC = sub2ind(size(Mw),ndxW,ndxC);

    if numwords(F)==1
      indS = sub2ind(size(Sc),F.words(indWC)',ndxC);
    else
      indS = sub2ind(size(Sc),F.words(indWC),ndxC);
    end;
    
    Mw(indWC) = Sc(indS);
  else
    % VECTORIZE THIS
    for w=1:numwords(F)
      for c=1:wordlen(F)
        
        Mw(w,c) = S{A(c)}(F.words(w,c));
        
      end;
    end;
  end;
  
  % (log) product of all messages (for each word/row) 
  Mw1 = sum(Mw,2);
  
  % product except the node at a given (column) index
  Mw = Mw1(:,ones(1,wordlen(F))) - Mw;
  
  % Add word node to factor message for word case
  Mw = S{F.wordnode}(2) + Mw;

  % Add factor for word case
  Mw = Mw + F.wordtable(:,ones(1,wordlen(F)));

  %% Character messages
  
  mx = max(maxv(Mc),maxv(Mw));
  %mx = 0;
  
  Mc = exp(Mc-mx);
  Mw = exp(Mw-mx);
  
  if any(isinf(Mc(:))) || any(isinf(Mw(:)))
    error('Infinite message after exp');
  end;
	
  if any(Mc(:)==0) || any(Mw(:)==0)
    warning('Zeros after exp');
  end;

  % Character messages: non-word + word parts
  if 1
    
    chList = 1:F.numLabels;
    chList = chList(ones(1,numwords(F)),:);

    for c=1:wordlen(F)
      M{c} = Mc(ones(1,F.numLabels),c) + ... % non-word
             (Mw(:,c)'*(F.words(:,c*ones(1,F.numLabels))==chList))'; % word
    end;
    
  else
    % VECTORIZE THIS
    
    for c=1:wordlen(F)
      M{c} = Mc(ones(1,F.numLabels),c); % non-word
    end;
    
    for w=1:numwords(F)
      for c=1:wordlen(F)
        
        M{c}(F.words(w,c)) = M{c}(F.words(w,c)) + Mw(w,c);
        
      end;
    end;
  end;
  
  % Normalize messages and return to log space
  M = cellfun(@(x) log(x)-log(sum(x)), M,'UniformOutput',false);
  
  
  %% Word variable message
  
  % log-sum-exp over words
  % (log) multiply factor for word case
  mw = Mw1 + F.wordtable; 
  mx = max(mw);
  mm = [Mw0 ; mx+log(sum(exp(mw-mx)))];
  
  % Normalize message
  mx = max(mm);
  M{end+1} = mm - mx - log(sum(exp(mm-mx)));
  


  
