function [Waux,Error,ErrorMin]=mrlLamark(netp,X,d,Xv,dv,Xt,dt);

%
% Evolutionary Programing for training MRL-NNs (Porto-Fogel Approach with Lamarckian Ideas).
%
% 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 set;
%              dv - Validation Target;
%              Xt - Input test samples;
%              dt - Target test;
%
% Output: W - Cell with the synaptic weight vectors;
%     Error - Vector with the MSE by epoch;
%

[N,K]=size(X);

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

% Population Size...
mu=50;
lambda=1;

% Initialization;
for i=1:mu
  W{i,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{i,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)];
end

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

it=0; Terminate=0;

% Environment Adaptation and Fitness Evaluation...
for i=1:mu
  clear Waux; Waux{1}=W{i,1}; Waux{2}=W{i,2};
   
  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=Waux{1}(:,3:N+2)+ones(J,1)*X(:,k)';
    alphah=sort(Ax',1);
    clear ind; ind=round(N-(N-1)./(1+exp(Waux{1}(:,2))));
    alphah=alphah(ind+[0:N:(J-1)*N]');
    Qh=(alphah*ones(1,N)==Ax);
    SQh=sum(Qh,2);
    betah=Waux{1}(:,N+3:2*N+3)*[X(:,k);1];
    neth=Waux{1}(:,1).*alphah+(1-Waux{1}(:,1)).*betah;
    Z=tanh(neth);
    % Output Layer...
    Az=Waux{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);
    SQo=sum(Qo,2);
    betao=Waux{2}(:,J+3:2*J+3)*[Z;1];
    neto=Waux{2}(:,1).*alphao+(1-Waux{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;
    
    Waux{1}=Waux{1}+rhoatual*dWatual{1};
    Waux{2}=Waux{2}+rhoatual*dWatual{2};
  end  
  E(i,1)=norm(d-mrlnn(Waux,X));
  W{i,1}=Waux{1}; W{i,2}=Waux{2};
end

while (it<=itmax)&(~Terminate)
  % Iteration Counter...
  it=it+1;
   
  % For each individual...
  for i=1:mu
    % Reproduction...
    for f=1:lambda
      W{i+f*mu,1}=W{i,1}+(E(i,1)/100)*randn(J,2*N+3);
      W{i+f*mu,2}=W{i,2}+(E(i,1)/100)*randn(L,2*J+3);
      % Environment Adaptation and Fitness Evaluation of the new individuals...
      clear Waux; Waux{1}=W{i+f*mu,1}; Waux{2}=W{i+mu,2};
      
      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=Waux{1}(:,3:N+2)+ones(J,1)*X(:,k)';
        alphah=sort(Ax',1);
        clear ind; ind=round(N-(N-1)./(1+exp(Waux{1}(:,2))));
        alphah=alphah(ind+[0:N:(J-1)*N]');
        Qh=(alphah*ones(1,N)==Ax);
        SQh=sum(Qh,2);
        betah=Waux{1}(:,N+3:2*N+3)*[X(:,k);1];
        neth=Waux{1}(:,1).*alphah+(1-Waux{1}(:,1)).*betah;
        Z=tanh(neth);
        % Output Layer...
        Az=Waux{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);
        SQo=sum(Qo,2);
        betao=Waux{2}(:,J+3:2*J+3)*[Z;1];
        neto=Waux{2}(:,1).*alphao+(1-Waux{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;
    
        Waux{1}=Waux{1}+rhoatual*dWatual{1};
        Waux{2}=Waux{2}+rhoatual*dWatual{2};
      end
  
      E(i+f*mu,it)=norm(d-mrlnn(Waux,X));
      W{i+f*mu,1}=Waux{1}; W{i+mu,2}=Waux{2};
    end
  end
  
  % Selection...
  
% Selection...
  %[aux,ind]=sort(E(:,it));
  
  clear wins;
  wins=zeros(size(E,1),1);
  for i=1:20*mu
    ind=ceil(2*mu*rand(2,1));
    if E(ind(1),it)<=E(ind(2),it)
      wins(ind(1))=wins(ind(1))+1;
    else
      wins(ind(2))=wins(ind(2))+1;
    end
  end
  clear ind;
  [wins,ind]=sort(-wins);
  
  clear Waux;
  for i=1:mu
    Waux{i,1}=W{ind(i),1};
    Waux{i,2}=W{ind(i),2};
    E(i,it+1)=E(ind(i),it);
  end
  clear W; W=Waux;
  
  % Evaluation of the Samples...
  clear Waux; Waux{1}=W{1,1}; Waux{2}=W{1,2};
  Y=mrlnn(Waux,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