function [J, grad] = costFunc(Theta_init, X, y, lambda)
%costFunc computes the cost and gradient of cross-entropy loss
    
    %Theta: d*q, X: m*d, y: m*q 
    m = size(X,1);
    d = size(X,2);
    q = size(y,2);
    Theta_init = reshape(Theta_init,d,[]);
    Theta = [Theta_init, zeros(d,1)];
    
    %objective function
    linR = X*Theta;%m*q
    eR = exp(linR);%m*q
    eR_sumrow = sum(eR,2);%m*1
    eR_sum_row_mtx = repmat(eR_sumrow,1,q);%m*q
    softR = eR./eR_sum_row_mtx;%m*q
    lnR = log(softR);%m*q
    yR = y.*lnR;%m*q
    %J = -sum(yR(:)) + lambda/2*sum(Theta(:).^2);
    J = -sum(yR(:)) + m*lambda/2*sum(Theta(:).^2);
    
    %gradient function
    grad = zeros(d,q-1);
    for tt=1:q-1
        I_mtx = zeros(m,q);%m*q
        I_mtx(:,tt) = ones(m,1);
        softRt_mtx = repmat(softR(:,tt),1,q);%m*q
        I_m_eR = I_mtx - softRt_mtx;%m*q
        yIR = y.*I_m_eR;%m*q
        yIR_sumrow = sum(yIR,2);%m*1
        yIR_repmat = repmat(yIR_sumrow,1,d);%m*d
        Xy = X.*yIR_repmat;%m*d
        grad_t = sum(Xy);%1*d
        %grad(:,tt) = -grad_t' + lambda*Theta(:,tt);%d*1
        grad(:,tt) = -grad_t' + m*lambda*Theta(:,tt);%d*1
    end
    grad = grad(:);        
end