function [ll,gradll]=volElemLoglikelihood(lambda,X,alphaNorms,blockSizes)
% [ll,gradll]=volElemLoglikelihood(X,lambda)
%
% This function computes the loglikelihood and the gradient of the
% loglikelihood function under the inverse volume element model with
% parameter vector lambda and data X.
%
%  X         A sparse matrix whose rows are examples - 
%                 TF representation of documents.
%  lambda    The model parameter as a row vector.
%
% Important Note: 
%    It is assumed that the number of dimensions 
%          size(X,2)=size(lambda,2) is even.

if size(lambda,2)==1, lambda=lambda';end
numExamples = size(X,1); n = size(X,2)-1; k=(n+1)/2;
if rem(k,1)~=0, disp('Error: k not integer!'); return;end
if size(X,2)~=size(lambda,2);disp('Error: Size mismatch!'); return;end

lamMat=  spones(X) * spdiags(lambda',0,n+1,n+1);
XdotLambda = full(sum(lamMat .* X ,2)) + alphaNorms;

%Compute the gradient (row vector) of the first loglikelihood term
gradll = full(sum( spdiags((1./XdotLambda),0,numExamples,numExamples) * X,1)) ...
    + sum(alphaNorms./XdotLambda,1);

if lambda==ones(size(lambda))*lambda(1),
    gradll = k*gradll  - numExamples * ones(size(lambda)) * ...
        (k/lambda(1))*exp(gammaln(5/2)-gammaln(3/2)+gammaln(3*k)-gammaln(3*k+1));
    Z = 0;
    return;
end

computeForwardBackwardVec(lambda,'forward',blockSizes);
computeForwardBackwardVec(lambda,'backward',blockSizes);
[Z,gradZ]=computeGradBackward(lambda,blockSizes);

ll = k * sum(log(XdotLambda),1) -  numExamples * log(Z);
gradll = k * gradll - numExamples * gradZ / Z;