  

 

  %%%%%%%%%%%%%%%%%%
  % MCRS Algorithm %
  %%%%%%%%%%%%%%%%%%

function [iterr,xmin,fmin] = mcrs2D(n, C0, C1, random, validate, costfunction, niter, eps_f)

  % mcrs2D(2, [-3;-3], [3;3],'random', 'validate', 'fitness', 1000,1e-2)
  % n - Number of parameters
  % C0,C1 - Search space X = {x in R^n: C0<=x<=C1}
  % random - Function which randomly generates any points in X
  % validate - Function which verifies the validity of x
  % if x does not belong to X: you must transform x to X
  % fitness - Function which evaluates the value of a fitness (cost) function
  % niter - Maximal number of iterations (evaluations of a cost functional)
  % eps_F - Precision in a stopping criteria 

  cesta = [cd '\'];
  global parameters;
  

    
  parameters(1,1).KP_chart = []; 
  parameters(1,1).iter = 0; 
  parameters(1,1).maxiter = niter; 
 


  % Setting used varibles and constants
  format short e
  m=max(5*n,n*n);   % Size of population
  alfa=4;  % alfa belongs to <4,8>
  if ~exist('niter','var'), niter=1e+3; end;
  if ~exist('eps_f','var'), eps_f=1.0e-10; end;
  buf=zeros(niter,n+3); indbuf=0; % Buffer with history of an optimization process
  StartTime=cputime; % Initial time

  P=zeros(m,n);  % Population
  f=zeros(m,1);  % Corresponding function values

  % Simplex
  is=zeros(n,1); % Indices of the first n points of simplex
  xt=zeros(n,1); % Center of gravity
  ip=0;          % Last point of simplex (reflexed point)

  xn=zeros(n,1); % New point
  fn=0;          % Function value

  % Other variables
  pk=zeros(m,1); % 
  vi=zeros(m,1); %
  w=zeros(n,1);  %

  iterr=0; % Number of iteration


  % Initialization of random generator
  rand('state',sum(100*clock));

  % Inititial population
  k=1;
  while k<=m
    OK=0;
    while ~OK
      w=feval(random,n,C0,C1);
      [f(k),OK]=feval(costfunction, w);
    end;  
    P(k,:)=w';
    %[k, cputime-StartTime, f(k), P(k,:)]
     
    disp(strcat(sprintf('Iter. %d, soltime: %1.4g, fmin: %1.5g, xmin: ', ...
      k, cputime-StartTime, f(k)), sprintf('%1.4e ,',P(k,:)) )); % Display
    
    %set(findobj('Tag','iterace'),'String',num2str(k));
    %drawnow;
    
    k=k+1;
  end;
  [ans,pk]=sort(f,1);
  iterr=m;


  % Improving population
  while f(pk(m))-f(pk(1))>eps_f
    vi=(1:m)';
    for k=1:n
      ni=floor((m-k+1)*rand+1);
      is(k)=vi(ni);
      vi(ni)=vi(m-k+1);
    end;
    ni=floor((m-n)*rand+1);
    ip=vi(ni);

    xt=sum(P(is,:)',2)/n;
  	Gama=alfa*rand;
    xn=xt+Gama*( xt - P(ip,:)' );
    xn=feval(validate,xn,C0,C1);
    [fn,OK]=feval(costfunction, xn);
    iterr=iterr+1;
    if OK
      if fn<f(pk(m))
        kam=1;
        for k=fliplr(1:m-1)
          if fn>=f(pk(k)),  kam=k+1; break;
     		  else
	          pom=pk(k);
	          pk(k)=pk(k+1);
	          pk(k+1)=pom;
          end;
        end;

        P(pk(kam),:)=xn';
        f(pk(kam))=fn;
        
      end;

      %if mod(iterr,100) == 0
        
        disp(strcat(sprintf('Iter. %d, soltime: %1.4g, fmin: %1.5g, xmin: ', ...
          iterr, cputime-StartTime, f(pk(1))), sprintf('%1.4e ,',P(pk(1),:)) )); % Display
    
        %set(findobj('Tag','iterace'),'String',num2str(iterr));
        %drawnow;
        
        indbuf=indbuf+1;
   	    buf(indbuf,:)=[iterr, cputime-StartTime, f(pk(1)), P(pk(1),:)];
      %end;
      if mod(iterr,100) == 0
        %save history.txt buf -ascii
      end;

    end; 
    
    if iterr==niter, break; end;

  end;

  buf=buf(1:indbuf,:);
  %save history.txt buf -ascii;
  
   beep;
   beep;
   beep;


disp('***** Optimization has finished! ***** ');

fmin = f(pk(1));
xmin = P(pk(1),:);

parameters(1,1).xmin = xmin;
parameters(1,1).fmin = fmin;
