function xmin=cmaes
  % ---------------------------------------------------------------
  % CMA-ES: Evolution Strategy with Covariance Matrix Adaptation for
  % nonlinear function minimization. To be used under the terms of the
  % GNU General Public License (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
  % References: See end of file. Last change: August, 13, 2003

  % ---------------------------------------------------------------  
  % YOU NEED TO SET strfitnessfct, xmean, coordist, and stopfitness
  % (and optionally stopeval, and mindx) to problem dependent values.
  strfitnessfct = 'felli'; %'frosen'; 'fcigar'; 'fsharpR'; 'ftwoaxes'; 'fsphere'; 'ftablet'; 'fplane'; 'frand';
  xmean = 0.50*ones(20,1); N=length(xmean); % object variable start point
  coordist = 0.3;   % expected initial distance to optimum per coordinate
  stopfitness = 1e-10; % minimization: stop if fitness < stopfitness
  stopeval = 900*(N+3)^2; % stop after stopeval number of function evaluations
  mindx = 0*coordist/1e12/N^.5; % minimal coordinate standard deviation
  
  % Strategy Internal Parameter Setting: Selection  
  lambda = 4+floor(3*log(N)); mu = floor(lambda/2); 
  weights = log(mu+1)-log(1:mu)'; % muXone array for weighted recomb.
  % lambda=10; mu=2; weights = ones(mu,1); % uncomment for (2_I,10)-ES
  weights = weights/sum(weights); % normalize weight array
  mueff=sum(weights)^2/sum(weights.^2); % variance-effective size of mu
  % Strategy Internal Parameter Setting: Adaptation
  cc = 4/(N+4); mucov = mueff; ccov = 2/(N+1.41)^2;
  ccov = max(ccov, (1/mucov) * ccov + (1-1/mucov) * ...
             min(1,(2*mucov-1)/((N+2)^2+mucov))); 
  cs = 10/(N+20); 
  damp = max(0.3,1-N*lambda/stopeval) * max(1,3*mueff/(N+10)) / cs + 1;
  
  % Initialize dynamic strategy parameters and constants
  sigma = coordist*mueff/N^.5; pc = zeros(N,1); ps = zeros(N,1); 
  B = eye(N); D = eye(N); BD = B*D; C = BD*BD';
  flginiphase = 1; % Initial phase flag
  randn('state', sum(100*clock)); startseed = randn('state'); % saved below
  chiN=N^0.5*(1-1/(4*N)+1/(21*N^2)); 
  % Initialize records and output
  arfitnesshist=Inf*ones(1,3); % arfitness(1)=feval(strfitnessfct, xmean)
  outt0=clock; outetime=0.5; outx=[];outy1=[];outy2=[];outy3=[];outy4=[];
  
  % Generation loop
  disp(['  ' strfitnessfct '(N=' num2str(N) '): (' num2str(mu) ',' ...
	num2str(lambda) ')-CMA-ES (w=[' num2str(100*weights','%3.0f') ']%)' ]);
  counteval = 0; flgstop = 0;
  
  while flgstop == 0 
    
    % Generate and evaluate lambda offspring
    arz = randn(N,lambda);
    arx = xmean*ones(1, lambda) + sigma * (BD * arz);              % Eq. (13)
    % You may handle constraints here and now. You may either
    % resample columns of arz and/or multiply them with a factor
    % between -1 and 1 (the latter will decrease the overall step
    % size) and recalculate arx accordingly. Do not change arx or arz
    % in any other way. You may also copy and alter (columns of) 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));               
      counteval = counteval+1;
    end
    
    % Sort by fitness and compute weighted mean in xmean
    [arfitness, arindex] = sort(arfitness); % minimization
    arfitnesshist(2:3) = arfitnesshist(1:2);
    arfitnesshist(1) = arfitness(1);
    xold = xmean; % for speed up of Eq. (14)
    xmean = arx(:,arindex(1:mu))*weights;
    zmean = arz(:,arindex(1:mu))*weights; % == sigma^-1*D^-1*B'*(xmean-xold)
    
    % Cumulation: update evolution paths
    pc = (1-cc)*pc + (sqrt(cc*(2-cc)*mueff)/sigma) * (xmean-xold); % Eq. (14)
    ps = (1-cs)*ps + (sqrt(cs*(2-cs)*mueff)) * (B*zmean);          % Eq. (16)
    % remove momentum in ps, if ps is large and fitness is getting worse.
    if sum(ps.^2)/N > 1.5 + 10*(2/N)^.5 & ...
         arfitnesshist(2)<arfitnesshist(1) & arfitnesshist(3)<arfitnesshist(1)
      ps = ps * sqrt(N*(1+max(0,log(sum(ps.^2)/N))) / sum(ps.^2));
    end

    % Adapt covariance matrix
    if ~flginiphase % do not adapt in the initial phase 
      C = (1-ccov) * C + ccov * ((1/mucov)*pc*pc' + ...            % Eq. (15)
          (1-1/mucov)*sigma^-2 * ...  % plus additional rank mu update 
          (arx(:,arindex(1:mu))-xold*ones(1,mu)) * ...
          diag(weights) * (arx(:,arindex(1:mu))-xold*ones(1,mu))');
    end
    % Adapt sigma
    sigma = sigma * exp((norm(ps)/chiN - 1)/damp);                 % Eq. (17)
    
    % Update B and D from C
    if ~flginiphase & ccov > 0 & mod(counteval/lambda, 1/ccov/N/3) < 1
      C=triu(C)+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); 
	warning([num2str(counteval/lambda) ': condition of C ' ...
		 'at upper limit' ]);
      end
      D = diag(sqrt(diag(D))); % D contains standard deviations now
      % D = D / prod(diag(D))^(1/N);  % C = C / prod(diag(D))^(2/N);
      BD = B*D; % for speed up only
    end % if mod

    % Adjust minimal coordinate axis deviations
    if any(sigma*sqrt(diag(C)) < mindx) ...
          | any(xmean == xmean + 0.2*sigma*sqrt(diag(C))) 
      C = C + ccov * diag(diag(C) .* ...
	  (xmean == xmean + 0.2*sigma*sqrt(diag(C))));
      sigma = sigma * exp(0.05+1/damp); 
      warning([num2str(counteval/lambda) ': coordinate axis std ' ...
	       'deviation at lower limit' ]);
      % flgstop=1;
    end
    % Adjust step size in case of (numerical) precision problem 
    if arfitness(1) == arfitness(1+min(mu,ceil(lambda/4))) ...
          | all(xmean == xmean ...
	    + 0.1*sigma*BD(:,1+floor(mod(counteval/lambda,N))))
      sigma = sigma * exp(0.2+1/damp); 
      warning([num2str(counteval/lambda) ...
          ': minimal main axis sigma increased to ' ...
	  num2str(sigma*min(diag(D))) ' due to an assumed problem' ]);
      % flgstop=1;
    end
    
    % End initial phase, if step size is not much too small
    if flginiphase & counteval/lambda > min(1/cs,1+N/mucov)
      if sum(ps.^2)/N/(1-(1-cs)^(counteval/lambda)) < 1.05 
        flginiphase = 0;
      end
    end

    % Set stop flag
    if arfitness(1) <= stopfitness flgstop = 1;
    elseif counteval >= stopeval flgstop = 1;
    end

    %%% Generate output
    if 1 & (counteval < 2*lambda | (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;
	xmin = arx(:, arindex(1));
	save 'variableswcmaes.mat'; % for inspection and possible restart
	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,'%03d') '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 % if etime(...
    end % if output

  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.
                             % Notice that xmean is expected to be even
                             % better.
  
% ---------------------------------------------------------------  
function outplot(outx, outy1, outy2, outy3, outy4);
  figure(324);
  dfit = outy1(:,1)-min(outy1(:,1)); dfit(find(dfit<1e-99)) = NaN;
  subplot(2,2,1);semilogy(outx,dfit,'-c');hold on;
  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(f) (blue), f-min(f) (cyan), Sigma (green), 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=fsphere(x)
  f=sum(x.^2);
  
function f=fschwefel(x)
  f = 0;
  for i = 1:size(x,1),
    f = f+sum(x(1:i))^2;
  end

function f=fcigar(x)
  f = x(1)^2 + 1e6*sum(x(2:end).^2);
  
function f=fcigtab(x)
  f = x(1)^2 + 1e8*x(end)^2 + 1e4*sum(x(2:(end-1)).^2);
  
function f=ftablet(x)
  f = 1e6*x(1)^2 + sum(x(2:end).^2);
  
function f=felli(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=felli100(x)
  N = size(x,1); if N < 2 error('dimension must be greater one'); end
  f=1e4.^((0:N-1)/(N-1)) * x.^2;

function f=fplane(x)
  f=x(1);

function f=ftwoaxes(x)
  f = sum(x(1:floor(end/2)).^2) + 1e6*sum(x(floor(1+end/2):end).^2);

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

function f=fsharpR(x)
  f = -x(1) + 100*norm(x(2:end));
  
function f=frosen(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=fdiffpow(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=frastrigin10(x)
  N = size(x,1); if N < 2 error('dimension must be greater one'); end
  scale=10.^((0:N-1)'/(N-1));
  f = 10*size(x,1) + sum((scale.*x).^2 - 10*cos(2*pi*(scale.*x)));

function f=frand(x)
  f=rand;

% ---------------------------------------------------------------  
%%% REFERENCES
%
% The equation numbers refer to 
% Hansen, N. and A. Ostermeier (2001). Completely Derandomized Self-Adaptation 
% in Evolution Strategies. Evolutionary Computation, 9(2), pp. 159-195.
% The modification of Eq. (15) corresponds to Eq. (11) in  
% Hansen, N., S.D. Müller and P. Koumoutsakos (2003). Reducing the Time 
% Complexity of the Derandomized Evolution Strategy with Covariance Matrix 
% Adaptation (CMA-ES). Evolutionary Computation, 11(1). 
%