tree_ga2.m

王小平《遗传算法——理论、应用与软件实现》随书光盘

%
%                 Tree Structured Genetic Algorithm
%
%      Note: Need to run ga_ini to initialise parameters for tree_ga.
%
            
%                      Ben.McKay & Mark Willis
%                             (c) 1995
%
%                 Dept. of Chemical Engineering               
%                       University of Sydney
%
%                               &
%
%            Dept. of Chemical and Process Engineering
%                University of Newcastle-Upon-Tyne
%
%       This package may not be modified, reproduced or sold
%              without the permission of the authors.
%

% Author             Date     Predecessor   Modification
% ======             ====     ===========   ============
% B.McKay M.Willis   17/5/95  tree_ga.m     user def. func.
%                    25/5/95                fitness proportionate selection
% B.McKay             7/8/95  tree_ga2.m    permutation operator
%		     24/10/95 tree_ga2.m    replaced strrep6.m with strep2.mex 
%					    routine
%
% Function Calls: pcincl.m, equ_gen3.m, opt_par.m, mutato3.m, x_over2.m,
% ==============  fit_fun2.m, strrem.m, strep2.m, fps.m, lin_seed, disp_res,
%                 permut.m
%    
% Last Modification: 3/10/95
% =================
%

a=0;
count=0;
results=[0 0 0];
clear s p r l;
% Set Random number seeds
rand('seed',sum(100*clock)) 
randn('seed',sum(100*clock)) 

%--------------------------------------------------------------------------------
% Initialise Graphics Screen
%-------------------------------------------------------------------------------- 

% Clear the screen
da_front
main=w1;
drawnow;

% Plot a couple of axes

ax1=axes(...
        'Units','pixels',...
        'Position',[60 270 520 130],...
        'Box','on',...
        'Color',[0 0 0],...
        'Visible','on');

ax2=axes(...
        'Units','pixels',...
        'Position',[60 110 230 100],...
        'Box','on',...
        'Color',[0 0 0],...
        'Visible','on');

ax3=axes(...
        'Units','pixels',...
        'Position',[350 110 230 100],...
        'Box','on',...
        'Color',[0 0 0],...
        'Visible','on'); 
drawnow;

%--------------------------------------------------------------------------------
% Initialise strings for ad_just2
%--------------------------------------------------------------------------------

% Maximum number of solutions in population prior to selection
%  => pop_size

y_front='y_est(:,1)=';
colons=';';
for i=2:pop_size,
   y_front=str2mat2(y_front,['y_est(:,' num2str(i) ')=']);
   colons=str2mat2(colons,';');
end
  
  
%--------------------------------------------------------------------------------  
% Set up the input (u) and output (y) variables
%--------------------------------------------------------------------------------

[D L]=size(output_var);
for i = 1 : L
   if output_var(i) == 1
      y=data(:,i);
      break;
   end
end


% routine to fill input data matrix

u=pcincl(data,include_var,output_var);
[u y uv vy]=psplit(u,y,tr_split);
[num_data,num_inp]=size(u);

stdy = std(y);

%--------------------------------------------------------------------------------  
% Seed population with a linear model
%--------------------------------------------------------------------------------  

lin_seed

%load c:\windows\desktop\data\weetabix\seed_pop
%seed_pop=seed_pop(1:50,:);

[seed_num,xxx]=size(seed_pop);
seed_flag=1;
%seed_num=0;
%seed_flag=0;

%--------------------------------------------------------------------------------  
% Initialise Inputs
%--------------------------------------------------------------------------------  

for i=1:num_inp,
         eval(['u' num2str(i) '= u(:,' num2str(i) ');']);
end

%--------------------------------------------------------------------------------  
% Generate initial population
%--------------------------------------------------------------------------------  

if seed_flag==0
   pop = equ_gen3(s_length,num_inp,equ_par,fun_set,num_fun_arg_set,user_def_fun);
   pop=opt_par(pop,u,y,num_inp,options,time_out);
else
   pop = seed_pop;
end

for i=seed_num+1:pop_size,

   x2 = equ_gen3(s_length,num_inp,equ_par,fun_set,num_fun_arg_set, user_def_fun);
   x2=opt_par(x2,u,y,num_inp,options,time_out);
   pop = str2mat2(pop, x2);
   
end


%--------------------------------------------------------------------------------  
% Initialise Plot
%--------------------------------------------------------------------------------  

theta=[0 0]';
text3=ga_text(main,210,410,200,15,'0 - Generations',[1 1 1],[0 0 0]);
text1=ga_text(main,20,35,600,15,'Best Tree from current generation:',[1 1 1],[1 0 0]);
text2=ga_text(main,20,20,600,15,'',[1 1 1],[1 0 0]);

text6=ga_text(main,20,2,145,15,'Fitness Value:',[1 1 1],[1 0 0]);
text7=ga_text(main,165,2,145,15,'',[1 1 1],[1 0 0]);
text8=ga_text(main,300,2,145,15,'RMS:',[1 1 1],[1 0 0]);
text9=ga_text(main,435,2,185,15,'',[1 1 1],[1 0 0]);


MAX_FIT=[];
run_pop=[];

%--------------------------------------------------------------------------------  
% Generation Loop
%--------------------------------------------------------------------------------  

for i=1:num_gen,
  [Q,R]=size(pop);
  disp(['Generation ' num2str(i)])

 
%
% Determine fitness of each equation in the population  
%
 
  
  [fit,pop,mod_size,str_len]=fit_fun2(ffunc,pop,pop_size,u,y,num_data,num_inp,stdy,y_front,colons,snip_length);
   
     
  [max_fit,max_fit_ind]=max(fit);
  [min_fit,min_fit_ind]=min(fit);

  set(text3,'string',[num2str(i) ' - Generations']);


%
% Prepare best string for evaluation, and calculate theta values
%

out=['y_e = ' pop(max_fit_ind,:) ';'];


% Remove '{', '}', '<1', '<2', '<3' and '>'  
out=strrem(out);      

out = strep2(out);

eval(out);
      
if ffunc=='ccf',
   dta=[y_e ones( num_data,1)];
   theta=inv(dta'*dta)*dta'*y;
   y_est1=y_e;
   y_p1=dta*theta;

   y_p=y_p1;
   
elseif ffunc=='mse',
   y_p = y_e;
end


count = count +1;
drawnow; 
a=fit(max_fit_ind);       
if count >= presults,

   %
   % Display results
   %

   disp_res
   count=0;
end


%
% Check for save interval
%
  if rem(i,dump_num) == 0,
    if i==dump_num,
       save_theta=theta';
       save_pop=pop(max_fit_ind,:);
       save_fit=a;
    else
       save_theta=[save_theta; theta'];
       save_pop=str2mat2(save_pop,pop(max_fit_ind,:));
       save_fit=[save_fit; a];
    end
  save save_pop save_theta save_pop save_fit
  end 

   % Determine size of Next Generation

   new_size=floor((1-p_old)*mod_size);

   %
   % Selecting and applying genetic opperators     
   %

  new_pop=[];
  j=0;
  while j < new_size,

     %
     % Select a genetic opperator
     %
     
     p_gen=rand(1);
     
     if p_gen < p_mutate,                % Mutation
        fps_num=fps(fit);
        
        deblank_pop = deblank(pop(fps_num,:));
        pop_mut = mutato3(deblank_pop,u,y,num_inp,p_c2u,mut_fun_wgt,...
  mut_fun_set,mut_num_fun_arg_set,fun_name,num_fun_arg,user_def_fun,equ_par);
        if strcmp(deblank(pop(fps_num,:)),pop_mut) == 0,
           pop_mut=opt_par(pop_mut,u,y,num_inp,options,time_out);
        end
        
        j=j+1;
        if j==1,
           new_pop=pop_mut;
        else   
           new_pop=str2mat2(new_pop,pop_mut);
        end
     elseif p_gen < p_mutate + p_cross,  % Cross Over
        x_fps_num=fps(fit);
        y_fps_num=fps(fit);
        
        [xnew,ynew]=x_over2(deblank(pop(x_fps_num,:)),deblank(pop(y_fps_num,:)));
        xnew=opt_par(xnew,u,y,num_inp,options,time_out);
        ynew=opt_par(ynew,u,y,num_inp,options,time_out);
        
        if j < new_size-1,
           j=j+2;
           if j==2,
              new_pop=xnew;
              new_pop=str2mat2(new_pop,ynew);
           else   
              new_pop=str2mat2(new_pop,xnew,ynew);
           end
        else
           j=j+1;
           new_pop=str2mat2(new_pop,xnew); 
        end
     elseif p_gen < p_mutate + p_cross + p_perm, % Permutation
        fps_num=fps(fit);
        
        deblank_pop = deblank(pop(fps_num,:));
        pop_perm = permut(deblank_pop,num_inp,fun_set,num_fun_arg_set,fun_name,...
                     num_fun_arg,user_def_fun);
 
        if strcmp(deblank(pop(fps_num,:)),pop_perm) == 0,
           pop_perm=opt_par(pop_perm,u,y,num_inp,options,time_out);
        end
        
        j=j+1;
        if j==1,
           new_pop=pop_perm;
        else   
           new_pop=str2mat2(new_pop,pop_perm);
        end
     else                                % Direct Reproduction
        fps_num=fps(fit);
     
        j=j+1;
        if j==1,
           new_pop=pop(fps_num,:)
        else
           new_pop=str2mat2(new_pop,pop(fps_num,:));
        end        
     end
 
  end


  % Sort old population

  sort_pop = pop;
  sort_fit = fit;
  for j=1:pop_size,
     [max_fit,max_fit_ind]=max(sort_fit);
     sort_pop(j,:) = pop(max_fit_ind,:);
     sort_fit(max_fit_ind) = -inf;
  end

  %
  % Add best member to run_pop
  %
  if i==1
	run_pop=sort_pop(1,:);
  else
	run_pop=str2mat(run_pop,sort_pop(1,:));
  end

  % Add surviving members of old population  
  
  keep_pop=sort_pop(1:pop_size-new_size,:);
  pop=str2mat(new_pop,keep_pop);
 
%
% Check for reinitialisation
%
  if rem(i,rein_num) == 0,
      best_pop=keep_pop(1,:);
      clear pop
      pop=best_pop; 
     for j=2:pop_size,
        x2 = equ_gen3(s_length,num_inp,equ_par,fun_set,num_fun_arg_set,user_def_fun)
        x2=opt_par(x2,u,y,num_inp,options,time_out);
        pop = str2mat2(pop, x2);
     end   
  end 
 
  clear new_pop keep_pop
 
%
% Shuffle population
%

pop=shuffle(pop);
 
end

%
%  Show analysis screen
%
da_gpev;
%
% End     
%