{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as pp\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATQAAAD8CAYAAAD5TVjyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJztnXecVNXZgJ+XXRYUjCyCCigCtmA0iqwoJpZYEVH0s7cQI/JhiUYxiSX2xBIJRiMWxJYvVlAEaYoUsQKLgoKIIF1QugLKLrv7fn/MvTgsuzuzc+89t8x5fr/5zdwy91yGmWdPec97RFWxWCyWJNAg7BuwWCwWv7BCs1gsicEKzWKxJAYrNIvFkhis0CwWS2KwQrNYLImhMOwbsFgsyUVEngZ6ACtV9cAajv8JuMjZLAQ6Ai1Vda2ILAI2AJVAhaqWZCzPxqFZLJagEJGjgY3Af2oSWrVzTwOuU9XjnO1FQImqrs62PNvktFgsgaGqk4G1WZ5+AfCil/KMNjlbtGih7dq1C7YQ3QSVy6CgLcgOgRZVqVtYV76MnQpb0qigaaBlVamy5IfVNC1sTItGOwVaFsCm8nKqVBERmhQVIQGXt658E+vKN9K2SUsKJdi/s+VVP/D9lm9oVtSGQmkUaFloOVQuhIJWID8Ltixg+vTpq1W1pZdrnHx8a12zpiy78maunQ1sTts1SFUH1bdMEdkR6AZcnbZbgbdERIEnsrmuUaG1a9eO0tLSwK6v5dPQdb2hwa+R5v9BCjz9v9bJ+vLlvLqkHxVVHfiftv1p2XjvwMrauGUz15Q+Q6Pvv+a+Qy7kmN0OCKysiqoqbhg3hhFzv9i6r+te7Xj81J40Kgzu6/Lx2oX8cfqztGrcjIFdegcq7Z/+7zYH/n+nVRvQdb+HLbORZvchjU8IrCwAEVns9RprvmtN6Yzsfqcisjmbvq0sOA14X1XTa3O/UtXlIrIrME5EvnBqfLWSmCbnTzJrbVBmwf8gXJnNMSyzaw/vyrDzLuLWo3/D5MWL6DtqOGUVFYGVfWjz9vyr8+9YsXk9V00dzOqyDYGV1ayoNWe1/SeFDRrz2pIbWLX5q8DKkgY7IcVPQ8NfoOuvQTe/HVhZMed8qjU3VXW587wSGAZ0yXSRRAgtH2R2r0GZ/enIX3Pt4Udy8G67c+khh3Lv8SfxjpVazlip1Y2I7AwcAwxP29dERHZyXwMnAbMyXSv2QgtPZg8YrZkda1BmV5Qcvs3xc39x0Fap/e9IK7Vc2F5q4wIrK0qIyIvAh8D+IrJMRC4Tkb4i0jfttDOBt1R1U9q+3YD3RGQmMBUYpapjM5ZnMmyjpKRE/exDy4eamclmZk0yS+fl2Z9x8/i3OKptO57oYfvUcmHbPrWHfe9TE5HpXvu06vM79aM8P4ltDc3KzDv1kRnAeU5N7d0ltqaWK7b5GSyxFJqVmXfqKzOXc63UPGOlFhyxE5qVmXdylZlL/khtfmBlWakFQ6yEZmXmHa8yc0mXWnJHP/9kBwpiRmyEZmXmHb9k5uJKzcap5c62UrvWSs0jsRCalZl3/JaZS7rUkt38tFKLA5EXWj7IzHTQrF8yczHfp9YrD6Rm+9RyIdJCyxeZhRk06xdba2pGpNYhD6RmBwpyIbJCy4cZACZldkPX4GTmYqXmHSs1b0RSaOHWzPYJrKwwpzNdeViwMnMxG3zbIQ/61KzU6kPkhJYPzcy4DgBki51R4B0rtdyIlNCszLwTtsxc8if41kotSkRGaFZm3omKzFzyI/jWSi1KREJo28rs/6zMciBqMnOxcWresVLLntCFtr3MWgRWlpVZONjmp3es1LIjVKFt38y0MqsvUZeZi5Wad6zUMhOa0GyfmXfiIjMXKzXvWKnVTShCszLzTtxk5mKl5h0rtdoxLjQrM+/EVWYuVmresamHasas0HSTlZlH4i4zFys179gJ7dtjVmgVixM9N9PKrH5YqXmnek0t3zErNGmU6LmZSZLZuo0/8s3aDaz5flPmkz1gpeadn6R2cGBlxAWzQivc2zYzc8C0zJ4aM5Xj//Q43W8ZzIl/GUT/IZMIcrlDKzXvSIOdkObPB3b9uBB6YK1fWJn5w+AxUxg44n2OOqg9t118It1K9ueFCZ8Yk5qZ1EMJlZok5uecM4n4BKzM/GHwmCk8OuIDunf5OQP6ns4ZvzqQv//+FC48rhMvTpyRQKklL59avhN7oVmZ+UO6zO7sdTIFDVJfDRGh39nHcNHxhxqRmvl8alZqSSLWQrMy84faZOYiIlx/1tHGamo2SaQlV2IrNCszf8gkMxe3pma6+WlTD8UbEXlaRFaKyKxajh8rIt+JyAzncVvasW4iMldE5ovIjdmUF0uhWZn5w1NjpmYlMxfTzU+77mcieBboluGcd1X1EOdxF4CIFAADgVOAA4ALRCTjjy92QrMy8wd3NDNbmbm4zc+fpPaOlVoO5IvUVHUysDaHt3YB5qvqAlUtB14CemZ6U2EOBYWGXTfTH7JtZtaGKzWA58d/DMAN5xyDiPh+r5CSGsCN49+i76jhPH5qTxoVBvPVdUc//zj9Oa6aOpiBXXrTotFOgZTlSu3VJf14bckNgX+ns2Vd+VKGLr7eZJFdRWQmsBy4QVVnA22ApWnnLAMy/jBiU0MLs2aWhHUzXbzKzGXbmpq5OLV3jGS+taOf9aCFiJSmPfrU8/0fA3up6sHAv4HXnf01/XXM/AVTVWOPzp07ay6sK/taB887Xx+fe4au/HF+TtfIlg3lP+qlHzyqR4y9RSd9MzvQsrZUVuq1Y0dq+4f666PTPgq0LFXVJ0d/pJ36DtBbnh6tFZWVvlyzqqpK+w+ZpJ36DtB/vDxBq6qqfLlubbw061Pt8FB/7TVsqG7esiXQsqavWaBHvXWbnjt5gK7a/H2gZfn1HQdK1eDvNJvygHbArGzKBhYBLYCuwJtp+28Cbsr0/qz+PIvIdSIyW0RmiciLItJYRNqLyBQRmSciL4tIUTbXqi+2z8wf6jsAkC3b96klKaRj2z61Nbam5jsisrs4fRUi0oVUq3ENMA3Y1/FMEXA+MCLT9TJ+q0WkDXANUKKqBwIFzsXvBx5U1X2BdcBluf2TasfKzB9yHQDIFtNSC2vu55VWavVGRF4EPgT2F5FlInKZiPQVkb7OKWcDs5w+tIeB853KXwVwNfAmMAd4RVN9a3WTRRXQ7ZxrTmoQYSRwMrAaKHTO2aZ66EdVNqnNzArDzczBo6dop74D9GYfm5m1sW3zc2Lgzc+XjTY/v4pF85MINjlNPjL+qVbVr4H+wBJgBfAdMB1YrymLQmoEok0thu7jdhiuWrUqo2Ah2TWzfiHVzO4KoGZWne0HCsyEdNhpUhaXbJqcxaTiP9oDrYEmpILdqlPjN1dVB6lqiaqWtGyZOXVQkmUWx9HM+rLtNClzo59WahbILmzjBGChqq5S1S3Aa8CRQDMRcYOB9iAVQ+IJKzN/CEtmLmHNKLBSs2TzTV8CHCEiOzqjEccDnwMTSXXoAfQChnu5ERs06w9hy8wl2QMFdkJ7VMmmD20KMJRUANxnznsGAX8BrheR+cAuwFO53kSYMrNBs8ERltRsPrX8JatvvKrerqo/V9UDVfUSVS3T1ByrLqq6j6qeo6pludyAnQHgD1GTmYvp1EPpUgt+7qetqUWNUL/1ts/MH1yZnRIxmbmYTj10Xto0KTuhPb8I7ZtvZeYP6TUzE6EZuRJW5lszcz+t1KJCKN9+KzN/CGo6U1Dk0zQpK7VwMP4LsDLzh6fGTA10OlNQmM6nZqWWXxj9FVTqFiszH4irzFzCSj1kpZZ8jP4S1pUvtTLzSNATzU2R/JCOcKSW7xj9NShqZeaBqIZm5IqVmj+kSy3fMfqLaNGovZVZjpiU2XvvfMH5pz/EmSf355xTH2TsyBmBlRXWQEESpXZx+ycDu35cMCo0CbC4fJCZiTizdyd9wd/++hrFzZtwQreDaLNnMf+8ZySjh38cWJnJXvfTnNQaFTQN7NpxId5tFod8kVnQcWbvTvqCv9/6Gvt3bE3/gZdw1XUn84+HLqZL13148P7RgUstjODbpEkt34m90PJlorlJmd3z4AU0adIIgKJGhdx+z9k/SW3EJ4HdQ7KzdFipmSDWQkvy3EyTQbO1ycxlG6ndN8pI89P0up9WaskgtkJLcjPTZJxZJpm5mK6pJTfzrZVakMRSaFZm/pCtzFzCq6nZ4FtLdsROaEmWmcmg2frKzMWV2mFd9zYyUGDj1Cz1IVZCS7rMotJnlomiRoXccc85VmoesVLzn9gIzcrMH7zKzGU7qRnoU7MhHZZMxEJo+SAzU0GzfsjMxZWaqT41G6dmyUTkhZYPMgszzswrpvvUTErNDhTEj0gLLV9kZmJuZhAyczHd/LTBt5baiKzQrMz8wZ2bGZTMXLZrfhqLU7NSs/xEJIWWL9OZTPWZ7ffzVoHKzGWb5qeVWs6YXCIvaUROaHbdTH9wm5n7/bwV9/7rwsBl5rJN89P4NKkkSc3cCu1BIiJPi8hKEZlVy/GLRORT5/GBiBycdmyRiHwmIjNEpDSb8iIltHyQmckUQPt3bG1UZi5JDumwUqs3zwLd6ji+EDhGVX8J3E1qEfN0fqOqh6hqSTaFFeZ0iwGQD31mplMAmWhm1oYrtTtuHsKD940CVbr3PDSQstyBAoAXJqTkecM5xyIigZR37i8OAuCm8W/xvyOH80SPnjQqDOan5C5m/Mfpz3LV1MEM7NKbFo12CqQsl8WbVtN3ij/JIlV1soi0q+P4B2mbHwF7eCkvEjW0fJBZnEMzciWM0c9k1tQiPVDQQkRK0x59PFzrMmBM2rYCb4nI9Kyvq6rGHp07d9bqbCj/US/94FE9YuwtOumb2dsd95MtlZV67diR2v6h/vrotI8CLUtVdfDoKdqp7wC95enRWlFZGWhZ706aoyf/+u96zeXP6MaNmwMtq76Ubd6iN1//op7Q9W4d9fr0QMuqqqrS/kMmaae+A/QfL0/UqqqqQMt7edan2uGh/vq714fq5i1bAi1r+poFetRbt+m5kwfo6s3f13gOUKoB/E5rI5vygHbArAzn/AaYA+yStq+187wrMBM4OlNZodbQklwzM5k14713Un1mpkYz60t4E9rNpR6avHgRfUeZq6ldOXUwa6JVU8sZEfklMBjoqapr3P2qutx5XgkMA7pkulZoQrMy84ewBwCyJbwJ7WalZrL5eWX0mp/1RkTaAq8Bl6jql2n7m4jITu5r4CSgxpHSdEIRWpJlFocUQGERrtRsn1oYiMiLwIfA/iKyTEQuE5G+ItLXOeU2YBfg0WrhGbsB74nITGAqMEpVx2Ys0Gt7u75t8/Q+s4kJ6zN7cvRHxvrMJk+Mbp9ZJsLrU5uQ2D61VU6fGhHsQzP5MFpDq1JN7BoASZqbGTThpfM2mE/NcPMz6jU1UxgV2pIfVid6OpPpoNk4ysxlu3TeNvg2J6pLLd8xKrTNleWJrZnlU9CsX5heoyAf4tTyHaMzBdo12TWRNbOkDQDMn7GQp295kbIfyigsKuSiW87il0cH8//mSu3Om4fy4P2jQYTup3cKpCxXaiLC8+NT8kzOjIL2PHPEFezDXYFcPy5k9QsUkWYiMlREvhCROSLSVUSai8g4EZnnPBdnuk7jgobe77gWkiwzk31m8z9ZyJ9PuIt50xcAsHTO19zS/R5mTpodWJl23U9/2Hun3QO7dmzIZuQAeA7o7bwuApoB/wBudPbdCNzv5+hJfbCjmf4w7+MFembzXnrhXn11+YJvVFV17Tfr9LJf/FF7NLlIZ0ycFWj5ZZu36E3Xv5DoGQW9hgU7+kmej3JmI7OfkZoRL9X2zwVaOa9bAXP9/KCyxcrMH2qSmYsrtVN3vDDBUjMX0hGk1KzQMgvtEFKBbc8Cn5CaotAEWF/tvHW1vL8PUAqUtm3bNusPKhuszPyhLpm5mJZa0uPUgpKaFVpmoZUAFcDhzvZDpPIWZSW0XD+oTFiZ+UM2MnOxNTV/CDL4Nt+Flk1P9jJgmapOcbaHAocC34pIKwDneWUW1/KFfBgAMDHRPDUAcCc77LQD/SfcQav2u9V5fvFuzXhg/O3s3n5Xbjk1+IEC09OkTId0mJjQnm9k/GWq6jfAUhHZ39l1PPA5MALo5ezrBQwP5A6rkS6zG7omM2jWRNrsbWQ28Q5adahbZi5JlloYcWomZhTkFdlU40j1o5UCnwKvA8WkJpSOB+Y5z839rMrWhG1m+kN9mpm1kR+jn/FrfpLnTU6jhXkRmpWZP/ghM5e1364PR2rDPw60rKqqKn3glYmxlFq+Cy0SKbgzEUY+s0TOAKjezMzQZ5aJ4l135oEJd7Bbu5b8tce95pqfBuZ+btv8TE4+taQTeaFVJjw5Y1xl5lK86848MP72xErNdJJIEzMKkkykhVZZVUW/hCZnND+dyX+ZubgDBaFILYFJIidbqeVMZIVWEYLMTIZmmFoDIGiZuaRLzfjoZxLzqVmp5UQkhZbkPrN0mRkPzQhIZi7FuzXjgQl3mA/pMJRPzZTUzrNSy5nICS3pawCYWtDEtMxc3D61UKRmNEuHGanZPrX6ESmhJb1mZrTP7MS7sp4B4Dfpwbem+tSSmM47XWp2RkF2REZo+VIzM9Zn1rRxSmZZzgDwG7f5aWqgYOu6nwmrqdlpUvUjEkKzNTN/yHU6U1CEFqdmdIk8O00qSoQuNFsz84ew+swyEUac2tbmZwKlZvvU6iZUoSU5a0YSgmb9wnSc2tbmp5Va3hGa0KzM/CHqMnMJLfjW+ECBoWlSVmo1EorQkr5uppVZzYQhNfPrfpqbJmWltj3GhWbXzfSHuMnMJYzRT/PrfgYvNRt8WwsmU3sc2rlzYlMAvTspnimAwsJ0PrWtaxQYSD1kcjWpl6qtUUDE0gcBT5PKZj2rluMCPAzMJ5Vv8dC0Y71I5VucB/TK5t6NCq3VfvsZk9lr735qTGafzVxiTGarvl6jZ+7yu1jLzCVdaos+XxpoWelSe3/y3EDLSpfak6OD/667+dT+MPqNKArtaFIp+2sTWndgjCO2I4Apzv7mwALnudh5XZzp3o2unL5pSzkn79Uu8GYmwLS5S2nWpHHgfWYAs2YsobKyirsfODfQZibAgpmL2bB2I7e+cn1smpm1UbxbM/768vVcftD1fDFlHnt13COwstzm51mn/JOZHy/iyKP2C6wst0/t/dmLmDZ3Kb1PCfb7fu4vDmLa8q+ZuHCBL9dbsG4dF7z6si/XUtXJItKujlN6Av9x5PiRs6h5K+BYYJyqrgUQkXFAN+DFusozKjSARgXmiixqWBi4zNJpWGTw39a4yFhZQVLUuKGxshoWFSANxEhZIkLjhga/DwUFxsqqRgsRKU3bHqSqg+rx/jbA0rTtZc6+2vbXiVGhNSooZPzCrxgz/0tO2Se4v5AAHVrtwtjSuTw5+iMu735EoGW1bd8SgAH3juTG286goDA4ibbeezcKGxbw2HXPcN+bt9K0WZPAygqash/LeOjKJxER9ty/daBlqSoDH3yTH38oZ692LQItC2Dou5/yxdKVnHfsIYGX9cHSJQz74nN+uevuTPfheh2Ki3nxrPOyOvclzl+tqiUeiqvpL4zWsb9OjI5y7vGzn3HI7q24ZsxIRs/7MtCyLu12GD0O78hjb3zIk6M/CrSsI4/aj95XHsektz/nvrtep7KiKrCy9tivNbcNvYGvZizipm5/Y+P6TYGVFSRlP5Zx2xn/4JO3P6PfU1dwQNf9M78pR1yZDR9ayjkXHsEpp3cKrCxIyeyeF8bz6wPbc93/HBVoWR8sXULvN4bRdudmDOx+WqBlBcQyYM+07T2A5XXsrxuvHYj17WzcUFamZ7/ygu7z8D919LxgO2crKiv11mfGaKe+A/SJkR8GWpaq6kv/976e0PVu/dutr2rFlmAHIt4fPlW7FZ2nVx9+o25YtzHQsvxm8w+b9c8n3aUnNjhHxz4zIdCyqqqq9OH+Y/SErnfr4w+PC3zUccg7M7RT3wH6h0eGaVm5/yujp/P+ksXaceC/9OT/PqurNm1S1cyd9Nk8/F4kBWhH7YMCp7LtoMBU/WlQYCGpAYFi53XGleWMC01VQ5PaoFHJktoHI6bFTmphyeyJf5uT2TUDg5fZe0sWbZXZakdmqtETGqlO/BXAFqfWdRnQF+jrHBdgIPAV8BlQkvbe35MK55gPXJrNvYciNNVtpTbqy2RJ7eX/fmClVgOJltnkmaHLTDV6QjP9CG0uZ9OiIp7peRaH7N6Ka8eOZMz84PrUCho04PbfnmSsT+3ci7oa61PreloJtw29gfmfLIx0n1r1PrOTf/ebwMpSTfWZjXi1lLMvOILLrzoekeBGN9P7zB64vAdFAY5uvr90MZe/8Tptd27G82eewy477hhYWXEk1Gwb6VK7ZkyypHbexUcaldqtQ/pFVmplP5Zx+5kPGJXZ8KEpmfW52ozMjjqoPf37BCuzD5Yu2Sqz/1qZ1Ujo+dDCkNppRxyQOKkdefphkaypuTL7eNynRmT2yADzMjNRM0sfzfzvmefQwsqsRkIXGpiX2m2XnBiO1O7Mr5qaaZm5zcxzLkxezczKLDsiITQIV2qDRgUvtcuvOp5J483U1KIgtTBk5saZJanPLF1mz1uZZSQyQoNwpNbj8I48PjJ4qZ17UdeU1AzU1NKbnzeefLdxqbkDAKabmaZlFnTN7P2li7eRme0zy0ykhAbhDRQYlZqBmprb/PxqxiKjUkv0aObkmUZlZgcA6k/khAbhSs1ESMfWmpqh5qepaVLpMrt+cMIGACbP5J4XJ4QiM9vMzJ5ICg2SHdJhMk4tvU8tyJpadZl1uzR4mRkdADAks+qhGVZm9SOyQoNkS83k6GfQNbWwZGYyaNaOZsaDSAsN8kRqMR79rB40G7TM0vvMTIZm2DizeBB5oUHypRbG6KcfUgsraNZknJkNzYgXWQtNRApE5BMRGelstxeRKSIyT0ReFpFAU6gmWWpb+9QMjn56lVpYQbNJG82sLjM7mumN+tTQrgXmpG3fDzyoqvsC60ilBQmUJEvN9DQpL1IzHWdmdG6mwdFMK7MAyCYlB6lskeOB44CRpHIYrQYKneNdgTf9TEtSFzafmj+4SSKv6vKXrFMPJToFUMjJGf2APE8flK3QhgKdSa3EMhJoAcxPO74ntWek7AOUAqVt27bN+oPKRJLzqUVVaukyG/O0zTSbK+n5zPyUmaoVWsYmp4j0AFaqavr6C1kvYKCqg1S1RFVLWrZsmam4rKmeTy3INQpCHSiISPCt6dAM4xPNbZxZMshkPOBeUqlzFwHfAD8AzxNikzOdvGh+/tVc5tuaamqJrpk5mWbj3MxMhzyvodXvZKfJ6bweApzvvH4cuNLPD6o+WKn5Q00Lr5iW2b//aV5mJtJmm5CZqhWaF6F1AKaSWsBgCNDIzw+qvpjuU7vt2bHm1ygwWFM7Qc7Wf/V9Qq8+4iZjAwBGZWZwQRNTMlO1QqtXZ4GqTgImOa8XAF3q8/5KDa4vyO1Tu3T4q1w7diQiPQJbzNhNPaSqPPbGhwCBLmZ87kVdUVUGPzoBhEAXM+56Wgl3vPYnRj35Nu8Nm0rDokL+9MxVnPjbYwIpD0LIZ5YWmmE6aDbI0Iz15T8Edu24YHTl9K82fssX333Nz3fOuKJ7TqRL7ZoxI3n4lGCldvtvTwLgsTc+RBX6nBqc1M67+EiAlNQIVmqHn9qZw0/tHMi1q5MuM2NBswmMM1uyaTVXTB0c2PXjgtGpTw0Qri59mi+++zqwMpKcT81k8K0JqsssSUGz6ckZg85n5spsS1VFYGXEBaNC26tJS5oUNDIuNVMhHSbyqZmc+xkkqiHMzUxgPrN0mT3WpXdg5cQFo0Jr2KCAx7pcblxqiV33c3w8pebWzMJYNzPJMtt7p90DKysuGM+20XrH4q1Su2raU8xJYPMzaamH/MR4M9OgzEwGzVqZ1Uwo6YNcqTUtbMwfEtqnlrQZBX6Q3sxMosxM5TOzMqud0PKhpdfUkjpQYKr5GYc+terNTLtuZm7EUWYi0k1E5orIfBG5sYbjD4rIDOfxpYisTztWmXZsRKayQk3wuFVqhY0TKTVTixlvs5pUBKUWZjMzSckZYyqzAmAgcApwAHCBiByQfo6qXqeqh6jqIcC/gdfSDv/oHlPV0zOVZzQOrSZa71jMY4f15oppg7m69GkeKfm9sTi1h7r1oPu+Nvg2SIwHzYbUzDQVZ2ZCZou/XcflA4b4dbkupDLzLAAQkZeAnsDntZx/AXB7roWF/43nJ6mFMfppMktHvsWpVe8zS2qmWZNxZhGsmbUQkdK0R59qx9sAS9O2lzn7tkNE9gLaAxPSdjd2rvuRiJyR8W5MzrPKNEfs601r9fSJ9+vxb9+lc9Yvq/Ncr4SVT+2JkcnKp1Yb6XMzk5acMch8ZtVZvHGVdp9wr5749t06//sVGc8nYnM5gXOAwWnblwD/ruXcv1Q/BrR2njuQyvizd13lRaKG5hLWQIHJODVTwbdh1tS0WtCskWZmAoNmF0e7ZpYty0glgHXZA1hey7nnAy+m71DV5c7zAlLzyDvVVVikhAY2Ts0vTK77mY4mOGjWZJzZ4k2ruTL+MgOYBuzrLKpUREpa241Wisj+QDHwYdq+YhFp5LxuAfyK2vvegAgKDWycml9sjVMzFHzrysz0aGYSQzMSIjNUtQK4GniT1CJLr6jqbBG5S0TSRy0vAF5ymrEuHYFSEZkJTATuU9U6hRapPrTqhNWnlrQkkVvzqQXYp2Z8QZOQMs2ujlifWXWIWB+a6Ucka2guNvjWH7bO/QyoT01DamaarpklKTQjqURaaBCu1JIW0hHENClXZsayZqSFZiQpOaMrswqttDLzQOSFBnb00y/8niZVXWY2OWNupMvs0cMuszLzQCyEBrb56RfbTJPyUFNLl5mxZmbCkzNamXknNkIDKzW/8NqnZryZmdDQDNtn5j+xEhpYqflFrn1qGkbQbAJTACUkaDZyxE5okPzgW5NZOuoTfOvKLOmjmTZoNr7EUmgQbvBt0KOft11yojGpZZv5trrMTAbNmhzNtEGz8Sa2QoNwRz9NSC09WZ1UAAANwklEQVSU5mcNNbV0mZnuM7P5zCz1IdZCg/wI6TC+8EpaTc140GxIKYBs0GwyiL3QIH8GCkzmU/v9hY/xxrDpPPC3EXbdTI9YmZkjEUKD/JCaqdRDfa4+gaoq5eEHxjBuzGece1FXu25mjpiUWZVWBnbtuCC6zeT2YOlwUEv99OM5NG3YIrAylv+wjiumPsmmyrJA03kDbCwv59LhrzLjmxWBpvMGqKyq4s7/vMXIKXO44rSugabzBti0qYyyzVsoLGzAz3YO7gcPdt1MP9hcuYHhS2/m/PaPTFfVEi/XKikp0dLS0qzOFRHP5fmJ0RpapW7h1aU3sHHL6sDKsH1q/tCkSSOa79I0UTJLatDs5soNDFv6F1aVzQ+sjLhgVGjNitqwqWJNYqWWtODboElq0GwYMltTtpDubXJeWyQxGBVawwY7cMYe91qpecB08G1QJDVoNkyZdWgabDdEHDA+KNB6xwN/ktqSflZqOWA6+NZvwgqaTVKcmZVZzYQyyrlVapVrbU0tR0wH3/pFmEGzSQnNsDKrndDCNrapqSVYaiaTREZdanbdTO9YmdVNqHFo+SC1pI1+5kqYQbO2mZk/ZBSaiOwpIhNFZI6IzBaRa539zUVknIjMc56Lc7mBfJCaqZpaVPvUbNCsd6zMsiObGloF0E9VOwJHAFeJyAHAjcB4Vd0XGO9s50Q+SM1ETS2KfWo2zsw7VmbZk1FoqrpCVT92Xm8gtbZeG6An8Jxz2nPAGV5uJB+klm99ajY0wztWZvWjXn1oItKO1FLsU4DdVHUFpKQH7FrLe/qISKmIlK5atarO6+eD1PKlTy3J+cyszKJL1kITkabAq8AfVfX7bN+nqoNUtURVS1q2bJnx/HyQWtKDb63MvGNllhtZCU1EGpKS2fOq+pqz+1sRaeUcbwWs9OumrNS8U71PLejUQy5JD5o1sW6mlVnuZDPKKcBTwBxVHZB2aATQy3ndCxju541ZqXmneuqhoKWW9KBZE+tmWpl5I5sa2q+AS4DjRGSG8+gO3AecKCLzgBOdbV+xUvOOKamZDJoNIzmjlVnuiEg3EZkrIvNFZLtoCBH5nYisSvNL77RjvZzQsHki0qv6e7dDVY09OnfurLnw9abPdODcHvrsV710Q/mqnK6RfVlr9fSJ9+vxb9+ln69fFmhZG8rK9OxXXtB9Hv6njvpybqBlVVRW6q3PjNFOfQfooFEf+nrtIe/M0E59B+gfHhmmZeVbfL12dd5bskg7DvyXnvzfZ3XVpk2BlrV44yrtPuFePfHtu3X+9ysCLevHiu/1hYVX6L+/6KZfbcj9/wcoVYO/00zlAQXAV0AHoAiYCRxQ7ZzfAY/U8N7mwALnudh5XVxXebHIWBtaTc3gEnlxHf1MatCsyXUzk1ozc+gCzFfVBapaDrxEKuQrG04GxqnqWlVdB4wDutX1huC+fT7jSu31ZTfx6tIbOGvP/oFlvnWldsXUJ7l62lM8cthldAwo860rtUuHv8o1Y0by8Ck9OGWfYDLfulITER5740MAT5lvTQ8AmJSZqaXmoiizZUvW0O+q//h1uTbA0vTLA4fXcN5ZInI08CVwnaoureW9df4QY1FDcwmjphbGup9xSD2U5NCMfJZZDrRw40ydR59qx2taiKJ63v83gHaq+kvgbX4K2M/mvdWOxqAPrTph9anNMdinNnpedPvUhkyeaazP7P0li22fWT0gen1oXYE307ZvAm6q4/wC4Dvn9QXAE2nHngAuqKu8WNXQXMLqU7tq2lOJqqnl0qcWVtpsm5wxtkwD9hWR9iJSBJxPKuRrK248q8PppKZXArwJnCQixU7yi5OcfbUSS6FBeM3PpIZ0ZBN8axcB9k6eyQxVrQCuJiWiOcArqjpbRO4SkdOd065xMvnMBK4hNeqJqq4F7iYlxWnAXc6+WjG6jF19lsfKluU/zOL1ZTfRpHCXQAcKUmXl7xJ5dkET75iQmR/Lytll7EIkH4Jvrx0bbpYO0+tmWplZciX2QoP8kZrJ5ued//cWM79azuDRU2w+M49YmZkjEUKD/JCayT614R/M5tL+L/PoGx/YZqYHrMzMYrYPrdPuOu3jhYjsEFgZ+dCnFmTwLaRCeRZ/u44VazfQsLABB+/dmoYFBYGVZ2XmHdUK9PvbKWh2j+1DM0bVenTdFaj+GFgRtqbmHRGh3e7N6XrAXpTst6eVWQ4Yl9l3/eDHIYGVERfMCq1gDyj/0EotR0xKzQQ2n5l3tsps8xhkp78EVk5cMCu0Bs2Qne+zUvNAUqRm85l5p7rMpMllgZUVF4wPCsgOZ1qpeSTuUgsjn9mWqgorszwglFFOKzXvmFxNyk/supnesTKrndDCNqzUvGMyTs0PrMy8Y2VWN6HGoVmpeScuzU8rM+9YmWUm9MBaKzXvRF1q6TKzWTNyw8osO0IXGlip+UFUpZY+nclmzcgNK7PsiYTQIESpLelnRmoJSz2UDTYFkHeszOpHZIQGIUmtcq2ZmtphvfNKalZm3rEyqz+REhokvPl5WO+8aH5Wn85kZVZ/rMxyI3JCg4RLLeFxajafmXeszHInkkKDPOhTS2Ccml030ztWZt6IrNAg4X1qCRv9tOtmesfKzDuRFhrY5qcfBC01u26md6zM/CHyQgMrNT8Iqk/N5jPzjpWZf8RCaGCl5gd+96nZdTO9Y2XmL7ERGlip+YFfzc+w8plZmVnqIlZCAys1P/AqNZuc0TtWZsEQO6GBlZof5Cq1MJIzWplZsiWWQgNXavcmXmpzIjRQEFYKICszS7bEVmgAssP/JL+mNu0pY1Kra6DAdJyZHQCw5EKshQbJb342LWxsVGo1NT9NhmbYoFmLFzwJTUS6ichcEZkvIjf6dVP1vo88kNofDPepDZ87h03l5by7eJENmvVIvssskydE5HoR+VxEPhWR8SKyV9qxShGZ4TxGZCwr15XTRaQA+BI4EVgGTAMuUNXPa3tPfVZkzgX9cRj63Y1Q1BUpfsyu0J4D7grt01cs37pvv11a2BRAOWJaZn6sZO7nyunZeEJEfgNMUdUfROQK4FhVPc85tlFVm2Z774XZnlgDXYD5qrrAKfgloCdQq9CCRnY4EwD97q9QPhMaBfdFdWtqry+7ia82vs/BxT0DLCtVU7ti6pMMWzaNmwIUWtOiIp4742w+XrGcOatX0rBBAT3370jxDsH9cQAYu3yGEZkBLN30iRGZAVAxDzZPjFXNbNnc5fT7ze1+XS6jJ1R1Ytr5HwEX51qYlxra2UA3Ve3tbF8CHK6qV1c7rw/Qx9k8EJiV682GQAsguDal/9j7DZY43O9eqtrSywVEZCypf2s2NAY2p20PUtVBadfKyhNp5z8CfKOqf3O2K4AZQAVwn6q+XtfNeKmhSQ37trOj848b5NxcqdfqsEns/QaLvd9ooqrdfLxcVp4AEJGLgRLgmLTdbVV1uYh0ACaIyGeq+lVthXkZFFgG7Jm2vQewvJZzLRZLfpKVJ0TkBOAW4HRVLXP3q+py53kBMAnoVFdhXoQ2DdhXRNqLSBFwPpBxFMJiseQVGT0hIp2AJ0jJbGXa/mIRaeS8bgH8igx99Dk3OVW1QkSuBt4ECoCnVXV2hrcNynA8atj7DRZ7vwmnNk+IyF1AqaqOAB4AmgJDRARgiaqeDnQEnhCRKlKVr/vqiqIAD4MCFovFEjViP1PAYrFYXKzQLBZLYjAitKhMkaoNEdlTRCaKyBwRmS0i1zr7m4vIOBGZ5zwXh32v6YhIgYh8IiIjne32IjLFud+XnU7YSCAizURkqIh84XzOXaP8+YrIdc53YZaIvCgijaP8+VpSBC40Z+rDQOAU4ADgAhE5IOhy60kF0E9VOwJHAFc593gjMF5V9wXGO9tR4lpgTtr2/cCDzv2uA6IUmv4QMFZVfw4cTOq+I/n5ikgb4BqgRFUPJNWZfT7R/nwtmKmhbZ36oKrlgDv1ITKo6gpV/dh5vYHUj60Nqft8zjntOeCMcO5we0RkD+BUYLCzLcBxwFDnlMjcr4j8DDgaeApAVctVdT0R/nxJRQDsICKFwI7ACiL6+Vp+woTQ2gBL07aXOfsiiYi0IxW8NwXYTVVXQEp6wK7h3dl2/Av4M1DlbO8CrFfVCmc7Sp9zB2AV8IzTRB4sIk2I6Oerql8D/YElpET2HTCd6H6+FgcTQst66kPYiEhT4FXgj6r6fdj3Uxsi0gNYqarT03fXcGpUPudC4FDgMVXtBGwiIs3LmnD68noC7YHWQBNSXSbVicrna3EwIbRYTJESkYakZPa8qr7m7P5WRFo5x1sBK2t7v2F+BZwuIotINeGPI1Vja+Y0kSBan/MyYJmqTnG2h5ISXFQ/3xOAhaq6SlW3AK8BRxLdz9fiYEJokZ8i5fQ/PQXMUdUBaYdGAL2c172A4abvrSZU9SZV3UNV25H6PCeo6kXAROBs57Qo3e83wFIR2d/ZdTypKSyR/HxJNTWPEJEdne+Ge7+R/HwtP2FkpoCIdCdVg3CnPvw98ELrgYj8GngX+Iyf+qRuJtWP9grQltSX/BxVXRvKTdaCiBwL3KCqPZyMBC8BzYFPgIvTJ/qGiYgcQmoAowhYAFxK6g9qJD9fEbkTOI/UCPgnQG9SfWaR/HwtKezUJ4vFkhjsTAGLxZIYrNAsFktisEKzWCyJwQrNYrEkBis0i8WSGKzQLBZLYrBCs1gsieH/AUvPBPtl/WHVAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f0725851198>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n = 100\n",
    "x = np.tile( np.linspace( -1.0, 1.0, n ), ( n, 1 ) )\n",
    "y = x.T\n",
    "\n",
    "pp.contour( np.abs( x ) + np.abs( y ) )\n",
    "# pp.contour( np.sqrt( x ** 2.0 + ( y / 1.5 ) ** 2.0 ) )\n",
    "pp.axis( 'square' )\n",
    "pp.colorbar()\n",
    "pt = [ -0.5, -0.5 ]\n",
    "pp.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercício\n",
    "\n",
    "Prove que se $f : \\mathbb R^n \\to \\mathbb R$ é diferenciável em $\\def\\vect{\\boldsymbol}\\vect x$, então o gradiente $\\nabla f( \\vect x )$ é ortogonal à superfície de nível\n",
    "$$\\large\n",
    "\\Gamma := \\bigl\\{ \\bigl( \\vect y, f( \\vect y ) \\bigr) \\in \\mathbb R \\times \\mathbb R^n : f( \\vect y ) = f( \\vect x ) \\bigr\\}\n",
    "$$\n",
    "no ponto $\\bigl( \\vect x, f( \\vect x ) \\bigr)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dica: note que se parametrizarmos a superfície de nível por\n",
    "$$\\large\n",
    "\\Gamma = \\{ \\vect\\gamma( \\vect t ) : \\gamma \\in \\mathbb R^{n - 1} \\}\n",
    "$$\n",
    "onde\n",
    "$\\vect \\gamma : \\mathbb R^{n - 1} \\to \\mathbb R\\times \\mathbb R^n$\n",
    "temos\n",
    "$$\\large\n",
    "f( \\vect \\gamma( \\vect t ) ) = c.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Subdiferencial e subgradientes\n",
    "\n",
    "Seja $f : \\mathbb R^n \\to \\mathbb R$ convexa. Definimos o subdiferencial de $f$ em $\\vect x \\in \\mathbb R^n$ como\n",
    "$$\\large\n",
    "\\partial f( \\vect x ) := \\{ \\vect v \\in \\mathbb R^n : f( \\vect y ) \\ge f( \\vect x ) + \\vect v^T( \\vect y - \\vect x ) \\quad\\forall \\vect y \\in \\mathbb R^n \\}.\n",
    "$$\n",
    "\n",
    "Se $\\vect v \\in \\partial f( \\vect x )$, dizemos que $\\vect v$ é um subgradiente de $f$ em $\\vect x$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Propriedade: Seja $f : \\mathbb R^n \\to \\mathbb R$ convexa, então, para todo $\\vect x \\in \\mathbb R^n$, $\\partial f( \\vect x ) \\neq \\emptyset$. Além disso, se $f$ é diferenciável em $\\vect x$, então $\\partial f( \\vect x ) = \\{ \\nabla f( \\vect x ) \\}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nós aqui iremos comparar o valor de\n",
    "$$\\large\n",
    "\\| \\vect x - \\vect y \\|_2^2\n",
    "$$\n",
    "com o valor de\n",
    "$$\\large\n",
    "\\| \\vect x - \\lambda \\vect v - \\vect y \\|_2^2\n",
    "$$\n",
    "onde $\\lambda \\ge 0$ e $\\vect v \\in \\partial f( \\vect x )$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\large\n",
    "\\begin{split}\n",
    "\\| \\vect x - \\lambda v - \\vect y \\|_2^2 &{} = \\| ( \\vect x - \\vect y ) - \\lambda v \\|_2^2 = \\| \\vect x - \\vect y \\|_2^2 - 2 \\lambda\\vect v^T( \\vect x - \\vect y ) + \\lambda^2\\| \\vect v \\|_2^2\\\\\n",
    " &{} = \\| \\vect x - \\vect y \\|_2^2 + 2 \\lambda\\vect v^T( \\vect y - \\vect x ) + \\lambda^2\\| \\vect v \\|_2^2\\\\\n",
    " &{} \\le \\| \\vect x - \\vect y \\|_2^2 + 2 \\lambda\\bigl( f( \\vect y ) - f( \\vect x ) \\bigr) + \\lambda^2\\| \\vect v \\|_2^2\n",
    "\\end{split}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Método do subgradiente\n",
    "\n",
    "1. $k \\leftarrow 0$\n",
    "\n",
    "2. Para $k \\in \\{ 0, 1, 2, \\dots \\}$\n",
    "\n",
    "    1.$$\\large\n",
    "        \\vect x^{( k + 1 )} = \\vect x^{( k )} - \\lambda_k\\vect v_k\n",
    "    $$\n",
    "    onde $\\vect v_k \\in \\partial f( \\vect x_k )$\n",
    "    \n",
    "    2. $k \\leftarrow k + 1$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Métodos de bundle são alternativas mais sofisticadas e que podem rodar bem em problemas médios/pequenos."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "O algoritmo converge se os subgradientes forem limitados e\n",
    "$$\\large\n",
    "\\sum_{k = 0}^\\infty \\lambda_k = \\infty, \\quad \\lambda_k \\to 0^+.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercício\n",
    "\n",
    "Qual o subdiferencial de $| x |$ em todos os pontos do domínio? Prove que é\n",
    "$$\\large\n",
    "\\partial | x | = \\begin{cases}\n",
    "                    \\{ -1 \\}  & \\text{se } x < 0\\\\\n",
    "                    [ -1, 1 ] & \\text{se } x = 0\\\\\n",
    "                    \\{ 1 \\}   & \\text{se } x > 0\n",
    "                 \\end{cases}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercício\n",
    "\n",
    "Prove que o máximo dentre funções convexas é uma função convexa.\n",
    "\n",
    "## Exercício\n",
    "\n",
    "Qual é o subdiferencial do máximo entre duas funções convexas?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\large\n",
    "\\| A\\vect x - \\vect b \\|_1\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}