{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Update 1**: A more complete and updated speed comparison can be found [here](http://rlhick.people.wm.edu/posts/comparing-the-speed-of-matlab-versus-pythonnumpy-partii.html).\n",
    "\n",
    "**Update 2**: Python and Matlab code edited on 4/5/2015.\n",
    "\n",
    "In this short note, we compare the speed of matlab and the scientific computing platform of python for a simple bootstrap of an ordinary least squares model.  Bottom line (with caveats): matlab is faster than python with this code.  One might be able to further optimize the python code below, but it isn't an obvious or easy process (see for example [advanced optimization techniques](http://scipy-lectures.github.io/advanced/optimizing/)).\n",
    "\n",
    "As an aside, this note demonstrates that even if one can't optimize python code significantly enough, it is possible to do computationally expensive calculations in matlab and return results to the ipython notebook. \n",
    "\n",
    "## Data Setup\n",
    "\n",
    "We will bootstrap the ordinary least squares model (ols) using 1000 replicates.  For generating the toy dataset, the true parameter values are\n",
    "$$\n",
    "\\beta=\\begin{bmatrix}\n",
    "10\\\\-.5\\\\.5\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "We perform the experiment for 3 different sample sizes (\\\\(n = \\begin{bmatrix}1,000 & 10,000 & 100,000 \\end{bmatrix}\\\\)). For each of the observations in the toy dataset, the independent variables are drawn from\n",
    "\n",
    "$$\n",
    "\\mu_x = \\begin{bmatrix} 10\\\\10 \\end{bmatrix}, \\sigma_x = \\begin{bmatrix} 4 & 0 \\\\ 0 & 4 \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "The dependent variable is constructed by drawing a vector of random normal variates from Normal(0,1).  Denoting this vector as \\\\(\\epsilon\\\\) calculate the dependent variable as\n",
    "$$\n",
    "\\mathbf{Y=Xb+\\epsilon}\n",
    "$$\n",
    "<!--TEASER_END-->\n",
    "\n",
    "## What is timed\n",
    "\n",
    "In the code below, the data generation steps described above are not timed.  For comparison purposes, we only include the steps required for the bootstrap process.  This includes:\n",
    "\n",
    "1. Sampling rows of data with replacement\n",
    "2. Calculating parameter estimates for the OLS model\n",
    "3. Storing results\n",
    "\n",
    "For these steps, we have done some basic code optimization in both python and matlab by preallocating arrays/matrices.\n",
    "\n",
    "It is also important to note that we use built-in functions for calculating parameter estimates (`ols` for matlab and `statsmodels.regression.linear_model.OLS` for python).  If these functions aren't optimized for speed, then these results will obviously be skewed towards that package having the more optimized linear regression model function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Starting MATLAB on ZMQ socket ipc:///tmp/pymatbridge\n",
      "Send 'exit' command to kill the server\n",
      ".MATLAB started and connected!\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "from timeit import timeit\n",
    "%load_ext pymatbridge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "reps, beta, n_array = 1000, [10,-.5,.5], [1000, 10000, 100000]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run the bootstrap in Matlab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%matlab -i reps,beta,n_array -o mat_time,store_beta\n",
    "\n",
    "mat_time = zeros(cols(n_array),2);\n",
    "    \n",
    "for i=1:3\n",
    "n=n_array(i);\n",
    "row_id =1:n;\n",
    "\n",
    "X = [normrnd(10,4,[n 2]) ones(n,1)];    \n",
    "    \n",
    "Y = X*beta' + normrnd(0,1,[n 1]);\n",
    "\n",
    "    store_beta = zeros(reps,3); \n",
    "    tic\n",
    "    for r = 1:reps\n",
    "       this_row = randsample(row_id,n,true);\n",
    "       store_beta(r,:) = (ols(Y(this_row),X(this_row,:)))'; \n",
    "    end\n",
    "    mat_time(i,:) = [n toc];\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Here are the Matlab results for n = 100,000\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>b0</th>\n",
       "      <th>b1</th>\n",
       "      <th>b2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td> 1000.000000</td>\n",
       "      <td> 1000.000000</td>\n",
       "      <td> 1000.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>    9.999913</td>\n",
       "      <td>   -0.501554</td>\n",
       "      <td>    0.518378</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>    0.000820</td>\n",
       "      <td>    0.000785</td>\n",
       "      <td>    0.011911</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>    9.996906</td>\n",
       "      <td>   -0.503905</td>\n",
       "      <td>    0.486200</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>    9.999368</td>\n",
       "      <td>   -0.502104</td>\n",
       "      <td>    0.510149</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>    9.999908</td>\n",
       "      <td>   -0.501557</td>\n",
       "      <td>    0.518545</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>   10.000502</td>\n",
       "      <td>   -0.501017</td>\n",
       "      <td>    0.526612</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>   10.002367</td>\n",
       "      <td>   -0.499192</td>\n",
       "      <td>    0.554077</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                b0           b1           b2\n",
       "count  1000.000000  1000.000000  1000.000000\n",
       "mean      9.999913    -0.501554     0.518378\n",
       "std       0.000820     0.000785     0.011911\n",
       "min       9.996906    -0.503905     0.486200\n",
       "25%       9.999368    -0.502104     0.510149\n",
       "50%       9.999908    -0.501557     0.518545\n",
       "75%      10.000502    -0.501017     0.526612\n",
       "max      10.002367    -0.499192     0.554077"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "matlab_results = pd.DataFrame(store_beta.copy(),columns=['b0','b1','b2'])\n",
    "print \"Here are the Matlab results for n = 100,000\"\n",
    "matlab_results.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run the bootstrap in python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def python_boot():\n",
    "    for r in np.arange(reps):\n",
    "        this_sample = np.random.choice(row_id, size=n, replace=True) # gives sampled row numbers\n",
    "        # Define data for this replicate:    \n",
    "        X_r = X[this_sample,:]\n",
    "        Y_r = Y[this_sample]\n",
    "        # Estimate model \n",
    "        store_beta[r,:] = sm.regression.linear_model.OLS(Y_r,X_r).fit(disp=0).params   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "python_time = np.zeros((len(n_array),2))\n",
    "count=0\n",
    "\n",
    "for n in n_array:\n",
    "    # Constructs toy data\n",
    "    row_id = range(0,n)\n",
    "    X1 = np.random.normal(10,4,(n,1))\n",
    "    X2 = np.random.normal(10,4,(n,1))\n",
    "    X=np.append(X1,X2,1)\n",
    "    X = np.append(X,np.tile(1,(n,1)),1)\n",
    "    error = np.random.randn(n,1)\n",
    "    Y = np.dot(X,beta).reshape((n,1)) + error\n",
    "    \n",
    "    store_beta = np.zeros((reps,X.shape[1]))\n",
    "    TimeIt = timeit(\"python_boot()\", setup=\"from __main__ import python_boot\", number=1)\n",
    "    python_time[count,:] = [n,TimeIt]\n",
    "    count=count+1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Here are the results for n = 100,000\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>b0</th>\n",
       "      <th>b1</th>\n",
       "      <th>b2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td> 1000.000000</td>\n",
       "      <td> 1000.000000</td>\n",
       "      <td> 1000.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>   10.001015</td>\n",
       "      <td>   -0.501155</td>\n",
       "      <td>    0.501976</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>    0.000806</td>\n",
       "      <td>    0.000773</td>\n",
       "      <td>    0.011775</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>    9.998609</td>\n",
       "      <td>   -0.503459</td>\n",
       "      <td>    0.461239</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>   10.000454</td>\n",
       "      <td>   -0.501663</td>\n",
       "      <td>    0.494047</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>   10.000983</td>\n",
       "      <td>   -0.501155</td>\n",
       "      <td>    0.502325</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>   10.001581</td>\n",
       "      <td>   -0.500628</td>\n",
       "      <td>    0.509737</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>   10.004093</td>\n",
       "      <td>   -0.498803</td>\n",
       "      <td>    0.546887</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                b0           b1           b2\n",
       "count  1000.000000  1000.000000  1000.000000\n",
       "mean     10.001015    -0.501155     0.501976\n",
       "std       0.000806     0.000773     0.011775\n",
       "min       9.998609    -0.503459     0.461239\n",
       "25%      10.000454    -0.501663     0.494047\n",
       "50%      10.000983    -0.501155     0.502325\n",
       "75%      10.001581    -0.500628     0.509737\n",
       "max      10.004093    -0.498803     0.546887"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results = pd.DataFrame(store_beta,columns=['b0','b1','b2'])\n",
    "print \"Here are the results for n = 100,000\"\n",
    "results.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Discussion\n",
    "\n",
    "Having run the bootstrap for \\\\(n = \\begin{bmatrix}1,000 & 10,000 & 100,000 \\end{bmatrix}\\\\), we see that \n",
    "\n",
    "1. Matlab outperforms python for any sample size.\n",
    "2. As the sample size increases, the gap between python and matlab increases\n",
    "\n",
    "This can be seen in the figure below.\n",
    "\n",
    "While running numerical analysis in python does take a performance hit, this short example shows the flexibility of the Ipython Notebook.  You can run code in matlab, transfer data in and out and retain all of the advantages of both numerical computing platforms.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes.AxesSubplot at 0x4f2fe50>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEPCAYAAABFpK+YAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xe4XFW5x/HvJLRQQgolhAAnJHRCUTooC+QiFxQpXhVE\nQRC5VyRIE0WEJaCoFCEIV9B7L71KV8QIcigCUiSI1BB6KCF0Akib+8e7J7PnnD2nzJk9a83s3+d5\nzpPZU9/Z52TNmne/+10gIiIiIiIiIiIiIiIiIiIiIiIiIiIRKwG3Am8Bc4FpyXUnAh+nfnYJFaCI\niDTfesDCwDjgXuDzwAnA7iGDEhGRxiwwgPvMSP79GPgAeDnZLuUSkYiI5GrYAO/3IfAiNuO/M7nu\nDOBN4Cpg6eaHJiIioa0I/A34AjARWAIYBZwFnB0uLBERGYzBpmsOxwb8I1PXTQHOBdbPuP/jwKTG\nQhMRKaRZwOSQAYwHNgIWSQK5F6vg+RSwJDAam/H/X53Hl1sQY8x86AAi4kMHEBkfOoDI+NABRCT3\ncbO/HP8I4HTgVeAm4GrgCuBbwNPAE8AY4LAcY2xnXaEDiEhX6AAi0xU6gMh0hQ6gSPqr6pkFbJhx\n/ddyiEVERDpA0VM9LnQAEXGhA4iMCx1AZFzoACLS9uNm278BEZEWC57jz8ur2JvTT+f8vErfXD+3\nF40LHUBkXOgAimQgZ+7mYTTFOPPXAd2BY2gVfbsTEaD+YKBBovPodyrSHB2b6hERkUA08OfLhQ4g\nIi50AJFxoQOIjAsdQJFo4G+dvYBrB3jfi4H9m/CaXdhaCiIi82ng7+2fwEdYU7q0W7HW1OsM4Dkm\nYwNud4MxVCpl+rMdtQvipH8+At7HeivFoDt0AJHpDh1AZLpDB1AkGvh7KwOPAXukrlsZWB74V5CI\n6rse+x0Ow5rlgVVqDQOGA88HiktEIqaBP9sF1Lal+Dpwfo/7fAq4G3gbeBY4NnXbncBiVGfea2W8\nxpnAbOybwa3A2j1unwDcnDz/dKxhXl+yymO7qE313AkcAtwDzAPOwXotXZG8zi3Asqn7T8Y+XF4D\nHgF27SeGvrghPLYTudABRMaFDqBINPBnexgbCCt9ir6KtZ5O+xD4LrYIzeewD4edk9s2xgbWrbGZ\n94MZr3F38vzjsIH/gtRtJeCbwA+xD4DnsEG6Gb4O7Aasin14dWOL6iwHvIK13gbryDo9+ZmAfRCe\nBqzWpDhEpEM1WMdfLjfnpyEPYDPbA7DF5bcA/prc9i71c/xHAycnlys5/rS9qH9wdwEsjVTJx18E\nnJK6fRT2zWFsH3GvjX3DSH+Yd/WI4w7sA6ViGrXfZHbBurCCLbhzf4/X+AXwozqvrzp+kebI/f9S\nqDN3+1EKfVZvGausuR9YnN6zfbAPgFOATwAjk+vOG8RreOzDYDz2eygDS1EdqJ9O3fd1bJnL8dis\nfCjSef951B63eAd7vwArYccNPu7x+N8O8fVFJDCleup7GVt45svYh0BPZwJ/xgbIYcCPqe7Pyie2\nq/Pcn8UOHn8WWBRYCEsdpT/wulKXR2MfLi8M7i0MyTPYMYFhPX6+1eDzueaE1TFc6AAi40IHUCQa\n+Pv2HewP8o2M2xbFZsj/AjYF9qE64L+c3L5cneddDBvo38LSOydhxwIqSsBXgM2xNM9JWApmbsPv\nJFtf36ymYwd6D0liWCKJJ+tAtYi0EQ38fXsaOwib5SDsJKvXsPx+Olf+JnACtnpZpaonXZt/NXAX\n8Ch2TOHB5DEVZSyl8jPswO4EYM8BxJuVG+wrX9jzfIH09jvAtsCnsbWTnwGOofH0YHeDj+tU3aED\niEx36ACi4FkodAjNoCZtxaHfqchQeHbAMxM1aWt7LnQAEXGhA4iMCx1AZFzoAILxrIbnOqwq8MBW\nvGSkVT0iIh3OMxIrj/4GltbdCc/7rXjpvMsmy3Veo9710r70OxUZCM8w7ETKn2Jnxh+B58XUPXL/\nv6QZv4hIq3g2xk6cLGMz/LtChKEcf75c6AAi4kIHEBkXOoDIuNAB5MqzHJ5zgCuxar/NQg360P/A\nX8L6yLyF1ZBPS64bCVyHtTC4ldrGXiIiAuBZGM/3sLLtF4HV8JyL73VGfEsNJI+0Hta0bDTwB6zV\nwIbYmaVTgSOAJYH9Mh6rHH9x6HcqkubZAWvr8ghwcFKqORBR/V9aBjuFfxNgBtVmZWOAOXUeozr+\n4tDvVASq5ZmeR/H8ewPPEE0d/4fY15R7scF/eapNxF7Fes0s2PTo2t9DwA6BY3DY18zQXOgAIuNC\nBxAZFzqAIfOMxHMi1s33RmAKnj8GjirTQKt6FsCWIrwM2KnOfep9NTkbeCq5/Dr2baHCJf92R7T9\nf1jjNbAePfdgjdpe7efxxwFrYsdBKrePwNol/yHHeHfE1gXI8hHwb1RX58rj9bWt7WJve24B9mQW\nJ/AOdzGFtZPyzIE+X+VyF5E6HBvgZgDrJteNpbNSPQ8AuyeXVwT+hnXi7M9x2EIlaXcA2zcvtH7t\ngK0GFkLMv1ORfHg2wXMXnjvw8xduGqrgqZ7xwEbYakyTgS8BfweuAQ7F8vvfB67KMcaQnsGWJVwH\n+CSW7krvs12xD8EdsIPc+2P969NlWqsl229hzdnSi59vh/X8fxOrjkov8nIn1giu3mOzZH3rctSm\nel7EWis/mjzvz7DU3U3YqmPXYJ1FKzYAbktinIGt2iVSbNXyzMuBXwGb4+s2dGw7k7DulO9gM8mj\nkusr5ZzvYQPWuDqPb9cZ/1eTyythA++vk+1/Qs3BmquwwRlszd1pPZ7roeQxa2Ifkn8FfpDcNhEb\naHfFFj85HNvHiyS335nEkvXYej5H7xm/o3bgfwHLPy4HrIGl327DBvgxyetOTe67DFbGuxuWttoW\n++AYk/Ha/f1OXT+3F40LHUBkXOgABsTKMw/HMxfPz/D9TsYaEXwFrlmQ+fXlTfJMYfgmvXHfUElU\nCVtJ6zxstavpVAfcc7C1Z/+Ipbi2pVrGWqL3jLsMnIp9AICt5LVNcvmL2AB8ebL9c2BfbJ3e65LH\nTqvz2KE6HvsAeAFbYP0f2LEMgEuB9ZPLu2EfOBcl29OxBeB3YHCrjYm0NxtLdgB+iZVnbjqI8szo\nxNmyobEBu1nK2OpYF2bcdj5wJDZD/wo2CL6UelxPbwKzU9vzqC5tOJ7a5RUBnsTSLhX1HjtUPZdf\nTG/3XH7x89Quv1im91q8A9HdwGM6WXfoACLTHTqAujyrYfX4E4ED8FwfOKIhi3Pgj9cL2Ax4V+zD\nIX0wd7DfUmZji5ykdVE72LdKvQ/aZ7BlJ3evc7tI5/IsiXXP3Av7lnxaq7pn5k29egbvbCz1syZ2\n4LdiDnYgN72CzkjquwL4DLAzdtD2MCy//5cmxjpUF2PppT2wbwGjsPLQFRt4Lte8sDqCCx1AZFzo\nAObzDMOzN5bSGQ2sheekThn0QQN/I67GDnpegR3crrgUGxzfhLrNl9JLGz6B5fl/jKVavoAdnH0v\n+6G9lkmsZ7DLL/a8Pf06L2ILwn8DWwJyJnDAAGIQaU+eTbACh32BHfHsg5+fzu0Y6sc/eCVs0N6H\nuGbnobXz71SKzjMeK23+DFaifkHARmrqxx+ZEjb7/QAN+iLtz7Mwdub7YcBvgdXxvBU2qPxp4B+c\nu4EJWEnnQDhirlZoLYf2RZpD+yPN0cr9YZWDn8PWuX0Y2ATP4y17/cA08A/OBqEDEJEh8qyO1eN3\nTHnmYCnHL82i36nEzcozjwL2xNa7/VWklTrK8YuIDIktbr4X8BPsrPi1OrFSJybt2KunmVzoAFpI\nvXoGx4UOIDIul2f1bIrnbjy349smVRu8V09eXqM4g39RvBY6AJH54irPjE6oHL+ISPNZeeZBWHnm\nb4CftGF5pnL8IiL9qpZn/hLraLtxkcozY1P0dI4LHUBEXOgAIuNCBxAZ1/AjPavjuR7PI3i2a15I\nwXRsjl9EZGjapzwzOsrxi0h7sfLMb2DrXF8HHNFh5ZnK8YuIzOfZFFuZ7gPg8/j5K8fJIKgtc75c\n6AAi4kIHEBkXOoDIuD5v9YzHcy7wO2zg30KDfuM04xeReFXLMw/FyjML0T0zb8rxi0h8epdnHlyg\n8kzl+EWkYKx75inYGtT74/lT2IA6j3L8+XKhA4iICx1AZFzoACLj8CyJ5yTgNmA6sI4G/Xz0N/CP\nB34PvAHMBqYm158IfJz62SWvAEWkw3mGsQ3bY4ubL4l1zzxZNfn56S+PNBHYCLgeWAFbbnArrMXp\nfcCF/TxeOX4Rqc+zGVal8z4wVZU6QAQ5/ieTH7BZ/wxg2WRbA7qINMa6Z/4cm0hWumcWvcVLywwm\nxz8JO9hye7J9BvAmcBWwdHPD6hgudAARcaEDiIwLHUAQnoXxfB94AHgOWAPP+Xi2DBxZoQx01j4W\nO9hyIHbgZSIwFxgO/AJYCEv/9FQGzgGeSrZfx741dCfbLvm3U7e/S7Heb1/blcuxxBN6u3I5lnjy\n3faUuJEfsBL7M5l7gEPwTEjdv3LfOOJt7XblcldyeU8iyKiMBO4AvlDn9ilYvj+LvrqJFJ1nDTx/\nwvMwns+GDqcNBO/OOQK4Fjtj7urU9Z8G7sdSRQdgs1oRkSrrnnk08DWq3TM/CBuUDMQ21JZtfgx8\nFTgPS9u8hvXOWKrO44s+43ehA4iICx1AZFzoAHLjGYZnHzwv4vktnmUG8CiXd1htJPiM/wayDwBf\nkEMsItLuasszd8Bzb+CIJIN69YjI0HmWx8ozHXA4cKHKMxsWvI5fRKQ+zyJY98xDgLOw7plvhw1K\n+qOBP1+OaulW0Tm0L9Ic7bw/rHvm57Humf/EFjefNYRndLTz/mgzGvhFZHA8a2DdM1cE/gvP9MAR\nySApxy8iA+MZhZVn7gH8BDhd5Zm5UI5fRALzDKe6uPm1WPfMOWGDkpgV/ai+Cx1ARFzoACLjQgcw\nIJ7N8dyL5694PpnjK7kcn7vdBK/jF5EiUnlmR1OOX0SqrDzzYKw889fA8SrPbDnl+EWkBaw8c0fg\nZKxl8kZDLM+UAiv6V0MXOoCIuNABRMaFDmA+6545Hc9DeLYNFIUL9LoxUo5fRHKi8szCUo5fpGis\nPHNv4FjgGuBIlWdGRTl+EWkiz+ZY98z3UPfMwhrMmrsyeC50ABFxoQOIjGvpq3mWx3MBcAlwErBF\nZIO+Cx1AkWjGL9LJquWZBwNnAvupPFOU4xfpRFae+QVsdv8Atri5yjPbg3L8IjJInjWx7pkTUPdM\nyaAcf75c6AAi4kIHEBnX9Gf0jMLzS+Bm4Dpg3TYa9F3oAIpEM36Rdte7PFPdM6VPyvGLtLPa8syp\nkVXqSGOU4xeRDNY98xfAlsD3gIvUPVMGSjn+fLnQAUTEhQ4gMq6hR3kWwXMEcD/wFLa4eSe0THah\nAyiS/gb+8cDvgTeA2cDU5PqR2MGjd4FbgWXzClBEsPJMz07Ag8CGWPfMH6omXxrRXx5pIrARcD2w\nAvAXYCvgy0AX9kFwBLAksF/G45XjFxmq2vLMA/H8OXBEkq/oxs3pwNbADGCd5LoxULeCoN2/foqE\nY+WZp+B5Gc+BeBYMHZK0RFTj5iTgUWAE8DI2y694HTL/KKN6AwG40AFExIUOIDKu7i2e4Xj2xfMi\nnrPwLN26sIJxoQOISDT9+McClwL7YHn9LPW+mpyNHYQC+4CYAXQn2y75t1O314ssHm3Hvr0Ba2M1\n+e9wJT/ifmZiE6044tN2HtuVy11EZCRwB9b3o2IGsG5yeSxK9YgMjWcCngvwPItn96TXjhRT8Bn/\nCOBa4DfA1anrrwEOBQ4Evg9clUt0Ip3OumcegnXP/G/UPVMisA3wcY+f3amWc76HlXOOq/P4os/4\nXegAIuJCBxCVEg7PTnhm4bkSz8qhQwrMhQ4gIsFn/DdQv9Z/+ybHIlIMnjV5mBOxb9T74bkhdEhS\nLOrVI9Iqtri5B74KHAecocXNJYN69Yi0PeueuQ9wDHasbE38/EodkY6jHL9UuNABBOHZAs/f8dyG\n5xOpW1yokCLlQgcQkeA5fhFphGcC1j3zU1j3zIs7oJGadAjl+EWaqXd55vF45oUNStqMcvwibcFO\nuNoJW9x8BrAhnifCBiUSRtG/2rrQAUTEhQ4gN5418fwZz4N4thngo1yeIbUhFzqAiCjHLxItz2is\nPHN3bL3b/1Z5prQD5fhFBqtannks1q7kSJVnShMpxy8SFc8W2OLm84Dt8NwXOCKRQdOau/lyoQOI\niAsdwJBY98wLgYuwMs1PD3HQd02Jq3O40AEUiWb8In2x8sxDgYOAM4B9VZ4p7U45fpEs1fLMk4H7\ngEPwPBk2KCkI5fhFWs6zFnAq1m78m3huDByRSFMpx58vFzqAiLjQAfTLMxrPqdjSeFcD6+U46Luc\nnrddudABFIlm/CJWnvlNrHvmlah7pnQ45fil2Dyfwsoz3wamqjxTIqAcv0guPCtgZZmbY90zL1H3\nTCkK5fjz5UIHEBEXOgDAyjM9R2KN1GYCa+CDtEx2LX692LnQARSJZvxSDLXdM+8DNlB5phSVcvzS\n+WrLMw9UeaZETjl+kYbVds88Buue+WHQmEQioBx/vlzoACLiWvZKnuF49gMeARbGyjNPi2zQd6ED\niIwLHUCRDGTgvwiYAzyQuu5E4OPUzy7ND02kAVaeeQ+wB9Y98z9Vky9SayB5pM2Ad4FzgSnJdSdg\nB8gu7OexyvFLa6g8UzpHFDn+24HJGddrQJfwPCOwxc0PAk7Heuuoe6ZIH4aS4z8DeBNbgWjp5oTT\ncVzoACLimvpsnhKeXYCHgPWx8syj2mjQd6EDiIwLHUCRNFrVcwZWJTEc+3p9ArBXnfueDTyVXH4d\nO3GmO9l2yb+dur1eZPF0xrZnLnAqM1mZWZzGnZwcVXza1vbgtiuXu2iRgaZrJmEz+ykZt03B8v/r\nZ9ymHL80j5Vn/hj4CtXFzWOq1BFphihy/GQE8WngfixVdAA2ixfJR233zCuw8sy5YYMS6Wy3UVu6\neSBwHpa2eQ34HbBUnccWvarChQ4gIq6hR3k+hec+PLfg56fOOoELHUBkXOgAIpL7uDmQGf8WeQch\n0ouVZ56AlRMfBlyq8kyR5lCvHomLlWceCnwXK8/8eRtV6og0QzQ5fpF8WffMnbHumfcCn8TPrwYT\nkTZS9K/mLnQAEXF1b/GsjedGPP/Es3XrQgrKhQ4gMi50ABGJIscvkg/PGKrlmT8Gfq3yTJH8Kccv\nrWflmftig/0VwI9Unikyn3L80mE8n8YWN38D+Cxe54CItJr68efLhQ4gGl38B56LgfOB4wFX8EHf\nhQ4gMi50AEWiGb/ky8ozD+MJDgFOAfbG807gqEQKTTl+yYeVZ+6CLdpzL3CoyjNFBkQ5fmlDninY\n4uZLA/vg+UvgiEQkRTn+fLnQAbSUZwye04AbsWqd9VODvgsWV5xc6AAi40IHUCSa8cvQWXnmtwAP\nXI66Z4pETTl+GZra8sypeO4PHJFIu1OOXyLlWRHrnrkJ1j3zMnXPFGkPyvHny4UOoOk8I/AcBdwH\nPAKsgR9Qy2SXe2ztxYUOIDIudABFohm/DIyVZ+6KlWfejbpnirQt5filf7XlmVPx3BQ4IpFOphy/\nBGTdM48BvoQ1VDtT3TNF2p9y/PlyoQNoiGc4nv8CHsb+RtbAc/oQB33XlNg6hwsdQGRc6ACKRDN+\nqVUtz3wd2FblmSKdRzl+MSrPFImFcvySs0r3TDgQOA34hrpninQ2Dfz5ckB34Biy9S7P/ASep3N8\nRUes+yIMh/ZHmkP7o2UGMvBfBHwGeAmYklw3ErgY2Aq4B/hicru0g2p55lLYDF/lmSIFMpA80mbA\nu8C5VAf+Y4AuYCpwBLAksF/GY5Xjj4nKM0XaQTTj5mTggdT2DGCd5PIYYE6dx+ngYAwq5ZmeOXjO\nwDM2dEgiUlfu42ajdfzLw/x88KvAQsCCTYmos7jQAeDZEvg78GXg3/B8G88rASJxAV4zZi50AJFx\noQMokmYe3K331eRsmN/T5XXs20J3su2Sfzt1e71gr+9ZkQc5m8dZk8lMxcozt6T2IFrr4tG2trVd\nb7tyuYvITKJ3qmfd5PJYlOqJh3XPPBrPK3g8nkVDhyQig5L7uDnQGX/P2fw1wKFY7ff3gauaGZQ0\noPXlmSLSwW4DPk79TMXKOa8D3gNuBcbVeWzRZ/yuJa/imYLnJjz/wLNVS15z8FzoACLjQgcQGRc6\ngIhEMePfos712zczEGmAyjNFpAHq1dOOPAtQXdz8MuCoQJU6ItJ86tUjPXgcdtbta8A2eP4RNiAR\naTfqx58v17Rn8qyE51LgHOA4YKs2G/Rd6AAi40IHEBkXOoAi0Yw/dlaOeRh2UP00YC91zxSRoVCO\nP1ZWnvlFrEf+XcBhKs8UKQTl+AvJsw6Wxx+LzfC7wwYkIp1EOf58uUHd2zMWz+nADVi1zic6aNB3\noQOIjAsdQGRc6ACKRDP+GPQuz1xD5Zkikhfl+EOz8sxpwCvAgW1WqSMiQ1JeAFgJWAVY1X5K+6Mc\nf4fyrIT11dkIOAS4XIubi3SicgkYT83gzqrJ9kRs9cLHkp+ZrYhIM/58OXquI2rlmd8DDsBm+icU\npDzToTVV0xzaH2mOtt8f5bFkD+6rAPOoDu6VAf4xYBaU3u35RGjG3yGq5ZknAn9D3TNF2lB5caqD\neXpwXxUbT9OD+1XMH+RLbwQJtw7N+FvByjOnAaOxPH532IBEpL7yQsDK9J65r4otNfs4vWfujwEv\nQ6kZ6drcx00N/HmytW2PxWb6RwO/UfdMkRiUhwMrkD24TwCeIXtwfw5KH+cdHEr1tCErz9yPWRzH\nJC4AVsfzauiwAnO0fQ63qRzaH2mOpu+PcglbKyQrLbMyMJfawX16cvlJKH3Q3FjiooG/2WwhlFOB\nV7iPQziP/w0dkkhnK48me3BfFVssKj24X5BcfhxKRSiqyKRUT7NUyzM3xJalVHmmSNOUFwUmkz24\nL0J2WmYmlF4LEu7QKMcfPSvPPBz4DjbTPwFPz/IsEelXeUGsrj1rcF8KeILMwZ0Xm3RQNRbK8UfL\nyjP/A+ueeSewPp5netzLoTxuhUP7Is1RyP1RHoYdPO0xuF+/Lja4z6Y6oD8IXJlsPwulj4KE3IE0\n8Deitjzz63huDhyRSETKJWBpsmfuk7DV49Iz95vgrFGw3SVQ+leYmItFqZ7BqJZn7oo1VFN5phRY\neUmqg3rPM1Y/Ah6ld1pmJpTeDhJu+1COPwqV8kyrxb8EOFrlmVIM5RHYLL3n4L4KsDjzB/OeB1dL\n6i7bOA38wVXLM+diZ90+MIhHOwqZx83k0L5Ic0SzP8oLAF1kz9yXBZ6kdnCvDPDPN/GgqiOa/RFc\n9Ad352KnMAO8DYwc4vPFo7Y88xDgCpVnSvsqD8M6RGYN7isBL1Ad3B8Frk22n4aS0pkdZqifKi8A\ny/Vxe/vN+FWeKW2rXMKW68wa3CcDb5I9c58FpfdCRCyZop/xd47a8sw7yC7PFIlAeQl6d4isDPbD\nqB3YL6d6UPXNIOFKdIb6qTIHWAhL89jsuFZ7zPg962Lxj8Ly+M0qz3Qob1nh0L5Ic/S5P8oLYwdV\nswb3UdR2iEyXRs5t05OZHPr7qIh+xr8l9ge4GtZ7+m56//LOBp5KLr8OzEjdxyX/htlelh3ZjL2B\nTYGjOY7H+LBmhw/19dbLNX5tt/n2IluDWxb++CqwKpy7FSwxAXZeChgP18+BebNh1zuAe+Ggh2HG\ns3DT5UmHyMjej7Yb3K5c7qJFmvmpciLwLDZzrohzxm/lmf8JHIXKMyVX5RJ2HCxr5r4y9q05a+b+\nVKd3iJS6op7xL4P9QT+K/RHvBOzTjKByZeWZ04CXgc8MsjxTpI7yGOovu/cutYP7uVSX3Stsh0gJ\nZygD/yjgfKxaYA4204+3dYGnC/tWsgGtK890KG9Z4Wj7fVFejNoOkenqmYWoHdyvodoh8vWMJ3Og\nSUeKo+3/PtrHUAb+x4ApzQokN73LM7+m8kypb/6ye1n93ccCs6gO7rcAv00uz2nTg6pSQJ175q6V\nZ34JqzS6HfieyjPFlIdRf9m9FbBjVfWW3VOHSMmbWjY0xMozpwFL0tzyTGkb5RJ2HCprcJ8EvEL2\n4P4klN4PEbFIQgP/oHiWBo4BdqG6uHnIGZpDecsKRy77ojyK+svuvU/24P44lOY1P5ZBcehvI82h\n/VERdVVPPDwLA1OxXP75wBoqz+wk5RHUX3ZvUWoH9z9ix3JmQkl/AyIZ2nvGb3n8nbBqnYeAQ/E8\nmtvrSY7KC1LtENlzcF+G7GX3HqPzlt0TUaqnLs/6wC+xSouD8fw5l9eRJioPA5Yne3BfEXie7DVV\nn9ZBVSkQpXp68SwHHAfsgOXx/yfiVbAchctblkvY2qk9Bvfr1sdO+Hud2sH95uTfJwq27J6jcH8b\nfXJof7RM+wz8nhHAQcDBwP8Cq+F5I2xQRVYeSf1l98rULrt3GVx8PWx/CZTeChSwiCTiT/VU6/F/\nDtyL1ePPGnpo0r/yImQvu7cqsAS9l92rzORfUd5dpGEFz/F7NsLy+COwPH53c8KSqvIC2ApMWYP7\nOKyzar1l9z4OELBIpyvowO+ZAPwU2AY4EjgncD1+oxxR5C3LJeovu9cFvEj24P5UE5fdc0SxL6Lh\n0P5Ic2h/VBTs4K5nMeBQrCb/11geXznhASv3teze29QO7n+luuyeeheJFEgcM37PMGB34HhsQDoc\nz9P5htauyouTfVB1FeyD/DF6z9xnQkkHwkXaQwFSPZ7NgFOS+x6E5/acY2oD5YXp3SGyMriPprZD\nZLo08mUdVBVpex088HtWwip1Ngd+AFyIp9MOFjrq5i3Lw7GTlrIG9+Wp7RCZHtyfa9ODqg7lcNMc\n2h9pDu2Pig7M8XuWAL6PLX04DdgHT+iGWTkZBny0HNlpmZWxVcDSeffpVDtEatk9EclF62b8nuHA\nnsCxwI3AEXiey/n1W6Q8muzBfVXgPbJn7o9r2T0RydAhqR7Pllg9/rtYHv+unF83B+VFqb/s3iJk\nD+4zofSLsrmcAAAFh0lEQVRakHBFpF11wMDvuQL4JPA94LIWrHM7BOWFgIlkz9yXotohsmcjsZfq\nHFR1KG9Z4dC+SHNof6Q5tD8qOiLHfw+wRzzr3M5fdi9rcF8BmE11QP8ncEVy+Vl1iBSRThC+nDOf\nl60su5c1uE8CXqP+sntF6hApIvHpgFRPLq9RXgBr8bsiNkvv+e8k4EPqL7v3dvNjEhFpiqgH/q2w\n9shjsbLMIzPu08AbKJeAMWQP6pXL47BSyGeBZ1L/Vi4/CaVXBvl+8uBQ3rLCoX2R5tD+SHNof1RE\nm+MvYYP+VKxV8o3AH4A7+n9oeVFs8O45mKf//Re9B/N/pK57HkrvNxh7K62H/pgrtC9qaX/U0v5o\noUYH/vWwlZSuTbZ/A+xC5sBfPo3aQX1x4DlqB/U7gcuq2x2zWMeo0AFERPuilvZHLe2PFmp04F8e\napqoPQlsUOe+jwM3UR3kX27TlgMiIh2hWeWcfeSjSqc26TXaUVfoACLSFTqAyHSFDiAyXaEDKJJG\nB/7nsFWbKiYm1/V0P8R8wlZL7Bk6gIhoX9TS/qil/WGiXVp2GJbe2RFL+zwMbBY0IhERyd3W2OD/\nFrZMooiIiIiISPsZD/weeAPr4zM1uX4kcB3WcfRWYNnUYw4CXknuv0vq+inAA8A84H+wlBnAgsC5\nwDvYcZHVc3gfzTQMK8G9Ndku8r5YGrgGi/cJYG2Kuz++i72vt4ArsdLtIu2Li4A5WOwVrXr/XwFe\nSJ7v2015NwU3EfgysCT2n3oOsBZwDPZLGAX8Ajgzuf8k4HmsN9Cm2C9jkeS2W7BfyljsZJXdkuv3\nBm7Azlg+APhTju+nGb6Dvfdbku0i74vLgVOxv4+J2H/sIu6PpbEeWKtgA/5VwIEUa19sBqxP7cDf\nive/BPASsDGwWvK8E5r5xsRWydoamAGsk1w3BvtAAPsUPzl1/2uAHZL7zE1dvwt20hrA1diBcoDh\nyf0Wa3bgTTIeO1N7U6oz/qLui3FYfAv3uL6I+2Mx7BvPathAdBWwK8XbF5OpHfjzfv+LAztjXYQr\npgH79xXksL5ulF4mYWWsd1B7EturwELJz3jsZLWKJ5P7jqe25PWp5HqS2yrP9RH21W9806NvjpOB\nw6FmfeSi7ovJ2Pu7GktvXI3N7Iq4P+YBR2CtzN/A/j4up5j7Iq0V77/ec9WlgX/gxgKXAvtA3bUF\n+mqs1AnnM2yHteq4h/6bSHX6vgA7D2Yd4BSsW+wb2Idilk7fH+OAE7BFl0ZjM9Lv1Llvp++L/gR/\n/xr4B2YkdoD3GOC25LrZVM82HAu8jzWXm03tyW0rY5/gL1Cbd0uf9JZ+ruHYp/XzTYy/WTYHvoXN\n5m5Ptu/G3kdXcp+i7AuwmF8Crgfexg7urUEx98e6WCPFf2AfgBdj6cCi/T/pOXDn/f5n9/FcMgQj\ngJuxAytpxwDnYTm5E4CzkusnY3+Mq2IHe16k9qDN/tgyjjcDuyfX743lzcdiVUPTc3gfzbYx1Rx/\nkffF/cC/Y3nt87B9UcT9sTKWv56CHei+CvghxdsXPXP8rXj/S2D7fhOs0kcHd5tgG2yGm/7ZnWqZ\n1nvYADgu9ZiDsXze89gBrop1sBzoO1hb68pXvgWwP45Yy9SybEK1qqfI+2Jj7Mz1t7CDcOkSxqLt\nj29jM803sLToYhRrX9xG7Tgxlda9/92wbwuvUj/FJiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIi\nItJc47BWub/C+utcR++2yiIdTU3apIjWxBa4GIedDr9T0GhEWmyB0AGIBPAU8Lvk8h+odjwUKQTN\n+KWI5qUuf4AmQFIwGvhFRApGA78UUc/FMoqw6pOIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIjk\n7f8BHQpiIH/bk3QAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x4f2fa50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Convert to pandas dataframe for plotting:\n",
    "matlab_data=pd.DataFrame(mat_time,columns=['n','Matlab Time'])\n",
    "python_data=pd.DataFrame(python_time,columns=['n','Python Time'])\n",
    "plot_data = pd.concat([matlab_data,python_data['Python Time']], axis=1)\n",
    "plot_data = plot_data.set_index('n')\n",
    "plot_data.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Matlab Time</th>\n",
       "      <th>Python Time</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>n</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1000  </th>\n",
       "      <td> 0.257269</td>\n",
       "      <td>  1.064270</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10000 </th>\n",
       "      <td> 0.808506</td>\n",
       "      <td>  3.776620</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>100000</th>\n",
       "      <td> 7.839279</td>\n",
       "      <td> 31.610367</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        Matlab Time  Python Time\n",
       "n                               \n",
       "1000       0.257269     1.064270\n",
       "10000      0.808506     3.776620\n",
       "100000     7.839279    31.610367"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot_data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}