notebooks/02.05-Computation-on-arrays-broadcasting.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computation on Arrays: Broadcasting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We saw in [Computation on NumPy Arrays: Universal Functions](02.03-Computation-on-arrays-ufuncs.ipynb) how NumPy's universal functions can be used to *vectorize* operations and thereby remove slow Python loops.\n",
    "This chapter discusses *broadcasting*: a set of rules by which NumPy lets you apply binary operations (e.g., addition, subtraction, multiplication, etc.) between arrays of different sizes and shapes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introducing Broadcasting\n",
    "\n",
    "Recall that for arrays of the same size, binary operations are performed on an element-by-element basis:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([5, 6, 7])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = np.array([0, 1, 2])\n",
    "b = np.array([5, 5, 5])\n",
    "a + b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Broadcasting allows these types of binary operations to be performed on arrays of different sizes—for example, we can just as easily add a scalar (think of it as a zero-dimensional array) to an array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([5, 6, 7])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a + 5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can think of this as an operation that stretches or duplicates the value `5` into the array `[5, 5, 5]`, and adds the results.\n",
    "\n",
    "We can similarly extend this idea to arrays of higher dimension. Observe the result when we add a one-dimensional array to a two-dimensional array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 1., 1.],\n",
       "       [1., 1., 1.],\n",
       "       [1., 1., 1.]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M = np.ones((3, 3))\n",
    "M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 2., 3.],\n",
       "       [1., 2., 3.],\n",
       "       [1., 2., 3.]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M + a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here the one-dimensional array `a` is stretched, or broadcasted, across the second dimension in order to match the shape of `M`.\n",
    "\n",
    "While these examples are relatively easy to understand, more complicated cases can involve broadcasting of both arrays. Consider the following example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 1 2]\n",
      "[[0]\n",
      " [1]\n",
      " [2]]\n"
     ]
    }
   ],
   "source": [
    "a = np.arange(3)\n",
    "b = np.arange(3)[:, np.newaxis]\n",
    "\n",
    "print(a)\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 1, 2],\n",
       "       [1, 2, 3],\n",
       "       [2, 3, 4]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a + b"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just as before we stretched or broadcasted one value to match the shape of the other, here we've stretched *both* `a` and `b` to match a common shape, and the result is a two-dimensional array!\n",
    "The geometry of these examples is visualized in the following figure. (Code to produce this plot can be found in the online [appendix](https://github.com/jakevdp/PythonDataScienceHandbook/blob/master/notebooks/06.00-Figure-Code.ipynb#Broadcasting), and is adapted from a source published in the [astroML](http://astroml.org) documentation. Used by permission.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Broadcasting Visual](images/02.05-broadcasting.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The light boxes represent the broadcasted values. This way of thinking about broadcasting may raise questions about its efficiency in terms of memory use, but worry not: NumPy broadcasting does not actually copy the broadcasted values in memory. Still, this can be a useful mental model as we think about broadcasting."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rules of Broadcasting\n",
    "\n",
    "Broadcasting in NumPy follows a strict set of rules to determine the interaction between the two arrays:\n",
    "\n",
    "- Rule 1: If the two arrays differ in their number of dimensions, the shape of the one with fewer dimensions is *padded* with ones on its leading (left) side.\n",
    "- Rule 2: If the shape of the two arrays does not match in any dimension, the array with shape equal to 1 in that dimension is stretched to match the other shape.\n",
    "- Rule 3: If in any dimension the sizes disagree and neither is equal to 1, an error is raised.\n",
    "\n",
    "To make these rules clear, let's consider a few examples in detail."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Broadcasting Example 1\n",
    "\n",
    "Suppose we want to add a two-dimensional array to a one-dimensional array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "M = np.ones((2, 3))\n",
    "a = np.arange(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's consider an operation on these two arrays, which have the following shapes:\n",
    "\n",
    "- `M.shape` is `(2, 3)`\n",
    "- `a.shape` is `(3,)`\n",
    "\n",
    "We see by rule 1 that the array `a` has fewer dimensions, so we pad it on the left with ones:\n",
    "\n",
    "- `M.shape` remains `(2, 3)`\n",
    "- `a.shape` becomes `(1, 3)`\n",
    "\n",
    "By rule 2, we now see that the first dimension disagrees, so we stretch this dimension to match:\n",
    "\n",
    "- `M.shape` remains `(2, 3)`\n",
    "- `a.shape` becomes `(2, 3)`\n",
    "\n",
    "The shapes now match, and we see that the final shape will be `(2, 3)`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 2., 3.],\n",
       "       [1., 2., 3.]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M + a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Broadcasting Example 2\n",
    "\n",
    "Now let's take a look at an example where both arrays need to be broadcast:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "a = np.arange(3).reshape((3, 1))\n",
    "b = np.arange(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, we'll start by determining the shapes of the arrays:\n",
    "\n",
    "- `a.shape` is `(3, 1)`\n",
    "- `b.shape` is `(3,)`\n",
    "\n",
    "Rule 1 says we must pad the shape of `b` with ones:\n",
    "\n",
    "- `a.shape` remains `(3, 1)`\n",
    "- `b.shape` becomes `(1, 3)`\n",
    "\n",
    "And rule 2 tells us that we must upgrade each of these ``1``s to match the corresponding size of the other array:\n",
    "\n",
    "- `a.shape` becomes `(3, 3)`\n",
    "- `b.shape` becomes `(3, 3)`\n",
    "\n",
    "Because the results match, these shapes are compatible. We can see this here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 1, 2],\n",
       "       [1, 2, 3],\n",
       "       [2, 3, 4]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a + b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Broadcasting Example 3\n",
    "\n",
    "Next, let's take a look at an example in which the two arrays are not compatible:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "M = np.ones((3, 2))\n",
    "a = np.arange(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is just a slightly different situation than in the first example: the matrix `M` is transposed.\n",
    "How does this affect the calculation? The shapes of the arrays are as follows:\n",
    "\n",
    "- `M.shape` is `(3, 2)`\n",
    "- `a.shape` is `(3,)`\n",
    "\n",
    "Again, rule 1 tells us that we must pad the shape of `a` with ones:\n",
    "\n",
    "- `M.shape` remains `(3, 2)`\n",
    "- `a.shape` becomes `(1, 3)`\n",
    "\n",
    "By rule 2, the first dimension of `a` is then stretched to match that of `M`:\n",
    "\n",
    "- `M.shape` remains `(3, 2)`\n",
    "- `a.shape` becomes `(3, 3)`\n",
    "\n",
    "Now we hit rule 3—the final shapes do not match, so these two arrays are incompatible, as we can observe by attempting this operation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "operands could not be broadcast together with shapes (3,2) (3,) ",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-13-8cac1d547906>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mM\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m: operands could not be broadcast together with shapes (3,2) (3,) "
     ]
    }
   ],
   "source": [
    "M + a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note the potential confusion here: you could imagine making `a` and `M` compatible by, say, padding `a`'s shape with ones on the right rather than the left.\n",
    "But this is not how the broadcasting rules work!\n",
    "That sort of flexibility might be useful in some cases, but it would lead to potential areas of ambiguity.\n",
    "If right-side padding is what you'd like, you can do this explicitly by reshaping the array (we'll use the `np.newaxis` keyword introduced in [The Basics of NumPy Arrays](02.02-The-Basics-Of-NumPy-Arrays.ipynb) for this):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3, 1)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a[:, np.newaxis].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 1.],\n",
       "       [2., 2.],\n",
       "       [3., 3.]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M + a[:, np.newaxis]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also notice that while we've been focusing on the `+` operator here, these broadcasting rules apply to *any* binary ufunc.\n",
    "For example, here is the `logaddexp(a, b)` function, which computes `log(exp(a) + exp(b))` with more precision than the naive approach:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.31326169, 1.31326169],\n",
       "       [1.69314718, 1.69314718],\n",
       "       [2.31326169, 2.31326169]])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.logaddexp(M, a[:, np.newaxis])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For more information on the many available universal functions, refer to [Computation on NumPy Arrays: Universal Functions](02.03-Computation-on-arrays-ufuncs.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Broadcasting in Practice"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Broadcasting operations form the core of many examples you'll see throughout this book.\n",
    "We'll now take a look at some instances of where they can be useful."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Centering an Array"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In [Computation on NumPy Arrays: Universal Functions](02.03-Computation-on-arrays-ufuncs.ipynb), we saw that ufuncs allow a NumPy user to remove the need to explicitly write slow Python loops. Broadcasting extends this ability.\n",
    "One commonly seen example in data science is subtracting the row-wise mean from an array of data.\n",
    "Imagine we have an array of 10 observations, each of which consists of 3 values.\n",
    "Using the standard convention (see [Data Representation in Scikit-Learn](05.02-Introducing-Scikit-Learn.ipynb#Data-Representation-in-Scikit-Learn)), we'll store this in a $10 \\times 3$ array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(seed=1701)\n",
    "X = rng.random((10, 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compute the mean of each column using the `mean` aggregate across the first dimension:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.38503638, 0.36991443, 0.63896043])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Xmean = X.mean(0)\n",
    "Xmean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now we can center the `X` array by subtracting the mean (this is a broadcasting operation):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "X_centered = X - Xmean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To double-check that we've done this correctly, we can check that the centered array has a mean near zero:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 4.99600361e-17, -4.44089210e-17,  0.00000000e+00])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_centered.mean(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To within machine precision, the mean is now zero."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting a Two-Dimensional Function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One place that broadcasting often comes in handy is in displaying images based on two-dimensional functions.\n",
    "If we want to define a function $z = f(x, y)$, broadcasting can be used to compute the function across the grid:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# x and y have 50 steps from 0 to 5\n",
    "x = np.linspace(0, 5, 50)\n",
    "y = np.linspace(0, 5, 50)[:, np.newaxis]\n",
    "\n",
    "z = np.sin(x) ** 10 + np.cos(10 + y * x) * np.cos(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll use Matplotlib to plot this two-dimensional array, shown in the following figure (these tools will be discussed in full in [Density and Contour Plots](04.04-Density-and-Contour-Plots.ipynb)):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATYAAAD8CAYAAAD9uIjPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA6HUlEQVR4nO2deZBdd3Xnv+fdt/VbepOszZIl2yg2mMUOisngTBYCxCEUpipMYqZCzBSUp6ZgQiaTBQ9VmQkJU8xMVchMDZPBBQ4QCEsgVDyJg+NgAyFgsE1kG29IyAJJ1t7qVq9vue/MH+9Jeuecn/rd7n7qVt93PlWvun/33eV3f/f2r+/9no2YGY7jOGkis9YdcBzH6Tc+sTmOkzp8YnMcJ3X4xOY4Turwic1xnNThE5vjOKkjm2QlIjoIYBpADKDJzHsuZaccx3FWQqKJrcPPMfOpS9YTx3GcPuGvoo7jpA5KEnlARM8DOAOAAXyEme8OrHMngDsBIKLcK8vZsfPfxZW82WejLNtD5Zpob8zNmG2yFIv2fMvud0LtuL6QE+1M3WzSPqvuZqTahZbZZLQ4L9rj2VnRLpId16Ya66nWkGw3ZBsAak35UB26XFEk+zeUbYh2ObInHUFuM6fGcqpm+0Kz8v9gblbug+bkNQQAbql1CgWzTrOSVW15ksMlOdYAMBbNiXaG5HGmY9v/ybpc1lyQx83IYQvSUu84mWJs1hkvyL7pvgJAjmS7pq7rZFwy25xtFC+sf3wKjal5MistgV/4uTKfnrD9D/HYE7X7mfnWlRxvNUn6KvpTzHyEiDYBeICInmXmr3ev0Jns7gaAkfxmfvWWt57/bupV280Oj94ir8kr9vxAtH9967fMNluiKdHeu3CVWecvX3ilaD//3FbRHjqsZi0A2QXZXhiXd1nrWvuH9UvXfU+03zb+TdG+Lmcnw+NxU7Tvm7lBtP/uhGwDwIGTG0S72bT9H67I/r30iqOi/RPDPzTbjERyIn58Vo7l3/7A9iX3SFW0t35T/sFGe/eZbVrzcnCjXbvMOqdevUm0T94ix+n1r5BjDQC/PP6oaJczclL92syLzTZf+tErRPvMc+OiPXTMvsCo+RLzV8h7Y+j6SbPN7dd8V7TfMvJds862SF7H/Q359/DXZ28y23zl2HXnf9/7rk+a75fKqYkY377f/m2GyG39wcYVH3AVSfQqysxHOj9PAPgSgJsvZaccx1kNGDG3En3WGz0nNiIqE1H13O8AXg/A/gt1HGddwQBa4ESfXhDRPUR0goiCcwO1+V9EtJ+IniCiH+/67g4i2tf53NGPc0vyKroZwJeI6Nz6f8HMX17KQQJyk3nEryvxYqEltTEAgHoLKwZEEa0vcV4eqJWzr3KtplqglAtmK2XMx7J/Cyz7H8PqTfq/iNaF4pb9P9NsqP3O2Ut2VvXvSGlUtK8tWWP25tykaF9XOibaP9h8hdnmiauk7nP2haJojx+Xr5TtHR0UTT5yzKwy+v2KaM9vlDrp10dfZLbZUjgr2reUvy/a1xXl6zgA/MTmMdH+Wk1ew7mMfNUGgOyMHFvOyZu50bDX40RD7mciLpp1Nivds5SRWlcpsvdPo+v+YH2TLpMW+vY09nEA/xvAxd6RfxHA7s7nVQD+FMCriGgcwH8GsAftufYxIrqXmc+spDM9JzZmPgDgFb3WcxxnfcFgNPr0msnMXyeiXYuschuAT3LbWvkwEY0S0VYAPwvgAWaeAAAiegDArQA+s5L+LMWPzXGcFMEA4gSvmR02ElG3xebukHfEIlwJ4FBX+3Bn2cWWrwif2BxngEmin3U4tZ4ijtxB13EGFAYQMyf69IEjAHZ0tbd3ll1s+YpYnSe2kPEgluJnI5ai/gIHjAeKHFnnwkpOiq5UkOu0CvaUW8qHyBoP7LGbLdnfBuv9WvE3UvsN9d/2Tf7voXn7vyhWhozjQ1K8PlyVojkAXF04IdpbspOifdPoIWie3yH9vs4el/stH5PfA0DhtNxvfMZqwrkD0qAwPr5TtI+NSeMCAPx95XrRruyQ/nK7C8fNNjeU5d/LxBZpDHkiZCQ6rRx91SqFjNWoZpvSCflsyxoPWizvj6qysF2RnTbblHIXDGOZkEVuGayiI8e9AN5NRJ9F23gwxcxHieh+AP+ViM7dTK8HcNdKD+avoo4zoDB4KRrbohDRZ9A2BGwkosNoWzpzAMDM/xfAfQDeAGA/gDkA/6bz3QQR/SGARzq7ev85Q8JK8InNcQYUZqDRp1pOzPzWHt8zgHdd5Lt7ANzTn5608YnNcQYWQtwnf7jLjVWZ2CggUmmJoBbLrtSNZgVzEXKkPWttwHc2L9eJ8zYIO6PkPE5gUmmoleoqcr4VOGe927zqPwV0E1ZaZDagsUGFss7mpS50cMRqXz9WHhHtLVkZh/vSocNmmx9tkft5cJd0pJ08YZMSbDq5WbRpzgaEx6dOi3Z5n9QIR8ess/Dxqgxd/EpBam7RJjuWG7PSqfflVam5NQMO0gfyMlZ3oS5vluqQ1VJzytm2oT3LYd0sqhl5v+/IyTEBgK2lC/3PZ5IFry8GA2iltPqmP7E5zgDjT2yO46SKtoOuT2yO46QIhpVU0sIa+rHJdkPpG7VAEHysVKp8wA+snJWaR6EgdazZQNLIVk3uN8m1Nn5sSkcJmdG10lIkGbCfDfhE6d2EEmVGdRWonZWX9ciw1NMA4PlhqVtpv68rs9bf7FUjB0T7wE6pPx0+tc1sUzoljz08ZQPlm4deEG1+QfZldJ9Nulgflr5tz5Vk3r1S1g7Uq8YOivZIVup9NwzbwPliJO8fnayyGNlEDFWd4C8BQyT1yR1KDwSA6ysX/P2+liQrZg8YZP6m0oI/sTnOANMKOCWnAZ/YHGdAcY3NcZwUQohdY3McJ020M+j6xLZsdLZcACb6tq6qMemMtADQUE6wmUAI75ASc4fysj2Tt6K+tlNwVq5Dgaf1pvpPp4PgQ6ccqR3pDMD5jHU41mR0wD6ArPF5lX2br9iKTfuU0+sNZRnQvi1gPLihIJ12f3qT3Obz1w6bbc6ckc62hTPW2TY/K0+gNSmdhbM/lAH7ADA2LMX2hqqE9njRFinRxpkbh+X5bMzaymiFirwmU6r6la72BQBb87L/ZQqVRpNEJK/Z5shOOC8f+tH534eC5daWBjMZx/K04E9sjjPAtFxjcxwnTbSNB/4q6jhOqnDjwcoIOOhmVHB3M5YDPBfbYHV9EaKAeFdS2kN3cj4AQME69bbyizvoJknql0SryKnHfu2gq51BAQCR0vsC4l2kguC172ZzyF7mH41KfeyZYelcuztvq0ntUIHy/6K8X7Sf3y4ddgHgH89cJ9pnztjresXZLaKdqcsTaE1ZZ9Xi8zJl11hFandx0eqKe7MylX5mhxxbHRQPABtzMuHjtrzUHkPJQkdVMerRjC243Ysi2Wu2K3fhnHUCheXgxgPHcVJJ7A66juOkCQYFUtqng3Q+hzqO05NzxoMkn14Q0a1E9Fyn0vt7A99/iIj2dj7fJ6LJru/iru/u7ce5pXO6dhynJwzqy6soEUUAPgzgdWjXBX2kU8396fPHYv4PXev/ewA3de1inplvXHFHulizDLrar1FXqaq1Ag66Kj9GFLBKaKfXUk4aE6Kg8UAei5WxIInxoKUsDqGSZYWMXEcbD7RzMQBQ1PvYUV31VznstnL25p0ZluL6U+NSwH9xWbYB4EqVceKanBTwXzv2NDRHr5XZPZ4/a2vh5qelE+/YvDw2HbaGDExMimblgKwE1ShZZ+GJosz4+2ROZgQJOUj/+LDMFrMlL4+7IbJOvVVlwBoN7DdST0E1lte+wfY+rXYZDEL3/nLok/HgZgD7mfkAAHQqUd0GwN4Qbd6KdrGXS4a/ijrOgMLc9jRI8kGnEnzX586uXSWu5k5EOwFcDeDBrsXFzj4fJqI39+Pc/FXUcQaUtvEgcUhVvyrB3w7gC8zikXQnMx8homsAPEhETzLzD1ZyEH9ic5wBpk/Gg6VUc78dwGe6FzDzkc7PAwC+Cqm/LYvVeWILOJVqv8ZYO+i2bMWjBbWsnLEVggpaY1OZVPN5q3fM53UUvGxmot71svXFD9UQipSDbkn1v5qzmVezObmnkKO4lnBys7q/dqNmSf6nPjQmK1DtHb7KbKOddl+inFVfWbTV449vlhrbX+wum3UmZ+Sxc/Ojol1pBkbztDx25pR0Hq4etPdPc0jqipN5qe19T2luADCak86123LyuFpPA4AtkexvgexTkb4XpltyP3MBjXau6+mqHzGeDOpXoslHAOwmoqvRntBuB/Cv9UpEdD2AMQDf6lo2BmCOmWtEtBHALQD++0o75K+ijjPA9CNWlJmbRPRuAPejnQH/HmZ+iojeD+BRZj7nwnE7gM92iief48UAPkJELbT/A3+w25q6XHxic5wBpV1XtD9qFDPfB+A+tez3Vfu/BLb7JoCX9aUTXfjE5jgDi1eCXxEhNzAdzN1sqsrqIT+2BBYc7RumNbZi3vqKLRSkHqOrr0cBjS2Jb5vZRj32l1Ug80jWBkvnlCbYzPU+bnZBrhM1rEYVF+RY1oelH9jeEWutf1Fpl2hrv7bNgXH6l6Xvi/ax7bZi1r2zLxftiXlZlSpTt5XsTd2qaRl4njsxrddAtSK11EZJts+UbN+eHJLJAXYUpe/etbmTZpuCujcqZAP/ayyv63RLto/HNoj/RHxBE6zxhPl+qbTL73miScdxUgQz9e1V9HIj8cTWCZt4FMARZn7jpeuS4zirhedjA94D4BkANlbFcZx1RzsfWzo1tkTTNRFtB/BLAD56abvjOM7qQUsJqVpXJH1i+xMAvwugerEVOrFjdwJAMVKrBZwNtfGg1VKVlWLlNAtggeWykEVHZzStKOOBrloFADNKoG+pgPxsyHiwjCBkU6VK9XVEp8IFUCrI/k5av1O0lP5LTdm3rHHYBYYi2ZdGWe7k1IgV0r9VvUa0ry7I6lEjRets/iLldPwLI0+adU5eXRHtf1z4MdGO6vakqSUNCkM/kv2nBes4WzgpswNUK/K4zXIg03BVHufJijSqXFc4arbZmZUGhSgTmBjU7TOnRPyDjY1mkwO1Ted/n2/90O5zibTdPQb0iY2I3gjgBDM/tth6zHw3M+9h5j35jLXoOI5zeXEuVjTJZ72R5IntFgBvIqI3ACgCGCaiTzHzr13arjmOc6lJa82DnmfFzHcx83Zm3oV2SMSDPqk5zvqnnbaIEn3WG5eNg26rKQcvqLGpku0hUTOnnF51xWxTtQpAPi+1rqbSqPJZ6+CazYTC3C8QCpvXDroFdb+MRKakOyoFGSh/uhioZJ9XO1LNaMEG/hdPy/3Eah+Nih3/p6syAeQ3ylIL26KqWAHAy3LynF6mAucB4A3jUnebvFa63z5ZtwH51JT9Y5Ka4NBR6bDb3kZelcKEuleO2/upNiYdl58d2yzaz1Vt4Pz1eamxjWV6J1FYUK97x5pW43x29sKxFwIO7MshrRrbkkaHmb+KdloRx3HWOe3sHul8FfXIA8cZUNohVT6xOY6TKvyJbWW0An5sSqJiFQS/0Oztxxay6ERK0NMFUspZ69+kfdtqGal35LNWoyr00NhCZEwleNkeDWhsYwW57GDJHrc5pMZFFW/R2hIAZOflOZUycpvmkPUdmyxL7eufKleL9pa8rdg+Ovy4aO/MWk3nJ1SCyrObpLvQbMP25fmm1raU5patQFOYkNc+EyvN7awdp+Jp5d93SvpoPr1BBskDwCuKPxLtzZEdlxwWd6GYiYtm2ZHZ0fO/h5JELIe0Rh74E5vjDCjnrKJpxCc2xxlg/FXUcZxU0ceaB5cd6ZyuHcfpCQNocibRpxdEdCsRPUdE+4novYHv305EJ4lob+fzzq7v7iCifZ3PHf04t8vGQRcqa+1CbLs2F8tMpPUkGXV7VK0CgCHltEuqw8WA8SCnjAdR0CW3R99U9aJqxgbBb1DGg6hsHYybZSmcN4fkjchRwMgyIx1/86dl/yv5QGUrZaQ4Ux4T7YdK0mEXADbmZCbbssqoCwDjKkj8J4cOiPbMFiuk/1VT3h8v8BWizZG9fypFeY/pal7mngSQ1X6+U3IMnj+7wWxzcET25fqAU/K4MtbkEjw4zTUuHLtfT1r9eBXt5Gr8MIDXoV0s+REiujdQlOVzzPxute042lXh96A91z7W2dYO2hLwJzbHGVS4/Sqa5NODmwHsZ+YDzFwH8FkAtyXsxS8AeICZJzqT2QMAbl32OXXwic1xBpRziSaTfABsJKJHuz53du3qSgDdPjuHO8s0v0xETxDRF4joXIHlpNsuCTceOM4As4RX2lPMvGcFh/p/AD7TKYz8bwF8AsBrVrC/RVmdiS2ksWk/U62xBRx0ayy722Dbfa116crwIQfdSl7qTboC1VDW6lo6uF4nuAwRkXxAzih9YzhjK8GP56XIUyrVzDrzFenQWi/LsYyH7DhlJ+U5Zs5KLc+qWsBwTjqnNtV+9xdlgDgAfDl3g2iXN9v+v7ookyZuiWTfXl3aZ7aZ3ip7+GX1B3qcbKLGVl5qmsWT6noELqHWh6N5uc3EnM09eKQutceJor2XR1RgfDlB0tFc1za0jESnmj4mmjwCYEdXe3tn2YVjMZ/uan4UF6q9HwHws2rbr660Q/4q6jgDCoPQbGUSfXrwCIDdRHQ1EeXRTm92b/cKRNQdKvImtOunAO3q8a8nojEiGgPw+s6yFeGvoo4zwPQjpIqZm0T0brQnpAjAPcz8FBG9H8CjzHwvgN8gojcBaAKYAPD2zrYTRPSHaE+OAPB+5pUXTfWJzXEGFe6f2wgz3wfgPrXs97t+vwvAXRfZ9h4A9/SlIx1WyY+tdzEXUokma7H1UZuLZTC0TjwJWL81rX2Vs1bjCfm29fq+kJG+bZE6oSTv+Dnjx2aPo/3AxktWe/nhsNS+6iNKi6zascxPyLGjGbnfzIStpD6k/M1G8jLQvJW3weqPF7aLdimgV+Y2yrG8sfCCaG8O6E0/VXlOtGtb5Dk/GPiDPZ6VhVm0RpibMZvAxJqr3TYa9k/odF2Oy2TL6nAtlgerKr+2bTnrxrVx6MI22QTJK3uR5mIu/sTmOAOMT2yO46QKBiHubRhYl/jE5jgDjOdjcxwnVXAfjQeXG2vooMuqLQe43rSC97w2HrA1HuRYCtHaeKBFfwCo5hZ30K2GDA6RdtC1+9XErA0Mi1eGB4DxSDrodgvI5zg8LB1CtfGgNmLHslCVCQXyU3IdnrHZfKMTUtCuqLeYVlYaMdrLpHD+7WiXWUcnFGiNyR2/JFBtfUskx+EnK/tFu7HFnvM31HU9lpfj1jhj/xyooZydy/Ia5gMifk1ZHGbZGlVi9UdRUoakHbnT0OyuXqh+9XDU+35LAvvE5jhOukhvPjaf2BxngPEnNsdxUgUzELd8Yls2yRx0ZbvetF2bVRpbLeCgW85IPUwHxZcCTrBlpZdllP4RdOrVQfBYetUqTSjZoK5ctaVoHWerFenAOjWmtMgNAY3trBy77FlZgYrmrFMsT0tdK2rJsR3OWNcBjsqifSYqmXX+ia4RbV1gZGHMXufd+WOivUFpbjeWZaUoAGhskuPwuEog+kLZVl9vLMj7MJOT51wdsvdGSMftxRDJa7YtULHs5aUL2X10Eobl4lZRx3FSBcNfRR3HSR1uPHAcJ4UEVKJUsEqV4O0i7bKlg+CbQT82qbXMtax/UCWQfLKbUELIIa2xKQGwEoU0NrlMB9+Hysy0eiQHDAW36AIvG/LWj21TRS47OyZ9xxbO2rSRubORakvtqzgd0NjmZSLM1pSscB4F/kqGaZPcB5XNOpOQy77JUnNbiK3G9uox2d9r8idEezTSVViAl5RkcL3WX0fyNtHnqXnb3242lazmqa9REv1VJyEdz9h7u1tXLJJNJrAc/FXUcZxU0baKeqyo4zgpw19FHcdJHWl9FU3nc6jjOD1hEJiTfXqRoBL8bxHR053ye18hop1d38VdFeLv1dsuh1UKgu/toKsrBDUbgQy6zd4ZdBuqOrwWiHMB50kb0C47k8R4kFcCcUS9b4ZWgurx2iixMWuNB5tLUsQ/OaKcYjfacZqfVdlj55TD7tyw2SZbk+PUOiUDtePTNuurNiiMsq1kRSz7O9mShoFHmzuhObtdGkT2jMsA/KsKNoi8SLL/u4pynaHICvITJdk37Ty8IWeNFFfl5X5D1cc0OkFC6P7Z0HXPZRPcO0nox5towkrw/wxgDzPPEdG/Q7tK1a92vptn5hv70JXz9HxiI6IiEX2HiB4noqeI6A/62QHHcdYIBrhFiT496FkJnpkfYuZz4RQPo11m75KR5FW0BuA1zPwKADcCuJWIfvJSdspxnNVhCa+i/agEf453APi7rnaxs8+HiejN/Tivnq+izMwAzr3/5DqflNpSHGewWIJVdKWV4AEARPRrAPYA+JmuxTuZ+QgRXQPgQSJ6kpl/sJLjJNLYOu/QjwF4EYAPM/O3A+vcCeBOAChGUu/Q1bSBQBC80thasX38nVfV4WfigllHB8brSvAhB92CisAvqCR+1chqJGUVhFxMkGhSa2qxuqtCqokN4rd63+aCdBI9VZVVkubGrbPnwoLUjqKafHjP1mxlpcrCqGhnarIv8cSk2aZ1WpaIzLTszTCidDdqyf5PNayD8XO1bXKdmlznpeMySB4Adg5J7auirquuCAYA41mpoen7ZyTgCHylqjA1Hrh/IhV8XlMJUnUbAOKubbgf9UDRN6toz0rwAEBErwXwPgA/w8znbx5mPtL5eYCIvgrgJgArmtgSWUWZOe6Ie9sB3ExELw2sczcz72HmPfmM/aNwHOcygwEwJfssTpJK8DcB+AiANzHzia7lY0RU6Py+EcAtALqNDstiSVZRZp4koocA3Argeys9uOM4a0s/HHQTVoL/HwAqAP6S2hbfHzHzmwC8GMBHiKiF9oPWB5U1dVn0nNiI6AoAjc6kNoS2Sfe/rfTAjuOsNYksnolIUAn+tRfZ7psAXtaXTnSR5IltK4BPdHS2DIDPM/Pf9LsjjuOsASk1Ayaxij6Btpi3fALPuxldpaop5b5WwEF3viENAzrbB2CddrXYGwWupHaC1dk9QoK9ztSbI12ByqKNBVogbgRusljtSfcVAMayMtvqlaUp0Z4PZKA9pLKnzDWkLpoJjH+mLkX9ckONbTMgeKsMILEyJgBApJxTh+Mt8rhNedx2/+Q5Ha9tFO3peWtwOL5BGrWuqZwS7a15OW4AMKaMBzpT7xVZeX4AsEUZFMYDN0NGOeDOsLyupwPGs+Pxhf432I7jkuH0hlR5rKjjDDKD+sTmOE6a8Sc2x3HSRn9CTi871q4SvA6C1/JMw/4nqanKVbNN66Crs+pqTSrkoKuX6aruWk8DAhl0jZZnhRXtoNtQA7PAVtfSQf0htAPx5oLUfZpV2xedYPCIso7Nxlajojir2jJQvtSyfyWRcsiNp60TrHbs1WdcbW2CJtOQellUU/dGIIj/6Xl5b0xulLri7hF7zj9WPi7am5SmNhoIcNeaWjWQDbfB8n6ZVuN0qDlqttlXv6A9zvNR8/2SOefHlkL8ic1xBhhPNOk4Tvrwic1xnNThr6LLJ1EleFMZ3upCNe171QwkUFS+bQsZ2c5EVgfSGpvW5UJ+bLpKkK7irv2UAKCh/LUW1LjMBipszbakjhjS3LTfnQ7cjobsOetq9xmVqeAwjZltZqA1KDX+NGq2KalxiI7a/reU7taalP5kmcD9U2kqv8G61NSihcBYLsj+H5rfINraTzLEiKrQviNr/clyaiwLZPerfRr1tT/UkH0DgKdmL2QCmm893rOvSQglqEgD/sTmOIMKE9CnkKrLDZ/YHGeQ8Sc2x3FSh09sjuOkDp/YVkAgayqZIHi1QtO++9frsru6ahUAzCsH3RorR1q2QeSRMQQoYwLZbbRDbjFBVapY3UWzykl2umUdRCdjWbFJGxMAa1DQ/Q9VtiqU5IDrCk25yDoyHyQpaM+ohKIcWZGco1HRLuXtLRcdk9esdVYaE3hWCvYAQCekaF9SYnwmtg66OnA+05DHPRWPmG2eVOp6NScdcrepbLkAsC06IdojCdK5NliudKpZNescnBk//3s97sOfrjvoOo6TRtwq6jhO+vCJzXGctOFPbCsgSZUqHQQfyKeIpkmOaDW2GaW7VbNSVymp6lIhdBB8KLmj1thyKv2LrkIEtAu0irbSxibjMjRTatlU3LtQjj5HXY0JAKrRvGiPqGSV5ax1Sh7KynHYl7tCtKdztv+tnLzFmkWrHZWHVHLQY1Jr5LM2cB4NeY50RganBzVPVf2KWN8/ViM8GUnd7cmCrI51ZcFqbLtyUv8bD+i6vdAJUwFZiUtXpF82fdoPEd0K4H+incPgo8z8QfV9AcAnAbwSwGkAv8rMBzvf3YV2rdEYwG8w8/0r7U+iKlWO46QQXsJnETplAz4M4BcBvATAW4noJWq1dwA4w8wvAvAhdOqmdNa7HcANaBeJ+j+d/a0In9gcZ5Dpw8QG4GYA+5n5ADPXAXwWwG1qndsAfKLz+xcA/Dy1y1XdBuCzzFxj5ucB7O/sb0X4xOY4Awy1kn0AbCSiR7s+d3bt5koAh7rahzvLEFqHmZsApgBsSLjtklmlRJMJguCV21Qm4MfWrMsn1IWm7b72besVFA8AlWjx+V1XYweAohIOc+rpORP4nxGr4i2zLPuifdYA4JQqZHKmYdfRAewjWamfFQIa4agqODKqgrt1GwBGlC5XUTrcM3lZ0R0AJovSnywu2mvWKMlzqpTlOoUTVlfMnFX9a8obiOasrpg/LfdbyqkCQjnbt7gor9GhyqhoP1vZara5riCTQG6JTpp1SqR9D+X3IV33kpDceHCKmfdcwp70FX9ic5wBhTj5pwdHAOzoam/vLAuuQ0RZACNoGxGSbLtkfGJznEGGKdlncR4BsJuIriaiPNrGgHvVOvcCuKPz+1sAPMjM3Fl+OxEViOhqALsBfGelp+V+bI4zyPTBj42Zm0T0bgD3o+3ucQ8zP0VE7wfwKDPfC+BjAP6ciPYDmEB78kNnvc8DeBpAE8C7mNnG8y0Rn9gcZ4Dpl4MuM98H4D617Pe7fl8A8K8usu0HAHygPz1ps4bGAxW03FRB8QHtlBvyzXmhYbs/05BB4rM52a5E1vFUB5HravGhylb6HT5kLDDHUe0FZTyYChgPTtalQ+tEPWA8UP2tq0DzgskwEDAWZGR7PGcD53UV9I056Tg7nrcGhyeK0qH18NC4WadRkQafRlVes/Kovc5DJ6VBITslryvV7TmTqqKVm5PrFKas+1TtjLyuM2ek8/CBDTbT7YGqrKq1O3fKrFPIymtWNokLrFPySOGCQSTqx4zE1oiXFvyJzXEGGQ+pchwndfjE5jhO2vAg+JUQeI/XslUSB11dHV4nngRsYPxsrBJPtgIOokpjiwMB7Jq8CrLWDrq66jsAqNyaJmlkKMBda2on5ytmnZYyx2snZa3BAUBJaY3VjHS+3UDSgRcANijH33GluWkNDgCuyEut6J+LO8w636/IYPqpqgymrw9b7WthVGpdQ6flOefPWo0tU5fXpJWV+plOfgoAWeXnG83IvpyetZrn0fqoaE8M2QSiW1juuJqR13BH7rTZZlflQnD9Y5E9P+cC/sTmOIOMP7E5jpMq3CrqOE4q8Se25ROsBB9rPzb5fSgGOFOXmkgz4Memq3lrvWkmtsVQaiqpX4t7+6RZPzaty9l96H+Oxo+taTW2MwtSwzkzZ9dpxFL3mVG+e83A+WQzUtTUySmHM4HklMq3aktG6nQbMi+Yba7IygSQW/OTZp1tpatE+/GK9H07OjJqtqmPynOsnZJjUJiw55ybXdx3Mi5YbVUPnX7CaQTuwcmGvEbTLXvNYqiEAqpa/I7slNnmhvKFEMr7+hAkT3DjgeM4aSSlE1vPRxMi2kFEDxHR00T0FBG9ZzU65jjOJaZ/2T0uO5I8sTUB/Edm/i4RVQE8RkQPMPPTl7hvjuNcagbVeMDMRwEc7fw+TUTPoJ3h0ic2x1nnrMensSQsSWMjol0AbgLw7cB3dwK4EwCKkapElCAI3jroBo6vHHTjWsBxU4m5s9phN2eNB7oiUJ2XXksiIpWNNUHmFX3c2abtW3dlIgCYmbXOns2G7O98Vv4brjft+WinXp0lOB8InNfVu3aoylYjGSu+VzPS0XdD5gdmnS3ZSdHeVpDtveXtZpt9w9Kp9+yYdFxeOG1v7fyUvEbK3zhIQxXeahWUASJjH3m08SnJ/VQg2d8tAQfcGwoXjAdD1LvaWiJSOrElTjRJRBUAXwTwm8x8Vn/PzHcz8x5m3pPP9C4R5zjOGtOnKlWXI4me2Igoh/ak9mlm/qtL2yXHcVaLgX0V7ZTI+hiAZ5j5jy99lxzHWTUGdWIDcAuAtwF4koj2dpb9p07GzGVjHXRVu2H1Gu2TGNftm7QOjJ8zDru2SpWpZKUqhOsg+eWig+sbrPTApq5MDszWVF/mbP9pXgXgq6E7u2AvczOWYxe3ZDt0znFZj7d0yNWaGwBUlfa4PXDHVTMnRFsH02/NTZpttpek7vZMdYtoHxobNdvMnZESSUYFtIfuuVZO3petYal9lYtW6yqr6l35QKJSjdZoS4Fqatu6HKRzfYqFWo2QKiIaB/A5ALsAHATwK8x8Rq1zI4A/BTCMdkX4DzDz5zrffRzAz6Bdsg8A3s7Mexc7ZhKr6DeABOkuHMdZX6yefvZeAF9h5g8S0Xs77d9T68wB+HVm3kdE29B2K7ufmSc73/8OM38h6QG9SpXjDCi0hM8K6a4C/wkAb9YrMPP3mXlf5/cXAJwAcIVeLyk+sTnOIJPcKrpYJfhebO74wwLAMQC2snYXRHQzgDyAbt+gDxDRE0T0ISKyflEKjxV1nAFmCVbRRSvBE9E/ANgS+Op93Q1mZqKLH5WItgL4cwB3MPM5BfAutCfEPIC70X6Nff9inV3DKlWyndEZdQMOulFdPhQ3A2Jvsy4F4YW6Mh4UrECvM37obB9a5AeWF4kSszYeqL7G9ji1mlxGc1bUz87oTLDy+9asfTCfq8llzyujy3zQyCLHbroinYWnC0eh2ZkVGjE2RHbktIHhGmWEqGYOmm1s1hCZDePZsv0be35YVpQ6NaMy9dbs+OvRHi5Jw8C2inHpxCZVvauUsZXRNDG3VNv+zXRXT6N+iWP9K7/32ot9R0THiWgrMx/tTFwnLrLeMIC/BfA+Zn64a9/nbqwaEf0ZgN/u1R9/FXWcQaWTaDLJZ4V0V4G/A8Bf6xU6FeS/BOCT2kjQmQzPuZ69GcD3eh3QJzbHGWRWJ/LggwBeR0T7ALy20wYR7SGij3bW+RUAPw3g7US0t/O5sfPdp4noSQBPAtgI4I96HdA1NscZYFYj8oCZTwP4+cDyRwG8s/P7pwB86iLbv2apx1y7KlUmCF476AZ0OaWpZWpWY9OB8Vpjm6lbg4oOPp9rSS0pFMSsNRCtkSQhVg/M9UAFrVgFsEfzgcyw02pclM9oyL84OyePVV+Q+z0c0JumF+Q4nR6T2X1PVIfNNtcPSSfeXXlbFX1LJAPli+qyVgOC667chGgXSXpvj0TWWXhTQWpfRyqjoj1ZtzHOTeW4PFqQkfM7S7IfALCzIM9xNJCNOFJOFPMsL9p0y57z6a6qZs1+vWwNcOSB4zgpZWBjRR3HSSmMwU006ThOOvFiLitkWVWqQn5sKgg+FLSsA+Mb2j8rUFVIB8ZPx9I/S1eTAoCV1wiyCQm1ngMArYZKjhjIL6hjz7NzSr8M/FeO83LssnNK75u1/n5Ts3LsnlBJL4+OWY3toPId212xLkxXF06K9jbl+xaqmJVTznpl5Su2LSf3EdpmJJJ62VRsNTadkHNEZae8KqAZXpuT5zgeuJkzJIXPOZWY9FCgmtqzta3nf19o2eMuC5/YHMdJG6GHjjTgE5vjDCrrNDtuEnxic5wBxjU2x3FSx2okmlwL1jAIvlcGXbsbvSxTt8aDjDIe9HLYBYDZhhRqdVD8XCtQ2UqJyq0+PNMzBzJfqXS4IYNJtCCPnZ+W7Sjg7KwP1ZyW45azsd1onJVjWZuSYvuxSWtwODUqK5Y9PzZu1tmhst/uGJLC/5WqahUAjGdlll3toBsiUn/FI8rqUopssLo2OIwqx98dudNmm23KwDCa6f1nNq3+Hg41Nph1npq/kDV4np/ruc9E+BOb4zipYp1WeU+CT2yOM8j4xOY4TppwB92VEvKV6RUE37TbaE0tIIkgVoHxrR4Ou0CgWrwKip8NaGwN5VzbUA6WObKR55G6izJJlFs1DKGCR5Fy2s3Ny/1m5+xxSI0vZ+W4xQXrLNyYlMsKZ+Q2tRE7tvUROQ7HR+1YHh8ZEe3nhjeJ9saKDJIHgM1DMqB9PC+1L10pCgCKSqTV+lkhIOxq7a6a0fqZdR7WiTOHyGqPNZZOu7MqmekLjTGzzb7pCyUAQklJl4PWutOCP7E5zqDifmyO46QRd/dwHCd9+BPbCgi8x1PcUu3Fg+IBm3wyWC1e6U2kipboYi8AMKd826ZN4smQH5vcTwNSr9GJBENEKmdMsHhPgqKOpjBOXe4nOxsIwl5QepL+zx3ZAxfy8pybJXn7NKoBja0qx782Yse/PqwCwkfkeB8Ylr5wAPDDqvTzKpWl1lUtWo2tmpfLhvNym7G8TU65rSCLxOjCLI3spNkmo66jrvIOwEwoWrMNBeSfmq+c/z2UMGE5pNV44DUPHGdQYbQNe0k+K4CIxonoASLa1/lpLSPt9eKuegf3di2/moi+TUT7iehzncIvi+ITm+MMMKtUpeq9AL7CzLsBfKXTDjHPzDd2Pm/qWv7fAHyImV8E4AyAd/Q6oE9sjjOgnPNjS/JZIbcB+ETn90+gXUIvWR/bJfdeA+BcSb5E2/vE5jiDStLX0Par6EYierTrc+cSjrS5q+jxMQCbL7JesbPvh4nozZ1lGwBMMp93/DsM4MpeB1wzq6j+L2CC4EMOurpKVSCbbKQddJWG3KxZ8dpUi2+qiucqoy4AzLJcp8GLV1oCZCVvwDqI5qOA921GOdKGdGh9LO3U27D7zcwr40Fdtqlpt4lIZd3NydunMGSlj7isxqlqkxBo40FtWDk/D4eMEiozsq5KX7H9p5JyyC3LG2i0Yo0HZ6qyEldLWXM2qGB8ALgiktltSyGv6mXQ6DIYcBKrUgKW8DR2ipn3XHQ/RP8AYEvgq/d1N5iZKWglAwDsZOYjRHQNgAc7tUSnLrLuori7h+MMMn2yijLzay/2HREdJ6KtzHy0U9Xd5odv7+NI5+cBIvoqgJsAfBHAKBFlO09t2wEc6dUffxV1nAFmlTS2ewHc0fn9DgB/bfpBNEZEhc7vGwHcAuBpZmYADwF4y2Lba3xic5xBhQHEnOyzMj4I4HVEtA/AazttENEeIvpoZ50XA3iUiB5HeyL7IDM/3fnu9wD8FhHtR1tz+1ivA65hEHwvB127jU6YqIO/AaClAuXjHg67gA2M10HxOvEkYJ12Gyy1Fl0pHrD/RXKqwnk+VM0oK8eJs4GkkeoqsnauzQT0GNU/qkmNjRcCGQaaAa/p7n1EgSr1BTlO2SE7lsWy1MeaVTW2w711uXqVVNve2o2KdiiW1/nYqNVSpzeovinHWF3pCgA2ZKTeWqVQAku5n6J6LKpENri+lLtwjbQT8HJZDQddZj4N4OcDyx8F8M7O798E8LKLbH8AwM1LOWbPJzYiuoeIThDR95ayY8dx1gGr4KC7FiR5Ff04gFsvcT8cx1kDVkljW3V6TmzM/HUAE6vQF8dxVhNewmed0TeNreOwdycAFCMVtJygEjw1pZakA97by1Q7oLHpZdqvTSeiBIBYBcbP1ntrbDr5pE5R2DJR5TauXCc+LAXKvEc56QPVCkTJxQV1jsqJrlWwvnuZnFxGWbuOhpXGxnXV3zjgr6W0JMrZW46UDpcvyQDwXNkGhBeH5bKm0svqgaSXtarylxtViTJnQoV+yqK9P7NRtMfyVmPTVeg3R4fMOuPqZtC+bluy1n3ryvLk+d+fCGWJWCIE+3eYFvo2sTHz3QDuBoCR/OZ0jpbjpAyvBO84TrpYp6+ZSfCJzXEGlvVp8UxCEnePzwD4FoDriOgwEfVMGeI4zvogrVbRnk9szPzWFR8liYNuj8rwQMhB1xoCWjW5jhbWQ9XjY+W0O28y6gaC4JWKryvDh9D/RXQFpErAeFAoSJF4rmjHpTkkj90oySNlKwHBXp+TzsIbumasVlLXkAPGA27UF20DAOalAE9zMhg9M22NN9G0DE7PVGQ7O2kNDnlVISs3J69zVA/9n5djN52XxoQDZVux/YdlaWDYnT9u1hlV4n9VOVHvyslAegC4vnxhPw/1wXgAILVPbP4q6jiDCrtV1HGcNJLOec0nNscZZNzdYwUEq033dNC1Dq4ZpYFEebvfuEe1+JAu11T7rTdUUHzTesXaIHi5jzjwrzCvEjWWVcWjcqC0fbkgNamZIatjNcuq/yognGLrfMuqL3lVCT4KOfUW5Dhk5mSgNs9bZ1WuKY1NO/UCtoqZamvHYABWl1N/oFFs759ef8StrL3OzZIKrp+W98bErNT2AOBEXTqoT7as3hfjrGhXM/LY27J2LF86dMHRdyjknb4cfGJzHCdVMGzZxZTgE5vjDCgE9ldRx3FSSCudj2xrN7HpRIe6MnyomEtT+7EFfN1yiyeazATyJ2rftkZDVSZPoLHpyvAtDiSNVG3txzYS0FWqBdnhk2W730ZFHltrkToQHQBaapya2vctUHQlO1dUbeVLpgvEACBVcT5TC2hDunCM9ocLPVVESgNUQfycIKg/ie+klrL0vaL1WMAmTVhgO5aaAsl1xjNWS702d7pr/T74saX4VdRTgzvOAEPMiT4rOkaCSvBE9HNdVeD3EtHCuRJ8RPRxInq+67sbex3TJzbHGWRWJ4Nuz0rwzPzQuSrwaBdIngPw912r/E5Xlfi9vQ7oE5vjDCwJJ7WVT2y3YWmV4N8C4O+Y2RZ6TYhPbI4zqCytStVqVII/x+0APqOWfYCIniCiD50r07cYa1elShsLejjsAkCmrpx4s9bZNsotXslKV5MPLWuqjLrzTSv+zrV0JXi5TWxy6gKRqt5dVsaDsaysbgQAYwX5T2uoZMX3WVUpvd5SGXRz9pwbZVXVfV47NgfGdkGeY7amgshr9jpHtdaibcBeV2NICjl462UqiJwDlbl0JuFmSbWLgW2U3YgjZXDIXBr1PUfW+FHqqnaV6VMs1BL0s9WoBI9OQeWXAbi/a/FdaE+IebST2f4egPcv1ll393CcQaZPfmz9qATf4VcAfImZz//X73raqxHRnwH47V798VdRxxlUGO0n3ySfldGzEnwXb4V6De1MhiAiQluf61kK1Cc2xxlYVs14kKQSPIhoF4AdAL6mtv80ET0J4EkAGwH8Ua8Drp3GppfpIPiAxsYNHSgf0IHUMh30Hqoer50wuSHn+4WmHaaFltSXtBNmC7aSd04FnhdUZaLRyBqBNhak7jZask68uiq6HpZ4yP7/yigNLdNUgfNWIjRjSc0EY6ucnRNVH1O+pxQqfqWeIpJkeVV5CtBUSUibZXs/1UdkO67Ie7BctCc9nJXXXjtih4hVEs+FgIP3NF+4D2PYvi6LVQipSlIJvtM+CODKwHqvWeoxXWNznEGFYYx4acEnNscZWNime08JPrE5ziDj2T0cx0kV56yiKWTtjAemSpWqeBQwHtgsuwFRXGX80GJ1KPGoqVylRPGFQAaH+VgaC+rcO5tEDnKdsupMNWMNA+M5aTzYOGSdeOsqQ+5MVo5To27734qVcK6cerkVEKdjbZXo0QZA8eIGh9AyndgiaDxoymtP6nbR7RD6kjVDFcCG5cEL4/IabavITLgAsDU/JdrVjDUkaWfteZb3wqlAxa+DjU3nf6/zafP9svAnNsdxUodPbI7jpApmm/suJfjE5jiDjD+x9RftYMkmKN7+J8koXUU77ALWaVc77AaD4LVDqFqnEajyVFMOug2ln8WBGyZSgdlF5bAbdNDNTYv25iGr6TSV52kpJx1CG61ABl1VuT6jPFyTuH/qM9T7BIBYHTsOaHdNNb5mm9j2P1b3Qktt0wpoeVo3VMOPTM7eT5WSzGC8bViO/+6qDXu8piCXbQikbc6oIPfplrwJDzaVZzCAJ+avOv/7fOuH5vtl4ROb4zjpoi9xoJclPrE5zqDCALuDruM4qcNDqlZA6HFXOxppv7aAH5t2cKJswI9N+bZpP7YoEI+sC/5ov6qQxlNryaFrcO+hzKhkKjmTeNI62Y1HM6J9RX7GrNNSGttszlbV0mTVWBZUeygQ0Z5TDmVZlWQxSlDyKA4klNHanB5b3QasH+FCjzYA1FQyA61NZgPObyMF6be2szQh2jeUjphtdmWlj9l4xp6z9mObVmNwqLHBbPPszIU8jjoJw7Jg9vJ7juOkEDceOI6TNtif2BzHSRd9SSJ5WeITm+MMKh4Ev0ISZNClphSmWXtPAkBGGwasEy9n9Tqq3QxkcNVB2ImMB9LBUgfBJwlUKZAc/pK2YsA67Y4HKllp8X04K/erRX4AKKkA/JGsPE4ocLusHE1Lqp0PRKvnAuekiZWIr52dQ4aZ2ZaswDYdF0V7Ki6ZbfQ6IaOEphLJc9yel8aDa/LWQXezslBVMkNmnQbLsZpTzsMnm1WzzdG54fO/11u9ky70ggFwSkOqEtU8IKJbieg5ItpPRKaKs+M46xDuJJpM8lln9PyXRUQRgA8DeB2AwwAeIaJ7mfnpS905x3EuLZzSV9EkT2w3A9jPzAeYuQ7gs2iXrHccZ72T0ic24h5WESJ6C4BbmfmdnfbbALyKmd+t1rsTwLmy9y9Fgtp/lwkbAZxa604sgfXU3/XUV2B99fc6ZrZC3BIgoi+jfc5JOMXMt67keKtJ34wHzHw32uXnQUSPMvOefu37UrKe+gqsr/6up74C66u/RPToSvexniaqpZLkVfQI2kVMz7G9s8xxHOeyJMnE9giA3UR0NRHlAdyOdsl6x3Gcy5Ker6LM3CSidwO4H0AE4B5mfqrHZnf3o3OrxHrqK7C++rue+gqsr/6up76uOj2NB47jOOuNRA66juM46wmf2BzHSR19ndjWU+gVEd1DRCeI6LL3tyOiHUT0EBE9TURPEdF71rpPi0FERSL6DhE93unvH6x1n3pBRBER/TMR/c1a96UXRHSQiJ4kor39cPtII33T2DqhV99HV+gVgLderqFXRPTTAGYAfJKZX7rW/VkMItoKYCszf5eIqgAeA/Dmy3hsCUCZmWeIKAfgGwDew8wPr3HXLgoR/RaAPQCGmfmNa92fxSCigwD2MPN6cSZedfr5xLauQq+Y+esAJnqueBnAzEeZ+bud36cBPAPgyrXt1cXhNudymOc6n8vWSkVE2wH8EoCPrnVfnP7Qz4ntSgCHutqHcRn/8a1XiGgXgJsAfHuNu7IonVe7vQBOAHiAmS/n/v4JgN8FEhRtuDxgAH9PRI91QhkdhRsP1hFEVAHwRQC/ycy2cvJlBDPHzHwj2pEqNxPRZfm6T0RvBHCCmR9b674sgZ9i5h8H8IsA3tWRVZwu+jmxeejVJaSjVX0RwKeZ+a/Wuj9JYeZJAA8BuFzjEm8B8KaObvVZAK8hok+tbZcWh5mPdH6eAPAltGUgp4t+TmweenWJ6IjxHwPwDDP/8Vr3pxdEdAURjXZ+H0LboPTsmnbqIjDzXcy8nZl3oX3PPsjMv7bG3booRFTuGJBARGUAr8f6yaSzavRtYmPmJoBzoVfPAPh8gtCrNYOIPgPgWwCuI6LDRPSOte7TItwC4G1oP03s7XzesNadWoStAB4ioifQ/of3ADNf9m4U64TNAL5BRI8D+A6Av2XmL69xny47PKTKcZzU4cYDx3FSh09sjuOkDp/YHMdJHT6xOY6TOnxicxwndfjE5jhO6vCJzXGc1PH/AePD6duM+VbXAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(z, origin='lower', extent=[0, 5, 0, 5])\n",
    "plt.colorbar();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is a compelling visualization of the two-dimensional function."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "jupytext": {
   "formats": "ipynb,md"
  },
  "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.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}