sparse_categorical_speed.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: pylab import has clobbered these variables: ['e', 'f']\n",
      "`%matplotlib` prevents importing * from pylab and numpy\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pysal\n",
    "import scipy.sparse as sp\n",
    "import itertools as iter\n",
    "from scipy.stats import f, chisqprob\n",
    "import numpy.linalg as la\n",
    "import pandas as pd\n",
    "from datetime import datetime as dt\n",
    "import matplotlib.pyplot as plt\n",
    "%pylab inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 282,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#OLD\n",
    "\"\"\"\n",
    "def spcategorical2(n_cat_ids):\n",
    "    '''\n",
    "    Returns a dummy matrix given an array of categorical variables.\n",
    "    Parameters\n",
    "    ----------\n",
    "    n_cat_ids    : array\n",
    "                   A 1d vector of the categorical labels for n observations.\n",
    "\n",
    "    Returns\n",
    "    --------\n",
    "    dummy        : array\n",
    "                   A sparse matrix of dummy (indicator/binary) variables for the\n",
    "                   categorical data.  \n",
    "\n",
    "    '''\n",
    "    if np.squeeze(n_cat_ids).ndim == 1:\n",
    "        cat_set = np.unique(n_cat_ids)\n",
    "        n = len(n_cat_ids)\n",
    "        C = len(cat_set)\n",
    "        row_map = dict((id, np.where(cat_set == id)[0]) for id in n_cat_ids)\n",
    "        indices = np.array([row_map[row] for row in n_cat_ids]).flatten()\n",
    "        indptr = np.zeros((n + 1, ), dtype=int)\n",
    "        indptr[:-1] = list(np.arange(n))\n",
    "        indptr[-1] = n \n",
    "        return sp.csr_matrix((np.ones(n), indices, indptr))\n",
    "    else:\n",
    "        raise IndexError(\"The index %s is not understood\" % col)\n",
    "\"\"\"\n",
    "def spcategorical2(n_cat_ids):\n",
    "    '''\n",
    "    Returns a dummy matrix given an array of categorical variables.\n",
    "    Parameters\n",
    "    ----------\n",
    "    n_cat_ids    : array\n",
    "                   A 1d vector of the categorical labels for n observations.\n",
    "\n",
    "    Returns\n",
    "    --------\n",
    "    dummy        : array\n",
    "                   A sparse matrix of dummy (indicator/binary) variables for the\n",
    "                   categorical data.  \n",
    "\n",
    "    '''\n",
    "    if np.squeeze(n_cat_ids).ndim == 1:\n",
    "        cat_set = np.unique(n_cat_ids)\n",
    "        n = len(n_cat_ids)\n",
    "        C = len(cat_set)\n",
    "        indices = n_cat_ids\n",
    "        indptr = np.arange(n+1, dtype=int) \n",
    "        return sp.csr_matrix((np.ones(n), indices, indptr))\n",
    "    else:\n",
    "        raise IndexError(\"The index %s is not understood\" % col)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 283,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def spcategorical1(data):\n",
    "    '''\n",
    "    Returns a dummy matrix given an array of categorical variables.\n",
    "    Parameters\n",
    "    ----------\n",
    "    data : array\n",
    "        A 1d vector of the categorical variable.\n",
    "\n",
    "    Returns\n",
    "    --------\n",
    "    dummy_matrix\n",
    "        A sparse matrix of dummy (indicator/binary) variables for the\n",
    "        categorical data.  \n",
    "\n",
    "    '''\n",
    "    if np.squeeze(data).ndim == 1:\n",
    "        tmp_arr = np.unique(data)\n",
    "        tmp_dummy = sp.csr_matrix((0, len(data)))\n",
    "        for each in tmp_arr[:, None]:\n",
    "            row = sp.csr_matrix((each == data).astype(float))\n",
    "            tmp_dummy = sp.vstack([tmp_dummy, row])\n",
    "        tmp_dummy = tmp_dummy.T\n",
    "        return tmp_dummy\n",
    "    else:\n",
    "        raise IndexError(\"The index %s is not understood\" % col)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 284,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def spcategorical1a(data):\n",
    "    '''\n",
    "    Returns a dummy matrix given an array of categorical variables.\n",
    "    Parameters\n",
    "    ----------\n",
    "    data : array\n",
    "        A 1d vector of the categorical variable.\n",
    "\n",
    "    Returns\n",
    "    --------\n",
    "    dummy_matrix\n",
    "        A sparse matrix of dummy (indicator/binary) variables for the\n",
    "        categorical data.  \n",
    "\n",
    "    '''\n",
    "    if np.squeeze(data).ndim == 1:\n",
    "        tmp_arr = np.unique(data)\n",
    "        n = len(data)\n",
    "        C = len(tmp_arr)\n",
    "        tmp_dummy = sp.dok_matrix((n, C))\n",
    "        for each in tmp_arr[:, None]:\n",
    "            row = (each == data).astype(float)\n",
    "            tmp_dummy[:,each[0]] = row.reshape((n,1))\n",
    "        return tmp_dummy.tocsr()\n",
    "    else:\n",
    "        raise IndexError(\"The index %s is not understood\" % col)\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 285,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def spcategorical1b(data):\n",
    "    '''\n",
    "    Returns a dummy matrix given an array of categorical variables.\n",
    "    Parameters\n",
    "    ----------\n",
    "    data : array\n",
    "        A 1d vector of the categorical variable.\n",
    "\n",
    "    Returns\n",
    "    --------\n",
    "    dummy_matrix\n",
    "        A sparse matrix of dummy (indicator/binary) variables for the\n",
    "        categorical data.  \n",
    "\n",
    "    '''\n",
    "    if np.squeeze(data).ndim == 1:\n",
    "        tmp_arr = np.unique(data)\n",
    "        n = len(data)\n",
    "        C = len(tmp_arr)\n",
    "        tmp_dummy = sp.lil_matrix((n, C))\n",
    "        for each in tmp_arr[:, None]:\n",
    "            row = (each == data).astype(float)\n",
    "            tmp_dummy[:,each[0]] = row.reshape((n,1))\n",
    "        return tmp_dummy.tocsr()\n",
    "    else:\n",
    "        raise IndexError(\"The index %s is not understood\" % col)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 286,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n",
      "True\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "n = 20\n",
    "o = np.tile(np.arange(n),n)\n",
    "print np.allclose(spcategorical1(o).toarray(), spcategorical2(o).toarray())\n",
    "print np.allclose(spcategorical1(o).toarray(), spcategorical1a(o).toarray())\n",
    "print np.allclose(spcategorical1(o).toarray(), spcategorical1b(o).toarray())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 287,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "spcat1 = []\n",
    "for n in np.arange(25,250,25):\n",
    "    o = np.tile(np.arange(n),n)\n",
    "    s = dt.now()\n",
    "    a = spcategorical1(np.array(o))\n",
    "    e = dt.now()\n",
    "    spcat1.append((e-s).total_seconds())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 293,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "spcat1a = []\n",
    "for n in np.arange(25,250,25):\n",
    "    o = np.tile(np.arange(n),n)\n",
    "    s = dt.now()\n",
    "    b = spcategorical1a(np.array(o))\n",
    "    e = dt.now()\n",
    "    spcat1a.append((e-s).total_seconds())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 180,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "spcat1b = []\n",
    "for n in np.arange(25,250,25):\n",
    "    o = np.tile(np.arange(n),n)\n",
    "    s = dt.now()\n",
    "    b = spcategorical1b(np.array(o))\n",
    "    e = dt.now()\n",
    "    spcat1b.append((e-s).total_seconds())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 288,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "spcat2 = []\n",
    "for n in np.arange(25,250,25):\n",
    "    o = np.tile(np.arange(n),n)\n",
    "    s = dt.now()\n",
    "    b = spcategorical2(np.array(o))\n",
    "    e = dt.now()\n",
    "    spcat2.append((e-s).total_seconds())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 289,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.00925,\n",
       " 0.016085,\n",
       " 0.023427,\n",
       " 0.036465,\n",
       " 0.053077,\n",
       " 0.071263,\n",
       " 0.105652,\n",
       " 0.125046,\n",
       " 0.158332]"
      ]
     },
     "execution_count": 289,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "spcat1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 294,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.017515,\n",
       " 0.118629,\n",
       " 0.371812,\n",
       " 0.926212,\n",
       " 2.062251,\n",
       " 3.920805,\n",
       " 6.131568,\n",
       " 9.242453,\n",
       " 13.739116]"
      ]
     },
     "execution_count": 294,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "spcat1a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 295,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.005511,\n",
       " 0.014289,\n",
       " 0.040431,\n",
       " 0.117726,\n",
       " 0.188659,\n",
       " 0.357239,\n",
       " 0.518341,\n",
       " 0.888993,\n",
       " 1.277546]"
      ]
     },
     "execution_count": 295,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "spcat1b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 292,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.00032,\n",
       " 0.000211,\n",
       " 0.000295,\n",
       " 0.000455,\n",
       " 0.000561,\n",
       " 0.000798,\n",
       " 0.001032,\n",
       " 0.001201,\n",
       " 0.001746]"
      ]
     },
     "execution_count": 292,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "spcat2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 296,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x1132d46d0>"
      ]
     },
     "execution_count": 296,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEZCAYAAAB1mUk3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4FFXa9/FvpwOEkBACQRZJCOAg44KKwKOIEtl0FARx\nARTEbRgdRRx0RtFRMo6PjgsuzzjozCsogguOqIgjDLKERUCUHQQRVHbZ95CQpOv941QnndAJTZLu\n6uX3ua6+urq6uuruSuXcVedUnwMiIiIiIiIiIiIiIiIiIiIiIiIiIiIip5QD3FXJz94L7AIOA6nV\nFZBElTXAFU4HIRLOOgMLgYPAPmAB0N6BOOYAd1biczWAXOC8Cpa5C1iHSRa/AP8BkiqxrVD6GfO9\nDgMHgK+A3wEuB2OqrGygADji83g4iNt7G/hrENcvEnXqYpJAf0whkwD0AM53IJbKJoN0wAO4y3m/\nCyYBXGC/TgUGE5xkUF4MlfET0NWeTgZ6Az8C46pxG6EyCngnhNt7GyUDkdPSHnPWWZ7bMWekf8ck\njXWUFFAAKcBYYAewDfMPGOfz/p3Ad8B+YDqQ4fNeD2C9vd6/U3E1US3gFWC7/XgZqAm0Bo5iksER\nYKafzz4MfFLBd3wbeAOYgTkLzykT56vAFuAQ8C3mSsorG/gImGC/fyfQ0V7uECYJjfZZ/hLMVdgB\nYAUmUZXHNxl4dQCKgHPs1zmU3me3A/N9XnswVWg/2N/tKaAVsAiz3z/AXFkBZGH+hn8EdmP+pn2B\na4ANmKvGR+1lGwPHgPo+22pnf85fQszG7KNTzc+0Y/YeQzl2zAvs+P8LNPBZ3ntVewDzNxoC/BY4\nAeRjjokp9rI/A93s6fKOJyjZDyMwVY87MPvV6xpgrR3PNuAhP99LJOIkA3sxBeLVnFzffjvm8n44\n5p/8ZkwhUs9+/xPgdaA20BD4Ghhqv9cHUwidjfnnfhyTWADSMP9M/ez1Pmhvp7wrg6cw//Rp9uMr\nex5Ac0oXIGV1xlS3ZAOXYQoCX2/bsXTGFAivULpAvRWzX+IwBcROSgqObEzBc539OgFT0N5qv04E\n/seePhOzr6+2X3e3X6eVE7e/ZACwGVNdBCdfTd3OycngE8xV0DmYAnI2ptCtiynUbrOXzcL8Df6M\n+Zvcbcf3LlDH/nwuZn+DqWq7x2dbL2MSpz/Z+E8Gozh1MvgBOAuzb+cAz9rvNcf83frb8dan5Orv\nLUqODy/f/VnR8ZSF2Q/Z9np/g0l8Kfb7OzHHEfa8i/x/ZZHI0wbzz7MV808wBTjDfu92zJmTr6+B\nQUAjIA/zT+o1EFPYAEyjdEEVh/mnysAUQAvLrHcr5SeDjZQUogA9Mf/ccHIB4s/VwGeYM8gjmLN1\n7/JvA+/5LFsHKMQU3v7sp6QaLRtTYPmaa88vW8g/wslVJdMpKYzLKi8ZLAJG2tOBJINLfV5/iznz\n93oRU4iDKQRzKWmTSLY/36HM572Jrz/mjB1MobmT8tuasjGJ6ID92A804dRXBnOAx3zevxdzXIHZ\nB5PL2d5bnFxN5Ls/KzqesjD7wfd42oW54gOTjIdikqlUUkX/rOKc9cAdmLr384CmmLNjr7LJYLO9\nTAamimEnJf/kb2CuEMCcub3q894+e/6ZmIJgW5n1bq0gxqb2dr222PMCNR1TiKVirlhux5z5Alhl\nYjmGKay8638YU9V1EPM9Uihd0Jf9Hndhqq/WAUuAa+35zYGbKNkfBzBnmI1P43sANLPjC9Qun+nj\nZV7nUbrtZB9mf3iX9fd57/JTMFcLmZgqP281WnkmYfZ/KuYsfmeA8f9SzvbTMW0olXGq42kfJil5\n5fps9wZMVdHPmBOBSyoZQ0yLdzoAOaXvgfGUVPXAyWfIzTEFwVbM2V4DSv/jeG3BnJ297+e9X2H+\nmb1cZV6XtQNT6KyzX2fY8ypjtv04t5xtJ2EKqx3A5Zgz6a6YKhUwBbHvHT0WpW0EbrGnb8C0KTTA\n7I8JlN63p6sDptDynpEfw1zJeJ1uYikb++nIA/6NuUpsQ8UNxBb+74I6iqlK8zqd+LdQcrbub3sV\nqcrx9C2mLcUNDAM+pHQbkwRAVwbh52xMPbi3wE/HVPUs8lnmDOABzFXATZh//C8wZ2wzgJcwVQpx\nmMZJ773cb2Au8b2NnSn257E/fy5wPeYk4QEqLgjex9Rle+t4n8R/HbQ/12GqNFIxBVJHTMPtYp9l\nrsGcpdfEJLBFmCuiZEyV0V77vSc5dfXAIEqujg5hCqYiYCLmjqCemIIkAVMlUV51FJQUoHWBXpj9\nMIGSxLQC0+5SG1OvHsjvNFzlTFfGO5iryuuo+O9R3nZWYI6XdMzxMdLPMuV99j1Mu8tNmGOoASVt\nBruAlhXEU9njqQamPSgF8zc9Yj/LaVIyCD9HMA2cX2PO0hYBqyh9h8TXmDP5PZiC8gZK7kC6DVNI\neu8Y+jclhfqnwHOYO1YOAauBq+z39mL+if9mT59FydmuP09jzshW2Y9v7XleFZ0JHsDcYbLBjmMC\n8DwlVywWpmAZhakeuAhToIOpXppuf/ZnTDXFljLbLbvtqzA/cDqCqY8fgLmC2oaponoMc9fNFsx+\nruj/YiqmkXQLpqAcjSl8vV7GNGDvwtSTTywTj7/9Uvb9ipY/1Rn2V5irwqVUXM3nbz+BuftrEuZv\n+g3m+1YUg+96tmCS+EOYv9tyoK393ljMScgB4GM/263K8TQI075wCHOVd2sFy4oDxmH+IVb7ee8h\nzAFb3897UrHbKd0gGY38NTZK4GZSud+HSAwL5pXBW5S+O8ArHdO4tdnPeyIQmb/oDRcdML8vmOR0\nIBJZgpkM5uP/x1MvAX8K4najXXmX99EkFr5jMIwHvsT8RuSYw7GIlJJJ6WqiPpTcQ/0TqiYSEQkL\noby1NBHTUNfDZ56qA0REwkAok0ErzJXCSvt1M8wdDx0xd3KULNiqlbVp06YQhiYiEhU2Ye4EPG2h\nvLV0Naa7hBb2YxslHWmVsmnTJizL0sOyGDVqlOMxhMtD+0L7Qvui4gfmpLtSgpkM3sf0ddMac7/z\nHWXeVwOhiEiYCGYyGIj5mX4tzO2kb5V5vyWn15+LiEjYySvMo9d7vcgtyHU6lCrRL5DDXFZWltMh\nhA3tixLaFyWc3hfvrHwHC4vEGomnXjiMhevdPJZd/yUiEraKPEW0+Ucbxl43liuaOz+cs8vlgkqW\n60oGIlJt6tevz4EDFQ3UJ9UhNTWV/ftPrmVXMhCRsOByudD/bvCVt5+rkgzUZiAiIkoGIiKiZCAi\nIigZiIgISgYiIkGRmZnJ7NmzS80bOnQobdq0we12M378eIci80/JQEQkCPzd8XPhhRcyZswY2rVr\n573zJ2woGYhIzHjuuedo1qwZdevWpU2bNsyePZvs7GxuvPFGBgwYQN26dbn44otZtWpV8We2bt1K\nv379OOOMM0hLS2PYsGGA6VCza9eupKWl0bBhQwYNGsShQ4cAGDx4MFu2bKF3794kJyfz4osvAvD7\n3/+erl27kpCQEPovfwpKBiISE77//nv+8Y9/8O2333L48GFmzJhBZmYmAJ999hk333wzBw4c4JZb\nbqFv374UFRVRVFREr169aNGiBZs3b2b79u0MGDCgeJ2PP/44O3fuZN26dWzdupXs7GwAJkyYQEZG\nBp9//jlHjhzh4YcfduAbnx4lAxEJKZer6o/KcLvd5Ofns3btWgoKCsjIyKBly5YAtG/fnn79+uF2\nuxkxYgR5eXksWrSIJUuWsHPnTl544QVq165NrVq1uOyyywBo1aoV3bp1o0aNGqSlpfGHP/yBuXPn\nVtduCjklAxEJKcuq+qMyzjrrLF555RWys7Np1KgRAwcOZOfOnQA0a9aseDmXy0WzZs3YsWMH27Zt\no3nz5sTFnVxU7tq1iwEDBtCsWTNSUlIYPHgw+/btq1xwYUDJQERixsCBA5k/fz6bN2/G5XLxyCOP\n4HK52Lp1a/EyHo+Hbdu2ceaZZ5Kens6WLVsoKio6aV2PPfYYbrebNWvWcOjQISZMmIDH4yl+P9wa\niE9FyUBEYsKGDRuYPXs2+fn51KpVi4SEBNxuNwBLly7lk08+obCwkFdeeYWEhAQuueQSOnToQJMm\nTXj00UfJzc0lLy+PhQsXAnD06FHq1KlD3bp12b59Oy+88EKp7TVq1Iiyw/cWFBSQl5eHx+PhxIkT\n5OXlhU1fTkoGIhIT8vPzGTlyJA0bNqRJkybs3buXZ555BoA+ffowadIk6tevz7vvvsvHH3+M2+3G\n7XYzdepUNm7cSEZGBunp6Xz44YcAjBo1imXLlpGSkkLv3r254YYbSl0NjBw5kqeffprU1FReeukl\nAHr06EFiYiKLFy9m6NChJCYmMn/+/NDvDD/C9TpGvZaKRKBI7LX0L3/5Cxs3bmTChAlOhxIw9Voq\nIlLNIi15BYuSgYjENJfLFXGNvcEQrntA1UQiESgSq4kikaqJREQkKJQMREQk6MlgHLALWO0z7wVg\nHbAS+BhICXIMIiJyCsFOBm8BV5eZNwM4F7gA2ACMDHIMIiJyCsFOBvOBA2XmfQl4f7P9NdAMERFx\nlNNtBncCXzgcg4hIzIt3cNuPAyeA9/y96e0XHCArK4usrKyQBCUiUh0yMzMZN24cXbt2LZ43dOhQ\n5s2bxw8//MC4ceMYMmRIlbaRk5NDTk5OFSM1nEoGtwPXAN3KW8A3GYiIRJryhr0cMGBAcW+pVVX2\nRPkvf/lLpdflRDXR1cAfgT5AngPbF5EYFc7DXi5ZsoRLL72U1NRUmjZtyrBhwygoKAjBXjGCnQze\nBxYCZwNbMW0EfweSMA3Jy4ExQY5BRCTsh72Mj4/n1VdfZd++fSxatIhZs2YxZkzoisdgVxMN9DNv\nXJC3KSJhzPWXqlePWKNOv8sL32EvGzRoQEZGRvF73mEvAUaMGMHo0aNZtGgRbre7eNhL72hnvsNe\ntmrVCqB42Munnnqq0t+pXbt2xdPNmzdn6NChzJ07l+HDh1d6nafDyQZkEYlBlSnIq4PvsJdr167l\nqquuKh5noLxhL10uV4XDXg4fPpwFCxZw5MgRPB4P9evXr3R8GzZsYMSIESxdupTc3FwKCwtp3759\npdd3upy+tVREJGTCedjLe++9l3POOYeNGzdy6NAh/vd//7fU+oJNyUBEYkK4D3t59OhRkpOTSUxM\nZP369bz++ush2CvhzxKRyBPO/7urVq2yOnbsaCUnJ1v169e3evfube3YscPKzs62brzxRqt///5W\ncnKy1a5dO2v58uXFn9uyZYvVt29fq0GDBlZaWpo1fPhwy7Isa+3atdbFF19sJSUlWRdddJE1evRo\nKz09vfhzU6ZMsTIyMqx69epZo0ePtizLsrp06WK5XC4rLi7OcrlclsvlsubOnWtZlmXNmzfPatOm\njZWUlGRdfvnl1pNPPmldfvnlfr9LefsZqHQdnMYzEJFqE4njGWjYS0PVRCIS0yIteQWLkoGIxDQN\ne2mE6x5QNZFIBIrEaqJIpGoiEREJCiUDERFRMhARESUDERFByUBERFAyEBEJC9nZ2QwePNix7SsZ\niIgEQWZmJrNnzy41b+jQobRp0wa328348eNLvef0bx2UDEREgqC8YS/HjBlDu3btTir8nf59hpKB\niMSMcB720uVykZeXV24cwaZkICIxIdyHvbQsiylTppwUR2FhYTB2x0mUDEQktFyuqj8qwXfYy4KC\nAjIyMmjZsiVQMuyl2+1mxIgR5OXlsWjRIpYsWVI87GXt2rWpVatWqWEvu3XrRo0aNYqHvZw7d26V\ndo2/OBYvXlyldQZKyUBEQsuyqv6oBN9hLxs1asTAgQPZuXMnUP6wl9u2batw2MsBAwbQrFkzUlJS\nGDx4MPv27avcPrH5i8MbY7ApGYhIzAjnYS8Bv3E0bdr0tNdTGUoGIhITwn3Yy4riCAUlAxGJCfn5\n+YwcOZKGDRvSpEkT9u7dyzPPPANAnz59mDRpEvXr1+fdd9/l448/xu1243a7mTp1Khs3biQjI4P0\n9HQ+/PBDAEaNGsWyZctISUmhd+/e3HDDDaWuBkaOHMnTTz9NamoqL730EgA9evQgMTGRxYsXM3To\nUBITE5k/fz5griT69u3rN45QCOavHMYB1wK7gfPtefWBSUBz4GfgZuCgn89qPAORCBSJ4xlo2Esj\nmFcGbwFXl5n3KPAl0BqYZb8WEXFMpCWvYAlmMpgPHCgz7zrA+xvs8UDfIG5fROSUNOylEew9kAlM\npaSa6ACQ6rPt/T6vfamaSCQCRWI1USQKRjVRfBVjqgrLfvjl/SUfQFZWFllZWcGPSEQkguTk5JCT\nk1Mt6wr1lcF6IAv4BWgCzAHa+PmcrgxEIpCuDEIj0hqQ/fkMGGJPDwE+DfH2RUTEj2BeGbwPdAHS\ngF3Ak8AU4EMgA91aKhJ1dGUQGsG4MgjXJnQlA5EIpGQQGtFQTSQiImFIyUBEJAjKDnu5YcMG+vTp\nwxlnnEGDBg24+uqr2bBhg4MRlqZkICISBGWrcg4dOkTfvn3ZsGEDu3btomPHjvTp08fBCEtTMhCR\nmOHksJcdOnTgjjvuoF69esTHx/Pggw/y/fffc+BA2Y4axJclIpEnnP93169fb6Wnp1s7d+60LMuy\nNm/ebG3atMkaNWqUVaNGDWvy5MlWYWGh9eKLL1otWrSwCgsLrcLCQqtt27bWiBEjrNzcXCsvL89a\nsGCBZVmWtXHjRmvmzJnWiRMnrD179lhXXHGF9eCDDxZvLzMz05o1a1a58XzyySdW06ZNK/VdytvP\nVPBD3lNx8hfIIhKDXNXwi1mrEj0S+A572aBBAzIyMorf8w43CTBixAhGjx7NokWLcLvdxcNeekc7\n8x32slWrVgDFw14+9dRTAcWybds27r///uKurcOBkoGIhFRlCvLq4Dvs5dq1a7nqqquKC+Pyhr10\nuVwVDns5fPhwFixYwJEjR/B4PNSvX/+UcezZs4eePXty33330b9//+r7glWkNgMRiRlOD3t54MAB\nevbsSd++fRk5cmRwvmQlKRmISExwetjLw4cPc9VVV9G5c+fiEdbCiZKBiMQEJ4e9HD16NJ9++inf\nfvstb731FsnJySQnJ1O3bl22bdvmyP4oS91RiEi1icTuKDTspaErAxGJaZGWvIJFyUBEYpqGvTTC\ndQ+omkgkAkViNVEkUjWRiIgEhZKBiIgoGYiIiLqjEJFqlJqaqsbYEEhNTa32dYbrX00NyCLC9sPb\nueifFzFj8AwubHyh0+GEPTUgi0jUsSyLu6fezf0d71ciCAElAxEJS2OXj2XPsT2M7BxeHbpFK1UT\niUjY2XxwM+3/X3tyhuRw7hnnOh1OxFA1kYhEDY/l4c7P7uThSx9WIgghp5LBSGAtsBp4D6jlUBwi\nEmZe/+Z1jp04xkOdHnI6lJgSSDJIAtz29NnAdUCNKmwzE/gt0A443173gCqsT0SixKb9mxiVM4rx\nfccTH6c730MpkGQwD3PmfibwX2Aw8HYVtnkYKAASMb9zSAS2V2F9IhIFPJaHO6bcweOXP87ZaWc7\nHU7MCSQZuIBcoB8wBrgJOK8K29wPjAa2ADuAg8DMKqxPRKLAq4tfBWD4JcMdjiQ2BXoddilwK3CX\n/boqbQ2tgAcx1UWHgH/b637Xd6Hs7Ozi6aysLLIcGkRbRILv+73f88yCZ1h812LiXLqvJVA5OTnk\n5ORUy7oCuQWpC/AQ8BXwHKYwHw48UMlt9gd6AHfbrwcDlwD3+SyjW0tFYkShp5DO4zpz2wW38fsO\nv3c6nIhWlVtLA7kymGs/vDZR+UQAsB54AqgN5AHdgSVVWJ+IRLAXF75InZp1uKf9PU6HEtMqSgZT\nfaYtSmcbC3NXUWWsBN4BvgU8wDLgX5Vcl4hEsDW71zB60Wi+/e23qh5yWEWXE1n28/VAY2CivfxA\nYBem3j9YVE0kEuUKigq4ZOwl3Nv+Xu5ud/epPyCnVJVqokA+tBS4OIB51UnJQCTKPTX3KRZvW8x/\nbvmPur2uJsFuM0jENBpvsl+3tOeJiFTK8p3LeW3Jayz/3XIlgjARSDL4AzAH+Ml+nQkMDVZAIhLd\n8gvzGfLpEEb3HM2Zdc90OhyxBZqSE4A2mIbj9UB+0CIyVE0kEqUen/U4a/as4dP+n+qqoJoFu5oI\nTD9CLezlL7DnvVOZDYpI7Ppm+ze8ufxNVt6zUokgzASSDCZi2glWAEU+85UMRCRgeYV5DPl0CP93\n9f/ROKmx0+FIGYGk5nXAOZgqolBRNZFIlPnjjD+y+dBmPrzpQ6dDiVrBriZaAzTBdConInLavtry\nFRNXT2T1vaudDkXKEUgyaAh8h+kywttwXJVfIItIDDl24hi3T7mdMdeMIS0xzelwpByBXE5k2c/e\nehuXPT3X79LVQ9VEIlHigWkPsP/4fib2m+h0KFEv2NVEOZjuKDpgksASYHdlNiYisSXn5xw+Xvex\nqociQCA9Q90MfI0Z1OZmTDK4KZhBiUjkO5J/hDun3Mm/ev+L1NqpTocjpxDI5cQqTDfT3quBhsAs\noG2wgkLVRCIR757P76GgqICxfcY6HUrMCHY1kQvY4/N6X2U3JiKxYcamGUzbOI1V96xyOhQJUCDJ\nYDrwX+A9TBLoD0wLZlAiErkO5h3k7s/uZux1Y0lJSHE6HAlQoGf4NwCX2dPzgU+CE04xVROJRKg7\nptxBgjuB13u97nQoMSfY1UQtgC+Ayfbr2pieS3+uzAZFJHp9vuFz5v48l1X3qnoo0gRyN9FHlO6T\nyGPPExEptv/4fn73+e94q89bJNVMcjocOU2BJAM3cMLndT5QIzjhiEikGjZtGDedcxNdMrs4HYpU\nQiDVRHuBPsAU+3Ufe56ICAAfr/uYb7Z/w4p7VjgdilRSIA0NZwHvAk3t19uAwcDGYAWFGpBFIsae\nY3to+0ZbJt88mU7pnZwOJ6ZVpQH5dD6UZC9/pDIbOk1KBiIRwLIsbvr3TbRMbcnzPZ53OpyYV5Vk\nEEibQWNgLKbR+AhmbIO7KrMxEYkuk9ZOYt3edTx15VNOhyJVFEgyeBuYQUk10Q/AH6q43XqY5LIO\n0z32JVVcn4iE2C9Hf2H49OGM7zuehPgEp8ORKgokGaQBkyi5vbQAKKzidl/F/Hbh15g+jtZVcX0i\nEkIey8Nvp/6Woe2G0r5pe6fDkWoQyN1ER4EGPq8vAQ5VYZspwOXAEPt1YRXXJyIhZFkWw74YxsG8\ngzzR5Qmnw5FqEkgyeAiYCrQEFmKuFKrShXULTMd3bwEXAEuB4UBuFdYpIiFgWRYjZ43k6+1fM+u2\nWdR013Q6JKkmFSWDjsBWTGHdBRiK6aPoS3t+VbbZDrgf+AZ4BXgUeNJ3oezs7OLprKwssrKyqrBJ\nEakOz8x/xnQ5cftcdUIXBnJycsjJyamWdVV0C9JyoBuwH7gC025wP3AR0Aa4sZLbbAwswlwhAHTG\nJINePsvo1lKRMPPK4lf4xzf/YN7t82iS3MTpcMSPYHVUF4dJBGC6rf4nprO6ycDKymzM9gvmyqI1\nsAEzcM7aKqxPRIJs7LKxvLz4ZSWCKFZRMnBj+iAqwBTYQwP8XCCGYX7VXBPYBNxRxfWJSJB8sOYD\nnsx5kjlD5tC8XnOnw5EgqahQfx+Yi+mHKBczjgHAr4CDVdzuSqBDFdchIkE29fupPDj9Qb4c/CWt\nG7R2OhwJolPVLV2KqeOfARyz57XGdE2xLIhxqc1AxGEzf5zJLZNv4T+3/IcOZ+rcLRKEqm+iUFIy\nEHHQV1u+ou+kvky+eTJXNL/C6XAkQMHum0hEYsiyncu4ftL1TLx+ohJBDFEyEJFi3+35jmvfu5Z/\n9vonV511ldPhSAgpGYgIAJv2b6LnhJ680OMFrv/19U6HIyGmZCAibDu8je4TuvPEFU8wqO0gp8MR\nBygZiMS4XUd30f2d7tzX4T5+1/53TocjDlEyEIlh+4/vp+fEnvQ/tz8Pd3rY6XDEQbq1VCRGHck/\nQvcJ3bks/TJG9xztvS1RIph+ZyAip+V4wXGuee8aWtdvzRu93lAiiBJKBiISsBNFJ+j7QV9Sa6fy\nTt93cMe5nQ5JqomSgYgEpNBTyMDJAykoKuDfN/2bGu4aTock1ShYXViLSBTxWB7u+uwuDuUdYurA\nqUoEUoqSgUgM8I5b/OOBH5l+63RqxddyOiQJM0oGIlGu7LjFdWrWcTokCUNKBiJRTuMWSyCUDESi\n2KuLX+XtlW8z7/Z5NEhs4HQ4EsaUDESi1NhlY3lp8Usat1gComQgEoU+WPMBT8x5gpzbczRusQRE\nyUAkykz9firDpw9n5uCZGrdYAqZkIBJFZv04i7s+u4v/3PIfzm90vtPhSARRMhCJEgu3LmTA5AFM\nvnmyBrCX06YurEWiwLKdy+j7QV+NWyyV5mQycAPLgakOxiAS8b7b8x3XvHsNb/R6Q+MWS6U5mQyG\nA98B6pFOpJJ8xy3u9+t+TocjEcypZNAMuAZ4k/DtOVUkrHnHLX788scZfMFgp8ORCOdUMngZ+CPg\ncWj7IhFt97HdxeMW39vhXqfDkSjgxN1EvYDdmPaCrPIWys7OLp7OysoiK6vcRUViyoHjB+g5QeMW\nC+Tk5JCTk1Mt63KiiuYZYDBQCCQAdYHJwG0+y2hwGxE/th7ayjXvXcPVra7m+R7Pa7hKKSWSRzrr\nAjwM9C4zX8lApIzVu1Zz7XvXMqzjMB7u9LASgZwk0kc6U6kvcgqzf5rNgI8G8OrVrzLw/IFOhyNR\nKFxPLXRlIGJ7b/V7PDj9QSbdOIkrW1zpdDgSxiL9ykBE/LAsixcWvsBrS15j1m2z1NeQBJWSgUgY\nKvIU8eD0B8nZnMPCuxbSrG4zp0OSKKdkIBJmjhcc59aPb+Vg3kHm3zGfegn1nA5JYoA6qhMJI/ty\n99F9QncS4hOYdus0JQIJGSUDkTDx04Gf6DSuE53TOzOx30RqxddyOiSJIUoGImFg6Y6ldH6rM8M6\nDuO5Hs8R59K/poSW2gxEHDZ943Ru++Q2/tnrn1z/6+udDkdilJKBiIPGLR/HY7Me49MBn9IpvZPT\n4UgMUzL/2IDQAAAQ+UlEQVQQcYBlWTw19ynGrxzP3Nvncnba2U6HJDFOyUAkxAo9hdz7+b0s+2UZ\nC+9aSOOkxk6HJKJkIBJKR08cpf9H/fFYHnKG5JBcK9npkEQA3U0kEjK7ju7iyvFX0qhOIz4b8JkS\ngYQVJQORENiwbwOdxnXi2l9dy9jrxlLDXcPpkERKUTWRSJAt2rqI6yddz9Ndn+budnc7HY6IX0oG\nIkE0Zf0U7p56N+P7jueaX13jdDgi5VIyEAmSMd+M4el5TzPt1mm0b9re6XBEKqRkIFLNPJaHx2Y9\nxsfrPmbBnQtomdrS6ZBETknJQKQanSg6wZ1T7mTTgU0svGshaYlpTockEhAlA5FqcijvEDd8eAPJ\ntZKZddssEmskOh2SSMB0a6lINdh+eDtXvH0FZzc4m49u+kiJQCKOkoFIFa3dvZZO4zpxy3m38No1\nr+GOczsdkshpUzWRSBXM/XkuN390M6N7jmZQ20FOhyNSaUoGIpU0ac0khk0bxvs3vE+3lt2cDkek\nSpxKBunAO8AZgAX8C/g/h2IROW0vLXqJlxe/zMzbZtK2UVunwxGpMpdD221sP1YAScBSoC+wzn7f\nsizLodBEylfkKeKhGQ/x5Y9fMu3WaWSkZDgdkkgxl8sFlSzXnboy+MV+ABzFJIGmlCQDkbCz/fB2\n7p92PweOH2DBHQtIrZ3qdEgi1SYc7ibKBC4CvnY4DhG/DuYd5LFZj9H2jba0rt+a/w76rxKBRB2n\nG5CTgI+A4ZgrhGLZ2dnF01lZWWRlZYUyLhHyC/MZ880Ynl3wLL1a92LF71aQnpLudFgixXJycsjJ\nyamWdTnVZgBQA/gcmAa8UuY9tRmIYzyWh/dWv8efZ/+Z8xudz7PdnuW8M85zOiyRU4rENgMXMBb4\njpMTgYgjLMtixqYZPDLzERLiExjfdzxdMrs4HZZISDh1ZdAZmAeswtxaCjASmG5P68pAQmrpjqX8\naeaf2HZ4G890fYZ+v+7nPcsSiRhVuTII16NdyUBCYtP+Tfx5zp+Z+/NcnuzyJHdddJeGpJSIVZVk\nEA53E4mE3O5ju3lg2gP8z5v/w7kNz2XDsA3c0/4eJQKJWU7fTSQSUkdPHOWlRS/x6tevMuj8Qay7\nbx0N6zR0OiwRxykZSEwoKCrgzWVv8td5fyUrM4sldy+hVf1WToclEjaUDCSqWZbF5HWTeWzWY2Sk\nZDB14FQubnqx02GJhB0lA4la8zbP409f/on8onxeu+Y1erTsoTuERMqhZCBRZ/Wu1YycNZK1e9by\n9JVPM/D8gcS5dK+EBNHBg1CvntNRVImSgUSNLYe2MCpnFF/88AUjO49k8s2TqRVfy+mwJJoUFcHG\njbByJaxYYZ5XroTjx2HnTqhZ0+kIK03JQCLegeMHeHbBs4xdPpZ729/Lhvs3kJKQ4nRYEumOHIHV\nq0sX/GvWQKNGcMEF5jF0qHlu3hwivApSyUAi1vGC47y25DWeX/g8/dr0Y/W9q2ma3NTpsCTSWBZs\n3Xry2f6OHXDuuaawv/BCGDwY2raFunWdjjgowjWV6RfIUq4iTxETVk3gyTlPcnHTi3m227O0SWvj\ndFgSCfLz4bvvTi74ExJKzvYvvNA8/+pXEB9Z58vqjkJigmVZfPHDFzw661FSaqXwfI/n6ZTeyemw\nJFzt2XNyob9xI7RqVbrQv+ACOOMMp6OtFkoGEtX25e5jzs9zeG3Ja+w+tpu/df8bvVv31m2iYhQV\nwQ8/lC74V6wwjbplC/1zzjFXAVFKyUCiSm5BLl9t+YqZP85k1k+z2LBvA5c3v5ybzrmJQW0HER8X\nWZfuUo127zaNur6P776Dxo1PLvgzMiK+Ufd0KRlIRCv0FLJ0x9Liwv+bHd9wYeML6daiG91bdqfj\nmR2p6Y7cW/akEnJzTSHvW+ivWgUFBXD++aUf550XtY26p0vJQCKKZVms37ueWT/NYuaPM5m7eS7p\nddPp3rI73Vt25/KMy0mulex0mBIKRUXw448nF/rbtkHr1qUL/bZtoWnTmDvbPx1KBhL2th/eXlz4\nz/ppFvFx8fRo2YNuLbrRtUVXGiU1cjpECbayVTyrVpmz/zPOOPlsv3VrqKHuxE+XkoGEnYN5B8n5\nOYdZP85i5k8z2X1sN11bdKV7i+50a9mNVqmt1AAcrcpW8axaZZ5PnDBn96riCRolA3FcXmEei7Yu\nYuaPM5n500y+2/MdndI7Fdf7X9j4QvUPFG3KVvF4C/2tW+Hss08+2z/zTFXxBJmSgYRckaeIFb+s\nKK72WbRtEec2PJfuLbvTrUU3Lk2/lIT46L2FL+pYljmj37vXPPbtq3ja+9yo0cln+6ricYySgQSd\nZVls3L+xuN5/zs9zaFSnUfGZf5fMLtRLiOxeG6OGZcGxY4EV6r7TcXGQlmYeDRqUfi5vXmKi099W\nfCgZSLXJL8xn6+GtbDm0hc0HN7P50GY2HdjEvM3zKPQUmjt+WnSna4uunFn3TKfDjQ3Hj5vCes+e\n0s8VFfDx8eUX5v4K9gYNVLBHASUDCdjh/MNsPrjZFPaHNhcX+N7pfcf30TS5Kc1TmtO8XnPznNKc\nyzIu4+wGZ6vRt6o8HtP3fdmC3fe57LyCAmjY0DzS0sxzgwYlz/4K+Nq1nf6m4oBITAZXA68AbuBN\n4Lky7ysZVIJlWew+tpvNhzaXOrP3LfRPFJ04qaBvXq85GSkZNE9pTtPkprjj3E5/lciRn19+Ie7v\nef9+SEoqKdS9z77TZZ+TktTwKgGJtGTgBr4HugPbgW+AgcA6n2WUDGw5OTlkZWUB5pe62w5vKy7Y\nyxb4Ww5toU6NOicV9M1T7MK+XnMa1G4QsWf3vvsiIJYFeXnmcfy4eXin/c0rb9p3Xm6uKdC9BXxe\nXsmZ+akKde+ZfDU0rp72vohi1bkvLAsKCyv3uPTSagmhSqqSDJzo5KUjsBH42X79AdCH0skgKlmW\nRV5hHkdPHOXoiaMcOXGkZDr/iN/5s9+aTcpPKWw+uJlfjv5Co6RGpQr59k3bc8M5N5CRkkFGSgZJ\nNZOc/prl83hKClTfwtXfs595OfPnk3XuuYEV2nl55qy9Zk3TMVnt2ubhnfY3r+z73uoW3/dr14b6\n9UsK97p1HTlrD7QA9HhMQVVQcHLh5W/e6S5bVGS2Ufa5vOnqnufxwLp1ObRsmVVujKfz8HhMc0tl\nHgsXRvYFnBPJ4Exgq8/rbcD/OBDHKRUUFfgtsAMpxMtbJj4unuRaySTVTCKpZhLJNX2mayWTVKNk\nOi0xjVaprXjwygfJSMmgWd1m1HDXMKcv3v+IoqKS/8oj+ViFuViFReZRUFg87V3WO794urBkvnc9\npeYVFWEVFprC9Vgu5B2H3OO4jnsL7Vxcx+3XeebZlXccV95x4o7n4soveXYVnMCqmYAnIRFPrdoU\n2c+emma6qGZtPLXMc5H9XFizNkU161BUM43dcY1YWf9KCmvUpjA+gcIatSmIr02B2552J5jX9rwC\ndwIe4ooLjdN+HAPPkZPn+xZE5c0LZJmqrGvXLpgw4dSFtmWZC5H4+JLnsg9/8wNZ1u02z3Fx5uF2\nl372N+12m/wcyGcqWo/v8wcfwJAhp449kEdcXGQX6FXhRDKIiPqfnk89Sm5yPcAFFrjsKy+XZR8p\nlsvMs+KAFFyk+Mw3kiwXyUATzGpc3r7WioDj4DpudoUHiyPAkeKtm/kuiti08RjPvr8EF1/jwoML\ncGEVL2EF/ChZvmQtLiyXy/98+wEuPC4XHlc8HncqnqQ0ipLceHDjccXhcZnpIpf9GreZ5yr9vmW/\nLv5Hc5Vcy9oh4PI33+f52C8/8NF57UrPt5d3uTy4XLlAbjnvn2JegMtb9rP3Ufbz5U57v1uZaX/L\nep/jXBAH1PCz7JGxu2h29ypccSXz4rzrjzs5xqoqtB95PvPCpSr3h5Rf2N9kZfkLeIOvQHV8k/+2\nbRuxVbDgTJvBJUA2phEZYCTgoXQj8kagVWjDEhGJeJuAs5wOIlDxmIAzgZrACuDXTgYkIiLO+A3m\njqKNmCsDERERERGR0q4G1gM/AI84HIsTfgZWAcuBJfa8+sCXwAZgBhCNnQCNA3YBq33mVfS9R2KO\nkfVAzxDFGCr+9kU25s675fbjNz7vRfO+SAfmAGuBNcAD9vxYPDbK2xfZROGx4cZUHWUCNYjN9oSf\nMAe6r+eBP9nTjwB/C2lEoXE5cBGlC8Dyvvc5mGOjBuZY2Yi58SZa+NsXo4ARfpaN9n3RGLjQnk7C\nVC//mtg8NsrbF9VybITbTvL9QVoBJT9IizVl7/K6DhhvT48H+oY2nJCYDxwoM6+8790HeB9zjPyM\nOWY6Bj/EkPG3L8D/3X/Rvi9+wRRoAEcxP049k9g8NsrbF1ANx0a4JQN/P0iLta4xLWAm8C3wW3te\nI0y1AfZzrIwRWd73boo5Nrxi5TgZBqwExlJSLRJL+yITc8X0NTo2MjH7YrH9usrHRrglg/D4FYuz\nLsP8kX8D3IepMvBlEZv76VTfO9r3yetAC0w1wU5gdAXLRuO+SAImA8Px/X2mEWvHRhLwEWZfHKWa\njo1wSwbbMY0kXumUzmyxYKf9vAf4BHNZtwtTXwjmB827HYjLCeV977LHSTN7XjTbTUmh9yYll/ux\nsC9qYBLBBOBTe16sHhvefTGRkn0RlcdGrP8gLRFItqfrAF9h7gB4npI7qx4lOhuQwfzdyzYg+/ve\n3oaxmpgzok2E79gclZVJ6X3RxGf6D8B79nS07wsX8A7wcpn5sXhslLcvovbYiOUfpLXA/PFWYG4d\n837/+ph2hGi+tfR9YAdwAtNudAcVf+/HMMfIeuCqkEYafGX3xZ2YQmAVpl74U0q3G0XzvuiM6a5m\nBSW3Tl5NbB4b/vbFb4jdY0NERERERERERERERERERERERERERESkrMcxv8dYibnnOtgdkuUAF5/G\n8pdg+o5ZDnyH6V0SoDex2T27iEi1uxRYiPl5PpgfITUpf/FqMQdodxrLfw+cb0+7iK1f00uECbe+\niUQC1RjYi+meF2A/Jf06PYEZGGg18E+fz+QALwHfYLr/7YDp/2kD8Fd7mUzMrzUnYs7m/w3U9rP9\nnphktBT4ENN9SFkNMd0Og+k3Zp09fTvwd3va99ekuZiOCetgBrj5GliG6a5ZRET8qIMpQL8H/gFc\n4fNeqs/0O0Ave3oO8Kw9/QCmy4dGmL5bttqfy8T85P9Se7mxwEM+n28HpAFzKUkSj2ASUFlPYJLU\nx8BQoJY9fwglycCrt73OeOAZ4FZ7fj37Oyb6Wb9ItdGVgUSqY5j6+6GYHl4nYQpZgK6YuvpV9vQ5\nPp/7zH5eYz92YfoA+pGSHh63Aovs6YmYPmG8XJi2gHMwVwbLgduADD8x/hVoj+k75xZgus86fP0K\n0/HazUAh5qrjUXvdczBJJB2RIIp3OgCRKvBgzqbnYqqEhmBGxxuDOYPfjmm0TfD5TL7PZ/N95nso\n+X/w7fPdhf8+4L/EFPCn8iPwBvD/MEmr7JCmSZhEdjclg7UA9MOMXSsSEroykEjVGnNG7XURZmi/\nBEzhvQ9T0N5UiXVnYM7+wRT4833eszBXHZcBrex5dcrE4nVtmXgLOXk4y3HAW5juyr3+S8lg52C+\nm0hQ6cpAIlUSpt69HqaQ/QFTZXQIcxa+BtN4+3U5n69odKzvMaPMjQPWYkaS8rUX0wj8PiXtAI9z\n8pn8IEyDda4d460+27UwSecGTCK50/7MXZjqpVcw1VxxmKsLNSKLiIRQJqUHlRGJCaomEjlZtI2Z\nKyIiIiIiIiIiIiIiIiIiIiIiIiIiIhKo/w9hiDZdsfygSwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x113aec490>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.arange(25, 250, 25)\n",
    "plt.plot(x, spcat1, x, spcat1a, x, spcat1b, x, spcat2)\n",
    "plt.legend(('spcat1', 'spcat1a', 'spcat1b', 'spcat2'))\n",
    "plt.title('Speed of Sparse Dummy Functions')\n",
    "plt.xlabel('Sample Size')\n",
    "plt.ylabel('Seconds')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### It is obvious that spcat1 and spcat2 are the fastest functions for creating sparse dummies. The differnece between these two functions is that spcat1 creates the dummies row by row (in sparse format) and then stacks the rows using sparse.hstack() while spcat2 builds an index for the non-zero dummy variable entries and then uses this index to instanitate the entire sparse dummy matrix, somewhat simularly to how the regimes are built in pysal. The slower functions spcat1a and spcat1b are riffs on spcat1 in that they also work row by row but they first instantiate an empty sparse matrix (either lil or dok) and then assign the rows to the empty matrix. The scipy doucmentation suggests that assignment for these two types of sparce matrices is faster than csr and csc and therefore is better for incrementally building a sparse matrix. However, it appears that avoiding assignment alother may be the best route if building a sparse matrix row by row. Therefrow csr is the best sparse data structure for row by row contruction in this context. Now lets lets row by row building (spcat1) vs entire instantiation using an non-zero-value index (spcat2) when the number of dummy variables starts to become much larger."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 311,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "spcat1 = []\n",
    "for n in np.arange(100,1000,100):\n",
    "    o = np.tile(np.arange(n),n)\n",
    "    s = dt.now()\n",
    "    b = spcategorical1(np.array(o))\n",
    "    e = dt.now()\n",
    "    spcat1.append((e-s).total_seconds())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 312,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "spcat2 = []\n",
    "for n in np.arange(100,1000,100):\n",
    "    o = np.tile(np.arange(n),n)\n",
    "    s = dt.now()\n",
    "    b = spcategorical2(np.array(o))\n",
    "    e = dt.now()\n",
    "    spcat2.append((e-s).total_seconds())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 314,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.039641,\n",
       " 0.127073,\n",
       " 0.313178,\n",
       " 0.692129,\n",
       " 1.286246,\n",
       " 2.14461,\n",
       " 3.737387,\n",
       " 5.416489,\n",
       " 7.583021]"
      ]
     },
     "execution_count": 314,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "spcat1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 315,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.001299,\n",
       " 0.00189,\n",
       " 0.004869,\n",
       " 0.009591,\n",
       " 0.010163,\n",
       " 0.012483,\n",
       " 0.0204,\n",
       " 0.026372,\n",
       " 0.032667]"
      ]
     },
     "execution_count": 315,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "spcat2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 316,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x124185d90>"
      ]
     },
     "execution_count": 316,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEZCAYAAACQK04eAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4U2Xax/FvKKvsOyKFIoq7IgKvjKgRFXFBGECBERF1\nZBw3GNRRXF46zuiMC4rjMo6OoqLijqivCzoQBMEFARFkEZR9XyxILdD2vH/cJ026p22SkzS/z3Xl\n6snJSc6d0+Q+T57zLCAiIiIiIiIiIiIiIiIiIiIiIiIiIiIeCwBXV/K5fwS2AnuAptEKSKqVJcAZ\nXgchkih6AXOBn4GdwBygmwdxzASuqsTzagHZwPFlbHM1sAw7MWwB/g9oUIl9xdMa7H3tAXYDnwN/\nAHwexlRZmcBBYG/Y7ZYY7u954K8xfH2RpNYIS/hDsIRSFzgXOMGDWCqb+NOBfCCtlMfPxJL9Se79\npsDlxCbxlxZDZfwE9HaXGwL9gB+B56K4j3gZD7wYx/09jxK/SKm6YaXJ0ozESpqPYSeIZYSSEUBj\n4FlgE7AB+7LVCHv8KuB7YBfwEdA+7LFzgeXu6z5G2VU9dYCJwEb39ghQG+gM/IIl/r3ApyU89xZg\nahnv8XngKWA6VroOFInzUWAdkAXMx34hBWUCbwKT3cevAnq422VhJ5wJYdufiv262g0swk5KpQlP\n/EHdgTzgWPd+gMLHbCQwO+x+PlYN9oP73u4BOgHzsOP+KvaLCcCP/Q9vBbZh/9MBwAXASuzX4O3u\ntm2AfUCzsH11dZ9X0skvEztG5a3PcGMOfoYCbsxz3Pg/BpqHbR/8tbob+x9dAVwDHAD2Y5+Jae62\na4Cz3eXSPk8QOg5jserDTdhxDboAWOrGswG4uYT3JZLQGgI7sOTXl+L14yOxn+ijsS/0pVjCaOI+\nPhX4F1APaAl8CYxyH+uPJZyjsC/yndhJBKAF9sUZ6L7uGHc/pZX478G+4C3c2+fuOoAOFE4WRfXC\nqkwygdOwL324591YemFf/okUTp6XYcelBpYMNhNKEplYkrnYvV8XS6qXufcPAf7HXT4MO9Z93fvn\nuPdblBJ3SYkfYC1W5QPFfyWNpHjin4r9ujkWS4YzsATbCEtgI9xt/dj/4C7sf/J7N76Xgfru87Ox\n4w1WXXZt2L4ewU6SJcmk5MQ/nvIT/w/AEdixnQn83X2sA/Z/G+LG24zQr7pJhD4fQeHHs6zPkx87\nDpnu656PneQau49vxj5HuOtOLvktiyS2o7EvynrsAz8NaOU+NhIrEYX7EhgOtAZysC9k0DAssQB8\nSOGkVAP7ArXHks3cIq+7ntIT/ypCCROgD/ZFhuLJoiR9gXexkuFerBQe3P554JWwbesDuViiLsku\nQlVhmVhyCjfLXV80od9G8eqOjwgl3qJKS/zzgHHuciSJv2fY/flYiT7oISxhgyW8bELXEBq6z+9e\n5PnBk9wQrCQOliA3U/q1oUzspLPbve0CDqX8Ev9M4I6wx/+Ifa7AjsFbpexvEsWresKPZ1mfJz92\nHMI/T1uxX3JgJ95R2IlTqqCsL6zE3nLgSqyu/HigLVbqDSqa+Ne627THqgk2E/pCP4WV/MFKZI+G\nPbbTXX8Y9qXfUOR115cRY1t3v0Hr3HWR+ghLWE2xXyIjsRItgFMkln1YYgq+/i1YddXP2PtoTOGk\nXvR9XI1VQS0DvgIudNd3AC4hdDx2YyXHNhV4HwDt3PgitTVs+dci93MofK1jJ3Y8gtuW9Pzg9tOw\nXwEZWLVdsCqsNK9hx78pVjrfHGH8W0rZfzp2zaMyyvs87cROQEHZYfsdhFX3rMFO+qdWMoaUV9Pr\nAKTACuAFQtU1ULzk2wH70q/HSnHNKfwlCVqHlbqmlPDYkdgXN8hX5H5Rm7AEs8y9395dVxkz3Ntx\npey7AZaYNgGnYyXk3li1CFjSDW9Z41DYKuB37vIg7BpAc+x4TKbwsa2o7liCCpa092G/UIIqehIp\nGntF5ABvYL/+jqbsi7cOJbdG+gWrDguqSPzrCJXCS9pfWaryeZqPXftIA24EXqfwNSGJkEr83jkK\nq7cOJvd0rLpmXtg2rYCbsNL9JdiX/AOsJDYdeBirFqiBXTgMtpV+CvuZHrwQ2dh9Pu7zjwN+i534\nb6LsL/0UrO45WCf7v5RcZ1ySi7FqiaZY8umBXVT9ImybC7DSd23sZDUP+6XTEKv22eE+9r+U/xN/\nOKFfPVlYEsoDXsJa5vTBkkZdrFqhtColCCXLRsBF2HGYTOgktAi7TlIPqwePpB+Er5TlyngR+7V4\nMWX/P0rbzyLs85KOfT7GlbBNac99BbtOcgn2GWpOqI5/K3B4GfFU9vNUC7t+0xj7n+51/0olKPF7\nZy928fFLrPQ1D1hM4ZYKX2Il9O1YUhxEqCXQCCwhBlvuvEEogb8D3I+1HMkCvgPOcx/bgX1h/+Eu\nH0GoFFuSv2ElrcXubb67LqisEt5urKXHSjeOycADhH6JOFgSGY/9xD8ZS95gVUQfuc9dg1U1rCuy\n36L7Pg/rLLQXqz8fiv0y2oBVM92BtX5Zhx3nsj7/72EXMNdhSXEClmiDHsEuLm/F6rVfKhJPScel\n6ONlbV9eyflz7NfeN5RdVVfScQJrhfUa9j/9Gnu/ZcUQ/jrrsBP2zdj/bSFwovvYs1iBYzfwdgn7\nrcrnaTh2PSAL+/V2WRnbiofGYSWk77AveNFWHVK6kRS+WFgdlXQhUCL3KZXrfyEpLpYl/gystNcV\na4mRhpXARIKSsSdsouiOfbde8zoQST6xTPx7sCaKh2D1gIdQvJWKlK60n+jVSSq8x1h4AfgE64Ox\nz+NYRIoZhdW3biPyC4IiIpKkOmEXHptjJf6p6GKMiIjnYtmOvxvWQzTYeeht4DdYN3QAOnXq5Kxe\nvTqGIYiIVEursRZ5lRLLOv7lWM+6ethFvHOwXwAFVq9ejeM4CX8bP3685zEoTsWpOBVj8IbVqFRa\nLBP/t1gnk2CbXYCnY7g/ERGJQKyHbHjAvYmISIJQz90I+P1+r0OIiOKMLsUZXckQZzLEGA1ed6Bx\n3PoqERGJkM/ngyrkb43OKSIx0axZM3bvLmuSOSlP06ZN2bWrIiOBR0YlfhGJCZ/Ph77fVVPaMaxq\niV91/CIiKUaJX0QkxSjxi4gkkYMHq/4aSvwiIknk3nur/hpK/CIiMZSRkcGMGTMKrRs1ahRHH300\naWlpvPDCCxG/1ldfwb/+VfWYlPhFRGKopJY5Xbp04cknn6Rr167BFjrlys6GESPgsceqHpMSv4ik\nnPvvv5927drRqFEjjj76aGbMmEFmZiaDBw9m6NChNGrUiFNOOYXFixcXPGf9+vUMHDiQVq1a0aJF\nC2688UbABpvs3bs3LVq0oGXLlgwfPpysrCwALr/8ctatW0e/fv1o2LAhDz30EADXXXcdvXv3pm7d\nuhHHfNtt0K0bXHpp1d+/Er+IpJQVK1bwxBNPMH/+fPbs2cP06dPJyMgA4N133+XSSy9l9+7d/O53\nv2PAgAHk5eWRl5fHRRddRMeOHVm7di0bN25k6NDQTLJ33nknmzdvZtmyZaxfv57MzEwAJk+eTPv2\n7Xn//ffZu3cvt9xyS6Vinj4dpk2Dxx+v6rs3Svwi4gmfLzq3ikpLS2P//v0sXbqUgwcP0r59ew4/\n/HAAunXrxsCBA0lLS2Ps2LHk5OQwb948vvrqKzZv3syDDz5IvXr1qFOnDqeddhoAnTp14uyzz6ZW\nrVq0aNGCP/3pT8yaNStqx2nXLrj6apg0CZo0ic5rKvGLiCccJzq3ijriiCOYOHEimZmZtG7dmmHD\nhrF582YA2rVrV7Cdz+ejXbt2bNq0iQ0bNtChQwdq1CieMrdu3crQoUNp164djRs35vLLL2fnzp3F\ntqus66+HgQPh7LOj9pJK/CKSeoYNG8bs2bNZu3YtPp+P2267DZ/Px/r16wu2yc/PZ8OGDRx22GGk\np6ezbt068vLyir3WHXfcQVpaGkuWLCErK4vJkyeTn59f8HikF29Ls2gR/OMfVXqJYpT4RSSlrFy5\nkhkzZrB//37q1KlD3bp1SUtLA+Cbb75h6tSp5ObmMnHiROrWrcupp55K9+7dOfTQQ7n99tvJzs4m\nJyeHuXPnAvDLL79Qv359GjVqxMaNG3nwwQcL7a9169YUnWL24MGD5OTkkJ+fz4EDB8jJySl1XKOX\nXoJ69WJwIDzkiEj1lKjf78WLFzs9evRwGjZs6DRr1szp16+fs2nTJiczM9MZPHiwM2TIEKdhw4ZO\n165dnYULFxY8b926dc6AAQOc5s2bOy1atHBGjx7tOI7jLF261DnllFOcBg0aOCeffLIzYcIEJz09\nveB506ZNc9q3b+80adLEmTBhguM4jnPmmWc6Pp/PqVGjhuPz+Ryfz+fMmjWrWKylHUOgSqPfaXRO\nEYmJZBud8y9/+QurVq1i8uTJXodSQKNziojEUDKdpKoq1on/KGBh2C0LuCnG+xQRqTCfz1flC7HJ\nIp7vsgawEegBBC+dq6pHpJpKtqqeRFQdqnrOAVYTSvoiIuKBeCb+ocArcdyfiIiUIF5VPbWxap5j\nge1h61XVI1JNqaqn6mJV1VOzCjFVxPnANxRO+gAFgxkB+P1+/H5/nEISEUkOgUCAQCAQtdeLV4n/\nVeBDoOiMAyrxi1RTKvFXXaxK/PFI/PWBtUBHYG+Rx5T4RaopJf6qS+ZWPfuAFhRP+iIi1V7RqRdX\nrlxJ//79adWqFc2bN6dv376sXLkyrjGp566ISAwVLbVnZWUxYMAAVq5cydatW+nRowf9+/ePa0xK\n/CKScrycerF79+5ceeWVNGnShJo1azJmzBhWrFjB7t27PTkWXqjc8HoikvAS9fu9fPlyJz093dm8\nebPjOI6zdu1aZ/Xq1c748eOdWrVqOW+99ZaTm5vrPPTQQ07Hjh2d3NxcJzc31znxxBOdsWPHOtnZ\n2U5OTo4zZ84cx3EcZ9WqVc6nn37qHDhwwNm+fbtzxhlnOGPGjCnYX0ZGhvPf//631HimTp3qtG3b\ntsTHSjuGVHF0zng15xQRKcT3l+i0LXHGVywHhk+92Lx5c9q3b1/wWHDqRYCxY8cyYcIE5s2bR1pa\nWsHUi8FZuMKnXuzUqRNAwdSL99xzT0SxbNiwgRtuuIGHH364Qu+hqpT4RcQTFU3Y0RI+9eLSpUs5\n77zzChJvaVMv+ny+MqdeHD16NHPmzGHv3r3k5+fTrFmzcuPYvn07ffr04frrr2fIkCHRe4MRUB2/\niKQcr6de3L17N3369GHAgAGMGzcuNm+yDEr8IpJSvJ56cc+ePZx33nn06tWL++67L35vPIwSv4ik\nlP379zNu3DhatmzJoYceyo4dOwoScP/+/Xnttddo1qwZL7/8Mm+//TZpaWmkpaXx3nvvsWrVKtq3\nb096ejqvv/46AOPHj2fBggU0btyYfv36MWjQoEKl/HHjxvG3v/2Npk2bMmHCBN555x3mz5/PpEmT\naNiwIQ0bNqRRo0Zs2LAhbsfA61kH3AvUIlLdJFvPXU29KCKSYpLpJFVVSvwiImjqxXhSVY9INZVs\nVT2JSFU9IiISFUr8IiIpRolfRCTFaMgGEYmJpk2bpszF0lhp2rRpTF7X6/+KLu6KSLWwfDn06gVz\n50LnzrHdly7uioh47OBBuPxy+OtfY5/0o0GJX0Skiu69F1q0gGuv9TqSyMS6jr8J8B/gOGzigKuA\nL2K8TxGRuPnqK/jXv2DhQkiWSxqxTvyPAh8Ag9191Y/x/kRE4iY7G0aMgMceg7ZtvY4mcrE8PzUG\nFgKHl7GNLu6KSNK68UbYvRteeim++63qxd1Ylvg7AtuBScBJwDfAaCA7hvsUEYmL6dNh2jQIm489\nacQy8dcEugI3AF8DE4Hbgf8N3ygzM7Ng2e/34/f7YxiSiEjV7doFV18Nzz8PTZrEfn+BQIBAIBC1\n14tlVU8bYB5W8gfohSX+i8K2UVWPiCSdYcOgVSt49FFv9p/IVT1bgPVAZ2AlcA6wNIb7ExGJuSlT\nYNEiWLDA60gqL9aNj07CmnPWBlYDVwJZYY+rxC8iSWPDBujaFT78EE45xbs4qlri97rVqRK/iCSF\n/Hzo2xdOPx3uvtvbWDRkg4hIHDz5JOzZA+PGeR1J1anELyJSjngOwBYJlfhFRGIo2QZgi4QSv4hI\nGZJtALZIaCIWEZFSJOMAbJFQiV9EpATZ2VbFk2wDsEXC63OYLu6KSEK64Qb4+ef4D8AWiUTuuSsi\nkpSmT4d3303OAdgiocQvIhIm3gOweUFVPSIiYbwegC0SquoREYmS6jAAWyRU4hcRAVautN65Xg/A\nFgn13BURqaItW+D88+Hvf0/8pB8NSvwiktL27oULL4SRI+2ibipQVY+IpKwDB+Cii6BjR3jqqeTp\nnavx+EVEKsFxYMQIG2r5rbegZhI1dVGrHhGRShg3Dlavhk8/Ta6kHw0p9nZFRGz8nXfegc8/h0MO\n8Tqa+FPiF5GU8uabcP/9MGcONG/udTTeiEcd/xpgD5AHHAR6hD2mOn4RiZvPPoPBg20sni5dvI6m\n8pKhjt8B/MCuOOxLRKRES5bAJZdY79xkTvrREK92/F63HhKRFLZ+PVxwAUycCGef7XU03otH4neA\nT4H5wDVx2J+ISIHdu61X7ujRNgCbxKeq5zRgM9AS+ARYDswOPpiZmVmwod/vx+/3xyEkEUkFOTkw\nYAD06QM33+x1NJUXCAQIBAJRe714V8GMB34BJrj3dXFXRGIiLw+GDLE2+q+8AjWq0QA1iT5I2yFA\nQ3e5PtAH+C7G+xSRFOc4MGaMTarywgvVK+lHQ6yreloDU8P29TIwPcb7FJEU98ADMGsWzJ4Ndep4\nHU3i8bq1jap6RCSqJk+Gu++2XrmHHeZ1NLGhQdpERFwff2wDrwUCcMwxXkcTO8nQgUtEJOa++QaG\nD7cxeKpz0o8GXfIQkaT3449w8cXw9NNw2mleR5P4lPhFJKlt3w59+8Jdd8Fvf+t1NMkhksTfAEhz\nl48CLgZqxSwiEZEI7dtn0yZeein88Y9eR5M8Irk4sADoBTQFPge+Bg4Al0Vh/7q4KyKVkpsL/ftD\nq1bw3HPJM21iNMSjA5cPyAYGAk8ClwDHV3aHIiJV5Tjwhz/Y36efTq2kHw2RturpiZXwg3PQ69qA\niHhm/HhYvBhmzoRaqniusEgS/xhgHNYDdynQCZgZy6BERErz73/b2Dtz50KDBl5Hk5y8/oGkOn4R\nidi0aXYRd/Zs6NTJ62i8E8ueu++FLTtFtnWw1j1VpcQvIhGZO9cu5n74IXTr5nU03oplz93g0Mm/\nBdoAL7k7GgZsrewORUQqavlyGDjQxuFJ9aQfDZGcMb4BTolgXWWoxC8iZdq0CX7zG8jMhJEjvY4m\nMcSjOech2AXdoMPddSIiMZWVZXPljhqlpB9NkZwx+gJPAz+59zOAUcDHUdi/SvwiUqL9+y3pH300\nPP642uqHi9ewzHWBo7GLusuB/ZXdYRFK/CJSTH6+jbSZkwNvvAFpaeU/J5XEa1jmrkBHd/uT3HUv\nVnanIiJl+fOfYd06+OQTJf1YiCTxv4TV6y8C8sLWK/GLSNQ98gh88AHMmQP16nkdTfUUSeI/BTgW\nq+YREYmZV1+Fhx+2aRObNfM6muorklY9S4BDq7CPNGAhhTuEiYgUMmMG3HSTlfbbt/c6muotkhJ/\nS+B74CtCF3Ur0nN3tPv8hhWOTkRSwrffwtCh8PrrcMIJXkdT/UWS+DPdv8GqHh+RV/u0Ay4A7gXG\nVigyEUkJa9faZCqPPw5+v9fRpIZIEn8AG7KhO5bwvwK2Rfj6jwC3Ao0qE5yIVG87d9q0ibfcYrNo\nSXxEkvgvBR4EZrn3H8eS+RvlPO8i7ASxEPCXtlFmZmbBst/vx69TvkhK+PVXmyD9ootgzBivo0ls\ngUCAQCAQtdeLpAPAYuAcQqX8lsB/gRPLed59wOVALtYBrBHwFjAibBt14BJJQXl5MHgw1K8PL74I\nNTS1U4XEo+fud1iSD2boGsC3QEUuwZwJ3AL0K7JeiV8kxTgOXHcd/PCDteCpXdvriJJPPHrufoSN\ny/OKu6MhwIeV2JcyvIhw333wxRcwa5aSvlciPWMMAk5zl2dj0zBGg0r8Iilk0iS45x6bVOXQqvQO\nSnHxqOrpCGwBfnXv1wNaA2squ9MwSvwiKeKDD+Cqq6ykf9RRXkeT3OIxHv+bFB6jJ99dJyISka++\ngiuugKlTlfQTQSSJPw04EHZ/P1ArNuGISHWzfLnNlfvss9Czp9fRCESW+HcA/cPu93fXiYiUac0a\n6NMH/vEPa7MviSGSOqIjgJeBtu79DVj7/FVR2L/q+EWqqc2b4fTTrXPWDTd4HU31Eq8ZuAAauNvv\nrezOSqDEL1IN7dxp4+4MHQp33ul1NNVPPC7utgGexS7o7sXG5r+6sjsUkept7144/3y73XGH19FI\nSSJJ/M8D0wlV9fwA/ClWAYlI8vr1V+jXD7p2hfvv1wTpiSqSxN8CeI1Qk86D2Pg7IiIFDhyASy6B\nww6DJ55Q0k9kkST+X4DmYfdPBbJiE46IJKO8PBgxwgZbe/55TZCe6CIZq+dmbNrEw4G52C+AS2IZ\nlIgkD8eBa6+Fbdusd24t9fJJeGWV+Htgc+1+g42ueQeQA3wCrI99aCKS6BwHbr0VvvsOpk2DunW9\njkgiUVbi/zehOXZ7AncBTwC7gadjHJeIJIG//Q2mT7eSfkPNqp00yqrqqQHscpeHYCeCt9zbtzGO\nS0QS3KOP2iQqs2dDs2ZeRyMVUVaJP43QmDznADPDHovk2oCIVFOTJsGECfDpp9CmjdfRSEWVlcCn\nYPPs7gCysXH4AY4Efo5xXCKSoN5803rjzpwJHTp4HY1URnktbXtiPXenA/vcdZ2x4RsWRGH/GrJB\nJIl89JENr/zxx9Cli9fRpK54jtUTC0r8Ikli9mwYNAjeeQd+8xuvo0lt8RirR0RS3IIFlvRfeUVJ\nvzqIdeKvC3wJLAK+B/4e4/2JSJQtWwYXXghPPw3nnON1NBINsW6dkwOchV0crgnMAXq5f0Ukwf30\nk02k8sADMGCA19FItMSjqifb/VsbayK6q4xtRSRBbNpkJfxx4+Dyy72ORqIpHom/BlbVsxXrC/B9\nHPYpIlWwYwecey78/vdw3XVeRyPRFo+OWPlAF6Ax8DHgBwLBBzMzMws29Pv9+P3+OIQkIqXZswf6\n9rVx9ceN8zoaAQgEAgQCgai9Xrybc94N/Ao85N5Xc06RBJKdbTNnHXecxtRPZInenLMF0MRdrgec\nCyyM8T5FpBIOHIDBgyE9HR5/XEm/Oot1Vc+hwAvYCaYGMBn4b4z3KSIVlJcHw4dD7do2Dk8N9fCp\n1rw+p6uqR8RjjgPXXANr1sD772tM/WRQ1aoejbIpksIcB26+Gb7/3sbVV9JPDUr8Iinsnntgxgwb\nabNBA6+jkXhR4hdJURMn2tg7n30GTZt6HY3EkxK/SAp67jlL/J99Bq1bex2NxJsu7oqkmDfegNGj\nIRCAzp29jkYqQ+Pxi0jEPvwQRo6ETz6BE0/0OhqpLLXqEZGIfPaZzZ717rtK+qlO3TREUsD8+dYr\nd8oUOPVUr6MRrynxi1Rz338PF10EzzwDZ5/tdTSSCJT4RaqxH3+0iVQmTID+/b2ORhKFEr9INbVx\no02kctddcNllXkcjiUSJX6Qa2r7dJlK59lq7iYRTc06RaiYrC3r3tslU7r3X62gkFtSOX0QKZGfD\needBly7wz39qTP3qSolfRACbSKV/f2jVSmPqV3dK/CJCbi4MG2YTqrz+OtRU18xqTT13RVJcfj6M\nGmV1+++9p6Qv5dNHRCSJOQ6MHQsrVthEKnXqeB2RJAMlfpEklZNjbfRnzbKJVOrX9zoiSRaxvvyT\nDswElgJLgJtivD+Ras9x4M034ZhjYPVq+PhjaNLE66gkmcT64m4b97YIaAB8AwwAlrmP6+KuSAUs\nXAhjxsDPP9tEKmed5XVE4oWqXtyNdYl/C5b0AX7BEn7bGO9TpNrZsgWuvhrOP9+GX1iwQElfKi+e\nLX0zgJOBL+O4T5GklpMD998Pxx8PzZrZRdxRoyAtzevIJJnF6+JuA+BNYDRW8i+QmZlZsOz3+/H7\n/XEKSSRxOQ5MnQq33gonnABffAFHHOF1VOKVQCBAIBCI2uvFowNXLeB94ENgYpHHVMcvUsSiRVaP\nv2sXPPKIxtCX4hK9jt8HPAt8T/GkLyJhtm6Fa66xwdWGDbN6fCV9iYVYJ/7TgOHAWcBC99Y3xvsU\nSSr798MDD8Bxx0HjxrB8OfzhD+qBK7ET64/WHDTmv0iJHAfeeQduucUu3s6bB0ce6XVUkgpUphDx\nwLffwp/+ZBOm/PvfNlOWSLyoNC4SR9u2WXPMPn3gkkusQ5aSvsSbEr9IHOzfDw8+CMceCw0bWnv8\nP/5R9fjiDX3sRGLIcWDaNKvHP/ZYmDsXOnf2OipJdUr8IjGyeLHV42/dCk8+adU7IolAVT0iUbZt\nmzXHPPdcGDTIOmQp6UsiUeIXiZIDB2DCBGuPf8gh1h7/uutUjy+JRx9JkSpyHJvy8Oab4aijYPZs\nOPpor6MSKZ0Sv0gVfPed1eNv2gSPPw7nned1RCLlU1WPSCVs327NMc85B377W7uQq6QvyUKJX6QC\nDhyAhx+2ppl16sCyZXD99arHl+Sij6tIBBwH3n/f6vGPPFL1+JLclPhFyrFkCYwdC+vXwz//acMm\niyQzVfWIlGLHDmuO2bs39Otn9fhK+lIdKPGLFLFlC9x7LxxzjNXdL18ON94ItWp5HZlIdCjxiwB5\nefDhhzBwoCX8n36CWbOsaqdZM6+jE4mueMy5WxbNuSue2rABnnsOnn0WWra0IZOHDoVGjbyOTKR0\nVZ1zVxd3JeXk5sIHH8Azz8Dnn8OQITB1KnTt6nVkIvER68T/HHAhsA04Icb7EinTmjVWsp80CdLT\nrXT/6qtQv77XkYnEV6zr+CehydXFQwcPwttvW2ucbt1gzx6ry583D668UklfUlOsS/yzgYwY70Ok\nmFWr4D/f45EpAAALlUlEQVT/gRdesA5Xo0ZZdU69el5HJuI91fFLtbF/vyX3Z56xwdNGjICZM9XD\nVqQoJX5JesuXW7KfPBlOPNFK9wMG2Fg6IlKc54k/MzOzYNnv9+P3+z2LRZLHr7/Cm29awl+50urr\n586FI47wOjKR6AsEAgQCgai9Xjza8WcA71Fyqx6145cK+e47S/avvALdu8M119hwCupVK6mkqu34\nY92qZwowF+gMrAeujPH+pBrat8+aYPbsaa1zmjSB+fNDPW2V9EUqRj13JWEtWGCl+9deg169rHR/\n/vka+15EPXelWtmzB6ZMsYS/fTv8/vc2Kma7dl5HJlJ9qMQvnnMc+PprePppeOstGwb5mmvg3HMh\nLc3r6EQSj0r8krR+/hleeslK9/v2Wel+2TJo08bryESqN5X4Ja4cx5pdPv00TJtmF2uvuQbOOgtq\naJBwkYhUtcSvxC8x5Tg2tv3MmXYLBKBBA0v2I0bYUMgiUjFK/JJw1qyxBB9M9rm5VqL3++1vp07g\n8/qTJ5LElPjFc+vXh0rzM2dCdrYl+GCy79xZiV4kmpT4Je42bSpcdZOVFSrNn3WWDYqmRC8SO0r8\nEnNbtoRK84EA7NwJZ54ZSvbHHqsLsyLxpMQvUbdtmyX4YLLfuhXOOCOU6E84QYlexEtK/FJlO3bA\nrFmhRL9hA5x+eqiO/qST1JFKJJEo8UuF7d5tiT5YdbNmDZx2WqiOvksXjYcjksiU+KVcWVnw2Weh\nC7KrV9tIl8FEf8opSvQiyUSJX4rZswdmzw5V3axYAaeeGqqj795dQxmLJDMl/hSWlwebN1tVzdq1\nNknJzJmwdCn06BGqo+/RQ9MQilQnSvzVWF4ebNwYSuxr1oRua9dax6nmzSEjAzp0sPbzZ55ppfu6\ndT0NXURiSIk/ieXmWguakpL6mjWW9Fu2tKSekVH41qEDtG+vBC+SipT4E9jBg1YqLy2xb94MrVsX\nT+zB++npqqIRkeKU+D104IAl9pKS+po11vGpTZuSk3pGhs0qVbu2V9GLSLJK9MTfF5gIpAH/Ae4v\n8nhCJv4DB6wJZPC2c6cl9KIl9+3boW3b0hP7YYep9YyIRF8iJ/40YAVwDrAR+BoYBiwL2yaqid9x\nICcnlLD37CmcwMu7BbfPzYXGjUM3ny9Aly7+YlUybdsmVvv3QCCA3+/3OoxyKc7oUpzRkwwxQmJP\nvdgDWAWsce+/CvSncOIv4Dg2nG8kibmsm89XOGmXdEtPh+OPt+VGjYo/Xq9e4dElMzMDZGb6Y3io\noiNZPrSKM7oUZ/QkQ4zREMvEfxiwPuz+BuB/im50+OGhpF6rVuEEXFJS7tix9ITeqJFauUj8BX+1\nOjjF1oWvL2ldec+P5DX35+4nKyeLfCcfBwfHccpdznfycRwnouVove6qXav4aNVHBeui9Td8P1X9\nO2/9PB6e93CZ24W/78r+reprVFUsq3oGYXX817j3h2OJ/8awbZyeT/WmZhqk1Sw+4mN51UCRHIBo\nvMaaqWvoMKBDsS9wot3f9n/baHFBi4IPSHCb8A9NSevCX6cy6yqyL4CDMw5Ss3f5ZY5IqwEr8kWo\nyGvmz8ynxlk1Sk3YpfG5Xytf2M/GstaFr490Xfj6A/89QL1z6+HDh8/no4avRpnLNXw18Pl8ES2X\n91oVed0f3/6RIwYdUfC8aPwN30+xv5V4vS9f/pKew3tWbn/RiDnCv6O6jbKPQCXFMvGfCmRiyR9g\nHJBP4Qu8q4BOMYxBRKQ6Wg0c4XUQJamJBZcB1AYWAcd4GZCIiMTe+VjLnlVYiV9ERERERKqT54Ct\nwHdh65oBnwArgelAk7DHxgE/AMuBPnGKESAdmAksBZYAN7nrEynWusCXWJXZ98DfEzDGcGnAQuA9\n934ixrkGWIzF+ZW7LhHjbAK8iTWF/h5rJJFocR6FHcfgLQv7HiVanMH9LsXy0itAnQSMc7Qb3xJ3\nmQSMsVSnAydTOPE/APzZXb4N+Ie7fCyW1Gph1wVWAfGa2bUN0MVdboBVTx2TgLEe4v6tCXwB9ErA\nGIPGAi8D77r3EzHOn7AvU7hEjPMF4Cp3uSbQOEHjDKoBbMYKVIkWZwbwI5bsAV4DrkiwOI/HcmZd\nrAD1CdYIJpFiLFcGhRP/cqC1u9zGvQ92xrotbLuPsJZBXngH63GcqLEegvWEPo7EjLEd8ClwFqES\nfyLG+RPQvMi6RIuzMZaoikq0OMP1AWa7y4kWZzOsYNcUO4m+B5ybYHEOxoa4CboLS/hRi9GLs0Jr\nrPoH92/wjbTFOnkFbcA6gcVbBvYr5UsSL9Ya2Jl9K6GqqUSLEeAR4Fas+W5QIsbpYCeo+YT6myRa\nnB2B7cAkYAHwDFA/AeMMNxSY4i4nWpy7gAnAOmAT8DNWok6kOJdgtSXNsELeBVhhKmoxev1zwHFv\nZT0eTw2At7A6tb0lxOJ1rPlYlVQ74AysRF00Bq9jvAjYhtXzltZPJBHiBDgNO8mfD1yPfdmKxuF1\nnDWBrsCT7t99wO0lxOF1nEG1gX7AG6XE4XWcnYAxWAGvLfadH15CHF7GuRzr7zQd+BAr7OWVEEOl\nY/Qi8W/FfqYAHIolCbCB3NLDtmvnrouXWljSn4xV9UDixpoF/B9wCokX42+Ai7FqlClAb+yYJlqc\nYPXQYCXqqdj4UokW5wb39rV7/03sBLCFxIoz6HzgG+yYQuIdz27AXGAnkAu8DfQk8Y7nc26sZwK7\nsQu6iXYsy5RB8Yu7wfqo2yl+gaI29vN2NfGbL8AHvIhVUYRLpFhbELqKXw/4DDg7wWIs6kxCdfyJ\nFuchQEN3uT7wOVY3nWhxgv2vO7vLmW6MiRgn2GCMV4TdT7Q4T8KqUuq5+3sB+7WXaHG2cv+2x1pz\nBS/oJ1KMpZqC1aMdwAZsuxKrt/qUkpsk3YFdkV4OnBfHOHth1SiLCDVH65tgsZ6A1fEuwpog3uqu\nT6QYizqTUKueRIuzI3YsF2GJINjBMNHiBEtWXwPfYiXUxgkaZ31gB6ETKiRmnH8m1JzzBezXfqLF\n+Zkb4yJCVbqJFqOIiIiIiIiIiIiIiIiIiIiIiIiIiIhUb3di7e2/xfpZ9Ijx/gJYz+hInYqNmLoQ\nGzp5vLu+H4UH0BIRkQj0xLrZ13LvN8O6rMfSTGxohEitwDrZgfWa1DSjkrC8HqRNJBJtsB6hB937\nuwiNs3M3NonKd8C/w54TAB7GerwuA7pj4/GsBP7qbpOB9XR8CSulv4F15S+qD3bi+QZ4HeuhWlRL\nbLwXsAGylrnLI4HH3OXwnuHZ2KBw9bFxWb7EemZfXNIBEBFJNfWxZLkCeAIbmTSoadjyi9jooGAl\n9uAsZTdhQ4e0xsYzWe8+LwMbqqOnu92zwM1hz++KjZE0i9AJ4TbsZFPU3dgJ6W1gFKGJPq4glPiD\n+rmvWRO4D7jMXd/EfY+HIBJDKvFLMtiH1bePwkZ9DM6aBDb65xfY+EW9sQGrgoLjBC1xb1uxcaN+\nJDSa4Xpgnrv8EjZuU5APq7s/FivxLwRGYANnFfVXbDTF6cDvsMkwgq8R7khssK1LsdEh+2ADbi3E\nTjZ1KDzSokjU1fQ6AJEI5WOl5FlYtc4V2EiQwXHqN2IXVOuGPWd/2HP3h63PJ/TZDx+33EfJ45h/\ngiXz8vwIPIVNlrKd4tM6NsBOWr8nNKEGwEBsvlSRuFCJX5JBZ6ykHHQyNlF6XSxR78SS6iWVeO32\nhKap+x2hKQNxX/sLbMKWTu66+kViCbqwSLy52Djq4Z7DZtL6PGzdx1hVVNDJFYhdpFJU4pdk0ACr\nJ2+CJdQfsGqfLKx0vQS7sPplKc8va7aiFdh47M9hw+D+q8jjO7ALtFMI1dvfSfES+nDsYnK2G+Nl\nYft1sBPMIOykEZw4/WqsimgiVlVVA/vVoAu8IiIxkkHhSYJEUoKqeiTVxXteZxERERERERERERER\nEREREREREREREUlt/w9xcioONwNDbAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1184c2e50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.arange(100, 1000, 100)\n",
    "plt.plot(x, spcat1, x, spcat2)\n",
    "plt.legend(('spcat1', 'spcat2'))\n",
    "plt.title('Speed of Sparse Dummy Functions')\n",
    "plt.xlabel('Sample Size')\n",
    "plt.ylabel('Seconds')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### It is clear that the spcat2 function, which uses the non-zero-value index, is more efficient than the row by row method in spcat1. However, from additional testing, it was noticed that spcat2 does use more memory and for larger n (somewhere between 10k and 20k), my notebook ran out of memory. This is likely due to to the fact that the spcat2 function requires an ($n^2$,) array ones ones, an index array of shape ($n^2$,) and an index pointer array of shape ($n^2$,). Combined these must soak up a lot of memory, which makes very large n infeasible. Then the spcat1 function still may be useful in scenarios where there is a very large sample size and no way of accessing more memory. It will be very slow in compairson, but will make analysis feasible where it would otherwise be infeasible. As a reasult,  I will keep both functions as options."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}