Logistic regression initialization.ipynb

init_logreg.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Logistic regression initialization.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "authorship_tag": "ABX9TyPAUZSyCkrzK+/4xnXOkJpj",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/gist/avidale/a640f7a8e353d9efdd79385e277caef1/init_logreg.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TNhcEdNxBn14"
      },
      "source": [
        "The problem: logistic regression, one of the simplest and fastest classifiers, still trains relatively slowly, because it uses gradient descent for training.\r\n",
        "\r\n",
        "The speed with which gradient descent converges depends on initialization: if we start descent from a good initial point, we will have a shorter way to go.\r\n",
        "\r\n",
        "So what is a way to initialize a logistic regression which is close to an optimal solution and at the same time fast to converge? Our suggestion is that linear regression is a good candidate. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nDmx2NVpoC8j"
      },
      "source": [
        "Our idea is to initialize logistic regression $f$ with such parameters $\\alpha, \\beta$ that a linear regression $g$ trained on the same data is tangent to the logistic regression in the center of the dataset $\\bar{X}, \\bar{Y}$ (linear regression always goes through it, so that $g(\\bar{X})=\\bar{Y}$).\r\n",
        "\r\n",
        "Linear regression: $$g(X) = a + b X$$\r\n",
        "Logistic regression: $$f(X) = \\sigma(\\alpha + \\beta X)$$\r\n",
        "where $\\sigma(t) = \\frac{1}{1+e^{-t}}$.\r\n",
        "\r\n",
        "\"Tangent\" means \"values are equal and first derivatives are equal at this point\", so we want\r\n",
        "$$ g(\\bar{X}) = f(\\bar{X})$$\r\n",
        "$$ g'(\\bar{X}) = f'(\\bar{X})$$\r\n",
        "\r\n",
        "Or equivalently\r\n",
        "$$ \\bar{Y}= f(\\bar{X})$$\r\n",
        "$$ b = f(\\bar{X}) (1-f(\\bar{X}))\\beta$$\r\n",
        "\r\n",
        "From the last equation we can immediately express the slope coefficient\r\n",
        "$$\\beta = \\frac{1}{\\bar{Y}(1-\\bar{Y})}b$$\r\n",
        "We can use this value to calculate the intercept \r\n",
        "$$ \\alpha = \\sigma^{-1}(\\bar{Y}) - \\beta \\bar{X}$$\r\n",
        "\r\n",
        "The inverse logistic function $\\sigma^{-1}(t) = \\log\\frac{t}{1-t}$ is called logit. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "cAUoZjlwr3tp"
      },
      "source": [
        "We start with a small generared dataset. To train linear regression, we use Ridge, so that we don't have to worry about multicollinearity. \r\n",
        "\r\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "88UUviBKkMyq"
      },
      "source": [
        "import numpy as np\r\n",
        "import matplotlib.pyplot as plt\r\n",
        "\r\n",
        "from sklearn import datasets\r\n",
        "from sklearn.linear_model import LogisticRegression, Ridge\r\n",
        "from sklearn.model_selection import train_test_split\r\n",
        "\r\n",
        "def logit(x):\r\n",
        "    return np.log(x/(1-x))\r\n",
        "\r\n",
        "\r\n",
        "# Create dataset of classification task with many redundant and few\r\n",
        "# informative features\r\n",
        "X, y = datasets.make_classification(n_samples=100, n_features=1, n_informative=1, n_redundant=0, random_state=42, n_clusters_per_class=1)\r\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 279
        },
        "id": "U0Q7KW5Hkg6b",
        "outputId": "b9a91372-3deb-4274-b559-52759cd29020"
      },
      "source": [
        "plt.scatter(X[:,0], y, s=3, c='k')\r\n",
        "g = np.linspace(X[:,0].min(), X[:,0].max())\r\n",
        "linreg = Ridge(alpha=0.01).fit(X, y)\r\n",
        "plt.plot(g, linreg.predict(g[:, np.newaxis]))\r\n",
        "log_init = LogisticRegression()\r\n",
        "ym = y.mean()\r\n",
        "log_init.coef_ = linreg.coef_[np.newaxis, :] / (ym * (1-ym))\r\n",
        "log_init.intercept_ = logit(ym)[np.newaxis] - np.dot(X.mean(axis=0), log_init.coef_)\r\n",
        "log_init.classes_ = np.array([0, 1])\r\n",
        "plt.plot(g, log_init.predict_proba(g[:, np.newaxis])[:, 1])\r\n",
        "log_final = LogisticRegression().fit(X, y)\r\n",
        "plt.plot(g, log_final.predict_proba(g[:, np.newaxis])[:, 1])\r\n",
        "\r\n",
        "plt.xlabel('X'); plt.ylabel('y'); plt.legend([ 'linear model', 'initial logistic model', 'final logistic model', 'data']);"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd3xN9//A8dfnZsoQQcyExAySIGKvqFlFS+vbqhptUVWly2i/NaqqaH/VQbW0pTpoi6pS1KzaESOJvYLEyhCyx72f3x838k0IglxZ72d7Hvfecz7nc943ifu+53zOeR+ltUYIIYS4maGgAxBCCFE4SYIQQgiRK0kQQgghciUJQgghRK4kQQghhMiVdUEHkJ/Kly+vPT09CzoMIYQoMoKDg6O11m65LStWCcLT05O9e/cWdBhCCFFkKKXO3m6ZHGISQgiRK0kQQgghciUJQgghRK6K1RhEbtLT04mIiCAlJaWgQxGFjL29Pe7u7tjY2BR0KEIUSsU+QURERODs7IynpydKqYIORxQSWmtiYmKIiIjAy8uroMMRolAq9oeYUlJSKFeunCQHkYNSinLlysmepRB3UOwTBCDJQeRK/i6EuLMSkSCEEKK4CgqP5at/Tlmkb0kQD4GTkxMAFy5c4KmnnirgaO7fli1b6NGjxwO3EUI8uITUDCb+EUbfr3by8+5zJKVl5Ps2iv0gdWFSpUoVli5datFtZGRkYG0tv1YhirN/jkfxzvJQLlxLZnArT8Z0rYuDbf7/u5c9iIcoPDwcHx8fABYuXEifPn3o1q0btWvXZuzYsVnt/v77b1q2bIm/vz99+/YlISEBgClTptC0aVN8fHwYNmwYN+4GGBgYyGuvvUZAQACfffZZjm1OnjyZQYMG0bZtW6pXr87y5csZO3Ysvr6+dOvWjfT0dAA2btxI48aN8fX15YUXXiA1NRWAtWvX4u3tjb+/P8uXL8/qNzExkRdeeIFmzZrRuHFj/vjjD8v94IQQAFxNTOONXw8w6Ls9lLK1YunwVkzu1QBHO8t8KSxRXzXf+/MQhy9cz9c+61cpzaSeDe5r3QMHDrB//37s7OyoW7cur776KqVKlWLq1Kls2LABR0dHZsyYwSeffMLEiRMZOXIkEydOBGDAgAGsWrWKnj17ApCWlnbbOlSnTp1i8+bNHD58mJYtW7Js2TJmzpxJ7969Wb16Nd26dWPw4MFs3LiROnXqMHDgQObOncvw4cMZOnQomzZtolatWjz99NNZfX7wwQc88sgjfPfdd8TFxdGsWTM6dep0Xz8HIcSdaa1ZE3aJiX+EEZeUzqhHavHKI7Wws7ay6HZLVIIobDp27IiLiwsA9evX5+zZs8TFxXH48GFat24NmD/4W7ZsCcDmzZuZOXMmSUlJxMbG0qBBg6wEkf3D+2aPPvooNjY2+Pr6YjQa6datGwC+vr6Eh4dz7NgxvLy8qFOnDgCDBg1izpw5BAYG4uXlRe3atQF47rnnmDdvHmDey1m5ciUff/wxYD6d+Ny5c/n9IxKixLtyPYUJf4Sx7tBlfKu6sOiF5tSvUvqhbLtEJYj7/aZvKXZ2dlnPraysyMjIQGtN586dWbx4cY62KSkpjBgxgr179+Lh4cHkyZNznMPv6Oh41+0YDAZsbGyyTu80GAxkZNzfwJbWmmXLllG3bt0c8y9fvnxf/QkhctJa81twBFNXHSY1w8Tbj3rzYhsvrK0e3siAxbaklPpOKXVFKRV2m+WBSqlrSqkDmdPEbMu6KaWOKaVOKqXGWyrGwqhFixZs376dkydPAuZj/cePH89KBuXLlychISFfB7vr1q1LeHh41jZ/+OEH2rdvj7e3N+Hh4Zw6ZT6FLnvS6tq1K1988UXWOMj+/fvzLR4hSrrzsUkM+HYPY5eG4F25NGtGt+Wl9jUfanIAy+5BLARmA4vu0OZfrXWOcyKVUlbAHKAzEAEEKaVWaq0PWyrQwsTNzY2FCxfSr1+/rIHiqVOnUqdOHYYOHYqPjw+VKlWiadOm+bZNe3t7FixYQN++fcnIyKBp06YMHz4cOzs75s2bx2OPPYaDgwNt27YlPj4egAkTJvDaa6/h5+eHyWTCy8uLVatW5VtMQpRERpPm+x3hfLTuGFYGxdQnfHi2WTUMhoK5qFPd+AZokc6V8gRWaa19clkWCLyVS4JoCUzWWnfNfP02gNb6w7ttLyAgQN88UHvkyBHq1at3n+9AFHfy9yEKixOX4xm7LIT95+IIrOvGtN6+VClTyuLbVUoFa60DcltW0GMQLZVSB4ELmJPFIaAqcD5bmwig+e06UEoNA4YBVKtWzYKhCiFE/kvLMPHVP6eYvekkjnZWzHq6IU80qlooSsEUZILYB1TXWicopboDK4Da99qJ1noeMA/MexD5G6IQQlhOSEQcY5eGcPRSPD38KjO5VwPKO9ndfcWHpMAShNb6erbnfymlvlRKlQciAY9sTd0z5wkhRLGQkm5k1vrjzP/3NG7Odswb0IQuDSoVdFi3KLAEoZSqBFzWWmulVDPMZ1TFAHFAbaWUF+bE8AzwbEHFKYQQ+WnX6RjGLwshPCaJfs08GP9oPVxKFc6bVlksQSilFgOBQHmlVAQwCbAB0Fp/BTwFvKyUygCSgWe0ecQ8Qyk1ElgHWAHfZY5NCCFEkRWfks70NUf5afc5qpV14KchzWldq3xBh3VHFksQWut+d1k+G/NpsLkt+wv4yxJxCSHEw7bp6GX++3sYl6+nMKSNF290qWOR4nr5TYr1PQStWrW6a5shQ4Zw+LD5Uo9p06bd8/o3SorndX5eZI8pNwsXLuTChQt5bm8p2YsgPkgbIfJbbGIary3ZzwsL9+Jsb82yl1vxbo/6RSI5QMGf5loi7Nix465tvvnmm6zn06ZN45133rmn9S0he0y5WbhwIT4+PlSpUiVP7YUoKbTW/BlykckrDxGfks7ojrV5pUMtbK2L1nfyohVtEXXjW/yWLVsIDAzkqaeewtvbm/79++co2b13717Gjx9PcnIyjRo1on///jnWT0hIoGPHjvj7++Pr63tPJba11owZMwYfHx98fX355ZdfADCZTIwYMQJvb286d+5M9+7ds8p43IjJaDQyePDgrHVnzZrF0qVL2bt3L/3796dRo0YkJydntQdzmXB/f38aNmxIx44db4ln4cKFPPHEE3Tu3BlPT09mz57NJ598QuPGjWnRogWxsbGAueJtixYt8PPzo3fv3ly9ehWA4OBgGjZsSMOGDZkzZ05Wv0ajkTFjxtC0aVP8/Pz4+uuv8/6LEiIfXLqWwtBFwYxavB8P11L8+WobXu9cp8glByhpexBrxsOl0Pzts5IvPDo9z83379/PoUOHqFKlCq1bt2b79u20adMma/n06dOZPXs2Bw4cuGVde3t7fv/9d0qXLk10dDQtWrSgV69eebqgZvny5Rw4cICDBw8SHR1N06ZNadeuHdu3byc8PJzDhw9z5coV6tWrxwsvvJBj3QMHDhAZGUlYmLmsVlxcHGXKlGH27Nl8/PHHBATkvAgzKiqKoUOHsnXrVry8vLI+7G8WFhbG/v37SUlJoVatWsyYMYP9+/fz+uuvs2jRIl577TUGDhzIF198Qfv27Zk4cSLvvfcen376Kc8//zyzZ8+mXbt2jBkzJqvPb7/9FhcXF4KCgkhNTaV169Z06dKlUFx0JIo3rTVLgs4zbfUR0k0m/tu9Hi+08cKqgMpk5Ieil9KKuGbNmuHu7o7BYKBRo0aEh4fneV2tNe+88w5+fn506tSJyMjIPFdP3bZtG/369cPKyoqKFSvSvn17goKC2LZtG3379sVgMFCpUiU6dOhwy7o1atTg9OnTvPrqq6xdu5bSpe9canjXrl20a9cOLy8vAMqWLZtruw4dOuDs7IybmxsuLi5ZpctvlCG/du0acXFxtG/fHjCXId+6dStxcXHExcXRrl07wHxvjBv+/vtvFi1aRKNGjWjevDkxMTGcOHEiTz8jIe7X2ZhEnp2/m7eXh9KgamnWjm7H0HY1inRygJK2B3EP3/QtJbcS33n1008/ERUVRXBwMDY2Nnh6euYo+W0prq6uHDx4kHXr1vHVV1/x66+/8t133z1wv9l/FgaDIUdZ8gcpQ/7FF1/QtWvXHPPvJRELkVdGk2bB9jN8/PcxbAwGPuzjyzNNPYrNHqvsQRRCNjY2WbcCze7atWtUqFABGxsbNm/ezNmzZ/PcZ9u2bfnll18wGo1ERUWxdetWmjVrRuvWrVm2bBkmk4nLly+zZcuWW9aNjo7GZDLx5JNPMnXqVPbt2weAs7NzVnXX7Fq0aMHWrVs5c+YMwG0PMd2Ni4sLrq6u/Pvvv8D/ypCXKVOGMmXKsG3bNsCcOG/o2rUrc+fOzfr5HT9+nMTExPvavhB3cuxSPH3m7mDq6iO0rlmev99oR79m1YpNcoCStgdRRAwbNgw/Pz/8/f1zfPj179+fnj174uvrS0BAAN7e3nnus3fv3uzcuZOGDRuilGLmzJlUqlSJJ598ko0bN1K/fn08PDzw9/fPusvdDZGRkTz//POYTCYAPvzQXFh38ODBDB8+nFKlSrFz586s9m5ubsybN48+ffpgMpmoUKEC69evv6+fxffff8/w4cNJSkqiRo0aLFiwAIAFCxbwwgsvoJSiS5cuWe2HDBlCeHg4/v7+aK1xc3NjxYoV97VtIXKTlmHiyy0nmbP5JM72NnzerzE9/SoXq8Rwg0XLfT9sUu77/iQkJODk5ERMTAzNmjVj+/btVKpU+OrCWIL8fYh7cfC8ubjescvxPN6oCpN6NqCso21Bh/VACnO5b1EI9OjRg7i4ONLS0pgwYUKJSQ5C5FVympFP1h/j221nqOBsz7eDAuhYr2JBh2VxkiBEruMOQgiznadiGL88hLMxSTzbvBrjH/WmtH3hLK6X3yRBCCFELq6npPPhX0dZvOcc1cs5sHhoC1rWLFfQYT1UkiCEEOImG4+Yi+tdiU9hWLsavN6pDqVsrQo6rIdOEoQQQmSKSUjlvT8Ps/LgBbwrOfP1gCY09ChT0GEVGEkQQogST2vNyoMXmLzyEAmpGbzeqQ4vB9YskvWT8lPJfvcPyeeff069evXo378/K1euZPr0+7+iW8p6S1lvkb8uxCXz4vd7Gb3kANXLObJ6VFtGd6pd4pMDyB7EQ/Hll1+yYcMG3N3dAejVq1cBR5STlPUWJZHJpPl5zzmmrzmK0aSZ0KM+g1t5Fvn6SflJUqSFDR8+nNOnT/Poo48ya9YsFi5cyMiRIwHzlcijRo2iVatW1KhRI6vMtpT1lrLewrLORCfSb/4u3l0RRkMPF9a91o4Xi3jlVUsoUXsQM/bM4Gjs0Xzt07usN+Oajbvt8q+++oq1a9eyefNmypcvz8KFC3Msv3jxItu2bePo0aP06tWLp556Ssp6S1lvYSEZRhPfbjvDJ+uPY2ttYHofX54uRsX18pvsQRSwJ554AoPBQP369bNKd0tZbynrLfLfkYvX6TN3Bx+uOUq7Om5seKM9zxSz4nr5zWJ7EEqp74AewBWt9S2jhkqp/sA4QAHxwMta64OZy8Iz5xmBjNvVCblXd/qmX1Cyl7y+URdLynrfnZT1FnmVmmFkzqaTfLnlFGUcbJjzrD/dfStJYsgDS+5BLAS63WH5GaC91toXeB+Yd9PyDlrrRvmVHIoSKestZb1F/th37io9Pt/G55tO0qthFda/3p7HimnlVUuw2B6E1nqrUsrzDst3ZHu5C3C3VCxFjZT1lrLe4sEkpWXw8brjLNhxhsql7VnwfFM61K1Q0GEVORYt952ZIFbldojppnZvAd5a6yGZr88AVwENfK21vnnvIvu6w4BhANWqVWty87dtKef8PyW5rPftyN9H8bP9ZDTjl4dwPjaZgS2rM7abN052Jep8nHtSqMt9K6U6AC8CbbLNbqO1jlRKVQDWK6WOaq235rZ+ZvKYB+b7QVg84CJMynqL4uxacjrTVh/hl73n8SrvyK8vtaSZV+4nToi8KdAEoZTyA74BHtVax9yYr7WOzHy8opT6HWgG5JogRN5JWW9RXK07dIkJK8KISUzj5cCajO5YG3ubkldcL78VWIJQSlUDlgMDtNbHs813BAxa6/jM512AKQUUphCiEIuKT2XyykOsDr1Ivcql+XZQU3zdXe6+osgTS57muhgIBMorpSKASYANgNb6K2AiUA74MvOMghuns1YEfs+cZw38rLVea6k4hRBFj9aa3/dHMmXVYZJSjYzpWpdh7WpgYyWXduUnS57F1O8uy4cAQ3KZfxpoaKm4hBBFW2RcMu8sD+Wf41H4VyvDzKf8qFXBuaDDKpYKfJBaCCHywmTS/Lj7LDPWHEUDk3vWZ0BLKa5nSZIgHrLJkyfj5OTEW2+9levyFStWUKdOHerXr/+QIxOi8DoVlcD4ZSEEhV+lbe3yTOvti0dZh4IOq9iTBFHIrFixgh49ekiCEAJzcb15/57m0w0nsLc28NFTfjzVxF2uhH5IZETnIfjggw+oU6cObdq04dixYwDMnz+fpk2b0rBhQ5588kmSkpLYsWMHK1euZMyYMTRq1IhTp07l2k6IkuDQhWs88eV2Zq49xiN1K7Dhzfb0DZDKqw+TJAgLCw4OZsmSJRw4cIC//vqLoKAgAPr06UNQUBAHDx6kXr16fPvtt7Rq1YpevXrx0UcfceDAAWrWrJlrOyGKs5R0Ix+vO8bjs7dz6Voqc/v789WAJlRwti/o0EocSRC5mD9/Ph4eHsyfP/+B+/r333/p3bs3Dg4OlC5dOutucmFhYbRt2xZfX19++uknDh06lOv6eW0nRHEQfDaWxz7/l9mbT/JE46pseKMdj/pWLuiwSixJELmYMmUKERERTJliuevzBg8ezOzZswkNDWXSpEm3Leed13ZCFGWJqRlMXnmIp77aSUq6iUUvNOPjvg0p42Bb0KGVaJIgcjFx4kTc3d2ZOHHiA/fVrl07VqxYQXJyMvHx8fz5558AxMfHU7lyZdLT03OUsb65tPbt2glRXGw9HkWXWVv5fmc4g1p68vfr7WhXx62gwxLIWUy5Gjp0KEOHDs2Xvvz9/Xn66adp2LAhFSpUoGnTpgC8//77NG/eHDc3N5o3b56VFJ555hmGDh3K559/ztKlS2/bToii7lpSOlNXH+a34AhquDny20stCfCU4nqFiUXLfT9sAQEBeu/evTnmSTlncSfy91Ew1oZdZMIfh4hNTOOldjUYJcX1CkyhLvcthCg5rsSnMOmPQ6wJu0SDKqVZ+HxTGlSR4nqFlSQIIYTFaa1ZGhzB1NVHSE43MrZbXYa2leJ6hV2JSBBaa7m4RtyiOB1eLczOxybxzu+h/Hsimqaerkx/0o+abk4FHZbIg2KfIOzt7YmJiaFcuXKSJEQWrTUxMTHY28vFV5ZiMmkW7Qxn5rpjKOC9Xg0Y0KI6BimuV2QU+wTh7u5OREQEUVFRBR2KKGTs7e1xd3cv6DCKpZNX4hm3LJTgs1dpX8eND3r74O4qxfWKmmKfIGxsbPDy8iroMIQoEdKNJuZtPc1nG07gYGfFJ/9pSO/GVWXvvYgq9glCCPFwhEVeY8zSEI5cvM5jfpWZ3LMBbs52BR2WeACSIIQQDyQl3cinG04w/9/TlHW05avnmtDNp1JBhyXygSQIIcR923MmlvHLQjgdncjTAR68070eLg42BR2WyCeSIIQQ9ywhNYMZa47yw66zuLuW4scXm9OmdvmCDkvkM0kQQoh7svnYFf67PJSL11N4obUXb3Wtg4OtfJQURxa9jFEp9Z1S6opSKuw2y5VS6nOl1EmlVIhSyj/bskFKqROZ0yBLximEuLuriWm88csBnl8QhKOdNUuHt2Jiz/qSHIoxS/9mFwKzgUW3Wf4oUDtzag7MBZorpcoCk4AAQAPBSqmVWuurFo5XCHETrTV/hV5i0sow4pLSGfVILV55pBZ21lJcr7izaILQWm9VSnneocnjwCJtrnmwSylVRilVGQgE1mutYwGUUuuBbsBiS8Yriof58+czbtw4tNbMnDkTMN8EauLEiVll3OfPn3/LvHvpf+zYsaSmpqKUyirZYW9vz4wZMxg6dGhWG6VU1rzb9ZU9jjvF1aJFC3bv3p312sHBgT59+rBly5Zc2w8YMIDFixfTr18/fvjhh6y+AwMD2bJlS9bjzT+XcePGkZKSgr29Pe9+8BE/Hc0gppQ76ZdP8XITZ/b/sBjH7jn7HTduHPHx8RiNRvr3788PP/yQ4/0FBgayevXqrN/Jzdu7ef4NWmsyTBmkm9Kzpu9//J7PZ3/OiJEjeOLJJ1i6fCnzv53P8y8+T/ce3ckwZWDURowmIxk6A6PJaH6tjZi0ic1bNrPyz5U81vMxWrVuhUmb0GhM2pRjyj5Pa42J/z3X6Kx5aHK8Dg0NZffu3TRr3owGDRpkLdPoHO/rxn83XsP/+rnx/OZlN9Nac+LECQ7uOciYJmPy7TYFN1i83HdmgliltfbJZdkqYLrWelvm643AOMwJwl5rPTVz/gQgWWv9cS59DAOGAVSrVq3J2bNnLfNGRJHh4eFBREQEQNaV0hEREbi7u3P+/PkcbbLPu5/+b3ajv5tjuN02bo7jTnHldrGZlZUVRqMx1/bW1tYYjUasrKzIyMjI6vvGOrmtmz1uR9/OlOs0FJQVcdt+4nrQCtyrVuHixYu39Hsh+gJWDlYYHAzYONmwev1qEtMTGfnGSK4lX8PawRpswWBrwNnVmR69e5BiTGHj1o2k63SUjcK2lC3unu6kGlNJM6aRakwl3ZR+T78bS1MoDMqAQqGUyvUxMSERk8mEQRko7VI6ZxvMv8Mbv8sby7KeozD/f2u77I9Zy7QmIuI86ddSSfom7Z7/ljP7Kb7lvrXW84B5YL4fRAGHIwqBiRMnZn0rvXFXwBvfyrO3uXnevfR/uz2IG/3daKOUuuM2bo7jTnE1b978jnsQN+vXr1/WHkT2vnPbg8gez/ipH2Mf2A87T3cqOpzGKXY3B8r9S9VnKtOqUysiYiO4HH8Zl0outF3SFtf3XSmjyuTY9isbXzHH+IQDDphLbOg0jSnNRBnHMpyIO4G9lT3VPatz+vhpdJLGu6o39SrVw87KDjsrO2ysbMyPBhtsDDbYWtliY7Bhx7YdrFyxkiefeJIunbqwaf0mlvy8hOeefY6ej/XE2mCNlcEKa2V+NCgDVsoqa/r1l1/57LPPeH306wx4bgBWytxGKYUBQ9ZzK2V1y7y8mD9/PlOmZu4F9sv8Rm8yQXoipMZnm65DWiKkJkBagvl5jsckSE8yv05PNj9PT8p8njkZUwG4nKhYOfGjPMV3Lwp6D+JrYIvWenHm62OY9x4CgUCt9Uu5tbud3G4YJITIyWgyciXpChEJEUTER3A+/jyRCZFcTrzMydhI4tKiUYaMW9azMdjgau9KOftyuNq74mLrQmm70pS2NU/Ots6UtiuNk40TjjaOWZODjQMO1g5YG4r491GTCVLiIPnqrVNSrHlZyrXcp9R4yOUQ0S2UAWwcwdYRbB3Mz21K5XyefbLOfCxVBpoMvq+3VZj3IFYCI5VSSzAPUl/TWl9USq0DpimlXDPbdQHeLqgghSiKUo2pnLl2hpNxJzl59SSn4k4Rfj2cyITIHIdurJQV5e0rci3BgfiEiriX9uMJn/rUKleVig4VsxKCk41T8auplJ4MCVcgMQoSLpufJ1yBxCuQGA1JMZmP0eYkoI2378uuNNiXAXsX8+TqaX60Kw32pcHOOdtUGmydwM4pMxk4mx9tSkEh+hlbNEEopRZj3hsor5SKwHxmkg2A1vor4C+gO3ASSAKez1wWq5R6HwjK7GrKjQFrIcStktKTOBJ7hNCoUEKjQzl+9Tjn4s9h0iYArJU1ni6e1HatTYdqHfBw9sDdyZ1KDlVZuTeJOZvP4GhnxfSeDXi8UZXikQhSrsG1iMzpPMRfgusXIf5C5vML5m/9uSnlCo5u4FAOytWEas3Nzx3Kg0NZKFXW3CZrKgOG4ndWV7G/J7UQxdHFhIvsuriLkOgQQqNCORl3EmPmt9uqTlWpV7YeNcvUpJZrLWq51KJ66erYWOUsgRESEcfYpSEcvRRPD7/KTO7VgPJORai4Xmo8XA3POcWd+19SSL2es72yAqeK4FwJSlcB58rgXBGcKoFThcypojkJWNs+/PdTQArzISYhRB5cT7tO0MUgdl7cye6Luwm/Hg6As40zPuV9eNH3RfzK++FT3odypcrdsa/kNCOfbjjO/H9P4+Zsx/yBAXSuX/EhvIv7kJ4MMacg5mTOKfaM+bBPdvZloEw1KFsDvNqBi3vm5AGlq5oTQDH8lm9JkiCEKKTOx5/n7/C/2XR+E2HRYZi0iVLWpQioGMB/6v6H5pWbU6tMLQwq7wURdp2OYfyyEMJjkujXzIO3u9ejtH0hKK6XlghRxyDqKFw5DFeOml9fO5eznXMV8yGfej3Mx/izT6Vcb+1XPBBJEEIUIueun+Pvs3/zd/jfHIk9AoBPOR+G+Q2jReUW+JX3u+VQUV7Ep6Qzfc1Rftp9jmplHfh5SHNa1SqA4npam4/9XwqBiyHmx8thcPUsWWf5WNlC+brg0Qz8B5gTQrlaULameVBXPDSSIIQoYNdSr/HnqT/549QfHI09CoBfeT/eCniLztU7U8WpygP1v+noZf77exiXr6cwpI0Xb3apSynbh3So5fpFiNwLkcFwYT9cCjWfGQSAMn/4V24EjfqDmzdUqAeuXmAlH02FgfwWhCgAWmvCosP45dgvrA1fS6oxFd/yvowJGEPn6p2p7FT5gbcRm5jGlD8PseLABepWdGbuc01o5FHm7iver/RkiNwHEXvMCSEi2HzGEIDBBirWh7rdoXJDqORnfm3nbLl4xAOTBCHEQ5SUnsTqM6v57dhvHIk9QinrUvSq2Yv/1P0P3mW982UbWmv+DLnI5JWHiE9J57VOtRkRWAtb63wu3px8Fc7thnM7zdOF/WBMMy8rWwM8W0PVAKjaBCr5go19/m5fWJwkCCEegsT0RBYfXcz3h74nLjWOOq51mNBiAt29uuNkm3/H1S9dS+HdFWFsOHKZhh5lmPmkH3Ur5dO39NR4OGdXvlIAACAASURBVLsDTv8DZ/6By4cAbd47qNIYWrwM1VqZxw4cyubPNkWBkgQhhAUlpCXw89GfWXR4EddSr9GmahuG+g6lcYXG+XoxmtaaJUHnmbb6COkmE//tXo8X2nhhZXiAbRgzICIITm8xJ4SIIDBlgJUdVGsBHf4L1Vua9xBsSuXbexGFhyQIISwgPi2en478xA+Hf+B62nXaubdjuN9wfN18831bZ2MSGb8slJ2nY2hZoxzTn/SlejnH++ssMRpOboDj6+DURvPVyCjzHkKrUVCjPXg0l4RQQkiCECIfmbSJ5SeW89m+z4hLjSPQI5DhDYfToFyDfN+W0aRZsP0MH/99DBuDgQ/7+PJMU4972zPR2nzdwdHV5qQQGQxocKwA3j2gdhdzUpBrDEokSRBC5JPQqFCm7Z5GWEwY/hX8GddsHPXL1bfIto5dimfsshAOno+jU70KTH3Cl0oueRwE1hou7IMjf8LhlRB7yjy/ij8EjjcnhcqNwGDROxKLIkAShBAPKDYlls/3fc7yE8spV6ocH7b9kMe8HrNIwbu0DBNfbjnJnM0ncba34bNnGtGrYR6K62ltHkM49Ls5MVw7b65N5NUOWr4C3o+ZaxQJkY0kCCHuk0mb+O3Yb3y+/3OS0pMYWH8gwxsOz9ezkrI7cD6OcUtDOHY5nscbVWFSzwaUdbxLUbmoYxDyK4T+BnFnzQPMNR+BwLeh7qNytpG4I0kQQtyHy4mXeXf7u+y6uItmlZrxTvN3qFmmpkW2lZxm5P/+PsZ3289QwdmebwcF0LHeHYrrXb8AoUsh9FfzlcvKADUCzYePvHuY700gRB5IghDiHq0/u573dr5HmjGNiS0n8lTtpyx2/4Qdp6IZvyyUc7FJPNu8GuMf9c69uJ4xHY6vhX0/wMn1oE3mMYVu06FBH3NZayHukSQIIfIoMT2RD3d/yB+n/qBBuQZMbzsdTxdPi2zrWnI609ccYfGe81Qv58DioS1oWTOXMt7RJ2DfIji42HxXNOfK0OZ1aPgslK9lkdhEySEJQog8OHDlAG//+zYXEi8w1HcoLzd6GRuDZcpkrz98mXdXhBIVn8pL7WrwWqc6OYvrZaTBkZUQ9I25xIXBGup0A/+BULOjFLoT+Ub+koS4A601iw4v4pPgT6jsWJkFXRfgX9HfItuKTkhl8spDrAq5iHclZ+YPDMDPPVtxvesXIXghBC8w3z/Z1Qs6vQcN+8khJGERkiCEuI1UYypTdk5h5amVdK7emfdavYezbf5XH9Va88eBC7z35yESU4282bkOL7WvaS6upzWc2wV75pn3GkxGqN0Zmg0z7y3ItQrCgiRBCJGLqKQoXtv8GiHRIYxoNIKX/F66pzu35dWFuGT++3som49F0biaubhe7YrO5jpIYStgxxfmi9rsXaD5cGj6orlSqhAPgSQIIW4SFh3G6E2jiU+PZ1bgLDpV75Tv2zCZND/tOceMNUcxmjQTe9RnUCtPrDKSYPc82DnbfN1C2Zrw2Cfmw0i2DvkehxB3YtEEoZTqBnwGWAHfaK2n37R8FtAh86UDUEFrXSZzmREIzVx2Tmvdy5KxCgGw6vQqJm2fhJuDGz90+oG6Zevm+zbORCcyblkIe87E0qZWeT7s44uHbSJs+cA88Jx81VwQr+s088Vshod09zchbmKxBKGUsgLmAJ2BCCBIKbVSa334Rhut9evZ2r8KNM7WRbLWupGl4hMiO601cw/OZe7BuQRUDOCTwE9wtc/fAnUZRhPfbjvDJ+uPY2ttYOaTfvSta43aMRX2fgcZKeaSF61GQbXm+bptIe6HJfcgmgEntdanAZRSS4DHgcO3ad8PmGTBeITIlUmbmLFnBj8f/ZnetXozoeWEfD+F9fCF64xbFkJo5DU616/ItI6uuB34DNYuMt9jwe8/0PZNKF87X7crxIO4a4LI/Gb/o9b66j32XRU4n+11BJDr1yKlVHXAC9iUbba9UmovkAFM11qvuM26w4BhANWqVbvHEEVJl25KZ+L2iaw6vYrBDQbzRpM38vWq6NQMI7M3nWTullO4lLLhu8cr0CHqB9S3PwHaPLbQ9g0ZeBaFUl72ICpiPjy0D/gOWKe11vkcxzPAUq21Mdu86lrrSKVUDWCTUipUa33q5hW11vOAeQABAQH5HZcoxlIyUnjrn7f4J+IfRjUexRDfIfmaHILPXmXcshBOXklgkI89bzuvwn79D6CU+aK2Nq9BGflSIwqvuyYIrfW7SqkJQBfgeWC2UupX4NvcPrCziQQ8sr12z5yXm2eAV27abmTm42ml1BbM4xN32p4QeZaQlsCrm14l+HIw7zZ/l6e9n863vpPSMvh43XEW7DhDXed0/mm8jeonfgRTOjQeAO3GgEvVfNueEJaSpzEIrbVWSl0CLmE+5OMKLFVKrddaj73NakFAbaWUF+bE8Azw7M2NlFLemf3tzDbPFUjSWqcqpcoDrYGZeX9bQtxebEosw9cP58TVE0xvO53uNbrnW9/bTkQzfnkIcVdj+MZzJ49c/RV1JME8xhA4Xg4liSIlL2MQo4GBQDTwDTBGa52ulDIAJ4BcE4TWOkMpNRJYh/k01++01oeUUlOAvVrrlZlNnwGW3HTYqh7wtVLKBBgwj0HcbnBbiDyLTYnlxXUvcj7+PJ898hnt3NvlS7/XktP5YPVhlu8N59XS2xjhshSbS7Hm8tod/gsVLXNnOSEsKS97EGWBPlrrs9lnaq1NSqked1pRa/0X8NdN8ybe9HpyLuvtAPL/7u6iRLuWeo2X1r/E+fjzzOk4h+aV8+dU0nWHLjHh91CaJO9gd5mllEs5B55tofN7ULVJvmxDiIKQlzGI2556qrU+kr/hCGEZCWkJDF8/nFNxp/jikS/yJTlExZuL60WG/ct3jkvwsTkMTnWh9y9Qp6t5MFqIIkxKbYhiLyk9iREbR3A09iizOsyiddXWD9Sf1prf90cy/89/GGn8gcfsdqLtKkC3WdB4oJTbFsWG/CWLYi0lI4VXN73KwaiDfNTuIwI9Ah+ov8i4ZCYv3YNP+AJWWK/GxtYArcehWr0Kdvlf6VWIgiQJQhRbacY0Xtv8GkGXgpjWdhpdPLvcd18mk+bHXeGErv2OKeonKlvHYPJ5CkPn98DFPR+jFqLwkAQhiqUMUwZv/fMW2y9sZ0qrKfSoccfzKe7oVFQC85Ys46mo2Qw0HCfNzRd6/IShest8jFiIwkcShCh2tNa8v+t9Np/fzDvN36F37d731U+60cSiTftw+PdDPjRsJM2+LLrr59g2fk4qrIoSQRKEKHbmHJjD8hPLecnvJfp597uvPsIirrLx5/9jQOICXAxJpDQZhkPn/5pv3CNECSEJQhQrvx77la9DvqZP7T680uiVu69wk5R0I7+s/BO/g+8z2nCSWLcArPp+jkPFBhaIVojCTRKEKDY2ntvIB7s/oL17eya0mHDPhff2Hwvn7NK3eS5tDYm2riR2/ZKyAc/K9QyixJIEIYqFfZf3MW7rOHzK+TCz3UysDXn/005MSWfNktm0PzMLPxXPpboDqdp7CpQqY8GIhSj8JEGIIu9U3ClGbhpJZcfKzO44GwebvN+7eXfwXtTqN3jKdJBIx3qkPf0lVav7WzBaIYoOSRCiSLuSdIXhG4ZjZ2XH3E5z83yb0Lj4BHYsmsQjV77HqKw523wy1buOkrOThMhGEoQospLSkxi5cSTXU6+zsNtC3J3zdsHars1/UuGf8XQngqPlO+LZ/3Oql5WL3YS4mSQIUSQZTUbG/TuOY1eP8cUjX1CvXL27rhN15TJHfnyDdtdXcdlQgbNdFuLd4v6ukRCiJJAEIYqkT4I/Ycv5LYxvNv6u93TQWrNz1QJqBU+htY7jgMcAGvT/EJtSUjtJiDuRBCGKnF+P/cqiw4t41vtZ+tfrf8e2F86d4sLikbRK3sEZ6xqk9f6JRg0erJqrECWFJAhRpGyP3M603dNo596OsU1vd7dbMBqN7P7t//A98gk+ZBBc9zUa9/0vBhvbhxitEEWbJAhRZJy4eoI3/3mTWmVqMbPdTKxuc8bR2eMHSfxtBK3Swzhk34hy/ebSxFNu+SnEvZIEIYqE6ORoXtn4Cg7WDszuOBtHG8db2qSnpRL08xSanPmaVGVLcKOp+Pd6BWUwFEDEQhR9kiBEoZdqTGX0ptHEpcaxoNsCKjlWuqXNiQPbMPz5Kq2Mp9nv3I5q/efQpHK1AohWiOLDol+tlFLdlFLHlFInlVLjc1k+WCkVpZQ6kDkNybZskFLqROY0yJJxisJLa82E7RMIiQ5hWptpNCiXs2heSlICO79+Fa/fe+JijGV/yy9o/NaflJPkIMQDs9gehFLKCpgDdAYigCCl1Eqt9eGbmv6itR5507plgUlAAKCB4Mx1r1oqXlE4fR3yNWvOrGG0/2g6Ve+UY9nhnWtw/vsNWuoLBLk+Rp1Bn9HY1a2AIhWi+LHkIaZmwEmt9WkApdQS4HHg5gSRm67Aeq11bOa664FuwGILxSoKobXha5lzYA69avbiRZ8Xs+bHX4vl8A9v0jx6ORdURcI6LqJp28cLMFIhiidLHmKqCpzP9joic97NnlRKhSilliqlPO5xXZRSw5RSe5VSe6OiovIjblEIhEWH8e62d2lcoTGTWk7KKt19cPNvJM5qStOo39lV4WnKvBmEjyQHISyioE/v+BPw1Fr7AeuB7++1A631PK11gNY6wM1NDi8UB5cSLzFq0yjKlyrPpx0+xdbKlqvRlwia1ZeG/wwhxVCK4z2X0WLEPByc5A5vQliKJRNEJOCR7bV75rwsWusYrXVq5stvgCZ5XVcUT0npSYzaNIqkjCRmPzIbV9syBP/1DXp2UxrFbWSX+4tUHrsH74COBR2qEMWeJccggoDaSikvzB/uzwDPZm+glKqstb6Y+bIXcCTz+TpgmlLqRu3mLsDbFoxVFAImbeLtf9/m2NVjzH5kNi5JVhyc/xhNknZwwqo2cb3n0MKneUGHKUSJYbEEobXOUEqNxPxhbwV8p7U+pJSaAuzVWq8ERimlegEZQCwwOHPdWKXU+5iTDMCUGwPWovj6NPhTNp3fxNiAMdjv3oN9aF/qaiM7a71O02fewVrKZAjxUCmtdUHHkG8CAgL03r17CzoMcR+WHl/Kezvfo2flLjy7dzc+aaGE2TbE5em5eNRscPcOhBD3RSkVrLUOyG2ZXEktCtzOCzv5YNdUfKjAO9u/x4QNO30n07z3aAxWBX0ehRAllyQIUaBOx53m9U2jqZpqZP6FYI46tsG9/2xaVvUs6NCEKPEkQYgCc/nqWYb88R/sjUl8eDmZQ00/o9mjA7OueRBCFCxJEKJAHN32K5MPvcd1WxidVI+qL39F2XJyHYsQhYkkCPFQJcdGcurHkSyyOsghJ0dGuPVnwPO31HEUQhQCkiDEw2EycXrtF1TYM52/XO1Z4+TEyz4jeLnJywUdmRDiNiRBCItLOLuf2F9GUCPpMNNc6rK4TDLP1H2Gl/2HF3RoQog7kHMIheWkJnB2yRvYL3gEh8QIpnoNYknZFDpW68j4ZuNlMFqIQk72IET+05r4fUsxrhlP9Yxo1th0Ib7bQJYfnkSj8o2Y3nb6be8nLYQoPCRBiHylo44T9etoKkTt4LCuzrqGM/Br58eL6wfj7uzOF498gb21fUGHKYTIA0kQIn+kJZGwYTr2e+Zgr2352vllAp8dR3uXVJ776znsrez5qtNXuNhJeW4higpJEOLBaI3p0B8krX4bp+QLrDC1JbHdJIY8EkBC+nUGr32ZhPQEFnZbSBWnKgUdrRDiHkiCEPfvUhjJf46hVOQOzps8+NXtYwY/+yzVyzmSkJbA8PXDOXv9LF92+hLvst4FHa0Q4h5JghD3LjEG06apELyQVO3AxwyhTvdXmNjcC6UUSelJvLLxFY7GHmVWh1m0qNyioCMWQtwHSRAi74wZsPdbjBs/gLR4FmV05kCNl3j7ydZUcjEPPKcaUxm9eTQHog4wo90MAj0CCzZmIcR9kwQh7k5rOLYGvX4iKuYEu0w+fGr9IgP7duNTv8pZ1zOkm9J5a8tb7Lq4i6mtp9LNs1sBBy6EeBCSIMSdRQbD3xPg7HYiDFV5L+1NnP168nXPBpR1/N8d3owmI2//+zZbIrbwbvN3ebzW4wUYtBAiP0iCELm7ehY2vQ+hv5Fo7cr09OfZ7PgoUwY24hHvijmamrSJSTsmsS58HW8FvMXT3k8XUNBCiPwkCULklBgD2z6BPfMwYuAn677MTOjKEy28WdPNG2d7mxzNM0wZTNg+gVWnVzGi0QgGNRhUQIELIfKbJAhhlnINds6BnXPQ6UkEl+nKyIvdsS/nwTfP+tGiRrlbVkkzpjHmnzFsOr+JUY1HMcR3SAEELoSwFEkQJV1aEgTNh22zIPkql9wfZfTlbgRdcmNo+xq83qkO9ja31k1KSk/itc2vsfPiTsY3G0//ev0LIHghhCVZNEEopboBnwFWwDda6+k3LX8DGAJkAFHAC1rrs5nLjEBoZtNzWuteloy1xElPgf0/wNaPIeESaV4d+cT4H7467ox3JWdWDPLDz71MrqteT7vOKxteISQ6hKmtp8qAtBDFlMUShFLKCpgDdAYigCCl1Eqt9eFszfYDAVrrJKXUy8BM4MYIZ7LWupGl4iux0pJg3/ew/TOIv4iu1pJtjWYyakcpElONvNm5FsMDa2JjlXsl+JjkGIZvGM7JuJN83P5jOlfv/JDfgBDiYbHkHkQz4KTW+jSAUmoJ8DiQlSC01puztd8FPGfBeEq21ATY+x3s+BwSo8CzLdFdvmBMUGk2b4imcTVHZj7pR+2Kzrft4lLiJYatH8bFhIt88cgXtKna5iG+ASHEw2bJBFEVOJ/tdQTQ/A7tXwTWZHttr5Tai/nw03St9YrcVlJKDQOGAVSrVu2BAi6Wkq9C0LfmAejkWKjRAVPbMfx02Z0ZS49iNF1lYo/6DGrliZXh9jfwCYkKYdSmUaQaU5nbaS4BlQIe4psQQhSEQjFIrZR6DggA2mebXV1rHamUqgFsUkqFaq1P3byu1noeMA8gICBAP5SAi4K4c7BrLuxbBGkJUKsztB/HmVL1GbcshD1nwmhTqzwf9vHFo6zDHbtafXo1E7dPxM3BjW+6fEMt11oP6U0IIQqSJRNEJOCR7bV75rwclFKdgP8C7bXWqTfma60jMx9PK6W2AI2BWxKEuMnFEPNhpLDloBT4PAmtXiXDrQHfbDvDrPVbsbM2MPMpP/o2cb/jbT9N2sTs/bOZHzqfJhWbMCtwFq72rg/xzQghCpIlE0QQUFsp5YU5MTwDPJu9gVKqMfA10E1rfSXbfFcgSWudqpQqD7TGPIAtcmMywon1sHsunN4Ctk7Q4mXz5OLOkYvXGfvlDkIjr9GlfkWmPuFDhdJ3vqtbUnoS72x7h43nNtKndh/ebf4uNlY2d1xHCFG8WCxBaK0zlFIjgXWYT3P9Tmt9SCk1BdirtV4JfAQ4Ab9lfpO9cTprPeBrpZQJMGAegzic64ZKsqRYOPATBH0DV8PBuTJ0eg+aDIZSZUjNMDLn72N8ueUUZRxs+LK/P4/6VLrjXgOYB6Nf3fQqx68eZ2zTsTxX77m7riOEKH6U1sXnsH1AQIDeu3dvQYdheZdCYc88CPkNMpKhemtoNhS8e0Dmt/zgs1cZtyyEk1cS6ONflQmP1cc1W3G929l4diOTdk4iw5TBR+0+oq17W0u/GyFEAVJKBWutcz3rpFAMUos8SI2HQ7+bB50jgsC6FPj9B5oNg0o+Wc2S0jL4aN0xFu4Ip4pLKRY+35TAuhXu2n1yRjIfBX3Eb8d/o365+sxoOwNPF08LviEhRGEnCaIw09qcDPZ9D2G/Q3oilK8LXT6Axv2hVM4B420nohm/PISIq8kMbFmdsd28cbK7+6/4aOxRxm4dy5lrZ3je53lebfSqjDcIISRBFErXIiFsKRz4GaKOgo0j+PQB/4Hg3tR8dlL25knpfPDXYX7dG0GN8o78+lJLmnmVvetmTNrEj4d/5NN9n1LGrgzzOs+jZZWWlnpXQogiRhJEYZEcB0dWQsivEL4N0OZk0PNzc3Kwy/0K57Vhl5jwRxixiWmMCKzJqI61cy2ud7Nz18/x/q732XVxF4EegUxpNUVOYRVC5CAJoiClJcLJDRD6GxxfB8Y0KFsTAseDb18oV/O2q0bFpzJ55SFWh16kfuXSLBjcFJ+qLnfdZKoxlW9Dv+Xb0G+xsbJhQosJ9K3TV85SEkLcQhLEw5Zy3ZwMjvwBJzaYz0JydIOAF8GvL1Txv+UQUnZaa5bvi2TKqsMkpxsZ07Uuw9rVuG1xvey2RW5j2u5pnI8/z6NejzImYAxuDm75+e6EEMWIJIiHIeEKnPgbDq+E05vNewpOlcwDzfV6mU9Ttbr7ryLiahLv/B7G1uNRNKnuyown/ahVwemu611KvMTMoJmsP7sez9KefNPlG5pXvlNZLCGEkARhGSYTXDxgTgrH18GFfeb5LtWg6VCo3wvcm4Hh7t/6zd1pftx9lhlrjqKByT3rM7ClJ4Y7FNcDiE2JZWHYQpYcW4LWmtH+oxlYfyC2Vne/HkIIISRB5JfEaDjzD5zcBCfXQ8JlQIF7AHR4F+p0gUp+dzx8lJtTUQmMXxZCUPhV2tVxY1pvH9xd71xcLy4lju8Pf89PR34i1ZhKd6/ujGw8kqpOVR/gDQohShpJEPcrLRHO7jQfMjr9D1zOvPmdnQvU6gh1ukKtTuBY/r66TzeamP/vaT7dcIJSNlZ83LchT/pXveNg8rXUa/xw+Ad+PPIjSelJdPPqxvCGw6nhUuO+YhBClGySIPIqOQ7O74FzO8yJITIYTOlgZQsezeGRCVCjA1RumKfxhDsJi7zGuGUhHLpwne6+lZjcqwEVnG9fXO/s9bP8cuwXfj/xOwnpCXSp3oWXG74sZbmFEA9EEkRutDYXv4sMhvO7zQnhchigwWANVRqbK6XWCIRqLcH2zod88iol3cgXm07w1T+ncXWw5avn/OnmUznXtkaTka0RW1lybAk7LuzAWlnTqXonhvgOoW7ZuvkSjxCiZJMEAeaqqBf2QUSwOSlEBkNStHmZjSN4NIXAt6F6S6gakG8JIbu94bGMXRbC6ahE+jZx593H6uPicGu5i6ikKP449Qe/HfuNC4kXqFCqAiMajeCp2k/JKatCiHwlCSIjFf6vrvnUUxS4eUOdbuDexJwMKtR/4ENGd5KQmsFHa4+yaNdZqriUYtELzWhXJ+cHfUxyDBvObmBt+FqCLwej0TSt1JQ3A96kQ7UO2BikbpIQIv9JgrC2g8c+gTLVzIeO7Es/tE3/czyKd5aHcuFaMoNaejKma10cM4vrXU25yubzm1l7Zi17Lu3BqI14uXgxvOFwunl1k4FnIYTFSYIA8B/wUDcXl5TG+6uOsGxfBDXdHFk6vCV+7s4ciNrPjgs72HFhB0dijqDReDh78ILPC3T17Eod1zpSEkMI8dBIgnjI1oReZMIfh4hLTuLp1oq6nuF8d/I3grYFkZyRjLWyxs/Nj1cavUIb9zbUL1tfkoIQokBIgngItNaEXjrDe+vXEhYdglPlSJxtIvkrNp2/YsHD2YNeNXvRukprmlZqipPt3ctnCCGEpUmCyGdxKXGciDvB8avHORl3khNXT3Ak5jhppmQAHMrb0sDNh0ZuHWjo1hBfN18qONz9jm9CCPGwSYK4RxmmDKKTo7mQcIFz8ec4H3+e89fPcz7+POfiz3E97XpWWyeb0ui0SiRcbUQ1pxqM7dCZwBp+ctaREKJIkASRKSk9idiUWGJTYrmacpXYlFhiUmKITo7mcuJlLidd5nLiZaJTojFpU9Z6BmWgsmNlqjlXo5tnN6qVrkYNl5ocOGXP3I1RKBTvdq9H/2bV7lpcTwghChOLJgilVDfgM8AK+EZrPf2m5XbAIqAJEAM8rbUOz1z2NvAiYARGaa3XWSJGrTWtF7cmPj0+1+WONo5UcqhERceK1CpTi4qOFanoUJFKjpWo5lyNyk6Vc+wRnLwSz7hloQSfvUxgXTc+6O1L1TKlLBG6EEJYlMUShFLKCpgDdAYigCCl1Eqt9eFszV4ErmqtaymlngFmAE8rpeoDzwANgCrABqVUHa210QJx0r9+f+yt7ClrX5ZypcpR1r4sZe3L4mrvSinrvH24pxtNfP3PKT7feBIHOytmPd2QJxrdubieEEIUZnm7IcH9aQac1Fqf1lqnAUuAx29q8zjwfebzpUBHZf5EfRxYorVO1VqfAU5m9mcRtkG2TO4+megt0Rxbc4y2tdtSu2JtqrhVYf78+Vnt5s+fj4eHR9a8AQMGYDAYcKjqjeewuXz893ESj+8g7uc36ePvgcFgQCllbuPggJ2dHQaDAScnJ5RSWFtbU758eZRSWZOjoyPz589n/vz5lC1bFgcHBxwcHHB0dMx6tLGxwWAwMGDAgKw4rK2tadGiBWXLlsXV1TVH3KJkuvnvVYh7pbTWlulYqaeAbv/f3r3H2FHWYRz/PuyxsCvQ3V6AwrbdIlVQglBXtlCigIraP0BuhjXRolxEaOItJlXMSjAEUEKiQTCEQlAQUa4VuVhuMVFaKAj0BtKC7W5TCrZS7tAtP/84Q3PYztk97dmZOe0+n+RkZ8682Xl4z/T8dt55mYmIM5P1rwNdETG7os2SpE1fsr4S6AIuABZExA3J+3OBeyLilpT9nA2cDTBp0qRPrVq1apuzTpw4kb6+Ptrb2wHo6+vbsq29vZ3e3t6t2vX29lLarYU9pn+VPQ8/ic1vvMKG+Vfx1nMLtnn/A6XlSNPU1ER/fz+lUonNmz94clWZ20amgcerWRpJj0dEZ9q2LM8gchERV0dEZ0R0jh+/fTer6+npob29nZ6eHnp6emhtbaWlpYXW1lZ6enpS2y18fj1TZ89l9PRTeXPJd++HPwAAB1xJREFUA/zvxh9A31O0tbXR0dHxgd8viebmZkaNGrXlLAHKX/Bjx479QNuWlpYtOdra2mhubqa5uZmWlpYtP0ulEpLo7u4GoLu7m6amJrq6umhra9sqt41Mlcer2fbI8gziCOCCiPhisv5jgIi4uKLNfUmbRySVgBeB8cCcyraV7QbbZ2dnZyxatCiL/5wtXnt7E5fe+ww3LFjNxDHNXHLSIcw4YPseCmRmVrTBziCynMX0GDBV0hRgDeWLzl8b0GYeMAt4BDgFeDAiQtI84A+SLqd8kXoq8GiGWWvy0DMvcf7ti1n76tuccdQUfnjcR2kZ5ZnCZrZzyuzbLSL6Jc0G7qM8zfXaiFgq6UJgUUTMA+YCv5e0AthAuYiQtPsTsAzoB87LYgZTrTa88S4/v2sZt/9rDVP32p1bv3Mk0ya1FRXHzCwXmQ0xFWG4h5gigr8uXsvP7lzKxrc2ce4xB3DeMR9h11LTsO3DzKxIRQ0x7dDWvfo2P71jCfOXreOQ9tHccGYXB03I71kRZmZFc4EYICK4+bFeLrp7Oe/2v8dPZh7It2ZModS0w0/4MjPbJi4QFVavf5M5tz3NP1eup2vKGC49+RA6xn246FhmZoVwgQA2vxdc948XuOxvz1LaZRcuOvFguj/tm+uZ2cg24gvExjc3Meu6R3my9xWOPXAvLjrxYCaM9s31zMxGfIHYs7nE5LEtfHNGB8d/cl/fXM/MLDHiC4QkfnXaYUXHMDNrOJ6aY2ZmqVwgzMwslQuEmZmlcoEwM7NULhBmZpbKBcLMzFK5QJiZWSoXCDMzS7VTPQ9C0svAqqJzDGIc8N+iQwzC+erX6BkbPR80fsZGzwfblnFyRIxP27BTFYhGJ2lRtQdzNALnq1+jZ2z0fND4GRs9HwxfRg8xmZlZKhcIMzNL5QKRr6uLDjAE56tfo2ds9HzQ+BkbPR8MU0ZfgzAzs1Q+gzAzs1QuEGZmlsoFIkOSTpW0VNJ7kqpOOZP0JUnPSlohaU6O+cZImi/pueRnW5V2myU9mbzm5ZBr0P6QtKukm5PtCyV1ZJ1pOzKeLunlin47M+d810p6SdKSKtsl6ddJ/qclTWuwfEdL2ljRfz0555so6SFJy5J/w99NaVNYH9aYr/4+jAi/MnoBBwEfAx4GOqu0aQJWAvsDo4CngI/nlO8XwJxkeQ5waZV2r+fYZ0P2B3Au8Ntk+TTg5pw/11oyng5cUeCx9xlgGrCkyvaZwD2AgOnAwgbLdzRwV4H9NwGYlizvAfw75TMurA9rzFd3H/oMIkMRsTwinh2i2eHAioh4PiLeBf4InJB9Okj2c32yfD3wlZz2O5ha+qMy9y3A55Tvw8SL/MxqEhF/BzYM0uQE4HdRtgBolTQhn3Q15StURKyNiCeS5deA5cB+A5oV1oc15qubC0Tx9gN6K9b7yOCDrmLviFibLL8I7F2l3W6SFklaICnrIlJLf2xpExH9wEZgbMa5UvefqPaZnZwMPdwiaWI+0WpW5HFXqyMkPSXpHkmfKCpEMoR5GLBwwKaG6MNB8kGdfViqM9uIJ+l+YJ+UTedHxJ155xlosHyVKxERkqrNeZ4cEWsk7Q88KGlxRKwc7qw7mb8AN0XEO5K+TfmM59iCM+1InqB83L0uaSZwBzA17xCSdgduBb4XEa/mvf+hDJGv7j50gahTRHy+zl+xBqj867I9eW9YDJZP0jpJEyJibXJq/FKV37Em+fm8pIcp/7WSVYGopT/eb9MnqQSMBtZnlCfNkBkjojLPNZSv9zSSTI+7elV+2UXE3ZKulDQuInK7SZ6kD1H+8r0xIm5LaVJoHw6Vbzj60ENMxXsMmCppiqRRlC+6Zj5TKDEPmJUszwK2OuOR1CZp12R5HDADWJZhplr6ozL3KcCDkVyVy8mQGQeMRR9PeYy4kcwDvpHMxJkObKwYbiycpH3ev64k6XDK31W5/RGQ7HsusDwiLq/SrLA+rCXfsPRhXlfdR+ILOJHyuOQ7wDrgvuT9fYG7K9rNpDwLYSXloam88o0FHgCeA+4HxiTvdwLXJMtHAospz9RZDJyRQ66t+gO4EDg+Wd4N+DOwAngU2L+Az3aojBcDS5N+ewg4MOd8NwFrgU3JMXgGcA5wTrJdwG+S/IupMsuuwHyzK/pvAXBkzvmOAgJ4Gngyec1slD6sMV/dfehbbZiZWSoPMZmZWSoXCDMzS+UCYWZmqVwgzMwslQuEmZmlcoEwy0hyx80XJI1J1tuS9Y5ik5nVxgXCLCMR0QtcBVySvHUJcHVE/KewUGbbwP8fhFmGktshPA5cC5wFHBoRm4pNZVYb34vJLEMRsUnSj4B7geNcHGxH4iEms+x9mfJtJQ4uOojZtnCBMMuQpEOBL1B+4tj383woj1m9XCDMMpLcSfMqyvfqXw38Eris2FRmtXOBMMvOWcDqiJifrF8JHCTpswVmMquZZzGZmVkqn0GYmVkqFwgzM0vlAmFmZqlcIMzMLJULhJmZpXKBMDOzVC4QZmaW6v9DZwukyaXuggAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "vGolOtBPkztp",
        "outputId": "33a185db-0e74-4402-bd61-210f19bb07df"
      },
      "source": [
        "l1 = LogisticRegression().fit(X, y)\r\n",
        "print(l1.n_iter_, 'iterations with default initialization')\r\n",
        "l2 = LogisticRegression(warm_start=True)\r\n",
        "ym = y.mean()\r\n",
        "l2.coef_ = linreg.coef_[np.newaxis, :] / (ym * (1-ym))\r\n",
        "l2.intercept_ = logit(ym)[np.newaxis] - np.dot(X.mean(axis=0), l2.coef_)\r\n",
        "l2.classes_ = np.array([0, 1])\r\n",
        "l2.fit(X, y)\r\n",
        "print(l2.n_iter_, 'iterations with initialization from linear model')"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "[9] iterations with default initialization\n",
            "[6] iterations with initialization from linear model\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TQfWnaqSoTsA"
      },
      "source": [
        "On a larger dataset, it turns out that smart initialization helped to make training almost 3 times faster. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vOzDc1_CnTUF"
      },
      "source": [
        "X, y = datasets.make_classification(n_samples=10_000, n_features=1000, n_informative=200, n_redundant=400, random_state=42, n_clusters_per_class=10)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "PBV7L29EnbT8",
        "outputId": "78cc64ec-019b-409c-80ae-81fc6703eec4"
      },
      "source": [
        "l1 = LogisticRegression(max_iter=1000).fit(X, y)\r\n",
        "print(l1.n_iter_, 'iterations with default initialization')\r\n",
        "\r\n",
        "linreg = Ridge(alpha=0.01).fit(X, y)\r\n",
        "l2 = LogisticRegression(warm_start=True, max_iter=1000)\r\n",
        "ym = y.mean()\r\n",
        "l2.coef_ = linreg.coef_[np.newaxis, :] / (ym * (1-ym))\r\n",
        "l2.intercept_ = logit(ym)[np.newaxis] - np.dot(X.mean(axis=0), l2.coef_[0])\r\n",
        "l2.classes_ = np.array([0, 1])\r\n",
        "l2.fit(X, y)\r\n",
        "print(l2.n_iter_, 'iterations with initialization from linear model')"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "[292] iterations with default initialization\n",
            "[107] iterations with initialization from linear model\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 279
        },
        "id": "-CR493ReoR4G",
        "outputId": "f46fe834-1815-4941-99fc-be22ab37d802"
      },
      "source": [
        "plt.scatter(linreg.coef_, l2.coef_);\r\n",
        "plt.xlabel('linear regression coefficients'); plt.ylabel('logistic model coefficients');"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEGCAYAAABLgMOSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de5gcdZ3v8fcnw4ATlEwCAWEggoAgHBR1VlB0F+SOAllAwEUNinLco+uqaw5BORAuK8F4Wx9c16C4EVBu4jAIaww3LwhCwgAhSCQitwblkoSLCRCS7/mjqpNO091T09PV3TPzeT1PP1NVXV31q8xkvvO7fX+KCMzMzIZqXKsLYGZmI5MDiJmZ1cUBxMzM6uIAYmZmdXEAMTOzumzU6gI00xZbbBHbb799q4thZjaiLFy48OmImFx+fEwFkO23354FCxa0uhhmZiOKpIcrHXcTlpmZ1aWlAUTSIZKWSFoqaUaF9zeRdFn6/u8lbV/y3lsk3SppsaRFkl7TzLKbmY11LQsgkjqA7wCHArsBH5K0W9lpJwHLI2In4JvAeelnNwIuBj4VEbsD+wKrm1R0MzOjtTWQdwJLI+LBiHgZuBQ4suycI4G56faVwP6SBBwE3BMRdwNExDMRsaZJ5TYzM1obQHqAR0v2H0uPVTwnIl4BngU2B94EhKR5ku6U9H+bUF4zMysxUkdhbQS8B/g7YCVwg6SFEXFD+YmSTgZOBpgyZUpTC2lm1kp9AwVmz1vC4ytWsU13F9MP3oWpbyv/O71+rayBFIDtSva3TY9VPCft95gAPENSW/l1RDwdESuB64C3V7pJRMyJiN6I6J08+VXDmM3MRqW+gQKnXrWIwopVBFBYsYpTr1pE30D5r9n6tTKA3AHsLGkHSRsDxwP9Zef0A9PS7WOAGyPJPz8P2EPS+DSw/ANwX5PKbWbW9mbPW8Kq1Rt2Da9avYbZ85Y07B4ta8KKiFckfYYkGHQAF0bEYklnAQsioh/4AXCRpKXAMpIgQ0Qsl/QNkiAUwHURcW1LHsTMrA09vmLVkI7Xo6V9IBFxHUnzU+mx00u2XwQ+WOWzF5MM5TUzszLbdHdRqBAstunuatg9PBPdzGwUmn7wLnR1dmxwrKuzg+kH79Kwe4zUUVhmZlZDcbRVnqOwHEDMzEapqW/raWjAKOcAYmY2AuQ9p6MeDiBmZm2uOKejOCy3OKcDaGkQcSe6mVmba8acjnq4BmJm1kYqNVU1Y05HPRxAzMzaRLWmqu7xnSxf+eoVKxo5p6MebsIyM2sTM/sXV2yqiiD3OR31cAAxM2sDfQMFVqyqvC7es6tWc+5Re9DT3YWAnu4uzj1qD4/CMjMzanaIb9Pdlfucjnq4BmJm1gYq5a0qanVTVTUOIGZmLdY3UEBV3ps4vrPtah5FDiBmZi02e94SosJxAWccvnuzi5OZ+0DMzJqgViqSavM5gtbONB+MayBmZjkbbHnZavM5elo8z2MwDiBmZjkbLBVJM9buyIObsMzMctQ3UKg6wqrYdNWMtTvy4ABiZpaT0/oWccltj1R9v7Tpqh3neQzGAcTMrA6Drc/RN1DgktseqTi6CkZGE9VgHEDMzIaob6DA9CvuZvXaJDwUVqxi+hV3Axs2R1ULHkBbpCIZLneim5kN0cz+xeuCR9HqtcHM/sXr9mvNLO9JU5OMdA4gZmZDVC3pYenxDlWbW96+qUmGygHEzGwIinM3BrMmqjdgjYbaB7S4D0TSIcB/AB3A9yNiVtn7mwA/At4BPAMcFxEPlbw/BbgPmBkRX2tWuc1s7CjvLF/+t5eqnjtxfOe67Z7urorNWO0+OXAoWlYDkdQBfAc4FNgN+JCk3cpOOwlYHhE7Ad8Ezit7/xvA/+RdVjMbmyrNIF+5em3V89//lq3XbY/UyYFD0comrHcCSyPiwYh4GbgUOLLsnCOBuen2lcD+UtKwKGkq8GdgMWZmOag0g7yWm+5/at321Lf1tOUiUI3UyiasHuDRkv3HgL2qnRMRr0h6Fthc0ovAKcCBwBdr3UTSycDJAFOmTGlMyc1sTKiW5DDr+SNxcuBQjNRO9JnANyPihcFOjIg5EdEbEb2TJ0/Ov2RmNuL1DRTYZ9aNNedxVFItKeJo1coaSAHYrmR/2/RYpXMek7QRMIGkM30v4BhJXwW6gbWSXoyI8/MvtpmNVn0DBWb2L646TLeW0da/kUUrA8gdwM6SdiAJFMcD/1R2Tj8wDbgVOAa4MSICeG/xBEkzgRccPMxsOIod5kPp8+iQWBNBzwhJfthoLQsgaZ/GZ4B5JMN4L4yIxZLOAhZERD/wA+AiSUuBZSRBxsys4YbSYd7ZIWYf89YxFzDKtXQeSERcB1xXduz0ku0XgQ8Oco2ZuRTOzMaUwTrMOyTWRoyYVOvNMKQAImkc8NqIeC6n8piZtcQ2VSb+QdK/MdqG4DbCoKOwJP1Y0maSNgXuBe6TND3/opmZNVZxdNUOM65ln1k3bpCWZL9dJ1Mpe9XE8Z0OHlVkqYHsFhHPSTqBZNb3DGAhMDvXkpmZNVB5J3lxXfKiny4svGrY7j47TuKST76riaUcWbIEkE5JncBU4PyIWK0aWSbNzNrRYOuSV+pA/92fltE3UHDto4osEwm/BzwEbAr8WtIbgGfzLJSZWaNV698orFhVtQM9YF2AsVfLEkCuiYieiDgsnYPxCPDxnMtlZtYQfQMF9jzzl1Xf75BqziAfajqTsSRLAPlp6U4aRC7NpzhmZkNXrXP8hAtu5XOX3VVzZvmaCKYfvEvFDnQYe+lJhqJqH4ikXYHdgQmSjip5azPgNXkXzMwsi0qd45+/7C4+d9ldmT5fXF52wcPLuOS2RzboSB+L6UmGolYn+i7AB0hyTR1ecvx54JN5FsrMLKtKneNZkyCK9cvLnjN1D3rfMGmDxaM8YbC2qgEkIq4Grpb0roi4tYllMjOrqXSVwKFmzC0ScMLeUzYIEKM9/XqjZRnGu1TSl4DtS8+PCHekm1lT9Q0UOPOaxSxfOfRsuaUmju/kjMN3d7AYpiwB5GrgN8D1QPY0lWZmw1C+Fvl+u07mpwsLQ8qWW8nOW27K/C/s25hCjnFZAsj4iDgl95KYmaUqdYxffNsjDbn2yperr2luQ5NlGO/PJR2We0nMzFJDXYt8KDyvo3GyBJB/JQkiL0p6TtLzkpyN18xyk+cvec/raJxBm7Ai4nXNKIiZGcBpfYvqHlk1GM/raKxBA4iSzIknADtExNmStgO2jojbcy+dmY0ZSb/HPaxa3bg+iu6uTiRYsXK153XkIEsn+n8Ca4H3AWcDLwDfAf4ux3KZ2RiyvtN8+MGju6uTmUd4iG4zZAkge0XE2yUNAETEckkb51wuMxvlisN0q2XJrdemm2zk4NEkWQLIakkdpNkBJE0mqZGYmdXltL5Fr8o71SgeZdU8WUZhfRv4GbClpH8Hfgt8JddSmdmo1TdQ4OKcggd4lFUzZRmFdYmkhcD+JOljpkbEH3IvmZmNSjP7FzfkOuOUrOWxeu36UORRVs1VK537Zula6JOAJ4GflLw3KSKWNaOAZjaynda3iJ/8/lHWRCCyZ8qtpauzg3OP2gPA2XNbqFYN5Mck6dwXsuH3vPgz8Mbh3lzSIcB/AB3A9yNiVtn7mwA/At4BPAMcFxEPSToQmAVsDLwMTI+IG4dbHjNrrBMuuJVb/rT+b81GBI8OiXOP2mNdoHDAaJ1a6dw/kH7dIY8bpx3z3wEOBB4D7pDUHxH3lZx2ErA8InaSdDxwHnAc8DRweEQ8Lul/AfMA/xSZtZG+gcIGwaMRijUPB432MGgnuqR/lDShZL9b0tQG3PudwNKIeDAiXiZZJvfIsnOOBOam21cC+0tSRAxExOPp8cVAV1pbMbM20ai+jqKJ4zsdPNpMlmG8Z0TEz4o7EbFC0hlA3zDv3QM8WrL/GLBXtXMi4hVJzwKbk9RAio4G7oyIlyrdRNLJwMkAU6ZMGWaRzaxcedr1Yj9ErXXIh6KrcxznHvUWB442lCWAVKqlZPlc7iTtTtKsdVC1cyJiDjAHoLe3N6+Rg2ZjUt9AgelX3s3qNcl/rcKKVXzusrv4/OXZ1iOvxYs+tb8sgWCBpG+Q9FcAfJqkY324CsB2JfvbpscqnfOYpI2ACSSd6UjalmR+ykcj4k8NKI+ZDdGZ1yxeFzxKRQP+VBs4verfhdYmskwk/BeSkU6Xpa+XSILIcN0B7CxphzQ1yvFAf9k5/cC0dPsY4MaICEndwLXAjIi4pQFlMbM6DHdp2Wp6PBlwRMgykfBvwIxG3zjt0/gMyQiqDuDCiFgs6SxgQUT0Az8ALpK0FFhGEmQAPgPsBJwu6fT02EER8WSjy2lmG+obKDCzf3HD+jjKeTLgyKGoUteU9K2I+Jyka6gwfDsijsi7cI3W29sbCxYsaHUxzEasvoEC06+4e4PZ343QIVgbeDJgm5K0MCJ6y4/XqoH8KP36tXyKZGbtrDRbbofEmoh1XxulQ+JDe23HOVP3aNg1rXlqBZDZJPmvDouIU5pUHjNrA+vX50jWJS8GjUYFj28dt6drGaNArQCytaR3A0dIupQkhck6EXFnriUzs5aZPW/JuuDRaOM7xzl4jBK1AsjpwP8jGV77dTYMIEGyQqGZjSKliQ/z0NkhvnLUW3K5tjVfrQDyREQcKun0iDiraSUys6bLYz3yoq7Ocby4eq07yEehWgHk2yRZcKcCDiBmo1ReI6sgaa667+xDG35daw+1AshqSXOAHknfLn8zIj6bX7HMrFlm9i/OJXgAudRorH3UCiAfAA4ADqYxqUvMrE2UJkDMM0Gcl5cd3WqtB/I0cKmkP0TE3U0sk5nlqHyIbiPss+Mk7nzk2Q2u6Rnlo1+WXFirJN0g6V4ASW+RdFrO5TKznOQxRPeST76Lc4/ag57uLkSSy8prd4x+WbLxXgBMB74HEBH3SPoxcE6eBTOzfDy+YlVDr/fhvZN1dqa+rccBY4zJEkDGR8Tt0gbzCF/JqTxm1kClfR3d4zt5cfWahvZ57LPjJKchGcOyBJCnJe1ImlBR0jHAE7mWysyGrbyvYzip1/fZcRKLH39+XQZeL/ZkkC2AfJpkRb9dJRWAPwMn5FoqMxu2RvV1jO8cxyWffFcDSmSjTZb1QB4EDpC0KTAuIp7Pv1hmVq9Gr9fhuRxWzaABRNIE4Azg79P9XwFnRcSzOZfNzIYoj1nlnsth1WQZxnsh8DxwbPp6DvhhnoUys6HrGyjw+cvvavis8v12ndzQ69nokaUPZMeIOLpk/0xJd+VVIDMbmr6BAl+66h5W5tTUdNP9T+VyXRv5sgSQVZLeExG/BZC0D9DYgeRmVpfT+hZx8W2PDPs6XZ0dVTvcGz1vxEaPLAHkn4G5aV8IwHLgxNxKZGYVlc7p2Ka7i/12ncwlDQgePWma9eLyteXcB2LVZBmFdRfwVkmbpfvP5V4qM9tAeed4YcWqhtQ8BNwyY/3acOU5spzPymoZtBNd0lckdUfEcxHxnKSJkpzGxKyJ8kq5Xlq7mPq2HuezsiHJ0oR1aER8qbgTEcslHQY4oaJZkzRqTkepSrUL57OyocgyjLdD0ibFHUldwCY1zs9M0iGSlkhaKmlGhfc3kXRZ+v7vJW1f8t6p6fElkg5uRHnM2k3fQIF9Zt3Y8Ou6dmGNkKUGcglwg6Ti3I+PAXOHe2NJHcB3gAOBx4A7JPVHxH0lp50ELI+InSQdD5wHHCdpN+B4YHdgG+B6SW+KiMbmqDZror6BAmdes3hdzqquznG8vCZY04Cmqw6Jrx/7VgcMa6hBayARcR5J6vY3p6+zI+KrDbj3O4GlEfFgRLwMXAocWXbOkawPVlcC+ytJC3wkcGlEvBQRfwaWptczG5H6BgpMv/LuDRIerlq9tiHBo3Ocg4flI0sNhIj4BfCLBt+7B3i0ZP8xYK9q50TEK5KeBTZPj99W9tmK/zsknQycDDBlypSGFNys0WbPW8LqNY3vJO/u6mTmEc6aa/nIFEBGsoiYQ5JNmN7e3jyXfzYbkuK8jkpzL4arp7trg+G5ZnloZQApANuV7G+bHqt0zmOSNgImAM9k/KxZ2zqtbxGX3PZIQxd3KvLcDWuWLKOw8nIHsLOkHSRtTNIp3l92Tj8wLd0+BrgxIiI9fnw6SmsHYGfg9iaV22xY+gYKuQUPj66yZqpaA5G0CCr+jAuIiHjLcG6c9ml8BpgHdAAXRsRiSWcBCyKiH/gBcJGkpcAykiBDet7lwH0ky+t+2iOwbKSYPW9JQ4OHR1hZqyj5g77CG9Iban0wIh7OpUQ56u3tjQULFrS6GDYGleaxamTwEPDN4/Z08LBcSVoYEb3lx6vWQEoDRBpMdo6I69OJhKO+892sEZJ1ye/JZVU/ASfsPcXBw1omy4qEnyQZBjsJ2JGkw/q/gP3zLZrZyNY3UOALl99Fo1NYiSSH1fSDd3HwsJbKUpP4NMkkvd8DRMQDkrbMtVRmI1Sj1yMvNQ74hpurrI1kCSAvRcTLyQRwSIfTej6FWZk81iMvEg4e1n6yDOP9laQvAV2SDgSuAK7Jt1hmI8/seUtyCR7gjnJrT1kCyAzgKWAR8L+B63Aqd7NXyWPpVwEfdke5taksKxKuBS5IX2ZWojyD7nA8NOv9r1q21h3l1s7qmUgIwHAnEpqNZI3uLJ84vhPwgk42stSqgXwg/frp9OtF6dcP4050G8NO61vUkPXIizo7xBmH796w65k1y6ATCSUdGBFvK3nrFEl3kvSNmI0pjQ4eTrduI1mWYbyStE9E3JLuvJvWJmE0a4m+gULDgocE3zzWI6tsZMsSQE4CLpQ0gWRQyHLg47mWyqyNNHrdjs5xYvYHnfzQRr4so7AWAm9NAwgR8WzupTJrE0kuq0WsWt2YZM8Tx3dyxuFusrLRIUsurAnAGcDfp/u/As5yILGxYPa8JQ0JHj0ekmujUJYmrAuBe4Fj0/2PAD8EjsqrUGbNVG3uRd9AoSHNVt1dnV5e1kalquuBrDtBuisi9hzs2Ejg9UCsXKUmqs4O0TlOrGxACnb3d9hoUG09kCyjqVZJek/JhfYBGp+zwawFKjVRrV4TDQkePd1dDh42qmVpwvpnYG7JKKxlwIl5FsqsWfLIXwVJ8HCzlY12WUZh3UUyCmuzdP+53EtllqPSPo9xEmsGacYdqq7ODqYfvEtDr2nWjrKMwuoGPgpsD2xUXBckIj6ba8nMclDe59Ho4OHRVjaWZGnCug64jSSde+MXdjZrokYNyy3X2SFmH+P+DhtbsgSQ10TEF3IviVnOGjUsd3znOI56x7bcdP9TTrtuY1qWAHKRpE8CPwdeKh6MiGW5lcqsQRqdhmTnLTdl/hf2bci1zEa6LMN4XwZmA7cCC9PXsCZTSJokab6kB9KvE6ucNy095wFJ09Jj4yVdK+l+SYslzRpOWWz0Kq5RPpzg0ZH2+XVIfHjvKQ4eZiWy1ED+DdgpIp5u4H1nADdExCxJM9L9U0pPkDSJJIVKL8n6Iwsl9ZPUgr4WETdJ2hi4QdKhEfE/DSyfjVClI6xg+AvX/Oncw4ZfKLNRKksNZCmwssH3PRKYm27PBaZWOOdgYH5ELIuI5cB84JCIWBkRNwFExMvAncC2DS6fjUDFEVaFFasIhh88irUPM6ssSw3kb8Bdkm5iwz6Q4Qzj3Soinki3/wJsVeGcHuDRkv3H0mPrpEOMDwf+o9qNJJ0MnAwwZcqUYRTZ2lWj+zmKPrTXdg29ntlokyWA9KWvIZF0PfD6Cm99uXQnIkLSkP9YlLQR8BPg2xHxYLXzImIOMAeSXFhDvY+1t0anWwcYJ/invaZwztQ9GnZNs9Eoy0z0uYOdU+VzB1R7T9JfJW0dEU9I2hp4ssJpBWDfkv1tgZtL9ucAD0TEt+opn40OjZzX4UmAZkOTpQaSh35gGjAr/Xp1hXPmAV8pGaF1EHAqgKRzgAnAJ/IvqrWz4TZbOWiY1a9VAWQWcLmkk4CHSdcakdQLfCoiPhERyySdDdyRfuas9Ni2JM1g9wN3pqlVzo+I7zf9KazlOurMZeWZ42bD15IAEhHPAPtXOL6AklpFRFxIsqBV6TmPkWQFtjGub6BQdy4rBw+z4cuSTHE+8MGIWJHuTwQujYiD8y6cWSV9AwVm9i9mxarVdX2+p7vLwcOsAbLUQLYoBg+AiFguacscy2RW1XBHXTnVulnjZAkgayVNiYhHACS9geHP0TIbkuHM9RgnWBvuMDdrtCwB5MvAbyX9iqTv4b2kE/PMmqGeWsfE8Z2ccfjuDhZmOcoyD+QXkt4O7J0e+lyD82KZ1TTUuR4Tx3cycPpBOZbIzKBGAJG0a0TcnwYPgMfTr1PSJq078y+ejRWlSRC36e5iv10nr1tvYyjtpV2dHZxx+O65ldPM1qtVA/kCSVPV1yu8F8D7cimRjTnlTVSFFau4+LZHMn22u6sTCVasXO2FncyarGoAiYhiP8ehEfFi6XuSXpNrqWzM6Bso8G+X3z3k+RxdnR2ce9QeDhZmLZSlE/13wNszHDPLrN65HALXNMzaRK0+kNeTpE/vkvQ21s/+3gwY34Sy2ShV71yOnu4ubpnhllOzdlGrBnIwcCJJFtyvsz6APA98Kd9i2WhWTwZdTwA0az+1+kDmAnMlHR0RP21imWwUKR1dNSHt8F6+cmjNVp4AaNaesvSBbCtpM5KaxwUkfR8zIuKXuZbMRrzypqp6c1e52cqsPWVZE/3jEfEcyXocmwMfIUnHblZTIxZ76u7qbFBpzKzRsgSQYt/HYcCPImIxTqduGTw+zMWeOseJmUd4UqBZu8rShLVQ0i+BHYBTJb0OWJtvsWykKu3zGJdxsSfx6uyc3V2dzDzCuazM2lmWAHISsCfwYESslLQ58LF8i2UjUXmfR9bJgUHSUV5MY+IOc7ORYdBcWCTBA+CN6fKxZhXV2+fRIbmj3GwEci4sa4i+gUJda3VA9pqKmbWXQXNhRcR+zSuOtavSBZ060r6N4vwMgFOvWlT3tXu6uxpVTDNroixroh9V4fCzwKKIeLLxRbJ2URo0Sju6izWGwopVnHrVIjbZaFzVpqvODkHA6rWVaxmeYW42cmXtRH8XcFO6vy+wENhB0lkRcVFOZbMW6hsoMP2Ku9f94q/WyLRq9Zqa/R6zj3lr8jUdmdU9vpMIeHaV06+bjXRZAshGwJsj4q8AkrYCfgTsBfwacAAZhWb2L65aa8iqp7trXXBwkDAbfbJMJNyuGDxST6bHlgH15aawtjeUtCPjO8fR1dmxwTE3TZmNflkCyM2Sfi5pmqRpQH96bFNgRT03lTRJ0nxJD6RfJ1Y5b1p6zgPpvcvf75d0bz1lsMZZvTY4+h099HR3IZKahxd7Mhv9sjRhfRo4CnhPuj8X+GlEBFDvCK0ZwA0RMUvSjHT/lNITJE0CzgB6SZrgF0rqj4jl6ftHAS/UeX+rothxPhSr1wQ33f+U53KYjTGDBpCICEm/BV4m+UV+exo8huNIks54SALSzZQFEJL1SOanTWVImg8cAvxE0mtZP0/l8mGWZUwrTT3SPb6TF158pa6+j+HmvTKzkWfQJixJxwK3A8cAxwK/l3TMMO+7VUQ8kW7/Bdiqwjk9wKMl+4+lxwDOJpnguHKwG0k6WdICSQueeuqpYRR59CmmHimsWEWQrNNRb8f5Np7LYTbmZGnC+jLwd8U5H5ImA9cDV9b6kKTrgddXud46aQ0n828tSXsCO0bE5yVtP9j5ETEHmAPQ29vrKc8lGpFuHdxhbjZWZQkg48omDD5DhppLRBxQ7T1Jf5W0dUQ8IWlrkpFd5Qqsb+aCZGndm0nmpPRKeoik/FtKujki9sWGJGuzU0eNrLpeLdBs7MoyCusXkuZJOlHSicC1wHXDvG8/UBxVNQ24usI584CDJE1MR2kdBMyLiO9GxDYRsT1Jx/4fHTzqk6XZqauzgw/ttV3FYbrfOm5PbpnxPgcPszEqS01iOkkT0FvS15yIKO/wHqpZwIGSHgAOSPeR1Cvp++l9l5H0ddyRvs4qdqhbY0w/eJeaK4MVh+OeM3UPzj1qDw/TNbMNaPgDqkaO3t7eWLBgQauL0Va2n3FtxeMC/jzr/c0tjJm1JUkLI6K3/Hit9UCep3IKJJH0fW/WwPJZi3R3dVacdT7Ba5Gb2SBqpXN/XTMLYq1RbY0wrx1mZoPJ0oluo9iKlZVzXlU7bmZW5AAyxlUbieWJgWY2GAeQMaBvoMA+s25khxnXss+sG+kbKKx7b/rBuziTrpnVJctEQhthSvNbTejq5G8vv8LqNRuuIgjJGh3FobjF873Ik5ll5WG8o0wxv9VgKUp6urucPdfMMhnyMF5rb6W1jNJaQ9b8Vs6ea2bD5QAyApXXMkqbpbIGBneSm9lwuRN9BKpUy1i1eg2z5y3JnN/KneRmNlwOICNQtVrG4ytWVRxV1TlOTBzf6TxWZtZQbsIagbbp7qJQIYhs093lUVVm1jQOICNA+bKzL1XoJC9tliodnmtmlhcHkDbWN1BgZv/iDZIdLq+QYqS7q5OZR+zuoGFmTeUA0qayzucA2HSTjRw8zKzp3InepoayXrnndJhZKziAtKmhBAXP6TCzVnAAaVNZg4LndJhZqziAtKlK8znKeU6HmbWSO9HbVOl8jsKKVXRIrImgx/M6zKxNOIC0Mc/nMLN25iYsMzOriwOImZnVpSUBRNIkSfMlPZB+nVjlvGnpOQ9ImlZyfGNJcyT9UdL9ko5uXunNzAxaVwOZAdwQETsDN6T7G5A0CTgD2At4J3BGSaD5MvBkRLwJ2A34VVNKbWZm67SqE/1IYN90ey5wM3BK2TkHA/MjYhmApPnAIcBPgI8DuwJExFrg6dxLPIhqKwSamY1WrQogW0XEE+n2X4CtKpzTAzxasv8Y0COpO90/W9K+wJ+Az0TEXyvdSNLJwMkAU6ZMaUDRX63aCoELHl7GTfc/5aBiZqNSbk1Ykq6XdG+F15Gl50VEADGES28EbAv8LiLeDtwKfK3ayRExJyJ6I6J38uTJ9TzKoKqtEHjxbY9QWLGKYH1Q6Rso5FIGM7Nmy60GEhEHVHtP0l8lbR0RT0jaGniywmkF1jdzQRI0bjV4H+sAAAriSURBVAaeAVYCV6XHrwBOakSZ69E3UKi4uFMlxWVnXQsxs9GgVZ3o/UBxVNU04OoK58wDDpI0Me08PwiYl9ZYrmF9cNkfuC/f4lZWbLoaCmfONbPRolV9ILOAyyWdBDwMHAsgqRf4VER8IiKWSTobuCP9zFnFDnWSDveLJH0LeAr4WDMLX+wwz1rzKOXMuWY2WrQkgETEMyQ1h/LjC4BPlOxfCFxY4byHgb/Ps4xF5aOr9tt1Mj9dWMi8VkcpZ841s9HEubBqqDS66pLbHhm0x1+8elTAxPGdnHG4l501s9HDAaSGSqOrBgseXZ0dHP2OHg/fNbNRzwGkhqF2eDvVupmNJQ4gNWzT3VWxo7y8iaqrs8MLO5nZmONsvDVUWhWwq7ODE/aeQk93F8KrAprZ2OUaSA2lqwK6P8PMbEMOIIPwqoBmZpW5CcvMzOriAGJmZnVxADEzs7o4gJiZWV0cQMzMrC5KsqOPDZKeIsn+OxZsQRss9dtEY+l5x9Kzgp+3HbwhIl61It+YCiBjiaQFEdHb6nI0y1h63rH0rODnbWduwjIzs7o4gJiZWV0cQEavOa0uQJONpecdS88Kft625T4QMzOri2sgZmZWFwcQMzOriwPICCNpkqT5kh5Iv06sct609JwHJE0rOf7vkh6V9ELZ+ZtIukzSUkm/l7R9vk8yuAY86zskLUqf6duSlB6fKakg6a70dViznqkSSYdIWpKWc0aF96t+bySdmh5fIungrNdspZye96H0e32XpAXNeZLB1fuskjaXdJOkFySdX/aZij/XLRERfo2gF/BVYEa6PQM4r8I5k4AH068T0+2J6Xt7A1sDL5R95v8A/5VuHw9cNgqe9fb0eQX8D3Boenwm8MVWP19alg7gT8AbgY2Bu4HdsnxvgN3S8zcBdkiv05HlmqPpedP3HgK2aPXzNfBZNwXeA3wKOL/sMxV/rlvxcg1k5DkSmJtuzwWmVjjnYGB+RCyLiOXAfOAQgIi4LSKeGOS6VwL7t/Qvm0Tdzyppa2Cz9HkD+FGVz7faO4GlEfFgRLwMXEry3KWqfW+OBC6NiJci4s/A0vR6Wa7ZKnk8b7uq+1kj4m8R8VvgxdKT2+3n2gFk5NmqJAD8Bdiqwjk9wKMl+4+lx2pZ95mIeAV4Fth8eEUdtuE8a0+6XX686DOS7pF0YbWmsSbJ8r2q9r2p9exD/f43Sx7PCxDALyUtlHRyDuWux3CetdY1a/1cN5VXJGxDkq4HXl/hrS+X7kRESBrR47Bb9KzfBc4m+aVzNvB14OMNura1xnsioiBpS2C+pPsj4tetLtRo5wDShiLigGrvSfqrpK0j4om0OvtkhdMKwL4l+9sCNw9y2wKwHfCYpI2ACcAzQyl3PXJ81kK6XXq8kN7zryX3uAD4eb3lb4Div3vRunJWOKf8e1Prs4Nds1Vyed6IKH59UtLPSJqPWh1AhvOsta5Z8ee6FdyENfL0A8WRRtOAqyucMw84SNLEtHnmoPRY1useA9yYtrG2Ut3PmjZ9PSdp77T9/KPFz6fBqOgfgXvzeoAM7gB2lrSDpI1JOlL7y86p9r3pB45PR/LsAOxM0sGa5Zqt0vDnlbSppNcBSNqU5Gegld/TouE8a0W1fq5botUjFfwa2oukffQG4AHgemBSerwX+H7JeR8n6WRcCnys5PhXSdpN16ZfZ6bHXwNckZ5/O/DGUfCsvSS/SP4EnM/6zAsXAYuAe0j+A2/d4uc8DPhjWs4vp8fOAo4Y7HtD0tT3J2AJJaNxKl2zXV6Nfl6SUU53p6/F7fS8w3zWh4BlwAvp/9Xdav1ct+LlVCZmZlYXN2GZmVldHEDMzKwuDiBmZlYXBxAzM6uLA4iZmdXFAcSaRmkGYEnbSLqy1eVpNUmfkvTRVpejSNLkNCPsgKT3SvqgpD+kWWF7JX17kM9fJ6m7zntPlbRbfSW3VvEwXmsaSS9ExGubdK+NIsktNKT3Ml5bJP931tZdwDYk6XjggIj4RLr/C+CcSJL65X3v/wZ+HhFj/g+LkcQ1EGs6SdtLujfdPlHSVZJ+oWQ9j6+WnHeQpFsl3SnpCkmvTY+fLukOSfdKmlPMGizpZknfUrIexL+W3XOmpIsk3QJclP61/dP0OndI2ic9b7KStUcWS/q+pIclbZGWeYmkH5FM4tpO0vT0s/dIOjP9/KaSrpV0d1q+49LjsyTdl577tZIyfTHd3lPSben7P0tn1Ref6TxJt0v6o6T3Vvk3PUXJGhF3S5o1yDV3TP+9F0r6jaRdJe1JMsn0SCVrapxBkk78B5JmS9pX0s/Tz79W0g/T+90j6ej0+EOStki3P5yW+S5J35PUkR5/QcmaNHenZdtK0ruBI4DZ6fk7Svpsyb/XpfX+rFnOWj1T06+x8yJdgwTYHrg33T6RZA2PCSSzch8myQ20BUkuo03T804BTk+3J5Vc8yLg8HT7ZuA/q9x7JrAQ6Er3f0ySgA9gCvCHdPt84NR0+xCShItbpGVeC+ydvncQMIdkTYZxJPm0/h44Grig5L4TSGbUL2F9jb+7pExfTLfvAf4h3T4L+FbJM3093T4MuL7Csx0K/A4YX/rvU+OaNwA7p9t7kaTPKH4vzi+57s1Ab7q9L0kNAeC84rXS/eL6Kw+l/1ZvBq4BOtPj/wl8NN2Oku/XV4HT0u3/Bo4puebjwCal/15+td/LyRStHdwQEc8CSLoPeAPQTbKA0C1pBWNj4Nb0/P0k/V9gPMlCUotJfmEBXFbjPv0RsSrdPgDYTeuXPNksreG8hyQ/FhHxC0nLSz7/cETclm4flL4G0v3XkuRm+g3wdUnnkfzC/Y2SJHkvkvw1/3PKkjdKmkDyS/JX6aG5JOktiq5Kvy4kCWTlDgB+GBEr03Ivq3bN9BnfnW4XP79JhWvWcgBJXifS+y0ve39/4B3AHek9ulifCPNl1j//QuDAKve4B7hEUh/QN8TyWZM4gFg7eKlkew3Jz6VIFor6UOmJkl5D8hdtb0Q8KmkmSc2l6G817lP63jiS2kT5gj21yln6eQHnRsT3yk+S9HaS2sI5km6IiLMkvZPkF+sxwGeA99W6UZniv0/x32Y4xgErImLPYV6nFgFzI+LUCu+tjrRaQe3neT9Jje5w4MuS9ohh9FtZPtwHYu3qNmAfSTvBur6FN7E+WDyd/jV9TJ3X/yXwL8WdtA8A4Bbg2PTYQSTL5FYyD/h4Sb9Mj6QtJW0DrIyIi4HZwNvTcyZExHXA54G3ll4orX0tL+nf+AjwK7KbD3xM0vi0LJOqXTMingP+LOmD6bmS9NaKV619v08Xd/TqBbluAI5RsjZHcW37NwxyzeeBYkbdccB2EXETSdPlBJIanrUZ10CsLUXEU5JOBH4iqdjEclpE/FHJGh73kqxSeEedt/gs8B1J95D8P/g1yfrTZ6b3/AhJk9lfSH65bfALLCJ+KenNwK1preUF4MPATiSdwWuB1cA/k/xivDqtPQn4QoXyTAP+Kw0CDwIfy/ogaVPbnsACSS8D1wFfqnHNE4DvSjoN6CRZavXurPcDziH5t7uXpBZxJuub2YiI+9Jr/zINBqtJAs7DNa55KXCBpM+SNI/9IG2GE/DtiFgxhPJZk3gYr1mJNFitiYhXJL0L+G7OzT1mI5ZrIGYbmgJcnv7l/DLwyRaXx6xtuQZiZmZ1cSe6mZnVxQHEzMzq4gBiZmZ1cQAxM7O6OICYmVld/j/JPxCbEbFNEgAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5uz6NkjJodVv"
      },
      "source": [
        "What about real data? For some reason, coefficients for linear and logistic regression seem to have no clean relation to each other, but still the initialization makes an x4 speedup. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "5CMhS3jEDpzt"
      },
      "source": [
        "X, y = datasets.load_breast_cancer(return_X_y=True)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "1ZxLCixMsJkO",
        "outputId": "4b4350bf-fc68-4888-8b9d-08b00fecdf86"
      },
      "source": [
        "l1 = LogisticRegression(max_iter=10000).fit(X, y)\r\n",
        "print(l1.n_iter_, 'iterations with default initialization')\r\n",
        "\r\n",
        "linreg = Ridge(alpha=0.01).fit(X, y)\r\n",
        "l2 = LogisticRegression(warm_start=True, max_iter=10000)\r\n",
        "ym = y.mean()\r\n",
        "l2.coef_ = linreg.coef_[np.newaxis, :] / (ym * (1-ym))\r\n",
        "l2.intercept_ = logit(ym)[np.newaxis] - np.dot(X.mean(axis=0), l2.coef_[0])\r\n",
        "l2.classes_ = np.array([0, 1])\r\n",
        "l2.fit(X, y)\r\n",
        "print(l2.n_iter_, 'iterations with initialization from linear model')"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "[2127] iterations with default initialization\n",
            "[554] iterations with initialization from linear model\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 279
        },
        "id": "QEybmq7fsNih",
        "outputId": "92b4964c-5ff0-4415-df61-96cfa1cd3c15"
      },
      "source": [
        "plt.scatter(linreg.coef_, l2.coef_);\r\n",
        "plt.xlabel('linear regression coefficients'); plt.ylabel('logistic model coefficients');"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAfU0lEQVR4nO3df5gcVZ3v8feHEGFASNAEJQEJIhtFBcLOghi8FxEJIEJEZEFFUZTVBxa97mYlwoJydUWzeF2uqERFARFECQNCZPgpqCvKhAAJYNwsC5IJLoMwBGSUkHz3j6pOOpPpmpqZ7q6a7s/refqZ+jV1vj3p1LfrnFPnKCIwMzOrZYuiAzAzs3JzojAzs0xOFGZmlsmJwszMMjlRmJlZpi2LDqARpkyZEjNmzCg6DDOzcWPJkiVPRsTUofa1ZKKYMWMGPT09RYdhZjZuSHq01j5XPZmZWSYnCjMzy+REYWZmmZwozMwskxOFmZllasleT2aN1LW0lwXdK1jdP8C0yR3MmzOTubOmFx2WWcM4UZiNQNfSXuYvWsbA2nUA9PYPMH/RMgAnC2tZrnoyG4EF3Ss2JImKgbXrWNC9oqCIzBrPicJsBFb3D4xou1krcKIwG4FpkztGtN2sFThRmI3AvDkz6Zg4YZNtHRMnMG/OzIIiMms8N2abjUClwdq9nqydOFGYjdDcWdOdGKytuOrJzMwyOVGYmVkmJwozM8tUaKKQdLGkJyQtr7H/IEnPSLo3fZ3d7BjNzNpd0Y3Z3wO+BlyacczPI+LI5oRjZmaDFXpHERF3Ak8VGYOZmWUbD20UB0i6T9JPJb2+1kGSTpHUI6mnr6+vmfGZmbW0sieKe4BdI2Jv4P8DXbUOjIiFEdEZEZ1Tp05tWoBmZq2u1IkiItZExHPp8mJgoqQpBYdlZtZWSp0oJL1SktLl/Uji/WOxUZmZtZdCez1JugI4CJgiaRVwDjARICK+CRwLfFzSi8AAcHxEREHhmpm1pUITRUScMMz+r5F0nzUzs4KUuurJzMyK50RhZmaZnCjMzCyTE4WZmWVyojAzs0xOFGZmlqno0WPNSqdraa/nxDar4kRhVqVraS/zFy1jYO06AHr7B5i/aBmAk4W1LVc9mVVZ0L1iQ5KoGFi7jgXdKwqKyKx4ThRmVVb3D4xou1k7cKIwqzJtcseItpu1AycKsyrz5sykY+KETbZ1TJzAvDkzC4rIrHhuzDarUmmwdq8ns42cKMwGmTtruhODWRVXPZmZWSYnCjMzy+REYWZmmZwozMws04gShaQtJG3fqGDMzKx8hk0Ukn4gaXtJ2wLLgQclzWt8aGZmVgZ57ij2jIg1wFzgp8BuwIkNjcrMzEojT6KYKGkiSaK4LiLW1qtwSRdLekLS8hr7JekCSSsl3S9p33qVbWZm+eRJFBcBjwDbAndK2hV4pk7lfw84LGP/4cAe6esU4Bt1KtfMzHLKkyh+EhHTI+KIiAjg98CH61F4RNwJPJVxyNHApZG4C5gsaad6lG1mZvnkSRRXV6+kyeLKxoSzmenAY1Xrq9Jtm5F0iqQeST19fX1NCc7MrB3UHOtJ0muB1wOTJB1TtWt7YOtGBzZSEbEQWAjQ2dkZBYdjZtYysgYFnAkcCUwG3lm1/Vngo40MqkovsEvV+s7pNjMza5KaiSIirgWulXRARPyqiTFVuw44TdKVwP7AMxHxeEGxmJm1pTzDjK+U9BlgRvXxETHmBm1JVwAHAVMkrQLOASam5/8msBg4AlgJPA98aKxlmpnZyORJFNcCPwduAdYNc+yIRMQJw+wP4NR6lmlmZiOTJ1FsExGfbngkZmZWSnm6x14v6YiGR2JmZqWUJ1F8giRZ/FnSGknPSlrT6MDMzKwchq16iojtmhGIWbN0Le1lQfcKVvcPMG1yB/PmzPQc2WYZ8gwzLknvl/TP6foukvZrfGhm9de1tJf5i5bR2z9AAL39A8xftIyupX48x6yWPFVPXwcOAN6brj8HXNiwiMwaaEH3CgbWbtp5b2DtOhZ0rygoIrPyy9Praf+I2FfSUoCIeFrSSxocl1lDrO4fGNF2M8t3R7FW0gQgACRNBdY3NCqzBpk2uWNE280sX6K4ALgG2FHSF4BfAP/S0KjMGmTenJl0TJywybaOiROYN2dmQRGZlV+eXk+XS1oCvA0QMDciHmp4ZGYNUOnd5F5PZvkpGSVjiB3S9hGxRtLLhtofEVkTDhWqs7Mzenp6ig7DzGzckLQkIjqH2pd1R/EDkmHGl5C2T1TOl66/um4RmplZaWUNM35k+nO35oVjZmZlk+eBu3dJmlS1PlnS3MaGZWZmZZGn19M5EfFMZSUi+knmjTAzszaQJ1EMdUyeB/XMzKwF5EkUPZK+Imn39PUVkgZuMzNrA3kSxd8DLwA/TF9/wbPOmZm1jTwP3P0JOKMJsZiZWQnVTBSSvhoRn5T0EzZ9jgKAiDiqoZGZmVkpZN1RXJr+/NdmBGJmZuWUlSgWkIzvdEREfLpJ8ZiZWclkNWbvJOnNwFGSZknat/pVj8IlHSZphaSVkjZrB5F0kqQ+Sfemr4/Uo1wzM8sv647ibOCfgZ2B80nGeKoI4OCxFJzOcXEh8HZgFXC3pOsi4sFBh/4wIk4bS1lmZjZ6WYni8Yg4XNLZEXFuA8reD1gZEQ8DSLoSOBoYnCjMzKxAWVVPF6Q/GzWu03Tgsar1Vem2wd4t6X5JP5a0S62TSTpFUo+knr6+vnrHambWtrLuKNZKWghMl3TB4J0RcXrjwtrgJ8AVEfEXSX8HXEKNKq+IWAgshGQ+iibEZmbWFrISxZHAIcAcGjNkRy9QfYewc7ptg4j4Y9Xqt4EvNyAOMzPLkDUfxZPAlZIeioj7GlD23cAeknYjSRDHA++tPkDSThHxeLp6FOApWM3MmizPWE8Dkm6VtBxA0l6SzhprwRHxInAa0E2SAK6KiAcknSup8tT36ZIekHQfcDpw0ljLNTOzkak5Z/aGA6Q7gHnARRExK922PCLe0IT4RsVzZpuZjUzWnNl57ii2iYjfDNr24tjDMjOz8SBPonhS0u6kAwNKOhZ4PPtXzMysVeSZqe5Ukm6nr5XUC/wX8L6GRmVmZqWRZz6Kh4FDJG0LbBERzzY+LDMzK4thq54kTUqnP70DuF3S+ZImNT40MzMrgzxtFBcDzwLHpa81wHcbGZSZmZVHnjaK3SPi3VXrn5N0b6MCMjOzcsn7wN2BlRVJs4GBxoVkZmZlkueO4uPAJVXtEk/jJ6TNzNpGnl5P9wJ7S9o+XV/T8KjMzKw08vR6+hdJkyNiTUSskbSDpM83IzgzMytenjaKwyOiv7ISEU8DRzQuJDMzK5M8iWKCpK0qK5I6gK0yjjczsxaSpzH7cuBWSZVnJz5EMtOcWcvqWtrLgu4VrO4fYNrkDubNmcncWUPN1GvW+vI0Zn8pnQ/ikHTT/42I7saGZVacrqW9zF+0jIG16wDo7R9g/qJlAE4W1pby3FEQETcCNzY4FrNSWNC9YkOSqBhYu44F3SucKKwt5WmjMGsrq/uHfp601nazVudEYTbItMkdI9pu1uqcKMwGmTdnJh0TJ2yyrWPiBObNmVlQRGbFqtlGIWkZ6ax2g3cBERF7NSwqswJV2iHc68kskdWYfWTTorC2V7buqHNnTXdiMEvVTBQR8WhlWdKuwB4RcUv6wF2u3lLDkXQY8G/ABODbEXHeoP1bAZcCfw38EfjbiHikHmVbeZSxO+pZXcu44tePsS6CCRIn7L8Ln5/7xkJiMStanrGePgr8GLgo3bQz0DXWgiVNAC4EDgf2BE6QtOegw04Gno6I1wD/D/jSWMu18snqjlqEs7qW8f27fs+6SGpe10Xw/bt+z1ldywqJx6xoeRqzTwVmk8xsR0T8B7BjHcreD1gZEQ9HxAvAlcDRg445mo1Pgf8YeJsk1aHszXQt7WX2ebex2xk3MPu82+ha2tuIYmwIZeuOevmvfz/sdn9erJ3kqUL6S0S8ULk+S9qSoRu5R2o68FjV+ipg/1rHRMSLkp4BXg48WYfyNyiq6qNs9fJFxTVtcge9QySFrO6ojYwxany6K9vLWFVm9VPW/5dFynNHcYekzwAdkt4O/Aj4SWPDGjlJp0jqkdTT19c3ot8touqjcrHp7R8g2HixKfqbaRFxjbQ7atF/u7JVlVn9FP3ZKqs8ieIMoA9YBvwdsBg4qw5l9wK7VK3vnG4b8pj0TmYSSaP2ZiJiYUR0RkTn1KlTRxRIEVUfZb3YFBHX3FnT+eIxb2T65A4ETJ/cwRePeWPNb3FF/+3KVlVm9VP0Z6us8gwKuB74Vvqqp7uBPSTtRpIQjgfeO+iY64APAr8CjgVui6hVMTB6o6n6GKuyXmyKimsk3VEbHeMO20zk6efXDrkdivm8WHOU9f9l0WreUUhaJun+Wq+xFhwRLwKnAd3AQ8BVEfGApHMlHZUe9h3g5ZJWAp8iubupu3o+iZu3kbOsw0SUNa5qjY7xHXvtlLndT263rvHw+S9CVtXTkcA7SUaNvRF4X/r6KUn105hFxOKI+KuI2D0ivpBuOzsirkuX/xwR74mI10TEfhHxcD3KHWykVR+1jKR+s6wXm7LGVa3RMV5/3+OZ2+v1ebHyGQ+f/yJouJocSUsjYtagbfdExL4NjWwMOjs7o6enp+nlzj7vtiGrJKZP7uCXZxy82fay9q4oa1zVGhnjjDNuqLnvkfPeUZcyrLzGw+e/ESQtiYjOofbl6R4rSbMj4pfpypvxYIJDGmn9ZlmHiShrXNXGQ4w2Pvmztbk8ieJk4GJJk0gGBHwa+HBDoxqn3Mg5/mV1g6w0Zo/mnO34DdVax7B3BhGxJCL2BvYG9oqIfSLinsaHNv64fnP8y+oGec47Xz/i87lfvrWCPGM9TZL0FeBW4FZJ56d3FzaIGznHv6xukKP5d3S/fGsFeaqeLgaWA8el6ycC3wWOaVRQ45nrN8e3WtWH00dZfeh++dYK8jRK7x4R56SD9z0cEZ8DXt3owMyKUO/qQ/fLt1aQJ1EMSDqwsiJpNuCvQ9YSBj8gCdS1+tDtVtYK8lQ9fRy4pKrX01PASY0MyqwZao0C+8Vj3jjkcy+j4WlVrRUM+8DdhgOl7QEiYk1DI6qDoh64s/FlpA9ImrWyMT1wJ2ky8AFgBrBlZV6KiDi9jjGaNZ0bms3yyVP1tBi4i2SY8fWNDae1+cGrcvEDkmb55EkUW0fEpxoeSYsbyaxoTijNMW/OzE3+TcANzTY+NfqakafX02WSPippJ0kvq7zqFkGbyPvglZ/kbR4/IGmtoBnXjDx3FC8AC4Az2ThXduBnKUYkb314VkLxBaz+/ICkjXfNuGbkSRT/ALwmIp6sS4ltKm99uBtYy8HVfzZeNOOakafqaSXwfN1KbFN5H7zyk7zFc/WfjSfNuGbkSRR/Au6VdJGkCyqvukXQJvLWh/tJ3uJ5ID8bT5pxzchT9dSVvmyM8tSH+0ne4rn6z8aTZlwzhk0UEXFJ3UqzXNzAWiw/X2HjTaOvGZ7S1GwQV/+ZbSpP1ZNZW3H1n9mmCkkU6QN7PyQZP+oR4LiIeHqI49aRDB0C8PuIOKpZMVp7c/Wf2UZ5pkK9OR0YsLK+g6TuMZZ7BnBrROxBMsXqGTWOG0jn6N7HScLMrBh52iimRER/ZSX95r/jGMs9Gqg0kl8CzB3j+czMrEHyJIr1kl5VWZG0KxuH8hitV0TE4+nyH4BX1Dhua0k9ku6SlJlMJJ2SHtvT19c3xvDMzKwiTxvFmcAvJN1BMsPdW4BThvslSbcAr6xxvg0iIiTVSjy7RkSvpFcDt0laFhH/OdSBEbEQWAjJxEXDxWdmZvnkeY7iRkn7Am9KN30yz7hPEXFIrX2S/lvSThHxuKSdgCdqnKM3/fmwpJ8Bs4AhE4WZmTVGzaonSa9Nf+4LvApYnb5elW4bi+uAD6bLHwSuHaL8HSRtlS5PAWYDD46xXDMzG6GsO4pPkVQxnT/EvgDGMqnwecBVkk4GHgWOA5DUCXwsIj4CvA64SNJ6koR2XkQ4UZSAR1Y1ay+KyK7Ol7R1RPx5uG1l0tnZGT09PUWH0ZIGz9QHyVPLnvDHbHyTtCQiOofal6fX07/n3GYtpGtpL7PPu43dzriB2efdtmGIbY+satZ+alY9SXolMB3okDSLpMcTwPbANk2IzQqSNb+3R1Y1az9ZbRRzgJOAnUnaKSqJ4lngM40Ny4qUddfgkVXN2k/NqqeIuCQi3gqcFBEHR8Rb09dREbGoiTFak2XdNXhkVbP2k6eNYmdJ2yvxbUn3SDq04ZFZYbKmVsw7U5+ZtY48T2Z/OCL+TdIc4OXAicBlwE0NjcwKM2/OzCF7NlXuGjyyqll7yZMoKm0TRwCXRsQDkpT1Cza+lX0+Bj/HYdZceRLFEkk3AbsB8yVtB6xvbFhWtLLeNWT1yCpjvGatIE8bxckk80X8TUQ8D7wE+FBDozKrwc9xmDVf1nMUr42I3wL7pJte7RonK5qf4zBrvqLGejIbFT/HYdZ8NRNFRJyS/nxr88IxyzZcjywzq79hG7MlHTPE5meAZREx5DwSZo1S9h5ZZq0oT6+nk4EDgNvT9YOAJcBuks6NiMsaFJvZkMraI8usVeVJFFsCr4uI/waQ9ArgUmB/4E6Sh+/MzKxF5ekeu0slSaSeSLc9BaxtTFhmZlYWee4ofibpeuBH6fqx6bZtgf6GRWZmZqWQJ1GcChwDHJiuXwJcHcnUeO4RZWbW4oZNFBERkn4BvEDy/MRvYrj5U83MrGUM20Yh6TjgNyRVTscBv5Z0bKMDMzOzcshT9XQmyThPTwBImgrcAvy4kYGZmVk55EkUWwx6sO6P5OstZTYueRhzs03lueDfKKlb0kmSTgJuABaPpVBJ75H0gKT1kjozjjtM0gpJKyWdMZYyzfKoDGPe2z9AsHEY866lvUWHZlaYYRNFRMwDFgJ7pa+FEfHpMZa7nKQn1Z21DpA0AbgQOBzYEzhB0p5jLNcsk4cxN9tcnqonIuJq4Op6FRoRDwEMM2z5fsDKiHg4PfZK4GjgwXrFYTaYhzE321zNOwpJz0paM8TrWUlrmhDbdOCxqvVV6bYhSTpFUo+knr6+voYHZ62p1nDlHsbc2lnNRBER20XE9kO8touI7Yc7saRbJC0f4nV0fd/ChngXRkRnRHROnTq1EUVYG5g3ZyYdEydsss3DmFu7y1X1NBoRccgYT9EL7FK1vnO6zaxhPIy52eYalijq4G5gD0m7kSSI44H3FhuStQMPY262qUKeh5D0LkmrSOa5uEFSd7p9mqTFABHxInAa0A08BFwVEQ8UEa+ZWTsr5I4iIq4Brhli+2rgiKr1xYzxmQ0zMxsbP2FtZmaZnCjMzCyTE4WZmWVyojAzs0xOFGZmlsmJwszMMjlRmJlZJicKMzPL5ERhZmaZnCjMzCyTE4WZmWVyojAzs0xOFGZmlsmJwszMMjlRmJlZJicKMzPL5ERhZmaZnCjMzCyTE4WZmWVyojAzs0xOFGZmlqmQRCHpPZIekLReUmfGcY9IWibpXkk9zYzRzMwSWxZU7nLgGOCiHMe+NSKebHA8ZmZWQyGJIiIeApBURPFmZjYCZW+jCOAmSUsknZJ1oKRTJPVI6unr62tSeGZmra9hdxSSbgFeOcSuMyPi2pynOTAieiXtCNws6bcRcedQB0bEQmAhQGdnZ4wqaGsZXUt7WdC9gtX9A0yb3MG8OTOZO2t60WGZjUsNSxQRcUgdztGb/nxC0jXAfsCQicKsomtpL/MXLWNg7ToAevsHmL9oGYCThdkolLbqSdK2krarLAOHkjSCm2Va0L1iQ5KoGFi7jgXdKwqKyGx8K6p77LskrQIOAG6Q1J1unyZpcXrYK4BfSLoP+A1wQ0TcWES8Nr6s7h8Y0XYzy1ZUr6drgGuG2L4aOCJdfhjYu8mhWQuYNrmD3iGSwrTJHQVEYzb+lbbqyWy05s2ZScfECZts65g4gXlzZhYUkdn4VtQDd2YNU2mwdq8ns/pworCWNHfWdCcGszpx1ZOZmWVyojAzs0xOFGZmlsmJwszMMjlRmJlZJkW03vh5kvqAR4uOYwSmAK0050YrvR+/l/JqpfdThveya0RMHWpHSyaK8UZST0TUnOlvvGml9+P3Ul6t9H7K/l5c9WRmZpmcKMzMLJMTRTksLDqAOmul9+P3Ul6t9H5K/V7cRmFmZpl8R2FmZpmcKMzMLJMTRUlIWiDpt5Lul3SNpMlFxzQWkt4j6QFJ6yWVtttfFkmHSVohaaWkM4qOZ7QkXSzpCUnjfiphSbtIul3Sg+nn6xNFxzQWkraW9BtJ96Xv53NFxzQUJ4ryuBl4Q0TsBfwOmF9wPGO1HDgGuLPoQEZD0gTgQuBwYE/gBEl7FhvVqH0POKzoIOrkReAfImJP4E3AqeP43wXgL8DBEbE3sA9wmKQ3FRzTZpwoSiIiboqIF9PVu4Cdi4xnrCLioYhYUXQcY7AfsDIiHo6IF4ArgaMLjmlUIuJO4Kmi46iHiHg8Iu5Jl58FHgLG7cQjkXguXZ2YvkrXw8iJopw+DPy06CDa3HTgsar1VYzjC1IrkjQDmAX8uthIxkbSBEn3Ak8AN0dE6d6PZ7hrIkm3AK8cYteZEXFtesyZJLfXlzczttHI837MGkHSS4GrgU9GxJqi4xmLiFgH7JO2S14j6Q0RUar2JCeKJoqIQ7L2SzoJOBJ4W4yDB1yGez/jXC+wS9X6zuk2K5ikiSRJ4vKIWFR0PPUSEf2SbidpTypVonDVU0lIOgz4J+CoiHi+6HiMu4E9JO0m6SXA8cB1BcfU9iQJ+A7wUER8peh4xkrS1EoPR0kdwNuB3xYb1eacKMrja8B2wM2S7pX0zaIDGgtJ75K0CjgAuEFSd9ExjUTaseA0oJukwfSqiHig2KhGR9IVwK+AmZJWSTq56JjGYDZwInBw+v/kXklHFB3UGOwE3C7pfpIvJzdHxPUFx7QZD+FhZmaZfEdhZmaZnCjMzCyTE4WZmWVyojAzs0xOFGZmlsmJwupO0nPpz2mSflx0PEWT9DFJHyg6joq07/6vJS2V9JZ0pN+H0lFZOyVdMMzvLx7t6MaS5o7zQfzakrvHWt1Jei4iXtqksrasGkwx976c5xbJ/5H1ow6whCQdDxwSER9J128EPh8Rv2hC2d8Dro+Itv8CMZ74jsIaRtKMyhwIkk6StEjSjZL+Q9KXq447VNKvJN0j6UfpOD5IOlvS3ZKWS1qYXriR9DNJX5XUA3xiUJmflXSZpF8Cl6Xfnq9Oz3O3pNnpcVMl3ZzOAfBtSY9KmpLGvELSpSTDKOwiaV76u/dX5guQtK2kG9J5BJZL+tt0+3npXAn3S/rXqpj+MV3eR9Jd2jjvyA5V7+lL6dwEv5P0lhp/009LWpaWe94w59w9/XsvkfRzSa+VtA/wZeDo9GG1c4ADge8omRPlIEnXp7//UknfTcu7X9K70+2PSJqSLr8/jfleSRcpGZ4dSc9J+kIa512SXiHpzcBRwIL0+N0lnV7197pytJ81a7CI8Muvur6A59KfM4Dl6fJJwMPAJGBr4FGSsZSmkMxZsW163KeBs9Pll1Wd8zLgnenyz4Cv1yj7s8ASoCNd/wFwYLr8KpKhHyB5En5+unwYydDOU9KY1wNvSvcdSjLxvUi+WF0P/C/g3cC3qsqdBLwcWMHGO/XJVTH9Y7p8P/C/0+Vzga9Wvafz0+UjgFuGeG+HA/8ObFP998k4563AHuny/sBtVf8WX6s678+AznT5IJJv/ABfqpwrXd8h/flI+rd6HfATYGK6/evAB9LlqPr3+jJwVrr8PeDYqnOuBraq/nv5Vb6XBwW0Zro1Ip4BkPQgsCswmWRioF+mNwwvIRluAuCtkv4J2AZ4GfAAyYUJ4IcZ5VwXEQPp8iHAnum5AbZP71gOBN4FEBE3Snq66vcfjYi70uVD09fSdP2lwB7Az4HzJX2J5ML6c0lbAn8m+XZ+PUlS2UDSJJKL4R3ppkuAH1UdUhngbglJwhrsEOC7kY4FFhFP1Tpn+h7fnC5Xfn+rIc6Z5RCSMa5Iy3t60P63AX8N3J2W0UEyVDbAC2x8/0tIxjAayv3A5ZK6gK4RxmdN4kRhzfSXquV1JJ8/kYxvc0L1gZK2JvmG2hkRj0n6LMmdSMWfMsqp3rcFyd3BnwedPyvO6t8X8MWIuGjwQZL2Jfn2/3lJt0bEuZL2I7mAHksyVtTBWQUNUvn7VP42Y7EF0B8R+4zxPFkEXBIRQ83GuDbS2wSy3887SO7Q3gmcKemNMYZ2JWsMt1FY0e4CZkt6DWyo+/8rNiaFJ9Nvx8eO8vw3AX9fWUnr6AF+CRyXbjsU2KHG73cDH65qN5kuaUdJ04DnI+L7wAJg3/SYSRGxGPg/wN7VJ0rvpp6uan84EbiD/G4GPiRpmzSWl9U6ZyRzNPyXpPekx0rS3kOeNbu8UysrlbaPKrcCx0rasRKPpF2HOeezJINfImkLYJeIuJ2kynESyR2blYzvKKxQEdGnZB6OKyRVqkbOiojfSfoWSYPyH0hG1hyN04ELlYzOuSVJe8jHgM+lZZ5IUtX1B5KL2CYXqoi4SdLrgF+ldyHPAe8HXkPSKLseWAt8nOQCeG16NyTgU0PE80Hgm+nF/mHgQ3nfSFpFtg/QI+kFYDHwmYxzvg/4hqSzSKbYvBK4L295wOdJ/nbLSe4KPsfG6jEi4sH03DelF/21JInl0YxzXgl8S9LpJNVa30mrzwRcEBH9I4jPmsTdY60tpUlpXUS8KOkA4BsNrqYxG7d8R2Ht6lXAVek34ReAjxYcj1lp+Y7CzMwyuTHbzMwyOVGYmVkmJwozM8vkRGFmZpmcKMzMLNP/AFk/rt61y+oRAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "0SF3WRJlsjGx"
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}