nnCostFunction.m

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

% Reshape nn_params back into the parameters Theta1 and Theta2, the weight matrices
% for our 3 layer neural network
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));

% initialize return variables 
m = size(X, 1);
J = 0;
Theta1_grad = zeros(size(Theta1));
Theta2_grad = zeros(size(Theta2));

% transpose X and y
% calculating labels in vector format based on observed labels in y
X = X';
y = y';
yVec = zeros(num_labels, m);
for c = 1:num_labels
    yVec(c, :) = (y==c);
end


% loop over example and calculate and accumulate cost and gradients for
% each example
for i = 1:m
    % calculating h(x(i))
    a1 = X(:, i);
    a1 = [1; a1];
    
    a2 = sigmoid(Theta1 * a1);
    a2 = [1; a2];
    
    a3 = sigmoid(Theta2 * a2);
    h = a3;
    
    % calculating cost
    J = J + (sum(-yVec(:, i) .* log(h) - (1 - yVec(:, i)) .* log(1 - h)));

    % back propagation algorithm
    delta3 = h - yVec(:, i);
    delta2 = (Theta2' * delta3) .* a2 .* (1 - a2);
    
    Theta2_grad = Theta2_grad + (delta3 * a2');
    Theta1_grad = Theta1_grad + (delta2(2:end) * a1');
end



% -------------------------------------------------------------

% =========================================================================

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


end


% % vectorized implementation
% X = X';
% y = y';
% yVec = y == (1:num_labels)';
% 
% a1 = [ones(1, size(X, 2)); X];
% 
% a2 = sigmoid(Theta1 * a1);
% a2 = [ones(1, size(X, 2)); a2];
% 
% a3 = sigmoid(Theta2 * a2);
% h = a3;
% 
% J = sum(sum(-yVec .* log(h(:, :)) - (1 - yVec).* log(1-h(:, :))));
% 
% delta3 = h - yVec;
% delta2 = ((Theta2)' * delta3) .* a2 .* (1 - a2);
% 
% Theta2_grad = delta3 * (a2)';
% Theta1_grad = delta2(2:end, :) * (a1)';
%