Spiked-Covariance-Model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spiked Covariance Model\n",
    "\n",
    "This notboke\n",
    "\n",
    "### Demostrate of Edge Statistics of this Model\n",
    "(It is partially complete and could use a little work )\n",
    "\n",
    "\n",
    "#### Bulk Edge conventions:  Johnstone(statistics) vs Boucaud(physics)\n",
    "\n",
    "Depending on your point of view, stats or phys, we can express our uncertainty in $\\lambda^{+}$ in terms of $N$ or $M$. \n",
    "\n",
    "\n",
    "- Create an ensemble of Gaussian random matrices, each  with a rank 1 perturbation\n",
    "\n",
    "$$\\mathbf{W}^{rand}\\;(M\\times N)+2\\mathbf{\\Delta}_{0}\\; (rank 1) $$\n",
    "\n",
    " and form the Empirical Spectral Density (ESD), $\\rho_{emp}(\\lambda)$, of the $(M\\times M)$ correlation matrix\n",
    " \n",
    "$$\\mathbf{X}=\\dfrac{1}{N}\\mathbf{W}^{T}\\mathbf{W}^{T}$$\n",
    "\n",
    "\n",
    "- Demonstrate the fluctuations in both the MP bulk edge $\\lambda^{+}$\n",
    "\n",
    "- and the (Gaussian) shape of the fluctuations of the spikes $\\lambda_{k} > $\\lambda^{+}$\n",
    "\n",
    "\n",
    "#### References\n",
    "\n",
    "- Spiked Covariance Model of Johnstone\n",
    "\n",
    "- Applications to Finance, Bouchaud and co-workers\n",
    "\n",
    "- Sornette's model of Self Organization\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-29T16:54:12.450642Z",
     "start_time": "2018-08-29T16:54:09.932439Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/charlesmartin14/anaconda3/lib/python3.6/importlib/_bootstrap.py:219: RuntimeWarning: compiletime version 3.5 of module 'tensorflow.python.framework.fast_tensor_util' does not match runtime version 3.6\n",
      "  return f(*args, **kwds)\n",
      "/Users/charlesmartin14/.local/lib/python3.6/site-packages/h5py/__init__.py:34: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
      "  from ._conv import register_converters as _register_converters\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.19.1\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import tensorflow as tf\n",
    "import pickle, time\n",
    "from copy import deepcopy\n",
    "from shutil import copy\n",
    "from tqdm import tqdm_notebook as tqdm\n",
    "import pandas as pd\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "import import_ipynb\n",
    "\n",
    "import sklearn\n",
    "from sklearn.decomposition import TruncatedSVD\n",
    "from sklearn.random_projection import sparse_random_matrix\n",
    "print(sklearn.__version__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-29T16:54:12.901495Z",
     "start_time": "2018-08-29T16:54:12.453155Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "importing notebook from RMT_Util.ipynb\n",
      "3.6.5 |Anaconda custom (64-bit)| (default, Apr 26 2018, 08:42:37) \n",
      "[GCC 4.2.1 Compatible Clang 4.0.1 (tags/RELEASE_401/final)]\n",
      "numpy version 1.14.2\n",
      "scipy version 1.0.1\n",
      "tensforflow version 1.5.0-dev20171116\n",
      "keras version 2.1.5\n",
      "sklearn version 0.19.1\n",
      "powerlaw version 1.4.3\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Using TensorFlow backend.\n"
     ]
    }
   ],
   "source": [
    "import RMT_Util"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create Random Matrix  \n",
    "\n",
    "$$N\\times M=1000\\times 500$$\n",
    "\n",
    "$$\\mu=0,\\; \\sigma=\\dfrac{1}{\\sqrt{N}}$$\n",
    "\n",
    "#### will add rank 1 perturbation\n",
    "\n",
    "using Python API\n",
    "\n",
    "https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.normal.html\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-29T16:54:13.447370Z",
     "start_time": "2018-08-29T16:54:12.903252Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "% spikes outside bulk 0.00\n",
      "% sigma 0.9953\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "RMT_Util.ipynb:7: RuntimeWarning: divide by zero encountered in true_divide\n",
      "  \"## RMT Util\\n\",\n",
      "RMT_Util.ipynb:7: RuntimeWarning: invalid value encountered in sqrt\n",
      "  \"## RMT Util\\n\",\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.9953047167192505"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEICAYAAAC6fYRZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd8FHX+x/HXJyEhVGkBpAhIFRRB\ningiIDawcZ7oASqKKHj+sICeXVRsJyrq2SvqWRA7d+JhA8sJSFRAqdKJSBWQaiD5/v6YCS7LhmxC\nksls3s/HYx/ZnfnOzGd2sp/9zne++x1zziEiIokpKegARESk+CjJi4gkMCV5EZEEpiQvIpLAlORF\nRBKYkryISAJTkpcSZWZbzezQoOMQKSuU5PfDzJ42syciXqeY2bY8pnUxsxZm9r6ZrTOzX81skpm1\n3M/6XzSzLD/x5T5mRcwfbGbzzWyLma0xsw/MrErUslv8x49mdq+ZHVSA/ZtiZs7Mjoya/p4/vYf/\n+nYz2+XHt8nMvjazY/JYZzUze8HMVvtxLTSz63PnO+cqO+eWxBtjSTOz8n78v/n7MCKfsg+Z2Soz\n22hmT5hZSsT8w8zsMzPbbGaLzOysiHmN/fc48tjfGjF/tJmt9ONYbmY3R227nZl9a2bb/b/tYsSX\n6v//ZEZMOy5qm1v9OM725x/u/9+uN7OYP6Ixs35mNs//v19sZsf508+LWu92f90d/PnHm9lk//1Y\nFmO9fzKzb/z/m9lm1jVi3mlm9pX//7fazJ7N/Sz482N9lpIj5p/gvxfb/Rgaxbts6Dnn9MjjAZwH\nzIl4fQwwD/gxato2IAXoDAwGaviv7wTm72f9LwJ35TGvO7AGaO+/rgFcCFSJXhZIAzoBk4EfgUpx\n7t8UYAHwYMS0msBqYC3Qw592O/CK/zwFGA38AliMdY4FxgPV8SoRrYC+QR/LAhzze4Ev/fgP89+L\nXnmUvc0vWwNIB6YBd/jzygELgRFAMtDT/z9p4c9vDDigXB7rbpl7HIH6wBzgL/7rVGA5MBwoD1zp\nv06NWsfNwBdA5n72twewJWJbLf3/4T5eetin/En+trr4x7c+UD+PdV8ELM79P/E/HxcAQ4BlUWVr\nAOuBc/z363xgI1Ddnz8A6AVU9I/Nh8BTcX6WagGb/XWnAfcD0+JZNhEegQdQmh9AAyAHqOW/vs7/\nYC+NmvZJHsvX8D/INfOYv79/zGuB9/YT2z7LAlXwku+wOPdvCjASyASS/WnDgCf9aT38abfjJ3n/\ndRt/v2rFWOePwJ/3s00HNPOf1wT+DfwGzADuAr6KKns58JOfiO4EmgJT/WXG4yc2/4P/H2Cdnxz+\nAzQoxDH/GTg54vWdwLg8ymYA50S8HgCs9J8fDmwl4osQ+Ai403/emP0k+ajt1Ad+AK7zX5/sxxm5\n7hVEfBkBTfAqJL3Zf5IfC4yNMb0ZsZP818DgON/LycBtMaafyL5J/nQiKlT+tIV5bQv4C/BDnJ+l\nIcDXEa8rATuAVvktmwgPNdfsh3MuE6/Wcpw/qRteze3rqGlf5LGKbsBq59yGQmx+OnCKmd1hZsea\nWfk44t0CfBwRWzxWAXPxEgfAQODlvAr7cVyElzjWxygyDbjbzAaZWfN8tv04Xu22Lt5ZyoUxyvQC\nOuDVHK8DnsE7w2qIl0j7++WS8BJWI+AQvA/xYxFxP+Gf6sd6zPbLVAfqAXuazPznbfKI3/xH5OsG\nfpOZ5VH+8Khpy80s08zGmlmtvQqb3WBmW/G+cCsBr/mz2gCznZ+hfLOj4nwUuMl/H2IHb1YR6Au8\nlFeZqPLJQEcg3W9+yjSzx8ysQoyyjfD+//P8X4pehH3fs1jvV65ueGc3kS43r5n029zmJ18bIo6p\nc24b3hlGmziWDT0l+fx9DnQzsyS8081peIk+d9qxfpm9mFkDvCSWZ5uu79qohPMSgHPuS7zaylHA\nB8AGMxsTR1vhKrwziIJ4GRho3vWDas65qTHKnGtmm4CVeEn3z3ms6wrgVbwzgrl+MugdXcjfj7Px\nanrbnXNziZ1s7nPO/eacm4N3lvCRc26Jc24z3il7ewDn3Abn3Nv+urYAd+M1eeHPv9w5Vy2PR1u/\nWGX/7+aI7W/GO0OK5UPgKjNLN7O6eM0m4DUpzMdr8vq7eddtTvbjqeiXWY/XxNYI7/2s4r9vezjn\n/uFPPwr4V0RclaNi3CtO89r+yznn3s0j7lxn+3Hs8/+bhzp4zXV98SoS7fDe/1tilB0IfOmcWxrn\nur8G6plZf//9uhDvrK1idEEzOwmvQjAyYvI/geZAbeBW4EUzO9aft9/3K59lQ09JPn9f4NUajgCW\nOOe2A19FTKuAV+vew8zS8U7Nn3DOvZ7P+h+ISjh7arPOuQ+dc2fgJe0+eDXoS/JZX33g13h3zvcO\nXpvxFXjJJJbxfny1nXM9nXPfxirknNvhnLvHOdcBrzlmPPCmmUV/8aTjtVuvjJi2kn2tiXi+I8br\nyuDVSs27UL7czH7DO27VCngBbav/t2rEtKp4TUWx3A18D8zES1LvAbuAtc65XXhfhKfhtetfg/de\nZAI457Y65zKcc7udc2vwvhRPNrPIbeM83/v7ekdEnHuVy43TzCrhXTO5Io79vRB4OeqMYH9yzwoe\ndc794p/JjQFOjVF2IHGeIYD3JY33Pz4C7xj3Aj7Bf79ymVkXvDOavs65hRHLf+d/0e92zk3E+8L8\niz87z/crjmVDT0k+f18AR+J9WL/0p83Bay44DZjhnNuZW9g/5f8ImOCcu7soAnDO5TjnPgU+I+/T\nV8ysMl5755d5lclj/dvxaqV/I+8kX2DOud+Ae/CaGppEzV4H7Ma77pGr4QFs7hq8i4ZHO+eq4n0J\ng98EYGZP2b69SnIfc/x4N+Jd04jsbXQk+zYL5O7fDufcMOdcfefcocAG4FvnXLY/f7ZzrrtzrqZz\n7hTgUOCbPOLPTbSxmnnA+0Js6j+fA7Q1s8iybf3pzfHa+780s9V4X+AH+z1SGucWNrOGeBdd421O\nyX1/MiNijcmvBdcD3op33f76P3fOdXLO1cC7QNuSiPfLzNoDE4CL/c/DflfHH+/lHCKOqf9F2JQ8\njmvUsqGnJJ8P59wivJrFVfjJ06/5TPen7WmP92thk4D/OeduOJDtmlkf87qqVTdPZ7zT/WkxypY3\nr5vae3gXHcf603O76TWOY5M3Ad2dc8sOMO5bzayTed330vDeo014vXj28BPhO8Dtfi28FV7tr7Cq\n4NU0N/lnDbdFbe8y53XfjPWIbJt9GbjFf99bAZfiXZiLta/1zayef3y64J3q3xYxv62Zpfn7dy1w\ncO66zOxoM2tpZklmVhOvyWCKc26zP21o1LH/PyA3sU0BsoEr/WM/zJ/+GV6TVkO8ppR2eGd+a/zn\nkWdKF+BdjFwctU/mH7dU/3Wa7X09aCxwhZnV9is0V+Nd5I50IfC232wWue4kf90p/qbSzCw1Yn57\nv6mmKvAA3nWfSf68w4H/Alc45/4d41j0NbPK/jZOxuudM8Gf/S5wuJmd7W9/JN41jflxLBt+hb1i\nW5YewOt4vWxqRky7Du8b/5SIaRf607bhnSLmPg7JY70vAllRZdf787rhfajX451WLsTvXRG17BZ/\ne3OA+/Da1HPLHAcsA1Ly2P4U4JI85uXZuyaf9+oWvETzG16z0RTgTxHzI3vXpONdb8jtXXMf8Gms\nsv7rr4CLIl7fBTznP6/nb2ur/14NJc7eK1Hxlwde8GNaA4yImHdI5PH0j9EyYDvel9h5Ueu6H+9L\ndyvemVLkvvTH66W1De/s4WWgrj8vCS+h/RqxPzexd2+a9sC3eF9s3+F3tY2xPz2I0bsG75rBPj1X\n+KPXT+RjWcT8FOAJvC/u1XhfTmkR89P8eSfkEUv0uqdEfc42+483gNoR88bifQYjPyuR3Zu/9Jf7\nDe8ia7+obZ/o7/MO//+kcbzLhv2R239VEpCZ3QKsc849HXQs8TCz+/ASXaxeNiJSCEryEhi/OSQV\nrw94J2Ai3pnFe4EGJpJAygUdgJRpVfBO0evhdTd8EHg/0IhEEoxq8iIiCUy9a0REElhgzTW1atVy\njRs3DmrzIiKh9O233653zqXHWz6wJN+4cWMyMjKC2ryISCiZ2fKClFdzjYhIAlOSFxFJYEryIiIJ\nTEleRCSBKcmLiCQwJXkRkQSmJC8iksCU5EVEEpiSvIhIAku4USinTPnjeY8eQUUhIlI6qCYvIpLA\nlORFRBJY4iV556j6w/84aPaXkJMTdDQiIoFKuDb5Fg8Opdrsz3HJKfBuOrz9NtSoEXRYIiKBSKya\n/IQJVJv9ORnPfM+M52ZB27Zw+umwY0fQkYmIBCKxkvyDD7L04rvISasIycnw0ENQvz5cf33QkYmI\nBCJxkvySJTB/Puu7/vmPaUlJ8Mwz8N578PHHwcUmIhKQxEnykyZB7964cil7T69eHZ54AoYNg6ys\nYGITEQlI4iT5jz+Gk0+OPe/006FZM3j44ZKNSUQkYImT5KdPhz/9Ke/5Dz8Mo0fDmjUlF5OISMAS\nI8mvXQvbt0OjRnmXad4czjsP7r235OISEQlYYiT5WbOgXTsw23+5m26Cf/0LVqwombhERAKWGEl+\n5kwvyeenTh247DIYNar4YxIRKQUSI8nPnw+tWsVX9tpr4f33YcGC4o1JRKQUSIwkv2wZNGkSX9nq\n1WH4cLjjjmINSUSkNEiIsWt2zFvK7LVN2DFl7+lTpsQqDcntruDoB5oxa+yPdBp0eHGHJyISmPDX\n5Hfvpvz6n9lZ+5C4F8muWIWV515Lo5dUmxeRxBb+JP/zz2RVq41LLV+wxf78fxz0w1feRVsRkQQV\n/iS/dCk76zYu8GI5aRVZ2f96uP32Ig9JRKS0CH+SX7WK39MbFG7RM4ZCRob3EBFJQOFP8mvWsKt6\nnUItmlO+Atx4I9x2WxEHJSJSOoQ/ya9dS1a12oVf/pJL4IcfYOrUootJRKSUCH+SX7OGrELW5AEo\nXx5uuUW1eRFJSAmR5HdVP4CaPMCgQbBoEXz5ZdHEJCJSSsSV5M2sl5ktMLNFZnZDjPmHmNlkM/ve\nzGab2alFH2oe1q49sJo8QEoK3HorjBxZNDGJiJQS+SZ5M0sGHgd6A62B/mbWOqrYLcB451x7oB/w\nRFEHmqcDba7JdcEFkJkJn3124OsSESkl4qnJdwYWOeeWOOeygHFAn6gyDqjqPz8IWFV0Ie6Hc7B2\n7YE31wCUK+e1y48c6a1XRCQBxJPk6wMrI15n+tMi3Q6cb2aZwETgilgrMrMhZpZhZhnr1q0rRLhR\ntmyBlBSvK2RR6N8fNmyAjz4qmvWJiAQsniQf604c0VXd/sCLzrkGwKnAv8xsn3U7555xznV0znVM\nT08veLTRNm70RpUsKsnJ3i9gVZsXkQQRT5LPBBpGvG7Avs0xg4HxAM65qUAaUKsoAtyvzZuhWrWi\nXec553i3Evzgg6Jdr4hIAOJJ8jOA5mbWxMxS8S6sTogqswI4AcDMDsNL8kXQHpOPTZuKPsknJXlj\nzas2LyIJIN/x5J1zu81sGDAJSAZecM7NMbNRQIZzbgJwDfCsmQ3Ha8q5yLkSyJBFkORjjjlf7c90\n2HIny+96j/XHnUWPHvkvm1cZEZEgxXXTEOfcRLwLqpHTRkY8nwscW7ShxaE4avIASUksGzSKJs/f\nzPpj+5AIvxkTkbIp3NmruJI8sOGY08lJTSP987eKZf0iIiVBST4vZiwbNIrGL90O2dnFsw0RkWKm\nJL8fv3Y6hd2Vq8G4ccW2DRGR4qQkvz9mLB10p9fbZvfu4tuOiEgxUZLPbxNH9YQGDeC554p1OyIi\nxSHcSX7LFqhSpXi3YQYPPuiNa7NpU/FuS0SkiIU7yW/dCpUqFf922reHM86Au+4q/m2JiBShcCf5\nbdtKJsmDl+BffNG7uYiISEiEP8lXrlwy26pbF669Fv7+95LZnohIEQh/ki+pmjzA1VfDzJkweXLJ\nbVNE5ACEO8mXVJt8rrQ0GD3aS/bqUikiIRDeJO+cNyRwxYolu92+faFWLXjssZLdrohIIYQ3ye/Y\nAeXLezf6KElm8PjjcNddpK77uWS3LSJSQOFN8iXdHh+pVSu47DKaPXlNMNsXEYlTXEMNl0bTPtlK\nu3KVmTal+LcVc9z4m26iSrPDqZ7xMRs7nhR7XHryHmdeY9GLSEkIbU0+eec2stMCqskDVKzIoiv+\nSfOHLycpa2dwcYiI7Edok3zSjm1kVwgwyeONOb+tyeEc8to/Ao1DRCQvoU3ygdfkfYuufJR67z1O\npSU/BB2KiMg+Qp3kc0pBkv89vQFLL72XVvcNwrLVd15ESpfwJvkdW0tFTR7gl1MHs6tKdRqMfzDo\nUERE9hLaJJ+UtZPstBL+IVRezFh47bMcMu5+Kq6YH3Q0IiJ7hDrJu5TyQYexx866jVl60R20HH2x\n7gkrIqVGqJN8Tmpa0GHsZVWfv+GSU2ioZhsRKSXCm+R3/U5OKarJA5CUxLwbX6bh+Aeo/NP3QUcj\nIhLiJF8Ka/IAv9dtxKLLH+Kwu88jaef2oMMRkTIutEnedv1eKpM8wNqTzmNr03Y0ffq6oEMRkTIu\ntEk+KWtn6WuuifDT8CeoOfXfMHFi0KGISBkW4iRfemvyALsrV2PeDS/DJZfA6tVBhyMiZVSIk3zp\nrskDbG7XHYYMgQED1K1SRAIR7iQfQE1+ypQ/HnG59VZISoLbby/Z7YqIEOYkX4ovvO4lORlefRXG\njoX//jfoaESkjAlvks/aiUst3c01e9SpA6+9BhddBCtXBh2NiJQh4U3yYanJ5+rWDYYPh3PPhd9/\nDzoaESkjwpvkQ3DhdR9//zvUrQvDhoFzQUcjImVAaJN8af4xVJ6SkuDll2HqVOq/93jQ0YhIGRDa\nG3mHsiYPUKUKTJhAo45/Yluj1mw6qmfQEYlIAourJm9mvcxsgZktMrMb8ihzrpnNNbM5ZvZa0Ya5\nr9I6dk1cDj2Uube8Ruu7BpC2aknQ0YhIAss3yZtZMvA40BtoDfQ3s9ZRZZoDNwLHOufaAFcXQ6x7\nCd2F1yibjurJ8vNv5vBb+sDmzUGHIyIJKp6afGdgkXNuiXMuCxgH9IkqcynwuHNuI4Bzbm3Rhrmv\n0DbXRPj5rGFsbtsN+vaFXbuCDkdEElA8Sb4+ENm5O9OfFqkF0MLM/mdm08ysV6wVmdkQM8sws4x1\n69YVLmJf2GvyAJix6IpHoEIFb/gD9bgRkSIWT5K3GNOis1E5oDnQA+gPPGdm1fZZyLlnnHMdnXMd\n09PTCxrrXkL1Y6j9cMnl4PXXYc4cuOOOoMMRkQQTT5LPBBpGvG4ArIpR5n3n3C7n3FJgAV7SLx7O\nkbQri5xyqcW2iRJVqRL8+99e98qxY4OORkQSSDxJfgbQ3MyamFkq0A+YEFXmPeB4ADOrhdd8U3zd\nRrKyyCmX4vU7TxR16sCHH8KNN8J//hN0NCKSIPLNks653cAwYBIwDxjvnJtjZqPM7Ey/2CRgg5nN\nBSYDf3fObSiuoNm1C1cupdhWH5iWLWHCBLj4Yg03KSJFIq4fQznnJgITo6aNjHjugBH+o/jt3u21\nZSeizp3hjTe8MW4++AA6dQo6IhEJsXBmyl27cMnhqMkXqkJ+/PHw3HNknXIGMx/8lO1N2sS9/h49\n8i+XVxkRSTzhbNTetctrk09kZ57JosvHcOR1p5D28+KgoxGRkAptkk/INvkoa08cwLILbqXdNT2V\n6EWkUMLbXFMGkjzAL2cOBaDdiOOZNeYzdtRvFnBEIhIm4U3yIWmTLwq/nDkUkpI4ckRPZj34KcX5\nEwQRSSzhTPKJ3LsmD7+cfimYceSIntD5U2jRIuiQRCQEwpkpy8KF1xh+Oe0SnCXRqmdP76bghx8e\ndEgiUsqFNsmXlTb5aKtPvZhW7SvAiSfCe+8BXYIOSURKsfD2rilDbfL76N8fXngBzjiD6hkfBx2N\niJRi4U3yZbQmv8epp8K773LYPeeTPuXNoKMRkVIqtEm+LLbJ76NrV2aN/ohmj13FwROeDjoaESmF\nwtkmXwZ71+RlW7MjmfnIFxxxfW/S1iyHbncl1uicInJAwpkN1Fyzlx31m/H9Y19TbdbnMGAA7NwZ\ndEgiUkqEN8mX5QuvMeyqls6sMZ96txA88URYvz7okESkFAhvkldNfh85qWnerQSPOw6OOQZ++ino\nkEQkYKFN8rrwmoekJLj3XrjuOujaFT5WF0uRsiycVy8TrLmmKG8CtWddzS+lx/gW0K8fXHstHDUC\nzPa7vchx5uMpk28M+1mnxrQXKRnhrMmrd018uneH6dPh1Vdpde9Akn7fEXREIlLCwpnk1SYfv0MO\nga++wnJyaH9FV8qvWRF0RCJSgpTky4KKFZl38yusOXEAR/1fF6p/+0nQEYlICQltkteF1wIyI/Pc\na5h386u0+seFNH7xdsjODjoqESlmoU3yiXThtSRtan883z79LQfN/oIjrzuFlF/XBB2SiBSj8CZ5\n1eQLLatGXWbd/zGb2/yJjkM7UG3mlKBDEpFiEs4kr941By45mWUXj2L+dS9w2J39afTSKNi9O+io\nRKSIhTPJqyZfZDZ2OtlrvvnhS6/L5ZIlQYckIkUotEleF16LTlateswePQn69oWjj4aXXvLGwBGR\n0Atnks/OhqTkoKNILElJMHw4fPopPPAAre84l3K//Rp0VCJygEKb5J2SfPFo2xZmzCCrVn06XnIk\n1b+ZFHREInIAlORlX2lpLBr2MAuue4GWY4bQcvRgym3dFHRUIlII4UzyOTlK8iVgY8eTmPH8D+Sk\npNLx4iOoMW1i0CGJSAGFM8lnZ+sWdyUku1JVfhr+JPNveInm/xwGF10EGzcGHZaIxCmcmVLNNSVu\n01E9yXh+NlSpAkccAe+8ox44IiEQzl8UlfEkH+/480U5Tj1AdoXKTDn7UQ5qfi4tRgyl0tix8Nhj\n0KhRkW1LY9GLFK1w1uRzciC57Cb5oG1uexwZz86ELl2gQwe4/35s966gwxKRGMKZ5LOzcRbO0BOF\nS0mFm2/2bkryySd0GNqBqnOmBh2WiESJK1OaWS8zW2Bmi8zshv2U62tmzsw6Fl2IMZTx5ppSpWlT\n+O9/WXHeTbS57WxaPDiUlM3rg45KRHz5JnkzSwYeB3oDrYH+ZtY6RrkqwJXA9KIOch/Z2WquKU3M\nWNuzHzNenEtOahqdLmpN/XcexbI14JlI0OKpyXcGFjnnljjnsoBxQJ8Y5e4ERgM7izC+2FSTL5V2\nV67GoiseYeaYydT8egIdL2kHn+guVCJBiifJ1wdWRrzO9KftYWbtgYbOuf8UYWx5y8lRm3wptr1J\nG2bf/xFLB98FQ4bAWWdpdEuRgMSTKS3GtD0dpM0sCXgIuCbfFZkNMbMMM8tYt25d/FFGU3NN6WfG\n+q5/hrlzoVMn6NyZJs/eRPLWzUFHJlKmxJPkM4GGEa8bAKsiXlcBDgemmNkyoAswIdbFV+fcM865\njs65junp6YWPWs014ZGWBjfdBLNmkfrrLxw9sAX133kUsrKCjkykTIgnyc8AmptZEzNLBfoBE3Jn\nOuc2O+dqOecaO+caA9OAM51zGcUSMSjJh1H9+iy4fiyz7v+YGtMnQuvW8Oab+tWsSDHLN8k753YD\nw4BJwDxgvHNujpmNMrMzizvAmDRAWWhta9qWH+77EJ56Cu69F445Br74IuiwRBJWXMMaOOcmAhOj\npo3Mo2yPAw8rHxqgLPxOPBEyMmDcOBg40BvH/q67gLZBRyaSUMKZKdVckxiSkmDAAJg/H3r2hJNP\npvWoflRcMT/oyEQShpK8BC8tDa6+GhYvZkvz9rS7qhut7r0QFi8OOjKR0AtnktcAZYmpUiVW9r+e\n6a8sYke9pt5NxS+9FFasCDoykdAKZ5LXAGUJLbtSVZZfOBIWLoTataF9e5o/Mozya5TsRQpK48nL\nHgUdE744xpDfS40acPfdcNVVZF/5AB2HtGfdcX9h2s83sLN+UyC+MefjHYs+rzgKui6NiS+lSTir\nw/rFa9lSuzZLLhvN9JcXklWzHh0uP5rD7j6fisvmBh2ZSKkXziSvfvJl0u6DarJs0B1Me3Ux2xq3\nod2I46FvX/j++6BDEym1wpnk1SZfpmVXPogV593ItFeXQNeucPrpcNppHDTrC/2CViRKODOlmmsE\nyKlQaU/XS848k5b3D+aoy7uQ/vlb3v+IiIQ3yau5RvZIS4OhQ/nmpfmsGHADDd4cAy1awOOPw/bt\nQUcnEqhwJnm1yUssycmsP+4svn/sa3j5Ze+GJY0b03jsbaRsXBt0dCKBCGeS19g1kp9jj4V334Uv\nvyT119V0HtiSFmMu84ZQEClDwpkp1Vwj8WrZkoXXPM03Ly8gq3od6N4dTjkFPvjA++W0SIJTkpcy\nYVf12iwbdAcsXw7nnQcjR0LLlvDII/Dbb0GHJ1Jswpnk1SYvhZWW5g1tnJEBL74IX38NjRvDlVdS\nYeXCoKMTKXLhTPJqk5cDZea127/xBsyeDVWq0P7Krhxxw6nUmP6hmnIkYYQzU6q5RopSgwZw991M\nG7ecdd3PoclzN3H0+c1o+Pp96pUjoackL+LLKV+B1b0H8e0z3zH31nFUXLGAoy9o4d3Y5Av9mlbC\nKZxJXuPJS3EyY8thnVlw/QtMe32pN6790KHQpg313/4n5bZuCjpCkbiFM8lr7BopIburVIerroK5\nc+HJJ6k6dypH929Cy9GDYcaMoMMTyZfGk5dildcY7fGMRX8gZQ5kuzGZQffuzLu1Oykb13Lwhy+w\n48xz2V2lOr/0HkyL2wdA9eqFWnVBY4pnLPqCjl2vse4TVzirwxqgTAK0q3ptVgy4gemvLmbJkPuo\n9sOX0KSJ1//+s8/UM0dKlXDW5HNy1FwjwUtKYmPHk9jY8SRqH7EBXn0Vhg/n6PVbWd1rEKt7XQQ0\nCDpKKePCmynNgo5A5A81a8KVV8LMmcwd+Qap63+m4+C2cOqp8PbbkJUVdIRSRoUvyed2Y1OSl9LI\njC0tO/LT8CeZOj4T+veHRx+Fhg1p+sQ1VFr6Y9ARShkT3iQvUsrlpFWECy7wrmr+73/kpJTniOt7\n02FIB2/MnLX6oZUUv3AmedWSPvdTAAANcElEQVTiJWyaNWPppfcw7fVlLB46Gr77zruxyRlnkP75\nWyRl7Qw6QklQSvIiJSk5mU0dToCXXoLMTDjnHOpNeIpjzqlPizGXUfXHr3W2KkUqfL1rlOQlUVSu\nDAMHMuuQgZRfu5I6n7xKy/sHY9m7WXPyQNacfAE76zYOOkoJOdXkRUqB32s3ZMWAG5jx4lzm3fIa\nKZvWctRlnWh3VXfqfvC8hlKQQlOSFylNzNjSqhOLrnyUqW/+zMpzRlBz+kS69GtEm5FnU+vzt2Gn\n2u8lfuFL8qAkL2WCS0llQ9c+zBn1NtPGLWdDl9OoP+FJqFcPBg3yblSenR10mFLKhS/J66KUlEG7\nK1dj9akXM+vBT+DHH6FtW7jxRm8s/OHDvcHS9NmQGMKZ5FWTl7KsXr0/Evvnn8NBB3lj3rdoAbfd\nBgsWBB2hlCJK8iJh1qIF3H47LFwIr73m3ZS8Rw86DO1Ig/FjSF33c9ARSsCU5EUSgRl06gQPPQSZ\nmSwZ8g8qLf2RToOPgBNOgOeeg19/DTpKCYC5ONrxzKwX8AiQDDznnPtH1PwRwCXAbmAdcLFzbvn+\n1tmxY0eXkZFR8Ii3boU6dZjywbaCLytSSkSO2V7oMe7jWE9S1k5qTP2A2pPHUSPjI8p17wr9+kGf\nPlC1apHEEB2HFC8z+9Y51zHe8vnW5M0sGXgc6A20BvqbWeuoYt8DHZ1zbYG3gNHxh1xAqsmLxC0n\nNY313c9m7u1vegOmnXcevPkmNGwIf/kLjB8P27cHHaYUo3iaazoDi5xzS5xzWcA4oE9kAefcZOdc\n7n/KNIpzEG0leZFCya5YxbtAO2ECLFsGZ5wBzz8P9epx2J39qfnV+1jW70GHKUUsniRfH1gZ8TrT\nn5aXwcCHsWaY2RAzyzCzjHXr1sUfZSQleZEDV72619d+0iRYuJDNbbvR8K2H+FPfg2l53yBqfPNf\nbPeuoKOUIhBPko+VUWM25JvZ+UBH4P5Y851zzzjnOjrnOqanp8cf5d4rUZIXKUq1a7Oqz9+Y+fAU\nZjz/A9uaHkmjl+7gmL71aDHmMqrNnKIfXYVYPAOUZQINI143AFZFFzKzE4Gbge7OueI751OSFyk2\nWen1yex7NZl9rybtl6WkTxlP08eHk7pxDeu6n8Panv34rXUXfQZDJJ4kPwNobmZNgJ+BfsCAyAJm\n1h54GujlnCveOyEoyYuUiJ0HN2Fl/+tZ2f96KqxYQO3Jb9By9MUk/76DtT3OZW3Pfmxt3l6fx1Iu\n3yTvnNttZsOASXhdKF9wzs0xs1FAhnNuAl7zTGXgTfMO+Arn3JnFErGSvEiJ23FIS5ZfOJLlA2+l\n0pIfqD35DdrccQ4uKZm1x/8V0vtBmzZBhykxxDWevHNuIjAxatrIiOcnFnFc+wtGSV4kKGZsa9qW\npU3bsnTwXVRZkEHtyW9Ar15QrRr89a9eP/xmzYKOVHz6xauIFI4/LPLivz0Ay5fDk0/C6tVw7LHe\nr2/HjIGfNaxC0JTkReTAJSVB167w2GNeYr/nHm+0zCOOgOOPh2eegQ0bgo6yTFKSF5GiVa4cnHQS\nvPACrFoFV17pjX1/6KHeD7Bee80bnkRKhJK8iBSftDQ46yxv+ITMTDj3XHjlFahf32u7f/99+F2/\nsi1OSvIiUjKqVIELLoCJE2HxYm9UszFjvPHxL7kEPv1UP7oqBkryIlLyatWCyy7zbnoycya0agXX\nXefd6erqq2H6dN3pqogoyYtIsBo2hGuvhW+/9cY8rl4dBg70umHefDPMnx90hKEW13jyxaHQ48ln\nZsLRRzPlVXXNEklYzlH5p++p8+lr1P70NX5Pb8jqkweytmc/uvapGXOReMbEz2v8/bzGw8+rTDzL\nxruugiry8eRLHdXkRRKfGVtbHMXivz3AtDdWsHTQKA6a8zVdzmvqXch9913Iygo6ylBQkheRUs0l\nl2Nj51OYd8urTH1jhdcN85FHvB46w4bBN9+o/X4/lORFJDSyK1WFiy/22j5mzIA6dby7XR12GIe8\ncg/l167Mdx1ljZK8iIRT48Zw662wcCGMHUv5dSvpeGk7jri+N7W+eEc3PfEpyYtIuJnBMcfw0/An\nmTo+k7UnDKDB2w/T5a+HcOgzN1Ah86egIwyUkryIJIyc8hVYc/IFzHzkC2Y9NBnL3k37K47lyOHH\nU/uT12DnzqBDLHFK8iKSkLYf0orFf3uAqW+sZFWfy6k76UXvx1YjRsCiRUGHV2KU5EUkobnU8qzr\ncQ6z7//Iu1ibkgLHHMMRN55O9W8mQU5O0CEWKyV5ESk7mjSB++6D5ctZ1/Usmj59HZ0vOgwefRR+\n+y3o6IqFkryIlD0VK7L6tMFkPDeTBdc8C1984fXWufpqWLEi6OiKlJK8iJRdZmw+shu8+SbMmgWp\nqdC+PVx4IcyZE3R0RUJJXkQEvIHSRo/2Lsq2bAknnMDhN59J1R+/DjqyA6IkLyISqXp1uOkmWLqU\nXzv35rB7zufIESfAtGlBR1YoSvIiIrFUqMCqPn/jm5cXsLZnP++uVmec4Y1/HyJK8iIi++HKpfDL\n6Zd6wyecdBL07k3rUf1IW70s6NDiUi7oAApMSV6kTItn3PiCLhvPOqdMS4O2V5L0wmAavjmGDkM7\n8HOf/4PON0DFijHXdSDjxhcV1eRFRAogp0Illg+8lYxnZ1Lh55+8Wxe++WapHe5YSV5EpBB+r92Q\nebe+Dq+8ArfdBn37wpo1QYe1DyV5EZED0a0bfPcdtGgBRx5J+uTxQUe0FyV5EZEDlZYG994LEybQ\n5Pmbaf7w/2FZvwcdFaAkLyJSdDp35runZpC6YRXtr+oGq1YFHZGSvIhIUdpduRpzRr3DhmNOh65d\n4adgb1qiLpQiIkXNjOUDb6XJMXWhe3eYPBloGUgoSvIiIsXl0kuhXDno1YvU+/9HVq16JR6CmmtE\nRIrToEFw8cW0HvVXLHt3iW8+fEm+YUPv9l0iImFx883kpKbRcNz9Jb7p8CX5unVhwICgoxARiV9S\nEgtHPE3D8Q+Q8mvJ/mAqfEleRCSEdtY7lHXd+nLwxOdKdLtxJXkz62VmC8xskZndEGN+eTN7w58/\n3cwaF3WgIiJht/qUC6nz6eslus18k7yZJQOPA72B1kB/M2sdVWwwsNE51wx4CLivqAMVEQm7LYd1\npvzaFSV60/B4avKdgUXOuSXOuSxgHNAnqkwf4CX/+VvACWbqAiMiEskll2Pqmz9D1aolts14+snX\nB1ZGvM4Ejs6rjHNut5ltBmoC6yMLmdkQYIj/cquZLYixvVrRyyWIRN0vSNx9036FT6LuW+R+NSrI\ngvEk+Vg18uiBk+Mpg3PuGeCZ/W7MLMM51zGOuEIlUfcLEnfftF/hk6j7diD7FU9zTSbQMOJ1AyB6\n1J09ZcysHHAQ8GthAhIRkaITT5KfATQ3syZmlgr0AyZElZkAXOg/7wt85lwpvU2KiEgZkm9zjd/G\nPgyYBCQDLzjn5pjZKCDDOTcBeB74l5ktwqvB9zuAmPbbnBNiibpfkLj7pv0Kn0Tdt0Lvl6nCLSKS\nuPSLVxGRBKYkLyKSwAJL8ok6VEIc+3WRma0zs5n+45Ig4iwoM3vBzNaa2Y95zDcz+6e/37PN7KiS\njrEw4tivHma2OeJ4jSzpGAvDzBqa2WQzm2dmc8zsqhhlQnfM4tyvsB6zNDP7xsxm+ft2R4wyBc+L\nzrkSf+BdwF0MHAqkArOA1lFlLgee8p/3A94IItZi2K+LgMeCjrUQ+9YNOAr4MY/5pwIf4v1mogsw\nPeiYi2i/egD/CTrOQuzXwcBR/vMqwMIY/4uhO2Zx7ldYj5kBlf3nKcB0oEtUmQLnxaBq8ok6VEI8\n+xVKzrkv2P9vH/oALzvPNKCamR1cMtEVXhz7FUrOuV+cc9/5z7cA8/B+mR4pdMcszv0KJf84bPVf\npviP6J4xBc6LQSX5WEMlRB+ovYZKAHKHSijN4tkvgLP90+O3zKxhjPlhFO++h9Ex/in0h2bWJuhg\nCso/pW+PVzOMFOpjtp/9gpAeMzNLNrOZwFrgY+dcnscs3rwYVJIvsqESSpl4Yv430Ng51xb4hD++\nlcMujMcrHt8BjZxzRwKPAu8FHE+BmFll4G3gaudc9NCHoT1m+exXaI+Zcy7bOdcOb2SBzmZ2eFSR\nAh+zoJJ8og6VkO9+Oec2OOd+918+C3QoodiKWzzHNHScc7/lnkI75yYCKWZWK+Cw4mJmKXiJ8FXn\n3DsxioTymOW3X2E+Zrmcc5uAKUCvqFkFzotBJflEHSoh3/2KavM8E69NMRFMAAb6PTa6AJudc78E\nHdSBMrO6uW2eZtYZ7zOzIdio8ufH/Dwwzzk3Jo9ioTtm8exXiI9ZuplV859XAE4E5kcVK3BejGcU\nyiLnSn6ohBIR535daWZnArvx9uuiwAIuADN7Ha/XQi0zywRuw7swhHPuKWAiXm+NRcB2YFAwkRZM\nHPvVF/ibme0GdgD9QlDZADgWuAD4wW/jBbgJOARCfczi2a+wHrODgZfMu1FTEjDeOfefA82LGtZA\nRCSB6RevIiIJTEleRCSBKcmLiCQwJXkRkQSmJC8iksCU5EVEEpiSvIhIAvt/jMfPm4xYIMYAAAAA\nSUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "W = np.random.normal(0.0, 1.0/np.sqrt(1000), (1000,500))\n",
    "Q = np.max(W.shape)/np.min(W.shape)\n",
    "evals, _ = RMT_Util.eigenspectrum(W)\n",
    "RMT_Util.plot_ESD_and_fit(eigenvalues=evals, Q=Q)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Add a random Rank 1 perturbation\n",
    "\n",
    "#### check that we see a spike\n",
    "\n",
    "$$\\mathbf{W}=\\mathbf{W}^{rand}+2\\mathbf{\\Delta}$$\n",
    "\n",
    "$\\mathbf{\\Delta}$ is a 1-D random vector , same mean and variance as $\\mathbf{W}^{rand}$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-29T16:54:13.991683Z",
     "start_time": "2018-08-29T16:54:13.449606Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "% spikes outside bulk 0.00\n",
      "% sigma 1.1695\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "RMT_Util.ipynb:7: RuntimeWarning: divide by zero encountered in true_divide\n",
      "  \"## RMT Util\\n\",\n",
      "RMT_Util.ipynb:7: RuntimeWarning: invalid value encountered in sqrt\n",
      "  \"## RMT Util\\n\",\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "1.1694848374954956"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEICAYAAACktLTqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmYFNXVx/HvmWEGEJBFQBFQVlFA\nRQGVGBQD7guJGtwjviKuURMTRU3c9xh3XOO+a4xKFOOOuwKKGwiIiIAgsomAwAwz5/3j1mDb9DAN\nzHRNF7/P89TDdNWtqtNF9+lbt27dMndHRESSpSDuAEREpPopuYuIJJCSu4hIAim5i4gkkJK7iEgC\nKbmLiCSQkrvklJktMbMOccchknRK7mtgZneY2a0pr4vMbGkl83Yxs63M7Fkzm2tmC8zsRTPrsobt\n32dmJVHCq5g+SVl+vJlNNLPFZjbHzJ43s0Zp6y6Ops/N7Eoza7wW72+UmbmZbZ82/5lofr/o9UVm\nVhrF94OZvWtmfSrZZhMzu8fMvovimmxm51Qsd/eG7j412xhzzcwGRe/vJzMbVUXZVmY2wsxmRcer\nXYYyA8zso+gzMsPMBqUsOzD6f1sS7bNrJft5Ldp+nZR5PczsLTNbZGYzzeyCSta9MFp3QMq8Zmb2\nuJnNi6aHzWzjlOWvR5/hH83sEzMbmLLsvLTP6zIzKzez5mn7bRZt4+20+UPMbEq07v/MbPOUZamf\ns4pptYqAmR0bvachKfOamNn9ZvZ9NF2Uts60KNaK7b6UtryDmT0XfWbnmdk1mY5nPlFyX7M3gd1T\nXvcCpgO7pc0D+BBoAowAugCbAqOBZ6vYxzVRwquYtgcws92BK4Aj3L0RsA3wRIZ1GwEtgOOAXYB3\nzKzBWrzHycAfKl6Y2SbRduamlXvc3RtG+3ob+I+ZWYbtXQ80jOJtDBwEfLUW8cRtAXADcFUWZcuB\n/wGHZFoYJetHgPMJx6IH4XOCmXUGHgZOInxu/guMSE3gUbmjgF/MizxC+Hw2I3xGTzazg9LW7Qgc\nCsxOW/cyoCnQAehI+KxelLL8DKCVu28MDAUeMrNWAO5+RernFbgaGOXu89L2cTXwRVo8FZ/pgVHc\nXwOPpq33eNr3YWraNpoC5wLj09a7HtgIaAfsBBxjZsellTkwZbt7pWyzGHgZeA3YDGgDPESeU3Jf\nszeAbVJqJX2Bx4AGafPec/dSdx/t7ne7+wJ3LyV84LpECXNt9Y62Ow4g2ub97r44vaC7L3f3MYRE\nugkh0WfrYeAwMyuMXh8BPA2UZCocva/7CV+CTO+rN/CIuy9093J3n+ju/65YGNW4OkV/b2Jm/41q\niGPM7LLUml5U9hQz+zKqUV1qZh3N7L1onSeiLyZm1jSqec01s4XR323W4jhUvL9X3P0JYFYWZee4\n+63AmEqK/A24w91fcPeV7j7f3St+6PYG3nL3t919JSEZtialMmHhLOxC4OwM224HPOzuZdE23wa6\npZW5BTiH1f8v2wPPuPuP7r6I8P+9al13/zSKCcCBIqBtegDRj/sxhM9D6vw+QHfg3rRVDgSedPfx\n7l4CXArsFv0IZetK4CYg/cfkQEJl5yd3nwbcDfxfltscDMxy9+vcfWn0ffp0LWKqlZTc18DdZwLf\nEBI4hBr7W8C7afPerGQTuwHfufv8ddj9B8DeZnaxme1qZnWziHcxoQbSt6qyKWYBE4CKmswfgAcq\nKxzFMRiYmaG2BvA+cLmZHRfVTtdkOLCU8ENxbDSl2wfoSTibOBu4EziKkGy6E36MIHyW7wW2BLYA\nlhGSW0Xct1poUso01dQXeZdo35+Z2Wwze8jMmlWEFE2kve6eMu8K4DbguwzbvgH4g4VmwS5AH+CV\nVRsz+z1Q4u4jM6w7HDgg+kFsSjjzeCG1QPTjuJzwORwFjM2wnb6EWv9TKesVRts/jfDD8IvNZnjP\npL3nAy00aY43s5PTYtqJcKZ8e4ZYUrdX8Xf3tOUPRz/+L9kvmyJ3AaaZ2QtRk8woM9u2kn3kD3fX\ntIYJuI9QAy8Aviec+p2UMm8hsHuG9doA3xKaVda07eXADynT/SnL9yWcrv8ALAGuAwpT1r0swzav\nAl7O8r2NAoYARxNOj7sAk6NlM4F+0d8XEWp/P0TH4DWgZyXbrA+cR2h+KAWmAPumLHegE1AYLe+S\nsuwy4O20srumvP4QOCfl9T+BGyqJowewcD3+34cQmhuyKVsnirVd2vwSYBqwFaGp6ilCbRtga8IP\nWz+gGPg7oZnn3Gh5L+DjaNvtou3XSdn2r6JjuzJadnHKsobAl0D76PU0YEDK8s0JPwTl0fQyUJzh\nfRVFn8E/VfK+7wbuS5v3J+C26O/Baf+f/Qk17u2iz8kd0f6PiJZ3jWIrjN7f7JRlhYQfmD6pn92U\nbT8E/AdoFH2+vgJWpCzfNdrnRoRmne+AJtGyl6LP4r7R/8VfgamZjkk+Taq5V+1NQg18W2Cqu/9E\nOAWumFefULtZxcxaED4wt7p7eptiumvdvUnKtKr26uF0/kBC++RAwpdlSCXbqdCa0G68Nv4D/Ab4\nI/BgJWWeiOJr6e6/cfcPMxVy92Ue2mV7EpptngCeTKmxVmhBSFwzUubNYHVzUv5eluF1QwAz28jC\nBfBvzOxHwv9bk5TmpjgsA+5198nuvoRQE98PwN0nEs5UbiEkseaEM6iZZlYA3Aqc4T83j6wSHcv/\nAZcA9QhnMXub2SlRkYuBB93960riepJwraURsDEhEa7WxuyhqfGFaNvp7fn1gd+T0iQTXRw9nXCN\nYTXu/iqhmekpwhnxNGAxoSKBu09w91kempreBW4kXDMAOAX41N3fq+Q9nU443l8SrnM9WrHdaNvv\nRJ/Nn9z9SkJFpeIMdxnhR+gFD81F1xI+u9tUsq+8oORetTeB7YH9CU0yEC7mtI3mjXH35RWFo9Pc\nl4AR7n55dQTgoe36VUKNOf1UcxUzawgMSIkz2+3/RDgtP5nKk/tac/cfCQmtAaGdN9VcQq0ztV18\ntXbdtXAW4cxjZw8XAisuehuAmd2e1gsjdUq/OFddPmX1polV3P3f7t7d3TchJL0tCe33GxNq7o+b\n2Xf83KY/08z6Ei6Elrn7Ax7a8mcSrgXtF5XrD5xuocfSd4Tj+oT93Gtpe8K1gKXRj87tKetmUodw\n4TXVwYRKxKiUeTsBrYAJ0X5vBHaK4iiM3vNwd+/s7i0JSb4O8Hllh4ifm1r6A79LeU+/Av5pZrdE\n213g7ke5+2bu3o2Q20av4T2lbnuN/0/5Ssm9Cu4+hVBbPIMoaXo4l/sgmreqvd1Cd7IXgXfcfdj6\n7NfMBprZ4VG7qEXtjbsT2rTTy9Y1s57AM4Rmonuj+e2ski56GZxHaF6atp5x/93MeptZsZnVIxyj\nH4BJqeXcvYxwxnBRVOvempReO+ugEaEG9kNUs70wbX8n+S97YaROqy4mmllhFHcdoMDM6plZ0Rre\nbz2g4npI3eh1hXuB4yx0s9uIcHHzuZR1e0b7a0FoovhvVKNfRGie6BFNFYm3J+FzNzmsbkeaWYGZ\nbQYcBlR0o+1PqARUrD8LOJHQFg7hx2KImdWPauBDK9Y1s63NbN9oWZGZHU34oXwj7a0fCzwQfRcq\nvEBoQqrY7wXAOKCHu5dFx7J79HnegnD95EZ3Xxjte2Da5/10fu5tNphQk67Y9ljCGcr50bodLVyg\nLzSzfaP3dFm0bAsL162Koxj+SjhTeifa9kPALha6rRYCZxKaj37R2yfvxN0ulA8T4RSvHNgkZd7Z\nhF/7vVPmHRvNW0poI6+Ytqhku/cR2mVTy86Llu0GvEr4kC0mfKHPzrDu4mh/4wk9LpqklOlLOPUt\nqmT/o0hpt0xblt7m/lCWx+pvhJrYj/xcs/tVynIHOkV/twCej8qOieJ/NVPZ6PXbwOCU15cB/4r+\n3jza15LoWJ1IWjt1lvEPjtZLne5LWb4E6JsW4y+mtO1dTDhLmUs4K2qa9n4WR8fpDqBBJTG1S38v\nhGa0MYQfgu+Au4CNKll/Gr9sc29PuJYzP9r3/4DO0bJtCD8giwk/ymOA36VtrzXhrKtTpv2lHcvU\nNvcmhFry0ijmK4muIaV8z+ZHx3gicPoatv2Lzy4wiPAj9hPhWkXq97Jbyn7nE75XvdK2dzDhGsaP\n0ba71UQuyeVk0RuTBDKzvwFz3f2OuGPJhpldDWzmKdcdRGTdKLlLbKKmmGLgM0L/+JGE2tgzsQYm\nkgCZ7nwTyZVGhFPxzQldLP9J1Xf0ikgWVHMXEUkg9ZYREUmg2Jplmjdv7u3atYtr9yIieenDDz+c\n5+4tqioXW3Jv164dY8dmGq5CREQqY2bfZFNOzTIiIgmk5C4ikkBK7iIiCaTkLiKSQEruIiIJpOQu\nIpJASu4iIgmk5C4ikkBK7iIiCZT3o0KOGpV5fr9+uYxCRKR2Uc1dRCSBlNxFRBJIyV1EJIGU3EVE\nEkjJXUQkgZTcRUQSSMldRCSBlNxFRBJIyV1EJIGU3EVEEijvhx+oTOqwBNkORbAu64iI1EaquYuI\nJJCSu4hIAim5i4gkkJK7iEgCKbmLiCSQkruISAIpuYuIJJCSu4hIAim5i4gkkJK7iEgCKbmLiCSQ\nkruISAJlldzNbB8zm2RmU8xsWIblW5jZ62Y2zsw+NbP9qj9UERHJVpXJ3cwKgeHAvkBX4Agz65pW\n7G/AE+6+A3A4cGt1ByoiItnLpua+EzDF3ae6ewnwGDAwrYwDG0d/NwZmVV+IIiKytrIZz701MCPl\n9Uxg57QyFwEvmdkfgQbAgGqJTkRE1kk2NXfLMM/TXh8B3OfubYD9gAfNbLVtm9lQMxtrZmPnzp27\n9tGKiEhWsknuM4G2Ka/bsHqzy/HAEwDu/h5QD2ieviF3v9Pde7l7rxYtWqxbxCIiUqVskvsYoLOZ\ntTezYsIF0xFpZaYD/QHMbBtCclfVXEQkJlUmd3dfCZwGvAh8QegVM97MLjGzg6JiZwEnmNknwKPA\nYHdPb7oREZEcyeoB2e4+EhiZNu+ClL8nALtWb2giIrKudIeqiEgCZVVzz3ejRlW+rF+/XEUhIpI7\nqrmLiCSQkruISAIpuYuIJJCSu4hIAim5i4gkkJK7iEgCKbmLiCSQkruISAIpuYuIJJCSu4hIAim5\ni4gkUF6OLbOmsWJqYlupZTQWjYjkA9XcRUQSSMldRCSBlNxFRBJIyV1EJIGU3EVEEkjJXUQkgZTc\nRUQSSMldRCSBlNxFRBJIyV1EJIHycviB2k7DFYhI3BJdc687Zzr1p08C97hDERHJqWTW3MvL6Xzz\n6bR4/XHK69antHELJp5zL0s7bhd3ZCIiOZHImvvmz95Go4mj+eCRqbz/2DfMPPh0tj+rP03HvBR3\naCIiOZG45F5Qspx2D1zMpLPvoWyjRmDGnH2OZfylT7PN5UfR+LO34w5RRKTGJS65N3/jKZZ02oGl\n7bv/Yv6ibX/NF+c9RNeLfk/dOdNjik5EJDcSl9xbvvEkc/Y6JuOyhTvtzcxBZ9H977+joGR5jiMT\nEcmdRCV3W1lKk49HsaDXXpWWmTHoLJa3ak+HO4flMDIRkdxKVHJv9MVolm3egdKmLSsvZMaks+6k\n+Vv/odno/+UuOBGRHEpWcp80lh+32aXKcis3bsbEcx+gyzX/R9HC73MQmYhIbiUruU8Zx5LOO2RV\n9oce/fhu72Ppcu0Q3eQkIomTqOTe8MuPWNx5x6zLTxt8MXW/n0GrkXfXYFQiIrmX1R2qZrYPcCNQ\nCPzL3a/KUGYQcBHgwCfufmQ1xll1jGUr2WjGJH5q3y3rdbyomC/Oe4gef+rHwh57sLx1x1XLUseH\nSVVd48Zo/BkRqUlV1tzNrBAYDuwLdAWOMLOuaWU6A+cCu7p7N+DMGoh1jerOmU5J000pL663Vuv9\n1L4b048+n22uPAYrW1lD0YmI5FY2zTI7AVPcfaq7lwCPAQPTypwADHf3hQDunvOrlPVnfcWyzTtW\nXTCDmQefTnlxfdo+ek01RyUiEo9skntrYEbK65nRvFRbAVuZ2Ttm9n7UjLMaMxtqZmPNbOzcuXPX\nLeJK1J/1FcvXMblTUMDEYffR5qkbaDj5o2qNS0QkDtkkd8swL717SR2gM9APOAL4l5k1WW0l9zvd\nvZe792rRosXaxrpG9daj5g6womVbppx6A9tccTQFK5ZVY2QiIrmXTXKfCbRNed0GmJWhzLPuXuru\nXwOTCMk+Z+rP/pplrTqs1za+738ESztsR4e7zq2mqERE4pFNch8DdDaz9mZWDBwOjEgr8wywB4CZ\nNSc000ytzkCrUjx/NiXNN1+/jZgx+cxbaf7mUzT98JXqCUxEJAZVJnd3XwmcBrwIfAE84e7jzewS\nMzsoKvYiMN/MJgCvA3919/k1FXQmxQtmU7JJq/XezsqNmzHp7HvocvVx1Fm8sBoiExHJvaz6ubv7\nSGBk2rwLUv524M/RlHvuFC/4jpJmm1XL5hb22pN5v/4dnW84lS/+/ki1bFNEJJcScYdq4dIf8YJC\nyuo3rLZtTh16FQ2njKPlq49W2zZFRHIlEcm9bjU1yaQqr7cRE897kE63nEHduTOrddsiIjUtEQ/I\nrs4mmVSLu/Ti24NPZ+urBvPJP16Cgl/+FlY2RIGISNwSUXMvXvAdJU2rP7kDTD9yGAXLl9L66Vtq\nZPsiIjUhEcm9zqL5lDZuXiPb9sI6TDzvQdo9cAkbTZtQI/sQEaluiUjuRUsWsrJR0xrb/rLWnZg6\n5Aq2ueJorLSkxvYjIlJdEpHc6yyu2eQOMPuAE1jRcgs63RpPb08RkbWRmORe2rBmkztmTBx2H03H\nvEjLlx+u2X2JiKynZCT3Gm6WqbCyYRPGX/wUnYafSYOpn9X4/kRE1lUiknvRjwtyktwBlnbcjq9O\nvZ5uFx5C4ZJFOdmniMjaSkRyr7NkIaWNmuVsf3P2PJqFPfdk66sHQ3l5zvYrIpKtZCT3HFxQTTfl\n1OspXvg97e67KKf7FRHJRiKSe1EMyd2Livn80qfZ9OUHafmKBhcTkdol75O7rSzFSldU66Bh2Spt\n2pLPLx9Bp+FnsvH493K+fxGRyuT92DJ1li6irMHGYJmeBljzlnbYloln30u3Cw/ho1veY8VmW1Za\ntrKxaFLn9+tX9XwRkarkfc29cNmSWGrtqRb02Z8Zh/2Vbc87gDpLfog1FhERSEByL1i2lLJ68SZ3\ngJmHnskPO/yG7ucPpKBkedzhiMgGLu+Te6i5N4g7DDBjyqnXU7JJK7a57EgoK4s7IhHZgOV/cl++\nNPZmmVUKCvhi2P0U/rSYrW44BdzjjkhENlD5n9xrQZt7Ki+uy/hL/kOjyR/S7t4L4w5HRDZQyUju\n9WpBs0yKso0a8elVI2nxxpNw5ZVxhyMiG6D8T+61qVkmRWnTlnxy3Wtw771w7bVxhyMiG5i87+de\nuGwJ5bWs5l6hZJNW8NprsPvuUFQE258Rd0gisoFQzb2mtWkTEvwNN7D5M7fGHY2IbCASUXNf2aBx\n3GGs2ZZbwquv0nbX/hSULGPmoLPijkhEEi4RyX1F89Zxh1Gpn4cQ6EDdG99k+78MoM7SH5k2+KK1\nGjJhbYco0NAFIhu2vG+WKVi+tNb1lqnMipZtGXfjWzR/51k63nqW+sGLSI3J++Re2/q5V6W0aUs+\nvv51Np7wHl2uPQErWxl3SCKSQPmf3Fcso7xu/bjDWCsrGzXl02tfpu730+l2wcEULFsad0gikjB5\nn9wLSpbnXXIHKKvfkM+ufJ7SjTehx5/6UbRgTtwhiUiCJCO5F9WNO4x14nWKmHT2PczvcwA7ntaH\njaZPjDskEUmIvO8tU1CynPLienGHse7M+ObYC1mx6Zb0OHN32PKJcNOTiMh6yPuau5WuyO/kHvlu\nn8F8cf7DMGgQ3HqretKIyHrJ++Se9zX3FAt7DoB33w3JfehQrGRF3CGJSJ5KRnLP0zb3jDp2hPfe\ng3nz6PHnPSiePzvuiEQkDyUjuSek5r5Ko0bw1FMs6L03O568E40/ezvuiEQkz2SV3M1sHzObZGZT\nzGzYGsodamZuZr2qL8Q1K0hIm/tqCgr45tgLmfyn2+l24SG0ffRqKC+POyoRyRNV9pYxs0JgOLAn\nMBMYY2Yj3H1CWrlGwOnABzURaGWSVnNPHRMGYEGf/fnw9rF0veQwmnzyBhPPfYBRo5rXyL5qYmwa\njXEjEo9sau47AVPcfaq7lwCPAQMzlLsUuAZYXo3xrVl5OVa2Eq9TlLNdxmFFy7Z8fMMbLG3XnZ5D\nd1QzjYhUKZvk3hqYkfJ6ZjRvFTPbAWjr7s+taUNmNtTMxprZ2Llz5651sOkKSleEi6lrMbpivvI6\nRUw96Rq+PPNWul14KO3uuQBbWRp3WCJSS2WT3DNlzlWdsM2sALgeqHKQcne/0917uXuvFi1aZB9l\nJZLWJJON+X0OYOxd42g0eSw7nqq7WkUks2yS+0ygbcrrNsCslNeNgO7AKDObBuwCjMjFRdWCkuV4\nkrpBZqlkk1Z8duXzzN5/CDv88de0fvoW3fQkIr+QTXIfA3Q2s/ZmVgwcDoyoWOjui9y9ubu3c/d2\nwPvAQe4+tkYiTrEh1txXMWPWQSfx0S3vsulLD7Ld2fvAN9/EHZWI1BJVJnd3XwmcBrwIfAE84e7j\nzewSMzuopgNck8R2g1wLy9puxbhb3uGH7XeHnj3h5puhrCzusEQkZln1c3f3ke6+lbt3dPfLo3kX\nuPuIDGX75aLWDht4zT2FF9Zh+tHnwTvvwBNPQN++bDRtQtUrikhi5fUdqokbemB9dekCb7wBxxzD\nDmfsxpb3X4yVlsQdlYjEIK+Tu5WsoLxYyf0XCgrg5JNDj5pJY+l9/Hbw8stxRyUiOZbfyb1sJV6Y\n7BuY1tWKlm35/PIRfHXiNXDiiXDooTB9etxhiUiO5PXDOgrKSvHCvH4L1Sp9OAHMmL/rQbzZa0+2\nePRqWm+7A0XnnAVnnQV1Vz/jWW39DPMrG0JgfYYZ0BAFItUvv2vuK0sTP/RAdSivW59pgy/iw9vG\nwPvvw7bbwrPPqm+8SILlfXIvV3LP2vLNO8CIEXDTTXD++fQ4sx+NJo6JOywRqQH5ndzV5r5u9tkH\nPv6YOXsdQ/e/DWSby46i7ne6AUokSfI7ua9Um/s6q1OH2fsPYfSDk1nWpjO9TtyRDnecQ53FC+OO\nTESqQV4n94Iytbmvr7L6DZk2+CLG3P0ZdRYvYOejO7PlA5dSuPTHuEMTkfWQ18ldF1SrT0nzzZn8\nl7v46Nb3qT9zMjsf05m2j/2DguU/xR2aiKyD/E7uZSt1QbWaLWvdiYnnPcjH/3yNRhNHs/PRnWj9\nn5uxkhVxhyYiayG/k7va3GvMT+27MeGiJ/nsyudpOvYldjmqI22evB6WLo07NBHJQn4n97JS9Zap\nYUs678DnV/yXzy4fQePP34H27eGyy+CHH+IOTUTWIK+Te4Ha3HNmyVY7Mv7if4eByb78Ejp2hHPP\nhe+/jzs0Eckgr5P7hvBw7Fpnm23g/vth7FhYtAi23hpOO436306JOzIRSZHXDdZqc197lY0fs/br\nt4dBt1Lc/++0fvoWdji1Dz9235UZv/8z7N63yoeWr28cIrJm+V1z1/ADsSvZpBVfD7mc9x+dxoLe\ne9Pl2iHQuzc88giUlsYdnsgGK7+Tuy6o1hrl9Rswa+DJjL5/Ilx4Idx1F3ToQNtHr6Zo0by4wxPZ\n4OR5clebe61TUAAHHgivvw7PPkuDaRPY+ahObH3lsfDBBxqJUiRH8jq5F6jNvXbbcUcmnns/Hzw8\nhaXtu8MRR4Qmm3vu0Z2vIjUsr5O7hh/ID6WNmzPj8L/ClClw6aXw9NP0OWwLOg7/M/VnTI47PJFE\nyvvkXq429/xRUAD77gv//S8f3jGW8uJ67HBGX3qcsTubvviAavMi1Si/k7va3PPW8s3a8fUJV/De\n4zOYecgZtHz9cfoMahOe9zp6tNrmRdZTXjdYm56hmve8qJh5ux3MvN0Opnjut/zqy/vhyCOhfn34\nv/+DY46B5s3jDlMk7+R1zV3DDyRLSYvWcN55MHky3HwzjBsHnTrBIYfA00/DCo1MKZKtvE7uuqCa\nUAUF0K8fPPAATJsG++0HN94IrVvDSSfB229DeXncUYrUanndphGeoZrXbyGR1mdogdXXbQIdj4eL\njqfunOls+uojbHr0iRSs+InvBxzFnD2P5qcttl5tO/36rft+13Zdkdoor2vueDleUBh3FJIjKzbd\ngulHDmPMvZ8z/tKnKShZzvZ//g07ntSbNv++geK538YdokitkdfVXisvwy2/f59kHZixpFMPlnTq\nwVdDr6bpx6/T8pWH6f3gpSzdsitz+w2CrQ6BzTePO1KR2OR5ci+HQtXcN2iFhSzsOYCFPQcwubSE\nph++QstRT0D3C2HbbWHQoHBBdrPN4o5UJKfyu9rr5aq5yypeVMyCXfZj4rD7YPZs+Mtf4P33wxj0\ne+wBt90Gc+bEHaZITuR1ZrSystCzQiRd3bphALMHHwyJ/swzQy+bLl1gt93g+uvh66/jjlKkxuR3\nZtQFVclGvXowcCA8/DB89x2ccw5MmAC77AI9esDFF9Pgq091V6wkSt63uatZRtZKvXqw//5huv12\nePddeOYZuv/9t2DGvF1/C4W/hV/9StdzJK/ld2b0cjXLyLorLIS+feGf/+SDh7/i80ueZmWDjeGP\nfww9bYYMgREjYOnSuCMVWWt5nRnVFVKqjRlLO27HN8deCB9/HC7Edu8ON90ErVqF0SyHDw93zIrk\ngawyo5ntY2aTzGyKmQ3LsPzPZjbBzD41s1fNbMvqDzVDXOoKKTWlfftwEfaVV2DGDDj+eBgzBnba\nKST9YcPgrbdg5cq4IxXJqMrkbmaFwHBgX6ArcISZdU0rNg7o5e7bAf8GrqnuQDNSV0jJhcaN4dBD\n4b77wgXZu++GoiI44wzYdFM46qjwQPAFC+KOVGSVbC6o7gRMcfepAGb2GDAQmFBRwN1fTyn/PnB0\ndQZZGXWFlMpU7/g2qQqAnaH/zvS79FKYORNGjmTe8MdpcsJJLO2wHY0H7Q177w09ezLqrZ/PLDVm\njeRSNsm9NTAj5fVMYOc1lD8/CJIbAAAMMUlEQVQeeCHTAjMbCgwF2GKLLbIMcQ3UFVLi1qYNDB3K\n51sNpWDFMhp/+hbbz/4fHHcczJnDNtvvycLee7Og995Aq7ijlQ1INtVeyzAvY4dgMzsa6AX8I9Ny\nd7/T3Xu5e68WLVpkH2VlgakrpNQi5XXrs7D3XnDddTB+PIwbxw879qfZ+8/T+7husP32cPbZ8Oqr\nGptealw2NfeZQNuU122AWemFzGwAcD6wu7vn5JNr5WqWkVqsbVtm7z+E2fsPwcpWsnv90fDii3D+\n+eEmqt12C803++wTHkpimepRIusmm8w4BuhsZu3NrBg4HBiRWsDMdgDuAA5y9++rP8xKqFlG8oQX\n1gk3Rl18cehm+fXX4RGC48bB7rtDhw5wwgnw2GPwfe6+QpJcVSZ3d18JnAa8CHwBPOHu483sEjM7\nKCr2D6Ah8KSZfWxmIyrZXLWyct3EJHlqk03gsMPgnnvg22/huedCF8tHHoHOncOwCGedBS+8AEuW\nxB2t5KGshh9w95HAyLR5F6T8PaCa48qObmKSJDCDbt3CdMYZUFoa+tS/+ipcdRX8/vfQsyf07w8D\nBkDv3qErpsga5PfYMmqWkSQqKgpNOL/6Ffz972H4g7feCjdUnXJKaNLZbbeQ6AcMgK5d1V4vq8nr\n5I6aZWRD0KBBuOi6zz7h9dy58NproWZ/ww2wfHmo1e+xR5jat1eyl/xO7hpbRjZILVqE9vrDDguv\np04NtfpXXgk9cerW/TnR77EHVMc9JZJ38jozamwZEUJPm6FD4dFHw4NJXnwRdt4Znn8eevWCjh3D\n2DgPPRQu3soGIa9r7hpbJr+sz5AAtdHavp/KyqcOS7DeZcxg663DdPLJ4QEk48fD66/D00+HwdCa\nNQs1+n79wr96vmwi5XVy101MIlUwC10su3cP49SXl8Nnn4Vk//jjcOqpYfCziiacfv1Cs4/kvbxO\n7pSrt4zIWikoCMMgbL99qMWXlcEnn4Rkf//94UaqNm1+Tva77x765EveyevkbmqWEVk/hYWw445h\nOuusMD79uHEh2d91FwweHNr0K5L9brtBkyZxRy1ZyO/MqK6QItWrTp1wk9TZZ4e7Y+fPh9tuC001\nN98MbduGi7R//SuMHAmLF8cdsVQirzOjukKK1LCiIujTB847D15+GebNg+uvh0aN4B//CI8g3GUX\nOPdceOklPW+2FsnvZhl1hRTJrbp1w0PF+/aFCy4IN1C9/35oxrnsMvjoozAuTkUzTp8+UL9+3FFv\nkPI6uasrpEjM6tULPWz69QsjXv70E7z7bkj2f/sbfPppaMapSPY77xx+IKTG5XVyV1dIkVpmo41+\nHvMGwoiWb78dkv1f/gJffBESfEUf+969obg41pCTKs+Tu7pCitRqDRv+clycRYvCIGivvx763U+Z\nEppuKmr2PXuGi7qy3vL3KHr0pD8NkCSSPxo3hgMOCBPAggXw5psh2Z9wAnzzDfz61z8n+x49dF1t\nHeVtcrfyMlxNMiL5rVkz+O1vwwRhxMs33gjJ/p574LvvwsXb/v1hr72gSxdV6LKUt8kdd37ceqe4\noxBZTXWNObM++61sLJrU+bVSixZw6KFhgpDcR40KI15ee22Yt9desOeeIeE3bx5bqLVd3lZ9vU4R\n44a/F3cYIlKTNtsMDj8c/vWv0GTz8sth6ISHHw6jXfbuHfrgjxoFJSVxR1ur5G/NXUQ2LGahWaZL\nl3AxtqQk9LF/+WU455zQE6dv31Cz32uvMDLmBtyEk7c1dxHZwBUXh7FuLr0UPvgApk0LY+GMHx96\n57RrF0a9fOEFWLYs5mBzT8ldRJKhWbPwMPE77wyJfuRI2HJLuPLKMKzxQQeFZRvIA0uU3EUkecyg\nW7cwANqbb4Zkf/jhoSfOdtvBDjuEh4+PHv1zt+qEUXIXkeRr1gyOPDJciJ0zB266CUpL4Q9/CM+Y\nPfPMcHNVWVnckVYbJXcR2bDUqRMuvF51FUycGJ4526wZnHZaeFDJKafAa6+Fse3zmJK7iGzYunYN\nI1x+8klowtlii9Ccs/nmIeHnadONkruISIXOnWHYMBg7NnSz3HRTOOqo0K3y8stDX/s8oeQuIpJJ\nhw7houvkyeH5st9+GwY2+81v4MknQ5t9LWYe0+lGr169fOzYseu0bnXdri0ia1brhyvItRUr4Nln\nYfjwMKLlCSfA0KGhCacK1TUMhJl96O69qiqnmruISLbq1oVBg0KXyhdfDAOdde8OxxwDEybEHd0v\nKLmLiKyL7t1DDf7rr8NF2T32CAOeffxx3JEBSu4iIuuncePwgPCpU2HXXWHffeG442DWrFjDUnIX\nEakODRrAn/4EkyaFXjbbbgvXXQfl5bGEo+QuIlKdNt443CA1ejQ8/XToXTN9es7DUHIXEakJHTuG\nLjJ77w19+tBo0rr1DlxXSu4iIjWlsDC0xw8fzrbD9qPhpA9ztmsldxGRmvbb3zL5z3ew7d8GUmfR\n/JzsMqvkbmb7mNkkM5tiZsMyLK9rZo9Hyz8ws3bVHaiISD6b1/d3zO37O9rdf3FO9ldlcjezQmA4\nsC/QFTjCzLqmFTseWOjunYDrgaurO1ARkXw344hz2PSVh3IydEE2NfedgCnuPtXdS4DHgIFpZQYC\n90d//xvob7YBP7xQRCSDFS3aMOvAE+HHH2t8X9k8ILs1MCPl9Uxg58rKuPtKM1sEbALMSy1kZkOB\nodHLJWY2aQ37bZ6+fi2i2NZNbY4Nand8im3d1M7YHrkK1j22LbMplE1yz1QDTx9tLJsyuPudwJ1Z\n7BMzG5vN4DhxUGzrpjbHBrU7PsW2bjbk2LJplpkJtE153QZIv692VRkzqwM0BhZUR4AiIrL2sknu\nY4DOZtbezIqBw4ERaWVGAMdGfx8KvOZxjSUsIiJVN8tEbeinAS8ChcA97j7ezC4Bxrr7COBu4EEz\nm0KosR9eDbFl1XwTE8W2bmpzbFC741Ns62aDjS22h3WIiEjN0R2qIiIJpOQuIpJAsSf32jy0QRax\nDTazuWb2cTQNyWFs95jZ92b2eSXLzcxuimL/1Mx2rEWx9TOzRSnH7YIcxdXWzF43sy/MbLyZnZGh\nTCzHLcvYYjlu0b7rmdloM/skim+1e+jj+q5mGVuc39VCMxtnZs9lWFZzx8zdY5sIF2i/AjoAxcAn\nQNe0MqcAt0d/Hw48XotiGwzcEtOx2w3YEfi8kuX7AS8Q7kHYBfigFsXWD3guhmPWCtgx+rsRMDnD\n/2ksxy3L2GI5btG+DWgY/V0EfADsklYmru9qNrHF+V39M/BIpv+7mjxmcdfca/PQBtnEFht3f5M1\n30swEHjAg/eBJmbWqpbEFgt3n+3uH0V/Lwa+INxdnSqW45ZlbLGJjseS6GVRNKX3xojlu5plbLEw\nszbA/sC/KilSY8cs7uSeaWiD9A/0L4Y2ACqGNqgNsQEcEp2+/9vM2mZYHpds449Ln+g0+gUz65br\nnUenvzsQanmpYj9ua4gNYjxuUfPCx8D3wMvuXumxy/F3NZvYIJ7v6g3A2UBlz9qrsWMWd3KvtqEN\nakA2+/0v0M7dtwNe4edf4NogruOWjY+ALd19e+Bm4Jlc7tzMGgJPAWe6e/oITrEetypii/W4uXuZ\nu/cg3KW+k5l1TysS27HLIracf1fN7ADge3df0xM6auyYxZ3ca/PQBlXG5u7z3X1F9PIuoGcO4spW\nNsc2Fu7+Y8VptLuPBIrMrHku9m1mRYTk+bC7/ydDkdiOW1WxxXnc0uL4ARgF7JO2KPZhSCqLLabv\n6q7AQWY2jdCs+xszeyitTI0ds7iTe20e2qDK2NLaYg8itJPWFiOAP0S9P3YBFrn77LiDAjCzzSra\nFc1sJ8LnsMYfTxPt827gC3e/rpJisRy3bGKL67hF+2thZk2iv+sDA4CJacVi+a5mE1sc31V3P9fd\n27h7O0L+eM3dj04rVmPHLJtRIWuMxze0QXXFdrqZHQSsjGIbnIvYAMzsUULvieZmNhO4kHAhCXe/\nHRhJ6PkxBfgJOK4WxXYocLKZrQSWAYfn6Ad7V+AY4LOofRbgPGCLlNjiOm7ZxBbXcYPQm+d+Cw/v\nKQCecPfnasN3NcvYYvuupsvVMdPwAyIiCRR3s4yIiNQAJXcRkQRSchcRSSAldxGRBFJyFxFJICV3\nEZEEUnIXEUmg/wdmUhKet+qAlgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "R = np.random.normal(0.0, 1.0/np.sqrt(1000), (1000,500))\n",
    "D = np.zeros_like(R)\n",
    "D[0]= np.random.normal(0.0, 1.0/np.sqrt(1000), (500))\n",
    "W = R+2.0*D\n",
    "evals, _ = RMT_Util.eigenspectrum(W)\n",
    "RMT_Util.plot_ESD_and_fit(eigenvalues=evals, Q=Q)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-30T01:28:28.196818Z",
     "start_time": "2018-06-30T01:28:28.191719Z"
    }
   },
   "source": [
    "### Generate ensemble (500) of Matrices w/(the same, simple) Rank 1 perturbation\n",
    "\n",
    "$$\\mathbf{W}=\\mathbf{W}^{rand}+2\\Delta_{0}$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\\Delta_{0}=S_{min}=Q^{-1/4}$$\n",
    "\n",
    "is the minimum value needed to generate a spike\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-30T07:15:41.208952Z",
     "start_time": "2018-06-30T07:15:41.181728Z"
    }
   },
   "source": [
    "### Mean and Variance for Spiked Covariance\n",
    "\n",
    "#### Johnstone Conventions\n",
    "\n",
    "http://www.maths.usyd.edu.au/u/austms2013/talks/aust13Johnstone.pdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-29T16:54:13.998116Z",
     "start_time": "2018-08-29T16:54:13.993906Z"
    }
   },
   "outputs": [],
   "source": [
    "def mp_var_scale(Q):\n",
    "    g = 1.0/Q #gamma\n",
    "    numer = np.cbrt(np.cbrt(g))\n",
    "    denom = np.cbrt(1.0+np.sqrt(g))**4\n",
    "    return numer/denom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-29T16:54:14.004777Z",
     "start_time": "2018-08-29T16:54:14.000975Z"
    }
   },
   "outputs": [],
   "source": [
    "def rho(Q,S):\n",
    "    \"\"\"MP mean value of max eigenvalue\"\"\"\n",
    "    g = 1.0/Q\n",
    "    h = S*S\n",
    "    return (1.0+h)*(1.0+(g/h))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-29T16:54:14.012037Z",
     "start_time": "2018-08-29T16:54:14.006695Z"
    }
   },
   "outputs": [],
   "source": [
    "def tau_squared(Q,S):\n",
    "    \"\"\"MP variance of max eigenvalue - rho (mean)\"\"\"\n",
    "    h = S*S\n",
    "    h1 = (1.0+h)\n",
    "    h2 = h*h\n",
    "    g = 1.0/Q\n",
    "    return 2.0*h1*h1*(1.0-(g/(h2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### check that these methods match the Bouchaud equations\n",
    "\n",
    "ONly checked mean value so far...still do not understand variance scaling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-29T16:55:32.776489Z",
     "start_time": "2018-08-29T16:55:32.773907Z"
    }
   },
   "outputs": [],
   "source": [
    "#rho(Q,S), RMT_Util.spike_lmax(S,Q)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-30T02:35:38.847668Z",
     "start_time": "2018-06-30T02:35:38.844936Z"
    }
   },
   "source": [
    "#### Compute adjusted errors \n",
    "\n",
    "$$adjusted\\;errors\\;=\\;\\dfrac{\\lambda_{max}-\\rho(Q,S)}{\\tau(Q,S)}$$\n",
    "\n",
    "#### Check finite size scaling \n",
    "\n",
    "$$adjusted\\;errors\\times N^{1/2}\\rightarrow N(0,1)$$\n",
    "\n",
    "Where we have\n",
    "\n",
    "$$\\gamma=\\dfrac{1}{Q}$$\n",
    "\n",
    "and \n",
    "\n",
    "$$h=S=2\\times S_{min}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-29T17:00:24.202362Z",
     "start_time": "2018-08-29T16:55:32.952581Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "S_min = RMT_Util.spike_min(Q)\n",
    "evals = []\n",
    "max_evals = []\n",
    "max_errs = []\n",
    "adjusted_errs = []\n",
    "\n",
    "\n",
    "N, M = 1000, 500\n",
    "for ne in tqdm(range(1,1000)):\n",
    "    R = np.random.normal(0.0, 1.0/np.sqrt(N), (N,M))\n",
    "    D = np.zeros_like(R)\n",
    "\n",
    "    D[0,0]= 1.2*S_min\n",
    "    \n",
    "    W = R+D\n",
    "    Q = np.max(W.shape)/np.min(W.shape)\n",
    "    evals, _ = RMT_Util.eigenspectrum(W)\n",
    "    \n",
    "    #RMT_Util.plot_ESD_and_fit(eigenvalues=evals, Q=Q, sigma=1.0, alpha=0.25)\n",
    "    #sigma = RMT_Util.plot_ESD_and_fit(eigenvalues=evals, Q=Q, alpha=0.25, verbose=False)\n",
    "\n",
    "    S = np.max(np.abs(D))\n",
    "    emax = np.max(evals)\n",
    "    max_evals.append(emax)\n",
    "    max_errs.append(emax-RMT_Util.spike_lmax(S,Q))\n",
    "    adj_err= (emax-RMT_Util.spike_lmax(S,Q))/np.sqrt(tau_squared(Q,S))\n",
    "    adjusted_errs.append(adj_err)\n",
    "    \n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-29T17:00:24.744788Z",
     "start_time": "2018-08-29T17:00:24.204306Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'distrbution of adjusted errors in predicting spike at $2\\\\Delta=2Q^{-1/4}$  for $Q=2\\\\;\\\\;N,M=1000,500$')"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjEAAAEMCAYAAAA1Y+caAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xu0JVV94PHvj+aliAHkirya1gTJ\nCAE0HdQxGnwhIIqJREETMNFpzcjKY5KJmEwSFyEZnFm+CBkJEQI4EfCFooBKVGLMoNhgg6AiTdtI\n0yjNS3kopuE3f1RdLE7XefU9r337+1nrrHvOrl1Vv9qn6tzf2VX7VGQmkiRJpdlq2gFIkiRtDpMY\nSZJUJJMYSZJUJJMYSZJUJJMYSZJUJJMYSZJUJJMYSZJUJJMYSVpEIuLnIuKqiLg/Ig5olO8QERc0\nXh8XERumE6U0GiYxkrS4PAi8HPhoR/mLgC8ARMRWwDHArZMNTRqtzUpiIuKciDilfn5DRBw6imAi\nYm1EvGQUy6qXN7LYhlzvfhHx9Yi4LyJ+f4HLerSt69dT2aamzphKMgvt182ojqtZ28bOeEZ9nOux\nMvM/MrOth+Vw4DP189dRJTmPDLPsUX62SaOw4J6YzNw/M6/oV28SH1yd6xg0tjH4U+CKzNwxM08b\n5YJHsU1b8j+RKe4TQ1nIcTVr2ziNePrt4xGxXUScFRG31P+Qvx4RR7TUuyIi7omI7UYU16Dr3SEi\nTomIm+t634yINy9w9Usz83sRsQR4DXDhZixj5J9tY9rWgdu6y7w7R0RGxJUd5f8QEe8ZcN5bOsr3\niogHI+Lu4bdmk3WcGBErI+KhiDinZfouEXFRRDxQb//rBp3eb94eMV0RET+pT2PeHxE3jnudAFsP\nWnGcImLrzNw47ThGaB/ggr61tnBt7/uw+8Is7juzGJM2sTXVqZRfA74HHAl8OCJ+KTPXAkTEMuBA\n4FvAK4GPTGi9O1Od9vkm8GLgFuA/AxdFxMbMPCsinsKmp4sAjsnM73cWRsR+wHfql78FfDgzH4mI\nYePfrM+2bsfEINs67Loa+rZ1DwcD3weeERG7Z+btjfL/M8C8twE7R8QTM/NHdfnfAuvqaQu1HjgF\neBnwuJbpfw/8FNitjueSiLg2M28YYHq/eXs5MTM/0GXaeNaZmX0fwDOBa4D7qLL3C4BT6mlrgZfU\nz99G9QbdB9wIvLgu/yBVt+WPgfupsvm1df3rgIeodri1wNupduh7gH8Ctm/EkcAvNF6f04ij2zpe\n0qj/n4ArgHuBG4BXNqatBf6kjueH9XZu36U9ei3nC8DDwE/qOJ7eMv9JwM11O30T+PVB2rqlvbu2\nR7f3o62d6rp7AB8DNgDfBX5/0Jg6tq3rchrxt73vnWX93qvO+q37Xkt8zfYb+D1v1G/dP7vE1K8t\nBjqu6td7Ax+vl3UXcHqP93LgbQSeBXy9juEj9fRu723XNh6gbV7S5T34xbptju3XXoMeR93aZYDP\nueuAVzde/yVwFrAC+PQgy9icR8t6z6M67RMt23vNEMs9Bzigfv4HjTZ/J/C5eh0/BE4bcHmbfLYx\n5HHassyRbOvmtnWPen8EfIoqUXxzXbaE6nqjAwec90rgOY3jbHV9XLx7hNtzCnBOR9kOVAnB0xtl\nHwRO7Te937x9YrkCeFOXaWNZZ2b2T2KAbamy4z8CtqG6GOw/6PiwBfajynr3qMuXAT/fsUN3fpCt\novpwflyj7Pq6bBfg33nsP+V+/7Tb1jF/4G5T70R/Vm/Ti6g+/PZr1L2K6oN0F6pvX29paY+ey+n3\nZtbTf7Nez1bAa4EHgN37tXXLNvVK6rq+Hy3ttBVwNdWH9rbA04A1VFl+35gGWc4A7/ujZQO+V836\nPfe9jhib7TfQe94xb+v+2RJTz7bo164dcS4BrgXeQ3XAbw/8att7Ocw2NmL4gzqG36D6MGl7bwc5\nvnu1zSbHJdUH+/eAo/q11zDHUbd26fM5txvVP+dfbJStruPcpV72bl3m/TTVP/G2R8/kp3O9VD0d\nDwPP6rK9dw+4PZdSfVu/EngDVXK4XUu9lYO2UednG0Mepy3LGnpbR9nWfeqeB/w1Va/VZXXZ/lTJ\n2DYDzvsPjbb6Qr1NXwWOH9V20Z7EPBP4cUfZnwCf6je937wD7BsbgDupjv9DB4lpIevMzIGuiXkO\n1c763qwuGPso8LWWeg8D21F1v22TmWsz8+Y+yz4tM2/NzB83yk6vy+4G/gY4boAYB/Ec4AlU2d1P\nM/MLVDtOc/mnZeb6et2fourW2pzl9JSZH6nX80hmXgjcBBzC4G09iGHej18B5jLz5Hqb1gD/SPXt\neJiYei2nqe19b5YN+l7N19+cfa+5nH7veVOv/bMZU7+2GKZdD6H6Z/3fM/OBzPxJZn55wO3rtY3P\noeoxOq2O4eNUCU+bQdp4mGP3+cDFwAmZ+WkG33ce1eM4GkpEbAP8M3BuZn67LvtVqoTxi/X2fIHq\nYti2OI7KzJ26PI4aZr1USdOtmXlNyyx7Up2O6Cszj8zMPTLzuZl5DvCxzHyopd7yQZbXxbDHaaeh\nt3XEbd3LwVQJ2CXA8yNix7rs+sz8jwHnvQ7YPyJeSfXF4+PAL1H1fI5ku7p4AlUPW9MPgR0HmN5v\n3l7eRvXlY0/gTOBTEfHzY17nQEnMHsBtWadHtVs6K2XmauAPgXcAd0TEBRGxR59ltw3va5bdUq9/\nFPagOmCaV+PfQtXg85rnkh+katzNWU5PEXF8RKyKiHsj4l7gAGBXBmzrQQz5fuwD7DEfTx3Tn1F9\ncxkmpl7Laer3vg/Sxo/W38x9b94g73m3ODv3z+a0fm0xTLvuDdySm3+dTbdtbIuhdcjtgG08zLH7\nFuD/ZeYX69eD7juP6nEcDSyqocYfpOqBOrEx6QTgwsx8uH59fl02Ej3WO0f3ROXXqYdIDyszP7Q5\n8/Ux1HHaYizb2qlHW3ervx3VabJVmXkPVWJ/BFWPwapB56VKYg6mOmXy36h6M7ei6g0dp/uBJ3aU\nPZGql6zf9H7zdpWZX83M+zLzocw8l6o35shxrhMGS2JuB/aMx14BtrStYmZ+KDN/leoDKanOvT46\nuW2WlrK9O9azvvH6QeDxjddPGWB589YDe9c7dHP5w15ktaDlRMQ+VN8wTwSelJk7UXXDB0O0da1n\ne/R4Pzrb6Vbgux3Z/46ZeeSQMfVazmNCa5m3WTZIGz9mGX32vVHqtX92JgO92mLYdl0aEW0X4vfa\n5/tpi2HvbpUHaONebdPpLVTbND/SY9B9B+h7HMEA7VJv91lUidKr579hR8TjqEbvnN+o/kngFyLi\noJblXNYYkdH5uGzQ9da+C+zTse8TES8FlgPv6rddEzT0cdph6G0dcVt3cwDV9VRr6tefAF5FlcRs\n0ovSY95vAIcC12XmV+r5r2/7MjLsdvXxHWDriNi3UXYQ1TVL/ab3m3cYyc+Ox7Gtc5Ak5kpgI/D7\nEbF1RPwGLV22Uf1+wIvqTPQnVG/kw40qP6DqaurnrVENRduF6ptYcxjgKuB1EbEkIg6nuuq8qdc6\nvkp1XvtPI2KbqH634hUMf6X9QpezA9WbuwEgIn6HaseHAdu6oWt79Hk/OtvpKuBHEfG2iHhcvbwD\nIuJXhoyp13KGMVQbD7DvjVKv/bOpX1sM2663A6dGNRx1+4h4Xj1t0OOqzZVU7XRiHcPR3WIYsI0H\nbRuovmUdDrwgIk5l+H2n13EEg7XL+6m+Nb+i43THq4C7gWvrtt6+3tZLgeM7F5KZR2TmE7o82ob0\ndlsvVKcvAE6JiMdHNUz4t6gSqt/MzFn6cbqFfhYOva0jbutunglc2+ih/CRVj8JjemKi+l2nc7rN\nm5n3Ai/kZ70/B9OlJ2fY7aqP1+2prpdbUu+nW9fLeoDq1NXJ9efF84CjqXqjek7vN2+37Y6InSLi\nZfNxRMTrgRcAnx3FOnvKwS7mWs7PRjBcSGMEAz+7SO9Aqg+i+6g+AD5NfRFgXe9oqov47qW6aGct\n7Rckzo9wuBc4F3h8Rxw31Ov4INXOfsqg66C6MOtfqc63dY4K6qz7DuD/dmmPrsvJjovfusz/N3Ub\n3Qm8u17Wm/q1dWecvdqj1/vR2U512R71/N+nGl3ylY71dI2pY9u6LqetnXuUDfNe9dz3uq1rmPe8\n3/7ZZRv6tUXf46pRdynVN8K76v3mtB7v5cDbWMewiqpL9yNUHyZ/0bLt/Y7vgdumI75dqC5a/ut+\n7TXkcbRJu3TMO9+bND/SZv7xeqrRMtnl8X1aRtkM+ui13kadp1Ndu3QH1Sira4D9N3edo3zQ8dnG\nEMdpl+WNbVsHaesu850O/F1H2ao6vh0bZZ8H/ku/eRvT/oVqCPIotu0dLfvmOxrTd6H6vHigPg5e\n1zF/1+kDzNu23XNU1/TdVx9zXwFeOqp19npEvQAVIiK+B/xWZn5p2rFsaSJiLdUH+L9MO5ZxiYiv\nAmdk5j8NOd9aFnnbTENEvAZ4H9U/9gX/SNosK2lbI2JbquT7wBzsFNWiMIvbPRM/dqfBRMQcVca7\ndsqhaJGIiF+j+s2XO6l6IQ7kZz9NrynLzA9HxJOoTpUt6i8uJW1rZv6U6jTVFmUWt9skphD1tQGX\nU3VVfm/a8WjR2A/4MNWIpZupfvH19t6zaJIy8/3TjmFStqRt1Wh4OkmSJBVpwTeAlCRJmgaTGEmS\nVCSviVkEdt1111y2bNm0w5Ckolx99dV3ZubctOPQ5jOJWQSWLVvGypUrpx2GJBUlIjbrti6aHZ5O\nkiRJRTKJkSRJRTKJkSRJRTKJkSRJRTKJkSRJRTKJkSRJRTKJkSRJRTKJkSRJRTKJkSRJRfIXe8cg\nIs4GjgLuyMwD6rILgf3qKjsB92bmwS3zrgXuAx4GNmbm8okErUVh2UmXPPp87akvn2IkkjR+JjHj\ncQ5wOnDefEFmvnb+eUS8C/hhj/lfmJl3ji06SZIWAZOYMcjML0XEsrZpERHAa4AXTTImSZIWG6+J\nmbznAz/IzJu6TE/gcxFxdUSs6LaQiFgRESsjYuWGDRvGEqgkSbPMJGbyjgPO7zH9eZn5LOAI4K0R\n8YK2Spl5ZmYuz8zlc3PeSV6StOUxiZmgiNga+A3gwm51MnN9/fcO4CLgkMlEJ0lSWbwmZrJeAnw7\nM9e1TYyIHYCtMvO++vlhwMmTDFDq5IgnSbPKnpgxiIjzgSuB/SJiXUS8sZ50LB2nkiJij4i4tH65\nG/DliLgWuAq4JDM/M6m4JUkqiT0xY5CZx3Upf0NL2XrgyPr5GuCgsQYnSdIiYU+MJEkqkkmMJEkq\nkkmMJEkqkkmMJEkqkhf2SgUa1bDnWVuOJA3DnhhJklQkkxhJklQkkxhJklQkkxhJklQkkxhJklQk\nRydJW5jmSKJu5Y4wklQCe2IkSVKRTGIkSVKRTGIkSVKRTGIkSVKRTGIkSVKRTGIkSVKRTGIkSVKR\nTGIkSVKRTGIkSVKRTGLGICLOjog7IuL6Rtk7IuK2iFhVP47sMu/hEXFjRKyOiJMmF7UkSWUxiRmP\nc4DDW8rfk5kH149LOydGxBLg74EjgGcAx0XEM8YaqSRJhTKJGYPM/BJw92bMegiwOjPXZOZPgQuA\no0canCRJi4Q3gJysEyPieGAl8MeZeU/H9D2BWxuv1wHPbltQRKwAVgAsXbp0DKGqdN7QUdJiZ0/M\n5Lwf+HngYOB24F0tdaKlLNsWlplnZubyzFw+Nzc3uiglSSqEScyEZOYPMvPhzHwE+EeqU0ed1gF7\nN17vBayfRHySJJXGJGZCImL3xstfB65vqfY1YN+IeGpEbAscC1w8ifgkSSqN18SMQUScDxwK7BoR\n64C/Ag6NiIOpTg+tBd5c190D+EBmHpmZGyPiROCzwBLg7My8YQqbIEnSzDOJGYPMPK6l+KwuddcD\nRzZeXwpsMvxakiQ9lqeTJElSkeyJkSZoWsOem+uVpMXCnhhJklQkkxhJklQkkxhJklQkkxhJklQk\nkxhJklQkRydJ2iyjHPHkzSolbQ57YiRJUpFMYiRJUpFMYiRJUpFMYiRJUpFMYiRJUpEcnSRNySyP\nyPFeS5JKYE+MJEkqkkmMJEkqkkmMJEkqkkmMJEkqkkmMJEkqkkmMJEkqkknMGETE2RFxR0Rc3yj7\n3xHx7Yi4LiIuioidusy7NiK+ERGrImLl5KKWJKksJjHjcQ5weEfZ5cABmXkg8B3g7T3mf2FmHpyZ\ny8cUnyRJxTOJGYPM/BJwd0fZ5zJzY/3yK8BeEw9MkqRFxCRmOn4XuKzLtAQ+FxFXR8SKbguIiBUR\nsTIiVm7YsGEsQUqSNMtMYiYsIv4c2Aj8c5cqz8vMZwFHAG+NiBe0VcrMMzNzeWYun5ubG1O0kiTN\nLpOYCYqIE4CjgNdnZrbVycz19d87gIuAQyYXoSRJ5fAGkBMSEYcDbwN+LTMf7FJnB2CrzLyvfn4Y\ncPIEw9QMm4WbMo4yhoUsa5ZvnilpcuyJGYOIOB+4EtgvItZFxBuB04Edgcvr4dNn1HX3iIhL61l3\nA74cEdcCVwGXZOZnprAJkiTNPHtixiAzj2spPqtL3fXAkfXzNcBBYwxNkqRFw54YSZJUJJMYSZJU\nJJMYSZJUJJMYSZJUJC/slQo3C0OvJWka7ImRJElFMomRJElFMomRJElFMomRJElFMomRJElFcnSS\nNAYLuUFh6Tc37DVaqsTtkTS77ImRJElFMomRJElFMomRJElFMomRJElFMomRJElFMomRJElFcoi1\npInxZpWSRsmeGEmSVCSTGEmSVCSTmDGIiLMj4o6IuL5RtktEXB4RN9V/d+4y7wl1nZsi4oTJRS1J\nUllMYsbjHODwjrKTgM9n5r7A5+vXjxERuwB/BTwbOAT4q27JjiRJWzqTmDHIzC8Bd3cUHw2cWz8/\nF3hVy6wvAy7PzLsz8x7gcjZNhiRJEo5OmqTdMvN2gMy8PSKe3FJnT+DWxut1ddkmImIFsAJg6dKl\nIw5Vmp7Sb4ApaXLsiZkt0VKWbRUz88zMXJ6Zy+fm5sYcliRJs8ckZnJ+EBG7A9R/72ipsw7Yu/F6\nL2D9BGKTJKk4JjGTczEwP9roBOCTLXU+CxwWETvXF/QeVpdJkqQOJjFjEBHnA1cC+0XEuoh4I3Aq\n8NKIuAl4af2aiFgeER8AyMy7gb8GvlY/Tq7LJElSBy/sHYPMPK7LpBe31F0JvKnx+mzg7DGFJknS\nomFPjCRJKpI9MdKIjOPmhlv6DRMHGW7tkGxpy2VPjCRJKpJJjCRJKpJJjCRJKpJJjCRJKpJJjCRJ\nKpJJjCRJKpJJjCRJKpJJjCRJKpJJjCRJKpJJjCRJKpJJjCRJKpJJjCRJKpI3gJTGbJCbOG7pN3os\nnTehlKbDnhhJklQkkxhJklQkkxhJklQkkxhJklQkkxhJklQkRydNUETsB1zYKHoa8JeZ+d5GnUOB\nTwLfrYs+npknTyxIDcwRReVw9JC0OJnETFBm3ggcDBARS4DbgItaqv5bZh41ydgkSSqNp5Om58XA\nzZl5y7QDkSSpRCYx03MscH6Xac+NiGsj4rKI2H+SQUmSVAqTmCmIiG2BVwIfaZl8DbBPZh4E/B3w\niS7LWBERKyNi5YYNG8YXrCRJM8okZjqOAK7JzB90TsjMH2Xm/fXzS4FtImLXlnpnZubyzFw+Nzc3\n/oglSZoxJjHTcRxdTiVFxFMiIurnh1C9R3dNMDZJkorg6KQJi4jHAy8F3twoewtAZp4BHAP8XkRs\nBH4MHJuZOY1YtSmHVUvS7DCJmbDMfBB4UkfZGY3npwOnTzouSZJK4+kkSZJUJJMYSZJUJJMYSZJU\nJJMYSZJUJC/slfpwRNJs8H2Q1MmeGEmSVCSTGEmSVCSTGEmSVCSTGEmSVCSTGEmSVCSTGEmSVCSH\nWEu15hDetae+fIqRaHMN8h52q+P7L5XHnhhJklQkkxhJklQkkxhJklQkkxhJklQkkxhJklQkRydp\ni9btpoLebLB8434P3Uek6bMnRpIkFckkRpIkFckkZsIiYm1EfCMiVkXEypbpERGnRcTqiLguIp41\njTglSZp1XhMzHS/MzDu7TDsC2Ld+PBt4f/1XkiQ12BMze44GzsvKV4CdImL3aQclSdKsMYmZvAQ+\nFxFXR8SKlul7Arc2Xq+ryyRJUoOnkybveZm5PiKeDFweEd/OzC81pkfLPNlZUCdAKwCWLl06nkgl\nPYbDqqXZYk/MhGXm+vrvHcBFwCEdVdYBezde7wWsb1nOmZm5PDOXz83NjStcSZJmlknMBEXEDhGx\n4/xz4DDg+o5qFwPH16OUngP8MDNvn3CokiTNPE8nTdZuwEURAVXbfygzPxMRbwHIzDOAS4EjgdXA\ng8DvTClWSZJmmknMBGXmGuCglvIzGs8TeOsk45IkqUSeTpIkSUWyJ0aLVnMkydpTXz7FSDSrHG0k\nlc2eGEmSVCSTGEmSVCSTGEmSVCSTGEmSVCSTGEmSVCSTGEmSVCSHWGuL47BaSVoc7ImRJElFMomR\nJElFMomRJElFMomRJElFMomRJElFcnSSijTszR0dkaTN5b4jzS57YiRJUpFMYiRJUpFMYiRJUpFM\nYiRJUpFMYiRJUpEcnTRBEbE3cB7wFOAR4MzMfF9HnUOBTwLfrYs+npknTzLOWTLIKKRhRypJ09A5\nysl9VVo4k5jJ2gj8cWZeExE7AldHxOWZ+c2Oev+WmUdNIT5Jkorh6aQJyszbM/Oa+vl9wLeAPacb\nlSRJZTKJmZKIWAY8E/hqy+TnRsS1EXFZROw/0cAkSSqEp5OmICKeAHwM+MPM/FHH5GuAfTLz/og4\nEvgEsG/LMlYAKwCWLl065oglSZo99sRMWERsQ5XA/HNmfrxzemb+KDPvr59fCmwTEbu21DszM5dn\n5vK5ubmxxy1J0qwxiZmgiAjgLOBbmfnuLnWeUtcjIg6heo/umlyUkiSVwdNJk/U84LeBb0TEqrrs\nz4ClAJl5BnAM8HsRsRH4MXBsZuY0gp0Wb7inkg065H8cPw3gzw1oS2MSM0GZ+WUg+tQ5HTh9MhFJ\nklQuTydJkqQimcRIkqQimcRIkqQimcRIkqQieWGvFhVHNmkUZmU/GiSOQUYhOWpJi5U9MZIkqUgm\nMZIkqUgmMZIkqUgmMZIkqUgmMZIkqUgmMZIkqUgOsdbUOOxTi93mDNWeheHdwx6bHsuaFntiJElS\nkUxiJElSkUxiJElSkUxiJElSkUxiJElSkRydtIUbZFRB52iJQev1qz/IvNKWYCH7f7d5hy0f1aii\nhS5/2M8kR0Nt2eyJkSRJRTKJkSRJRTKJkSRJRTKJmbCIODwiboyI1RFxUsv07SLiwnr6VyNi2eSj\nlCRp9pnETFBELAH+HjgCeAZwXEQ8o6PaG4F7MvMXgPcA75xslJIklcEkZrIOAVZn5prM/ClwAXB0\nR52jgXPr5x8FXhwRMcEYJUkqQmTmtGPYYkTEMcDhmfmm+vVvA8/OzBMbda6v66yrX99c17mzY1kr\ngBX1y/2AG8cQ8q7AnX1rTZ9xjl4psRrnaJUSJ4wm1n0yc24UwWg6/J2YyWrrUenMIgepQ2aeCZw5\niqC6iYiVmbl8nOsYBeMcvVJiNc7RKiVOKCtWjY+nkyZrHbB34/VewPpudSJia+DngLsnEp0kSQUx\niZmsrwH7RsRTI2Jb4Fjg4o46FwMn1M+PAb6QnvOTJGkTnk6aoMzcGBEnAp8FlgBnZ+YNEXEysDIz\nLwbOAj4YEaupemCOnV7E4z1dNULGOXqlxGqco1VKnFBWrBoTL+yVJElF8nSSJEkqkkmMJEkqkknM\nFiYidomIyyPipvrvzl3qnVDXuSkiTqjLdoyIVY3HnRHx3nraGyJiQ2Pam6YZa11+RX2Lh/mYnlyX\nj/TWDgts08dHxCUR8e2IuCEiTm3UH0mbLuRWFxHx9rr8xoh42aDLnGScEfHSiLg6Ir5R/31RY57W\nfWCKsS6LiB834jmjMc8v19uwOiJOi1j4j1wuIM7Xdxzrj0TEwfW0kbfpAHG+ICKuiYiNUf3eVnNa\nt+N/5O2pGZSZPragB/C/gJPq5ycB72ypswuwpv67c/1855Z6VwMvqJ+/ATh9lmIFrgCWt8zzX4Ez\n6ufHAhdOK07g8cAL6zrbAv8GHDGqNqW6gPxm4Gn18q8FnjFIe1DdGuNaYDvgqfVylgyyzAnH+Uxg\nj/r5AcBtjXla94EpxroMuL7Lcq8Cnkv1W1GXze8H04izo84vAWvG1aYDxrkMOBA4Dzim33E1jvb0\nMZsPe2K2PM3bGpwLvKqlzsuAyzPz7sy8B7gcOLxZISL2BZ5M9U93pmPts9xR3Nphs+PMzAcz84sA\nWd2K4hqq3w8alYXc6uJo4ILMfCgzvwusrpc3yDInFmdmfj0z539v6QZg+4jYboHxjCXWbguMiN2B\nJ2bmlZmZVP+s2/ajacR5HHD+AmNZUJyZuTYzrwMe6Zi39bgaU3tqBpnEbHl2y8zbAeq/bV3BewK3\nNl6vq8uajqP61tYc3vbqiLguIj4aEXuzcKOI9Z/qLu+/aHw4PzpPZm4Efgg8acpxEhE7Aa8APt8o\nXmibDvJedmuPbvMOssxJxtn0auDrmflQo6xtH5hmrE+NiK9HxL9GxPMb9df1Weak45z3WjZNYkbZ\npgvZn3rto6NuT80gfydmEYqIfwGe0jLpzwddREtZ51j8Y4Hfbrz+FHB+Zj4UEW+h+nb3IvoYc6yv\nz8zbImJH4GN1vOf1mWcacc7/OvP5wGmZuaYu3qw2HWa9fep0K2/78rPQ32pY8C05ImJ/qru+H9aY\n3m0fmFastwNLM/OuiPhl4BN13EPvkwMYRZs+G3gwM69vTB91my5k24fdd7XImMQsQpn5km7TIuIH\nEbF7Zt5ed7ne0VJtHXBo4/VeVOfB55dxELB1Zl7dWOddjfr/SPXPZKqxZuZt9d/7IuJDVN3W5/Gz\nWzusiwFv7TDuNqX64a6bMvO9jXVuVpu2rHfQW110tkevefstc5JxEhF7ARcBx2fmzfMz9NgHphJr\n3XP5UB3T1VHd4PXpdf3macSpt2ntWDp6YcbQpoPE2WveQzvmvYLxtKdmkKeTtjzN2xqcAHyypc5n\ngcMiYueoRtocVpfN2+Qcef3Hieb7AAABgElEQVTPe94rgW9NM9aI2Doidq1j2wY4Cpj/NjnqWzss\nqE0j4hSqfx5/2JxhRG26kFtdXAwcW49geSqwL9XFkoMsc2Jx1qfhLgHenpn/Pl+5zz4wrVjnImJJ\nHdPTqNp0TX0a8r6IeE59euZ42vejicRZx7cV8JtU16hQl42jTReyP7UeV2NqT82iaV9Z7GOyD6rz\n3Z8Hbqr/7lKXLwc+0Kj3u1QXcq4GfqdjGWuAX+wo+59UF1VeC3yxc/qkYwV2oBo9dV0d1/uAJfW0\n7YGP1PWvAp42xTj3ourm/hawqn68aZRtChwJfIdqBMif12UnA6/s1x5Up8tuBm6kMbqjbZkjeL83\nK07gfwAPNNpvFdV1SV33gSnG+urGe3oN8IrGMpdTJQQ3A6dT/6L6NOKspx0KfKVjeWNp0wHi/BWq\n3pUHgLuAG3odV+NqTx+z9/C2A5IkqUieTpIkSUUyiZEkSUUyiZEkSUUyiZEkSUUyiZEkSUUyiZEk\nSUUyiZEkSUX6/zSZbhtO7uqPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(adjusted_errs,bins=100, density=True);\n",
    "plt.title(r\"distrbution of adjusted errors in predicting spike at $2\\Delta=2Q^{-1/4}$  for $Q=2\\;\\;N,M=1000,500$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Check that finite size scaling is $N^{1/2}$ and NOT $M^{1/2}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-29T17:00:24.752394Z",
     "start_time": "2018-08-29T17:00:24.746913Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.027236485069126987, 0.044721359549995794, 0.60902632082726)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.std(adjusted_errs), np.sqrt(1/M),np.std(adjusted_errs)* np.sqrt(M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-29T17:00:24.760310Z",
     "start_time": "2018-08-29T17:00:24.754771Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.027236485069126987, 0.03162277660168379, 0.8612932827560988)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.std(adjusted_errs), np.sqrt(1/N), np.std(adjusted_errs)* np.sqrt(N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bouchaud Scaling - at the EDGE he has $M^{-2/3}$ \n",
    "\n",
    "### BUT for the spikes it is $T^{-1/2}$\n",
    "\n",
    "see:  https://arxiv.org/pdf/0910.1205.pdf\n",
    "\n",
    "equations (66) for MP,  and (72) for spikes, with the text just below\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  },
  "widgets": {
   "state": {
    "19da624521684c0fa4c779bd52c45781": {
     "views": [
      {
       "cell_index": 15
      }
     ]
    }
   },
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}