function [LL,g,h]= LL_poiss_MATlike(Y,Xa,Xnb,W,EstimOpt,OptimOpt,b0)

% save tmp_LL_lc_MATlike;
% return

LLfun = @(B) LL_poiss(Y,Xa,Xnb, EstimOpt,B);

% if EstimOpt.Test == 1
%     [f,j] = LLfun(b0);
%             g = sum(j,1)'; ...
%     EstimOpt.NumGrad = 1;
%     LLfun2 = @(B) LL_poiss(Y,Xa,Xnb, EstimOpt,B);
%     f2 = LLfun2(b0);  
%             j2 = numdiff(LLfun2,f2,b0,1,EstimOpt.BActive);...
%             g2 = sum(j2,1)';  
%     [g,g2, abs(g-g2)]
%    % abs(j(:,4) - j2(:,4))
%    % abs(j-j2)
%    % pause;
%     EstimOpt.NumGrad = 0;
% end

if isequal(OptimOpt.GradObj,'on')
    if EstimOpt.NumGrad == 0
        [f,j] = LLfun(b0);
        j(:,EstimOpt.BActive ==0) = 0;
        f = f.*W;
        j = j.*W(:, ones(1, size(j,2)));
        g = sum(j,1)'; ...
        if isequal(OptimOpt.Hessian,'user-supplied') == 1
            h = j'*j;
        end
    else % => EstimOpt.NumGrad == 1 
        f = LLfun(b0);  
        j = numdiff(LLfun,f,b0,isequal(OptimOpt.FinDiffType,'central'),EstimOpt.BActive);...
        f = f.*W;
        j = j.*W(:, ones(1, size(j,2)));
        g = sum(j,1)';   
        if isequal(OptimOpt.Hessian,'user-supplied') == 1
            h = j'*j;
        end
    end
    
else % No gradient
    EstimOpt.NumGrad = 1;
    f = LLfun(b0);   
end
LL = sum(f);
