label theend;

const	   {* Change any of these parameters to match your needs *}

   VTR = 1.0e-25;      {* "Value To Reach". devec3 will stop its minimization
                           if either the maximum number of iterations "itermax"
                           is reached or the best parameter vector "bestmem"
                           has found a value f(bestmem,y) <= VTR. *}

   D = 3;              {* number of parameters of the objective function *}

   NP = 30; 	       {* number of population members *}

   itermax = 50000;     {* maximum number of iterations (generations) *}  

   F = 0.8;            {* DE-stepsize F ex [0, 2] *}

   CR = 0.5;           {* crossover probabililty constant ex [0, 1] *}

   strategy = 8;       {* 1 --> DE/best/1/exp           6 --> DE/best/1/bin
                          2 --> DE/rand/1/exp           7 --> DE/rand/1/bin
                          3 --> DE/rand-to-best/1/exp   8 --> DE/rand-to-best/1/bin
                          4 --> DE/best/2/exp           9 --> DE/best/2/bin
                          5 --> DE/rand/2/exp           else  DE/rand/2/bin *}

   XVmin : array [1..D] of double = (-1000,-10000,-1000);
   XVmax : array [1..D] of double = (1000,1000,1000);
                      {* XVmin,XVmax   vector of lower and bounds of initial population and  *}
                      {*               futhref iterations                                    *}
                      {*    	       the algorithm seems to work well only if [XVmin,XVmax]*}
                      {*    	       covers the region where the global minimum is expected*}

  XY = 200;            {* number of points; first element of *}
                      {* INPUT DATA  *}           
                      {* load x_dos.txt -ASCII  *}
                      {* load y_dos.txt -ASCII  *}


{%      The first four arguments are essential (though they have
%       default values, too). In particular, the algorithm seems to
%       work well only if [XVmin,XVmax] covers the region where the
%       global minimum is expected. DE is also somewhat sensitive to
%       the choice of the stepsize F. A good initial guess is to
%       choose F from interval [0.5, 1], e.g. 0.8. CR, the crossover
%       probability constant from interval [0, 1] helps to maintain
%       the diversity of the population and is rather uncritical. The
%       number of population members NP is also not very critical. A
%       good initial guess is 10*D. Depending on the difficulty of the
%       problem NP can be lower than 10*D or must be higher than 10*D
%       to achieve convergence.
%       If the parameters are correlated, high values of CR work better.
%       The reverse is true for no correlation.
%
% default values in case of missing input arguments:
% 	VTR = 1.e-6;
% 	D = ?; 
% 	NP = 10*?; 
% 	itermax = D*100; 
% 	F = 0.8; 
% 	CR = 0.5; 
% 	strategy = 7;
%
% Cost function:  	function result = f(x,y);
%                      	has to be defined by the user and is minimized
%			w.r. to  x(1:D).
%
% About devec3.m
% --------------
% Differential Evolution for MATLAB
% Copyright (C) 1996, 1997 R. Storn
% International Computer Science Institute (ICSI)
% 1947 Center Street, Suite 600
% Berkeley, CA 94704
% E-mail: storn@icsi.berkeley.edu
% WWW:    http://http.icsi.berkeley.edu/~storn
%
% devec is a vectorized variant of DE which, however, has a
% propertiy which differs from the original version of DE:
% 1) The random selection of vectors is performed by shuffling the
%    population array. Hence a certain vector can't be chosen twice
%    in the same term of the perturbation expression.
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 1, or (at your option)
% any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. A copy of the GNU 
% General Public License can be obtained from the 
% Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.}




type

    popsizematrix = array[1..NP,1..D] of double;    {* matrix of size NPxD *}
    indexarray = array[1..NP] of integer;             {* index array *}

var

   a1,a2,a3,a4,a5 :   indexarray;                        {* index arrays *}
   bestmem:           array[1..1,1..D] of double;        {* best population member ever *}
   bestmemiter:       array[1..1,1..D] of double;        {* temporary best member of iteration *}
   bestval:           double;                            {* best value of member in iteration *}
   bestvalev:         double;                            {* value of best member ever *}
   bm:                popsizematrix;
   cur_gen:           integer;                                          
   DC,odch_st :       double;
   FDX, FDY, datnum:  integer;
   i,j,k,n:           integer;
   ik,iP:             integer;
   ind:               array[1..4] of integer;
   jrand:             integer;
   memberval:         array[1..NP] of double;            {* array of members values *}
   mpo:               array[1..NP,1..D] of integer;
   mui:               array[1..NP,1..D] of integer;
   muitrans:          array[1..D,1..NP] of integer;
   muitratm:          array[1..D,1..NP] of integer;
   n_pop:             integer;
   oldrand:           array[1..55] of double;
  pm1,pm2,pm3,pm4,pm5:popsizematrix;                     {* population matrix 1, 2, 3, 4, 5 *}
   pop:               popsizematrix;                     {* pop is a matrix of size NPxD *}
   popold:            popsizematrix;                     {* toggle population *}
   randmatrix:        array[1..4] of double;
   randmatrixint:     array[1..4] of integer;
   randmatrixnp:      array[1..NP] of double;
   randmatrixnpint:   array[1..NP] of integer;
   randomseed:        double;
   res:               text;
   rot:               array[1..NP] of integer;           {* rotating index array (size NP) *}
   rotd:              array[1..D] of integer;            {* rotating index array (size D) *}
   rt:                array[1..NP] of integer;           {* another rotating index array *}   
   rtd:               array[1..D] of integer;            {* rotating index array for exponential crossover *}
   tempmata:          array[1..NP] of double;
   st:                integer;
   tempmemberval:     double;
   teval:             double;
   transmata:         array[1..NP] of integer;
   ui:                popsizematrix;
   valmin:            double;
   xs,ys:             text;   
xdata,ydata,ydata_obl:array[1..XY] of double;