%This program uses Genetic Algorithm Toolbox to find Thicknesses and 
%Shear Wave Velocities (vs) by forward solution of Rayleigh function using 
%experimental dispersion data. 
%Genetic Algorithm Toolbox by: Chipperfield, A., P. Fleming, H. Pohlheim, 
%and C. Fonseca (1994).
%Forward Solution of Rayleigh Function by: Rix, G. J., and C. G. Lai (1999)
%
%Authors: Morteza Zarrabi and Shahram Pezeshk
%Department of Civil Engineering 
%The University of Memphis
%Spring 2004 
%REVISED: April 2005 by Morteza Zarrabi and Shahram Pezeshk
%To save and insert the best solution of previous 
%generations in place of the worst solution of the next generations. 
%This has been done to keep and not lose the best solutions of each 
%generation. 
%REVISED: Summer 2005 by Morteza Zarrabi to calculate and save the best 
%solutions of the Shear Wave Velocity
%
%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 2 of the License, 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.
%
%You should have received a copy of the GNU General Public License
%along with this program; if not, write to the Free Software
%Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
%
%This program is associated with the Following publications: 
%1- Pezeshk, S., and M. Zarrabi (2005). A new inversion procedure for 
%spectral analysis of surface waves using agenetic algorithm, Bull. 
%Seism. Soc. Am. 95(5), 1801-1808. 
%2- Morteza Zarrabi (2006). Ph.D. Dissertation, the University of Memphis.
%For more information please visit http://umdrive.memphis.edu/mzarrabi/www
%
%--------------------------------------------------------------------------
clc;
clear all;
StartTime = clock

MaxGen  = 70;  % Maximum number of generations
popSize = 70;  % Number of populations in each generation

geneLength = [15,15,15 15]; % 4 Layers including half space

Max_Layer = 4;

% Parameter Max_Con is used to control convergence. The smaller the number, 
% the better solution will be obtained; however, the number of populations 
% will be much higher and will take longer to converge.
% 
Max_Con = 45; % Maximum value of best objective function



%********************* Assign lower and upper bounds **********************

for i = 1:Max_Layer     
    lowBound(i) = 400;    % m/s
    
    uppBound(i) = 600.0;  % m/s 
end

% ********************** Develope FieldD matrix ***************************

code     = zeros(1,length(geneLength)) ;  % code for decoding
scale    = zeros(1,length(geneLength)) ;  % scaling for decoding
lowBin   = ones(1,length(geneLength))  ;  % 
uppBin   = ones(1,length(geneLength))  ;

FieldD = [geneLength;lowBound;uppBound;code;scale;lowBin;uppBin] ;


%***************** Read experimental dispersion curve *********************
 
           load NumExample1.txt;
           vr_exp = NumExample1;   

%********** Generate the initial binary string population *****************

chromosome = crtbp(popSize,sum(geneLength));

%****************** Loop over the number of generations *******************
genNum  = 1;   %Initial number of generations

while genNum <= MaxGen 
    Generation = genNum
    
%*************** Decode chromosomes to unknown real values ****************
%
%*** Line('state',sum(100*clock)) added to the end of crtbp to link the ***
%*** initial random number generation to computer clock.                *** 

    unknowns = bs2rv(chromosome,FieldD); 
      
%*********************** Estimate shear wave velocities *******************
%
    for j = 1:popSize
        Generation = genNum
        Population = j
        for i = 1:Max_Layer     
           vs(i) =    unknowns(j,i);
           vsPop(j,i)= unknowns(j,i);    
        
        end
   
         velocity = Mlayer_NumEx1(vs); 
       
         vr_theo(:,j) = velocity; 
 
    end
                    
%************ Compute the Fitness of the Population and sort **************

   [i1,i2] = size(vr_exp);
   for j = 1:popSize 
       ObjV1(j) = 0.0;
       for i = 1:i1       
           ObjV1(j) = ObjV1(j) + (abs(vr_theo(i,j) - vr_exp(i)))^1;
       end
   end 
       
%********** Trasnpose the objective vector into a column vector************
    
    ObjV = ObjV1';
    
    
%******* Find the maximum objective function which is the worst case ******   
     [bad_pop,iworst] = max(ObjV);
 
    
%********* Change the maximum objective function with the best one   ******
%********* of the previous generation for generations greater than 1 ******
%     
      if genNum > 1 
          chromosome(iworst,:) = old_chromosome(best_index(genNum - 1),:); 
          ObjV(iworst) = ObjV_old(ibest);
          vr_theo(:,iworst) = old_vr_theo(:,ibest);
      end

%************** Replace the corresponding objective function **************
    
       ObjV_old = ObjV;
    
    
% ***** find the minimums of objective function for each generation *******
% *****(Best), and its mean                                         *******

      [Best(genNum,1),ibest] = min(ObjV);
      Best(genNum,2) = mean(ObjV);
     

% **********Save the best index number (ibest) in a vector*****************
%      
     best_index(genNum) = ibest;
%      
     
% *********** Save all the generated Shear Wave Velocities  ***************

     vsAll(:,:,genNum) = unknowns(:,:);
% 
% ******** save the best Shear Wave Velocity set for each generation ******
%        
     vsGen_best(:,genNum) = unknowns(ibest,:)'; 
%     
%**** Save the worst solution number ibest in a vector(in A_GA9 043005)****
%      
     worst_index(genNum) = iworst;
%     
%*** Save the best dispersion curve of each solution for each generation **
%    
     best_vr_theo(:,genNum) = vr_theo(:,ibest);
%  
%** Save the worst dispersion curve of each solution for each generation **
%   
     worst_vr_theo(:,genNum) = vr_theo(:,iworst); 
% 
%* Parameter vr_theo corresponding to the Best (ibest) for each generation*   
%    
     vr_theo_best(:,genNum)= vr_theo(:,ibest); 
% 
%******** Parameter vr_theo corresponding to the LAST BEST (ibest) ********
%     
    Last_vr_theo_best= vr_theo(:,ibest);
     
%************ Save the chromosome of the previous generation***************
   
     old_chromosome = chromosome;
   
%**** Save the theoretical phase velocities of the previous generation ****   
    
     old_vr_theo = vr_theo;
   

%************************ Ranking and fittness ****************************   
  
   selectPress   = 2.0 ;      % Ranking: 1.0 < selectPress < 2.0
   rankingMethod = 0    ;      % Ranking method: 0 = linear; 1 = nonlinear
 
   FitnV = ranking(ObjV,[selectPress,rankingMethod]);
    
%********** Perform selection from populations for Recombination **********
    
    SEL_F = 'sus';   
    SelCh = select(SEL_F,chromosome,FitnV);
    
%************************* Execute Recombination **************************
    
    XOV_F = 'xovdp';     
    xovRate = 0.6 ;       
    SelCh = recombin(XOV_F, SelCh, xovRate);
    
%************************** Mutate offspring ******************************
    
    MUT_F   = 'mut';          
    mutRate = 0.01 ;         
    
    SelCh   = mutate(MUT_F, SelCh, [], mutRate);

%************ Insert offspring in population replacing parents ************
    
    reinsOpt(1) = 1   ;      
    reinsOpt(2) = 1.0 ;      
   
    chromosome  = reins(chromosome, SelCh, 1, reinsOpt, ObjV);
     

    
%******* Terminate the loop if Best reaches to the limit (Max_Con) ********   
      if Best(genNum,1)<Max_Con
          genNum
          genNum = MaxGen
      end   
%      
%*************** Increment the generation and go to the next **************
      
      genNum = genNum + 1;    

 end   %End the program   

EndTime = clock

 figure(1)
 plot(Best(:,1));
 xlabel('Number of Generations')
 ylabel('Error Function')
 title('Convergence History for Example 1 (without thickness estimation)')
 %
 figure(2)
 freq = linspace(5,100,50);
 plot(freq,vr_exp,'xr');
 hold on
 plot(freq,vr_theo(:,ibest),'k');
 xlabel('Frequency (Hz)')
 ylabel('Phase Velocity (m/s)')
 title('Dispersion Curve for Example 1 (without thickness estimation)')
 legend ('Experimental Phase Velocity', 'Theoretical Phase Velocity') 
%

 figure(3)
 thk = [5 10 10]; %The simulated thicknesses
 vs_sim = [500 400 500 600]; %The simulated shear wave velocities
 vs_est = vsGen_best(:,length(vsGen_best));
 X = sum(thk)+17;
 thk_cum=[X, sum(thk), sum(thk), sum(thk(1)+ thk(2)), sum(thk(1)+ ...
          thk(2)),thk(1),thk(1),0];
 vs_sim_r = [vs_sim(4), vs_sim(4), vs_sim(3), vs_sim(3), vs_sim(2), ... 
             vs_sim(2), vs_sim(1), vs_sim(1)];
 vs_est_r = [vs_est(4), vs_est(4), vs_est(3), vs_est(3), vs_est(2), ...
             vs_est(2), vs_est(1), vs_est(1)];
 plot(vs_sim_r, thk_cum, 'k');
 hold on
 plot(vs_est_r, thk_cum, '--r');
 axis([0, vs_sim(4)+200, 0, X+3])
 axis ij
 xlabel('Shear Wave Velocity (m/s)')
 ylabel('Depth (m)')
 title('Soil Profile for Example 1 (without thickness estimation)')
 legend ('Simulated Shear Wave Velocity Profile', ...
         'Estimated Shear Wave Velocity Profile',3)