{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 602 Satisfiable Modulo Theories(SMT)\n",
    "\n",
    "## PREPARATION\n",
    "pip install z3-solver\n",
    "\n",
    "## Example: Tonkoke\n",
    "http://mechlog.jp/materials/601_Linear_programming.html\n",
    "\n",
    "\n",
    "$$ \\exists x_1 \\exists x_2 \\exists x_3 \\varphi(y,x_1,x_2, x_3) $$\n",
    "where\n",
    "$$ \\varphi(y,x_1,x_2, x_3) = (y=15x_1+18x_2+30x_3 \\wedge \\\\\n",
    "                                2x_1+x_2+x_3 \\le 60 \\wedge \\\\\n",
    "                                x_1+2x_2+x_3 \\le 60 \\wedge \\\\\n",
    "                                x_3 \\le 30 \\wedge \\\\\n",
    "                                x_1, x_2, x_3 \\ge 0) $$\n",
    "\n",
    "### To parametric Optimization problem\n",
    "\n",
    "$$ \\varphi(y,x_1,x_2, x_3) = (y=15x_1+18x_2+30x_3 \\wedge \\\\\n",
    "                                2x_1+x_2+x_3 \\le 60 \\wedge \\\\\n",
    "                                x_1+2x_2+x_3 \\le 60 + {\\color{red}{\\theta_1}} \\wedge \\\\\n",
    "                                x_3 \\le 30 \\wedge \\\\\n",
    "                                x_1, x_2, x_3 \\ge 0) $$\n",
    "\n",
    "\n",
    "## Developler of SAT/SMT solver named Z3:\n",
    "https://stackoverflow.com/users/1096362/nikolaj-bjorner\n",
    "\n",
    "## Description example of Z3:\n",
    "\n",
    "e.g. ForAll([x], Implies(x >= 25, Exists([u,v], And(u >= 0, v >= 0, x == 3*u + 5*v))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import z3\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_vars():\n",
    "    x1, x2, x3, y, theta1 = z3.Reals('x1 x2 x3 y theta1')\n",
    "    return x1, x2, x3, y, theta1\n",
    "\n",
    "def check_feasibility(funcs):\n",
    "    solver = z3.Solver()\n",
    "    solver.add(funcs)\n",
    "    return solver.check(), solver.model()\n",
    "\n",
    "def maximize(problem, obj_func):\n",
    "    model = z3.Optimize()\n",
    "    model.add(problem)\n",
    "    model.maximize(obj_func)\n",
    "    return model.check(), model.model()\n",
    "\n",
    "def eliminate(problem):\n",
    "    return z3.simplify(z3.Tactic('qe')(problem).as_expr())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem definition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_problem_with_sensitivity_parameters(x1, x2, x3, y, theta1):\n",
    "    return z3.Exists([x1, x2, x3], z3.And(15*x1 + 18*x2 + 30*x3 == y,\n",
    "                                           2*x1 + x2 + x3 <= 60,\n",
    "                                           x1 + 2*x2 + x3 <= 60 + theta1,\n",
    "                                                       x3 <= 30,\n",
    "                                           x1 >=0, x2 >=0, x3 >=0\n",
    "                                          )\n",
    "                     )\n",
    "\n",
    "def make_problem(x1, x2, x3, y):\n",
    "    return z3.Exists([x1, x2, x3], z3.And(15*x1 + 18*x2 + 30*x3 == y,\n",
    "                                           2*x1 + x2 + x3 <= 60,\n",
    "                                           x1 + 2*x2 + x3 <= 60,\n",
    "                                                       x3 <= 30,\n",
    "                                           x1 >=0, x2 >=0, x3 >=0\n",
    "                                          )\n",
    "                     )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Execution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=== Evaluate whether feasible problem ===\n",
      "(sat, [y = 0])\n",
      "(sat, [y = 0, theta1 = -60])\n",
      "=== Optimize the problem ===\n",
      "(sat, [y = 1230])\n",
      "=== Do QE to the problem ===\n",
      "Or(And(Not(y >= 1440),\n",
      "       Not(650/7 <= 5/63*y + -5/7*theta1),\n",
      "       y >= 900),\n",
      "   And(y <= 1800,\n",
      "       1/42*y + -5/7*theta1 <= 300/7,\n",
      "       Not(y >= 900),\n",
      "       y >= 0),\n",
      "   And(Not(y >= 1800),\n",
      "       Not(300/7 <= 1/42*y + -5/7*theta1),\n",
      "       y <= 900,\n",
      "       Not(y <= 0)),\n",
      "   And(Not(y >= 1440),\n",
      "       5/36*y + -5/4*theta1 <= 325/2,\n",
      "       Not(y <= 900)),\n",
      "   And(Not(3075/14 <= 5/28*y + -5/4*theta1),\n",
      "       Not(-325/2 <= -5/36*y + 5/4*theta1),\n",
      "       Not(225/2 <= 1/12*y + -5/4*theta1)),\n",
      "   And(5/28*y + -5/4*theta1 <= 3075/14,\n",
      "       Not(-325/2 <= -5/36*y + 5/4*theta1),\n",
      "       1/12*y + -5/4*theta1 <= 225/2),\n",
      "   And(y <= 1440, 5/36*y + -5/4*theta1 <= 325/2, y >= 900),\n",
      "   And(Not(-750/7 <= -5/84*y + -5/7*theta1),\n",
      "       Not(y >= 1440),\n",
      "       Not(y >= 1800),\n",
      "       y >= 1080),\n",
      "   And(Not(y >= 1080),\n",
      "       Not(300/7 <= 5/63*y + -5/7*theta1),\n",
      "       Not(y <= 0)),\n",
      "   And(5/84*y + 5/7*theta1 <= 750/7,\n",
      "       Not(650/7 <= 5/63*y + -5/7*theta1),\n",
      "       Not(300/7 <= 1/42*y + -5/7*theta1),\n",
      "       -5/63*y + 5/7*theta1 <= -300/7),\n",
      "   And(Not(3075/14 <= 5/28*y + -5/4*theta1),\n",
      "       -5/36*y + 5/4*theta1 <= -325/2,\n",
      "       Not(225/2 <= 1/12*y + -5/4*theta1)),\n",
      "   And(Not(y >= 1440),\n",
      "       Not(325/2 <= 5/36*y + -5/4*theta1),\n",
      "       y >= 900))\n"
     ]
    }
   ],
   "source": [
    "x1, x2, x3, y, theta1 = make_vars()\n",
    "original_problem = make_problem(x1, x2, x3, y)\n",
    "sensitivity_analysis = make_problem_with_sensitivity_parameters(x1, x2, x3, y, theta1)\n",
    "\n",
    "print('=== Evaluate whether feasible problem ===')\n",
    "print(check_feasibility(original_problem))\n",
    "print(check_feasibility(sensitivity_analysis))\n",
    "\n",
    "print('=== Optimize the problem ===')\n",
    "print(maximize(original_problem, y))\n",
    "\n",
    "print('=== Do QE to the problem ===')\n",
    "print(eliminate(sensitivity_analysis))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbIAAAD5CAYAAABCm7AHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xt8TVf6+PHPk4SEuhOiiGsiV4qUDnUp1aItLa2qJIhbiDvt1GWGkjHVX2Vavv2ipYqvSxnTTtG7ThUzoxWlQlCXqrtQlwhpIjnr98fZ6WQiISE5Jyd53q/XeWWftdfe+zm5nCdr7XXWEmMMSimllKtyc3YASiml1L3QRKaUUsqlaSJTSinl0jSRKaWUcmmayJRSSrk0TWRKKaVcmsMSmYgsFZEkEdmXrewVETktInusR49s+6aIyBEROSQij2cr72aVHRGRyY6KXymlVPEkjvocmYh0AFKAFcaYEKvsFSDFGDM3R90gYA3QGrgf2Az4W7t/BLoCp4CdwAvGmERHvAallFLFj4ejLmSM2SoiDfJZvRfwvjEmDfhJRI5gT2oAR4wxxwBE5H2rriYypZQqpYrDPbLRIrLX6nqsapXVAU5mq3PKKsurXCmlVCnlsBZZHhYCsYCxvsYBgwHJpa4h98Sba9+oiAwHhgPcd999rQICAgojXqWUg+zateuiMca7APVrenh4LAFCKB7/pKvCYQP2ZWRkDG3VqlVSbhWcmsiMMeeztkVkMbDJenoKqJetal3gjLWdV3nOc78DvAMQFhZm4uPjCylqpZQjiMjPBanv4eGxxMfHJ9Db2/uym5ubTiJbQthsNrlw4ULQuXPnlgA9c6vj1P9aRKR2tqfPAFkjGjcA/UTEU0QaAn7Ad9gHd/iJSEMRKQv0s+oqpVSIt7d3siaxksXNzc14e3tfxd7SzpXDWmQisgboBNQQkVPADKCTiDyAvXvwOBANYIzZLyLrsA/iyABGGWMyrfOMBj4H3IGlxpj9jnoNSqlizU2TWMlk/VzzbHg5rEVmjHnBGFPbGFPGGFPXGPOuMSbSGBNqjGlmjOlpjDmbrf5sY0xjY0xTY8yn2co/Mcb4W/tmOyp+pZS6E3d391YBAQFBTZo0CW7atGnQK6+8UiszM/OOx0VHR9dt0qRJcHR0dN27uW758uVbABw6dKjsokWLqt3NOXIzefJkn+zPW7RoUSwHG+gNUaWUKiSenp62gwcPJh45cmT/P/7xjx+/+OKLyi+++OL9dzpu1apV3gkJCYlvv/32qXu5/uHDhz3Xrl1baIls/vz52W//sHv37oOFde7CpIlMKVUqLdq5qNr9cfeHus10a3V/3P2hi3YWXksGoE6dOhlLliw5/t5779W02WxkZGQQHR1dNyQkJNDf3z/o9ddfrwHQuXPnJqmpqW4tWrQIXLx4cdXVq1dXbtasWUBgYGBQ27Zt/U+ePOkBMHHixPunT59eK+v8fn5+wYcOHSqb/ZrTpk2rEx8fXyEgICBo5syZNbPvs9lsREdH1/Xz8wv29/cPWrx4cVWATZs2VQwLC2vatWvXxo0bNw7u37+/b2ZmJjExMXXS0tLcAgICgnr27NkQ/tPy27RpU8UHH3ywaY8ePRo1aNAgJCYmps7ChQurhYaGBvr7+wft37/fE6BPnz4N3nvvvayPVRX4+Pxy9vB7pZRyuEU7F1Wb8MWE+r9m/OoGcDblbNkJX0yoDzDiwRGXCus6QUFB6TabjdOnT3usXbu2SuXKlTP37dt3IDU1VR588MGAp556Kvkf//jHkfLly7c4ePBgIsCFCxfc+/Xrd9DNzY2//OUvNWbNmuWzePHifLXUZs+efTouLq7W119/fSTnvhUrVlRJSEgod+DAgf1nz571aN26deBjjz2WApCQkHDf7t279/n7+6d36NDBb8WKFVUXLFhwetmyZTWz4srp4MGD5davX3+sZs2aGfXr1w/19PS8mJCQcCA2NrZmXFxczaVLl57M7bjCOj47bZEppUqdWVtn1clKYll+zfjVbdbWWYU+wULWNICbN2+utG7duuoBAQFBLVq0CLx8+bJHYmKiV876P/30U9n27dv7+fv7B82fP9/n4MGD5Qojjm3btlXs27fvJQ8PD+rVq5fRpk2blO3bt5cHCA0NvR4UFJTu4eFB3759L23btq3Cnc4XGhp6vX79+jfLlStnfH1907p3734VoHnz5qknTpwoW9THZ6ctMqVUqXMu5Vyub5R5ld+txMTEsu7u7tSpUyfDGCNxcXEn+vTpk3y7Y0aPHu07bty4c+Hh4Vc3bdpUcdasWfcDeHh4GJvN9lu9tLS03CaOyNPt5tUVkds+z42np+dvJ3Rzc8PLy8tkbWdmZkpWzFmDXWw2Gzdv3pSCHJ9f2iJTSpU6PhV80gtSfjfOnDnjMWzYsPpRUVFJbm5udO3a9erChQu9sxLQ3r17PZOTk295D7527Zq7r6/vTYBly5ZVzypv0KBB2p49e+4D2L59e/nTp0/fch+pcuXKmSkpKe65xdOxY8dr69evr5aRkcGZM2c8vvvuuwrt27e/DvauxYMHD5bNzMxk/fr11dq3b38N7ImooAkzu/r166fv2rWrPMCqVauqZGRk3PW5bkcTmVKq1JneYfppLw8vW/YyLw8v2/QO00/fy3mzBkc0adIk+JFHHvHv0qVL8ty5c88ATJgw4WJAQMCvoaGhgX5+fsHDhg2rn72FkmXatGlnXnjhhcatWrVqWr169Yys8gEDBly+fPmye0BAQNBbb73lXb9+/V9zHtu6detUDw8P07Rp01sGe0RGRl4JDg5ODQwMDO7UqZP/zJkzT/n6+mYAPPDAAymTJk2q6+/vH+zr65sWGRl5BSA8PPxCYGDgb4M9CmrMmDEX/vWvf1UMDQ0N3LFjx33lypWz3fmognPYMi7OpFNUKeV6RGSXMSYsv/V/+OGH482bN7+Y3/qLdi6qNmvrrDrnUs6V9angkz69w/TThTnQw1Vs2rSpYl4DRIqTH374oUbz5s0b5LZP75EppUqlEQ+OuFQaE1dJpF2LSilVij355JPXintr7E40kSmllHJpmsiUUg6xe/dupk6detth4ErdDU1kSqkilZmZyZw5c2jTpg3Lly3j/Pnzdz5IqQLQRKaUKjI//fQTnTp1YsqUKTzTpg0JW7fi4+Nz5wOVKgAdtaiUKnTGGJYvX87YsWMRm43/mzCB8E6dkKpV73ywUgWkiUwpVaguXrxIdHQ0H3zwAZ1CQ1k+fjy+3t7ODkuVYNq1qJQqNJ9++imhISFs2riR16Oi+Co2VpNYLjIyMoiKiqrXpEmTYH9//6DExMRCneOxoC5evOg+Z86cu/pB5WexzXs5f35oIlNK3bMbN24watQoevToQQ0vL3bOncuLzzyDm5u+xeRm6tSptRs1apR25MiR/dHR0UlvvvlmzTsfVXR++eUX93fffTfXGGw2G7db5To/i23e7vyFQX/LlFL3ZOfOnbR44AEWLlzIpKefZmdcHM0a3tXUfC5v3Lhx98fGxv72hj1mzJg6f/rTn/7rDTw5Odnt448/rvLHP/4xCaBx48Zpx44dy9dCkm+99VZ1f3//oKZNmwY9/fTTv32TX3nllVp+fn7Bfn5+wbNmzaoJcOjQobKNGjUK7tevX/0mTZoEt2vXzi8lJUWSk5PdOnXq1KRp06ZBfn5+wYsXL646adKkuidPnvQMCAgIio6Orpt1bEREhG9wcHDQ0aNHyz766KONg4ODA5s0aRI8d+7cGlnXzlosM6/rAeQ8/718j3PjsHtkIrIUeBJIMsaE5Nj3IvA64G2MuSj2NQTmAT2AG8AgY8z3Vt2BwB+sQ/9kjFnuqNeglPqPjIwMXn31VWbNmkXtqlX5KjaWR5o1c3ZYThUTE3PxmWeeafzHP/4xKTMzk7///e9Vd+7ceSB7nQ0bNlQ6e/Zs2YCAgCCAq1evuj/88MPX7nTu+Ph4r7lz59b+97//fbB27doZ58+fdwfYtm1b+dWrV1fftWvXAWMMrVq1CuzSpcu1GjVqZJ44ccJr5cqVx9q2bftzjx49Gq1YsaJq+fLlbT4+Pje3bNlyBOytpQ4dOlx/8skny2Utonno0KGyx48f91q8ePHxlStXngBYtWrV8Vq1amWmpKRIixYtgiIiIi77+Pj8V1Mtt+vFxMRciouLO5X9/IXNkYM9lgFvASuyF4pIPaArcCJbcXfAz3q0ARYCbUSkGjADCAMMsEtENhhjLhd59Eqp3xw9epTIiAj+vWMH4R078lZ0NFUq3HEtRocZPHhwvX379pUvzHOGhITcuNOqxU2bNk2vUqVKxj//+c9yZ8+eLRMcHHwj55v97t27y02ePPnM73//+wsAzz//fP3Q0NDUxMTEsq+88krt5ORk988+++xYznN//vnnlZ566qnLtWvXzgCoVatWJsCWLVsq9OjR40qlSpVsAE888cTlr7/+uuJzzz13pU6dOmlt27ZNBWjRosWN48ePe0ZERFyaNm1avZEjR9bp1avX1W7duqVcvHjxlqVfateund6lS5frWc9fe+21Wh9//HEVgHPnzpXZv3+/l4+Pz/Xsx+R2vfx8b++Vw7oWjTFbgdwm6HwD+D32xJSlF7DC2O0AqohIbeBx4EtjzCUreX0JdCvi0JVSFmMMixcvpnnz5hzYt4/3X3qJlZMmFask5mxRUVEXlyxZUuO9996rERUV9curr77qHRAQEBQQEBB0/PjxMpcvX/YoX768DeDmzZts3bq10rPPPnslKCgofd26dT/ndV5jDCJyy7Qot5sppWzZsr/tdHd3NxkZGdKsWbO077//PjE0NDR12rRpdV588cXauR2bFSPYZ8j/5ptvKsbHxx88dOhQYmBgYGpqauot+SO36+UZXCFy6vB7EekJnDbG/JBjRdI6QPb/fE5ZZXmVK6WKWFJSEsOGDWPDhg10ad6cZePGUbdGjTsf6AR3ajkVpcjIyCuzZ8+uk5GRIX369Dnm4eHBlClTLmTt9/f3/3XHjh33jR49+pdZs2bV6ty589WAgIA7LujZrVu35GeffbbJ1KlTz/v4+GSeP3/evVatWpmdO3dOGTx4cIPY2Nhzxhg++eSTqsuWLbulRZfl+PHjZWrWrJkRExNzqWLFirbly5dXr1y5ctL169fzbNhcuXLFvXLlypkVK1a07d692+uHH364ryDfk8qVK2fe7vz3ymmJTETKA9OAx3LbnUuZuU15bucfDgwH8PX1vcsolVIAGzduZMjgwSQnJ/Pm0KGMefJJHZGYBy8vL9O2bdvkKlWqZHp43PoWO2TIkEuPPvqon6+vb0jLli2vr1y58nh+zhsWFvbrpEmTzrZv3z7Azc3NhISE3Pjb3/52/OGHH77Rv3//X1q2bBkIEBkZeaFdu3aphw4dynVI/65du8pNmTKlrpubGx4eHmbBggU/+/j4ZLZq1SrFz88vuHPnzlcnTpyYlP2YPn36XH3nnXe8/f39gxo3bvxr8+bNr+d27rzkPP/bb799qiDH34lDF9YUkQbAJmNMiIiEAl9hH8wBUBc4A7QGZgJbjDFrrOMOAZ2yHsaYaKv87ez18qILayp1d1JSUpg4cSKLFy/mgUaNWDlxIsH38o9h27ZQvXq+qhb1wppFJTMzk+Dg4KC//vWvR0NDQ9Pye9y5c+fcJ06cWGfbtm2VIiIiLr766qvnijJOV1MsF9Y0xiQAvw1LFZHjQJg1anEDMFpE3sc+2OOqMeasiHwO/FlEsua5eQyY4uDQlSoVduzYQWREBEePHWNynz7M7N+fsmXKODusYm3Xrl1evXr18uvevfvlgiQxsLdaVq9efeLONVVOjhx+vwZ7i6qGiJwCZhhj3s2j+ifYh94fwd5iiwIwxlwSkVhgp1VvljFGV3hVqhDdvHmT2NhYZs+eTT1vb775859pHxzs7LBcQqtWrX49depUgrPjKG0clsiMMS/cYX+DbNsGGJVHvaXA0kINTikFwKFDh4iMiGBnfDwDO3dm/vDhVCpfqKPYlSp0OmmwUgpjDIsWLWLSpEmUL1OG9ZMn06dtW2eHpVS+aCJTqpQ7d+4cg6Oi+PSzz+jWsiVLx46ldrVqzg5LqXzTRKZUKfbhhx8ybOhQbly/zv+OGMHI7t3J8ZlOpYo9/SCIUqVQcnIygwcPpnfv3jSsXp3v33iDmB49NIkpl6QtMqVKme3btxMZEcGJkyf54/PP88fnn6dMLh/cVcpVaItMqVIiPT2dKVOm0KFDB9zT09k+Zw6zwsM1iRUid3f3VlnzKgYEBATlNbvG3Xj++efr79q1ywv+s3RKTn369Gnw3nvvVc1tX37kdd7CsmrVqspTp071Kezz6m+wUqVAYmIiEeHh7N6zh6GPPcYbQ4ZQoVw5Z4dVtFavrlyo5+vf/+qdqnh6etqKaqmStWvX5jmhsKsIDw+/Ctzx+1hQ2iJTqgSz2WzMnz+fVq1acer4cT6aNo3Fo0eX/CRWjGRkZBAdHV03JCQk0N/fP+j111+vAXD16lW33/3ud/5BQUGB/v7+QStXrqwC9oU3cy58CdC6deumW7du/e1DfcOGDasbFBQU+Lvf/c7/zJkztzRKtm3bVv7BBx9sGhwcHPjwww/7/fzzz7dMy3Lw4MGyDzzwQEBISEjguHHj7s8qt9lsREdH1/Xz8wv29/cPyoph06ZNFR988MGmPXr0aNSgQYOQmJiYOgsXLqwWGhoa6O/vH7R//35PgNWrV1du1qxZQGBgYFDbtm39T5486QEwf/786gMGDPAFe+tx0KBB9Vq0aBFQt27d0HtpSWoiU6qEOn36NN0ef5xx48bxaGgoCfPm0bNNG2eHVaKlpaW5ZXUrdu3atTHAm2++WaNy5cqZ+/btO/DDDz8cWL58uffBgwfLli9f3vbxxx8fSUxMPPDNN9/8OHXq1Lo2m40PPvigko+Pz81Dhw4lHj58eH/v3r2Tc14nNTXVrWXLljcSExMPtGvX7trkyZPvzxGHjB071vejjz46un///gMDBw68+OKLL96yUkhMTIzv0KFDL+zbt++Aj4/PzazyFStWVElISCh34MCB/V999dWP06dPr5uVCA8ePFhu4cKFJw8cOLB//fr11X/88UevhISEA5GRkRfj4uJqAnTt2jVlz549Bw8cOJD47LPPXpo1a1au3Ynnz58vEx8ff/Cjjz46PGPGjLteyUS7FpUqgdatW8eIESNIS03l7ZgYhj3+uI5IdIDcuhY3b95c6eDBg+U3bNhQFeDatWvuiYmJXg0bNrw5fvz4ujt27Kjg5uZGUlJS2VOnTnm0bNkyNefClzmv4+bmxtChQy8BDB48+JfevXs3yb5/7969nocPHy7XuXNnf7C3sLy9vW/mPM/3339f4dNPPz0KEB0d/UtsbGxdgG3btlXs27fvJQ8PD+rVq5fRpk2blO3bt5evXLmyLTQ09Hr9+vVvAvj6+qZ17979KkDz5s1Tv/nmm4oAP/30U9mnn3667oULF8qkp6e71atXL9d5J3v27HnF3d2dVq1a/frLL7/c9USemsiUKkGuXLnCmDFjWLlyJW2aNmXlhAk0uf/+Ox+oiowxRuLi4k706dPnv1pW8+fPr/7LL794JCQkHPD09DR16tQJTU1Ndcta+PJvf/tb5WnTptXZvHlz8ty5c8/e7ho5/0kxxkiTJk1S9+zZc/BO8bm5uRVosU5PT8/fdrq5ueHl5WWytjMzMwVg9OjRvuPGjTsXHh5+ddOmTRVnzZqV6y9h1rF3uuYdX8NdH6mUKla2bNlCs9BQ1qxZw8z+/dk+Z44msWKga9euVxcuXOidlpYmYG8tJScnu129etW9Ro0aNz09Pc3GjRsrnjlzpizYF76sWLGiLSYm5tL48ePP79mz55bJLm02G1n3lJYtW1a9devW17Lvb9as2a+XLl3y2Lx5831g72qMj4/3ynmeli1bpixevLgawOLFi39bX6djx47X1q9fXy0jI4MzZ854fPfddxXat2+f7zXIrl275u7r63szK778Hne3tEWmlItLS0vjD3/4A3FxcTS5/37+9dprtPb3d3ZYyjJhwoSLx48f9wwNDQ00xki1atVufvLJJ0eHDh16qXv37k1CQkICg4ODbzRs2PBXyH3hy5znLFeunG3//v3lgoODfSpWrJj5wQcf/NeK0F5eXub9998/OnbsWN9r1665Z2ZmysiRI8+HhYX9mr3eggULTvTr16/RggULavXs2fNyVnlkZOSVf/3rXxUCAwODRcTMnDnzlK+vb8bevXvz9ZqnTZt25oUXXmhcq1at9LCwsOsnTpzwvKtvXj45dGFNZ9GFNVVJlZCQQHj//iTs28fI7t15PSqK+7xu+ce7+CgFC2uqonG7hTW1a1EpF2Sz2YiLiyMsLIwLZ87wyYwZLBg5sngnMaWKiHYtKuViTpw4wcCBA9myZQvPPPQQ74weTY1KlZwdllJOo4lMKRdhjGH16tWMGjWKzPR0lo4dy6AuXXRYvSr1NJEp5QIuXbpETEwMa9eupV1QEP83fjwNfQp9yjpXZ7PZbJLbcHLl2mw2mwC2vPZrIlOqmNu8eTODBg7k/Pnz/Dkykt/37o27u7uzwyqO9l24cCHI29v7qiazksNms8mFCxcqA/vyquOwRCYiS4EngSRjTIhVFgv0wp5pk4BBxpgzYu8rmQf0AG5Y5d9bxwwE/mCd9k/GmOWOeg1KOVJqaipTpkxh3rx5BNarx4a5c2nZuLGzwyq2MjIyhp47d27JuXPnQtCBbCWJDdiXkZExNK8KDht+LyIdgBRgRbZEVskYk2xtjwWCjDEjRKQHMAZ7ImsDzDPGtBGRakA8EAYYYBfQyhhz+dYr/ocOv1euZvfu3USEh5N44ABjn3qKOQMGUM6zSD+K4xhFOPxelV4O+6/FGLMVuJSjLPuULfdhT05gb6WtMHY7gCoiUht4HPjSGHPJSl5fAt2KPnqlHCMzM5M5c+bQpk0briQl8cXMmcwbNqxkJDGliojT75GJyGxgAPY1ah6xiusAJ7NVO2WV5VWulMs7fvw4AwYMYNu2bTzXrh2LYmKoVrGis8NSqthzej+yMWaaMaYesAoYbRXnNp7Y3Kb8FiIyXETiRST+woULhROsUkXAGMOyZcto1qwZP3z/Pf83YQJrf/97TWJK5ZPTE1k2q4E+1vYpoF62fXWBM7cpv4Ux5h1jTJgxJszb27sIwlXq3l28eJFnn32WqKgoWtavz95584h45BH9bJhSBeDURCYiftme9gSylhzYAAwQu4eAq8aYs8DnwGMiUlVEqgKPWWVKuZxPP/2U0JAQNm3cyOtRUXwVG0v9mjWdHZZSLseRw+/XAJ2AGiJyCpgB9BCRptiHV/4MjLCqf4J9xOIR7MPvowCMMZesIfs7rXqzjDH/NYBEqeLuxo0bvPTSSyxYsICQ+vX5fO5cmjVs6OywlHJZDktkxpgXcil+N4+6BhiVx76lwNJCDE0ph9m5cycR4eEcPnKESU8/zZ8iIvAqW9bZYSnl0orTPTKlSqyMjAxiY2Np27YtqVeu8FVsLHMHD9YkplQhcPrwe6VKuqNHjxIZEcG/d+wgvGNH3oqOpkqFCs4OS6kSQxOZUkXEGMOSJUuYMGECZUR4/6WXeL59e2eHpVSJo4lMqSKQlJTEsGHD2LBhA12aN2fZuHHUrVHD2WEpVSJpIlOqkG3cuJEhgweTnJzMm0OHMubJJ3Fz09vRShUV/etSqpCkpKQwfPhwevbsSZ1Kldj1l78wrmdPTWJKFTFtkSlVCHbs2EFkRARHjx3j5T59mNW/P2XLlHF2WEqVCvqvolL34ObNm0yfPp127dpxMyWFb/78Z+YMHKhJTCkH0haZUnfp0KFDREZEsDM+noGdOzN/+HAqlS/v7LCUKnU0kSlVQMYYFi1axKRJkyjn4cH6yZPp07ats8NSqtTSRKZUAZw7d47BUVF8+tlndGvZkqVjx1K7WjVnh6VUqaaJTKl8+vDDDxk2dCg3rl/nrehoYnr00OVWlCoGdLCHckmrElbR4M0GuM10o8GbDViVsKrIrpWcnMzgwYPp3bs3DapV4/s33mDUE09oElOqmNAWmXI5qxJWMXzjcG7cvAHAz1d/ZvjG4QCEh4YX6rW2b99OZEQEJ06e5A99+zK9Xz/KeOifjVLFibbIlMuZ9tW035JYlhs3bzDtq2mFdo309HSmTJlChw4dcE9PZ/ucOcRGRGgSU6oY0r9K5XJOXD1RoPKCSkxMJCI8nN179jD0scd4Y8gQKpQrVyjnVkoVPm2RKZfjW9m3QOX5ZbPZmD9/Pq1ateLU8eN8NG0ai0eP1iSmVDGniUy5nNldZlO+zH9/8Lh8mfLM7jL7rs95+vRpuj3+OOPGjePR0FAS5s2jZ5s29xqqUsoBtGtRuZysAR3TvprGiasn8K3sy+wus+96oMe6desYMWIEaampvB0Tw7DHH9cRiUq5EIclMhFZCjwJJBljQqyy14GngHTgKBBljLli7ZsCDAEygbHGmM+t8m7APMAdWGKMmeOo16CKj/DQ8HseoXjlyhXGjBnDypUradO0Kf83YQJ+999fSBEqpRzFkV2Ly4BuOcq+BEKMMc2AH4EpACISBPQDgq1jFoiIu4i4A/8LdAeCgBesukoVyJYtW2gWGsqaNWuY2b8/2+fM0SSmlIvKdyITkc0i0vxuL2SM2QpcylH2hTEmw3q6A6hrbfcC3jfGpBljfgKOAK2txxFjzDFjTDrwvlVXqXxJS0vjpZdeonPnzngZw79ee43p/frh4e7u7NCUUnepIF2LvwfeEJGfganGmLOFHMtgYK21XQd7YstyyioDOJmjXO/Iq3xJSEggvH9/EvbtY2T37rweFcV9Xl7ODkspdY/y3SIzxnxvjOkMbAI+E5EZIlIo45JFZBqQAWTNM5TbnXZzm/LczjlcROJFJP7ChQuFEaZyUTabjbi4OMLCwkg6fZqPp09nwciRmsSUKiEKdI9M7EO5DgELgTHAYRGJvJcARGQg9kEg4caYrKR0CqiXrVpd4Mxtym9hjHnHGBNmjAnz9va+lxCVCztx4gRdunThxRdfpEfLliTMn0+PsDBnh6WUKkRkMzK/AAAYkUlEQVQFuUe2HTgNvIG9m28Q0AloLSLv3M3FrRGILwM9jTHZ5xzaAPQTEU8RaQj4Ad8BOwE/EWkoImWxDwjZcDfXViWbMYZVq1bRrFkz4r/9lqVjx/LBlCl4V67s7NCUUoWsIPfIRgD7s7WasowRkQN3OlhE1mBPfDVE5BQwA/soRU/gS+tzOzuMMSOMMftFZB2QiL3LcZQxJtM6z2jgc+zD75caY/YX4DWoUuDSpUvExMSwdu1a2gUFsWL8eBr5+Dg7LKVUEZFb89JdnESkkTHmWCHEUyTCwsJMfHy8s8NQDrB582YGDRzI+fPnmdW/P7/v3Rt3HZFYfLRtC9Wr56uqiOwyxmg/sLqjQvkcWXFOYqp0SE1NZfz48XTt2pVK7u58O3cuU557TpOYUqWATlGlXN7u3buJCA8n8cABxj71FHMGDKCcp6ezw1JKOUhBBnu8lp8ypRwlMzOTOXPm0KZNG64kJfH5zJnMGzZMk5hSpUxBuha75lLWvbACUaogjh8/ziOPPMKUKVN4unVrEubP57EWLZwdllLKCe7YtSgiI4EYoJGI7M22qyLwz6IKTKncGGNYvnw5Y8eORWw2/m/CBMI7ddLZ6pUqxfJzj2w18CnwKjA5W/k1Y8yl3A9RqvBdvHiR6OhoPvjgAzqGhLB8/Hjq16zp7LCUUk52x0RmjLkKXAVeKPpwlMrdp59+yuCoKC5dusTrUVFM6NlTRyQqpYB83COzZvRARK6JSHK2xzURSS76EFVpduPGDUaNGkWPHj2o4eXFzrlzefGZZzSJKaV+k5+uxZ+sr380xrxZlMEold3OnTuJCA/nx8OHmdirF7MjI/EqW9bZYSmlipn8jFpsKSL1gSgRqSoi1bI/ijpAVfpkZGQQGxtL27ZtSb1yha9iY4kbMkSTmFIqV/lpkb0NfAY0Ar7Psc9Y5UoViqNHjxIZEcG/d+ygf8eO/G90NFUqVHB2WEqpYiw/gz3mA/NFZKExZqQDYlKlkDGGJUuWMGHCBMqIsObFF+nXoYOzw1JKuYB8T1FljBkpIs2B9lbRVmPM3tsdo1R+JCUlMWzYMDZs2ECX5s1ZNm4cdWvUcHZYSikXUZApqsZiX8G5pvVYJSJjiiowVTps3LiR0JAQPv/sM94cOpQvZs7UJKaUKpCCTBo8FGhjjLkOv82z+G/gf4oiMFWypaSkMHHiRBYvXswDjRrxjxkzCPb1dXZYSikXVJBEJkBmtueZVplSBbJjxw4iIyI4euwYL/fpw8z+/fEsU8bZYSmlXFRBEtl7wLci8qH1/Gng3cIPSZVUN2/eJDY2ltmzZ1PP25tv/vxn2gcHOzsspZSLy1ciE/uMrH8FtgAPY2+JRRljdhddaKokOXToEJEREeyMj2dg587MHz6cSuXLOzsspVQJkK9EZowxIvJ3Y0wrbv0smVJ5MsawaNEiJk2aRDkPD9ZPnkyftm2dHZZSqgQpyHpkO0Tkwbu9kIgsFZEkEdmXrew5EdkvIjYRCctRf4qIHBGRQyLyeLbyblbZERHJPhu/KmbOnTvHE088QUxMDB0DA9n3P/+jSUwpVegKco/sEWCEiBwHrmPvXjTGmGb5PH4Z8BawIlvZPqA39tlDfiMiQUA/IBi4H9gsIv7W7v/FvsjnKWCniGwwxiQW4HUoB/jwww8ZNnQo11NSeCs6mpgePXTNMKVUkShIIrun1aCNMVtFpEGOsgNAbm9wvYD3jTFpwE8icgRobe07Yow5Zh33vlVXE1kxkZyczPjx43nvvfdo1aQJK//0JwLq1nV2WEqpEiw/K0R7ASOAJkAC8K4xJqOI46oD7Mj2/JRVBnAyR3mbIo5F5dP27duJjIjgxMmT/KFvX6b360cZj4L8r6SUUgWXn3tky4Ew7EmsOxBXpBHZ5dYHZW5TfusJRIaLSLyIxF+4cKFQg1P/LT09nalTp9KxY0fc0tPZ9uqrxEZEaBJTSjlEft5pgowxoQAi8i7wXdGGBNhbWvWyPa8LnLG28yr/L8aYd4B3AMLCwnJNdureJSYmEhEezu49exj62GP8ZfBgKuqweqWUA+WnRXYza8MBXYpZNgD9RMRTRBoCftgT6E7AT0QaikhZ7ANCNjgoJpWNzWZj/vz5tGrVilPHj/P3qVNZPHq0JjGllMPlp0XWXESSrW0BylnPs0YtVsrPhURkDdAJqCEip4AZwCXsczV6Ax+LyB5jzOPGmP0isg77II4MYJQxJtM6z2jgc8AdWGqM2Z/P16oKyenTp4kaNIgvN2/myQcfZMno0dSqWtXZYSmlSikxpuT3uoWFhZn4+Hhnh1EirFu3jhEjRpCWmsobgwcz7PHHdVi9yr+2baF69XxVFZFdxpiwO9dUpV1BPhCtSrErV64QGRnJ888/j3/Nmux5802Gd+umSUwp5XQ6rEzd0ZYtWxgQGcmZs2eZ2b8/U597Dg93d2eHpZRSgLbI1G2kpaXx0ksv0blzZ7yM4V+vvcb0fv00iSmlihVtkalcJSQkEN6/Pwn79jGye3dej4riPi8vZ4ellFK30BaZ+i82m424uDjCwsJIOn2aj6dPZ8HIkZrElFLFlrbI1G9OnDjBoEGD+Prrr3n6oYd4Z9QovCtXdnZYSil1W5rIFMYYVq9ezahRo8hMT2fp2LEM6tJFRyQqpVyCJrJS7tKlS8TExLB27VraBQWxYvx4Gvn4ODsspZTKN01kpdjmzZsZNHAg58+fZ3ZEBC/36YO7jkhUSrkYHexRCqWmpjJ+/Hi6du1KJXd3vp07l6l9+2oSU0q5JG2RlTK7d+8mIjycxAMHGPPkk7w2cCDlPD2dHZZSSt01bZGVEpmZmcyZM4c2bdpwJSmJz2fOZP7w4ZrElFIuT1tkpcDx48cZMGAA27Zt47l27VgUE0O1ihWdHZZSShUKTWQlmDGG5cuXM3bsWMRmY8WECUR06qTD6pVSJYomshLq4sWLREdH88EHH9AxJITl48dTv2ZNZ4ellFKFThNZCfTZZ58RNWgQv/zyC/9v0CAm9uqlIxKVUiWWDvYoQW7cuMGoUaPo3r07Nby82BkXx0u9e2sSU0qVaNoiKyF27txJRHg4Px4+zMRevZgdGYlX2bLODksppYqctshcXEZGBrGxsbRt25Ybly/zVWwscUOGaBJTSpUaDktkIrJURJJEZF+2smoi8qWIHLa+VrXKRUTmi8gREdkrIi2zHTPQqn9YRAY6Kv7i6OjRo3Ro357p06fTt1079s6fT+fmzZ0dllJKOZQjW2TLgG45yiYDXxlj/ICvrOcA3QE/6zEcWAj2xAfMANoArYEZWcmvNDHGsHjxYpo3b86BfftY8+KLrJo0iaoVKjg7NKWUcjiHJTJjzFbgUo7iXsBya3s58HS28hXGbgdQRURqA48DXxpjLhljLgNfcmtyLNGSkpJ4+umnGT58OG0aNyZh/nz6dejg7LCUUsppnD3Yo5Yx5iyAMeasiGR90KkOcDJbvVNWWV7lpcLGjRsZOmQIV69e5Y0hQxj71FO4ueltTqVU6VZc3wVzm3rC3Kb81hOIDBeReBGJv3DhQqEG52gpKSkMHz6cnj17UrtiReLj4hjfq5cmMaWUwvmJ7LzVZYj1NckqPwXUy1avLnDmNuW3MMa8Y4wJM8aEeXt7F3rgjrJjxw5aPPAAS5Ys4eU+ffj29dcJqV/f2WEppVSx4exEtgHIGnk4EPgoW/kAa/TiQ8BVqwvyc+AxEalqDfJ4zCorcW7evMmMGTN4+OGHuZmSwpbZs5kzcCCeZco4OzSllCpWHHaPTETWAJ2AGiJyCvvowznAOhEZApwAnrOqfwL0AI4AN4AoAGPMJRGJBXZa9WYZY3IOIHF5hw4dIjIigp3x8Qzs3Jl5w4ZR+b77nB2WUkoVSw5LZMaYF/LY1SWXugYYlcd5lgJLCzG0YsMYw6JFi5g0aRLlPDxYP3kyfdq2dXZYSilVrDl71KKynDt3jsGDB/Ppp5/yeMuWLB0zhvurV3d2WEopVexpIisGPvzwQ4YNHcr1lBTeio4mpkcPXTNMKaXyydmDPUq15ORkBg8eTO/evWlQrRq733yTUU88oUlMKaUKQFtkTrJ9+3YiIyI4cfIkf+jbl+n9+lHGQ38cSilVUNoic7D09HSmTp1Kx44dcUtPZ9urrxIbEaFJTCml7pK+ezpQYmIiEeHh7N6zh6GPPcZfBg+mYvnyzg5LKaVcmiYyB7DZbLz11lu8/PLLVPD05O9Tp9LroYecHZZSSpUImsiK2OnTp4kaNIgvN2/mibAw3h0zhlpVS93KM0opVWQ0kRWhdevWMWLECNJSU3k7JoZhjz+uIxKVUqqQ6WCPInDlyhUiIyN5/vnn8a9Zkz1vvsnwbt00iSmlVBHQFlkh27JlCwMiIzlz9iyvvPAC0/r2xcPd3dlhKaVUiaUtskKSlpbGSy+9ROfOnfEyhn+99hozXnhBk5hSShUxbZEVgoSEBML79ydh3z5Gdu/O61FR3Ofl5eywlFKqVNAW2T2w2WzExcURFhZG0unTfDx9OgtGjtQkppRSDqQtsrt04sQJBg0axNdff83TDz3EO6NG4V25srPDUkqpUkcTWQEZY1i9ejWjRo0iMz2dpWPHMqhLFx2RqJRSTqKJrAAuX77MyJEjWbt2Le2CglgxfjyNfHycHZZSSpVqmsjyafPmzQwaOJDz588zOyKCl/v0wV1HJCqllNPpYI87SE1NZfz48XTt2pVK7u58O3cuU/v21SSmlFLFRLFIZCIyTkT2ich+ERlvlVUTkS9F5LD1tapVLiIyX0SOiMheEWlZVHHt3r2bsFatmDdvHmOefJJdf/kLLRs3LqrLKaWUugtOT2QiEgIMA1oDzYEnRcQPmAx8ZYzxA76yngN0B/ysx3BgYWHHlJmZyZw5c2jTpg2Xz5/ns1deYf7w4ZTz9CzsSymllLpHxeEeWSCwwxhzA0BEvgGeAXoBnaw6y4EtwMtW+QpjjAF2iEgVEaltjDlbGMEcP36cAQMGsG3bNp5r146FI0dSvVKlwji1UkqpIlAcEtk+YLaIVAdSgR5APFArKzkZY86KSE2rfh3gZLbjT1ll95zI3n//fYYPH47YbKyYMIGITp10WL1SShVzTk9kxpgDIvIa8CWQAvwAZNzmkNwyi7mlkshw7F2P+Pr65isWd3d3WrRowYply6hfv36+jlFKFYD+Y6iKgNh76IoPEfkz9lbWOKCT1RqrDWwxxjQVkbet7TVW/UNZ9fI6Z1hYmImPj8/X9W02G25uTr91qFSpJyK7jDFhzo5DFX/F4h07q9tQRHyB3sAaYAMw0KoyEPjI2t4ADLBGLz4EXC2s+2OAJjGllHIxTu9atPzNukd2ExhljLksInOAdSIyBDgBPGfV/QT7fbQjwA0gyhkBK6WUKh6KRSIzxrTPpewXoEsu5QYY5Yi4lFJKFX/aj6aUUsqlaSJTSinl0jSRKaWUcmmayJRSSrk0TWRKKaVcmiYypZRSLk0TmVJKKZemiUwppZRL00SmlFLKpWkiU0op5dI0kSmllHJpmsiUUkq5NE1kSimlXJomMqWUUi5NE5lSSimXpolMKaWUS9NEppRSyqVpIlNKKeXSNJEppZRyacUikYnIBBHZLyL7RGSNiHiJSEMR+VZEDovIWhEpa9X1tJ4fsfY3cG70SimlnMnpiUxE6gBjgTBjTAjgDvQDXgPeMMb4AZeBIdYhQ4DLxpgmwBtWPaWUUqWU0xOZxQMoJyIeQHngLNAZWG/tXw48bW33sp5j7e8iIuLAWJVSShUjTk9kxpjTwFzgBPYEdhXYBVwxxmRY1U4BdaztOsBJ69gMq351R8aslFKq+PBwdgAiUhV7K6shcAX4K9A9l6om65Db7Mt+3uHAcOtpiogcymdINYCL+azrCBrP7RW3eKD4xeSq8dQv6kBUyeD0RAY8CvxkjLkAICIfAG2BKiLiYbW66gJnrPqngHrAKasrsjJwKedJjTHvAO8UNBgRiTfGhN3VKykCGs/tFbd4oPjFpPGoks7pXYvYuxQfEpHy1r2uLkAi8DXwrFVnIPCRtb3Beo61/x/GmFtaZEoppUoHpycyY8y32AdtfA8kYI/pHeBlYKKIHMF+D+xd65B3gepW+URgssODVkopVWwUh65FjDEzgBk5io8BrXOp+yvwXBGGU+DuyCKm8dxecYsHil9MGo8q0UR75ZRSSrkyp3ctKqWUUvdCE5lFRGJFZK+I7BGRL0TkfqtcRGS+NSXWXhFp6aB4XheRg9Y1PxSRKtn2TbHiOSQijzsonuesacRsIhKWY5/D47Gu28265hERcfi9UhFZKiJJIrIvW1k1EfnSmlrtS+vjJY6Kp56IfC0iB6yf1ThnxmRNNfediPxgxTPTKs91+jml7poxRh/27tVK2bbHAous7R7Ap9g/v/YQ8K2D4nkM8LC2XwNes7aDgB8AT+yfvTsKuDsgnkCgKbAF+3RiODked+tajYCyVgxBDv6d6QC0BPZlK/t/wGRre3LWz81B8dQGWlrbFYEfrZ+PU2Ky/mYqWNtlgG+tv6F1QD+rfBEw0pE/N32UvIe2yCzGmORsT+/jPx+y7gWsMHY7sH++rbYD4vnC/Gdmkx3YP0uXFc/7xpg0Y8xPwBFyGRRTBPEcMMbk9qFyp8RjXeOIMeaYMSYdeN+KxWGMMVu59TOM2adQyz61miPiOWuM+d7avgYcwD4TjlNisv5mUqynZayHIe/p55S6K5rIshGR2SJyEggHplvFv02JZck+XZajDMbeKiwu8WTnrHiK2/chSy1jzFmwJxagpjOCsFaFaIG9FeS0mETEXUT2AEnAl9hb0XlNP6fUXSlViUxENltLxeR89AIwxkwzxtQDVgGjsw7L5VSFMtTzTvFYdaYBGVZMTo8nt8OKKp47cNZ1iz0RqQD8DRifo6fB4YwxmcaYB7D3KLTG3kV9SzXHRqVKmmLxOTJHMcY8ms+qq4GPsX+2LWtKrCzZp8sq0nhEZCDwJNDFGJP1x+60ePJQZPEU0+veyXkRqW2MOWt1QSc58uIiUgZ7EltljPmgOMQEYIy5IiJbsN8jy2v6OaXuSqlqkd2OiPhle9oTOGhtbwAGWKMXHwKuZnXTFHE83bDPbtLTGHMj264NQD+xLzDaEPADvivqeG7DWfHsBPysEXBlsa9ht8EB172T7FOoZZ9archZU7y9CxwwxvzF2TGJiHfWaFsRKYd9XtUD5D39nFJ3x9mjTYrLA/t/sfuAvcBGoI5VLsD/Yu/bTyDbiL0ijucI9ntAe6zHomz7plnxHAK6OyieZ7C3gtKA88DnzozHum4P7CPzjgLTnPA7swb70kM3re/NEOzTqX0FHLa+VnNgPA9j76bbm+33poezYgKaAbutePYB063yRtj/2TmCfbULT0f/7PRRsh46s4dSSimXpl2LSimlXJomMqWUUi5NE5lSSimXpolMKaWUS9NEppRSyqVpIlNKKeXSNJEppZRyaZrI1D0RkUxrDbd9IvJXESlfDGKqIiIxhXAedxGZZ62llSAijQojPqVU4dJEpu5VqjHmAWNMCJAOjMjPQdaUX0X1+1cFKFAiyyOeKcAxY0wwML+g51RKOYYmMlWYtgFNAETk7yKyy2rNDLfKGlirFy8Avgfq3abeQRFZYrX0VonIoyLyT2tV4d/WOxORCGsV4j0i8raIuANzgMZW2et51cstnmznvQ94xhgzzyr6Keu1KaWKF52iSt0TEUkxxlQQEQ/s81V+ZoxZKCLVjDGXrMlidwIdsa9afAxoa+yLlHKbekewr6e13yr/Aftchj2BKGPM0yISiH31497GmJtWQtoBbAU2Wa1E7lDvv+LJ9rp6AYv5z8zs1YDNxpjBhfwtVErdo1K1jIsqEuWshRPB3iJ719oeKyLPWNv1sM+Kfw74OUfSyKveT8aYBAAR2Q98ZYwxIpIANLDqdwFaATvtE79TDvsSJVtzxHi7ejnjyfIA9kluF1kxLAH2WvfJpgGVjTHP5nKcUsrBNJGpe5Vq7Asn/kZEOmFfsuN3xpgb1jpUXtbu6/msl5btlLZsz2385/dWgOXGmCk5rt8gR4y3q3ed3FXF3p2I1dp8DJhtjPkJGCIi6/M4TinlYHqPTBWFysBlKzkFYF9M8V7q5eUr4FkRqQn2bkoRqQ9cw949ead6t/NjtngmAB9bSUwpVcxoIlNF4TPAQ0T2ArHY70fdS71cGWMSgT8AX1jn+BKobYz5BfinNVDk9bzq3eH0a4CWInIE+7paEwsSm1LKcXSwh1IFICLVgdlAV2CJMeZVJ4ekVKmniUwppZRL065FpZRSLk0TmVJKKZemiUwppZRL00SmlFLKpWkiU0op5dI0kSmllHJpmsiUUkq5NE1kSimlXJomMqWUUi7t/wMhRbdqCG1dRgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(0, 1230, 'o', color='green', label='Default optimum')\n",
    "\n",
    "y_fill = np.arange(900, 1440, 1)#900 <= y <= 1440\n",
    "x_fill = (5/36*y_fill - 325/2)*(4/5)# theta1 >= (5/36*y - 325/2)*4/5\n",
    "ax.plot(x_fill, y_fill, color='black', label=r'y-$\\theta_1$ constraint')\n",
    "plt.fill_between(x_fill, y_fill, 900, where=y_fill<=1440, facecolor='red',\n",
    "                 alpha=0.3, label='Feasible domain')\n",
    "ax.set_ylim(800,1500)\n",
    "ax.set_xlabel(r'Parameter $ \\theta_1 $')\n",
    "ax.set_ylabel(r'Profit $ y $')\n",
    "\n",
    "plt.subplots_adjust(left = 0.20, right = 0.75, bottom=0.2)\n",
    "plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', \n",
    "           borderaxespad=0, fontsize=10)\n",
    "plt.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}