function Current_Label = LBP(adj_list, node_pot, beta,N_Label,Index)
% format Longe
%% Definition
Massage=ones(N_Label,size(adj_list,2));
Old_Massage=zeros(N_Label,size(adj_list,2));

% adj_list = adjmatrix2list(adj_mat);

%% Loopy Blief Propagation
while 1
    %% checking convergence criteria
    if any(any(isnan(Massage)))
       error('Potential function is NaN'); 
    end
    if all(all(abs(Old_Massage-Massage)<1.0000e-05))
      break;
    end
    Old_Massage=Massage;
    Massage=Massage.*0;
    for i=1 : size(Massage,2)
        tmp_Massage=zeros(N_Label);
        Ne_Massage= Old_Massage(:,(adj_list(2,:)==adj_list(1,i)) & (adj_list(1,:)~=adj_list(2,i)));
        if ~isempty(Ne_Massage)
         MaxMassage= max(Ne_Massage,[],2).*0;
            for source=1 : N_Label
                Source_Pot = node_pot(ismember(Index,adj_list(1,i)),source);
                 Income_Massage=exp(sum(Ne_Massage(source,:)-MaxMassage(source)));
%                 Income_Massage = sum(Ne_Massage(source,:));
                for target= 1 : N_Label
                     CurrentEdge_Pot=exp(-beta*(target~=source));
%                      CurrentEdge_Pot=log(exp(-beta*(target~=source)));
                    tmp_Massage(target,source)=Source_Pot.*CurrentEdge_Pot.*Income_Massage;
%                     tmp_Massage(target,source)=Source_Pot + CurrentEdge_Pot + exp(sum(Income_Massage));
                end
            end
    %         tmp_Massage=sum(tmp_Massage,2);
            tmp_Massage=max(tmp_Massage,[],2);
            if any(any(tmp_Massage<=0))
                error('Log Zero');
            end
    %         Massage(:,i)=MaxMassage+log((1./Z).*tmp_Massage);

            Massage(:,i)=MaxMassage+log(tmp_Massage);
        else
            Massage(:,i) = Old_Massage(:,i);
        end
    end
    
%     Massage(1,:)=Massage(1,:)./sum(Massage(1,:),2);
%     Massage(2,:)=Massage(2,:)./sum(Massage(2,:),2);
    Min=repmat(min(Massage),N_Label,1);
    Max=repmat(max(Massage),N_Label,1);
    Massage=((Massage-Min.*(Min~=Max)))./((Max-Min)+((Min==Max).*Massage));
    
    Z=sum(Massage);
    Massage=Massage./repmat(Z,N_Label,1);
end
%% Computing MAP
N_node=size(unique(adj_list),1);
map=zeros(N_Label,N_node);
for i=1 : N_node
    map(:,i)= node_pot(i,:)'.*prod(Massage(:,(adj_list(2,:)== Index(i))),2);
end
[Max,Current_Label] = max(map); 

return

% function adj_list = adjmatrix2list(adj_mat)
% adj_list=[];
% for i=1:size(adj_mat,1)   
%     tmp1 = find( and( adj_mat(i,:)>0,  adj_mat(i,:)~=Inf) );
%     tmp2=[i.*ones(1,size(tmp1,2));tmp1];
%     adj_list=[adj_list,tmp2];
% end

