mogan-2d.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW0AAAEGCAYAAACuMsS7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnXmYnFWV/z/f7qyEbHQnEbJDIhAQAgTMCEhEQRAEITMICsKIgjMwoziguPxQQQYccZsRwSDIqoisEcMmJiwqS1gCBAwmIUuTkJBAFrYk3XV+f9xbyZtKVVdVd63d5/M89fT73nvf+56urvfbp849916ZGY7jOE590FBtAxzHcZzCcdF2HMepI1y0Hcdx6ggXbcdxnDrCRdtxHKeOcNF2HMepI1y0nU6jwF8lfbTC9z1f0gpJb0nav5L3dkqHpN9JOr3adtQLFRdtSbMkbZC0XtJaSQsl3SBpv0rb0hWQNEnSE5LekbRA0skFXHOipOckvS3pNUnfytHud5JM0kF5ujwBaDWzBxPXjpP0p3iPFkn/lcemayVtigKcfv17O+1HAP8NfMTMtjezJ/PYWBCSpsTf+S1J6yStlvQXSV+V1Lud68bE60aUwo5uxneA/5bUt9qG1APV8rQvMrP+ZjYQ+AiwGHhM0nFVsqdgJPWstg1pJA0E7gFuAwYDXwKulPRP7VxzCvAT4KvAQGA8MD1Lu+OBpgJN+QpwVeLaRuAPwEvAEOAY4OuSPp2nn+uiAKdfv2in7RggZWYvFmjjNrTzt2yL9x8ADAe+C3wemCWpV0fv19Uo1bNgZn8H5gMnlaK/Lo+ZVfQFzAK+naX8aqAFUI7rTgTmAOuA5cAvgX6J+u2By4CFwHpgLnBQrOsJfBOYF+sWAFNj3bXArzLutQg4OR6fRvhAnRftmxvL/zve663Y31cy+hgD/D7augb4C0EEfwDcldH2UGBt8vcp8L38V2BJ8j0DbgB+naN9A/Aq8KU8/TbF321nwNLvY462w2KbHRNlHwHeAbZPlF0EzGynn23+Du20/TTwbrzvW8CChN3Xx/f8NeA6YIeMv+sFwEzgbeDELH1PIXxryCwfD7wHfDGHTWujPW9Hm/5fLB8N3AWsApYCPwX6tvO7tfc53g74WexnFXAnMCrj2foR4Z94+nN+bKzbIdo/MeN+DwEXxOMehOfk5cRndr+Mv9FNwK+BN4ArYvlRwIvx976b4BTMyvg8XR3tfh24BRiWYcd3gekd0ZTu9qr8DXOL9mHxQ79bjuuOBPYgCM+4+CG5JFH/O+ARYCyg+JCNi3U/iO33inUjgL0SH8R8ot0aP4h9ge1i+cnATrG/Qwki8vFYt1186C4neLM9gH8C+gO7ABvZWuR+A1yZOF+T5zUqtvspcGeG7ecAT+d4D3eL7/G5wN+BlfEhG5fR7rfAOfE4n2gfCbyRUfYV4NmMsuMy22XUXwu8GcXgZeCHJEQ/S/spZIgrcC/Bwx8cX38E/pjxd10K7BP/btuIZ7Z+E3V/AW7OUTcmvlcjEmU9gBeIDgbBa38SuLyd36u9z/EvgcdjP/2AXxEcmcbEs7UKOJDwnJwTPy/pz+wtwE8T99qZ8NkeHc//O/a/M9AInB77G5z4G20k/NNsJHzOdwE2ELzkHsBHCY7VrHiN4u/zK8KzsB1BwB/M+L2nAi2V1qN6fFX+hrlFe/f4oT+wwH7OBp6Ix0PjtXtkaSeCB3BUjn6uJb9ovwv0zmPPrcD/xOMTCN5ejxxt/wScH48Hx/73a6//HP1cTQgpJMv+FZifo/1B8X16PorCdsAvCGGMHrHNp+KD2xDP84n2Z4BFGWX/D3goo+wj5BDDWL8fwWtvIPxzfgL4bTvtpyT7I/wDNWB8omxXEt8C4t/1gjzv6ZRcdhIE9YEcdWPYVrQ/RBC05DfCj8e/9zbfKPN8jhvidYclyrYniOg/JZ6tyxP1/WJ/e8fzIwmebs94fmH69yE8J+uBD2fc9/nEs3At8OeM+m8DD2eU3cAW0Z5E+NbVO1HflOW9Ogx4p9hnoDu+ail7JD2AszpbpaTDJD0i6XVJ6wje85BYPSb+fDnLpUMIH95sdYWy3Mw2ZNjzn5Kel/SmpDXAJzPsWWhmrTn6+yXBi4Hgsb9kZk91wK71BO8lySCCp5OrPcDPzOwVM3uH8HV4N+D9knYA/hf4gpmlCrThTWBAJ+3CzJ4ysxVmljKzuQQv8Z/bG/zLYGT8+UqibEFGHQTh7igjyPH5bMemlWb2doZNfdjyWUkyJv7M9TnuQ/gGB4CZvUX4tpT8/ZYn6tP37R9/3k8Q+U9KEvA54JpY10z4J/AHSWvSL4LXnRxcXZRh13DCmFSS5PlYoDewItHnAkKoZlSi3QDCtywnD7Uk2p8mxFvnZVbEwZ87gZsJoYEBwNcJ3gFs+SCNz9Lv64Q4Y7Y6CF54v8S9ehA8niRbCZikAwn/NM4Ems1sEOFredKesXFALht3Av0lHUIQ76uSlRkZFNle6Q/7HMJX/ST7xPJszGNLLDgTI4SPdgJmSlolaVWsu1vSpTn6fAYYLOl9ibI5hH8C/RJl7dmVjfR7rnZbbWFp/DkmUbZzRl2y36KQNI7wbeDPOZpk63cpMFTSdhk2vUcIO2SyKP7M9TneQBDBtE3bEz6rS7O03wYzayPE/E8jhPQGAnfE6lWE5+RjZjYo8epnZsm/febv+Sohbp8kKcaLY787ZPTb18z+mmi3J+Gz5OSj0q49GeERgpfwPUJsbWqOa7YH2oDPxPMJhJjsokSb38e+xxAe9HFsiQX+kBBb3DPWDQc+EOtOI3graY/gB8AmMgYiM+w5kiD24wn/+I4ifDCvjfX9CA/g/xIejEZgMtA/0celwNPxuoEdfC8HER7m84BehHjiW8SvyzmuuRx4Lr7vvYH/i+9NYzwfkfEyQjx6UDt9PpZ+v+J5IyHk8jPCOMBEYAVZBv4S15yYvkd8X/8K3NZO+ylsG9O+jzDoN4gQdvoDMCNRvyhpZyH9xvfko/E9e5wcYbL4e7YRUhDTZT0IA4lXEEJRO8U+ftHO/dv7HE8D/hb72Y7wje05to5pfzujv63CW8D7CZ/v+zPtAC4mDEyOTzx3Hwd2iufXsm0oMR3TPiH+3aewdUy7AXg4fs6aYtmQzM8C8Chwerl0pyu9Kn/D8MHaQPgKvY7wdfYm4IA8151ByN54izD6fwFbi3Z/wsDc4tgmOereK7afH+vmA8fHup6EQZI3Y/9nkSV7JMOWhvggpgfOfg3cSBTt2GZnghezkjAY9AhbZzLsQvBaft3J93N/Qvz3XcJX55Mz6ucC30yc9ybEsd8geFfTgbHt9N9uTDu2OZFEtkAsGwc8SIhnLgPOzai/h60HX2dFm96On4kfAwPauecUthXtIfHv8Brhn8QNhG9C6frNf9c8/Vr8nKyLNv0N+BrQJ8+134z3XgN8K5aNJfzzWBU/X/9LHBjM0Ud7n+N+BPF7NfG3G5PxHrYr2rHsL7F8UkZ5D0Iq6ItsydK6gxh7JkeGDyGl8yW2ZI9cDtyXqN8hli0iPPcLM/72u8a/V86smlp4EUJJK4EXctQr/n3nE/6Z7puoOxX4R3ydmijfjzBuMD9emzV7LvlSvNCpMDF0sAI43Lb+mlh3xPjoXwlClSt84HQTJP0WWG9mZxTR/kEz+1V5Lesckj5M+Md0vZntmaX+E8B/AJ8APkgYO/pgHCuaTRiUNeApQuLBm5KeAL5M+LY6A/hfM7unPTt6lPB3cgokitxXCAOQdS3YABb+8+ec0ON0bSR9khDeWE8IFU4lhFUKwszqYlKNmT0saUw7TY4lCLoRJgsOkrQj4dvbA2b2BoCkB4AjJM0ifJv8Wyy/npC95aJdS0gaSvh6uBL4lyqb4zil4BBCiLAPYbLXl8xsZnVNqgrD2XpQuCWWtVfekqW8XSoq2s3NzTZmzJhK3rLm2G+/zUusjAVmT5o0qYrWOE7nSXymIcy3uHrSpElXd7bfp556apWZZUuNLJhhwyfaxg3r87Zbs3rhXEJWT5ppZjatyNtly3SyDpS3S0VFu42BjN7za5W8peM4dcpTT306M/+7aDZuWM+hR1+St93t1336PTPrrAfVwtY58yMIg/AthBBJsnxWLB+RpX271FKetuM4Tj0zHfhcXKp4MrDWzJYTUlEPlzRY0mDgcEJ2zXJgvaTJiclOd+W7ScGedpwoMht41cyOlnQtIZa1NjY5zcyeLeIXdBzHqRtilssUoFlSC2FJ2Z4AZnYlIfvjE4T0vXcIS0pgZm9Iuoiw7gzAhelBSeDfCKmUfQkDkO0OQkJx4ZEvE3Ixk1OWzzOzW4vow3Ecpy7Jl+USs0bOylF3DVuWDEiWzyZM+iuYgkQ7Lux+FGHG1FeLuYHjOMWzXd9Gjjn8fQxt6o0Knchfh5jBytUbmH7/a7zzblu1zakLCvW0f0qYEdY/o/xiSRcQZr6dbxmLKjmO0zGOOfx97Pb+4fTp2x91YdU2M3ZoCtkdN9/1apWtqQ/yDkRKOpqwUlnmKnTfIKwOtz9hmurXc1x/hqTZkmZveC/nIm+O4yQY2tS7yws2gCT69O3P0KZCF3N0CskeORA4RtIiwip7h0q60cyWW2ADIbH+gGwXm9k0M5tkZpN698lcwdPpTqxe2cS853dj9cpCdzHrvkh0ecFOI6lLh4BKTV7RNrNvmNkIMxtDWBjoz2Z2cpyemZ6S/SnCSnGOk5XVK5t49P5DmPv0njx6/yEu3HXAruOG8MmjDuETRxzIGV/4DOvWrc1/UQ6mHDyRN94oZilyJxedydO+SdLzhBWqmoHvl8YkpyuyasUQ2toagAbaUg2sWtGpiW5OBejTpy9/+ONDzLj3LwwaNIgbb6jp9Zy6DUXNiDSzWYSZPJjZoWWwx+miNA97ncbGFG0paGxI0Tzs9Wqb5BTBxH32Z97f524+v2ra/3HPH+9k48aNHHb4UXz5nPMB+LczT2b58mVs2PAep552JieedGq1TO6y+IJRTkVoGrqagw5/iFUrhtA87HWahvpX5VLzzNM9ePyxHnxwciv77Jtrp7viaWtr429/fZh/OeFkAB55ZCaLFy3ktjv/hJlx5hc/yxNP/JUDDvgQl/zg/xg0aDDvvfcux3/qY3z8iE8yePAOJbPFcdF2KkjT0NUu1mXimad78LmT+7NxI/TqBdffuL7Twv3ee+/yyaMO4dWWJey5594ceNAUAP7yyEwefWQmxxwdzt95520Wv7KQAw74ENdfO40H7v8jAMuXvcqiRQuqKtqbejXy2ujM7UrrGxdtx+kCPP5YDzZuhFRKbNpkPP5Yj06LdjqmvX7dOs74wknceMOvOPW0M4N3/W9f4aTPnJZhw6P89S8Pcctt99K373Z89qRj2LjBp26UGl8wynG6AB+c3EqvXtDYaPTsGc5LRf8BA/j2dy7h6qsuZ9OmTRz04UO59fc38fbbbwHw2mvLWL3qddavX8eAgYPo23c7Fix4mWefmV0yG5wtuKftOF2AffZt5fob15clpg2wxx57sdvue/LHu2/nU8d9mgXzX+aEqUcAsF2/flz24ys5+MMf5be/uZajjzyYsTuPY+I+vlZ8OXDRdpwuwj77llas57ywZKvzab/6zebj0/71TE771zO3uebqX9+Sta9Zj/gCoKXCwyOO4zh1hIu24zhOHeGi7TiOU0e4aDuO4xSApCMkzZM0X9L5WepHS3pQ0nOSZsV9CNJ1P5D0Qnx9OlF+U+zzBUnXSOqZzw4XbcdxnDzE7RYvB44EJgAnSZqQ0ewy4Hoz2wu4ELgkXnsUsC8wEfggcJ6k9JKnNxGWuP4AYcuxL+SzxUXbcRwnPwcA881soZltJCxTfWxGmwmEDWEAZibqJwAPmVmrmb0NzAGOADCzGXGJawOeYOvd2bPiou04TlbSS7OmXy0tS3K2ffyxR/ni6e1uoVjvDAeWJs5bYlmSOcDUeHwc0F9SUyw/UtJ2kpqBjwAjkxfGsMgpwL35DPE8bcdxspKexl7PNPQy+o0pKHe9WVJyCuc0M5uWOM+2TYNlnJ8L/FzSacDDwKtAq5ndL2l/4K/A68DfgEyjfgE8bGaP5DPUPW3HcQqmpWUJJ51wFMd+8iMc+8mP8PRTT2zT5rk5T3PM0VNYsmQR77zzNud/7T84/tiPcszRU/jTAzOqYHVBrErvsBVf0zLqW9jaOx4BLEs2MLNlZna8me0DfCuWrY0/LzaziWZ2GOEfwD/S10n6DjCEAjdNd0/bqRlWr2zKu3RrIW26K32efpK+jz3Ku5MP4r199+90f+lV/gBGjhzFL668gaamZq694TZ69+7DolcWcM6Xv8gd0/+8+Zqnn3qCC793Plf+8kZ2Gj6CH/3wIv7pQwdz6f/8H+vWrWXqpw7jQwcewnbb9eu0fRXmSWC8pLEED/pE4DPJBjH08YaZpQh76F4TyxuBQWa2WtJewF7A/bHuC8DHgY/G6/Liou3UBOntyNraGmhsTHHQ4Q9tI8qFtOmu9Hn6SXY6+Ti0cSPWqxfLbryj08KdLTzSuqmV7333a7z04gs0NDay6JUFm+sWLHiZb3/rHH593a0MG7YjAI8+MpMHH7yXX111OQAbNrzHsmUtjBu3a6dsqzRm1irpbOA+oBG4xszmSroQmG1m04EpwCWSjBAeOSte3hN4JO75uQ442czS4ZErgcXA32L97WZ2YXu2uGg7VWf1yiZemjMhsR1Z2J4sU5C33rIse5vuSt/HHkUbN6JUG2zaSN/HHi2Jt53Jr6+5gubmofxhxsOkUin23H2nzXVDhwxjw4b3eHHu85tF24Cf/+Jadt55fMltqTRmNgOYkVF2QeL4VuDWLNe9R8ggydZn0RrsMW2nqqS955XLhhFCfamc25GltyxDudt0V96dfBDWqxfW2Ij17MW7kw8qy33Wr1/HkKHDaGho4M47fkdbW9vmuv4DBnLV1Tfz48u+z+OPPQrAwQd/hBuuu4qQ0QZz5z5XFru6E+5pO1Ul6T1DiqE7rWD3vV/M6kH7lmW5eW/f/Vl24x0ljWln47Mnf56z//007plxF5MnH7RNbLp5yFB+edVvOP3zJ3DJD/6Xs/7jXL5/0bc4+siDMYzhw0dx1dW/LYtt3QUXbaeqNA97nYaGFKkUNDSk2H3vFwGY9/xuWYXZtyzLzXv77l9Ssc5cmhVgzNhduPueLVlp534tRAc+OPkgPhi9+52Gj+Ce+/66uc33L/5xyWxyXLSdmkCbX2vfHMDzT+7jg42OkwOPaTtVZdWKIaRSQbBTJpYtGbFlsLGtgVUrhlTbRMepKVy0naqSObg4cPCbBK/bANGzl28M6zhJCg6PxATx2cCrZnZ0TDK/GdgBeBo4JS6k4jgFkzm4GDxrIz0wuWlj7ypbWB3MwMyIubtdGjPDMieEOzkpJqb9ZeAlIL2k4A+An5jZzZKuBE4HriixfU6dUszMxczBxcbGFG0punVa38rVG9ihaT19+vbv0sJtZrz37npWri7PN6pevdsYOXp9WfquFgWJdlzM+yjgYuCrCp+iQ9kyjfM64Lu4aDt0buaip/UFpt//GgBDm3rThTUbs/APKv37Ovkp1NP+KfA1oH88bwLWJKZiZlum0OmmdHbmoqf1wTvvtnHzXa9W2wynBsk7ECnpaGClmT2VLM7SNGtUStIZkmZLmr3hvXUdNNOpZVavbGLe87uxemUT4DMXHaecFOJpHwgcI+kTQB9CTPunwCBJPaK3vc0yhWniEofTAAY37+LDDV2MXKGQjoQ4Fs4by7IlI9hpVAs77/pKmS13nPokr2ib2TcIywwiaQpwrpl9VtLvgX8mZJCcCtxVRjudGiVXKKTYEMfCeWN59rFJAKxc9j4AF27HyUJn8rS/ThiUnE+IcV9dGpOceqJUoZBlS9Jb4ynj3HGcJEVNYzezWcCseLyQsNml040pVbbHTqNaoodtm88dx9kWX3vE6TSlyPZIh0Lai2n7rjVONZF0BPAzwiYIvzKzSzPqRxN2qxkCvEHY7KAl1v0PIW26AXgA+LLZlilFkqYDO5vZnvnscNF2aoadd30lZxy73LvW+D8Epz3ijPDLgcMIKc5PSppuZi8mml0GXG9m10k6FLgEOEXShwgJHXvFdo8ChxCjFpKOB94q1BZfe8SpC7Ye8CztQlLpfwhzn96TR+8/ZHPqouMkOACYb2YL43IdNwPHZrSZADwYj2cm6o2QedcL6E3YfmwFgKTtCRv6fr9QQ1y0naqTmeedjXLmfpfzH4LTZRgOLE2cZ5tQOAeYGo+PA/pLajKzvxFEfHl83WdmL8V2FwE/At4p1BAPjzhVpdCwRzmnt6f/IXT39U66In0aYdwOBW1y3ixpduJ8WpxjkqaQCYXnAj+XdBphY99XgVZJ44DdCfNZAB6Q9GHCJr/jzOwcSWMKMRJctJ0qU8yU93JNb/f1ThxglZlNaqe+BRiZON9mQqGZLQOOh81hj6lmtlbSGcBjZvZWrLsHmAysB/aTtIigxUMlzTKzKe0Z6uERp6pUesp7ZigmfQ6w6wf+7oLt5OJJYLyksZJ6AScC05MNJDVLSmvqNwiZJABLgEMk9ZDUkzAI+ZKZXWFmO5nZGOAg4OV8gg3uaTs1wKhdFsWfi8sqmpmhmA/s/4xvbeYUhJm1SjobuI+Q8neNmc2VdCEw28ymA1OASyQZITxyVrz8VsKqqM8TQir3mtkfOmqLi7ZTNTJFdNQuizvUR6FhjcxQzFZbm3VgNUKne2FmM4AZGWUXJI5vJQh05nVtwJl5+l4E5M3RBhdtp4p0dgnXYnO3MwccdxrVwuoVQ3wA0qkrXLSdqtHZrI1iRT/bgOPAwet8ANKpK1y0narR2ayNjoh+ZgaKb7jg1Bsu2k5V6Yxoeqqe0x1x0XbqGveUne6G52k7juPUES7ajuM4dYSHRxzH6bL0aYTdBha09kjd4J624zhOHeGi7TiOU0e4aDuO49QRLtqO4zh1hIu24zhOHeGi7TiOU0e4aDuO49QReUVbUh9JT0iaI2mupO/F8mslvSLp2fiaWH5zHcdxqoOkIyTNkzRf0vlZ6kdLelDSc5JmSRoRyz+S0MlnJb0n6VOxTpIulvSypJck/Wc+OwqZXLMBONTM3opb5Twa9zgDOC8u/O04jtNlkdQIXA4cRtgv8klJ083sxUSzy4Drzew6SYcClwCnmNlMYGLsZwdgPnB/vOY0wt6Tu5lZStLQfLbk9bQt8FY87RlfmbsQO47jdGUOAOab2UIz2wjcDByb0WYC8GA8npmlHuCfgXvM7J14/m/AhWaWAjCzlfkMKSimLalR0rPASuABM3s8Vl0cvwr8RFLvQvpyHMepQ4YDSxPnLbEsyRxgajw+DugvqSmjzYnAbxPnuwCfljRb0j2SxuczpKC1R+IeZxMlDQLukLQnYbfh14BewDTg68CFmdfG7ePPAOjbr7mQ2zmO45SEXo0pxg7YWEjTZkmzE+fTzGxa4lxZrsmMOJwL/FzSaYSNfV8FWjd3IO0IfICwOXCa3sB7ZjZJ0vGEHdwPbs/QohaMMrM1kmYBR5jZZbF4g6RfR4OzXTONIOoMbt7FwyqO49Qiq8xsUjv1LYTYc5oRwLJkAzNbBhwPIGl7YKqZrU00OQG4w8w2ZfR7Wzy+A/h1PkMLyR4ZEj1sJPUFPgb8Pf7XQJKATwEv5OvLcRynTnkSGC9prKRehDDH9GQDSc2S0pr6DYLXnOQktg6NANwJHBqPDwFezmdIIZ72jsB1cfS0AbjFzO6W9GdJQwhfG54FvlRAX47jOHWHmbVKOpsQ2mgErjGzuZIuBGab2XRgCnCJJCOER85KXy9pDMFTfyij60uBmySdA7wFfCGfLXlF28yeA/bJUn5oluaO4zhdEjObAczIKLsgcXwrkDUF2swWse3AJWa2BjiqGDt8RqTjOE4d4aLtOI5TR7hoO47j1BEu2o7jOHWEi7bjOE4d4aLtOI5TR7hoO47j1BFFTWN3HMepJ3o1wvB+2ZYNqV/c03Ycx6kjXLQdx3HqCBdtx3GcOsJF23Ecp45w0XYcx6kjXLQdx3HqCBdtx3GcOsJF23EcpwAkHSFpnqT5ks7PUj9a0oNxs/NZkkYk6kZJul/SS5JejJsiIOmjkp6W9KykRyWNy2eHi7bjOE4e4s5dlwNHAhOAkyRNyGh2GXC9me1F2OT8kkTd9cAPzWx34ABgZSy/AvismU0EfgN8O58tLtqO4zj5OQCYb2YLzWwjcDNwbEabCcCD8Xhmuj6Kew8zewDAzN4ys3diOwMGxOOBZGwWnA0XbcdxnPwMB5YmzlvYdvuwOcDUeHwc0F9SE/B+YI2k2yU9I+mH0XOHsCfkDEktwCmEPSPbxdcecRyny9KongzoObSQps2SZifOp5nZtMR5tgVMLOP8XODnkk4jbOz7KtBK0NmDCXvtLgF+B5wGXA2cA3zCzB6XdB7wY/Js7uui7TiOA6vMbFI79S2E3dTTjCAjlGFmy4DjASRtD0w1s7XRi37GzBbGujuByZKmA3ub2eOxi98B9+Yz1MMjjuM4+XkSGC9prKRewInA9GQDSc2S0pr6DeCaxLWDJQ2J54cCLwJvAgMlvT+WHwa8lM8Q97Qdx3HyYGatks4G7gMagWvMbK6kC4HZZjYdmAJcIskI4ZGz4rVtks4FHpQk4CngqtjnF4HbJKUIIv75fLa4aDuO4xSAmc0AZmSUXZA4vhW4Nce1DwB7ZSm/A7ijGDvyhkck9ZH0hKQ5kuZK+l4sHyvpcUn/kPS7+JXBcRzHKSOFxLQ3AIea2d7AROAISZOBHwA/MbPxBLf+9PKZ6TiO40ABom2Bt+Jpz/gyQjA9/VXgOuBTZbHQcRzH2UxB2SOSGiU9S5h6+QCwAFhjZq2xSbZEc8dxHKfEFCTaZtYW58aPIEzn3D1bs2zXSjpD0mxJsze8t67jljqO4zjF5Wmb2RpgFjAZGCQpnX2yTaJ54pppZjbJzCb17jMgWxPHcRynQArJHhkiaVA87gt8jJAAPhP459jsVOCuchnpOI7jBArJ094RuC4ucNIA3GJmd0t6EbhZ0veBZwjz6B3HcZwykle0zew5wkInmeULCfFtx3FaLXNIAAAeuUlEQVScmkQpo/fG1vwN6whfe8RxHKeOqKhob+rVyGujB1bylo7jOF2KqnjaLtyO4zgdo2rhkddGD3TxdhzHKZKqx7RdvB3HcQqn6qKdxoXbcRwnPzUj2uBet+M4tYukIyTNkzRf0vlZ6kdLelDSc5JmSRqRqBsl6X5JL0l6UdKYWF70Etc1JdppXLwdx6kl4uTCy4EjgQnASZImZDS7DLjezPYCLgQuSdRdD/zQzHYnzG9ZGcuLXuK6JkU7jYu34zg1wgHAfDNbaGYbgZuBYzPaTAAejMcz0/VR3HvE3Wsws7fM7J249VjRS1zXtGincfF2HKfKDAeWJs6zLUc9B5gaj48D+ktqAt4PrJF0u6RnJP0weu5NdGCJ67oQ7TQu3o7jlInm9BLS8XVGRr2yXJO5HPW5wCGSngEOAV4FWgnLhRwc6/cHdgZOK7DPbajLjX1fGz2Q9y1eW20zHMepddo2Ym8uzd8OVpnZpHbqW4CRifNtlqM2s2XA8QCStgemmtlaSS3AM3G9JiTdSVje+hriEtfR2865xHWSuvK0k7jX7ThOBXkSGB+zPXoBJwLTkw0kNUtKa+o3CKKcvnawpCHx/FDgRTMzOrDEdd2KdhoXb8dxyk30hM8G7iPsJ3CLmc2VdKGkY2KzKcA8SS8Dw4CL47VthNDIg5KeJ4RFrorXfB34qqT5hBh33iWu6zI8ko20cHvYxHGccmBmM4AZGWUXJI5vZUsmSOa1DwB7ZSkveonruve0M3Gv23GcrkyXE23wkInjOF2XLinaaVy8HcfpanRp0U7j4u04TlehW4h2Ghdux3HqnW4l2uBet+M49U23E+00LtyO49Qj3Va0wb1ux3Hqjy4zuaYz+MQcx+mibNoEy1+rthUlJa+nLWmkpJlxx4W5kr4cy78r6VVJz8bXJ8pvbnnpiOe9emUT857fjdUrm8pkleM4zhYK8bRbgf8ys6cl9QeekvRArPuJmV1WPvOqQ6Zw5/LAV69s4tH7D6GtrYHGxhQHHf4QTUNXV8JEx3G6KXlF28yWA8vj8XpJL1HAQt3ZaOiVd6nYmiSX992ydCRtqQaggbYUrFoxxEXbcZyyUtRAZNyMch/g8Vh0dtzE8hpJgwvpo9+Y1vyN6oSBo9+godFAKRoajYa939kcYvEBTsdxykHBA5FxUe/bgK+Y2TpJVwAXEXZauAj4EfD5LNedAZwB0KcpLCebFu63F9X3OGj/kWvZ45QnWbt4BwaOfoP+I7cOoxQaZnEcxymUglRTUk+CYN9kZrcDmNmKRP1VwN3ZrjWzacA0gAFjxm8VH+kK4t1/5NptxDoXSRF3AXccpyMUkj0iwsLcL5nZjxPlOyaaHQe80FEjulLIpFCSYRQPpThO7SPpCEnzJM2XdH6W+tGSHowh41mSRiTq2hKZdtMT5TfFPl+IYeae+ewoxMU9EDgFeF7Ss7Hsm8BJkiYSwiOLgDML6CsnXcHr7gweSnGc2iXunn45cBhhv8gnJU03sxcTzS4Drjez6yQdClxC0E6Ad81sYpaubwJOjse/Ab4AXNGeLYVkjzxK9l2DZ2Qp6zTdXbzT+IQfx6kpDgDmJzbnvRk4FkiK9gTgnHg8E7gzX6dxNxxin08QNvdtl5qdxt5vTGu3DJtk4iEUx6kJhgPJbd1b2Db1eQ4wNR4fB/SXlJ5110fSbEmPSfpUZucxLHIKcG8+Q2pWtNO4eG/BxdtxykZzFNX064yM+mzRhsyJJ+cCh0h6BjgEeJUwORFglJlNAj4D/FTSLhnX/gJ42MweyWdo3cQgPGyyhddGD/SwiQOEWbmrVgyhedjrPrErGxtbsVcLWntkVRTVXLQAIxPnI4BlyQZmtgw4HjanSE81s7WJOsxsoaRZhPkuC2Lb7wBDKHBcsOY97Uzc6w641+2kl1GY+/SePHr/Ib7+TXl5EhgvaaykXsCJwPRkA0nNktKa+g3gmlg+WFLvdBtCcseL8fwLwMeBk8wsVYghdSfa4CGTJC7etUu5FxNbtWIIbW3pZRQaWLViSFnu44CZtQJnA/cBLwG3mNlcSRdKOiY2mwLMk/QyMAy4OJbvDsyWNIcwQHlpIuvkytj2bzEd8IJ8ttR1rKHfmFYPl0Q826S2WL2yiUfum0Iq1UBDQ4qDPz6r5OGL5mGv09iYoi0FjQ0pmoe9XtL+na2JmR4zMsouSBzfCtya5bq/Ah/I0WfRAlb3ileuWPf6pQNzTk+vZVy8a4MlC0aTSjUAIpVqYMmC0SUX7aahqzno8Ic8pt3NqHvRTlOIeBcqxOuXDmTuDfuTahMNjcYepzxZV8INPmW+u9A0dLWLdTejy4h2mlziXYwQr128A6k2gTWQakuxdvEOdSfaSVzAK8+oXRazeP6YzeGRUbssrrZJThehoqLdq3cbI0evZ+ni/mW/V2a8uxghTi+5mmoLS64OHP1G2e2tFD5dvjI0DV3NwR/30IVTeqriaY8cvR6g7OKd9LqLEeL0kqsrn9uprPbVAu6Flw8PXTjloKrhkcqKd/trX2fj9TnDSbWJ1+cMr8u4drG4F17/+GSbrk9NxLQrId4ht3t1wcLb1eLaHcG98MrSWcH1PUu7BzUh2mkqEe9OTsppL9OkK8e1O4ILeHkpheAuWTA6MdnG9yztqtSUaEPlQiaQb0p88eGUclMrueMu4KVn69mNWwtuIR746pVNLJ4/lrCukdEgn2wDYBtaaV3wZrXNKCk1J9ppKine2djUuj3vrhvEsMmvM2jcms3l1ZqBWau54x4HLw25ZjcW6oGvWjGEVEqkRXv0uEXuZXdRala001RDvNfMH8RTP5xMqlU09DD2O++xzcKd6Z1XSsTrJcbuIt4xcs1ubM8DT5Ip+p4X3nWpedFOU6n8boA35zWRat0ikG/Oa9rK205SKRGv1xi7h1IKJ1uKYKHri/iU9u5D3Yg2VM7rHrzrahp6bBHIwbtu/QCsmT+IN+c1MXjX1duIeVLE18wfxIrHhpQkBp3OHa+FmHZHcS+8eIoRY88L7x5UVLT7NJamn3J73YPGrWG/8x7LKszthU6SZGvXs8dbnbKr/8i1dSnWuXARLwwXYydJxdfTHrdDinE7FLTWd7uMHL1+s+ddDgaNW8PYoxYA8Mofd2HN/EFAZuhEvDkv+1rJ2drlWwd8/dKBtDw6lvVLS78+djn7LhW+Nrjj5Kdq4ZFxO6SY/0bn/2eUM2SSzVvOFzpJ0167bLni5cwOqdXMk1z48rK58RmPTlVj2mmPu1TiXWrhzjYgOfaoBTlDJ0naC7EkSS9sVc7skHrJPMnExXtrkul/DQ0pRo9bxKhdFrt4VwhJRwA/AxqBX5nZpRn1owlbjA0B3gBONrOWRP0Awq43d5jZ2RnXTgd2NrM989lRE9uNVTtksmb+oK1CIGnS3jINW3vL6dBJLiFOU2i7fmNaGTb5dRoaDVT67JB05kk5+q4EHjIJJNP/UqlGXnl5l4L2hiz3tmfdAUmNwOXAkcAE4CRJEzKaXQZcb2Z7ARcCl2TUXwQ8lKXv44GCB7zyetqSRgLXA+8DUsA0M/uZpB2A3wFjgEXACWbWqalHpfK8iwmZtDewWKi3XAoGjVvDfl8L9+o7YE1JPeGulHnSnb3uzel/bRAm0Wjz3pDtzZT09UhKwgHAfDNbCCDpZuBY4ga9kQnAOfF4JnBnukLSfoS9IO8FJiXKtwe+CpwB3FKIIYWoYyvwX2a2OzAZOCv+hzkfeNDMxgMPxvOSUEnPO9/AYnveci4PvaOk7/W+g1eXfOPi/iPXMuKgV+pSsJN0Z687nf439v0LaWhIgVJ594bMtfmve9/b0CxpduJ1Rkb9cGBp4rwlliWZA0yNx8cB/SU1xR3afwScl+W+F8W6dwo1NK+nbWbLgeXxeL2kl6KxxxJ2Hwa4DpgFfL3QGxdCqQcrYVvvu9CBxUwKTf3rDOXa/7Le6c5edzr9b9QuiwsakMw2Oac7ed+2oY0N/1hXSNNVZjapnXpl6z7j/Fzg55JOAx4GXiU4vf8OzDCzpdKWbiRNBMaZ2TmSxhRiJBQ5EBk73gd4HBgWBR0zWy5paDF9FUopByuBbTzvkaPXM+x977B8zlB23HslG7bfWFA/xcya7CyFrkzY3XDxzi+02SbnzHt+N18NsHhagJGJ8xHAsmQDM1sGHA+bwx5TzWytpH8CDpb078D2QC9JbwGLgf0kLSJo8VBJs8xsSnuGFKwA0YjbgK+Y2brkf4w8151BiNcwZKcmdhuY4u9rixfgUnnd2Ri6x2qG7pH7Q5stNt5RD72zVGvtk1rmtdEDu6VwF0qmwOfyvj2VsF2eBMZLGkvwoE8EPpNsIKkZeMPMUsA3CJkkmNlnE21OAyaZWTqcfEUsHwPcnU+woUDRltSTINg3mdntsXiFpB2jl70jsDLbtWY2DZgGMO4DOxvAbgOD91yseJfa6y6UTO986eL+FR2kbI9yiHitLAFbDPXmdXdGJDsrsJneN9BtwiUdxcxaJZ0N3EdI+bvGzOZKuhCYbWbTCeHiSyQZITxyVjlsKSR7RMDVwEtm9uNE1XTgVODS+POuYm9eb+KdJi3iI0evh48uzeqJt7c+STnpbCil3ibiZFItr7sYIe1MTLnQa/PZk/S+PVxSGGY2A5iRUXZB4vhW4NY8fVwLXJulfBGQN0cbCvO0DwROAZ6X9Gws+yZBrG+RdDqwBPiXQm6YjVoMmRRD5kBnJQYpC6EjAl6vE3GSVFq4ixXhXMutFiL8hSzVWqw9ha4k6NQGhWSPPEr2kVOAj5bKkHr1ujMZOXo9b/xlxJZBytYUC+4azy7H/qNqIRToftusVVK4C13zOk1nMjoKEdhi7fFlXeuLmhvF6ox4V1u4V85tYvmcofQZsIHGXinaNgEp8cbcZta83FQ1jzuT9uLgXWEiTppKCXexnmoxGR0L541l2ZIR7DSqhZ13faUgge2I5+wrCdYPFRXtXo2FT5jpiHgnJ+R0RsDT4rvj3ivbzSrJvOae8z5CalMDDT1TTP73p3nlkZEse2pYRdICO0OmF96VloCthHB3xFPNltHR0JAilWLz/o4L543l2cdC6vDKZe8D2CzcnfWcPVukfqm4pz12wEZeWder4PadiXcnKVTEM8X3yB/OLEi4l88ZSmpTA5ZqILUJ3lvXm30/9wIrnh9CahM09Eyx2yFL2VD0b1JZss3ErPe0wkoJd+fFT4kXLFsyIlFuLFsygp13faXT9nSnyTVdkao8jWMHhAkshYp3R0MmSdqbFp8U9EzxXT5naEGivePeK2nomdos0Gkv/cgfzszqtVdrw+KO4Lnh5Se5MW/KxKoVQ9hpVEv0sMPEu51GtbTbRzH38myR+qWqT19HxLszwp2LpKD3OngFz964x1biWwi5BDrXxJ1q7zbfGepxhmatT8DJFodOC2kypl2uezn1g8wyp8+XjwkTx9hvHrwgZ30xYZNyiHeaJXOaWTh7KDtPCoL95CPDiopvd5R6FPA09SLetSzclYwz10NM+/brPv1UnvVA8rLX4EF295SD8rYbfecfO32vSlFTT1oxnne5vG6AUXuvYtTeq1gyp5mrzzyU1k0N9OiZ4vRf/plRe68qW5ZKewtb1TrpzRxqnVr2uCuZweHZIvVLRXPkehW4sW9avPOx28DU5nh3OVg4eyitMb7d2trAwtlhTaz00rGlWD42F+llZcu9F2YpybcHZq3QnZd3deqfirtGw/uFkfFX324/LJPP6573dDMvPDaMPSevYLd9VwGlD5nsPGklPXqmaG2FHj1Sm8MlSSo1uaeevHBfUra2qIdQiFM4VXuqihHvTOGe93Qz3zn5Y7RubKBHrxTfu/FP7Lrvqm287s6K+Ki9V3H6L/+8Ob6dDpkkzyHEwFti2caR5Z9BWC8CXsshk1oOk5QST+/relT9iRreT0V73S88NozWjQ2kUg20bgrnu0ZvO0l7oZNCBT0d3wayxriBisa9M8m2AmEtUcted3cQbk/v63rUxJNUrNe95+QV9OiVonUT9OiZYs/JK4q+ZzZBzyfkW8e42RzjziwbtfeqmlhGtpYEvFbFu6sLt6f3dT0q+gQ1qme79YWI99gBGxk7ZRnc+KfNMe1sXnZHyBdeyRXjbi/uXc0FrWrRC6/lkEkpqZU4si8G1fWo+NMzoGfwTtdtyj1ppZCQyRFTlpVMrHORFPG/r23IGuMGspZlUgurEdaKiNeacJfa2661OLKn95UGSUcAPyNsgvArM7s0o340YbeaIcAbwMlm1hLLb4/X9QT+z8yujNf0An5O2EAhBXzLzG5rz46qPTn5xLtQrztNMRNzOkJawHf78Er+niHMybh3Pjq6Jko5qKaI11q4pJTC7XHkroekRuBy4DDCfpFPSppuZi8mml0GXG9m10k6FLiEsBfBcuBDZrYhbtv4Qrx2GfAtYKWZvT/u2r5DPluq/sQM6Dm001435M7tLoeYl2ItlDSlWpmwFFRjan2ted2lwOPIXZIDgPlmthBA0s3AsUBStCcA58TjmcCdAGaWFKfebD0/5vPAbrFdCsjr/dXE01IKrzsX2cS8VEJeSvGG2vHCKy3etSLcpfK2PY7cJRkOLE2ctwAfzGgzB5hKCKEcB/SX1GRmqyWNBP4IjAPOM7NlkgbF6y6SNAVYAJxtZu1mVlT0SVHK6L2xlQ29st+2VF53Pkot5KUW7zTV9sIrKd5dUbhdrKtP60ZYubSgZ6dZ0uzE+bS4KXmabLt3ZYrRucDP447rDxN2bW8FMLOlwF6SdgLulHQr0AaMAP5iZl+V9FVCiOWU9gytylPSe2OIZ2YT73J63e1Rivh4ucQbqivgI0evr4nME8cpI6vyLBjVAoxMnI8AliUbxBj18QAxdj3VzNZmtpE0FzgYuA14B7gjVv8eOD2foVUNoqbFOxtp8c7F8H7a/Co1Ywds3PzqCOVeEyW59kk51z9JUok1UOph3RKn2/IkMF7S2JjxcSIwPdlAUnMcTAT4BiGTBEkjJPWNx4MJm6XPs7DE6h8ImSMQ9txNxsizUvXvo53xutPkEu5Sh1KK9cAzUwbLRSVj4eUOmdRCmKSrT7hxisfMWiWdDdxHSN27xszmSroQmG1m0wnie4kkI4RHzoqX7w78KJYLuMzMno91XwdukPRT4HXgX/PZUtmno20j9uZSNHjkNlX5xDufcGcjm5h3RsiL3bQhSaUEHCoj4uUMmdSCcDtOJmY2A5iRUXZB4vhW4NYs1z0A7JWjz8XAh4uxoypPhr0ZBmFziXdnvO58ZAp5ZzJSOhv7TlNJES+lgJfT6662cLu37dQqVY1pp8U7k94bW3PGuwf0HJo33l0Mydh4sfHxzsS9k6Rj4OWMg6cpRxy8Xtb7LpZ8626vXtnEvOd3Y/XKpgpZ5DgFeNqSrgGOJsza2TOWfRf4IiEGA/DN+NWhfTZt2qaoI1435B6oLJUnXowH3lnPO0mlvPBST6svR7ik2t52e9TaVHWn+1DIE3stcESW8p+Y2cT4yi/YaZa/Fl4ZtOd1F0PaE0++OkJHvO9Sed5Jyu2Fl9LrLofHXe2Mklze9tZT1RtYtWJIZQ1zui153Rgze1jSmJLfeflrsOP7tr5XDq+7vUHKQsgm3MV45MV635nCXeoZmElKOZW+FF53V/S4s8W3faq6Uy068yScLelzwGzgv8zszaJ7SHvcWcS7vQwT6LiAp0kKeaEC3tGJPeVc2KqU4ZRxO6RqVrhrDZ+q7lSLjj6hVwC7ABMJK1j9KFdDSWdImi1p9uur1mOvbhsayRUuyRUygS2DlZmvjlBsKKUzk3qSE3dKHUqBzk/sKdVAZalDJbUYJmkauppdP/B3F2ynonTIXU0uaCLpKuDudtpOA6YBTBo3xIDNwq3hCQ+7SK87F9mEuxivvJjUwlJMqS+XF97ZKfWl8rpLSS2GSZzaZuMGY8krG6ptRknp0FMpacfE6XHACx3pp1Redz464o13xPPu7JT6cnjhnfG8O+t1d8VUwHxpgI5TbgpJ+fstYXpms6QW4DvAFEkTCatcLQLOLORmtqGV1gVv0mOXwVvKcnndGR435M4wKcYTh+Ji48VO6inF5J00pRzQ3G1gqiped6nj29X2tsE9bqe6FJI9clKW4qs7c9NM4YYg3oWES7KRTcwLFfJCBbyjMzKTIt7ZtVA6K+KdCZm4cDtObVC1oGXrgjdpXbB1wknOcEmWkEk+0iGVYkIrhYRPOpP/XeqVCTsaRulMuKSjdLVQiYdJnGpRUXfFNrSx4R/r6D1+wOaygsIlkFu4C/DEYVtvvD1PvNC88I6kDaYpRxilGM+7o153rQxQ1oK37WESpxpU5enb8I91bPjHus3nBXvd2Uh74slXARTigXdk4LIUXnhn0wmLoSNed0c97q7mbYN73E7lqahot27a+nZJ4QayCnf6VRRFiHghIZSO5ICXchp9R6bTF0OlFquC0gp3tXO307hwO5Wk4p72m8t7b3WeTbgzxRu2FvCixbxIAc9FRyfxdFbA0xQj3uX2umslvu3C7VQKSUdImidpvqTzs9SPlvSgpOckzZI0IpZPlPQ3SXNj3acT13xU0tOSnpX0qKRx+eyoSnjkzeW9txLvzHAJbOt156JoES9AwIsJnRQr4qVczKoQalW4HaeekNQIXA4cCUwATpI0IaPZZcD1ZrYXcCFwSSx/B/icme1BWHzvp4md2K8APmtmE4HfAN/OZ0tlwyMZOyMX6nXn8r6zUZSIl0C805TCCy9WyMvpdRdKLcS33dt2KsABwHwzW2hmG4GbgWMz2kwAHozHM9P1Zvaymf0jHi8DVgLpZSENSGdmDCRjs+BsVMXTzhTufF53mkwRL0TMCxLwAsW73AKephgBL4fXXWx824V7Cy7cXZbhQFIAWmJZkjnA1Hh8HNBf0lY7ZEg6AOgFLIhFXwBmxImLpwCX5jNEYUPgyiDpdWBxGbpuBlaVod/OUIs2QW3a5TYVRnezabSZdWqhckn3EmzMRx/gvcT5tLhuUrqffwE+bmZfiOenAAeY2X8k2uwE/BwYS9jYdyqwh5mtjfU7ArOAU83ssVh2O/ADM3tc0nnArul75KKyedqd/APkQtJsM5tUjr47Si3aBLVpl9tUGG5T8ZhZtg1cOkILkJzcMYKMUEYMfRwPIGl7YGpCsAcAfwS+nRDsIcDeZvZ47OJ3wL35DKn+LAnHcZza50lgvKSxknoBJwLTkw0kNUtKa+o3gGtieS/gDsIg5e8Tl7wJDJT0/nh+GPBSPkN8AQfHcZw8mFmrpLOB+4BG4BozmyvpQmC2mU0nLKx3iSQjhEfOipefAHwYaJJ0Wiw7zcyelfRF4DZJKYKIfz6fLRWNaZcLSWck40+1QC3aBLVpl9tUGG6TA11EtB3HcboLHtN2HMepI+petCUtkvR8nAY6u0o2XCNppaQXEmU7SHpA0j/iz8Ht9VEhm74r6dX4Xj0r6RMVtmmkpJmSXopTer8cy6v2XrVjU7Xfqz6SnpA0J9r1vVg+VtLj8b36XRzkqrZN10p6JfFeTayUTd2Rug+PSFoETDKzquWvSvow8BZhdHjPWPY/wBtmdmlcp2CwmX29yjZ9F3jLzC6rlB0ZNu0I7GhmT0vqDzwFfAo4jSq9V+3YdALVfa8E9DOztyT1BB4Fvgx8FbjdzG6WdCUwx8yuqLJNXwLuNrNbK2FHd6fuPe1awMweBt7IKD4WuC4eX0cQgmrbVFXMbLmZPR2P1xPSm4ZTxfeqHZuqigXeiqc948uAQ4G0OFb6vcplk1NBuoJoG3C/pKcknVFtYxIMM7PlEIQB6NwSf6XjbIWVxq6pdMgmiaQxwD7A49TIe5VhE1T5vZLUKOlZwloVDxCmPq8xs/Sc/WxTqStqU2JiyMXxvfqJpN7tdOF0kq4g2gea2b6E1bfOimEBJztXALsAE4HlwI+qYUScLXYb8BUzy77QTIXJYlPV3ysza4urv40gLFi0e7Zm1bRJ0p6EiSS7AfsDOwAVCwN2R+petOPUUcxsJWHW0QHVtWgzK2K8NB03LW4/sjJgZiviQ5cCrqIK71WMhd4G3GRmt8fiqr5X2WyqhfcqjZmtIaxZMRkYJCk9KW6bqdRVsOmIGGIyM9sA/JraeQa7JHUt2pL6xcEjJPUDDgdeaP+qijEdODUenwrcVUVbgM2CmOY4KvxexYGsq4GXzOzHiaqqvVe5bKqB92qI4prLkvoCHyPE22cC/xybVfq9ymbT3xP/cEWIsdfKM9glqevsEUk7E7xrCFPyf2NmF1fBjt8SprA2AyuA7wB3ArcAo4AlwL+YWcUGBnPYNIXwdd+ARcCZ6VhyhWw6CHgEeB5Ir+f6TUIMuSrvVTs2nUR136u9CAONjQTn6hYzuzB+5m8mhCGeAU6OHm41bfozYX1oAc8CX0oMWDolpq5F23Ecp7tR1+ERx3Gc7oaLtuM4Th3hou04jlNHuGg7juPUES7ajuM4dYSLtuM4Th3hou04jlNHuGg7juPUEf8fzgAuMEpWocAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW0AAAEGCAYAAACuMsS7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnXmYVMXV/z9fBgYUlR2DAmKUGBMTQYlCjICoiJq4xbxioq9GDTE/F7KZqEnUaPTVJMYlJhoExahxieIao6IIuACKLAJixAWVJSIgICozMHN+f1S1XJqe6Z6ZnunpnvN5nvv0vVV1656+y7nnnqo6JTPDcRzHKQ5aFVoAx3EcJ3dcaTuO4xQRrrQdx3GKCFfajuM4RYQrbcdxnCLClbbjOE4R4UrbaTAKvCDp4CY+7vmS3pe0XtLXmvLYTv6QdI+k0wstR7HQ5Epb0mRJFZI+krRW0luSbpe0b1PLUgpIGiDpRUmfSHpT0klZym8r6RZJH0paI2mcpG3SypwnaamkjyU9JenzWcT4H2CTmT2dqGP3uO/HkpZI+lkWucZL2hgVcGr5f7WU7wlcARxkZtuZ2UtZZMwJSUMlWTz+OkmrJD0v6aeS2tayX5+4X898yNHCuBi4Iv0+dGrAzJp0ASYDv05s7wL8DtgIHNvU8tRD/jaFliEhSwfgA+CXQFvgUGA9MKiWfW4GXgB2BLrH9RsT+d8DVgD7ANsC1wPzgbJa6pwGnJzYLgMWAn+OdewT6zyhljrGA2Pr8N+/QXhR1JR/Szzm/LpcS2Bosl6gXTyvBnwMzAUeTuTfCfwnLgb0qeM1LANmA4/WkD8YmAVsAo5Py6sC5sQlKZOAy4HX43U4N5F+PfAG8AqwT0zvF6/hgph+QqKug+Px5wDPAbsnzx9wfPzfA2LafgmZ5iafaeAn8RjzgbuAdjH92XjfvgssAx7Mw7PxeDz+AuCm1P0L7B3/6zzgEWCHpn5uG/zfmvyAaUo7kT4OWAKohv1GxouwDlgO/A1on8jfDvgj8BbwUbxY30jcXBfGB+sj4E3g2zFvPGnKAlgMnBTXT403+XlRvgUx/Yp4rPWxvh+n1dEH+GeUdQ3wPNAFuAp4KK3sMGBt8v/keC6/H290JdJuB26tofw2wKfAwYm0g4FPEg/QFOCytPP6CTCkhjp3jA9tj0TaQXGf7RJplwHP1PJftroOtZQ9If4PS53/mN4F+Hs856viQ7kw7bpeBDxDUMAjM9Q9lAwvg1h+A/CDtPQjCMpwbZSnIsr0m5i/S7xGa4H3gGuBbRL7/xT4B1Fps/V9/DpwSvxf3wWui/WsJCjy3mnP1tXATIIR9CZwNOHl3BmoJCheAQOBGfF6Xwf0BVoD/xfrXRvv2XeAPWP902KdtwKrgbHAVOA1Nj8Lj8X/OBnoQXhxdico6U0EI+Ne4CHg1ITslwAPA/cD/5sHPbND/FWsc2Tcfol4LwOnkbjXi2VpTj7tu4GdgT1qyF9LuGk7AgfG5deJ/HHA/gQltANwDPDfmPc74CTgOzFvCLCoDrL1AXYi3Ngp3+mrBGtve+AHwP9JOgyCCwKYRLhhvwh0BX5OeGjGAIdL6pGo/wzgLjP7OO6/JsvSO+63NzDL4h0YmRXTM7EHwXJ8Oa38NsAXEnV+lm9m6wnnqqY69wE+NLPlibS9gdfjvrnIleLbklZLel3SHyRtl6mQmd0DHA5UWXCN7Baz7gQ6AV8CdgfaA0l3RWvgAsKLZhbBCs0VI5yXLfz2ZvZYPP+p//a3KNNlkloD/yK8YE4gKMoDCEo55eI5kqD8UqTfx0cSlGc1wYAYGJddokyPSCpL7H8K4evmS8ANwG3AejNbTVD26ywwHegW5fmTmS0CLiUYEIsIFvMthGdy51h3OeFZeDzuu4HwMvlC/A8dCS+N02L5dlHG++PvfwlW/fp43GUJuecB+8bjPxjPT/voyntJ0mxJR5MjZrYurraOcqeekT0ILxqAicC3c62z2dDUbwlqtrT3JJzYA3Ks52zgxbjePe775QzlRLhJjqyhnvFkt7Q/Bdpmkec+4Pdx/X8I1l7rGso+BZwf1zvF+vetx7kcB9yWlvZ94I0ayh8Yz1PSMm8V01JfJVUEP3FyvymZrlnM+y6wOC3tN8CUtLSDqN2dsS9BmbYCvgy8SHiR1VR+KFu6MXaK/6NvIm0Yia+AeJ6vjev7A5Oy1ZtI30SwLtcAx2TI3z0e67hE2tcJlvcdwNCYdliUQ/Ge2Tce89Es9/H4WNehaTJVE9wNxxCerb8QvjJ+RXgxGfCtWH4G8CHRLQS8DcxIPCcfAaMIL7NWMf3NmL4k/vcpMb0/QRn/OqYPSMj1GOHLZD3BWPmE4L4bHdNWRrl6JvY5NP6/+xJpV7D5OexI+OrI+WsUeCL+33+w2T3yAnB0XP8p8FFdn7tCL83J0k5ZRKsyZUo6VNKzkj6QtI7gZugWs/vE39cz7NqNYHFlysuV5WZWkSbPuZLmpRr0gG+lyfOWmW2qob6/AanW8pMIn/Av11C2Nj4i+LWTdCS4kGoqT9o+qfV1iTJ1qfNDgkXYELkws5fN7H0zqzazBQT/5/G1Nf6l0Sv+vp1IeyeVF632tsCxkuYQrkEPAEnHSZovaT7h078sbj+RqKs3QZk9C1wraTe25LL4+2Ks8zCCshDBhTI2HvcqggX6XWBF2nXvE38z3avtCBbjW2kyLSP4bK+NZZbH/7nBzPaJ5X4Xf1MulW9JEvA5gpsCwtfgdsCNhGdxdbyvdwH+bmY9CRb45yS1Aq4BfkawwjekyTqL4Ib4GvDDKM8KwhdGFcH6rSZ8oabYgaDI70qkDQfOj+dtcvx/vSX1Sl2vDMtOqZ3N7DDCNW5LeIFD+Ao4S9LLhK/kyvQT3dxpTkr7BGApwe+8BZLKCZ9MdxN8eDsQGt8UiyyOv30z1PsB4a2fKQ/Cm7994litCRZPkuo0eQ4gPHw/BLqaWUeC/zQpz65pn61JHgS2lzSEoLxvTqt/fZYl5R6ZS7B4kvSP6Zn4D+EB2yeR1p9g+aUUxdxkflR2fWupczbQSdLnEmlzgS9Iap9Iq02uTKTOuWottZn34m+fRFrvRF6rWOevzKxfXPYEMLMJZraXme1F+FKpituHJeralmAVP0JQIJ+dd0kXE3zGn2FmTxAUsxFcJGeYWT/gF4RrsBdwlKTFhPt6GKHdBDLfqxsIfupdE2nrCPfqy1Gm7WP6EoIVnGK3RPoLhK/HYYSXwD9iXkXi/Gwf7+m+wDtmdlYsszgeb/so/2RCW1N34GFJA2K53vEcLCS4NT8lKOi/m1mHWPf32WzkQFDwreK5SiFC21PqevU2s4Vm9l7qemVYki4XzGwDwVd+dNx+zcyGm9m+hBfEmxQbTW3as3XvkV7AbwkWwLdr2Gc7whv6u3H7S4TGj8WJMv+MdfchXOzdiS3dwB8In5B7xbydga/EvFMJVsCuhDfyVYSHY4uGyDR5Dico+76EG+1IwothfMxvT7jBrydYnGUEP+T2iTquJFgkHwMd6nkuOxJeSucRHsCDya33yHOEB617XL8pkf894H2CUtqGYMEtoPbeI9NT5ytup3qPXBfr6Bfr3KrhL7HPSKBjXO9LUC7311J+KGluDMLn8EPxvHQiuKE+SuRvAK6P6wL2zlZvvCeOJvhcZxDcMIuAL8X8M6KsnUhzLREsygWEe/WwuO8M4K8ZjplqiMx4HxPcIxMJjYE7xeVmQm+P7lGmFwnuiiuB02J9RmyMJdynUwj39wzg/ZheDjwNPBnz+8b0DgSlm3KdPQ+8nSb7bgRlf3687iMJL5PJBCt9WTw/98bzsC1BWU8BzknU8wbwfFrdVxD88orb/XN8LrZjs0usNXAPcHbc7h5/WxH88ac1hp5rzKXpDxguZgXhE3od4XP2TmC/LPuNIlgK6wmt/xexpdLenqBg3ollkr1HymP5N2LeG0TfI6FnyVjCZ/4S4Cwy9B5Jk6UV4TPyQ4Kf81aC33J8oszngQcIL4Q1hM/qzhlu9lsbeD6/RnhYPyV8Op+Ulr8AuDCxvS2hgWlNXMaR6M0Qy/wiPmyfEB7m3bLIMBKYnJa2e9z3k1jXz9Py/82WL4vJ8Vx+HO+JP1FLdywyK+1u8Tr8l6CgPyEoqCWEL5olBGt/LqEh+aIa6rV4n6yLMs0juB1eieunJ8pvIlhrc+L//Cie11/F/F0Jlv7aePzrgW0zHDOltLcnvARXsLln0op4XlbFZWk8xrr4P+bF/zeZoLQ7EizWefG/nBLrFsHnnep5c3JMPymepzmx7g3x2MujLAvjOfsvcE+GczYvXrP1hC+v1VG2WQRfe+d43DWEF1sl4cuubdx/j5h2VFq92xDcWPMIRlfGbpEZ5NmR4J55hXD//5nYvkTwq78elyupobdac15SbzCniYmug/eB4Wb2QqHlaQjRP/oCQVFNKrQ8TmGRdBfhC2dUHco/bWZjsxZ2XGkXgqjkLiT0QvDh105RI+lbBDfbRwQXzD+Bw8zsmYIKVqK0LrQALQ1J3QlujBWEfuOOU+wMIbgI2xEGEp3pCrvxaFJLu2vXrtanT58mO57jOMXLyy+/vNLMumUvWTOf26mfVVZ8lLXch6vfesLMRjTkWE1Fk1rau/TpxAszLm3KQzqOU6S0bX3kO9lL1U5lxUccfMQVWcvdd8fIrg09VlPRnPppO47jOFnIWWlLKovj/x+N2+MlvS1pTlz6NZ6YjuM4DtTNPTKa0F8zOWT5PDO7L78iOY7jODWRk9JORCO7nBBkxXGcRmTTRrF8WWcqKtpsjk9Xigjatt1Ij51W07pNKf/R/JGrpX0tYZTc9mnpl0u6iDDy7XxLC6rklBbTp3Vk6pQuDB6yioGD1hRanJJm+bLO7LBDDzp33oHQrb80MTNWr1rL8mXQa5eMseKcNLL6tCV9k62jkUGIS/xFwjDqzoQATpn2HyVppqSZKz9Y21B5nQIxfVpHRgzfn0su7suI4fszfVrHQotU0lRUtCl5hQ0gic5dOoQvCicncmmIPIC0aGSS7jCz5RaoIHSs3y/TzmY2xswGmNmArt3So3U6xcLUKV2orBRVVa2orBRTp3QptEiljVHyCjuFpNJ2AeWZrErbzC4ws55m1ocQGGiSmZ2UmnklDsk+hhDQxSlRBg9ZRXm5UVZWTXm5MXiIf8qWOm1ad2af/t/gq18ZxFFHncCaNfV3iX1+16+wcqXfM/mgIf2075Q0jxCBqyubA607JcjAQWt4/MkZXPLbRTz+5Az3abcAttlmG2bNfo5X5k2jc+dO/PUvHs+pOVCnEZFmNpkQ/hEzG1ZrYafkGDhojSvrFsrAgfsx75UFn23/8Q/X889/PkBFRQXHHPNNLvnthQAce+x3WfLeUjZsqOCcc89k1KhTCyRx6eIBoxynRJg2rYwpk1szZOgmBg2qylu9VVVVTJo0hdNOOxmAJ5+cxKJFbzJ9xiTMjKOPPpGpU59n8OADGDfuL3Tu3IlPP/2U/fcbxre/fRRdunTOcoTGY2N5K97vlT4bXnHjSttxSoBp08o49JD2VFZCeXlbJj71cYMV96effso+/b/B4sXvsu++/Tj00IMAmPjkJCZOnMS++xwIwPr1H7No0ZsMHnwAf77+Jh588FEA3ntvKYsWvVlQpV2KeOwRxykBpkxuTWUlVFWJysqw3VBSPu23F8+jsnIjf/1LmMrUzPjl+T9l1uznmDX7OV5fNJvTT/9fJk9+lqefnsLzL0xk9pzn6d//K2zYkD7nr9NQXGk7TgkwZOgmysuhrMwoLw/b+aJDhw5ce92VXH31DWzcuJHhhx3M+FvvYP369QAsXbqMFSs+YO3adXTs1IFtt92W1157nenTZ+ZNBmcz7h5xnBJg0KAqJj71caP4tAH699+br+79Ze6++35OPnkkry38Dwd8fTgA223Xnr/fPoYRIw7hb3+7lX57f50v7NGXgQMHZKnVqQ+utB2nRBg0qCqvynrdR0u32H744Xs+Wz939I84d/SPttrnsccyx4976+15eZOrpePuEcdxnByQ9BNJCyTNl3SXpHZp+dckQlW/LmlNIu+quN98SSck0iXp8lh+oaRzs8nhlrbTLPHgVE5zQtLOwLnAl8zsU0n3EkaIj0+VMbOfJMqfA/SP60cC+wD9gLbAFEn/NrN1wKlAL+CLZlYd55CtFVfaTrMjFZyqslKUl+/uIzCd5kJrYBtJG4FtgWW1lD0RuDiufwmYYmabgE2S5gIjgHuBHwHfNbNqADNbkU0Id484zQ4PTuU0N8xsKfBHwmzzy4G1ZvZkprKSdgF2BSbFpLnA4ZK2ldQVOIhgXQPsBpwQI6H+W1LfbLK40naaHR6cyikAXVMhpOMyKpkpqRNwNEEZ7wS0l3RSDXWNBO4zsyqAqNwfA14A7gKmAak+mW2BDWY2ALgZuCWboO4ecZodqeBU7tN2Gky5oGdOtunKqDhr4hDgbTP7AEDSBODrwB0Zyo4EzkommNnlhJm/kPQPYFHMWgLcH9cfIIS5rhW3tJ28M31aR35/5W4+UUKRkwrNmloWL36nxrKTJz/Lt751Qo35JcC7wMDo4hBwMGHO3C2QtAfQiWBNp9LKJHWJ618FvgqkXCsPAqnge0OA17MJ4pa2k1fy0YjoDZHNg9QwdgfMbIak+4BZBNfGbGCMpEuBmWb2cCx6InC3mSWndWgDPBsntVgHnBQbJQGuJIS5/gmwHjgjmyxuaTtZqYvlnGxErKgQl13at84WtzdENl8WL36HIYMPZ8C+gxmw72BeeGHGVmVeemkW++5zIG+9tZiPP/6Y0087i/33O4h99zmQhx76VwGkzg9mdrGZfdHM9jKzk82swswuSihszOwSMzs/bb8NZvaluAw0szmJvDVmdqSZfcXMBpnZ3GxyuKXt1Epdrd7QiLg7FRXVVFeLSU935fnnOtfJWk7VUVnZsIbIFtfXe9qLMPk5GPoNGJRx9r86kYryB9Bn112YMOFOunfvxhNPPkC7du1YtOhNvvfd03nxpcmf7fPCCzMYfe4veODBf9C7dy9+deGlHDRsMONu+Qtr1qxh4P4Hc8ghQ2nfvn2D5WupuNJ2amVLq7eaqVO61KoAU42Il13al0lPd6W6Wp9Zy7kqznw0RLY4F8u0F9EhRxNjs2JPPdRgxZ3JPbJx40bOOec85s6ZT1lZK15//c3P8l5b+B/O/OGPefyJCey0Uw8AJk6cxCOP/Js/Xf1nADZsqODdd5ew5557NEi2lowrbadW6mP1Dhy0ht9ctIjnn+sclWbdreWGzpJT15dN0TP5OaisRFVVWGVl2M6DtZ3Otdf8lR27d2f2nOeorq5m2212/CyvR4/PsWHDBmbPfuUzpW1m/PO+v7PHHlm7Hzs54j5tp1bqOzdkU80pWZO/vcX19R76jWBhl5VBeXnYbgTWrl1Hjx470qpVK26//W6qqjYHqOrQsQOPPHovv/7VpUye/CwAw4cfzA1/HkOqXW727KwuWycLbmk7Wamv1dvYc0rW5gJpcX29B+0XXCJ59Gln4kf/7wy+c/zJ3HffQwwdeuBWvukdd+zOQw/fzZFHfIex427g1785j5/8+AL67X0AZsYufXrzyCP31FC7kwuutJ2iJZsLpMVNRDxov7wq6/TQrAB9++7GnLkvfLZ9xf+F8BpDhx7I0KFh+rHevXsxb/70z8rc9Ldr8yaT40rbKWLy1cvEKV3alFexc6/1hRYjr7jSdurF9GkdueP2nQE46eSlOVu0+eyG1+JcII6DK22nHkyf1pHhh+xPRUVox75tfE8mPp29sbExuuG1OBeI0+LJufdIHD8/W9KjcXtXSTMkLZJ0j6TyxhPTaU4EX3IrQIDYuDG3UYt33L4zGza0avSRjiUR+0Sw5Ujo0sXMwq3k5ERduvyNZssAKVcB15hZX+BD4PR8CuY0X4IvuRowwGjTJrs/efq0jtw2vidBDxmtWzeOD3rszT05eNhALr7oC4wYvn/RKu62bTeyetXaklfcZsbqVWtp23ZjoUUpGnJyj0jqCRxJCC340xjlahjw3VjkNuAS4MZGkNFpZgwctIYnn5pRJ5/21CldqKoKljlUc9iIFXl3a0yf1pHR5+zFpk3hOBUVdRuJ2ZzosdNqli+DD1auCu/GUkXhBdVjp9WFlqRoyNWnfS3wC2D7uN0FWJOIVLUE2DnPsjnNmLr6kgcPWUVZ2e5UVRkgnni8O2Nv7snqVW3z1og4dUoXqqshvBiMVmUUbY+S1m2MXrsUp+xO45LVPSLpm8AKM3s5mZyhaEZ7QNKo1GwQKz9YW08xnWJn4KA1nHLqEkJ0yuAHH33OXlxycd+8uTEGD1lF27ZGq1bVtG5tXHf9/KK0sh2nNnLxaR8AHCVpMXA3wS1yLdBRUspS70kNk1ya2RgzG2BmA7p265AHkZ1i5aSTl9KuXTVlZdW0KoPqavLaKJnqAvjbSxfx9DPTOeMHS/IgteM0L7K6R8zsAuACAElDgZ+b2fck/RM4nqDITwEeakQ5nRIgpVTvuH1n3n+/nCce786mTfkdGONdAJ1SpyH9tH8J3C3pd4RZHMblRySn1Lnj9p5UVoqyMuO009+r0+Acx2np1Elpm9lkYHJcfwtonKg0TsmSjBcC1fTuvcEVttNolLc2enevKLQYecVDszpNSosLmeo4ecaHsTtNiscLcZyG4UrbaXKSjYUtbh5Hx2kgrrSdgtHi5nF0nDzgPm2nYGw5iUHjBZBynFLClbZTMLxR0nHqjittp2A01eS/jtNQJO0haU5iWSfpx2llhkpamyhzUSJvtKT5khYk95N0iaSliX2OyCaL+7SdguIjGJ1iwMz+A/SDMLcAsBR4IEPRZ83sm8kESXsBPyCMa6kEHpf0LzNbFItcY2Z/zFUWt7Qdx3HqxsHAm2b2To7l9wSmm9knMTLqFODY+h7clbZTb0pihhjHqTsjgbtqyBskaa6kf0v6ckybDwyW1EXStsARQK/EPmdLekXSLZI6ZTu4u0eceuHd9ZxioG0Z7LZDTrNIdJU0M7E9xszGpBeK0yoeRQyil8YsYBczWx990w8Cfc1soaSrgInAemAukJqL4EbgMkJo68uAq4HTahPULW2nXnh3PafEWJkKIR2XrRR25HBglpm9n55hZuvMbH1cfwxoI6lr3B5nZvuY2WBgNbAopr9vZlVmVg3cTA7xnFxpO/XCu+s5LZQTqcE1IulzcSpGJO1H0K+r4nb3+NsbOC5Vh6QeiSqOJbhSasXdI0698BgiTksj+qMPBX6YSDsTwMxuIswv8CNJm4BPgZG2eWbm+yV1ATYCZ5nZhzH995L6Edwji5N114Qrbafe1NZdb+zNPXlgQg+OPW65zyDjlARm9glhftxk2k2J9RuAG2rY98Aa0k+uqxyutJ0Gkx70aezNPTnrR18B4KmJXQFccTtOnnCl7TSITL1IHpiQctOFWdEfmNDDlbbj5AlviHQaRKZeJMcetzzmBnfe5m3HcRqKW9pOgwi9SHansnLLXiRt2lSzcWMr2rSpZq+91hdYSscpHdzSdhpEpqBPU6d0obpagKiu9j7cjpNP3NJ2Gkx6L5KarG/HcRqOK20n73gfbqe50LaVsev2VYUWI6+40nYaBQ+56jiNg/u0nWaLRxF0nK3JamlLagdMBdrG8veZ2cWSxgNDgLWx6KlmNqexBHVaFh5F0HEyk4t7pAIYFsMNtgGek/TvmHeemd3XeOI5LZUt+39XM3VKF1fajkMO7hELpDratolLTgFqHae+eBRBx8lMTj5tSWWS5gArgIlmNiNmXR5nXLhGUttGk9Jpcfikv46TmZx6j5hZFdBPUkfggThR5QXAf4FyYAzwS+DS9H0ljQJGAfTu3S1PYjstAe+B4jhbU6feI2a2BpgMjDCz5dF1UgHcSg0zLpjZmNRsEF27dWiwwI7jOC2ZrEpbUrdoYSNpG+AQ4LXUjAtxpoZjyGHGBcdxHKdh5OIe6QHcJqmMoOTvNbNHJU2S1I0Qf3MOcGYjyuk4juOQg9I2s1eA/hnShzWKRI7jOHmivMzo2b6y0GLkFR8R6TiOU0S40nYcxykiXGk7juMUEa60HcdxighX2o7jOEWEK23HcZwckNRR0n2SXpO0UNKgGsp9TVKVpOMTaadIWhSXUxLpl0t6T1LOE6m60nYcx8mN64DHzeyLwN7AwvQCcTzLVcATibTOwMXA/oSR4xdL6hSzH6GG0eQ14UrbcRwnC5J2AAYD4wDMrDKG9UjnHOB+QnC9FIcRAu2tNrMPgYnAiFjPdDNbXhdZXGk7juNk5/PAB8CtkmZLGiupfbKApJ2BY4Gb0vbdGXgvsb0kptULV9qO4zjQVdLMxDIqLb81sA9wo5n1Bz4Gzk8rcy3wyxgVNYkyHK/ecxL4xL6O45QsrVuJ7tuU51J0pZkNqCV/CbAkMZfAfWyttAcAd4cYenQFjpC0Ke47NFGuJyFaar1wpe04jpMFM/tv7OWxh5n9BzgYeDWtzK6p9TiH7qNm9mBsiLwi0fg4nDAfQb1w94jjOE5unAPcKekVoB9BEZ8pqdYIp2a2GrgMeCkul8Y0JP1e0hJgW0lLJF2STQi3tB3HcXLAzOYQXCBJ0hsdU2VPTdu+BbglQ7lfAL+oixxuaTuO4xQRrrQdx3GKCFfajuM4RYQrbcdxnCLClbbjOE4R4UrbcRyniHCl7TiOU0R4P23HcUqWMrVm+zZdCy1GXnFL23Ecp4hwpe04jlNEZFXaktpJelHSXEkLJP02pu8qaUacPuceSTmF0nIcx3HqTy6WdgUwzMz2JgRJGSFpIGFKnWvMrC/wIXB644npOI7jQA5K2wKpSSfbxMWAYYSYsgC3Acc0ioSO4zjOZ+Tk05ZUJmkOYd6zicCbwBoz2xSLNGj6HMdxHCc3clLaZlZlZv0IMy7sB+yZqVimfSWNSk3hs/KDtfWX1HEcx6lb75E4+/BkYCDQUVKqn3dPYFkN+4wxswFmNqBrtw4NkdVxHKfFk0vvkW6SOsb1bYBDgIXAM8DxsdgpwEONJaTjOI4TyGVEZA/gNkllBCV/r5k9KulVwiSWvwNmA+MaUU7HcVoQ5WVHFFqEZktWpW1mrwBn9q/JAAASeklEQVT9M6S/RfBvO47j5JXKqsfyUo8MyjdW5aWu5oKPiHQcxykiXGk7juMUEa60HcdxighX2o7jOEWEK23HcZwiwpW24zhODkhaLGmepDmSZmbI/56kV+LygqS9E3kjJP1H0huSzk+knx3TTFJOszX4zDWO4zi5c5CZrawh721giJl9KOlwYAywfxzj8hfgUEKcppckPWxmrwLPA48SRprnhCttx3GcPGBmLyQ2pxPCe0AYz/JGHNuCpLuBo4FXzWx2TMv5OO4ecRzHyQ0DnpT0sqRRWcqeDvw7ru8MvJfIa1BUVLe0HccpXTZVYmuW5lKya5qfeoyZjUkrc4CZLZPUHZgo6TUzm5pekaSDCEr7G6mkDMfLGBU1F1xpO47jwEozG1BbATNbFn9XSHqA4PbYQmlL+iowFjjczFbF5CVAr0SxGqOi5oK7RxzHcbIgqb2k7VPrwHBgflqZ3sAE4GQzez2R9RLQN86rWw6MBB6uryyutB3HcbKzI/CcpLnAi8C/zOxxSWdKOjOWuQjoAvw12S0wzvB1NvAEIaz1vWa2AEDSuZKWEKzvVySNzSaIu0ccx3GyEHt+7J0h/abE+hnAGTXs/xiwVehCM7seuL4usril7TiOU0S40nYcxykiXGk7juMUEa60HcdxighX2o7jOEWEK23HcZwiwrv8OY5TumzaCB/8t9BS5BW3tB3HcYoIV9qO4zhFhCttx3GcIiKr0pbUS9IzkhZKWiBpdEy/RNLSOMZ+jqQjGl9cx3Gclk0uDZGbgJ+Z2awY5eplSRNj3jVm9sfGE89xHMdJklVpm9lyYHlc/0jSQhow64LjOI5Tf+rk05bUB+gPzIhJZ8eZh2+R1CnPsjmO4zhp5Ky0JW0H3A/82MzWATcCuwH9CJb41TXsN0rSTEkzV36wNg8iO47jtFxyUtqS2hAU9p1mNgHAzN43syozqwZuJky9sxVmNsbMBpjZgK7dOuRLbsdxnBZJLr1HBIwDFprZnxLpPRLFjiVt6h3HcYqH6dM68vsrd2P6tI6FFsXJQi69Rw4ATgbmSZoT0y4ETpTUjzCr8GLgh40ioeM4jcr0aR0ZMXx/KitFefnuPP7kDAYOWlNosfLDxk3Y+6uylysicuk98hyZp4Dfauocx3GKj6lTulBZKaqqWlFZWc3UKV1KR2mXID4i0nFaOIOHrKK83Cgrq6a83Bg8pLQs01LDo/w5Tgtn4KA1PP7kDKZO6cLgIavcym7muNJ2nCJm+rSOeVG2AwetcWVdJLjSdpwipaQbEJ0acZ+24xQpWzYgiqlTuhRaJKcJcKXtOEWKNyA2PZLKJM2W9GiGvLaS7pH0hqQZMexHKu+CmP4fSYcl0kfEtDcknZ+LDO4ecZwixRsQC8JoYCGwQ4a804EPzWx3SSOBq4ATJH0JGAl8GdgJeErSF+I+fwEOBZYAL0l62MxerU0AV9qO00jkq5GwNrwBsemQ1BM4Ergc+GmGIkcDl8T1+4Ab4ojyo4G7zawCeFvSG2wO+/GGmb0V6787lnWlXWw0xcPuNC75aCT0+6DZcS3wC2D7GvJ3Bt4DMLNNktYCXWL69ES5JWwOb/1eWvr+2YRwpd3M8B4BpUFDRxn6fZAfrLKK6ndz8vV3lTQzsT3GzMakNiR9E1hhZi9LGlpDHZlGjlst6ZnaFC2boN4Q2czwHgGlQUMbCf0+aHJWpqKRxmVMWv4BwFGSFgN3A8Mk3ZFWZgnQC0BSa6ADsDqZHukJLKslvVZcaReA2iKqeY+A0iDVSHjJbxfVy0r2+6B5YWYXmFlPM+tDaFScZGYnpRV7GDglrh8fy1hMHxl7l+wK9AVeBF4C+kraVVJ5rPfhbLK4e6SJyfbZ6z0CSoeGNBL6fVAcSLoUmGlmDxNCWN8eGxpXE5QwZrZA0r2EBsZNwFlmVhX3Pxt4AigDbjGzBdmO6Uq7icnF1+k9AhxonPvAGzcbjplNBibH9YsS6RuA79Swz+WEXifp6Y9Rx4ipLUppN4cbNnz27k5lpX/2Ok2LN26WBi1GaTeXG9Y/e51CUei42c3BaCoFWozSLvQNm8TdHy2bsTf35IEJPTj2uOWc8YMlTXbcQn7lNRejqRRoMUo71xvWrQGnMRl7c0/O+tFXAHhqYleAJlPchfzKa05GU7HTYpR2LjdsrtZAroq9oS8Af4EUH9mu2QMTUvNhCzAemNCjSa3tQn3leVtO/mgxShuy37C5WAN1UewN+Rz0z8niI5drduxxy6OFbZ9ttwS8LSd/tCilnY1crIFcP/Ma+jnon5PFRy7XLGVVF8KnXWgKYeVb5SY2vr22SY/Z2LjSTpCLNZDrZ15DPwf9c7L4yPWanfGDJfVW1u4yc1xpp5HNGsj1M6+hn4P+OVl8NPY1c5eZA66060Wun3kN/Rz0roHFR2NeM3eZOZBDwChJvSQ9I2mhpAWSRsf0zpImSloUfzs1vrilRW2BoxwnHQ8i5UBulvYm4GdmNkvS9sDLkiYCpwJPm9mVcW6z84FfNp6opYV/6jp1JZ/uF/eNFy9ZlbaZLQeWx/WPJC0kzLpwNDA0FruNEEDFlXaO+KeuUx/y4X5xg6G4qVM87Ti7cH9gBrBjVOgpxd4938I1NU3prvBPXadQ+AQLxU3ODZGStgPuB35sZuvCfJU57TcKGAXQu3e3+sjYJDS19eG9Q5z60lDXhncnLW5yUtqS2hAU9p1mNiEmvy+ph5ktl9QDWJFp3zhtzxiAfQf0zTr/WaEohLvCe4c4dWX6tI4cMmwgGzeKNm2MpyZN9+6kLYxceo+IMCPDQjP7UyIrObXOKcBD+Rev/tTV1eHuCqcYuPqPu7JxowCxcaO4+o+71quegYPW8Ivz33SFXYTkYmkfAJwMzJM0J6ZdCFwJ3CvpdOBdapixoRDUx9Xh1ofTlGRyceTi9li2rF2t286WVFfCx0s2FVqMvJJL75HnyDwFPMDB+RUnP9TX1eHuCqcpyGRUADkZGt8/7T1mvtSRVMCp75/2XlOK7jQDSnJEZL4bWpoqFKvTMshkVAA5GRotOeCUEyhJpZ1PV8fYm3sy+ty9qK4SbdtWN1ooVqflUJNRkauh0ZCAUyncwCheSlJpQ/4GIYw+Zy82bQoNPxUVapRQrP4AtSxqMiqaqk3FDYzipmSVdj6YOqUL1dWQmmWkVRl5D8XqD1DLJJNR0VRtKj4at7hxpV0Lg4esom3b3amoqKZVK7ju+vl5C8Wasq7ffbedP0AlSnP9gmqMwTXN9b/mC0m9gL8DnwOqgTFmdl1amfOA78XN1sCeQDczWx0D7f2AYAHebGbXxn32Bm4CtgMWA98zs3W1yeJKuxbqqohztZSS1nVZmdG6tQHeP7yUaM5fUPnu3tqc/2seyRg4z8xeTRUwsz8AfwCQ9C3gJ1Fh70VQ2PsBlcDjkv5lZouAscDPzWyKpNOA84Df1CaIK+0sNMYna/LzFKo57fT36N17Q8laKS2R5u6CyOd93dz/az6oJXDeqzXsciJwV1zfE5huZp8ASJoCHAv8HtgDmBrLTQSeIIvSrlPAKCc/pI++POnkpT46rcRoSSNsW9J/ha0C52XK3xYYQQj9ATAfGCypS8w7AuiVyDsqrn8nkV4jRWtpF7MPzUdflj4t6RqXyH/tKmlmYntMjJu0BemB82qo61vA82a2GsDMFkq6imBJrwfmEtwtAKcB10u6iBAapDKboEWptEvBh5Zvt0sxv8RKlZY0wra5/teqjWLtf8tzKbrSzAbUVqCGwHmZGMlm1wgAZjaOEMMJSVcAS2L6a8DwmP4F4Mhsghal0m4JPrS6UAovMcdpztQSOC+9XAdgCHBSWnp3M1shqTdwHDAoLb0V8GtCT5JaKUqfdqn60Oo7CYMHtXecRicVOG+YpDlxOULSmZLOTJQ7FnjSzD5O2/9+Sa8CjwBnmdmHMf1ESa8DrwHLgFuzCVKUlnaJ+NC2oCHWsge1Lz7cnVVcZAmclyw3HhifIf3AGspfB1yXKa8milJpQ/P1odWXhrh8SvElVsq4O8tpCEWrtEuNhlrLpfYSK2W8TcZpCK60mwluLbccki/o1q2Nd99tx/RpHf2aOzlRlA2RpYpPAdUySL2gTzv9PczglnG9GDF8/zo3QDstE1fajlMABg5aQ+/eG6iq8l4/Tt1wpe04BaJUu646jYv7tB2nQHg7hlMfXGk7TgHxXj+Ny6aN8MGyQkuRX9w94jiOU0S40nYcxykiXGkXAfWNSeI4TumR1act6Rbgm8AKM9srpl1CmD7ng1jsQjN7rLGEbMn4kGfHcZLkYmmPJ8zCkM41ZtYvLq6wGwmP4Oc4TpKsStvMpgKrm0AWJwPel9dxnCQN6fJ3tqT/BWYSZin+MNsOTt3xvryO4ySpr9K+EbgMsPh7NWGus62QNAoYBdC7d7d6Hq5l4315HcdJUa/eI2b2vplVmVk1cDOwXy1lx5jZADMb0LVbh/rK6TiO41BPS1tSDzNbHjePJUwD7ziO06yorKxmyTsVhRYjr+TS5e8uYChhivklwMXAUEn9CO6RxcAPG1FGx3EcJ5JVaZvZiRmSxzWCLI7jOE4WfESk4zhOEeFK23Ecp4hwpe04jlNEuNJ2HMcpIlxpO47j5ICkWyStkJSxi7OkoZLWSpoTl4sSeaMlzZe0QNKPE+l7S5omaZ6kRyTtkE0OV9qO4zi5MZ7MwfOSPJsIpHcpgKS9CFFR9wP2Br4pqW8sPxY438y+AjwAnJdNCFfajuM4OdCA4Hl7AtPN7BMz2wRMIQxKBNgDmBrXJwLfzlaZK23HcZz8MUjSXEn/lvTlmDYfGCypi6RtgSOAXom8o+L6dxLpNSIzy7fQNR9M+gB4pxGq7gqsbIR6G0JzlAmap1wuU260NJl2MbMGRZmT9DhBxmy0AzYktseY2ZgM9fUBHk1NCJOWtwNQbWbrJR0BXGdmfWPe6cBZwHrgVeBTM/uJpC8C1wNdgIeBc82s1qD5Taq0GwtJM81sQKHlSNIcZYLmKZfLlBsuU+GpTWlnKLsYGGBmK9PSrwCWmNlf09K/ANxhZjUG4AN3jziO4+QFSZ+TpLi+H0G/rorb3eNvb+A44K609FbAr4Gbsh2nIZMgOI7jtBhqCJ7XBsDMbgKOB34kaRPwKTDSNrsy7pfUBdgInJWYNOZESWfF9QnArdnkKBWlvZXvqRnQHGWC5imXy5QbLlMBqSF4XjL/BuCGGvIOrCH9OuC6ushREj5tx3GcloL7tB3HcYqIolfakhbHIaBzJM0skAxbDW+V1FnSREmL4m+nZiDTJZKWJobZHtHEMvWS9IykhXE47+iYXrBzVYtMhT5X7SS9GPv8LpD025i+q6QZ8VzdI6m8Gcg0XtLbiXPVr6lkaokUvXukpm41TSzDYEL/y7+nugJJ+j2w2syulHQ+0MnMfllgmS4B1pvZH5tKjjSZegA9zGyWpO2Bl4FjgFMp0LmqRab/obDnSkD72Oe3DfAcMBr4KTDBzO6WdBMw18xuLLBMZxK6wd3XFHK0dIre0m4O1DC89Wjgtrh+G0ERFFqmgmJmy81sVlz/CFgI7EwBz1UtMhUUC6yPm23iYsAwIKUcm/pc1SST04SUgtI24ElJL0saVWhhEuyYmvw4/nYvsDwpzpb0SnSfNKnLJkkcpNAfmEEzOVdpMkGBz5WkMklzgBWEuBRvAmti/AqAJTTxCyZdJjNLnavL47m6RlLbppSppVEKSvsAM9sHOBw4K7oFnMzcCOwG9AOWA1cXQghJ2wH3Az82s3WFkCGdDDIV/FyZWZWZ9QN6EiLE7ZmpWCFlUohgdwHwReBrQGegydyALZGiV9pmtiz+riCENqx1CGgT8n70l6b8pisKLA9m9n586KqBmynAuYq+0PuBO81sQkwu6LnKJFNzOFcpzGwNMBkYCHSUlBpf0RNYVmCZRkQXk5lZBWFwSHN5BkuSolbaktrHxiMktQeGE6JmNQceBk6J66cADxVQFuAzhZjiWJr4XMWGrHHAQjP7UyKrYOeqJpmawbnqJqljXN8GOITgb3+GMPIOmv5cZZLptcQLVwQfe3N5BkuSou49IunzBOsawujOf5jZ5QWQ47PhrcD7hOGtDwL3Ar2Bd4HvmFmTNQzWINNQwue+AYuBH6Z8yU0k0zeAZ4F5QHVMvpDgQy7IuapFphMp7Ln6KqGhsYxgXN1rZpfGe/5ughtiNnBStHALKdMkoBsgYA5wZqLB0skzRa20HcdxWhpF7R5xHMdpabjSdhzHKSJcaTuO4xQRrrQdx3GKCFfajuM4RYQrbcdxnCLClbbjOE4R4UrbcRyniPj/RM/r5NlkxVQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-1-b1d432bd1fa3>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m    100\u001b[0m     \u001b[1;32mfor\u001b[0m \u001b[0m_\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m3\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    101\u001b[0m         \u001b[0mopt_D\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mzero_grad\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 102\u001b[1;33m         \u001b[0mD_loss\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mbackward\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mretain_graph\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mTrue\u001b[0m\u001b[1;33m)\u001b[0m    \u001b[1;31m# reusing computational graph\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    103\u001b[0m         \u001b[0mopt_D\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mstep\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    104\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\AppData\\Local\\Continuum\\anaconda3\\lib\\site-packages\\torch\\tensor.py\u001b[0m in \u001b[0;36mbackward\u001b[1;34m(self, gradient, retain_graph, create_graph)\u001b[0m\n\u001b[0;32m     91\u001b[0m                 \u001b[0mproducts\u001b[0m\u001b[1;33m.\u001b[0m \u001b[0mDefaults\u001b[0m \u001b[0mto\u001b[0m\u001b[0;31m \u001b[0m\u001b[0;31m`\u001b[0m\u001b[0;31m`\u001b[0m\u001b[1;32mFalse\u001b[0m\u001b[0;31m`\u001b[0m\u001b[0;31m`\u001b[0m\u001b[1;33m.\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     92\u001b[0m         \"\"\"\n\u001b[1;32m---> 93\u001b[1;33m         \u001b[0mtorch\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mautograd\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mbackward\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mgradient\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mretain_graph\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcreate_graph\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     94\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     95\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0mregister_hook\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mhook\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\AppData\\Local\\Continuum\\anaconda3\\lib\\site-packages\\torch\\autograd\\__init__.py\u001b[0m in \u001b[0;36mbackward\u001b[1;34m(tensors, grad_tensors, retain_graph, create_graph, grad_variables)\u001b[0m\n\u001b[0;32m     87\u001b[0m     Variable._execution_engine.run_backward(\n\u001b[0;32m     88\u001b[0m         \u001b[0mtensors\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mgrad_tensors\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mretain_graph\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcreate_graph\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 89\u001b[1;33m         allow_unreachable=True)  # allow_unreachable flag\n\u001b[0m\u001b[0;32m     90\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     91\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "import torch\n",
    "import torch.nn as nn\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline  \n",
    "\n",
    "# reproducible\n",
    "torch.manual_seed(2019)\n",
    "np.random.seed(2019)\n",
    "\n",
    "# data\n",
    "data = pd.read_csv('datasets/toy.csv',header=None)\n",
    "data=data.values\n",
    "\n",
    "# hyper-parameters\n",
    "BATCH_SIZE = 64\n",
    "LR_G = 0.001 # learning rate for generator\n",
    "LR_D = 0.001 # learning rate for discriminator\n",
    "random_neuron = 2 # number of neurons as input of generator\n",
    "components = 2 # total number of points that G generated\n",
    "\n",
    "def create_G():\n",
    "    G = nn.Sequential(\n",
    "    nn.Linear(random_neuron, 128),\n",
    "    nn.ReLU(),\n",
    "    nn.Linear(128, 128),\n",
    "    nn.ReLU(),\n",
    "    nn.Linear(128, components),\n",
    "    )\n",
    "    return G\n",
    "\n",
    "number_G = 20\n",
    "G_group = {}\n",
    "for k in range(number_G):\n",
    "    G_group['G'+str(k)] = create_G()\n",
    "\n",
    "# # define G/generator\n",
    "# G = nn.Sequential(\n",
    "#     nn.Linear(random_neuron, 128),\n",
    "#     nn.ReLU(),\n",
    "#     nn.Linear(128, 128),\n",
    "#     nn.ReLU(),\n",
    "#     nn.Linear(128, components),\n",
    "# )\n",
    "\n",
    "# define D/discriminator\n",
    "D = nn.Sequential(\n",
    "    nn.Linear(components, 128),\n",
    "    nn.ReLU(),\n",
    "    nn.Linear(128, 128),\n",
    "    nn.ReLU(),\n",
    "    nn.Linear(128, 1),\n",
    "    nn.Sigmoid(),\n",
    ")\n",
    "\n",
    "# optimization\n",
    "opt_D = torch.optim.Adam(D.parameters(), lr=LR_D)\n",
    "# opt_G = torch.optim.Adam(G.parameters(), lr=LR_G)\n",
    "opt_G_group = {}\n",
    "for k in range(number_G):\n",
    "    opt_G_group['G'+str(k)] = torch.optim.Adam(G_group['G'+str(k)].parameters(), lr=LR_G)\n",
    "\n",
    "plt.ion() \n",
    "\n",
    "D_loss_history = []\n",
    "G_loss_history = []\n",
    "\n",
    "for step in range(10000):\n",
    "    \n",
    "\n",
    "    # random samples of real data\n",
    "    idx = np.random.choice(len(data), BATCH_SIZE)\n",
    "    #idx = range(len(data))\n",
    "    selected_real = data[idx]\n",
    "    selected_real = torch.from_numpy(selected_real).float()\n",
    "    \n",
    "    D_loss_group = {}\n",
    "    G_loss_group = {}\n",
    "    D_loss = 0\n",
    "    G_loss = 0\n",
    "    for k in range(number_G):\n",
    "    \n",
    "        # random noises\n",
    "        G_noise = torch.randn(BATCH_SIZE, random_neuron)\n",
    "        G_data = G_group['G'+str(k)](G_noise)\n",
    "\n",
    "        prob_real = D(selected_real) # D try to increase this prob\n",
    "        prob_fake = D(G_data) # D try to decrease this prob\n",
    "\n",
    "        D_loss_group['G'+str(k)] = - torch.mean(torch.log(prob_real) + torch.log(1. - prob_fake))\n",
    "        G_loss_group['G'+str(k)] = torch.mean(torch.log(1. - prob_fake))\n",
    "    \n",
    "        D_loss = D_loss + D_loss_group['G'+str(k)] / number_G\n",
    "        G_loss = G_loss + G_loss_group['G'+str(k)] / number_G\n",
    "    \n",
    "#         D_loss_history.append(D_loss)\n",
    "#         G_loss_history.append(G_loss)\n",
    "\n",
    "    for _ in range(3):\n",
    "        opt_D.zero_grad()\n",
    "        D_loss.backward(retain_graph=True)    # reusing computational graph\n",
    "        opt_D.step()\n",
    "        \n",
    "    for k in range(number_G):\n",
    "        opt_G_group['G'+str(k)].zero_grad()\n",
    "        G_loss_group['G'+str(k)].backward()\n",
    "        opt_G_group['G'+str(k)].step()\n",
    "\n",
    "    if step % 500 == 0:\n",
    "        plt.cla()\n",
    "        plt.plot(selected_real.data.numpy()[:,0],selected_real.data.numpy()[:,1],'.b', label='Real')\n",
    "        plt.plot(G_data.data.numpy()[:,0],G_data.data.numpy()[:,1],'.r', label='Fake')\n",
    "        plt.xlim(round(min(data[:,0])) - 1, round(max(data[:,0])) + 1)\n",
    "        plt.ylim(round(min(data[:,1])) - 1, round(max(data[:,1])) + 1)\n",
    "\n",
    "        plt.text(round(min(data[:,0])) - 1, round(max(data[:,1])) + 1.5, 'D accuracy=%.2f (0.5 for D to converge)' % prob_real.data.numpy().mean(), fontdict={'size': 13})\n",
    "        plt.legend(loc='upper right', fontsize=10);\n",
    "\n",
    "        # visualize the background\n",
    "        x = np.arange(round(min(data[:,0])) - 1, round(max(data[:,0])) + 1, 0.1)\n",
    "        y = np.arange(round(min(data[:,1])) - 1, round(max(data[:,1])) + 1, 0.1)\n",
    "        xx, yy = np.meshgrid(x, y, sparse=True)\n",
    "\n",
    "        a = xx + yy - yy\n",
    "        b = xx + yy - xx\n",
    "        mem = np.concatenate([a.reshape([a.shape[0],a.shape[1],1]),b.reshape([a.shape[0],a.shape[1],1])],axis=2)\n",
    "\n",
    "        z = D(torch.from_numpy(mem).float()).data.numpy()\n",
    "        h = plt.contourf(x,y,z.reshape(z.shape[0],z.shape[1]), 20, cmap='Spectral')\n",
    "        plt.colorbar();\n",
    "        plt.show()\n",
    "        \n",
    "plt.ioff()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# visualize the background\n",
    "x = np.arange(round(min(data[:,0])) - 10, round(max(data[:,0])) + 10, 0.1)\n",
    "y = np.arange(round(min(data[:,1])) - 10, round(max(data[:,1])) + 10, 0.1)\n",
    "xx, yy = np.meshgrid(x, y, sparse=True)\n",
    "\n",
    "a = xx + yy - yy\n",
    "b = xx + yy - xx\n",
    "mem = np.concatenate([a.reshape([a.shape[0],a.shape[1],1]),b.reshape([a.shape[0],a.shape[1],1])],axis=2)\n",
    "\n",
    "z = D(torch.from_numpy(mem).float()).data.numpy()\n",
    "plt.plot(data[:,0],data[:,1],'.b', label='Real')\n",
    "h = plt.contourf(x,y,z.reshape(z.shape[0],z.shape[1]), 20, cmap='Spectral')\n",
    "plt.colorbar();\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}