function [B,G,F,M] = asyncbp(E,H,S,A,M,C,T)
%    ASYNCBP Loopy BP using a random asynchronous schedule
%
% [B,G,F,M] = ASYNCBP(E,H,S) where E is an adjacency matrix, H is the
% set of node energies (log compatibilities) and S is the set of
% edge energies.
%
% [B,G,F,M] = ASYNCBP(E,H,S,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,G,F,M] = ASYNCBP(E,H,S,A,M) where M are the set of messages to
% begin with.
%  
% [B,G,F,M] = ASYNCBP(E,H,S,A,M,C) where C==1 will return beliefs
% in logspace.
%
% [B,G,F,M] = ASYNCBP(E,H,S,A,M,C,T) where T is a struct of length N
%   containing the list of neighbors for each node.
  
%	$Id: asyncbp.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/>.
  
  N = size(E,1);
  
  if size(E,2)~=N
	error('Adjacency matrix must be square');
  end;
  
  if (size(H,2)~=N)
	error('Incorrect dimension of node compatibilities');
  end;

  if nargin<4
	A = 0;
  end;
  
  if nargin<6
	C = 0;
  end;
  
  if nargin<7
	T= {};
  end;
  
  L = size(H,1);

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

  if nargin<5 || isempty(M)
	% Initialize Messages
	M = cell(N,N);
	
	% Find all edges
	[a,b] = find(E);
	
	% Initialize (log) messages: ln 1 = 0
	for n=1:length(a)
	  M{a(n),b(n)} = zeros(L,1);
	  M{b(n),a(n)} = zeros(L,1);
	end;
  else
	
	mp = zeros(L,1);
	
	for n=1:N
	  if isempty(T)
		nbr = neighbors(E,n);
		T{n} = nbr;
	  else
		nbr = T{n};
	  end;
	  
	  mp(:) = 0;
	  for k=nbr(:)'
		mp = mp + M{k,n};
	  end;
	  
	  msgProd(:,n) = mp;
	end;
	
  end;

  MAX_ITER = 5*N;
  MSG_TOL = 1e-3;

  [E0,E1] = find(E);
  
  nEdg = length(E0);
  
  for iter=1:MAX_ITER
	
	% Random shuffling of edges
	pEdg = randperm(nEdg);
	
	converging = 1;
	l2diff = 0;
	
	for d=0:1
	  for e=pEdg(:)'
		if d
		  n0 = E0(e);
		  n1 = E1(e);
		else
		  n0 = E1(e);
		  n1 = E0(e);
		end;
		
		% Product (sum log) of messages coming into sender n0
		prod = msgProd(:,n0) - M{n1,n0};
		
		%throw caution to the wind for a 5% speedup:
		%if any(isnan(prod))
		%error('NaN in prod');
		%end;
		
		% We want to eliminate the last dimension, which will be the sender
		%
		% psi  LxL
		
		psi = S{E0(e),E1(e)};
		
		if d % sender is E0
		  psi = psi';
		end;
		
		% Local sender compatability
		phiMsgs = (H(:,n0)+prod)'; % 1xL
		
		table = psi + phiMsgs(ones(L,1),:); % LxL

		% find the largest entry for each label 
		% (of recipient)
		mx = max(table,[],2); % Lx1

		if A
		  % MAX PRODUCT
		  mNew = mx;
		  
		else
		  % SUM PRODUCT
		  
		  % We need to eliminate (marginalize) the sender from
		  % the last column, a sum in log space.
		  
		  % step 2: for accuracy, remove largest factors before exponentiating
		  intermediateTable = exp( table - mx(:,ones(1,L))); % LxL
		  
		  % step 3: marginalize/sum out child/sender label
		  marg = sum(intermediateTable,2); % Lx1
		  
		  % step 4: put it all together, noting that we haven't
		  % exponentiated the largest factors
		  mNew = mx + log(marg);
		end;
		
		% normalize
		mx = max(mNew);
		Z = mx + log(sum(exp(mNew-mx)));
		
		mNew = mNew - Z;
		
		
		if 0 % still throwing caution to the wind for another 10% speedup
		if any(isnan(mNew))
		  if marg==0
			sst = 'marg == 0';
		  elseif isnan(mx)
			sst = 'mx == NaN';
		  else
			sst = '';
		  end;
		  
		  psi
		  phiMsgs
		  table
		  
		  error(sprintf('New message is NaN: mx=%f\tmarg=%f\tln(marg)=%f',...
						mx,marg,log(marg)));
		elseif any(imag(mNew))
		  psi
		  phiMsgs
		  table
		  
		  error(sprintf('New message is complex: mx=%f\tmarg=%f\tln(marg)=%f',...
						mx,marg,log(marg)));
		end;
		end;
		
		mOld = M{n0,n1};
		
		l2diff = l2diff + sum((mNew-mOld).^2);
		
		if converging && any(abs(mNew - mOld)>MSG_TOL)
		  converging = 0;
		end;
		
		msgProd(:,n1) = msgProd(:,n1) - mOld + mNew;
		
		M{n0,n1} = mNew;
	  end;
	  
	end;
	
	%fprintf('%d\t%f\n',iter,l2diff);
	
	if converging
	  %disp(['asyncbp converged in ' num2str(iter) ' iterations']);
	  F = 1;
	  break;	  
	end;
  end;
  
  if iter==MAX_ITER
	disp(sprintf('async(%d) failed to converge',A));
	F = -1;
  end;
	
  if (A || nargout==1)
	B = beliefs(E,H,S,M,A,C,T,msgProd);
	G = {};
  else
	[B,G] = beliefs(E,H,S,M,A,C,T,msgProd);
  end;
