hmw2_2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Homework 2 - Problem 2\n",
    "### Optimization theory with applications\n",
    "\n",
    "\n",
    "#### Andrea Senacheribbe\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The minimisation problem we are considering is:\n",
    "\n",
    "$$ \\min_{x}{f(x)} \\quad \\textit{st} \\quad \\sum_{i=1}^{m} x_i \\le 100 \\quad \\textit{and} \\quad x_i \\ge 0, \\forall i \\ge 1 $$\n",
    "\n",
    "where\n",
    "$f(x_1, x_2, ..., x_m)= \\sum_{i=1}^{m} a_i(x_i-b_i)^2 + 3$ $\\quad$\n",
    "\n",
    "---\n",
    "\n",
    "We can derive the Lagrangian, by rewriting the inequalities in the form $h_i(x)\\le 0$ and introducing the multipliers $u$:\n",
    "\n",
    "$$\\mathcal{L} (x,u)= f(x) + u_0\\left(\\sum_i{x_i} -100\\right) + u_{1...m}^T(-x)$$\n",
    "Note that $u$ is a vector with $m+1$ elements: $u_0$ is associated with the constraint on the summation, while  $u_{1...m}$ with the constraints on the individual $x_i$. Notice that the index $i$ in $x_i$ starts from 1.\n",
    "\n",
    "\n",
    "The dual problem becomes therefore\n",
    "\n",
    "\n",
    "\n",
    "$$\\max_{u}{g(u)} \\quad \\textit{st} \\quad u_i\\ge0, \\forall i$$\n",
    "\n",
    "where $g(u)= \\min_{x}{\\mathcal{L} (x,u)}$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1) KKT conditions\n",
    "\n",
    "We can write the KKT conditions as follow\n",
    "\n",
    "1) stationarity \n",
    "$$ \\frac{\\partial \\mathcal{L}}{\\partial x_i} = 2 a_i (x_i - b_i) +u_0 - u_i = 0$$\n",
    "\n",
    "\n",
    "2) complementary slackness\n",
    "$$ u_0\\left(\\sum_i{x_i} -100\\right) = 0$$\n",
    "$$ u_i x_i =0,\\, \\forall i \\ge 1$$\n",
    "\n",
    "\n",
    "3) primal feasibility\n",
    "$$\\sum_{i} x_i \\le 100$$\n",
    "$$-x_i \\le 0, \\forall i \\ge 1$$\n",
    "\n",
    "\n",
    "\n",
    "4) dual feasibility\n",
    "$$u_i\\ge0,\\, \\forall i$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Some insights on the problem using the KKT\n",
    "From 1) and 2) for a generic $i\\ge 1$ we can state the following\n",
    "\n",
    "$\\textit{if} \\quad  x_i\\neq 0 \\implies u_i = 0 $\n",
    "\n",
    "$\\quad  2 a_i (x_i - b_i) + u_0 = 0 \\iff x_i = b_i - \\frac{u_0}{2 a_i}$\n",
    "\n",
    "Since $ x_i > 0 $, we have $ b_i - \\frac{ u_{0} }{2 a_i} = x_i >  0 \\implies u_{0} < 2 a_i b_i $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "From 1) and 4) we can write\n",
    "\n",
    "$ 2 a_i (x_i - b_i) +u_{0} = u_i \\ge 0 \\implies  2 a_i (x_i - b_i) +u_{0} \\ge 0 \\implies  x_i \\ge b_i - \\frac{u_{0}}{2 a_i}$\n",
    "\n",
    "We can use this result to derive the complementary case of above, using again 2)\n",
    "\n",
    "$\\textit{if} \\quad  u_i\\neq 0 \\implies x_i = 0$\n",
    "\n",
    "$\\quad 0= x_i \\ge b_i - \\frac{u_0}{2 a_i} \\implies u_{0} \\ge 2 a_i b_i  $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "Overall we can write then\n",
    "\n",
    "$x_i=\\begin{cases}\n",
    "    b_i - \\frac{u_{0}}{2 a_i} & u_{0} < 2 a_i b_i\\\\\n",
    "    0 & u_0 \\ge 2 a_i b_i\n",
    "  \\end{cases}$\n",
    "  \n",
    "or equivalenty $x_i=max(0, b_i - \\frac{u_{0}}{2 a_i})$. To have $x$ feasible we should also impose $\\sum_{i} x_i \\le 100$.\n",
    "\n",
    "With this method we didn't end with a complete analytic solution, but we can at least understand the actual form of the solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = lambda x,a,b: np.sum(a*((x-b)**2))+3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((500,), (500,))"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m=500\n",
    "a=np.random.randint(1,100+1, size=m)\n",
    "b=np.random.randint(1,100+1, size=m)\n",
    "a.shape, b.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2) Dual ascent\n",
    "\n",
    "The dual ascent algorithm allows to solve a minimisation problem with equalities and inequalities constraints, by maximising the dual of $f$.\n",
    "\n",
    "Our optimisation problem can be reformulated as\n",
    "\n",
    "$$\\min_{x}{f(x)} \\quad \\textit{st} \\quad  Cx \\le d $$\n",
    "\n",
    "where $C$ is a matrix $(m+1, m)$, associated with the constrain reported above.\n",
    "\n",
    "In particular the first row of C is full of ones to represent the constraint $\\sum_{i} x_i \\le 100$, while the other $m$ rows have $-1$ on the diagonal, and represent the constraints $-x_i \\le 0, \\forall i \\ge 1$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  1.,  1., ...,  1.,  1.,  1.],\n",
       "       [-1., -0., -0., ..., -0., -0., -0.],\n",
       "       [-0., -1., -0., ..., -0., -0., -0.],\n",
       "       ...,\n",
       "       [-0., -0., -0., ..., -1., -0., -0.],\n",
       "       [-0., -0., -0., ..., -0., -1., -0.],\n",
       "       [-0., -0., -0., ..., -0., -0., -1.]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C=np.concatenate([np.ones((1,m)),-np.eye(m)],axis=0)\n",
    "C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(501, 500)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "$d$ is instead formed by putting $100$ in the first row (for the copuling constraint on the summation) and $0$ elsewhere. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([100.,   0.,   0.,   0.])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d=np.zeros(m+1, dtype=np.float64)\n",
    "d[0]=100\n",
    "d[:4] # printing first 4 rows"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "To compare the result with our algorithm, we firstly solve the problem using ```scipy.optimize``` "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(True, 'Optimization terminated successfully.', 80630285.25494224)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.optimize import LinearConstraint, minimize\n",
    "\n",
    "# here optimise sqrt(f(x)), since grows less rapidly and allow better convergence of the algorithm\n",
    "# x that minimises this problem will be the same of our problem\n",
    "min_sol=minimize(lambda x: np.sqrt(f(x,a,b)),np.ones(m), constraints=LinearConstraint(C,-np.inf, d))\n",
    "f_opt=f(min_sol.x, a,b)\n",
    "min_sol.success, min_sol.message, f_opt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "Essentially the dual ascent algorithm works by doing gradient ascent on the dual variable $u$ and updating at each iteration the value $x$ by solving a minimisation problem.\n",
    "\n",
    "$$\n",
    "x^{k+1}= argmin_{x}{(f(x) + (u ^k)^T C x)}$$\n",
    "and\n",
    "$$\n",
    "u^{k+1}=max(u^k+t_k ( A x^{k+1} -b),0)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "For our case we can derive a closed form for the argmin, by simply computing the gradient and setting it equal to 0:\n",
    "\n",
    "$$ \n",
    "\\nabla f(x) + (u ^k)^T C  =0\n",
    "$$\n",
    "therefore\n",
    "$$\n",
    "2 a_i (x_i - b_i) +(u ^k)^T C_i =0 \\implies x_i=b_i - \\frac{(u ^k)^T C_i}{2a_i}\n",
    "$$\n",
    "\n",
    "To help the algorithm we can clip to zero all the negative values of $x_i$ (recall that $x_i$ should be non negative).\n",
    "Notice that this step is not part of the dual ascend algorithm (indeed the algorithm works even when not clipping the negative values to 0), but from our tests this provides faster and better convergence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "argmin_f=lambda u: np.maximum(b-(u@C)/(2*a),0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dual_ascent(u0, argmin, fun, t, C, d, eps=1e-5, max_iter=10000):\n",
    "    assert(u0.shape[0]==C.shape[0]==d.shape[0])\n",
    "    n_var,n_con=C.shape[1], u0.size\n",
    "    \n",
    "    x=np.empty(n_var, dtype=np.float64)\n",
    "    u=np.empty(n_con, dtype=np.float64)\n",
    "    f=np.empty(max_iter, dtype=np.float64)\n",
    "    \n",
    "    u=u0.copy()\n",
    "    \n",
    "    for k in range(max_iter-1):\n",
    "        \n",
    "        # solve argmin problem\n",
    "        x=argmin(u)\n",
    "        f[k]=fun(x)\n",
    "        \n",
    "        if  np.abs(f[k]-f[k-1]) < eps:\n",
    "            print(\"converged after %d iterations\" % k)\n",
    "            break\n",
    "\n",
    "        #update u\n",
    "\n",
    "        u=np.maximum(u+t*(C@x-d),0)\n",
    "        \n",
    "    else: #reached max_iter-1\n",
    "        k=max_iter-1\n",
    "        x=argmin(u)\n",
    "        f[k]=fun(x)\n",
    "\n",
    "        print(\"stopped since reached max_iter\")\n",
    "        \n",
    "    \n",
    "    return x, f[:k+1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(501,)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u0=np.ones(m+1)\n",
    "u0.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "t=0.15"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "converged after 1977 iterations\n"
     ]
    }
   ],
   "source": [
    "x_da,f_da=dual_ascent(u0, argmin_f ,lambda x: f(x,a,b), t, C, d, max_iter=5000, eps=1e-4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAt8AAAFzCAYAAADmPGXMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeXSf1WHn//fVZlnebVnyIukrG2/YxuDdYl8DIYQtQMBg0WwkaTNtJ53p6TJtOr/fdJrOTNomTaY0ISSRwUBYw5KwBAgYkHdjYxtjjK2vJC/yvq+S7vxhpXEIOBZIer6S3q9zdKzvI+mrj+55JH18de/zhBgjkiRJktpfVtIBJEmSpO7C8i1JkiR1EMu3JEmS1EEs35IkSVIHsXxLkiRJHcTyLUmSJHWQnKQDdJTCwsJYXl6edAxJkiR1cUuXLt0RYxz8QW/rNuW7vLycJUuWJB1DkiRJXVwIIf1hb3PZiSRJktRBLN+SJElSB7F8S5IkSR3E8i1JkiR1EMu3JEmS1EEs35IkSVIHsXxLkiRJHcTyLUmSJHUQy7ckSZLUQSzfkiRJUgexfEuSJEkdxPLdzpbU7GJl/Z6kY0iSJCkDWL7b2V8/vorvvrQ+6RiSJEnKAJbvdtavZy77jhxPOoYkSZIygOW7nfXtmcPew41Jx5AkSVIGsHy3s775uew77My3JEmSLN/trm9Py7ckSZJOsHy3s0G98th/tJEjx5uSjiJJkqSEWb7bWdmgAgBqdx1KOIkkSZKSZvluZ6lBvQBI77R8S5IkdXeW73ZW3jLznd55MOEkkiRJSprlu531L8ijX89cNuywfEuSJHV3lu8OMH5oX1Zv2pt0DEmSJCXM8t0BJpX04+0t+znW2Jx0FEmSJCXI8t0Bzirpx7GmZtY17E86iiRJkhJk+e4AZ5f0B2BF/Z6Ek0iSJClJlu8OUDKgJ4N65bF4466ko0iSJClBlu8OEELg3FGFvP7eTmKMSceRJElSQizfHeSCUYVs33+UdQ0Hko4iSZKkhFi+O8h5owsBeGXdtoSTSJIkKSmW7w4yvH9PJg7vy9MrtyQdRZIkSQmxfHega88exsr6vWz0bpeSJEndkuW7A10zaRgAT63YnHASSZIkJcHy3YGG9e9JxchBPLS4jsYm73YpSZLU3Vi+O9jnzx/Bpj2HeXb11qSjSJIkqYNZvjvYZeOKKB9UwA/mb/Sa35IkSd2M5buDZWUFvnD+CFbU7WFpenfScSRJktSBOm35DiFcEEK4O4RwTwjhjaTztMZnppbQr2cu3391Q9JRJEmS1IEyqnyHEO4NIWwLIax63/GrQgjvhBDWhxD+AiDGOD/G+BXgaeAnSeT9qArycrjz3HKeX9PAyvo9SceRJElSB8mo8g38GLjq5AMhhGzge8AngfHAbSGE8Se9y2xgXkcFbCtfumAEA3vl8Y/Prk06iiRJkjpIRpXvGOOrwK73HZ4BrI8xbogxHgMeBK4DCCGUAXtjjPs/6PlCCHeFEJaEEJZs3769PaO3Wp/8XL52ySheX7+TX73jLeclSZK6g4wq3x9iOFB30uP6lmMAXwB+9GEfGGP8foxxWoxx2uDBg9sx4kdz+6wyRhT24u+eXM2R401Jx5EkSVI76wzl+0PFGL8RY+xUmy1P1iMnm//vugnU7Dzk5ktJkqRuoDOU701A6UmPS1qOdQkXjB7MNZOG8t2X15PeeTDpOJIkSWpHnaF8LwZGhxBGhBDygFuBJxPO1Kb+5prx5GVn8V8fWUlTszfekSRJ6qoyqnyHEB4AqoGxIYT6EMIXYoyNwNeA54C3gZ/GGFcnmbOtFffN579fO4FFG3fx76++l3QcSZIktZOcpAOcLMZ424cc/znw8w6O06FunDKcl97Zxj89v46zhvfjgtGZt0FUkiRJH09GzXx3ZyEEvnnjWYwq6s1X71vGms37ko4kSZKkNmb5ziB98nP50eem0yc/h8/9eBG1Ow8lHUmSJEltyPKdYYb268mPPzeDo43NfPb71Wzc4RVQJEmSugrLdwYaO6QPD3xp1okC/u/VrN92IOlIkiRJagOW7wx15tC+PHjXLJoj3Pr9at7Zuj/pSJIkSfqYLN8ZbExxHx68axbZWYGb7n6D19fvSDqSJEmSPgbLd4YbVdSbR796LsP69eTOexfx08V1SUeSJEnSR2T57gRKBhTw8FcrqDhjEH/+6Er+93NrafZOmJIkSZ2O5buT6Jufy71/MJ3bZpTyvZff448fXM6R401Jx5IkSVIrZNQdLnVqudlZ/M8bzqJ8UC+++exaNu44yN13TKV0YEHS0SRJknQanPnuZEIIfPmiM/jhndOo3XWIa7/7GvPf3Z50LEmSJJ0Gy3cndem4Yp762vkU9cnnznsX8W+/eo8YXQcuSZKUySzfnVh5YS8e+8Nz+eRZQ/nHZ9fypaql7Dl0LOlYkiRJ+hCW706uV48cvnvbZL7x6fG8sm4bn/rOayyv3Z10LEmSJH0Ay3cXEELgc+eN4JGvnEsIcPPd1dwzf4PLUCRJkjKM5bsLObu0P8/88QVcdmYR/+OZt/lS1VL2HjqedCxJkiS1sHx3Mf165nL3HVP522tOLEO5+jvzXYYiSZKUISzfXVAIgc+fP4KHv3IuALf8ezU/fG2jy1AkSZISZvnuws4p7c/P//gCLh5bxP//9Bq+PNdlKJIkSUmyfHdx/Qpy+f6cqfzNNeN5ae02PvWv83mzbk/SsSRJkroly3c3EELgC+eP4OGvVBAj3Hz3G9wzfwPNzS5DkSRJ6kiW725kctkAnvnj87lozImrodx09xusa9ifdCxJkqRuw/LdzfQvyOMHlVP5p1vOZsOOg3zqO/P55xfWcbSxKelokiRJXZ7luxsKIXDjlBJ++fWLuPqsoXz7xXe55juvsTTtJQklSZLak+W7Gyvs3YNv3zqZH/3BdA4ebeSmu9/gb3+2in1HvCKKJElSe7B8i0vGFfH81y/izopy5i5Ic/m3XuHplZu9LrgkSVIbs3wLgN49cvi7ayfwxB+ex+A+PfjavOXc+aPFpHceTDqaJElSl2H51m85u7Q/P/uj8/jGp8ezLL2bT/zzq/zri++6IVOSJKkNWL71O3Kys/jceSP45dcv4vIzi/nWC+u4+tvzqX5vZ9LRJEmSOjXLtz7UkH75fO/2Kfzoc9M51tTMbT9YwNcfepMdB44mHU2SJKlTsnzr97pkbBHP/+lFfO2SUTy1cjOXfesV5i2s9Q6ZkiRJrWT51mnpmZfNf7lyLL/4kwsYN6QPf/X4W9x09xus2bwv6WiSJEmdhuVbrTKqqA8P3jWLb918NjU7D/Hp777G3z+zhoNHG5OOJkmSlPEs32q1EAKfmVrCS392EbdMK+EH8zdy+T+9wnOrt3ptcEmSpFOwfOsj61+Qxz/cOIlHv1pBv565fHnuUr5UtYT63YeSjiZJkpSRLN/62KamBvLUfzqfv7p6HK+v38kV//Qqd7/yHsebmpOOJkmSlFEs32oTudlZ3HXhGfzyzy7i/NGFfPMXa7nmO6+xYIPXBpckSfo1y7fa1PD+PflB5TR+UDmNA0cbufX7C/ij+5dRt8ulKJIkSTlJB1DXdMX4Ys4fVcj3X93Av72ynhfebuDLF47kqxefQUGep50kSeqenPlWu+mZl82fXD6al/7sYj45cQj/+tJ6Lvk/v+Lx5fXeoEeSJHVLlm+1u2H9e/LtWyfz6FcrKO6bz39+aAXX/9/XeeO9HUlHkyRJ6lCWb3WYqamBPPGH5/Gtm89mx/6jzP7BQj7/48W8s3V/0tEkSZI6ROguN0WZNm1aXLJkSdIx1OLI8SZ+8kYN3315PQePNnLT1BK+fsVYhvTLTzqaJEnSxxJCWBpjnPaBb7N8K0m7Dx7jey+vp6o6TVYWfOH8EXz5ojPom5+bdDRJkqSPxPKN5TvT1e06xLeef4cn3tzMwF55/PGlo5g9M0VejiujJElS53Kq8m2zUUYoHVjAv9w6mae+dj7jhvTh755awxX//AqPL6+nySujSJKkLsLyrYxyVkk/7v/iTH78uekU5OXwnx9awZX/8irPrNzi5QklSVKnZ/lWxgkhcPHYIp75T+fzb7dPIQB/NG8ZV39nPs+v3kp3WSolSZK6Hsu3MlZWVuCTZw3l2T+9kG/feg5HG5u5a+5Srvve67z8zjZLuCRJ6nTccKlOo7GpmceWb+I7L75L/e7DTCnrz599YiznnjGIEELS8SRJkgCvdgJYvruSY43NPLy0ju++tJ4te48wc8RA/uwTY5kxYmDS0SRJkizfYPnuio4cb+KBRbV87+X32HHgKBeMLuTrV4xhctmApKNJkqRurEteajCEcHEIYX4I4e4QwsVJ51HHy8/N5nPnjWD+n1/CX109jtWb93HD/32DL/x4Mas27U06niRJ0u/IqPIdQrg3hLAthLDqfcevCiG8E0JYH0L4i5bDETgA5AP1HZ1VmaNnXjZ3XXgGr/75JfzXK8eyuGYX1/zra3xl7lLe2bo/6XiSJEn/IaOWnYQQLuREoa6KMU5sOZYNrAOu4ETJXgzcBqyNMTaHEIqBf4ox3n6q53bZSfex9/BxfvjaRu59bSMHjjZyxfhivnLRGUxNuRxFkiS1v06z7CTG+Cqw632HZwDrY4wbYozHgAeB62KMzS1v3w30+KDnCyHcFUJYEkJYsn379nbLrczSr2cuX79iDPP//BL++NJRLNq4i8/82xvccnc1L6/1EoWSJCk5GTXzDRBCKAeePmnm+ybgqhjjF1sezwFmAi8BVwL9gX+LMf7qVM/rzHf3dfBoIw8truOe+RvYvPcIY4v78OWLRnLNpGHk5WTU/z8lSVIXcKqZ75yODtNWYoyPAY8lnUOZr1ePHD5//gjmVKR4asVm/v2VDXz9pyv4X8++w53nljN7Rhn9CnKTjilJkrqBzlC+NwGlJz0uaTkmtUpudhY3Tinh+nOG88q72/nh/I3847Nr+deX3uWWaaV8/rwRlA0qSDqmJEnqwjpD+V4MjA4hjOBE6b4VmJ1sJHVmWVmBS8YWccnYItZs3sc9r23g/oVpqqpruHLCEL54wQimprxhjyRJansZteY7hPAAcDFQCDQA34gx/jCEcDXwL0A2cG+M8e9b+9yu+dapNOw7wk/eqOH+hbXsPXycyWX9+dIFI/nE+GJysl0XLkmSTp93uMTyrdNz6Fgjjyyt54evbSS98xAlA3ry+fNGcMv0Unr36Ax/KJIkSUmzfGP5Vus0NUd++XYD98zfwOKa3fTJz2H2jDL+4LxyhvbrmXQ8SZKUwSzfWL710b1Zt4d75m/gF6u2EoBPTRrKly4YycTh/ZKOJkmSMpDlG8u3Pr66XYf48Rs1PLS4jgNHG5k1ciBfPH8kl44rIisrJB1PkiRlCMs3lm+1nX1HjvPQojp+9PpGNu89wsjCXnz+/BF8ZkoJPfOyk44nSZISZvnG8q22d7ypmV+s2so98zewsn4vAwpymTMrxZyKcgb36ZF0PEmSlBDLN5ZvtZ8YI4s27uKe1zbyy7cbyM3K4vrJw/jiBSMZU9wn6XiSJKmDdcnby0uZIoTAzJGDmDlyEBu2H+BHr9fw8NI6frqkngvHDOZLF4zg/FGFhOC6cEmSujtnvqV2sPvgMe5fmOYn1Wm27z/KuCF9+ML5I7j2nGH0yHFduCRJXZnLTrB8KxlHG5t48s3N/PC1jazdup/BfXpwZ0WK22emGNArL+l4kiSpHVi+sXwrWTFGXlu/gx/M38ir67aTn5vFZ6aUcOe55a4LlySpi3HNt5SwEAIXjB7MBaMH887W/fzwtQ08vLSe+xfWMnPEQCoryvnEhGJys7OSjipJktqRM99SQnYdPMZPl9Rx34I09bsPU9y3B7NnpLhtRilFffOTjidJkj4il51g+VbmamqOvLx2G1UL0ry6bjs5WYGrJg6hsqKc6eUDvEqKJEmdjMtOpAyWnRW4fHwxl48vZuOOg9y3IM3DS+p4euUWxg3pw5yKFNefM5xePfx2lSSps3PmW8pAh4818bM3N1FVnWbNln306ZHDZ6aWMKcixRmDeycdT5IknYLLTrB8q3OKMbKsdjdV1Wl+/tYWjjdFzh9VyJyKFJeNKyLHDZqSJGUcyzeWb3V+2/cf5aHFtdy/sJYte48wvH9PZs8s49bppQzq3SPpeJIkqYXlG8u3uo7GpmZ++XYDVdVp3nhvJ3nZWXxq0lDmVKSYXNrfDZqSJCXM8o3lW13T+m37mVud5tFlmzhwtJGJw/tSOauca88ZRn6ut7GXJCkJlm8s3+raDhxt5PHlm5hbXcO6hgP0L8jllmml3DEzRdmggqTjSZLUrVi+sXyre4gxsnDjLuZWp3l29VaaY+TiMYOprCjnojGDycpySYokSe3N8o3lW91Pw74jzFtYy7xFtWzff5SygQXcMauMW6aV0r8gL+l4kiR1WZZvLN/qvo41NvPc6q3MrU6zqGYXPXKyuPbsYVRWlHNWSb+k40mS1OVYvrF8SwBrt+6jqjrNE8s3cehYE+eU9qeyIsWnJg2lR44bNCVJaguWbyzf0sn2HTnOo0vrmbsgzYbtBxnUK4/PTi/l9lkphvfvmXQ8SZI6Ncs3lm/pg8QYeX39Tqqqa/jl2w0AXHZmMZUVKc47o9ANmpIkfQSnKt85HR1GUuYIIXD+6ELOH13Ipj2HmbcwzYOL6nhhTQMjC3txx6wUn5laQr+euUlHlSSpS3DmW9JvOdrYxM/f2kJVdZrltXvomZvN9ZOHU1mR4syhfZOOJ0lSxvvYy05CCGWn+bn2xBj3tSZcR7F8S623atNeqqpr+Nmbmzna2Mz08gHMqSjnqglDyMvJSjqeJEkZqS3K98tABE61ADQCP44xVn2klO3M8i19dHsOHePhJSc2aNbuOsTgPj24bUYZs2eUMaRfftLxJEnKKG64xPIttYXm5sgr67ZTVV3Dr9ZtJysErpxQzJxZ5cwaOZAQ3KApSZIbLiW1iayswCXjirhkXBG1Ow9x38I0P11Sx8/f2sqY4t7MmZXihikl9O7hjxZJkj7I7535DiH0ijEeDCH0jjEe6KBcbc6Zb6l9HDnexJMrNjO3Os1bm/bSu0cON045sUFzVFGfpONJktThPtaykxBCCXA9sD7G+Gw75OsQlm+pfcUYebNuD3Or0zy9cgvHmpqpGDmIyooUV4wvJifbDZqSpO7h45bvO4H/BNwLPBJj3Nb2Eduf5VvqODsPHOWhJXXcv6CWTXsOM7RfPrNnlHHrjDIG9+mRdDxJktrVxy3fZwK5QGmM8Zl2yNchLN9Sx2tqjrz4dgNzF6SZ/+4OcrMDn5w4lMqKFFNTA9ygKUnqkrzaCZZvKWkbth9g7oI0jyytZ/+RRs4c2pfKihTXnTOMgjw3aEqSuo62uM73t2OMfxJC6BljPNzmCTuA5VvKDIeONfLE8s1UVdewdut++uTncMu0Uu6YlWJEYa+k40mS9LG1RfleHmOc3PJEU9s8YQewfEuZJcbIkvRuqqrT/OKtLTQ2Ry4cM5jKWSkuGVdEdpZLUiRJnVNblO//A5wHlAF/A6wAVsUYj7Zl0PZk+ZYy17b9R3hwUR33L0zTsO8oJQN6cvvMFJ+dXsrAXnlJx5MkqVXaZM13COEM4GXgJ8BZwATgGCdK+GfbKGu7sXxLme94UzMvrGmgqrqGBRt2kZeTxTWThlJZUc45pf2TjidJ0mlpsw2XIYQxMcZ1Jz3uDUyMMS74+DHbl+Vb6lzWNexnbnWax5bVc/BYE5NK+lFZUc41k4aSn5uddDxJkj6UVzvB8i11VvuPHOfx5Zuoqk6zftsBBhTkcsv0Uu6YmaJ0YEHS8SRJ+h2WbyzfUmcXY6R6w07mVqd5fk0DzTFy6dgi5lSkuHD0YLLcoClJyhCnKt8f6eK6IYRPxxif+nixJOn0hRA494xCzj2jkC17D/PAwlrmLarjxR8tpnxQAXfMSnHz1FL6FeQmHVWSpA/1kWa+QwgrY4yT2iFPu3HmW+p6jjU284tVW5hbnWZJejf5uVlcd/Zw5lSkmDi8X9LxJEndVJvPfAP+fVdS4vJysrjunOFcd85w1mzex9wFNTyxfDMPLaljSll/7jy3nKsmDqFHjhs0JUmZwZlvSV3K3sPHeWRpPfctSLNxx0EKe+dx6/QyZs8sY1j/nknHkyR1A22+4dLyLSnTNTdHXlu/g6rqGl5cu40AXDG+mMqKcs49YxAh+Ac8SVL7aI9lJ5KU0bKyAheOGcyFYwZTt+sQ9y+s5aHFtTy3uoEzBvdizqwUn5laQp98N2hKkjrOR535fiHGeEU75Gk3znxLOnK8iWdWbqFqQZoVdXsoyMvmhsnDqawoZ+yQPknHkyR1EV3yOt8hhDOBPwEKgRdjjP92qve3fEs62cr6PVRVp3lyxWaONTYzY8RA7qwo5xMTisnNzko6niSpE+s05TuEcC9wDbAtxjjxpONXAd8GsoF7YozfPOltWUBVjPGOUz235VvSB9l18Bg/XVLHfQvS1O8+THHfHtw2o4zZM8oo6pufdDxJUifUmcr3hcABTpTpiS3HsoF1wBVAPbAYuC3GuCaEcC3wVWBujHHeqZ7b8i3pVJqaI796ZxtV1WleWbednKzAlROHUDkrxYwRA92gKUk6be1xh8tewJEYY9PHSvY+McZXQwjl7zs8A1gfY9zQ8rkfBK4D1sQYnwSeDCE8A5yyfEvSqWRnBS47s5jLziymZsdB7luQ5qdL6nhm5RbGDenDHbNS3DB5OL16uE9dkvTRndbMd8vSjluB24HpwFGgB7ADeAb49xjj+jYJdKJ8P33SzPdNwFUxxi+2PJ4DzAQeAW5sybEyxvi9D3iuu4C7AMrKyqam0+m2iCipmzh8rIknV2yiqjrN6s376NMjh89MLeGOWSlGFfVOOp4kKUN97GUnIYRXgF8CPwNWxRibW44PBC4BZgOPxxjva4Ow5ZxG+Y4xfq01z+uyE0kfVYyRZbV7mFtdwzNvbeF4U+S8UYOorCjnsnFF5LhBU5J0krZYdnJNjHH/+w/GGHcBjwKPhhDa62K5m4DSkx6XtByTpA4RQmBqagBTUwP460+N56HFtdy/sJYvz13KsH753D4rxWenl1LYu0fSUSVJGe50p2v+6tevhBDO+6B3iDEeb5NEv2sxMDqEMCKEkMeJ5S9PttPnkqRTGtynB1+7dDTz//wS7r5jKiMG9+J/P/cO5/7DS/zpg8tZmt5NJm1klyRlltOd+X6uZf30QeBM4PX2CBNCeAC4GCgMIdQD34gx/jCE8DXgOU5cavDeGOPq9vj8knS6crKzuGriEK6aOIT12w5w34I0jy6t54k3NzNhWF/urCjn02cPo2dedtJRJUkZ5Peu+Q4h/BDYC0wGFsQY/7IjgrU113xLam8Hjzby+PJNzK1O807Dfvr1zOWWaSc2aKYG9Uo6niSpg7TFhsuewBROXOnkzBjjl9s2YvuzfEvqKDFGFm3cRdWCNM+t2kpTjFw0ZjCVFSkuGlNEdpbXDJekrqzT3GSnPVm+JSWhYd8R5i2s5YFFtWzbf5TSgT25Y2aKW6aVMqBXXtLxJEntoM3Kdwjh2zHGP2mzZB3I8i0pScebmnlu9VaqqtMs2riLHjlZfPrsYVRWpJhU0j/peJKkNtSWd7jcH0J4Crg1xngwhHAl8Lcxxg+8Aook6YTc7CyumTSMayYNY+3WfcytTvP48k08srSec0r7U1mR4uqzhpKf6wZNSerKWr3sJIQwG/jPwDHgAPA/Yozz2yFbm3LmW1Km2XfkOI8tradqQZoN2w8ysFcen51eyu0zyygZUJB0PEnSR9SWy04uA/4bEIChwLUxxnfaJGU7s3xLylQxRt54bydV1TW8sKYBgEvHFVNZkeL8UYVkuUFTkjqVtlx28tfA38QYXwshnAU8FEL4eozxpY+dUpK6qRAC540q5LxRhWzac5h5C9M8uKiOX77dwIjCXtwxK8VNU0vo17O9biQsSeooH+tqJyGEocCjMcZz2y5S+3DmW1JncrSxiV+8tZWq6hqW1e6hZ242108expxZ5Ywf1jfpeJKkU2iL63yH+CHvGELoGWM8fKr3yQSWb0md1apNe5lbneZnKzZx5Hgz08sHMKeinKsmDCEvJyvpeJKk92mL8v0y8Bjwsxhj7UnH84DzgTuBl2OMP26TxO3A8i2ps9tz6BgPL6nnvoVp0jsPUdi7B7NnlDJ7Zooh/fKTjidJatEW5fv7wJvAFzix0XIPkA9kA88D/zfGuLzNErcDy7ekrqK5OfLKu9uZW53m5Xe2kRUCnxhfzJyKFBUjBxGCGzQlKUltUb6XxxgnhxCWATOAwcDhGOOeto3afizfkrqi2p2HuH9hmoeW1LHn0HFGF/VmTkWKG6eU0LtHa/fUS5LaQluU7/8DnAeUAX8DrABWxRiPtmXQ9mT5ltSVHTnexFMrNjN3QZqV9XvplZfNjVNKqKxIMbq4T9LxJKlbaZPrfIcQzgBeBn4CnAVM4MSNdlbFGD/bRlnbjeVbUnfxZt0eqqpreHrFFo41NVMxchCVFSmuGF9MTrYbNCWpvbXlTXbGxBjXnfS4NzAxxrjg48dsX5ZvSd3NzgNHeWhJHfcvqGXTnsMM6ZvP7Jll3DqjlKI+btCUpPbSFstOyk7zc+2JMe5rTbiOYvmW1F01NUdeWruNquoa5r+7g9zswFUTh1JZkWJaaoAbNCWpjbXFHS5/AkRO3Fb+w0Tgx0BVq9JJktpVdlbgivHFXDG+mA3bD3DfgloeXlrHUys2M25IHyoryrl+8jAK8tygKUnt7WPd4bIzceZbkn7j0LFGfvbmZqqq07y9ZR998nO4eWopcypSjCjslXQ8SerU2mzNd2dm+Zak3xVjZGl6N1XVaX6xagvHmyIXjC6ksqKcS8cVkZ3lkhRJai3LN5ZvSfp9tu0/woOL6pi3sJat+44wvH9Pbp9VxmenlTKod4+k40lSp2H5xvItSaersamZF9Y0UFWdpnrDTvJysrjmrKFUnlvOOaX9k44nSRnP8o3lW5I+incb9jN3QZpHl9Zz8FgTk0r6MWdWik+fPYz83Oyk40lSRrJ8Y4w+bIYAAB73SURBVPmWpI/jwNFGHl9WT1V1mne3HaB/QS6fnVbKHbNSlA4sSDqeJGUUyzeWb0lqCzFGFmzYxdwFNTy3uoHmGLlkbBFzKlJcNHowWW7QlCTLN1i+Jamtbdl7mAcW1jJvUR07DhwlNaiAO2amuHlaCf0L8pKOJ0mJsXxj+Zak9nKssZlnV29lbnUNi2t2k5+bxbVnD6OyopyJw/slHU+SOpzlG8u3JHWENZv3MXdBmieWb+Lw8SamlPWnsqKcT541hB45btCU1D1YvrF8S1JH2nv4OI8sree+BWk27jjIoF553DqjlNkzUwzv3zPpeJLUrizfWL4lKQnNzZHX1u+gqjrNS2sbALj8zGIqK8o5b9QgQnCDpqSu51TlO6ejw0iSuo+srMCFYwZz4ZjB1O8+xP0La3locR3Pr2lg5OBezJmV4jNTS+ibn5t0VEnqEM58S5I61JHjTfz8rS1UVad5s24PBXnZXD95OJUVKcYN6Zt0PEn62Fx2guVbkjLRW/V7qaqu4ckVmzna2MyMEQOprEhx5YQh5GZnJR1Pkj4SyzeWb0nKZLsPHuOnS+q4b2Gaul2HKerTg9tmlDF7ZhnFffOTjidJrWL5xvItSZ1BU3PklXXbqKpO88q67WSHwJUThjCnIsXMEQPdoCmpU3DDpSSpU8jOClw6rphLxxWT3nmQ+xak+emSep55awtji/twR0WKGycPp1cPf31J6pyc+ZYkZbTDx5p4asVmqhbUsGrTPnr3yOEzU4Yzp6KcUUW9k44nSb/DZSdYviWps4sxsrxuD1Vv1PDzt7ZyrKmZ80YNYs6sci4/s4gcN2hKyhCWbyzfktSV7DhwlIcW13H/gjSb9x5hWL98Zs8s49YZZRT27pF0PEndnOUby7ckdUWNTc28uHYbc6vTvLZ+B7nZgavPGkplRTlTyvq7QVNSItxwKUnqknKys7hywhCunDCE97YfYG51mkeX1vOzNzczYVhfKitSXHv2cHrmZScdVZIAZ74lSV3MwaONPPHmJuZWp1m7dT/9euZy89QS7piVorywV9LxJHUDLjvB8i1J3U2MkcU1u/lJdQ3PrdpKY3PkojGDqaxIcfHYIrKzXJIiqX247ESS1O2EEJgxYiAzRgxk274jzFtUy7yFtXzhJ0soHdiT22em+Oy0Ugb0yks6qqRuxJlvSVK3cbypmedXN1BVXcPCjbvIy8ni2rOHUVmRYlJJ/6TjSeoiXHaC5VuS9Nve2bqfuQtqeGzZJg4da+Ls0v5UzkrxqUlDyc91g6akj87yjeVbkvTB9h85zmPLNlFVXcN72w8ysFcet0wr5faZZZQOLEg6nqROyPKN5VuSdGoxRqrf28lPqmt4YU0DEbhsXBFzKsq5YFQhWW7QlHSa3HApSdLvEULg3FGFnDuqkM17DjNvYS0PLq7ll28vYkRhL26fWcbNU0vpV5CbdFRJnZgz35IkfYijjU08u2orVdVplqZ3k5+bxQ2ThzNnVjnjh/VNOp6kDOWyEyzfkqSPZ9WmvcytTvOzFZs4cryZaakBzKlI8cmJQ8nLyUo6nqQMYvnG8i1Jaht7Dx3n4aV1zF2QJr3zEIW9e3DbjFJmzyxjaL+eSceTlAEs31i+JUltq7k58uq725lbneald7aRFQJXnFlMZUWKijMGEYIbNKXuqktuuAwhjAT+GugXY7wp6TySpO4lKytw8dgiLh5bRN2uQ9y3MM1PF9fx7OqtjCrqzZxZKW6cMpw++W7QlPQbGTXzHUK4F7gG2BZjnHjS8auAbwPZwD0xxm+e9LZHTqd8O/MtSWpvR4438fTKLcytrmFF/V565WVz45QS5lSkGFPcJ+l4kjpIp1l2EkK4EDgAVP26fIcQsoF1wBVAPbAYuC3GuKbl7ZZvSVLGebNuD1XVNTy9cgvHGpuZNXIglRXlXDG+mNxsN2hKXVmnWXYSY3w1hFD+vsMzgPUxxg0AIYQHgeuANR2bTpKk03dOaX/OKT2H//ap8Ty0uI77FqT5w/uXUdy3B7NnpLhtRilFffOTjimpg3WG/3oPB+pOelwPDA8hDAoh3A1MDiH85Qd9YAjhrhDCkhDCku3bt3dEVkmSfsvAXnl89eIzePXPL+GeymmMHdKXf/7lOs795kt8bd4yFm3cRSb9FVpS+8qome/WiDHuBL7ye97n+8D34cSyk47IJUnSB8nOClw+vpjLxxezccdB7luQ5uEldTy9cgvjhvShsqKc6ycPoyCv0/5qlnQaOsPM9yag9KTHJS3HJEnqlEYU9uJvrhnPgr+6jG/eeBZZIfBXj7/FzP/5Iv/9qdVs2H4g6YiS2kln+O/1YmB0CGEEJ0r3rcDsZCNJkvTxFeTlcOuMMj47vZSl6d1UVae5b0GaH71ewwWjC5kzK8VlZxaTneU1w6WuIqPKdwjhAeBioDCEUA98I8b4wxDC14DnOHGpwXtjjKsTjClJUpsKITCtfCDTygeyff94HlxUy7xFtdw1dynD+/dk9swybp1eyqDePZKOKuljyqhLDbYnLzUoSepMGpua+eXbDVRVp3njvZ3kZWdxzaShzKlIcU5pf++gKWWwTnOd7/Zk+ZYkdVbrt+1nbnWaR5dt4sDRRs4a3o85FSmuPXsY+bnZSceT9D6WbyzfkqTO78DRRh5fvom51TWsazhA/4JcbplWyh0zU5QNKkg6nqQWlm8s35KkriPGyIINu5i7oIbnVjfQHCMXjxlMZUU5F40ZTJYbNKVEWb6xfEuSuqate48wb1EtDyyqZfv+o5QNLOCOWWXcMq2U/gV5SceTuiXLN5ZvSVLXdqyxmedWb2VudZpFNbvokZPFdecMo7KinInD+yUdT+pWLN9YviVJ3cfbW/ZRVZ3mieWbOHy8icll/amsSHH1WUPpkeMGTam9Wb6xfEuSup+9h4/z6NJ67luQZsOOgwzqlcdnp5dy+6wUw/v3TDqe1GVZvrF8S5K6r+bmyOvv7aCqOs2LbzcAcNmZxVRWpDh/VKHXDJfa2KnKd0bd4VKSJLW9rKzABaMHc8HowdTvPsS8hbU8uLiOF9Y0MLKwF3fMSnHTtBL65ucmHVXq8pz5liSpGzra2MTP39pCVXWa5bV7KMjL5vrJw6msSDFuSN+k40mdmstOsHxLkvRh3qrfS1V1DU+u2MzRxmZmlA9kTkWKKycMIS8nK+l4Uqdj+cbyLUnS77P74DEeXlrHfQtqqd11iMF9enDbjDJmzyhjSL/8pONJnYblG8u3JEmnq7k58sq67VRV1/CrddvJCoErJxQzZ1Y5s0YOdIOm9Hu44VKSJJ22rKzAJeOKuGRcEemdB7l/YS0PLa7j529tZUxxb+bMSnHDlBJ697BGSK3lzLckSfq9jhxv4skVm6mqrmHVpn307pHDZ6YMZ05FilFFfZKOJ2UUl51g+ZYkqS3EGFlet4e51WmeWbmFY03NnHvGICorUlx+ZjE52W7QlCzfWL4lSWprOw4c5aHFdcxbWMumPYcZ2i+f2TPKuHVGGYP79Eg6npQYyzeWb0mS2ktTc+TFtxuYuyDN/Hd3kJsduPqsoVRWpJhSNsANmup23HApSZLaTXZW4BMThvCJCUN4b/sB7luQ5pGl9fzszc2MH9qXyooU150znJ552UlHlRLnzLckSWpzh4418sTyExs0127dT9/8HG6eVsqcWSnKC3slHU9qVy47wfItSVISYowsrtlNVXUNz67aSmNz5MIxg6mcleKScUVkZ7kkRV2Py04kSVIiQgjMGDGQGSMGsm3fER5YVMe8RWm+WLWEkgE9uWNWilumlTKwV17SUaUO4cy3JEnqUMebmnlhTQNV1TUs2LCLvJwsPj1pGJUVKc4u7Z90POljc9kJlm9JkjLRuob9zK1O89iyeg4ea+Lskn7MqSjnmklDyc91g6Y6J8s3lm9JkjLZ/iPHeWzZJqqqa3hv+0EGFORyy/RS7piZonRgQdLxpFaxfGP5liSpM4gxUv3eTqqq07zwdgPNMXLp2CLmVKS4cPRgstygqU7ADZeSJKlTCCFw7qhCzh1VyJa9h5m3sJYHFtXy4o+2UT6ogDtmpbh5ain9CnKTjip9JM58S5KkjHassZlfrNpCVXWapend5Odmcf05w5lTkWLCsH5Jx5N+h8tOsHxLktQVrN68l7nVaZ54cxNHjjczNTWAyooUn5w4lLycrKTjSYDlG7B8S5LUlew9dJyHl9Zx34I0NTsPUdg7j1unlzF7ZhnD+vdMOp66Ocs3lm9Jkrqi5ubI/PU7mFtdw4trt5EVApefWURlRTnnnjGIENygqY7nhktJktQlZWUFLhozmIvGDKZu1yHuX1jLQ4treW51A2cM7sWcWSlumFJCv55u0FRmcOZbkiR1KUeON/HMyi1ULUizom4PPXOzufbsYdw+q4xJJd5BU+3PZSdYviVJ6o7eqt/LvEVpnli+mcPHmzhreD/umFXGp88eRkGeCwDUPizfWL4lSerO9h05zhPLN3HfgjTrGg7Qp0cON04Zzu2zUowp7pN0PHUxlm8s35Ik6cQdNJekd3P/gjQ/f2srx5qamVE+kNtnlXHVxCH0yMlOOqK6AMs3lm9JkvTbdh08xsNL6pi3qJb0zkMM7JXHzdNKuH1GirJBBUnHUydm+cbyLUmSPlhzc+T193Zw/4JaXni7gabmyIVjBnP7zDIuG1dETrY371HrWL6xfEuSpN9v694jPLS4jgcW1bJ13xGG9M3n1hml3Dq9jCH98pOOp07C8o3lW5Iknb7GpmZeWruN+xfW8uq72//j5j23z0xx/qhCsrK8eY8+nDfZkSRJaoWc7Cw+MWEIn5gwhNqdh5i3qJaHl9Tx3OoGygYWMHtmGTdPLWFQ7x5JR1Un48y3JEnSaTja2MSzq7Zy/8JaFm3cRV52Fp88awi3z0wxvXyAt7LXf3DZCZZvSZLUdt5t2M/9C2t5dFk9+480Mqa4N7fPTHHDlOH0zfdW9t2d5RvLtyRJanuHjjXy9Iot3Lcwzcr6vfTMzeaaSUO5dUYZU8r6OxveTVm+sXxLkqT2tbJ+Dw8squPJNzdx8FgTY4p7c9uMMm6YPJz+BXlJx1MHsnxj+ZYkSR3jwNFGnl6xmQcW1bKifi95OVl86qyh3Dq9lBkjBjob3g1YvrF8S5Kkjrd6814eXFTHE8s3sf9oIyMH9+K26WXcOGW4V0rpwizfWL4lSVJyDh9r4pm3tvDAolqWpneTmx24csIQbptRRsXIQV43vIuxfGP5liRJmWFdw34eWFTLY8s2sffwcVKDCvjs9FJumlpCUR/votkVWL6xfEuSpMxy5HgTz63eyryFtSzcuIucrMBlZxZx24wyLhg9mGxnwzstyzeWb0mSlLk2bD/AQ4vreHhpPbsOHmN4/57cMq2UW6aXMLRfz6TjqZUs31i+JUlS5jvW2MwLaxp4YFEtr63fQQhwwejB3DKthMvPLCY/NzvpiDoNXbJ8hxBGAn8N9Isx3vT73t/yLUmSOpPanYd4ZGkdjyytZ/PeI/Trmct15wzjlmmlTBjW10sWZrCMK98hhHuBa4BtMcaJJx2/Cvg2kA3cE2P85mk81yOWb0mS1FU1NUfeeG8HDy+p59nVWznW2My4IX24eVop158zzEsWZqBMLN8XAgeAql+X7xBCNrAOuAKoBxYDt3GiiP/D+57i8zHGbS0fZ/mWJEndwt5Dx3ly5WYeXlLHyvq95GYHLhtXzC3TS7hw9GBysrOSjihOXb5zOjoMQIzx1RBC+fsOzwDWxxg3AIQQHgSuizH+AydmyVsthHAXcBdAWVnZR84rSZKUCfoV5DJnVoo5s1K8s3U/Dy+p4/Hlm3h29VaK+vTgxikl3DythDMG9046qj5EYmu+W8r30yfNfN8EXBVj/GLL4znAzBjj1z7k4wcBf8+JmfJ7Wkr6h3LmW5IkdUXHGpt5+Z1tPLykjpff2U5Tc2RKWX9unlbK1WcNpV/P3KQjdjsZN/PdFmKMO4GvJJ1DkiQpSXk5WVw5YQhXThjCtv1HeHzZJh5eWs9fPvYW33hyNVecWcwNk4dz0djB5LosJXGZVL43AaUnPS5pOSZJkqTTUNQnny9fdAZ3XTiSFfV7eWL5Jp5csZln3trCwF55fHrSUG6YUsLZJf28WkpCMmnZSQ4nNlxexonSvRiYHWNc3Rafz2UnkiSpOzre1Myr67bz2PJNvLCmgWONzYws7MUNk4dz/eThlA4sSDpil5OJVzt5ALgYKAQagG/EGH8YQrga+BdOXOHk3hjj37fV57R8S5Kk7m7fkeP84q0tPLZsEws37gJgevkAbphcwqfOGkq/AteHt4WMK99JsHxLkiT9Rv3uQ/zszc08tqye97YfJC87i0vHFXHN2UO5bFwxPfO8m+ZHZfnG8i1JkvRBYoys2rSPx5bX8/TKLWzff5SCvGwuP7OYayYN5aKxg+mRYxFvDcs3lm9JkqTfp6k5smjjLp5auZlfvLWF3YeO0yc/hysnDOGaSUM5b1ShV0w5DZZvLN+SJEmtcbypmTfe28lTKzbz3Oqt7D/SyICCXK6aOJRPnz2UmSMGkZ3lFVM+iOUby7ckSdJHdbSxiVfX7eDplZt5YU0Dh441Udg7jyvGD+GqiUOoGDmIvBxnxH/N8o3lW5IkqS0cPtbES2u38YtVW3h57TYOHmuiT34Ol40r4qqJQ7hwzGAK8jLpVjIdz/KN5VuSJKmtHTnexBvv7eDZVVt5YU0Duw8dJz83iwtHD+aqiUO4bFxxt7x8YZe8vbwkSZKSlZ+bzaXjirl0XDGNTc0sqtnF86sbeHbVVp5f00BOVmDWyEFcOq6IS8cVUV7YK+nIiXPmW5IkSW2quTmyctNenl21lV++3cD6bQcAGFnY6z+K+LTygV12nbjLTrB8S5IkJaV25yFeWtvAi2u3sXDDLo41NdO7Rw4XjinkkrFFXDy2iMF9eiQds81YvrF8S5IkZYKDRxt5ff0OXlq7jZfWbmPb/qMAjB/alwtGF3L+6EKmlw8kP7fz3tjH8o3lW5IkKdPEGFm9eR+/emcb89/dwbLa3RxviuTlZDG9fADnjxrMBaMLGT+0L1md6Jrilm8s35IkSZnu4NFGFtXs4rV3d/Dauzt4p2E/AAN75VFxxiBmjRjI9BEDGVPUJ6PLuFc7kSRJUsbr1SOHS8YWccnYIgC27TvC6+/tYP67O3hj/U6eWbkFgH49c5lePpCZLWV8wrC+nea29858S5IkKePFGKnffZiFG3exeOMuFtXsYuOOgwAU5GUzNTWAqakBTCkbwNml/enXM7nrizvzLUmSpE4thEDpwAJKBxZw09QS4MTM+OKa3SzauJOFG3fx7RffJUYIAUYN7s3XrxjDJ88amnDy32b5liRJUqdU1DefT00ayqcmnSjY+48cZ2X9Xpald7OsdjcFPTKv6mZeIkmSJOkj6JOfy3mjCjlvVGHSUT5U51iZLkmSJHUBlm9JkiSpg1i+JUmSpA5i+ZYkSZI6iOVbkiRJ6iCWb0mSJKmDWL4lSZKkDmL5liRJkjqI5VuSJEnqIJZvSZIkqYNYviVJkqQOYvmWJEmSOojlW5IkSeogIcaYdIYOEULYDqQT+vSFwI6EPndn5Zi1nmPWeo5Z6zlmreeYtZ5j1nqOWeu155ilYoyDP+gN3aZ8JymEsCTGOC3pHJ2JY9Z6jlnrOWat55i1nmPWeo5Z6zlmrZfUmLnsRJIkSeoglm9JkiSpg1i+O8b3kw7QCTlmreeYtZ5j1nqOWes5Zq3nmLWeY9Z6iYyZa74lSZKkDuLMtyRJktRBLN/tKIRwVQjhnRDC+hDCXySdJ1OEEEpDCC+HENaEEFaHEP6k5fjfhRA2hRDebHm5+qSP+cuWcXwnhHBlcumTE0KoCSG81TI2S1qODQwhvBBCeLfl3wEtx0MI4TstY7YyhDAl2fQdL4Qw9qRz6c0Qwr4Qwp96nv22EMK9IYRtIYRVJx1r9XkVQriz5f3fDSHcmcTX0lE+ZMz+dwhhbcu4PB5C6N9yvDyEcPik8+3ukz5masv39PqWcQ1JfD0d4UPGrNXfi93p9+qHjNlDJ41XTQjhzZbjnmecsl9k1s+0GKMv7fACZAPvASOBPGAFMD7pXJnwAgwFprS83gdYB4wH/g74Lx/w/uNbxq8HMKJlXLOT/joSGLcaoPB9x/4X8Bctr/8F8I8tr18N/AIIwCxgYdL5Ex67bGArkPI8+52v+0JgCrDqo55XwEBgQ8u/A1peH5D019bBY/YJIKfl9X88aczKT36/9z3PopZxDC3j+smkv7YOHrNWfS92t9+rHzRm73v7t4C/9Tz7ra/1w/pFRv1Mc+a7/cwA1scYN8QYjwEPAtclnCkjxBi3xBiXtby+H3gbGH6KD7kOeDDGeDTGuBFYz4nx1Ymx+UnL6z8Brj/peFU8YQHQP4QwNImAGeIy4L0Y46lutNUtz7MY46vArvcdbu15dSXwQoxxV4xxN/ACcFX7p0/GB41ZjPH5GGNjy8MFQMmpnqNl3PrGGBfEE7/tq/jNOHc5H3KefZgP+17sVr9XTzVmLbPXtwAPnOo5uuF59mH9IqN+plm+289woO6kx/WcumB2SyGEcmAysLDl0Nda/vRz76//LIRj+WsReD6EsDSEcFfLseIY45aW17cCxS2vO2a/7VZ++5eU59mptfa8cux+2+c5MZv2ayNCCMtDCK+EEC5oOTacE+P0a911zFrzveh59hsXAA0xxndPOuZ5dpL39YuM+plm+VZiQgi9gUeBP40x7gP+DTgDOAfYwok/qek3zo8xTgE+CfxRCOHCk9/YMqvh5YveJ4SQB1wLPNxyyPOsFTyvWieE8NdAI3B/y6EtQFmMcTLwdWBeCKFvUvkyjN+LH91t/PaEgufZST6gX/yHTPiZZvluP5uA0pMel7QcExBCyOXEN8b9McbHAGKMDTHGphhjM/ADfvMnf8cSiDFuavl3G/A4J8an4dfLSVr+3dby7o7Zb3wSWBZjbADPs9PU2vPKsQNCCH8AXAPc3vILnpalEztbXl/KiTXLYzgxPicvTel2Y/YRvhc9z4AQQg5wI/DQr495nv3GB/ULMuxnmuW7/SwGRocQRrTMvN0KPJlwpozQslbth8DbMcZ/Oun4yWuSb/h/7d3vi1VFHMfx9wcJExHphw8iAhOsoCgLBS0LH4iEWKQFRoH9AjOwIAiR/AcEHwUF0SMhxB4U5T7KSrDMCCXbXX+kpAUSUoSFlkKaTA/ObN3d1l1X5e6P+37BgTlzZ+acM8w998u9M/cAfSu8u4Ank0xOciswm2YBScdIMjXJtL40zeKuAzR907cK+xlgW013AavqSu75wKmWn9w6Tb9viBxnl2Sk42o7sCTJdXXqwJKa1zGSPAysAx4tpZxtyZ+RZFJNz6IZVz/UfjudZH69J67iv37uCJfxXvRztbEYOFxK+Xc6ieOscbH4grF2T7taKzfdBl11u5Rmpe0xYMNon89Y2YCFND/59ALddVsKvAvsr/ldwE0tdTbUfjzCBF6pPUSfzaJZ2d8DHOwbT8ANwA7ge+Az4PqaH+Ct2mf7gbmjfQ2j1G9TgZPA9JY8x1n/PtpK85P1eZp5jS9czriimed8tG7PjfZ1jUKfHaWZI9p3T3u7ln28vme7gX3AIy3tzKUJOI8Bb1IffDcRt4v02Yjfi530uTpYn9X8zcCaAWUdZ2XI+GJM3dN8wqUkSZLUJk47kSRJktrE4FuSJElqE4NvSZIkqU0MviVJkqQ2MfiWJEmS2sTgW5LGmSQzkxwYvmS/Ol+11H3qKp/P64MdS5L0fwbfktQBSin31+RMYETBd32i3lD6Bd8tx5IkDWDwLUnjWJJZSb5NMm+Ycn/W5EbgwSTdSV5NMinJpiR7k/QmebGWX5RkV5Iu4FDN+yjJN0kOJlld8zYCU2p7W1qPVZ8atynJgST7k6xsaXtnkveTHE6ypT6ZTpImvOG+zZAkjVFJbgfeA54tpfRcYrX1wGullGW1jdU0j1Sel2QysDvJJ7XsfcBdpZQf6/7zpZTfkkwB9ib5oJSyPsnaUsqcQY61ApgD3APcWOt8UV+7F7gTOAHsBh4AvhzB5UvSuGTwLUnj0wxgG7CilHLoCtpZAtyd5Im6Px2YDZwD9rQE3gCvJFle07fUcieHaHshsLWUcgH4JcnnwDzgdG37J4Ak3TTTYQy+JU14Bt+SND6dAo7TBLhXEnwHeLmUsr1fZrIIODNgfzGwoJRyNslO4NorOO5fLekL+HkkqUM451uSxqdzwHJgVd+/lyS5OcmOYer9AUxr2d8OvJTkmtrGbUmmDlJvOvB7DbzvAOa3vHa+r/4Au4CVdV75DOAhYM+lXJwkTVR+0yBJ41Qp5UySZcCndZHjCeDvYar1AheS9ACbgTdopnzsq4sefwUeG6Tex8CaJN8BR4CvW157B+hNsq+U8nRL/ofAAqAHKMC6UsrPNXiXpI6UUspon4Mk6SpIshY4XkrpGu1zkSQNzuBbkiRJahPnfEuSJEltYvAtSZIktYnBtyRJktQmBt+SJElSmxh8S5IkSW1i8C1JkiS1icG3JEmS1Cb/AAVu/EKxrxDdAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(12,6))\n",
    "\n",
    "plt.semilogy(np.abs(f_da-f_opt))\n",
    "\n",
    "plt.ylabel(\"$|f(x^{k}) - f^*|$\")\n",
    "plt.xlabel(\"k, iteration\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that all the constraints are respected:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sum(x_i) =\t\t 100.00000 \n",
      "# negatives x_i =\t 0\n"
     ]
    }
   ],
   "source": [
    "print(\"sum(x_i) =\\t\\t %.5f\" % np.sum(x_da), \"\\n# negatives x_i =\\t\",np.sum(x_da<0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analysis\n",
    "\n",
    "In the plot above, we report $|f(x^{k}) - f^*|$, that is the distance of the value of the function from the optimal computed before using ```scipy.optimize```.\n",
    "Notice that it is not meaningful to plot  $f(x^{k})$. The step of the dual ascent is not always a descent on $f(x)$ since it takes care also of the constraints (the function can increase at some point if it doesn't respect the constraints).\n",
    "\n",
    "\n",
    "By a trial and error proceedure, we concluded that $t=0.15$ yields the best convergence. As we can notice from the plot, the convergence is exponential (linear in log scale).\n",
    "\n",
    "For stopping condition, we can no longer use the check on the norm of the gradient (as we did in Problem 1), since now this may not converge to 0. Indeed for the gradient ascent in the variable $u$,  $\\nabla g(u)=C x -d$ and this will tend to be $Cx -d \\le 0$ (since it's a constraint of our minimisation problem).\n",
    "\n",
    "A reasonable choice seems to check the difference of the function we want to minimise ($f(x)$) with respect its value in the previous iteration: if the function is not changing we stop the execution. Therefore the stopping condition is $|f(x_{k})- f(x_{k-1})|<\\epsilon$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3) Parallelisation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solution can be parallelised using dual decomposition. Indeed since the function to minimise can be decomposed in $m$ terms $f_i(x_i)$ that depends only on one variable $x_i$ (recall $f(x_1, x_2, ..., x_m)= \\sum_{i=1}^{m} a_i(x_i-b_i)^2 + 3$ ).\n",
    "\n",
    "This allows to split the argmin problem like follows:\n",
    "$$ x^{k+1}= argmin_{x}{(f(x) + (u ^k)^T C x)} \\iff x^{k+1}_i= argmin_{x_i}{(f_i(x_i) + (u ^k)^T C_i x_i)} $$\n",
    "\n",
    "In this way, we don't reduce the number of iterations required, because after each iteration we still need to gather all the $x_i$ and compute the updated $u^{k}$.\n",
    "\n",
    "Anyway this decomposition can speed up the execution time, since different cores/processors can run smaller argmin problems in parallel and so each iteration lasts less time."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "project-optim",
   "language": "python",
   "name": "project-optim"
  },
  "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}