function [J grad] = nnCostFunction(nn_params, ...
                                   input_layer_size, ...
                                   hidden_layer_size, ...
                                   num_labels, ...
                                   X, y, lambda)

Theta1 = reshape(nn_params(1:hidden_layer_size * (input_layer_size + 1)), ...
                 hidden_layer_size, (input_layer_size + 1));

Theta2 = reshape(nn_params((1 + (hidden_layer_size * (input_layer_size + 1))):end), ...
                 num_labels, (hidden_layer_size + 1));

% Setup some useful variables
m = size(X, 1);
J = 0;
Theta1_grad = zeros(size(Theta1));
Theta2_grad = zeros(size(Theta2));

a1 = [ones(size(X,1),1) X];
z2 = a1 * Theta1';
a2 = sigmoid(z2);
a2 = [ones(size(a2,1),1) a2];
z3 = a2 * Theta2';
result = sigmoid(z3);

new_y=zeros(m,num_labels);
for i=1:m,
    new_y(i,y(i))=1;
end
helperT1 = 2 : size(Theta1,2);
helperT2 = 2 : size(Theta2,2);
regularization = sum(sum(Theta1(:,helperT1).*Theta1(:,helperT1))) + sum(sum(Theta2(:,helperT2).*Theta2(:,helperT2)));
regularization = regularization * lambda /(2*m);
J = 1/m * sum(sum(-new_y.*log(result)-(ones(size(y))-new_y).*log(ones(size(result))-result))) + regularization;


for t = 1:m
    a1 = [1; X(t,:)'];
    z2 = Theta1 * a1;
    a2 = [1; sigmoid(z2)];
    z3 = Theta2 * a2;
    a3 = sigmoid(z3);
    yk = ( [1:num_labels] == y(t) )';
    delta_3 = a3 - yk;

    delta_2 = (Theta2' * delta_3) .* [1; sigmoidGradient(z2)];
    delta_2 = delta_2(2:end);
    Theta1_grad = Theta1_grad + delta_2 * a1';
    Theta2_grad = Theta2_grad + delta_3 * a2';
end

Theta1_grad = (1/m) * Theta1_grad + lambda * [zeros(size(Theta1, 1), 1) Theta1(:,2:end)] / m;
Theta2_grad = (1/m) * Theta2_grad + lambda * [zeros(size(Theta2, 1), 1) Theta2(:,2:end)] / m;

% Unroll gradients
grad = [Theta1_grad(:) ; Theta2_grad(:)];

end