ch04.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "第4章 追跡アルゴリズムの性能保証\n",
    "=="
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "定理4.1（OMPの最適解保証：二つの直交行列の場合）\n",
    "--\n",
    "$\\mathbf{\\Psi}$,$\\mathbf{\\Phi}$を$n×n$の直交行列とする．線形連立方程式$\\mathbf{Ax}=[\\mathbf{\\Psi},\\mathbf{\\Phi}]\\mathbf{x}=\\mathbf{b}$に対して，もし，解$\\mathbf{x}$が存在して、最初の$n$個の要素の中に$k_p$個の非ゼロ要素があり，\n",
    "$$max(k_p,k_q)<\\frac{1}{2\\mu(\\mathbf{A})}$$\n",
    "が成り立つならば，しきい値パラメータ$\\epsilon_0=0$で実行されるOMPは$k_0=k_p+K_q$回の反復で最適解を与える．"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "定理4.2（基底追跡の最適解保証：二つの直交行列の場合）\n",
    "--\n",
    "$\\mathbf{\\Psi}$,$\\mathbf{\\Phi}$を$n×n$の直交行列とする．線形連立方程式$\\mathbf{Ax}=[\\mathbf{\\Psi},\\mathbf{\\Phi}]\\mathbf{x}=\\mathbf{b}$に対して，もし，解$\\mathbf{x}$が存在して、最初の$n$個の要素の中に$k_p$個の非ゼロ要素があり，\n",
    "$$2\\mu(\\mathbf{A})^2k_pk_q+\\mu(\\mathbf{A})k_p-1<0$$\n",
    "が成り立つならば，得られた解は$(P_1)$の一意解であり，${P_0}$の一意解でもある．"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "線形計画問題(LP)\n",
    "--\n",
    "\n",
    "BP（基底追跡）の解$\\mathbf{x}$より$L_0$ノルムが大きくて$L_1$ノルムが小さい解$\\mathbf{y}$の集合\n",
    "\n",
    "$C = \\{\\mathbf{y}\\,|\\,\n",
    "\\mathbf{y}\\neq \\mathbf{x},\\,\n",
    "\\|\\mathbf{y}\\|_1\\leq\\|\\mathbf{x}\\|_1,\\,\n",
    "\\|\\mathbf{y}\\|_0\\geq\\|\\mathbf{x}\\|_0,\\,\n",
    "\\mathbf{A}(\\mathbf{y}-\\mathbf{x})=\\mathbf{0}\\}\n",
    "$\n",
    "\n",
    "が与えられた条件\n",
    "\n",
    "$2\\mu(\\mathbf{A})^2k_pk_q+\\mu(\\mathbf{A})k_p-1<0$<br>\n",
    "\n",
    "もしくはもっと解釈が容易な条件\n",
    "\n",
    "$\\|x\\|_0=k_p+k_q<\\frac{\\sqrt{2}-0.5}{\\mu(\\mathbf{A})}$\n",
    "\n",
    "の元で空となることを証明する"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "     fun: 50.0\n",
       " message: 'Optimization failed. Unable to find a feasible starting point.'\n",
       "     nit: 100\n",
       "  status: 2\n",
       " success: False\n",
       "       x: nan"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.optimize import linprog\n",
    "\n",
    "# パラメータ\n",
    "n, kp, kq, mu = 50, 9, 7, 0.1\n",
    "\n",
    "one_p = [1.] * kp + [0.] * (n - kp)\n",
    "one_q = [1.] * kq + [0.] * (n - kq)\n",
    "C = np.array(one_p + one_q)\n",
    "A_ub = np.diag([1. - mu] * n + [mu - 1.] * n)\n",
    "b_ub = np.zeros(2 * n)\n",
    "A_eq = np.eye(2 * n)\n",
    "b_eq = np.ones(2 * n)\n",
    "\n",
    "# 最大化したいので負をとる\n",
    "C *= -1\n",
    "A_ub *= -1\n",
    "b_ub *= -1\n",
    "A_eq *= -1\n",
    "b_eq *= -1\n",
    "\n",
    "linprog(C, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "nanだと…解集合が空であることを示したいからこれでいいのか…"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "定理4.3（OMPの最適解保証）\n",
    "--\n",
    "線形連立方程式$\\mathbf{Ax}=\\mathbf{b}$（$n<m$である$\\mathbf{A}\\in\\mathbb{R}^{n×m}$はフルランクの行列）に対して，もし，解$\\mathbf{x}$が存在して，\n",
    "$$\\|\\mathbf{x}\\|_0<\\frac{1}{2}\\left(1+\\frac{1}{\\mu(\\mathbf{A})}\\right)$$\n",
    "を満たせば，しきい値パラメータ$\\epsilon_0=0$で実行されるOMP(OGA)は最適解を与える．"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "定理4.4（しきい値アルゴリズムの最適解保証）\n",
    "--\n",
    "線形連立方程式$\\mathbf{Ax}=\\mathbf{b}$（$n<m$である$\\mathbf{A}\\in\\mathbb{R}^{n×m}$はフルランクの行列）に対して，解$\\mathbf{x}$が存在して<br>\n",
    "（その非ゼロ要素の最小値は$|x_{min}|$，最大値は$|x_{max}|$とする）\n",
    "$$\\|\\mathbf{x}\\|_0<\\frac{1}{2}\\left(1+\\frac{1}{\\mu(\\mathbf{A})}\\frac{|x_{min}|}{|x_{max}|}\\right)$$\n",
    "を満たせば，しきい値パラメータ$\\epsilon_0=0$で実行されるしきい値アルゴリズムは最適解を与える．"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "定理4.5（基底追跡の最適解保証）\n",
    "--\n",
    "線形連立方程式$\\mathbf{Ax}=\\mathbf{b}$（$n<m$である$\\mathbf{A}\\in\\mathbb{R}^{n×m}$はフルランクの行列）に対して，解$\\mathbf{x}$が存在して，\n",
    "$$\\|\\mathbf{x}\\|_0<\\frac{1}{2}\\left(1+\\frac{1}{\\mu(\\mathbf{A})}\\right)$$\n",
    "を満たせば，得られた解は$(P_1)$の一意解であり，$(P_0)$の一意解でもある．"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "驚くべきことに，OMPの上界は基底追跡の成功も保証する．ただし，これはいつでも二つのアルゴリズムは似たように振る舞うということではない．実際，二つの直交行列の場合について，これらの二つのアルゴリズムの違いをすでに議論している．一般の場合のこれらのアルゴリズムについての定理は，最悪の場合の振る舞いは，二つのアルゴリズムは同じであり，これらの上界はタイトであることを意味している．"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "しかし，上記の結果はかなり弱いものである．$n$個の要素のうち非ゼロ要素が$\\sqrt{n}$個よりも少ないという，かなりスパースなときにしか成功しない．それほどスパースな状態になる問題はほとんどない．"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Troppの厳密性復元条件\n",
    "==\n",
    "定義4.1（厳密性復元条件）\n",
    "--\n",
    "サポート$S$と行列$\\mathbf{A}$が与えられたとき，厳密復元条件（exact recovery condtion; ERC）は次式で与えられる．\n",
    "$$ERC(A,S):max_{i\\notin S}\\|\\mathbf{A}_{S}^{+}\\mathbf{a}_i\\|_{1}<1$$\n",
    "\n",
    "$\\mathbf{A}_{S}^{+}$は$\\mathbf{A}_{S}$の疑似逆行列"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "定理4.6（ERCと追跡アルゴリズムの性能）\n",
    "--\n",
    "線形連立方程式$\\mathbf{Ax}=\\mathbf{b}$の解である，サポート$S$を持つスパースなベクトル$\\mathbf{x}$に対して，もしERCが満たされれば，OMPとBPは$\\mathbf{x}$を求めることに成功する．"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "この条件は，相互コヒーレンスを用いた追跡アルゴリズムの成功条件よりも強いが，構成的ではない．つまり，サポートが与えられなければ確認することはできないのである．"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "（参考）Python SciPy: 線形計画問題を解く\n",
    "--\n",
    "http://org-technology.com/posts/scipy-linear-programming.html\n",
    "\n",
    "SciPy では関数 linprog で線形計画問題を解くことができます。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "     fun: -10.0\n",
       " message: 'Optimization terminated successfully.'\n",
       "     nit: 3\n",
       "   slack: array([ 0.,  0.,  1.])\n",
       "  status: 0\n",
       " success: True\n",
       "       x: array([ 2.,  2.])"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.optimize import linprog\n",
    "\n",
    "c = [-3, -2]\n",
    "A = [[2, 1], [1, 1]]\n",
    "b = [6, 4]\n",
    "x0_bounds = (None, None)\n",
    "x1_bounds = (1, None)\n",
    "res = linprog(c, A_ub=A, b_ub=b, bounds=(x0_bounds, x1_bounds))\n",
    "\n",
    "res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x114c78310>"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW8AAAEACAYAAAB8nvebAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4FGXXwOHfEyCEDlJFQIqgUiJSQhMIJRIFBGkiEOor\noiCCUpSO6Ov7gV1BBOldirQoEDoWEFEEFUWkqYCUAIEAqef7Y7JLKtkku9lNcu7r2kt2Znaes2M4\nTGbOc8aICEoppbIWL3cHoJRSKu00eSulVBakyVsppbIgTd5KKZUFafJWSqksyKHkbYzJG+/PuY0x\n+VwXklJKqdQYR0oFjTHbsBK9ALmAhSIyx8WxKaWUSkFuB7eLFZFWLo1EKaWUwxy95l3AGPOxMWan\nMeZDY0wBl0allFLqjlJN3saYgkBl4C0R8QeOAK+5OC6llFJ34Og17zwiEhX3ZwN8LSKNXR2cUkqp\n5Dl6zdvY/iAiYoy5lWQDY7RJilJKpYOImNS3SsiRyyYlga+NMUXi3vcG9ia37fPPC9WqCb/+Koh4\n5itoaBAlp5bk+3++d3ssd3pNnDjR7TFkhxg1To3T01/plWryFpELwBvAZmPMdsA/7n0S06fD6NHQ\nvDmsW5fumFyqcrHKzG4/m3bL2vHrhV/dHY5SSqWLQ5dNRGQtsNaRbfv3h5o14dixDMXlUh0e6EBY\nRBhtFrdhd9/dVCpWyd0hKaVUmjh6zTtN/Pyslyfy9/cHIOihIK5GXCVgUQB7+u3h7kJ3uzewRGxx\nerKsECNonM6mcXoGh6pNHNqRMeKsfWWm13e/zopfVrCr7y7uyneXu8NRSuUwxhgkHTcsMzV5//sv\nlC7tlOGcRkQYGTKSPaf3sDVoK4XyFnJ3SMrNKlasyKlTp9wdhspm7r33Xk6ePJlkuccn7+hoqFUL\nnnoKJkwALw/qZygiDNwwkONXjhPcIxif3D7uDkm5UdxfJneHobKZlH6uPD55A5w7B126QPHisHAh\nFCnilKGdIiY2hh5renAr+haruq4iT6487g5JuYkmb+UKzk7emXr+W6YMbN8O5cpBgwbw22+ZOfqd\n5fLKxaInFxEVE0X/9f2JlVh3h6SUUily2w3LuXNh8mQ4cgTy53dKCE5xI+oGgYsD8S3ty4ePfYjV\nDUDlJHrmrVwhS595x9e/Pxw65FmJGyB/nvxseHoD3/79LeN3jHd3OEqlSd26dfnmm29SXD9o0CCC\ng4OdOubGjRt57rnnnLpPgIkTJ1KtWjWaNGni9H3fyeTJk/n0008zdcz0yPGlgim5EH6BZvOb0b92\nf0Y2GenucFQmyspn3v7+/nz44YfUqlXLZWNs2LCBDz74gJCQEJeNERYWxn333cfZs2cREXLndsmU\nFABefvllihcvzpgxY1w2BmSjM29PV7JASUKCQpjx/QxmH5jt7nCUcogxhly5crls/zdv3qRAgQIU\nKVKE6OhoIiIiXDLGpUuXqFChAsYYoqKinD6GzY0bNyhYsCCFCxcmIiKCmJgYl43ldE5sriLOsHSp\nyNq1TtmVU/xx6Q8p+3ZZWX54ubtDUZnEWT/LrnL48GEJCAiQunXrSv369WXmzJn2df7+/vLFF19I\n27ZtpX79+uLn5yebN2+2r+/UqZOsXr1aRETOnTsnnTp1El9fX/Hz85O18f7iXbhwQbp16ya1a9eW\nOnXqyKeffioiIp07d5Zq1apJsWLFpHbt2vLJJ5/I+vXrpV27dnL58mUpUqSIhIWF2fezY8cOqV27\ntoiIHD16VFq3bi0PPfSQNGnSRL766qtkv9+QIUOkZs2aUqBAAalTp468+OKLEhYWJnnz5k2w3Qcf\nfCCDBg0SEZFhw4bJRx99JIMGDZKGDRvKww8/LCEhIfZtf/jhB2nWrJnUrl1bGjVqJPv27ZPw8HB5\n8MEHpWzZslKhQgWpWbOmHDx4UF566SV5++23RUTkxo0bMmzYMHn44YelXr160qdPH7l06ZKIiKxd\nu1a6d+8uH330kTRt2lR8fX1l6tSpKf5/S+nnKm552nNuej6U7I6c9AO/b59IuXIiEyaIxMQ4ZZcZ\ndujcISk1rZQEHw12dygqE3hy8g4LC5OqVavKgQMHREQkPDxcAgICZM2aNSJiJe/GjRvLP//8IyIi\nJ0+elMqVK8uZM2dERKR3796yceNGERFp1aqVPfGHhoaKn5+fHDlyREREAgMDZf78+SIicu3aNald\nu7Zs2rRJRER27twpXbt2tce0fft26dKli4iIDBgwQJYtW2ZfN3z4cPnkk08kKipKqlevLl988YWI\niPz1119Ss2ZNuXDhQrLf8+TJk1K/fn37+9jYWClevHiCbebOnSsvv/yyiIhMmDBBKlWqJBs2bBAR\nkb1798qDDz5oP2ZVqlSRw4cPi4jIkSNHpEyZMnLx4kUREZk0aZJMnz7dvt8JEybIRx99JCIigwcP\nlldffdW+btasWdK+fXv7977nnnvklVdeERGR69evS4kSJeT8+fPJfidnJ2+Pu2zi5wf791slhR07\nwtWr7o4IapWuxbru6+i7ti+7T+12dzjKA0yaBMYkfU2a5Nj2KW2XmuDgYAICAqhTpw4A+fPnZ8qU\nKcyeffvS3gsvvEDZsmUBa1bf008/zfr16wHw8vLCGMM///zDiRMnePbZZwEoVqwYzz33HEuXLuXc\nuXP89ddf9OnTB4CCBQuyfft2WrdunWxMtn0C9O/fnyVLlgAQGxvLhg0bePrpp9m3bx933303jz32\nGADlypWjc+fO9rhSY4zBK9HMPi8vL/syLy8v6tevT7t27QBo0KABly5dIiIiguDgYAIDA6lZsyYA\nDzzwAIcOHaJ48eJ3/D6xsbGsXr2a1167/eCwZ555hmPHjnH+/Hn72FOmTAGgQIEC+Pv7c/DgQYe+\nU0a57i5ABpQpA9u2wfDhVjJftw4eeMC9MTUs15BlnZfR5bMufNnzS+qWrevegJRbTZqUtgSc1u1T\ncvr0aSpXrpxgWeXKlTl9+rT9/X333Zdgfbly5Thz5kyCZf/88w+hoaE0btzYfhMtIiKC+vXrc/r0\naSpVSthps1ixYg7F17hxY06dOkVoaCgHDx6kcePGFCpUiH/++Yeff/6Zxo2tB3CJCDdu3KBr166O\nffFk2OK2adasWYL3BQoUICYmJtnvU7JkyVT3f/HiRUqUKJHkZmmlSpXsx7tRo0YJ1hcsWNCl1+jj\n88jkDeDtbfUHnzfPmlrvCVpVbsWs9rNou7Qt2/tsp3rJ6u4OSeUwFSpUYM+ePQmWHT9+nAoVKtjf\nHz16lHr16tnfnzx5MknCL1euHA8//DDbt2+3L4uNjcXLy4szZ85w4sSJBNtfuXIFHx8ffHxSbx0R\nFBTEihUr+P777xkwYIB9vMcff5y5c+cmGc9R0dHRCT6zc+dOSpUqZV+f0pyMcuXKJSmfPH/+PCVL\nlrzjPI4SJUpw6dIloqOjEyToEydOUKFCBY4cOeLWeSAed9kksX79rP7gnqLjAx2ZGjCVNovbcOLy\nidQ/oJQTtW3blpCQEPuv5tevX2fChAkMHDgQsBLi22+/zT///APAsWPHWLlyJR06dEiwn7Jly5I3\nb162bNkCQGRkJL169WLevHmULVuWe+65h3nz5tnXBQUFsWbNGgDy5MnDlStXAOtsPTY2NsFZcK9e\nvVi6dCmHDx+2nw37+fnx66+/2uMODw+nRYsW7NixI9nvGRsbm6Tyo1q1auzatQuALVu2sGfPHmJj\nY+3b2/5sExMTQ2xsLG3btmXTpk0cOnQIgLNnz9K0aVN78zHb94mJibH/AyEieHl50alTJyZPnmzf\n5+zZs6latSqlSpVKdszklrmKx555e7LeD/UmLCLMY3uBq+yrUKFCrFmzhuHDh3P58mWMMTzzzDN0\n7NgRsJLi6NGj6dWrF1evXsXb25u5c+dSOpl2ngsWLGDIkCGMGTOGmzdv0r17d/r27QvAkiVLGDJk\nCB9++CFeXl706tWLHj16AODr60toaCi+vr4899xzlCtXLkHJ4N13303evHkTXMbInTs3S5cuZciQ\nIfz7779ERkYyePBgWrRokez3vHnzJteuXUuw7K233uLFF1/E29ubBx98kHHjxvHtt98CVl14/kQz\n/sLDwwkLC6Ns2bKsWrWKoUOHcu3aNby9vfn444+pWLEiAIGBgXTp0oWVK1cye/ZswsPD7d9n2rRp\nvPLKK9StWxdjDDVr1mTBggX2McPDw5MdMzNk2Uk6585Z18bdSXuBZ09ZeZJOavr27Uu3bt14/PHH\n3R1KjqOTdIA//gBfX/c/J3Ns07EEVgnksSWPcS3iWuofUMoNYmJi+PLLLxER/vrrL4oWLerukJQT\nZMnkXbUqbNwIQ4bAxImQSZeYkjDGMDVgKg+VfoiOKzpyK/qWewJR6g7Cw8OZMWMGVapUIU+ePDRo\n0MDdISknyLKXTcC6dNK1KxQrBosWua8/uK0XeER0BCu7rtRe4Flcdr5sotxHL5vEY6sHL18e4uYT\nuIWtF3hkTKT2AldKZYosfeYdX1gYFC7stuEB7QWeXeiZt3IFPfNOgbsTN2gvcKVU5sk2ydtTFPEp\nwqaem1h9ZDXTvp7m7nCUUtlUtk3eItZT6t3xnEztBa6UcrVsm7yNgYoVoVkzcLBxmVOVK1yOLb22\nMGnXJFb8vCLzA1DZVmhoKD179sTX1xc/Pz969+5NaGiou8NKIiYmhhEjRuDr60uNGjUYO3as3ktw\nomybvMF6TuaGDTB4sNXRLbPrwasWr8qmnpsYumkoX/zxReYOrrKtnj17Ur16dQ4dOsR3331H8+bN\neeGFF9wdVhIvvPAC0dHR/PDDDxw8eJADBw4wY8YMd4eVbWSbapM7sdWDFy0KS5Zk/s3NvX/v5Yll\nT7Cq2yqa3dss9Q8ot/LkapPr169TvXp1jh49au/wFxoaSr169Th+/Libo7vt0qVLNG7cmJ9//pk8\neax5D7///juXL1+mYcOGbo7OPbTaJB1s9eCNGoELH++Xovi9wL8/833mB6CyjYIFC3L69OkErVkX\nLVqEr68vCxcuJDAw0L589OjRjByZ8OHZu3btol+/fgmWjRo1imHDhiVY1q9fPz744IMEy2JjY5kx\nY4b94QN38tVXX9GyZUt74ga4//77c2zidoUccebtKdb9to5BwYPY1nub9gL3YJ585p3YBx98wNtv\nv822bdu477776NixI0FBQdSoUYPu3buzb98+8ubNa99+165dLFiwIEFf7ePHjxMYGMjRo0cBK0lX\nqlSJAwcOUKJECQD27dvH0KFDadSoEa+99hpr1qzhk08+sR+r+P9dv349K1as4PLly3h7e7Nhwwai\noqJo1aoV48aNI1++fJl7kDyEs8+8tSVsJurwQAfCIsJos7gNu/vuplKxSql/SHkkMznjE7BkYvr/\ngbh58ybPPfccP//8M/v27aNMXIvNjz/+mJYtW1KiRAlmzpyZIHHbx02UQCpXrkyVKlXYs2cPTZs2\nZffu3dSpU8f+MIKxY8fyxx9/MHfuXGrUqAFY3Qlt7WOTc+PGDT799FPGjh3L7t27iYmJ4ZVXXmHE\niBFMnz493d9bxZOeB18m98KDH9p6J9evi1y5krljfrjvQ6nyfhU5E3YmcwdWDvH0n+WIiAhp3ry5\ntGnTRsLDw5Os7927t/2J7SLWw3wbNmwojRo1kurVq0upUqWkUaNG0qhRI/nkk09EROTzzz+X3r17\ni4jIwIED5fPPPxcRkRo1akjnzp3THONnn30mAQEBCZbdunVL7rnnnjTvK7tI6eeK7PL0+Mw2b55I\ntWoiv/6aueO+vut1qTmjply6cSlzB1ap8vSf5ZEjR0r79u0lKioqybq9e/dKgwYNpE+fPvLpp58m\nWb9z507p169fkuVRUVFSrVo1uXDhglSrVk0iIyNFROTKlSvy/PPPS7NmzeSHH36wbz9nzhypV6+e\n1K9f3/6yvT9//rz8/PPPEhgYmGCMK1euaPJOebkm7/SYM0ekZEmRdesyb8zY2FgZsXmE+M32k7Bb\nYZk3sEqVJ/8s37hxQ/Lnzy9//PFHknU3b94UX19fOXjwoFy6dEnuu+8++euvvxJsk1LyFhGZNGmS\ntG/fXl544YUk6w4cOCANGzaUwYMHS2hoqEOxdurUSTZs2CAiIjExMfLKK6/I888/79BnsyNnJ2+9\nYRnnu++gc2cYMMCamZmG56Kmm4gwcMNAjl85TnCPYHxyp/5wV+V6nnzDcu/evfj7++Pr62u/QRgb\nG8vNmzdp0qQJPj4+9iqR+fPns2TJEkJCQuyfT+6Gpc3ff/9NxYoV+e6776hTp06S9SLC7NmzOX/+\nPOPGjUs11nPnzvHCCy9w/PhxIiMjadGiBf/73/+SPK4sp3D2DUuHk7cxpi4QLCLJPnwsqydvuF0P\n3rUrDB2aOWNqL3DP48nJ2/Zw3FyJal5v3brl0JPd7+Ts2bM89thj9ocEK+dyS/I2xhQBVgAFRKRp\nCttk+eQNEBlpzcTM4N+DtI0ZE0nH5R0pnr84CzouwMvkiPJ7j+XJydtVRIRJkyZRuHBhXn75ZXeH\nky25a5LOJ8BYINs/ZcDbO3MTN4B3Lm9WdVvFqSunGPrl0ByXOJR7HTp0iAcffJADBw545DR7lbxU\nz7yNMaOASyIyxxizW0SSnd+dXc683enqrau0XNiSx+97nCktU5/FplwjJ555K9fL1Ek6xpjWQFUR\nmWpbdKftx7ZqRZ5GjSB3bvz9/fH3909rPB4pLAzGjYMpU1z7nExbL/Bm85tRxKcIIxqPcN1gSim3\n2LlzJzt37szwfu545m2MWQmUBSKwEvfDwG/AOBHZmmhbOd+gJoVPX0QmTsLnP/9xTyMRF4iMhOHD\nrf4oa9fCAw+4dry/w/6m6bymjHlkDM/Ufca1g6kk9MxbuYLbqk3iBvlKRB5JYZ1ER3/O2ZXfUGDK\nanwiINf//g/vzp2t5trZwNy58MorMHs2dOjg2rH+uPQH/gv8eefRd3iq5lOuHUwloMlbuYK7k/fv\nInJ/CutExHrqQVRkFOc+3UbRqWvJdVcpvN9+h9wtWqQ1No9kqwfv39/qEe7Kf5cO/3uY1otaM6/D\nPB6v+rjrBlIJaPJWruDWlrApJe7E8njnofzzgeT55V0uPulLVPeuhLdoQcyPP6Y1Po/j5wf790Ox\nYq7/haJW6Vqs676Ovmv7svvUbtcOppTKUlxaUOxTIB8Vxnch5tf3uVSnBNGtWnC905PE/vmnK4d1\nuTJlIFH7Y5eJ3wv8wJkDmTOoUmm0aNEiqlat6u4wcpRMmQ1SsHhhKrzdm5s/vcPlYpHEPFyb6wMG\nIP/+mxnDZ3mtKrdiVvtZtF3aliMXjrg7HKUSuHjxIq+99lqOnfbuLpk6la9o+ZKUnzOIq3v/x9Wr\nx4muVpXwESOsWrxs4OxZ1z0ns+MDHZkWMI02i9tw8spJ1wyiVDoMHjyYIUOG4O3t7e5QchS3zMMu\nUb0CZVcOJ3TLBMIO7SGy4r3cfPNNiIhwRzhOM3o0dOwIV6+6Zv9BDwUxqskoWi9szdlrZ10ziPJo\nO3bsoFatWsTGnSV8/PHHdO3a1b7+1KlTtEhUHDBjxgw6duyYYNnkyZN56aWXkux/xYoVaZpluWrV\nKm7cuMGTTz6JySZVZVmF25poGGMo3eB+Sm8aw8Vlw7m+bikRlSpya9YsiIlxV1gZ8umnUL48NGgA\nv/3mmjGG+A2hX+1+PLr4UUJvhrpmEOWxWrRoQatWrXj33Xe5cOECH374ITNnzkywTeIkGhQUxP79\n+7l586Z92ZdffkmfPn3s73/77TcCAgLYuHEjY8eOZfPmzTRq1IjGjRsn+e/hw4cBuHDhAmPHjmXW\nrFlaneMO6ekjm9zL2tX6dL+ioz+Xv5aNlNDqleXGfVUkYtUqkdjYZPvfero5c0RKlBBZu9Y1+7f1\nAm8wu4Fci7jmmkFyMBzp5w0Zf6XT9evXpWbNmtKmTRv57LPPEqw7ceKE+Pv7J/nMoEGDZOHChSIi\ncvz4cXn44YdFRCQ8PFxeffVV8fPzk6+//jpNcTz11FOyYMECEbGe1lO/fv30fJ0cI6WfKzyjn/f6\nDO8nKjKKc7O3UnTaOnIVL23ViGfBafb79kGXLvDFF1CrlvP3L9oL3GWyQp33xIkTmTlzJn///Tci\nQvPmzTHGcOvWLY4fP0716tYDrtu0acPEiRP58ccfGTZsGLt27eLNN98kX758DBs2jHbt2hEaGsqO\nHTuSfd5lSpYvX86iRYsIDg4G4OTJk3Tt2pX9+/e75PtmB26dpJNKYE5J3ja3wm9y/r1gSn4QTGwN\nX3zeeYdctWs7bf+Z4do1KFTIdfuPiY2h55qe3Iq+xapuq8jtpc+TdgZPT94nT56kbdu2tGvXDm9v\nb6ZMud3E7NSpU/Tr14/t27cn+VyjRo1YuHAhTz31FJs3b6ZkyZJERETwxhtv8MUXXzB16lRatmwJ\nWJdVJkyYkOASjMQ9/OHTTz9l0qRJHDt2jIIFCwJWP/GjR49StWpVgoKCkr2entPlmORtc/1SGKH/\nXUvpeSFEtWhJ/qnT8KpSxenjZFXaC9z5PD15BwQEMGzYMAICAqhfvz7z5s2zP/nm1KlT9O3blx07\ndiT53IIFC1i9ejUA69cn/Lv6559/MnjwYIoVK8Zbb73FPffck6aY9Mw7dW6dYekOCWrEi0ZojXgi\n8XuBv/jlix6ddFTGTZ8+nXz58tG2bVu8vb1577336N27NxHxKrVSqvro1q0bX331FX379k2yrkqV\nKmzatIkuXbowderUpB9OxbVr17h27VqaP6fSz+PPvBO7+OtpIiasolTIASKfGUiBCROgcGGXj+ss\nP/4IGzZYLWad+ZxM7QXuPJ5+5p1eERER1KxZk19++UVrst0gx515J5bVa8TvvhtCQpxfD27rBb7q\nyCre+uYt5+1YZRtLliwhMDBQE3c2keXOvOOLjY3l3JaD5Jm0yuojPmkyPgMGeHwfcVt/8K1brf7g\nDz7ovH1rL/CMy25n3qGhoTRq1IjSpUuzdu1a7rrrLneHlCPluBuWjoiJieHsZ99Q4PXV+EQaq494\np04e30fc1h98zhxo3955+9Ve4BmT3ZK38gyavO/A3kd82rrbfcQ9vEZ83z44d875D3c49O8hAhYF\naC/wdNDkrVxBk7cDskONuDPs/Xsv7Ze1Z3W31TS7N9nnRqtkaPJWrpDjb1g6wqdAPiqMtfqIh9Yp\nTnRL/2zRRzytGpZryPLOy7UXuFLZULZM3jYFixem/Fu2PuJZq0b877+ds5/4vcB/vfCrc3aqlHK7\nbJ28baw+4s9lmT7iERHQooX1jExn9AeP3wv8xOUTGd+hylaioqKIjY0lJiaGmCza0TMnyhHJ28Ze\nI755vFUjXqkiN/73P4+rEc+bF/bsgW3bnFcPHvRQEK80eYWARQHaC1wlsGLFCurWrUuxYsX49ttv\n3R2OclCOSt4Q10e84QNWH/ElLxK+dolH9hEvU8ZK3uXLWw89PuKEp58N9husvcBVEr169eLHH3+k\nXr16FChQwN3hKAfluORt4+XlRdnAOtz19RQuvB3Ezff/j5v3VyNy1Sqr27IH8PaG6dOtJ/S0b++c\nXxDGNB1DmypteGzJY1yL0F4U6jYR0dmXWUiOTd42uXLlotzTTSn44zQuvhRIxMhh3KhTh+hkurK5\nS//+cPCgdTklo4wxTAuYhm8pXzqu6Mit6FsZ36nKNN26deP9999PsMzf35+HH36YX375JcP710eZ\nZR05Pnnb5PHOQ/nnA8nz8ztc7OJLVPeuhLdoQczBg+4ODYC4tslOYYxhZruZlMhfgu6ruhMVE+W8\nnSuXev755xO0c7169Sp//vknH3zwAX379mX06NFcv37dvv7NN99M9nFmTZs2dUf4yomy5SQdZ7h+\nKYzLb66l1Nzs20c8MiaSDss7UCJ/Ce0FHo8jk3TMzp0ZHkfSOfu3Ro0afPHFF9x7773Mnz+fn376\niXfffZfY2FhmzpzJrFmzGDduHF26dEnTflu0aMH06dPtT+FRzqUzLDPZlb8ucG3Sasqs3E1Et6co\n8MYbmNKl3R0WALNnQ6lSGZtafyPqBoGLA/Et7cuHj32ovzbj+TMs33//fS5fvsykSZNo06YN//d/\n/0fteDOIf/nlFxo3bsw777zDgAEDHN6vJm/X0hmWmcyqER9k1Yhf+dOqER850iNqxB96CIYMgYkT\n018Pnj9PfjY8vYFv//6W8TvGOzdA5RK9e/dm+fLlnD17losXL9oTd3R0NO+99x49evRg4cKFDBgw\ngNdff5369evj5+dnf9WvX59GjRq5+VuojNIz7zQQEc7v+x0Zv5K7fjhG9MhR5B8+3Dl3EtPp3Dno\n2hWKFYNFi6BIkfTt50L4BZrNb0b/2v0Z2WSkc4PMYjz9zBugb9++hIaG4u/vz0svvcSOHTsYMWIE\nbdq0Yfz48eTLly/N+2zevDnTp0+nZs2aLohY6Zm3G9lqxEttHpuwRvyTT9xWI564Hvz339O3n5IF\nShISFMKM72cw+8Bs5wapnG7QoEF8+eWX9OrVC4A33niDZcuW8d///jddiRsgPDycMA/4jVI5Rs+8\nM+B2H/E1+ESI1Ue8c2e39RGfPx+aNoWM3Fc9FnqM5vOb5+he4FnhzHv//v1MnjyZjRs3ujsU5SC9\nYemBsmIf8Ts5/O9hWi9qnWN7gWeF5N2nTx+eeOIJOnfu7O5QlIP0sokHSlAj3qkWUU91IbxlS4+p\nEU+rWqVrsa77Ovqu7cvuU7vdHY6KZ+PGjdSoUQMvLy9N3Dmcnnm7gL1GfM4Wolq2cmuNuIjVXrZ8\n+bR/dtvxbTy9+mm+7PkldcvWdX5wHiornHmrrEfPvLOABH3Ei96y+oj/5z9u6SP+889Qty6sW5f2\nz2ovcKU8lyZvFypaodTtPuJuqhGvVQs2bkx/Pbj2AlfKMzmUvI0xvY0xe4wxe40xW4wxZV0dWHaS\noI/4T7szvY+4nx/s3w/bt6evP3jQQ0GMbjJae4Er5UFSveZtjCkAfAAMFJEYY0wg0EZEhifaTq95\nOyA2NpbnAyb1AAAcCklEQVRzIT+RZ+JKCp++iEychM9//gO5crl87MhIGD4cLl2C5cvT/vk3dr/B\n8l+Ws6vvLu7Kd5fzA/QQFStW5NSpU+4OQ2Uz9957LydPnkyyPFNKBY3V+GI0cF1EPkq0TpN3GsTE\nxHB25TcUmLIGn0isGvFOnTKlRjw8HNLTc19EGBUyit2nd7M1aCuF8hZyfnBK5TAuT97GmI+AQGA/\n0EdEIhOt1+SdDlmtRlxEGLhhIMevHCe4RzA+uX3cHZJSWVqmTdIxxnQBWorI84mWa/LOgFvhNzn/\nbjAlPwwmtoYvPu+8Q654neI8SUxsDD3W9OBW9C1WdV1Fnlx53B2SUlmWy5K3MSY3kF9EwuIt+0pE\nHkm0nUyc2N3+3t+/Fv7+tdIaT453/VIYof9dS+m5mVcjLgIvvQQDB8KDDzr2mciYSDou70jx/MW1\nF7hSabBz5052xusHP3nyZJcl78eBl4C2IhJhjGkCDBORrom20zNvJ7py+jzXJq/JtD7ic+daz8r8\n9FPH+4NrL3ClMs6ll02MMa8CPYELQBjwvIj8k2gbTd4ucPHX00RMWEWpkANEDnyWAuPHQ+HCLhnr\nu++gc2frmZkTJ4KXAyfTV29dpeXCljx232O83vJ1l8SlVHamjamyuX/3/U7s+JUU//4oMSNHke+l\nl1zSR9zWH7xoUVi6FAo5UFCivcCVSj+dHp/NlW5wP6U3jeHisuFcX7fU6iM+a5bT+4jb+oM3bw7e\n3o59RnuBK5X59Mw7C7pdI74an0iTqTXid6K9wJVKO71skgPZa8SnriVX8dIeUSNu6wU+94m5tK3W\n1q2xKJUV6GWTHMjeR/yXd7n4ZFwf8RYtiPnxR5eMd+VK6n1RbL3A+63rp73AlXIhTd7ZgE+BfFQY\n14WYX9/nUp0SRLf053qnJ4n980+njrNiBTRoAL/9duftGpZryLLOy+jyWRcOnDng1BiUUhZN3tlI\nweKFqfB2XB/xYhFWH/EBA5zWR/zZZ2HUKOs5man1B7f1Am+3rB1HLhxxyvhKqds0eWdDCfqIXz1u\n9REfMcIpfcT794fgYKs/+KRJd+4P3vGBjkxtPZU2i9tw8srJDI+tlLpNk3c2Zu8jvmUCYYf2EFnx\nXm6++WaG+4jb+oNv2wZz5tx526CHghjVZBStF7bWXuBKOZFWm+QQsbGxnNtykDyTVll9xCdNxmfA\ngAz1EY+MtKoT8zjQlyqn9AJXKq20VFA5JCYmhrOffUOB1zO3j7j2AlcqeZq8VZq4o0Zce4ErlZTW\neas0SVAj3imuRrxlywzXiF+6BIMGJV8PboxhZruZlMhfgu6ruhMVE5WhsZTKyTR553A+BfJRYaxV\nIx5apzjRrVpkqEa8UCHrMnpK9eC5vHKx6MlFRMZE0n99f2IljY+zV0oBmrxVnILFC1P+rYzXiHt7\nw/TpMHIkNGsG65O5kuady5tV3VZx6sophn45FGddulMqJ9HkrRIoWr7k7RrxK3+mu0Z8wADYuBEG\nD7bqwRPn5/x58rPh6Q18+/e3jN8x3nlfQKkcQpO3SlaJ6hUou+olQjeP59pPu9NVI26rB7/nnuSL\nWYr4FGFTz02sPrKaaV9Pc2L0SmV/Wm2iUuWKGvH4/g77m6bzmvLqI68ysO5Ap+xTqaxCSwWVy93u\nI+78GnHtBa5yKk3eKtPYa8SnrSPXXaXSVSN++jSUK5fwOZm2XuDzOszj8aqPOzdopTyU1nmrTGOv\nEf/5nds14mnsIz5qFHTsmLAe3NYLvO/avtoLXKlUaPJW6Ra/RjytfcQXLoTy5ZPWg2svcKUco8lb\nZVjCPuKRDtWI2+rBR42y6sHj9wdvVbkVs9vP1l7gSt2BJm/lNFYf8UEJ+4iPHAnXrqX4mf79rXrw\noUPh6NHbyzs80EF7gSt1B3rDUrmEiHD+u6PEjl9J8QN/ED1yFPmHD4e8eZPd/uZNyJcv6fKPvvuI\n9/a+x55+e7i70N0ujlqpzKfVJsojxcbGci7kJ/JMXJnuGnHtBa6yM03eyqMlqBGPEKtGvHNnh2rE\nRYSRISPZc3qP9gJX2Y4mb5UlJKgRL1bSqhFv0SLZbb/6ynrU2vjxYIzwzIZnOHHlhPYCV9mK1nmr\nLCFBjXhnX6K6d7VqxA8eTLLtfffB1q3w5JMQFmb4pN0nFM9XnKdWPaW9wFWOp8lbuUWSGvEWzbne\nuVOCGvEyZawz73LlrHrwP47mYnGnxdoLXCk0eSs3S1AjXjSuj/h//mOvEY9fD960KWz+wpvV3VZr\nL3CV42nyVh4hQY24rY/4yJH2PuL9+0NwsFVOqL3AldIblsoDiQjn9/1O7ITVFD9wNMUa8QvhF2g2\nvxn9a/dnZJORbopWqYzRG5Yq2zDGULrhA5Te9CoXlw4jfO0SIipV5NasWRATY9+uZIGShASFMOP7\nGcw+MNt9ASvlBnrmrTxeajXiO346Rq9t2gtcZU1a562yPXuN+NS15CpeGu+33yHSz58aNSCwz2HW\nFNRe4Crr0eStcoxb4Tc5/24wJT8MJraGL9fHvEOXybXhnn38Vqc9q59aRbN7m7k7TKUc4tJr3saY\n9saYb+JeIcaYe9MeolLO4VMgHxXGxdWI176Lot38CS7eCf/cJfHZuIwnl2kvcJX9pXrmbYzxAfYA\nrUQkzBhTDxguIj0Tbadn3sotrvx1gWuTVlNm5W4O1+5Bf9OQc+1fZXuf7VQvWd3d4Sl1Ry67bGKM\nKQHUFJGdce+rA+NF5OlE22nyVm518dfTRExYRcktBzjY5RF63vUDbU51IdfFXPjc48NzU56jYqWK\n7g5TqQQy5Zq3MaYqMBt4XkR+TbROk7dyO1uN+KkRS/nom4I8LWPIRz5ucpPlVZYzOWSyJnDlUVxe\n522M6Qa8D3RPnLiV8hS2GvHV91a3J26AfOSj+5/d+Xj8x26OUCnnyO3IRsaY/sCjQAcRSbGd26RJ\nS+1/9vevhb9/rQwHqFR63Dpj7InbJh/5uHXoHxBxqI+4Uq6wc+dOdu7cmeH9OHLNuzgQAtSVO2ys\nl02UJxndawX+S3omSOA3ucmmvCN4r14xzP/+B4884sYIlbK48rJJc6AksM0Ys90Ys9MY83maI1Qq\nEz03pQXLq8zlJjcBK3G/U+It5ncsyTvh/Ynt0Qvat4fDh90cqVLp48iZt0l8xm2MyZP48omeeStP\nc/LEeT4ev4NbZ8CnLPSf+Ah9N3/EmT8foFDwQnZ0m03JT9+ERx+FyZOhUiV3h6xyIJ1hqZQDrkfc\n5JGZY/C52pjVAxZwT+Hr8Pbb8OGH0LMnjBsHpUq5O0yVg2hXQaUcUDBvPrYNmMzlgjuZ/NUzUKgQ\nTJoER45YNzEffBAmTrT3EVfKU2nyVjlO8YKF2dH/NTaf3siodSOshaVKwfvvw/ffw4kTULUqvPce\nRES4N1ilUqDJW+VIZYsWZ3vfKSw6Mp83N78BWK3CT1AJFi6EkBDr6cf33w8LFiToI66UJ9DkrXKs\nKiXLsCVoIm//MI0Zu6fz44/g5wfr1gG+vrBxIyxaBLNmQe3asGGDVSOulAfQG5Yqx9t34ncCl73G\nR60/pio96dzZembmxIng5YWVsDduhFdfhaJFQWvElRPpDUul0qlBpftZ2Wkkg0Oe41ze9ezfD9u3\nQ8eOcPUq1o3M9u3hp59g4EDo1QvatYNDh9wdusrBNHkrBbR+oDbz2g2h98ZeHLm6g23boHx5GD48\n3ka5ckHv3vD77xAQYNWHBwVZNziVymR62USpeObv28bL2+az6ekQ6leqz61b4OOTwsbXrmmNuMow\nvWyilBP0bdCKCY278vjyQH7555eUEzckrBEHrRFXmUqTt1KJvOj/BC88HMiji1tx/MLx1D9QqhR8\n8IHWiKtMpclbqWRMCOxO9/sb0Xp+C85eOWtfHhsLzz57+2Q7gUrxasS3bbNqxBcu1Bpx5RKavJVK\nwVsd+tGi/P20mutP6PVQwCodbNAAmjWLqwdPjq+vVRO+aBF88onWiCuX0BuWSt2BiNB14VROXg1n\n17NfUSBvAQC++46k9eDJ78CqER8zBooU0RpxlYR2FVTKRaJjonl87utEROdjy3+2kTdPXgDOnYOu\nXa15O8uWQcGCd9hJTAwsWQITJkCtWvDGG9YZusrxtNpEKRfJnSs3G/qNIVou8+SCDsTEWtewy5Sx\nLm23aQP58qWyk/g14q1bW3XiWiOuMkCTt1IOyJvbm80DJnIm/Bg9FnfH9hurtzcMGWLlZsd2lBde\nfBH++AOqVIF69WDoUDh/3nXBq2xJk7dSDrL1Aj94cR+DVj6bsZ0VLqw14ipDNHkrlQa2XuCbTq3n\nlfWjUtzu/Pm4viipSVwjXq2a1Vdca8RVKjR5K5VGZYsWZ1uf11jw61x7L/DEVqywSgp/+83Bndpq\nxLdsud1HXGvE1R1o8lYqHe4rdXeCXuCJvfACjBoFTZveoR48OfFrxLWPuLoDLRVUKgP2nviNx5ZN\n4aPWH9PTr2eS9fv2QZcuMGCAVSWYYj14crRGPEfQUkGl3KBhpQfsvcDX/5T05KVBA9i/3yopXLYs\njTu39RE/eFD7iKsk9MxbKSdYc/Ab+gdP5/Mu62hxf4sk6yMjIXfuNJ55JxYRATNnwn//a/USf+01\n61q5ytL0zFspN+pUuzHvte5Pl9VPsv/E/iTrvb0zmLghYY145cpaI57DafJWykn6NmjFxCbd7L3A\nXaZwYZg82aoRN8aqEZ80yXo4hMoxNHkr5URDm7d3uBf4mTPQr5+D9eDJKVXKqgk/cACOH7f6iGuN\neI6hyVspJ4vfC/zc1XMpbleiBOTPn8Z68ORUrHi7j7itRnzBAq0Rz+Y0eSvlArZe4C3nNLf3Ak/M\n2xumT7fqwe/YH9xRtWpZNeGLF8Ps2fDQQ7B+vdaIZ1NabaKUi6TUCzw5tnrwZ56x6sGdMLjWiGcR\nWm2ilIcxxrC818vclc+Lx+e2ISIq5WvRtnrwKlWcNnjSGvH27eHwYScNoNxNk7dSLpRSL/DklCkD\nPZNO0syYxH3EW7fWPuLZhCZvpVwsfi/wnoufxlmXKtMWhPYRz240eSuVCWy9wH+8uJdBKwem6bPH\njllPrXeK+H3EbTXi2kc8S9LkrVQmud0LfMMde4EnNnIkdOyYgXrw5NhqxG19xKtWhffe0xrxLMSh\n5G2MyW0s3q4OSKnsrGzR4mzvO+WOvcATW7ECypd3Qj14cmx9xLdutbpnaY14luHomXcHYBewz4Wx\nKJUjVClZ5o69wBOLXw+e5v7gjkpcI659xD1emuq8jTFfi0iTFNZpnbdSabDvxO8ELn2N+U8uokP1\nDg595rvv4KmnYMcOa2KlS9hqxF99FYoW1RpxF8usOm/9XUopJ2lQ6X4mVx/EM8HPsPvUboc+4+dn\nXTpxWeKG2zXiP/1kzRrSGnGPpDcslXKj6oUqM7/9fLp81oUDZw449Jm8eV0clE2uXNCnz+0a8YAA\nrRH3IGm9bLJbRJqlsE7o0+f2gtq1rZdSSqnbDh60XjYLFqTrsolTk/c17SesVJp4eXmRP39+ABb9\ntIgx28ewu+9uKhVL2xNyNm+Gb79Nx3My0+P8eXjjDevm5pAh8PLLVv24Spf0XvNOa/L+SkSSvXNh\n3bDUO9NKZcRH333Ee3vfY0+/Pdxd6G6HP3fuHHTtCsWKWQ+eL1LEhUHanDxp/WuxebN1c/O55zLx\nmk724fIblsaYfED+tA6glHLcEL8h9Kvdj0cXP8qlG5cc/lyZMlaZtsvqwZNj6yMev0Z84UKtEc8k\nTm4Jq2feSmWUiDAqZBS7T+9ma9BWCuUtlKbPz5ljnQjPnw+PP+6aGJP11VcwerQ11f6//7WedG/S\nfEKZ42TKZZNUAtDkrZSTiAgDNwzk+JXjBPcIxie3T5o+/913EB0NjRu7KMCUiFiTe8aM0RpxB2ny\nViqbiYmNoceaHkRER7Cq2ypye+V2d0iOi4mBJUusa+K1alln4rVquTsqj6QPY1Aqm8nllYtFTy4i\nMiaS/uv6EyvOai2YCRL3EQ8IsN5rjbjTaPJWyoN55/JmVbdVnLxykqFfDs1wL/Dff3dSYI6K30e8\ncmXtI+5EmryV8nD58+Rnw9Mb+Pbvbxm/Y3y69xMWBo8+arXvdlp/cEcVKqR9xJ1Mk7dSWUARnyJs\n6rmJ1UdWM+3raenaR+HC1oOOt293QX9wR9n6iB84oH3EM0iTt1JZRMkCJQkJCmHG9zOYfWB2uvbh\nlnrw5NhqxENCtEY8nbTaRKks5o9Lf+C/wJ+3H32b7jW7p3s/c+bAjBnWU+tdPqU+NXv2wCuv5Mga\ncS0VVCoHOfTvIQIWBTCvwzwer5r+mTiRkdbDHjxC/D7iRYpYNeJNm7o7KpfT5K1UDrP37708sewJ\nVnVbRbN7k+0XlzXFrxGvWdM6E/f1dXdULqN13krlMA3LNWRZ52Vp6gWeJcSvEQ8IsEpktI94Epq8\nlcrCWlVuxaz2s2i7tC1HLhxxyj6nTHHRczLTKn6NeJUqVo34iy9qjXgcTd5KZXEdH+jItIBpPLr4\nUU5czvjZaZs2Vptut9SDJyd+jThojXgcTd5KZQNBDwXxSpNXCFgUwNlrZzO0Lz8/qwJl+3bo0MFN\n9eDJsdWIf/+9dQmlWjXrfQ6tEdfkrVQ2MdhvsL0XeOjN0Azty1YPXqGCVQ9+7JiTgnSGSpWsmvAt\nW6w68RxaI67VJkplI7Ze4HtO72Fr760U9C6Y4X0uWQKBgVC8uBMCdAVbjfjVq1ZlSvv2WapGXEsF\nlVLA7V7gJ66cYGOPjWnuBZ4l2WrEx4yx+gBkoRpxTd5KKbss3Qs8IxLXiL/5psf3Edc6b6WUnat7\ngUdFedh1cJvENeLZuI+4Jm+lsiln9wKPb/9+aNTIQ+rBk5NcH/FsViOuyVupbCx+L/Bx28c5bb+N\nG0NwsIfVgycnuRrxSZOyRY24Jm+lsjlbL/A1v61Jdy/w5MSvB3dbf3BHxe8jfvx4tqgR1+StVA7g\njF7gyYnfH3zsWKft1nXi9xHfujVL14hrtYlSOcix0GM0n9+cdx59h6dqPuXUfXtUe1lHffWVVSN+\n5YrbasS1VFAp5ZDD/x6m9aLWGe4Fnm3ErxG39RF/5JFMG15LBZVSDqlVuhbruq+j79q+7D61293h\nuJ8x1hn3wYMwcCD06mW9P3zY3ZHdkSZvpXKgzOgFHh1t5cEjzulU63rxa8Rbt7ZeHlwjrslbqRwq\nfi/wXy/86vT9584NLVtC8+YeXA+enPg14pUqeWyNuCZvpXIwWy/wNovbOKUXeGL9+1uXkz2+Hjw5\nhQvD5Mke20dck7dSOVzQQ0GMbjLaKb3Ak5O4HvzmTacP4Vrxa8RPnICqVeG999xeI67JWynFEL8h\n9K3d1ym9wJNjqwd/4gnwyapNDm014lu3Wl/m/vthwQK31YhrqaBSCrBayY4MGWn1Ag/aSqG8hdwd\nkmez1Yjb+oi3a5euGnGt81ZKZZitF/jxK8cJ7hGcM3qBZ4StRvzVV6Fo0XTViGvyVko5RUxsDE+v\nfprImEhWdl1Jnlx5XDre339b/aOKFHHpMK4Vv494rVrWmbiDfcR1ko5SyilyeeVicafFVi/w9c7v\nBZ7Y8uXWTc0sUw+eHDfUiDuUvI0xjY0xm40xu40xnxljirksIqWU29l6gZ+6csrpvcATGzECRo+G\nZs2yWD14cpLrIz50qEtqxFNN3saYgsD7QBcRaQbMA952eiSZZOfOne4OwSFZIc6sECNonOkVvxf4\n+B3j7ctdEacr6sHdejwLF77dR9wYl9SIO3Lm/SiwXkSuAYjIl0BZp0WQyTztL0hKskKcWSFG0Dgz\nwtYLfPWR1fZe4K6Ks0GD2/Xgn3+e8f15xPF0YY24I8m7IvB7omU3jDEFMjy6Usrjxe8FPuvALJeO\nVaaMlbw7dXLpMJkvuRrxDPYRdyR5R8S94ssPuPYWtFLKY5QrXI6QoBAm75rMz+d/dulYefJkekvt\nzFOrFmzYAIsXw6xZULt2uneVaqmgMaYVUE9E/i/esp9E5KFE22mdoFJKpUN6SgVzO7DNAeB1Y8ws\nEblsjGkHfOuMwZVSSqWPQ5N0jDEtgP8C4cC/wLMict3FsSmllEqB02ZYZiXGGG8gSkTEGOMF+IjI\nDTfEkRuIAfKISGRmj++o1OI0xuQVkYh42+YRkazWO85j6PF0nux8LNM0wzK1yTrGmGrGmGBjzE5j\nzBfGmHudG65z4gRmAl8bY7YD24GXMj9KADoAu4B9ya30oMlRd4wT+MIYsyPueG4DemRaZPEYY9ob\nY76Je4Uk/vnzlOOZWpx4zvHsbYzZY4zZa4zZYowpm2i9249najHiIcfSxhhT1xhzLpnlac+dIuLQ\nCygI7AcKxb1/DJgbb73B+stdNu79w8AXju7fWa/U4oxbtgYomNmx3SHmr9PzPTwhzrjlIR5wDH3i\njlfhuPf1gCWedjxTi9ODjmcBYA6QK+59IPCuJx3P1GL0lGMZL5YiwCZgT6Ll6cqdaTnzTm2yzsPA\nzyJyJm79j0C0MSazSwodmVRUFBgX9y/yQmNMyUyOMbHkij09cXJUSkWpBYwxH8edNXzopjkABYGR\nImKbwnaDhL9ZesrxTC1O8IDjKSLhIjJARGKMMQaoDfwZbxO3H08HYgQPOJbxfAKMBRLPH01X7kxL\n8q7InSfrJLf+OFA1DWM4Q3JxJJ5UdB+wRkRaACuwLqN4mopkgclRce0TKgNviYg/cAR4LbPjEJGL\nIrIzLqaqwAxgSrxNKuIBxzO1OD3leMaL5yPgD+AhIP4MnYp4wPGElGP0pGNpjBmF9VvAAawz7fgq\nko7cmZbkndpkneTW+wDeaRjDGRyZVFRFRL4DEJFg4O5Mii0tssTkKLGqjsqLiO2M52OgkbviMcZ0\nw+rF011E4j9V16OOZ0pxetrxFJEhInIfsBp4L94qjzmeKcXoKcfSGNMaqCoic2yLEm2SrtyZluT9\nG/BAomV3i8iVO6y/HzidhjGcIbU4IfmD52kc+R6ewn48xbpod8stQRjTH+gEdBCRxDeFPOZ4phIn\neMDxNMbkNsYUjhfHKsA33iZuP54OxAgecCyBZ4HqxpjtxpgdQK24G6yt49anK3emJXkfADra7ign\nM1nnBFDbdpfUGPMwcF1EnP9AvAzEGVcmuNd2V9oY4w/8k8kxJpbc/4fUjrc7JIkz7n7B18aYInHv\newN7MzswY0xxYAjwtIhEJbOJRxzP1OL0lOOJdU17jTEmb1wcTYD4Tyf2hON5xxg95ViKSFcRaSIi\nLeMu1f4sIg1FZGvcJunKnY7MsLQFcMUYMwar9MY+WccYMxr4UUS2GGOeAxYaY2KAMOCZNH/TDHIw\nzqHAZ8aYKOAiMDiz47QxxuTD+nXT9v5dYJ6IHErue7gpzNTifAPYbIy5AZwEXnBDiM2BksC2uJtX\nAFew/mLM96Dj6Uicbj+eIvKFMeYh4IAx5gLW3+fnPenn08EY3X4sk1ESIKO5M0dO0lHZjzFJn8Nn\njMmTwlm422SVOJXn0+StlFJZkD7DUimlsiBN3koplQVp8lZKqSxIk7dSSmVBmryVUioL0uStlFJZ\nkCZvpZTKgv4fu5z8fntIV+sAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x114b5bad0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(0, 4, 100)\n",
    "\n",
    "y0 = (-3*x+10)/2\n",
    "y1 = -2*x+6\n",
    "y2 = -x+4\n",
    "y3 = np.ones_like(x)\n",
    "\n",
    "y4 = np.minimum(y1, y2)\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(x, y0, \"--\", label=\"objective function\")\n",
    "plt.plot(x, y1, label=\"2x+y<=6\")\n",
    "plt.plot(x, y2, label=\"x+y<=4\")\n",
    "plt.plot(x, y3, label=\"y>=1\")\n",
    "plt.plot(res.x[0], res.x[1], \"o\")\n",
    "plt.fill_between(x, y3, y4, where=y4>y3, facecolor='yellow', alpha=0.3)\n",
    "plt.ylim(0, 5)\n",
    "plt.legend(loc=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PuLP で解くと以下のようになります。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[2.0, 2.0]"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pulp\n",
    "\n",
    "x = pulp.LpVariable(\"x\")\n",
    "y = pulp.LpVariable(\"y\")\n",
    "problem = pulp.LpProblem(\"A simple maximization objective\", pulp.LpMaximize)\n",
    "problem += 3*x + 2*y, \"The objective function\"\n",
    "problem += 2*x + y <= 6, \"1st constraint\"\n",
    "problem += x + y <= 4, \"2nd constraint\"\n",
    "problem += y >= 1, \"3rd constraint\"\n",
    "problem.solve()\n",
    "\n",
    "[x.value(), y.value()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}