function [psi, phi] = slvcrf_feature(model, X, Y, G)

EI = G.edgeIdx;
D = G.D;
A = G.A;
% energy = func(w, X, Y)
% psi = d
% w1 
psi1 = X*Y; 
% w2
psi2 = 0;
for i = 1 : size(EI,1)
    if Y(EI(i,1))~=Y(EI(i,2))
        psi2 = psi2 - .5;
    end
end
% w3
psi3 = 0;
for i = 1 : size(EI,1)
    if Y(EI(i,1))~=Y(EI(i,2))
        psi3 = psi3 - .5*(2-D(EI(i,1),EI(i,2)));
    end
end
psi = [psi1;psi2;psi3];

if nargout > 1
    num_nodes = size(A,1);
    w = mat2cell(model.w,model.s,1);
    dim_fea = length(w{1});
    phi = zeros(dim_fea, num_nodes);
    for i = 1 : size(A,1)
        phi(:,i) = w{1}*Y(i);
        ni = find(A(i,:)); % neighbors of node i
        for j = 1 : length(ni)
            if Y(i)~=Y(ni(j))
                phi(:,i) = phi(:,i) + w{3}*(X(:,i)-X(:,ni(j)));
            end
        end
    end
end