usage_example.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import xarray as xr\n",
    "from scipy import signal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from xrbispectral import calculate_bispectral"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dask.diagnostics import ProgressBar"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A signal containing frequencies $f_1$, $f_2$ and $f_1+f_2$ which are all in phase is created and its complex spectrum is estimated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 2**16\n",
    "x = np.arange(N)        # time\n",
    "NFFT = 32\n",
    "w = 2*np.pi*np.fft.fftfreq(NFFT) # angular frequencies\n",
    "i, j = 3, 9                      # selected frequency indices\n",
    "y = np.sin(w[i]*x) + np.sin(w[j]*x) + np.sin(w[i+j]*x) # signal\n",
    "# the spectrum of the signal is calculated for several windows (for averaging)\n",
    "f, t, s = signal.spectrogram(y, fs=1.0, window='boxcar', nperseg=NFFT, \n",
    "                             scaling='spectrum',\n",
    "                             mode='complex', return_onesided=False, noverlap=0)\n",
    "xs = xr.DataArray(s, coords=[('f', f), ('t', t)], name='sample_spectrum')\n",
    "xs.attrs['base_fs'] = xs.f.isel(f=[i,j]).values\n",
    "xs = xs.sortby('f')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The power spectral density contains these frequencies as expected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO2de7Rj11nYf5+kqztz77WdOJ44wY9MIKbFUEpgMM+SQhNwgDpdJAHTUhJegYIJbQos80qpYS1KWDwW4EUSCjQ8TR48HDANSXi1tEk9hpDgBJOJSeKJQzy2k9iSZqQr6esfR1s6V3MknX3OPkfaW/u31qy5ko509iedfb79Pfb3iaoSiUQike2lse4BRCKRSGS9REUQiUQiW05UBJFIJLLlREUQiUQiW05UBJFIJLLltNY9gCJcccUVevLkyXUPIxKJRLzinnvueVhVT8w/76UiOHnyJKdPn173MCKRSMQrROQDWc9H11AkEolsOVERRCKRyJYTFUEkEolsOVERRCKRyJYTFUEkEolsOZUrAhG5UUTuE5EzInJrxusvFpFzIvKOyb9vrnpMkUgkEplRafqoiDSB24HnAGeBu0XkTlV999yhv62qt1Q5lkgkEolkU7VFcANwRlXvV9UBcAfwvIrP6SWHozGvvfsBxuN6y4Lf/f5Hue8fH6/1nOvif733HB94pLvuYdTCH9/7jzz02IV1D6MWfv8dH+Lj5w9rPeeZhx7nbfc/Uus5q6RqRXAV8EDq8dnJc/M8X0TeKSKvF5Frsj5IRF4iIqdF5PS5c+eqGOta+cszD/O9b3gnf/3AR2s97w/87rv4qTffV+s518V33fEOXv0X9697GJUzGI75tl+/h9/8fx9c91Aq5yOPXeC77ngHb/ybB2s978/9yRm+9/XvrPWcVVK1IpCM5+aXvG8ETqrqpwNvAV6T9UGq+mpVPaWqp06cuGiHtPc8dmF45P/aznt+yOM1n3MdqCqPnT/cCll7gyFjZStkffzC4eT/uufN4fTcIVC1IjgLpFf4VwNHVLeqPqKq/cnDXwQ+q+IxbSTd/vDI/3Wet+5zroP+cMxwrFsha2dN19I66PRHwDrmzYju5NwhULUiuBu4TkSeLiJt4GbgzvQBIvLU1MObgPdUPKaNZB2KQFXpDobTG0fImO91O2RNblDbIet6ftdOf8hgNGYwHNd63qqoNGtIVYcicgvwJqAJ/LKq3isitwGnVfVO4KUichMwBB4FXlzlmDaV7nRlU98q48LhmLFCbxDOymYRRsZtkLU7SG6KWyFr38haryLoDWbnbbfatZ67CiqvPqqqdwF3zT338tTf3wd8X9Xj2HTM5K3TIuhs0Sp5m9wlW2X9TOdNvUqvk7K6nrDnvyKIO4s3hOlNucaVTdodpVpv2mrdbNXNcYuUXmdNbrDZdxyG1RUVwYawjslrJs9YEzdRyGyTRbCuAOo6WMe8GY2V84dhxWGiItgQzMqiV+MKI+1D7tbsY62baYzgcFT7pr26Mf7r7hbECHr9+mVNxyPqjk1URVQEG8I6XBfpVVToq0fzvaoyXc2FSrR+qiXtDgrlO46KYEOYBr1qXGGklU4oJu4itknpzTJpwrd+1ulSTf4OY1ERFcGGsI7g01ETN4wLehFH3WBhy5q+hnqBWz/rWEBF11CkMrprMHHTq5nQLYLOFloE83+HiJHvwuGY4aiehIcQLemoCDaEdZi423jDgHAm7yLSq+PgZe3Xb+nFGEGkEkypB4jB4qrYJougE+CNahHr+F2PzpswXG9REWwAptRDqyH0BqPaNnd1B0NajaRAbCgX9CJ6/dFM1sBjBL3+Fv2uqWu4Ln+9WbS1GhKMoo2KYAMwq5oTl+wyHCv9mgpZdfsjrjjYnfwdxgW9iO5gyIlLtkPWTn+bZB1NZa0rg6ebmq+h7L+JimADMBfWk2uevEmdlB3azUatpS3WQac/rP37XRfdQUrWwH/X7hp+V6NwrjjYjemjEXeYyfrkS48B9aVy9gZD9ndb7O82a93RvA56/REnLkm+3+DdJf3R9FoKOVhsSj3ULWuvP2Sv3eSSY63pzmbfiYpgAzA3pidPTdz6VjaJImgFv0ru9IdcerzFsZ1G8KvktPUTsoI3MYGprDXGCMy8CUXRRkWwAcxcQ8eOPK7jvPvtJvvtcC7oRXQHQ/bbreBlHY7G9IfjaewnZFlnCyhjEdSj9Dr90WTeNINZVERFsAGYyXrlpfVO3m5/5hoK5YJexEzWsK0fc3O85FiLvXYzaFnn502tC6jptRSGxRUVwQZgTNora44RdPtDDgK7oLMYDMccjpSD3WbwshqFPv1dA06VNfPmioNdRKjNX28UwUFAi4qoCDYAY9KeqDFGkGxiG7G/2wzqgs7CyJZM3rBXyUdlDft3NfPk4Jhx+dW0gBrMFlD9YX2lLaokKoINYB3po/3hmNFY2Wu32Gtvxw1j38gasBtsKutuM3jXkLHskt+1Plm7/RF77eT7TY/DZ6Ii2AC6gyHtVoPL9naAelxDZtIcmFVy0C6EyQ1jC1bJU1nbYWW1ZNFNKb2D3VZte2GMS/VgN2n5HsLCIiqCDcBcWLutJjtNqWXyTldTqQBqqH2L06vk/d1mECu4RXTmXEMhlxefj4fUHSPYN4ogAGUbFcEGYExNoDY3zcxdkgRQ6yxtUTdpv3nobrCjsobuGkpk25vKWr3SG48nsbV2sqiAMFJ0oyLYADoTiwCS1U0tFsFgdsPYn/o6/b+gs+imYgQHu0mMIFTr5yJ3SaC/KcySLPZ26pPVNPpJ5o2xCPy3uqIi2ABMqQegtnIP6ZWjOXeobgQT/zAuhLEmFV9DZF7WUH9TSNJF99tNGg2ZyFqDIsiYNzFGEHGCKfUATHK/64sRpINeoa4ej66SwzHns+j2h4jA8Z3m9FoK1vo5soCqJ310mrKaDhYHcC1FRbABmFIPQG0lEKb+1XaTvYAu6Cw6czECCFvW/XYLEWG/3UQ1XEvvyAKqpniIWUAl8yYcl2pUBBtAr78G19AgvbKZXNCB3jB6gyHNhrDbagRlzmfR64+mQcyQslqySLJ3ZrKePxwxGldr/RydN+Za8n/eREWwAaSDxXXlfmfFCMK9YSRZHiKSMuf9n7xZdFLukpBuVFl0J9YPzGStOk6QnjfHd5qIhDFvoiJYM+lSD8A0q6VqOv0RO02h3WpMJ1OofvOjijYccz6LbkrWvdCzweZiBFC9gk+7GRP3WxiZWVERrJl0qQeobx9Bt581ify/oLPo9ofTOIiRNYTJm0V30jQF2IIkgNERlypUL+tsI2Zz+n8I8yYqgjWTLvWQ/N/kcKQMKt7cZerzw+yiDjWomFhcRxVBXU1M6qbbHx1xM0LIsg6n8S1zLVctay+1/8b8H4LrrXJFICI3ish9InJGRG5dctwLRERF5FTVY9ok0qUe0v9XvcpIuxDqLG2xDtI3jIOpG8z/yZtFlrskWFn76cVMPdZPuoAhEEztqkoVgYg0gduB5wLXA18rItdnHHcJ8FLg7VWOZxNJl3pI/q/ngu72R9P0N6jPJbUOEnfJxPW2BTGCvXb48RBT6mHvothPtUqv2x9ybKdBsyEAwZTxqNoiuAE4o6r3q+oAuAN4XsZxPwK8ArhQ8Xg2jixTM3m+4gt6MLMIwKxsAl05pmTdaTZot8LtW5y4hsJPHzWlHuqWtTsYBTlvqlYEVwEPpB6fnTw3RUSeCVyjqn+w7INE5CUiclpETp87d879SNdEOgsh+b+uoNfMrDbnDfGGASaoOLN+QjHn5xmNlfOH6U1W4abK9ubmTV0lodNJFub8ISwqqlYEkvHcdMeHiDSAnwb+86oPUtVXq+opVT114sQJh0NcL+lSD+n/q48RHHUNhXJBZ9G5aPKGWYo6vdkJoNkQju2Eaf3M++rrjK2lF1ChuFSrVgRngWtSj68GHkw9vgT4NODPROT9wOcCd25TwDhd6iH5v54LOp1bD/WVtqibw9GYwXB81PoJVNbZtXTUdRGmrEeTLPZ2jCVd/T6Co9ZlM4jvt2pFcDdwnYg8XUTawM3AneZFVf24ql6hqidV9STwNuAmVT1d8bg2hvlVXF27QdMVT6G+0hZ105u7YZi/Q0ypnM9xT/6ur2FLnczKqCeyNhrCXrtZuay9VCpycv4WFw7HlZe2qJpKFYGqDoFbgDcB7wFeq6r3ishtInJTlef2he6CGEGVFkF/OOJwpEctgkBXjp2poj16cwwxpXJ+TwpQa1P3OsmUtQb35rybMZR2la3Vh5RDVe8C7pp77uULjv2XVY9n00iXeoB68qHTFRQN+4E2dc9yl+y3mzz4sfPrGlJlZMoaaBJAZ8HvWrXSS1cKTp+/2x9y6bGdSs9dJXFn8ZqZd9HstpIc5SpdF/NWiPk7RNfQopVjmO6So4kHEG4SwHySBTDtvV0lvf68a6ie/QtVExXBmunMZSGYOvJVXljzcYnk7yaDSWA1JOaDihByAPWo3zz5O4yslnl6g/plTQpEXrz/BvzfqxEVwZrpzmXvQPU3qvlMJQh381En8+bYpDsYBde5q5MZIwgzVTbLNVR15d7zhyPGenRREUqjo6gI1sx8Pj/AXsUrm06WWR1oKep043rDXrvFaKz0g7N+JjfHmt0l66DbH3J8pzkt9QCm3EN1Sm++HAyEU+E1KoI1M29qQvUVDed3Zab/Dq0C6XwJDwjHnJ/HXDMmpx5mq+TQrJ/uXBonVL9jPDsVOYzKvVERrJn5nYqQ+OurtQguXiXXVdqibjKtn0C7lJmMlsaRVXKLscKFw/Csn/05S7pq62e+HEz6b9/nTVQEaybTNVTxtvVFQcX0a6HQ7Q9pCBzbmV3qxrT3ffLOk27AYzgIVMFnLaD220nsZ1zR5q4sN2Mo8yYqgjWT5RqqOuhlXAj7GTGC0HbcmgY8IrNVcqgNW+YrY0I4N6p55kukQOp3PazG0usNLt65bdxwvjeniYpgzcxXM4Tqi6J1+0NaDWG3Nfv5Z0Evvy/oebK/3zDM+XkWuUvA/52v8/QGo4WyVrVHJCsrq9Ew6d5+f79REayRrFIPUH25B9PX9ugqOcwmJvMlqCEdLA5L6c3vSYFwS1F3Mt1g1Sr4rI2YUH2WXx1ERbBGsko9QDJ5B8Mxh6NqAnyd/mIXQmir5PnaMDD7vn2fvPMssi7NayHR7Q+nbUcNs9+1GqWXlWQBYWxQjIpgjSxaYcxM3Kp8nRffMOoobbEOeoOsrKyQ3SX1rpLXRbdfv6xZMQLz2Pf00dxF50TkCcDXAyfT71PVl7of1naQVeoheTzJ9BgMuWzPfSGrLLO6jtIW66DTH3HVE9pHngs7gHrx5kQIKzBuSj0sjBFUJGu3P2S31aDVPLp+3gugv4VN9dG7SPoFvAsIKyl5TWSVekgeV3uj6mbcMCDMUtRZAdR2q8FOU4IMjO/NWz/THePhyHr+cIRqliVdbapslpsRkoXcRx7zu926jSI4pqovq2wkW0hWBcX046oUQW8w4oqD3YueD7FhS5YbDMKTdTzWTNdQiDGCrI1d6cdVWbVZmUrmvL67hmxiBL8mIt8iIk8VkcvNv8pGtgWsihFUGfSaVz7mvCGtHGGJrAGY82lM7vy8pddqNththdW3uNfPlrVq11BWVpYZh+/Xko1FMAB+AvgBZg3oFfhE14PaFrLykqF6EzfZgXrxyqbq0hZ1MxyNuXA4XjB5/U/5S7NoUWGeC0nWrMqjUH3hxKxKwWYcvn+/NorgZcAzVPXhqgazbSyKEexXHiO42IWQjKPFI51eJedcB90FWR4AexVv2qubRamNUP0GxbrJajYE0GwIx3YalcbWLttrX/S8cQ2Nx3qkzpNP2LiG7gXCuUtsAFmlHtKPqzBxB8Mxg9H4ohxsqL60Rd1kVR41BCdrRmVMQ2husO6K37Uq92ZSwiPbkobqSlvUgY1FMALeISJ/CvTNkzF9tDhZpR6g2nIP5uY4nz4K4a4cF90cfc/0SJPVgMcQWmB82nWuXa+sWYXu4GiWX5bryAdsRv17k38RR2SVeoCkUmZDqnENzeISGZMosJVjZ8kNIzTXUFZlTMP+bouPnz+se0iVsUzBV+mvX5Y+al6/spIzV09uRaCqr6lyINtIVrVIMJu7qnFd9Ba4o8xzprTFTtP/TedZDXgMobmGlrtLmjz4sfN1D6kyFqWPQnUZPKq6NH0UqqsEUAc2O4v/gVm20BRVjVlDBcmqDWOoKtNj2SRKX9CX7fmvCBZlZUF4mTSL9qRAYiVUVZFzHUwXMwtcQ492B87P2R+OGY11wbzxv+eDjWvoVOrvY8ALgbiPoASLTE2ozl+/zIVQdWmLulkVVDwcKf3hiN3WxTcU38hqNmQIbcf4olIPkMj6wUfd57QsXVQE0MA+97JPVR9J/fuQqv4M8CUVji14ssofGKqavMtuGFWXtqibpTGCiitV1s2i3HqYLCoGo2D6Fi/aJAhU1hugu/T79b+IoY1r6DNTDxskFsIlzke0RSwq9QATc76CC2uZCyG0pu7LYgTpwnOX71+cG+4bvcGQ4ztNmhl57Pu7LUZjpT8cc2wnDOtnuUu1Cks6ezdz8pz/PR9sXEM/mfp7CPwD8NVuh7NdLF3Z7Lb4UAUBPrNqWbqy8fiCTtPtDxGB4xk3v9BKUXcWbBKEo66LIBTBYHTRJkyDSQJQ1Yuy8cqdc0mmUgD1nHIpAhFpAK9U1d+ueDxbxaJSD2BiBFWmj2al3vkf9ErT6Y/Y22lm7vYMrTnNKjdjcsyIJx3UOapqWJavv9duoZpUKM1a7BRlqeut4tIWdZArRqCqY+A7Kh7L1tHNqBZpqGpjTK8/oiHJXoV5zOQKZfPRosqjEIY5nyarAY/hIICsljTLXENVydpb4lJtNoTjO02v541NjuCbReS7ReSaWH20PIejMYNhdqkHqK79nclUyjKbQ2vYssr1Btslq883qjSL9t9AdTn9y5IszHl9rtxrowi+kcQq+Avgnsm/06veJCI3ish9InJGRG7NeP3bRORdIvIOEfnfInK9xZi8ZdnuSEjMzQuHY4aO+xYv2iYP1Za2WAfLV47+m/NpkkKC2TepvQBcF2nMjvwsquq9vcylmpzX78q9Nk60T1HVI8VZROTYsjeISBO4HXgOcBa4W0TuVNV3pw77TVV95eT4m4CfAm60GJeXLKuMmX6+dzjiUoe7fBftjoRZaYuQVo6Lbhjmed8bihi6gyHX7u5lvhaaG2zp/puKUqB7S5IszHl9njc2d5j/k/O5NDcAZ1T1flUdAHcAz0sfoKqPpR7uk7F7OURWWgQVuS6WuRBMaYuQVo6rXAhBybroJhVAVotBVVf8rhNZHd+UO/0R7WaDdiv7llmVK7cuVloEIvIU4CrguIg8EzDO5UuB7CXIjKuAB1KPzwKfk3GO7yDpd9BmwSY1EXkJ8BKAa6+9dtWwN55lpR7Sz7uevFl9befPG8INA5a7hnZbDVoNCUjWHOmjHq9YDf3hmLGyMNuuKutnWVYWJON5pOO+tEVd5HENfRnwYuBqkr0ERhE8Bnz/ivdmJfJm1Su6HbhdRP4t8IPAizKOeTXwaoBTp055bzUsK/WQPG+yH9xe0J3+kKufuFh/h1SVs7PEby4i7FW0C7VuVJXuIE/6qP+yrvLV761xAfXBR/xt17JSEUyqjr5GRJ6vqm+w/PyzwDWpx1cDDy45/g7gFyzP4SXTmuorJq/rYmG9Bc01DCFV5VyWUglGVv+V3vnDEaqLrct2q0G72QgiCWDVAuqgosB4d7C818BBRdWC68ImRvBZIvIE80BEnigiP7riPXcD14nI00WkDdwM3Jk+QESuSz38CuC9FmPylkXt9gxVZbUsc5dAMsFCWDmOx6ZscPhusFVuxuQ1v/PcDbMF1PIYgeskgGVZWWY8PlvSNorguar6MfNAVT8KfPmyN6jqELgFeBPwHuC1qnqviNw2yRACuEVE7hWRd5DECS5yC4XIsi3r6efdB71WKALP86EN5ntbtooLpSrnsjo4hr1AkgBm8yZb1lazwW7Lfd/i1fOmOS1t4SM26aNNEdlV1T6AiBwHsiumpVDVu4C75p57eerv77IYQzDM2u0tjxG4XGUMR2P6w/FSd0koK0ezIlwUVAQjawBKb0n5A8PBVlk/7pVebzDkKZcuzpbf362mtEVd2Iz414G3isivkAR8vxGIXcsK0u0PF5Z6gGoCfKv2LpjzhnTDWGoRtFs80vE3wGdY5WaEcPpR55fVtUt1tZsRkusuaEWgqq8QkXcCzybJBvoRVX1TZSMLnGWlHoBJL2PHiiDHasr3fGjDqqAiBCTrCjejeS0EWU3piEUbBWES53Js6SX7b5adM2XBe1ic31Z1vQcYqupbRGRPRC5R1cerGFjoLNsUA+nNXe4u6DyKIF3aIqsDlC/kdSGEYf2sjhHst1t85LELC1/3hTyWnms3mNnElsci8PV6yj3TReRbgNcDr5o8dRXwe1UMahvoLSl/YNhru/XXm1XS0htGqrSFz/RWpOfCZM9EADGCXo4Yge9ZLYY8i5k9x4pgMBozXNCv2OB7UyebJd93AF9AspEMVX0v8OQqBrUNLCv1YHDtusjjLvF9ZWPI4y45aLcYDMccOi7sVzd5rJ+D3WYQrqHOYEi71WBnibXqWtZZYsfy2Br4u3vbRhH0J/WCABCRFltSF6gKVpma4N51kdddYsbnM7mCxYHImudGtVdRf4u66fVHS+UEUwCubpdqNZUA6sJGEfy5iHw/Sc2h5wCvA95YzbDCJ6mMuVwR7LXdui56eYKKgTR1zxVUnBYo81zWwTCpnbR0ldzicKT0h37LmncB5dQiyBmMB/eVAOrCRhHcCpwD3gV8K8negB+sYlDbQHdFFgK4D3p1cvjNQ1kld7bIDZbHzRiKgs8l6yR91NXmrlwWgefVbG3SR8ci8hrg7SQuofvU1210G8A6XEN5UyrB3wvaYJqXZPUrNvg+eQ15rqV0MbbL99t1DKsS8iRZ7O+2GGtSqfTYzvJj85AvK8tvRWuTNfQVwPuAnwV+HjgjIs+tamChk29l4z59VGSVu8TvoJehu6RfscH3TA9DZ8VmJ0jJ6vnvuqrUA7hfzOSxCKalLTz9fm32Efwk8MWqegZARD4J+EPgj6oYWMiYUg+rYgT7rtNH+yP224s3sZlzmmN9ppsjqLgXiKxJldXVq2TwX+l1+0OeetnSxojTedXtD7niYGUVnFznhOWWNPhdxsMmRvCQUQIT7gcecjyerSBPqYfk9ST7YTx25+vMc05zrM/kcZeEYhHkk9XvrBaDnaz1WQTmdV+vJRuL4F4RuQt4LUmM4IUkPYi/CkBVf6eC8QVJnnop6de7gyGXHNspfd7Oivr8AMd33Je2WAd5XAihuMFWNRuC2SrZ16wWQ3eQI33UZPA4ygbLu3Dbaze9VbQ2iuAY8BHgWZPH54DLgX9NohiiIsjJtBH2ygDfrLa6C0XQy3FzbDSEvR3/d9z2BiOuOFgeFJ01Ovdf1lUB1BCSAPKUeoCZ0nNpEbQaQntFyZUDj/dq2GQNfUOVA9km8mQhJK/PLugrHZx3VXMNg88mrqHbH/K0Jy1fJR/badDYNuvHY1n7w9WlHsC9y88on2WxNUi+44/1/OxbbJM19AoRuVREdkTkrSLysIh8XZWDC5W8waf9ttsLutNf7RqCMKpy5snKEhHvq3KaVfIqWaeBcY8tPePqWe0amljSjiy9Tn+08vs15/X1WrIJFn+pqj4GfCVJL+JPBr6nklEFTp5SD+nXnZm4OVIqzXl9XjlCvqAi+J3pAXDhcMx4Sb9iw26rQashXsuaN2hbRfpoLku67W9hPxtFYJzUXw78lqo+WsF4toI8pR6S192ubFY11zC4Lm1RN+Ox0jtcHVQEU+HVX1lXtW40GOvHZ0WQp34UHE0fdUF3kK/ZzP6uvw3sbRTBG0Xk74BTJJ3KTgD+FzhfA3lKPSSvu81q6fZX55uD/6vk84cjNMcqGfx3g+V1MybH+JvVAqmWnCt+13arQbvZoONw3uR1DbksbVEnuRWBqt4KfB5wSlUPgR7wPPP6pBBdJAe26aMublSjsXL+MJ9F4PvKMa8LwRzjs6x53YzmGF+zWiBfPw3D/m7TsSWdL8lirIm7zjesWlCp6kdVdTT5u6uq/5h6+cedjixgepNSD8dX1EExAT4XF7S5AeRe2XjtLslncUHiRvBZ1p6FrL4Hxm0U/F7bnYLv5th/A36X8XDZi3B5blVkSidHqQeYmfsuJu+0Zn0eReBwEq0DG3fJQQWNzuvExiLw3eWXp6KswaXLL2/igessvzpxqQj8c4ytibxZCI2GJIFbBxfW7IaRb+XosrRF3eQNKoL/rqG8bkbwPzDes3KDuZM1b5LFvuPSFnXib3dyj+nkTOMEd5kIdqtkf01csHMhBBMs3gZZ1+AGGwzHDEbjnHEJf3equ1QE73f4WUHTy7mxC5JMDxcXVp4uS4Z0aQsfsY0R9Idjhp72Lc7TptLgu/XT6Q/ZaQq7rbw5/eVlnZaDyZk+Cn4uoGx2Fu+JyA+JyC9OHl8nIl9pXlfVr6pigCGSNwsB3E3ebs6UVfC/Lo1d1pDfO26tAqi7bhYV6yJpNmRhSTt0qeZyM25JjOBXgD5JCikku4t/1PmItoA85Q8Mrkxcq5ujxxc02LtL0u/xjc5gSLvVYGdFQTSAg3aLwWjMYOiv9ZN33hw4ynyzSrIwiwoPryUbRfBJqvoK4BBAVc8TM4UK0cu5UxFMcxp3rqFcQcXpBe3n6tGMey9Hm0KzOcnX/PpejgY8hn3PZc2bZAHJ7+pic1d3Wik4vyXt47yxUQQDETnOJDto0qGsX8moAidPa0GDO9fQ9qySu4Mhx3YatPKskj1v2JI3tRECcPlZJFkc7LYYjpV+SevHLivL33lj04/gvwD/E7hGRH4D+ALgxVUMKnSSLev5VjauMj06Fqtkn4NeYOl683jygp2s3icBWCZZQCJrmQb2Ntl2rktb1IlNP4I3i8hfAZ9L4hL6LlV9uLKRBYpNqQdwaxHstZs0Gqu9ed6vHC1Wya4rvNaNzSrZe1n7I05ckq8Hcbr/wuX7yxsULSj/1Y0AABxPSURBVGPWOyTvd+znBsWVtrOIfKb5BzwN+DDwIHDt5LlV779RRO4TkTMicmvG6y8TkXeLyDsnfQ6eVkQQX5hWHrWJERyW39zVs7hhuCxtsQ66/ZFVdgn47Ddf3Z3M4LvLL08DHoMrpdeziBGY8/o4b/J8qz+55DUFvmTRiyLSBG4HnkOSZXS3iNypqu9OHfbXJIXseiLyH4BXAF+TY1xeYpOFYI5TxcqKyKJjE1R03Oqvbmxcb/sBxAieetmxXMdOm9P4+rvmrPkD7jqy2aSPQjJ3fJw3K6VT1S8u8fk3AGdU9X4AEbmDpGLpVBGo6p+mjn8bEHTXM5tSD8lxswu6jCKweb/L0hbroDvI7w7wfZVcJFjsY1YLTDKkLNJHofz+kG5/SLMh7Lby5dUkBRv9u5ZsNpQdm7hxfkdE3iAi/1FEVi1FrgIeSD0+O3luEd8E/NGC879ERE6LyOlz587lHfbGYe0acnhB2ygSn6ty2gQVj+80EZnVsfGN7sA+fdTHG5VNqQdwl8FjXG+rCkQakpief/PGJn30V4FPBX4O+HngeuDXVrwn69vLdHZP+h+fAn4i63VVfbWqnlLVUydOnMg96E3DplokuMtq6Q7yZ5eA31U5bXZui8jEnPdv8pp+xduQPmqT/gzuZM3blCZ9Xh/njY2v4Z+o6j9PPf5TEfmbFe85C1yTenw1SaD5CCLybOAHgGepatB7E7qWWQjuLugRe0/Kn0bnc10aW+vH10yP/nDMcKy5Zd1tNWiIn0kAXWtLepIE4GABZW1Je3gt2VgEfy0in2seiMjnAH+54j13A9eJyNNFpA3cDNyZPkBEngm8CrhJVR+yGI+XdEvECMpgk29uzuvjylFVra2f/d2Wl7nfNpudYNa32Mff1T7Jwo1L1WbzJySWtI/fr41F8DnA14vIByePrwXeIyLvAlRVP33+Dao6FJFbgDcBTeCXVfVeEbkNOK2qd5K4gg6A1038cB9U1ZuKi7TZ2FQBTY5zc0H3bFfJ7SYPdwalzrkOLhyOGWu+apGG/XbLyxiB2RiWN30U/HVd2CZZtJsNWg0pnz7az9fn22B6eahq7rjCJmCjCG4scgJVvQu4a+65l6f+fnaRz/UVW1+nC4tgPFaroKI57wce6RU+57qYpfvZyOpnVU7b1EZIlIaPwWLbeWOsHxeW9DX7e7mP30+Vtiizo7lubHYWf0BEnkji82+lnv+rKgYWKjalHsCNIugd2pnV4G8TE9sbBiSyPvixC1UNqTKKyuqj0rPNtgM3slq7GVN7NYJUBCLyIyS1hd7HLPNn6YayyMX0LEo9wExhlLmgbVr8Tc/bbnlZk6Zr0UjEkMjqodKzaMBj8DUJwLbUA+BkL0zPYuc2pHeqj3hSqTPXi41r6KtJSlH75zjeIGyzEFrNBsd2GqXM+SIuhIPJxhjffJ22WVlgAuMeKr0CCn5/t8WjXf9cfrZJFsmx5du82iZZ+Jqia5M19LfAE6oayLbQsWiuYSjrppnW57dc2aj6V6myyA3D1z0T0wCqVWDczxiB7f4bKD9vhqMx/eHYzpL2dKe6zR3px0hSSP+WVB+CkDN8qsCmuYahrDlfxCJwVdqiborKev5wxGisNHO67DYB2/RR8LcoWm9gV+oBksXAuceLb0uyTVmFdH+LcBXBa4AfB94F+NnrbgOw6btqSDaplIgRWKasJsf62ct3Vi3SLn3UvPeSYzuVjKsKpumjVtaPr0kASdabjZuybAG42Sa2YjECn7C5Iz2sqj9b2Ui2hO5gyJMvyVct0lDWdWGbgw3+NmyZBhVt3CWpYmw+KYJOf8hOU9ht5f9d99qtZEfyaJyrg9umYFOC2lA2RlAoBuNp5V6bb/YeEfkxkp3BaddQTB+1oNsfsX+F/QX9aLd4jL6YievnBW0mr80qed9Tc76I2y5t6V123B9FYNNPw1DWDVbWpeoTNt/sMyf/f27quZg+aknXcqciJKuMBx4tnulRxDXka1P37mBIu9Vgx2K1m3YN+UTiLrFPPEjeO+Sy4z5ZP/b9OPbbTQajMYPhmLZFbMFQZOe2UbTBuoZK9iWITCi6iisTIyiSXeJrU3fbapHgbwvHookH5r0+YdNsyJCWtd2yb1dZJFNpt9Vkp1m+tEXdWM0YEfkKklLUUye3qt7melChMi31UMTXWSp9dMixnYZVRoy/N4z8JagNvjZssd2TAn67wS63KPUAR92bTyzQt7hIVhb4WYHUpjHNK0laSH4nSZ+BF5L0MI7kxJR6sF3ZHEyCXqrF+hYX2bvgqyLoWDSlMUz95h7Kav27Tt1g/im9otdwUVmLBIvBz8wsG8fZ56vq1wMfVdX/CnweR3sNRFZgSj0USR8da1JZs9B5C6wcXZS2WAdFg4rgX+cu2/IH4LMbzN7S2ytp/RQp4WGO922vho0iMFW5eiLyCcAQeLr7IYVLkSyE5PiSF3SBvQsuSlusg0JBRZ+tn22StURgvAjd/pCGJO1MbUjavPr1/doogjeKyBNI+gf8FfAPwG9VMqpAKZLGmT6+6AXdKRBoAz9N3CJBRWP9eBcYL+Qu8W+j4OEk88c+a6hcNphRPra1tnycNzbf7N8BI1V9g4hcD3wm8HvVDCtMZk1pCmY/FLyge4MRlxcIliV52H5d0L0CK8dGQ9hvNz2U1d76KbtKXge9gguoWbC4mNIr8v1C+dIW68DGIvghVX1cRL4QeA7wP4BfqGRUgTINPtkGM6e7fItd0EVcCJCYuL6tkgvL6qBSZZ0MhmMGo7H1npTjO01E/FIEpo2otaVXMgmgMxhabUw0+NgO1EYRmDvCVwCvVNXfB+yXmVtMkbzk5PhyF3SRTWzgX1XOpF+xfVARjDnvj9IrmtEiIqVr8NTNdLd4wRhBmdiaresNkoWbT4sKsFMEHxKRV5H0JbhLRHYt37/1FKmVnz6++AVd1MT164LuD8eMxlrYnPdJ6RVdVCTv8SurpWg+/24r2TtTNEbQLeBmBD+b/9jcyL+apAn9jar6MeBy4HsqGVWgzCpj2pq4xYNeySq5xMrGowu6qOsN/NsEZHLji96oOh4p+KJJFiIy6VJWdB9BUeuyyeFIGQz9KdJsU2KiB/xO6vGHgQ9XMahQKVLqAWaVNIu4Ls4fjlC1N6vBv6buRW8YkKw2H3rcn77FRSrKGnxT8J2pa6jezLciO7dhNteKlrZYB9G1UyPd/pDjO03r5idlYgSzvQvFgl4+3jCKy+qT0ivmLgH/XEO9QRlZi/ejLtqUycfKvVER1EiROkOQbO7abRXb3FU09Q7Kl7aomyJVVg2+BcbLyepnsLhonKtoEkC3QGkWc07wq4xHVAQ1UqRapKHo6rxMULFsaYu66RTMLjHv8UkRmJtb4WCmRzGCTsEkC5j0aC7wu47GyvlD+xIe4Gdhv6gIaqRoFgIU99eXCaD61n+1aFYWmJvjiPHYD+unWyJG4JvSM6Ueju3Y366KLqC6Jd1R4NdejagIaqRItUhD0dzvoruZk/f4dUGXuTkapWcqxG46ZSy9A9+SAAbFSj3AzL1pfc4y7igP27xGRVAjvcGo0E5FKB70KrNKnmY/eOJGmDUbLy6rL2UmeoMhzYawW6Dz1v5ui/OHI0YeWT9FbshA4fRR856imUrgVz2nqAhqpGj5Ayge9Jr18C2e/eDL6rHMKs63TI+kTWWz0Cp53zcF3y++gCoaGC+TlVW2tMU6iIqgRrr94XRPgC1Fs1qmKZUF9xGAPxd0pz+i3WwU6k+775nSK+VmNFktnshadEMkJLIOhmMOR3YJD9u0qICoCGqlaDVDSFZxRdwW0wbcBevvgD8rx6QpTVHXmynP7JOsxRMPwJ8bVbkki2JKz7h1iiigsqUt1kHlikBEbhSR+0TkjIjcmvH6F4nIX4nIUEReUPV41oUp9VAmfbSoidtuNdhp2v/Ue54FizsFGvAYfAvwdfqjQu4+8FPWwvNm4uO3LakxK3Rnf96ypS3WQaWKQESawO3Ac4Hrga+d9DJI80HgxcBvVjmWdXP+cMRYi5maMEkfHYysN3eVcSGUKW2xDopWiwT/WjgWacBj8DEbrExszXyGDUW7CRp827RXtUVwA3BGVe9X1QFwB/C89AGq+n5VfSfgx66lgpRJ9zPvG42VvmUhqzKb2HwLehUtEgZ+BsaLukt8y2op4wYr2oinTIzAvM+XeQPVK4KrgAdSj89OnrNGRF4iIqdF5PS5c+ecDK5OpqUeCpiayfsKXtCDUeEbxs4k8OqL37xokTCYKT1f/LplZPUvCaC4pWdcO7YK3ihJ237FBrNB0ReqVgRZuW2FkpdV9dWqekpVT504caLksOrHhUUABS7oEmY1TDbkeHLDKBVUbPvmGirhN/fIDTYcjblwOC7kq4fisppmTg3LApHT8xYsbbEuqlYEZ4FrUo+vBh6s+JwbSZm85OR9xTI9yioCn0pRF23AA9BsCMd3/Jm8ZfekgB8WQZnsnfT7iriGys0bfxZQUL0iuBu4TkSeLiJt4GbgzorPuZFMG4mUzf22zX4YjAoHFcGv2vVJvnkJWT0x5w9HYwbDceE9KXs7JlV282UtU2U1/b5i86akJe2JmxEqVgSqOgRuIels9h7gtap6r4jcJiI3AYjIZ4vIWeCFwKtE5N4qx7QuZk1pCgZuC7ouuiVSKsGvSpXd/rBwSiX4067SxJuKytpomPTGzZe1TBonpPdM2LtUi+5mhuKlLdZF8VmTE1W9C7hr7rmXp/6+m8RlFDRlsxCKZrWUCbRBMt6Pnz8s/P666A9HHI60nKyeWD+daWXMktaPD7KWqJUFSbC3IcXSR4vGmyCmj0YWUD5YbJ/poaql0kfBn4Yt3ZJZWeDP5C27qADjutj8FWuvpKwikij4AhvKyi6gipS2WBdREdTErNl4yfRRiwu6PxyX2sQGiUvKh4qcZYrrGfZ2m150lSrTY8Lgixus7MYuSH5XW1mTSsHlFAH4U88pKoKa6PaH7LYatAqUeoBimR4dBzcMb1bJJRqJGIqW8aibqfVTUsF7IeugXIwAivWj7pTYuQ3FS1usi6gIaqKsr77datBuNqyCXi5cCEVLW9SNE3eJLzGCEg14DL7sDykbIzDvLbaPoLxF4MN3DFER1EavYOP6NPu7Tas0uFlTmnKrqSKlLerGlaw+mPI9R9aPD26wsjECmFTutZg347GWnq9F9y+si6gIaiKpjFn8JgX25vzMrC43iWDzL+hZmmFZ62fojfVT7ndt+uEamoyxaKkHSH5XG0vatCstY3EVLW2xLqIiqImyWQhgb86XzVRKv3fTL2gXQcX93RZjTSrFbjIu3CU+pY+WKfUA9rK6can6U8YDoiKojbJb1sG+3EPZshbJe/1oYrJNk7fbH9IQOLZTfPoa19B4w/sWl6k8arDt9+1iURFdQ5FMuoPiRcIMtrt8pxVPS5m4xbbo143JiS+VXTJ576bHCbqDJJBZpF+xwSj43sZbP+UtaVs32HTntoNg8abPG0NUBDVRNgsB7He+ukgf9WmV3GoIuwX6FRt8krXsKnnPo9hPmVIPkPyuFw7HDHNu7nKRlVW0tMW6iIqgJspUizTY5kO72oGafNZmX9Dm5lhulezLzbG8delLg/Vuv3g/DYNtIx4XLtWipS3WRVQENaCapKOVDxY3rVxD3cGI9qS5TFF8aepetlokpM35DVd6AwfuEk92vjqVNec13C1Z8RSKl7ZYF1ER1EB/OGY01tIm7l6B7IfSZrVPLoSS6bnT3aBeyFrebw6eyFraDWZXp2tWt6rkeT0p4wFREdSCiywE8/7DkdIf5jdxS8clPHGXuHK9gQ+yutic6I+sZTYJQtoNZucacpLcseEWlyEqghpwUSQseb/dJhUXGRdFSlusAxf7NHwKFpe9OU4VwYa7LnoDd4uZvMUTze9f1urypU4XREVQC10HaZwwq6yZdxXXc5CympzXrrTFOugNRs5cQ5seI+gNyrtLfEgCcFHqAez7UfcGQ47vNGmW2MRmzrvp88YQFUENuAg+QTr7If/Kpuw5IbmgN31l48L6aTUb7LYaHrhLystqYkebLOts3pS1fuwSHly43sx5N92SNkRFUAMuSj2k358/6FXerAY/KlW6yK2HzTfnh6MxFw7HTvakwGa7wVyU24ZiMYKyrjfwp4wHREVQCz0HtWGS99vFCFyY1WCqnm72yqbrTNbNrsrpoiAaQLMhHN/ZbJefix4TYB8jcFHWwpx3k7/fNFER1EDZBtwG292gnZJtKg2b3rDlcDRmMByXalNp2NvwqpwuNgkaNt114aLKKswql1rNGweWtC8VXiEqglpwmT6a/rxlzPoVu4kRbLKJ6/LmuOluMLeKYLNldVHqAaDRkMlNOa9ryE2ShW1pi3USFUENuFrZ2MQI+sMxw7GWVj7mvJucXeJK0YIPN8fyDXgMm57V4sqlCnZuGpfxJshf2mKdREVQAy5KPUBqh2SOC8v4uV24S/YtS1vUjZG17C5qmLXm3FR6jhYVsPmBcVfZduYz8sradbB3wZwT/KhAGhVBDXQd+ep3Ww1aDcm1Yp1aIVuxSnZ4w9hwN5hL62fPsr9F3bionmvYtyj30HWUPmpb2mKdREVQA65MTRHJfVN2ecOwLW1RNy6qRRo2PTDuepW8yTcpV6UewCj41devqk4K3ZU/p23a6jqJiqAGXFRQNBzstnK6hlyukje7YYurImGQfL+9wWhj+xa72qUOcLDh1TGNrK7cYHlkPX84QtWdooX8aavrJCqCGuj2y5c/MOy185m4nb67GMHebv5spXXgcuW4t9tkNFb6w83M9HBVtwo23zXU7bsp9QD5K/d2XLpUPdi0Z4iKoAZclXqA/K4L1ymVsLkFyly6Sza9YUu3P0Sk/J4UmK2SN9b6cbSxC5Isqzwumq7LrCxPenlAVAS14KIypiFvnrvrlErY3KCXU1k3vP9CZ9Kxq0wnNsP+bgvVxB2yiXQdlKA25E2VdWlxxRhB5AiuSj1A/nIPPacWgV1pi7rp9Uc0S/YrNsyU3obKOnCTgQabX3bbVZIFzEqHjMfLrR/XiQcQYwSRCcmWdXcrm1yuIZNb76TswqavkpPuZG5WyZttzrsqfwD2/S3qxqmsOX9X87qLGIFtaYt1UrkiEJEbReQ+ETkjIrdmvL4rIr89ef3tInKy6jHVictSD5A/5a/TH9JytEr2wW/uyvW2batk85mbSLcC62eV0nO5c9u2tMU6qVQRiEgTuB14LnA98LUicv3cYd8EfFRVnwH8NPDjVY6pbgajpNSDU0WQ0zW0v+vOlwyb27DFpevtYGrOb6asXUfNhiDdnGYzFUHP0cYuyJ/w4NKlaj7Hh53FbqRdzA3AGVW9H0BE7gCeB7w7dczzgB+e/P164OdFRLSCVIZX/vn7eMM9Z11/7FJGEzHcuYaaDIZjnv1Tf86yW/xHHrvgbJVs3Es/9yfv5dff9gEnn+mSBz92nmc8+cDJZxlZb/uDe/mZt/y9k890yQcf7fGFz7jCyWcZWV/22r9xlt7skg8+2uOzT17u5LOMe/MbfuXupVbyx84fHjm+LPu7Lf7wnR/mng981MnnAfzwTZ/KFzi6BgxVK4KrgAdSj88Cn7PoGFUdisjHgScBD6cPEpGXAC8BuPbaawsN5oqDXa670s0Nw4Z/dtVlfMk/vdLJZ934aU/h7x/qMBovz3O/7soDPu8Tn+TknMd2mnznlzyD953rOPk811x35QFfev1TnHzWJ1x2nBd//kkeevyCk89zzXVXHvCCz7rayWd9ylMv5ebPvobHLhw6+TzXfPKVl/B8R7J+9skn8vzPvJrzh6tX59dcvselx9zcGr/1iz6Rv3jvOSefZXBlraSRKnOIReSFwJep6jdPHv974AZV/c7UMfdOjjk7efy+yTGPLPrcU6dO6enTpysbdyQSiYSIiNyjqqfmn686WHwWuCb1+GrgwUXHiEgLuAx4tOJxRSKRSGRC1YrgbuA6EXm6iLSBm4E75465E3jR5O8XAH9SRXwgEolEItlUGiOY+PxvAd4ENIFfVtV7ReQ24LSq3gn8EvBrInKGxBK4ucoxRSKRSOQoVQeLUdW7gLvmnnt56u8LwAurHkckEolEsok7iyORSGTLiYogEolEtpyoCCKRSGTLiYogEolEtpxKN5RVhYicAzav1kE+rmBu1/SWsK1yQ5Q9yr45PE1VT8w/6aUi8BkROZ21sy90tlVuiLJH2Tef6BqKRCKRLScqgkgkEtlyoiKon1evewBrYlvlhij7tuKN7DFGEIlEIltOtAgikUhky4mKIBKJRLacqAgqRkQuF5E3i8h7J/8/ccmxl4rIh0Tk5+scYxXkkVtEPkNE/q+I3Csi7xSRr1nHWF0hIjeKyH0ickZEbs14fVdEfnvy+ttF5GT9o6yGHLK/TETePfmd3yoiT1vHOKtgleyp414gIioiG5dSGhVB9dwKvFVVrwPeOnm8iB8B/ryWUVVPHrl7wNer6qcCNwI/IyJPqHGMzhCRJnA78FzgeuBrReT6ucO+Cfioqj4D+Gngx+sdZTXklP2vgVOq+ukkvclfUe8oqyGn7IjIJcBLgbfXO8J8REVQPc8DXjP5+zXAv8k6SEQ+C7gS+OOaxlU1K+VW1b9X1fdO/n4QeAi4aNejJ9wAnFHV+1V1ANxB8h2kSX8nrwf+lYhIjWOsipWyq+qfqmpv8vBtJN0KQyDP7w7JIu8VwEY2w46KoHquVNUPA0z+f/L8ASLSAH4S+J6ax1YlK+VOIyI3AG3gfTWMrQquAh5IPT47eS7zGFUdAh8HnlTL6Kolj+xpvgn4o0pHVB8rZReRZwLXqOof1DkwGypvTLMNiMhbgKdkvPQDOT/i24G7VPUBnxaIDuQ2n/NU4NeAF6nq2MXY1kDWDzefm53nGB/JLZeIfB1wCnhWpSOqj6WyTxZ5Pw28uK4BFSEqAgeo6rMXvSYiHxGRp6rqhyc3vIcyDvs84F+IyLcDB0BbRDqquiyesHYcyI2IXAr8IfCDqvq2ioZaB2eBa1KPrwYeXHDMWRFpAZeRtGf1nTyyIyLPJlkkPEtV+zWNrWpWyX4J8GnAn00WeU8B7hSRm1T1dG2jXEF0DVXPncCLJn+/CPj9+QNU9d+p6rWqehL4buBXN10J5GCl3CLSBn6XRN7X1Ti2KrgbuE5Enj6R62aS7yBN+jt5AfAnGsaOzpWyT9wjrwJuUtXMRYGnLJVdVT+uqleo6snJ/H4byXewMUoAoiKog/8GPEdE3gs8Z/IYETklIv99rSOrljxyfzXwRcCLReQdk3+fsZ7hlmPi878FeBPwHuC1qnqviNwmIjdNDvsl4EkicgZ4GcszyLwhp+w/QWLtvm7yO88rSS/JKfvGE0tMRCKRyJYTLYJIJBLZcqIiiEQikS0nKoJIJBLZcqIiiEQikS0nKoJIJBLZcqIiiEQcIiIvFZH3iMhvrHsskUheYvpoJOIQEfk74Lmq+g/rHkskkpdoEUQijhCRVwKfSFJC4D+tezyRSF6iRRCJOERE3k9Sd//hdY8lEslLtAgikUhky4mKIBKJRLacqAgikUhky4mKIBKJRLacGCyORCKRLSdaBJFIJLLlREUQiUQiW05UBJFIJLLlREUQiUQiW05UBJFIJLLlREUQiUQiW05UBJFIJLLl/H9cyjfJ6h+LjwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "psd = np.abs(xs).mean(dim='t').plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The bispectrum and bicoherence are created as Dask arrays and are evaluated (possibly in a parallelized way) on demand."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.Dataset>\n",
       "Dimensions:      (f: 15, f1: 15)\n",
       "Coordinates:\n",
       "  * f1           (f1) float64 0.03125 0.0625 0.09375 ... 0.4062 0.4375 0.4688\n",
       "  * f            (f) float64 0.03125 0.0625 0.09375 ... 0.4062 0.4375 0.4688\n",
       "Data variables:\n",
       "    bispectrum   (f1, f) complex128 dask.array<shape=(15, 15), chunksize=(1, 15)>\n",
       "    n1_norm      (f1) float64 6.56e-25 2.097e-25 0.25 ... 3.232e-25 2.228e-25\n",
       "    n_norm       (f) float64 6.56e-25 2.097e-25 0.25 ... 3.232e-25 2.228e-25\n",
       "    v_norm       (f1, f) float64 dask.array<shape=(15, 15), chunksize=(1, 15)>\n",
       "    norm         (f1, f) float64 dask.array<shape=(15, 15), chunksize=(1, 15)>\n",
       "    bicoherence  (f1, f) float64 dask.array<shape=(15, 15), chunksize=(1, 15)>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bisp_ds = calculate_bispectral(xs, freq_dim='f', avg_dim='t')\n",
    "bisp_ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[########################################] | 100% Completed |  0.2s\n"
     ]
    }
   ],
   "source": [
    "with ProgressBar():\n",
    "    bicoh = bisp_ds.bicoherence.compute()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.QuadMesh at 0x7f1498d5cc18>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEKCAYAAAAB0GKPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAep0lEQVR4nO3de5Ad5X3m8e+DjPAWl+IyiqMghCClMlFsr8CDoNZVWGZ9EWQLZENs5HW4hI0WAmtiBRdQ8ZpdMOvYEBNTZgFtkEE2AWNtMKpYWPYSyZQTQySBEEgOIAtsBpEFmZuzBJBmfvtHv0c0wzndPTPdMzpnnk9VF6fvb0PVPLyXflsRgZmZWZm9JroAZmbWHRwYZmZWiQPDzMwqcWCYmVklDgwzM6vEgWFmZpU4MMzMxomkZZKek/Roh/1HSfqppNclXTxs3wJJj0naKunS3PYjJD0g6QlJ35E0tanyOzDMzMbPLcCCgv0vAJ8FrslvlDQFuB44CZgDLJI0J+3+CnBtRMwGXgTOrbnMuzkwzMzGSUTcRxYKnfY/FxHrgJ3Dds0DtkbEtoh4A7gDOFWSgBOBFem4W4GF9Zc8846mLtyUvr6+mDVr1kQXw8y6wIYNG3ZExLSxXONjH9o3fvXCYLX7bXp9M/BabtPSiFg6lvsnhwJP59YHgOOAQ4CXImJXbvuhNdyvra4LjFmzZrF+/fqJLoaZdQFJvxjrNX71wiD/uHpmpWOnTH/itYjoH+s921CbbVGwvRFdFxhmZuMpgCGGJroYA8BhufUZwHZgB3CgpHekWkZreyPch2FmViAIdsZgpaVB64DZaUTUVOAMYGVks8euAU5Px50F3N1UIVzDMDMrUVcNQ9LtwHygT9IAcDmwN0BE3CjpN4H1wAHAkKQ/AeZExCuSLgRWA1OAZRGxOV32EuAOSV8CHgJurqWwbTgwzMwKBMFgTZ+BiIhFJfv/maxZqd2+VcCqNtu3kY2iapwDw8ysxFBz/chdxYFhZlYggEEHBuDAMDMr5RpGxoFhZlYggJ3+lDXgwDAzKxSEm6QSB4aZWZGAQecF4MAwMyuUvelt4MAwMyshBttO2TT5ODDMzApknd4ODHBgmJkVyt7DcGCAA8PMrNSQaxiAA8PMrJBrGG9yYJiZFQjEoL8EATgwzMxKuUkq48AwMysQiDdiykQXY4/gwDAzK5C9uOcmKXBgmJmVcqd3xoFhZlYgQgyGaxiA61lmZmWGUKWljKRlkp6T9GiH/ZJ0naStkjZJOiZt/5CkjbnlNUkL075bJD2Z2ze31ofPcQ3DzKxA1uld25/KW4BvAMs77D8JmJ2W44AbgOMiYg0wF0DSwcBW4Ie58z4fESvqKmQnrmGYmRVodXpXWUqvFXEf8ELBIacCyyNzP3CgpOnDjjkduCciXh3tM42WA8PMrMRgqNJSg0OBp3PrA2lb3hnA7cO2XZWasK6VtE8dBWnHgWFmVqD1pneVBeiTtD63LB7h7dqlzu7PN6XaxnuB1bn9lwFHAccCBwOXjPCelbkPw8ysxFD1UVI7IqJ/DLcaAA7Lrc8AtufWPwncFRE7Wxsi4tn083VJ3wQuHsP9C7mGYWZWIJt8sHINY6xWAmem0VLHAy/nAgFgEcOao1p9HJIELATajsCqg2sYZmYFArGzpqlBJN0OzCdruhoALgf2BoiIG4FVwMlko6BeBc7JnTuLrPbx42GXvU3SNLLmrI3AebUUtg0HhplZgQhqe3EvIhaV7A/ggg77nuLtHeBExIm1FK4CB4aZWaFqL+VNBg4MM7MCQX01jG7nwDAzK+EPKGUa/bcgaYGkx9K8KJcWHHe6pJA0luFoZma1C8RQVFt6XWM1DElTgOuBj5CNLV4naWVEbBl23P7AZ4EHmiqLmdloBbCzvrmkulqTNYx5wNaI2BYRbwB3kM2TMtyVwFeB1xosi5nZKInBikuvazIwSudEkXQ0cFhE/G3RhSQtbr1q//zzz9dfUjOzDoLsTe8qS69r8gnL5kTZC7gW+NOyC0XE0ojoj4j+adOm1VhEM7NyrmFkmmyYK5sTZX/gPcDa7I12fhNYKemUiFjfYLnMzCqL0KSoPVTRZGCsA2ZLOgJ4hmxK3k+3dkbEy0Bfa13SWuBih4WZ7UmyTu96pgbpdo0FRkTsknQh2TS8U4BlEbFZ0hXA+ohY2dS9x8uRX/9abdfadtGS2q5lb+f/VjZ6/qZ3S6NjxSJiFdlkWvltX+xw7Pwmy2JmNhpZp3fv909U4cHFZmYl/KZ3xoFhZlag9aa3OTDMzEoNuYYBODDMzApFwM4hBwY4MMzMCmVNUg4McGCYmZWaDG9xV+HAMDMr4GG1b3I9y8yskGqbfFDSMknPSXq0w35Jui59Q2iTpGNy+wYlbUzLytz2IyQ9IOkJSd+RNLWWx27DgWFmVmIofde7bKngFmBBwf6TgNlpWQzckNv3rxExNy2n5LZ/Bbg2ImYDLwLnjuTZRsKBYWZWIBslNaXSUn6tuA94oeCQU4HlkbkfOFDS9E4HK5u59URgRdp0K7Cw8sONkAPDzKzACD/R2tf6dk9aFo/wdkXfEXpnuub9klqhcAjwUkTsanN87dzpbWZWomJzE8COiOgfw62KviM0MyK2SzoS+DtJjwCvFBxfO9cwzMwKtEZJVaxhjFXH7whFROuf24C1wNHADrJmq3cMP74JDgwzsxLj+InWlcCZabTU8cDLEfGspIMk7QMgqQ/4ALAlIgJYA5yezj8LuLuOgrTjJqkx8HcRuof/W9loRYhdNb3pLel2YD5ZX8cAcDmwd3afuJHscxAnA1uBV4Fz0qm/A9wkaYjsf/T/PCK2pH2XAHdI+hLwEHBzLYVtw4FhZlairhf3ImJRyf4ALmiz/R+A93Y4Zxswr5YClnBgmJkV8Jveb3JgmJmVcGBkHBhmZgX8AaU3OTDMzEqM4D2MnubAMDMrEAG7/AElwIFhZlbKTVIZB4aZWQH3YbzJgWFmViIcGIADw8yslDu9Mw4MM7MCEe7DaHFgmJkVEoMeJQU4MMzMSrkPI+PAMDMr4Lmk3uTAMDMrElk/hjkwzMxKeZRUxoFhZlYg3Om9mwPDzKyEm6QyDgwzsxIeJZVxPcvMrEBEFhhVljKSlkl6TtKjHfZL0nWStkraJOmYtH2upJ9K2py2fyp3zi2SnpS0MS1za3v4YRwYZmYlhkKVlgpuARYU7D8JmJ2WxcANafurwJkR8bvp/L+UdGDuvM9HxNy0bBzp81XlJikzsxJ19WFExH2SZhUcciqwPCICuF/SgZKmR8TjuWtsl/QcMA14qZ6SVeMahplZgUAMDe1VaQH6JK3PLYtHeLtDgadz6wNp226S5gFTgZ/nNl+VmqqulbRP0Q0kvUvSzZLuSetzJJ1bpXCNBoakBZIeS+1xl7bZf56kR1K7208kzWmyPGZmoxEVF2BHRPTnlqUjvFW7dq3d9RtJ04FvAedExFDafBlwFHAscDBwSck9bgFWA7+V1h8H/qRK4RoLDElTgOvJ2uTmAIvaBMJfR8R7I2Iu8FXga02Vx8xsVGrs9K5gADgstz4D2A4g6QDg+8AXIuL+3cWLeDYyrwPfBOaV3KMvIu4EhtL5u4DBKoVrsoYxD9gaEdsi4g3gDrL2ud0i4pXc6r7kktTMbI8xgirGGK0EzkyjpY4HXo6IZyVNBe4i69/4bv6EVOtAkoCFQNsRWDn/T9IhrRK37lOlcE12erdriztu+EGSLgCWkLXJndjuQqkdcDHAzJkzay+omVmRut7DkHQ7MJ+sr2MAuBzYO7tH3AisAk4GtpKNjDonnfpJ4ATgEElnp21npxFRt0maRtactRE4r6QYS8iC6bcl/T1Z5/npVcrfZGAUtsXt3hBxPXC9pE8DXwDOanPMUmApQH9/v2shZsmRX6+vFXfbRUtqu1YvCWBoqJ7AiIhFJfsDuKDN9m8D3+5wTtv/0S64x4OSPgi8m+zv9GMRsbPKuU02SXVsi+vgDrLqlJnZniOAULWlC6RWnf0iYnNEPArsJ+mPq5zbZGCsA2ZLOiK1v51BVg3aTdLs3OrvAU80WB4zs1GJqLZ0iT+KiN3vb0TEi8AfVTmxsSapiNgl6UKy4VtTgGURsVnSFcD6iFgJXCjpw8BO4EXaNEeZmU247gmDKvaSpNT81RrROrXKiY2+6R0Rq8g6cfLbvpj7fVGT9zczG7vahszuKVYDd0q6kSwKzwN+UOVETw1iZlamt2oYlwD/GTifrNP7h8BfVTnRgWFmViQgaholtSdIb4jfwJsTG1bmwDAzK9U7gSHpA8B/Aw4nywCRjeg9suxcB4aZWZneapK6GfgcsIGKU4K0ODDMzMr0VmC8HBH3jOZEB4aZWZHWi3u9Y42kq4G/AV5vbYyIB8tOdGCYmZXoopfyqmjN6def2xZ0mMsvz4FhZlamt0ZJfWi05/qLe2ZmJRTVlm6wx35xz8ys61X9FkaXBAZj+OKem6TMupinJB8P3TMTbUV9EXGnpMtg97x/lYbXOjDMzMp0T+2hij3yi3tmZr1haKILUKs98ot7Zmbdr4few5C0F/BOYI/74p6ZWU+oa5SUpGWSnpP0aIf9knSdpK2SNkk6JrfvLElPpOWs3Pb3S3oknXOdpI7pliYe/IuI2NX64l7VsIBRBoak/UZznplZV6pvlNQtwIKC/ScBs9OymDSjrKSDgcvJXrqbB1wu6aB0zg3p2NZ5RdcH+KGk04qCpZPR1jC2jPI8M7NJKyLuA14oOORUYHlk7gcOlDQd+Bjwo4h4IX1S9UfAgrTvgIj4afqC3nJgYUkxlgDfBd6Q9IqkX0t6pUr5O/ZhSOo0Xk+AaxhmNmmM4KW8Pknrc+tLI2LpCG51KPB0bn0gbSvaPtBme0cRsf8IyvMWRZ3e/wO4GtjVZp/7PsxscghGMjXIjojoLz+so3Y3ilFs73yDrCnqPwJHRMSVkg4DpkfEP5YVrigwHgS+FxEb2tzwP5Vd2MysZ4zfexgDwGG59RnA9rR9/rDta9P2GW2OL/I/yQYKnwhcCfwLcD1wbFnhimoKzwC/kHRRm31jSVAzs64yjnNJrQTOTKOljif7dsWzZFN5fFTSQamz+6PA6rTv15KOTzWHM4G7S+5xXERcALwGkPpEplYpXFENYw6wL/CHkpbz1qpP5WFYZmZdr6YahqTbyWoKfZIGyEY+7Q0QETcCq4CTga3Aq8A5ad8Lkq4E1qVLXRERrc7z88lGX/0b4J60FNkpaQpvvuk9jYqvJhYFxk3AD4AjyT7llw+MSNvNzHpfTYEREYtK9gdwQYd9y4BlbbavB94zgmJcB9wF/Iakq8je8v5ClRM7BkZEXAdcJ+mGiDh/BIUxM+sZ3TR1eRURcZukDcC/J6sILIyIn1U5t3RqEIeFmU16PfQBpeQJ4BVSBkiaGRG/LDvJc0mZmZXopRqGpP9C1nfyf4FBslpGAO8rO9eBYWZWpocCA7gIeHdE/GqkJzowzMyK9FgfBtkb45W+fzGcA8PMrEwPBEZuuqdtwFpJ3wdeb+2PiK+VXcOBYWZWQr3xAaXWHFK/TMtUKr6w1+LAMDObBCLiv+fXJe2fbY5/qXoNTyJoZlamvu9hTDhJ75H0EPAosFnSBkm/W+Vc1zDMzIr0Xqf3UmBJRKwBkDQf+F/Avys70YFhZlamtwJj31ZYAETEWkn7VjnRgWFmVqa3AmObpP8KfCutfwZ4ssqJ7sMwMysgslFSVZYu8YfANOBv0tJHmhW3jGsYZmZFeqwPI33/4rOjObfRGoakBZIek7RV0qVt9i+RtEXSJkn3Sjq8yfKYmY1Kb42S+pGkA3PrB0laXeXcxgIjfaDjeuAkso8xLZI0Z9hhDwH9EfE+YAXw1abKY2Y2aj0UGEBfRLzUWkk1jt+ocmKTNYx5wNaI2BYRbwB3AKfmD4iINRHxalq9n7d+m9bMbI8wjp9oHQ9Dkma2VlLLTqXSN9mHcSjZJFctA8BxBcefS4dPC0paDCwGmDlzZrtDzMya0z1hUMWfAT+R9OO0fgLp72uZJmsY7b440vZfu6TPAP3A1e32R8TSiOiPiP5p06bVWEQzsxJR3yipCv26h6f+3E2S1kqakbZ/SNLG3PKapIVp3y2Snsztm1v4OBE/AI4BvgPcCbw/Iir1YTRZwxgADsutzwC2Dz9I0ofJEu+DEfH68P1mZhOuhhpGrl/3I2R/H9dJWhkRW3KHXQMsj4hbJZ0IfBn4g/Si3dx0nYOBrcAPc+d9PiJWlNz/qIj4J0nHpE2tv8cz0xf3Hix7hiYDYx0wW9IRwDPAGcCn8wdIOhq4CVgQEc81WBYzs1GrqX9id78ugKRWv24+MOYAn0u/1wDfa3Od04F7cv2/VS0ha3r6C94aga0v7p1YdoHGmqQiYhdwIbAa+BlwZ0RslnSFpFPSYVcD+wHfTVWplU2Vx8xs1KqPkuqTtD635PsG2vXrHjrsTg8Dp6XfHwf2l3TIsGPOAG4ftu2q1Ix1raR92j5CRKssJwPfJ/uI0kvAyrStVKMv7kXEKmDVsG1fzP3+cJP3NzMbs5ENmd0REf0d9lXp170Y+Iaks4H7yFpndu2+gDQdeC/Z/4i3XAb8M9m3LZYClwBXFJTxVuAV4Lq0vghYDnyy4BzAb3qbmRUStTVJlfbrRsR24BMAkvYDTouI/OdUPwncFRE7c+c8m36+LumbZKFT5N0R8W9z62skPVzlATyXlJlZiZrew9jdrytpKlnT0lua4SX1SWr9Xb4MWDbsGosY1hyVah1IErCQ7DsXRR6SdHzu/OOAvy8tPa5hmJmVq6GGERG7JLX6dacAy1r9usD6iFgJzAe+LCnImqQuaJ0vaRZZDeXHwy59m6RpZJWhjcB57e4v6ZH0JHsDZ0r6ZVo/nLd2vHfkwDAzK1PTi3sV+nVXkE2T1O7cp3h7JzkRUTq6KfkPlQvagQPDzKxId0370VFE/GKs13BgmJmV6YHAqIMDw8ysRBd9HKlRDgwzsxK90CRVBweGmVmR7vrWRaMcGGZmZRwYgAPDzKxQjW96dz0HhplZCQ05McCBYWZWzH0YuzkwzMxKuEkq48AwMyvjwAAcGGZmpVzDyDgwzMzKODAAB4aZWbHw1CAtDgwzswJ+D+NNDgwzszLhxAAHhplZKdcwMg4MM7MifnFvt73KDzEzm9w0VG0pvY60QNJjkrZKurTN/sMl3Stpk6S1kmbk9g1K2piWlbntR0h6QNITkr4jaWpdzz2cA8PMrEQdgSFpCnA9cBIwB1gkac6ww64BlkfE+4ArgC/n9v1rRMxNyym57V8Bro2I2cCLwLljetgCDgwzsyJB1uldZSk2D9gaEdsi4g3gDuDUYcfMAe5Nv9e02f8WkgScCKxIm24FFlZ/uJFxYJiZlVBUW4A+Setzy+LcZQ4Fns6tD6RteQ8Dp6XfHwf2l3RIWn9nuub9klqhcAjwUkTsKrhmbdzpbWZWpnqn946I6O+wTxWufDHwDUlnA/cBzwCtMJgZEdslHQn8naRHgFfGVNoRcmCYmRWo8cW9AeCw3PoMYHv+gIjYDnwCQNJ+wGkR8XJuHxGxTdJa4GjgfwMHSnpHqmW87Zp1cpOUmVmRCDRUbSmxDpidRjVNBc4AVuYPkNQnqfV3+TJgWdp+kKR9WscAHwC2RESQ9XWcns45C7i7hqduy4FhZlYmKi5Fl8hqABcCq4GfAXdGxGZJV0hqjXqaDzwm6XHgXcBVafvvAOslPUwWEH8eEVvSvkuAJZK2kvVp3DzWx+3ETVJmZiXqetM7IlYBq4Zt+2Lu9wreHPGUP+YfgPd2uOY2shFYjXNgmJkVCcDf9AYcGGZm5ZwXgAPDzKyUJx/MODDMzEpUGAE1KTgwzMyKeLba3RwYZl3syK9/rbZrbbtoSW3X6iXZi3tODHBgmJmV8ze9gYZf3Ksw9/sJkh6UtEvS6e2uYWY20RRRael1jQVGxbnffwmcDfx1U+UwMxuTqm95935eNNoktXvudwBJrbnfW6+zExFPpX2u8JnZHqrSPFGTQpNNUlXmfq9E0uLW/PLPP/98LYUzM6usng8odb0mA6PK3O+VRMTSiOiPiP5p06aNsVhmZiMQ9X3Tu9s12SRVOve7mVlXmAS1hyqarGGUzv1uZtYV3OkNNBgYVeZ+l3SspAHg94GbJG1uqjxmZqOloaFKS69r9MW9CnO/ryNrqjIz2zMFfnEv8ZveZmYFxOR4Ka8KB4aZWRkHBuBvepuZlavpPYwK0yUdLuleSZskrZU0I22fK+mnkjanfZ/KnXOLpCclbUzL3FqfPcc1DDOzIjX1YeSmS/oI2WsH6yStjIgtucOuAZZHxK2STgS+DPwB8CpwZkQ8Iem3gA2SVkfES+m8z6fvgTfKgWHWxTwl+fioaQRU6XRJZPPufS79XgN8DyAiHm8dEBHbJT0HTANeYhy5ScrMrFDF5qjyJqkq0yU9DJyWfn8c2F/SIfkDJM0DpgI/z22+KjVVXStpn5E+YVUODDOzIsFIAqOvNe9dWhbnrlRluqSLgQ9Kegj4IPAMsGv3BaTpwLeAcyKiVe25DDgKOBY4GLikhqduy01SZmZlqrdI7YiI/g77SqdLiojtwCcAJO0HnBYRL6f1A4DvA1+IiPtz5zybfr4u6ZtkodMI1zDMzErU9AGl0umSJPVJav1dvgxYlrZPBe4i6xD/7rBzpqd/ClgIPDrGx+3IgWFmVqaGPowq0yUB84HHJD0OvAu4Km3/JHACcHab4bO3SXoEeAToA75U45O/hZukzMyKRMBgPXODVJguaQXwtuGxEfFt4NsdrnliLYWrwIFhZlbGb3oDDgwzs3IODMCBYWZWLAB/0xtwYJiZlQgIz28ODgwzs2JBbZ3e3c6BYWZWxn0YgAPDzKycAwNwYJiZlaj2rYvJwIFhZlYkgHqmN+96DgwzszKuYQAODDOzEvVNDdLtHBhmZkUCwu9hAA4MM7NyftMbcGCYmZVzHwbgwDAzKxbhUVKJA8PMrIxrGIADw8ysRBCDgxNdiD2CA8PMrIinN9/NgWFmVsbDagHYa6ILYGa2JwsghqLSUkbSAkmPSdoq6dI2+w+XdK+kTZLWSpqR23eWpCfSclZu+/slPZKueZ0k1fXswzkwzMyKRPqAUpWlgKQpwPXAScAcYJGkOcMOuwZYHhHvA64AvpzOPRi4HDgOmAdcLumgdM4NwGJgdloW1PHY7TgwzMxKxOBgpaXEPGBrRGyLiDeAO4BThx0zB7g3/V6T2/8x4EcR8UJEvAj8CFggaTpwQET8NCICWA4sHPsTt9d1fRgbNmzYIekXE12OMegDdkx0IRo2GZ4RJsdzdvszHj7WC/yaF1f/n1jRV/Hwd0pan1tfGhFL0+9Dgadz+wbIagx5DwOnAV8HPg7sL+mQDucempaBNtsb0XWBERHTJroMYyFpfUT0T3Q5mjQZnhEmx3NOhmcsExF1NfG061sY3vFxMfANSWcD9wHPALsKzq1yzdp0XWCYmXWpAeCw3PoMYHv+gIjYDnwCQNJ+wGkR8bKkAWD+sHPXpmvOGLb9Ldesk/swzMzGxzpgtqQjJE0FzgBW5g+Q1Cep9Xf5MmBZ+r0a+Kikg1Jn90eB1RHxLPBrScen0VFnAnc39QAOjPG3tPyQrjcZnhEmx3NOhmccFxGxC7iQ7I//z4A7I2KzpCsknZIOmw88Julx4F3AVencF4AryUJnHXBF2gZwPvBXwFbg58A9TT2DwnOkmJlZBa5hmJlZJQ4MMzOrxIHRkApTAJwg6UFJuySdPhFlHKsKz7hE0pY0zcG9ksY8Jn68VXjG89K0DBsl/aTNm7t7vLJnzB13uqSQNKmH2U5m7sNoQJoC4HHgI2TD3tYBiyJiS+6YWcABZOOuV0bEivEv6ehVfMYPAQ9ExKuSzgfmR8SnJqTAo1DxGQ+IiFfS71OAP65x3H7jqjxjOm5/4PvAVODCiFg//FrW+1zDaEbpFAAR8VREbAK6dRrMKs+4JiJeTav389bx4t2gyjO+klvdlwZfmmpIlekqIBuh81XgtfEsnO1ZHBjN6PQafy8Z6TOeS4PD/RpS6RklXSDp52R/UD87TmWrS+kzSjoaOCwi/nY8C2Z7HgdGM8b1df0JUvkZJX0G6AeubrRE9av0jBFxfUT8NnAJ8IXGS1WvwmdML5FdC/zpuJXI9lgOjGaUTgHQAyo9o6QPA38GnBIRr49T2eoy0v+Od9DgTKENKXvG/YH3AGslPQUcD6x0x/fk5MBoRukUAD2gyjQHRwM3kYXFcxNQxrGq8oyzc6u/BzwxjuWrQ+EzRsTLEdEXEbMiYhZZX9Qp7vSenBwYDagyBYCkY9OEYr8P3CRp88SVeOQqTnNwNbAf8N007LSrQrPiM14oabOkjcAS4KwOl9sjVXxGM8DDas3MrCLXMMzMrBIHhpmZVeLAMDOzShwYZmZWiQPDzMwqcWBY15P0WUk/k3TbRJfFrJd5WK11PUn/BJwUEU9OdFnMeplrGNbVJN0IHEk2XcXnJro8Zr3MNQzremmOo/6I2DHRZTHrZa5hmJlZJQ4MMzOrxIFhZmaVODDMzKwSd3qbmVklrmGYmVklDgwzM6vEgWFmZpU4MMzMrBIHhpmZVeLAMDOzShwYZmZWyf8HchGEiLFAnm4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "bicoh.plot()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (base)",
   "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": 4
}