function [J,gradJ] = funcosto3(nn_params)
%
% Implements the neural network cost function for a two laye neural network 
% which performs classification
%
%   [J,grad] = NNCOSTFUNCTON(nn_params,hidden_layer_size,num_labels,X,y,lambda) 
%
%   computes the cost and gradient of the neural network. 
%   The parameters for the neural network are "unrolled" into the vector
%   nn_params and need to be converted back into the weight matrices. 
% 
%   The returned parameter grad should be a "unrolled" vector of the
%   partial derivatives of the neural network.
%
% Material del curso "Machine Learning" del Prof. Andrew NG, Coursera, Stanford University.

global n_ils  n_hls  n_lab  X  y  lambda

input_layer_size =n_ils;
hidden_layer_size=n_hls;
num_labels       =n_lab;


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

% Setup some useful variables
m = size(X, 1);
         
% You need to return the following variables correctly 
% J = 0;
% Theta1_grad = zeros(size(Theta1));
% Theta2_grad = zeros(size(Theta2));
%-------------------------
% Part 1: Feedforward the neural network and return the cost in the
%         variable J. After implementing Part 1, you can verify that your
%         cost function computation is correct by verifying the cost
%         computed in ex4.m
%

% Add ones to the X data matrix
X1 = [ones(m, 1) X];  % 1 fila por ejemplo (es A1')
%-------------
Z2=Theta1*X1';        % 1 col por ejemplo
A2=[ones(1,m); sigmoid(Z2)];
%-----
Z3=Theta2*A2;         % 1 col por ejemplo
A3=sigmoid(Z3);
%-----------
J=0;
for t=1:m,
  h=A3(:,t);
  dt=y(t);   % numero 1 a 10
  yvec=zeros(num_labels,1);
  yvec(dt)=1;
  
  J=J+sum(-yvec.*log(h)-(1-yvec).*log(1-h))/m;  % sum en k (labels)
end

% Part 2: Implement the backpropagation algorithm to compute the gradients
%         Theta1_grad and Theta2_grad. You should return the partial derivatives of
%         the cost function with respect to Theta1 and Theta2 in Theta1_grad and
%         Theta2_grad, respectively. After implementing Part 2, you can check
%         that your implementation is correct by running checkNNGradients
%
%         Note: The vector y passed into the function is a vector of labels
%               containing values from 1..K. You need to map this vector into a 
%               binary vector of 1's and 0's to be used with the neural network
%               cost function.
%
%         Hint: We recommend implementing backpropagation using a for-loop
%               over the training examples if you are implementing it for the 
%               first time.
%

grad=zeros(size(nn_params));

Delta1 = zeros(size(Theta1));
Delta2 = zeros(size(Theta2));
for t=1:m,
   %Step1: Compute activations for example t
   a1 = [1; X(t,:)'];  % ejemplo t
   z2=Theta1*a1;        
   a2=[1; sigmoid(z2)];
   z3=Theta2*a2;     
   a3=sigmoid(z3);   % valor de h para ejemplo t
   %------
   % Step 2: Compute delta3 in output layer
   dt=y(t);   % numero 1 a 10
   yvec=zeros(num_labels,1);
   yvec(dt)=1;
   delta3=a3-yvec; 
   %--------
   % Step 3: Compute delta2 in hiden layer
   aux=Theta2'*delta3;  % 1er  coord corresponde a a20=1
   delta2=aux(2:end).*sigmoidGradient(z2);
  
   %[f,g] = sigmoid(z2); delta2=aux(2:end).*g;   
   %--------
   % Step 4: Accumulate gradient from example t
   Delta2=Delta2 + delta3*a2';
   Delta1=Delta1 + delta2*a1';
   %------
   % Step5: Unregularized gradient
   Theta1_grad= Delta1/m;
   Theta2_grad= Delta2/m;
end

% Part 3: Implement regularization with the cost function and gradients.
%
%         Hint: You can implement this around the code for
%               backpropagation. That is, you can compute the gradients for
%               the regularization separately and then add them to Theta1_grad
%               and Theta2_grad from Part 2.
%

aux1=Theta1(:,2:end);
J=J+lambda*sum(sum(aux1.^2))/(2*m);

aux2=Theta2(:,2:end);
J=J+lambda*sum(sum(aux2.^2))/(2*m);
% ------------------------
Theta1_grad(:,2:end)=Theta1_grad(:,2:end) + lambda*Theta1(:,2:end)/m;

Theta2_grad(:,2:end)=Theta2_grad(:,2:end) + lambda*Theta2(:,2:end)/m;
% =========================================================================
% Unroll gradients
gradJ = [Theta1_grad(:) ; Theta2_grad(:)];