function M = spmaxprodmessages(F,S)
% SPMAXPRODMESSAGES Calculate max-product messages from factor sparsely 
% 
% M = SPMAXPRODMESSAGES(F,S) were F is the wordfactor and S contains
% the variable to factor messages

    
% Jerod Weinman
% jerod@acm.org
    
% $Id: spmaxprodmessages.m 1 2009-06-17 20:04:43Z 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);
  
  % Maximize each dimension

  Mc = zeros(1,length(A));
  
  for a=1:length(A)

	Mc(a) = max(S{A(a)});
	
  end;
  
  % (log) product of all message max
  Mw0 = sum(Mc);

  % product 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));
  
  Vwx = find(Vw);
  W = F.words(Vwx,:);

  % Calculate word message products
  Mw = zeros(size(F.words));

  if 1
    Sc = cell2mat(S(A)');
    
    [ndxW,ndxC] = ndgrid(Vwx,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: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;
  
  % (log) multiply factors for word case
  Mw = Mw + F.wordtable(:,ones(1,wordlen(F)));
  Mw1 = Mw1 + F.wordtable;
  
  % Add word node to factor message for word case
  Mw = S{F.wordnode}(2) + Mw;


  %% 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;

  
  if 1

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

  else
    % Initialize character messages with non-word part
    for a=1:length(A)
      M{a} = Mc(ones(1,F.numLabels),a);
    end;
    
    % Compare word part to character messages
    for w=1:size(W,1)
      for c=1:wordlen(F)
        
        M{c}(F.words(w,c)) = max( M{c}(W(w,c)), Mw(w,c));
        
      end;
    end;
  end;
  
  %% Word variable message

  if size(W,1) == 0 % no valid words
	M{end+1} = [Mw0 ; -Inf];
  else
	M{end+1} = [Mw0 ; max(Mw1)];
  end;
  
  % Normalize messages
  for m=1:length(M)
	mm= M{m};
	mx = max(mm);
	M{m} = mm - mx - log(sum(exp(mm-mx)));
  end;
  

