function [xopt,stats,options,bestf,fgen,lgen] = genetic(fun, ...          x0,options,vlb,vub,bits,P1,P2,P3,P4,P5,P6,P7P,P8,P9,P10)%GENETIC tries to maximize a function using a simple genetic algorithm.%	X=GENETIC('FUN',X0,OPTIONS,VLB,VUB) uses a simple (haploid) %       genetic algorithm to find a maximum of the fitness function %       FUN (usually an M-file: FUN.M).  The user may define all or%       part of an initial population X0 (or supply an empty argument%       in which case an initial population will be chosen randomly%       between the lower and upper bounds VLB and VUB.  Use OPTIONS%       to specify optional parameters such as population size and%       maximum number of generations produced.  Type HELP GOPTIONS%       for more information.%%       The default algorithm uses a fixed population size, OPTIONS(11),%       and no generational overlap.  Three genetic operations: %       reproduction, crossover, and mutation are performed during %       procreation.  The probability that an individual of the %       population will reproduce is proportional to its fitness.%       Individuals chosen for reproduction are mated at random.%       Mating produces two offspring (re: constant population size.)%       Crossover in mating occurs with probability Pc=OPTIONS(12)%       and the crossover index is randomly selected.  Each feature%       of the offspring can mutate independently with probability%       Pm=OPTIONS(13).  Default options are OPTIONS(11:13)=[30 1 0].%       The default maximum generations OPTIONS(14) is 100.%%	X=GENETIC('FUN',X0,OPTIONS,VLB,VUB,BITS) allows the user to %	define the number of BITS used to code non-binary parameters%	as binary strings.  Note: length(BITS) must equal length(VLB).%%	X=GENETIC('FUN',X0,OPTIONS,VLB,VUB,BITS,P1,P2,...) allows up %	to ten arguments, P1, P2, ... to be passed directly to FUN.%	F=FUN(X,P1,P2,...).%%	[X,STATS,OPTIONS,BESTF,FGEN,LGEN]=GENETIC(<ARGS>)%          STATS   - [max min mean std] for each generation%          OPTIONS - options used%          BESTF   - Fitness of indivadual X (i.e.: best fitness)%          FGEN    - first generation population%          LGEN    - last generation population%%       The algorithm implemented here is taken from Genetic%       Algorithms in Search, Optimization, and Machine Learning,%       David E. Goldberg, Addison-Wiley Publishing Company, Inc.,%       1989.%%	Copyright (c) 1993 by the MathWorks, Inc.%	Andrew Potvin 1-10-93.% OPT_STOP used by user to halt optimization prematurely% OPT_STEP will be true during the evaluation of the last individual's% cost function.  Often used to determine when to update graphics.global OPT_STOP OPT_STEPOPT_STOP = 0;% Argument and error checkingif nargin<4,   error('No population bounds given.')elseif (size(vlb,1)~=1) | (size(vub,1)~=1),   % Remark: this will change if algorithm accomodates matrix variables   error('VLB and VUB must be row vectors')elseif (size(vlb,2)~=size(vub,2)),   error('VLB and VUB must have the same number of columns.')elseif (size(vub,2)~=size(x0,2)) & (size(x0,1)>0),   error('X0 must all have the same number of columns as VLB and VUB.')elseif any(vlb>vub),   error('Some lower bounds greater than upper bounds')else   x0_row = size(x0,1);   for i=1:x0_row,      if any(x0(x0_row,:)<vlb) | any(x0(x0_row,:)>vub),         error('Some initial population not within bounds.')      end % if initial pop not within bounds   end % for initial popend % if nargin<4   if nargin<6,   bits = [];elseif (size(bits,1)~=1) | (size(bits,2)~=size(vlb,2)),   % Remark: this will change if algorithm accomodates matrix variables   error('BITS must have one row and length(VLB) columns')elseif any(bits~=round(bits)) | any(bits<1),   error('BITS must be a vector of integers >0')end % if nargin<6% Form string to call for function evaluationif ~( any(fun<48) | any(fun>122) | any((fun>90) & (fun<97)) | ...    any((fun>57) & (fun<65)) ),    % Only alphanumeric implies must be a function   evalstr = [fun '(x'];   for i=1:nargin-6,      evalstr = [evalstr,',P',int2str(i)];   end   evalstr = [evalstr, ')'];else   evalstr = fun;end% Determine all options% Remark: add another options index for type of termination criterionif size(options,1)>1,   error('OPTIONS must be a row vector')else   options = foptions(options);   if options(11)==0,      % Default size_pop      options(11) = 30;   end   if options(12)==0,      % Default Pc      options(12) = 1;   end   if options(14)==0,      % Default max_gen      options(14) = 100;   endendPRINTING = options(1);terminate = options(2);size_pop = options(11);Pc = options(12);Pm = options(13);max_gen  = options(14);% Ensure valid options: e.q. Pc,Pm,size_pop,max_gen>0, Pc,Pm<1if any([Pc Pm size_pop max_gen]<0) | any([Pc Pm]>1),   error('Some Pc,Pm,size_pop,max_gen<0 or Pc,Pm>1')end% Encode fitness (cost) function if necessaryENCODED = any(any(([vlb; vub; x0]~=0) & ([vlb; vub; x0]~=1))) |  ...          ~isempty(bits);if ENCODED,   [fgen,lchrom] = encode(x0,vlb,vub,bits);else   fgen = x0;   lchrom = size(vlb,2);end% Display warning if odd number in initial populationif rem(size_pop,2)==1,   disp('Warning: pop_size should be even.  Adding 1 to population.')   size_pop = size_pop +1;end% Form random initial population if not enough supplied by userif size(fgen,1)<size_pop,   fgen = [fgen; (rand(size_pop-size(fgen,1),lchrom)<0.5)];endxopt = vlb;bestf = -Inf;new_gen = fgen;% Header displayif PRINTING>=1,   if ENCODED,      disp('Cost function encoding as binary successful.')      disp('')      fgen = decode(fgen,vlb,vub,bits);   end   disp('                   Fitness statistics')   disp('Generation Maximum      Minimum         Mean    Std. dev.')end% Set up main loopSTOP_FLAG = 0;for generation = 1:max_gen+1,   options(10) = generation-1;   old_gen = new_gen;   % Get fitness of each string in population   % Decode first if necessary   if ENCODED,      x_pop = decode(old_gen,vlb,vub,bits);   else      x_pop = old_gen;   end   for i=1:size_pop,      x = x_pop(i,:);      if i==size_pop,         OPT_STEP = 1;      else         OPT_STEP = 0;      end      fitness(i) = eval(evalstr);   end   [max_fit,INDEX] = max(fitness);   stats = [stats; max_fit min(fitness) mean(fitness) std(fitness)];   if max_fit>bestf,      bestf = max_fit;      xopt = x_pop(INDEX(1),:);   else      % Remark: may want to regenerate to guarantee cost decrease      % Remark: be careful not to get stuck in infinite loop   end   % Display if necessary   % Remark: consider alternate printing options   if PRINTING>=1,      disp([sprintf('%5.0f %12.6g %12.6g ',generation-1, ...                     stats(generation,1),stats(generation,2)), ...            sprintf('%12.6g %12.6g  ',stats(generation,3), ...                    stats(generation,4))]);   end   % Check for termination   % Remark: add more termination options   if terminate>0,      if generation>5,         if ( (stats(generation,1)-stats(generation-5,1))/ ...               stats(generation,1)<terminate ) & ...            (stats(generation,1)>stats(generation-5,1)),            STOP_FLAG = 1;         end      elseif generation>1,         if ( (stats(generation,3)-stats(generation-1,3))/ ...               stats(generation,3)<terminate ) & ...            (stats(generation,3)>stats(generation-1,3)),            STOP_FLAG = 1;         end      end   end   if STOP_FLAG | OPT_STOP,      fprintf('\n')      if STOP_FLAG,         disp('Genetic algorithm converged.')         else         disp('Genetic algorithm terminated by user.')      end      return   end   % Reproduce   new_gen = reproduc(old_gen,fitness);   % Mate   new_gen = mate(new_gen);   % Crossover   new_gen = xover(new_gen,Pc);   % Mutate   new_gen = mutate(new_gen,Pm);   if ENCODED,      lgen = decode(new_gen,vlb,vub,bits);   else      lgen = new_gen;   endend % for max_gen% Maximum number of generations reached without terminationif PRINTING>=1,   fprintf('\n')   disp('Maximum number of generations reached without termination')   disp('criterion met.  Either increase maximum generations')   disp('or ease termination citerion.')end% end genetic