train_energy_model.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "import os\n",
    "import glob\n",
    "%matplotlib inline\n",
    "from scipy.stats import zscore\n",
    "import imp\n",
    "# modules\n",
    "import decoders, utils, tuning, learning, mainfigs, suppfigs\n",
    "\n",
    "### WHERE YOU DOWNLOADED THE FIGSHARE\n",
    "dataroot = 'Z:/releases/stringer-et-al-2019a/fs125'\n",
    "\n",
    "# file list\n",
    "db = np.load(os.path.join(dataroot, 'database.npy'), allow_pickle=True)\n",
    "fs = []\n",
    "for di in db:\n",
    "    mname = di['mouse_name']\n",
    "    datexp = di['date']\n",
    "    blk = di['block']\n",
    "    stype = di['expt']\n",
    "    \n",
    "    fname = '%s_%s_%s_%s.npy'%(stype, mname, datexp, blk)\n",
    "    fs.append(os.path.join(dataroot, fname))\n",
    "\n",
    "### WHERE YOU WANT TO SAVE THE OUTPUTS OF THE ANALYSIS SCRIPTS AND THE FIGURES (if save_figure=True)\n",
    "saveroot = 'D:\\DATA\\stringer-etal2019'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = fs[0]\n",
    "dat = np.load(f, allow_pickle=True).item()\n",
    "sresp, istim, itrain, itest = utils.compile_resp(dat, npc=32)\n",
    "\n",
    "x = sresp[:, itrain]\n",
    "y = istim[itrain]\n",
    "NN, NT = np.shape(sresp)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "th_range = np.arange(0, 2*np.pi, 2*np.pi/360)\n",
    "rcos = np.cos(th_range)\n",
    "rsin = np.sin(th_range)\n",
    "\n",
    "def log_prob(a):\n",
    "    atx = a[0].T @ x\n",
    "    btx = a[1].T @ x\n",
    "    \n",
    "    Z = atx * rcos[:, np.newaxis] + btx * rsin[:, np.newaxis]\n",
    "    Zmax = np.max(Z, axis=0)\n",
    "    \n",
    "    Z = np.exp(Z-Zmax)\n",
    "    Zsum = np.sum(Z,axis=0)\n",
    "\n",
    "    logL = np.mean(atx * np.cos(y) + btx * np.sin(y) - np.log(Zsum) - Zmax)\n",
    "\n",
    "    Zcos = rcos.T @ Z / Zsum\n",
    "    Zsin = rsin.T @ Z / Zsum\n",
    "    \n",
    "    da = (x @ (np.cos(y) - Zcos))#/x.shape[1]\n",
    "    db = (x @ (np.sin(y) - Zsin))#/x.shape[1]\n",
    "    \n",
    "    return logL, np.stack((da,db))\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-5.886107559374568\n",
      "-0.4387410879499554\n",
      "-0.16913275074296877\n",
      "-0.15764865463476466\n",
      "-0.15310199531616198\n",
      "-0.14904430211714798\n"
     ]
    }
   ],
   "source": [
    "a = np.random.randn(2, NN)/1e7 # initializdr with very small values\n",
    "\n",
    "eps0 = 0.05 # learning rate\n",
    "niter = 501\n",
    "lam = .0 # regularization parameter, makes fitting unstable\n",
    "\n",
    "logL = np.zeros(niter,)\n",
    "pa = np.zeros(a.shape)\n",
    "\n",
    "for it in range(niter):    \n",
    "    logL[it], da = log_prob(a)\n",
    "    \n",
    "    pa = .95 * pa + .05 * (da - lam * a)\n",
    "    \n",
    "    if it<20:\n",
    "        eps = eps0/(20-it)\n",
    "    a += eps * pa    \n",
    "    \n",
    "    if it%100==0:\n",
    "        print(logL[it])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x27b8e382780>]"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAG55JREFUeJzt3Xl0nHd97/H3d0aWbFm25U22E0txYjuxTeqYoISEEMjiNCFA2LdeGtrm1C0tve257aWE3EuhnJ7TBmgovZRbp6W0BQqX03ADoYXEISGESxI7ux3HSxw7sWVbkhct1j7zvX/MM9JIGS3WPNIzz6PP6xwdzbPMM9/fWP48v/k9y5i7IyIiyZGKugAREQmXgl1EJGEU7CIiCaNgFxFJGAW7iEjCKNhFRBJGwS4ikjAKdhGRhFGwi4gkTEUUL3rTTTf5j3/84yheWkQkzmwiK4XSYzezm8xsj5ntN7NPjbd+a2trGC8rIiJFlBzsZpYGvgq8DdgAfMTMNpS6XRERmZwweuyXA/vd/YC79wHfAd4VwnZFRGQSwgj2c4FXC6YPB/NERCQCYQR7scH819wL2My2mNkOM9vR0tISwsuKiEgxYQT7YaC+YHol0DRyJXff6u6N7t64dOnSEF5WRESKCSPYtwNrzex8M6sEPgz8IITtiojIJJR8Hru7D5jZJ4CfAGng6+6+q+TKRERkUkK5QMnd/wP4jzC2JeXB3enuz9DZM8CZvgy9Axl6+7P0ZbL0DWTpHcgEv3PTWXcyWci4k806mayTdR+cn/s9+tcwjnbVhY2ywMwwg5QZqeC3mZE2SKVyj1MFy3PThevmlxUuD+alJrZ+OjX68lRq+PNTKSM9YlnaDEsx9NgItjn0GiKTEcmVpzL9BjJZjrb1cPhUN80dPZzo7KO1szf46eNUV18uxHsH6OjN/R4jh2Ua5HcWxUI/v6PI75AGl42yo8jtiIbvfNKpIs8fsSxdsKMrtmMadVkqeP0iO7uRy4a1LdhZDj6vYOdZuGzYjrXgfRlrWb59I5fl36sk7XwV7AnT1tXP3uYO9jd3su94J/uaOzh0ooum090MjEjqipSxuKaSJTVVLKyuZPn82cytqqAm/zO7IphOU1WRpqoiRWVFisp0iqpZaSrTQ9Pp9NB/9HT+P8/gf5Sh+cX+I4z2hepj7Vey7riDe/BpwB0PPhnkfnLbzfrQPB98nH9+wfIsI9YZexuZ7Gu3l806DsM+rWRH1JR7Xv4TzWjLCucH62WHf/oZbZkH70Xh/PGWZbO5Hf/IZYVtzARtH2xbUH+m4H0cb1kSjbqDGmPn+8UPXMKVqxdPaV0K9pg7fKqLHQdP8cTBk+w4eJK9xzsHl82ZlWZ13Vwuqa/lHRtXUL+omvqF1SybX8WSmioWzJlFKlUOPY5yqEGmWrEdWn56tGWFO9Jiy4Z2okM715E7oeHrjdgJTcPON79TzT9eOHfWlL/XCvaYcXeeefU09z13lIdebOZA6xkA5lVVcOl5C3nnxnO4+NwFrKmr4dzaOWUS3CK5YYuKtP4ep4OCPSZeaunk3qePcO+zTRw60UVlRYorL1jMR684jzdesIh1y+eTVoiLCAr2subuPLy3hb/Zto9nXj1NyuBNq5fwiWvXcOPFy5k/e+o/0olI/CjYy9TjB07wxfv3sP3gKeoXzeF/vmMD79y4grr5s6MuTUTKnIK9zDx/uI0v3L+HR/a2UDevis+/+2I+1FhPZYW+7EpEJkbBXib6M1n+4ke7+cb/O0ht9Sw+ffM6fv2KVcypTEddmojEjIK9DDS39/D7336K7QdP8ZtXreK/3XAh8zR+LiKTpGCP2JOHTvLxbz5FR88AX/nI67nlknOiLklEYk7BHhF355uPHeLP73uBc2rn8C+3Xc665fOjLktEEkDBHoFM1vnzH+7in395iOvW1XHXhzaxYI6GXkQkHAr2CPzZD3byzcde4bevPp/b37ZeV4eKSKgU7NPsezte5ZuPvcKWt1zAp29eH3U5IpJAOjl6Gv1sbwuf/v7zvGn1Yj5540VRlyMiCaVgnyZ7j3fw8W8+yZq6eXzto2+gIq23XkSmRknpYmYfMLNdZpY1s8awikqa9p5+fvdfn2RuVQXf+M3LdKBURKZUqd3GncB7gUdCqCWR3J3//r1nOXSyi6/+2qUs071eRGSKlXTw1N13Q/l+PVQ5uPvnB/jJruP8j7ev5/LzF0VdjojMABronUIHWjr5qx/v4W0XL+e2N58fdTkiMkOM22M3s23A8iKL7nD3eyf6Qma2BdgC0NDQMOEC4+yubfuoTKf4/Lsv1qcaEZk24wa7u28O44XcfSuwFaCxsTGZ32xbYOeRNn74bBO/d81qltRURV2OiMwgGoqZAgOZLLff8zxLair5nbesjrocEZlhSj3d8T1mdhi4EviRmf0knLLi7V8fO8TzR9r47C2vY0G1Tm0UkelV6lkx3we+H1ItidDZO8Df/nQ/V61ZzNt/ZUXU5YjIDKShmJD9w88PcPJMH5+8cZ0OmIpIJBTsITp5po9/+PnL3PS65VxSXxt1OSIyQynYQ/R3D+2nq2+AP7nxwqhLEZEZTMEekqNt3fzLY4d476UrWVM3L+pyRGQGU7CH5O5HXsbd+aPNa6MuRURmOAV7CPozWe595gg3bFjGyoXVUZcjIjOcgj0Ej+xt4cSZPt7z+pVRlyIiomAPwz1PH2Fh9SzeeuHSqEsREVGwl6qlo5f7dx3j3a8/l8oKvZ0iEj0lUYm+88Qr9Gecj15xXtSliIgACvaSZLLOt594hTevWcLqpTVRlyMiAijYS7L94EmOtvXwwcvqoy5FRGSQgr0E9z3XxOxZKa5fVxd1KSIigxTskzSQyfKfzx/j+nXLmFtV0k0yRURCpWCfpEf3t3LiTB/v2Khb84pIeVGwT9J3t7/KormVXLdewzAiUl4U7JPQ0tHLAy8c532XnktVRTrqckREhin1q/G+YGYvmtlzZvZ9M5sRNyG/56nDDGSdD+lsGBEpQ6X22B8ALnb3jcBe4PbSSypv7s53t7/KZasW6va8IlKWSgp2d7/f3QeCyceAxN8Fa+/xTg60nuG9lya+qSISU2GOsf8W8J8hbq8sPbK3BYBrLtINv0SkPI17AraZbQOWF1l0h7vfG6xzBzAAfGuM7WwBtgA0NDRMqthy8LO9LVy4rIYVC+ZEXYqISFHjBru7bx5ruZl9DHgHcL27+xjb2QpsBWhsbBx1vXLW1TfAEy+f5NYrdcMvESlfJV0yaWY3AX8KvNXdu8IpqXw9duAEfZksb9F910WkjJU6xv6/gHnAA2b2jJn97xBqKlvbdjcztzLNGy9YFHUpIiKjKqnH7u5rwiqk3Lk7P93dzNVrl+qiJBEpa7rydIJ2NbVzrL2H63ULAREpcwr2CXp4TzNmcK1u0SsiZU7BPkGPv3ySi5bNY0lNVdSliIiMScE+Af2ZLE8eOsXl5+ugqYiUPwX7BOxqaqerL6NgF5FYULBPwPaXTwJw+SoFu4iUPwX7BPzywAlWLa6mbv7sqEsRERmXgn0c3X0ZfrG/lWsu0tkwIhIPCvZxPLq/ld6BLJvXL4u6FBGRCVGwj+PB3ceZV1WhA6ciEhsK9jFks8623c285aKlVFborRKReFBajeG5I220dvayWbcREJEYUbCP4aEXm0kZXHOhgl1E4kPBPoadR9pYU1fDwrmVUZciIjJhCvYx7GpqZ8OK+VGXISJyVhTsozjR2cux9h42nKNgF5F4UbCPYvfRDgBed86CiCsRETk7CvZR7GpqA2C9hmJEJGZKCnYz+7yZPRd83+n9ZnZOWIVF7YWj7axYMJtFOnAqIjFTao/9C+6+0d03AfcBnwmhpsi5O9tfPsklK2ujLkVE5KyVFOzu3l4wORfw0sopDy+1dNLU1sPVFy6JuhQRkbNWUeoGzOwvgFuBNuDaMdbbAmwBaGhoKPVlp9TDe1oAeMvapRFXIiJy9sbtsZvZNjPbWeTnXQDufoe71wPfAj4x2nbcfau7N7p749Kl5R2Yj+xr5YKlc6lfVB11KSIiZ23cHru7b57gtr4N/Aj4s5Iqilh/JssTL5/gw5eV96cKEZHRlHpWzNqCyVuAF0srJ3p7j3fQ05/l9Q06cCoi8VTqGPtfmtlFQBY4BPxu6SVF67nDufPXdUaMiMRVScHu7u8Lq5By8dzh08yfXcF5izW+LiLxpCtPR3jucBsbV9ZiZlGXIiIyKQr2Aj39GfYc62DjSt0fRkTiS8Fe4NF9rQxknU31Gl8XkfhSsBf4+0de4tzaOVy7Tt+YJCLxpWAPPHnoJNsPnuK3rz6fWWm9LSISX0qwwKP7TmAG72+sj7oUEZGSKNgD+5o7qF9YTU1VybfPERGJlII9sL+5k7V1NVGXISJSMgU7MJDJcqDlDGuWKdhFJP4U7MChk130ZbKsrZsXdSkiIiVTsAP7jncCaChGRBJBwQ7sb+4AYLWCXUQSQMEO7GpqZ+XCOTojRkQSYcYHeybr/PLACa64YHHUpYiIhGLGB/sLTe2c7urnzWv0xdUikgwzPtgf3d8KwFUKdhFJiBkf7L/Y38q65fNYOq8q6lJEREIRSrCb2Z+YmZtZrLq97s7zR9q49LyFUZciIhKakoPdzOqBG4BXSi9nep0400dbdz+rl+o0RxFJjjB67HcBnwQ8hG1NqwMtZwBYvXRuxJWIiISnpGA3s1uAI+7+bEj1TKuXWnJXnKrHLiJJMu4VOWa2DVheZNEdwKeBX53IC5nZFmALQENDw1mUOHW2v3yShdWzOLd2TtSliIiEZtxgd/fNxeab2a8A5wPPmhnASuApM7vc3Y8V2c5WYCtAY2Nj5MM23378Fe55+gi3XHIOqZRFXY6ISGgmfQ29uz8PDH45qJkdBBrdvTWEuqbc3z28H4APNK6MuBIRkXAl8jz21s5ebvjrn/Eb//RE0eWvnOji8KluPvvODVy9duk0VyciMrVCu+uVu68Ka1ul2nusg33Nnexr7qQ/k33Nl1PvbGoDoHHVoijKExGZUonssZ/u7h98fPR0z2uWHznVDUD9ouppq0lEZLokM9i7hoL90Mkzr1l+5HQ3NVUVzJ+t2/SKSPIkMthPdfUNPn7lZNdrlh853c25tXMIzuYREUmURAZ7W3c/lRUp0ikrOhTTdLqbc2pnR1CZiMjUS2Swn+7qY1F1JTVVFXT2Dgxb5u68cqJL4+sikliJG2R+9WQXrZ191FbPIt1jtPf0D1t++FQ3Hb0DrF8xP6IKRUSmVqKCvaOnn6vvfAiAN56fO5Wxs2d4j31XUzsAGxTsIpJQiQr2/c2dg4+XzKsi6/6aoZgXjraTMrho+bzpLk9EZFokaoy9MNh/7fKGomPsLzS1s3ppDbNnpae7PBGRaZGsYA9uw3v3rY1ctWYJ82bPomPEUMzuo+1sOEfDMCKSXIkK9oOtZ1hTV8MNG5YBUDO7Yliwt3X1c+R0tw6cikiiJSrYj7X1sGLB0Pnp86oq6OwdOivm1VO5i5VWLdY3JolIciUr2Nt7WD5/KNhrqiro6c/Sn8my51gHe451ALB8gS5OEpHkSsxZMQOZLC0dvcN77MG9YNq7+7nxy48Mzi8MfxGRpElMj72ls5esw7KCYF9QPQsYOncdIJ0yls6rmvb6RESmS2KC/Whb7p4whb3xZfNyj7/35OHBeUtrqkjrq/BEJMESE+w/3nmMdMp43TkLBufle+8/fLZpcN7imsppr01EZDqVFOxm9lkzO2JmzwQ/N4dV2Nn6wTNNbF5fN+zAaLGx9D+4bu10liUiMu3COHh6l7t/MYTtTFpHTz/H2nu4tf68YfPnVg01b9fnbqSnP8PiGo2vi0iyJeKsmAMtuW9JWr20ZtR15lZVDAt6EZGkCiPpPmFmtwI7gD9291MhbPOsvBTcSqBYsN/5vo1UViTmUIKIyLjGDXYz2wYsL7LoDuBrwOcBD35/CfitUbazBdgC0NDQMMlyizscfDl1Q5Evz/jgZfWhvpaISLkbN9jdffNENmRmdwP3jbGdrcBWgMbGRp9ogRNxpm+AyoqUeuYiIpR+VsyKgsn3ADtLK2dyuvsyVFfqNrwiIlD6GPudZraJ3FDMQeB3Sq5oErr6MlTr/uoiIkCJwe7uvx5WIaXo7sswRz12EREgIVeenukb0KmMIiKBRAR7V1+GORqKEREBEhLsOngqIjIkEcHe1TdAdaWGYkREICHBroOnIiJDEhHsXf0aihERyUtGsKvHLiIyKPbBPpDJ0jeQpXqWxthFRCABwd7VnwFgbpV67CIikIBg7+7LBbuGYkREcmIf7F1BsOvgqYhITgKCfQCAORpjFxEBEhDs3eqxi4gME/tg11CMiMhwiQl2HTwVEclJQLDnxth1rxgRkZwEBHtwHrt67CIiQAjBbmZ/YGZ7zGyXmd0ZRlGj+dL9e3j7V34+bJ7OYxcRGa6k8QszuxZ4F7DR3XvNrC6csopr6+7nyOnuYfOGDp5qKEZEBErvsX8c+Et37wVw9+bSSxpdOmVksj5sXlf/AJUVKdIpm8qXFhGJjVKD/ULgajN73Mx+ZmaXhVHUaCqKBLu+PUlEZLhxxy/MbBuwvMiiO4LnLwSuAC4D/o+ZXeDuPnJlM9sCbAFoaGiYVLGplDEwssfel6Fa33cqIjJo3GB3982jLTOzjwP3BEH+hJllgSVAS5HtbAW2AjQ2Nr4m+CdUbMrIFumx68CpiMiQUodi/i9wHYCZXQhUAq2lFjWatOV67IUfCNp7+qmp0oFTEZG8UhPx68DXzWwn0Ad8rNgwTFjSqdx+KOuQDo6VtnT0snJh9VS9pIhI7JQU7O7eB3w0pFrGVRGkeSbrg2fBNHf08obzFk5XCSIiZS9WV56mbCjYAfoGspw800fdvNlRliUiUlZiFewVQS99IJsFoKWzF4Bl86siq0lEpNzEKtjzwy9BrtPc3gNAnYJdRGRQLIM932Nv7sj12DUUIyIyJJbBnglOvGnr6gdg4dzKyGoSESk3sQr2/Bh7/uDp6e4+AGrnzIqsJhGRchOrYE/lh2IyQbB39VORMt0rRkSkQKyCPd9jz3q+x95PbfUszHRnRxGRvFgF+9DB06Ex9gUahhERGSaWwV44xl5brQOnIiKFYhXsrzl42tWvA6ciIiPEKthH3lLgdFc/C6oV7CIihWIV7PmbgOXH2Nu7+5k/W8EuIlIoVsGev21vJrgn+5m+Ad2LXURkhHgFe8FQTO9AlqxDdZXOYRcRKRSr7m7+rJi7f36A8xblvlxjbmWsmiAiMuVilYr5MfYHXjg+OE9XnYqIDFdSsJvZd4GLgsla4LS7byq5qlGkilxhOldj7CIiw5T61Xgfyj82sy8BbSVXNIb8eeyF1GMXERkulO6u5W7W8kHgujC2N5p0kWBXj11EZLiwzoq5Gjju7vtC2l5RRYNdB09FRIYZNxXNbBuwvMiiO9z93uDxR4B/G2c7W4AtAA0NDWdZZk6xoZi5Ot1RRGSYcYPd3TePtdzMKoD3Am8YZztbga0AjY2NfhY1DkoVHWNXj11EpFAYQzGbgRfd/XAI2xqTeuwiIuMLI9g/zDjDMGEpNsY+u0LBLiJSqORxDHf/jRDqmJDCYL+kvpZNKxcUHZ4REZnJYjVAXRjs16+r479evzbCakREylOsbgJWkRoqt7IiVqWLiEybWKVjuuCWArPSsSpdRGTaxCod0+mhYK9Ma2xdRKSYeAW7euwiIuOKVToWHjxVsIuIFBerdCy8QEkHT0VEiotVOqbUYxcRGVds07GyQgdPRUSKiW2wq8cuIlJcbNOxUsEuIlJU7NIxf8bjLB08FREpKnbpuLC6ElCPXURkNLFLx9rqWQCkTAdPRUSKiV2wn1s7B4Du/kzElYiIlKdY3bYX4EsfvIR/+sVBNtXXRl2KiEhZil2w182bzZ/etC7qMkREylbshmJERGRsJQW7mW0ys8fM7Bkz22Fml4dVmIiITE6pPfY7gc+5+ybgM8G0iIhEqNRgd2B+8HgB0FTi9kREpETm7pN/stl64CeAkdtJvMndD42y7hZgC0BDQ8MbDh0qupqIiIxuQhfwjBvsZrYNWF5k0R3A9cDP3P3fzeyDwBZ33zzeizY2NvqOHTsmUp+IiAwJJ9jHfLJZG1Dr7m5mBrS5+/zxnqdgFxGZlAkFe6lj7E3AW4PH1wH7StyeiIiUqNQe+5uBvyF3oVMP8Hvu/uQEntcCTHaQfQnQOsnnxpXaPDOozTNDKW1udfebxluppGCPgpntcPfGqOuYTmrzzKA2zwzT0WZdeSoikjAKdhGRhIljsG+NuoAIqM0zg9o8M0x5m2M3xi4iImOLY49dRETGEKtgN7ObzGyPme03s09FXU9YzOzrZtZsZjsL5i0yswfMbF/we2Ew38zsK8F78JyZXRpd5ZNjZvVm9pCZ7TazXWb2h8H8xLYZwMxmm9kTZvZs0O7PBfPPN7PHg3Z/18wqg/lVwfT+YPmqKOufLDNLm9nTZnZfMJ3o9gKY2UEzez5/59tg3rT9fccm2M0sDXwVeBuwAfiImW2ItqrQfAMYeW7qp4AH3X0t8GAwDbn2rw1+tgBfm6YawzQA/LG7rweuAH4/+LdMcpsBeoHr3P0SYBNwk5ldAfwVcFfQ7lPAbcH6twGn3H0NcFewXhz9IbC7YDrp7c271t03FZzaOH1/3+4eix/gSuAnBdO3A7dHXVeI7VsF7CyY3gOsCB6vAPYEj/8e+Eix9eL6A9wL3DDD2lwNPAW8kdzFKhXB/MG/c3I32LsyeFwRrGdR136W7VwZhNh1wH3kLolPbHsL2n0QWDJi3rT9fcemxw6cC7xaMH04mJdUy9z9KEDwuy6Yn6j3Ifi4/XrgcWZAm4NhiWeAZuAB4CXgtLsPBKsUtm2w3cHyNmDx9FZcsi8DnwSywfRikt3ePAfuN7MngzvbwjT+fcfpO0+L3fxmJp7Sk5j3wcxqgH8H/sjd23P3kSu+apF5sWyzu2eATWZWC3wfWF9steB3rNttZu8Amt39STO7Jj+7yKqJaO8IV7l7k5nVAQ+Y2YtjrBt6u+PUYz8M1BdMryTZX+xx3MxWAAS/m4P5iXgfzGwWuVD/lrvfE8xOdJsLuftp4GFyxxhqzSzfySps22C7g+ULgJPTW2lJrgJuMbODwHfIDcd8meS2d5C7NwW/m8ntwC9nGv++4xTs24G1wRH1SuDDwA8irmkq/QD4WPD4Y+TGofPzbw2OpF9B7lbJR6MocLIs1zX/R2C3u/91waLEthnAzJYGPXXMbA6wmdxBxYeA9werjWx3/v14P/BTDwZh48Ddb3f3le6+itz/15+6+38hoe3NM7O5ZjYv/xj4VWAn0/n3HfVBhrM8IHEzsJfcuOQdUdcTYrv+DTgK9JPbe99GbmzxQXK3Qn4QWBSsa+TODnoJeB5ojLr+SbT3zeQ+aj4HPBP83JzkNgft2Ag8HbR7J/CZYP4FwBPAfuB7QFUwf3YwvT9YfkHUbSih7dcA982E9gbtezb42ZXPqun8+9aVpyIiCROnoRgREZkABbuISMIo2EVEEkbBLiKSMAp2EZGEUbCLiCSMgl1EJGEU7CIiCfP/ATZG2Fim+TW+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(logL[:it])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.149023215611123"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dx = a[0].T @ sresp[:, itest]\n",
    "dy = a[1].T @ sresp[:, itest]\n",
    "\n",
    "apred = np.angle(dx + 1j * dy)\n",
    "apred[apred<0] = apred[apred<0] + 2*np.pi\n",
    "\n",
    "nangle = 2*np.pi\n",
    "error = istim[itest] - apred\n",
    "error = np.remainder(error, nangle)\n",
    "error[error > nangle/2] = error[error > nangle/2] - nangle\n",
    "\n",
    "np.median(np.abs(error)) * 180/np.pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'plt' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-1-807a70c7e323>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0merror\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m: name 'plt' is not defined"
     ]
    }
   ],
   "source": [
    "plt.plot(error)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}