{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "In this note, I extend a [previous post](http://rlhick.people.wm.edu/posts/comparing-the-speed-of-matlab-versus-pythonnumpy.html) on comparing run-time speeds of various econometrics packages by\n", "\n", "1. Adding Stata to the original comparison of Matlab and Python\n", "2. Calculating runtime speeds by\n", " * Comparing full OLS estimation functions for each package\n", " - Stata: reg\n", " - Matlab: fitlm\n", " - Python: regression.linear_model.OLS from the statsmodels module.\n", " * Comparing the runtimes for calculations using linear algebra code for the OLS model: \\\$$(x'x)^{-1}x'y \\\$$ \n", "3. Since Stata and Matlab automatically parralelize some calculations, we parallelize the python code using the Parallel module.\n", "\n", "\n", "\n", "In addition to the above, I attempted to do some optimization using the Numba python module, that has been [shown to yield remarkable speedups](http://jakevdp.github.io/blog/2012/08/24/numba-vs-cython/), but saw no performance improvements for my code. \n", "\n", "\n", "The computational problem considered here is a fairly large bootstrap of a simple OLS model and is described in detail in the [previous post](http://rlhick.people.wm.edu/posts/comparing-the-speed-of-matlab-versus-pythonnumpy.html).\n", "\n", "**tl;dr**\n", "\n", "\n", "Time consuming econometric problems are best performed in Python or Matlab. Based on this comparison, Stata is dramatically slower (particularly when Parallel processing in either Python or Matlab).\n", "\n", "\n", "Matlab is the fastest platform when code avoids the use of certain Matlab functions (like fitlm). While slower, Python compares favorably to Matlab, particularly with the ability to use more than 12 processing cores when running jobs in parallel." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "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", "from numba import double\n", "from numba.decorators import jit, autojit\n", "from joblib import Parallel, delayed" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For boostrapping standard errors, we will consider 1,000 bootstrap replicate draws. We will explore several sample sizes (\\\$$n=\\begin{bmatrix}1000& 10,000& 100,000\\end{bmatrix}\\\$$) for the underlying dependent and independent variables. The true parameters are\n", "$$\n", "\\beta = \\begin{bmatrix} -.5 \\\\ .5 \\\\ 10\\end{bmatrix}\n", "$$" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [], "source": [ "reps, beta, n_array = 1000, [-.5,.5,10], [1000, 10000, 100000]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# A comparison of Bootstrapping with OLS Functions\n", "\n", "The first comparison we will perform uses the following functions:\n", "\n", "| | function|\n", "|:--:|:--:|:---|\n", "|Stata|reg|\n", "|Matlab|fitlm|\n", "|Python|regression.linear_model.OLS|\n", "\n", "\n", "**What is calculated**\n", "\n", "It is important to note several features of these OLS functions. The Stata reg command only calculate robust standard errors by request [need to verify this], whereas fitlm and regression.linear_model.OLS calculate several variants of robust standard errors, and all other factors equal should run slower due to these additional calculations.\n", "\n", "\n", "**Multiple Threads and Parallel Processing**\n", "\n", "In Stata and Matlab, the reg and fitlm are automatically multi-threaded without any user intervention. Consequently, all other factors equal python should run slower as by default regression.linear_model.OLS is not multithreaded. These comments are based on my observing cpu load using the unix top command. Python never extends much beyond 100%, whereas Stata and Matlab extend to the 200% to 300% range.\n", "\n", "### First consider the bootstrap in Matlab" ] }, { "cell_type": "code", "execution_count": 3, "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": [ "%load_ext pymatbridge" ] }, { "cell_type": "code", "execution_count": 4, "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:cols(n_array)\n", "n=n_array(i);\n", "row_id =1: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", " est = fitlm(X(this_row,:),Y(this_row),'linear','Intercept',false);\n", " store_beta(r,:) = (est.Coefficients.Estimate)'; \n", " end\n", " mat_time(i,:) = [n toc];\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here are the results for \\\$$N=100,000\\\$$." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
b1b2constant
count 1000.000000 1000.000000 1000.000000
mean -0.500783 0.499098 10.012360
std 0.000778 0.000803 0.011445
min -0.503204 0.496346 9.974214
1% -0.502515 0.497165 9.984511
2.5% -0.502226 0.497509 9.991038
50% -0.500783 0.499143 10.012116
95% -0.499493 0.500333 10.031592
97.5% -0.499266 0.500623 10.035487
99% -0.499011 0.500946 10.038212
max -0.498060 0.501551 10.046288
\n", "
" ], "text/plain": [ " b1 b2 constant\n", "count 1000.000000 1000.000000 1000.000000\n", "mean -0.500783 0.499098 10.012360\n", "std 0.000778 0.000803 0.011445\n", "min -0.503204 0.496346 9.974214\n", "1% -0.502515 0.497165 9.984511\n", "2.5% -0.502226 0.497509 9.991038\n", "50% -0.500783 0.499143 10.012116\n", "95% -0.499493 0.500333 10.031592\n", "97.5% -0.499266 0.500623 10.035487\n", "99% -0.499011 0.500946 10.038212\n", "max -0.498060 0.501551 10.046288" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "matlab_results = pd.DataFrame(store_beta.copy(),columns=['b1','b2','constant'])\n", "matlab_results.describe(percentiles=[.01,.025,.5,.95,.975,.99])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The same Bootstrap in Python\n", "\n", "Here is the python function implementing each replicate of the bootstrap. It samples with replacement from the data, calculates the OLS estimates, and saves them in a numpy matrix." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def python_boot(arg_reps,arg_row_id,arg_n,arg_X,arg_Y):\n", " store_beta = np.zeros((reps,X.shape[1]))\n", " for r in range(arg_reps):\n", " this_sample = np.random.choice(arg_row_id, size=arg_n, replace=True) # gives sampled row numbers\n", " # Define data for this replicate: \n", " X_r = arg_X[this_sample,:]\n", " Y_r = arg_Y[this_sample]\n", " # Estimate model \n", " store_beta[r,:] = sm.regression.linear_model.OLS(Y_r,X_r).fit(disp=0).params\n", " return store_beta" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Matlab employs a just in time compiler to translate code to machine binary executables. This substantially increases speed and is seemless from the user perspective since since it is performed automatically in the background when a script is run. The python [Numba Project](http://numba.pydata.org/) has developed a similar just in time compiler, with very minimal addtional coding required. One only needs to add @jit before functions you would like to compile, as shown below:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [], "source": [ "@jit(nogil=True)\n", "# rewriting python_boot to make function args explicit:\n", "def python_boot_numba(arg_reps,arg_row_id,arg_n,arg_X,arg_Y):\n", " store_beta = np.zeros((1000,3))\n", " for r in range(arg_reps):\n", " this_sample = np.random.choice(arg_row_id, size=arg_n, replace=True) # gives sampled row numbers\n", " # Define data for this replicate: \n", " X_r = arg_X[this_sample,:]\n", " Y_r = arg_Y[this_sample]\n", " # Estimate model \n", " store_beta[r,:] = sm.regression.linear_model.OLS(Y_r,X_r).fit(disp=0).params\n", " return store_beta" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number of observations= 1000\n", "1 loops, best of 3: 836 ms per loop\n", "1 loops, best of 3: 841 ms per loop\n", "Number of observations= 10000\n", "1 loops, best of 3: 3.58 s per loop\n", "1 loops, best of 3: 3.56 s per loop\n", "Number of observations= 100000\n", "1 loops, best of 3: 30.7 s per loop\n", "1 loops, best of 3: 32.5 s per loop\n" ] } ], "source": [ "for n in [1000,10000,100000]:\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", " print \"Number of observations= \",n\n", " %timeit python_boot(reps,row_id,n,X,Y)\n", " %timeit python_boot_numba(reps,row_id,n,X,Y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The numba speed (the second entry for each value of n) up actually is very small at best, exactly as predicted by the numba project's documentation since we don't have \"native\" python code (we call numpy functions which can't be compiled in optimal ways). Also, it looks like run times scale linearly.\n", "\n", "Next, is a printout of the results for \\\$$N=100,000 \\\$$." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
b1b2constant
count 1000.000000 1000.000000 1000.000000
mean -0.501051 0.499463 10.015594
std 0.000815 0.000770 0.011765
min -0.503515 0.496934 9.973254
1% -0.502901 0.497677 9.987255
2.5% -0.502635 0.498047 9.991328
50% -0.501082 0.499447 10.015637
95% -0.499684 0.500709 10.034549
97.5% -0.499387 0.501007 10.038958
99% -0.499134 0.501479 10.042510
max -0.498420 0.502184 10.048871
\n", "
" ], "text/plain": [ " b1 b2 constant\n", "count 1000.000000 1000.000000 1000.000000\n", "mean -0.501051 0.499463 10.015594\n", "std 0.000815 0.000770 0.011765\n", "min -0.503515 0.496934 9.973254\n", "1% -0.502901 0.497677 9.987255\n", "2.5% -0.502635 0.498047 9.991328\n", "50% -0.501082 0.499447 10.015637\n", "95% -0.499684 0.500709 10.034549\n", "97.5% -0.499387 0.501007 10.038958\n", "99% -0.499134 0.501479 10.042510\n", "max -0.498420 0.502184 10.048871" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "store_beta = python_boot(reps,row_id,100000,X,Y)\n", "results = pd.DataFrame(store_beta,columns=['b1','b2','constant'])\n", "results.describe(percentiles=[.01,.025,.5,.95,.975,.99])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Next we'll perform the same bootstrap in Stata. \n", "\n", "This is run in Stata 12.1 MP (2 cores). Admittedly, this is a fairly old version of stata, so perhaps newer ones are faster.\n", "\n", "\n", "clear\n", "set more off\n", "\n", "*before looping, generate variables so we can use replace in the loops\n", "set obs 10\n", "gen x1 = 0\n", "gen x2 = 0\n", "gen xb = 0\n", "gen error =0\n", "gen y = 0\n", "\n", "foreach num of numlist 1000 10000 100000 {\n", "set obs num'\n", "*generate x's\n", "replace x1 = rnormal(10,4)\n", "replace x2 = rnormal(10,4) \n", "\n", "replace xb = 10 + .5*x1 -.5*x2\n", "replace error = rnormal()\n", "\n", "replace y = xb + e\n", "\n", "timer clear\n", "timer on 1\n", "quietly bootstrap _b, reps(1000): quietly reg y x1 x2\n", "timer off 1\n", "quietly timer list\n", "di \"Time: \" r(t1)\n", "}\n", "\n", "estat bootstrap, percentile\n", "\n", "\n", "Produced these results\n", "\n", "\n", "obs was 10, now 1000\n", "Time: 3.445\n", "obs was 1000, now 10000\n", "Time: 10.346\n", "obs was 10000, now 100000\n", "Time: 91.113\n", "\n", ". \n", ". estat bootstrap, percentile\n", "\n", "Linear regression Number of obs = 100000\n", " Replications = 1000\n", "\n", "------------------------------------------------------------------------------\n", " | Observed Bootstrap\n", " y | Coef. Bias Std. Err. [95% Conf. Interval]\n", "-------------+----------------------------------------------------------------\n", " x1 | .50049835 4.44e-06 .00081187 .4990234 .5021256 (P)\n", " x2 | -.49990621 .0000552 .0008086 -.5013727 -.4982354 (P)\n", " _cons | 9.9916064 -.0004185 .01208307 9.967021 10.01425 (P)\n", "------------------------------------------------------------------------------\n", "(P) percentile confidence interval\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Discussion\n", "\n", "The following chart shows the performance of each statistical package using native OLS functions" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEPCAYAAABCyrPIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXm8XPP5x98ni0Q2JEKSBlcSQmSVImj1W3tRbfFDUJSW\nFrWrUmpE1VZborW0FKlS+74vEUu0kkgIsS8h+yISe/H8/njOdM7cO3PvzL0z8z1zzvN+ve7rzpk5\ny+d75t7ne77P8/0+DxiGYRiGYRiGYRiGYRiGYRiGYRiGYRiGYRiGYRiGYRhGm7kRWAS8FHmvB3A/\n8BnwFLB25LPjgKXAXGCPGmk0DMMwasBWwCjyO4RxwPXA6sD5wJXh+wOBecAGwJbAfGDVmik1DMMw\nqs4g8juEGcDw8HVPdAQBOjq4KLLf3cCuVVdnGIZhVIR2rTjmW8B74etlwCrhTz9gTmS/d8J9DcMw\njDqgNR1CIYIKnccwDMPwRIcS9pFG23OBBmAm0Av4EvgifH+9yH4DgIeKnPNNNOZgGIZhlMZbqAvf\nK41jCOOAiWj84ALgqsh+84AN0WD0AqBzkXM27mTSRsa3gBiR8S0gZmR8C4gZGd8CYoR3u/k08E3k\n52hy004/R6ed9onsfzwaV5gH7NnMeb03zDPX+hYQI671LSBmXOtbQMy41reAGFF1u9mSy+g7Rd7f\npcj7F5E/08gwDMMwmiXtIwTnW0CMcL4FxAznW0DMcL4FxIjE2s3ENswwDKNKVN1uVmraaaVYhjba\nfpLzs4zmcS18njacbwExw/kWkCZKmXZaS9YgHWsaHDDJs4ZaYaNBwzCapZiRMOORPOw7NYzKkDqX\nkWEYhuEJ6xD84HwLiBHOt4CY4XwLiBnOt4A0YR2Cfw4G7ilx35uAIytwzQZgZQXOYxhGgrAOoXRm\nAV8D6zZ6/yl0FffwJkc0ZRBqiCe1UkN25k5L7Ez+CvPoz9do/qnurdRQaSb5FhAzJvkWEDMm+RaQ\nJqxDKB0BXgcOiLw3AE3x/YUXRcV5EP1u2wHDwvc6hNvt0dQihmEYeViHUB43AD+NbB8I/KPRPt8F\nngc+Bt4Hzop89hzQldyT+iYFrnElmjl2JTr6GNro8/7Ak+H5H0brUDRHoWm8DeS7jJ4DTgCmAp8A\n16HJC28PrzOZ/FKpg9BO50PgVZrPW9USrg3HJhHnW0DMcL4FpAnrEMpjNmogNwu390fLiUb5CjgW\n6A3shnYaPwk/2wI1uNuiT+ovF7jG8+H5+6Adwg2RzwLg58Dv0I7hA9R4V4IDgbFottrvokP1vwB9\n0TrZJ4f7dUY7oodDDT8FJgCDK6TDMIyU0cp1CCKV+WkVL6FPwr8GxqOJ/54JP/uM4jGEM8gl/MvG\nEKIcTPGgcgfUHZX1998IXBL5fHV0pNGrGd1D0RFJtPNvaKRjCtrRZBlP/shnD+CJ8PWP0FoYUc4H\nTi9yfVuHYBiVoer/S3FbqdwCge9VzILO9JkJdKPp6AC0Y7gE2BRNFQ5aP6JUMmgn0Q/9fgRYk5wB\nfy+y73JgRbjv0jKuUYhoXOET8uMin6LtBS2CNAztZKL8rY3XNwzDM+YyKp/FwDRgH7RzaMyVwCOo\n4WwHnEnuPmd7eFfk3DuhQeudgC5oreqvyI8DNERer4F2OvPLa0KbmIPGHNo1+jmsledzlZGVGJxv\nATHD+RaQJqxDaB1HoX+oHxX4rAv6RP0FsCVwKLmOYHH4ed8i5+2KdgArUTfRhWisIUsA7AtsjbqL\nLkRdOUta3ZLCNDcSexgNMJ8Qauge6ikUIDcMo40I9JOmk1eqgnUIreM9NPhbiOPQxWMfovGD6Be5\nAi07+mdys4yiawvuAv4DvIbGLF4Oj8kiqGvmXDSg3B84qAS9hXyPzfkjG693iG5/CuwIbIPWxp6D\nllVtrftxUiuPSyqTfAuIGZN8C/CFwCoCJwEvojMWE4slt0sP9p0aRpkI7CjwqsD9ojP/wreTSdo7\nBOdbQA1p6Tt1tRBRRzjfAmKG8y2glgg0CNwu8JbADyXffWvZTg3DMJKOwKqiLuZpwHRgkwDuCVLy\nkJz2EUKasO/UMIogEAj8SOAdgVtFZyc2s3sysQ4hPdh3ahgFEBgs8KDAbIHtSzskmaS9Q3C+BdQQ\niyGUh/MtIGY43wIqjUB3gfMElgicILreqMRDq4vFEAzDMGpA6B7aD82J1hcYFsCFgaajTzVpHyGk\nCftOjdQjMFzgSYEXRBdytvI0ycQ6hPRg36mRWgTWEBgvsEjgV5KfeaAVp0smae8QXgF29azBoauh\nq43FEMrD+RYQM5xvAa1BoJ3AoQILBK4QTVBZgdNWF4shlM4scmUoF6F1CnqWcNwf0HoBUUothdkW\nfknxMppfoikBhhU72DCM1iGwOZoA8lBg1wB+GVQ+31iiqMcRwktoQAi0rvK/0cymLVGoQ5gC7FI5\naS2yK/5yocT5OzWMiiGwlsDfBOYJHCiVf+C2EUJMmYOWlxwOjAYWkH8v9wRmoIb4VDTZ3Tdo4ros\ng8PtlWhSu2jR+53Rmgsr0Kpp0eI7z6EJ9IodW4hC2Usd+S6jBWgK69fC856L1ot+Aq0SdzeaqTXL\nt4GnQ40z0CprhpE6BDqIFs7KJqPcOIDrg6Y1Q4wi1OsIYf/w9XqoQb4i3J4F/CCy752o0QatqTy+\n0bleCY8ZgrqdngFOCT9bHzXAe6JFaU5Gn+47h58/F2opdGwxdqPpCMGR3yHMBx5Dp8NtjBbfeRo1\n/D3D6x4d7rsWOgQeC6yKZj9dQGEXmsUQysP5FhAznG8BzSGwjcCLAo9J9VPAW8W0PDIVuiGZZvP9\nFyNAK59NRKuTPUzOEF+H1hZ+AC1nuSNweOS4xtcT4FK0YwCtvJZdqbgXaphvC7fPA36B1mG+Pzx2\nfJFj28o5aMcwH5iMpt2dGn52MzAqfD0W7YhuDLcfBp5ER0TlVIczjLpEdPR8AVpK9wTg1iTkHaq3\nDsFnCU1Bq5n9s8Bn/wBOQ5/o90WN48LIcY1ZAcyNbH9CrkRlP/LLZAK8g/4BZil2bFtpXEYzut24\njOYPyR8SC01rLZfCpFYck2Qm+RYQMyb5FhBFoBNwLPAb4HLgF4H+rySC+uoQ4st89Il5T7TTiAaR\ny31qmIsWn4nSQH4nUCuKdcBz0PKh+xX53DASh2hs71LgdWCLQAtEJQoLKleOa1EX0hA04JxlERpA\njuYr6dHMeW4HtgN+ggaLT0LjB49XUGtbuQl1Ux2AjhpWB3ZAZ1+Vi6ucrETgfAuIGc63AIH1ReOC\nlwEnBPDDJHYGYB1CJbkLDbbeDnweef9m1GiuIH+WUZTouoS30TjCmajL5kdoUPjzwoeWvKah3DKa\njT+PXmcBsBPwM7SU5xvoLAvDSAwCXUT/D6ei/7tDA7jXs6zYcizqxlgJ3IEavR5o4PMzdLrk2kWO\nrcdZRi0RoL7+bX0LiRn1/J0aKSRMQreHwLsCN0vrRr7VILb/S73RIvIboB3BncAxaLH161EXwvkU\nX7iVtA4hAA5BfYtGPvX6nRopRGAjgYcFXpb4PdzF9n+pK+raGIz6ue8ktxgru4iqJ+o/L0TSOoSp\nqBtlhxL3d9WTEjtsHUJ5ON8CYoarxUXCGgUXhDUKjhXoWIvrlkms7eO+wH/RqYfZIOpiYLXIPssp\nfGOT1iGUi/MtoIZYh1AezreAmOGqefLQPbS/wFyBawX6VPN6bSS2C9P6oIsyRqNz5q8Hjiqyb7Gp\ni9cC74avl6Ojiywu/D0podvZ9+Kix+f2pJjp8b09KWZ6fG9Pqtb5Re3OhPtgrbvh7KvgL1XQ35bt\n7OsGYs5OwH2R7bFo9s8ZwIjwvV6kx2VkFMe+UyNWCPQU+LPAQoHD2lijoJbENrndG8BmaPrk1YB9\n0FQKdwMnovGD36KxBaMpzreAGOF8C4gZzreAmOEqdSKB9qJpYGajnoshAVwVwNeVukaaOQKdg/4R\nOte+K7lpp5+j006L+ePSPkJwvgXUEIshlIfzLSBmuEqcRGCMwFSBZySXk6veSKx9THuHkCbsOzW8\nIbC2wDVh0PinUjymWQ8k9n/JOoTq0YAuFowL9p0aNUego8AxAosF/iTNp4upFxL7v1SPHcIYNCX0\nCmAZWisgm4Rue8qrT7w/rTfaO1O8NObXaLbUOGEuo/JwvgXEDFfuAQLfF5gl8IhobY+kENtpp2mj\nM5rD5Hy09GVHtIP4yoOWB8lNBtgE7Yg6YNWZjJQj0B/4E7AlcDxwexJqFKSBehshbIQa3FUKfNaN\npk/pXdCSks+j1c/eRyunZVnSaP8h6JPMk2iQfgE6J7rQ9aIMDc8RnS3WQP7o4zm0gMdUNG/7degs\nsNtDbZPJzzk1CO10PgReRVegt4W4fqdGQhDoJHCKwFKBcZJf6jVJJPZ/qd46hE6oUb8RrYa2RqPP\nt6Opy2hLYGu0xOQIdAHfT8LPBtLUZTQEXd/RLfx8OrkynMUotUOYiead+haacuRF1M3VHU1MeFG4\nb+fw8+PRWWOboRlXB7egozni+p0aCUBgF4E3BO4SGOBbT5VJ7P9SqzoEAanETys1r4PmQ38d+BJd\nmNc//KyUGMIZ5AxvKTGEg8ivq1CIUjqEKcDPI9vj0QpvWfYAnghf/4imVc/OB05vQUdzWAyhPJxv\nATHDFXpTYKDA3QKvS3498yQT24VpXgggqMRPKy//PpqeY0O0hGR74G/N7D8cLWqzHDXaZ6Crt4ux\nDnAPWq/5G+DvLexfDo1LY86PbDcujTmMfBfYiaE2w/COaI2Cs9D6BM8CwwKtZW5UgLrqEGLEfOBq\ncpldC/XcVwKPoEa2HVpoI3u/nyuw//lobqcNw/1+Ru2/nzmotnaNfg6r4jUnVfHc9cgk3wJixiT4\nXxK6vdBVxoOAEQGcG8AXPsUlDesQSqMBNeiDUT/7IOCX5Az7ItQ/H32i74I+fX+BxhMOJddxLA4/\nH9Bo/y/CY4aifvxq0NwI6WE0wHwCWtOiOxoH2aRKWgyjRUTja4+go+yDAhgbaJYEo8JYh1AaK9Cq\nSQ+iLqDJqAsp++T8EhpTeIfcLKPjgCPR2TpnkO+33xTNFjuV3Cyj09Hg9DLgr8A/qU5pTGn0eXT7\nUzRovg1aM3YOWvSomtOTXRXPXY843wLigkCPa+Bf6Oy7u4FRgY2gEkm9zTKqNM63gBpiQeXycL4F\n+EagncCBAvNugPtEa5UbCbaPae8Q0oR9p0bJCIwKE9A9L7CFbz0xI7H/S9YhpAf7To0WEeglcLlo\njYKfi7mzC5HY/6W0dwjOt4AaYi6j8nC+BdQS0RoFh4cdwWWiq+ijOB+6YorlMjIMI5kIbAVMIJzM\nEDRdFGmkhLSPENKEfadGHgJ9BK4TrVGwn9R3jYJaktj/JesQ0oN9pwbwvxoFxwksEThfdJ2LUTqJ\n/V8q1rBl5ObF208yfpbRPK6Fz9OG8y2gGghsJ/CywEOi2YNLxVVLUx2Sug4hLTjfAmKE8y0gZjjf\nAiqJwLoCNwu8I/BjKd895Kqhq05JrN1MbMMMwwCBzgK/E3UPZUTTwBttI7F2M7ENM4y0I7CbwJsC\ndwis71tPgkis3Uxsw0rE+RYQI5xvATHD+RbQWgQGCdwr8JposadK4Cp0niRQdbtpqwENw2gTAl0F\nzkaz/05GaxQ85FmWUUekfYRgGHWPaI2CvQXmCNwgmgLeqB6JtZuJbZhhpAGBTQQeF5gpmi7dqD6J\ntZuJbViJON8CYoTzLSBmON8CmkNgNYGLBRYLHCXVT3/jqnz+esJiCIZh+Ee0RsFBwKtoDe4hAVwW\nwFeepRkJIO0jBMOoGwRGC0wR+LfAZr71pJjE2s3ENswwkoLAmgJXCiwQOETMo+CbxNrNxDasRJxv\nATHC+RYQM5xvAaI1Co4QWCRwqcDqHuU4j9eOG1W3m1YPwTCM/yGwNXAZsALYPoAXPUsyUkDaRwiG\nESsE+gpcL/CBwL6tSEJnVJ/E2s3ENsww6gmBVQRODJPQnSM6g8iIDdIdZDeQS0iw3Uxsw0rE+RYQ\nI5xvATHD1epCAjsIzBZ4QGDDWl23TJxvAbVFOoJ8ByQD8jTIxyCPgZxCgu1mYhtWIs63gBjhfAuI\nGa7aFxBYT+BWgbcFdo+5e8j5FlBdJADZBOQYkHtBPgKZBnIeyA4gXaI7e5NZAr2Bu9EC2W8DQ4Ee\nwP3AZ8BTwNpFjo11wwwjiYQ1Ck4XWBr+thoFXpBvgRwEMhFkPsg7IFeB7A2yZnMH1kxiK7gNuBRY\nDc15vjYwDrgenaZ2PnBlkWNj3TDDSBJhErrdBd4SuE2gwbemdCGrgewOMh5kNshSkJtBDgMZUM6J\nqiaxjfQBlgCdGr0/Axgevu4JLCpyfGwbViOcbwExwvkWEDNcJU8msKHA/WGsYIdKnrtGON8CykdW\nAdkGZBzIsyArQR4BORlkNEj71p64ojIryHeA6cCDwErgLnRUsBgdMWRZDnQscHxsG1YjnG8BMcL5\nFhAzXCVOItAtnDW0JJxFtEolzusB51tAy0gAMgzkOJD7QFaATAU5B2Q7kEq55mJrNx2a1GpndJra\n9cA5FO4QCv0hxrZhhlHPhO6hfQXeF5go0M+3pmQi64AcDHIDyAKQN0GuANkLpFe1Llql8/6P1q5U\n/gBYiI4QAG4EDg/fbwBmAr2AL8OfQlwLvBu+Xo66myaF2y78bdu2bdslbgssBSbcB/3vg/Mvhwlx\n0lff2xt1g9kdgO3hgd3h0dVg+weBR2Dne+ChBVW4fvZ1A3XATOAHQHdgIhpQHhe+7glcAFxV5Ni0\njxCcbwExwvkWEDNcuQcIrB7mHFoU5iBqrY86jjg/l5VOIN8DOQvkuTAO8BDISSCjQHwk+ou13dwC\nmI3GEG5BXUfZaaefo9NO+xQ5NtYNqwHOt4AY4XwLiBmu1B3DGgWHCMwPs5I2N2WxXnG1uYy0AxkB\ncgLIA2Ec4D8gfwTZFqRzbXQ0S2LtZmIbZhi1QGCzsD7BFIFv+9ZTn8i6IIeA/BNkIcgbIH8B2QOk\np291BUis3Uxswwyjmgj0FvhrOCo42GoUlIOsERr7v4C8DrII5EaQQ0EafKsrgcTazcQ2rEScbwEx\nwvkWEDNcoTcFOggcGcYJLpb82XxJxrX+UOkE8n2Qs0H+HbqBHgjdQiM8xQHaQmxnGRmGUSMEvovW\nKFgKbBvALM+SYoq0QxfG7gBsD2wFvAw8CpwMTIHgC3/6jGKkfYRgGC0i0E/ghnBNwd4xT0LnCWkA\n+TnITSCLQV4D+TPIj0F8VnqrBom1m4ltmGG0lbBGwUnhKuOzBbr61hQfpCfIniCXh4vBFoaLw36m\nQeJEk1i7mdiGlYjzLSBGON8C4sTxmmbiVYH7BAb51uOf3juG6R/OAXk+jAPcD3J8mC4iTaOmxNrN\nxDasRJxvATHC+RYQBwQaBG5/AOYK7OZbjz+kHcimIL/RhHCPfoomiBuHJoyr15xMlSCxdjOxDTOM\nchBYVeCMsEbBaQJxWABVY2QAmgr6ZpAlaIroCSA/AknLbKpSSKzdTGzDDKMUwiR0PxZ4R+AWgaT7\nvyNIL5D/A7kS5C20SMxEtGhMf9/qYkxi7WZiG1YizreAGOF8C6g1AoMFHhR4RWC7Rh87H5qqi6wK\nsj1aFnIaWibyXpBjQYa2EAdwtVJZB9g6BMNICqKJIE8DDgXOBi4L4L9+VVUDaQ+MQtcCbI/mPXsR\nXQ9wDPAfCIplQTZSSNpHCEaKCN1D+wl8IHCdFE/6WKdIADIQ5HCQW9ESkS+DXAryQ5AevhUmhMTa\nzcQ2zDCiCAwXmCwwXXTlbEKQ3iD7gPwVLRI/D+R6kANBrChPdUis3Uxsw0rE+RYQI5xvAdVAYA2B\nCQILBQ4vo0aBq6au1iNdQHYEOR/khTAOcDfI0SBDqrgewFXpvPWIxRAMo54Is48eAvwBuAMYEmgO\nojpD2gOjycUBNgdeQOMAR6FxgATGPwwfpH2EYCQQgc0F/iPwrMCmvvWUhwQgG4D8CuQ2kGUgs0Au\nBtkVpLtvhUZy7WZiG2akD4G1BK4WmCdwYP3UKJC1QMaCXA3yHshckGtBDgDp61ud0YTE2s3ENqxE\nnG8BMcL5FtBawhoFRwssFriwQjUKXAXOUQTpCrIzyJ9AZoAsB7kL5CiQjWKaF8j5FhAjLIZgGHFE\n4HvABGAR8L0AXvEsqQDSgVwcYAe01OY0NA7wK+B5CL7yp88wlLSPEIw6RaC/wI0C7wnsFa8aBRKA\nDAY5EuQOkA9BXgS5CGQXkG6+FRptIrF2M7ENM5KJQCeBk8MaBWfFp0aBrA2yH8g1IHNA3g9f7weS\nsAVwqSexdjOxDSsR51tAjHC+BbSEwM4CrwvcLTCwypdzLajpFj7tXxQ+/X8YjgaODEcHMRqxVATn\nW0CMsBiCYfhCYABwMTAEOCaA+z2o6ABsRm49wKbAVDQO8AtgmsUBjHon7SMEI8YIdBE4M3QPnSLQ\nqYZXD8IZP78OZwAtD2cE/SmcIRQTV5XhgcTazcQ2zKhfwiR0ewi8K/AvgXVqdOW+4dz/a0E+CGMB\nV4Psq2sFDANIsN1MbMNKxPkWECOcbwEAAhsJPCwwS+D7Vb5a93D178XhauBl4ergX8EP909gHKAt\nON8CYkRi7WZiG1YizreAGOF8Xlygh8AF4eKyYwQ6VuEqHUG2BjkD5CmQj0EeBzkVZPMwb1AWV/nr\n1zXOt4AYkVi7mdiGGfVB6B46QLSo/TUCa1fw7EGYAfRoNCPocpDpaKbQHUG6VO5aRopIrN1MbMOM\n+CMwUuApgakCYyp01m+htQCuR2sDvIvWCtgHpHdlrmGknMTazcQ2rEScbwExwtXqQgI9Bf4c1ig4\nrIwaBYXO1gOtBnYpWh1sKcgtaNWwgW2IA7jWa0okzreAGGHrEAyjrYSG/1DgLOBWYOMAlpV5lo5o\nbeAd0PUAw4F/o+sBDgJegODryqk2jPSQ9hGCUSMExoSuoacFRpZxZAAyFORYkHvRCmHTQM4D2R5k\n1eqpNoyCJNZuJrZhRjwQWFvg72HQ+IDSktBJf5CDQCaCzAd5C+RKkP8DWbP6qg2jWRJrNxPbsBJx\nvgXECFfJkwl0DKePLg6nk/ZoZu/VQH4EMgFkNsgSkJtBDgMZUEldZeA8XTeuON8CYoTFEAyjVMIF\nZROAecA2AcxutMcq6KyibH2AocAUNA6wPzADgm9qKNkwDGyEYFQQgXXCVBPvCvwk5x6SAGQ4yPEg\n94OsAHke5ByQ7UA6+1VuGGURe7vZDngOeCrc7oFmhPwsfK/YYp/YN8yIP2GNglPDJHRnCnQBWQfk\nZyA3gCwEeRPkcpA9QXr61mwYbaDqdrOtxcCPAF4nJ/REYAnQFx2Kj2vj+ZOK8y0gRrjWHCSwCzDr\nczptvRe3/C5A1gyQF4DpwE7AE8AYCAZB8CsIboOgzKmmXnC+BcQM51uAURr9gMeALcmNEGag87MB\neqL1ZguR9hGC8y0gRrhydn6NDTaeS99nFtJ76d7cNBtkJchDICeCjARp60OOb5xvATHD+RYQI2Jt\nN29Ci3ZvQa5DWAysFtlnOYWThcW6YUackHYgI3qz8LeXc/ibS+j5zdmcMqc7H50L8n2LAxgpIrZ2\nc2fgivD1GJrvEFYpcHxsG2bEAVkP5FCQGwO+XnQA189bTK+Vr7HB5InsP8S3OsPwRGynnW4NHBb+\nZHke+ABoAGYCvYAvw59CXAu8G75ejrqbJoXbLvyd1O1jSVd7m9t2MKQb7DEKzuoL7ACPrAnLpu0N\nT1zDIes8yaf9zoTzLmPppfCGb73V3s6+jose39vZ13HRU8vt7OsG6oioy2gcMBGNH1wAXFXkmLSP\nEJxvAX6RziDbgvwR7pkdTgd9AOQEkBHvsu5qAheGi8t+LelaL+N8C4gZzreAGFEXdnMMMDl8nZ12\n+jnaSfQpckxdNMyoFNIOZBTISWEAeCXIFJCzQL4H0glAoJ3AgQLzBP4mYOUjDSNHYu1mYhtmZJEG\nkJ+D/AtkMcirIJeB/Bhk9SZ7w6YCzwj8R2BzD4INI+4k1m4mtmEl4nwLqDzSC2QvkCvCxWALw8Vh\nB+tiscIMh90FLhdYIHCotH1tTL3jfAuIGc63gBgR26CykXqkMzq5IJsXaEPUTfgo8GdgllapLHK0\n1ij4xSPwR+AGtEbBh1WXbRhG7Ej7CKEOkXYgo0FOBnkkjAM8CzIOZJswcVxpZ4KtBKYLTBYYUU3V\nhpEgEms3E9uwZCEDwlTQN4epoV8BGQ+yu6aOLvNs0EfgOoEPBMaWVqPAMIyQxNrNxDasRJxvAYWR\nNUH2BrkK5G20SMxEtGhM/1afVWsUHB8moTtPoHvkY9d23YnC+RYQM5xvATHCYghGNZFVge+gcYDt\ngUHoFOJHgUuBV5qLA5R0BdgOrVHwPrB1AK+1SbJhGIkj7SMET0h7kG+DnALyWBgHeBokA/IdtJB8\nZa4E6wrcIvCOwI/NPWQYbSaxdjOxDYsXEoAMAvklyK0gS0FeBrkUZDeQZspLtvKK0Fngd6F76AwB\nK0ZvGJUhsXYzsQ0rEVe9U0tvkH1A/gryLshckOtAfgrSr3rXBYHdBN4UuF1Kz7/iqiipHnG+BcQM\n51tAjLAYgtES0gX4Lrk4wADgSTQOcBHwalvjAC0qgA2AS4CBwBEBPFzN6xmGkSzSPkJoA9IeZHOQ\nU0GeAPkYZDLI70G2qmQcoEUl0FXg7NA9dJIUTnVuGEZlSKzdTGzDKo8EIBuCHAFyO8gykJdALgbZ\nFaR7y+eosCIdcuwtMEfgH6LV8wzDqC6JtZuJbViJuOY/lrVAxoJcDfIeyAcgfwfZH6RvTRQWUwZD\nBR4XmCnqqmorrgLnSBLOt4CY4XwLiBEWQ0gH0hU1rjugcYD10GIZj6J1JV6rdhygJUQr4WWA/YEz\ngSsD+MqnJsMwkkHKRwjSAWQMyGkgk8I4wJMgp4fvx6ajDmsUHCwwX+Aqgd6+NRlGSkms3Uxswwoj\nAchgkCNB7gRZDjIT5EKQH4B0862wEAKjBaYIPCewmW89hpFaMqyBuYzqGemDpm3ITgcV1AV0M2xx\nHfznDp/SPbC+AAAR9UlEQVTqmkNgTeBsYHfgVOC6AL6p0uUcuVqyht2PxjjSdj8y9ARGN/qpycjc\nOoSKId2AbcjVB+gPPIF2AucAb0TiAM6HwpYIaxQcjsYK/onWKFjuVZRhJJkMvWhq/HsBLwDTgDuB\n04E3qEHMzld+GfF47QohHdBSj9kRwKbA82gH8CgwDYK6CbqKJrmbAHwE/DqAlzxLMoxkkaE3TY3/\n6uSMf/bnDTIFR+RVt5vWIZSMBMBG5DqA7wHvkOsAnobgE3/6WodAX+B8dNRyInBzkLoYj2FUmAxr\n0dT49wCmk2/83ypi/AthHYJfpC+5DmB7dMj2CNoBPA7Bolae2OHZLxquKj4a+C3wV+DsAD72IMWR\nNh9x8zjsfkRxxP1+ZOhDzuhvGv7uRlPj/3YZxr8QVbebFkPIQ7qjT/7ZDqAfGgd4BDgLeMv3eoBK\nIBrjGI+OcLYK4HXPkgyjPsjQl6ZP/quSM/o3ACegxr/ubEXKRwjSkVwcYAdgJPBvcm6g6RB87U9f\nZRFd8HYRMAo4BrjX3EOGUYAMAfpAGH3qHw10Iv+pfxrwbo2Mv7mMKnzZABhCbgSwDfAW+XGAT2uv\nq7qENQlOQl1ElwIXBPC5X1WGERPU+H+Lpk/+7ckZ/az75z2PT/7WIbThEtkefmT4MwKdSfMFuTjA\nExAsrq6Ogjhq4BcNq5T9EE1NPR04IYD3qn3dMnHE3UdcWxx2P6I4Knk/1Pj3p6nxD2j65P9+zNw+\nFkMoDekIbIwa/ZGR398AM8Kfu4BTIHjLl8paIrAhOhpoAA4LtAM0jPSgxn9d8g3/pqhhzRr9q8Lf\nH8TM+HuhDkcIsjo5g581/huhT74zyXUAM4EFSQgCl4Po7IbTgJ+jC+ImBPClX1WGUWXU+K9H0yf/\n/9L0yX9enRr/NLuMJADWp6nx7wW8SL7xn5VE3385hO6hfdDsqE8AJwcw368qw6gCavzXp+mT/xc0\nNv4Z5vmSWQXS0iFIZ2Ao+cZ/OLCC3NN+1vi/DUG18urUCkcF/aICw9BVxqsBRwXwTKXOXQMc5jOP\n4rD7kSPAcQbv09T4f0q+8Z9OJvEPQEmOIchJ5J76B6C5OrLG/w79HSz1py/+iC57PxMYC5wBXBVA\nYqbJGikjQzu0LndumufbbI7m08oa/ovQJ/+F3nQmGI8jBLmE3JP/bAi+8KSl7hBoBxwM/BFNfnVa\nAEu8ijKMclDjP4j8J/9RaC6t/KmeGVqbESBppMVlZJRKWJfgMnQG1VGB/uMYRnxR478h+Qu8RgEf\n0tTt42MaeL1gHUJCcZTpJw4rlf0R2A3NPzSxijUKaonDfOZRHPV8PzK0R41/9Ml/JDqCjS7wmk6m\npFGto57vR2VJcgzBKAXR7+iXwO+BfwAbBTqsNgy/qPEfTFPjv4jcU/841Pgv8yXTKB0bIcQYge+i\n7qElwNEBvOxZkpFWMnRA1/tEZ/qMABbS1O3zoS+ZCcdcRmlENK/K+WiHcAJwqyWhM2qGGv+NyX/y\nHw7MI9/4v0DGKurVkNh2CP3QJd/fRXPon4emU+4B3AR8H5gK7AUFp4elvUNwFPCLhjUKjgV+A1wB\nnBNA3RXdKROH+YijOGp5PzJ0RBM+RrN6DgPm0tT4+3BVOuzvI0tsYwidgInA/sA6wOPAY+hK2SVo\nFa5TUf/h4W2XmXwEdkI71TeAMQG86VmSkTTU+G9C/pP/UGAOOcN/C2r8V/iSafijUr3Nw8C56KKR\nA9HUEj2BV4G1Cuyf9hHC/xBdgn8R+lR2bAD3epZkJIEMq9DU+G+C5vyKPvnPIMNKXzKNsojtCCHK\nQDSp1BTU951Nr7wMdYF0RBNMGRHCGgW/AX6NdghjrUaB0SoydEKf9KPGfwjwNrmpnjegxt9HmVSj\nTmhrh9ALuBk4FPisyD7FerRrgXfD18vRFcuTwm0X/k7ctkCwC/ztIdhlJ3gK2DTQ1B1j4qDPw3b2\ndVz0+N7Ovi78eYZO3MnBrMWGbEV3YDTvsAlfMY8NmAxM4x5e5A3eYgUPRo7vQK5mdi3b09bt7Ou4\n6KnldvZ1AzWiLcOPHsBDqKvorvC9GcBBaEqKXsBszGX0P0TnbI9/ADb8ARwaaOwl7TgsaBjFkb0f\nGTqjs3uiUz03QuNLUbfPTDIkNduvw/4+ssR2ltGqwIPAdcA1kffHoT7xY4BT0OybhxU4PlUdgkB3\n4HTgZ+hq48sCc6MZUTKsSr7xH42u+H2dfOP/Ipmio3Ej2cS2Q9geDSRHOQANiN4EbAs8D/wfsKDA\n8anoEMIaBWPRNQWPAr8NCt8PI01k6IIu6opO9dwAeI2mxt/iSkaW2HYIbSXxHYLoP/wEoCuahG5K\n5GOHDYOzOJJ8L9T4jyT/yX8g6k6NGv9ZofF3JPl+lI/D7keWuphlZEQQWAM4Cx0d/R74m9UoSAkZ\nupF78s/+DABeQY3+s+hDwiwyWLp3I3bYCKFChDUKDgH+ANwOnB6AFfhJKmr8R5Fv/BuAWeQyemaf\n/K2mtVEJzGVUDwhsjiah+y/qHnrBsySjkmToQVO3z7qo8Y+6fV4mY5MFjKphHUKcEZ1Sew6wM1qj\n4B8lJqFzmF80iyNO90KNf7SQy2igP7r6Pvrk/0qVjL8jTvfDPw67H1kshhBHwhoFR6BTSa8DNg6w\n3C91R4bVaGr8+6HGfxo6k+4cYDYZvvIl0zBqhY0QykT0iWUCOn306EBnixhxJ8Ma5Ix/9ndfdBFl\n1O3zqhl/I6aYyyguiLoNLgC2Ao4HbrcaBTElQ0+aPvmvha6kjxZwf5WMzQAz6gbrEHwjmur7OOBE\n4C/AuQFtThPgML9oFkdb7kWGXuSndhgNrEm+8Z8GvF4nxt9hfxtRHHY/slgMwSeiweLxqFto80Cz\nRxq+yLAm+U/9o9F1Hy+gRv8ONK7zOhm+8SXTMOoVGyEUQHQx0cVoCuFjArjfs6T0kaE3TY3/6uTP\n9JkGvGnG30gJ5jKqJQJd0OmjRwB/Ai4OsBWlVSfD2jQ1/t1pavzfMuNvpBhzGdWCMAndT9BCNVOA\nkQF8UMVLOtLqF83Qh6jhf4ut0L/DrNH/J3AC8DaZVAbtHWn92yiMw+5HzUh9hyCwMRon6AMcHNgf\nX+XI0I+mT/6dyD35/4Mp3MRAbkqp8TeMWJFal5FogZ/fowV9zgL+EmDzz1tFhgAKGv+O5Lt8pgHv\nmfE3jFZhMYQqXDgA9gfOQyu+nRLAQh9a6hI1/v3Jn+Y5GmhPU+M/x4y/YVQM6xAqfNGRaBK6TmgS\nun/XWkOIox5cU2r816Hpk7+Qv8BrGvB+K42/ox7uRe1w2P2I4rD7kcWCypVANC3xacBu4e9rAmy2\nSh5q/Nej6Qrfr8kZ/yvC33Ptyd8wkkeiRwiiKYpPRYvVXA5cGMCH1b5u7FHj30C+4d8U+JL8p/5p\nwDwz/oYRC2yE0BrCvEOnAPsCVwGDA1jiV5Un1PivT1Pj/zk5oz8emEaG+b5kGobhn0R1CKLZK09B\ng8ZXAxsFsNivqoI4quEXVeM/kKbG/xNyxv8S1PgvqPj1W4fDfMRRHHY/ojjsftSMRHQIAmsDJwMH\nA38HhiR+5lCGdhQ2/ivIGf8/AdPJJPxeGIZREeo6hhBWLDsJrWU8ETgvIIFuDzX+G5Af8B0FfET+\nNM/pZFjkS6ZhGFXFYgiFEE1vfCLwCzTVwfAA5vpVVSHU+G9I/pP/SGAZOcN/Lur2SWdcxDCMqlBX\nHYJATzTPzS+Bm9GcQ+/7VdUqHDCJDO2BweQ/+Y9EA+BZ4/8H9Ml/qR+pVcdhPuIoDrsfURx2P2pG\nXXQIojnvjwOOBG4DNg3gPb+qykSN/0bAaGawOyP5AzACjXVkp3qOQ43/Mn9CDcNIK7GOIQisBhwL\n/Bq4C/hDAO9UWVvbydCBrPHP/YxA4xtRn/8LZGxdhGEYJZHOGEKYeO5o4BjgPmBMAG/6VVUENf5D\nyJ/pMxyYR87w34Ea/+W+ZBqGYbRErDoEgW7oaOA4NPHc1gG87ldVhAwdyTf+o4FhaBwja/xvQY3/\nimbO5DC/aBaH3YsoDrsfURx2P2pGLDoEga5ofOAE4HHge4HWMfZHhlWATch/8h8KzCFn/P+FGv+V\nvmQahmFUCq8xhLBk5a/QtQSTgTMDeLnmatT4DyX/yX8TNF4R9fnPNONvGIYnkhtDEA0W/wYtWblj\nAC/W5MIZOtHU+A8B3iZn+CcCM8jwSU00GYZhxACfLiMH7BLAjKpdIWf8v40a/m+js3/eJDfV8zr0\nyb+Wxt9hftEsDrsXURx2P6I47H7UDG8dQgA/rugJ890+2Q5gCDnjPxW4BniRDJ9W9NqGYRgJINbr\nEIqis302If/Jfwjq859KrgOYacbfMIyEYCU0G031zHYAQ4F3yRn+aZjP3zCMZFO3HcL3UfdML7T4\nymmNPi/cMF3ktTH5T/7D0Kme0Sf/GWT4uEraa4HD/KJZHHYvojjsfkRx2P3IUpezjAK0MzgaNeCP\noauNp+Ttpbl9Nib/yX848AE5w38zyZznPxL7I89i9yIfux/52P2oIdXoEEYCy4F7wu2/AnvQuEPQ\nXP7zyD3530rLK3yTwuq+BcQIuxf52P3Ix+5HDalGh/At8jORvoOOAJrul+GjKlzfMAzDaAXtanCN\nwj6vdHcGDb4FxIgG3wJiRoNvATGjwbeANFGNEcIHwHqR7fXD96LMRAMkaeYg3wJihN2LfOx+5GP3\nQ3nLt4DW0A51E+2Ouo9mA1t5VWQYhmF4Y1u0U1gJ/NGzFsMwDMMwDMMwjOrTD7gXnT47F113AVrt\n7X7gM+ApYO3IMccBS8P994i8Pwx4CfgEuJpcsL8jcD3wKRp32agK7agk7YDn0HZDuu9Fb+BuVO/b\n6Or6tN6PY9F2rUSrB3YjXffiRmARqj1Lrdq/L1q6dylwREVaYxRkfWAftO7zUPQL3wQYh345qwPn\nA1eG+w9E111sAGyJfkmdw88mo19WL3QRztjw/UOAR4GeaCW5h6rYnkpwFNr2yeF2mu/FbcCl6N/H\n+ug/fBrvR2/gQ7Rt3YA70dK4aboXWwGjyO8QatH+7sBCYAtgcHje/pVsmFGch9HYyQx0xTXoF7Qo\nfH0ccFFk/7uBXcN9lkTe3wMtwwlwFxqcB2gf7te10sIrRD90VfqW5EYIab0XfVB9nRq9n8b70RUd\nIQ1GDdSdwJ6k714MIr9DqHb7uwE/AW6PHDMerUxZlFqsQ0gDA9GptlPIX5i3DFgl/OmH5mTK8k64\nbz/yp+W+G75P+Fn2XF+jQ8h+FVdfGS4CTga+ibyX1nsxCG3fXaib5C70STCN9+MT4FRgFupe/QYd\nPaXxXkSpRfuLnaso1iG0nV5ozqVDUX9gIZpLSJWE9Rg7o+lKptJy8q2k3wvQ9T3DgUuAvqghPLnI\nvkm/H32AC9BcZWugT7BHFdk36feiJby33zqEttEDDSyPA54O35tLbnVlL+BL4Ivw/eiCvQFojz+f\nfL9edCFf9Fzt0d59XgX1V4qtgcPQp79nw+3n0XY0hPuk5V6Aal4IPAh8jAYVNyad92MEWh73RbRj\nvAl1K6bt/6SxQa92++c2cy6jCqwKPIkGdKKMQ2sy90SfjK4K3x+E/pFuiAaZFpAfLDoSWDM8537h\n+4egfvle6Cymh6vQjkqzBbkYQprvxUzgB6jffCJ6L9J4Pwag/vFhaID9TuB3pO9eNI4h1KL93dF7\nPwadeWRB5SqyPfpEHP3Zj9x0ss9Rw9gncszxqL9wHhpYyzIc9bF+iqYOzw4dO6B/NHGdTleIMeRm\nGaX5XmyBrtJfiQb/olMt03Y/jkCfTD9C3atdSde9eJp8O3E0tWv/WHR0sYzirjrDMAzDMAzDMAzD\nMAzDMAzDMAzDMAzDMAzDMAzDMAzDMAyjZQahc8KXtbDfc+F+u1ZdkWF4xlJXGGnmGXSlaHOMQXPP\nG0bisQ7BMNTo3w/ciq4svs6vHMPwg3UIhqFsgxa0WQ9NEbCFXzmGUXusQzAM5Vk0p8wyNFFYg1c1\nhuEB6xAMQ/kk8vq/aMIww0gV1iEYhmEYgHUIhpGlcQGTNFToMow8bFhsGLrWYI/I9imNPm+pLKhh\nGIZRx5SzMO1rYJeqKzIMwzAMwzAMwzAMwzAMwzAMwzAMwzAMwzAMwzAMwzAMw4g//w+uJCk1ijCl\naAAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "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([[1000,.836],[10000,3.58],[100000,30.7]],columns=['n','Python Time'])\n", "stata_data=pd.DataFrame([[1000,3.445],[10000,10.346],[100000,91.113]],columns=['n','Stata Time'])\n", "plot_data = pd.concat([matlab_data,python_data['Python Time'],stata_data['Stata Time']], axis=1)\n", "#plot_data = plot_data.set_index('n')\n", "plot_data.plot(x=['n'],y=['Matlab Time','Python Time','Stata Time'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Having run the bootstrap for \\\$$n = \\begin{bmatrix}1,000 & 10,000 & 100,000 \\end{bmatrix}\\\$$, we see that \n", "\n", "1. Python outperforms Matlab and Stata for any sample size.\n", "2. As the sample size increases, the gap between python and matlab is constant, whereas for larger \\\$$n\\\$$, Stata's performance relative to either package deteriorates rapidly.\n", "\n", "\n", "Because we are relying on the \"canned\" OLS functions, the comparison above may be capturing the relative inefficiency of these functions rather than the underlying speed of the statistical platform. " ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "# A comparison of Bootstrapping (no functions)\n", "\n", "We will perform the exact same analysis as before with slight modifications to the functions for calculating the OLS estimates using linear algebra code for each package (\\\$$(x'x)^{-1}x'y\\\$$). For the sake of brevity, I won't show results, but instead just focus on runtimes.\n", "\n", "The Matlab code is:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%%matlab -i reps,beta,n_array -o mat_time_la,store_beta\n", "\n", "mat_time_la = zeros(cols(n_array),2);\n", " \n", "for i=1:cols(n_array)\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,:) = (inv(X(this_row,:)'*X(this_row,:))*X(this_row,:)'*Y(this_row))'; \n", " end\n", " mat_time_la(i,:) = [n toc];\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The python code is:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number of observations= 1000\n", "1 loops, best of 3: 309 ms per loop\n", "Number of observations= 10000\n", "1 loops, best of 3: 2.12 s per loop\n", "Number of observations= 100000\n", "1 loops, best of 3: 21.7 s per loop\n" ] } ], "source": [ "def python_boot_la(arg_reps,arg_row_id,arg_n,arg_X,arg_Y):\n", " store_beta = np.zeros((reps,X.shape[1]))\n", " for r in range(arg_reps):\n", " this_sample = np.random.choice(arg_row_id, size=arg_n, replace=True) # gives sampled row numbers\n", " # Define data for this replicate: \n", " X_r = arg_X[this_sample,:]\n", " Y_r = arg_Y[this_sample]\n", " # Estimate model \n", " store_beta[r,:] = (np.linalg.inv(np.dot(X_r.T,X_r)).dot(np.dot(X_r.T,Y_r))).T\n", " return store_beta\n", "\n", "for n in [1000,10000,100000]:\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", " print \"Number of observations= \",n\n", " %timeit python_boot_la(reps,row_id,n,X,Y)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEPCAYAAABCyrPIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXncXOP1wL83+45ESCJ4ZbEEIdKitHpHbbXW8mttRWlp\nVa1VpdSgaitqaS1tNaRKFRW1hprXVoqQEGtEBNmEkMQWxfn9ce7t3Jn3zvvOO++ducuc7+czn8yd\nu51787znPM85z3MOGIZhGIZhGIZhGIZhGIZhGIZhGIZhGIZhGIZhGIZhGEaXuQF4G3gu8Nsg4C7g\nY+BhYPXAvuOAd4F5wF4NktEwDMNoAFsBEyg1CGcC1wErA+cDV3m/jwbmA2OBrwALgL4Nk9QwDMOo\nO2MoNQjTgfHe98HoCAJ0dHBR4LjbgV3qLp1hGIYRCd1qOGcNYK73fQnQy/uMAN4IHDfHO9YwDMNI\nAbUYhDCciK5jGIZhxESPKo6Rsu15QAswAxgCfAqs8H5fO3DcKODeCtd8FY05GIZhGNUxG3Xhx0p5\nDOFMYDIaP7gAuDpw3HxgXTQYvRDoU+Ga5Uam2cjHLUCCyMctQMLIxy1AwsjHLUCCiF1vPgJ8Efgc\nTXHa6SfotNNhgeOPR+MK84G927lu7A8WM5PiFiBBTIpbgIQxKW4BEsakuAVIEHXXmx25jL5a4fed\nK/x+EaUzjQzDMAyjXZp9hODGLUCCcOMWIGG4cQuQMNy4BUgQmdWbmX0wwzCMOlF3vRnVtNOoWII+\ntH3sU+tnCenGjVuAhOHGLYBRf6STv2cNN24BEoQb8fXS3obcuAVIGG7cAiSItLftijS7QTDqh7Uh\nI6vUvW0nzWVkGIZhxIQZhHhw4xYgQbhxC5Aw3LgFSBhu3AI0E2YQ4ucQ4J9VHnsj8OMI7tkCLI/g\nOoZhZAgzCNUzE/gcWKvs94fRVdzj25zRljGoIm6tUQZ/Jk1H7ETpCvPg53M0/9TAGmWImta4BUgY\nrXELkDBa4xagmTCDUD0CvAIcGPhtFJrie0UsElXmHvT/thuwsfdbD2+7O5paxDAMowQzCJ3jeuC7\nge2DgL+UHfM14EngA+BN4KzAvseB/hR76huG3OMqNHPscnT0sVHZ/pHAg971p6J1KNojLDV5C6Uu\no8eBE4CngA+Ba9Hkhbd693mI0lKpY1Cj8x7wEu3nreoItwvnZhE3bgEShhu3AM2EGYTO8SKqIL/s\nbR+AlhMN8hlwLDAU2BU1Gnt6+7ZAFe62aE/9+ZB7POldfxhqEK4P7HOA7wO/QA3DW6jyjoKDgP3Q\nbLVfQ4fqvweGo3WyT/KO64MaoqmeDN8FLgPWi0gOwzCajBrXIYhE86mJ59Ce8E+AS9HEf496+z6m\ncgzhdIoJ//wYQpBDqBxU7oG6o3x//w3AbwP7V0ZHGkPakXsjdEQSNP4tZXI8hhoan0spHfnsBRS8\n73ugtTCCnA+c1o4MjcTWIRhZpe5tu5oCOQnCibsym6AzfWYAA2g7OgA1DL8FNkNThYPWj6iWPGok\nRqD/PwKsSlGBzw0c+z6wzDv23U7cI4xgXOFDSuMiH6HPC1oEaWPUyAT5YxfvbxhGzJjLqPMsBqYB\n30GNQzlXAfehirMbcAbF9+xbeLfCtXdEg9Y7Av3QWtWfURoHaAl8XwU1Ogs69whd4g005tCt7HN4\njddzoxErM7hxC5Aw3LgFaCbMINTGUWhDXRqyrx/ao14BfAU4jKIhWOztH17huv1RA7AcdRNdiMYa\nfBxgX2Br1F10IerKeafmJwmnvZHYVDTAfIInw0BPnrAAuWEYXURghLSdvFIXzCDUxlw0+BvGceji\nsffQ+EHwP3IZWnb0dxRnGQXXFkwBngBeRmMWz3vn+AjqmjkXDSiPBA6uQt4w32N7/sjy9Q7B7Y+A\nHYBt0NrYb6BlVWt1P7bWeF5WaY1bgITRGrcAcSHQS+BE4Fl0xmJmseR2Rr2wNmSkHoEdBF4SuEt0\n5p/3czZpdoPgxi1AgnAjvl7a25AbtwAJw41bgEYi0CJwq8Bsgd2k1H1r2U4NwzCyjkBfURfzNOBp\nYEMH/umkv4NTFc0+QjDqh7UhIzUIOAJ7CMwRuFl0dmI7h2cTMwhGvbA2ZKQCgfUE7hF4UWC76k7J\nJs1uENy4BUgQbsTXS3sbcuMWIGG4cQsQNQIDBc4TeEfgBNH1RlWeWl8shmAYhtEAPPfQ/mhOtOHA\nxg5c6Gg6+qam2UcIRv2wNmQkDoHxAg8KPCO6kLPGy2QTMwhGvbA2ZCQGgVUELhV4W+BHUpp5oIbL\nZZNmNwgvALvELIOLroaOGzfi66W9DblxC5Aw3LgFqAWBbgKHCSwUuFI0QWUEl60vFkOonpkUy1C+\njdYpGFzFeb9C6wUEqbYUZlf4IZXLaH6KpgTYuNLJhmHUhsDmaALIw4BdHPihE32+sUyRxhHCc2hA\nCLSu8n/QzKYdEWYQHgN2jk60DtmFjOdCCZDkNmRkGIHVBP4oMF/gIIm+w20jhITyBlpecjwwEVhI\n6bvcG5iOKuJT0GR3X6CJ63zW87aXo0ntgkXvd0JrLixDq6YFi+88jibQq3RuGGHZS11KXUYL0RTW\nL3vXPRetF11Aq8TdjmZq9fkS8Ign43S0ypphNB0CPUQLZ/nJKDdw4Dqnbc0QowJpHSEc4H1fG1XI\nV3rbM4FvBo69DVXaoDWVLy271gveOeNQt9OjwMnevnVQBbw3WpTmJLR338fb/7gnS9i5ldiVtiME\nl1KDsAD4FzodbgO0+M4jqOIf7N33aO/Y1dAh8H5AXzT76UKqc6GV49ZwTnskuQ1Vgxu3AAnDjVuA\n9hDYRuBZgX9J/VPAW8W0EvIRvZB8u/n+K+Gglc8mo9XJplJUxNeitYXvRstZ7gAcETiv/H4CXIIa\nBtDKa/5KxX1QxXyLt30e8AO0DvNd3rmXVji3q5yDGoYFwENo2t2nvH03ARO87/uhhugGb3sq8CA6\nIupMdTjDSCWio+cL0FK6JwA3ZyHvUNoMQpwlNAWtZvbXkH1/AU5Fe/T7ospxUeC8cpYB8wLbH1Is\nUTmC0jKZAHPQBuhT6dyuUl5GM7hdXkZzN0qHxELbWsvV0FrDOVmmNW4BEkZr3AIEEegNHAv8DLgC\n+IGjfyuZIF0GIbksQHvMe6NGIxhE7myvYR5afCZIC6VGoFFUMsBvoOVD96+w3zAyh2hs7xLgFWAL\nRwtEZQoLKkfHJNSFNA4NOPu8jQaQg/lKBrVznVuBbwB7osHiE9H4wQMRytpVbkTdVAeio4aVge3R\n2VedxY1OrEzgxi1AwnDjFkBgHdG44OXACQ7slkVjAGYQomQKGmy9Ffgk8PtNqNJcRuksoyDBdQmv\noXGEM1CXzR5oUPiT8FOrXtPQ2TKa5fuD91kI7Ah8Dy3lOQudZWEYmUGgn+jf4VPo3+5GDtwRs1iJ\n5VjUjbEc+Aeq9Aahgc+P0emSq1c4N42zjDrCQX3928YtSJOT5jZkJAAvCd1eAq8L3CS1jXzrQWLb\n9lC0iPxY1BDcBhyDFlu/DnUhnE/lhVtZMwgOcCjqWzTiJa1tyEgAAusLTBV4XpLXuUts2+6PujbW\nQ/3ct1FcjOUvohqM+s/DyJpBeAp1o2xf5fFu/URJHW7E10trG/Jx4xYgYbiNuIlXo+ACr0bBsQI9\nG3HfTpLotr0v8F906qEfRF0MrBQ45n3CX2zWDEJnceMWIEG4EV8v7W3IjVuAhOHW8+Kee+gAgXkC\nkwSG1fN+XSSxC9OGoYsyJqJz5q8DjqpwbKWpi5OA173v76OjCx/X+7c1o9v+b0mRJ87t1oTJE/d2\na8LkiXu7tV7XF9U7l90Jq90OZ18Nv6+D/F3Z9r+3kHB2BO4MbO+HZv+cDmzi/TaE5nEZGcnB2pDR\nLgKDBX4nsEjg8C7WKGgkiU1uNwv4Mpo+eSXgO2gqhduBn6Lxg5+jsQWjLW7cAiQIN24BEoYbtwAJ\nw43qQgLdRdPAvIh6LsY5cLUDn0d1j2bmSHQO+lJ0rn1/itNOP0GnnVbyxzX7CMGNW4AE4UZ8vbS3\nITduARKGG8VFBLYUeErgUSnm5EobaW/bFWl2g2DUD2tDxv8QWF3gGi9o/F2pHNNMA5lt22YQ6kcL\nuliwWbE2ZCDQU+AYgcUCv5H208Wkhcy27TQahC3RlNDLgCVorQA/Cd12dK4+8QHUrrR3onJpzM/R\nbKlpwo34ekluQ9Xgxi1AwnA7e4JATmCmwH2itT2yQmKnnTYbfdAcJuejpS97ogbisxhkuYfiZIAN\nUUPUA6vOZDQ5AiOB3wBfAY4Hbs1CjYJmIG0jhPVRhdsrZN8A2vbS+6ElJZ9Eq5+9iVZO83mn7Phx\naE/mQTRIvxCdEx12vyAbedcIzhZroXT08ThawOMpNG/7tegssFs92R6iNOfUGNTovAe8hK5ATxNJ\nbUNGnRDoLXCywLsCZ0ppqdcskdm2nTaD0BtV6jeg1dBWKdv/Ddq6jL4CbI2WmNwEXcC3p7dvNG1d\nRuPQ9R0DvP1PUyzDWYlqDcIMNO/UGmjKkWdRN9dANDHhRd6xfbz9x6Ozxr6MZlxdrwM5kkRS25BR\nBwR2FpglMEVgVNzy1JnMtu2aDIKARPGpUeY10XzorwCfogvzRnr7qokhnE5R8VYTQziY0roKYVRj\nEB4Dvh/YvhSt8OazF1Dwvu9B26pn5wOndSBHV3Ajvl7a/2jcuAVIGG7YjwKjBW4XeEVK65lnmcQu\nTIsFB5woPjXe/k00Pce6aAnJ7sAf2zl+PFrU5n1UaZ+Ort6uxJrAP9F6zV8Af+7g+M5QXhpzQWC7\nvDTmxpS6wH7qyWYYsSNao+AstD7Bv4GNHa1lbkRAqgxCglgA/IliZtcwy30VcB+qZLuhhTb89/14\nyPHno7md1vWO+x6N//95A5WtW9nn8Dres7WO104jrXELkDBa4X9J6PZBVxmPATZx4FwHVsQpXNYw\ng1AdLahCXw/1s48BfkhRsb+N+ueDPfp+aO97BRpPOIyi4Vjs7R9VdvwK75yNUD9+PWhvhDQVDTCf\ngNa0GIjGQTaskyyG0SGi8bX70FH2wQ7s52iWBCNizCBUxzK0atI9qAvoIdSF5Pecn0NjCnMozjI6\nDvgxOlvndEr99puh2WKfojjL6DQ0OL0E+APwV+pTGlPK9ge3P0KD5tugNWPfQIse1XN6slvHa6cR\nN24BkoLAoGvgb+jsu9uBCY6NoDJJ2mYZRY0btwAJwo34emlvQ27cAsSNQDeBgwTmXw93itYqN9Lf\ntivS7AbBqB/WhlKMwAQvAd2TAlvELU/CyGzbNoNg1AtrQylEYIjAFaI1Cr4v5s4OI7Ntu9kNghu3\nAAnCjfh6aW9DbtwCNBLRGgVHeIbgctFV9EHcOORKKJbLyDCMbCKwFXAZ3mQGp+2iSKNJaPYRglE/\nrA0lHIFhAteK1ijYX9Jdo6CRZLZtm0Ew6oW1oYQiWqPgOIF3BM4XXediVE9m23azGwQ3bgEShBvx\n9dLehty4BagHAt8QeF7gXtHswdXi1kumFGIxBMMw0ovogs7foJlzjwOmWI0Co5xmHyEY9cPaUAIQ\n6CPwC889lBdNA290jcy2bTMIRQ5Bs5xWw41oOoyu0kJj6y5/F7gyguv0QpObtbdytRnbUKIQ2FXg\nVYF/CKwTtzwZIrNtO40GYSaad2itst8fRtNEj29zRlvGoIrYDfx2CNUbhBuAI6s4Lkl1l7sBsygq\nhlU9OXwZ3JBz1vGO+X3IvhMp1pUII8ltqBrcuAWoFYExAncIvCxa7CkK3IiukwWsHkKCELQ4zoGB\n30ahWU6TloLXr7vcDa1vABov6obWcZhf4bx6sB2aDXZOJ845CK3q9h3alhH9q7e/ZyTSGV1GoL/A\n2Wj234fQGgX3xiyWUQNmEDrH9aj7w+cgSrOYQvu1lB9HS1M+gPbUw9JKXwXMQ0cSD6OpsIOMRLM/\nfoCmq+6otx82x7uFxtVd3oViRbYwWkPk/S6aIXYxsFvZ/nloTeqt2rlmmmmNW4Bq8WoUfBt147Wg\nNQrOd7SiYFS0RngtowPSNcuoUIhmyJTL1boQ5kVUQX4ZVfoHoOX7Tgwc8xlwLFoTeV00be/TaO3i\nLYDplM6//nLZPZ5Eay8sBX6BGqFNvH0OWg5zL9SFdRGqvLev8XmCHIQWIPkINUStaE2Gg4HrgJO8\n7T6oIbocNQTjgCmePC+HXHcT2q8sV85XUbfS3Wia8IOBW8qOeQGYgBpGIwZEOzOXoTVADnS002AY\nNZHGGMJzqAL8CVqX+KvAo96+j6kcQzidos+7szGEHqg7yjcgNwC/DexfGR1ptFdqM+66y88Cuwe2\nO4oh/BGY5H0fh/Y2h5Yd82fUaIaR5DZUDW7cArSHwEoCFwssFjhK6t+pdOt8/TRhMYSEIehMn32A\nQ9GeczmdraVcTh4tpfmp9+mBKlGfuYHv76PFe6IIEter7vKywLkd0Rd9tzd52y+go479y44bgI6g\njAbh1Sg4GHURDgDGOXC5oyNiIyOYQeg8i4FpaMDzxpD97dVS9i18a4Vr74gGrXdEq671Qv/ggi6u\nlsD3VYBBlCrvetPZussz0ZFRJVoD3/dEn+cOisZmQ1QRBVnXu24WaY1bgHIEJqKj4SOB3R34gaN/\nB42gtUH3MTCDUCtHoUPZsF5qZ2spB+mPGoDlqJvoQnRWkI8D7IvWOV7Z219Ag6xREmXd5amoe62a\nexwMXEKpoVkTHZH4wfXBqFF8NOR8I0IEVhXt4NyJlnX9iqMxLiOjmEGojblU/sNor5byMrSW8nSK\ns4yEosGYAjyBukmeA573zvER1Md+LlpkfCRte89hhPke2/NHStn+4HZn6y7fgcYCVi/7/U1K3U6n\nADng6rLj5qEB5oO87f8DbkbdWlnEjVsAr0bBkajL7hNgfQeucfT/qdG4MdzTaDBpDCpHiRu3AA3m\naNSIheF24jrd0IB2e6tf096G3DhvLrC1wDMCD0p1iy3rjRu3AAki7W27Is1uEJqN7sDwiK4zrINj\nrA3VgMBwgesE3hLY12oUJJLMtm0zCEa9sDbUCQR6CfzUS0J3jlQ/I8xoCDIQZFeQ35Lhtt3sBsGN\nW4AE4UZ8vbS3IbdRNxLYXuBFgbtFZ24lETduARqL9AT5Kkge5BGQD0D+BXIy6W/bFTGDYPi4EV8v\n7W3IrfcNBNYWuFngNYHdE+4ecuMWoL6IA7IhyDEgd4AsBZkGch7I9iD9ggfHJmYVDEXTMnwEvIZO\nCxwE3IWu3H2YtjNLfJrdIBj1w9pQBbwaBacJvOv9azUKYkHWADkYZDLIApA5IFeDfBtk1fZObJiI\nNXALOmd8JXTWx+ro9MPr0Pnp56NzmMMwg2DUC2tDZXhJ6HYXmC1wi5QubjTqjqwEsjvIpSAvgrwL\nchPI4SCV1iSFXqhuInaRYehiqN5lv0+nOFVtMJr2OIxmNwhu3AIkCDfi66W9DblRXkxgXYG7vFhB\nFEkQG40btwCdR3qBbANyJsi/QZaD3AdyEshEkO4dXyP8wpGKGSFfRTN43oOuqp2CjgoWoyMGn/cJ\nz1vf7AbhBTQtdJy46OK3RnEa8PMKcnSGoWg+nfLOiE/a25AbxUUEBnizht7xZhGV15VIC27cAnSM\nOCAbgxwHcifIMpCnQM4B+QZIVK65xLZtF02xsBM6Te064BzCDUJYQ0yjQZhJcVXt22ha6sFVnPcr\nNE1wkMeAnSOVri0/pHLVtCjz1VfDIHTFcX9v+0vAfzs4x0Vl/VnIvt+hi93CSHIbqjuee2hfgTcF\nJktjq+M1EbImyCEg14MsBHkV5EqQfUA6k8yyUzet03X/R62pa98CFqEjBNC0zEd4v7egq0mHUMzY\nGcYkNKsnqOGYHtjnev+2Jmi7H5p47q9oUZA8agSPqOL8Nbzf/O1BaH6eu+oo70sUU5P8HK1lsFpg\nv0vj3t8v0VGRn25iIqUzW8LOPwlNnX0Qms4juP95tAbFpQ2SPxXbAu8Cl90JI++E868odkQSIV+6\nt9cfAC/2ALaDu3eH+1eC7e4B7oOd/gn3LqzD/f3vLaSAGWhxmIHAZDSgfKb3fTCas6c8L41PGkcI\nz1GahvkktKc/EVhIaV6ovVEDtwulPXNfsT2P5jx6gqLLLVg0Zyf0/S5DZ2sFUwg83sG5YeyK5g4K\n4lLqMlqIZix92bvuuaghK6BFgW5HjaLPl4BHPBmno5XiKjGF0noLwRGCG3J8f++6W6CdhYll+3t4\nMpXXt4Zkt6FqcDt7gsDKApcIvC1wpJQmREw7bjy3ld4gXwc5C+RxLw5wL8iJIBNA4sgDl+i2vQVa\nQWw58HfUdeRPO/0EVWSV0gyk1SAc4H1fG1XIV3rbM1Hj6HMbqrRBS2heSikveOeMQ43no8DJ3r51\nUGW3N/pOT0KVeR9v/+OeLGHnVqIag7AA+BeaYmIDVBE/girvwd59fTfNauikgv3QqYs7oAalkgvt\ndUoznnZkEL6Lti3QinDl7w80hvWtkN+T3Iaqwa32QNEaBYcKLBC4SkrrZmQFtzG3kW4gm4CcAHK3\nFwd4AuTXINuC9On4GnUnsS4jgP+giqOcuvnGC0RTQjNHTSU0HXT0Mxkdmk+lqIivRZXY3airbAfU\nleSfV36/pcA1qGEAjcFs533fB1XMftnI84AfANuixlZQBRl2blc5BzUMC9CSiM+idZZBi9ZM8L7v\nhxqiG7ztqWg5y13Q91POYEqztgZpDfntYOBv3ve/oe7F4yktxrIUnciQNVqrOUi09Orl6MhzN6f4\n/5Q1Wut3aVkL/dvZDvgG2kbvQ1N9HwDOkvrdO5mkqqZyjYo8KoRiDKGcvwCnoj36fVHluChwXhjz\nAt8/pJhDZgSlVdEA5qDum47O7SrlVdOC2+VV03ajNB2y0La0ps8yigHljlgT7RX6o5H7UBfILqjr\nyWcAlY1MZhGdZfVrdNR3MnBdTGmpU4isgqZY943Aymjn637gFHBej0+2ZGD1EKJhAdpj3hs1GsFe\ncphBGNTOtebRNojUQqkRaBSVDPAbaLW4boFPd1RRhTETGFthn1u2/V3vev6srhXoH+7BZXKNIZtV\n09ywHwV6iNbZeB51Ka7vwKQmMAZu7adKb5AcyNkg/0E7Wj8AZqN1NYaBsx84fzJjoJhBiI5JaI9t\nHHBr4Pe3gfWofh74rejwdU80WHwiGj94ICpBI+BGtId1INpTXxld9BQW5IXKVdPCOBiNvwSNzdao\nK9KPUYxHRwev1CB76hAN2E9DOxzbOnCcYzWlQ5BuIJt6gd970TjXOWin7CRgKDjfBOdCcGaAk3Vj\n2mnMIETHFDTYeisaVPe5iaJ7w59lVO7qEIojidfQOMIZqMtmD9Q98AnhBM9tj7BjOjpPyr772wvR\nus/fQ6cazwJ+0s51JqMun6BR7I72bh+gOAvrJNRldG3Z+Y+jyn9fb3tftHJcFmn1vwiMEF3v8lfg\nbOAbTjZHRe3R2v5uaQH5PsiNqJv2b+iI+gpgTXC2BOdUcFrBWVFXSY2aSeMso45wUF//tnELklAu\nQhfLdZWBqHGo5HZLcxsC8GsUnOitMj5bqo+/NAEyGGRvkCu8xWCLvMVh3/OCxFkm9W27ElkzCA5w\nKNW7MNz6iZJYeqMB0XLcGq6T2oyQHXG8ppl4SeBO0ThJkzN0By/9wzkgT3rTQe8COd5LF5Hk1N1R\nk+q23R5ZMwhPoW6UapOHufUTJXW4EV8vlW1IoEXg1rthnqiLsEmRbiCbgfxME8Ld/xGaIO5MNGFc\nWnMyRUEq23Y1ZM0gGMkhVW1IoK/A6V6NglOluACxiZBRaCrom0DeQVNEXwayB8hKHZ/fNKSqbXcG\nMwhGvUhFG/KS0H1LYI7A36XyDK0MIkNA/g/kKpDZaJGYyWjRmJFxS5dgUtG2a6HZDYIbtwAJwo34\neolvQwLrCdwj8ILoFOMgbhwy1RfpC7IdWhZyGlom8g6QY0E26iAO4DZKyhSQ6NQVhmF0AtEZUqcC\nh6HTSC93Ok4DnkKkO5rmxF8RvAWaBuV+4BjgCXAanYLdSDCVLN0SivPd7WOfWj6Jyz8j6h7aX+At\ngWulctLHlCIOyGiQI0BuRktEPg9yCchuIO2tzDeqR+IWoF5k9sEMI4jAeIGHBJ4W2CpueaJDhoJ8\nB+QPaJH4+SDXgRwEYkV56kNm9WZmH6xK3LgFSBBu3ALUA4FVBC4TWCRwhFRfo8Ctp1y1I/1AdgA5\nH+QZNA5wO8jRIOOo33oAt07XTSN115sWQzCMCBFNB3MoWjr1H8A4R9OlpwzpjhYm8uMAmwPPoHGA\no9A4QAbjH0YcNPsIwcggApsLPCHwb4HN4panc4gDMhbkRyC3gCwBmQlyMcguIB1V5TPqT2b1ZmYf\nzGg+BFYT+JPAfIGDJDVJI2U1kP1A/gQyF2QeyCSQA0GGxy2d0YbM6s3MPliVuHELkCDcuAWoFdEa\nBUcLLBa4UCCKVbVuBNeogPQH2QnkNyDTQd4HmQJyFMj6dYwDdAU3bgEShMUQDCOJCHwduAytd/F1\np1jSNEFID4pxgO3RWtbT0DjAj4Anwfms8vmG0RiafYRgpBSBkQI3CMwV2EcqV5WLAXFA1gP5Mcg/\nQN4DeRbkIpCdQaIqtWrEQ2b1ZmYfzMgmAr0FThKtUXCWJKZGgawOsj/INSBvgLzpfd8fJGML4Jqe\nzOrNzD5YlbhxC5Ag3LgF6AiBnQReEbhdYHSdb+d2IM0Ar7d/kdf7f88bDfzYGx0kaMQSCW7cAiQI\niyEYRlwIjAIuRutkH+PAXTFI0QP4MsX1AJuh9TfuRwvGT7M4gJF2mn2EYCQYgX4CZ3juoZNFq7Q1\n6u6ON+PnJ94MoPe9GUG/8WYIJcRVZcRAZvVmZh/MSC9eErq9BF4X+JvAmg2683Bv7v8kkLe8WMCf\nQPbVtQKGAWRYb2b2warEjVuABOHGLQCAwPoCUwVmCuTqfLeB3urfi73VwEu81cE/gt0OyGAcoCu4\ncQuQIDKrNzP7YFXixi1AgnDjvLnAIIELvMVlxwj0rMNdeoJsDXI6yMMgH4A8AHIKyOZe3iAfN/r7\npxo3bgFiyiDQAAAe8ElEQVQSRGb1ZmYfzEgHnnvoQNGi9tcIrB7h1R0vA+jRaEbQ90GeRjOF7gDS\nL7p7GU1EZvVmZh/MSD4Cmwo8LPCUwJYRXXUNtBbAdWhtgNfRWgHfARkazT2MJiezejOzD1YlbtwC\nJAi3UTcSGCzwO69GweGdqFEQdrVBaDWwS9DqYO+C/B2tGja6C3EAt3aZMokbtwAJwtYhGEZX8RT/\nYcBZwM3ABk6nS21KT7Q28PboeoDxwH/Q9QAHA8+A83l0UhtG89DsIwSjQQhs6bmGHhHYtBNnOiAb\ngRwLcgdaIWwayHkg24H0rZ/UhhFKZvVmZh/MSAYCqwv82QsaH1hdEjoZCXIwyGSQBSCzQa4C+T+Q\nVesvtWG0S2b1ZmYfrErcuAVIEG6UFxPo6U0fXexNJx3UztErgewBchnIiyDvgNwEcjjIqCjl6gRu\nTPdNKm7cAiQIiyEYRrV4C8ouA+YD2zjwYtkRvdBZRX59gI2Ax9A4wAHAdHC+aKDIhmFgIwQjQgTW\n9FJNvC6wZ9E9JA7IeJDjQe4CWQbyJMg5IN8A6ROv5IbRKRKvN7sBjwMPe9uD0IyQH3u/VVrsk/gH\nM5KPV6PgFC8J3RkC/UDWBPkeyPUgi0BeBbkCZG+QwXHLbBhdoO56s6vFwI8EXqEo6E+Bd4Dh6FD8\nzC5eP6u4cQuQINxaThLYGZj5Cb233oe//8JBVnWQZ4CngR2BArAlOGPA+RE4t4DTyammseDGLUDC\ncOMWwKiOEcC/gK9QHCFMR+dnAwxG682G0ewjBDduARKE25mDX2bsBvMY/ugihr77bW58EWQ5yL0g\nPwXZFKSrnZy4ceMWIGG4cQuQIBKtN29Ei3ZvQdEgLAZWChzzPuHJwhL9YEaSkG4gmwxl0c+v4IhX\n32HwF2dz8hsDWXouSM7iAEYTkVi9uRNwpfd9S9o3CL1Czk/sgxlJQNYGOQzkBofP3z6Q6+YvZsjy\nlxn70GQOGBe3dIYRE4mddro1cLj38XkSeAtoAWYAQ4BPvU8Yk4DXve/vo+6mVm/b9f7N6vaxNNfz\ntrftwrgBsNcEOGs4sD3ctyosmfZtKFzDoWs+yEcjzoDzLufdS2BW3PLWe9v/nhR54t72vydFnkZu\n+99bSBFBl9GZwGQ0fnABcHWFc5p9hODGLUC8SB+QbUF+Df980ZsOejfICSCbvM5aKwlc6C0u+4k0\n13oZN24BEoYbtwAJIhV6c0vgIe+7P+30E9RIDKtwTioezIgK6QYyAeRELwC8HOQxkLNAvg7SG0Cg\nm8BBAvMF/ihg5SMNo0hm9WZmH8zwkRaQ74P8DWQxyEsgl4N8C2TlNkfDZgKPCjwhsHkMAhtG0sms\n3szsg1WJG7cA0SNDQPYBudJbDLbIWxx2iC4WC2c87C5whcBCgcOk62tj0o4btwAJw41bgASR2KCy\n0fRIH3RygZ8XaF3UTXg/8DtgplaprHC21ij4wX3wa+B6tEbBe3UX2zCMxNHsI4QUIt1AJoKcBHKf\nFwf4N8iZINt4ieOquxJsJfC0wEMCm9RTasPIEJnVm5l9sGwho7xU0Dd5qaFfALkUZHdNHd3Jq8Ew\ngWsF3hLYr7oaBYZheGRWb2b2warEjVuAcGRVkG+DXA3yGlokZjJaNGZkzVfVGgXHe0nozhMYGNjt\ndl3uTOHGLUDCcOMWIEFYDMGoJ9IX+CoaB9gOGINOIb4fuAR4ob04QFV3gG+gNQreBLZ24OUuiWwY\nRuZo9hFCTEh3kC+BnAzyLy8O8AhIHuSraCH5aO4Eawn8XWCOwLfMPWQYXSazejOzD5YsxAEZA/JD\nkJtB3gV5HuQSkF1B2ikvWeMdoY/ALzz30OkCVozeMKIhs3ozsw9WJW79Li1DQb4D8geQ10HmgVwL\n8l2QEfW7LwjsKvCqwK1Sff4Vt44ipRE3bgEShhu3AAnCYghGR0g/4GsU4wCjgAfROMBFwEtdjQN0\nKAGMBX4LjAaOdGBqPe9nGEa2aPYRQheQ7iCbg5wCUgD5AOQhkF+CbBVlHKBDSaC/wNmee+hECU91\nbhhGNGRWb2b2waJHHJB1QY4EuRVkCchzIBeD7AIysONrRCyRDjm+LfCGwF9Eq+cZhlFfMqs3M/tg\nVeK2v1tWA9kP5E8gc0HeAvkzyAEgwxsiYSXJYCOBBwRmiLqquoobwTWyhBu3AAnDjVuABGExhOZA\n+qPKdXs0DrA2WizjfrSuxMv1jgN0hGglvDxwAHAGcJUDn8Upk2EY2aDJRwjSA2RLkFNBWr04wIMg\np3m/J8ZQezUKDhFYIHC1wNC4ZTKMJiWzejOzDxaOOCDrgfwY5DaQ90FmgFwI8k2QAXFLGIbARIHH\nBB4X+HLc8hhG05JnFcxllGZkGJq2wZ8OKqgL6CbY4lp44h9xStceAqsCZwO7A6cA1zrwRZ1u51Ks\nJWvY+yjHpdneR57BwMSyT0NG5mYQIkMGANtQrA8wEiigRuAcYFYgDuDGIWFHeDUKjkBjBX9FaxS8\nH6tQhpFl8gyhrfIfAjwDTANuA04DZtGAmF1c+WUkxntHhPRASz36I4DNgCdRA3A/MA2c1ARdRZPc\nXQYsBX7iwHMxi2QY2SLPUNoq/5UpKn//M4t86Ii87nrTDELViAOsT9EAfB2YQ9EAPALOh/HJVxsC\nw4Hz0VHLT4GbnKaL8RhGxORZjbbKfxDwNKXKf3YF5R+GGYR4keEUDcB26JDtPtQAPADO2zVe2CVm\nv6i3qvho4OfAH4CzHfggBlFcms1H3D4u9j6CuCT9feQZRlHpb+b9O4C2yv+1Tij/MOquNy2GUIIM\nRHv+vgEYgcYB7gPOAmbHvR4gCkRjHJeiI5ytHHglZpEMIx3kGU7bnn9fikr/euAEVPmnTlc0+QhB\nelKMA2wPbAr8h6Ib6GlwPo9PvmgRXfB2ETABOAa4w9xDhhFCHgftEAZ7/ROB3pT2+qcBrzdI+ZvL\nKOLbOsA4iiOAbYDZlMYBPmq8XPXFq0lwIuoiugS4wIFP4pXKMBKCKv81aNvz705R6fvun7kx9vzN\nIHThFr6F39T7bILOpFlBMQ5QAGdxfeUIxaUBflGvStluaGrqp4ETHJhb7/t2Epek+4gbi4u9jyAu\nUb4PVf4jaav8Hdr2/N9MmNvHYgjVIT2BDVClv2ng3y+A6d5nCnAyOLPjkrKRCKyLjgZagMMdNYCG\n0Tyo8l+LUsW/GapYfaV/tffvWwlT/rGQwhGCrExR4fvKf3205zuDogGYASzMQhC4M4jObjgV+D66\nIO4yBz6NVyrDqDOq/Nembc//v7Tt+c9PqfJvZpeROMA6tFX+Q4BnKVX+M7Po++8MnnvoO2h21AJw\nkgML4pXKMOqAKv91aNvzX0G58s8zPy4x60CzGATpA2xEqfIfDyyj2Nv3lf9r4NQrr06jcInQLyqw\nMbrKeCXgKAcejeraDcDFfOZBXOx9FHFwOZ03aav8P6JU+T9NPvMdoCzHEOREir3+UWiuDl/5/0P/\ndd6NT77kI7rs/QxgP+B04GoHMjNN1mgy8nRD63IXp3m+xuZoPi1f8V+E9vwXxSZnholxhCC/pdjz\nfxGcFTHJkjoEugGHAL9Gk1+d6sA7sQplGJ1Blf8YSnv+E9BcWqVTPfPUmhEgazSLy8ioFq8uweXo\nDKqjHP3DMYzkosp/XUoXeE0A3qOt2yeOaeBpwQxCRnHppJ/Yq1T2a2BXNP/Q5DrWKGgkLuYzD+KS\n5veRpzuq/IM9/03REWxwgdfT5Ksa1bqk+X1ES5ZjCEY1iP4f/RD4JfAXYH1Hh9WGES+q/NejrfJ/\nm2Kv/0xU+S+JS0yjemyEkGAEvoa6h94Bjnbg+ZhFMpqVPD3Q9T7BmT6bAIto6/Z5Ly4xM465jJoR\n0bwq56MG4QTgZktCZzQMVf4bUNrzHw/Mp1T5P0PeKuo1ggKFPjlyH5NQl9EIdMn319Ac+ueh6ZQH\nATcCOeApYB+w6WEhuIT4Rb0aBccCPwOuRFNOpK7oTidxMR9xEJdGvo88PdGEj8GsnhsD8ygq/ltQ\n5R+Hq9KlSdpHgUJvdAr+WO8zJvB9WCNkqNUg9AYmAwcAawIPAP9CV8q+g1bhOgX1Hx7RdTGzj8CO\nqFGdBWzpwKsxi2RkDVX+G1La898IeIOi8v87qvyXxSVmlilQ6IUq/aCy95X/CPT/Ypb3eR7NwTYL\nTc3z33rLF9XwYypwLrpo5CA0tcRg4CVgtZDjzWXkIboE/yK0V3asA3fELJKRBfL0oq3y3xBVLEG3\nz3TyLI9LzCziKf0W2ir8sag7+E2KSv/VwPe5OXLtKf1UzDIajSaVegx9WD+98hLUBdKTBli2tOHV\nKPgZ8BPUIOxnNQqMmsjTG+3pB5X/OOA1ilM9r0eVfxxlUjNHgUJPSpV+sMc/EniLorJ/BbjT+/56\nB0o/VrpqEIYANwGHAR9XOKaSRZsEvO59fx9dsdzqbbvev5nbFnB2hj/eCzvvCA8Dmzk6hNwyCfLF\nsO1/T4o8cW/738P35+nNbRzCaqzLVgwEJjKHDfmM+YzlIWAa/+RZZjGbZdwTOL8HxZrZjXyerm77\n3xt+/wKFR4CWq7hqryEMWWMf9nGAsU/wxPhneGboBCbMA2ZNZeqHy1k+b2/2vhyYtRu7rfkBH3wW\ncv1ZnZTH/95Cg+jK8GMQcC/qKpri/TYdOBhNSTEEeBFzGf0P0Tnbl94N634TDnM09tLsuDRJ0LBK\nXPz3kacPOrsnONVzfbTnGXT7zCBPVrP9utSxfRQo9EBrJoS5d9YCFhLu3pmTI9fodDuJnXbaF7gH\nuBa4JvD7mahP/BjgZDT75uEh5zeVQRAYCJwGfA9dbXy5Y240I0ievpQq/4noit9XKFX+z5KvOBo3\nQihQ6E6p0g+6d9ZGZ0IGlb2v/F/LkUuSGzexBmE7NJAc5EA0IHojsC3wJPB/qIUtpykMglejYD90\nTcH9wM+d8PdhNBN5+qGLuoJTPccCL9NW+SdJISUWT+mvSVuFPxZ1uSymrcKfhSr9tBjYxBqErpJ5\ngyD6B38Z0B9NQvdYYLeLuUl8XLL8LlT5b0ppz3806k4NKv+ZnvJ3yfL76Dwu3vsoUOiGBmzD3Dvr\nAO8S7t6ZnSKl3x6pmGVkBBBYBTgLHR39Evij1ShoEvIMoNjz9z+jgBdQpf9vtJMwkzyW7j0ET+mv\ngafw7+bu3Df55rHe9ih09mJQ4T9KUelnNY7SMGyEEBFejYJDgV8BtwKnOdpjMbKIKv8JlCr/FmAm\nxYyefs/faloH8JT+CNq6dnyl/z7hPv1Xc+SyvnK/PcxllAYENkeT0P0XdQ89E7NIRpTkGURbt89a\nqPIPun2eJ2+TBQAKFBw0Y0GYe2c0sJxwn/6rOXK2ViIcMwhJRnRK7TnATmiNgr9UmYTOxfzEPi5J\neheq/IOFXCaifutnKe35v1An5e+SpPfRDp7SH0b47J0x6LqHcn++796pNjWGS0reRwOwGEIS8WoU\nHIlOJb0W2MDBcr+kjjwr0Vb5j0CV/zR0Jt05wIvk+SwuMePEU/qrEz57Zwxa7D6o8P9O0b1jdTtS\nho0QOoloj+UydPro0Y7OFjGSTp5VKCp//9/h6CLKoNvnpWZT/p7SH0q4e2cM8CmV3TuW/rpxmMso\nKYi6DS4AtgKOB261GgUJJc9g2vb8V0NX0gcLuL9EvjlmgHlKf1Uqu3c+I3zK5qs5clbwJhmYQYgb\n0VTfxwE/BX4PnOvQ5TQBLuYX9XHpyrvIM4TS1A4TUcUXVP7TgFdSovxdanwfntIfQmX3zheE+/Rf\nzZFLaolLF/tb8bEYQpyIBosvRd1CmzuaPdKIizyrUtrrn4iu+3gGVfr/QOM6r5Dni7jErDcFCkMI\nV/hjUYURVPZ34hmAHDmbBm20i40QQhCdC30xmkL4GAfuilmk5iPPUNoq/5UpnekzDXg1i8q/QGEV\nwt07Y4HuhLt3ZgHv5siZKzObmMuokQj0Q6ePHgn8BrjYwVaU1p08q9NW+Q+krfKfnSXlX6CwMpXd\nO72o4N4BFpvSb0rMIDQCLwndnmihmseAEx0tcFEvXJrVL5pnGEHFP5utGE0PShX/NOA18ukP2hco\nrERl904fygK4F3PxgOM47ibgbVP6QDP/rbTFYgj1RmADNE4wDDjEscYXHXlG0Lbn35tiz/8vPMaN\njObGNCv/AoWBVHbv9CMwTRNtX3/wtheFKH33dm5f1BjJDaOUph0hiBb4+SVa0Ocs4PcOzTX/PDLy\nOBCq/HvStuc/N43Kv0BhAJXdOwPxFmPR1r2zwHr6RkSYy6gON3aAA4Dz0IpvJztaIMOoBlX+Iymd\n5jkRDXSWK/830qT8CxT6U9m9sxIwm/AFWvNN6RsNwAxCxDfdFE1C1xtNQvefRsvg4ZIG15Qq/zVp\n2/MXShd4TQPerFH5uzTwXRQo9KN0FW5Q+a9CqdIP9vjn58g1IqDtkoa20Thc7H34WAwhCkTTEp8K\n7Or9e41DdmarRIIq/7Vpu8L3c4rK/0rv33lJ7vkXKPRFlX1Yb38Iup7EV/hPAn/1tuc1SOkbRiLJ\n9AhBNEXxKWixmiuACx2wZfiq/FsoVfyboTlrgr3+acD8JCr/AoU+aBrlMPfOUGAO4e6dt3Lk0rBi\n2TDKMZdRjRcfCZwM7AtcjRqCd+p1v0Sjyn8d2ir/Tyj3+edZEJeYYXhKfxTh7p3VgNcJd++8aUrf\nyCBmEDp50eGoITgA+BNwgaPFtZOGSz38oqr8R9NW+X9IW+W/MPL718AgBm0/hSlvEe7eGQbMJXyB\n1ps5clmcFeZiPvMgLvY+fCyGUA2i+dpPAg4B/gyMy/zMoTzdCFf+yygq/t8AT5OP910UKPRCRylt\n3DvP8MwaaE/fV/jPA7d523MzqvQNI5GkeoTgVSw7Ea1lPBk4zyFZbo9IUOU/ltKA7wRgKaU9/6fJ\n83YcIhYo9KRU6Qd7/GsAbxLu05+bI2dlJ43moFDoia5bGeR9gt/b387ltsJGCG0RTW/8U+AH6AyR\n8Q7Mi1eqiFDlvy6lPf9NgSUUFf+5qNunoXERT+m3EF4cfSSa7iOo7O/y/n3dlL6RWgqFbkB/alHi\nbff1REfxy9C60ssqbL9Vtr0ceLjej5qqEYLAYOAE4IfATcCvHe15pg0XaCVPd2A9Snv+m6IB8PKe\nf0NSFxco9ECnn4bN3lkTmE/4itw5OXKf1nBLF/MRB3Gx9xHEpZb3USg46HqjWpR2+XZ/tAZKR0q8\nvX3+9sfkal7EaDEEANEFQ8cBPwZuATZzNNiYHlT5rw9MZDq7sym/AjZBYx3+VM8zUeVf12IlntJf\ni3D3zlpoedCgwr/P+3dOjpxlfzXqR6HQg6JCHshFF23I8cf3oXO9cn/7C6pT0os72P8BueaYtZbo\nEYJouoBjgZ8AU4BfOTq/PNnk6YGv/IufTdD4RrDn/wz5+qyLKFDoTlHpl7t41kYNUdiUzddM6Rud\nQnvj/alNaZdv96E6Jd7R9nJymWvHzTlC8BLPHQ0cg1Z82tJRpZU8VPmPo3Smz3jUteIr/n+gyj/S\nguSe0h9JuHunBe35BJV9q7c9O0fukyhlMVJIodCbaJT4QHRdSzVKu9w3Xn7sR11wqRhdJFEjBIEB\n6GjgODTx3FkOvNJg2SqTpyelyn8isDEaxyjv+S9r50ouVfpFCxS6Uar0g739dYB3CZ+9MztH7uPO\nPF5MuJjPPIhLe++jUOiO/p10VYkPQv8Gu9IL979/QK5u04NdrH34NMcIQXS4+WM0YPwA8HVH6xjH\nR55ewIaU9vw3At6gqPj/hir/5V25laf01yB89s4odIZRUOH/m6LS/6gr9zYagLpU+lKNkp46dQN2\n2OFQKgc9+wIf0LGSfhtN1NdeAHSF9caNILGOELySlT9C1xI8BJzh6MKkxqLKfyNKe/4bovGKYM9/\nRq3K31P6wwl374wG3id8Re7sHLkPa300owsUCr2ofjZKR/s+JZpZKh+SswR8TUp2U1eIuoV+hpas\nPMOBZxty5zy9aav8x6EZMIPKfzp5OqWICxQcSpV+sMc/Gv2jDpuy+WqO3AddfTQDf8540KVS2yIg\n/XQnugCnrcMwukqmDcIUIO/A9Lrdpaj8v4Qq/i+hs39epTSr54xqlb+n9IcRnntnDJo3KMyn/2qO\nnD+6cDG/qI9LofAgOrukK/5w/9OP4pzxWpR48PsnMbhUXKxtBHGx9+GT3RiCA9+K9IKlbh/fAIyj\nqPyfAq4BniVPu353T+mvRrh7ZwzwMaXK/u8UlX57weRsUVyGX7sSnzZtMOoX/5zqlPaidvb5AU5z\nqRhGDSRqllHV6GyfDSnt+Y9Dff5PUTQAMyopf0/pDyXcvTMGWEG4e2dWjtzSmmWPm+Iy/K70xP3v\nvYlmlspycjWtcjaMZiK7LqOq71061dM3ABuhGTJ9xR/q8w8Eckd5n3Kl/xmV3TvJKaRTXIbflRwq\n/vYAdITTVSW+jK4twzcMo3Ok1iDkUPfMEOBStGxlkPAH00VeG1Da898YneoZ7PlPJ88HAAUKA9D5\n+KNCPi3o7J3XvE9JvdwcuSXRPG4FdM54W8U8adKWHHLIW1SvxAd5V1xK12epLE/YMnwX8xEHcbH3\nEcTF3odPKmMIDmoMjkYV+L/Q1caPlRyluX02oLTnPx5dyegr/pvGLBgz4w9X/WEQRSW/IzpV1d8e\nSFHhz0F7+VP97U7P09feeD9qm15Yvu0vwy9V2itWrIYGtH0l7aeRqKzUs7cM32dT7A8+iL2PUux9\nNJB6GIRN0V75P73tPwB7UW4QtLc7H3hq6NKhM3eYscMf93hyj0+GLh/qu3h2QLOaro2uxn0t8Lnb\n/z7pYBZdewh9UVdI8DMCWBcdQfifYtKs9hX+CqpznczrYP+HFVwqeW68MV/tC804K8ctQMKw91GK\nvY8GUg+DsAalmUjnoCOAEm7+7dRH+67ouUavT9nFEfb8sD9vL12JJU+PYen8EXwwd22WzxrLk7PG\n8p+P+v9P4Y9CRxED0MCo/+8KdLrnBxU+wX1z6TjAaVW6DMNoOhox7TTU53XN93tuvGA4S98aycx3\nVuU96RaqvNtT7MXfkuUTr4aWuAVIEC1xC5AwWuIWIGG0xC1AM1GPAMWmaF3jCd72CWjN458FjpmO\npoM2DMMwqmM2OjsyVXRD3US7o+6jF4GtYpXIMAzDiI1tUaOwHPh1zLIYhmEYhmEYhmHUnxHAHej0\n2XnougvQqat3oauBH0bjJz7HoVNo56FTcX02Bp5Dg+Z/Qt1uAD2B69BkbTPQxHxJphvwOPrc0Nzv\nYihwOyrva+jq+mZ9H8eiz7UcrR7oZ6NtlndxA1qf4rnAb416/n3R0r3vAkdG8jRGKOsA30HrPm+E\n/odvCJyJ/uesDJwPXOUdPxpddzEW+Ar6n9TH2/cQ+p81BF2Es5/3+6HA/cBgtJLcvXV8nig4Cn32\nh7ztZn4XtwCXoO1jHfQPvhnfx1DgPfTZBgC3oaVxm+ldbIVOsgkahEY8/0B00esWwHredUdG+WBG\nZaaisZPp6FoJ0P+gt73vxwEXBY6/HdjFO+adwO97odlTQVOE7+597+4d1z9qwSNiBLoq/SsURwjN\n+i6GofL1Lvu9Gd9Hf3SEtB6qoG4D9qb53sUYSg1CvZ9/ALAncGvgnEvRypQV6dbeTqNqRqMrqh+j\ndGHeEqCX9xmB5mTymeMdOwJN1+Hzuvc73j7/Wp+jQ8gRkUsfDRcBJwHB1NPN+i7GoM83BXWTTEF7\ngs34Pj4ETgFmou7VL9DRUzO+iyCNeP5K16qIGYSuMwS4CTgM9QeG0d56jyxkC90JTVfyFB2vbcn6\nuwBd8Dke+C2abXcpaizDyPr7GAZcgOYqWwXtwR5V4disv4uOiP35zSB0jUFoYPlM4BHvt3kUV1cO\nQWvprvB+Xztw7ijU4i+g1K+3DsWeQPBa3VHrPj9C+aNia+BwtPf3b2/7SfQ5WrxjmuVdgMq8CLgH\nXU1/A5rIsRnfxyZoedxnUcN4I+pWbLa/k3KFXu/nn9fOtYw60Bd4EA3oBDkTmIz6/C4ArvZ+H4M2\n0nXRINNCSoNFPwZW9a65v/f7oahffgg6i2lqHZ4jaragGENo5ncxA/gm6jefjL6LZnwfo1D/+MZo\ngP024Bc037sojyE04vkHou9+S3TmkQWV68h2aI84+Nmf4nSyT1DFOCxwzvGov3A+GljzGY/6WD9C\nU4f7Q8ceaKNJ6nS6MLakOMuomd/FFugq/eVo8C841bLZ3seRaM90Kepe7U9zvYtHKNUTR9O4598P\nHV0sobKrzjAMwzAMwzAMwzAMwzAMwzAMwzAMwzAMwzAMwzAMwzAMw+iYMeic8CUdHPe4d9wudZfI\nMGLGUlcYzcyj6ErR9tgSzT1vGJnHDIJhqNK/C7gZXVl8bbziGEY8mEEwDGUbtKDN2miKgC3iFccw\nGo8ZBMNQ/o3mlFmCJgpriVUaw4gBMwiGoXwY+P5fNGGYYTQVZhAMwzAMwAyCYfiUFzBphgpdhlGC\nDYsNQ9ca7BXYPrlsf0dlQQ3DMIwU05mFaZ8DO9ddIsMwDMMwDMMwDMMwDMMwDMMwDMMwDMMwDMMw\nDMMwDMMwjOTz/2xHqT7qzDQbAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# add results for plotting:\n", "# Convert to pandas dataframe for plotting:\n", "matlab_la_data=pd.DataFrame(mat_time_la,columns=['n','Matlab Time (LA)'])\n", "python_la_data=pd.DataFrame([[1000,.309],[10000,2.12],[100000,21.7]],columns=['n','Python Time (LA)'])\n", "plot_data_2 = pd.concat([plot_data,python_la_data['Python Time (LA)'],matlab_la_data['Matlab Time (LA)']], axis=1)\n", "#plot_data = plot_data.set_index('n')\n", "plot_data_2.plot(x=['n'],y=['Matlab Time','Python Time','Stata Time','Matlab Time (LA)','Python Time (LA)'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Discussion\n", "\n", "The linear algebra model run times for both Python and Matlab are denoted by LA. We add them to the previous figure.\n", "\n", "\n", "The python results are very similar, showing that the statsmodels OLS function is highly optimized. On the other hand, Matlab shows significant speed improvements and demonstrates how native linear algebra code is preferred for speed. For this example, Matlab is roughly three times faster than python.\n", "\n", "\n", "Stata was dropped from the comparison because of lack of support in Stata's linear algebra environment (Mata) for sampling with replacement for large \\\$$N\\\$$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Parallelizing the Bootstrap\n", "\n", "All of the results above are run using default settings with respect to multi-threading or using multiple processing cores. Matlab and Stata automatically take advantage of multiple cores, whereas Python doesn't. The following comparison manually creates worker pools in both Matlab and Python. The current version of Matlab requires the license for the Parallel Computing Toolbox that [supports 12 workers](http://www.mathworks.com/help/distcomp/parpool.html) and to get more, one would need to purchase and configure the Matlab Distributed Computer Server and the price is conditional on the number of nodes (or roughly speaking, cores) one wants to use. To get any multi-core support in Stata, you must purchase the MP version of the program.\n", "\n", "\n", "Here is the Matlab code starting a worker pool and running the bootstrap code:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Starting parallel pool (parpool) using the 'local' profile ... connected to 12 workers.\n", "Parallel pool using the 'local' profile is shutting down.\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%%matlab -i reps,beta,n_array -o mat_time_la_par,store_beta\n", "\n", "if isempty(gcp('nocreate'))\n", "parpool;\n", "end\n", "\n", "mat_time_la_par = zeros(cols(n_array),2);\n", " \n", "for i=1:3\n", "n=n_array(i);\n", "row_id = 1: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", " parfor 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_la_par(i,:) = [n toc];\n", "end\n", "\n", "delete(gcp('nocreate'));" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following runs the bootstrap in parallel in Python. Note, when passing the n_jobs parameter to the Parallel procedure, one is not arbitrarily restricted due to licensing limits." ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1000\n", "1 loops, best of 3: 409 ms per loop\n", "1 loops, best of 3: 429 ms per loop\n", "1 loops, best of 3: 445 ms per loop\n", "10000\n", "1 loops, best of 3: 488 ms per loop\n", "1 loops, best of 3: 531 ms per loop\n", "1 loops, best of 3: 469 ms per loop\n", "100000\n", "1 loops, best of 3: 2.91 s per loop\n", "1 loops, best of 3: 1.92 s per loop\n", "1 loops, best of 3: 1.95 s per loop\n" ] } ], "source": [ "def python_boot_parallel(arg_reps):\n", " np.random.seed(arg_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", " beta = (np.linalg.inv(np.dot(X_r.T,X_r)).dot(np.dot(X_r.T,Y_r))).T\n", " return beta\n", "\n", "for n in [1000,10000,100000]:\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", " print n\n", " %timeit Parallel(n_jobs=10)(delayed(python_boot_parallel) (arg_reps) for arg_reps in range(reps))\n", " %timeit Parallel(n_jobs=20)(delayed(python_boot_parallel) (arg_reps) for arg_reps in range(reps))\n", " %timeit Parallel(n_jobs=25)(delayed(python_boot_parallel) (arg_reps) for arg_reps in range(reps))" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEPCAYAAABFpK+YAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXmYHFXVxn892RMSsrCGaEIIBAIJCAFk0wofmyIqIApI\nAAFFUZFFBFHhgqgogoqCgMiugAuL8CnwKRkRIawmrCJ7JASSkD0hhCT1/fFW09U93T09M119q7vO\n73n6ma6qrqq3K51Tt9577zlgGIZhGIZhGIZhGIZhGIZhGIZhGIZhGIZhGIZhGIZhpJiRwJ3AYmA2\ncGK0/sfA2tjrIC/qDMMwjLqzKfAZYF1gG2AusDVwAXC4R12GYRhGN+ndyfaXoxeo1T8D2DBaziUl\nyjAMw0gHmwH/AfqjFv9iYAlwG7C+R12GYRhGAowAHgN2j5Y3BQYDQ4ErgGv8yDIMwzC6Si12zRDg\nbuB84PYy2ycC1wEfKLPtBfSkYBiGYdTGi8A4nwIGAH8HjilZ/yHU4TsMtfivrrB/mJy0psD5FpAi\nnG8BKcP5FpAynG8BKSLxuNnWyfbdgD2AKykM3fws8HngVeAlYDhwWoIam5kxvgWkiDG+BaSMMb4F\npIwxvgVkic5G9fyV8jeH3ySgxTAMw2gBsm71BL4FpIjAt4CUEfgWkDIC3wJSRNPHzab/AoZhGA0m\n8bjZmdWTFAtQx7BhpIWFqL+qUQRAewPPl3YC7Hq0DJXuXFl5Egh8C0gRgW8BndDo32TQ4POlncC3\ngBTR9PEx64HfaB7sN2mkBe/DOQ3DMIwWwwJ/sgS+BaSIwLeAlBH4FpAyAt8CsoQFfsMwDKOumMdf\n4Gjgjho/exPw5TqccwywtA7HqZWpwGV1OE5f4Flggzocq1ay+Js00ol5/B54ClgDvL9k/T9QyopJ\nNRxjHD0LuCG1/ePvR3EltPhrDbAKZVFtBG3AWcAPo+X1Ih0jq+yzafSZS0vWrwKuAs6os0bDMLDA\nX44Q1R04IrZuLLAJ8E4XjxXUSVMl7kL/hm0oSypobkYb0At4PeHzx9kLVWh7ucL2oMy6I4EnUJW3\nviXbfhtt71MnfWkj8C0gZQS+BWQJC/zl+Q2yLfIcCdxQ8pk9gEeAZcB/ge/Gtk0HBgH3opb31mXO\ncTmqY7wUPU1sU7J9FMqMugy4h+otZyifYnsMxU8e04FTgUeB5cC1aNLSLdF57qNQYQ305HIXmtz0\nb+DgKuffH5jWicZSvVOBs4F5wAEl22cD84Fdu3BMwzBSQDc9/jCsz6tbPIkC3CPAjtG651EQfJuC\n1bMLyl46ANgWZSs9MNq2GR2tnqMp9viPQ8F8EPB9YGZs202o9bw7KnZzFfB/nejeBtkm8Zv5GDoG\n/pnA5ugJ5iXU4t4LWUK3AhdFn+0fbT8l0rgjeoIYX+H87RQ/JXVm9ewBLEIt/XOAP5X5zC3ASRX2\nrzfm8RtpIasefy5Xn1ePuA61SHdHQfiFku0PAv9EN4OZKDjvkf8CNRz/ShRIlyNvfEsKfnyIrI77\nUXA8BdgTVULrCSHwc3Qjmw3ciQL/X9EN4noKBXX2jdZdFGl8BD31fLrCsYejUpy1chQq27kKuBn1\nV5SW8FyMpfYwjLqT0sDvnRC1uj+FitBcV+Yzk5CVswi1bM+mY2AOqpzDAa+gwLcKefPrxba/Gnu/\nCAXVzuyeWoj7/suBObHlFcA60fvRqN8g3mH8deB9FY67JLZvOYLY+wHo2v4uWn4GeA44vGSfdVDw\nb0UC3wJSRuBbQJawwF+ZeajO8GfQTaCUy5H9Mhpdx3MoXM/OHtX2RbbIvsBAZHespvhJYUzs/TBU\nAjMepJNmFrKG2kpeX6jw+aeovVzcgej73EnhprI1egqIs0V0XMMw6ogF/up8BbVEyrU6B6IW8jvI\n7z+WQsCfF22fVeG4g1CgX4rsnQvRKJw8OeBQ1IcwNNo+DXV21pNqltQ9qKP31EjD4EhPuY7q/Od3\nr3KO9ti6o4CfUXxDeR96wsh3cg9HN79/Vv8KTUu7bwEpo923gCxhgb86ryJvuxwno0lWC5HNEx/1\nswS4AI2eyY/qiY/Nvx14GNkbTwJPU+yPh6gP4HzgNTTCp7Q1XI5yTxrVnj5K5wvEl1cA+6D6yi+g\nm9i5VE7lfScwgeJRQaART3G76ExgCqrVHGc28Bc0ggrgEOAPyI4yDKOJyPrM3cC3gAZzIrpZlSPo\nwnHaUIf5pj0V1AUsLbNfAt8CUkTLFmIxWpNLqE+ahRzq/3ijDscyDKPBZL3FbzQP9ps00kJWx/Eb\nhmEYSWGBP1kC3wJSROBbQMoIfAtIGYFvAVnCAr9hGIZRV8zjN5oF+00a/nFsiI3qMQzDyACOfsDX\ngG804nRm9STLMyhdsU8CNEmsUXyH8gVUgi4eZ32UCrpfTwWllMC3gJQR+BbgBUcOxyfRJM7dURaA\npqcZrZ6nKMwynYty8w+vYb/zUObLOE8DH62ruo58kcpVuFYlfO5ShqAZuIOi5cnAu9H7oMI+AdJa\nrqVzCZoU1ghsApdfAt8CGo5jIo6/4ngaxz6xLTac0wMhSqDWhgLXOOAH3TxWV9IUd5fLKOS7OQAF\n3vxyaVWrpDkU5RQql2ahvcI+R6HU0EeW2XYDcHxdlKWPdt8CUka7bwENw7EejktROvRbgW1x3NNI\nCRb4qzMLFQOZBOyAZpLGr9nBwAxk55yJcvesRXl48oyPlpeiHD3xGrj7odQES1AVrng93+koH1Cl\nfctRLulaQLHV8wbKsPlcdNzzUVGWaagK159Qgrk8k1FdgCXRd92DyuxP1/4DD0LX8HiUj2iHku2P\noOynpfWPDaP5cPTBcRKygN8FtsJxCY7VnpXVnWa0ep4EPhu9H40C72XR8lPAR2KfvQ0FZ1DpxYtL\njvVMtM8EZBf9E/hmtG1TFGgPRnnnT0cJzfpH26dHWsrtW4mPRceIE1Ac+OcAfwM2BrZCuf7vRwF+\neHTevL2yAcoIehjKob8PunFUsr5eoThDZ2dWz1Tg2ej9tXS8fgCPA5+scL56YlaPXwLfAhLF8VEc\n/8ZxF46tOvl0Rkf1uDp9cVdTJaxScqgS1fXAWyjdcD7gXouC1V9Q0ZV9KFgROTq2uEOUfviZaPk6\nVOYQVIjkb8Afo+UfAp9Hlbb+HO17cYV9e8oP0A1gDqqz+wTKJAoqjpKvwnUYuuHcGC3fg+oA74+u\nTyndqcJ1c/T+ZuAaVG0s3gJajNJCG0bzoSB/ETAW/bb/XLf41gPSGvh7WjaxJ+Q9/t+W2XYD8G3U\nQj8UBcE3Y/uVsgR57nmWU6hSNZLiKlsALyPbJU+lfXtKaRWu+HJpFa4DkH2VJ6S4PnCcJRQ6dktp\nL1l+H2rl5Z8u/g/VJNgf2Vp51qExfSWNpt23gJTR7ltAXXEMR+naDwe+B1yKa/hgi4qkM/Cnlzmo\nBXwwujnER/F09S4+G+W6jzOG4mDfKCrdaGeh6mOlJREr8RQq5P5gDZ+divpL4hW2QvQUkA/8OdS5\nblW4jObA0Ru5AGeh/sEJOOb5FdUR69ztOtcg62cC+ofNMxd15MZH0gypcpxbgP9BZQgHA6chf//e\nOmrtKTche+kI1PIeCuxN5c7WSlW4oKOHexTqH4lX4doNDX/N9yFMQq39/3RLfboJfAtIGYFvAT3G\nsTcaAHEQsDeOL6Ux6IMF/u5wO+r0vAVYGVv/Owq2xMNl9oPiClcvIZ//HGS1fAJ1zq4sv2uHalmV\n6GoVrtLt8fO8gfLifw5VAnse+GqV41yPrJr4za8XsorupTC/4HRk9Vxbsv90FOQPjZYPRZXIDCO9\nODbH8Sfgl8gK3gvHE55V9YiRqKTeYmRB5P3YIagD8m00DLG03F6eZhzV0xk55MXv6VtISrkITSrr\nKYPRTaDaU1M9aebfpOEDx7o4foxjPo5vRGkX6oH33+KmwGeAdVER7Lmofuy5aJTJUOBHwOUV9m+1\nwJ8DjqE1rYd60Q+lW6jHcdarw3FqpVl/k0ajcfTC8QUcb+C4MkqsVk9S91u8B7V0Z1CYbDQc3RDK\n0WqB/1Fkf+xd4+eD5KQ0HYFvAZ1g4/j9EvgWUBOOAMcMHPfh2D6hs6RqHP9maHjfg2jIYX4o4gLk\n6fahMFmnVZnsW4BhGB5wjAUuQLPLTwP+kIbx+N2l1sA/AnVeHot8/XJUGhJ4DZrRCZolOiO2LYj+\ntrfocn5dWvT4XG5PmR7fy+0p0+N7uT1lerQ8lAGcxIeAz/MIt/F3vsCy9/Lq1Ot8+fdjSBFDUCv/\nE7F1M4Bto/cjyI7VY7Qu9ps0CjjacHwOx+s4rsUxsoFn9/5bHIBmpx5Tsv5cNHRvOHr8uaLC/lkP\n/IFvASki8C2gE8zj90vgW8B7OHbH8SiOB3Ds5EGBd49/N5SNcQ8K46mPAH6MJve8jjIoHpKUQMMw\njIbgeD8apbgrKiZ0YzP7+NVIOidOWOEcldYbhi/sN5lVHIPQpMIvozQsF+DK1pRoFIn/Fm3mbuM4\nGrijxs/ehH6EaaMd5SmCrn2fa4BTu3G+rpzDMLqGfPwjUInPccAHcDjPQb8hWODvyFPAGjrmo/kH\nSjcwqcMeHRmHipwE3dRQa3oG0CSnfCqE1WiY7bndPG89dZXuN7bM+vOpXDbyP+iGcUB3hDYBgW8B\nKSNo6NkcOwMPoALnn8FxOI5ZDdXgEQv8HQlR0Dkitm4smrvwjhdFtTEK9dkcAnwF5dLvCklnai33\n6HoGhQRtX0UFYfLLWySsx8gijk1wXI9ybf0S2BnHA55VNRwL/OX5DUobnOdIlIs/zh6oY3sZqnr1\n3di26Sgv/b3o6WHrMue4HOU/WoqeJrYp2T4KjahahmZM1zqc7OHo/BOBLaNjLEYzji+lkECtP2pZ\nH4XmWTyAJuH9EQ3PXYzyMb2vxvOOA+4CFqJH54NLtr/Yyf7lCtkcTcHqyes9AqWLXohuFlsDj6Hk\neFeXHGNfVMFrCbqplPt38EW7bwEpoz3RozsG4Pg2qiUxCxgfDdNc28meLYkF/vI8iwLujtHyZ1Fu\nojirgZNQXpqPoZvDgdG2nVGBkzaUnfLpMud4JDr+Rijw/ya2LQccB3wL3QBeo2Mmy1Jy0fl2Bj6I\nqmq1Ad9HTyu7RetL+w4+HunYJzrGnShAjkY3gF90cl5QUL4neo1CN82fozTV9eajaA7JAWgo8YXo\n6WYr4MMUrKGtUeWw01ESwctQZtU+CWgy0oojh+PT6P/0dsCOOL6FY5lnZS1Nt8bxhxDW49VNzU+i\n1upXUenD3VHxFdCs5Uoe/9koMyWU9/iPpnJHZW9kI+WLqd8I/DS2fSh6chhRZt9Sj38WlT3+oyjU\nEMi3oKu1gkejlBx5pqFc41D8fT5Bx6pcPwK+E72/Gj1tVOMr6AYYJ36OvN5xse1PoBtknovRvwPA\nT1DZyzgPUblYvI3j90tQ9yM6tsfxjyi3zofrfvzk8D6O3ws5/8PqQjSyZibKsV/a2gfdAH4KbE8h\ndXC5OrSVcCiwjUT/DiEK4kuj7fGyjIuQXTES1QEuxyiKSyiCbJpL0bjkYdG60uAaP08fFDAPRK3k\nNgpDy6r9GEcja6n0sTmJXPrVykYup1D6cTQq0h6vHxBSu3VlNCuOjVC5w4+ixsfVONb4FZUuzOqp\nzDzkHX8G3QRKuRzViR2NruM5FK5nPki2Vzj2vsir3hcYiHz31RTf8MbE3g9DN5c5XfsK/Aj591tE\n2j5Hx3/zeLA+BiWh2hWlRR5Lee+9lFmoX6Gt5PWF2Gc68/jrzSw0aiiupxflayn7oN23gJTR3uMj\nOPrhOB2NzFsAbBmlTbagX4IF/up8BT2CLi6zbSAqTP4OsAtKYJcP+POi7eWGMIJapatR634w8ql7\nxbbnUPWp3ZDNcyGyWeZ3Uf/ASN8K1Hl8So2fX4r6Ln5Q43nuQU8Ip0Z6B0faG92ZGr9J/Rr4PPAR\n9L1GoCeAdRusyUga+fgHAs+gRssuOE7Dlf1/a2CBvzNeRZ2w5TgZdZQuRL5yfNTPEtTxOIPCqJ74\nGPjb0eib51CfwtPRPnlCZJOcjzp2RyF/vhKVbJjvoLq+C4BfodZuaZnFOFciu+S/aBTMXVWOHf8+\nK1Dn8IeAFyj0M8StxEo3wXLHq7SulhKS+c88iW6eZwNvolbgoTUco1EEvgWkjKBbezkmAX9Do+qO\nx/EJHM/XUZfRDSxJm5En8C2gE6xz1y9Blz7tWB/HZTjexHECLp39ld0k8d+i5eoxDGG/yWbA0RdZ\nsN9EQ6DPwbHQr6i6k/hvsZXukoZhtCqOHLA/6u96EdgDx7/9ijIqYVaPkSfwLaATzOrxS1Bxi2MC\njrtxPIvjI42T5I1sjuM3DMPAMRwNkz4UOA+4FNfydb0bgnn8hiHsN5kWHH2AL6JRab8HzsZ1eShz\nM2Mev2EYGcKxL5o9PhvYE8dTnhW1JDaOP1kC3wJSROBbQMoIfAtIFeM5AscdKCngGcA+FvSTw1r8\nhmH4wzEUOIuXOBZNwvoULtV1L4wayPqonuloCFor8QpKTAdKNPfzGvdrp2OO/lo4uwvn6AlZ+U2m\nA0cvHF/E8QaOK3Bs6FtSirBRPR54CpgQvZ+PErF9leL0xOU4D+WBKc0G2YiAMhmlgIBC+cWfApck\ncK6upFAo3a/c528CPl1hn7+iVBBGK+HYE/0+FwD74ZjhWVHmMI+/IyHKnNmGAuo4ak9WVsqQzj9S\nN9YgzQOAE1FtgF27eIwkGwJDK6w/lEL2zAtRzqP8cisH/cC3gIbjGIvjFpRA7xxgSizoB950ZRAL\n/NWZhQqXTELpit+g+JodjBKx7Q+ciZK2raXQ+gZVoXoYZby8nUKxFYD9UM7/JShPfrzIy3SUCK7S\nvpVYA/wFzW7cBhWSqVQickv0VHMCymt/I7ABSs62AOX+v5HKQbuUySi52xJ0XSoVPalEuRTQjoLV\nk9d7PEq89gaqMfwhVO5xIao4Fmcqytq4GF2XUV3UZPQUxxAc56Pf8qPAVjj+iDN7zRcW+MuTDz6j\nUXCfiXLzz0c59PNMRSUR/xcVfvgFuqY7RduXonTNR0fHWg/lGQHYFPgDymI5EpU8/F9UaSrPMRX2\nrUZvVIBiM1Shag2VS0SC7KmJqHTh51F66F+hbJrboFoA59Rw3vwN4xKUovkbaAx2vgDMohqOUY7S\n4LAuugGMQ9/rCnSD3BOlxz6Bwg10X/S0djSwMao/fGM3ddSbdt8CEkc+/rHoprwRMAnH93GsLPPp\n9oZqyzjp9PinTatPS2DKlO5MgsihSlrXoxbvPSghFCjIT0UtxxHIijg+tl/p+UJU/u+ZaPk6YK/o\n/adQOtk/Rss/RIF3T1TkPESlBMvtW45eFMovvopy708v+cxM4CrUEr81WtcGnAbv1SBdFNO0CAXO\nn1Q5b57DUInKfGC9BwXaj9G1ymS1cAaqG3ATCvyXoCeW14H7UG3VJ9AkoAspPIGdj24W70NPP0ZS\nOPZAPv5K4OM4HvWsyIiRzsDfvYBdL/Ief7lKTTcA30blGA9Fge3N2H6lDEETUfIsj/YFtfJfLfn8\ny6gwep5K+5ZjDeULiXdWInIRFBWeXheVa9wLPWXkqC1IjkaFzuMVvUIKtXhrtYs6YxEUDfcrV34x\nf51Go8B/YYmmNAT+gFZs5TpGo8pvu6CnvptrtHQCWvF6pBSzerrGHNSqPRjdHOIBtKtPKbMpLq9I\ntDy7wyd7RrUSkdCxTu4ZqIN4O9Qw2JPafiezUAu8tNRhqedejXp7vrPQE1mppgfqfB7DsQ6O7wKP\no6fULXHcZD5+OrHA33WuQdbPBNTxm2cu6sjtG1sXr6pVyi2oOtaBqNP2NOTv31tHrVC9RGSlz+fL\nL45B+VIqEX8yuwk9JRyBWtxDgb2B90fba/H4yz3pddeuA7gM+BYqA9kf9T18thvHS4J23wLqgqMN\nx1Tk448FtsNxDo4VXTxSe921GRWxwN91bkcdmbdAUSfV71DAW0LxqJ448bHsLyGf/xxkVXwC+eHl\nOr5K9620vRzVSkSW2+8CNPJlHvAn1AFdS/nFN1Bn6udQucjn0ZyGno7172r5xfhn7gJOR9bVW2h0\n055d0GNUw/FB9PT0VeAQHJ/FebfQjBTQijN3c8iLryWABMlKaSoC3wI6wfLx14pjFI4bcLyG40hc\nXRqQQR2O0SrYzN2UkUMt2nepvyVjGOnGMRD4OvA14JfIx19WfScjjVg+/q7xKLJBpqIOU6N1aNbf\nZPKo7OFn0JDjh4Bv4HjFq6bWJvHfogV+wxD2myyHYzIaDjwQOAnHfZ4VZYHEf4vWuZssgW8BKSLw\nLSBlBL4FVMWxMY6rgTuAq4EdEw76QYLHNkowj98wjAKO/mgk2Kkomdp4XNVhyUaLciMao/5kbN2P\n0cSf/OugCvu24qgeozXJ9m/SkcNxMI6XcNyKYzPfkjJMKkb1/BxNwb4utq5aWoNaWEjW/6MZaWOh\nbwHecGyH8jGNAI7D2Yi1VqcWj/8BNIuzlJ50PgynkNSslV9TUqAhLa+0X4vhNJagwefriGMDHFcA\ndwM3A9t7DPqBp/Nmkp54/Jeisbz3oqyS8+qiyDCMZHH0RbNtz0D5psbjup0222hCcjV+bhxK4zsx\nWt4U5abvhWygvijneSk2RM4w0oIjh9KCXIhSapyC4zm/oowyJB43a23xl/rxL8fe/5xi/7+Ua+C9\nyR6LUGWm9mg5iP7asi3bcpLLjq15lmvozQZszvE47oq2b5wKfdlezr8fQ8oYR/Gong+hvO3DUCGM\nqyvsl/UO3MC3gBQR+BaQMoKGnMUxAscvcMzFcSKubM2GNBD4FpAiUjGq534KRbvXojG+k1HRjRBV\nkfpiIuoMw+geCvBfQoWDbkZ1bt/yK8pIC0n77+bxG0ajcewHXITSY5+M42nPioyukRqP3zCMtOMY\njwL+5mjm7Z1WAcsoh+XqSZbAt4AUEfgWkDKCuh3JMQzHT5Atey+wDY47mizoB74FZAlr8RtGs+Lo\njebQOOA2YGscc71qMpoC8/gNoxlx/A9KlzwfpUue6VmRUT/M4zcMI4ZjHEqSOAlVw7q1ySwdIwWY\nx58sgW8BKSLwLSBlBF36tGMIjh8B06PXBBy3tFDQD3wLyBLW4jeMNOPoheo8fxf4CzARxxy/ooxm\nxzx+w0grjg8hH38F8DUcj3lWZDQG8/gNI3M4xgAXADsB3wB+10KWjpECzONPlsC3gBQR+BaQMoIO\naxzr4Pge8BjwBLAljpszEvQD3wKyhLX4DcM3jjZgKvB9NAFrWxyv+RVltDLm8RuGTxy7Ih8/RD7+\ndM+KDP+Yx28YLYnjfcAPUYrzM4Df4ljrV5SRFczjT5bAt4AUEfgWkAocA3GczUs8BbyAyh7eYEHf\nfh+NxFr8htEIVPbwUNTKf4D7+QLXcbNnVUZGMY/fMJLGsSPy8fujvDr/8KzISDfm8RtG0+IYCfwA\n2BtVwrrGLB0jDZjHnyyBbwEpIvAtoGE4+uM4E9WpnoPG419VEvQDL9rSS+BbQJawFr9h1Av5+Aej\nWbf/AnbC8aJfUYbREfP4DaMeOD6AfPyhyMef5lmR0byYx28YqcaxIXAecABwNnAljjV+RRlGdczj\nT5bAt4AUEfgWUFcc/XCcBjwNLEE+/uVdCPpBYtqak8C3gCxhLX7D6Ary8T+OqmD9G9gVx3/8ijKM\nrmEev2HUimMi8BNgY+BkHPd4VmS0JubxG4Z3HOsB5wKfiv5ehmO1X1GG0X3M40+WwLeAFBH4FtBl\nHH1wnAQ8C6xGPv4v6hT0gzoco5UIfAvIEtbiN4xyOD4KXAS8AnwYxzN+BRlG/TCP3zDiOLZCAX8s\ncArw54xUwDLSg3n8htEQHMMABxyOKmFdgmOVV02GkRDm8SdL4FtAigh8CyiLozeOL6Ohmf2ACTh+\n0oCgHyR8/GYj8C0gS1iL38gujr3R8Mw3gb1xPOFZkWE0BPP4jezh2BxNwNoGOBW43Xx8I0WYx28Y\ndcOxLsqL/zngR8CncbzjV5RhNB7z+JMl8C0gRQTezuzohePzyMcfDmyD40eeg37g8dxpJPAtIEtY\ni99obRwBSpe8FNgfx+N+BRmGf8zjN1oTx6aoIMpk4DTgD+bjG01C4nGzFqvnRmAuKiOXZwjwZ+Bt\n4B/AhvWXZhjdwDEYx/eBR4EZwFY4fm9B3zAK1BL4fw7sW7Lu68B8lKXwQZS4yuhI4FtAiggSPbqj\nDcfRyMcfBUzCcR6OtxM9b/cJfAtIGYFvAVmiFo//AWBcybqPA0cCi4Dz0X+24+srzTBqxLEb8vFX\nAwfheMizIsNINbX6SOOAW4GJ0fK8aN3iaHkRsD7wbsl+5vEbyeF4P/BDYHfgDOC3ZukYLUBTjeO3\nAG80Bscg4BvAV4BfAMfhWO5XlGE0D7UG/tJW1GxgDDATGAGsil7luAaltgU9GcwA2qPlIPrbqssn\nka3vW205/777x8sRMIW9gKOA+7mWL/Eyc+G9oJ+k/nov59+nRY/v5fz7tOhp5HL+/RhSxjiKR/Wc\nC1yPJsNcAFxRYb+sP3YHvgWkiKBHezt2wvEgjkcjT7/ZCXwLSBmBbwEpIhVx835gbex1IoXhnCvR\ncM6NKuybii9gNDGOTXBch2M2jqNxNtvcaHkSj5s2gctIJ44BKIHaSeiJ8gc4lvoVZRgNoak6d42O\nBBT8vKwTUMu1cOSAQ1AStUeBHXG8nKQwTwTYbyNOgF2PhmGB30gPju3RePwhwNE4CwSGkQRm9Rj+\ncWwEfA/YHzgL+DWONX5FGYY3zOoxWhhHP+ThnwZcDYzHvTcp0DCMhLDAnywB5lvmCchfC/n4n0RV\nsJ4CdsHxvC9hngiw30acALseDcMCv9FYHJOQj78B8EUc/+dZkWFkDvP4jcbgWB9N/Dso+ns5jtV+\nRRlGKjGP32hyHH2BLwNnAr9F+fEX+BVlGEaSZH3mbuBbgDccORz743gOx1/YiiN9S0oZgW8BKSPw\nLSBFJB43rcVv1B/HBOAilHTqZBx/xv5jG0ZqMI/fqB+O4YADDkPj8i/BdajRYBhGdczjN5oAR2/g\ni2jy1R9Eaz1WAAAUXUlEQVSQjz/fryjDaEbCPq3QVjaPv9Vx7IPjaRx/xb1Xoa0cQaMkNQmBbwEp\nI/AtoPGEOQhHQ/hpCC+E8H4Il2Eev5FaHFsAFwJboSyaf7Kyh4ZRjXAIsCOwc+wF8FD0OhslJlyU\ntBLz+I2u4RgKfAdVwfohcDGOd/yKMoy0EfYGtqE4yI9GFfkeir1mQa60wWQev5ESHL2A44BzgDuA\nrXG86VeUYaSBMAeMojjIbw/8l0KA/wXwJORSMdjBWvzJEtAK+UccU1CahUXASTj+1Y2jBLTCtagf\nAXY94gQ0zfUIByPLZicKgb43xS35RyDXXcvGWvyGRxxjUSK1D6AMmn80H9/IFmFvYGuKW/ObUrBs\nbgJOBl4tY9mkFmvxGx1xDAa+haydi4CLcKz0K8owGkFYzrKZTXFr/omELRtr8RsNRIXMjwbOA+4B\nJuF43asmw0iMcB1gMsWBvi+FAP89ZNks9CaxSWmaR5+ECHwLqBnH7jgew/EAjh0TOEOQwDGbmcC3\ngJQRJHv4sBeEEyE8DsJfQfgEhMshfBDCn0J4GIRjo45a39g4fiNhHKNRYfNdgDOAG83HN5qfcBOK\nW/I7AHMotOavAGZCbpU3iR4xjz+rOAahQH8CcDFwAY4VfkUZRncIB9HRsulPx1E2zZIOPPG4aYE/\na8jHPxw4H7gPOB3Hf/2KMoxaCXuh2eLxID8OeJLiQP9SM42yKcECf5MTkKaxyY6dgZ8BvYCv4Xig\ngWcPSNO18E+AXY84AWWvR7gxxUF+MvAmxUF+JuRaafa4jeox6oBjE9TC3xNVwroex1q/ogyjlHAQ\nGj4ZD/SDgIdRgL9A73NveZPYIliLv5VxDAS+DnwNuAz4AY5lfkUZBkDYRkfLZnPgKQqB/iHghSa2\nbLqLtfiNbuDIAZ9Go3UeBibjeNmvKCPbhBvR0bKZRyHAXwXMaDHLJrVY4E+WgEb7uI4dUF6ddYAj\ncfy9oeevTIB52nECWvZ6hAPpaNkMptCSv1Dvc/FiPQFYltdGYYG/VXBsBHwf+AhKm3w1jjV+RRmt\nT9gGjKc4yG8JPI2C/B3At4HnM2jZpBbz+JsdR3/gJOTlXwWch2OJX1FG6xJuQHGQ3xFYQPEom39B\nznI7dR/z+I0KyMc/EI10eBL4II4X/IoyWotwAMWWzU7AUAqWzc/0NzfPm0QjlWT90S5I5KiObXFM\nw/Ekjv9J5Bz1J/AtIGUEvgUUE7ZBuCWER0F4KYSPQbgCwkchvATCIyEcH1k7SRAkdNxmxHL1GDEc\nGwDfBT4JOOBXOFZ71WQ0KeH6dLRsFlGwa25Als3b3iQaiWEefzPg6At8FeXWuQE4F4elijVqJOyP\niunEA/1w4BEKgf5hyFkpzXRgKRsyjXz8j6Hhb88Dp+L4t19RRroJ29BEqHhZwAnAcxR3wD4HOZu9\nnU5SHzfnA2ujV7mRJObxdxfH1jjuwfEsjv3qJ8kbgW8BKSOoz2HC9SDcH8JzIbwbwgUQvgLhzRCe\nAuFu0bj6tBP4FpAiUu/xvwsk1dmTTRwjgHPQzNvzgF/iSLLMm9E0hP3oaNmsR8GyuQRZNm94k2g0\nBT19nJgDbFxle+ofWVKDow/wJTTZ5XfA2TgsGVVmCXMo3XA8yG8N/Idiy+bfZtm0HKkfx98LjQRY\nhsb0XtBjRVlEVs5FqKjznjie8qzIaDjhCIp9+Z2A5RQC/O+AxyG33JtEo2XoaeD/MPACmrJ9G3rk\nbC/5zDXAK9H7RcCM2GeC6G+rLp9Ete+7JVOZxAnocf0UzmEpIetRwLf+ei7n36dFj8fl8EE49nMw\ncSvYaAIcOgbYAG57Ad76Nxx7OXAM5ManQ6/9PhJezr8fQxPyY5T+N4517pbDMRTHRTjm4Tg1Gq7Z\n6gS+BfghzEE4DsLPQngxhA+pyPedL0B4OYTHQLh1VFkqywS+BaSIxONmT3ykDZC//xwaPnYrcCwU\nZYM0jz+OozdwHJp8dTvwHRxzvWoy6kw4nI6WzdsU+/KPmWVjVCHVHv9QNJloHDAXefxpSQGcPhx7\nonTJC4D9cMzwrMjoMWFfYFuKO2A3Ah5DAf5XwHGQe92bRMMog03gSpYAx2uo03tb4DTgFlwmLbCA\nps4/H+aATSkO8pOAFyluzT8DuVrSYQc09fWoOwF2PfKkusVvVMMxhMf5ArAPmnl7GA5LVds0hMPo\naNmsohDgvwk8CjkrZWk0HdbirzeOXsDRKJna3cCZOOZ41WR0QtgXtd7jrfmRFCyb6JWb7U2ikSUs\nV09T4dgD9XWsAE7C8ahnRUYHwhwaNhcP8tsCL9HRsrHMp4YPLPA3BY4xqLD5zsDpwM2Rjx9gvmWe\nAC/XIhyKUg7HA/1qioP8o5Bb2mBhAfbbiBNg1yOPefypxrEOSpX8JdTSPxrHCr+iskzYB5hIcZB/\nH/A4CvDXAicAr1n9VyPLWIu/OzjagCNQcfN24Ixo9I7RMMIcMJriztft0CzxKL88DwFPmWVjNBlm\n9aQOxy5oPD7A13BM9yknO4TrUmzZ7IR+X6WWjRWaN5odC/ypwTEK+CHKT/RN4Dc4OsuKGGC+ZZ6A\nmq9F2JuOls37gX9RHOj/28SWTYD9NuIE2PXIYx6/dxwD0cSrE4FLgeNx2NjtuhHmkA8fD/IfAGZR\nCPAXI8vG6hIYRh2wFn8lVPbwUNTKfxD4Bo5X/YpqBcIhwGSKA30bxS35RyC32JtEw/CLWT1ecOyI\nfPz+aDz+PzwralLC3sA2FAf50ShVdazIN682sWVjGPXGrJ6G4tgY+AGwL/At4FocteRdqURAZnzL\nMAeMoqNl8xrwEFy4AE6dCjxplg2Qqd9GTQTY9WgYFvgBHP2BU6LXr4HxuLLF4433CAfT0bLpTaEl\n/11k2SyKdgjg6497EGoYRgnZtnrk4x+Msmf+CzgNx4t+RaWRsBeq9xoP8mOBmRR786+YZWMYPcY8\n/sRwbId8/GHAyTju9awoRYSbUBzkd0D1gPOToh4CnoDcKm8SDaN1scBfdxwbAOcBnwDOAq7soY9f\njYDU+5bhOhQsm3wa4n50HGWzsIcnCkj9tWgoAXY94gTY9chjnbt1Q3VtT0S5da5DPv6i6ju1GmEv\nYALFrfnNgCdQgP8DmrPwslk2htG6tH6LXz7+AagYynPAqTie86qpYYQj6WjZvEFxa36mWTaGkSrM\n6ukRjm2AnwCbIB//bm9aEicchAJ7PNAPoNiXfxhyC7xJNAyjFizwdwvHesA5wCFoWOFlOHyMHQ9I\nxLcMewFbUhzkNweepDgz5YspsmwCzMONE2DXI06AXY885vF3CY3HPwH5+DcDW+F4y6+oehBuTEfL\nZi6FlvyvkWXzjjeJhmE0Da3R4nf0BqYCDqUDOBPH04mfNxHCgXS0bAbR0bJpgRuaYRhlMKunKuq4\n/STwPWAeKojyYGLnqzthGx0tmy2ApynugH0hRZaNYRh5pk3LoZxeg4CBsb8Dy6yrbduUKdthVk8F\nHFOA84G+KNXC3VGd2zQRUORbhhtSHOR3BOZTCPDXADMgt7KhKhtDgHm4cQLsesQJSOJ6TJvWh54G\n4urbBgCrgBXA8pK/5dYtB5Yhq7bcthWoVGiiNF/gd+yASh5uBnwHFTbvrCCKB8KBcMJEuHR7CoF+\nCAXL5id6n5vnUaRh+GPatDbyreVTT92QCy+cQP2DdBudB+LSbW9U2Vb6922mTElqAmhiNI/V49gC\njdDZI/r7axwpGX8etgHjKW7Nb0lHy+Z5s2yMpmHatL4k31peSdday13d9i5TpjTb/zkb1YNjE5Ra\n4SDgIuAYHMv9igo3oKNls4BCgL8O+FeLWjZGGlBreQD1DcSl63IoeHYlAL9eZVvpfm8zZUoKn9aT\nJdR17YP+/Upf/Rsx/j29gd8xHDgdOA64EqVY8DD5KByA8srHA/1Q4BEU5H+GLJu5ZXYOMB83T0BW\nroU6/Kq3lq+/fgemTn217LbagnR/4G261hJehAJzbS3oKVMaOfclwNPvI5Qd1J/ygTip11r071fu\nlTjpC/yOQSinzinAH4FJOGY35uRhG5oIFQ/yE4BnUZD/Cxoy+h/IZa6l0jJMm9aL4tZyElbGWqq1\ndkeNWgfYuGTbAmpvLa9sxdZyJ63her/ywb4v8A6VA3G112LUJ9Cl/XKwuvplSJb0ePyOPqh1/x3g\nPuAsHP9JThpAuD7FWSl3Qq2ieFnAxyHXkLuwQby1nISnnH/fj0IQrZeXXLytsa3lxIhaw/UMsj1t\nDdf7tRJYmWtAsO0CGRjH72gDPoM6bF9Ak68SGM4U9qejZTOcgmWTnxj1Zv3P3UKotZyUp5z/u5pk\nO/xWNmGHX7413Jf6BdhaXn1QcKxHgK1HazgLtHDg1+Sr/VCN25XAN3FMq9Npc5S3bJ6jeJTNcwlb\nNgGN9C3VWu5Hsq3lvnQnAP/+95twyCEzq+wXby03xX/8nrSGr4ItjtEcjq62mtdQxyBbw+udBrWG\nA7LSB9Q5LTqqx7ErCvjrA2cCt3d/8lXYBw2l3A7YNvq7A7CUQoC/CVk2K3oqvUdMm9ab5FvLpZNJ\nagnS86psK/37TjdbywGXXtrejf1qosbWcL394W63hldqEs+LXdxvpbWGjXrQ2Ba/0iR/DwXns4Hr\nu1b9KhyKgns+wG8LbAX8F+XomRn9fRxyb3RJaeWp1/UM0r1JylMuDI9LRWAIId+B2ih/uNbWcD1f\nq1LmDRutQYtYPY4xKE1y3tq5DEeVMe5hjqGrxtJn7WR6h9uTC7clxzb0DocxePULDFv1AhuvfIVx\nS//L5IWvM+zdfNDuRyFIdDU4D0Q9+0l1+C0HVvnwlqPWcD/q2wHX2as39e2Aq8UbbroZlIZRhlQH\n/inAVcAI4GLg22U+Ex72za+vXDigV78VAwYvXTlw6MIV/QesXt6//9q3+/VjZd++uXd69+n1Ln37\nrwl7DQhzuf7k6EPvsBdtIazOrWVtbhU5VtIrXEavcBk53iHqiYey7+PBuytBOomp1wFlfMtQQbGR\n44b7A+/SWG+4tDVc9lpkmAC7HnEC7HrkSa3Hn0NB/0TgMeBvwP9Cx8yY0zdY+MS4t0c+PHbekrXr\nzp07rP+8cOMBi1dvNGjZyhFD3l4xdNiqxf1HrpqzeuNVc9v6r1oV9no3XLjOqhVzBqx55822MJyP\nOsDmIx96fsnrrZwCfbeJtYaHUOeAey68/yx5uaXbetH9gPpWN/ZZmYLW8HbYf+w4dj2KsevRQLob\n+LdD493viJZ/hVIqdAj8L5985bSXC358P+TB/5WCJ/9svoBIKLtlvTKv9YGJ5baFCmylN4jVdM0b\nrqU1XO61PDpn2e03w9SzNLO3dNu7GfSGh/oWkDLsehRj16OBdDfwbwK8Glt+GZhc4bPLgEtQoH+t\nWpKynGyXWdGrU6LW+hCKbxDrU2hR12JdJNka3isHzyR0bMMwjG5Rr+GcVfyo3Hl1Oke5k4ZoyvRi\nNDQubYzxLSBFjPEtIGWM8S0gZYzxLSBLdDfwvwaMji1vGq0rZSbZszRKOcq3gBRh16IYux7F2PUQ\naWzEApqx+DLwcWT7PAvs6lWRYRiGkTh7ouC/FFXEMgzDMAzDMAzDaD5GAneiDuTZaO4CaDTRn9HI\noH8AG8b2ORmNt5+NhrnmmQg8iYaB/hpZZqCcL9ehUU0zUdnGNNMGTEffG7J9LdYH/oT0vgRsQ3av\nx0noey0FbgXWIVvX4kZURP3J2LpGff9DgTnR8U6oy7fJOJuitNHrov/Uc4GtgXPRP8JQ4EfA5dHn\nN0NVjjYHdkH/GP2jbfehf5QRaLLKYdH6Y9DchuHAV4G7E/w+9eAr6LvfFy1n+Vr8Ec3VWBf9VjYk\nm9djfWAh+m7rALcBXyNb12JXlPo9Hvgb8f0HA2+iLMTjo+OOqucXM+Ae1L8xA5gUrRuObgigu/hF\nsc//Cdg/+sz82PqDgN9H729HHeWgOQjzUR6hNDISzdTehUKLP6vXYiOkr1/J+ixej0HoiWc8CkS3\nAQeTvWsxjuLAn/T3Xwc4ELglts/FwJeriWyrttHowGZoGOuDFE9iW4BSAvdFgTE+Ae3l6LMjKR7y\n+kq0nmhb/lhr0KPfyLqrrw8XoVrI8ToGWb0W49D3ux3ZG7ejll0Wr8dylGL9KWSLrkVPQ1m8FnEa\n8f0rHasiFvhrZwTwO+BY5NeVo1pipVaYz7AfStXxKJ0nkWr1awGaBzMJ+Cmqn7sY3RTL0erXYyPg\nAlQLYxhqkX6lwmdb/Vp0hvfvb4G/NoagDt5zgfujdbMpzDYcgQqgvBOtj09uG4vu4HMo9t3ik97i\nx+qF7tav11F/vdgN+AJqzT0QLT+CvseY6DNZuRYgzW8Cd6HUJDei+hBZvB7bAk9Er8Wo+NEuZO//\nSWngTvr7z65yLKMHDAD+jjpW4pwLXI88uQuAK6L149CPcQvU2fMGxZ02X0Z5hf4OHB6tPwb55iPQ\nqKF7Evge9WZnCh5/lq/FTOAjyNe+Hl2LLF6Psci/nog6um8DvkX2rkWpx9+I7z8YXfsPopE+1rlb\nB/ZCLdz463AKw7RWogC4UWyfU5Cf9zrq4MozCXmgK1Ba6/wjX2/040jrMLVyfJDCqJ4sX4ud0cz1\npagTLj6EMWvX4wTU0lyMbNFBZOta3E9xnDiRxn3/w9DTwgIqW2yGYRiGYRiGYRiGYRiGYRiGYRiG\nYRiGYRiGYRiGYRiGYRiNYBwaU72gk89Njz63f+KKDMMzlrLByAL/RDMnq/FBlPvcMFoeC/xGlvgg\nmkX5BzTT9lq/cgzDDxb4jazxIVQ4ZTSaGr+zXzmG0Xgs8BtZ4wGUM2UBSng1xqsaw/CABX4jayyP\nvX8XJb4yjExhgd8wDCNjWOA3skZpoYwsVHwyjCLsMdfIEtNR8eo83yzZ3lk5ScMwDKMJ6MoErjXA\nRxNXZBiGYRiGYRiGYRiGYRiGYRiGYRiGYRiGYRiGYRiGYRhGlvl/5QItqJhLYJUAAAAASUVORK5C\nYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# use n_jobs=20\n", "python_para_data=pd.DataFrame([[1000,.429],[10000,.531],[100000,1.92]],columns=['n','Python Parallel Time'])\n", "matlab_para_data=pd.DataFrame(mat_time_la_par,columns=['n','Matlab Parallel Time'])\n", "\n", "plot_data_3 = pd.concat([plot_data_2,python_para_data['Python Parallel Time'],matlab_para_data['Matlab Parallel Time']], axis=1)\n", "plot_data_3.plot(x=['n'],y=['Matlab Time (LA)','Python Time (LA)','Matlab Parallel Time','Python Parallel Time'])" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "**Discussion**\n", "\n", "Both Matlab and Python show dramatic improvements when bootstrap replicates are distributed across multiple processor cores. While Matlab is the fastest for this example, Python's parallel performance is impressive. In terms of percentage gains, Python shows the largest percentage improvements in run times when the linear algebra code is distributed over multiple processors. \n", "\n", "\n", "It is notable that Matlab's Parallel Toolbox is limited to 12 workers, whereas in Python there is no limit to the number of workers.\n", "\n", "\n", "The full table of results is shown below." ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
nMatlab TimePython TimeStata TimePython Time (LA)Matlab Time (LA)Python Parallel TimeMatlab Parallel Time
0 1000 24.366323 0.836 3.445 0.309 0.163328 0.429 0.479305
1 10000 26.902861 3.580 10.346 2.120 0.796194 0.531 0.211717
2 100000 63.204921 30.700 91.113 21.700 8.114325 1.920 1.406549
\n", "
" ], "text/plain": [ " n Matlab Time Python Time Stata Time Python Time (LA) \\\n", "0 1000 24.366323 0.836 3.445 0.309 \n", "1 10000 26.902861 3.580 10.346 2.120 \n", "2 100000 63.204921 30.700 91.113 21.700 \n", "\n", " Matlab Time (LA) Python Parallel Time Matlab Parallel Time \n", "0 0.163328 0.429 0.479305 \n", "1 0.796194 0.531 0.211717 \n", "2 8.114325 1.920 1.406549 " ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot_data_3.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Detailed info on machine this was run on:\n", "\n", "\n", "Architecture: x86_64\n", "CPU op-mode(s): 32-bit, 64-bit\n", "Byte Order: Little Endian\n", "CPU(s): 28\n", "On-line CPU(s) list: 0-27\n", "Thread(s) per core: 2\n", "Core(s) per socket: 14\n", "Socket(s): 1\n", "NUMA node(s): 1\n", "Vendor ID: GenuineIntel\n", "CPU family: 6\n", "Model: 63\n", "Model name: Intel(R) Xeon(R) CPU E5-2697 v3 @ 2.60GHz\n", "Stepping: 2\n", "CPU MHz: 1201.281\n", "CPU max MHz: 3600.0000\n", "CPU min MHz: 1200.0000\n", "BogoMIPS: 5187.74\n", "Virtualization: VT-x\n", "L1d cache: 32K\n", "L1i cache: 32K\n", "L2 cache: 256K\n", "L3 cache: 35840K\n", "NUMA node0 CPU(s): 0-27\n", "`" ] } ], "metadata": { "hide_input": false, "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 }