**4. Optimization (Python vs C++)**
=========================================================

4.1. Differential Evolution Algorithm
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
In this section we implement a powerful optimization algorithm which is called **Differential Evolution** (Storn & Price, 1997). This algorithm is one of the most popular population-based algorithms that is very efficient in optimization of complex and multi-dimensional problems. We explain how it works and then implement it in Python:

Lets say we have a 2-dimensional problem meaning that there are 2 parameters that we want to find in a way that the cost function is minimized. First, we need to set up a list of population. The population list include vectors :math:`\vec{x_1} , ... , \vec{x_n}`. Each vector has two components corresponding to the unknowns of interest. The number of the population should be set by the user:


.. math:: 
   :name: eq.93 

   Populations:\begin{cases}
     \vec{x_1}=(a_1,b_1) \\
     . \\
     . \\
     . \\
     \vec{x_n} = (a_n,b_n)
   \end{cases}

Then we need to define a trial vector :math:`\vec{u}` for each member in the population list. The components of the trial vector are determined based on the below equation: 

.. math:: 
   :name: eq.94

   u_{nj}=x_{kj}+F \times (x_{lj} - x_{mj}) \quad    n \neq k \neq l \neq m 

.. note::

   In the above equation the :math:`F` is called **Mutation Factor** that should be in the range of [0,2]

It should be noted that we need to set a limit for each parameter for the optimization. In the trial vector, if one of the components corresponding to a parameter is higher than the maximum limit for that parameter, it will be pushed to the maximum value of the limit. Similarly, if it is less than the minimum of the limit it should be pushed back to the minimum value of the limit that has been set for that parameter. 

.. math:: 
   :name: eq.95

   check:\begin{cases}
     If \quad u>upper\,bound \Rightarrow u=upper\,bound \\
     If \quad u<lower\,bound \Rightarrow u=ower\,bound \\
   \end{cases}

Now, we check in the trial vector to know if we should keep that particular component or replace it with the corresponding component from the member in the population list that the trial vector has been produced based on. In this regard, we should define a cross-over value (:math:`C.R` )which is a number we should set between 0 and 1. A random number should be produced 2 times to determine each component of the trial vector. If the :math:`C.R` is larger than the produced random number, then the component is kept inside the trial vector and if it is less than the produced random number it should be replaced with the component of the corresponding member in the population list: 

.. math:: 
   :name: eq.96

    X_{nj}=\begin{cases}
    u_{nj}  \quad if \quad urand < C.R \\
     x_{nj}  \quad if \quad urand  > C.R
   \end{cases}

In the above, the :math:`\vec{X}` is called **Child** and corresponding member in the population list is called **Parent**. 
Now we should decide if we should keep the child or the parent. This should be done by replacing the component of the **Child Vector** and **Parent Vector** in the function that we want to minimize and evaluate the function. We should choose the one resulting in the smaller value of the function. Based on that, only one of the **Child** or **Parent** will survive in the population list. 

  

.. math:: 
   :name: eq.97

     Population-Update:\begin{cases}
     If \quad f(X_1)<f(x_1) \Rightarrow  Replace \quad x_1 \quad with \quad X_1 \quad in \quad Population \\
     Elif \quad f(X_1)>f(x_1) \quad keep \quad x_1 \quad in \quad Population \\
     . \\
     . \\
     . \\
          If \quad f(X_n)<f(x_n) \Rightarrow  Replace \quad x_n \quad with \quad X_n \quad in \quad Population \\
     Elif \quad f(X_n)>f(x_n) \Rightarrow Keep \quad x_n \quad in \quad Population 
   \end{cases}


After the above step, the population list has been updated. 
Finally, we should compare all members in the population list to see which one leads to a smaller value of the function that we want to minimize. Then we should select that particular member as the best vector containing the best optimum parameters. 


.. math:: 
   :name: eq.98

    \textbf f = [f_1,...,f_n] \Rightarrow Select \quad \text{min}[\textbf f]

Now we are done with the first generation. The same procedure should be repeated based on the updated population list. 


4.2. Python Implementation
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
The python code is implemented for a simple function :math:`f = (parameter.1)^3+ (parameter.2)^3` with the limit of :math:`[-3,3]` for both parameters to find the minimum of the function :math:`f`

.. code-block:: python

        import random
        import numpy as np

        def urand():
            return random.random()


        def func(X):

            #print ("aaa")
            return pow(X[0],3)+pow(X[1],3)

        up_bound = [3,3]
        low_bound = [-3,-3]

        dim= 2
        NP = 10
        cr = 0.90  # crossover probability
        F5 = 0.50  # Scaling factor

        Iteration_Number=20

        parent_start = []
        parent_val_start = []

        def start():

            for j in range(NP):
                x = []
                for k in range (dim):
                    x.append(random.uniform(low_bound[k],up_bound[k]))
                parent_start.append(x)
                parent_val_start.append(func(x))

        start()


        all_ind = []
        for j in range (len(parent_start)):
            all_ind.append(j)

        def DE():

            generation = 0
            iteration = []
            param1 = []
            param2 = []
            cost_func = []

            for z in range(0, Iteration_Number):

                generation=generation+1

                for i in range (len(parent_start)):

                    trials = []

                    p1 = random.choice([e for k,e in zip(all_ind,all_ind) if k != i])
                    p2 = random.choice([e for k,e in zip(all_ind,all_ind) if k != i and k != p1])
                    p_target = random.choice([e for k,e in zip(all_ind,all_ind) if k != i and k != p1 and k != p2])

                    for k in range(0, dim):
                        if urand() <= cr:
                            trials.append((parent_start[p_target])[k] + F5 * ((parent_start[p1])[k] - (parent_start[p2])[k]))
                        else:
                            trials.append((parent_start[i])[k])

                    for s in range(dim):

                        if trials[s] > up_bound[s]:
                            trials[s]=up_bound[s]
                        elif trials[s] < low_bound[s]:
                            trials[s]=low_bound[s]
                        elif low_bound[s] <= trials[s] <= up_bound[s]:
                            trials[s]=trials[s]

                    child_val = func(trials)

                    if child_val <= parent_val_start[i]:
                        parent_start[i] = trials
                        parent_val_start[i] = child_val

                soretd_val = sorted(parent_val_start,reverse=False)
                ind = parent_val_start.index(soretd_val[0])

                iteration.append(generation)
                param1.append(parent_start[ind][0])
                param2.append(parent_start[ind][1])
                cost_func.append(soretd_val[0])

            col_format = "{:<30}" * 4 + "\n"
            with open("RESULTS-Python.txt", 'w') as of:
                for x in zip(iteration, param1, param2,cost_func):
                    of.write(col_format.format(*x))

        DE()


The results are saved in a text file (RESULTS-Python.txt). In this file the columns 1-4 corresponds to the number of generation, best value found for the first parameter, best value found for the second parameter and the minimum value of the objective function respectively.



.. code-block:: python

        1                             -3                            0.13575177113689574           -26.997498292598483           
        2                             -3                            -0.32372692704626427          -27.033926298141374           
        3                             -2.8804034576214943           -3                            -50.897912723155216           
        4                             -2.8914790672098825           -3                            -51.17464792180381            
        5                             -2.91210613749475             -3                            -51.69571468276654            
        6                             -2.91210613749475             -3                            -51.69571468276654            
        7                             -2.923664426142372            -3                            -51.99093876783297            
        8                             -2.9752511297830933           -3.0                          -53.33727790449086            
        9                             -2.9752511297830933           -3.0                          -53.33727790449086            
        10                            -2.9778505435169755           -3.0                          -53.40636919427237            
        11                            -2.9792931692582787           -3.0                          -53.444765647001205           
        12                            -2.9832801955402197           -3.0                          -53.55107657228403            
        13                            -2.9832801955402197           -3.0                          -53.55107657228403            
        14                            -2.9909008960491836           -3.0                          -53.75506858321391            
        15                            -2.9979827647658714           -3.0                          -53.945571263611825           
        16                            -2.9979827647658714           -3.0                          -53.945571263611825           
        17                            -2.9979827647658714           -3.0                          -53.945571263611825           
        18                            -2.9979827647658714           -3.0                          -53.945571263611825           
        19                            -3                            -3.0                          -54.0                         
        20                            -3                            -3.0                          -54.0                    


4.3. C++ Implementation
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. code-block:: c
 
        #include <random>
        #include <vector>
        #include <fstream>
        #include <iostream>

        using namespace std;


        float func(float x, float y){

            return pow(x,3)+pow(y,3);
        }

        vector<float> up_bound {3,3};
        vector<float> low_bound {-3,-3};


        int dim= 2;
        int NP = 10;
        float cr = 0.90;  // crossover probability
        float F5 = 0.50;  // Scaling factor

        int Iteration_Number=20 ;


        vector<vector<float>> parent_start;
        vector<float> parent_val_start;


        random_device rd;
        default_random_engine eng(rd());


        void start() {

            for (int j = 0; j < NP; j++ ){
                
                vector< float > x;

                for (int k = 0; k < dim; k++ ){

                    uniform_real_distribution<> distr{low_bound[k], up_bound[k]};

                    x.push_back(distr(eng));
                }
                parent_start.push_back(x);
                parent_val_start.push_back(func(x[0],x[1]));
            }
        };

        int find_index (int size, float target, vector<float> val_array ){

            int iterator;

            for (int q = 0 ; q < size; q++){

                if (val_array[q] == target){

                    iterator = q;

                    }
                }
                return iterator;
            };

        int generation = 0;
        vector<int> iteration;
        vector<float> param1;
        vector<float> param2;
        vector<float> cost_func;


        uniform_real_distribution<> distr_cr {0, 1};

        int main(){
        srand(unsigned(time(NULL)));

            start();


            for (int z = 0; z < Iteration_Number; z++){

                generation = generation+1;

                        for (int i = 0; i < parent_start.size(); i++){

                            vector<float> trials;

                            int p1;
                            int p2;
                            int p_target;

                            do {

                                p1 = rand() % parent_start.size();
                                p2 = rand() % parent_start.size();
                                p_target = rand() % parent_start.size();
                                
                            //cout << "p1 :" << p1 << endl;

                            } while (p1 == i || p1 == p2 || p1 == p_target || p2 == p_target);


                                for (int k = 0; k < dim; k++){

                                    float rand_0_1 = distr_cr (eng);

                                    if (rand_0_1 <= cr){

                                        trials.push_back((parent_start[p_target])[k] + F5 * ((parent_start[p1])[k] - (parent_start[p2])[k]));
                                    }

                                    else {

                                        trials.push_back((parent_start[i])[k]);
                                    }
                                };

                                for (int s = 0; s < dim; s++){

                                    if (trials[s] > up_bound[s]){

                                        trials[s]=up_bound[s];
                                    }

                                    else if (trials[s] < low_bound[s]){

                                        trials[s]=low_bound[s];
                                    }

                                    else if (low_bound[s] <= trials[s] <= up_bound[s]){

                                        trials[s]=trials[s];
                                    }
                                };

                                float child_val = func(trials[0],trials[1]);


                                if (child_val <= parent_val_start[i]){

                                    parent_start[i] = trials;
                                    parent_val_start[i] = child_val;
                                }
                        };


                        vector<float> soretd_val = parent_val_start;
                        sort(soretd_val.begin(), soretd_val.end());

                        int ind = find_index (parent_val_start.size(),soretd_val[0],parent_val_start);

                        iteration.push_back(generation);
                        param1.push_back(parent_start[ind][0]);
                        param2.push_back(parent_start[ind][1]);
                        cost_func.push_back(soretd_val[0]);
            }

            string filename("Results-CPP.txt");
            fstream file_out;
            file_out.open(filename, std::ios::out);

            for (int r=0 ; r < iteration.size(); r++){

                file_out << iteration[r] << "       " << param1[r]<< "       " << param2[r]<< "       " << cost_func[r]<< endl;     
        
            };

            file_out.close();

            return 0;
        };

The results are saved in a text file (RESULTS-CPP.txt). In this file the columns 1-4 corresponds to the number of generation, best value found for the first parameter, best value found for the second parameter and the minimum value of the objective function respectively.

.. code-block:: c

        1       1.18459         -3            -25.3377
        2       -1.54879       -2.90724       -28.2871
        3       -1.83923       -2.95362       -31.9886
        4       -1.83923       -2.95362       -31.9886
        5       -2.30088       -2.9091        -36.8003
        6       -2.56735       -2.93374       -42.1723
        7       -2.56735       -2.93374       -42.1723
        8       -2.56735       -2.93374       -42.1723
        9       -2.6831        -2.90672       -43.8746
        10       -2.73581       -2.88452       -44.4773
        11       -2.81982       -2.93374       -47.6718
        12       -2.90422       -2.90321       -48.9656
        13       -2.91029       -2.93327       -49.8878
        14       -2.91029       -2.96987       -50.8442
        15       -3             -2.96256       -53.0016
        16       -3             -2.96256       -53.0016
        17       -3             -3              -54
        18       -3             -3              -54
        19       -3             -3              -54
        20       -3             -3              -54

Like the previous Implementation in Python, the C++ code minimizes the objective function to -54 where the best value for each parameter was found equal to -3.  

.. note::

   The limits of the parameter that we want to find as well as the number of populations, :math:`C.R` , :math:`F` and number of iterations should be set in way a that it leads to best possible results. It usually needs some try and error and good understanding of the physics of the problem.