function [B,F,M] = syncbp(G,A,M,C)
%    SYNCBP Loopy BP using  a synchronous schedule
%
% [B,F,M] = SYNCBP(G) where G is a list/flat graphical model struct.
%
% [B,F,M] = SYNCBP(G,A) where A==0 (default) indicates use of
% the sum-product algorithm for marginals. A==1 indicates use of
% max-product to find a maximizing configuration, which will be
% returned in B.
%
% [B,F,M] = SYNCBP(G,A,M) where M are the set of messages to
% begin with.
%  
% [B,F,M] = SYNCBP(G,A,M,C) where C==1 will return beliefs
% in logspace.
    
% Jerod Weinman
% jerod@acm.org

%	$Id: syncbp.m 1 2009-05-07 13:41:00Z weinman $
 
% Jerod Weinman
% jerod@acm.org
% Copyright 2004, 2008, 2009
 
%    This file is part of Matlab GridCRF.
%
%    Matlab GridCRF 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.
%
%    Matlab GridCRF 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 Matlab GridCRF.  If not, see <http://www.gnu.org/licenses/>.
  
% Damping coefficient for message updates
%   alpha of zero would mean ignore updates and 
%   alpha of one would mean use no damping
    
  alpha = 0.75;
  %alpha = 1;
  
  N = G.numNodes;
  L = G.numLabels;
  
  if nargin<2
	A = 0;
  end;
  
  if nargin<4
	C = 0;
  end;
  
%  L = size(G.node,2);

  E0 = G.adjmat(:,1);
  E1 = G.adjmat(:,2);
  
    
  nEdg = length(E0);

  % Product of all messages into node n
  msgProd = zeros(N,L);

  if nargin<3 || isempty(M)
	% Initialize Messages
    M = zeros(nEdg,L,2); 

  else
	
    for n=1:N
      msgProd(n,:) = sum(M(find(E0==n),:,2),1) + sum(M(find(E1==n),:,1),1);
    end;
	
  end;

  MAX_ITER = min(5*N,5000);
  MSG_TOL = 1e-3;

  % E0->E1: 1
  % E1->E0: 2
  Mnew = zeros(size(E0,1),L,2); 
  
  
  Un = G.node;
  Ue0 = reshape(G.edge,[size(G.edge,1) L L]); % (E0,E1)
  Ue1 = permute(Ue0,[1 3 2]); % (E1,E0)
  
  for iter=1:MAX_ITER
    
    %tic
    
	converging = 1;
	l2diff = 0;

    % E0 -> E1
    
    % Message Product * Local Term
    mp = Un(E0,:) + msgProd(E0,:) - M(:,:,2);

    % Message "Summands"
    ms = Ue0 + mp(:,:,ones(1,L));

    % Calculate max
    mx = max(ms,[],2);

    if A==0 % SUM-PRODUCT      
      Mnew(:,:,1) = mx + log(sum( exp(ms - mx(:,ones(1,L),:)), 2));
    else    % MAX-PRODUCT
      Mnew(:,:,1) = mx;
    end;
    
    %%
    %% E1 -> E0
    %%

    % Message Product * Local Term
    mp = Un(E1,:) + msgProd(E1,:) - M(:,:,1);

    % Message "Summands"
    ms = Ue1 + mp(:,:,ones(1,L));

    % Calculate max
    mx = max(ms,[],2);
    
    if A==0 % SUM-PRODUCT      
      Mnew(:,:,2) = mx + log(sum( exp(ms - mx(:,ones(1,L),:)), 2));
    else    % MAX-PRODUCT
      Mnew(:,:,2) = mx;
    end;
    

    % Normalize messages
    mx = max(Mnew,[],2);    
    Z = mx + log(sum(exp(Mnew-mx(:,ones(1,L),:)),2));
    Mnew = Mnew - Z(:,ones(1,L),:);

    % Convergence checks/progress
    l2diff = l2diff + sumv( (M-Mnew).^2 );
    
    if converging && anyv(abs(M - Mnew)>MSG_TOL)
      converging = 0;
    end;

    % Assign new messages
    M = M*(1-alpha) + alpha*Mnew;
    

    msgProd = zeros(N,L);

    % Update message products

    if 1 % MEXed!
      msgProd = calcMsgProd(E0,E1,M,N);
    else % Non-vectorized:
      for n=1:N
        msgProd(n,:) = sum(M(find(E0==n),:,2),1) + sum(M(find(E1==n),:,1),1);
      end;
	end;
	%fprintf('%d\t%f\n',iter,l2diff);
	
	if converging
	  %disp(['syncbp converged in ' num2str(iter) ' iterations']);
	  F = 1;
	  break;
    else
      %fprintf('L2 = %f\n',l2diff)
	end;

    %toc
  end;
  
  %keyboard
  
  if iter==MAX_ITER
	disp(sprintf('syncbp(%d) failed to converge',A));
	F = -1;
  end;	


  if A==0
    M = permute(M,[2 1 3]);
    B = beliefsf(G,M,A,C,msgProd');
  else
    B = argmax(msgProd,2);
  end;
