ex3_combinatorial_auction_valuations.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "import lsp \n",
    "\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Combinatorial Auctions Example\n",
    "\n",
    "In this example we can use the identity $\\hat{s}_B = (-1)^{|B|}(\\mbox{DSFT}' s)_B$ in order to utilize our existing implementations of fast Fourier transform for set functions [1].\n",
    "\n",
    "[1] Pueschel, M.; A Discrete Signal Processing Framework for Set Functions; https://acl.inf.ethz.ch/publications/icassp18-setfcts.pdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cartesian_product(*arrays):\n",
    "    la = len(arrays)\n",
    "    dtype = np.result_type(*arrays)\n",
    "    arr = np.empty([len(a) for a in arrays] + [la], dtype=dtype)\n",
    "    for i, a in enumerate(np.ix_(*arrays)):\n",
    "        arr[..., i] = a\n",
    "    return arr \n",
    "\n",
    "def submatrix_liF(rows, columns):\n",
    "    prod = cartesian_product(rows, columns)\n",
    "    intersections = prod[:, :, 0] & prod[:, :, 1]\n",
    "    subsets = (prod[:, :, 0] & prod[:, :, 1]) == prod[:, :, 1]\n",
    "    cards = []\n",
    "    for row in intersections:\n",
    "        cards.append([lsp.pypopcount(entry) for entry in row])\n",
    "    cards = np.asarray(cards)\n",
    "    return (subsets).astype(np.float64)\n",
    "\n",
    "def lattice_spectrum(s, N):\n",
    "    popc = lsp.popcount(np.arange(2**N).astype(np.uint32)).astype(np.int32)\n",
    "    s_hat = lsp.fdsft3(s)\n",
    "    return s_hat * (-1)**popc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.load('./data/GSVM.npy')\n",
    "data = data[0]\n",
    "N = int(np.log2(data.shape[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### For an example valuation function, sort subsets (bundles) by their cardinality and plot the valuation and its spectrum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "popc = lsp.popcount(np.arange(2**N))\n",
    "s = data[0]\n",
    "s_hat = lattice_spectrum(s, N)\n",
    "example = []\n",
    "example_hat = []\n",
    "sizes = []\n",
    "sets = []\n",
    "for size in range(N+1):\n",
    "    example += [s[popc == size]]\n",
    "    example_hat += [s_hat[popc == size]]\n",
    "    sets += [np.where(popc == size)[0]]\n",
    "    sizes += (popc == size).sum()*[size]\n",
    "example = np.concatenate(example)\n",
    "example_hat = np.concatenate(example_hat)\n",
    "sets = np.concatenate(sets)\n",
    "sizes = np.asarray(sizes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deZgU1b3/8fd3FsAFZUcUcFBBxWgQJ4i7xhW8CZqYG0yiRGPITUxM8ktyM5qNmwSjZjHxalQUE8yNcYkxEkEUERUUIQMiCogMiIIgi8jiMsAw5/dHV/dU93RP93R39VLzeT3PPNN1ajs11fPt06fOYs45REQkXCqKnQEREck/BXcRkRBScBcRCSEFdxGREFJwFxEJoapiZwCgV69erqamptjZEBEpKwsXLtzinOudbF1JBPeamhrq6+uLnQ0RkbJiZm+mWqdqGRGREFJwFxEJIQV3EZEQUnAXEQkhBXcRkRBScBcRCSEFdxGREFJwFxEpkt8/9TpzVm4O5Nhpg7uZdTGzBWb2spktNbP/8dIHmdl8M1tpZg+YWScvvbO33OCtrwkk5yIiZe6Ps1fxwqp3Azl2JiX3XcAnnXMfB4YBF5jZSOBG4Gbn3GDgPeAr3vZfAd5zzh0B3OxtJyIiBZQ2uLuI973Fau/HAZ8E/u6lTwEu8l6P8Zbx1p9tZpa3HIuISFoZ1bmbWaWZLQY2ATOBVcA251yTt8k64BDv9SHAWgBv/XagZ5JjjjezejOr37w5mDonEZGOKqPg7pzb65wbBvQHRgBHJ9vM+52slN5qolbn3CTnXK1zrrZ376SDmomISJba1VrGObcNeAYYCXQzs+iokv2B9d7rdcAAAG/9gcDWfGRWREQyk0lrmd5m1s17vQ9wDrAcmA1c4m02DnjUez3VW8Zb/7RzrlXJXUREgpPJeO79gClmVknkw+BB59xjZrYMuN/Mfgm8BEz2tp8M/MXMGoiU2McGkG8REWlD2uDunFsCHJ8kfTWR+vfE9Ebgc3nJnYiIZEU9VEVEQkjBXUQkhBTcRURCSMFdRCSEFNxFREJIwV1EpEhc6877eaPgLiJSREGNqqjgLiISQgruIiJF4Jxjz17HW1s/DOT4Cu4iIkUw9eXIWIuPLdkQyPEV3EVEimBHY1P6jXKg4C4iEkIK7iIiIaTgLiISQgruIiIhpOAuIlIEQXVeilJwFxEJIQV3EZEi+PE/Xw30+AruIiIFtvDN9wI/h4K7iEiBNe7ZG/g5FNxFRAos6IepoOAuIlJwZsGHdwV3EZECK0BsTx/czWyAmc02s+VmttTMvu2lTzCzt81ssfcz2rfPtWbWYGYrzOz8IC9ARKTcVBQguldlsE0T8D3n3CIz6wosNLOZ3rqbnXO/8W9sZkOBscAxwMHAU2Y2xDkX/BMEEZEyUBIld+fcBufcIu/1TmA5cEgbu4wB7nfO7XLOvQE0ACPykVkRkTCoKEBwz6TkHmNmNcDxwHzgFOCbZnY5UE+kdP8ekcD/om+3dST5MDCz8cB4gIEDB2aRdRGR8lG/ZiuX3DGvYOfL+IGqme0PPAx8xzm3A7gdOBwYBmwAfhvdNMnurab4ds5Ncs7VOudqe/fu3e6Mi4iUk0IGdsgwuJtZNZHA/lfn3D8AnHMbnXN7nXPNwF20VL2sAwb4du8PrM9flkVEJJ1MWssYMBlY7pz7nS+9n2+zi4HoQAlTgbFm1tnMBgGDgQX5y7KIiKSTSZ37KcBlwCtmtthLuw641MyGEalyWQN8DcA5t9TMHgSWEWlpc7VayoiIFFba4O6cm0vyevTpbewzEZiYQ75ERCQH6qEqIhJCCu4iIiGk4C4iEkLt6sQkIiLxJj23iuunv9bmNj84/8gC5aaFSu4iIjlIF9gBfv3EigLkJJ6Cu4hICCm4i4iEkIK7iEgIKbiLiISQgruISAgpuIuIhJCCu4hICKkTk4h0CM45Bl2bcrzDNl15yiDuef6NuLTXfnEBU15Yk4ecBUPBXUQ6hPsWvJX1vomBHeBLd8+n/s33cslSoFQtIyIdwge7mvJ6vG0f7cnr8fJNwV1EOgTXaibn3FQkm+WihCi4i4hkwZLOYVQ6FNxFpEPIc8EdK+3YruAuIhJGCu4i0iHkv869tIvuCu4iIlmoKPHoqXbuIhKIo37yOI17mpOu++8LjuSmGSs4cJ9qtntNCt/41Wg+detcXn17BwCdqirY3RTZv1NlBbv3NvPYt05lb7NjzG3Ps3/nKt73mje+/stR/HbmCu58djU/vOAorjptEIN/9HjsuAA3zkg/qUZ7lPoDVQV3EQlEqsAOcNOMyMxE231txec2bIkFdiAW2AF27428/o//nRtLe9/Xbv2av73EjKXvAJEgvnT99ti65Rt2Ul1Z2oE4CAruIlIS9uxN/WGQzvsJHZR2NLYsNzU3Y5b/OpSyb+duZgPMbLaZLTezpWb2bS+9h5nNNLOV3u/uXrqZ2S1m1mBmS8xseNAXISLlz3J4QJm4ayHibi75LYRMPs6agO85544GRgJXm9lQoA6Y5ZwbDMzylgFGAYO9n/HA7XnPtYiETi6hsq1Am+9WMlFlX3J3zm1wzi3yXu8ElgOHAGOAKd5mU4CLvNdjgHtdxItANzPrl/eci0io5FISTgy0ifE8iEL2ore25f+gedSuiigzqwGOB+YDfZ1zGyDyAQD08TY7BFjr222dl5Z4rPFmVm9m9Zs3b25/zkUkVHIpCSfuWuKF6oLIOLib2f7Aw8B3nHM72to0SVqrL0bOuUnOuVrnXG3v3r0zzYaIhFQuTQtLvf67GDIK7mZWTSSw/9U59w8veWO0usX7vclLXwcM8O3eH1ifn+yKSFjlUnJva9+AqtxLXiatZQyYDCx3zv3Ot2oqMM57PQ541Jd+uddqZiSwPVp9IyJt2/L+LmrqpvH9h16mpm5a3E/jnr3cOOM1auqm8avHl8fSgdjr5Rt28MaWD+LWNe7ZS03dNH7x2LJY+nsf7GbitMjyF+56EedcbF20SaL/GP87ayU1ddOYOG0ZAIvXbqOmbhpL1kXqnevXbKWmblqsffm8Ve+2/+JzKny3vXOpdzgKQiYl91OAy4BPmtli72c0cANwrpmtBM71lgGmA6uBBuAu4Bv5z7ZIOP3a69zz94XrWq07/uczuf2ZVQDc+ezqWPrmnbtir0f9YQ5n/eaZuP2+/n8LAZg8t2U2oeN/MZO75kSWX1j1Li+tbXk4+L9PN8Ttv3FHI7+d+TpAbJ+LbnsegE/fGvl9yR3zALjwlkgno0vvejH9xSbIJQCXesuVYkjbick5N5fUH4tnJ9neAVfnmC+RDundD3anXPfRnr1J012atn47G9PPQNS0t+UYuxLOs7e5MBUbuVXLtNUUsmNWzJT40DciklZISq157cQUkr9JLhTcRcpcPuqT2wqGhSr3bvsw9beWdFZv/qDNdR0x2Cu4i5SU9ofSsASuuQ1bst53xcadcctb3m95DvF/899kw/bGrI9drhTcRcpcutieyceFpVwoXJ31vfPezNux/KNLvrZhJ+PuWZC3Y5cLBXeRElIKz/7C1mww1YPosFNwFykhzVlE93z0zvQfwiWU9UvhA0faT8Fdyk60c82Oxj1c87eXqKmbxlVT6mlubumIE22+5++IM2HqUmrqpvGbJyJtyeeu3EJN3TSWrY98hZ+9YhM1ddN43au/fe71zTy1bCMA331gMTV107jz2Ug780cXv01N3TTeevdDAB5euI6aummsey+yPHPZRuasjIyZFO0gdMavZwOws3EPv3tyBU1eZ6Hpr2zgxdVZdPrxhKucLfmi4C5l67gJTzL15cjIFk8t38hzK1sGoHuofm3cth/ubuLPL6wB4NbZkU46X5o8H4DRt8wB4Io//RuA825+DoDL71nAVffWA/DIS28D8KvHI1O1ffv+xQCc7gXs7z30MgCn3hhZ/uq99Vw2OVLP+8aWSEuON70PgptmrOCWpxv415JI3r/x10WMnRTp9JNNITk/D1TN90ofF2Gg4C6h4e+IkzirT4H64SSVWNUSrQPe09Q6U9nkM9/BOPHDQtUy5UnBXSRgicEx3+Xiturp9za7jHqY+gP6ojffi1v3xrup25BL6VJwl9BosyNOEYufqc6c+OASssvnsg2pR+Beu/XD2DOFtuzxTUY9/42t7GpqaWGS2IwwMY+Jy83F/JokMQru0iEUM9y0Krm38SE0Z2X7O/J88e75Kdd9+U8L2J3BxNOfnxQ/0NeRP56RcttB105vc/mw6+KXpTgU3CWcyqDbZiG+TKzxHuJKx6PgLqGUGNqL+VAwsfol+gBUlRcSJAV3CY02A3gxg3s7qmVE8kXBvQMZ/KPp1NRN44ml7wAw6blV1NRN4+1tHwFw9V8XUVM3jQ3bI8uPLVnPXK8OONoR57ybnwXgvQ9285snVsRaYjzy0jrmJ+mIE+1E1NzscM5x2+wG1m6NVBW8sm47981/q9U+M5dt5OnXNub56osn1YeOmhhKkBTcO5A9Xjvwr/0lMjPP9dMjHXJOueFpAKa9EpkN8aRfRZa/ed9LsY4+87zA/frG9wH42dSl3Dq7gdmvRabO/e4DL7d6KOc3p2EL7+xo5NdPrODLf4q0vvjUrXO57pFXWm371XvrufLP9TlcaZK22iVUCRLNWynlScJHwV0ykthWutHriNPUnL4lBkSax0UP8dHu4Ady8k89B/BBDufMpWlfsmaC0T/ZwwlT6TVl0KpFJFMK7pKRxF6QudQb51JezbQd+O+fWhm3fLE352c2Hn/1naz3XbZhR1wJfdorG3jAGxph0Vvb4rb1z3Eqkqu0c6iKQOpgnnG9seWnZ+YxP3si5broODBR0QHDADYllOT969ItX33foqz3fWbFZn7tDVQGkcHJUm0bHbdGJB9Ucpes5NKcL5cHiR8WoEonn+7/d/wD4+gAZCJBU3CXrLS3Wsay2CcMKhIuWi1kpFAU3CUjqeJypsEqHxNKlKOOedVSCtIGdzO7x8w2mdmrvrQJZva2mS32fkb71l1rZg1mtsLMzg8q41JgCVEql+Z8HakJYEf9UJPiy6Tk/mfggiTpNzvnhnk/0wHMbCgwFjjG2+ePZlaZr8yWG/8sQI179jJx2jI+3N0EwJNL34l1JvL7y7w1LFm3rVW634ur3+XvCc3oCi2bMcQ74iQQiu1SLGmDu3PuOWBrhscbA9zvnNvlnHsDaABG5JC/spXYZO8v897krjlvcMczkWnaxv9lYawzkd9PHl3Kp29tu9ne2Ekv8n1v5p9CSQzM0dK3f5Z5iFz33xeuazVZxvuNTbHXG3fEt1zZ1bSXuSu3xHqu+k19eT3v72ripbfeY3kbQ9uWqsQ6d5FCyaXO/ZtmtsSrtunupR0C+Oc3W+eltWJm482s3szqN2/enGyTUNnj9VzZVYYdVdZ7wxNEvfr2dqa/EvnWcYc3p2jUfQve4vsPvRz7EIu6+r5FKUuxtb94ii9Nns+Zv3kmLv3+BW9xzd9e4tp/vMLFf3yBUX+Yk9uFFIFCuxRLtsH9duBwYBiwAfitl57svZy0gtU5N8k5V+ucq+3du3eW2ShdrWffiVVSF0V72mYnLs9avpFL72oZWuCfCc35/Nv+6JHIo5l3P9jd6pgnXj8r6T47d0VK9Xu9Ca6j6v4RGZrgX948qeVIBXcplqyCu3Nuo3Nur3OuGbiLlqqXdcAA36b9gfL9z8yjcv4nv29B/GTTd6snZcY64nMGKQ1ZBXcz6+dbvBiItqSZCow1s85mNggYDCxI3L8j6zjtRESkmNIOP2BmfwPOBHqZ2TrgZ8CZZjaMSKxaA3wNwDm31MweBJYBTcDVzrny6lKYJ4lBPFp+K+ZcntkqxzyXio7U7FNKS9rg7py7NEny5Da2nwhMzCVTYRRrF67/dRHxOeqgroEcVz1UC0R1rx2T7rukU1UZzHtEwZ3IjEMzl+V35p9UVRmvvbOzVdritdvY1RRfe+Wco37NVm+fHWz/aE/Scyxbv4OdjXvYtKORNVs+iK2rX7M1p3HI48+Tl8OISBJB/X8puBOZceir9+Y2809b7nx2Fe/saARgbsOWuHUzXt3ARbc9z4Spy+LSp7ywhkvumMes5Ru54Pdz+Pyd81od9w+zVjL6ljmMu2cBI66fFWsn/sKqLVxyxzzueG5Vq32ksFZsbP1hLuIXVHDXeO4B2LijMa5N94yl7/CSb2IGf1vu//q/yFjhy9Zvj0uf8K9IsI/2Yo2W+P3bRCekSJz04Qt3RabGa/CmxMuVHgqKlB+V3AOQWPXyyrrtWR+rKU9VKyLSsSi4F0BzmVdaZ5P96ByrItK2oKKDgnsBFKvwna/TZnOcN3wPd0UktV77dwrkuAruAVDjNxHJ1AH7VAdyXAV3SSubHqqJQ/6KSApqClkYm3Y00tzs2NG4hw+80Qqdc2z0mjL6JUsLo1Wb21/FktiCR0QKS8Hd54WGLYy4fhZ3Prea4yY8yYiJTwHwYP1aTrx+VtwMSbNXbOLE62cxa3nrzk/lPAKkiIRDhw/u/nbjX7g70j78xhmvAfDB7kiLjx8+HBlXfKWv3fgVf/o3AC+vzVMJNYBPBA34JdJxqRNTjvIVPht3p286+Ozr6Wes+uMzDbHX/1y8nurKCh4q8nyrIpLazz41NJDjdviSexCyKTBn0k193D3ph8a/acaKuGUFdpHS1rtr50COq+DeDslitmo+RKQUKbjnSOOuiEguLKAWGAruAVC4F5FiU3DPkaplRKQUhTa4/+Chl+OaORZSJg8+RUSCFMqmkJt37spbK5EXVrVMrnH3nNV8tLuJCz7WL5b2p+fXcHC3fTj2kANZun6Hut2LSEkIZXBftmFH3o4VnfgCIuO0/+TRpTy86O1Y2kd79vLjf74aWz663wF5O7eISLZCWy0TpC3v70q5rkkldxEpAQruWVj33kcp1+1sbCpgTkREklNwz7N3OshIkSJS2tIGdzO7x8w2mdmrvrQeZjbTzFZ6v7t76WZmt5hZg5ktMbPhQWZeRESSy6Tk/mfggoS0OmCWc24wMMtbBhgFDPZ+xgO35yebwbhl1spiZ0FEJBBpW8s4554zs5qE5DHAmd7rKcAzwA+99HtdZKzZF82sm5n1c85tyFeG26t+zVYa9zQzsMe+7HWO/TpVsvXD3Rx10AH8bubrxcqWiEigsm0K2TcasJ1zG8ysj5d+CLDWt906L61VcDez8URK9wwcODDLbCS37cPdsdeX3DEvbl1VhdHU7Fhzw4V5PaeISCnJ9wPVZCPgJO2g75yb5Jyrdc7V9u7dO6+Z+LCNsdGbmjVegIiEX7bBfaOZ9QPwfm/y0tcBA3zb9QfWZ5+94Kg9uoiEWbbBfSowzns9DnjUl36512pmJLC9GPXtmQzmNW/1u8FnRESkSNLWuZvZ34g8PO1lZuuAnwE3AA+a2VeAt4DPeZtPB0YDDcCHwBUB5Dmtpub0pXKN5igiYZZJa5lLU6w6O8m2Drg610zl6qePLi12FkREiioUPVSdc2xqZ89QFdxFJMxCEdwnz32DEdfPomFT+kmmRURKxZz/PiuwY4ciuD/fEBlz/c13P8x4H02oISLFNqDHvoEdOxTBfc/eSCXLwjffK3JORERKQyiC++sbI9Uxf3xmVZFzIiKS2uRxtfTu2rkg5wpFcLdk/WJFRErM2Uf35SunDirIuUIR3CsU3UWkTFQUKFyFIrj7/1aL124rWj5EpGP6f+cOyXhbSzoEV/6FI7j7Su4X3fZ8EXMiIh3RNWcPznjbQlU0hCK4i4iUCytQdA9FcK8IxVWISEegOvd2WLv1o2JnQUQkI4Vq/hGK4C4iUi4qClR0L/vgvv2jPcXOgohIxlRyz9DMZRuLnQURkczpgaqIiGSr7IP7+42qlhERSVT2wX3Cv5YVOwsiIhlTnbuIiGRNwV1EpIA0/ICISAhp4DARkTJw42ePbdf2nxl+SEA5iafgLiKSg89/YmC7tu9SXRlQTuLlFNzNbI2ZvWJmi82s3kvrYWYzzWyl97t7frIqIpKbQs2CVAryUXI/yzk3zDlX6y3XAbOcc4OBWd6yiEjRnTq4V7GzUDBBVMuMAaZ4r6cAFwVwDgAaNu0M6tAiEkLD+ndLuW5gj30LmJPg5RrcHfCkmS00s/FeWl/n3AYA73efZDua2Xgzqzez+s2bN2d18ne278pqPxHpmLrv1ynluv+s7V/AnMDlJx0a6PGrctz/FOfcejPrA8w0s9cy3dE5NwmYBFBbW+tyzIeIdCBdO1exc1dTXo/pChiF1txwYeDnyKnk7pxb7/3eBDwCjAA2mlk/AO/3plwzmUqhOgOISPuNGNQjq/0O671f2m2u/uQRWR07SP/v3CGM/cQAbvzssZx1ZG9G1PRg304tLWNqeha22ifrkruZ7QdUOOd2eq/PA34OTAXGATd4vx/NR0ZFJHuHdNuHt7cVdsaycit7VVXmVkvtnyTb3zyypm4aAMccfGBOx2+vXKpl+gKPeJO9VgH3OedmmNm/gQfN7CvAW8Dncs9mcuX25hEplrOP7sO9894s6DkrsvxqncleuVShXH/xsfTp2pnnV23h7fc+4klvTogvjRzI5p27qDC4e+4bAPx8zDH89NGl2Z+siLIO7s651cDHk6S/C5ydS6ZEJL+KURA6uNs+We3Xv/u+rNr8QZ5z0+ILJ0ZK1ecM7cvarR/GgnuX6kp++qmhQEtwHz6wfLvplHUP1camvcXOgkhZsCI8oPrFRcfwh7HD+PLJNbG06declna/Wy49nurKtvNbVaB5SMtZWQf3++avLXYWRApu/87t/8JdjMYH+3aqYsywQ7j4+JaxVIYefEDa/Q7cp5oLPtavzW2a89S0JcyNMso6uO/QLEzSAQ0/tP1VBYUaibBQ1HY6vbIO7rv2qFpGOh6XRam1mCXUIAJxIduk54sr8EdSWQf3veV4h0WKoJjl9mw+jNLtk69qmTAr6+DetFc3WCQTFWX2ALJQ/9nFeNBcKGUd3Pc2K7hLx5NNobWoJfcgjlmgkns5x/6yDu7d9009CJBIWGVVd1vMOvcA4nC5HLOYyjq4H9ytS7GzIFIWittaJv9RM5BvAyFrg1PWwT1ct0IkM1lVy5Rb9UKaa1TJPb2yDu6qcpeOqNzq3IOQr1J2If8uJxya3SiZ2Srr4F6ohyoi5a6o7dzLpJRdmaRFUe/9O+ft+FeeUpO3Y2WivIN7sTMgUgTZlFqLWecexP9pJkMwnDGkd9zyhcelHtKgU2UF1UmG/O1zQMtzvZs+exwAp/nmYX36e2ekzUdUoZtdlnVwV3SXjiiTUuvR/eLHcBlyUNeAcpO5oxLy8BlvzJmrzzo8lrb4p+fGbXPjZ49tdZyvnjaI/bzg/umPHxxLX3Bd/GC0d152QtzyrZce3+pYB+5TDcDXzzy81bpLR0RGjxznTYc3/NDI/Kv/WTsgts1hvfdvtV+iRT85l4U/PiftdvlW1sFdvdQknV55/FqdC38QytRNl0RKiv4S55IJ58Ve+yeHeO0XF8Tte3S/lkB652UncIQXhAb1apnl6JFvnBy3z6rrR8ctr/YtL/pJfNBNXE4Myv71Q/pE8vKD848E4Msn13Bk3658cWQkaH7uhJZg2c1r3nyZF1DPGNKHqgrj52OOiW1TW9MjVnq+wlfV0eeALuxT3TLzUSdfSXz4wG5JS877da5i5cRRfOecwXHpKyeOYuJFHwNgwqePoWHiKI7o05XXfzmKT7XzXvbYrxM9i/A+LOvgrtheWG19rS2kU47oGXudWCJ641fxAerfP2opzb3+y1Fx61b8Mj4gJq73L6+cmHpdsvWzfF/Xxww7mC+eOJBUJnhjiENkcogpV47gwuP6MepjBzFiUA9+eP5RsfUHdKnmf8YcwylH9OSq0wbF0rv4gtqXRg7kutFHx5Zreu7H4L77c8aQ3vxh7LBY+kEHduHBr50UW06sc/b3aj1wn2rOP6ZvbLnHfp04ok9LqbVbQp+THr6JqA/ct5o1N1zI2UdH9p/w6WN44runc8Kh3Vlzw4XU9Go9rd7Iw3qy5oYLOejALjRcP5rLT6rhdK+apXNVBQN67MuaGy7k+ITx1h//dmRI4UN77huX/7aqRKorK1qtr66siO1vZrFZmjpVlU/ILJ+cJlHq7VLPObpPsbMQJ3FS3nTLM77TMvb2f51xOJ+P+zoa/w/5vXOHZHROvye/e3rS9HuvHJFyH4CbP98SoLpUVzL7+2fGlhP/Sf3Lnaoq+MzwluFnO1dVcvLhPePW+/mXqysr4gJWp6qKuIeUifW1h/kCVt8DuvCJmtQtJS7yDYlbVVHBGUN6c9sXhtO1SzUPfu0kBibMvXnUQQfw16tGckCX6rj0754TuQc99u3U6htLdWUFU64cwXH9u3GQrx453Tyn/mu84GMHxa2LVleUIhX8yj24l/gN/O3nhqVc5w9I7ZVYOg3KUQfF19ue5AuEiS72Bc22nJjBpMkf79+tzfV9unaJm3h4YI/MJx5OnMfylCN6pdiytZMTrn/0sam/ySR+yLQ1tot/7s58PHML6t+i1P/foAzb8weorIN7vurcn/3BmXk5jt/Iw3pw4L7VKdcPSvJVFOCykelLQ7k8da9NGAt8cJ/0D4Qg/bekTPPUc/8QDRmR4dsvXZNd//pcYlMQgc1/yHII7tKirIN7vgYOO7Rn8kCbi2ybnpXZ4H3tlkmAaG91W7H+ZEFUCwYRoEu9+jJT6tfSPuUd3HWvy1ox/ldzKxnH752v/PuPW6ozJjlXOh8R2XxzLc2/arDKO7g3Nxc7CylVlPVfNjl/MEv8Z8nnP097g2ax6lk7QkHSH0hLpeTcEQN1Nso6BJVwbA+9xH/zYv7bt6ckl8985qss69r61CwhpRHaJVNlHdw1zV77JcbBjONiSP7UuZQ+E/9UQQxcV6qxPSS3v0MJLLib2QVmtsLMGsysLohzNJfwsJClWncalIw/I3y3LFXJt1TvauIHYaafE+m2i6tzL7G2fHG5KdUbk4ES+7MWRCDB3cwqgduAUcBQ4FIzG9r2Xu1XyiX3jvhm6nhK9/0nElTJfQTQ4Jxb7ZzbDdwPjMn3SSpLOIJ2rqpMv1X2MrMAAAchSURBVFESQXdv9ndTB+LG4mhLVWX837pLwvVVZHgvOlf7OuykKO9ncqRk+a6uTL9nVUJb08Tl9pwz03tcneae+nOQyTWkPI/XGaoq4Wl+4t+5i3cPMvl2uY/XWcy51sMTVBa41UBn7++YbGjeqOj7sEt1fN6y/X8sZ0HdnUOAtb7ldV5ajJmNN7N6M6vfvHlzVif54aij0m+UpUwGB2orLtz6hfgR6B7++sk88Z3k3e2jAyp17VLF973Xfh8fkLrHZrT7eGLvyUS/+syxXr6GA3DRsMj1TfG6+ke7kj8wfmTcfid4nZ6+e86QuA+eB//rpNg1du1SRd8DWrq73315LdAyzse5QyNjilw3+ihu9IZNBRjSd/9YZ67oIFkA3fdr3dEpmo/bvzg8dv66UUfFRgf88YVH89i3ToudB2DetZ8EIsMo/OQ/Il8cL/XGeIkODzDu5Bqg5W8894dnxZ330atP4ZfeAFI/9o4RHdvmBm/EwlFet/zo9UZFt4sOCxA157/P4s9XfCK2vJ9v+NoLU/R6nXbNqXFj0ETd8+Va7vhS5G9yxSk1fO30wxh/+mFx2wzpG99R7d4rT+QH5x8Zu2cTPjWUadecCsDEiyPXGu3Y98g3Tua60UfRqaqCT3vvmZ7e/bnkhP5AS4e8xFEZ2+O+q07k959P3aMb4IbPHsfVZx3OSYfFv9cf+9apsYHF+nffh++dO4TJ4yJ/3+hwE4kjRObq8W+fxk//I++VEXllQTRvMrPPAec7567yli8DRjjnvpVs+9raWldfX5/3fIiIhJmZLXTO1SZbF1TJfR0wwLfcH1gf0LlERCRBUMH938BgMxtkZp2AscDUgM4lIiIJ0s9VlQXnXJOZfRN4AqgE7nHOLQ3iXCIi0logwR3AOTcdmB7U8UVEJLWy7qEqIiLJKbiLiISQgruISAgpuIuIhFAgnZjanQmzzcCbWe7eC9iSx+yUqo5yndBxrlXXGS7FuM5DnXO9k60oieCeCzOrT9VDK0w6ynVCx7lWXWe4lNp1qlpGRCSEFNxFREIoDMF9UrEzUCAd5Tqh41yrrjNcSuo6y77OXUREWgtDyV1ERBIouIuIhFBZB/dCTMIdBDNbY2avmNliM6v30nqY2UwzW+n97u6lm5nd4l3jEjMb7jvOOG/7lWY2zpd+gnf8Bm/fgsxHaGb3mNkmM3vVlxb4daU6R4Gvc4KZve3d08VmNtq37lovzyvM7HxfetL3rzdU9nzveh7whs3GzDp7yw3e+pqAr3OAmc02s+VmttTMvu2lh+qetnGd5X1PnXNl+UNkKOFVwGFAJ+BlYGix85Vh3tcAvRLSbgLqvNd1wI3e69HA40Sm2hwJzPfSewCrvd/dvdfdvXULgJO8fR4HRhXouk4HhgOvFvK6Up2jwNc5Afh+km2Heu/NzsAg7z1b2db7F3gQGOu9vgP4uvf6G8Ad3uuxwAMBX2c/YLj3uivwunc9obqnbVxnWd/TwP/hA7whJwFP+JavBa4tdr4yzPsaWgf3FUA/35tthff6TuDSxO2AS4E7fel3emn9gNd86XHbFeDaaogPeoFfV6pzFPg6UwWCuPclkTkOTkr1/vWC3BagKvF9Ht3Xe13lbWcFvLePAueG9Z4muc6yvqflXC2TdhLuEuaAJ81soZmN99L6Ouc2AHi/+3jpqa6zrfR1SdKLpRDXleochfZNrzriHl81QnuvsyewzTnXlJAedyxv/XZv+8B51QXHA/MJ8T1NuE4o43tazsE9WT1yubTrPMU5NxwYBVxtZqe3sW2q62xveqkJ23XdDhwODAM2AL/10vN5nUX5G5jZ/sDDwHecczva2jRJWtnc0yTXWdb3tJyDe9lOwu2cW+/93gQ8AowANppZPwDv9yZv81TX2VZ6/yTpxVKI60p1joJxzm10zu11zjUDdxG5p9D+69wCdDOzqoT0uGN56w8Etub/alqYWTWRgPdX59w/vOTQ3dNk11nu97Scg3tZTsJtZvuZWdfoa+A84FUieY+2IhhHpN4PL/1yryXCSGC79zX1CeA8M+vufV08j0g93gZgp5mN9FoeXO47VjEU4rpSnaNgooHIczGRewqRvI31WkUMAgYTeYiY9P3rIpWvs4FLvP0T/2bR67wEeNrbPqhrMmAysNw59zvfqlDd01TXWfb3tFAPKQJ68DGayJPtVcCPip2fDPN8GJGn6C8DS6P5JlLPNgtY6f3u4aUbcJt3ja8Atb5jXQk0eD9X+NJribwRVwG3UqCHbsDfiHx93UOkRPKVQlxXqnMU+Dr/4l3HEiL/sP182//Iy/MKfC2XUr1/vffIAu/6HwI6e+ldvOUGb/1hAV/nqUSqCJYAi72f0WG7p21cZ1nfUw0/ICISQuVcLSMiIikouIuIhJCCu4hICCm4i4iEkIK7iEgIKbiLiISQgruISAj9fww2kpApifh7AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAATVUlEQVR4nO3dfbBd1X3e8e+DhHBLqHnRhYgXWdghtMQphNxiGFrXTsyb6sZNhtZSOzWJ6ShxzEzcaTOBesb2pO3Ubmp3muCxrBgVO+MQNxOTaGphoLYTwgzBXDECCYOMTOUgS0bC2IIU23ClX/+4W83R1TnS1Tn3RVr9fmbOnL3XXnuvte4+erTvOvvck6pCktSukxa6A5KkuWXQS1LjDHpJapxBL0mNM+glqXGLF7oD/SxdurRWrFix0N2QpBPGpk2bnq+qsX7bjsugX7FiBRMTEwvdDUk6YST55qBtTt1IUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktS4poL+T7ft4dkXXl7obkjScaWpoP/F//4Ib/vYny10NyTpuNJU0AP8cPLAQndBko4rzQW9JOlQBr0kNc6gl6TGGfSS1Lij/pniJOuBtwN7quqNXdnngIu7KqcD36uqy/rsuwN4CdgPTFbV+Cz1W5I0QzP5e/R3ArcDnzlYUFXvPLic5KPAviPs/9aqen7YDkqSRnPUoK+qB5Ks6LctSYB/BvzM7HZLkjRbRp2j/wfAc1X19IDtBdyXZFOSNSO2JUkawqhfJbgauOsI26+uql1JzgbuT/JUVT3Qr2L3H8EagOXLl4/YLUnSQUNf0SdZDPwC8LlBdapqV/e8B7gbuOIIdddV1XhVjY+N9f1+W0nSEEaZunkb8FRV7ey3McmpSU47uAxcC2wdoT1J0hCOGvRJ7gIeAi5OsjPJzd2mVUybtklybpKN3eo5wINJHgO+Cnyhqr44e12XJM3ETO66WT2g/Bf7lO0CVnbLzwCXjtg/SdKI/GSsJDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxM/mGqfVJ9iTZ2lP2oSTfSrK5e6wcsO/1SbYl2Z7k1tnsuCRpZmZyRX8ncH2f8v9aVZd1j43TNyZZBHwcuAG4BFid5JJROitJOnZHDfqqegB4YYhjXwFsr6pnquoV4A+AdwxxHEnSCEaZo78lyePd1M4ZfbafBzzbs76zK5MkzaNhg/4TwBuAy4DdwEf71Emfshp0wCRrkkwkmdi7d++Q3ZIkTTdU0FfVc1W1v6oOAL/L1DTNdDuBC3rWzwd2HeGY66pqvKrGx8bGhumWJKmPoYI+ybKe1Z8Htvap9ghwUZILkywBVgEbhmlPkjS8xUerkOQu4C3A0iQ7gQ8Cb0lyGVNTMTuAX+7qngt8qqpWVtVkkluAe4FFwPqqemJORiFJGuioQV9Vq/sU3zGg7i5gZc/6RuCwWy8lSfPHT8ZKUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhp31KBPsj7JniRbe8p+K8lTSR5PcneS0wfsuyPJliSbk0zMZsclSTMzkyv6O4Hrp5XdD7yxqv4u8HXgtiPs/9aquqyqxofroiRpFEcN+qp6AHhhWtl9VTXZrf4FcP4c9E2SNAtmY47+3cA9A7YVcF+STUnWzEJbkqRjdNQvBz+SJO8HJoHPDqhydVXtSnI2cH+Sp7rfEPodaw2wBmD58uWjdEuS1GPoK/okNwFvB/5FVVW/OlW1q3veA9wNXDHoeFW1rqrGq2p8bGxs2G5JkqYZKuiTXA/8BvBzVfXygDqnJjnt4DJwLbC1X11J0tyZye2VdwEPARcn2ZnkZuB24DSmpmM2J1nb1T03ycZu13OAB5M8BnwV+EJVfXFORiFJGuioc/RVtbpP8R0D6u4CVnbLzwCXjtQ7SdLI/GSsJDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxMwr6JOuT7EmytafszCT3J3m6ez5jwL43dXWe7r5nVpI0j2Z6RX8ncP20sluBL1XVRcCXuvVDJDkT+CDwJqa+GPyDg/5DkCTNjRkFfVU9ALwwrfgdwKe75U8D/6TPrtcB91fVC1X1XeB+Dv8PQ5I0h0aZoz+nqnYDdM9n96lzHvBsz/rOrkySNE/m+s3Y9CmrvhWTNUkmkkzs3bt3jrslSf//GCXon0uyDKB73tOnzk7ggp7184Fd/Q5WVeuqaryqxsfGxkboliSp1yhBvwE4eBfNTcCf9KlzL3BtkjO6N2Gv7cokSfNkprdX3gU8BFycZGeSm4EPA9ckeRq4plsnyXiSTwFU1QvAvwce6R6/2ZVJkubJ4plUqqrVAzb9bJ+6E8C/6llfD6wfqneSpJH5yVhJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklq3NBBn+TiJJt7Hi8med+0Om9Jsq+nzgdG77Ik6VjM6KsE+6mqbcBlAEkWAd8C7u5T9c+r6u3DtiNJGs1sTd38LPCNqvrmLB1PkjRLZivoVwF3Ddh2VZLHktyT5CcGHSDJmiQTSSb27t07S92SJI0c9EmWAD8H/GGfzY8Cr6uqS4HfAf540HGqal1VjVfV+NjY2KjdkiR1ZuOK/gbg0ap6bvqGqnqxqv6qW94InJxk6Sy0KUmaodkI+tUMmLZJ8qNJ0i1f0bX3nVloU5I0Q0PfdQOQ5G8C1wC/3FP2KwBVtRa4EXhPkkng+8CqqqpR2pQkHZuRgr6qXgbOmla2tmf5duD2UdqQJI3GT8ZKUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekho3G18OviPJliSbk0z02Z4kv51ke5LHk1w+apuSpJkb6Rumery1qp4fsO0G4KLu8SbgE92zJGkezMfUzTuAz9SUvwBOT7JsHtqVJDE7QV/AfUk2JVnTZ/t5wLM96zu7skMkWZNkIsnE3r17Z6FbkiSYnaC/uqouZ2qK5r1J3jxte/rsU4cVVK2rqvGqGh8bG5uFbkmSYBaCvqp2dc97gLuBK6ZV2Qlc0LN+PrBr1HYlSTMzUtAnOTXJaQeXgWuBrdOqbQDe1d19cyWwr6p2j9KuJGnmRr3r5hzg7iQHj/X7VfXFJL8CUFVrgY3ASmA78DLwSyO2KUk6BiMFfVU9A1zap3xtz3IB7x2lHUnS8PxkrCQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcUMHfZILknwlyZNJnkjya33qvCXJviSbu8cHRuuuJOlYjfINU5PAv6mqR7vvjd2U5P6q+tq0en9eVW8foR1J0giGvqKvqt1V9Wi3/BLwJHDebHVMkjQ7ZmWOPskK4KeAh/tsvirJY0nuSfITRzjGmiQTSSb27t07G92SJDELQZ/kR4A/At5XVS9O2/wo8LqquhT4HeCPBx2nqtZV1XhVjY+NjY3aLUlSZ6SgT3IyUyH/2ar6/PTtVfViVf1Vt7wRODnJ0lHalCQdm1HuuglwB/BkVX1sQJ0f7eqR5Iquve8M26Yk6diNctfN1cC/BLYk2dyV/TtgOUBVrQVuBN6TZBL4PrCqqmqENiVJx2jooK+qB4Ecpc7twO3DtiFJGp2fjJWkxhn0ktQ4g16SGmfQS1LjDHpJOg78h//5Nf7ef/xfc3LsUW6vlCTNkk89+L/n7Nhe0UtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6QFNrHjhTk9vkEvSQvsxrUPzenxDXpJapwfmJKkBfDyK5NM7Pguu/d9f87bMuglaQGs+cwmHtz+/Ly0NVLQJ7ke+G/AIuBTVfXhadtPAT4D/DRTXyH4zqraMUqbknS8qiouvG3jQnfjMEMHfZJFwMeBa4CdwCNJNlTV13qq3Qx8t6p+LMkq4CPAO0fpsDTdxi27eeO5r+Vv/Y3FPPqX3+XMU0/hwrNO5eVXJ/nTbXu5+g1LOf3Uk3ll8gCnLD78banua43/er1PG9OqkD61pteZyfbpx5leZ/FJ6dkW9h8oFp0Uqor9B4rFi6bGU1UkU+WTB4qTu/JXJg+wpBvzgQPFqwcOcMriRUfc59X9B1iUcFLX9g8n97Nk0Ul9607uP8CB4pA2Jg/UIev7q/j9h/+SD2544sg/oOPEr193Mb917zaWvfY17N73AwB+4/q/zUe++BQAJy8K11xyDhu3fJsP/8JPcuvnt/DP37Sct/2ds3n3nRML2fWBMuxXuCa5CvhQVV3Xrd8GUFX/qafOvV2dh5IsBr4NjB3te2PHx8drYuLYf2Arbv3CMe8jSceTHR/+R0Ptl2RTVY332zbK1M15wLM96zuBNw2qU1WTSfYBZwGHTUwlWQOsAVi+fPlQHbry9Wfyw8kDvPmisaH214nnz76+l83Pfm/G9U9ZfBK/ft3F/299+iVHcfg1yOF1Dne048zkemr69c89W7/NE7tePPqO0lGMEvT9flGd/nKeSZ2pwqp1wDqYuqIfpkN/sOaqYXbTCexfX/PjC92FOXPLz1y00F1o2sYtu/nVzz56WPlZpy7hO//nlUPKTnvNYl76weTAY/3jS8/l6ede4qlvvzTr/ZwNowT9TuCCnvXzgV0D6uzspm5eC8ztR8AkaQZW/uSyoadJRnX7l5/mv9z39Xlrb5QPTD0CXJTkwiRLgFXAhml1NgA3dcs3Al8+2vy8JLXuvW/9sUPWL1x66py2N/QVfTfnfgtwL1O3V66vqieS/CYwUVUbgDuA30uynakr+VWz0WlJOpElOey3ibm8mWSk++iraiOwcVrZB3qWfwD801HakCSNxr91I0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4/ziEUk6Dnz63Vfw0g9enZNjG/SSdBz4hz8+d39116kbSWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuNyPH6Fa5K9wDeH3H0p8Pwsdud45Tjb4jjbshDjfF1V9f3U1XEZ9KNIMlFV4wvdj7nmONviONtyvI3TqRtJapxBL0mNazHo1y10B+aJ42yL42zLcTXO5uboJUmHavGKXpLUw6CXpMY1E/RJrk+yLcn2JLcudH9mKsmOJFuSbE4y0ZWdmeT+JE93z2d05Uny290YH09yec9xburqP53kpp7yn+6Ov73bN/M0rvVJ9iTZ2lM25+Ma1MY8j/NDSb7VndPNSVb2bLut6/O2JNf1lPd9/Sa5MMnD3Xg+l2RJV35Kt769275ijsd5QZKvJHkyyRNJfq0rb+qcHmGcJ/Y5raoT/gEsAr4BvB5YAjwGXLLQ/Zph33cAS6eV/Wfg1m75VuAj3fJK4B4gwJXAw135mcAz3fMZ3fIZ3bavAld1+9wD3DBP43ozcDmwdT7HNaiNeR7nh4B/26fuJd1r8xTgwu41u+hIr1/gfwCruuW1wHu65V8F1nbLq4DPzfE4lwGXd8unAV/vxtPUOT3COE/oczrn/+Dn49G9OO7tWb8NuG2h+zXDvu/g8KDfBizreeFt65Y/CayeXg9YDXyyp/yTXdky4Kme8kPqzcPYVnBoAM75uAa1Mc/jHBQKh7wugXu7127f128XeM8Di6e/zg/u2y0v7uplHs/tnwDXtHpO+4zzhD6nrUzdnAc827O+sys7ERRwX5JNSdZ0ZedU1W6A7vnsrnzQOI9UvrNP+UKZj3ENamO+3dJNWazvmWo41nGeBXyvqianlR9yrG77vq7+nOumFH4KeJiGz+m0ccIJfE5bCfp+884nyn2jV1fV5cANwHuTvPkIdQeN81jLjzetjesTwBuAy4DdwEe78tkc54L8DJL8CPBHwPuq6sUjVe1TdsKc0z7jPKHPaStBvxO4oGf9fGDXAvXlmFTVru55D3A3cAXwXJJlAN3znq76oHEeqfz8PuULZT7GNaiNeVNVz1XV/qo6APwuU+cUjn2czwOnJ1k8rfyQY3XbXwu8MPuj+WtJTmYq/D5bVZ/vips7p/3GeaKf01aC/hHgou7d7CVMvZGxYYH7dFRJTk1y2sFl4FpgK1N9P3g3wk1MzRPSlb+ru6PhSmBf96vsvcC1Sc7ofqW8lql5v93AS0mu7O5geFfPsRbCfIxrUBvz5mAodX6eqXMKU31b1d1dcSFwEVNvQPZ9/dbUZO1XgBu7/af/zA6O80bgy139uRpTgDuAJ6vqYz2bmjqng8Z5wp/T+XpTYx7eNFnJ1Dvk3wDev9D9mWGfX8/Uu/GPAU8c7DdT83JfAp7uns/sygN8vBvjFmC851jvBrZ3j1/qKR9n6kX5DeB25ukNO+Aupn7FfZWpK5Wb52Ncg9qY53H+XjeOx5n6x7usp/77uz5vo+cOqEGv3+418tVu/H8InNKVv6Zb395tf/0cj/PvMzWN8DiwuXusbO2cHmGcJ/Q59U8gSFLjWpm6kSQNYNBLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxv1fRotFr0psSdQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(example)\n",
    "plt.show()\n",
    "plt.plot(example_hat)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "fourier_support = np.where(popc <= 2)[0]\n",
    "errs = []\n",
    "for s in data:\n",
    "    s /= np.linalg.norm(s)\n",
    "    M = submatrix_liF(fourier_support, fourier_support)\n",
    "    s_spectrum, _, _, _ = np.linalg.lstsq(M, s[fourier_support], rcond=-1)\n",
    "    s_hat = np.zeros(2**N)\n",
    "    s_hat[fourier_support] = s_spectrum\n",
    "    s_approx = lsp.fidsft3(s_hat * (-1)**popc)\n",
    "    errs.append(np.linalg.norm(s - s_approx))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean relative of reconstruction error 4.314737525604432e-06\n",
      "standard deviation of reconstruction error 1.7453649107869637e-06\n"
     ]
    }
   ],
   "source": [
    "print(\"mean relative of reconstruction error\", np.mean(errs))\n",
    "print(\"standard deviation of reconstruction error\", np.std(errs))"
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}