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

%
% Evolutionary Programing for training MLP-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=10;
lambda=1;

% Initialization;
for i=1:mu
  W{i,1}=(3/sqrt(N+1))*(2*rand(J,N+1)-1);
  W{i,2}=(3/sqrt(J+1))*(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; rhoatual=0.01;

% Environment adaptation and fitness evaluation...
for i=1:mu
  clear Waux; Waux{1}=W{i,1}; Waux{2}=W{i,2};
  % Network evaluation...
  % Hidden Layer;
  neth=Waux{1}*[X;ones(1,K)];
  Z=tanh(neth);
  % Output Layer;
  neto=Waux{2}*[Z;ones(1,K)];
  Y=tanh(neto);
  
  for k=1:K 
    % Backpropagate the Error...
    gamma=(d(:,k)-Y(:,k)).*((1+Y(:,k)).*(1-Y(:,k)));
    % Updating the heights of the hidden layer...
    dWatual{1}=((Waux{2}(:,1:J)'*gamma).*((1+Z(:,k)).*(1-Z(:,k))))*[X(:,k);1]';
    Waux{1}=Waux{1}+rhoatual*dWatual{1};
    % Updating the heights of the output layer...
    dWatual{2}=gamma*[Z(:,k);1]';
    Waux{2}=Waux{2}+rhoatual*dWatual{2};
  end
  E(i,1)=norm(d-mlpnn(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,N+1);
      W{i+f*mu,2}=W{i,2}+(E(i,1)/100)*randn(L,J+1);
      % Environment Adaptation and Fitness Evaluation of the new individuals...
      clear Waux; Waux{1}=W{i+f*mu,1}; Waux{2}=W{i+mu,2};
      % Network evaluation...
      % Hidden Layer;
      neth=Waux{1}*[X;ones(1,K)];
      Z=tanh(neth);
      % Output Layer;
      neto=Waux{2}*[Z;ones(1,K)];
      Y=tanh(neto);
  
      for k=1:K 
        % Backpropagate the Error...
        gamma=(d(:,k)-Y(:,k)).*((1+Y(:,k)).*(1-Y(:,k)));
        % Updating the heights of the hidden layer...
        dWatual{1}=((Waux{2}(:,1:J)'*gamma).*((1+Z(:,k)).*(1-Z(:,k))))*[X(:,k);1]';
        Waux{1}=Waux{1}+rhoatual*dWatual{1};
        % Updating the heights of the output layer...
        dWatual{2}=gamma*[Z(:,k);1]';
        Waux{2}=Waux{2}+rhoatual*dWatual{2};
      end
      E(i+f*mu,it)=norm(d-mlpnn(Waux,X));
      W{i+f*mu,1}=Waux{1}; W{i+f*mu,2}=Waux{2};
    end
  end
  
  % 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=mlpnn(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