function [Lhat, belief, edgeIdx] = crf_infer_lbp(model, X, G, L, bp_pars)

if ~exist('bp_pars','var')
    bp_pars.maximize = 0;
    bp_pars.max_iters = 200;
    bp_pars.tol = 1e-3;
end

w = mat2cell(model.w, model.s, 1);
w1 = w{1};
w2 = w{2};
w3 = w{3};

[dim_feature, num_nodes] = size(X);

% node potentials
node_energy = -w1' * X;
if isempty(L)
    node_potentials_p = exp(-node_energy);
    node_potentials_n = exp(node_energy);
else
    node_potentials_p = exp(-node_energy+double(L'~=1));
    node_potentials_n = exp(node_energy+double(L'~=-1));
end
node_potentials = [node_potentials_p' node_potentials_n'];

V = G.V; E = G.E; edgeIdx = G.edgeIdx; num_edges = G.num_edges;

edge_energy = zeros(1, num_edges);
for i = 1 : num_edges
    diff = X(:,edgeIdx(i,1)) - X(:,edgeIdx(i,2));
    edge_energy(i) = 0.5*(w2 - w3*sum(diff.^2));
%     edge_energy(i) = 0.5*(w2+w3*2 - w3*sum(diff.^2));
end

edge_potentials_pp = exp(zeros(1,num_edges));
edge_potentials_pp = exp(edge_energy);
edge_potentials_pn = exp(-edge_energy);
edge_potentials_np = edge_potentials_pn;
edge_potentials_nn = edge_potentials_pp;

edge_potentials = zeros(2, 2, num_edges);
edge_potentials(1, 1, :) = edge_potentials_pp;
edge_potentials(1, 2, :) = edge_potentials_pn;
edge_potentials(2, 1, :) = edge_potentials_np;
edge_potentials(2, 2, :) = edge_potentials_nn;

node_potentials = double(node_potentials);
edge_potentials = double(edge_potentials);

[belief, num_iter, msgs, edge_belief] = BP_General_C(...
        edge_potentials, node_potentials, bp_pars.maximize, bp_pars.max_iters,...
        bp_pars.tol, int32(V), int32(E));

[junk, Lhat] = max(belief, [], 2);
Lhat = -Lhat*2 + 3;
    
    
    
    
    
    
    
    
    
    