
% ------------------------------------------------------------------------
% Script: genetic_algorithm_2d
% ------------------------------------------------------------------------
% Aim:
% To run the GA to determine the dislocation core structure
% ------------------------------------------------------------------------

global Kx Kz;
global Ax Dx Wx dx;
global Az Dz Wz dz;

%   edge or screw dislocation

    % "theta" is the angle between the dislocation line and its Burgers 
    %         vector for a general mixed dislocation in an isotropic solid

    if mistype=='e'   % for edge dislocation
        theta=pi/2;
    end
    
    if mistype=='s'   % for screw dislocation
        theta=0;
    end
    
    Kx=shear_modulus*(sin(theta)^2/(1-poisson_ratio)+cos(theta)^2);
    Kz=shear_modulus*(sin(theta+pi/2)^2/(1-poisson_ratio)+cos(theta+pi/2)^2);
    
% Lower bound (lb) and Upper bound (ub)
    
    if Nmis>=1

        % % To define some matrix

        lb=zeros(6*Nmis-4,1);
        ub=zeros(6*Nmis-4,1);
        InitialPopulationMatrix=zeros(PopulationSize,6*Nmis-4);

        % % For "x" variables
        for i=1:Nmis-1
            lb(i,1)=Xcoef_range(1);
            ub(i,1)=Xcoef_range(2);
            InitialPopulationMatrix(:,i)=(rand(PopulationSize,1)*0.2+0.9)*Xcoef_range(3);
        end
        for i=1:Nmis-1
            lb(Nmis-1+i,1)=Xdist_range(1);
            ub(Nmis-1+i,1)=Xdist_range(2);
            InitialPopulationMatrix(:,Nmis-1+i)=(rand(PopulationSize,1)*0.2+0.9)*Xdist_range(3);
        end
        for i=1:Nmis
            lb(2*Nmis-2+i,:)=Xwid_range(1);
            ub(2*Nmis-2+i,:)=Xwid_range(2);
            InitialPopulationMatrix(:,2*Nmis-2+i)=(rand(PopulationSize,1)*0.2+0.9)*Xwid_range(3);
        end

        % % For "z" variables
        for i=1:Nmis-1
            lb(3*Nmis-2+i,:)=Xcoef_range(1);
            ub(3*Nmis-2+i,:)=Xcoef_range(2);
            InitialPopulationMatrix(:,3*Nmis-2+i)=(rand(PopulationSize,1)*0.2+0.9)*Xcoef_range(3);
        end
        for i=1:Nmis-1
            lb(4*Nmis-3+i,:)=Xdist_range(1);
            ub(4*Nmis-3+i,:)=Xdist_range(2);
            InitialPopulationMatrix(:,4*Nmis-3+i)=(rand(PopulationSize,1)*0.2+0.9)*Xdist_range(3);
        end
        for i=1:Nmis
            lb(5*Nmis-4+i,:)=Xwid_range(1);
            ub(5*Nmis-4+i,:)=Xwid_range(2);
            InitialPopulationMatrix(:,5*Nmis-4+i)=(rand(PopulationSize,1)*0.2+0.9)*Xwid_range(3);
        end    

    elseif Nmis==0
        lb=[Xdist_range(1);Xwid_range(1);Xdist_range(1);Xwid_range(1)];
        ub=[Xdist_range(2);Xwid_range(2);Xdist_range(2);Xwid_range(2)];
        InitialPopulationMatrix=[(rand(PopulationSize,1)*0.2+0.9)*Xdist_range(3),...
                                 (rand(PopulationSize,1)*0.2+0.9)*Xwid_range(3),...
                                 (rand(PopulationSize,1)*0.2+0.9)*Xdist_range(3),...
                                 (rand(PopulationSize,1)*0.2+0.9)*Xwid_range(3)];
    elseif Nmis==-1
        lb=[Xcoef_range(1);Xdist_range(1);Xwid_range(1);dx_range(1);...
            Xcoef_range(1);Xdist_range(1);Xwid_range(1);dx_range(1)];
        ub=[Xcoef_range(2);Xdist_range(2);Xwid_range(2);dx_range(2);...
            Xcoef_range(2);Xdist_range(2);Xwid_range(2);dx_range(2)];
        InitialPopulationMatrix=[(rand(PopulationSize,1)*0.2+0.9)*Xcoef_range(3),...
                                 (rand(PopulationSize,1)*0.2+0.9)*Xdist_range(3),...
                                 (rand(PopulationSize,1)*0.2+0.9)*Xwid_range(3),...
                                 (rand(PopulationSize,1)*0.2+0.9)*dx_range(3),...
                                 (rand(PopulationSize,1)*0.2+0.9)*Xcoef_range(3),...
                                 (rand(PopulationSize,1)*0.2+0.9)*Xdist_range(3),...
                                 (rand(PopulationSize,1)*0.2+0.9)*Xwid_range(3),...
                                 (rand(PopulationSize,1)*0.2+0.9)*dx_range(3)];
    end

% The number of variables in the problem

    if Nmis>=1
        nvars=6*Nmis-4;
    elseif Nmis==0
        nvars=4;
    elseif Nmis==-1
        nvars=8;
    end

% Setting options setting

    % % Size of the population and Maximum number of iterations before the algorithm halts

    % PopulationSize=[];
    % MaxIterations=[];
    
    % % Modify options setting
    
    % options = optimoptions('ga','MaxGenerations',MaxIterations,'PopulationSize',PopulationSize,'InitialPopulationMatrix',InitialPopulationMatrix,'Display','iter','FunctionTolerance',1E-25);
    % options = optimoptions(options,'PlotFcn', {[]},'UseParallel',false);
    
    options = gaoptimset('Generations',MaxIterations,'PopulationSize',PopulationSize,'InitialPopulation',InitialPopulationMatrix,'Display','iter','TolFun',1E-25);
    options = gaoptimset(options,'PlotFcn', {[]},'UseParallel',false);
  
%% To run the GA

    fprintf('        >> The output of Genetic Algorithm is shown as follows:\n');
    
    % % Classical P-N model (mode=1)
    if pnmode==1
        ga_result=evalc('[variable_ga,fval] = ga(@minfunc,nvars,[],[],[],[],lb,ub,[],options);');
    end

    % % Semidiscrete Variational P-N model (mode=2)
    if pnmode==2
        ga_result=evalc('[variable_ga,fval] = ga(@minfunc_svpn,nvars,[],[],[],[],lb,ub,[],options);');
    end
    
    fprintf('%s',strrep(['           ',ga_result],char(10),[10,'           ']));
    
    % % To output the results of the various energies of dislocation (1e-10*J/m)\n')
    
    run output_Edis.m
    
%% To output the results of variational constants in disregistry function

    if Nmis>=1

        if Nmis==1
            Ax(Nmis)=1;
            Az(Nmis)=1;
            Dx(Nmis)=0;
            Dz(Nmis)=0;
        else 
            Ax(Nmis)=1-sum(Ax(1:(Nmis-1)));
            Az(Nmis)=0-sum(Az(1:(Nmis-1)));
            Dx(Nmis)=0-sum(Dx(1:(Nmis-1)));
            Dz(Nmis)=0-sum(Dz(1:(Nmis-1)));
        end

    end

    fprintf('\n');
    fprintf('    >> The results of variational constants in disregistry function (A) \n');
    fprintf('\n');

    if Nmis>=1

        fprintf('        >>   NO.   Xcoef   Xdist   Xwid \n');

        for i=1:Nmis
            fprintf('           %3i %10.4f %10.4f %10.4f \n',i, Ax(i), Dx(i), Wx(i));
        end

        fprintf('\n');

        fprintf('        >>   NO.   Zcoef   Zdist   Zwid \n');

        for i=1:Nmis
            fprintf('           %3i %10.4f %10.4f %10.4f \n',i, Az(i), Dz(i), Wz(i));
        end

    elseif Nmis==0

        fprintf('        >>   Xdist   Xwid   Zdist   Zwid \n');
        fprintf('           %10.4f %10.4f %10.4f %10.4f \n', Dx, Wx, Dz, Wz);

    elseif Nmis==-1

        fprintf('        >>   Xcoef   Xdist   Xwid   dx \n');
        fprintf('           %10.4f %10.4f %10.4f %10.4f \n', Ax, Dx, Wx, dx); 

        fprintf('\n');

        fprintf('        >>   Zcoef   Zdist   Zwid   dz \n');
        fprintf('           %10.4f %10.4f %10.4f %10.4f \n', Az, Dz, Wz, dz);        

    end

    fprintf('\n');
   
