function [M,N] = spsumprodmessages(F,S)
% SPSUMPRODMESSAGES Calculate sum-product messages from factor sparsely 
% 
% M = SPSUMPRODMESSAGES(F,S) were F is the wordfactor and S contains
% the variable to factor messages. M returns the factor to variable messages.

% Jerod Weinman
% jerod@acm.org

% $Id: spsumprodmessages.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/>.

  doDisp = 0;
  
  % Arguments / Node names
  A = chararguments(F);
  
  % Marginalize each dimension

  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;
  
  %% Filter word list
  
  Ec = cell(1,wordlen(F));
  
  % Find eliminated characters
  for c=1:wordlen(F)
	Ec{c} = isinf(S{A(c)});
  end;
  
  % Find valid words
  Vw = true(numwords(F),1);
  
  if 1
    for c=1:wordlen(F)
      Vw = Vw & ~Ec{c}(F.words(:,c));
    end;
    
  else
    % VECTORIZE THIS!
    
    for w=1:numwords(F)
      for c=1:wordlen(F)
        if Ec{c}(F.words(w,c)) % if this word/slot's character is eliminated
          Vw(w) = 0;           % the word is not valid
          break;               % go to the next word
        end;
      end;
    end;
  end;
  
  %fprintf('Validated V(%d) = %d/%d\n',F.wordnode,sum(Vw),numwords(F));
  
  if doDisp
    defchars;
    
    for c=1:wordlen(F)
      fprintf('Ch %d: %s\n',c,sprintf('%s ', chars(find(~Ec{c}))));
    end;
    
    if sum(Vw)~=numwords(F)
      disp(chars(F.words(find(Vw),:)));
    end;
  end;
  
    
  N = sum(Vw);
  
  W = F.words(find(Vw),:);
  
  % Calculate word message products
  Mw = zeros(size(W));

  if 1
    Sc = cell2mat(S(A)');

    [ndxW,ndxC] = ndgrid(1:size(W,1),1:wordlen(F));
    ndxW = ndxW(:);
    ndxC = ndxC(:);
    
    indWC = sub2ind(size(Mw),ndxW,ndxC);
    if size(W,1)==1
      indS = sub2ind(size(Sc),W(indWC)',ndxC);
    else
      indS = sub2ind(size(Sc),W(indWC),ndxC);
    end;
    
    Mw(indWC) = Sc(indS);
  else
    
    % VECTORIZE THIS
    
    for w=1:size(W,1)
      for c=1:wordlen(F)
        
        Mw(w,c) = S{A(c)}(W(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(find(Vw),ones(1,wordlen(F)));

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

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

    for c=1:wordlen(F)
      M{c} = Mc(ones(1,F.numLabels),c) + ...
             (Mw(:,c)'*(W(:,c*ones(1,F.numLabels))==chList))';
    end;
    
  else

    for c=1:wordlen(F)
      M{c} = Mc(ones(1,F.numLabels),c); % non-word
    end;

    % VECTORIZE THIS
    for w=1:size(W,1)
      for c=1:wordlen(F)
        
        M{c}(W(w,c)) = M{c}(W(w,c)) + Mw(w,c);
        
      end;
    end;
  end;
  
  % Normalize messages and return to log space
  M2 = cellfun(@(x) log(x)-log(sum(x)), M,'UniformOutput',false);
  
  
  %% Word variable message
  
  if size(W,1)==0 % no valid words
	mm = [Mw0 ; -Inf];
  else
    % log-sum-exp over words
    % (log) multiply factor for word case
    mw = Mw1 + F.wordtable(Vw); 
    mx = max(mw);
    mm = [Mw0 ; mx+log(sum(exp(mw-mx)))];
  end;
  
  % Normalize message
  mx = max(mm);
  M2{end+1} = mm - mx - log(sum(exp(mm-mx)));
  


  for v=1:length(M)
	if any(isnan(M2{v}))
	  error('NaN in new factor->variable message');
	end;
  end;
  
  M = M2;
