function [evolPesos,Ep,evolSal] = bp_train(patrones,salida,validacion,NeuOcultas,funsActivacion,dersFunsActivacion,alpha,beta,errorCuadraticoMedioMaximo,iterMax,Pesos0)
%function [evolPesos,Ep,evolSal] = bp_train(patrones,salida,validacion,NeuOcultas,funsActivacion,dersFunsActivacion,alpha,beta,errorCuadraticoMedioMaximo,iterMax,Pesos0)
%
% Entrenamiento de una red hacia delante (feedforward) con algoritmo de
% aprendizaje Back Propagation.
%
% [evolPesos,Ep] = bp_train(patrones,salida,NeuOcultas,funsActivacion,dersFunsActivacion,alpha,beta,errorCuadraticoMedioMaximo,iterMax)
%
% patrones:     (NexP) Ne entradas y P patrones
% salida:       (NsxP) Ns salidas y P patrones (salida deseada)
% validacion:   (PxNe+Ns) P patrones con Ne entradas y Ns salidas
% NeuOculta:    (1xO) vector con el número de neuronas en cada capa
%               oculta. Ej:  [10 3], dos capas ocultas, con 10 y    
%               3 neuronas respectivamente.
% FunsActivacion: ((O+1)xM) funciones de activación para la 1a capa
%                 oculta, para la segunda, ..., para la de salida
% dersFunsActivacion: ((O+1)xM) derivadas de las funciones de activación 
%                 para la 1a capa oculta, para la segunda, ..., para la 
%                 de salida
% alpha:        Velocidad de aprendizaje.
% beta:         Momento
% errorCuadraticoMedioMaximo: Error máximo permitido 
% iterMax:      Número máximo de iteraciones 

% APRENDIZAJE AUTOMÁTICO (Gonzalo Pajares y Jesús M. de la Cruz, Eds.)
% Editorial: RA-MA
% Capítulo 4: Redes Neuronales Artificiales
% Autores: José Antonio López Orozco


[numEntradas, numPatrones] = size(patrones);
[numSalidas, numPatrones] = size(salida);

NOval=0;
if isempty(validacion), NOval=1; errorCuadraticoMedioValidacion=0; end;

% Topología de la red deseada : nº de neuronas por capas
% el número de capas es variable, admite tantas como se deseen
% (las conexiones son siempre de una capa a la siguiente)
%topologia = [numEntradas numNeuronasOcultas numSalidas];

topologia = [numEntradas NeuOcultas numSalidas];
numCapas = length(topologia);

rangoPesos = 1;

maxNeuronas = max(topologia);

% variaciones en pesos en cada etapa y para cada capa
% de forma que deltasPesos(x, y, zz) es el incremento en el peso que va
% de la neurona x de la capa zz a la neurona y de la zz + 1
deltasPesos(maxNeuronas+1, maxNeuronas+1, numCapas) = 0;

% cada columna es la neurona a 1 y el patrón corresp.
% tiene filas : neuronas de entrada + 1 ,  columnas : tantas como patrones
redSuma(1:topologia(1), 1:numPatrones, 1) = patrones; 
red(1:topologia(1), 1:numPatrones, 1) = patrones; 

if nargin==10,
   if sum(size(Pesos0))==2,  rand('seed',Pesos0); Pesos0=Pesos0*2-1;  end;   % Pesos aleatorios entre [-1,1]
end;

% damos valores aleatorios a los pesos, entre -rangoPesos/2 y rangoPesos/2
for capa = 1:numCapas - 1
   neuronasCapa1 = topologia(capa);
   neuronasCapa2 = topologia(capa + 1);
   pesos(1:neuronasCapa2, 1:(1+neuronasCapa1), capa) = rand(neuronasCapa2,  1+neuronasCapa1) * rangoPesos - (rangoPesos / 2);
end   

if nargin==10,
   if sum(size(Pesos0))>2,  pesos=Pesos0;  end;   
end;

iteracion = 1;
errorCuadraticoMedio = errorCuadraticoMedioMaximo + 1;

while (iteracion <= iterMax) & (errorCuadraticoMedio > errorCuadraticoMedioMaximo)
   
   % simulamos la red para todos los patrones
   
	for capa = 1:(numCapas - 1)
   
   	neuronasCapa1 = topologia(capa);
   	neuronasCapa2 = topologia(capa + 1);
   
   	% pesosCapa es la matriz bidimensional con los pesos de las neuronas de
      % la capa "capa" a las de la capa siguiente, incluyendo los pesos de la
      % neurona "siempre a 1" para calcular Theta
      pesosCapa = pesos(1:neuronasCapa2, 1:(1+neuronasCapa1), capa);
      % entradas es una matriz bidimensional con las salidas de la capa anterior y 
      % las neuronas "siempre a 1", y son todas las salidas para todos los patrones
      entradas = [ ones(1, numPatrones); red(1:neuronasCapa1, 1:numPatrones, capa) ];
      % sumamos, para todos los patrones, la multiplicación de la etapa anterior por los pesos
      redSuma(1:neuronasCapa2, 1:numPatrones, capa + 1) = pesosCapa * entradas;
      % lo mismo de antes pero aplicando ahora la función de activación
      vectAct = 0;
      funActivacion = deblank(funsActivacion(capa,:));
      for x = 1:neuronasCapa2
         for y = 1:numPatrones
            vectAct (x,y) = feval(funActivacion, redSuma(x,y,capa+1));
         end
      end
      red(1:neuronasCapa2, 1:numPatrones, capa + 1) = vectAct;
      %red(1:neuronasCapa2, 1:numPatrones, capa + 1) = feval(funActivacion,(redSuma(1:neuronasCapa2, 1:numPatrones, capa + 1)));
    end
   
   %  
    
   % calculamos el error de la última capa :
   
   % error : diferencia entre las salidas deseadas (de todos los patrones) y las obtenidas
   %         error(x,y,numCapas - 1) es el error de la neurona x de la capa de salida para el patrón y
   error(1:numSalidas, 1:numPatrones, numCapas - 1) = red(1:numSalidas, 1:numPatrones, numCapas) - salida;
   % error cuadrático medio de todos los patrones
   errorCuadraticoMedio = sum(sum(error(1:numSalidas, 1:numPatrones, numCapas - 1).^2));
   % error2 : error anterior (diferencia) multiplicada por la derivada de la función de activación
   %          luego error de las neuronas de la capa "z" (3º) para todos los patrones
   vectAct = 0;
   derFunActivacion = deblank(dersFunsActivacion(numCapas - 1,:));
   for x = 1:numSalidas
      for y = 1:numPatrones
         vectAct (x,y) = feval(derFunActivacion,redSuma(x,y,numCapas),red(x,y,numCapas));
      end
   end
   error2(1:numSalidas, 1:numPatrones, numCapas - 1) = error(1:numSalidas, 1:numPatrones, numCapas - 1) .* vectAct;
   %error2(1:numSalidas, 1:numPatrones, numCapas - 1) = error(1:numSalidas, 1:numPatrones, numCapas - 1) .* (feval(derFunActivacion, redSuma(1:numSalidas, 1:numPatrones, numCapas), red(1:numSalidas, 1:numPatrones, numCapas)));
   
   % registramos la evolución de la salida de la red, de los errores, etc.
   evolSal(iteracion,:) = red(1:numSalidas,1:numPatrones,numCapas);
   evolError(iteracion,:) = error2(1:numSalidas, 1:numPatrones, numCapas - 1);
   
   % propagamos el error hacia las capas anteriores, mapeándolo sobre toda
   % la red simultáneamente para todos los patrones
	capa = numCapas - 1;
	while capa > 1; 
   	neuronasCapa1 = topologia(capa);
   	neuronasCapa2 = topologia(capa + 1);
      
      % pesosCapa : pesos de las conexiones entre la capa "capa" y la siguiente
      %             pesosCapa(x,y) : peso de la neurona x de la capa "capa" a la y de la capa siguiente
      pesosCapa = pesos(1:neuronasCapa2, 1:(1+neuronasCapa1), capa);
      error(1:(1+neuronasCapa1), 1:numPatrones, capa - 1) = pesosCapa' * error2(1:neuronasCapa2, 1:numPatrones, capa);
      vectAct = 0;
      derFunActivacion = deblank(dersFunsActivacion(capa - 1,:));
      for x = 1:neuronasCapa1
         for y = 1:numPatrones
            vectAct (x,y) = feval(derFunActivacion, redSuma(x,y,capa), red(x,y,capa));   
         end
      end
	   error2(1:neuronasCapa1, 1:numPatrones, capa - 1) = error(2:(neuronasCapa1 + 1), 1:numPatrones, capa - 1) .* vectAct;   
      %error2(1:neuronasCapa1, 1:numPatrones, capa - 1) = error(2:(neuronasCapa1 + 1), 1:numPatrones, capa - 1) .* (feval(derFunActivacion, redSuma(1:neuronasCapa1, 1:numPatrones, capa), red(1:neuronasCapa1, 1:numPatrones, capa)));   
   	capa = capa - 1;
   end
   
	% actualizamos los pesos
	for capa = 1:(numCapas - 1)
   
   	neuronasCapa1 = topologia(capa);
   	neuronasCapa2 = topologia(capa + 1);
      
      % entradas : salida de las neuronas de esta capa para todos los patrones
   	entradas = [ ones(1, numPatrones); red(1:neuronasCapa1, 1:numPatrones, capa) ];
      % varPesos : la salida de las neuronas de la capa multiplicada por el error
      %            de la capa siguiente, multiplicada por alpha
      %				 varPesos(x,y) = variación en el enlace de la neurona y de la capa "capa" a
      %				 la neurona x de la capa siguiente; la última neurona en "y" es la que
      %				 está siempre a 1 para simular theta
      %				 la calculamos de una sola vez para todos los patrones como la suma de la
      %				 que sería para cada uno de ellos por separado
      varPesos = 0;
   	varPesos(1:neuronasCapa2,  1:(neuronasCapa1+1)) = -alpha * (error2(1:neuronasCapa2, 1:numPatrones, capa) * entradas');
      % lo que se modificaron los pesos en la etapa anterior
      deltaPesos = 0;
      deltaPesos = deltasPesos(1:neuronasCapa2,  1:(neuronasCapa1+1),  capa);
      % modificación para esta etapa
   	deltaPesos = varPesos + beta * deltaPesos;
   	deltasPesos(1:neuronasCapa2,  1:(neuronasCapa1+1),  capa) = deltaPesos;
      
      % actualizamos los pesos
 	 	  pesosCapa = pesos(1:neuronasCapa2, 1:(1+neuronasCapa1), capa);
  		pesosCapa = pesosCapa + deltaPesos;
  		pesos(1:neuronasCapa2, 1:(1+neuronasCapa1), capa) = pesosCapa;
     
	end
   
   tam = prod(size(pesos));
   evolPesos(iteracion,:) = reshape(pesos,1,tam);

   salidaRed=Simredmcapa(topologia,pesos,funsActivacion,patrones);
   errorCuadraticoMedioRed = sum(sum((salida-salidaRed).^2))/numPatrones;
         
   if NOval==0,
   % Cálculo de la salida del conjunto de validación:
   patronValidacion=validacion(:,1:numEntradas);
   salidaValidacion=validacion(:,numEntradas+1:end);
   [numEntradaValidacion,aux]=size(patronValidacion);
   salidaRed=Simredmcapa(topologia,pesos,funsActivacion,patronValidacion);
   errorCuadraticoMedioValidacion = sum(sum((salidaValidacion-salidaRed').^2))/numEntradaValidacion;
   fprintf('iteración %d :  e.c.m. : %f -- e.c.m. Val.: %f\n', iteracion, errorCuadraticoMedioRed,errorCuadraticoMedioValidacion);
   else
   fprintf('iteración %d :  e.c.m. : %f\n', iteracion, errorCuadraticoMedioRed);    
   end;
   
    Ep(iteracion,:)=[errorCuadraticoMedioRed,errorCuadraticoMedioValidacion];
    iteracion = iteracion + 1;

end
