hybridOptimisation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "thirty-engine",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Question 3a: GP (10%)\n",
    "\n",
    "# Packages\n",
    "import numpy as np\n",
    "import random as rand\n",
    "from copy import deepcopy\n",
    "from inspect import signature\n",
    "\n",
    "\n",
    "# Let numpy throw errors\n",
    "np.seterr(all='raise')\n",
    "\n",
    "# Constants\n",
    "b = 1        # Range of datapoints\n",
    "D = 40       # Number of datapoints\n",
    "trials = 50  # Number of times PSO is run\n",
    "XOR = 0.9    # Crossover rate\n",
    "MUR = 0.1    # Mutation rate\n",
    "PMR = 0.02   # Point mutation rate\n",
    "GEN = 40     # Maximum allowed number of generations\n",
    "MINDEPTH = 2 # Minimum initial tree depth\n",
    "MAXDEPTH = 5 # Maximum initial tree depth\n",
    "\n",
    "# Functions\n",
    "def epsilon_close(x, y):\n",
    "    return abs(x - y) < 0.01\n",
    "\n",
    "def target(x):\n",
    "    return x**4 + x**3 + x**2 + x\n",
    "\n",
    "# Data\n",
    "X = np.array([rand.uniform(-b, b) for _ in range(D)])\n",
    "Y = np.array([target(x) for x in X])\n",
    "\n",
    "# Node values\n",
    "def add(x, y): return x + y\n",
    "def mul(x, y): return x * y\n",
    "F = [add, mul]\n",
    "T = ['x']\n",
    "\n",
    "# Genetic program class\n",
    "class GP: \n",
    "    \n",
    "    # Constructor: \n",
    "    def __init__(self, v = None, l = None, r = None): \n",
    "        self.v = v        # node value\n",
    "        self.l = l        # left arg\n",
    "        self.r = r        # right arg\n",
    "        \n",
    "    # String repr. of tree root node\n",
    "    def to_str(self): \n",
    "        if (self.v in F): return self.v.__name__\n",
    "        else: return str(self.v)\n",
    "    \n",
    "    # Print whole tree\n",
    "    def tree_print(self, prefix=\"\"): # Print program tree\n",
    "        print(\"%s%s\" % (prefix, self.to_str()))        \n",
    "        if self.l: self.l.tree_print(prefix + \"  \")\n",
    "        if self.r: self.r.tree_print(prefix + \"  \")\n",
    "\n",
    "    # Run the program on input\n",
    "    def apply(self, x):\n",
    "        if self.v in F:  \n",
    "            try:\n",
    "                return self.v(self.l.apply(x), self.r.apply(x))\n",
    "            except FloatingPointError:\n",
    "                return float(\"inf\")\n",
    "        elif self.v == 'x': return x\n",
    "        else: return self.v\n",
    "    \n",
    "    # Generate random tree\n",
    "    def random_tree(self, grow, max_depth, depth = 0):\n",
    "        if depth < MINDEPTH or (depth < max_depth and not grow): \n",
    "            self.v = rand.choice(F)\n",
    "        elif depth >= max_depth:   \n",
    "            self.v = 'x'\n",
    "        else: # intermediate depth, grow\n",
    "            if rand.uniform(0, 1) > 0.5:\n",
    "                self.v = rand.choice(T)\n",
    "            else: self.v = rand.choice(F)\n",
    "        if self.v in F:\n",
    "            self.l = GP()\n",
    "            self.l.random_tree(grow, max_depth, depth = depth + 1)            \n",
    "            self.r = GP()\n",
    "            self.r.random_tree(grow, max_depth, depth = depth + 1)\n",
    "    \n",
    "    # Mutate\n",
    "    def mutate(self): \n",
    "        if rand.uniform(0, 1) < MUR:\n",
    "            self.random_tree(grow = True, max_depth = 2)\n",
    "        elif self.l: self.l.mutate()\n",
    "        elif self.r: self.r.mutate() \n",
    "            \n",
    "    # Point mutate\n",
    "    def point_mutate(self):\n",
    "        if rand.uniform(0, 1) < PMR:\n",
    "            if self.v in T: self.v = rand.choice(T)\n",
    "            else: self.v = rand.choice(F)\n",
    "        if self.l: self.l.point_mutate()\n",
    "        if self.r: self.r.point_mutate()\n",
    "            \n",
    "     # Tree depth\n",
    "    def depth(self, curr_depth=0):\n",
    "        below = curr_depth\n",
    "        if self.l: below = max(below, self.l.depth(curr_depth + 1))\n",
    "        if self.r: below = max(below, self.r.depth(curr_depth + 1))\n",
    "        return below\n",
    "    \n",
    "    # Tree size in nodes\n",
    "    def size(self): \n",
    "        if self.v in T: return 1\n",
    "        l = self.l.size() if self.l else 0\n",
    "        r = self.r.size() if self.r else 0\n",
    "        return 1 + l + r\n",
    "\n",
    "    # Builds an identical tree\n",
    "    def build_subtree(self): \n",
    "        t = GP()\n",
    "        t.v = self.v\n",
    "        if self.l: t.l = self.l.build_subtree()\n",
    "        if self.r: t.r = self.r.build_subtree()\n",
    "        return t\n",
    "    \n",
    "    # Scan tree\n",
    "    def scan_tree(self, count, second): \n",
    "        count[0] -= 1 # note: count is list, so it's passed \"by reference\"\n",
    "        if count[0] <= 1: \n",
    "            if not second: # return subtree rooted here\n",
    "                return self.build_subtree()\n",
    "            else: # glue subtree here\n",
    "                self.v = second.v\n",
    "                self.l = second.l\n",
    "                self.r = second.r\n",
    "        else:  \n",
    "            ret = None              \n",
    "            if self.l and count[0] > 1: ret = self.l.scan_tree(count, second)  \n",
    "            if self.r and count[0] > 1: ret = self.r.scan_tree(count, second)  \n",
    "            return \n",
    "        \n",
    "    # Crossover two trees at random\n",
    "    def crossover(self, other):\n",
    "        if rand.random() < XOR:\n",
    "            second = other.scan_tree([rand.randint(1, other.size())], None) # 2nd random subtree\n",
    "            self.scan_tree([rand.randint(1, self.size())], second) # 2nd subtree \"glued\" inside 1st tree\n",
    "\n",
    "# Fitness\n",
    "def fitness(gp):\n",
    "    return 1 / (1 + np.sum([abs(gp.apply(X[i]) - Y[i]) for i in range(D)])) \n",
    "    \n",
    "# Selection\n",
    "def max_select(gps): \n",
    "    fitnesses = np.array([fitness(gp) for gp in gps])\n",
    "    return deepcopy(gps[np.argmax(fitnesses)])\n",
    "\n",
    "# Roulette-wheel selection\n",
    "def roulette_select(gps, fitnesses):\n",
    "    sum_fitness = np.sum(fitnesses)\n",
    "    if sum_fitness == 0.0:\n",
    "        return deepcopy(rand.choice(gps))\n",
    "    else:\n",
    "        norm_fitnesses = fitnesses / sum_fitness\n",
    "        choice = rand.choices(gps, weights=norm_fitnesses, k=1)[0]\n",
    "        return deepcopy(choice)\n",
    "    \n",
    "# Generation initial population\n",
    "def init_population(POP): \n",
    "    pop = []\n",
    "    for md in range(3, MAXDEPTH+1):\n",
    "        for i in range(int(POP/6)):\n",
    "            t = GP()\n",
    "            t.random_tree(grow = True, max_depth = md) # Grow method\n",
    "            pop.append(t)\n",
    "        for i in range(int(POP/6)):\n",
    "            t = GP()\n",
    "            t.random_tree(grow = False, max_depth = md) # Full method\n",
    "            pop.append(t)\n",
    "    return pop\n",
    "\n",
    "\n",
    "\n",
    "# Records\n",
    "res = 15 # Resolution\n",
    "POPs = list(range(60, 601, (600 - 60) // res))\n",
    "F_g = np.zeros((len(POPs), trials))       # Fitnesses per population size and trial\n",
    "G_c = np.zeros((len(POPs), trials))       # Generation until convergence per population size and trial\n",
    "\n",
    "for k, POP in enumerate(POPs):\n",
    "    for trial in range(trials):\n",
    "        # 1. Generate an initial population of GPs\n",
    "        GPs   = init_population(POP) # Programs\n",
    "        F_p   = np.zeros(POP)     # Fitnesses\n",
    "        f_g   = 0.0               # Trial-best fitness\n",
    "        gp_g  = None              # Trial-best program \n",
    "        converged = False         # Did this run converge?\n",
    "        \n",
    "        # 2. Run \n",
    "        for gen in range(GEN):\n",
    "            # Evaluate fitnesses, replace GB\n",
    "            for i in range(POP):\n",
    "                f_i = fitness(GPs[i])\n",
    "                F_p[i] = f_i\n",
    "                if f_i > f_g:   # Global best\n",
    "                    gp_g = deepcopy(GPs[i])\n",
    "                    f_g = f_i\n",
    "\n",
    "            # Check for termination criterion\n",
    "            if epsilon_close(f_g, 1.0): \n",
    "                converged = True\n",
    "                G_c[k][trial] += gen\n",
    "                break\n",
    "\n",
    "            # Create new population\n",
    "            GP_new = []\n",
    "            for i in range(POP):\n",
    "                gp_new = deepcopy(GPs[i])\n",
    "                mate   = roulette_select(GPs, F_p)\n",
    "                gp_new.crossover(mate)    # Crossover\n",
    "                gp_new.mutate()           # Mutation\n",
    "                gp_new.point_mutate()     # Point-mutation\n",
    "                GP_new.append(gp_new)     # Add to new population\n",
    "\n",
    "            GPs = GP_new # Replace old population\n",
    "        \n",
    "        # Not solved, took time.\n",
    "        if not converged:\n",
    "            G_c[k][trial] += GEN\n",
    "        \n",
    "        # Record trial-best fitness\n",
    "        F_g[k][trial] = f_g\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "second-latex",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAeG0lEQVR4nO3debwcVZ338c83gbDIEiVBMZAEFMHgCANXFnEwgrKNY5gREYwbohkUnBj1kSiOxpnhUVAfRmUzKIIjmmETkEHQRwU3thv2CMEYAoksSVSWEAwk+c0fdS5UOn3vbXKr+lZ3fd+vV79ud53q0+fUvbe+Xae6TisiMDOz+hox3A0wM7Ph5SAwM6s5B4GZWc05CMzMas5BYGZWcw4CM7OacxBYV5K0maQfSXpc0sWSpkr6yXC3y6yKHASGpEWS3tyw7P2Sft2w7GhJN0l6StLSdP8jkpTKz5f0jKQVkv4s6aeSdk1lsyQ9m8oek/RbSfulslGSvippSSq/X9LpTdpzl6SVkh6RdLak0QN060jgpcA2EfGOiLgwIg7O1ReSXjmU7VYmSddJ+mvaHiskzW8oP0jSvWl7/ELShDa0abKkJQXUU+ltX0cOAmuJpE8AXwO+DLyMbCd7PLA/MCq36mkRsQWwPbAUOD9X9t+pbCzwa+CyFCKfBnqAvYEtgTcBtzW89qnA/wG2BvYFJgA/lZR/7bwJwH0RsXrDez3sToyILdJtl76FksYAlwH/CrwE6AX+u5UKlRmW/3tJGw3H61oLIsK3mt+ARcCbG5a9H/h1ur818BTw9kHqOR/4j9zjvwdWpPuzgO/lynYDAhgDXAV8rJ86twJWAEc1LN+CLGg+0OQ5XwCeAZ5Nzz2uoT+/TK/9VCp/JzAZWAJ8ItX7MHBsrs5NgK8ADwKPAucAm6Wyvj48BvwZ+BUwIpWdBPwReBKYDxzU4u/kOuCD/ZRNA36be/wi4Glg1wHqOgX4TVrvlcCuwE9Te+fnty9wOPC71OY/Ap/MvcbatM1WAC8nC+8bUt8fBs4ARuXqCuAE4PfA/f1s+363n2/tufmIwFqxH9mO8IpWnyBpC2AquXf2ubJNyHbMSyJiOXAj8PE0zPQ3fUNNyeuBTcneAT8nIlYAPwbe0lh/RHwe+L+kI5CI+HZD+QHp7u6pvO/d9MvIQm8cWXicKenFqexU4FXAHmQ70nHA51LZJ8hCZCzZkdJngJC0C3Ai8LqI2BI4hCx0W/VFScsl/UbS5Nzy3YA7cv15CvhDWt6f95AFyJbAMrIQ+D6wLXAMcJakvud/G/jn1ObXAD9Pr3EY8FA8f5TyELAGmEG2M98POAj4SMNrHwHsA0zqZ9s33X6Dbx4rioPA+lyexu4fk/QYcFaubAywPHLDLGmM/zFJT0s6ILfuJ9PzF5C9a39/ruyoVLYY2ItsBwHwRbId7VSyYY4/Snpff6+d83AqL8qzwL9FxLMRcTXZO9ZdUjB9CJgREX+OiCfJgubo3PO2Ayak5/4qsrfDa8gCdJKkjSNiUUT8ocW2nATsRBY4s4EfSXpFKtsCeLxh/cfJdvL9OT8i5qXteCiwKCK+ExGrI+JW4FKy8yp9/ZkkaauI+Esqbyoi5kbEjameRcA3gTc2rPbFtN2e7qea/raftYmDwPocERGj+26s+67uT8CY/BhvRLw+rfcn1v07+kqq42UR8baGHd9FqWzbiDgwIuamutZExJkRsT8wmmwY4zxJrwaWN752znapvCh/agiclWQ73bHA5sDcXFBek5ZDdt5kAfATSQslzUz9WgB8jGxYbKmkOZJe3kpDIuKmiHgyIlZFxAVkwzqHp+IVZENmeVuRDeX0Z3Hu/gRgn4bgn0p2RATw9vRaD0i6vu+kfjOSXiXpqnQC/wmygGwM58VNnprXdPtZ+zgIrBU3AKuAKWW/UEQ8HRFnAn8BJuVe+5/y60l6EdlQxc/KbhNZ2DwN7JYLy60jO/FN2mF/IiJ2Av6BbJjroFT2/Yh4A9nON8iOfDZEAH1DZvOA3fsK0rZ4RVo+0PP7LAauzwd/Gqb5cGrzLRExhWzY6HLgoiZ19DkbuBfYOSK2IhvWUcM6A767H2j7WXs4CGxQEfEY2QnYsyQdKWkLSSMk7UF2EnFIJH0sfTRxM0kbpWGhLYHbIuLx9NrfkHSopI0lTQQuJhtX/q8NfNlHyYZeBhURa4FzgdMlbZvaPE7SIen+WyW9Mg0hPUE2JLRG0i6SDkznRP5KFiZrBns9SaMlHSJp07Q9pgIHANemVX4IvEbS2yVtSnau4s6IuLfFvl8FvErSe9L23FjS6yS9On2Ud6qkrSPi2Vx/INtm20jaOlfXlmmdFco+KvzhFl5/nW3f3/ZrsS9WAAeBtSQiTgM+DnyK7FM1j5KNB58E/HaI1T8NfBV4hOzd9wlkn1BamHvtz5B9aucJ4Cayd7UHRcSqDXzNWcAFaWjkqBbWP4ls+OLGNATy/4G+j3TunB6vIDuCOSsiriM7P/Cl1KdHyN5hfwYg7Wz7ewe/MfAfZCd1lwMfJRu6mw8QEcvIhm9OITty2ofnz1cMKp3jODg956HUtlNTeyE7sbwo9fN44N3pefcCPwAWpu32crJPFL2LbFjqXFr7GOss1t32/W0/axP5nIyZWb35iMDMrOYcBGZmNecgMDOrOQeBmVnNddwkUGPGjImJEycOdzPMzDrK3Llzl0fE2GZlHRcEEydOpLe3d7ibYWbWUSQ90F+Zh4bMzGrOQWBmVnMOAjOzmnMQmJnVnIPAzKzmHARmZjXnIDAzqzkHgZlZzTkIzMxqzkFgZlZzDgIzs5pzEJiZ1ZyDwMys5hwEZmY15yBoQUQM+Lju9VSxTVWrp4ptct/aV0/RdRWttCCQdJ6kpZLu7qdckr4uaYGkOyXtWVZbhmLWrFnMmDHjuV9aRDBjxgxmzZrleirapqrVU8U2uW/tq6foukoREaXcgAOAPYG7+yk/HPgxIGBf4KZW6t1rr72iXdauXRvTp08PIKZPn970cZ3rqWKbqlZPFdvkvrWvnqLrGgqgN/rbX/dXUMQNmDhAEHwTOCb3eD6w3WB1tjMIItb9JfbdNuSX1631VLFNVaunim1y39pXT9F1baiBgkBZeTkkTQSuiojXNCm7CvhSRPw6Pf4ZcFJErPc9lJKmAdMAxo8fv9cDD/T7jWuliAhGjHh+FG3t2rVIcj0VblPV6qlim9y39tVTdF0bQtLciOhpVjacJ4ubbYGmqRQRsyOiJyJ6xo5t+t3LpYnIxvLy8mN9da+nim2qWj1VbJP71r56iq6rFP0dKhRxo8OHhqo23li1eqrYpqrVU8U2uW/tq6fouoaCAYaGNiopX1pxJXCipDnAPsDjEfHwMLZnPZIYPXo006dP5/TTT0cSp59+OgCjR49u+bCuW+upYpuqVk8V2+S+ta+eousqS2nnCCT9AJgMjAEeBT4PbAwQEeco6/0ZwKHASuDYaHJ+oFFPT0/09g66WqEiYp1fVuPjutdTxTZVrZ4qtsl9a189Rde1IQY6R1DaEUFEHDNIeQAnlPX6RWr8ZW3oL69b66lim6pWTxXb5L61r56i6yqaryw2M6s5B4GZWc05CMzMas5BYGZWcw4CM7OacxCYmdWcg8DMrOYcBGZmNecgMDOrOQeBmVnNOQjMzGrOQWBmVnMOAjOzmnMQmJnVnIPAzKzmHARmZjXnIDAzqzkHgZlZzTkIzMxqzkFgZlZzDgIzs5pzEJiZ1ZyDwMys5hwEZmY15yAwM6s5B4GZWc05CMzMas5BYGZWcw4CM7OacxCYmdWcg8DMrOZKDQJJh0qaL2mBpJlNyreW9CNJd0iaJ+nYMttjZmbrKy0IJI0EzgQOAyYBx0ia1LDaCcDvImJ3YDLwVUmjymqTmZmtr8wjgr2BBRGxMCKeAeYAUxrWCWBLSQK2AP4MrC6xTWZm1qDMIBgHLM49XpKW5Z0BvBp4CLgLmB4RaxsrkjRNUq+k3mXLlpXVXjOzWiozCNRkWTQ8PgS4HXg5sAdwhqSt1ntSxOyI6ImInrFjxxbdTjOztpg8ObtVTZlBsATYIfd4e7J3/nnHApdFZgFwP7BriW0yM7MGZQbBLcDOknZMJ4CPBq5sWOdB4CAASS8FdgEWltgmMzNrsFFZFUfEakknAtcCI4HzImKepONT+TnAvwPnS7qLbCjppIhYXlabzMxsfaUFAUBEXA1c3bDsnNz9h4CDy2yDmZkNzFcWm5nVnIPAzCqlqp+s6WYOArMO5J2lFclBYGZWcw4CM7OacxCYmdWcg8Csxoo819Ct5y26tV95DgIzs5pzEJiZ1ZyDwMys5hwEZmY15yAwM6s5B4GZWc05CMzMas5BYGZWcw4CM7OacxCYmdWcg8DMrOYcBGZmNecgMDOrOQeBmVnNOQjMzGrOQWBmVnMOAjPrSnX4QpmiOAjM2sg7J6siB4GZWc05CMzMas5BYGZWcw4CM7OacxCYmdXcoEEgabqkrZT5tqRbJR3cSuWSDpU0X9ICSTP7WWeypNslzZN0/QvtgJmZDU0rRwQfiIgngIOBscCxwJcGe5KkkcCZwGHAJOAYSZMa1hkNnAW8LSJ2A97xglpv1gb+yKd1u1aCQOnn4cB3IuKO3LKB7A0siIiFEfEMMAeY0rDOu4DLIuJBgIhY2lqzzcysKK0EwVxJPyELgmslbQmsbeF544DFucdL0rK8VwEvlnSdpLmS3tusIknTJPVK6l22bFkLL21mZq3aqIV1jgP2ABZGxEpJLyEbHhpMs6OGaPL6ewEHAZsBN0i6MSLuW+dJEbOB2QA9PT2NdZiZ2RC0ckSwHzA/Ih6T9G7gs8DjLTxvCbBD7vH2wENN1rkmIp6KiOXAL4HdW6jbzMwK0koQnA2slLQ78CngAeC7LTzvFmBnSTtKGgUcDVzZsM4VwN9J2kjS5sA+wD0tt97MzIaslSBYHRFBdqL3axHxNWDLwZ4UEauBE4FryXbuF0XEPEnHSzo+rXMPcA1wJ3Az8K2IuHvDumJmZhuilXMET0r6NPAesnfvI4GNW6k8Iq4Grm5Ydk7D4y8DX26tuWZmVrRWjgjeCawiu57gEbJP/njHbevx5+3NOtOgQZB2/pcCm6RFy4EfltkoMzNrn1ammPgQcAnwzbRoHHB5iW0yM7M2amVo6ARgf+AJgIj4PbBtmY0yM7P2aSUIVqUpIgCQtBHrXxhmVjk+Z2HWmlaC4HpJnwE2k/QW4GLgR+U2y8zM2qWVIJgJLAPuAv6Z7OOgny2zUWZm1j6DXkcQEWuBc9PNzMy6zKBBIGl/YBYwIa0vICJip3KbZmZm7dDKlcXfBmYAc4E15TbHzMzarZUgeDwiflx6S8zMbFi0EgS/kPRl4DKyqSYAiIhbS2uVmZm1TStBsE/62ZNbFsCBxTfHzMzaraVvKIuIhfkFknyi2ErTdxHYddcNZyvM6qOV6wguabLs4qIbYmZmw6PfIwJJuwK7AVtL+qdc0VbApmU3zNrD777NbKChoV2AtwKjgX/ILX8S+FCJbTIzszbqNwgi4grgCkn7RcQNbWyTmZm10UBDQ5+KiNOAd0k6prE8Iv6l1JaZmVlbDDQ0dBJwGvAH4C/taY6ZmbXbQEHwqKQJwLHAm9rUHjMza7OBguBs4BpgJ6A3t1xkF5T5WgIzsy4w0MnibwDfkHR2RHy4jW0yM7M2GvSCMoeAmVl3a+XKYjMzG6oLL2TOjRP5+fUjYOJEuPDC4W7Rc1qZa8jMzIbiwgth2jRetmpl9viBB2DatOz+1KnD167ERwRmZmU7+WRYuXLdZStXZssrwEFgZla2Bx98YcvbzEFgZtVR4XH0IRk//oUtbzMHgZlVw3Pj6A8wgnh+HL0bwuCUU2Dzzdddtvnm2fIKcBCYWTVUfBx9SKZOhdmzeWSTCaxFMGECzJ5diRPFUHIQSDpU0nxJCyTNHGC910laI+nIMttjZhVW8XH0IZs6laP3XcSBb1wLixZVJgSgxCCQNBI4EzgMmAQcI2lSP+udClxbVlvMukpR4+hVG4+v4jh61bZRSco8ItgbWBARCyPiGWAOMKXJeh8FLgWWltgWs+5Q1Dh60ePxRewwqzaO3s3nLBqUGQTjgMW5x0vSsudIGgf8I3DOQBVJmiapV1LvsmXLCm+oWccoahy9yPH4onaYVRtH7+ZzFg3KDAI1WRYNj/8TOCki1gxUUUTMjoieiOgZO3ZsUe17QSZPfv77fc2GTVHj6EWOxxe5w6zSOHq3n7PIKTMIlgA75B5vDzzUsE4PMEfSIuBI4CxJR5TYJrPOVtQ4epHj8d26w6ziOYuSlBkEtwA7S9pR0ijgaODK/AoRsWNETIyIicAlwEci4vIS22TW2YoaRy9yPL5bd5hVO2dRotKCICJWAyeSfRroHuCiiJgn6XhJx5f1umZdrahx9CLH47t1h1m1cxYlKnX20Yi4Gri6YVnTE8MR8f4y22LWNaZO5ehzs53RdddVox6AR447mW1XPciICeOzEOiGHWZR26jiPA21mQ1dTXaY3cpTTJiZ1ZyDwMys5hwEbVbU9Qi+rsFsADWZGqIoPkdgZt2l4l8LWUU+IjCz7lKjqSGK4iAwaxcPV7RHt17pXCIHgVk71Ggmy2HXrVc6l8hBYNYOHq5on2690rlEDgKzdvBwRfvUaGqIojgIzAZS1Li+hyvaq0rTWXcAB4FZf4oc1/dwhVWYg8CsPwV/4YqHK6yqfEGZWX+KHtf3xGxWUT4iMOuPx/WtJhwE1p2KOMnrcX2rCQeBFaNKV80WdZLX4/pWEw4CG7qqXTVb8ElefwzRup2DwIaualfN+uItsxfEQWBDV7Udr0/ymr0gDgIbuiJ3vD7Ja9Z2DgIbuqJ2vD7JazYsHAQ2dEXteH2S12xYOAjqrMiPfBax463auQazmnAQ1FXVPvIJPslrNkwcBHVVtY98gk/ymg0TB0FdVXEYxid5zYaFZx+tq/Hjs+GgZsuHk2foNGs7HxHUlYdhzCxxENSVh2HMLCk1CCQdKmm+pAWSZjYpnyrpznT7raTdy2yPNfBn7c2MEoNA0kjgTOAwYBJwjKRJDavdD7wxIl4L/Dswu6z2mJlZc2UeEewNLIiIhRHxDDAHmJJfISJ+GxF/SQ9vBLYvsT1mZtZEmUEwDlice7wkLevPccCPmxVImiapV1LvsmXLCmxii6r0pStmZgUrMwjUZFk0XVF6E1kQnNSsPCJmR0RPRPSMHTu2wCa2oMgrcB0oZlZBZQbBEmCH3OPtgYcaV5L0WuBbwJSI+FOJ7dkwRV2B60Axs4oqMwhuAXaWtKOkUcDRwJX5FSSNBy4D3hMR95XYlg1X1BW4VQwUMzNKDIKIWA2cCFwL3ANcFBHzJB0v6fi02ueAbYCzJN0uqbes9mywoiZCq1qgmJklpU4xERFXA1c3LDsnd/+DwAfLbMOQnXJK9o47v/PdkCtwi5rSoYpzBJlZR/OVxYMp6grcoqZ08FTNZlYwB0ErirgCt2qBYmaWePbRdipiZs0UHI8cdzLbrnqQERPGZyHg6SHMbAM5CDqRp2o2swJ5aMjMrOYcBGZmNecgMDOrOQeBmVnNOQjMzGrOQWBmVnMOAjOzmnMQmJnVnIPAzKzmHARmZjXnIDAzqzkHgZlZzTkIzMxqzkFgZlZznobazKxNqjptvI8IzMxqzkFgZlZzDgIzs5rzOQIzq5SqjqN3Mx8RmJnVnI8IzGrM774HV4dt5CAws0LUYYfZrTw0ZGZWcz4iMOtAfvdtRXIQmFlXcli2zkNDZmY15yAwM6u5UoNA0qGS5ktaIGlmk3JJ+noqv1PSnmW2x8zM1lfaOQJJI4EzgbcAS4BbJF0ZEb/LrXYYsHO67QOcnX52LY9bmlnVlHmyeG9gQUQsBJA0B5gC5INgCvDdiAjgRkmjJW0XEQ+X2K6u4EAZXBW3URXbZFZmEIwDFuceL2H9d/vN1hkHrBMEkqYB09LDVZLuLraprZMKrW4MsLzQGjdQAf16ri9FbaOCt3Wr+v2dDFN7hqIyf18FcF+GbkJ/BWUGQbN/m9iAdYiI2cBsAEm9EdEz9OYNP/elerqlH+C+VFUV+1LmyeIlwA65x9sDD23AOmZmVqIyg+AWYGdJO0oaBRwNXNmwzpXAe9Onh/YFHvf5ATOz9iptaCgiVks6EbgWGAmcFxHzJB2fys8BrgYOBxYAK4FjW6h6dklNHg7uS/V0Sz/AfamqyvVF2Qd2zMysrnxlsZlZzTkIzMxqrqOCYLApK6pE0nmSluaveZD0Ekk/lfT79PPFubJPp37Nl3TI8LS6OUk7SPqFpHskzZM0PS3vuP5I2lTSzZLuSH35QlrecX2B7Ap+SbdJuio97tR+LJJ0l6TbJfWmZZ3al9GSLpF0b/qf2a/yfYmIjriRnXD+A7ATMAq4A5g03O0aoL0HAHsCd+eWnQbMTPdnAqem+5NSfzYBdkz9HDncfci1eztgz3R/S+C+1OaO6w/ZtStbpPsbAzcB+3ZiX1L7Pg58H7iqw//GFgFjGpZ1al8uAD6Y7o8CRle9L510RPDclBUR8QzQN2VFJUXEL4E/NyyeQvZHQvp5RG75nIhYFRH3k32Kau92tLMVEfFwRNya7j8J3EN2BXjH9ScyK9LDjdMt6MC+SNoe+HvgW7nFHdePAXRcXyRtRfYm8NsAEfFMRDxGxfvSSUHQ33QUneSlka6TSD+3Tcs7pm+SJgJ/S/ZOuiP7k4ZTbgeWAj+NiE7ty38CnwLW5pZ1Yj8gC+OfSJqbppSBzuzLTsAy4DtpyO5bkl5ExfvSSUHQ0nQUHaoj+iZpC+BS4GMR8cRAqzZZVpn+RMSaiNiD7Er2vSW9ZoDVK9kXSW8FlkbE3Faf0mTZsPcjZ/+I2JNsRuITJB0wwLpV7stGZEPCZ0fE3wJPkQ0F9acSfemkIOiG6SgelbQdQPq5NC2vfN8kbUwWAhdGxGVpccf2ByAdsl8HHErn9WV/4G2SFpENkx4o6Xt0Xj8AiIiH0s+lwA/Jhkc6sS9LgCXpKBPgErJgqHRfOikIWpmyouquBN6X7r8PuCK3/GhJm0jakez7GW4ehvY1JUlkY573RMT/yxV1XH8kjZU0Ot3fDHgzcC8d1peI+HREbB8RE8n+F34eEe+mw/oBIOlFkrbsuw8cDNxNB/YlIh4BFkvaJS06iGzq/Wr3ZbjPsL+QG9l0FPeRnVk/ebjbM0hbf0A2nfazZKl/HLAN8DPg9+nnS3Lrn5z6NR84bLjb39CXN5Adrt4J3J5uh3dif4DXArelvtwNfC4t77i+5No3mec/NdRx/SAbV78j3eb1/W93Yl9S2/YAetPf2OXAi6veF08xYWZWc500NGRmZiVwEJiZ1ZyDwMys5hwEZmY15yAwM6s5B4FZPyRNVG722AHWeVfucY+kr5fYpn+T9Oay6rd68sdHzfqR5lW6KiL6nYJC0mTgkxHx1jY1y6xwPiKwjpTeid8r6QJJd6b53zdPZQelCb/uUva9EJuk5Ysknars+whulvTKtPx8SUfm6l7Rz+v9StKt6fb6VPQl4O/SPPozJE3Wut8NcHlq342SXpuWz0rtuk7SQkn/0uT1RqZ23Z36MSPf1nTkcXu63SUpUvkrJF2TJm/7laRdi9zu1p0cBNbJdgFmR8RrgSeAj0jaFDgfeGdE/A3ZJGAfzj3niYjYGziDbPbOVi0F3hLZxGjvBPqGf2YCv4qIPSLi9IbnfAG4LbXvM8B3c2W7AoeQzanz+TSXU94ewLiIeE3qx3fyhRHRm15zD+Aa4CupaDbw0YjYC/gkcNYL6KPVlIPAOtniiPhNuv89sqkwdgHuj4j70vILyOaH7/OD3M/9XsBrbQycK+ku4GKyLxQZzBuA/wKIiJ8D20jaOpX9T2Rz0C8nC5mXNjx3IbCTpG9IOpQs6NYj6SiySc1mptlhXw9cnKbZ/ibZlwqZDWij4W6A2RA0nuAKmk/r299z+u6vJr0pShPsjWryvBnAo8Duad2/ttC+gaYYXpVbtoaG/8WI+Iuk3cmOGk4AjgI+sE7l0m5kRx0HRMQaSSOAx9JRglnLfERgnWy8pL539ccAvyabSXRi3/g/8B7g+txz3pn7eUO6vwjYK92fQvbuv9HWwMMRsTbVOTItf5Ls6zub+SUwFZ47qbw8Bv4eh+dIGgOMiIhLgX8le9efL9+abPrp90bEMoBU9/2S3pHWUQoTswH5iMA62T3A+yR9k2xWx7Mj4q+SjiUbHtmIbPryc3LP2UTSTWRvgo5Jy84FrpB0M9nMkE81ea2zgEvTTvYXuXXuBFZLuoPs3MRtuefMIvumqjuBlTw/DXErxqXn9r1Z+3RD+RHABLLhKgDSkcBU4GxJnyULtDlks3qa9csfH7WO1MpHO5s8ZxHQk8blzSzx0JCZWc35iMDMrOZ8RGBmVnMOAjOzmnMQmJnVnIPAzKzmHARmZjX3vw3R2GU0BCE1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAki0lEQVR4nO3df7wcVX3/8debEMOvhB/mAuFnBBEKiFFuEUVt5IcGi4KtVZAiUmu0/vgCxSqgLdHqF/QrIvpVJAiCEqCKIpTytVIkIirgjSA/jIhAEEpILtAIAQQhn+8f51wyLHfv3Zvs7M7uvJ+Pxz525szMmXNmZz8zc2b2rCICMzOrj3W6XQAzM+ssB34zs5px4DczqxkHfjOzmnHgNzOrGQd+M7OaceC3jpB0uKQfdrsca0PSbZJm5+F5ks7vbonM1kzfBn5JSyTt35D2bknXNqQdKul6SY9JWp6HPyBJefq5kp6StFLSw5KulLRLnjZP0p/ytBWSfibpVXnaCySdKum+PP1uSaeNUp5bJD0u6QFJZ0japNQN0wGSZkoKSeuOpEXEgoh4QzfLNRH5c/90MS0idouIhV0qUq3l78oz+bs08ppdmL6ZpEvy9/geSe/sULmeF2fWII/n7Wtl69vA3wpJxwGnA/8H2BLYAng/sA/wgsKsn4uIjYBtgOXAuYVp/5anDQDXAt/LB40TgEFgL2Aq8HrgxoZ1fxb4J2BjYG9ge+BKScV1V46kSd0ug6254gG5x/w8IjYqvBYWpn0FeIr0HT4cOEPSbq1k2s3t0bXvUkT05QtYAuzfkPZu4No8vDHwGPDX4+RzLvDpwvhfAivz8Dzg/MK03YAApgOXA8c0yXMasBJ4e0P6RqQDy981We6FwL8DjwC/AD49Up88fRfgSuBh4PZi/rkeXwH+A3gUuB7YcQLLngFckbfZ/nk73JjLci8wrzD/7/N2WJlfrypu+zzPq3Md/pDfX12YthD4V+Cnuaw/BKbnaesB5wMPASvysls02V4BvHi0zxKYDdwHHJe3+VLgqDxtLvAnUiBZCfx74z7V+NmPsu6DgZvy9rkTmJPTtwIuy9v5d8B7C8vMA74NfDPX+zZgME87Hri4YR2nA18q7M9n53r8d943JhX2+58Cp+X1fppy96XdCssuA07M6evketyZP79vA5u1+H1+d7F8DdM2zJ/VSwpp3wJOGSOvxu0xBfg8ad9dBnwNWD/PP/J9XpHn/0muy7eAVcATeT/5aJ7/O8ADpH37GmC3Mb5Lzfa1j+XP8dG8/fdra3xsZ2ZVejF+4J8DPA2sO04+57I6WGwEXAD8pPBFPT8PTyFdOdybxz+Rd6IPAC8FVMiz6bqB84ALm5TlovzaANiVFHBH6rNhHj8KWBd4BfDgyE6X6/Ew6QpkXWABcNEElv0D6UpoHVLwnZ3rtQ6wR/6yHJLnn0kKuus22fabAf8DHJHXd1gef2GevpAUHF4CrJ/HT8nT3kcKWBsAk4A9gWlNttd4gf9p4FPAZOBNwOPApo3zjrZPMUbgz9v4D8ABeftsDeySp/0Y+GrehrOAYfKXOuf5x1yWScDJwHV52va5fNPy+CRSkN87j38fODN/lpsDNwDvK2z7p4EP5+29PuXtS1NzuY7LdZwKvDJPOwa4jnTlPCWXd9R9fZRt+m5SoHwQ+C3wz+T9C3g58ETD/B8hB9EmeTVujy+SDsib5TL/O3Bynv9k0oFgcn69lvx9ZvQ483c5jyk535sa9sHG79K5PPfkcue8/bcqfJ92bGU7tRwf25lZlV75A1lJOkqPvB4v7Nx/CzzQsMzP8nxPAK8rfFB/zOkP5J1jx8IX9ak8bTnwI2DPwhfzg6QziyeB+4Ejm627UIZTgCtHSZ9EOjPYuZD27Fka8A7yAakw/UzgpEI9vl6Y9ibgNxNY9pvjbO8vAqcVdtSxAv8RwA0Ny/8ceHceXgh8ojDtA8APYvWX6mfAHi3sA+MF/icayric1YH02Xkb9qlWAv+ZI9uiIX1b4BlgaiHtZODcQp7/VZi2K4WARmpKfFcePgC4Mw9vkfex9QvzHgZcXdj2v+/QvnQYcGOT7bKYwpkrMCOXY8yTrzzvDsCLSMHypcCvgRPytNfy/O/ye4GFTfJq3B4iHVSKVy2vAu7Ow58CLi3uS6PtE03WtUneDzdu9l1q3NeAF+d9cX9g8njbZk1e/d7Gf0hEbDLyIgWQEQ8B0xtuQL46z/cQz73/8fmcx5YR8ZaIuLMw7dt52uYRsW9ELMp5PRMRX4mIfUgf/meAcyT9GemsZXqTtsUZeXqjAdLZyb2FtOLw9sAr803mFZJWkNo6tyzM80Bh+HHSFUyryxbXhaRXSrpa0rCkP5DujUwfpdyj2Qq4pyHtHtKZ8Xhl/Rbwn8BFku6X9DlJk1tcb6OHIuLpJutZG9uSrlgabQU8HBGPFtLGq/d6hf3kAlJgBXhnHof0+U0GlhY+vzNJZ/4jip9fmftSs7qP5HtJIc/FpAPhFk3mf1ZE3BURd0fEqoi4hRSM35YnryQ1nxZNIzWTNNO4PTYAFhXK9oOcDulK/nfADyXdJen4ZplKmiTpFEl3SnqEdGCA53437n3+kqtFxO9IV0fzgOWSLpK01VjLTFS/B/6x/Jx0lnRw2SuKiCci4iuk5oxdC+v+q+J8kjYEDgSuGiWbYdLl6TaFtG0Lw/cCPy4e6CLdAPuHForYyrLRsMwFpKufbSNiY9KlsJrM2+h+UhAo2o7UpjmmiPhTRHwyInYl3Sc4CHhXk9kfJ32hR2zZZL5RVzWBeRvdC+w4Svr9wGaSphbSWqp39h1gtqRtgLeyOvDfS9qfphc+v2kRUby5WaxP2fvSaHUfmXZgQ77rRUSr9S8KVu9vvwXWlbRTYfrLSPdIxlp+xIOkq7/dCuXaONJDG0TEoxFxXETsALwZ+EdJ+42SD6QD8sGks/WNSVe/FMo62jLP29ci4oKIeA3pexKkB0HapraBPyJWAJ8EvirpbZI2krSOpFmkNs61IukYSbMlrS9pXUlHktr9boyIP+R1f1nSHEmTJc0kfbHvI53VNpb3GeB7wDxJG+RHSosB73LgJZKOyPlNlvTn+QpjPGuy7FTS2esfJe1F2uFHDJNueu3QZNkr8vrembfNO0gHxMvHK6ik10t6aX4a4hFSU8EzTWa/CXhnPgubA/zFePkXLBuj/OM5GzhK0n55n9pa0i4RcS+pmepkSetJ2gN4D6mNfFwRMUxqBvsGqRlicU5fSroBfqqkaXmdO0oatb4d2Je2zPv/FElTJb0yT/sa8BlJ2wNIGpDU0omXpAMlbZGHdyG18V+a6/NYrs+nJG0oaR9S8H3e92g0EbEKOAs4TdLmeR1bS3pjHj5I0ovz03qPkPa3kX2ucT+ZSjoIP0Q66fjfLRThOXlI2lnSvpKmkJqZn6D5Pr5Gahv4ASLic8A/Ah8ltaktI10if4z0BV0bTwCnki6JHyS19/91RNxVWPeJpCcJHiE9GXEvqQ30ySZ5foh0FvEAaae+kLSTkZsP3gAcSjqzfIB0ljBlvIKu4bIfIH3RHgX+hfSExkh+j5Oatn6aL533bljfQ6Qz9eNIX5CPAgdFxGhNXI22BC4mbbPFpJulzX5IdTTpDG0Fqani+y3kP+JsYNdc/oksR0TcQLoxehrpRt6PWX2FcxjpLPB+4BJSu/mVE8j+AtLZ5AUN6e8iPYL8a9KV5cWkZsNmytyXDiBt9weAO0iPMkN6CukyUpPJo6QbvSMHBZSezX9tk6z3A26W9BjpxOF7PDeofoB0k3Z5rss/RMRYZ/yNPkZqzrkuN9H8F+kmK8BOeXwl6Wr9q7H6UdKTgU/k/eQjpCey7iFdxf0613E8jfvaFNK9vgdJ23BzUqxom5E709aDJH0W2DIijux2Way3eV+ql1qf8fcaSbtI2kPJXqRmgku6XS7rPd6X6q1Xf8FXV1NJl7FbkS5pTyW3c5pNkPelGnNTj5lZzbipx8ysZnqiqWf69Okxc+bMjq7z9tvT+847jz1fJ/Pq13zamVe/5tPOvKqWTzvz6td81tSiRYsejIiBxvSeCPwzZ85kaGioo+ucPTu9L1xYnbz6NZ925tWv+bQzr6rl0868+jWfNSWp8RfygJt6zMxqx4HfzKxmHPjNzGrGgd/MrGYc+M3MasaB38ysZhz4zcxqxoHfzKxmHPjNzGrGgd/MrGYc+M3MasaB38ysZhz4zcxqxoHfzKxmHPjNzGqm9MAvaZKkGyVdnsc3k3SlpDvy+6Zll8HMzFbrxBn/0cDiwvjxwFURsRNwVR43M7MOKTXwS9oG+Evg64Xkg4Hz8vB5wCFllsHMzJ6r7DP+LwIfBVYV0raIiKUA+X3z0RaUNFfSkKSh4eHhkotpZlYfpQV+SQcByyNi0ZosHxHzI2IwIgYHBp73X8FmZraGyvyz9X2At0h6E7AeME3S+cAySTMiYqmkGcDyEstgZmYNSjvjj4gTImKbiJgJHAr8KCL+FrgMODLPdiRwaVllMDOz5+vGc/ynAAdIugM4II/3rwULuOi6mfzox+vAzJmwYEF38zGz2iuzqedZEbEQWJiHHwL268R6u27BApg7ly2ffDyN33MPzJ2bhg8/vPP5mJnhX+6W6+Mfh8cff27a44+n9G7kY2aGA3+5fv/7iaWXnY+ZGQ785dpuu4mll52PmRkO/KNr143Uz3wGNtjguWkbbJDSu5GPmRkO/M/37I3Ue1iHWH0jdU2C/+GHw/z5PDBle1Yh2H57mD9/4jdk25WPmRkdeqqnp4x1I3VNAu3hh3PoWWm5hQvXolztysfMas9n/I18I9XM+pwDfyPfSDWzPufA38g3Us2szznwN+r3G6lV7PqhimUy62MO/KM5/HAO3XsJ+/7FKliypK+CftueWKpimXwAMWuJA3+dVLHrh3aVqYoHNbOKcuCvkyo+sdSuMlXxoGZWUQ78dVLFJ5baVaYqHtTMKsqBv06q+MRSu8pUxYOaWUU58NdJFZ9YaleZqnhQM6uo0rpskLQecA0wJa/n4og4SdI84L3AcJ71xIi4oqxyWIMqdv3QjjLlA8UD7/k4mz/5e9bZfrsU9PvliSyzNiqzr54ngX0jYqWkycC1kv5fnnZaRHy+xHVbHVXxoGZWQWX+2XpExMo8Ojm/oqz1WYf5mXmznlVqG7+kSZJuApYDV0bE9XnShyTdLOkcSZs2WXaupCFJQ8PDw6PNYt3iZ+bNelqpgT8inomIWcA2wF6SdgfOAHYEZgFLgVObLDs/IgYjYnBgYKDMYtpE9fsz876asT7Xkad6ImIFsBCYExHL8gFhFXAWsFcnymBt1M/PzPtqxmqgtMAvaUDSJnl4fWB/4DeSZhRmeytwa1llsJL08zPz/X41Y0a5Z/wzgKsl3Qz8gtTGfznwOUm35PTXA8eWWAYrQz8/M9/Oqxk3GVlFlfY4Z0TcDLx8lPQjylqndUg/PzO/3XapeWe09Il4tskoXz2MNBlBf2ynKlqwgIuuS/skM/tonyyBf7lra6Zfu65u19WMm4w6y/dmJsSB36yoXV1I9PMN8CrygXZCyvzlrllvascvgNvVZGSt8YF2QnzGb1aGfr4BXkX9/KRZCRz4zcpQxZ5Q+5kPtBPiwG9Wln69AV5F7TzQ1uAxXLfxm1l/aMe9mZo8huszfjOzETV5OsiBvwMWLnT/8NZn+rU5pCZPBznwm9nE9POPpWrydJADv5lNTD83h9Tk6SDf3O0hbi6ySujn5pB+7oeqwGf8ZnXSjrb5fm8OqcFjuA78Zr2gHQG7XW3z7W4O6dcbxRXmwG9Wde0K2O1qm2/zj6X69kZxhZX5D1zrSbpB0q8k3Sbpkzl9M0lXSrojv4/6Z+tmlrUrYLezbb5dzSH9fKO4wso8438S2DciXkb6Y/U5kvYGjgeuioidgKvyuJk1066AXcW2+X6+UVxhpQX+SFbm0cn5FcDBwHk5/TzgkLLKYNYX2hWwq/ioYhUPRu1S4XsXpbbxS5ok6SZgOek/d68HtoiIpQD5ffMyy2DW89oVsKvYY2gVD0btUPF7F6UG/oh4JiJmAdsAe0navdVlJc2VNCRpaHh4uLQymlVeOwN21R5VrOLBqB0qfu+iIz/giogVkhYCc4BlkmZExFJJM0hXA6MtMx+YDzA4OBidKKdZZbWj58mq6se6VfzeRZlP9QxI2iQPrw/sD/wGuAw4Ms92JHBpWWUwM+uKit+7KLOpZwZwtaSbgV+Q2vgvB04BDpB0B3BAHjcz6x8Vv3dRWlNPRNwMvHyU9IeA/cpar5lZ11W8zx930mZmVoYK37twlw1mZjXjM/4mqnaENjNrF5/xm5nVjAO/mVnNOPCbmdWMA7+ZWc1MOPBLWkfStDIKY2Zm5Wsp8Eu6QNI0SRsCvwZul/RP5RbNzMzK0OoZ/64R8Qip7/wrgO2AI8oqlJmZlafVwD9Z0mRS4L80Iv5E+lMVMzPrMa0G/jOBJcCGwDWStgceKatQZmZWnpZ+uRsRXwK+VEi6R9LryymSmZmVqaXAL2kK8NfAzIZlPlVCmczMrESt9tVzKfAHYBHwZHnFMTOzsrUa+LeJiDmllsTMzDqi1Zu7P5P00olkLGlbSVdLWizpNklH5/R5kv5b0k359aYJl9r6zsKF7hHVrFNaPeN/DfBuSXeTmnoERETsMcYyTwPHRcQvJU0FFkm6Mk87LSI+v8alNmvCBw+z8bUa+A+caMYRsRRYmocflbQY2Hqi+Vj7OTia1VtLTT0RcQ+wCfDm/Nokp7VE0kzS/+9en5M+JOlmSedI2rTJMnMlDUkaGh4ebnVVZmY2jlb76jkaWABsnl/nS/pwi8tuBHwXOCZ3+3AGsCMwi3RFcOpoy0XE/IgYjIjBgYGBVlZlZmYtaLWp5z3AKyPiMQBJnwV+Dnx5rIVyNw/fBRZExPcAImJZYfpZwOVrUG4zM1tDrT7VI+CZwvgzOa35ApKAs4HFEfGFQvqMwmxvBW5tsQxmZtYGrZ7xfwO4XtIlefwQUlAfyz6kHjxvkXRTTjsROEzSLFInb0uA97VeXDMzW1ut9tXzBUkLSY91CjgqIm4cZ5lrGf2q4IqJFtLMzNpnzMAvaVpEPCJpM9LZ+ZLCtM0i4uFyi2dmZu023hn/BcBBpD56iv3vK4/vUFK5zMysJGMG/og4KL+/qDPFMTOzsrX6HP9VraSZmVn1jdfGvx6wATA9/8J25GbtNGCrkstmZmYlGK+N/33AMaQgv4jVgf8R4CvlFcvMzMoyXhv/6cDpkj4cEWP+StfMzHpDq8/xf1nS7sCuwHqF9G+WVTAzMytHq/+5exIwmxT4ryB103wt4MBvZtZjWu2r523AfsADEXEU8DJgSmmlMjOz0rQa+J+IiFXA05KmAcvxj7fMzHpSq520DUnaBDiL9HTPSuCGsgplZmblGTfw5+6VT46IFcDXJP0AmBYRN5ddODMza79xA39EhKTvA3vm8SUll8l6hP+716w3tdrGf52kPy+1JGZm1hGttvG/Hni/pCXAY+TeOSNij7IKZtZNvpqxftZq4D9wohlL2pb0nP+WwCpgfkScnvv2/zdgJql//7dHxP9MNH8zM1szLTX1RMQ9wLbAvnn48RaWfRo4LiL+DNgb+KCkXYHjgasiYifgqjxuZmYd0mq3zCcBHwNOyEmTgfPHWiYilkbEL/Pwo8BiYGvgYOC8PNt5pP/vNTOzDmn15u5bgbeQ2veJiPuBqa2uRNJM4OXA9cAWEbE057MU2LzJMnMlDUkaGh4ebnVVZmY2jlYD/1MREeS/X5S0YasrkLQR8F3gmIh4pNXlImJ+RAxGxODAwECri5mZ2ThaDfzflnQmsImk9wL/RfoV75gkTSYF/QUR8b2cvEzSjDx9Bqn7BzMz65BWu2X+vKQDSH/AsjPwLxFx5VjL5F/8ng0sjogvFCZdBhwJnJLfL12Tgpv1Aj8WalXU6uOc5EA/ZrBvsA9wBHCLpJty2omkgP9tSe8Bfg/8zQTyNDOztdRqf/x/BXyWdCNWrP4B17Rmy0TEtaz+q8ZG+02wnGZm1iatnvF/DnhzRCwuszBmZmvDTWutafXm7jIHfTOz/jCR/vj/Dfg+8ORIYuFJHTMz6xGtBv5ppG4a3lBIC8CB38ysx7T6OOdRZRfEzMw6o9W+el4i6SpJt+bxPSR9otyimZlZGVq9uXsWqYO2PwHkv108tKxCmZl108KF/f2EUKuBf4OIaPxz9afbXRgzMytfq4H/QUk7srqTtrcBS0srlZmZlabVp3o+CMwHdpH038DdwOGllcrMzErTauA/BLgCuJp0lfAYsL+kRRFxUzlFMzOzMrTa1DMIvB/YFNgEmAvMBs6S9NFSSmZmZqVo9Yz/hcArImIlPPtXjBcDrwMWkfryMTOzHtBq4N8OeKow/idg+4h4QtKTTZYxsz7Wz4879rtWA/8FwHWSRv405c3AhfkvGH9dSsnMzKwUrXbZ8K+SrgBeQ+pj//0RMZQn++keM7MeMpF/4FpEas9viaRzgIOA5RGxe06bB7wXGM6znRgRV7RcWjMzW2utPtWzJs4F5oySflpEzMovB30zsw4rLfBHxDXAw2Xlb2Zma6bMM/5mPiTpZknnSNq02UyS5koakjQ0PDzcbDYzM5ugTgf+M4AdgVmkvn5ObTZjRMyPiMGIGBwYGOhQ8czM+l/LN3fbISKWjQxLOgu4vJPrN+tl/fzcfD/XrYo6esYvaUZh9K3ArZ1cv5mZlXjGL+lCUn8+0yXdB5wEzJY0i9S98xLgfWWt38xG57NrKy3wR8RhoySfXdb6zMysNd14qsfMzLrIgd/MrGYc+M3Maqajj3OamdVJVW+k+4zfzKxmHPjNzGrGgd/MrGYc+M3MasaB38ysZhz4zcxqxoHfzKxmHPjNzGrGgd/MrGYc+M3MasaB38ysZkoL/PnP1JdLurWQtpmkKyXdkd+b/tm6mZmVo8wz/nOBOQ1pxwNXRcROwFV53MzMOqi0wB8R1wAPNyQfDJyXh88DDilr/WZmNrpOt/FvERFLAfL75s1mlDRX0pCkoeHh4Y4V0Mys31X25m5EzI+IwYgYHBgY6HZxzMz6RqcD/zJJMwDy+/IOr9/MrPY6HfgvA47Mw0cCl3Z4/WZmtVfm45wXAj8HdpZ0n6T3AKcAB0i6Azggj5uZWQeV9p+7EXFYk0n7lbVOMzMbX2Vv7pqZWTkc+M3MasaB38ysZhz4zcxqxoHfzKxmHPjNzGrGgd/MrGYc+M3MasaB38ysZhz4zcxqxoHfzKxmHPjNzGrGgd/MrGYc+M3MasaB38ysZhz4zcxqprQ/YhmLpCXAo8AzwNMRMdiNcpiZ1VFXAn/2+oh4sIvrNzOrJTf1mJnVTLcCfwA/lLRI0tzRZpA0V9KQpKHh4eEOF8/MrH91K/DvExGvAA4EPijpdY0zRMT8iBiMiMGBgYHOl9DMrE91JfBHxP35fTlwCbBXN8phZlZHHQ/8kjaUNHVkGHgDcGuny2FmVlfdeKpnC+ASSSPrvyAiftCFcpiZ1VLHA39E3AW8rNPrNTOzxI9zmpnVjAO/mVnNOPCbmdWMA7+ZWc048JuZ1YwDv5lZzTjwm5nVjAO/mVnNOPCbmdWMA7+ZWc048PeIiBhzvJf1c93MqsiBvwfMmzePY4899tmAGBEce+yxzJs3r7sFa4N+rptZVTnwV1xEsGLFCk4//fRnA+Sxxx7L6aefzooVK3r67Lif62ZWaRFR+deee+4ZdbZq1ao4+uijg/SXlQHE0UcfHatWrep20dZaP9fNrNuAoRglpip64KxqcHAwhoaGul2MrooI1lln9QXaqlWryP9p0PP6uW5m3SRpUUQMNqa7qWcUjQfDtTk4tiOvyE0gRcV28U6Xp9359HPd2pFPFcvkunUun1KMdhlQ9guYA9wO/A44frz5O9nUc9JJJz2nqWGkKeKkk07qSl7FppCRvBrHO1medubTz3Wr2n5UxXyqWKaq5bO2aNLU042gPwm4E9gBeAHwK2DXsZbpVOBvZyCqWlBrV3naWa9+rVsV96Oq5VPFMlUtn3aoUuB/FfCfhfETgBPGWqaTZ/ztvNnY7rzGGu9kedp9Q7Yf61bF/ahq+VSxTFXLZ201C/wdv7kr6W3AnIj4+zx+BPDKiPhQw3xzgbl5dHfg1o4WFPYsDC9qY173AA+uZX5rqx11mw5s34Z82m1N6jad538m7fr8y9qPmuU1Wl3WJJ92lWdt8mqlLq3k067ydDOfNbV9RAw0Jnb8z9aB0R7XeN7RJyLmA/MBJA3FKHeme1G/1CXXY2a3y9EO/fKZgOtSRVWsRzee6rkP2LYwvg1wfxfKYWZWS90I/L8AdpL0IkkvAA4FLutCOczMaqnjTT0R8bSkDwH/SXrC55yIuG2cxeaXX7KO6Ze69Es9wHWpqn6pS+Xq0RO/3DUzs/bxL3fNzGrGgd/MrGYqHfglzZF0u6TfSTq+2+UZj6RzJC2XdGshbTNJV0q6I79vWph2Qq7b7ZLe2J1Sj07StpKulrRY0m2Sjs7pPVUfSetJukHSr3I9PpnTe6oeRZImSbpR0uV5vCfrImmJpFsk3SRpKKf1al02kXSxpN/k78yrKl2X0X7VVYUXa9C1Q7dfwOuAVwC3FtI+R+6PCDge+Gwe3jXXaQrwolzXSd2uQ6HcM4BX5OGpwG9zmXuqPqTfjWyUhycD1wN791o9Gur0j8AFwOU9vo8tAaY3pPVqXc4D/j4PvwDYpMp1qfIZ/17A7yLiroh4CrgIOLjLZRpTRFwDPNyQfDBppyC/H1JIvyginoyIu0kd1u3ViXK2IiKWRsQv8/CjwGJga3qsPpGszKOT8yvosXqMkLQN8JfA1wvJPVmXJnquLpKmkU76zgaIiKciYgUVrkuVA//WwL2F8ftyWq/ZIiKWQgqmwOY5vWfqJ2km8HLS2XLP1Sc3jdwELAeujIierEf2ReCjwKpCWq/WJYAfSlqUu2iB3qzLDsAw8I3cBPd1SRtS4bpUOfC31LVDD+uJ+knaCPgucExEPDLWrKOkVaI+EfFMRMwi/Up8L0m7jzF7Zesh6SBgeUS02udLZeuS7RMRrwAOBD4o6XVjzFvluqxLauI9IyJeDjxGatppput1qXLg75euHZZJmgGQ35fn9MrXT9JkUtBfEBHfy8k9W598+b2Q9H8QvViPfYC3SFpCavrcV9L59GZdiIj78/ty4BJSc0cv1uU+4L58JQlwMelAUNm6VDnw90vXDpcBR+bhI4FLC+mHSpoi6UXATsANXSjfqCSJ1Ga5OCK+UJjUU/WRNCBpkzy8PrA/8Bt6rB4AEXFCRGwTqXO8Q4EfRcTf0oN1kbShpKkjw8AbSD3w9lxdIuIB4F5JO+ek/YBfU+W6dPtu+Dh3yt9EeprkTuDj3S5PC+W9EFgK/Il0VH8P8ELgKuCO/L5ZYf6P57rdDhzY7fI31OU1pMvPm4Gb8utNvVYfYA/gxlyPW4F/yek9VY9R6jWb1U/19FxdSO3iv8qv20a+371Yl1y2WcBQ3s++D2xa5bq4ywYzs5qpclOPmZmVwIHfzKxmHPjNzGrGgd/MrGYc+M3MasaB36xA0kwVelcdY553FsYHJX2pxDJ9StL+ZeVv9ePHOc0Kcr9El0dE024dJM0GPhIRB3WoWGZt5TN+6xn5TPs3ks6TdHPu/3yDPG2/3EHWLUr/izAlpy+R9FmlPvlvkPTinH6upLcV8l7ZZH0/kfTL/Hp1nnQK8Nrcj/yxkmbruX3jfz+X7zpJe+T0eblcCyXdJel/jbK+Sblct+Z6HFssa76yuCm/bpEUefqOkn6QOzv7iaRd2rndrf848Fuv2RmYHxF7AI8AH5C0HnAu8I6IeCmp06x/KCzzSETsBfxfUu+WrVoOHBCpI7F3ACPNOccDP4mIWRFxWsMynwRuzOU7EfhmYdouwBtJfdKclPtCKpoFbB0Ru+d6fKM4MSKG8jpnAT8APp8nzQc+HBF7Ah8BvjqBOloNOfBbr7k3In6ah88ndS2xM3B3RPw2p59H6h99xIWF91dNYF2TgbMk3QJ8h/QHGuN5DfAtgIj4EfBCSRvnaf8RqQ/2B0kHlS0alr0L2EHSlyXNIR3YnkfS20mdgB2fe099NfCd3PX0maQ/0TFrat1uF8BsghpvSgWjd3PbbJmR4afJJz65Q7oXjLLcscAy4GV53j+2UL6xutx9spD2DA3fv4j4H0kvI10VfBB4O/B3z8lc2o10VfG6iHhG0jrAinwVYNYSn/Fbr9lO0shZ+2HAtaTeNmeOtN8DRwA/LizzjsL7z/PwEmDPPHww6ey+0cbA0ohYlfOclNMfJf0d5WiuAQ6HZ28CPxhj/4/BsyRNB9aJiO8C/0w6qy9O35jUHfO7ImIYIOd9t6S/yfMoHzzMmvIZv/WaxcCRks4k9Xp4RkT8UdJRpOaOdUlden+tsMwUSdeTTnQOy2lnAZdKuoHUc+Jjo6zrq8B3c1C9ujDPzcDTkn5FurdwY2GZeaR/YroZeJzV3fK2Yuu87MgJ2QkN0w8Btic1PwGQz/QPB86Q9AnSAewiUq+XZqPy45zWM1p51HKUZZYAg7ld3cxwU4+ZWe34jN/MrGZ8xm9mVjMO/GZmNePAb2ZWMw78ZmY148BvZlYz/x/UW2fuhN37FgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "F_g_best = np.max(F_g, axis=1)\n",
    "F_g_avg  = np.mean(F_g,axis=1)\n",
    "F_g_std  = np.std(F_g, axis=1)\n",
    "\n",
    "# Plot. X = population. Y = fitness\n",
    "fig, ax = plt.subplots()\n",
    "title = \"PSO fitness. \" + str(trials) + \" restarts\"\n",
    "ax.set_title(title)\n",
    "ax.set_xlabel(\"population size\")\n",
    "ax.set_ylabel(\"fitness\")\n",
    "ax.set_xlim(0, 600 * 1.1)\n",
    "ax.set_ylim(0, 1 * 1.1)\n",
    "plt.scatter(POPs, F_g_best, c='#000000',marker='x')\n",
    "plt.errorbar(POPs, F_g_avg, F_g_std, ecolor='blue', c='#FF0000', linestyle='None', marker='o')\n",
    "plt.show()\n",
    "\n",
    "G_c_best = np.min(G_c, axis=1)\n",
    "G_c_avg  = np.mean(G_c,axis=1)\n",
    "G_c_std  = np.std(G_c, axis=1)\n",
    "\n",
    "# Plot. X = population. Y = generations\n",
    "fig, ax = plt.subplots()\n",
    "title = \"PSO generations until convergence. \" + str(trials) + \" restarts\"\n",
    "ax.set_title(title)\n",
    "ax.set_xlabel(\"population size\")\n",
    "ax.set_ylabel(\"generations\")\n",
    "ax.set_xlim(0, 600 * 1.1)\n",
    "ax.set_ylim(0, 40 * 1.1)\n",
    "plt.scatter(POPs, G_c_best, c='#000000',marker='x')\n",
    "plt.errorbar(POPs, G_c_avg, G_c_std, ecolor='blue', c='#FF0000', linestyle='None', marker='o')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "female-property",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}