MSR_Cog_Learn_IZH.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "rIB1Lpb1vb8y"
   },
   "source": [
    "# MSR function written in JUlia"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "id": "v_kDZ2GEtUnc"
   },
   "outputs": [],
   "source": [
    "# Import the required modules\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from collections import Counter\n",
    "import pandas as pd\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "import h5py\n",
    "file_name = \"../../Project_MSR/Codes/Ryan_Original/src/data.hdf5\"\n",
    "f1 = h5py.File(file_name,'r+')\n",
    "unfil_spikes =np.array(f1[\"Spikes_unfil\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def opt_rel(signal,dt,Tnum_points):\n",
    "    total_time=np.shape(signal)[0] * dt\n",
    "    Tbin=np.logspace(np.log10(0.01),np.log10(total_time), Tnum_points)\n",
    "    df = pd.DataFrame(columns=('Time', 'Hs','Hk','Hs+Hk'))\n",
    "\n",
    "    for count, tt in enumerate(Tbin):\n",
    "\n",
    "        mapping_ks=np.histogram(np.where(signal==1.0)[0] * dt, bins=np.arange(0.0, total_time + tt, tt))[0]\n",
    "\n",
    "        ks_counts = np.asarray(Counter(mapping_ks).most_common())\n",
    "        positive_values = np.where(ks_counts[:,0]>0)[0]\n",
    "        kq, mq = ks_counts[:,0][positive_values], ks_counts[:,1][positive_values]\n",
    "        assert np.sum(kq*mq)==np.sum(mapping_ks)\n",
    "        M = float(np.sum(kq*mq))\n",
    "\n",
    "        hk=-np.sum(((kq*mq)/M)*np.log2((kq*mq)/M))/np.log2(M)\n",
    "        hs=-np.sum(((kq*mq)/M)*np.log2(kq/M))/np.log2(M)\n",
    "        df.loc[count]=[tt,abs(hs),abs(hk),abs(hs+hk)]\n",
    "\n",
    "    return df,df.loc[df[\"Hs+Hk\"].idxmax()]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "s1,s2=opt_rel(unfil_spikes[8],0.01,100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAadElEQVR4nO3df2xd9XnH8c+DcTqHbjMt3g8upMk6lowqBRcX0LIfDVMboFrjQsWPolbrKiG00Q3GMsxWUVA74QpVZVPpMoRQNRWVoEK9dKFLJ6VbN2jWOEogDSVVBiLEdGpaEioSA/7x7I97r3Nzfb73Htv3nHt+vF8Sanzvsf09te/j732+z/f5mrsLAJB/p3V7AACAziCgA0BBENABoCAI6ABQEAR0ACiI07v1jc866yxfuXJlt749AOTS7t27f+ruA1HPdS2gr1y5UuPj49369gCQS2b2Yug5Ui4AUBAEdAAoCAI6ABQEAR0ACoKADgAF0bUqF6BIxvZM6N7tB/TysUmd3d+nTRtWa3iw0u1hoWQI6MASje2Z0B2P79Pk1IwkaeLYpO54fJ8kEdSRKgI6EFNoFn7v9gNzwbxucmpG924/kEpA590B6gjoQAyhWfj4i69o4thk5Oe8HHg8jXFJvDsoIxZFgRhCs/CHdx4Kfs7Z/X1JD6vluwOUDzN0IIbQbDt03ldfb482bVg993FzWmT9mgF957kjwTRJuzRK/fluvjtA9hDQgRjO7u8LBs8o91y1di4AR6VFvtows29Ok7RLozQ/HxovyoeUCxDDpg2r1dfbc8pjFri20t93ymw6Ki3SbHJqRrc9+rRWjWzTbY8+3TKN0u7rNb87QHkQ0IEYhgcruueqtar098lUDdo3XLpiXpCPCqZx0x8z7vLa/0aZODapdaM7Wr5TqPT3nfLuAOVCygWIaXiwMi9QDr3jbW1LBheargkxqW0wf3LksiV/H+QXAR2IELe2OyrIN9u0YXXbnHc7pvACrHTqOwPq0suLgA7o1CDYv7xXr70+ranZaghdam13/XNCVS6nmUWmWXrMNOvedoZfaQja1KWXm3kgX5e0oaEh58QiZEGcqhEpuZRG1Pfv6+05JRceyp03jynudcgvM9vt7kNRz7EoitKLU4UiJVfbHbXg2rywGVVls5AF2NDjY3smtG50h1aNbNO60R0a2zOxtJtBV5FyQWmEcstxA3WStd3tcvFRaZuFLMC6qrP3xs8hPVM8BHSUQqvgFacKJQu13UtdgG0O2N1uKobOI+WCUmgVvKLSGb09pv6+3mAKJKsa0zdRGjcoLTQ9g+xjho5SaBW84qYz8qI+k181si2y1LH+/0XonQltA/KLgI5SaBe84qQz8iZ0z/3Le+eqYZrr23t7TMffmNaqkW25/8NWRqRcUApxq0SKJJRKeu316blA7zrZk+bM5b2SS8cmp+Q6mXOn8iU/COgohTilgUUTdc9nLDt9bsNUndeeWx7xHL3V84WNRUCJhPLq9Vl6KBqYRAomI9hYBEBSeMHz7P6+louh9RTMLVv26sK7v00aJqMI6CgMdj2212otIeq5KMcmp8itZxRVLigEdj3GE6dEs/5cq2QsG5CyiRw6CoGmVJ3X7jANqfr/bxFq9/OEHDoKj12PnRcnBTNRm8lT4pgNBHQUQqvFPixOvezxzOW9sa6nxLH7YgV0M7vczA6Y2UEzG2lx3XvNbMbMPtK5IQLthRb71q8ZYKF0CYYHK9pz5wd037UXqsdCx2KfxDui7mob0M2sR9L9kq6QdL6k683s/MB1n5e0vdODBNpVsERtorn6oooe2z1BWqADhgcrmo2x3sY7ou6KU+VysaSD7v68JJnZI5I2Snq26bpPSXpM0ns7OkKUXtwKluZ+LOtGd9AetoPatRkueiuFPIiTcqlIeqnh48O1x+aYWUXShyVtbvWFzOxGMxs3s/EjR44sdKwoqVatb1thobSzotJajX1g3nL6abp1y15SW10UZ4YelThrfu91n6Tb3X3GWuTZ3P0BSQ9I1bLFmGNEyS02MNMetrNCNeySIt9Bjb/4ytxB2JQ1piNOQD8s6dyGj8+R9HLTNUOSHqkF87MkXWlm0+4+1olBonwaj4s7zUwzEfnbehvYUMCIOr2HtMDSRLUZDqW2vrrz0NzHbPRKR5yAvkvSeWa2StKEpOskfbTxAndfVf+3mX1F0r8SzLFYzTnzqGBebwN79MSUpOiAUbSDK7IqbgqL9YvktQ3o7j5tZjerWr3SI+khd99vZjfVnm+ZNwcWKipnLkk9Zpp119n9fTr+xrSOTU6d8nxUwCjiwRVZE+dM1jrWL5IVq5eLuz8h6YmmxyIDubv/8dKHhTILvehn3fXC6AclVdvALuRzkZxWB1M3Y/0iWewURebE2fXJztDsqO8BaLfxiPWL5BHQkTlxjosr45FyWTY8WNEXrrkgWNZYPyFKEjt3E0T7XGROnMVMFjyzp93P5NNj+/TwzkNzNc9UvnQe7XMBJG5sz4Ru3bI32GO9wh/k2GifC6Cr7t1+oOWBGfTZ6YyupVyeP3Jc1/7T97r17QGkKE5Z4+TUjP7668/oa98/1PZaRGOGDiBxy3rihZo3Z2YTHkmxkUNHVzVu8Wdhs7iad/+GcGRge61y6FS5oGs42Lk8mitgli/r0fE35wf39WsG0h5aoZByQdcsti0u8ml4sKInRy7TC6MfVP/yZZHXfOc52movBTN0pKY5vRJaKGP7fvHRqz4ZzNCRinp6pfE4uNBGcbbvFx+tG5JBQEcqotIrrvmnp7B9vxxo3ZAMUi5IReittKta2UCVS7nQuiEZBHSkIpQzp0ytvJp71Y/tmWh5AhXaI+WCVPAWG61ErbHcsmWvLrz727QDWABm6EgFb7HRSuiUqmOTU+xNWAACOlLDcXAIaVWuyFmk8ZFyAdB17coVJ45NknqJgYAOoOui1lia0V63PQI6gK6rn0t65vLe4DW0hWiPgA4gE4YHK9pz5wd037UXBq+hNUBrBHQAmTI8WFGF1gCLQkAHkDmhnPrxN6bJo7dA2SKAzKmXKN79zf06emJq7nHq0ltjho4lq2/ZXjWyTetGdzCDQkcMD1a0fNn8Oefk1Ixue/Rpfs8iMEPHgjT3NF+/ZkCP7Z7g1CEkIrQIOuOuW7bs1V1b9+uuD72L37UazhRFbFHnQpqqHROb9Zhp1p0t/liSdaM7ggeh1PX19uieq9aW5nes1ZmipFwQW6ineZQZ97kmS2wIwWLF2XBEffpJBHTEttgaYF5wWKz6hqMeC51vVUV9ehUBHbGFaoBbv9SqeMFhsYYHK/rCNRe0nKlTn15FQEdsoZ7mN1y6QpX+PpkUnEnxgsNStGoNQF/9k6hyQWxxeppHLZzygkMn1Nsv1yutJo5NqsfslJReWRZGQ2IFdDO7XNLfS+qR9KC7jzY9v1HSZyXNSpqWdIu7/3eHx4oMaNfTnIMskLT671LjxIFy2aq2ZYtm1iPpR5LeL+mwpF2Srnf3Zxuueauk4+7uZvZuSY+6+5pWX5eyxWxrrjcnKCNLQuWMZTijtlXZYpwZ+sWSDrr787Uv9oikjZLmArq7v9Zw/RkKV7MhB5rTJsx+kDWhRfayL77HWRStSHqp4ePDtcdOYWYfNrPnJG2T9CdRX8jMbjSzcTMbP3LkyGLGixRE1ZtTeogsCS2yl33xPU5AjypbmDcDd/dv1NIsw6rm0+d/kvsD7j7k7kMDAwMLGijSw+wHWRequCr74nucgH5Y0rkNH58j6eXQxe7+XUnvNLOzljg2dAmzH2RdvYyxXi5b6e/T1RdVdO/2A6VuEhcnoO+SdJ6ZrTKzZZKuk7S18QIz+02zagGymb1H0jJJP+v0YJEOZj/Ig+HBip4cuUwvjH5Qmzas1mO7JzRxbHKu5cStW/bq02P7uj3MVLVdFHX3aTO7WdJ2VcsWH3L3/WZ2U+35zZKulvRxM5uSNCnpWu9W1y8sGaWHyJtQn6GHdx7S0DveVprfXbotAsi9VSPbgqV1RStlXGrZIgqOmnPk3dn9fcE2u2VazKeXS8nVa84bc4+0u0XebNqwOtgkrkyL+QT0kqPmHEUwPFjRDZeuiAzqZTpYmoBectScoyg+N7xWX7z2wnkdGesHS5chqBPQS46acxRJq4Oly/Cuk4BectSco2jK/K6TgF5yUTvuynTgLoqnzO86KVtE2x7nQJ5s2rB63iEr0snF0SL/rhPQARRKPWDf/c39Onpiau7x+uJo4zVFQ8qlJMb2TGjd6I5SNy5CeZR1cZQZeglwYAXKqIyLo8zQS4DNQyijMi6OEtALbmzPBD0uUEpRJbmm6jvUoqYdSbkUWD3VElLkmQrQ2AZ64tikTCePWitq2pEZeoFFpVrq2DyEMqgfglHp75vXXreIaUcCeoG1SqmweQhlUpYFUgJ6gYVSKpX+PoI5SqUsC6QE9AKjTwtQFfVakIrXWpdF0QLjbFCgqiy7RzlTFEBprBvdEVnGm6dzRzlTtAQ4FxRoL7QIGtqrkTfk0AuAc0GBeEKLoCYV4vVCQC8AtvYD8YQOk3apEK8XAnoBlKXGFliq4cHKvA1GdUV4vRDQC6AsNbZAJ1QK/HohoOdYvcd5vU9FI+rNgWhFbtpFlUtONfc4d2mu+VCFKhcgqMhNu5ih51TUQmg9mD85clkufxmBtBS1aRcBPadYCAWWrmivIwJ6TrEQCixd0V5HBPScovEWsHRFex2xKJpTNN4Clq5oryOacwFATR56IrVqzkXKBQAU3RPp1i179emx8Lm8WRMroJvZ5WZ2wMwOmtlIxPM3mNkztf+eMrMLOj9UAEhOqBT44Z2HcrPZqG0O3cx6JN0v6f2SDkvaZWZb3f3ZhstekPQH7n7UzK6Q9ICkS5IYcNnk4S0gUAShUsV64648vO7izNAvlnTQ3Z939zclPSJpY+MF7v6Uux+tfbhT0jmdHWY50RYXSE+rUsW81KXHCegVSS81fHy49ljIJyV9K+oJM7vRzMbNbPzIkSPxR1lStMUF0hNqrSvlpy49TkAPtQ+ef6HZelUD+u1Rz7v7A+4+5O5DAwMD8UdZUkXbxQZk2fBgRTdcuiLXje7iBPTDks5t+PgcSS83X2Rm75b0oKSN7v6zzgyv3Iq2iw3Ius8Nr9UXr71Qlf4+maq9ke65am0u8udSvIC+S9J5ZrbKzJZJuk7S1sYLzGyFpMclfczdf9T5YZZHvSXuqpFtOvHmtHpPO3W+kKfZApBH9cZdX7z2QknSrVv25qatbtsqF3efNrObJW2X1CPpIXffb2Y31Z7fLOlOSW+X9GUzk6TpUOE7wppb4h49MaXeHlN/X69enZyiygVISfNrMS9tddkpmiH1wyqa1VviAkhHll+L7BTNCRZBgWwIveYmjk1mOvVCQM8QFkGBbGj1msvyXhACeoYUrZUnkFdRr8W6LO8FoX1uhhStlSeQV/XX3C1b9kY+n9U0KAE9Y4YHKwRwIAOGBytzB0k3y2oalJQLAATkLQ3KDB0AAvKWBqUOHQBiyEor61Z16MzQAaCNvOwcJaCnLCt/5QHE16qVdZZevwT0FOXlrzyAU+VlFzdVLikZ2zOh2x59mgMrgBzKyy5uAnoK6jPzmcACdNb+ygM4VV7KF0m5pCAq/9Yoa3/lAZwqL+WLzNBT0GoGnsW/8gDmy8PBF8zQE1SvaAlV+veY5ep4K6Dssl7YwAw9IfUffFQfCKk6M//CNRdk4pcAQDytyhezgBl6QlrlzSsZzb8BaC3r5YvM0BMS+gGbpCdHLiOYAzkUKmA4zSwTuXQCekLyUrcKIL7QwRcz7pk4yYiAnpC81K0CiG94sKJ7rlqrHrN5z2Uhl05AT0j9B1/p75OpmjenogXIv+HBimYzukmQRdEEcfoQUExn9/dFVrD1L+/twmhOYoYOAAu0acNq9fbMT7u89vp0V/PoBHQAWKDhwYrOWDY/wTE1613NoxPQAWARjk1ORT4e2kyYBgI6ACxCVKVLq8fTwKLoInHyEFBuoXbYocfTwAx9ERr7tLhONujp9qYCAOmpBDYJntnFShcC+iJkvUEPgORlsdKFgL4IWW/QAyB5rSpdbnv06a4EdQL6ItCnBYAkvRqodOlWbxcC+iLQpwWA1HoS1400bKyAbmaXm9kBMztoZiMRz68xs++Z2Rtm9ledH2a20KcFgBTuvliXdhq2bdmimfVIul/S+yUdlrTLzLa6+7MNl70i6c8lDScxyCyiTwuAegz4y0f3ajaiWjHt3i5xZugXSzro7s+7+5uSHpG0sfECd/+Ju++SFJ1QAoCCGh6s6Jd+ITpwp12SHmdjUUXSSw0fH5Z0yWK+mZndKOlGSVqxYsVivkTq2EAEoJ3Q4mjo8aTEmaFH7WNd1N8dd3/A3YfcfWhgYGAxXyJVbCACEEdocTSLKZfDks5t+PgcSS8nM5zsGNszodsefZoNRADaysomozgBfZek88xslZktk3SdpK3JDqu76jPzUE8GNhABaJSVdrptc+juPm1mN0vaLqlH0kPuvt/Mbqo9v9nMfk3SuKRfkjRrZrdIOt/df57c0JNz19b982bmjdhABKBZKF+eZjvdWN0W3f0JSU80Pba54d//p2oqJvfG9kwE+xxLbCACEC10LJ2pGlfSKKZgp2iTVm+PeszYQAQg0qYNq4MVJGmlXQjoTVrlx79wzQUEcwCRhgcrwfK/tNbdCOhNQvnxM5f3EswBtBTqkZ5W+SIBvWZsz4TWje7QxLHJeW+b+np79Jk/eldXxgUgP0Lli6+emEqlfJGArlM3EEnVnFf9R0LjLQBxDQ9WdPpp8wP6rKrVc0njTFFFn0DkqgbzJ0cu686gAOTS5NRs5OOtquc6hRm6OIEIQDEQ0MUJRAA6p9Uh0Unn0Qno4gQiAJ3TqoAi6Xp0Aro4gQhA57SKG0mncVkUreEEIgCdUgm0AUg6jcsMHQA6LCqNa5LWr0n2HIhSBfT65qFVI9u0bnQHB1UASMTwYEVXX3TqO36XtGXXS4nGndIEdE4fApCmbc/8eN5jUzOuv/3GvsS+Z2kCetTmIU4fApCUoyeiNxIdf3MmsYlkaQI6m4cAZEVSE8nSBHQ2DwFIU39feINRUhPJ0gR0Ng8BSNNdHwpvMEqqnW5pAjqbhwCkaXiwor7e6BAbOH9+yUq1sYjNQwDS9Hqg82LoQOmlKs0MHQDSFlqjI+UCADkTOsHotdenEyldJKADQEKGBys6Y9n8zPbUrCdSukhAB4AEhU4qimretVQEdAAoCAI6AHRJp/PoBHQASFCPzV8Uret0Hp2ADgAJuv6Sc4PPdboFAAEdABL0ueG1CmwY7Xg9OgEdABJ2xluiA3enWwAQ0AEgYaHSxdDji0VAB4CEhRZGWy2YLgYBHQASNhPIrYQeXywCOgAkrBJo0hV6fLFiBXQzu9zMDpjZQTMbiXjezOwfas8/Y2bv6egoa8b2TGjd6A6tGtmmdaM7OOAZQC5EHbBjktavGejo92kb0M2sR9L9kq6QdL6k683s/KbLrpB0Xu2/GyX9Y0dHqWowv+PxfZo4NilXtQ/CHY/vI6gDyLzhwYquvqiixoy5S3ps90RHY1icGfrFkg66+/Pu/qakRyRtbLpmo6R/9qqdkvrN7Nc7NkpVd1RNTs2c8tjk1Exih60CQCd957kjas6YdzqGxQnoFUkvNXx8uPbYQq+Rmd1oZuNmNn7kyJEFDTS0oyqpw1YBoJPSiGFxAnpUXU3zH5o418jdH3D3IXcfGhhYWO4odPJH6HEAyJI0YlicgH5YUmMzgnMkvbyIa5YkalGhr7dHmzas7uS3AYBEpBHD4gT0XZLOM7NVZrZM0nWStjZds1XSx2vVLpdKetXdf9yxUaq6qHDPVWtV6e+TqVruc89Vazn0GUAupBHDzGMUtpvZlZLuk9Qj6SF3/zszu0mS3H2zmZmkL0m6XNIJSZ9w9/FWX3NoaMjHx1teAgBoYma73X0o6rn5h91FcPcnJD3R9Njmhn+7pD9byiABAEvDTlEAKAgCOgAUBAEdAAqCgA4ABRGryiWRb2x2RNKLi/z0syT9tIPDyQPuuRy453JYyj2/w90jd2Z2LaAvhZmNh8p2iop7LgfuuRySumdSLgBQEAR0ACiIvAb0B7o9gC7gnsuBey6HRO45lzl0AMB8eZ2hAwCaENABoCAyHdCzcjh1mmLc8w21e33GzJ4yswu6Mc5OanfPDde918xmzOwjaY4vCXHu2czeZ2Z7zWy/mf1n2mPstBi/279sZt80s6dr9/yJboyzU8zsITP7iZn9IPB85+OXu2fyP1Vb9f6vpN+QtEzS05LOb7rmSknfUvXEpEsl/U+3x53CPf+OpDNr/76iDPfccN0OVbt+fqTb407h59wv6VlJK2of/0q3x53CPf+NpM/X/j0g6RVJy7o99iXc8+9Leo+kHwSe73j8yvIMPROHU6es7T27+1PufrT24U5VT4fKszg/Z0n6lKTHJP0kzcElJM49f1TS4+5+SJLcPe/3HeeeXdIv1s5XeKuqAX063WF2jrt/V9V7COl4/MpyQO/Y4dQ5stD7+aSqf+HzrO09m1lF0oclbVYxxPk5/5akM83sP8xst5l9PLXRJSPOPX9J0m+renzlPkl/4e6z6QyvKzoev2IdcNElHTucOkdi34+ZrVc1oP9uoiNKXpx7vk/S7e4+U5285V6cez5d0kWS/lBSn6TvmdlOd/9R0oNLSJx73iBpr6TLJL1T0r+b2X+5+88THlu3dDx+ZTmgZ+Jw6pTFuh8ze7ekByVd4e4/S2lsSYlzz0OSHqkF87MkXWlm0+4+lsoIOy/u7/ZP3f24pONm9l1JF0jKa0CPc8+fkDTq1QTzQTN7QdIaSd9PZ4ip63j8ynLKJROHU6es7T2b2QpJj0v6WI5na43a3rO7r3L3le6+UtLXJf1pjoO5FO93+18k/Z6ZnW5myyVdIumHKY+zk+Lc8yFV35HIzH5V0mpJz6c6ynR1PH5ldobu7tNmdrOk7Tp5OPX+xsOpVa14uFLSQdUOp+7WeDsh5j3fKentkr5cm7FOe4471cW850KJc8/u/kMz+zdJz0ialfSgu0eWv+VBzJ/zZyV9xcz2qZqOuN3dc9tW18y+Jul9ks4ys8OSPiOpV0oufrH1HwAKIsspFwDAAhDQAaAgCOgAUBAEdAAoCAI6ABQEAR0ACoKADgAF8f9SbHp/peLZYAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(s1.Hs.values,s1.Hk.values)\n",
    "plt.axhline(s2.Hk)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Time     1.843933\n",
       "Hs       0.798917\n",
       "Hk       0.371052\n",
       "Hs+Hk    1.169968\n",
       "Name: 44, dtype: float64"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def opt_relv(signal,dt,Tnum_points):\n",
    "    total_time=np.shape(signal)[0] * dt\n",
    "    Tbin=np.logspace(np.log10(0.01),np.log10(total_time), Tnum_points)\n",
    "    df = pd.DataFrame(columns=('Time', 'Hs','Hk','Hs+Hk'))\n",
    "\n",
    "    for count, tt in enumerate(Tbin):\n",
    "\n",
    "        mapping_ks=np.histogram(np.where(signal==1.0)[0] * dt, bins=np.arange(0.0, total_time + tt, tt))[0]\n",
    "\n",
    "        ks_counts = np.asarray(Counter(mapping_ks).most_common())\n",
    "        positive_values = np.where(ks_counts[:,0]>0)[0]\n",
    "        kq, mq = ks_counts[:,0][positive_values], ks_counts[:,1][positive_values]\n",
    "        assert np.sum(kq*mq)==np.sum(mapping_ks)\n",
    "        M = float(np.sum(kq*mq))\n",
    "\n",
    "        hk=-np.sum(((kq*mq)/M)*np.log2((kq*mq)/M))/np.log2(M)\n",
    "        hs=-np.sum(((kq*mq)/M)*np.log2(kq/M))/np.log2(M)\n",
    "        df.loc[count]=[tt,abs(hs),abs(hk),abs(hs+hk)]\n",
    "\n",
    "    return df.loc[df[\"Hs+Hk\"].idxmax()]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "fzzabiedvsZe"
   },
   "source": [
    "# Seting up Izhikevich  Parameters\n",
    "[Simple Model of Spiking Neurons; IEEE Transactions on Neural Networks (2003) 14:1569- 1572 ](https://www.izhikevich.org/publications/spikes.htm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "id": "-1GIvLQRteYE"
   },
   "outputs": [],
   "source": [
    "\n",
    "pars={#                             a            b       c        d         I      V0  \n",
    "    \"tonic spiking\" :              [0.02,      0.2,     -65,      6,       14,   -70.0], \n",
    "    #------------------------------------------------------------------------------------#\n",
    "    \"phasic spiking\" :             [0.02,      0.25,    -65,      6,      0.5,   -64.0],\n",
    "    #------------------------------------------------------------------------------------#\n",
    "    \"tonic bursting\" :             [0.02,      0.2,     -50,      2,       15,   -70.0], \n",
    "    #------------------------------------------------------------------------------------#\n",
    "    \"phasic bursting\" :             [0.02,      0.25,    -55,     0.05,    0.6,   -64.0], \n",
    "    #------------------------------------------------------------------------------------#\n",
    "    \"mixed mode\" :                 [0.02,      0.2,     -55,     4,        10,   -70.0],\n",
    "    #------------------------------------------------------------------------------------#\n",
    "    \"spike frequency adaptation\" : [0.01,      0.2,     -65,     8,        30,   -70.0],\n",
    "    #------------------------------------------------------------------------------------#\n",
    "    \"Class 1\" :                    [0.02,     -0.1,     -55,     6,         0,   -60.0],\n",
    "    #------------------------------------------------------------------------------------#\n",
    "    \"Class 2\" :                    [0.2,       0.26,    -65,     0,         0,   -64.0],\n",
    "    #------------------------------------------------------------------------------------#\n",
    "    \"spike latency\" :              [0.02,      0.2,     -65,     6,         7,   -70.0],\n",
    "    #------------------------------------------------------------------------------------#\n",
    "    \"subthreshold oscillations\" :  [0.05,      0.26,    -60,     0,         0,   -62.0],\n",
    "    #------------------------------------------------------------------------------------#\n",
    "    \"resonator\" :                  [0.1,       0.26,    -60,    -1,         0,   -62.0],\n",
    "    #------------------------------------------------------------------------------------#\n",
    "    \"integrator\" :                 [0.02,      -0.1,    -55,     6,         0,   -60.0],\n",
    "    #------------------------------------------------------------------------------------#\n",
    "    \"rebound spike\" :              [0.03,      0.25,    -60,     4,         0,   -64.0],\n",
    "    #------------------------------------------------------------------------------------#\n",
    "    \"rebound burst\" :              [0.03,      0.25,    -52,     0,         0,   -64.0],\n",
    "    #------------------------------------------------------------------------------------#\n",
    "    \"threshold variability\" :      [0.03,      0.25,    -60,     4,         0,   -64.0],\n",
    "    #------------------------------------------------------------------------------------#\n",
    "    \"bistability\" :                [1,         1.5,     -60,     0,         0,   -61.0],\n",
    "    #------------------------------------------------------------------------------------#\n",
    "    \"DAP\" :                        [1,         0.2,     -60,    -21,        0,   -70.0],\n",
    "    #------------------------------------------------------------------------------------#\n",
    "    \"accomodation\" :               [0.02,      1,       -55,     4,         0,   -65.0],\n",
    "    #------------------------------------------------------------------------------------#\n",
    "    \"inhibition-induced spiking\" : [-0.02,     -1,      -60,     8,        80,   -63.8],\n",
    "    #------------------------------------------------------------------------------------#\n",
    "    \"inhibition-induced bursting\": [-0.026,    -1,      -45,     0,        80,   -63.8]\n",
    "    #------------------------------------------------------------------------------------#\n",
    "    };\n",
    "name=[\n",
    "    \"tonic spiking\",#1\n",
    "    \"phasic spiking\",#2\n",
    "    \"tonic bursting\",#3\n",
    "    \"phasic bursting\",#4\n",
    "    \"mixed mode\",#5\n",
    "    \"spike frequency adaptation\",#6\n",
    "    \"Class 1\",#7\n",
    "    \"Class 2\",#8\n",
    "    \"spike latency\",#9\n",
    "    \"subthreshold oscillations\",#10\n",
    "    \"resonator\",#11\n",
    "    \"integrator\",#12\n",
    "    \"rebound spike\",#13\n",
    "    \"rebound burst\",#14\n",
    "    \"threshold variability\",#15\n",
    "    \"bistability\",#16\n",
    "    \"DAP\",#17\n",
    "    \"accomodation\",#18\n",
    "    \"inhibition-induced spiking\",#19\n",
    "    \"inhibition-induced bursting\"#20\n",
    "     ];"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "uBO60vHwwVoT"
   },
   "source": [
    "# Importing MNIST DataSet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "id": "ADzgG2Votu6D"
   },
   "outputs": [],
   "source": [
    "from keras.datasets import mnist\n",
    "\n",
    "(train_images, train_labels), (test_images, test_labels) = mnist.load_data()\n",
    "\n",
    "train_images = train_images.reshape((60000, 28 * 28))\n",
    "train_images = train_images.astype('float32') / 255\n",
    "test_images = test_images.reshape((10000, 28 * 28))\n",
    "test_images = test_images.astype('float32') / 255"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "VfxPDELqOAiw",
    "outputId": "3e66d23e-7ba6-4bf7-83bc-0df228afb784"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(6742, 784) (5851, 784)\n"
     ]
    }
   ],
   "source": [
    "D0=train_images[np.where(train_labels==0)]\n",
    "D1=train_images[np.where(train_labels==1)]\n",
    "D2=train_images[np.where(train_labels==2)]\n",
    "D3=train_images[np.where(train_labels==3)]\n",
    "D4=train_images[np.where(train_labels==4)]\n",
    "D5=train_images[np.where(train_labels==5)]\n",
    "D6=train_images[np.where(train_labels==6)]\n",
    "D7=train_images[np.where(train_labels==7)]\n",
    "D8=train_images[np.where(train_labels==8)]\n",
    "D9=train_images[np.where(train_labels==9)]\n",
    "print(np.shape(D1),np.shape(D8))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "LGn-NcdPOCTw",
    "outputId": "c48d2b32-e48b-4d60-8ddb-1ae6c0f8c5fd"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7fe85020f1c0>"
      ]
     },
     "execution_count": 90,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAS4AAAD8CAYAAADJwUnTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAVz0lEQVR4nO3de6xV5ZnH8e8PxFslox1QCcLAGBRxtOIcpal1vLaCrWLj1IFOW2vqUCcyHZI21fjHaDqZxlHLVKKVHBXFaJUmYsUGq0ZGbeNoOaByEWlPEPEI5TLWqr0Z5Jk/9sbZ7HP22vvs29rv4fdJds5Z61nrXQ8r8PCud71rbUUEZmYpGZZ3AmZmg+XCZWbJceEys+S4cJlZcly4zCw5LlxmlhwXLjNrGUmLJO2QtK5CXJIWSOqVtEbSqbW068JlZq10LzA9Iz4DmFT8zAHuqKVRFy4za5mIeA54O2OTmcB9UfACcLikMdXaPaBZCdZi1KhRMWHChHYe0my/snnzZnbt2qVG2pA0mMdp1gN/KlnujojuQew/FnizZLmvuG5b1k4NFS5J04FbgeHAXRFxY9b2EyZMoKenp5FDmlmGrq6udh/yTxHRyEEHKrJVC2fdl4qShgO3U7hGnQLMljSl3vbMrHNIqunTBH3AuJLlY4Ct1XZqZIzrdKA3IjZFxAfAQxSuV80sccOGDavp0wTLgK8W7y5+EvhdRGReJkJjl4oDXZtOK99I0hwKdwsYP358A4czs3ZpUm8KSQ8CZwOjJPUB1wMjACJiIbAcuBDoBf4AXFFLu40UrpquTYsDdd0AXV1dfoeOWYdr4mUgETG7SjyAqwfbbiOFq65rUzPrfM0qXK3SyEXqSmCSpImSDgRmUbheNbPEtXFwvi5197giYrekucATFKZDLIqI9U3LzMxy0+k9robmcUXEcgqDa2Y2REhq1h3DlmnrzHkzS8OQ7nGZ2dDkwmVmyXHhMrPkuHCZWVI8OG9mSXKPy8yS48JlZslx4TKzpOT9OE8tXLjMrB8XLjNLju8qmlly3OMys6R4jMvMkuTCZWbJceEys+R4cN7MkuIxLjNLkguXmSXHhcvMkuPCZWbJceEys6T4RYJmliT3uMwsOS5cZpYcFy4zS4onoJpZkly4zCw5vqtoZslxj6tN3n///cz4kiVLMuMHHXRQZnz16tUVY++9917mvvfff39m/JxzzsmMjx07NjPeSkcffXRmfObMmZnxrq6uZqZjbTDkx7gkbQbeAz4EdkeE/5aaDQGdXriacSF7TkSc4qJlNnTs7XVV+9TY1nRJGyX1Srp2gPhfSHpM0iuS1ku6olqbQ+ZS0cyap1mD85KGA7cDnwH6gJWSlkXEqyWbXQ28GhEXSRoNbJT0QER8UDG/BvMK4ElJqyTNqZD4HEk9knp27tzZ4OHMrNVq7W3V2OM6HeiNiE3FQvQQUD4wGsBIFRo8DHgb2J3VaKM9rjMiYqukI4GnJL0WEc/tk1FEN9AN0NXVFQ0ez8zaYBBjXKMk9ZQsdxf/ze81FnizZLkPmFbWxm3AMmArMBL4h4jYk3XQhgpXRGwt/twh6REK1fW57L3MrNMNonDtqjK+PVBD5R2YC4CXgXOBYyl0gn4eEe9WarTuS0VJH5M0cu/vwGeBdfW2Z2ado4mXin3AuJLlYyj0rEpdASyNgl7gdWByVqON9LiOAh4pJn8A8KOI+FkD7TXku9/9bmb85ptvblMmg/f444/nnULdvve972XGTzzxxIqxWbNmZe47e/bszPjEiRMz41a/Jk6HWAlMkjQReAuYBXypbJstwHnAzyUdBRwPbMpqtO7CFRGbgE/Uu7+ZdaZmvkgwInZLmgs8AQwHFkXEeklXFeMLgX8H7pW0lsKl5TURsSurXU+HMLN+mjkBNSKWA8vL1i0s+X0rhaGmmrlwmVk/nT5z3oXLzPpx4TKzpAz5h6zNbGhy4WqThx9+OLdjjxo1KjN+0kkntSmT/iZPzpwOw2uvvZYZf+eddzLjL730UmZ87dq1dcUATj755My4p0O0jl8kaGbJcY/LzJLiMS4zS5ILl5klx4XLzJLjwXkzS4rHuMwsSS5cbfLkk09mxjdu3JgZP/744+s+9qGHHpoZHzNmTN1t563aV69Vm6P2xhtv1H3sxx57LDP++c9/vu62LZsLl5klx4XLzJLjwmVmSWnmiwRbxYXLzPpxj8vMkuPCZWbJceEys6R4AmobHXvssQ3FbWDV5lI1Mk/r4IMPzoxfeeWVdbdtjXHhMrPk+K6imSXFl4pmliQXLjNLjguXmSXHhcvMkuJHfswsSe5xWa4++OCDzPg3v/nNzPjixYubmc4+nn/++cz41KlTW3Zsy9bphatqf1DSIkk7JK0rWfdxSU9J+nXx5xGtTdPM2mnvlIhqn7zUciF7LzC9bN21wNMRMQl4urhsZkNE8oUrIp4D3i5bPRPYew2xGLikuWmZWV5qLVp5Fq56x7iOiohtABGxTdKRlTaUNAeYAzB+/Pg6D2dm7dTpdxVbnl1EdEdEV0R0jR49utWHM7Mm6PQeV72Fa7ukMQDFnzual5KZ5a2ZhUvSdEkbJfVKGnA8XNLZkl6WtF7Ss9XarLdwLQMuL/5+OfBone2YWYdp5hiXpOHA7cAMYAowW9KUsm0OB34IXBwRJwJfrNZu1TEuSQ8CZwOjJPUB1wM3Aj+W9HVgSy0HstZZsWJFxdj999+fue8999zT0LEPPPDAzPiCBQsqxk444YSGjm2t08TLwNOB3ojYVGz3IQo3914t2eZLwNKI2AIQEVWv4KoWroiYXSF0XrV9zSxNgxicHyWpp2S5OyK6S5bHAm+WLPcB08raOA4YIekZYCRwa0Tcl3VQz5w3s34G0ePaFRFdWU0NsC7Klg8A/pZCZ+gQ4H8kvRARv6rUqAuXme2jyXcM+4BxJcvHAFsH2GZXRPwe+L2k54BPABULV2dP1jCzXDTxruJKYJKkiZIOBGZRuLlX6lHgTEkHSDqUwqXkhqxG3eMys36a1eOKiN2S5gJPAMOBRRGxXtJVxfjCiNgg6WfAGmAPcFdErKvcqguXmQ2gmZNLI2I5sLxs3cKy5ZuBm2tt04UrAb/85S8z4xdccEHF2O7du5udzj6q/QUfN25cxdjw4cObnY41gV8kaGZJ6vT3cblwmVk/LlxmlhwXLjNLjguXmSUl71fW1MKFy8z68V1FM0uOe1zWsCVLlmTGWz1XK8uf//znzPjnPve5irHTTjstc9+LLrooM37JJZdkxk866aTMuFXmwmVmSfEYl5klyYXLzJLjwXkzS457XGaWFI9xmVmSXLjMLDkuXNawSy+9NDO+YUPlt9z29PRUjAHs3LmzrpyaYeXKlQ3Fb7jhhsz4vHnzKsauueaazH2PPPLIzPhQ58JlZknxiwTNLEnucZlZcly4zCw5LlxmlhwXLjNLiiegmlmSfFfRGvapT30qM758+fKKsS1btmTuu2vXrsz49u3bM+NLly7NjN99990VYxGRuW81e/bsyYzPnz+/Ymz16tWZ+z799NOZ8U7/h92oTu9xVT37khZJ2iFpXcm6GyS9Jenl4ufC1qZpZu2093Kx2icvtfy3cS8wfYD1/xURpxQ/lf/LN7Ok1Fq08ixcVS8VI+I5SRPakIuZdYjkLxUzzJW0pngpeUSljSTNkdQjqSfP5+LMrHbDhg2r6ZNbfnXudwdwLHAKsA34fqUNI6I7Iroiomv06NF1Hs7M2in5S8WBRMRHt5ok3Qn8tGkZmVmu8i5KtairxyVpTMniF4B1lbY1s/Qk3+OS9CBwNjBKUh9wPXC2pFOAADYD32hditaI8ePHNxSvZsaMGZnxs846q2Lstttuy9z3xRdfrCunWjzzzDOZ8VtuuSUz/p3vfKeJ2XSeTu9x1XJXcfYAqyvPKjSz5CVfuMxs/5LCiwQ7Ozszy0Uzx7gkTZe0UVKvpGsztjtN0oeS/r5amy5cZtZPswqXpOHA7cAMYAowW9KUCtv9J/BELfm5cJlZP03scZ0O9EbEpoj4AHgImDnAdv8CPAzsqKVRFy4z62cQhWvU3idjip85ZU2NBd4sWe4rris91lgK06oW1pqfB+etpb785S9XjM2aNStz3/PPPz8z/uyzz9aVUy16e3tb1nanG+QcrV0R0ZXV3ADryt9n9APgmoj4sNbjunCZWT9NvKvYB4wrWT4G2Fq2TRfw0N4eHHChpN0R8ZNKjbpwmVk/TZzHtRKYJGki8BYwC/hS6QYRMbHkuPcCP80qWuDCZWYDaFbhiojdkuZSuFs4HFgUEeslXVWM1zyuVcqFy8z20eznEIsvGl1etm7AghURX6ulTRcuM+vHj/yYWXI6/ZEfFy4z20fer6yphQuX5eaAA7L/+p166qmZ8VbO4zruuONa1nYKXLjMLDkuXGaWHBcuM0uOC5eZJSWFFwm6cJlZP+5xmVlyXLjMLDkuXMa2bdsy43feeWdmfPLkyZnxyy67bNA5dYIPP/wwM/7KK6+07NgjRozIjE+bNq1lx+50noBqZkny4LyZJcc9LjNLjguXmSXFY1xmliQXLjNLjguXmSUn+buKksYB9wFHA3uA7oi4VdLHgSXABGAzcFlE/LZ1qXau3/zmN5nx6dOnZ8bXrFmTGX/nnXcGm1LH2L59e8XY/PnzM/ddsWJFs9P5yAknnJAZP/PMM1t27E6XwhhXLWV1N/CtiDgB+CRwtaQpwLXA0xExCXi6uGxmQ8Agvsk6F1ULV0Rsi4jVxd/fAzZQ+ArtmcDi4maLgUtalKOZtVmnF65BjXFJmgBMBV4EjoqIbVAobpKObH56ZpaHTr9UrLlwSToMeBiYFxHv1voHkzQHmAMwfvz4enI0szbr9MJV060DSSMoFK0HImJpcfV2SWOK8THAjoH2jYjuiOiKiK7Ro0c3I2cza6G9LxKs5ZOXqkdWofTeDWyIiNLbQMuAy4u/Xw482vz0zCwPQ2GM6wzgK8BaSS8X110H3Aj8WNLXgS3AF1uSYQLmzZuXGa823aGa119/PTN+/PHHV4wdcsghDR37j3/8Y2b8pptuyoxnTXl4991368qpViNHjqwYW7BgQUuPnbpOv1SsWrgi4hdApT/Fec1Nx8w6QfKFy8z2L3lfBtbChcvM+kn+kR8z2/+4x2VmyXHhMrOkeIzLzJLkwrUfOO+87FkhS5Ysaaj9qVOn1h0//PDDGzp2tVfqvPTSSw2134iseVoAjzzySMXYWWed1ex0hpRmFi5J04FbgeHAXRFxY1n8H4FriovvA/8cEZnfTefCZWb9NOuuoqThwO3AZ4A+YKWkZRHxaslmrwNnRcRvJc0AuoHML7Z04TKzfTR5jOt0oDciNhXbfojCK7E+KlwR8XzJ9i8Ax1Rr1IXLzPoZROEaJamnZLk7IrpLlscCb5Ys95Hdm/o68Hi1g7pwmVk/gyhcuyKiK6upAdZFhWOeQ6FwfbraQV24zKyfJl4q9gHjSpaPAbYOcLyTgbuAGRHxv9Ua7ex5/WaWiya+1mYlMEnSREkHArMovBKr9FjjgaXAVyLiV7U06h6Xme1j74sEmyEidkuaCzxBYTrEoohYL+mqYnwh8G/AXwI/LBbD3VUuP124muH888/PjM+ePTsz/uCDDzZ0/DznUjVixIgRmfFq7zm79NJLM+PTpmXeUbcMzZzHFRHLgeVl6xaW/H4lcOVg2nThMrN+PHPezJLjwmVmSfFD1maWJL9I0MyS4x6XmSXHhcvMkuIxrv3ExIkTM+P33HNPZvziiy/OjK9YsSIzftxxx1WMLVu2rGKsFpMnT25o/3PPPbdiLOv7IKH6e8isdVy4zCw5LlxmlhzfVTSzpHiMy8yS5MJlZslx4TKz5LhwmVlyki9cksYB9wFHA3sovAz/Vkk3AP8E7Cxuel3xvTtW5qCDDsqMz5o1q6F4lm9/+9t172v7p2a+SLBVaulx7Qa+FRGrJY0EVkl6qhj7r4i4pXXpmVkeku9xRcQ2YFvx9/ckbaDwlUNmNkR1euEaVH9Q0gRgKvBicdVcSWskLZJ0RIV95kjqkdSzc+fOgTYxsw7TxC/LaImaC5ekw4CHgXkR8S5wB3AscAqFHtn3B9ovIrojoisiukaPHt14xmbWUrUWrTwLV013FSWNoFC0HoiIpQARsb0kfifw05ZkaGZt1+mD81WzU6Gs3g1siIj5JevHlGz2BWBd89MzszwMhR7XGcBXgLWSXi6uuw6YLekUCl+nvRn4RgvyM7McdPrgfC13FX8BDPSn8JwtsyEo795ULTxz3sz6ceEys+S4cJlZUobKIz9mtp9xj8vMkuPCZWbJceEys+S4cJlZUjyPy8yS5LuKZpYc97jMLDmdXrg6uz9oZm3X7PdxSZouaaOkXknXDhCXpAXF+BpJp1Zr04XLzPppVuGSNBy4HZgBTKHwVpkpZZvNACYVP3MovKQ0kwuXmfUzbNiwmj41OB3ojYhNEfEB8BAws2ybmcB9UfACcHjZ+/76aesY16pVq3ZJeqNk1ShgVztzGIROza1T8wLnVq9m5vZXjTawatWqJySNqnHzgyX1lCx3R0R3yfJY4M2S5T5gWlkbA20zluKX9AykrYUrIvZ56byknojoamcOterU3Do1L3Bu9eq03CJiehObG+h6MurYZh++VDSzVuoDxpUsHwNsrWObfbhwmVkrrQQmSZoo6UBgFrCsbJtlwFeLdxc/Cfyu+H2uFeU9j6u7+ia56dTcOjUvcG716uTcGhIRuyXNBZ4AhgOLImK9pKuK8YUUXgN/IdAL/AG4olq7isi8lDQz6zi+VDSz5LhwmVlycilc1R4ByJOkzZLWSnq5bH5KHrkskrRD0rqSdR+X9JSkXxd/HtFBud0g6a3iuXtZ0oU55TZO0n9L2iBpvaR/La7P9dxl5NUR5y0lbR/jKj4C8CvgMxRug64EZkfEq21NpAJJm4GuiMh9sqKkvwPepzCr+G+K624C3o6IG4tF/4iIuKZDcrsBeD8ibml3PmW5jQHGRMRqSSOBVcAlwNfI8dxl5HUZHXDeUpJHj6uWRwAMiIjngLfLVs8EFhd/X0zhL37bVcitI0TEtohYXfz9PWADhZnYuZ67jLxskPIoXJWm93eKAJ6UtErSnLyTGcBRe+e4FH8emXM+5eYWn/BflNdlbClJE4CpwIt00Lkryws67Lx1ujwK16Cn97fZGRFxKoUn1q8uXhJZbe4AjgVOofCc2ffzTEbSYcDDwLyIeDfPXEoNkFdHnbcU5FG4Bj29v50iYmvx5w7gEQqXtp1k+94n54s/d+Scz0ciYntEfBgRe4A7yfHcSRpBoTg8EBFLi6tzP3cD5dVJ5y0VeRSuWh4ByIWkjxUHTZH0MeCzwLrsvdpuGXB58ffLgUdzzGUfZa8i+QI5nTsVXhR1N7AhIuaXhHI9d5Xy6pTzlpJcZs4Xb/f+gP9/BOA/2p7EACT9NYVeFhQeh/pRnrlJehA4m8JrT7YD1wM/AX4MjAe2AF+MiLYPklfI7WwKlzsBbAa+Ue2Zsxbl9mng58BaYE9x9XUUxpNyO3cZec2mA85bSvzIj5klxzPnzSw5LlxmlhwXLjNLjguXmSXHhcvMkuPCZWbJceEys+T8H98F994PeeQLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(D3[2,:].reshape(28,28),cmap=\"binary\")\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7fe85014f610>"
      ]
     },
     "execution_count": 91,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAS4AAAD8CAYAAADJwUnTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAVeklEQVR4nO3df6wV5Z3H8fcH1LhbiVKpQBBW1rKurF2kIrax6w+aumi7Ys2WhRp/tV3qRnQb20brH2qyMXHtLzVrJbdK1MSKJljFikVr1x/NKgEs5adQQlm8QmXVRm2N2qvf/ePMbQ/nnjPn3HvOuTPP9fNKTu6Z+c488/WkfPvMM8/MKCIwM0vJqKITMDMbLBcuM0uOC5eZJceFy8yS48JlZslx4TKz5LhwmVnXSFoqaZ+kTQ3iknSLpB2SNkj6eCvtunCZWTfdCczNiZ8JTMs+i4DbWmnUhcvMuiYingZey9lkHnB3VDwHHCZpYrN2D+hUgq2Q5Gn6Zl0WEWpn/0H+O90MvF213BMRPYPYfxLwYtVyb7Zub95ObRUuSXOBm4HRwO0RcUM77ZlZct6OiFlt7F+vyDYtnEM+VZQ0GriVyjnqdGChpOlDbc/MykNSS58O6AUmVy0fCexptlM7Y1yzgR0RsTMi3gWWUTlfNbPEjRo1qqVPB6wALsiuLn4CeD0ick8Tob1TxXrnpifVbiRpEZWrBWaWiA71ppB0L3AaME5SL3AtcCBARCwBVgJnATuAt4CLW2m3ncLV0rlpNlDXAx6cN0tBB08DiYiFTeIBXDrYdtspXEM6NzWz8utU4eqWdk5S1wDTJE2VdBCwgMr5qpklbhgH54dkyD2uiOiTtBhYRWU6xNKI2NyxzMysMGXvcbU1jysiVlIZXDOzEUJSp64Yds2wzpw3szSM6B6XmY1MLlxmlhwXLjNLjguXmSXFg/NmliT3uMwsOS5cZpYcFy4zS0rRt/O0woXLzAZw4TKz5Piqopklxz0uM0uKx7jMLEkuXGaWHBcuM0uOB+fNLCke4zKzJLlwmVlyXLjMLDkuXGaWHBcuM0uKHyRoZklyj8vMkuPCZWbJceEys6R4AqqZJcmFy8yS46uKZpYc97iGySGHHJIbnz9/fm78nXfeyY2fcMIJDWNjxozJ3fe8887LjT/55JO58Zdeeik33k2//e1vc+MPPfRQbnzt2rWdTMeGwYgf45K0C3gTeA/oi4hZnUjKzIpV9sLViRPZ0yPieBcts5Gjv9fV7NNiW3MlbZO0Q9JVdeKHSnpY0q8kbZZ0cbM2R8ypopl1TqcG5yWNBm4FPgP0AmskrYiILVWbXQpsiYh/kvQRYJukeyLi3Yb5tZlXAI9JWidpUYPEF0laK8mDHWYJaLW31WKPazawIyJ2ZoVoGTCvZpsAxqjS4CHAa0BfXqPt9rhOjog9ko4AHpf0QkQ8vV9GET1AD4CkaPN4ZjYMBjHGNa6mU9KT/ZvvNwl4sWq5Fzippo3/AlYAe4AxwL9ExPt5B22rcEXEnuzvPkk/plJdn87fy8zKbhCF65Um49v1GqrtwPwjsB6YAxxNpRP0TES80ajRIZ8qSvqQpDH934EzgE1Dbc/MyqODp4q9wOSq5SOp9KyqXQw8EBU7gN8Af5vXaDs9rvHAj7PkDwB+FBE/baO9tlxzzTW58W984xvDlMngzZ07t+gUhuxb3/pWbnzLli0NY/fee2/uvs3iu3btyo3b0HVwOsQaYJqkqcBLwALgizXb7AY+DTwjaTxwDLAzr9EhF66I2AnMGOr+ZlZOnXyQYET0SVoMrAJGA0sjYrOkS7L4EuA/gDslbaRyanllRLyS166nQ5jZAJ2cgBoRK4GVNeuWVH3fQ2WoqWUuXGY2QNlnzrtwmdkALlxmlpQRf5O1mY1MLlzD5Nxzzy3s2K+++mpufMOGDcOUyUDbtm3LjR9zzDG58cMOOyw3PnPmzNz4cccd1zB2/fXX5+7b7HfzdIju8YMEzSw57nGZWVI8xmVmSXLhMrPkuHCZWXI8OG9mSfEYl5klyYVrmJxxRv49ms3mKzWb75Tnrbfeyo03e8VXmTV79drGjRtz41OmTBnysc8+++zc+COPPDLkti2fC5eZJceFy8yS48JlZknp5IMEu8WFy8wGcI/LzJLjwmVmyXHhMrOkeALqMNq5M/dtRk3jVt/nPve53Hg787Teeeed3Pjtt98+5LatPS5cZpYcX1U0s6T4VNHMkuTCZWbJceEys+S4cJlZUnzLj5klyT0uK9RBBx2UG7/lllty4xdccEEn09nPJz/5ydz4+vXru3Zsy1f2wtW0PyhpqaR9kjZVrfuwpMcl/Tr7O7a7aZrZcOqfEtHsU5RWTmTvBObWrLsKeCIipgFPZMtmNkIkX7gi4mngtZrV84C7su93Aed0Ni0zK0qrRavIwjXUMa7xEbEXICL2Sjqi0YaSFgGLhngcMyvAB/6qYkT0AD0AkqLbxzOz9iU/ON/Ay5ImAmR/93UuJTMrWidPFSXNlbRN0g5JdcfDJZ0mab2kzZKeatbmUAvXCuDC7PuFwENDbMfMSqaTY1ySRgO3AmcC04GFkqbXbHMY8APg7Ij4O+ALzdpteqoo6V7gNGCcpF7gWuAG4H5JXwZ2t3Ig657TTz+9Yez888/P3feiiy5q69h//OMfc+OXX355w9gLL7zQ1rGtezp4qjgb2BERO7N2l1G5uLelapsvAg9ExG6AiGh6Bte0cEXEwgahTzfb18zSNIjB+XGS1lYt92Tj2v0mAS9WLfcCJ9W08TfAgZKeBMYAN0fE3XkH9cx5MxtgED2uVyJiVl5TddbVXqQ7ADiBSmfoL4BnJT0XEdsbNerCZWb76fAcrV5gctXykcCeOtu8EhF/AP4g6WlgBtCwcJV7soaZFaKDVxXXANMkTZV0ELCAysW9ag8B/yDpAEl/SeVUcmteo+5xmdkAnepxRUSfpMXAKmA0sDQiNku6JIsviYitkn4KbADeB26PiE2NW3XhMrM6OjkBNSJWAitr1i2pWf428O1W23ThSsDs2bNz44899ljD2OjRozudzn4i8m+G2L17d8PYe++91+l0rAP8IEEzS1LZb/lx4TKzAVy4zCw5LlxmlhwXLjNLStEPCWyFC5eZDeCrimaWHPe4rG3z58/PjXd7rlaeZq8/e+SRRxrG1q5d2zAG8PDDD+fGH3zwwdz4xo0bc+PWmAuXmSXFY1xmliQXLjNLjgfnzSw57nGZWVI8xmVmSXLhMrPkuHBZ25YvX54bP/bYYxvGTjzxxNx9x40bN6ScOmHWrLx3LDSPX3vttbnxm266qWHsxhtvzN13374P9juOXbjMLCl+kKCZJck9LjNLjguXmSXHhcvMkuPCZWZJ8QRUM0uSrypa25599tnc+Gc/+9mGsSlTpuTue/jhh+fGJ0yYkBs/99xzc+Nf+tKXGsba/X/1Zv+4rrjiioaxE044IXffOXPm5MabvU8ydWXvcTUtq5KWStonaVPVuuskvSRpffY5q7tpmtlw6j9dbPYpSiv9wTuBuXXWfz8ijs8+K+vEzSxBrRatIgtX01PFiHha0lHDkIuZlUTyp4o5FkvakJ1Kjm20kaRFktZKyn/AuJmVxqhRo1r6FJbfEPe7DTgaOB7YC3y30YYR0RMRsyIi/45ZMyuN5E8V64mIl/u/S/oh8JOOZWRmhSq6KLViSD0uSROrFj8PbGq0rZmlJ/kel6R7gdOAcZJ6gWuB0yQdDwSwC/hq91K0duzevbuteDOPPvpobvzJJ59sGLvsssty9509e/ZQUmrJqaeemhv/5je/mRtv9jyv1JW9x9XKVcWFdVbf0YVczKwkki9cZvbBksKDBMudnZkVopNjXJLmStomaYekq3K2O1HSe5L+uVmbLlxmNkCnCpek0cCtwJnAdGChpOkNtvtPYFUr+blwmdkAHexxzQZ2RMTOiHgXWAbMq7PdZcByoKW3lLhwmdkAgyhc4/rvjMk+i2qamgS8WLXcm62rPtYkKtOqlrSanwfnravuueeehrH77rsvd9+f/exnufFTTjllSDm14qMf/WjX2i67Qc7ReqXJXTH1Gqp9JtBNwJUR8V6rx3XhMrMBOnhVsReYXLV8JLCnZptZwLL+HhxwlqS+iHiwUaMuXGY2QAfnca0BpkmaCrwELAC+WL1BREytOu6dwE/yiha4cJlZHZ0qXBHRJ2kxlauFo4GlEbFZ0iVZvOVxrWouXGa2n07fh5g9aHRlzbq6BSsiLmqlTRcuMxvAt/yYWXLKfsuPC5eZ7afoR9a0woXLCtPX15cbX7duXW68m/O4tm/f3rW2U+DCZWbJceEys+S4cJlZcly4zCwpKTxI0IXLzAZwj8vMkuPCZWbJceEyJkyYkBtftKj22Wv7e+GFF3Lj999//6BzKoNm4ygzZszo2rGbzSFbvXp1145ddp6AamZJ8uC8mSXHPS4zS44Ll5klxWNcZpYkFy4zS44Ll5klJ/mripImA3cDE4D3gZ6IuFnSh4H7gKOAXcD8iPhd91Itr2bztFatyn+r+Mc+9rHc+NixYwedU1mMHz++YeyKK67I3XfOnDmdTudPtm7dmht/5plnunbsskthjKuVstoHfD0ijgU+AVwqaTpwFfBEREwDnsiWzWwEGMSbrAvRtHBFxN6IeD77/iawlcortOcBd2Wb3QWc06UczWyYlb1wDWqMS9JRwExgNTA+IvZCpbhJOqLz6ZlZEcp+qthy4ZJ0CLAc+FpEvNHqf5ikRUD+zXhmViojonBJOpBK0bonIh7IVr8saWLW25oI7Ku3b0T0AD1ZO9GBnM2si1J4kGDT7FQpvXcAWyPie1WhFcCF2fcLgYc6n56ZFWEkjHGdDJwPbJS0Plt3NXADcL+kLwO7gS90JcME3HTTTbnxZtMdmpk6dWpuPO+xN2+//XZbxz744INz41deeWVuPG/Kw5gxY4aUU79m/3DefPPNhrHLL7+8rWOPdMmfKkbEL4BG/xWf7mw6ZlYGyRcuM/tgKfo0sBUuXGY2QNkH5124zGwA97jMLDkuXGaWFI9xmVmSXLg+AJ544onc+Pz589tq//nnn8+N//KXv2wYe/3119s69qGHHpobnzlzZlvttyNvnhbAOeec0zD21FNPdTibkaWThUvSXOBmYDRwe0TcUBM/D+ifEPh74N8i4ld5bbpwmdkAnbqqKGk0cCvwGaAXWCNpRURsqdrsN8CpEfE7SWdSuUXwpLx2XbjMbD8dHuOaDeyIiJ1Z28uoPBLrT4UrIv6navvngCObNerCZWYDDKJwjZO0tmq5J3uwQr9JwItVy73k96a+DDza7KAuXGY2wCAK1ysRMSuvqTrr6j4lRtLpVArXp5od1IXLzAbo4KliLzC5avlIYE+d4/09cDtwZkS82qzRcs/rN7NCdPCxNmuAaZKmSjoIWEDlkVjVx5oCPACcHxHbW2nUPS4z208nHyQYEX2SFgOrqEyHWBoRmyVdksWXANcAhwM/yIphX5PTTxeuTnj88cdz48uWLcuNL1iwoK3jFzmXqh19fX258WbPOVu+fHlufPXq1YNNyTKdnMcVESuBlTXrllR9/wrwlcG06cJlZgN45ryZJceFy8yS4puszSxJfpCgmSXHPS4zS44Ll5klxWNcHxC7du3KjV988cW58RUrVuTG58yZkxvfvr3xZOOzzz47d99m8t7Z2Iqf//znQ257/fr1bR3bhs6Fy8yS48JlZsnxVUUzS4rHuMwsSS5cZpYcFy4zS44Ll5klp+yFSxF1H//85w2kycDdwATgfSoPw79Z0nXAvwL/l216dfbcnby28g9mZm2LiLaqzowZM+LRR5u+rwKASZMmrWv20L9uaKXH1Qd8PSKelzQGWCep/8l534+I73QvPTMrQtl7XE0LV0TsBfZm39+UtJXKK4fMbIQqe+Ea1CwzSUcBM4H+Z+IulrRB0lJJYxvss0jS2pp3r5lZiXXwZRld0XLhknQIsBz4WkS8AdwGHA0cT6VH9t16+0VET0TMKuI82MwGr9WiVWThaumqoqQDqRSteyLiAYCIeLkq/kPgJ13J0MyGXdlv+WmanSpl9Q5ga0R8r2r9xKrNPg9s6nx6ZlaEkdDjOhk4H9goaX227mpgoaTjqbxOexfw1S7kZ2YFKPvgfCtXFX8B1PuvyJ2zZWZpKro31QrPnDezAVy4zCw5LlxmlhRJpb+q6MJlZgO4x2VmyXHhMrPkuHCZWXJcuMwsKZ7HZWZJ8lVFM0uOe1xmlpyyF65y9wfNbNh1+nlckuZK2iZph6Sr6sQl6ZYsvkHSx5u16cJlZgN0qnBJGg3cCpwJTKfyVJnpNZudCUzLPouoPKQ0lwuXmQ0watSolj4tmA3siIidEfEusAyYV7PNPODuqHgOOKzmeX8DDPcY1yvA/1Ytj8vWlVFZcytrXuDchqqTuf1Vuw2sW7dulaRxLW5+cM37JHoioqdqeRLwYtVyL3BSTRv1tplE9pKeeoa1cEXER6qXJa0t67Poy5pbWfMC5zZUZcstIuZ2sLl655O171dtZZv9+FTRzLqpF5hctXwksGcI2+zHhcvMumkNME3SVEkHAQuAFTXbrAAuyK4ufgJ4PXufa0NFz+Pqab5JYcqaW1nzAuc2VGXOrS0R0SdpMbAKGA0sjYjNki7J4kuoPAb+LGAH8BZwcbN2FZF7KmlmVjo+VTSz5LhwmVlyCilczW4BKJKkXZI2SlpfMz+liFyWStonaVPVug9LelzSr7O/Y0uU23WSXsp+u/WSzioot8mS/lvSVkmbJf17tr7Q3y4nr1L8bikZ9jGu7BaA7cBnqFwGXQMsjIgtw5pIA5J2AbMiovDJipJOAX5PZVbxcdm6G4HXIuKGrOiPjYgrS5LbdcDvI+I7w51PTW4TgYkR8bykMcA64BzgIgr87XLymk8JfreUFNHjauUWAAMi4mngtZrV84C7su93Ufkf/rBrkFspRMTeiHg++/4msJXKTOxCf7ucvGyQiihcjab3l0UAj0laJ2lR0cnUMb5/jkv294iC86m1OLvDf2lRp7HVJB0FzARWU6LfriYvKNnvVnZFFK5BT+8fZidHxMep3LF+aXZKZK25DTgaOJ7KfWbfLTIZSYcAy4GvRcQbReZSrU5epfrdUlBE4Rr09P7hFBF7sr/7gB9TObUtk5f775zP/u4rOJ8/iYiXI+K9iHgf+CEF/naSDqRSHO6JiAey1YX/dvXyKtPvlooiClcrtwAUQtKHskFTJH0IOAPYlL/XsFsBXJh9vxB4qMBc9lPzKJLPU9Bvp8qDou4AtkbE96pChf52jfIqy++WkkJmzmeXe2/iz7cAXD/sSdQh6a+p9LKgcjvUj4rMTdK9wGlUHnvyMnAt8CBwPzAF2A18ISKGfZC8QW6nUTndCWAX8NVm95x1KbdPAc8AG4H3s9VXUxlPKuy3y8lrISX43VLiW37MLDmeOW9myXHhMrPkuHCZWXJcuMwsOS5cZpYcFy4zS44Ll5kl5/8B9UH4qUr3FakAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(abs(D3[2,:]-1).reshape(28,28),cmap=\"binary\")\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fe84f83cb80>"
      ]
     },
     "execution_count": 123,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABGMUlEQVR4nO19e7hdVXXvb5xz8k5IAgkh5IQk0shDBcQQsGgJpQpoLd5evKIWH/dyqV8V23q/Vlu/295+vUKtrfVZkGupcrVSL1JFRcEq+ICChDchEAIkEELI+52ck3POuH/svdeea6055hpz7bX3Xnvt+eMLZ6+5xhpzrNdYY/7mmHMSMyMgICAgoPcx0G0DAgICAgKKQXDoAQEBARVBcOgBAQEBFUFw6AEBAQEVQXDoAQEBARXBULcqnjdvHi9durRb1QcEBAT0JB544IHtzDzftq9rDn3p0qVYvXp1t6oPCAgI6EkQ0UZpX6BcAgICAiqC4NADAgICKoLg0AMCAgIqguDQAwICAiqC4NADAgICKoJMh05ENxDRViJ6XNhPRPR5IlpPRI8S0ZnFmxkQEBAQkAVNhP5VABc59l8MYHn935UArm3drICAgIAAX2TmoTPzz4loqUPkEgA3cm0e3nuJaA4RLWTml4oysl/AzPjaPRuw88AoAGDB7Kk4fs40PLZpD9599gmYN3NKJDsyNo6v3r0B5500HycfdxQA4Jlt+/HdhzdjkAjvPGsxjps9FQCw5+ARfONXG3HMjMmYYOC4o6bioed3YXjudGw/MIL3nL0Es6dNAgC8tOcQ/vX+F7D74BGsWDoXF73qOPzz3RuwaddBLJo7De//9WWYPGSPA/aPjOHG/9iAl/ccxonHzsTl5ywBEaXkjoxP4Kt3b8Cug6MYGZvAZWctxvIFs6w6H3p+F+58civ2HDqCuTMm4x0rFmPRnGlWnf9893OYOmkQB0bGMTY+gctfvwRzpk9Oye49fARfv3cj3vqahRgZm8D3H9mM4aOn47+sWJySfXnvYdz0qxcwffIgPnDuUvxozRas27IPv3nKApyxeE5K/p5ntuPeZ3YAAM476Vjc99wOLJozDZecsch6fnsO1Ww5cf5M7B8ZwzEzJuP8k4+1ygLAN+7biEOj4yAinD48GyuWHi3K3r1+O55+eR9Gxiaw6qRjcdJxzWv8y6e341fP7cCZS+ZixdKj8bV7NmDr3sOYM30yLn3dMBYfPR1A7Z5+7Z4NGDkyjvNOOhaHRsfxwMZdOPaoKZgyNIDfPXM40rlt3wj+770bsf/wGMYnJjB72iTsPTyGo6YO4TXDc/CmUxdEsht3HMAPHnsJs6ZOwjNb9+OUhbNwzIwpWPvSXvzeOUswd0bzvj3+4h7c8cTLeNtpC7F8wSzc++wO3LN+O04bnoPfMnQCwM4Do/jGvRtxZHwCB0fHMThAOGXhUXj7a+3Xv2ooYmDRIgAvGNub6mUph05EV6IWxeOEE04ooOpqYdOuQ/hf33siVjZAwAQDR8+YjN87Z0lU/vDzu3HND5/EDx/fgu986FwAwNfu2YAb/6M25mDqpAH8/nknAgB++tTL+NsfPRUdO2f6JOw+eCTaXjx3Ot52+vEAgFsefBGf/fenAQBfvWcDvvfhN+CTt62NZF+3ZC5et8TuRH759PZYPW86dQEWzk473yc2743pfPzFPfjX33+9VefnfvI07npqm3E9CB+5YHlK7rEX9+Dq256MlS2cMw2Xvm44JfuLdTU7n9qyD1OGBvCt1ZsAAG877XhMmzwYk7314c34h39fBwB4/YnH4GM3P4oDo+N44qW9+Mr7zkrp/tQPn8Qjm/YAAG68d2N0nSWHftdTW/Hp25+KlW34m7daZXceGMUn/q3JfM6cMoTH/+pCqywAvOcr90W/b1+zBbf8wbnR9tW3rcUTL+3FifNn4GMXnRyzgQF89E2vBFD7KDT2PfTCbry46xCe3X4gkr3kjEUYHKh9tL/94CZ8/idPW205fvbUmEP/uzvW4XuPbLbKLpg9NfZxvfZnz+AHj76EbftGcM3vvgZ/+6Mn8eDzu7FozrSUQ//xE1vw9z9el9LZLw69iE7RdAhWeybShczXM/MKZl4xf7515GpfY6K+2Mg/vPN0/PXbX10vi+9rYLy+vXXv4WbZBGPWlKHY/lp5vJ4jY/ECU/fEhL0eSZekp2GPDUmd2/ePiDqTOiSdSbulMrP+fYfHYueTtD9p6/gER9uuc7vg5GOxbN6M1HW22u2xwEyyzv0jY+pjt+6LX+NGvROctsG8bo3fs6YOYYIZRybkc3Kdb/Ke7zl0RJBM37fGduPvuPBOAO7nsx9QhEPfBMBsqw4DsH96A5xwvdupfRZZBqLPqymfXJUqpYpd+9i57bJREm1lkSzpUFs5C9LmOZgyVh2Ja9PYFu1wXMss+UxZlUalrsZ5MKfvm+WaEOy2alc88ztP+7GRXfUCqz0FXqNeRBEO/VYA761nu5wDYE/gz1sDgazNHt2x2QirDrYXRLID7FXY+kICyodMDp2IvglgFYB5RLQJwF8CmAQAzHwdgNsAvAXAegAHAXygXcZWHa73PyvKrsnYX7x0cJ/UZURkqWjNw0Z1dKT3dOmWiRR168pcptiizfi1aW5pWh+a6+Hl9Av8QDRsY4vaeOuu9pdIioi19XnYJrQYkq0j2/Wt0kc0DzRZLu/K2M8APlSYRQG1CC9nQKQ5rt8f+vajFqJX6TqH+Lw3EEaKlghOfjola5eiaL9M5ro59CRnrqnXT7YV3ljmrm3RdXb9poyGQwfLshp9LvkiZTN1GXS0q1XWuP5E1FJE3NI9T0bmcSo9caxUf4W+rg4Eh15S5GXRVVxnfzzbXUODQ6/SdQ4Rem8gOPQSwc2hJ2Xt0RJZ5F2ceUo2pdZ9rCzpyDIRNVhktVG/4timbJwXz5I39zf5XInLNzNENBx6C/0JLcDkoTXPg8yht7/fJBWZG/x/5sHu4sohOPQSgogCh97DINQpikpd5xCj9wKCQy8RnPy0QrZWRil5n8yVzHqcHLqbf88q18jKueWWMoVsTMJ6TROMshcXXIyMUbuHdIYuI5dblYdOfq2gvHJmncljU9kuPhy6vvqeRnDoJQQhfzyki9D75fHuDpp56NW5ziE+7w0Eh14qOPhpTR56LMtFlk1HQDKJ7hGgW+wRyn1444ysG0nOKSucbpaOGofuoTtW7kH+C2gLh84WE2LnXNsggt/Nl1UqhIWnLhWZ+9zz6nxcXQgOvYQIeei9DQKJnYi9ivxjlwM6ieDQSwQ/btkenVDEoXtkbwi/bcf65aELHLbbHLdOUVBfj3S+WZwsw+CexQwevywLH1680O+DJQp31VPj0PURcVou/3mmR4jKdRfxzPUygkMvEaIOKI/ZXJIPaq7IXmZcuo6WzCn4ZKz0hAU1Dr1aEa1qjqC2WxGQheDQewSaSJWhy0P3qzcZvTl4/lTOuqQztzleeeiyDrb/9qnPt1y0RdhhlS3OZTajXFvUnW7dSWmY2meriFZZMzOHY+WaesoWqLQLwaGXCOYgjrwBXp5Z8eJpaiV78jvVEafRx7q0xVb6QAICWkFw6D2ClKO18of23S1FxI46XPW7ZFtqMXhy1HYdwm9btOrNiQvlbehPaAVmlOu6b825XPxaMHnlTNuS2+ln0SPLpWyBSpsQHHqJEL086GzeryNrsetojS4q0BDonW/V+HMgtDh6BcGh9wg0IyZjZRnccP569bKidEsGScWth+heHLoYuftF4t36gJqZIq6+j4gGBAktGGV9rWS5JGxJZr1o7ClboNIuBIdeIpgcupYLT0rl4W8lCqIMaK0DtfgsFxU87l+vQHM6ZTrjil1+NYJD7xGouGwp2mzBsaVHarqyXBLbnQvQPTl0eydwdh66WzaPfV556IVy6M36nRx6FKFLLZgOZLkkIvLknC42e/rUnweHXiY0H9AW1hQl/4fZyaF3OWRvh/Nvt8JO94F0Ar3W4ug1e4tCcOg9Ak0GiZTZUqRTdHPoMhcbl2vBnoKzRTLz9YW+CLE+b85dUqTXkQdmlOvqn2n8IlFWXaGHbfZtc4ZIc9uGpDsPHHpAx2GmiOWeywXkHZ248tC7/SK0woMXnuWijdBbmM++tOixhbAqd/2VCA69R6AZhSmtltNaHro+JNNz6MU76bxZFLHDPDh03wSeYka4eggrdbHFhngLpslJs8UGdR66l232fpsmhy7rNLNy4vWX6XPTPgSHXiKYHVCtzLbYCoee8t/djtC7dKxVnw+HXrEIUdPq6/azEkPFrr8WwaH3CNKRUfrtETn0lrJcktv6EL0I3jjLnqzyrPptfLEsaz8uLu95ci3kZ7cCM6dbk7lE5kE5jCokhTTJnVvNabYo4vW3Xn0vIDj0EoJIP9ti6lgUm4fe7YUBypWHrtNXydkWVRx6ebxm1VpIWgSH3qPwiyaLq8c5l4uyQ7UlewoI0aX+BXsAao/gi8tD16M9H9e0G2bLg0SgOoee7MvR2dRKX4GZM2/q8uPQ+wPBoZcIxXDoOaJ7szM1d15ae1CqF7GvOfRsmTLRGlW7/loEh94jSPf8W2SE/UVmuTgjdGWHamucvsBde7h+H17ch29Pymvs88tDLw6xnG7HNzzZmvFpseWRS9Zpbkc6klx6TLaGNIdeoq9NGxEceolg5qHnBUX/86nX+F2yLJeWpi0o2HavPPRiq+46VFkuHbBDizBSNKDUSDsTWzRpRpDZ0WSeel26UtGbGE0XZ09WeVb92Ry6JOvXUigkf75Aj2kGu65WWHwuF8u8L+r68p+nmZFj1pnVjxSvvz8QHHqJEL08rczW1+pcLsl9+awoDK116BZrvVZbrZVUrQhRdTpdb841Ua2rr0dw6D0CDXfJwkZrw+eT0ZusS03XtOG99zpFYSRVVsSXtf6oy45iLkVxF86cDyUddafPk6jBoeufB1t9KlnJ1ohLT3DqsWNb76voZagcOhFdRERPEdF6Ivq4Zf9sIvoeET1CRGuI6APFm1p9NDt0WphtETny0GNZLnabuoXW8tCLswPwuBY5Wkllh+aJ7PazEkPVboASmQ6diAYBfAnAxQBOBfAuIjo1IfYhAE8w8+kAVgH4eyKaXLCtfQ2Vo3VQJ22tN9qnjNZasK6Q2RZj+ty6fXLW3eWtR47Fcuj1KBfu1p9JA9ZGlSb1aOvzsE1oMaQ5dFs/UgEG9DA0EfpKAOuZ+VlmHgVwE4BLEjIMYBbViN+ZAHYCGCvU0j5A9NK3Mttijjz0OAFRLha9tQWli4ZOY23Gy8Ir7yp6jELv1wBd5dAXAXjB2N5ULzPxRQCnANgM4DEAf8jME0lFRHQlEa0motXbtm3LaXJ/QpMPHo8mZRrFq14tL+4h21pevKQzXxaFFK3bBFydx5l1epa3Kpupy8gccUXd0TNFlBnNuyv0sE141jm53YLOqkLj0G0fu+TVuRDAwwCOB3AGgC8S0VGpg5ivZ+YVzLxi/vz5nqZWH00OPf9cIPk4dPtv06ZuoRc59DwzXpYdqgi9609LEyEPXcYmAIuN7WHUInETHwBwC9ewHsBzAE4uxsQAQOYVJRlphGOr9brcWpqsKT7joIg843jrxX2kGMGLnLjfOXc9D51Zl4dOyMyIcdeX/zxTkbnB/6ePLf6Z6yVoHPr9AJYT0bJ6R+dlAG5NyDwP4AIAIKIFAE4C8GyRhvYDmi9Pfg4213zozn1d5tBbomeKtV2rr9ZKqlaE2HMcerUuvxpDWQLMPEZEHwZwO4BBADcw8xoi+mB9/3UA/hrAV4noMdTu/ceYeXsb7e47aLhLL25YW68Hh57c2Y5Re0XPiZJ1naS+CG9OvABnV+iaogYv7WqFRTRg5CHzcXItOXuO/4h0ZbwDVhUVR6ZDBwBmvg3AbYmy64zfmwG8uVjT+hG1x866anzSWVqP5mhNUZ8OPBcFMaGoV9qn7cD00unhMTWdslnUlO9HktneSirkQ1QoVdX06M5rzM1nMrsD1VG/1lDYno+4I09SL7Z6wuRcAQEBTvSLU7ChXymMXkNw6CWC2QGVfIHSUZQ9mmxEh3kXic6KiP3SFoWoVG+OOppvG+Ui6pDOja156DIV0KVOUYNycVFlzYiXkO4+lVs76fpa6BQ1bXXImWXJPox++RQHhx7gRLc7RcuEPg7Qq5eHWVEEh14imNFQ8g3SOJNIhLK5YUl3VkTu0tWOaFp7rLVPQZFamDXhlnQdszl0XYTYrcExUb8i28YGp8+zyaHb+e0s+7xaZQ5bY3+txzY5/1h5n3yMg0MPcKJfXgQN+vlaBA69NxAceong5tBlntMss62Wk+WIXFGWPsuktYwUUVLKzlDU5euAC5mcC/V7l+LQi+hP8BDO1NWMcl2tsoZctEi0yyaP/hW3bbKtsb9CPxKQ/f5UFcGhBzjRz5kdSfTzlQgRem8gOPQSoRkNZXOA9oiDo1GKuSerStXj3nbZ6BNNizoVLRObnFy7fL6ZHHoG396Ut7SSFLZ0Es0o1933EfXriEP/7b9bs81uTzoP3XZsAxlpYhVFcOgBbvTJi6BBP1+K6i17XU0Eh14imFkqWXm0Modej6Zi5RlZLq56UhGyPpOhiGHY2jx4Hw5d6jPIjJaFaD0pZO0D8WhZiNUXyqE363fdN5OTrvHt8vOQ+ZwpTyArDz05ctR5cHRMfyA49AAnAoXeRL90rAX0LoJDLxGaOb/ZmSrSqMZopKgiI8O2P4tDd+uRuVipvkydqW19BKaZPyVuiy1rwiOCr8vY5rMvZlKx4j4o5kLLrv6ZJodeC9E10bxcp9I2Ydu1OLQpa+vA7ZfAJDj0ACf65UXQoJ8vRWDQewPBoZcI0Sg3TR667Xg2Zlu0RFlZ9dqktVkmVr2F8MbKLBdHTnK6fruM1OrRyjbkbfeviK9BoRy68Tel1tJXEPXLZETHmjqzjXM3R139Q40WUrru/vgcB4ce4ESI0JvoF6dgQ4jQewPBoZcJEYeeHeHZRzW2n0N3z+WSbaPGHnf9+Y+11Z/VkpFz1qVzE2ZbFG3xaK2oJT2Usa7vozG2oTUOXZnlImwn89BtdTLYulpUvwQmwaEHuNEvb4IC/XwlQoTeGwgOvUQwMwqyZuuToslGdB9nxX04a3c97gwDZZaL2pq0sE8eupyTbM9csXPodoFMDl05W6Zfn0SBWS7RX1vfQ/o8xTx0IaffVWembcI9Z8NqSafMofcHgkMPcCIE6E3096UIMXovIDj0EsGMhjJn65MyMuoH5+U2s0aGOiN0dTSdnzd2rRSUdazNrqyWjMyhy7ptfSA++fMSivyguPPQzd+1LarPsd8ah660Tci0snPo6VahPQ+9Pz7HwaFXEEXOjNcfr4EO/eIUbOil+Lyf550JDr1EMFdbyTPbYi3DoimR/iXU68oayG4YiPuKSPPTtExqcrqytFCGDum3mMHjl2VRRF5/HphstGusQZxDd8+2mF2nTlrk0CMuXa6/iOvcywgOvYIoMj7pk/dAhX5xCjb0VNTbQ6YWjeDQSwQzGsoz22Lj2NT+DE/kjnjcnLrTpiKiUp3K7AyVmKy9zyArynddp2S5dtV5v1ZMcV+UZtSrm+O8yaHLz0P2bItK2yRbIy5dvmf9PPgLCA69kgiry7QHfc2h99Az1UOmFo7g0EuEprtQzLZoO56bD7MmI8OmuzUOPZ1xoJHT2lbb1meLaPjU7NGfkhJ7Mbg5Wjdep9RasJXpZfMiinbhbpVFrcZ65pQmmm/ZNqEOO4eev4+liggOvYIoku/s9yasiX5xClb0UNjbS62JohEceolgzmyXOduiwBk3uFuflXhco/20vLhNtpDMDm3Ub1GaxXOnfmdy6OZ1knVb53LxCPR9ZPPCzBxR5aELNnQkD72RM6+os5i+it5FcOgVRMhDbw/6OUK3pWKWFT2VkVMwgkMvEaKMAmTPBSKNarRz6B7ZB5kcugf/LWaZqFXoo36PelwjDZ2ysXLp3NjKoasqsNSTIZobWh66yaHXzzllg9y6S0uqQ3RnuTvTqv39D2WGyqET0UVE9BQRrSeijwsyq4joYSJaQ0Q/K9bMAB8Um4feJ2+CAv3iFGzooQC9p2wtGkNZAkQ0COBLAN4EYBOA+4noVmZ+wpCZA+AfAVzEzM8T0bFtsrfaaERDlOZgU6IC3xvloWfIWqpN/bYd69KlXl3IbY67flHQViREawIXnpXLrufQLXnoXi2LhhZbeUFw8tDpeqh+iGZUqVhlzgA9lYfu0NmJ/ocyQxOhrwSwnpmfZeZRADcBuCQh824AtzDz8wDAzFuLNTPAB73Ed/YS+jtC751nqncsLR4ah74IwAvG9qZ6mYlXAphLRHcR0QNE9F6bIiK6kohWE9Hqbdu25bO4wnDP5ZKd7cHgiHv3yUM3hbNGhrp0qTMOWvCMch66NbwWdAi/hX6JLB3Jcmseukd/QgaFXAjiWVCuqLv5TLaah661X5pzPSp23LN2PHO9BI1D18wXPwTgdQDeCuBCAP+TiF6ZOoj5emZewcwr5s+f721s1ZHnmSviMfWhXDqNlpx/gXYAwMRENZ2C6hIrwt5uPysN9FJromhkcuioReSLje1hAJstMtuZ+QCAA0T0cwCnA1hXiJV9But86EkhW3aEQb365KG74OPg07JSNF2cPc26/HUk18jMipbjv+WWgq0PROZ2W28V5IGrBWdrwRAIbLHWL2NJJyzx4pYAPXDoCWgi9PsBLCeiZUQ0GcBlAG5NyHwXwBuJaIiIpgM4G8DaYk2tPop6YX0DFOdAjW5H6K0cW7Dx3b4W7YLmvHTPVDkuUP/G54oInZnHiOjDAG4HMAjgBmZeQ0QfrO+/jpnXEtGPADwKYALAV5j58XYaXmUQLGuKKvjpWG5EBjeshbR6jFVWyaG3I0T3UmmOPDT7D2TR+v7sVk9zLIBytkVrq8CjnyAnXK0NWzZPk0PXPw+uOn3kzBaV+dcuq++rqCI0lAuY+TYAtyXKrktsfxrAp4szrf9Q1DPnHaE70vG6/SK0RBcVbHu3r0W7UFSEXprr08chehgpWkLkncsFbM9Db+VL4ceh60L0IlsMUXkuDj3Jx1r4bOFj5+RqbffPoz9BPJdCOXTdRzzOobf2wdeO0E3TfkaLCu57Jqsty9emvQgOvUQoivPVzGVhOhzni9blsKs1518wShOCFgzFefVShN7HAXpw6GVF5nzoQjTZzEN3c8Pah17iM62yqQC9WD6zsWqOVWeiLpcDMrMmpNG1huLUfqdzY1h6QPw4dIfqwsCQz9vGr0crFrXyvVfKSrnl0X3Lumc2nSX52LQbwaGXCJ3k0M1c3fgLHEe3XwSzfp/Iq9GJp5XV7G/SD9WAbSK3bGkZnZz3x2VNP+ehB4deQqjWFLUcJ0UutmheH6HrQzLtxyDva09E6hxt10ttrk0Zm0PeopsNXc1olZwjVq19IB79CZ3KQ7fNnZ+sJ/qIUU0ypccry8VoNbr6YoSHvcmhy/esHWMfegmqLJeAzqCwPHSNjMChe/jvjsCs3jtCV77GIULXy7qQ51k5asoArjp7LpbMmRQRVTOnjGPt2uYwlr9442wcGT8Kk4cGsHbtWnzpLQuiDu1Nzz2NzcbDfNHwOM47dh5AhHFjZO/hbc9j7c7eil+nTp2K4eFhTJo0SX1McOglhGbFGzmaTO+3yVI0f54bPhx6sqKiI00iBy9vkZWQvDZNWVu0zFYO3cWJ2+azF4+w3kc5+i8SMg+d7n+R+i9yZ7nUNV919lyceeLxGJo+K4q6582cguPnTItkB7bsw8jYOKZPHsSvHTsL4y/uqfd/ME5aeBSGBpuO+sXdh7D74CgIhLGJiaj8lQtmYeqkQb2xXQYzY8eOHdi0aROWLVumPq63PlmVRzEvrIpDNCP02Ausc8qdQpxD18fGjUWNtbKa/WYKXxVgm8jNhSyxPI/KkjmTYs68CFTh7hARjjnmGBw+fNjruODQSwj7ijf2nv9YGRtZLhmyeTl0dx66e1vSqYYjyyVNostqzKwJ5qas1OpJzY/jzLYRxhF49Ce0m0Pn5pepboPjHtc3ah/I1gywqK23RrNG1XKs3N1K9CouNfJ85IJDLxG6xaGbT7vUIdU9NA3w5dCLkk0yMlWIAAH7NBEuZEboHXxWstpU/Yrg0EsITYQnrXLf5NDdGQVa2iDl33OM8MuS08LJXVtkJTTz0NmZj92QtXPo0rmxvQ9EtMXFYOvKtXqb+2p/dXnozeexdQ7d/UyKx9kMy5ZWFXcbn/3sZ3Hw4MHC9AWHXiIU9syp8tDt9WqdcqfQGoeus13PoXNsu9dhG4TmQuaC2mX1mm3A+Pi4c1uL4ND7AGRZs0jzzsWzCMzf6YPVHHrSwbtkFTx/lg4XnCNFkxS6k0Pn6JhYFovAoSf3u0esNvpAdPfPHvXK0b8Wmr4O6byTkXTjabS2YLzy0JXHdfi7cOONN+K0007D6aefjssvvxzvf//7cfPNN0f7Z86cCQC46667cP755+Pd7343XvOa16S2x8fH8Sd/8ic466yzcNppp+HLX/5ydNyqVatw6aWX4uSTT8Z73vMeMDM+//nPY/PmzTj//PNx/vnnF3IuIW2xROgsh26MFDVfYKVT7hTM+v3z0PWymv2c2O51JM8rE5kken5bAOD//OJZPLftACYNDmDyUDPWPDg6DmbGwABh2qRBHBgdi+qbPmUodj9GxyYwNtH4YDOWzZ+B//7GV4h1rlmzBp/85Cdx9913Y968edi5cyc++tGPivK/+tWv8Pjjj2PZsmW46667YtvXX389Zs+ejfvvvx8jIyM499xz8eY3vxkA8NBDD2HNmjU4/vjjce655+Luu+/GRz7yEXzmM5/BnXfeiXnz5uW/cAaCQy8h8s62WIsO08SonUPXwWfuDrVszi9X7dx00aszQ6DBoXNc1spnc3pUIhFl5qEnL7A2f14qc5X7ykbUkTTaMva7sQKT/Zxz56G34vxZ7k9x9BxYS3/605/i0ksvjRzq0Ucf7ax65cqVsbxwc/uOO+7Ao48+GkX3e/bswdNPP43Jkydj5cqVGB4eBgCcccYZ2LBhA97whjc468qD4NBLhKI4SO88dMeL1m1e1HSy3hG60nR9hM6x7V5H8ryykM2ht4ZGJH3MjMlYNHd6VP7klr0YHZvAtEmDWL5gFh7btAcDBIwz49TEwKJNuw5i76FaBG8OLBJtNj7YDQwNDWGifiwzY3R0NNo3Y8aMmKy5zcz4whe+gAsvvDAmc9ddd2HKlCnR9uDgIMbGxjJty4PAoZcQaQZd6WiZ1S9pezj07JIsHU5kcNdJWQkc/eV4Hrokm+SaXfnwAECWcaJeHLqgN2dGiahG5NDjv13UjF+rwaT2vA5Mw3HPfHDBBRfgW9/6Fnbs2AEA2LlzJ5YuXYoHHngAAPDd734XR44cUem68MILce2110by69atw4EDB5zHzJo1C/v27WvhDOIIEXqJ0DUO3bShGBMKQ2scujbLxW9/1SJ0LbKez25nREXw6EB51atehU984hM477zzMDg4iNe+9rX41Kc+hUsuuQQrV67EBRdckIrKJVxxxRXYsGEDzjzzTDAz5s+fj+985zvOY6688kpcfPHFWLhwIe68806d0Q4Eh15CqGZbVGZk1H6nhbWD0NIrwjgiPyWHnj8P3c1dJ2UlNNembHLEol2c5tglLr+Z1mi5f2LLQriRWlkBbg699lfqO0h+4KM8dMc5exvl6ovJLDeaVQXgfe97H973vvfFyu69997o9zXXXAMAWLVqFVatWhWVJ7cHBgZw9dVX4+qrr47pSsp98YtfjH5fddVVuOqqqwo4i7oNhWkKaBlFxTeeFHoGh95dxPLQPd5hV2qhTVaznxPbvY6iF4no9rPSQEVuTy4Eh15KZHOwUsYBGb9dstp5IvwCMDnSi+vM9+qT49hkqevsmiNFkxyxrReAU9dU6nA189RT90/MzrGUechKcOehxzt3s/PQa/+1zqGnbciDPPF5WT427UZw6CVCHkdne7BVa4qa9TrawqXKcvEIjV30jE1Ws99c4KIKkBa4EJHJofvbUOuYzvOMVeMeuJDnugSHXkLkz0Nn60tqzUNXc+ju7YQBKtm8nwhXX1cqD92hh40fsWhPiLqtEbpDr2Y+exeK6HtwOWnzPLL0MmoXQObQPWyyUHsbdx/B2MG9uqwcVSU+wuVFYz70qVOneh0XOkUrCJ2zNkeKwvob6P77YdZfJg69LAkdraAMHPoX7tuFqwAsmbM9alnunzKIfdMnRzJb9hzG2ARj0iBhfOdUvLzrEAYHgPEJYGDPVAwONG/groOjOHyklkNurljEu6bERp/2AhorFvkgOPQSopaHngrRE5v2aMnOoadlc0forsgvtd06FxyHftEKV4ze5NAbzf1Gq0aSTWbBEJjTg1aaWTA+ef423l6QVeqs6dVICSNFEwxcY1yErvfCYZMlD33vyAQ++fMdMbnLzlqMv/nPp0Tb//Wan+ClPYdx0oJZuP2PX4uLP/4DzJs5Bdv3j+C+P78AC45qRrF/evMj+Pm6nRibYGzfPxKVf/dD5+KUxXPUtvYqeuuTVXF0MuqTHI7evXQILWS5aG3XRugNfRWh0J1L79mQOVK0gw+w6x5U5f7kQXDoJURt7ox4WSr6lfheSzO6Uxx6OmddkNNVnYKbQ0/LSojNtgh7q8aU9efQffLQLWVSJk9BDjOLQ09G0s08dFmXT72139kcf3Lb7ECVM3QEnXozexrBoZcInZzLhWIcOlt/17YLMSk38nLoQBs4dI5v9zpc0wbbUKY8dGeE3jkzSofg0EsI2wOpic7M2CU7D11ni890utoO1fx56PKiFUk7dSsWoc4Ry+l7bNkvLUBdXB66YLdQrtWbtEU672QkbVv3M5dNyuPSz5zRoopaF4Ltks5uRyYdQnDoJUJRz5zGV8fz0O2/gRJE6LE8dN9jdXIhQi9IYSf7gJwcekVuUA4Eh15C2PPQE9u27Ahma0eXnUNXjhRN1eEjW+wbrhn92ZTNPj9mVnDo5kjSph327JRGp2k6RJe5XXurwG6vvVyrN6kneV62ehr9MmInuheH7n4mpX2xkb31ssCh26Fy6ER0ERE9RUTriejjDrmziGiciC4tzsT+QWEReguRrNdAog4gzqH7nZh6tkXtSFGOb/c6nJOS5UAnRxVX5R4UjUyHTkSDAL4E4GIApwJ4FxGdKsh9CsDtRRvZb2jMnWFC62jtkUtaOD+Hrs9OUMvlzJRwqXBz6BwdE19+z11f1nzepmzq/om2WMoKYNG9WlKOexy1UITrmXedUzeHbt9mTk8TkHUutrqrDE2EvhLAemZ+lplHAdwE4BKL3FUAvg1ga4H29RWKeuZUsy0aMk5H1e0I3XSSXebQG9eiKsFh8rxaRUfHUYQ8dCs0Dn0RgBeM7U31sghEtAjAfwJwnUsREV1JRKuJaPW2bdt8be0bqDh0y3FeeejK5C4/Dl3OlnDKqSyR5xOx1aWZyyXKQ486B61X1T59rsNo6/0Ts3O0hb4cumOfMaLVJpvm0NMtRk09bgP1LT0rhy7dM1Ftf4ToGoduzaJLbH8WwMeYedyliJmvZ+YVzLxi/vz5ShP7B0V1IqpmWzQj9NgLnHTK3X0RTAeu/Qg1j9UhS2+UIhcNaqlGCJg8r1bR2ZHO8j3o5whdM5fLJgCLje1hAJsTMisA3FTvqJgH4C1ENMbM3ynCyH5EVoQnj2pMv6T2CF0Hn3fUNeLQKaf0BK5JsdJ56Jq5XORWjSlrnZzLodcWz/px6IKsUG7X64iA63/FCD32m+sjX6V6fGyy15G2T7CIjWsspFz6XOcqQuPQ7wewnIiWAXgRwGUA3m0KMPOyxm8i+iqA7wdn7o/OcujGSFGHXLdfhHhHY/5jXcjSaxv6XwW40jXzoKNZLs59VblD/sh06Mw8RkQfRi17ZRDADcy8hog+WN/v5M0D/FHL+3VnSYijGm0cui3LRWlLml9tPdNC0x9gA5E822KuHG3mWgRKA6IhDCO9L2aH5fpHeei2FpZsTaqk7Rx67W/yvGwHMwNw5aH7ZN6Yne/uED2+GePQG9dY4PR9+ioqCNX0ucx8G4DbEmVWR87M72/drD5FQU+dylmLHHpcrPscehPeEbryguojdI5t9zqS59UqOsuhO/ZV5QblQBgpWkLYZluUopakiI1bLJJDd4/wy+b5beVqR+DIckmJOvPQ638bv6NWjS1aZmO/aYesl2Br9uv6EyQ7XOVavSlTRA490R8Bx/XMzaFnc/zJ7VoeumEULM+Sou4qIzj0EiEPB2k9QjPbYkzGaAqnvhzeJhWLFjh0re3aCB2c2O5xWGaJaAmdfFTcHHr/Ijj0EoKQfig1uduNmfGS+22y+gg9u15pn5ytkdSpjbpdHHoiqnRF6NExjX4HeQh8nENv8reu7BQfDt2qpwAOPd2iS3+0pUWik2msReWhx55Jr5Yep47X2O7SWVUEh14i5Hnm8tIpch566zYVCfOFHfAMvTyYHCcGEs17XzvKiuR5tYrycOgWyrJPEBx6CaGK8KxZFqhlJFAiIrM6feVIUQ8Hr+1QzfvRIJCaI3WdX5NDb0Sqsh01Cj0ewdfsEPj2+v6skb7NcuE+WmVt9un4dlt0nDwvqyzsz2NSlwbq2RaF7UaLCnDnoWtGQlYVwaGXCHkeury5v+oIvcuvQiwP3TdC1w5Yyhopmug0rUr0lzyv1tG5Z6Wfc81dCA69hCALi67qzWcd/+4DHwefd46WVrhhqdjtdDmq15xD3hoBm78dHz5TtqbPff+cesSWja1VoNPrGqXpalmx8Ew1delvXla/jmxPs65Uqyp1LFuN7TZ12CkEh14iFMWhaxAfKZruMGtudxe2jrA8x7YC28CiKkAcWJQT/eI0y4zg0EsI+2yL2b35Nb6T0lkYOTtObfX4cOhiNJ3zzSfIEWGaQ5cRRXzcvGaSXbXMofhxNTvcNrQy22IRI0LTt8KMuhtcv72++KPD9RWYpJGZalNT2TOiXOo4s0VVK5Pumcyh98fXJjj0EqGwPHQFtBx6t2P0jqwpmrU/kbBdkQDdulxhK+jkk1KRW1A4gkMvIaw8uILLbkSTyUjW9qKpVyzScPfCvirMttio1zxOnA/diB5140TzR7jNMh3fbuPQxUyRxAfezaH7QCctrUrUaFEBjnvGwr3vjwA9OPQyoaMcuvGKutiZbr8HMQ69XfOhZ6hNDtaqSoaFbRBaK+jk4J2qtJKKRnDoZQSlo4yszISGTIPuzOIr9RG6nkTXiub9aEjzkNvqcnLoUZZLfbZFIR+7UZZcHUdaOSnKwIDHSFFrS8tD1i7qvMbpXG57RNz47c5D1zvxrCwhW/3JguQiI7aWa8hDDygF8j10+TsYIw2OEL3bL0Kso9H7WD2V49zf0JfTjrIieV6toqMcelVuQsEIDr2EsK54o3C0Tb4znuVifdGUb4QPh57cq42m1UGeMA+5tS7XXC4WjljQUouKk3O9SHO5GJFvej57ZXPFIeuqM6s8/jzEQ/R0dhInfrZhLhfXkcIzV2tR1SDdM7NFlSzvBwSHXibkeOryc+iGjhLnoZvwjtDVijNGiib0VSU4LDxCDyF61xEceglh5yzlrIWmRCNnOOGkLbJ5OXSfPHQtF+wzx7komajMzaE3D2FD1o9Dl/XaOHSPAN1vVKkUzbtYdA8OHWAnh+7zNdBz6MlnjqOqzLEANj21Flfa2JCHHtBx5Hnk8j6msUde+aJ1Gqkpcb0V6MS0HHpuO0qKos8jBOjdR3DoJYIZfaQ4WEtvfq08HolHeegZHLq0SLSDTo3Va4NWVh3JJ8qd86EntonsMyKaehkc49DFfon6tZqIInRhTVFjvz4P3cbbC7KCfRphex66kCli6WOQA3Tz+cvqYTGpPYecyKE3N6R7xuDAoQf0H9RD/0vUVG1XUJY1N4tt6H8V0HSK5bnHWlRlPp2iERx6iRCL8JJ5zCnhdDnXCeEkx5uVhx6P8t2cuR/3Kcm5t6Vy3zz0LA6/tj89V0tS1jZS1Mmh28YReHHoUstC73hd1zidh56UTbT6nHnoapMyx0ZIOmMtqrptLg7dqlNvZk8jOPQ+Rd4FLrqJdo3Q1I8U5dh2r8M1oKrsCAG6HcGhlwhmlkRqn5CJkI5M6znDGRy6WYmTQ/dIY1TnrDMnWgiSmD7sTs/lks0vN7Imkg47YWp0rdi4Qc48dIV9yWMs5qnK9XnoJn8dN9S5pihqI2klikObuZLW65ATsqCYzWtsv2cMe19LWFM0oNJQc+gleg/aFpT1+0jRNt/jtkTTVbkJBSM49BIhNtIwg0NvzhEd57+juVwyQnT9EnR2G632pwqkaDr5PspyJvw4dEeWS3QMx2yRomVKeHSRTzYE9HO5WMo8QnRtHrqTQ3fUL428tOrNYKqznklb/eY2G4e589DddlYZwaH3KdQ8cInehG5x6A1UlUNv9z1uS4BejVtQOIJDLxGa0Ueas3RFLebxjVGKSR40iViEbrFBOtYvD12Wi+XBK3lg93zoHrKxkYfGbIuC4mTnoTTC0qSm03O5uG2xKkqVChy/Va+8HT1nEoee+G0d+RrpNVuIdhmXDVY5YbvGocc/quk67XnoZQpM2ong0PsUeRe46CbaFZXpI/T22tFpSGmLxddT/AWryj0oGsGhlwhRtEPpZqoUtcRoyTrfWeOP4+VJxBa4cJDolgBIhD4PPT5ntUwZx/c4eXGB/7brbf4wOXL7qE2jXyIjOjQdfjvWFG0lI8Y21kBqmaRGH5NutkW/nHR9S8+8YZmtC4Hz78XBU3kQHHqfIvcCF11E2yJ0bU5+m+3oNFwfskLraYvOityEgqFy6ER0ERE9RUTriejjlv3vIaJH6//uIaLTize1f6CJ8OwceiNn2JPbFH5rtl075RGP8kjVpJymrFZ1Ipp3zftiZk0IEV3SVlsKfDqyrddtcTfidWsh6q7VKV07By+u7Ouo7atH84qRoplZLkoOPbnXmocujnK1O/sSxSVtRaZDJ6JBAF8CcDGAUwG8i4hOTYg9B+A8Zj4NwF8DuL5oQwOKhZbXLNOL0O35O5oRejWiQ2lyruLraa/+gCY0EfpKAOuZ+VlmHgVwE4BLTAFmvoeZd9U37wUwXKyZ/QFzpGFWk1KcbdGSiWbNcokrS9lg2ZWqL22Te9sslxapdsE1P0tSiSNvxcia4HruvuzYGraSkTZDxvFx2WY0q58P3cLbF8GhK25GcsCU9VijhSLU5LQvLqmTdeehN66xNHI1zLaYhUUAXjC2N9XLJPw3AD+07SCiK4loNRGt3rZtm97KgMKhjppK9Ca0j0PXoaocervvcTv47orcgsKhcejqrE4iOh81h/4x235mvp6ZVzDzivnz5+ut7BOY0UfWSEMxD53S/LE9yyVdb/J3Ur+kq7nPzvOn5eIGqPPQLfZEsoltdR465IyJhmySQ8+KbGt528k8dI9+ggLy0NPHpo8RzyP2PAgtDkv9WaZo+XapccHM0YabQ8/WWVUMKWQ2AVhsbA8D2JwUIqLTAHwFwMXMvKMY8wLahcChG3qVclnN/V5Dpzj0doTT1bgDxUMTod8PYDkRLSOiyQAuA3CrKUBEJwC4BcDlzLyueDP7A3EOPbFPEeHVflMqkrUdGYvQnRy6O2J37ZOj6UQeupijHS8nl2yi3DnvS3QM6hyxm0OP9CUjW6G/wcY5e/HiHrIiN59q0aVbYZTYth3bGEmryUPPCoNtrQSrXKqlF29RAcY9s9lu+eCWKf22nciM0Jl5jIg+DOB2AIMAbmDmNUT0wfr+6wD8BYBjAPxj/UKPMfOK9pkd0CrCSFFDr1IumTLX6+jYSNEe0VkFaCgXMPNtAG5LlF1n/L4CwBXFmtZ/iHW6ZXLojaglHotHfGdmHjpZ90tcvVuX28a0YCIPXanPmWuR4tvt636aFTbWv3E5NnOkZHO3FB02It90iK7Nn3fL6o632mbV7x7x2vit59DdXwbt2AiZQ09/VAOHHkcYKdqn8J2/pAxoX2Ss7U9ocOjtsqOzkKbPbVc9BWtth9KeR3DoJUJ8pCFZ90XblnLmJv+exXzLWS5Jm5Ka9NkJLrnY+SkjUlceekrWyaHXWzeMWN6yRFFTo3LDjsbxMVluCuhnW5Tt0whr+XZbNoqYT58c22CvImVrdh66DmJGFxtrigr3LOShB/QlejHKDBx6sejUeYQ89M4hOPQSwTXSUMpEsPOdmjx0iUOXMx4kXdpjTTkdhx7f4zPbojOa5+bfOIduJ9GjPPQoO8TtTgiW+ydm8sj2qWQFG9I8tBlJx88jm0NXrikq2GKVdTxIrnNKc+ju1qfeumogOPQ+he8qPWVA12dbrGiE3u573I7rFWZbtCM49BLBDDay8pijlzAR9TTYdxcvDuiyTGz1uuCKCpM6dXnocbhXLEqG6CRHesl+BwWHbnLyYoZF5PD1sy16XV9rA0K6xnJrKes80s+UTHHEovnMk3E/kzY9nDA8ZXvqWPvAr8ChB1QavnOAVxm+c8NXJTp0Lr1XaD1t0FmNW1A4gkMvIRo8uAlbJFIrj0c9zdV1TFkriZ7SlfxtO9YZhSn59gYvKxwmHp9cickt64r8uX5MLWvCtViyNQ9doCrMEZhZc/GYR9nqtEraslyUWuORdP2HwEMnOXQ489B1UXesXriuh+MYoz6Z/2frB6QfAhMgOPS+hTqzo61WlASeF6MqwWFz6H+b62kPiR5gQXDoJYIZfaSeVyETwZaHDlAqskki/pKZUZYj5UHQJR3riqYpsS0ojIM0PG1d1PHCmxlCpi3WUZv1iC/OoQvRodFpqmXRffLQrRy/GM27Dq4/ZwIRneStXRy6xhabCc7xDMKzW5u/vvZbat+ZfSI+tlUFwaH3KUKE3oR+5skmpVIFdCxCb4fOQKJbERx6iRCL8PLkoTNHc6lbOVMD8TU9ZVltbrn9WEGuEfaZ25KcgfQIWNlOzajSKA/dYW8j4jPrFjMsjP1aDr2lqBuuSNfO75t6pFzuZCTdnjx0h1zMVjtH7xqta11TtC9Ck+DQ+xYhQm/Cd+bJqsSGnYrQ23HBQoBuR3DoJYKZfJA1l0tD2hbBJDM8sjj0VFaDo14t9+mS1XLoqcwVxwQtKbqdskeVcj2x2eXYarIU09fk0JOtF7kPRNuycMqqC92tJW1fQOO3LWunqUsZdiP7mXTZ2ihPflTTrSS2fkAChx5Qaegj9Oq/Cb6rN1UlOHR1BrejnkJ1hhDdiuDQSwRzpGFqPm1BNpm8YOdvbSS6fX9Lc7kIPL8Nqjz05DEu2WQ076y/0bpBLG/ZzlEbHLrZBLLYEuPYk+MIvDh0qWWRP5q39qlII14TfLttBaaULoctThuydMb0G7MtCnqSrT+tbVVBcOh9Ct/5SyoNz2CvciNF23yPQx565xAceongHmloj36TEVJjLvWsl1Say8UnykrtUzqG1JzVIg8c30GOPPRURgzJ9qTz0BtD4C0RMNcjeKPurLlc6uJO+2zHNGX10PY/WNcUFUd/IiZLjqGiPqM/4/sdfTFI29o4vkl7Cf0YkPLQ+yEyCQ69b6Hm0PvgRfBdvakqwWGTtmgzhx6yXDqG4NBLhPhIw8S+1LYQTVJjVGM8skliwOSwHSG6xN3bkJaVoulElouYhx6HH4fuyHIxjmEAA0LEHelKcOjRpUvJNyJf6noeumZNUdeMheZvJ4eusMXW+ZqPQ2/WNyDYDqPF1Y8IDr2EcA3kyDxWXYdOrl3BW9uzFDzUaxxALG0xQ96WdlpG+HDordyudtzrXri+3UBw6CWCDz8dj2I4dnxqtkWLZi2Hnj7WwX06bCxCLrkSU5ZVst4oRo/x+ZKOpOOTFobwab2kbcmWtrfKJL3ydtZ5ZFti6tJF3UllkqhrLViwoh8DYU3RgBLCfCZ9Ahy9bHezXNqTm5xPv+aa2Ra4cMm2554VC3GBi5Scdbo4fT0FyCavUeDQ7QgOvUwQ3iwbd2zjRM0MANv+mE5tHrrORKu0GJXmlCNLmWmXenSmcUyDRIkKLDakHZ8wwlKoryYr703Z7cWh62Wbx8S7d119H5mZK8o607LCs54I0c3BXBy7+/bMpAbnn667P0L04NBLhsbDGIs2PcIRraw+y0VdtRfak/lAxu9ibTFdX2aE7hnPdmvUozpCR2v3y+teiDoosZ3fniojOPQSQeQVgdRbZxuhF0UhlB1lyXNzuEu02QkuWVeOtNMYR245YFmpJqN+RjyyFzn0hOOT5PN8/PK0LOLH+/DtHNUJmPW6Wkz6sDubb8+WTWVoGW2oeB66XRHD3mEaOPSAroCivzmjTbVcdzn0drDocQ7dp1WjkWnSWP3Hobd6tzzuhSAbOHQdgkMvEeQ86PRamrY5o+N8o122qdNebxYv7JMRoc0v10akrhGwtUg6/pZn56E3oj85fY9j+5PyyVhWvs5yI4TTdvtw6C3ISotEe43+tDyHOlm7TCpDq/FMR+9A4h4k6+CQ5RJQIjRe7rwZG1phPYfenjeh3dN7tINDV0foMO6jpn6FTDsgLXBhk+sYhy7IJotDhG5HcOglgjjSzsYdWzn0hrw9UyCu02RQ05ylpRpRV3Ofjm/Pn4fuqju9PysPPbpeifKUbCLKboyy9bo2ruyclN16XtwHzZZJ7a94HrHnIUOnh0m2sRNJpDK0GuX1fqFUP4bQokxdU72ZPY3g0EsGm8/y4oMLHkMnpYG1ijwqs47x7Xcwm/NZ9bomBkvJe2YodXvUY3aWSwfz0MUIPZnlImvt57nSVQ6diC4ioqeIaD0Rfdyyn4jo8/X9jxLRmcWbWn04eUVH5GwMfLTvF3Ta6s3KQHG9+1puXBsvuTI1bLK+UZnm3KRosXZ8Up+rLrlca3frHHqjZVLnoQXKRRNJR/sVtlhlBRlK1GmODGXLcdrRuv0wyRygcOhENAjgSwAuBnAqgHcR0akJsYsBLK//uxLAtQXb2Tew5qF7Hq/jhJVZLp71m3a0sj+Xztjv7ArUnDhRTXvU3Fe0FRr3MdOKLnLoQsdiSq5++rnr8eLQhSwXD539G58DQwqZlQDWM/OzAEBENwG4BMAThswlAG7k2mfwXiKaQ0QLmfmlog3+2bpt+N/ffyJbsAexff+ItZyIsO7l/XjTZ34WlW3Zczj6ffHnfo4BIoyOT0Rlm3YdjORf2HXQorP5+6Hnd0eypl4A2D8yhkEiTNQjnNvXbMEjL+y22vni7kOx7a/fuxE/fCz9CLy89zCmTR6Mlf3WZ36WehHHJtKu5uV9I7Hr0MCWvYdTzuDxF/dYZTftqtm5f2QMQNOxfenOZ/Av9z0fk31h50G86vjZNfnRuny9mnd++T8wONCsc2RsAkk0ZD/1oydx7V3PWG0xdQDA//jWI5ieuD4AsPPAaKrsvTfch0mD6bjs8Nh4quyiz/0CZOxr2Na4Dg3sODAaXbcXdx/Crx07M6Wrgc/8eB1u+OVzAICNO9LPGQBMGaqdy599+1HMmDJkrbMBArBm896o/vEoQq9lOF3xtdV1uZrxH//2Y7Fr9fzOgzjnFcdg2qT49fv07U/h+p8/K55Hp/HOsxbjije+onC9Goe+CMALxvYmAGcrZBYBiL3NRHQlahE8TjjhBF9bAQAzpwxh+QL5AetlLF8wEycfdxSA2kvwB6tOxLZ9Izj31+bhjie2pGTnzZyCXQePYHyi5kheedwsvOXVC7H70CiOGM59+YKZOHbWVOwfGcPQAGF8gnH5OUuw++ARzJo6hK37Dsdk33ba8Xhg4y7sHxnD3sNHcOYJczEyNoHt+0fw8t64w0/atHD2NJxw9HTsPDCKp7fuE+VWLDkaJx83Cz97ehte3HUo+mAkcfrwbFz06uOwdd8Ilh87CwdGxq1UzPIFM3HW0qPx0p7DeMW8GZg9bRKOmmZ/vBvXY9u+ERAB71gxjMGB9AepIfuO1w3j9OHZePD5XZg7fTLec/YJ2L5/JHaNG1ixdC5WLJ2LE+fPwNqX9uHcE4/Bms17seOA/WO9fMFM/ObJC/DAxp0Yn6h1+h0YtTs7AFh6zAzMmT4JG3YcxKHRcYxYHHcDK5ceg0mDhMlDA7HnpLHvva9fgr2HjuDA6BhWnXQsNmw/gNOG5+B7j2yOrvHyBTNxwckLcM6Jx+DtZxyPwYEB/MYr5+GXT28HUdwxL18wE8Nzp+OUhbPw3PaDeOdZi/H1ezfi8nOW4HP//jT2jRyJZF953CzMnzkFI2MTmDllEIeOjOPIGOOcE4/Gj594OXYerz5+Nt78qgX40eNbMMGMFUvm4r2vX4L9I0dSH4blC2bid05fhFMWzsKf3fIY3rFiGA8/vxvbhGCpW5g3c0pb9FJ2yhK9A8CFzHxFfftyACuZ+SpD5gcArmHmX9a3fwLgT5n5AUnvihUrePXq1QWcQkBAQED/gIgeYOYVtn2aTtFNABYb28MANueQCQgICAhoIzQO/X4Ay4loGRFNBnAZgFsTMrcCeG892+UcAHvawZ8HBAQEBMjI5NCZeYyIPgzgdgCDAG5g5jVE9MH6/usA3AbgLQDWAzgI4APtMzkgICAgwAZNpyiY+TbUnLZZdp3xmwF8qFjTAgICAgJ8EEaKBgQEBFQEwaEHBAQEVATBoQcEBARUBMGhBwQEBFQEmQOL2lYx0TYAG3MePg/A9gLNKRLBtnwItuVDWW0rq11A79u2hJnn23Z0zaG3AiJaLY2U6jaCbfkQbMuHstpWVruAatsWKJeAgICAiiA49ICAgICKoFcd+vXdNsCBYFs+BNvyoay2ldUuoMK29SSHHhAQEBCQRq9G6AEBAQEBCQSHHhAQEFAR9JxDz1qwulsgohuIaCsRPd5tW5IgosVEdCcRrSWiNUT0h922CQCIaCoR/YqIHqnb9VfdtikJIhokooeI6PvdtsUEEW0goseI6GEiKtVKMfUlKG8moifrz9zru20TABDRSfXr1fi3l4j+qNt2AQAR/XH9HXiciL5JRFNz6eklDr2+YPU6AG9CbVGN+wG8i5m7vsgoEf0GgP2ora366m7bY4KIFgJYyMwPEtEsAA8AeHu3rxvVFgGdwcz7iWgSgF8C+ENmvrebdpkgoo8CWAHgKGb+7W7b0wARbQCwgplLN0CGiL4G4BfM/JX6GgrTmXl3l82Koe5LXgRwNjPnHeBYlC2LUHv2T2XmQ0T0LQC3MfNXfXX1WoQeLVjNzKMAGgtWdx3M/HMAO7tthw3M/BIzP1j/vQ/AWtTWfO0quIb99c1J9X+liTCIaBjAWwF8pdu29AqI6CgAvwHgnwCAmUfL5szruADAM9125gaGAEwjoiEA05Fzxbdec+jSYtQBShDRUgCvBXBfl00BEFEaDwPYCuDHzFwKu+r4LIA/BZBeDbr7YAB3ENED9cXXy4JXANgG4J/rVNVXiGhGt42y4DIA3+y2EQDAzC8C+DsAzwN4CbUV3+7Io6vXHDpZykoT0ZUdRDQTwLcB/BEz7+22PQDAzOPMfAZq69CuJKJS0FVE9NsAtroWOu8yzmXmMwFcDOBDdcqvDBgCcCaAa5n5tQAOAChNXxcA1Gmg3wHw/7ptCwAQ0VzUmIZlAI4HMIOIfi+Prl5z6GEx6pyoc9TfBvANZr6l2/YkUW+W3wXgou5aEuFcAL9T56pvAvCbRPT17prUBDNvrv/dCuDfUKMjy4BNADYZLa2bUXPwZcLFAB5k5pe7bUgdvwXgOWbexsxHANwC4NfzKOo1h65ZsDoggXrn4z8BWMvMn+m2PQ0Q0XwimlP/PQ21B/vJrhpVBzP/GTMPM/NS1J6znzJzrqipaBDRjHrnNup0xpsBlCK7ipm3AHiBiE6qF10AoOtJCwm8CyWhW+p4HsA5RDS9/q5egFo/lzdUa4qWBdKC1V02CwBARN8EsArAPCLaBOAvmfmfumtVhHMBXA7gsTpfDQB/Xl8rtptYCOBr9YyDAQDfYuZSpQeWFAsA/Fvt3ccQgH9h5h9116QYrgLwjXrQ9SxKtGg8EU1HLUvu97ttSwPMfB8R3QzgQQBjAB5CzikAeiptMSAgICBARq9RLgEBAQEBAoJDDwgICKgIgkMPCAgIqAiCQw8ICAioCIJDDwgICKgIgkMPCAgIqAiCQw8ICAioCP4/Qa5SAs4xRN4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "imgP=D3[2,:];\n",
    "imgN=abs(D3[2,:]-1);\n",
    "curr=np.zeros((np.shape(imgP)[0]),dtype=\"int\")\n",
    "curr[np.where(abs(np.diff(imgN))>0.001)]=1\n",
    "\n",
    "niter=np.shape(current)[0]\n",
    "dt = 0.01;\n",
    "ti = 0.0; tf = niter*dt;\n",
    "t=np.linspace(ti,tf,niter)\n",
    "\n",
    "# plt.figure(figsize=(50,30))\n",
    "# #plt.plot(imgN,label=\"Negative\")\n",
    "# #plt.plot(imgP,label=\"Positive\")\n",
    "plt.plot(t,curr,label=\"current\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7fe84f7f9f40>"
      ]
     },
     "execution_count": 124,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAS4AAAD8CAYAAADJwUnTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAATDElEQVR4nO3df4wc5X3H8ffHBtQWrGLqg1r+UbvI/eFWxdArIJG2ThGpjdQapNJiKuIg0gMJV0HKHyD+KEhVJdoEEqISrAMsjESwkDDBRW4chNKaKAX5jFywcYGT48Jhyz9KFUiiCBm+/WPG6Xp3b3dud3ZnnvPnJa1uZ3b2medmzx8/88wzzyoiMDNLyZyqK2BmNlMOLjNLjoPLzJLj4DKz5Di4zCw5Di4zS46Dy8wGRtJmScck7ZvmdUn6hqRJSa9LuqxIuQ4uMxukJ4A1HV5fC6zIH2PAI0UKdXCZ2cBExC7ggw6brAOejMwrwPmSFnYr96yyKljEggULYtmyZcPcpdkZ5dChQ5w4cUL9lCFpJrfT7Ad+1rA8HhHjM3j/IuC9huWpfN2RTm/qK7gkrQEeAuYCj0XE/Z22X7ZsGRMTE/3s0sw6GB0dHfYufxYR/ey0Xch2Dc6eTxUlzQUeJjtHXQmsl7Sy1/LMrD4kFXqUYApY0rC8GDjc7U399HFdDkxGxMGI+BjYSna+amaJmzNnTqFHCbYDn8+vLl4J/CgiOp4mQn+niu3OTa9o3kjSGNnVApYuXdrH7sxsWEpqTSHpaWA1sEDSFHAvcDZARGwCdgDXApPAT4FbipTbT3AVOjfNO+rGAUZHRz2HjlnNlXgaSESs7/J6AHfMtNx+gqunc1Mzq7+ygmtQ+jlJ3Q2skLRc0jnAjWTnq2aWuCF2zvek5xZXRJyUtBHYSTYcYnNE7C+tZmZWmbq3uPoaxxURO8g618xslpBU1hXDgRnqyHkzS8OsbnGZ2ezk4DKz5Di4zCw5Di4zS4o7580sSW5xmVlyHFxmlhwHl5klperbeYpwcJlZCweXmSXHVxXNLDlucZlZUtzHZWZJcnCZWXIcXGaWHHfOm1lS3MdlZklycJlZchxcZpYcB5eZJcfBZWZJ8USCZpYkt7jMLDkOLjNLjoPLzJLiAahmliQHl5klx1cVzSw5bnGVqO4Hc1AiYmBlD/qYDrLuNhizvo9L0iHgI+AT4GREjJZRKTOrVt2Dq4wT2c9GxCqHltnscarV1e1RsKw1kt6SNCnp7jav/7Kkf5H0n5L2S7qlW5lJnSqa2XCU1TkvaS7wMHANMAXslrQ9It5s2OwO4M2I+DNJI8Bbkp6KiI+nrV+f9Qrgu5L2SBqbpuJjkiYkTRw/frzP3ZnZoBVtbRVscV0OTEbEwTyItgLrmrYJYJ6yAs8DPgBOdiq03xbXVRFxWNKFwIuS/isidp1Wo4hxYBxgdHTUPbVmCZhBH9cCSRMNy+P5v/lTFgHvNSxPAVc0lfHPwHbgMDAP+KuI+LTTTvsKrog4nP88Juk5snTd1fldZlZ3MwiuE136t9sV1NyA+VNgL/AnwMVkjaCXI+LD6Qrt+VRR0rmS5p16DnwO2NdreWZWHyWeKk4BSxqWF5O1rBrdAmyLzCTwQ+C3OhXaT4vrIuC5vPJnAd+KiO/0Ud5AeTxRe92OS7+XxQd5Wd2f6eCU+LntBlZIWg68D9wI3NS0zbvA1cDLki4CfhM42KnQnoMrIg4Cl/T6fjOrpzInEoyIk5I2AjuBucDmiNgv6fb89U3A3wNPSHqD7NTyrog40alcD4cwsxZltpQjYgewo2ndpobnh8m6mgpzcJlZi7qPnHdwmVkLB5eZJWXW32RtZrOTg2uW6OeDHPRl+yr/yPr53er+j+NM5okEzSw5df9PxcFlZqdxH5eZJcnBZWbJcXCZWXLcOW9mSXEfl5klycFVE90+iJSnSKlyLNVsPq5nMgeXmSXHwWVmyXFwmVlSypxIcFAcXGbWwi0uM0uOg8vMkuPgMrOkeABqyaocE+TxSO0N+uvNrBp1/9ySCi4zGw5fVTSzpPhU0cyS5OAys+Q4uMwsOQ4uM0uKb/kxsyS5xWUDV/c/sl557Fx16v431bU9KGmzpGOS9jWsu0DSi5LeyX/OH2w1zWyYTg2J6PaoSpET2SeANU3r7gZeiogVwEv5spnNEskHV0TsAj5oWr0O2JI/3wJcV261zKwqRUOryuDqtY/roog4AhARRyRdON2GksaAMYClS5f2uDszG6a6X1UceO0iYjwiRiNidGRkZNC7M7MS1L3F1WtwHZW0ECD/eay8KplZ1coMLklrJL0laVJS2/5wSasl7ZW0X9K/dyuz1+DaDmzIn28Anu+xHDOrmTL7uCTNBR4G1gIrgfWSVjZtcz7wTeDPI+J3gBu6ldu1j0vS08BqYIGkKeBe4H7gGUm3Au8W2ZH1bpBN8n7n06r7eB/rTYmf6+XAZEQczMvdSnZx782GbW4CtkXEuwAR0fUMrmtwRcT6aV66utt7zSxNM+icXyBpomF5PCLGG5YXAe81LE8BVzSV8RvA2ZL+DZgHPBQRT3baqUfOm1mLGbS4TkTEaKei2qxrbuafBfw+WWPoF4H/kPRKRLw9XaEOLjM7TclXDKeAJQ3Li4HDbbY5ERE/AX4iaRdwCTBtcNV7sIaZVaLEq4q7gRWSlks6B7iR7OJeo+eBP5R0lqRfIjuVPNCpULe4zKxFWS2uiDgpaSOwE5gLbI6I/ZJuz1/fFBEHJH0HeB34FHgsIvZNX6qDy8zaKPNqcUTsAHY0rdvUtPwV4CtFy3Rw1UCVwx3MmnkiQTNLUt3H5zm4zKyFg8vMkuPgMrPkOLjMLClVT1lThIPLzFr4qqKZJcctLhu4QY7V6rfsfv4B9PuPx2PYeufgMrOkuI/LzJLk4DKz5Lhz3syS4xaXmSXFfVxmliQHl5klx8FlA9fPH9mgxzr1U36//3jqfFzqzsFlZknxRIJmliS3uMwsOQ4uM0uOg8vMkuPgMrOkeACqmSXJVxWtq5TnvOpHt9+7zsdlto/zqnuLq2usStos6ZikfQ3r7pP0vqS9+ePawVbTzIbp1Olit0dVirQHnwDWtFn/tYhYlT92tHndzBJUNLSqDK6up4oRsUvSsiHUxcxqIvlTxQ42Sno9P5WcP91GksYkTUiaOH78eB+7M7NhmTNnTqFHZfXr8X2PABcDq4AjwAPTbRgR4xExGhGjIyMjPe7OzIYp+VPFdiLi6Knnkh4FXiitRmZWqapDqYieWlySFjYsXg/sm25bM0tP8i0uSU8Dq4EFkqaAe4HVklYBARwCbhtcFa2bKscU1XmsVKf3171FUbW6H58iVxXXt1n9+ADqYmY1kXxwmdmZJYWJBOtdOzOrRJl9XJLWSHpL0qSkuzts9weSPpH0F93KdHCZWYuygkvSXOBhYC2wElgvaeU02/0jsLNI/RxcZtaixBbX5cBkRByMiI+BrcC6Ntv9LfAscKxIoQ4uM2sxg+BacOrOmPwx1lTUIuC9huWpfF3jvhaRDavaVLR+7py3vnQbslDn4RLW3gzHaJ2IiNFOxbVZ1/zBfR24KyI+KbpfB5eZtSjxquIUsKRheTFwuGmbUWDrqRYccK2kkxHx7ekKdXCZWYsSx3HtBlZIWg68D9wI3NS4QUQsb9jvE8ALnUILHFxm1kZZwRURJyVtJLtaOBfYHBH7Jd2ev164X6uRg8vMTlP2fYj5RKM7mta1DayI+EKRMh1cZtbCt/yYWXLqfsuPg8vMTlP1lDVFOLhsoDy1TJrq/tk4uMyshYPLzJLj4DKz5Di4zCwpKUwk6OAysxZucZlZchxcZpYcB9cZoOoPuc7zUlV5bKr+XFLlAahmliR3zptZctziMrPkOLjMLCnu4zKzJDm4zCw5Di4zS07dryp2rZ2kJZK+J+mApP2SvpSvv0DSi5LeyX/OH3x1q9Pnt/nOWiV80/G0IqLjY5Cq3HfVin4ZbJV/+0Vi9STw5Yj4beBK4A5JK4G7gZciYgXwUr5sZrNA8sEVEUci4rX8+UfAAbKv0F4HbMk32wJcN6A6mtmQ1T24ZtTHJWkZcCnwKnBRRByBLNwkXVh+9cysCnXvAikcXJLOA54F7oyID4v+YpLGgDGApUuX9lJHMxuyugdXoUsHks4mC62nImJbvvqopIX56wuBY+3eGxHjETEaEaMjIyNl1NnMBujURIJFHlUpclVRwOPAgYh4sOGl7cCG/PkG4Pnyq2dmVZgNfVxXATcDb0jam6+7B7gfeEbSrcC7wA0DqeEs0O3yeb9/AJ3eP+h996POdTvT1f3Ydw2uiPg+MN1vcXW51TGzOkg+uMzszFL1aWARDi4za1H3W34cXGbWwi0uM0uOg8vMkuI+LjNLUt2Dq949cGeIlKdQ6Vb3uv5eda5bHZQ5AFXSGklvSZqU1DKLjKS/lvR6/viBpEu6lekWl5m1KOuqoqS5wMPANcAUsFvS9oh4s2GzHwJ/HBH/K2ktMA5c0alcB5eZnabkPq7LgcmIOJiXvZVsSqyfB1dE/KBh+1eAxd0KdXCZWYsZBNcCSRMNy+MRMd6wvAh4r2F5is6tqVuBf+22UweXmbWYQXCdiIjRTkW1Wde2E1HSZ8mC6zPddurgMrMWJZ4qTgFLGpYXA4fb7O/3gMeAtRHxP90K9VVFM2tR4lXF3cAKScslnQPcSDYlVuO+lgLbgJsj4u0ihbrFZWanOTWRYBki4qSkjcBOYC6wOSL2S7o9f30T8HfArwDfzMPwZJfTTwfXMHT7n6nbuKFBjisa9JilOs/3ZdMr83OLiB3AjqZ1mxqefxH44kzKdHCZWYu6j5x3cJlZCweXmSXFN1mbWZI8kaCZJcctLjNLjoPLzJLiPq5ZpNOYoEF+L2K/+v3uwiq/+9DjsKrj4DKz5Di4zCw5vqpoZklxH5eZJcnBZWbJcXCZWXIcXGaWnOSDS9IS4EngV4FPySbDf0jSfcDfAMfzTe/J590546Q83qjfuqf8u1t7ZU4kOChFWlwngS9HxGuS5gF7JL2Yv/a1iPjq4KpnZlVIvsUVEUeAI/nzjyQdIPvKITObpeoeXDNqD0paBlwKvJqv2ph/bfZmSfOnec+YpAlJE8ePH2+3iZnVTIlfljEQhYNL0nnAs8CdEfEh8AhwMbCKrEX2QLv3RcR4RIxGxOjIyEj/NTazgSoaWlUGV6GripLOJgutpyJiG0BEHG14/VHghYHU0MyGru6d811rpyxWHwcORMSDDesXNmx2PbCv/OqZWRVmQ4vrKuBm4A1Je/N19wDrJa0i+zrtQ8BtA6ifmVWg7p3zRa4qfh9o91uckWO2zGa7qltTRXjkvJm1cHCZWXIcXGaWlNlyy4+ZnWHc4jKz5Di4zCw5Di4zS46Dy8yS4nFcZpYkX1U0s+S4xWVmyal7cNW7PWhmQ1f2fFyS1kh6S9KkpLvbvC5J38hff13SZd3KdHCZWYuygkvSXOBhYC2wkmxWmZVNm60FVuSPMbJJSjtycJlZizlz5hR6FHA5MBkRByPiY2ArsK5pm3XAk5F5BTi/ab6/FkPt49qzZ88JSf/dsGoBcGKYdZiButatrvUC161XZdbt1/otYM+ePTslLSi4+S9ImmhYHo+I8YblRcB7DctTwBVNZbTbZhH5l/S0M9TgiojTJp2XNBERo8OsQ1F1rVtd6wWuW6/qVreIWFNice3OJ5u/jLPINqfxqaKZDdIUsKRheTFwuIdtTuPgMrNB2g2skLRc0jnAjcD2pm22A5/Pry5eCfwo/z7XaVU9jmu8+yaVqWvd6lovcN16Vee69SUiTkraCOwE5gKbI2K/pNvz1zeRTQN/LTAJ/BS4pVu5iuh4KmlmVjs+VTSz5Di4zCw5lQRXt1sAqiTpkKQ3JO1tGp9SRV02SzomaV/DugskvSjpnfzn/BrV7T5J7+fHbq+kayuq2xJJ35N0QNJ+SV/K11d67DrUqxbHLSVD7+PKbwF4G7iG7DLobmB9RLw51IpMQ9IhYDQiKh+sKOmPgB+TjSr+3XzdPwEfRMT9eejPj4i7alK3+4AfR8RXh12fprotBBZGxGuS5gF7gOuAL1DhsetQr7+kBsctJVW0uIrcAmBAROwCPmhavQ7Ykj/fQvaHP3TT1K0WIuJIRLyWP/8IOEA2ErvSY9ehXjZDVQTXdMP76yKA70raI2ms6sq0cdGpMS75zwsrrk+zjfkd/purOo1tJGkZcCnwKjU6dk31gpodt7qrIrhmPLx/yK6KiMvI7li/Iz8lsmIeAS4GVpHdZ/ZAlZWRdB7wLHBnRHxYZV0atalXrY5bCqoIrhkP7x+miDic/zwGPEd2alsnR0/dOZ//PFZxfX4uIo5GxCcR8SnwKBUeO0lnk4XDUxGxLV9d+bFrV686HbdUVBFcRW4BqISkc/NOUySdC3wO2Nf5XUO3HdiQP98APF9hXU7TNBXJ9VR07JRNFPU4cCAiHmx4qdJjN1296nLcUlLJyPn8cu/X+f9bAP5h6JVoQ9Kvk7WyILsd6ltV1k3S08BqsmlPjgL3At8GngGWAu8CN0TE0DvJp6nbarLTnQAOAbd1u+dsQHX7DPAy8Abwab76HrL+pMqOXYd6racGxy0lvuXHzJLjkfNmlhwHl5klx8FlZslxcJlZchxcZpYcB5eZJcfBZWbJ+T/0nQ2VGtm54AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(curr.reshape(28,28),cmap=\"binary\")\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Time     0.088182\n",
       "Hs       0.726915\n",
       "Hk       0.383324\n",
       "Hs+Hk    1.110239\n",
       "Name: 65, dtype: float64"
      ]
     },
     "execution_count": 125,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "opt_relv(curr,dt,200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 143,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cal_Orel(D):\n",
    "    no_of_images=np.shape(D)[0]\n",
    "    rel=[]\n",
    "    \n",
    "    for ig in range(no_of_images):\n",
    "        img=D[ig,:];\n",
    " \n",
    "        spk=np.zeros((np.shape(img)[0]),dtype=\"int\")\n",
    "        spk[np.where(abs(np.diff(img))>0.001)]=1\n",
    "\n",
    "\n",
    "        rel.append(opt_relv(spk,dt,200).Hk)\n",
    "    return rel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 146,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(\"0.dat\",cal_Orel(D0),delimiter=\"\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 2 µs, sys: 0 ns, total: 2 µs\n",
      "Wall time: 3.81 µs\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "#lMSRD0=cal_Orel(D0);\n",
    "# lMSRD1=cal_msr(D1,nature);\n",
    "# lMSRD2=cal_msr(D2,nature);\n",
    "# lMSRD3=cal_msr(D3,nature);\n",
    "# lMSRD4=cal_msr(D4,nature);\n",
    "# lMSRD5=cal_msr(D5,nature);\n",
    "# lMSRD6=cal_msr(D6,nature);\n",
    "# lMSRD7=cal_msr(D7,nature);\n",
    "# lMSRD8=cal_msr(D8,nature);\n",
    "# lMSRD9=cal_msr(D9,nature);"
   ]
  },
  {
   "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": []
  },
  {
   "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": "markdown",
   "metadata": {
    "id": "eNwZ1nQvwJU5"
   },
   "source": [
    "# Solving ODE using Euler's Method "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {
    "id": "OUxuE4Thtg4R"
   },
   "outputs": [],
   "source": [
    "def IZH_NeuModel(Indx,niter,dt,current):\n",
    "    a,b,c,d,I,v=pars[name[Indx-1]];\n",
    "    u=b*v;\n",
    "    \n",
    "    if niter!=np.shape(current)[0]:raise Exception(\"Sorry, Check current array!!\") \n",
    "    \n",
    "    state=[[v,u]] #Initial condition\n",
    "    \n",
    "    for t in range(1,niter+1):\n",
    "        dv,du=((0.04*(v**2)) + (5*v) + 140 - u + current[t-1] ,a * ((b*v) - u) )\n",
    "        v+=dt*dv;u+=dt*du\n",
    "    \n",
    "        if(v>=30.0):\n",
    "            v=c\n",
    "            u+=d\n",
    "        state.append(np.array([v,u]))\n",
    "\n",
    "    return np.array(state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tonic spiking\n",
      "CPU times: user 15.8 ms, sys: 0 ns, total: 15.8 ms\n",
      "Wall time: 14.7 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fe84fbe80d0>]"
      ]
     },
     "execution_count": 110,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAApzUlEQVR4nO3dd3xc9ZX38c9RsyW5SJYsW7ZsuXcMBmGMTceAARNIW0hClpZ1CiRhE9iEAMmm8GSX8JDyZJMsSwnJQgiETgBTDIRmjGzcLbnbsmR1q/eZ8/wxI3lkj6wyc3VHM+f9evllaWY0czD3fnXu7/7u74qqYowxJjrFuV2AMcYY51jIG2NMFLOQN8aYKGYhb4wxUcxC3hhjoliC2wUEyszM1ClTprhdhjHGDCnr16+vVNWxwZ6LqJCfMmUK+fn5bpdhjDFDiogc6Ok5G64xxpgoZiFvjDFRzELeGGOimIW8McZEMQt5Y4yJYhbyxhgTxSzkjTEmilnIG0e9uaOMkppmt8swJmZZyBvHeL3KTY/mc/UDH7pdiolQ6w8c4ddv7HK7jKgWUVe8muhS39oBQFldq8uVmEj12d9/AMC3l890uZLoZZ28cUxdczsAo5MTXa7EmNhlIW8cU+sP+VHD7YDRGLdYyBvHdHbyo6yTN8Y1FvLGMXUtnZ28hbwxbrGQN46ptU7eGNdZyBvH2Ji8Me6zkDeOqWv2TaFMHWYhb47X1uF1u4SYYCFvHNPZyRsTjG0fg8NC3jjGdmJzIrXNbW6XEBMs5I1jOmfXqKrLlZhIVNNkTcBgsJA3jrFO3pzIEQv5QRFyyIvIJBF5S0R2iMg2Efm2//ExIvK6iOzy/50eerlmKOkMeWvkTTA1TTZcMxjC0cl3AN9V1bnAEuBmEZkHfB94U1VnAm/6vzcxpHN2jTHB2JHe4Ag55FX1sKpu8H9dD+wAJgJXAo/6X/YocFWon2WGDlXtWtbAGnkTjI3JD46wjsmLyBRgEfARME5VD4PvFwGQ1cPPrBKRfBHJr6ioCGc5xkUt7V7aPDYP2vTsiA3XDIqwhbyIjACeBm5V1bq+/pyqPqCqeaqaN3bs2HCVY1zWObMGbEzeBFdjwzWDIiwhLyKJ+AL+MVV9xv9wmYhk+5/PBsrD8VlmaLBDcdObWttGBkU4ZtcI8BCwQ1XvD3jqBeA6/9fXAc+H+llm6KhuPHoorjYqb4KoCbgYyq6lcE44FhVZBnwZ2CIiG/2P/QD4D+BJEbkJOAh8PgyfZYaIwPFW239NMEcarZMfDCGHvKq+B0gPT18Y6vuboSmwkzcmmMAplKogPaWICYld8WocccQf8smJ8S5XYiJRu8dLQ2sHCXGW7E6zkDeOqG5qY+SwBJISbBMzx+s8MZ+W4ruhjI3oOcf2QOOII41tpKcmAXZSzRyvczgvI3WYy5VEPwt544jqpnbSU5NsnNUEVdXQCkDGCGsEnGYhbxxxpLGNMXYobnpQ5e/kM0dYJ++06Lkv2623wsaNbldh/KoXrWJmXRGSNg2e/Qh+eYvbJZkIUjVuEUxdTsarL0B2Hrp8OWiML4Nxyinwq1+F/W2tkzeOOJKQzJj2ZgC0xxm2JlZVJaYQp17SO1rcLiXqRU8n78BvQDMwLe0emu5+lfQbrkXe2wdnXwVX3eV2WSaCVD6zhTHbS4m76UZ4bSf6+htgM7EcYf+qJuw6Z06M6ZxdY6Py5hhVDa1kpA5D7My84yzkTdh1hnx6SpIN1JigqhvbumbWgDUCTrKQN2HXuW5NVydv+685RlVjGxkBM2tsG3GOhbwJu6qGoyFvR+MmmMqGVjJSk3p/oQmZhbwJu0r/hS5jR/o6NWvSTKDWDg/1LR1kWBMwKCzkTdhV1LeSlBDHqOEJ9LxAqYlVXUsajBiG2PbhOAt5E3YV9a2MHeGbOSFi462mu87hvG4nXm0bcYyFvAm7iobWrqEaH9uDzVGdw3mZI2y4ZjBYyJuwq6hv7VqTxPZhc6zyOl/IZ40c3vWYTaF0joW8CbvKYzp5OxQ3gcrqfEsZZI0aZk3AILCQN2HV4fFS1djWFfJ2OG6OVV7fSlpKIsMSjt41zBoB51jIm7CqbmxDFevkTY/K6loY5x+qsSbAeRbyJqzK6/1z5LvG5G0vNt2V1beSNar7OvLWBzjHQt6E1bEXQoGdVDPdlde1dJ10tSbAeRbyJqwqju3kbR82AbxepaK+lXHHdvI2pucYC3kTVp3DNZkj7UIXc7zqpjY6vMq4UTYmP1gs5E1YldW1MGp4AilJvvvR2D5sAnVNnxxpY/KDxULehFVpbQvZo5O7PWY7sOnUdSHUqOG9vNKEi4W8CavSuhbGjz66A9udf0yg8voeOnnrBBxjIW/C6nBtC9mju3dptgObTiU1LYgQMCbvbwJsG3GMhbwJm7YOL5UNrd06eWMCldQ0kzVyGEl20+5BY//SJmzK61tQhfGjAodrbJ68OaqktpkJaUfP2XQO5tk24hwLeRM2pbW+8Vbr5E1PSmpauoe8nbJxnOMhLyIrRKRQRHaLyPed/jzjnlL/9LjA2TUi2HirAXwXPBXXNDMxLTnIcy4UFCMcDXkRiQf+C7gUmAd8QUTmOfmZxj09dfK2/xqAqsY22jq8TAicfeViPbEiweH3XwzsVtW9ACLyBHAlsD3cH7R2LVx9NbS1QUYG1NbC6NFH/66u9nUL9pxzzxWV5xCXG0fbdxPA36x1NCax5uFscr4LHk/k1m7P9f+5/r6+oiqBmrgz8ZzV0W3fbSkezaIFCYxIjYz/Ljeey8iAtDT47/+G2bPDHI6q6tgf4HPAgwHffxn47TGvWQXkA/mTJ0/WgZozR9X3z2d/3P5z771H/79MvWyX6/XYn8j6c/YFbV3bxx/f36cJY+pcrylS/lx22cDyD8gPzNXAP0538sGOxrTbN6oPAA8A5OXlaZDX98kjj1gn7/Zze0qamXxKLTfcML7r/8v400sZ2T6Squ3jrJOPsuf6+/ri8naaEuu5774R3fbdjMs2Mey9pYxMjYuI/y43O/n77x9oAvbM6ZA/BEwK+D4HKHHig5YsgQMHnHhn01en/vQ9ViwYT2bm0ZBPTG3nvBtL+PU141yszESCn760i8c/Osjpp1zS9ZgIDJ9Yy/ot7WSMGHaCnzYD5fTsmo+BmSIyVUSSgGuAFxz+TOOCxtYOqhvbmJSe0u3xYIdyJjYdrG4iJz056FIXAz6EN71ytJNX1Q4RuQVYDcQDD6vqNic/07jj0JFmAHLSbXqcCa6ouoncDGsCBpvTwzWo6svAy05/jnHXoSNNwPEhbwuUGQBV5WB1E0unZ/bw/CAXFEPsilcTFkc7+ZTjnrP911Q2tNHU5mHymGOO9KwJcJyFvAmLQ0eaGJ4YR+aIpG6P2y5sAA5WNwKQm5Ea9Hlbu8Y5FvImLA4daSYnPSX4STU7Fo95B6p8w3mTxvQwJm+biGMs5E1YFB1pCromibXyBnwza0Rg0rHDNcZxFvImZKrK/sompmYefyguWJNm4GBVE9mjhjMsIb7b43bPEOdZyJuQVTS00tDaETTkjQHYX9V43FCNGRwW8iZk+yt9461TgnXyvruGmBimquypaGR61ojjnhP/eJ6dtnGOhbwJ2b7KBgCm9dDJ28yJ2FbV2EZtczszxgYJeTtn4zgLeROyvZWNJMXHdbvjTyfbh82ecl8TEKyT72SNgHMs5E3I9lc2Mjkjhfi44JFuh+KxbXeFP+THBj8xb5xlIW9Ctq+ysceTrnY4bvaUN5KcGM+E0T1Pn7RGwDkW8iYkXq+yvyr49MlOtgPHtj0VDUwbm0pckCM9awKcZyFvQlJS20xbh7fnTt4OyGPe7vIGpgc56RrI+gDnWMibkOyr9K1JcsJO3nbhmNXc5qG4ppkZPZx0tSbAeRbyJiSdMyd6mj5ph+OxbU/XSddeOnkb03OMhbwJSWFZA+kpiYwd2fOt22z/jV27/U1AT528NfLOs5A3ISksrWPWuJF2cxAT1PbDdSTFxzEtyPTJQNYIOMdC3gyYqrKzrIHZ40ee+HWDVI+JPDsO1zFz3AgS44NHjbUGzrOQNwNWUttCQ2vHCUPeOvzYpapsL6ljXvYot0uJaRbyZsB2ltYDMHvcCUIeOxSPVRX1rVQ1tjH3BCHf2QTYNuIcC3kzYAX+kJ95gpA3sWv74ToA5k2wTt5NFvJmwHaW1ZM9ejijkxN7fI2vUbM2LRZ1hvzc8Sfo5P1/27UUzrGQNwNWWFrf60lXE7t2HK5nYloyo1N6bgKM8yzkzYC0dXjZXd77zBoRG2+NVdtLapmb3fv2AbaNOMlC3gzIzrJ62jxeTpo4utfX2v4be+pa2tlb2cjCnLQTvs4mXznPQt4MyOZDtQAsnJh2wtfZ2iSxaXNRLapwyqS0Pr3eGgHnWMibAdlSXMvo5EQmjel5jfBOti5J7NlYdASAk3sJeWsCnGchbwZkS3ENJ00c3evFTnY4Hps2FtUwbWzqCWdeBbJGwDkW8qbfWjs8FJbWs6AP4/Fgh+KxRlXZWFTTp6EaawKcZyFv+m1naQPtHmVhTu8hb/tw7Dl0pJnKhjYW9XE8HqwRcJKFvOm3zcU1AH2aWQM2PS7WfFJUA8Apk9LdLcQAFvJmADYcqCEjNYmc9N5PutrxeOzJ319NSlI8c3qZIx/IGgHnhBTyIvILESkQkc0i8qyIpAU8d4eI7BaRQhG5JORKTcRYf6Ca03LT+7zCpO2/sWXt3ipOy03vcXnhQEe3IdtKnBJqJ/86sEBVFwI7gTsARGQecA0wH1gB/E5E4kP8LBMBKupb2V/VRN6Uvh2KWx8fW6ob29hZ1sCSaRlul2L8Qgp5VX1NVTv8364FcvxfXwk8oaqtqroP2A0sDuWzTGRYf8A3//m03DF9er1vWQPr0mLFun1VACyZ1sftw/+3bSLOCeeY/I3AK/6vJwJFAc8d8j92HBFZJSL5IpJfUVERxnKME/L3V5OUEMeCibZ8rDne2r3VDE+M46ReroQ2gyehtxeIyBvA+CBP3amqz/tfcyfQATzW+WNBXh/0d7WqPgA8AJCXl2e/zyNc/oEjnJwzmmEJfRt9s+Ga2PLRPt/5mqSEvvWPXQuUOVhTrOs15FV1+YmeF5HrgJXAhXr0uPwQMCngZTlAyUCLNJGhuc3DtpJabjprmtulmAhU1dBKQWkd31k+y+1STIBQZ9esAL4HfEpVmwKeegG4RkSGichUYCawLpTPMu77eH817R7t83gr+GZP2HhrbHh3VyWqcO7ssX3+mc61a2wbcU6vnXwvfgsMA173T4Vaq6pfU9VtIvIksB3fMM7NquoJ8bOMy97fU0livLB4at9D3sSOd3ZWMCY1iQUT+naRnBkcIYW8qs44wXP3APeE8v4msnywu4pFk9NJSer7ZiPYrd1igder/GNnBefMzCQuru9nYo6Oyds24hS74tX0SU1TG1tLalk2PbPfP2uH4tFvW0kdVY1t/RqqATsxPxgs5E2frN1bhSosm9G/i1xsVYPY8HZhOQBnz+xfyHeyRsA5FvKmT97dVUlqUnyvN4EIxnbg6Ld6eymLJqeROWJYv37OmgDnWcibXqkqawrKOWtmZp/WIwlkd/6JfkXVTWwtruPSBcEup+kbawScYyFverX9cB2Ha1u4cO64Af28nVSLbqu3lQKwYn72AH7amgCnWcibXr25oxwROH92Vv9/2PbhqPfq1lLmZY9ickbKgN/DGgHnWMibXr25o4yTc9IYO7J/462d7FA8epXXtbD+4BFWDHCoxsbknWchb06ovK6FTYdqWT53AF081shHu9XbSlElpPF4sEbASRby5oTWFPimxg10PF7EFp+KZs98UszscSOZkTViQD9vTYDzLOTNCb28tZSc9GTmjO/7rdxMbNhT0cAnB2v47GkT+3yXMDP4LORNj6oaWnl/dyVXnDxhwDuxYK18tHp2QzFxAledEvRWEX3SuV3ZcI1zLORNj17ZWorHq1yxcILbpZgI4/Uqz35SzNkzx5I1arjb5ZgTsJA3PXpxUwnTx6YyN3vgQzW+MXlr06LN2r1VFNc089nTcnp/8QkcvY23bSNOsZA3QZXWtrBuf3VIQzUmej227iCjhidw8byBnZA3g8dC3gT10uYSVGFliEM1vht5h6koExFKa1tYvbWUq0+fxPDEvt0GsiddSw3bNuIYC3lzHFXlyfwiTp6UNuCpcd3eLww1mcjx+LqDeFS5dkmu26WYPrCQN8f5pKiGnWUNXHP6pN5f3AtboCy6tHV4efyjg5w/O4vcjNSQ389u5O08C3lznL+uKyIlKZ4rTg7PrBq1Y/Go8fctJVQ2tPLlM8PTxVsT4DwLedNNQ2sHL24u4YqFExgxLNRbANvaJNHE61X+8PZeZo0bwbkDvDlIT6wRcI6FvOnmhY0lNLV5+KcwDNV0st03OrxZUE5hWT1fP296v+7jekLWBDjOQt50UVUeeX8f8yeM4tTJaW6XYyKIqvJfb+0mJz3ZkYvjrBFwjoW86fLe7kp2lTdw47KpYZ0bb0fiQ98He6rYWFTDV8+dTkI/7w52ItbIO89C3nR56L19ZI4YxsqTB3KHn+DsQqqhT1W599UCJowezudDvMK1589w5G0NFvLGb3d5A28XVvDlJbkMSwjtApdAgh2KD3Wvbi1l06Fabr1oVsgXPx3raBNgW4lTLOQNAH94Zw/DEuL40pLJbpdiIkiHx8svXitkRtYIPrNo4KtNGvdYyBuKqpt49pNivrB4MpkjBnaLv56IYMfiQ9jj6w6yt6KR2y6eHdax+E5dfbxtIo6xkDf87u09xIvwtXOnu12KiSBVDa3ct7qQpdMzuGS+LUQ2VFnIx7iSmmb+tr6Ifzo9h/Gjw78uuI3JD133vlpIU5uHn1w537ET6LasgfMs5GPcL1/fiWBdvOnuk4NH+Gt+ETeeNZUZWXbrx6HMQj6GFZTW8fSGQ/zzmbnkpKc48hkiYuOtQ0xrh4c7ntnCuFHD+NaFMx39rM61a2wbcU7oi5OYIes/XylgxLAEbrlghtulmAjymzd3UVBaz8PX54Vl/SLjrrB08iJym4ioiGQGPHaHiOwWkUIRuSQcn2PC54M9lbxVWME3zp9BWkqSY5/jG5O3Nm2o2FRUw+/f3sPnT8vhgjnOn2w9etMQ20acEvKvaRGZBFwEHAx4bB5wDTAfmAC8ISKzVNUT6ueZ0HV4vPz0pR1MGD2c65dOcfzzbP8dGpraOvjOkxsZN2o4d62cNyifaddDOy8cnfwvgX+j+wnyK4EnVLVVVfcBu4HFYfgsEwZ/+vAAOw7XcffKeWG/gvFYtqrB0PHD57ext7KR+z5/MqOTEwf1s60PcE5IIS8inwKKVXXTMU9NBIoCvj/kfyzYe6wSkXwRya+oqAilHNMHZXUt3P/6Ts6dNZYVC8YPymdaJx/5nsov4m/rD/HNC2aybEZm7z8QLtYEOK7X4RoReQMIlgZ3Aj8ALg72Y0EeC7qrq+oDwAMAeXl5FgcO+9nfd9Dm8fLjTzk397k724sj3Y7Dddz9/FaWTBvDtx2eTdMTawSc02vIq+ryYI+LyEnAVGCTPyxygA0ishhf5x5414kcoCTkak1IXttWyoubSvjX5bOYkhn6/Tn7yvbfyFVR38pXHs1ndHIiv7lmEfHhuhlIH9nt/5w34OEaVd2iqlmqOkVVp+AL9lNVtRR4AbhGRIaJyFRgJrAuLBWbATnS2MYPnt3KvOxRfP28wbvwycbkI1dLu4ev/jmfqsZWHvzn08kaFf4rnvvKZmA5x5FJsKq6TUSeBLYDHcDNNrPGXT98YRu1zW38+abFJCUM3jVwvvXJbAeONF6v8r2nN7PhYA2/+9KpnJQz2pU6bKVh54Ut5P3dfOD39wD3hOv9zcA9v7GYFzeVcNvFs5ibPcrtcozLVJWfvLSd5zeWcPsls7nspPDdJMZEHlvWIMrtq2zkB89sIS833ZX1aWy4JvL86o1d/PGD/dy4bCrfGMShu2CskXeehXwUa2n3cPNjG0hMiOM3X1jkyHrgZmj5n3/s5ddv7uLzp+Vw1+Vz7faMMcAWpohi9/x9B9sP1/Hw9XlMSEt2pQbBFiiLFL9ds4v7XtvJ5Sdl8/PPnETcIM+kCabzl4xtI86xkI9ST35cxJ/XHuBfzp46KGuQmMilqty7upDfv72HTy+ayC8+t9CO6mKIhXwUyt9fzZ3PbeGsGZl8b8UcV2sRselxbvJ4lR+/uI0/fXiAL54xmZ9duSAiOvhOR28aYtuIUyzko0xxTTNf+9/1TExL5rdftHH4WNbY2sG3/vIJbxaUs+qcadxx6Rwbg49BFvJRpLa5nZv++DGt7V6eWHW6o0sI95WIjbe6obS2hZse/Zgdh+v46VUL+PKSXLdLCspu5O08C/ko0dLu4V8ezWdPRQOPXL+YGVkj3C7JuOTj/dXc8vgGGlo6eOj60zl/dpbbJRkXWchHgQ6Pl2/95RM+PlDNb65ZxFkzB3EVwV4IYqOtg0RVeei9ffz8lQImpSfzxxsWR/zFb3Yjb+dZyA9xXq9y13NbeW17Gf9+xTyuOHmC2yUdx5Y1cF5dSzvff3ozL28p5eJ547jvn05m1PDBXRN+YOwcgdMs5Icwr1e5+/mtPPFxEd+8YAbXL5vqdknHs33YcR/uqeK2pzZRWtfCHZfOYdU504bcCVZrBJxjIT9EdQb8Yx8d5BvnTec7F81yu6Qe2e7rjJZ2D/etLuTB9/YxJSOFJ796JqflprtdVr8Msd9FQ5KF/BB0bMDffsnsiO3cIrOqoe+Tg0f43tOb2VnWwLVLJvODy+aSkjR0d2drBJwzdLeKGNXu8fK9pzfzzIbiiA948F+2bntw2NQ2t/OL1QU89tFBskYO45Ebhvbsma4t17YRx1jIDyHNbR5ufnwDawrK+c5Fs/jmBTMiOuBN+KgqL20+zE9e2k5VQyvXnTmF7148i5FD4uSqcZOF/BBR09TGjX/8mI1FNdzz6QV86YzIvLjlWII1aaHacqiWn/59O+v2VXPSxNE8fN3prt3kI9y6FiizrcQxFvJDwIGqRm56NJ+DVU381xdP5VK7yUNMOFzbzC9WF/LMhmIyUpP42VUL+MLiyYN+H1YztFnIR7iP9lbxtf9dj1fh0RsXc+b0DLdL6hffsgbWpfXHkcY2Hnh3L4+8vw+vwtfPm843zpselUMztqyB8yzkI9iT+UXc+ewWJo1J4aHrTmdqZqrbJRkH1TS18eC7+3jk/X00tXv41MkTuO3i2Uwak+J2aWYIs5CPQO0eL//5SgEPvrePs2dm8tsvnsro5KHZxdmYfO9qm9t55P19PPTuPupbO7h8YTa3XjiTmeNGul2a47qWNbCNxDEW8hGmvK6FWx7/hHX7q7nuzFzuXjnPlguOUsU1zTz83j6eWHeQxjYPK+aP59aLZjJnfGSvN2OGFgv5CLJ2bxW3PP4Jja0d/PqaU7jylIlulxQyEbv937G2l9TxP+/u5cVNJShwxcJsVp0znXkTYi/chc7ZNcYpFvIRQFX573/s5RerC8nNSOHxfzmDWTFwqB5LPF5lTUE5f/pwP+/uqiQlKZ7rlk7hxrOmMtGl+++a2GAh77KK+lZue2oT7+ys4PKF2fznZxcyYlj0/G/xjcnHbp9WUd/Kk/lFPP7RQYprmhk/aji3XzKba8/IZXTK0DzPEk5Hx+RjdxtxWvSkyRC0pqCM25/aTENrBz+5cj5fXpIblVewxtr+q6rkHzjCnz88wCtbD9PuUZbNyODulfNYPjfLzrGYQWUh74KWdg//5+Ud/OnDA8wZP5K/rFoSvcMz0fc7q0eltS08veEQT68/xN7KRkYOT+DaJbl86Yxcu1NXL2KsDxhUFvKDbGtxLf/6143sKm/gK2dN5fYVsxmWEO92WY6K5k6+pd3D69vLeGr9Id7bVYFXYfHUMXztvOmsXJg9pFeGHAxReOAacWwLHCStHR5+u2Y3v3t7DxmpSfzpxsWcM2us22U5TqKwlfd6fcMxz28s5sVNJdS1dDAxLZlbzp/BZ0/LITfDLlrrr2huBNxmIT8INh+q4fanNlNYVs9nTp3ID1fOIy0lye2yBkW0dGqqypbiWl7cVMJLmw9zuLaF4YlxrJg/ns/nTeLMaRnE2Zoy/RaNTUCksZB3UGuHh9+8uYs/vLOXzBFJPHx9HhfMGed2WaYfCkvreXFTCS9uLuFAVROJ8cK5s8by/UvnsHzuOFKjaCaUu6yVd4ptoQ7ZVFTDbU9tYld5A58/LYe7Vs4bsksThEIYWtPjVJXth+tYva2MV7ceZmdZA3ECS6dncvN5M7hk/nib+hhGtqyB8yzkw6y+pZ3/+9pO/vThfrJGDh/yd+6JBR6vsuHgEVZvLWX19lKKqpsRgbzcdH5y5XwuXZDN2JHD3C7TmAEJOeRF5JvALUAH8HdV/Tf/43cANwEe4FuqujrUz4pkqsorW0v58YvbKK9v5dozcrl9xWxGReHysP0RqXf/a+3w8MGeKl7bVsrr28uobGgjKT6OpTMy+MZ5M1g+d5wF+yDo6uTdLSOqhRTyInI+cCWwUFVbRSTL//g84BpgPjABeENEZqmqJ9SCI1FRdRM/fH4rbxVWMC97FH+49jQWTU53uyxzjNrmdt7dVcFr28p4q6Cc+tYOUpPiOW9OFpfMH8/5s8dG5ZrtJraF2sl/HfgPVW0FUNVy/+NXAk/4H98nIruBxcCHIX5eRGn3ePmfd/fymzd3ESfCXZfP5fqlU+yKxgCCewuUqSp7KhpYU1DOmzvKyT9wBI9XGZOaxKUnjeeS+eNZNiOT4YnRfZ1CJOtaoMxaeceEGvKzgLNF5B6gBbhNVT8GJgJrA153yP/YcURkFbAKYPLkySGWM3g+3l/Nnc9uYWdZA5fMH8ePrpjPBFtoynUt7R4+2lfNWwXlvFlQRlF1MwBzxo9k1TnTuHBOFqdMSrNfxCZm9BryIvIGMD7IU3f6fz4dWAKcDjwpItMIfjF70N/VqvoA8ABAXl5exP8+P9LYxn+8UsBf84uYmJbMg/+cx/J5Ni2yJ74xeWf/t5bWtvBWYTlrCsp5b1clze0ehiXEsWxGJl89Zzrnz8mylR4j1NEx+Yjf9YesXkNeVZf39JyIfB14Rn1z5NaJiBfIxNe5Twp4aQ5QEmKtrlJVntlQzD0v76C2uZ2vnjONby+faZetu8DjVTYdquGtAl+wbyupA2BiWjKfOy2HC+Zkceb0DBuGMYbQh2ueAy4A3haRWUASUAm8ADwuIvfjO/E6E1gX4me5Znd5A3c9t4W1e6s5dXIa93z6JOZmx94NHgbCdyPv0N+nrqWdf+ysYE1BOe8UVlDV2EacwGm56XxvxRwumJPFrHEjonIVz2hmN/J2Xqgh/zDwsIhsBdqA6/xd/TYReRLYjm9q5c1DcWZNS7uH3721m9+/s4fkxHh+/pmTuDpvkl2+PggCT5quKSgnf/8ROrxKWkoi580ay/lzsjh31tiYWR7CmIEKKeRVtQ24tofn7gHuCeX93fTurgrufm4r+6ua+PSiidx5+VwyR9i86f6TPo+2trR7+HBvFW8VlPNWYflxJ00vmJPFosnpxNsv2ahh8+SdZwPKxyivb+FnL+3ghU0lTM1M5bGvnMGyGZlulzWknehQvLim2RfqBeW8v6eSlnYvwxPjOMtOmhoTFhbyfl6v8ti6g9z7agGt7V5uXT6Tr5073U7ehejYIfIOj5cNB2tY4w/2wrJ6ACaNSebqvEmcPyeLJdPspGns6Jwnb728UyzkgW0ltdz57FY2FtWwdHoGP7tqAdPG2p18wkHwXTT2zIZDrCko5x87K6hr6SAhTjh9yhh+cJnvpOn0sXbSNBZFw/9yVUXVN+Tk9X/t9f/SUgWPKl5VvF7F41W8/ud9Xyter+81qUnxZI0aHvb6YjrkG1s7+OXrO3nkg/2kJSfyy6tP5qpTJlrYhFGcCLXN7XznyU1kjkji4vnjuWBOFmfNzIz5dX2Mb/sA+P7TW7j7ua0AXfufyNHZN12PdX3f+Q7BXtv5jBzzvY9XffPyveoLYVXfLH3VzseOD+3O1+MP6MDnwmXlwmx++8VTw/eGfjEb8q9tK+XfX9hGSW0LX1g8me+vmGNLyDrg2iW5TEhLZun0DE6aONpmJplucsek8N2LZlHd1NYtMDuHb7Tre//faMDX3Z8j8LmA13f/HuLE98tFxPfLQ/D9Iojr+lqO+d7/yyjgsa6fD/L6zm3c97wQH+d7Lk6E+DghLk6IE4iXzq99r8lJTwnPP+oxYi7kK+pb+dELW3l5Sylzxo/k/31xEafljnG7rKg1e/xIZo+P0puUm5DFxQnfvHCm22VEtZgJeVXluY3F/PjF7TS1erj9ktmsOmcaibaGiTEmisVEyJfUNHPns1t4q7CCUyence/nFjIjy7pLY0z0i+qQ93qVv3x8kJ+/XIDHq/xw5TyuWzrFLqYxxsSMqA35kppmbntqEx/sqWLZjAx+/umFTM5w5sSGMcZEqqgM+ec3FnPXc1vxeJWff+Ykrjl9kk2LNMbEpKgK+dqmdu56fisvbirh1Mlp/PLqU8jNSHW7LGOMcU3UhPz7uyu57alNVNS38t2LZvH186bb3X+MMTEvKlLw/d2VfOnBj0hOiueZbyzlmxfOtIB3QWEhnHsuLFoEkybBggUwZQrccgtUVrpdnYkEa9dCbi5kZ/u2j0mTYOlS37ZjnCGRtDBQXl6e5ufn9/vnPF7lkff38aUzcklOsoWt3HL55fDyy8Gfu/deuP32wa3HRJ65c6Gg4PjHL7sM/v73wa8nWojIelXNC/ZcVAzXxMcJXzl7mttlxLz774eGBqir83Xuo0f7vl+5Em64we3qTCR45BG4+mpoa4OMDKit9XXz99/vdmXRKypC3kSG2bPhnXfcrsJEsiVL4MABt6uILTZwbYwxUcxC3hhjopiFvDHGRDELeWOMiWIW8sYYE8Us5I0xJopZyBtjTBSLqCteRaQCCGUWbSYQiRfQR2pdYLUNlNXWf5FaFwz92nJVdWywJyIq5EMlIvk9XdrrpkitC6y2gbLa+i9S64Lors2Ga4wxJopZyBtjTBSLtpB/wO0CehCpdYHVNlBWW/9Fal0QxbVF1Zi8McaY7qKtkzfGGBPAQt4YY6JYVIS8iKwQkUIR2S0i33e7nk4i8rCIlIvIVrdrOZaITBKRt0Rkh4hsE5Fvu11TJxEZLiLrRGSTv7Yfu11TIBGJF5FPROQlt2sJJCL7RWSLiGwUkf7fYs1BIpImIn8TkQL/Nnem2zUBiMhs/79X5586EbnV7bo6ici/+veBrSLyFxEZ3u/3GOpj8iISD+wELgIOAR8DX1DV7a4WBojIOUAD8CdVXeB2PYFEJBvIVtUNIjISWA9cFSH/bgKkqmqDiCQC7wHfVtW1LpcGgIh8B8gDRqnqSrfr6SQi+4E8VY24i3pE5FHgXVV9UESSgBRVrXG5rG78WVIMnKGqrt/aREQm4tv256lqs4g8Cbysqn/sz/tEQye/GNitqntVtQ14ArjS5ZoAUNV/ANVu1xGMqh5W1Q3+r+uBHcBEd6vyUZ8G/7eJ/j8R0Y2ISA5wOfCg27UMFSIyCjgHeAhAVdsiLeD9LgT2RELAB0gAkkUkAUgBSvr7BtEQ8hOBooDvDxEhYTVUiMgUYBHwkculdPEPiWwEyoHXVTVSavsV8G+A1+U6glHgNRFZLyKr3C4mwDSgAnjEP8z1oIikul1UENcAf3G7iE6qWgzcBxwEDgO1qvpaf98nGkJegjwWEV3fUCAiI4CngVtVtc7tejqpqkdVTwFygMUi4vpwl4isBMpVdb3btfRgmaqeClwK3OwfLowECcCpwO9VdRHQCETMuTMA/xDSp4Cn3K6lk4ik4xuVmApMAFJF5Nr+vk80hPwhYFLA9zkM4JAmFvnHu58GHlPVZ9yuJxj/Yf3bwAp3KwFgGfAp/9j3E8AFIvK/7pZ0lKqW+P8uB57FN5QZCQ4BhwKOxv6GL/QjyaXABlUtc7uQAMuBfapaoartwDPA0v6+STSE/MfATBGZ6v9tfA3wgss1RTz/yc2HgB2qer/b9QQSkbEikub/Ohnfxl7galGAqt6hqjmqOgXfdrZGVfvdWTlBRFL9J9DxD4VcDETErC5VLQWKRGS2/6ELAddP8B/jC0TQUI3fQWCJiKT499cL8Z0765eEsJc1yFS1Q0RuAVYD8cDDqrrN5bIAEJG/AOcBmSJyCPiRqj7kblVdlgFfBrb4x74BfqCqL7tXUpds4FH/bIc44ElVjajpihFoHPCsLwtIAB5X1VfdLambbwKP+RuxvcANLtfTRURS8M3O+6rbtQRS1Y9E5G/ABqAD+IQBLHEw5KdQGmOM6Vk0DNcYY4zpgYW8McZEMQt5Y4yJYhbyxhgTxSzkjTEmilnIG2NMFLOQN8aYKPb/AbBMo//HfL5TAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "img=D3[2,:];\n",
    "#img=abs(D3[2,:]-1);\n",
    "##############################################\n",
    "nature=1\n",
    "_,_,_,_,I0,_=pars[name[nature-1]];\n",
    "print(name[nature-1])\n",
    "##########################################################\n",
    "# Designing Input Current\n",
    "current=I0*np.ones((np.shape(imgP)[0]),dtype=\"int\")\n",
    "#current[np.where(abs(np.diff(imgN))>0.001)]=255\n",
    "###############################################################\n",
    "niter=np.shape(current)[0]\n",
    "dt = 0.01;\n",
    "ti = 0.0; tf = niter*dt;\n",
    "t=np.linspace(ti,tf,niter)\n",
    "###########################################################\n",
    "#Solving the IZH Equation\n",
    "sol=IZH_NeuModel(nature,niter,dt,current)\n",
    "vsol=sol[1:,0]\n",
    "Spikes=np.zeros(np.shape(vsol),dtype=int)\n",
    "Spikes[np.where(vsol>0.0)]=1\n",
    "\n",
    "plt.plot(t,current,\"r-\",label=\"current\")\n",
    "plt.plot(t,vsol,label=\"Mem Pot\")\n",
    "plt.plot(t,Spikes*(-70.0),\"bo\",ms=1,label=\"Spikes\")"
   ]
  },
  {
   "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": []
  },
  {
   "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": []
  },
  {
   "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": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Older Investigation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Encoding "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "id": "j-N_Vgm8tldC"
   },
   "outputs": [],
   "source": [
    "def lenc(fw,r):\n",
    "    return [r]*(fw) #[length]*(width) of the pulse\n",
    "\n",
    "def wenc(fl,r):\n",
    "    return [fl]*(r) #[length]*(width) of the pulse\n",
    "\n",
    "from functools import reduce\n",
    "from math import sqrt\n",
    "def factors(n):\n",
    "        step = 2 if n%2 else 1\n",
    "        return set(reduce(list.__add__,\n",
    "                    ([i, n//i] for i in range(1, int(sqrt(n))+1, step) if n % i == 0)))\n",
    "\n",
    "\n",
    "def arenc(r):\n",
    "        a=np.sort(list(factors(r)))\n",
    "        N=np.shape(a)[0]\n",
    "        if N==1:\n",
    "            return [r]\n",
    "        elif N==2:\n",
    "            return [a[0]]*(a[-1]) #[length]*(width) of the pulse\n",
    "        elif N==3:\n",
    "            return [a[0]]*(a[-1]) #[length]*(width) of the pulse\n",
    "        elif N==4:\n",
    "            return [a[1]]*(a[-2]) #[length]*(width) of the pulse\n",
    "        else:\n",
    "            s=np.random.randint(1,(N/2))\n",
    "            if a[s]*a[-s-1]!=r:raise Exception(\"Sorry, Area not conserved!!\") \n",
    "            return [a[s]]*(a[-s-1]) #[length]*(width) of the pulse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Q6qfl_AYOXVt"
   },
   "source": [
    "# Testing "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "WVAInsn2tr3A",
    "outputId": "c6ea0112-2703-4b98-c1f2-7af6c69e14ee"
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "img=np.diff(np.array(255*D1[2,:],dtype=int));\n",
    "##############################################\n",
    "nature=1\n",
    "_,_,_,_,I0,_=pars[name[nature-1]];\n",
    "print(name[nature-1])\n",
    "##########################################################\n",
    "# Designing Input Current\n",
    "curr=[]\n",
    "for i in range(np.shape(img)[0]):\n",
    "    if img[i]!=0:\n",
    "        curr.append(wenc(I0,img[i]))    \n",
    "    else:\n",
    "        curr.append([0.0])\n",
    "    \n",
    "current=np.concatenate(curr, axis=None) #Width ENcoding\n",
    "###############################################################\n",
    "niter=np.shape(current)[0]\n",
    "dt = 0.01;\n",
    "ti = 0.0; tf = niter*dt;\n",
    "t=np.linspace(ti,tf,niter)\n",
    "###########################################################\n",
    "#Solving the IZH Equation\n",
    "sol=IZH_NeuModel(nature,niter,dt,current)\n",
    "vsol=sol[1:,0]\n",
    "Spikes=np.zeros(np.shape(vsol),dtype=int)\n",
    "Spikes[np.where(vsol>0.0)]=1\n",
    "\n",
    "plt.plot(t,current,\"r-\",label=\"current\")\n",
    "plt.plot(t,vsol,label=\"Mem Pot\")\n",
    "plt.plot(t,Spikes*(-70.0),\"bo\",ms=1,label=\"Spikes\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 283
    },
    "id": "JitNK5Bktywp",
    "outputId": "adba1ce5-2427-4257-9e5c-fde6d35bfff3"
   },
   "outputs": [],
   "source": [
    "opt_relv(Spikes,dt,200)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "oXD4_PHSld6n"
   },
   "source": [
    "## Testing Inverse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "img=np.diff(np.array(255*abs(D1[2,:]-1),dtype=int));\n",
    "##############################################\n",
    "nature=1\n",
    "_,_,_,_,I0,_=pars[name[nature-1]];\n",
    "print(name[nature-1])\n",
    "##########################################################\n",
    "# Designing Input Current\n",
    "curr=[]\n",
    "for i in range(np.shape(img)[0]):\n",
    "    if img[i]!=0:\n",
    "        curr.append(wenc(I0,img[i]))    \n",
    "    else:\n",
    "        curr.append([0.0])\n",
    "    \n",
    "current=np.concatenate(curr, axis=None) #Width ENcoding\n",
    "###############################################################\n",
    "niter=np.shape(current)[0]\n",
    "dt = 0.01;\n",
    "ti = 0.0; tf = niter*dt;\n",
    "t=np.linspace(ti,tf,niter)\n",
    "###########################################################\n",
    "#Solving the IZH Equation\n",
    "sol=IZH_NeuModel(nature,niter,dt,current)\n",
    "vsol=sol[1:,0]\n",
    "Spikes=np.zeros(np.shape(vsol),dtype=int)\n",
    "Spikes[np.where(vsol>0.0)]=1\n",
    "\n",
    "plt.plot(t,current,\"r-\",label=\"current\")\n",
    "plt.plot(t,vsol,label=\"Mem Pot\")\n",
    "plt.plot(t,Spikes*(-70.0),\"bo\",ms=1,label=\"Spikes\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "opt_relv(Spikes,dt,200)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "qsNt8JGkUOWR"
   },
   "source": [
    "# Bluk Function for width encoding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "ubUdQT3EUS7d"
   },
   "outputs": [],
   "source": [
    "def enc(f,r):\n",
    "    return [f]*(r) #[length]*(width) of the pulse\n",
    "\n",
    "def cal_Orel(D,nature):\n",
    "    no_of_images=np.shape(D)[0]\n",
    "    rel=[]\n",
    "    _,_,_,_,I0,_=pars[name[nature-1]];\n",
    "    \n",
    "    for ig in range(no_of_images):\n",
    "        #Scalling the image from 0.0 1.0 to 0 to 255 \n",
    "        img=np.array(255*D[ig,:],dtype=int);\n",
    "\n",
    "        #Encoding \n",
    "        curr=[]\n",
    "        for i in range(np.shape(img)[0]):\n",
    "            if img[i]!=0:\n",
    "                curr.append(enc(I0,img[i]))    \n",
    "            else:\n",
    "                curr.append([0.0])\n",
    "        current=np.concatenate(curr, axis=None)\n",
    "\n",
    "        #declaring time\n",
    "        niter=np.shape(current)[0]\n",
    "        dt = 0.01;\n",
    "        #ti = 0.0; tf = niter*dt;\n",
    "        #t=np.linspace(ti,tf,niter)\n",
    "\n",
    "        sol=IZH_NeuModel(nature,niter,dt,current)\n",
    "        Spikes=np.zeros(np.shape(sol[1:,0]),dtype=int)\n",
    "        Spikes[np.where(sol[1:,0]>0.0)]=1\n",
    "\n",
    "        #plt.plot(t,current,\"r-\",label=\"current\")\n",
    "        #plt.axhline(30.0)\n",
    "        #plt.plot(t,vsol,label=\"Mem Pot\")\n",
    "        #plt.plot(t,Spikes*(-70.0),\"bo\",ms=1,label=\"Spikes\")\n",
    "\n",
    "        rel.append(opt_relv(Spikes,dt,200).Hk)\n",
    "    return rel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "D8rhzoKcwduR"
   },
   "source": [
    "# Solving IZH Neuronal model with MNIST images as current Input"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "TyE14cgft7Ag"
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "nature=1\n",
    "lMSRD0=cal_Orel(D0,nature);\n",
    "# lMSRD1=cal_msr(D1,nature);\n",
    "# lMSRD2=cal_msr(D2,nature);\n",
    "# lMSRD3=cal_msr(D3,nature);\n",
    "# lMSRD4=cal_msr(D4,nature);\n",
    "# lMSRD5=cal_msr(D5,nature);\n",
    "# lMSRD6=cal_msr(D6,nature);\n",
    "# lMSRD7=cal_msr(D7,nature);\n",
    "# lMSRD8=cal_msr(D8,nature);\n",
    "# lMSRD9=cal_msr(D9,nature);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "a8GXDOx3mQY7",
    "outputId": "44420eed-6600-4ac1-dd5f-04a7d3445ae5"
   },
   "outputs": [],
   "source": [
    "for s in [lMSRD0,lMSRD1,lMSRD2,lMSRD3,lMSRD4,lMSRD5,lMSRD6,lMSRD7,lMSRD8,lMSRD9]:\n",
    "  print(np.mean(s),np.std(s))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "U7iXAQZGiq6Q",
    "outputId": "276c750e-143b-43f4-e94e-fbd6ce75cdd7"
   },
   "outputs": [],
   "source": [
    "import csv\n",
    "a=[lMSRD0,lMSRD1,lMSRD2,lMSRD3,lMSRD4,lMSRD5,lMSRD6,lMSRD7,lMSRD8,lMSRD9]\n",
    "with open(\"lMSR.csv\", \"w\", newline=\"\") as f:\n",
    "    writer = csv.writer(f)\n",
    "    writer.writerows(a)\n",
    "\n",
    "from google.colab import files\n",
    "files.download('lMSR.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "AgOepRrCUIQf"
   },
   "outputs": [],
   "source": [
    "#About 15 mins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "background_save": true
    },
    "id": "eRhcD3iOUXx0"
   },
   "outputs": [],
   "source": [
    "nature=2\n",
    "wMSRD0=cal_msr(D0,nature);\n",
    "wMSRD1=cal_msr(D1,nature);\n",
    "wMSRD2=cal_msr(D2,nature);\n",
    "wMSRD3=cal_msr(D3,nature);\n",
    "wMSRD4=cal_msr(D4,nature);\n",
    "wMSRD5=cal_msr(D5,nature);\n",
    "wMSRD6=cal_msr(D6,nature);\n",
    "wMSRD7=cal_msr(D7,nature);\n",
    "wMSRD8=cal_msr(D8,nature);\n",
    "wMSRD9=cal_msr(D9,nature);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "background_save": true,
     "base_uri": "https://localhost:8080/"
    },
    "id": "LMSATsrNU12c",
    "outputId": "03bfd762-8385-4f76-f535-65c62710354c"
   },
   "outputs": [],
   "source": [
    "for s in [wMSRD0,wMSRD1,wMSRD2,wMSRD3,wMSRD4,wMSRD5,wMSRD6,wMSRD7,wMSRD8,wMSRD9]:\n",
    "  print(np.mean(s),np.std(s))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "background_save": true,
     "base_uri": "https://localhost:8080/"
    },
    "id": "WSCw_oUOji3e",
    "outputId": "05a71d1f-9e10-42a9-d670-14405168c635"
   },
   "outputs": [],
   "source": [
    "import csv\n",
    "a=[wMSRD0,wMSRD1,wMSRD2,wMSRD3,wMSRD4,wMSRD5,wMSRD6,wMSRD7,wMSRD8,wMSRD9]\n",
    "with open(\"wMSR.csv\", \"w\", newline=\"\") as f:\n",
    "    writer = csv.writer(f)\n",
    "    writer.writerows(a)\n",
    "\n",
    "from google.colab import files\n",
    "files.download('wMSR.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "background_save": true
    },
    "id": "oA4FK7-2UlZQ"
   },
   "outputs": [],
   "source": [
    "#About 25 mins each"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "0KskHQSvgz7Q"
   },
   "source": [
    "# Bluk Function for Area encoding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "background_save": true
    },
    "id": "i76YHO2ccrqk"
   },
   "outputs": [],
   "source": [
    "from functools import reduce\n",
    "from math import sqrt\n",
    "def factors(n):\n",
    "        step = 2 if n%2 else 1\n",
    "        return set(reduce(list.__add__,\n",
    "                    ([i, n//i] for i in range(1, int(sqrt(n))+1, step) if n % i == 0)))\n",
    "\n",
    "\n",
    "def enc(r):\n",
    "        a=np.sort(list(factors(r)))\n",
    "        N=np.shape(a)[0]\n",
    "        if N==1:\n",
    "            return [r]\n",
    "        elif N==2:\n",
    "            return [a[0]]*(a[-1]) #[length]*(width) of the pulse\n",
    "        elif N==3:\n",
    "            return [a[0]]*(a[-1]) #[length]*(width) of the pulse\n",
    "        elif N==4:\n",
    "            return [a[1]]*(a[-2]) #[length]*(width) of the pulse\n",
    "        else:\n",
    "            s=np.random.randint(1,(N/2))\n",
    "            if a[s]*a[-s-1]!=r:raise Exception(\"Sorry, Area not conserved!!\") \n",
    "            return [a[s]]*(a[-s-1]) #[length]*(width) of the pulse\n",
    "\n",
    "def cal_msr(D,nature):\n",
    "    no_of_images=np.shape(D)[0]\n",
    "    msr=[]\n",
    "    \n",
    "    for ig in range(no_of_images):\n",
    "        #Scalling the image from 0.0 1.0 to 0 to 255 \n",
    "        img=np.array(255*D[ig,:],dtype=int);\n",
    "\n",
    "        #Encoding \n",
    "        curr=[]\n",
    "        for i in range(np.shape(img)[0]):\n",
    "            if img[i]!=0:\n",
    "                curr.append(enc(img[i]))    \n",
    "            else:\n",
    "                curr.append([0.0])\n",
    "        current=np.concatenate(curr, axis=None)\n",
    "\n",
    "        #declaring time\n",
    "        niter=np.shape(current)[0]\n",
    "        dt = 0.01;\n",
    "        #ti = 0.0; tf = niter*dt;\n",
    "        #t=np.linspace(ti,tf,niter)\n",
    "\n",
    "        sol=IZH_NeuModel(nature,niter,dt,current)\n",
    "        Spikes=np.zeros(np.shape(sol[1:,0]),dtype=int)\n",
    "        Spikes[np.where(sol[1:,0]>0.0)]=1\n",
    "\n",
    "        #plt.plot(t,current,\"r-\",label=\"current\")\n",
    "        #plt.axhline(30.0)\n",
    "        #plt.plot(t,vsol,label=\"Mem Pot\")\n",
    "        #plt.plot(t,Spikes*(-70.0),\"bo\",ms=1,label=\"Spikes\")\n",
    "\n",
    "        msr.append(Main.calculate_MSR(Spikes))\n",
    "    return msr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "background_save": true
    },
    "id": "Cg_t1mBAcsmA"
   },
   "outputs": [],
   "source": [
    "nature=2\n",
    "aMSRD0=cal_msr(D0,nature);\n",
    "aMSRD1=cal_msr(D1,nature);\n",
    "aMSRD2=cal_msr(D2,nature);\n",
    "aMSRD3=cal_msr(D3,nature);\n",
    "aMSRD4=cal_msr(D4,nature);\n",
    "aMSRD5=cal_msr(D5,nature);\n",
    "aMSRD6=cal_msr(D6,nature);\n",
    "aMSRD7=cal_msr(D7,nature);\n",
    "aMSRD8=cal_msr(D8,nature);\n",
    "aMSRD9=cal_msr(D9,nature);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "background_save": true
    },
    "id": "PMnNUSmYhaAe",
    "outputId": "7ffbb31a-89dd-446e-e081-0c33c08c44a9"
   },
   "outputs": [],
   "source": [
    "for s in [aMSRD0,aMSRD1,aMSRD2,aMSRD3,aMSRD4,aMSRD5,aMSRD6,aMSRD7,aMSRD8,aMSRD9]:\n",
    "  print(np.mean(s),np.std(s))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "background_save": true
    },
    "id": "-6hXzdcsjkp0",
    "outputId": "d3edc690-b666-4b30-b878-ed6490b38af5"
   },
   "outputs": [],
   "source": [
    "import csv\n",
    "a=[aMSRD0,aMSRD1,aMSRD2,aMSRD3,aMSRD4,aMSRD5,aMSRD6,aMSRD7,aMSRD8,aMSRD9]\n",
    "with open(\"aMSR.csv\", \"w\", newline=\"\") as f:\n",
    "    writer = csv.writer(f)\n",
    "    writer.writerows(a)\n",
    "\n",
    "from google.colab import files\n",
    "files.download('aMSR.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "background_save": true
    },
    "id": "3hu_jMcZj623"
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "background_save": true
    },
    "id": "BAbUra-L17dM"
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "background_save": true
    },
    "id": "pbtjHHyA17iB"
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "background_save": true
    },
    "id": "W3oWB2TG17lz"
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "ALByISjtorQ9",
    "outputId": "7cc349b8-5a1b-4e31-b32e-d65e9d8d2baa"
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "df=pd.read_excel('1_aMSR.xlsx',header=None)\n",
    "arX=[];arNS=[];arY=[];arYerr=[];\n",
    "for index, row in df.iterrows():\n",
    "  b=row.to_numpy(dtype=float)\n",
    "  print(index+1,np.nanmedian(b),np.nanstd(b))\n",
    "  arX.append(index+1)\n",
    "  arY.append(np.nanmedian(b))\n",
    "  arYerr.append(np.nanstd(b))\n",
    "  arNS.append(b.size - np.isnan(b).sum())\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "fAfHHH9h69_R",
    "outputId": "13f41c8f-bb37-4f0f-fd9b-181c494540c8"
   },
   "outputs": [],
   "source": [
    "\n",
    "df=pd.read_excel('1_lMSR.xlsx',header=None)\n",
    "lX=[];lNS=[];lY=[];lYerr=[];\n",
    "for index, row in df.iterrows():\n",
    "  b=row.to_numpy(dtype=float)\n",
    "  print(index+1,np.nanmedian(b),np.nanstd(b))\n",
    "  lX.append(index+1)\n",
    "  lY.append(np.nanmedian(b))\n",
    "  lYerr.append(np.nanstd(b))\n",
    "  lNS.append(b.size - np.isnan(b).sum())\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "eKKxUcvS6_Ip",
    "outputId": "c5b22c8f-5723-475f-ca47-6c7b9f1a4182"
   },
   "outputs": [],
   "source": [
    "\n",
    "df=pd.read_excel('1_wMSR.xlsx',header=None)\n",
    "wX=[];wNS=[];wY=[];wYerr=[];\n",
    "for index, row in df.iterrows():\n",
    "  b=row.to_numpy(dtype=float)\n",
    "  print(index+1,np.nanmedian(b),np.nanstd(b))\n",
    "  wX.append(index+1)\n",
    "  wY.append(np.nanmedian(b))\n",
    "  wYerr.append(np.nanstd(b))\n",
    "  wNS.append(b.size - np.isnan(b).sum())\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 499
    },
    "id": "HOUu03K6rF-f",
    "outputId": "3881b3e1-20e6-432c-bdad-9f9d0174dd97"
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(30,8))\n",
    "\n",
    "plt.subplot(1,3,1)\n",
    "#plt.ylim([0.14,0.26])\n",
    "plt.title(\"Length Encoding\")\n",
    "plt.scatter(lX,lY,c=lNS)\n",
    "plt.errorbar(lX,lY,yerr=lYerr)\n",
    "plt.colorbar()\n",
    "plt.xticks(lX)\n",
    "plt.grid()\n",
    "\n",
    "plt.subplot(1,3,2)\n",
    "#plt.ylim([0.14,0.26])\n",
    "plt.title(\"Width Encoding\")\n",
    "plt.scatter(wX,wY,c=wNS)\n",
    "plt.errorbar(wX,wY,yerr=wYerr)\n",
    "plt.colorbar()\n",
    "plt.xticks(wX)\n",
    "plt.grid()\n",
    "\n",
    "plt.subplot(1,3,3)\n",
    "#plt.ylim([0.1,0.26])\n",
    "plt.title(\"Area Encoding\")\n",
    "plt.scatter(arX,arY,c=arNS)\n",
    "plt.errorbar(arX,arY,yerr=arYerr)\n",
    "plt.colorbar()\n",
    "plt.xticks(X)\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "XZsjRW09rFpZ"
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "hpYYsI35rFlp"
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "M8AHAlhB19-5"
   },
   "source": [
    "# Bluk Function for amplitude encoding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "background_save": true
    },
    "id": "zoWWG94E19_V"
   },
   "outputs": [],
   "source": [
    "def enc(r):\n",
    "    return [r]*(10) #[length]*(width) of the pulse\n",
    "\n",
    "def cal_msr(D,nature):\n",
    "    no_of_images=np.shape(D)[0]\n",
    "    msr=[]\n",
    "    \n",
    "    for ig in range(no_of_images):\n",
    "        #Scalling the image from 0.0 1.0 to 0 to 255 \n",
    "        img=np.array(255*D[ig,:],dtype=int);\n",
    "\n",
    "        #Encoding \n",
    "        curr=[]\n",
    "        for i in range(np.shape(img)[0]):\n",
    "            if img[i]!=0:\n",
    "                curr.append(enc(img[i]))    \n",
    "            else:\n",
    "                curr.append([0.0])\n",
    "        current=np.concatenate(curr, axis=None)\n",
    "\n",
    "        #declaring time\n",
    "        niter=np.shape(current)[0]\n",
    "        dt = 0.01;\n",
    "        #ti = 0.0; tf = niter*dt;\n",
    "        #t=np.linspace(ti,tf,niter)\n",
    "\n",
    "        sol=IZH_NeuModel(nature,niter,dt,current)\n",
    "        Spikes=np.zeros(np.shape(sol[1:,0]),dtype=int)\n",
    "        Spikes[np.where(sol[1:,0]>0.0)]=1\n",
    "\n",
    "        #plt.plot(t,current,\"r-\",label=\"current\")\n",
    "        #plt.axhline(30.0)\n",
    "        #plt.plot(t,vsol,label=\"Mem Pot\")\n",
    "        #plt.plot(t,Spikes*(-70.0),\"bo\",ms=1,label=\"Spikes\")\n",
    "\n",
    "        msr.append(Main.calculate_MSR(Spikes))\n",
    "    return msr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "4dAwh7yo19_X"
   },
   "source": [
    "# Solving IZH Neuronal model with MNIST images as current Input"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "background_save": true
    },
    "id": "IZDLaS0H19_X"
   },
   "outputs": [],
   "source": [
    "nature=1\n",
    "lMSRD0=cal_msr(D0,nature);\n",
    "lMSRD1=cal_msr(D1,nature);\n",
    "lMSRD2=cal_msr(D2,nature);\n",
    "lMSRD3=cal_msr(D3,nature);\n",
    "lMSRD4=cal_msr(D4,nature);\n",
    "lMSRD5=cal_msr(D5,nature);\n",
    "lMSRD6=cal_msr(D6,nature);\n",
    "lMSRD7=cal_msr(D7,nature);\n",
    "lMSRD8=cal_msr(D8,nature);\n",
    "lMSRD9=cal_msr(D9,nature);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "background_save": true
    },
    "id": "ttBJvIE919_Y",
    "outputId": "4f560a32-7999-436d-a03a-566441b73355"
   },
   "outputs": [],
   "source": [
    "for s in [lMSRD0,lMSRD1,lMSRD2,lMSRD3,lMSRD4,lMSRD5,lMSRD6,lMSRD7,lMSRD8,lMSRD9]:\n",
    "  print(np.mean(s),np.std(s))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "background_save": true
    },
    "id": "ZDbcpmaD19_b",
    "outputId": "c93d4e33-05ba-4d09-afd9-39e46774825f"
   },
   "outputs": [],
   "source": [
    "import csv\n",
    "a=[lMSRD0,lMSRD1,lMSRD2,lMSRD3,lMSRD4,lMSRD5,lMSRD6,lMSRD7,lMSRD8,lMSRD9]\n",
    "with open(\"1_lMSR.csv\", \"w\", newline=\"\") as f:\n",
    "    writer = csv.writer(f)\n",
    "    writer.writerows(a)\n",
    "\n",
    "from google.colab import files\n",
    "files.download('1_lMSR.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "background_save": true
    },
    "id": "77zb6XOA19_c"
   },
   "outputs": [],
   "source": [
    "#About 15 mins"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "UOyOK9Db19_d"
   },
   "source": [
    "# Bluk Function for width encoding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "background_save": true
    },
    "id": "dRy6dQEp19_e"
   },
   "outputs": [],
   "source": [
    "def enc(r):\n",
    "    return [10]*(r) #[length]*(width) of the pulse\n",
    "\n",
    "def cal_msr(D,nature):\n",
    "    no_of_images=np.shape(D)[0]\n",
    "    msr=[]\n",
    "    \n",
    "    for ig in range(no_of_images):\n",
    "        #Scalling the image from 0.0 1.0 to 0 to 255 \n",
    "        img=np.array(255*D[ig,:],dtype=int);\n",
    "\n",
    "        #Encoding \n",
    "        curr=[]\n",
    "        for i in range(np.shape(img)[0]):\n",
    "            if img[i]!=0:\n",
    "                curr.append(enc(img[i]))    \n",
    "            else:\n",
    "                curr.append([0.0])\n",
    "        current=np.concatenate(curr, axis=None)\n",
    "\n",
    "        #declaring time\n",
    "        niter=np.shape(current)[0]\n",
    "        dt = 0.01;\n",
    "        #ti = 0.0; tf = niter*dt;\n",
    "        #t=np.linspace(ti,tf,niter)\n",
    "\n",
    "        sol=IZH_NeuModel(nature,niter,dt,current)\n",
    "        Spikes=np.zeros(np.shape(sol[1:,0]),dtype=int)\n",
    "        Spikes[np.where(sol[1:,0]>0.0)]=1\n",
    "\n",
    "        #plt.plot(t,current,\"r-\",label=\"current\")\n",
    "        #plt.axhline(30.0)\n",
    "        #plt.plot(t,vsol,label=\"Mem Pot\")\n",
    "        #plt.plot(t,Spikes*(-70.0),\"bo\",ms=1,label=\"Spikes\")\n",
    "\n",
    "        msr.append(Main.calculate_MSR(Spikes))\n",
    "    return msr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "background_save": true
    },
    "id": "8_h3ryqg19_f"
   },
   "outputs": [],
   "source": [
    "nature=1\n",
    "wMSRD0=cal_msr(D0,nature);\n",
    "wMSRD1=cal_msr(D1,nature);\n",
    "wMSRD2=cal_msr(D2,nature);\n",
    "wMSRD3=cal_msr(D3,nature);\n",
    "wMSRD4=cal_msr(D4,nature);\n",
    "wMSRD5=cal_msr(D5,nature);\n",
    "wMSRD6=cal_msr(D6,nature);\n",
    "wMSRD7=cal_msr(D7,nature);\n",
    "wMSRD8=cal_msr(D8,nature);\n",
    "wMSRD9=cal_msr(D9,nature);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "background_save": true
    },
    "id": "FVx7mTq419_g",
    "outputId": "7552bdfe-be2f-4726-ab73-a922e2c29893"
   },
   "outputs": [],
   "source": [
    "for s in [wMSRD0,wMSRD1,wMSRD2,wMSRD3,wMSRD4,wMSRD5,wMSRD6,wMSRD7,wMSRD8,wMSRD9]:\n",
    "  print(np.mean(s),np.std(s))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "background_save": true
    },
    "id": "yPOZul3h19_h",
    "outputId": "56f30c58-4073-4af5-fd48-5e850d45d936"
   },
   "outputs": [],
   "source": [
    "import csv\n",
    "a=[wMSRD0,wMSRD1,wMSRD2,wMSRD3,wMSRD4,wMSRD5,wMSRD6,wMSRD7,wMSRD8,wMSRD9]\n",
    "with open(\"1_wMSR.csv\", \"w\", newline=\"\") as f:\n",
    "    writer = csv.writer(f)\n",
    "    writer.writerows(a)\n",
    "\n",
    "from google.colab import files\n",
    "files.download('1_wMSR.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "background_save": true
    },
    "id": "RuAbxARg19_h"
   },
   "outputs": [],
   "source": [
    "#About 25 mins each"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "xoiu9hAl19_i"
   },
   "source": [
    "# Bluk Function for Area encoding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "background_save": true
    },
    "id": "g4Wch-Bv19_i"
   },
   "outputs": [],
   "source": [
    "from functools import reduce\n",
    "from math import sqrt\n",
    "def factors(n):\n",
    "        step = 2 if n%2 else 1\n",
    "        return set(reduce(list.__add__,\n",
    "                    ([i, n//i] for i in range(1, int(sqrt(n))+1, step) if n % i == 0)))\n",
    "\n",
    "\n",
    "def enc(r):\n",
    "        a=np.sort(list(factors(r)))\n",
    "        N=np.shape(a)[0]\n",
    "        if N==1:\n",
    "            return [r]\n",
    "        elif N==2:\n",
    "            return [a[0]]*(a[-1]) #[length]*(width) of the pulse\n",
    "        elif N==3:\n",
    "            return [a[0]]*(a[-1]) #[length]*(width) of the pulse\n",
    "        elif N==4:\n",
    "            return [a[1]]*(a[-2]) #[length]*(width) of the pulse\n",
    "        else:\n",
    "            s=np.random.randint(1,(N/2))\n",
    "            if a[s]*a[-s-1]!=r:raise Exception(\"Sorry, Area not conserved!!\") \n",
    "            return [a[s]]*(a[-s-1]) #[length]*(width) of the pulse\n",
    "\n",
    "def cal_msr(D,nature):\n",
    "    no_of_images=np.shape(D)[0]\n",
    "    msr=[]\n",
    "    \n",
    "    for ig in range(no_of_images):\n",
    "        #Scalling the image from 0.0 1.0 to 0 to 255 \n",
    "        img=np.array(255*D[ig,:],dtype=int);\n",
    "\n",
    "        #Encoding \n",
    "        curr=[]\n",
    "        for i in range(np.shape(img)[0]):\n",
    "            if img[i]!=0:\n",
    "                curr.append(enc(img[i]))    \n",
    "            else:\n",
    "                curr.append([0.0])\n",
    "        current=np.concatenate(curr, axis=None)\n",
    "\n",
    "        #declaring time\n",
    "        niter=np.shape(current)[0]\n",
    "        dt = 0.01;\n",
    "        #ti = 0.0; tf = niter*dt;\n",
    "        #t=np.linspace(ti,tf,niter)\n",
    "\n",
    "        sol=IZH_NeuModel(nature,niter,dt,current)\n",
    "        Spikes=np.zeros(np.shape(sol[1:,0]),dtype=int)\n",
    "        Spikes[np.where(sol[1:,0]>0.0)]=1\n",
    "\n",
    "        #plt.plot(t,current,\"r-\",label=\"current\")\n",
    "        #plt.axhline(30.0)\n",
    "        #plt.plot(t,vsol,label=\"Mem Pot\")\n",
    "        #plt.plot(t,Spikes*(-70.0),\"bo\",ms=1,label=\"Spikes\")\n",
    "\n",
    "        msr.append(Main.calculate_MSR(Spikes))\n",
    "    return msr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "tS4BzpMN19_j"
   },
   "outputs": [],
   "source": [
    "nature=1\n",
    "aMSRD0=cal_msr(D0,nature);\n",
    "aMSRD1=cal_msr(D1,nature);\n",
    "aMSRD2=cal_msr(D2,nature);\n",
    "aMSRD3=cal_msr(D3,nature);\n",
    "aMSRD4=cal_msr(D4,nature);\n",
    "aMSRD5=cal_msr(D5,nature);\n",
    "aMSRD6=cal_msr(D6,nature);\n",
    "aMSRD7=cal_msr(D7,nature);\n",
    "aMSRD8=cal_msr(D8,nature);\n",
    "aMSRD9=cal_msr(D9,nature);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "X0wy9Ziy19_k"
   },
   "outputs": [],
   "source": [
    "for s in [aMSRD0,aMSRD1,aMSRD2,aMSRD3,aMSRD4,aMSRD5,aMSRD6,aMSRD7,aMSRD8,aMSRD9]:\n",
    "  print(np.mean(s),np.std(s))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "M58VAW9W19_l"
   },
   "outputs": [],
   "source": [
    "import csv\n",
    "a=[aMSRD0,aMSRD1,aMSRD2,aMSRD3,aMSRD4,aMSRD5,aMSRD6,aMSRD7,aMSRD8,aMSRD9]\n",
    "with open(\"1_aMSR.csv\", \"w\", newline=\"\") as f:\n",
    "    writer = csv.writer(f)\n",
    "    writer.writerows(a)\n",
    "\n",
    "from google.colab import files\n",
    "files.download('1_aMSR.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "0kxaIhfC19_l"
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [
    "T-6BRn8yu4n9",
    "wDF_HcqkuyQn",
    "rIB1Lpb1vb8y",
    "fzzabiedvsZe",
    "uBO60vHwwVoT",
    "eNwZ1nQvwJU5"
   ],
   "name": "IZH_Soumik.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python [conda env:work1]",
   "language": "python",
   "name": "conda-env-work1-py"
  },
  "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.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}