function [W,Error,ErrorMin,rho_vet]=mrlbackpro(netp,X,d,Xv,dv,Xt,dt);

%
% Back-propagation algorithm for training MRL-NNs.
%
% Inputs:    netp - A vector [m,n], where m is the number of hidden units and n is the number of output units;  
%               X - Training samples. Each pattern is a column of X;
%               d - Target or Desired System output (row vector);
% (Optional):  Xv - Validation samples;
%              dv - Validation Target;
%              Xt - Test samples;
%              dt - Test Target;
%
% Output: W - Cell with the synaptic weight vectors;
%     Error - Vector with the MSE by epoch;
%

J=netp(1); L=netp(2); 
if length(netp)>2
  itmax=netp(3);
else
  itmax=100;
end 
clear netp;

[N,K]=size(X);

% Initialization;
W{1}=[0.2*rand(J,1)+0.4,rand(J,1)-0.5,0.2*rand(J,N)-0.1,1/sqrt(N)*(2*rand(J,N+1)+1)];
W{2}=[0.2*rand(L,1)+0.4,rand(L,1)-0.5,0.2*rand(L,J)-0.1,1/sqrt(J)*(2*rand(L,J+1)+1)];
clear Waux;

rho_vet(:,1)=[W{1}(:,2);W{2}(:,2)];

Kv=0; Kt=0;
ErrorMin=inf;
if nargin>3
  X=[X,Xv];
  Kv=size(Xv,2);
  clear Xv;
  ErrorMin=ones(2,1)*inf;
  if nargin>5
    X=[X,Xt];
    Kt=size(Xt,2);
    clear Xt;
    ErrorMin=ones(3,1)*inf;
  end
end

% Backprop parameters
rhoatual=0.01;
rhoanterior=rhoatual*0.8;

it=0; Terminate=0;

dWanterior{1}=zeros(J,2*N+3);
dWanterior{2}=zeros(L,2*J+3);

while (it<=itmax)&(~Terminate)
  % Iteration Counter...
  it=it+1;
  
  Kind=randperm(K);
  for i=1:K
    % Select k;
    k=Kind(i);
  
    % Network evaluation...
    clear Ax; clear Az;
    clear alphah; clear alphao;
    clear Qh; clear Qo;
    clear SQh; clear SQo;
    clear neto; clear neth;
    clear Z; clear Y;
  
    % Hidden Layer...
    Ax=W{1}(:,3:N+2)+ones(J,1)*X(:,k)';
    alphah=sort(Ax',1);
    clear ind; ind=round(N-(N-1)./(1+exp(W{1}(:,2))));
    alphah=alphah(ind+[0:N:(J-1)*N]');
    Qh=(alphah*ones(1,N)==Ax);
    %Qh=exp(-0.5*((alphah*ones(1,N)-Ax)/0.05).^2);
    SQh=sum(Qh,2);
    betah=W{1}(:,N+3:2*N+3)*[X(:,k);1];
    neth=W{1}(:,1).*alphah+(1-W{1}(:,1)).*betah;
    Z=tanh(neth);
    % Output Layer...
    Az=W{2}(:,3:J+2)+ones(L,1)*Z';
    alphao=sort(Az',1);
    clear ind; ind=round(J-(J-1)./(1+exp(W{2}(:,2))));
    alphao=alphao(ind+[0:J:(L-1)*J]');
    Qo=(alphao*ones(1,J)==Az);
    %Qo=exp(-0.5*((alphao*ones(1,J)==Az)/0.05).^2);
    SQo=sum(Qo,2);
    betao=W{2}(:,J+3:2*J+3)*[Z;1];
    neto=W{2}(:,1).*alphao+(1-W{2}(:,1)).*betao;
    Y=tanh(neto);
    clear Az; clear Ax;
    
    % Back-propagate the error...
    % Output Layer...
    deltao=(d(:,k)-Y).*((1+Y).*(1-Y));
    gammao=[alphao-betao,W{2}(:,1).*(1-SQo/J),diag(W{2}(:,1)./SQo)*Qo,(1-W{2}(:,1)).*[Z',1]];
    % Hidden Layer...
    deltah=((gammao(:,3:J+2)+diag(1-W{2}(:,1))*W{2}(:,J+3:2*J+2))'*deltao).*((1+Z).*(1-Z));
    gammah=[alphah-betah,W{1}(:,1).*(1-SQh/N),diag(W{1}(:,1)./SQh)*Qh,(1-W{1}(:,1))*[X(:,k)',1]];
  
    % Updating...
    dWatual{1}=diag(deltah)*gammah;
    dWatual{2}=diag(deltao)*gammao;
    
    W{1}=W{1}+rhoatual*dWatual{1}+rhoanterior*dWanterior{1};
    W{2}=W{2}+rhoatual*dWatual{2}+rhoanterior*dWanterior{2};
    
    dWanterior{1}=dWatual{1};
    dWanterior{2}=dWatual{2};
  end
  rho_vet(:,it+1)=[W{1}(:,2);W{2}(:,2)];
  
  % Computing the Error...
  Y=mrlnn(W,X);
  Error(it,1)=norm(d-Y(:,1:K));
  disp(sprintf('Iteration: %d, Training Error = %2.4f (%2.4f)', it, Error(it,1), sum(abs(ceil(d)-ceil(Y(:,1:K))))))
end  

rho_vet(1:J,:)=round(N-(N-1)./(1+exp(rho_vet(1:J,:))));
rho_vet(J+1,:)=round(J-(J-1)./(1+exp(rho_vet(J+1,:))));