function xmin=wcmaes
  % CMA-ES: Evolution Strategy with Covariance Matrix Adaptation 
  % and weighted recombination for nonlinear function minimization.
  % To be used under the terms of the GNU General Public License,
  % see http://www.gnu.org/copyleft/gpl.html . 
  % Author: Nikolaus Hansen, 2001.
  % e-mail: hansen@bionik.tu-berlin.de
  % URL:http://www.bionik.tu-berlin.de/user/niko
  % Last change: February,7,2003

  % Set dimension, fitness fct, start values, stop criteria
  N=10; strfitnessfct = 'elli'; % 'rosen'; 'cigar'; 'sphere'; 'tablet'; 'sharpR';
  maxeval = 300*(N+3)^2; stopfitness = 1e-10; % stop criteria
  xmeanw = 0.099*ones(N, 1); % object parameter start point
  sigma = 0.25e-0; minsigma = 1e-15; % initial step size, minimal step size
  flginiphase = 1; % Initial phase
  
  % Parameter setting: selection,
  lambda = 4+floor(3*log(N)); mu = floor(lambda/2); 
  arweights = log((lambda+1)/2)-log(1:mu)'; % muXone array for weighted recomb.
  % lambda=10; mu=2; arweights = ones(mu,1); % uncomment for (2_I,10)-ES
  % parameter setting: adaptation
  cc = 4/(N+4); ccov = 2/(N+2^0.5)^2;
  cs = 4/(N+4); damp = (1-min(0.7,N*lambda/maxeval)) / cs + 1;

  % Initialize dynamic strategy parameters and constants
  B = eye(N); D = eye(N); BD = B*D; C = BD*transpose(BD);
  pc = zeros(N,1); ps = zeros(N,1); 
  cw=sum(arweights)/norm(arweights); chiN=N^0.5*(1-1/(4*N)+1/(21*N^2));
  % Initialize fitness for output and output
  arfitness(1)=feval(strfitnessfct, xmeanw); arx(:,1)=xmeanw; arindex=1:mu; 
  outt0=clock; outetime=0.5; outx=[];outy1=[];outy2=[];outy3=[];outy4=[];
  
  % Generation loop
  disp(['   (' num2str(mu) ',' num2str(lambda) ')-CMA-ES (w=[' num2str(arweights','%5.2f') '])' ]);
  counteval = 0; flgstop = 0;
  while 1 % if flgstop break;
    % Set stop flag
    if arfitness(1) <= stopfitness flgstop = 1;
    elseif counteval >= maxeval flgstop = 1;
    end
    
    %%% Generate output
    if 1 & (mod(counteval/lambda, 1) < 1 | flgstop)
      outx = [outx counteval];
      outy1=[outy1; [arfitness(1) sigma cond(D)]];
      outy2=[outy2;(arx(:,arindex(1)))'];
      outy3=[outy3; sigma*sqrt(diag(C))'];
      outy4=[outy4;sort(diag(D))']; 
      if etime(clock,outt0) > 10*outetime | flgstop | ...
	    etime(clock, outt0) <= 0 % bug fix for etime
	outt0 = clock;
	if 1 % graphic output
	  outplot(outx,outy1,outy2,outy3,outy4); % outplot defined below
	elseif 1 % save data for later graphic output with outplot
	  if ~exist('sfile', 'var') 
	    spath = '';
	    sfile = [spath 'cmaN' num2str(N, '%02d') 'on' strfitnessfct '.mat'];
	    if exist(sfile, 'file') error(['File ' sfile ' already exists.']);
	    end
	  end
	  save(sfile,'counteval','arfitness','arx','outx','outy1', ...
	      'outy2','outy3','outy4');
	end
	outetime = etime(clock, outt0); 
      end
    end % output
    if flgstop break;
    end

    % Generate and evaluate lambda offspring
    arz = randn(N,lambda);
    arx = xmeanw*ones(1, lambda) + sigma * (BD * arz);
    % You may handle constraints here and now. 
    % You may either resample columns of arz and/or multiply
    % them with a factor < 1 (the latter will result in a 
    % decreased overall step size) and recalculate arx
    % accordingly. Do not change arx or arz in any other
    % way. You may also use an altered arx only for the
    % evaluation of the fitness function, if you leave 
    % arx and arz unchanged for the algorithm.  
    for k=1:lambda
      arfitness(k) = feval(strfitnessfct, arx(:,k));             % Eq. (13)
      counteval = counteval+1;
    end
    
    % Sort by fitness and compute weighted mean in xmeanw
    [arfitness, arindex] = sort(arfitness); % minimization
    xold = xmeanw; % for speed up of Eq. (14)
    xmeanw = arx(:,arindex(1:mu))*arweights/sum(arweights);
    zmeanw = arz(:,arindex(1:mu))*arweights/sum(arweights);
    
    % Adapt covariance matrix
    pc = (1-cc)*pc + (sqrt(cc*(2-cc))*cw/sigma) * (xmeanw-xold); % Eq. (14)
    if ~flginiphase % do not adapt in the initial phase 
      C = (1-ccov)*C + ccov*pc*transpose(pc);                    % Eq. (15)
    end
    % adapt sigma
    ps = (1-cs)*ps + (sqrt(cs*(2-cs))*cw) * (B*zmeanw);          % Eq. (16)
    sigma = sigma * exp((norm(ps)-chiN)/chiN/damp);              % Eq. (17)
    
    % Update B and D from C
    if mod(counteval/lambda, 1/ccov/N/5) < 1
      C=triu(C)+transpose(triu(C,1)); % enforce symmetry
      [B,D] = eig(C);
      % limit condition of C to 1e14 + 1
      if max(diag(D)) > 1e14*min(diag(D)) 
        tmp = max(diag(D))/1e14 - min(diag(D));
        C = C + tmp*eye(N); D = D + tmp*eye(N); 
      end
      D = diag(sqrt(diag(D))); % D contains standard deviations now
      BD = B*D; % for speed up only
    end % if mod

    % Adjust minimal step size
    if sigma*min(diag(D)) < minsigma ...
          | arfitness(1) == arfitness(min(mu+1,lambda)) ...
          | xmeanw == xmeanw ...
              + 0.2*sigma*BD(:,1+floor(mod(counteval/lambda,N)))
      sigma = 1.4*sigma; 
      warning(['minimal axis sigma increased to ' num2str(sigma*min(diag(D))) ...
	     ' due to an assumed problem' ]);
      % flgstop=1;
    end

    % Test for end of initial phase
    if flginiphase & counteval/lambda > 2/cs
      if (norm(ps)-chiN)/chiN < 0.05 % step size is not much too small
	flginiphase = 0;
      end
    end
    
  end % while, end generation loop

  disp([num2str(counteval) ': ' num2str(arfitness(1))]);
  if exist('sfile', 'var') disp(['Results saved in ' sfile]); 
  end
  xmin = arx(:, arindex(1)); % return best point of last generation 

% ---------------------------------------------------------------  
function outplot(outx, outy1, outy2, outy3, outy4);
  figure(324);
  subplot(2,2,1);semilogy(outx,abs(outy1(:,1))+1e-99,'.-b');hold on;
  subplot(2,2,1);semilogy(outx,(outy1(:,2)),'-g'); hold on;
  subplot(2,2,1);semilogy(outx,(outy1(:,3)),'-r'); hold off;
  title('abs(Function Value), Sigma (gr.), Axis Ratio (red)');grid on;zoom on; 
  subplot(2,2,2); plot(outx,outy2,'-');
  title('Object Parameters');grid on;zoom on; 
  subplot(2,2,3); semilogy(outx, outy3, '-'); 
  title('Standard Deviations in All Coordinates');grid on;zoom on; 
  subplot(2,2,4); semilogy(outx, outy4, '-');
  title('Scaling (All Main Axes)');grid on;zoom on; drawnow;
  
% ---------------------------------------------------------------  
function f=sphere(x)
  f=sum(x.^2);
  
function f=schwefel(x)
  f = 0;
  for i = 1:size(x,1)
    f = f+sum(x(1:i))^2;
  end

function f=cigar(x)
  f = x(1)^2 + 1e6*sum(x(2:end).^2);
  
function f=tablet(x)
  f = 1e6*x(1)^2 + sum(x(2:end).^2);
  
function f=elli(x)
  N = size(x,1); if N < 2 error('dimension must be greater one'); end
  f=1e6.^((0:N-1)/(N-1)) * x.^2;

function f=parabR(x)
  f = -x(1) + 100*sum(x(2:end).^2);

function f=sharpR(x)
  f = -x(1) + 100*norm(x(2:end));
  
function f=rosen(x)
  if size(x,1) < 2 error('dimension must be greater one'); end
  f = 100*sum((x(1:end-1).^2 - x(2:end)).^2) + sum((x(1:end-1)-1).^2);

function f=diffpow(x)
  N = size(x,1); if N < 2 error('dimension must be greater one'); end
  f=sum(abs(x).^(2+10*(0:N-1)'/(N-1)));
  
function f=rastrigin10(x)
  N = size(x,1); if N < 2 error('dimension must be greater one'); end
  scale=10.^((0:N-1)'/(N-1));
  q = 10*size(x,1) + sum((scale.*x).^2 - 10*cos(2*pi*(scale.*x)));

function f=fct_rand(x)
  f=rand;