Advanced Numerical Linear Algebra/week_1_introduction_to_matrices.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Why are we here?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note: Future lessons have a lot more code than this one**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Why study Numerical Linear Algebra?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    **Key Question of this course**: How can we do matrix computations with acceptable speed and acceptable accuracy?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A list of the [Top 10 Algorithms](http://www.cs.fsu.edu/~lacher/courses/COT4401/notes/cise_v2_i1/index.html) of science and engineering during the 20th century includes: the **matrix decompositions** approach to linear algebra. It also includes the QR algorithm, which we'll cover, and Krylov iterative methods which we'll see an example of. (See here for [another take](https://nickhigham.wordpress.com/2016/03/29/the-top-10-algorithms-in-applied-mathematics/))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    There are 4 things to keep in mind when choosing or designing an algorithm for matrix computations:\n",
    "    - Memory Use\n",
    "    - Speed\n",
    "    - Accuracy\n",
    "    - Scalability/Parallelization\n",
    "\n",
    "    Often there will be trade-offs between these categories."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Motivation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    Matrices are everywhere-- anything that can be put in an Excel spreadsheet is a matrix, and language and pictures can be represented as matrices as well.  Knowing what options there are for matrix algorithms, and how to navigate compromises, can make enormous differences to your solutions. For instance, an approximate matrix computation can often be thousands of times faster than an exact one.\n",
    "\n",
    "    It's not just about knowing the contents of existing libraries, but knowing how they work too. That's because often you can make variations to an algorithm that aren't supported by your library, giving you the performance or accuracy that you need. In addition, this field is moving very quickly at the moment, particularly in areas related to **deep learning**, **recommendation systems**, **approximate algorithms**, and **graph analytics**, so you'll often find there's recent results that could make big differences in your project, but aren't in your library.\n",
    "\n",
    "    Knowing how the algorithms really work helps to both debug and accelerate your solution. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Matrix Computations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    There are two key types of matrix computation, which get combined in many different ways. These are:\n",
    "    - Matrix and tensor products\n",
    "    - Matrix decompositions\n",
    "\n",
    "    So basically we're going to be combining matrices, and pulling them apart again!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Matrix and Tensor Products"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Matrix-Vector Products:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true,
    "hidden": true
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "#### Matrix-Matrix Products"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<img src=\"images/shop.png\" alt=\"floating point\" style=\"width: 100%\"/>(Source: [Several Simple Real-world Applications of Linear Algebra Tools](https://www.mff.cuni.cz/veda/konference/wds/proc/pdf06/WDS06_106_m8_Ulrychova.pdf))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "#### Image Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Images can be represented by matrices.\n",
    "\n",
    "(Source: [Adam Geitgey\n",
    "](https://medium.com/@ageitgey/machine-learning-is-fun-part-3-deep-learning-and-convolutional-neural-networks-f40359318721))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "#### Convolution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "    *Convolutions* are the heart of convolutional neural networks (CNNs), a type of deep learning, responsible for the huge advances in image recognitionin the last few years.  They are now increasingly being used for speech as well, such as Facebook AI's results (https://code.facebook.com/posts/1978007565818999/a-novel-approach-to-neural-machine-translation/) for speech translation which are 9x faster than RNNs (the current most popular approach for speech translation).\n",
    "\n",
    "    Computers are now more accurate than people at classifying images.\n",
    "\n",
    "\n",
    "\n",
    "    You can think of a convolution as a special kind of matrix product\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### Matrix Decompositions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "    We will be talking about Matrix Decompositions every day of this course, and will cover the below examples in future lessons:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "    - **Topic Modeling** (NMF and SVD.  SVD uses QR)  A group of documents can be represented by a term-document matrix\n",
    "\n",
    "  (source: [Introduction to Information Retrieval](http://player.slideplayer.com/15/4528582/#))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "    - **Background removal** (robust PCA, which uses truncated SVD)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "    - **Google's PageRank Algorithm** (eigen decomposition)\n",
    "\n",
    "\n",
    "  (source: [What is in PageRank?](http://computationalculture.net/article/what_is_in_pagerank))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "    - List of other decompositions and some applications \n",
    "  [matrix factorization jungle](https://sites.google.com/site/igorcarron2/matrixfactorizations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Accuracy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Floating Point Arithmetic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "        To understand accuracy, we first need to look at **how** computers (which are finite and discrete) store numbers (which are infinite and continuous)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Take a moment to look at the function $f$ below.  Before you try running it, write on paper what the output would be of $x_1 = f(\\frac{1}{10})$.  Now, (still on paper) plug that back into $f$ and calculate $x_2 = f(x_1)$.  Keep going for 10 iterations.\n",
    "\n",
    "This example is taken from page 107 of *Numerical Methods*, by Greenbaum and Chartier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    if x <= 1/2:\n",
    "        return 2 * x\n",
    "    if x > 1/2:\n",
    "        return 2*x - 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Only after you've written down what you think the answer should be, run the code below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.1\n",
      "0.2\n",
      "0.4\n",
      "0.8\n",
      "0.6000000000000001\n",
      "0.20000000000000018\n",
      "0.40000000000000036\n",
      "0.8000000000000007\n",
      "0.6000000000000014\n",
      "0.20000000000000284\n",
      "0.4000000000000057\n",
      "0.8000000000000114\n",
      "0.6000000000000227\n",
      "0.20000000000004547\n",
      "0.40000000000009095\n",
      "0.8000000000001819\n",
      "0.6000000000003638\n",
      "0.2000000000007276\n",
      "0.4000000000014552\n",
      "0.8000000000029104\n",
      "0.6000000000058208\n",
      "0.20000000001164153\n",
      "0.40000000002328306\n",
      "0.8000000000465661\n",
      "0.6000000000931323\n",
      "0.20000000018626451\n",
      "0.40000000037252903\n",
      "0.8000000007450581\n",
      "0.6000000014901161\n",
      "0.20000000298023224\n",
      "0.4000000059604645\n",
      "0.800000011920929\n",
      "0.6000000238418579\n",
      "0.20000004768371582\n",
      "0.40000009536743164\n",
      "0.8000001907348633\n",
      "0.6000003814697266\n",
      "0.20000076293945312\n",
      "0.40000152587890625\n",
      "0.8000030517578125\n",
      "0.600006103515625\n",
      "0.20001220703125\n",
      "0.4000244140625\n",
      "0.800048828125\n",
      "0.60009765625\n",
      "0.2001953125\n",
      "0.400390625\n",
      "0.80078125\n",
      "0.6015625\n",
      "0.203125\n",
      "0.40625\n",
      "0.8125\n",
      "0.625\n",
      "0.25\n",
      "0.5\n",
      "1.0\n",
      "1.0\n",
      "1.0\n",
      "1.0\n",
      "1.0\n",
      "1.0\n",
      "1.0\n",
      "1.0\n",
      "1.0\n",
      "1.0\n",
      "1.0\n",
      "1.0\n",
      "1.0\n",
      "1.0\n",
      "1.0\n",
      "1.0\n",
      "1.0\n",
      "1.0\n",
      "1.0\n",
      "1.0\n",
      "1.0\n",
      "1.0\n",
      "1.0\n",
      "1.0\n",
      "1.0\n"
     ]
    }
   ],
   "source": [
    "x = 1/10\n",
    "for i in range(80):\n",
    "    print(x)\n",
    "    x = f(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What went wrong?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Problem: math is continuous & infinite, but computers are discrete & finite"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    Two Limitations of computer representations of numbers:\n",
    "    1. they can't be arbitrarily large or small\n",
    "    2. there must be gaps between them\n",
    "\n",
    "    The reason we need to care about accuracy, is because computers can't store infinitely accurate numbers.  It's possible to create calculations that give very wrong answers (particularly when repeating an operation many times, since each operation could multiply the error)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How computers store numbers:\n",
    "\n",
    "\n",
    "The *mantissa* can also be referred to as the *significand*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "IEEE Double precision arithmetic:\n",
    "- Numbers can be as large as $1.79 \\times 10^{308}$ and as small as $2.23 \\times 10^{-308}$.\n",
    "- The interval $[1,2]$ is represented by discrete subset: \n",
    "$$1, \\: 1+2^{-52}, \\: 1+2 \\times 2^{-52},\\: 1+3 \\times 2^{-52},\\: \\ldots, 2$$\n",
    "\n",
    "- The interval $[2,4]$ is represented:\n",
    "$$2, \\: 2+2^{-51}, \\: 2+2 \\times 2^{-51},\\: 2+3 \\times 2^{-51},\\: \\ldots, 4$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Floats and doubles are not equidistant:\n",
    "\n",
    "\n",
    "Source: [What you never wanted to know about floating point but will be forced to find out](http://www.volkerschatz.com/science/float.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Machine Epsilon**\n",
    "\n",
    "Half the distance between 1 and the next larger number. This can vary by computer.  IEEE standards for double precision specify $$ \\varepsilon_{machine} = 2^{-53} \\approx 1.11 \\times 10^{-16}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Two important properties of Floating Point Arithmetic**:\n",
    "\n",
    "- The difference between a real number $x$ and its closest floating point approximation $fl(x)$ is always smaller than $\\varepsilon_{machine}$ in relative terms.  For some $\\varepsilon$, where $\\lvert \\varepsilon \\rvert \\leq \\varepsilon_{machine}$, $$fl(x)=x \\cdot (1 + \\varepsilon)$$\n",
    "\n",
    "- Where * is any operation ($+, -, \\times, \\div$), and $\\circledast$ is its floating point analogue,\n",
    "    $$ x \\circledast y = (x * y)(1 + \\varepsilon)$$\n",
    "for some $\\varepsilon$, where $\\lvert \\varepsilon \\rvert \\leq \\varepsilon_{machine}$\n",
    "That is, every operation of floating point arithmetic is exact up to a relative error of size at most $\\varepsilon_{machine}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### History"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    Floating point arithmetic may seem like a clear choice in hindsight, but there have been many, many ways of storing numbers:\n",
    "    - fixed-point arithmetic\n",
    "    - logarithmic and semilogarithmic number systems\n",
    "    - continued-fractions\n",
    "    - rational numbers\n",
    "    - possibly infinite strings of rational numbers\n",
    "    - level-index number systems\n",
    "    - fixed-slash and floating-slash number systems\n",
    "    - 2-adic numbers\n",
    "\n",
    "For references, see [Chapter 1](https://perso.ens-lyon.fr/jean-michel.muller/chapitre1.pdf) (which is free) of the [Handbook of Floating-Point Arithmetic](http://www.springer.com/gp/book/9780817647049).  Yes, there is an entire 16 chapter book on floating point!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    Timeline History of Floating Point Arithmetic:\n",
    "    - ~1600 BC: Babylonian radix-60 system was earliest floating-point system (Donald Knuth).  Represented the significand of a radix-60 floating-point representation (if ratio of two numbers is a power of 60, represented the same)\n",
    "    - 1630 Slide rule.  Manipulate only significands (radix-10)\n",
    "    - 1914 Leonardo Torres y Quevedo described an electromechanical implementation of Babbage's Analytical Engine with Floating Point Arithmetic.\n",
    "    - 1941 First real, modern implementation.  Konrad Zuse's Z3 computer.  Used radix-2, with 14 bit significand, 7 bit exponents, and 1 sign bit.\n",
    "    - 1985 IEEE 754-1985 Standard for Binary Floating-Point Arithmetic released.  Has increased accuracy, reliability, and portability.  [William Kahan](https://people.eecs.berkeley.edu/~wkahan/) played leading role."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    \"Many different ways of approximating real numbers on computers have been introduced.. And yet, floating-point arithmetic is **by far the most widely used** way of representing real numbers in modern computers. Simulating an infinite, continuous set (the real numbers) with a finite set (the “machine numbers”) is not a straightforward task: **clever compromises must be found between, speed, accuracy, dynamic range, ease of use and implementation, and memory**. It appears that floating-point arithmetic, with adequately chosen parameters (radix, precision, extremal exponents, etc.), is a very good compromise for most numerical applications.\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    Although a radix value of 2 (binary) seems like the pretty clear winner now for computers, a variety of other radix values have been used at various point:\n",
    "\n",
    "    - radix-8 used by early machines PDP-10, Burroughs 570 and 6700\n",
    "    - radix-16 IBM 360\n",
    "    - radix-10 financial calculations, pocket calculators, Maple\n",
    "    - radix-3 Russian SETUN computer (1958).  Benefits: minimizes beta x p (symbols x digits), for a fixed largest representable number beta^p - 1.  Rounding = truncation\n",
    "    - radix-2 most common.  Reasons: easy to implement.  Studies have shown (with implicit leading bit) this gives better worst-case or average accuracy than all other radices."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### Conditioning and Stability"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "    Since we can not represent numbers exactly on a computer (due to the finiteness of our storage, and the gaps between numbers in floating point architecture), it becomes important to know *how small perturbations in the input to a problem impact the output*.\n",
    "\n",
    "    **\"A stable algorithm gives nearly the right answer to nearly the right question.\"** --Trefethen\n",
    "\n",
    "    **Conditioning**: perturbation behavior of a mathematical problem (e.g. least squares)\n",
    "\n",
    "    **Stability**: perturbation behavior of an algorithm used to solve that problem on a computer (e.g. least squares algorithms, householder, back substitution, gaussian elimination)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Example: Eigenvalues of a Matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[   1. 1000.]\n",
      " [   0.    1.]]\n",
      "[[1.e+00 1.e+03]\n",
      " [1.e-03 1.e+00]]\n"
     ]
    }
   ],
   "source": [
    "import scipy.linalg as la \n",
    "\n",
    "A = np.array([[1., 1000], [0, 1]])\n",
    "B = np.array([[1, 1000], [0.001, 1]])\n",
    "\n",
    "print(A)\n",
    "\n",
    "print(B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true,
    "hidden": true
   },
   "outputs": [],
   "source": [
    "np.set_printoptions(suppress=True, precision=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true,
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([1.+0.j, 1.+0.j]), array([2.+0.j, 0.+0.j]))"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "wA, vrA = la.eig(A)\n",
    "wB, vrB = la.eig(B)\n",
    "\n",
    "wA, wB"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "hidden": true
   },
   "source": [
    "**Reminder: Two properties of Floating Point Arithmetic**\n",
    "\n",
    "- The difference between a real number $x$ and its closest floating point approximation $fl(x)$ is always smaller than $\\varepsilon_{machine}$ in relative terms.\n",
    "\n",
    "- Every operation $+, -, \\times, \\div$ of floating point arithmetic is exact up to a relative error of size at most $\\varepsilon_{machine}$ \n",
    "\n",
    "Examples we'll see:\n",
    "- Classical vs Modified Gram-Schmidt accuracy\n",
    "- Gram-Schmidt vs. Householder (2 different ways of computing QR factorization), how orthogonal the answer is\n",
    "- Condition of a system of equations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### Approximation accuracy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "    It's rare that we need to do highly accurate matrix computations at scale. In fact, often we're doing some kind of machine learning, and less accurate approaches can prevent overfitting.\n",
    "\n",
    "    If we accept some decrease in accuracy, then we can often increase speed by orders of magnitude (and/or decrease memory use) by using approximate algorithms. These algorithms typically give a correct answer with some probability. By rerunning the algorithm multiple times you can generally increase that probability multiplicatively!\n",
    "\n",
    "    **Example**: A **bloom filter** allows searching for set membership with 1% false positives, using <10 bits per element. This often represents reductions in memory use of thousands of times. \n",
    "\n",
    "\n",
    "\n",
    "    The false positives can be easily handled by having a second (exact) stage check all returned items - for rare items this can be very effective. For instance, many web browsers use a bloom filter to create a set of blocked pages (e.g. pages with viruses), since blocked web pages are only a small fraction of the whole web. A false positive can be handled here by taking anything returned by the bloom filter and checking against a web service with the full exact list.  \n",
    "(See this [bloom filter tutorial](https://llimllib.github.io/bloomfilter-tutorial/) for more details)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Expensive Errors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*The below examples are from Greenbaum & Chartier.*\n",
    "\n",
    "European Space Agency spent 10 years and $7 billion on the Ariane 5 Rocket.\n",
    "\n",
    "What can happen when you try to fit a 64 bit number into a 16 bit space (integer overflow):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/jpeg": "/9j/4AAQSkZJRgABAQAAAQABAAD/2wCEAAUDBAgICAgICAgICAgICAgICAgICAgICAkICAYHBwcICAgHChALBwgaCQcHGCEYGhERHx8fBwsgIBYeIBwSEx4BBQUFCAcIDAkJDRIMDQwSEhISEhISEhISEhISEhISEhISEhISEhISEhISEhISEhISEhISHhISEhISEhISEhISEv/AABEIAWgB4AMBIgACEQEDEQH/xAAdAAACAwEBAQEBAAAAAAAAAAACAwABBAUGCAcJ/8QAUhAAAQMCBAMEBQgGBwUFCQEAAQACEQMhBBIxQQVRYQZxgZETIqGx8AcyVcHR1OHxCBRCUqSlFRcYI2KS0zNERXKTJTQ1Q3UWU1RjZISUldIk/8QAGAEBAQEBAQAAAAAAAAAAAAAAAAECAwT/xAAiEQEBAQEAAgIDAQEBAQAAAAAAARECEiEDMRNBUWFxgSL/2gAMAwEAAhEDEQA/APjJRRRBFFFEEUUUQRRRRBFFFEEUUUQRRRRBFFFEEUUUQRRRRBFFFEEUUUQRRRRBFFFEEUUUQRRRRBFFFEEUUUQRRRRBFFFEEUUUQRRRRBFFFEEUUUQRRRRBFFFEEUUUQRRRRBFFFEEUUUQfo7PkP7UOcWjhnrC8fr3Dhprribo/6iu1Vv8AsrXT/wD28O+8r7ZwbAQHOBz5b5gJbOostTGDujbkvV+GOfk+Hx8gnaz6J/juG/elY+QPtb9E/wAdw370vueg2ei1Uqcn2KfhjWvhFv6P3a4/8I/j+GfekbP0ee2B04P/ADDhf3tfetNsbdFtpMBIssfjhr4AH6OnbL6G/mHCvvaYz9G3to4SODW/9R4T9eMX9CaLCbbLTRpQRayv44a/ngP0ae230L/MuE/fEQ/Rm7b/AEJ/MuEffF/Rqk0zcRy7k1gv9a42Y0/nEf0Yu3Gn9CfzLg/1YxF/Zg7c/QfT/wAT4P8AfV/SCgeVu9aGGPKVGsfzY/sv9uvoP+Z8H++qv7L/AG6if6D/AJnwf76v6VCoLoqZt0RMfzUH6L3br6D/AJnwf76r/su9uvoP+Z8H++r+lxVBylpj+aX9l3t39B/zPg/31Cf0YO3X0H/MuD/fV/S3NdC56pj+aR/Rj7cD/gf8y4R98VH9GXtwP+CfzLhH1Yxf0me5Ie625jTkrF8X83Xfo09thrwX+ZcJ++IP7N3bT6G/mPCff+tr+jFZ0iIjms1VxNo0W/CHi/ne79G/toP+DD/9lwn74qP6OPbP6G/mPCtv/u1/QlzjrB+ClVH3nTVbnxxmv58H9Hbtj9D/AMw4X97VH9Hjtj9EfzDhf3tf0Df7ktzGq/ijfi/n+f0ee2H0R/H8L+9Kh+j12wP/AAj+P4X96X3+T7ELWztron4oTl8Cj9Hftj9D/wAw4X97Vn9HTtl9D/zDhX3tf0DpNAsNXWHerey1/NZvxyL4R/Ps/o6dson+hrf+ocL+9qh+jr2x+h/5hwv72v6CMBM6/OVZNL5T77lZ8IzeX8/P7O/bH6H/AJhwv72hd+j12wGvCNLf+IcL+9L+gT2Eam2v1bJD9b3v7IXSfFF8XwJ/Z57YTH9EfzDhf3tT+zz2wt/2Rrp/2hwv70vvZl/HW231KvS7ck/FDxfBP9nztfp/RH8fwz70qP6Pva/6I/j+Gfel97vdqRdJfon4oni+DXfo/drhrwn+P4Z96Q/1A9rfon+P4Z96X3ZUcDE6z7Ekkyb2JT8UMfDDvkE7WDXhP8fw370l/wBRPar6K/juG/eV9zucb/urM+esdwUnxxl8Q/1E9qvov+O4b95Qn5De1P0X/HcO+8r7YrOkQdNkh17b7LV+KQfFh+Q/tR9GfxvDvvKE/In2n+jP43h/3lfZxfHfoe9Iq2N99FJ8UHxv/Up2m+jf43h/3lUPkV7TfRun/wBZw/7wvsN52AS33m/gp+KI+Pv6mu0v0b/GYD7wh/qd7R/R38XgPvC+vq/LlokuGsnvU/HDXyJ/VD2i+j/4vA/66r+qPtD9H/xeC/119ZmIPek1dN/Ba/FFfKB+SftAP9w/isF/roD8lfHh/uH8Tg/9dfVJHxus9RxU/FB8uO+S7jo1wP8AE4P/AFks/Jpxsf7if/yML9VZfUFfnHislSTf2yl+OD5nd8nPGhc4J3/Ww3+qln5P+MD/AHJ3/VofVUX0pUnbN5pDp0k+e6z+MfNx7C8W/wDg3f8AUo/6iF3Yjio1wjv+pR//ALX0RVnmfNIcDuT5qeC4+i8IyS4nd7p74WgM6JWHacz4gtzCO/LdaxTO1/FevXllUxmi10dQhYzTn5p9KldY6uukp1Bstk87LdTZEWlKw1DbkZ6LfRpkfgsxYNgNrLVk306KmMtfZNAJAcfCdVi9KOm0jx9iaWHnYXgeSWwApjWDa1oPfMrlWoJp3P4hNpg7adVGUzFrqgXDqstrc6LbSownXW+nJW6nZCGQY6abIezWm+u/1K3G5hLLJ3vqraZHK6EMBPRKqH1fco4xEJTwT4JFA8lKqg622TXbCdUqoLHv0C3IuslV0Tfn71neCeS1VGgg2OqQ/wBXpaFoZnNMRB87eKQ9hkbTpF4T67p0CTiRb3de/mtzpislUGbKg345poAtOvsRFnLXpaAr5NB9FZEykf2SiE6m+n4oCTtzlL22dJnu5c+itzHjUkA6QOSS07g669Fbau0yeunVZvSznTGhw5GTZA8Dl+SJ1SI5TIKRUcTGwH5pEsG9vIrNUBlObfSY/BJqDS569y1rLI8kEgaHVSbDqreAbwRtyS4Me7orKKLhoNSUovIm/wCShIBGp6oKjh175v3JahTxJ1sCDbvQVT15+9XUfy96Q65vskqWrc+AZus9R5EHojeTokv/AAVZKe/2pbzAJgdTv4InHRKde/groS8mxFvs5pdRx+OSZUZ1j60stnUdEC3tJ03SjTMfinkeXLZZ3gSbKJS6oEWv36pM+1Pe2yQR0hSIz1fFJdZPq9SkvWljLUGvl9azvJ0Wmo3n4pL+ntWVZag5rJWp6i+kra88/wAlnrMvqrRiqjvSaoun1R5C3XvWd9vdBWRlqCPNZnla64sstQQSmLK+mKbR6oEC3K+kT3p1Fulvjql4Sk4TJmdOYW7D0CdPGdVrrqPIOhTB21+LLZTpRFp2Q4ejsfgLdSYQbBc7WuRMp/vDyWqjSjuQlPolZtdIJlOU1tMxeY77JlDyTxSsdYJWNWECkLTbqEdOmLkRBTMvId07Ii02iLahZaAwkG0pu55BLGtz5WVhonUwi6IE3kRGl1QfJkBBUbYoc0W840RRCdvbeyqoqJMeKF5I5ePxqihe6LRc6fggBO8hU8k6n7fBDczr4qwR5S80zrzVOn8UhwvO3etxIXXqnqPj3JNWTF5PsTXtmZ35qhTtb43RrYQ5tiN5VGjNjNltY0G5VVaon1QDa6eSubUZASQ+8xYWiVrxzyNWwCuW98Sdvi6nk3hrqsm9gLQqdVnS34LBUqHz5IRWI6d6s6i5rc5/f4WUa4mO7RYTW/xJjXbzJAidlG76je11vqQVH6cllZUPO6jqvx15rTlWh9WNNR5JdSuTp8DVJc/Qb7lLL41+N1WR1DKU7TdC6tFvbKTUqC9ygt7+WoQ1XmPylU8zOk8wkF/cVUQu84SKjeuuoTajeYHO2sLO8xzKrFUec30Wd7jomvFrd6UVqBYMEW6Jbjp8bonNO53lLImAO6fyVAudPhp+PNLcT+SMj2FJcdZWoBcZkT9qXU8EbhZZ3W9qylR5gbykvMpjzBSnn46KRCX8/Kb2SahHNOPn37dEl4I2F1vmhFaVnqbp9Qfks1Qws1dKeFlrR5BaHFZ6rd+atVles+IbF99+5aqjfiUl+3VYGGqkVButlQSYWWppGv1JSPqmlRtcLXSpAR8eapogSdk1iy8x9KlF9eXcntAS2d/VMaDYxbT61Mbkw1nLZOpsQMF4Gy10h4LKjoifBamVNvJJpiO9NidNlzrQkDu+Dv16I2+1A436hFAmMEnRWATqiAhAFRsfWlZe/qE/qUpz9xEct0aDB6/gge1EXxpv+yUsVOnX4KKB406aISw5Y5o84B5cjyVVHk9VYpD9B01WY3m2u61VHWSWnVaZKFP3bom07W0VwnU6U2JUtXACkIidUdPDR+wTJseiaaUjYRz0S8PxXC7VmToRm3CzW5A4jBF4h0CNF5nieGNMmbTp1ScV2mrekPreq1+gGt0XFcd6d2c+rAsO8LM9ukuORVqb9I5bpNOtrI9v2p+KYADdc6sQNVvwZ8pG+m8DkURxTQSOa4r6hHcoKpN4W+ecZveu2MU3VQ4kSuGzEcgjGIPPRakZ12P1kKqlaVyXVtyVPTRBGm6LrpFw5ylF0wkmuHQRorkeSJ5Cc481QfEpebb2KMf4IyIH4+NEB7kJd5JbqnkqKc/8kl0o3QlOdcX2VlAudzSs/LY/V7UVWpzSS4HvGncqyokhJcic/fZLJVUNTbVLnx1TCRbkNUioTJ5bIKdyS536K/eqeb+CIQ7dLfoFbjZJefeqFVHbJBCeW6pNQRpqiUtwhIe2dh4rRFo3S3Ka1GR7Gjbe/ks1UaX30WurN40WaobdUtGSoDJ8llqtieZWyoTB71lrqfofWAafzRtBmTdNp0xrPfKY1m9ljXHnkbGzEJ7Bsgo31stI2hTWotjPsWigw+G6BkaEx3LS3WLxCxauHNCNqodytvJZaR56oQ6NkwgIQEEAVE3mEXsUbCkaUDF0tzeQsjkdUOvRVcZ67CIIAS3ag3jktfulJqmfNE0qoCRdo6DolE2vaE17xmtyQlu096rTO4bFKe0jS60VAgbGk3TQovHK6bTdG5vYoH2uhc7NI6H3J9o8l8oHHzSNOnRqOaQ2Kkb96/N28ZeysHxebjmvQ9uOFPo1nPkuY4znN41suHwbhratQ5xmbElPFryb8PxX0hzfta2tqtf9JGPWfcLiYuiKbyGi0x4DRZX1tV0nEYtegfxM85CA41pF7nmuC2vG/mocR3LeRHYdjBz81nOI0PnC5j8QD8X8FRryJnzS8jrnESNYU/WBGq5Irz0TDWkWOmvJDXU/W46j61oZWBHzvjl3LhvrRHfCKjiInlNx0/NDXWFfJJFxy5dQtVOs1zZGbrdckVfHfwS6eKyEQTDrEbApIrtmpurFTcLIHyL6j2o6b9tCVMGguQVHJZPjt9at52zRv+CCFxHX3pTjvp13VuIHM9Up/wASqI907JThvItaNwm6ckp8HX2LTJbjYzISyB7ET/GEs76oAfCB1j0RIHiUA1Cl1DZMiUt4SBDuaS4SnvCU4bc1oJc0+aSQeQWipa06JTiAs2pShaZWd6bVekPBuVPtqEPHks1R21vBPqGVnduoM9S9llqNMrU9tz3LNWPxsiPqtte15up6fbSAsNCq43W+hQLiJG11K48io4hx3vomsxbx5rVQwgEIxhddNVn02LC4jN0ldGjWMwuZTokarXh325rKz7dpkKySk4OoXCeSeFzranSfJRqsqKKtyp0fkrKoBUUR4KZUGYi2qE1SqJUgJFVw23uUZfrKzOqydERTxN48kDXRKt5KQ43ke1DTnGxSXkD8Ep1W8JbnjWCtYaKrVsliv6piNPeYVVKgNuiy1HAfHRXFcbtpiabaBa4Zi6zB11uvA8JcTXYGkgEwRzBXsO2+Z1FoYyb3PL7F4fBYg06rdocJKW4ldLi2GknlJ9i4VanuOq9bjmy0EEEOkg964FWhdv8AikLcvpXIIPelPGq3PpxMdyQ+jfUaT4qyssubnspTI3mFC2/sQVFb0sMa/wAPFFTfrfVLAnT3BDedN1NRqDrwjzX02j67rPGlt9lowzoBsCQZg62MFNbkacEdZ2Em99YhC62ZvO7eaGm4tc12UZb5o5HvRvf6zSOceHLvTRq4e9zmHcsMXOo+tac8O3uuXZtQOEi+k2vzWh1aNjGbyCsSuiKk/WhDtUllWCDaHe9R7yCD593NEPdzBug+LqZ+SF5lWIkoZF9JlWZSnmPj3qiO8bpb0RQlAsoXKyUB2RFPPJKeEZSnPCFpZSqhHjKJz1ne+471SKfqkG/4Jj3Xsk1HctVmhTzsEio6+qaXQkPI5KGlvO3NZ6oRVCkPcqpdUrNXTKjlnquKU19W4XDe/wBi6GFpwSb627lno1csXWzDVhGqxXHlups0+NkynTmddbJLH6XWyg/SDdYdIW7CjzTqGFIjSE6k2+m1+i0sELNrUVRp5bDRNaVRIS31FhTELnDcwsOJqnnCQCT3c1odP0ojW6F+JA6LnsOXdSQdborWytmnL9nvUzCNYIOvVZxFjEqOqWi3kgbWvJGu/RZag9qs1Iknfw9iW582F0Slk+xC6fzROjldZ3uMrSJWdusbz0066pzpNys9Vp0MStQKcOXW2u3NJeeenjyTqkx3ckioSdSYHNaajHj6jGU3F8FsG06r884lVbmJAaBcxvK6/bXGPbVAmGtEj3LxdcuqOJEnzXO9/osev4FixWZk/wDdWVV2AQY0d9aDs/w5+HYxxcf74A5YiPsWvEwc0iZd5XWv0SubXpDMfKOuq59VonTp+C7NZgzm9wR7lixtP5sbyo05TqcyCNO5Jq0rrS9u41lIqkG6upgaTOfmjbRu7bdLD9L3mPNOz+sJ0dYe5AbacAc9fBPFJskxMR7UNEg5h0MDuunYW7iP3mDKOo1SqzPYRpcA+r156og0WOsO7k+oAWtEXvH/ADD4KytMkjWbgK6GYunYHk4/ghZU9U31F9/ctNOnmaQ9wAyZgdZPJZabcpIgQ4SElDMFUmRoYkTudPctDrjvt4rnMBEXOuvXl5LW18kfuPHk8Gy1rNPw7zoeUeSM9Fm0JvGaP8wTWW1uDorKmCsDcd0FVVIN9FTiTpsgkFaQQS6h+PrVlxQv1QpZ5e3mqnyRP5JTzARA1J3KQACjc9Ie66CnG6XU+O9U915SKtZQU9Z6jkVV3ks76t+QUsRVUpNQi91VRw2KTVddBdQhZqpn42RVH69yyud7kaSq8WHJZahnojqnVZ33VI+sqVYbkeSdSfcQVhzx4rXhzpaVhynLr0KkRK34eoFxw/QAx7loo1TdZ6jcd70kDQo2VgVzaFfSb96dSqjW3mudajYTcH2fYh9OCY8ikio3aeisg67e1FwNZpcQCdELhEtjU6pnpR5bHXw6JOIqxt1QKe/KkHEg9EvFVLjlv9iyVB3lXDHTpV+XrIi52unvWClm31RmuSbnRMMac5J57SVTqgWSpiI6pRrbpi43mr6tkio/uWZtUmUh9SFZDGt9QAEzrsEqu+3MTqsr624JCW+vY3t71pLBvqXOq5/F8W+nTLmesR7LarQ+vItZZa7pBBEg2uqSvy7jeOqVXmbmbg9673Y2jQexzz61UatJ2XL7UYH0dUkHU37pW3sJGesG7jXpy7lwk/8Apq/T01VskHouZXuagi0D2Ls1gYHOdVzvRyXf8seMmV2tYc3Mc5MaQUBaTMjY+9an04zSOnkbIXgB4gasM/HNRqVxjRbMExAzD3QufWokECNS65XVLb2G512ScVTAIP7XMrn7dHIcyZt0CtjyQzoYW99NoYIs4uiO9ZcJTdBa4Qc5IPPkrOkw2k7K8SmB5a8f4TfuddAaZIzg62vsd1eIF37WHmBMremNIJBcJn9oLHiIY7ONdVqZSd88QYDSTyBsbbocSxsAiPVdlcdTfQwluGaZnkZv+UgJOJbuOfKRB2R0wRIBHRpH7O0KFsw0Wm99suquoN7AG/OaQ8elAaZhzPVynkYCAG9vmvHpGj/FoR3oJAzSPmuDidodyCItsWi7mHOw8wRok9oY5oc251kjaHjTuWzgeNZSe2pUotrjK8GnVMMLnMyjT4ssFB2o1Dhmb3j5w80TT6xG5Ae20afOgbJc/bU1oefnEazI80p03sIiyIdwQVHcvYunN9MdX2suj7Ep56FCXwJSK1WVqVk4myzVKsdyB9UrPUfJWUFWqefJA6prISqlRJfUTQdSoPD61lq1JVVHSNVnJ6poZUqrK5+3JW8pLz1SUxHPWd71HO/JJqOTyTEqPlIqvVvclPciheUio/ZE9yS42Upj61osMSfxWrCtKy0n36rbSeIupWTwy3tTWPj7NkplXxHL61ootzDMYtsPxUakMZUn3ptKpGqgpjWItZE9ukLnVaGu1WqjUae8LC1h2UcHzv4KK01nNJtr7fySK3VDkPU9d0BpE2lCfbPVbJht1baDhc+S20qIbEaonj4KrVYnM5hA9ka2Wuo4d/ck1b6jzSEYKl9jayU8xAuFoqU728b/AGLNWZIkzCqpUI2KQQbz4ImvAEbe1LfVzWK2gTH2DdBUGvJG0gGNeZ+xLqQCBeOqM6S7vSzP2J1SAlVDYnYb7KXrCPCdtix1SWu3hwOnhCrsPScKzzbIG3IRdpAPSGACHG/RN7J08tX9oDLYbFcp76bv09UXRrfcc1ka351lrzXnokmp6pC645sZZAMnUzz1SK4AM8hput7cvqicrDDXu1j1rmyOvh6OaoA/OQQKRykZxEEidAs7+nSc/twKlM5mi17/AJpNfDmSHDeI6dF0atOH7erM/hzQVGn1dzMk9NolZbYMZhqP6u5/pD6RroFPITLA0lxzDeViwjMzMzheARF9dO8LsYkNHpG+uGFrhfLmuNossGGLWk04PqBoHdE7brn9VfKWSMbqVgI0JBPfe3mmOpwSCNGgR36LeMLmJBgT82XDnqZV1cjRVBbmqEtGbZrWiSR/ildNTxc6kYa6JJc0t6t7oT8eBUYCGMZNJvzNT6Pdw3KHDNklskAzd2xOkxomsw803P8AVhvqkgwSCYIbN5Ut1cxhaHEWmWASBErQ4eqfrkcpv4oTR9G+Hf3YNPOCTdw5QN10qhFanRZTYWGlRd6QveGl5BJzBp1Fx5K3o4n3WbHcGfQa1znU3sdDC+m8OALmCo1hH70ELIwH+7eASRmpmNTsO4x7l2K3D3jC0XkS1zS9j6IJyuY4l7a0/t5fZC1disTSoYsem9GKdem5pe8A5Q9shw/dcuU+XJa69fFNjz2EwVR/pQ3IPQevd+U5ZGbLPf7F6bs5hqoq18NhnYOr6XDgvqVSMsOHrCkX3a/Vcrj1GizF1G0XmrRzhnpIn1XQTduo18lq4jUwFD0jKA/W2VKAa2q+WPoVJvlDNRCnXVv1+1nMjg17EiQYLgO4OgGd0mqY+xdDheMo0vTenw4xGei5lPM8g06hMteI2XIdU5T4mV6uerufx5u5M3+hqVEiq7fqpUqc0h9RdNcll/1pLn2VOqpD3pqLqvnRKqO6KVKnSEhzupWVxHuskOcre8T7Emo+9u9UwL6lkDn+xA959qU52qm4VHuSHOujcbmEs81EA5yW53O6MoHmFrVwk3nml5iAmPbugPTZQfXNDuWqmBFwsNM2Wqm9Kw1UWRvC2Ug2QOaxU38+SfQdrA0vKjTe11u4omt3Wak6T8arQ1y51WhtO9kxrdBzS6Rk38E5l/P6lFGaYIj2qhRgbHlKtolMpOGkaI0zubAvCSQStxbImIukVR3IVhcyOncsznk7krbWB5LI6nGysZJfbkCszn6j3LW+m7u6fGiQ5utrqtMlRvIpDzYxFtVsc0biFndIOjbLTOslSpeIlUSNNtU19KDZLcyIVQuoV5ntTxwUQ6m0w46x7l6So1eF7dUgKkkXiZ59Fz+T6a5+3Ew+IfVqBvzpI+Lr33D8D6MA2zR9V1+ccEw9apWaKZiHAz0lfqFOQGzqGmT4LHx+17qqhj7EFuWys03EEgfN2QAkGRbKNDou7BTZ8Npteb6IDykzNr2AlXiOp3n1dAeSqqw26i0X8+SzY3Ec1mSdXAls8jrKRVd6QtkHMcobHSw9y34GowCDSNTNd7wYOWZIA5whx2GBc6rSHoqTcrm+kqDMQXbABco6eP7JwXDa2IeabGgEAkmpDQ0cjO65FKlFepECAJ/dkZhae5dijUdJy5iXk5w506uEZlkpeq+sHgFnpWgiBIyiWj2lRbJhdDDseHekDjAluT2XKR+pPysqBoLC4gS4TLbmwM6BasdT9aWnKyA8R6oHNp5mUp2Im7oED1QGwXTqCO6VZKep6OpYN9alXrFjWUqRaDlj1nEeq0jXWFzajWuacurYMHctHrADdPp1wxlSn6MFzwRnD3SJMtOWYJgHzS6dCcrMzQZAYf2Q55vnesTferOpZP3XOdLqjPnOdIDbCejb2AXSqY01C52JDajyW5spyEAerksLaBMxvBzT/WZl/wCrmmw1qTgaIL4kzM7rDji4uEmk+GBoNP5paP3ifnFa56nVxq85P+rrY4jM2mX06RJcKYeXNacuV0ncrCwj1ffMm5mL6hdjhfFqdKk5lbC0sT6sUs393khxcdPnIWegxtZ16eCcfSVg53+zyxDcO1o/aOvgpJeb9emc859q4SytiXjCUn02DEEMPpPmerLhmIHqjXTkldouE/qcB2IpVKmZzXUqU+pl0vvK5tb0lMw6WHURLd7ETcaLJVrSS7e5J+cZPetzm75b9pOs58bKp9RIc/ZVUclPK9Dj0lR6z1HK3EpTwozFVOiS9XUOqXmQ5U8pTzf4hW96BxRSXnpHNAfqRuNo5pLrKoU4mfcgTHtUyjXZZpSCEJEJrmwluGiIU9qAjbomPCFxVUg9UuodSmuSaghXR9W0n2aFppVBzXLpVLLTh6lxdRwjr03z+Kex50gQLrmtdJB30+xPZUvBt1R1jp0qgOtlppPEa3XNbUBEeS0YYiyxWo6THea0UXHdYqL1oYdLLCxua6SBp47Jrmzppyjfms1N3insdKLiOok3IBSKjSIkSTotWaL+GntQPbbrNu5GmOq2PBIyW6rZVA6+KyuF9IhAh3TVZ6zJ5z8WT8Q4ja2/NZC8+3xRNA9h0ss1VpCfVudD1lZqxjRbiMztDIhKEJ1U+33pNQjcqlZ8Q6A5xk5QTC/N+1GLdUcJJzGY6DkvcdqsQW0Ja4NnnqQvzHG1QSTM6i64/N1+muY73Yemcz3R8wXK9eCYlcDsDRb6B75nM4bctQvQmjrddPj+mL9+iiSBHraweUJbzaJ+1Ma+AQTAdYjmmZBpFtgVtGR7bCDB6oTTvYmN/JaXUvj7FQpkkjciBFt1LFJw82IcA4gt9b61TAHTmaDsHToeic9kHL5j8Uv0c30uo1OsUajiwwwZmH1njl/i5lZano2H+7BOZvr+luCTu3kFqqjK4ta4lsattPTqsFQ6n1r6jbks+K3puwHBnVGF5xGGosYDHpX3vfK1ouSuXVdD/nEwYzatIGpE3AVOMQLX5/YU2tSgF4Z6RgpkkzES4Akc7qSWX7b8tjLULCCSAHAnQkEj9m26Q6qIiSDEE7E8k7E4osYWS1wqgT6slo2aHbLIMY/0Xoi2mBmzSG+sfFaZvtqwjazmYkUafpmBgNcNa45GtcC2IsHylcU4bUwzKdR5pllZmdr2OzNP+G+jvsTeHcS/VWF9OtUY+q4MfSpmA5gd6wfY2WYvq4irWqspF4IqVi0kZQAMrnNzQNAPJcc6nTvfG8/2hwXomVJxNN1SmaZORjw18keo8kdVXCcGzE+lb6T0dRjHPpgtL2lrbuBcNHWXK9M0G8kWMCRMtuiGNezMKTnU2v2Bg5IgtJ/dXW89X6rjOpL7M4jxCriKnpKr31HhoYHOn5jQA3XT8Via8AknTcIHzMyUD3wt885HO9b9rM6dY6JM6xstFXFg0m0vRsGV5f6QD+8IOgJWOo5JUqn3uNUp5I1nx0RGI6pVRaZoSlPRu70p7rSihdM6ygqH2oyEOX42QKLPtS3ck8goCzzQKa0jVDCY4wkucs0oXG/clONoTCSB3pSIWgeUxyzvN+iIF5KRUKc53lz5pFR3hf2KkfTdKoIHTVaKbrhcqnV0OxWunVv8QtOcjqsq8imseRqZlc9lVNbWM2i9lHSOnTrRy+xa6FWwM2XJpx1lPpviLws1qOzSrkn61sbVO5XCp4rKrqcQy6LGLHpKeIHNMdii2xjwK8sziB11sm4PGZnXGvMq4uvU0sW06+WyJtQm4INog/Gq4dGsT0Gi0sqgaOus9RrXUa05ZJvySagPwVkdjIGpASK+KJBFhdZkS1oq1gOp0gLJvob3Wc12MPU7zZaG1RGxVRnqucgqu2hPriSCICU+8rUoyuE+Cy1hZ1pgE381qcRyHKFy+NcQbSpvJIBLYDbfEK300/POPcZe8va8mA4gDYdAvJYmoS6GyS4iANdbBdHizXOc9+gJJhaexGD9JimOdENk33IvHevLL5dNdeo/Qez+G9FhqLC3Kcsu8b+a2RrdJNQeJvf2aIm1CNrL2SenLRPpSNghLT0smB5KWdeXRTQYcDqgfGwuLyPJCe+Ty+1GwGE0VWdYCL7lKsIB2TnG07rO5wBkp9qTWEnRZ679uqa6oFmq1GA7nu/FEIqDW+/j11Saj/W5iIyknKQdZG6Ks761jeRMnZLNINgpvHrVsgvAiXSNB3LJVdLB6xLgYPq2A2PUqqriZ6keASXG31bDuRot7okExoRG9wTdDjcQwvPos9OkRDWucSbgZm22mfNLqO2G+2yUZ35z4rUiahJ/dMCL6gDTwCFrS4gNvYzbYFF+sPax7A4hr4zDnBmehSS8j1gYMRa0+Sns2Aqnvn6kqoVHG3xogqOstJgC+58PYhUcVHnRNTAOshcZ/BUXEoXk/Wka0t7Z0PmhnaPJG8jcIHdyIoAKnO0AV6JRsfBESo6DslZjJKJzvd9aW5ZtA1CllyOqdkgn8FEVUd+SVsfjyROdulOcqKft7Ut5GnX2IX1EGbc8lTFF1h3pFZ47rRdG51vrWeu6RooPoajXad5Gy1U64mN/YvH0scToLrdQrPmZ6Qt6zj19Gv3HuWlj1wsJV0W6jXBKVOXXD48URJ2Kw0aomNlspkac1itylurOCsPcd040Z6o6VGPsWdUVIEBaaDY9afBU1rSRpbWVobh5i8CduUaKtSCD3EWMJrXG9yOt1ZaAbZuiMuGkmVnpaFz7TOugWIuPziSdikcVxLKYGYxyXmeJdo2wW0tZ+Cskj09TFspm7p9vvS3doqTbFwkneNIX51iuIPeTLzM3HVINSSJmOeq1hj9NbxxjzDXtJGg5rdgcb6TX2EL8ww+a2XbQru8DxTqZgyRumJHtcXWDWOcS0Q1xkCNB1X4zxfidRz3FznO9Z0E6ROy7XHuPVK1R7DLabXOaI3XDdh2Vh6rXnnB0815/k6125kcXE4/qNfgL1vye4Jr82Iky2WluwJvK8xh+EB1YMqF7KZd86NphfrHB+F08LSbSp/NsZjUkdFr4udY7q8mhjwVPmNbJ1WxS517l6Ncl4c2Vki/xCW/TdZcRWLT3jRRWwER71VTEMb0EbrE2qS3Q96y40EgCTHVFaqvEGDT3rFV4iSbLEaQEb+C0UsM06hMFtqHXVCRJmE4sGmw0QlsCVqMkuCRWpgplSoEsvF1RirNj7FkrVJ2jqttdqw1aZPco1GdsRb8UDimuYeXigIWolIIQgBNefsSaiUwt5hILvyTnjRJKzFC4wl1jumEoKhvdVATZA42vr9SupzQEhVA1HBAXclTz8dEsu89kFud1lA/6lRdCUXLOiybJedBUfslOeoT2Iv5pbnadxQF6VVPirDBmpPck1HaKidksm/JUxbnJT5nVFUcNUkn2+5NFVHeSS/nmVPfCW9/VQx+w4djY1v0Wuk+PxWN78tykVsf0BK0z6d2hiIOoXRw9bQyvEtx7wtWH4k8C5Pcl6PT3OGrj97fbXwXSo4sDffdfn9DiDjED2wurg8aS4SFNHr3cTAjc6StNCu9+9lwsG9hGi6YrtYBchZ8VvTr0mDmmPqlgytPSfauVS4g217afimUsay/rDXzQnbSceR6riUjGccFJrgLktS61RpEAgH2lYMZgg67ip461OnmOPcRqVXxLv+XZcgB0n1j4r01XgpLhqRqY5JGK4Q+ZDTB5qY1OnCc5NwoJ5rps4QReCtVLAOkWt3b9VU1lwrnW5aLuUcS0AA6peH4cSBIJttz28UY4e8GS0jkUo8d2spmnWcQSWkErif0nUZEAwb+xe84xw5r2OLyRkY6bexfneLGcwCeQ2kdV5u+crpKc3j1R7xmuBGvKV+v8CxdOphqbmGWkAeIEFfjHDOC1azjkDnQQHBouATqF+08G4ezDYalQYcwYBJPzsxEmVv4v8Z6NfU6JDnp9YdfBc7EErs5ml5vCQReT3JRrPGk96pz385HJUaSZEAxCU9vO5SQ683RGrAmfBFXkbr4InN6pLqnqxzuqNTT3oVocbpL3aoHP6oc+2quIVUCyvZmkbdFrf1SXWPxooMtSnobwUNRo/O6OpF+miBz+9FjO9p5WSKjU+q74Cy1jC0FPSKgt8c06o5JJVCnj3/UszjC0PWWtpZZiqNTpBSnut4qzKU5qWso6JSn396j3IVNAFLc9E4FKV2hVQlAXQmP1S3i/gsxIS93vS3lMekPbJnZaaxHmySDrKY5pjn0QPHgqFjdA5MCWUQqoY2SJTamspL3lRAOWesPanvduk1J9bflOyRX7C9oJ1Kp2EYdwJWJ+IMdUsF1jEq1z9N1bDs2IMboKbG6QT1SSXnaB8ckym0jn9Szi61sIAgfHituGfl1WSky0mPrWbE4o6clb6V2KnFclmn2qN40Tq620ry4eXaCSmDD1TtCm1qe3phxPeVY4g4mWuj6/Ncihg6nL48Vpp4Codvbus+y5HbwvF3AwIJjddXC4yR6+a/kvNUcG9h0M7rcapywZlbhuvV0scGiwvC1UHB4BcF4J9etaCY0ldbh2LrgDNpCmD2lEUtI9i2MoUtgPBedwGJcYJA81024jS47lmxXXw9CmCDAA3TMTgaVTSD3rk/rBGpkdCjHEWt3AhMVze3tXDYXBVBU/2lZuRmQiQZi87L8OxNYAktudvcv3fjWCwGPaP1yl6SxDXSQQdQR1X4lxzhVKnXqNomoKQeWsbUMvDQdzvuuPetw3s3xipha3pW2P7Q2cOvNfrHZPjbMYyq40y19LLmj5pnvvK/H6lJgBFpMAX5XX6l8lOCb+ouOcGoXy8DYAeqDKzxWrHosTSBMjQ8/jRYK1KbQuzUw9tDqOXsWTGU40IXZycKowxHvSSwgGAulXjcLBjGB1pLea3ozPqaXlA5yJ1MAQFmFN53ViQ0HmilCGGLnwKIUfG8KWpUzHnoqnRNbhndw70L6UQJGi1qFPSnFMqpFQ+PP4Kis7zqkVahTH3NtFmrmbIsKquPckVHSjrk7rM47puNI8pTn3Ctzo96RVfJU8kVUfv1jwSPSK3lLOi1qF138kg5vD2+KZUdHekmr171ki9LlC4/G6VUqckl7xurqmVHxuge6Uqo4aIc6aGOJS3E66qpHPvS3PhSM4t5CCRoUmrWSjVVaaHnXuSXlUahMJT6isi6soXvGiU55idEp7lUE93dZIqv5R1VvJWZxnXT45Ii31PYkvfKqoRtslOcFFfq1Oo0SSfanUsUzmF5unWE9/NdDDVhAnKFrzcPCR26WLYeSp3EG6Na5x25SuX6Vk8+7609uOA9UNF91J1pI6mHLnXn1h7E12Gb+0665dDE5raRumVg6JLxCarqYbB0xBHiui0Mju1K81h8aAYkuXRpVpFrdTurKjq18bTpg8tZ1/JZH8aaPmG/gudVoPqTdsRsVzcVg6rf2fV5j4us9X+NzmPQ0OM1HHUW16jonVOOGwDR3ryDKz2bSOadTqO1IsseVXMe2wPEWvEECd/wAE48VYDlAsB7V5TCvj5wIPVdTDYmmxvpCP+ULULXpcPjC9u4OpSW8VsQDJFgvM4zjkyGWvdZ28UiPi6WmvV1uL1dBa3kkNx9YEF89y5LOI1HQWAIcZUxFQh5IETACzqya9SHV8TAmAId864HgvL9oOGsmpVZXa91N3rsNtYi+6v0tZjBNQAusR03WGrxCmwPE3gAHLNwZ8dVxvW+nSRw61MzmiYOk6L2PyZcafhqzgSPR1WHMDMAg2K8XX40KjyxzZcNCBAudOq/RODUMP+rtgsDwyTNrmySe1tx6B3aB7qm7hsAfVi90dXiJf850ToBqvN/qhu1j2l3RC3h1YAOdUPcFrXPXXHEpdlM96aagIHVcOnRdHrG+ZbqRDRC3IH1agBWPE4kyQPNBicWxupnmudXxwJtYDRWqaeIVGOuQRP1rfT4hHzrXmy8ljcYS+Y7uSx4jiLySZPjos7SY91W4k0AnM2B1XPp8blxGWbrydKo9+uaI5rs4Wg0gRc8uSsS49AcWCAUl+IJ5LM14yxyWfF4sCABtCUns2tiWibwsVXEysD3OJ13uEiu6pKkp9Oi+oDzSS6VlDzEHdCKnIpapr3JL952UNQG0yfwS3LUiBzSluMH2IH1iNRPXkkOrHuGqIa+8Trus9YCbJb6qF1VrkVb4SnDqpUcluPNDVOeOuiW98hLrk67IRboqYt1UBJ9J5KnpDyQq0Kd5UzhZDVOhS3udoNVKmtdSqLjNdJdiG85WSox6zvadymrJrpOqgwhdfRZKdhPJMOLLbQO9XV8TXU3cj5JdSi7WCp/ST9pSn8QdzM7ptTxCaXSEt9M6FPbjXDZvegfiyT6zbKa1zHpqVED9qStuGpwATGi/P8L8pPDXfOZiGcjlBnvutzPlI4XGuI2/8v8UzXDHuWDW8WhaWOaAIbnMarw9Dt7wtxPr1W6Xcy2nQpuN7d4BlImi81Khs0EZW3GqXIc869o1xJ1yjcWspVDZ+efBfnPY7tNh2Gu/GYkhzneo0BzmwT0C9XT7YcKGmLphwscwdrz0Uns9z9PR4drQJvO2ZU99Q6u7vwXFo9suGkwMfSJjkfsW/D8YwNWIxlFx2l4aPbotWek92/TeypUGjiegRZ6p1JA5FZa3FcHTicZhgCY/2zTfuWyi+lUgtrUng6FtRuUjz1TP9XYJtIshxgjUnVOPFXAZWhnflSyWRDa1MRr/eMMebrrj9qsacBhX4mmadaoMoYwva7MSRJIaeqhz7ejp497gMzGgAwocXTzAvDTFoCxNxAfhm1akA+hFWowFohxEwL6ryeE7aYdutCqCXaGDN48U1cr3YxlHem0J9Grh3bALFQwwqMZUDHEVG5riIkSJ5plagWhmUBs6rUY/8bMbxClQa2W5pMNhXQx9GpAIPncLhcY4Y6vSs+HtktDzvtcafgvLMFTC4hgxtU0mBv+0Ev9JyIIs1c/ltl+vTrxZn+v0muaDoBcQb2m4BXkO2FXKCaNstpOpXl+M9qjTql1CsaoiBIseS8/xDtbXrEB5AGtvtXCx2kes7HPL8R/eGd4i0c5XtncaoU3kawL30uvzHhfFWNaajXBrosZg22HVdHFU6tXCDG1c7KecNinq6T+0SFedh1P6/XOFcQo1WtcCLibouI8Wy2YQSd9l4bslxMYhj/R0yynSIptn1s1pOmi7n6y9wyCmYzROX5vcusrh9enQw1fEEkuIid10v1oRLjfKb7HuC8zUp1s4DS8A7qVX1BkDgT36DuUnpHbq4ymCZN4GnclVqtOJmZHKfyXLdUklpY7vi2nvS2PGax0tBELerrVUqUDFr7g/YlVY1DWQlvretdovvCyYwveBDh5rMXR08dlI0IldGljCI07xvdceiwMjMGX67rfnZESPBJcMHjOJANN4KxUar6nrAyhfhWGZBKugMgswgJR0KR3IEqqlQHYLKJy5vYkOed1nFBjHvOllle95FzKLFVXfs5Y3vdL9O3SIG5VFNrPZePyRHFk9xQfrDNjOyz4qvlHqid0XW3PI96zVxrdc6pi6k6EA7bJf6zUOrZC1amNjnNA1SfS2MWukmoQfWaEFQyZHkNJUVRrOzEXSnYx4dG19U62u6zYmlPx0VxBjGaEj8kwV2lc51C9yUbWAaErRrVUd1Sw/mlF3wUp9TZMDatUNSH1t90D6kyEh5AuUwNfXnVA+Dus9QbrNVL4VWc62WG6X6RpWL9Ye2xCH0w6A6lFmtbhyKU9u+aEDX5roKoO35I1pjK0amU2nXadVzyHjqlOquGyYPywflfRS/M+aHMo0rmz6MNV3M+ZRCqevmUkqwixoZi3iwe4DlmPkp+sO/eJJ3JJlZ4VtClKc2u7SYtCJtdw/aPeCUgjuUgKGGuqzqZ77p1HiFRggVKgGwa9wHsKyhSFRtbxF+73nnL3377ptDjNZhJbVeJEXcXC3R265oVQh09Ng+1+KY3KHAt3DnPM9NdE6l2yeILqDXls5XSRE72XkwFZKzkZnVn09iz5Q+Kt+bjsQwQGtZLXBrRoBIsulg/lX43TFsWH8hUpsOncF+dCeimYrftH6i75auNZSM2FBIgO9FJHW9iuZxb5VeLYmm+nUfRh9PIS2kJ/5gSLOXgxNvrCKPgKE+9dap2jxbgGl5I2gAeCbV49mblyFpA1nUrhEXnl7lSOj3XYftRg8M4vxeGfWfY0yauWkzb1mn5116UfKPSrTh31nswoe6pBFi4izWhv7K/IQia889t7+HcnUi/f2/buxfyrcPwFF9GpTxFUvrF+emGhuUnqJXdq/Lpw7MQ3CYstB2cwT7F85Enr5+5FmPNJ6c7Jr6IZ8ufDyb4bGdIcw/UmM+WrhryM9HGNA3LWmPIL50k81XpHaSfjqnoyPpF/yw8JIcc+IMCQz0REnvCw/10cMifQ4qeWVp8R0Xz96R3ND6V32Jp4x9ED5ZuFer/d4rr6gtZQfLDwiwyYoTqTRB9y+eS92ioPg9eYRJzP8AX0K/5V+CuqMzDEmZ9c0sobGkjkumz5ReDPeKYxbQSJDiwtZ3ElfNDX79/t6KjVJ85+Oia1OI+oW9tOFEw3H0BGvr69eiOp214bEjH4Yc/XXy4ag5ewKi4HYeX2q7UvMj6gPbDh5H/fsOZv8A7QSl1O1fDx87G4c8vXt7F8yZhMwFWYKGPpX/ANrOHuOUYygTE6i3ikDtNw92mMw5B2LwL6X8l835hy9uiheOXsT2ZX0ZU49gAf8AvVBt9qgI6X2TKnF8LocTRmJj0rfC6+bweg6z7IV5vjXpa6pj6EPGMJJJxVG3/wA1qB3HsJtiqHT+9bovn22qou6dPD7U+1kr9+rcfwoucVh/F/2JFXtNgw2RiqEdHjXu5L8HzfBUHhcIuP27/wBr8FnLP1ijMTmkx5qP7UYMf71R6+svxE/HNTMOS15RMftY7T4Mm2Kof5lmr9q8EyxxLCel/cvx23L2BQu6DwU8kx+vDtZgjpiWaX19iQ3tdgnEg1wOsGF+UFyolPJbzj9dZ2gwp0xFKOroPtVDj+GJy+mpTt64+tfkR+LR7lPKVfJMfsTcewj1XMN4kOB8kD8YO/u/NfkVOoRYEgdCUfpXbPf/AJinksj9Tr4hp1t8d6zvqNjv3+CvzP0jt3P8yqdWd+87p6xTyi7X6cMUG8lBjuUGy/MRVf8Avu/zFF+s1B/5jv8AMU8k1+j1MeRsEv8ApDp7Cvzs4l/77/8AMVbcXVH/AJjv8yu/w2kR0VJipc5TFABWoVbQiorVFSVBFUqxdSyCvNQeKuFC1XRYcrzIS3REWaIIorjZVTaSYAkxdSCiVSlQHcR+aIBXQIHVFKp0aKNaoLNlWZWWoQ1BQRtIQyraqCcoVaoqLilJPJQORIgSVMxVOCtoQQOKkgqnFVTRVlSVVRQFEHBVOlVdU8qxKjlZ0QTzTGocgCJypRFW0jMDEgESOfRU+CSQIBNhyChCiaqFRWokRSiitKBJUVFXy70FqFUdVaaBIVAoxuhhBahKofWpCCAqw5VlVQqCN0LkaEhQCrVwrcgBRFCmUoGZVWVENFAFnQIaihW0SrhFAWqsiMqJpgWhVCMhQN7kAwoNR3ooVhv2prUgOaJqtzbK2676clC+1dyFwOs6Jnn5KvjRNYLyTeZRZSiI+IUCaBLQryqyPiFUXRftcIcqIHv8lRCF5wuFaYAqcNLexCQF1cFGR3+SkIUvIpCMKwmr4lqJh+LKR8QmplLIUATMqgHxCGUoqwETh3qAdPYqsgYUREKDxRMA4KwEeX4hWWxspq/9Kyq4R5fj4CpgM/bKqUEKJhG6gHRDAQhIRxqry9/kiEnZWUzJ0Pkpk6HyV0Kgqwm5e/yVZOhTQLhCpMy96rJ0KmgFE4N6exXl6exNXCGhQpjx3+SHJ8QrqBKFNy9FeUdfJNNJARBMySiDensU0JVOCa9vRU1vTZXVwuFaYW9EOU/BT7Cc5UznmhUWsc9FnPNT0h5oVExdH6Q8/cpnPxCBRMNH6Q8/cqznmhUTE0Wc81C880KiYavMVec8yhUVNFnPMqsx5qlEF5jzUzHmqUQXmPNTMeapRAWc81WY81SiLq8x5q8x5oVENFnPMqsxVKIi8x5qZjzVKIurzHmpmPNUoiavMeamY81SiGizHmqzHmqUQ1eY81JVKIavMeavMeZQqICznmfNTOeZ81SpFEXHmVA48yhURNFnPM+amY8yhUQ0WY8ypmPMoVENFmPMqZjzKFRDRZjzKmY8yhUQ1eY81Mx5qlENXmPNSTzVKILk81Mx5lUogvMeamY81SiAg88yoXnmUKiGrzHmpKpRDUUUUQRRRRBFFFEEUUUQRRRRBFFFEEUUUQRRRRBFFFEEUUUQRRRRBFFFEEUUUQRRRRBFFFEEUUUQRRRRBFFFEEUUUQRRRRBFFFEEUUUQRRRRBFFFEEUUUQRRRRBFFFEEUUUQRRRRBFFFEH//2Q==",
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"400\"\n",
       "            height=\"300\"\n",
       "            src=\"https://www.youtube.com/embed/PK_yguLapgA\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.YouTubeVideo at 0x7fb330414cc0>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import YouTubeVideo\n",
    "YouTubeVideo(\"PK_yguLapgA\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is a floating point error that cost Intel $475 million:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[1994 NYTimes article about Intel Pentium Error](http://www.nytimes.com/1994/11/24/business/company-news-flaw-undermines-accuracy-of-pentium-chips.html)\n",
    "![article](images/pentium_nytimes.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Resources**: See Lecture 13 of Trefethen & Bau and Chapter 5 of Greenbaum & Chartier for more on Floating Point Arithmetic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Memory Use"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sparse vs Dense"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    Above we covered how *numbers* are stored, now let's talk about how *matrices* are stored. A key way to save memory (and computation) is not to store all of your matrix. Instead, just store the non-zero elements. This is called **sparse** storage, and it is well suited to sparse matrices, that is, matrices where most elements are zero.\n",
    "\n",
    "\n",
    "\n",
    "    Here is an example of the matrix from a finite element problem, which shows up in engineering (for instance, when modeling the air-flow around a plane). In this example, the non-zero elements are black and the zero elements are white:\n",
    "\n",
    "[Source](https://commons.wikimedia.org/w/index.php?curid=2245335)\n",
    "\n",
    "    There are also special types of structured matrix, such as diagonal, tri-diagonal, hessenberg, and triangular, which each display particular patterns of sparsity, which can be leveraged to reduce memory and computation.\n",
    "\n",
    "    The opposite of a sparse matrix is a **dense** matrix, along with dense storage, which simply refers to a matrix containing mostly non-zeros, in which every element is stored explicitly.  Since sparse matrices are helpful and common, numerical linear algebra focuses on maintaining sparsity through as many operations in a computation as possible."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Speed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    Speed differences come from a number of areas, particularly:\n",
    "    - Computational complexity\n",
    "    - Vectorization\n",
    "    - Scaling to multiple cores and nodes\n",
    "    - Locality"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computational complexity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you are unfamiliar with computational complexity and $\\mathcal{O}$ notation, you can read about it [on Interview Cake](https://www.interviewcake.com/article/java/big-o-notation-time-and-space-complexity) and [practice on Codecademy](https://www.codecademy.com/courses/big-o/0/3). Algorithms are generally expressed in terms of computation complexity with respect to the number of rows and number of columns in the matrix. E.g. you may find an algorithm described as $\\mathcal{O(n^2m)}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### Vectorization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "    Modern CPUs and GPUs can apply an operation to multiple elements at once on a single core. For instance, take the exponent of 4 floats in a vector in a single step. This is called SIMD. You will not be explicitly writing SIMD code (which tends to require assembly language or special C \"intrinsics\"), but instead will use vectorized operations in libraries like numpy, which in turn rely on specially tuned vectorized low level linear algebra APIs (in particular, BLAS, and LAPACK)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Matrix Computation Packages: BLAS and LAPACK "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "[BLAS (Basic Linear Algebra Subprograms)](http://www.netlib.org/blas/): specification for low-level matrix and vector arithmetic operations. These are the standard building blocks for performing basic vector and matrix operations.  BLAS originated as a Fortran library in 1979.  Examples of BLAS libraries include: AMD Core Math Library (ACML), ATLAS, Intel Math Kernel Library (MKL), and OpenBLAS.\n",
    "\n",
    "[LAPACK](http://www.netlib.org/lapack/) is written in Fortran, provides routines for solving systems of linear equations, eigenvalue problems, and singular value problems.  Matrix factorizations (LU, Cholesky, QR, SVD, Schur).  Dense and banded matrices are handled, but not general sparse matrices.  Real and complex, single and double precision.\n",
    "\n",
    "1970s and 1980s: EISPACK (eigenvalue routines) and LINPACK (linear equations and linear least-squares routines) libraries\n",
    "\n",
    "**LAPACK original goal**: make LINAPCK and EISPACK run efficiently on shared-memory vector and parallel processors and exploit cache on modern cache-based architectures (initially released in 1992).  EISPACK and LINPACK ignore multi-layered memory hierarchies and spend too much time moving data around.\n",
    "\n",
    "LAPACK uses highly optimized block operations implementations (which much be implemented on each machine) LAPACK written so as much of the computation as possible is performed by BLAS."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### Locality"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "    Using slower ways to access data (e.g. over the internet) can be up to a billion times slower than faster ways (e.g. from a register). But there's much less fast storage than slow storage. So once we have data in fast storage, we want to do any computation required at that time, rather than having to load it multiple times each time we need it. In addition, for most types of storage its much faster to access data items that are stored next to each other, so we should try to always use any data stored nearby that we know we'll need soon. These two issues are known as locality."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Speed of different types of memory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Here are some *numbers everyone should know* (from the legendary [Jeff Dean](http://static.googleusercontent.com/media/research.google.com/en/us/people/jeff/stanford-295-talk.pdf)):\n",
    "- L1 cache reference 0.5 ns\n",
    "- L2 cache reference 7 ns\n",
    "- Main memory reference/RAM 100 ns\n",
    "- Send 2K bytes over 1 Gbps network 20,000 ns\n",
    "- Read 1 MB sequentially from memory 250,000 ns\n",
    "- Round trip within same datacenter 500,000 ns\n",
    "- Disk seek 10,000,000 ns\n",
    "- Read 1 MB sequentially from network 10,000,000 ns\n",
    "- Read 1 MB sequentially from disk 30,000,000 ns\n",
    "- Send packet CA->Netherlands->CA 150,000,000 ns\n",
    "\n",
    "And here is an updated, interactive [version](https://people.eecs.berkeley.edu/~rcs/research/interactive_latency.html), which includes a timeline of how these numbers have changed.\n",
    "\n",
    "**Key take-away**: Each successive memory type is (at least) an order of magnitude worse than the one before it.  Disk seeks are **very slow**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "This video has a great example of showing several ways you could compute the blur of a photo, with various trade-offs. Don't worry about the C code that appears, just focus on the red and green moving pictures of matrix computation.\n",
    "\n",
    "Although the video is about a new language called Halide, it is a good illustration the issues it raises are universal.  Watch minutes 1-13:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "image/jpeg": "/9j/4AAQSkZJRgABAQAAAQABAAD/2wCEABALDA4MChAODQ4SERATGCgaGBYWGDEjJR0oOjM9PDkzODdASFxOQERXRTc4UG1RV19iZ2hnPk1xeXBkeFxlZ2MBERISGBUYLxoaL2NCOEJjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY//AABEIAWgB4AMBIgACEQEDEQH/xAAbAAEAAwEBAQEAAAAAAAAAAAAAAwQFAgEGB//EADsQAAEDAgMECAYCAgICAwEBAAEAAgMEERIhMQUiMlETFBVBQ1JhkSMzQmJxoSSBU7EGNMHRFpLwVLL/xAAaAQEAAgMBAAAAAAAAAAAAAAAAAQQCAwUG/8QAMBEBAAIBAQYFBAEDBQAAAAAAAAECAxEEEhMUUVIhMTKRoQUVIjNBcYHwIyRCYfH/2gAMAwEAAhEDEQA/APz9ERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREH1lJs6hfSwufTxEmIEkvOql7M2fZ/8aHIm2+VPQYup09sfyG8lPv2fx6u8q5NsltZ8fl52+W+9P5T7qPZmz7N/jQ5kX3yvezNn4v+tDw+c81e37M49W8uSb+P6+H7eajiW6/LDjX7p91DszZ9nfxocibb5TszZ9m/xocyL75V7ftJx6u5Jv2Zx6t5ck4luvyca/dPuzKvZ1Cymmc2niBERIIedVL2Zs+zP40OZF98qxXYuqT8fyXclZwu6Nh39W8lM5LaebKc192Pyn3/AKM7szZ+L/rQ8PnPNOzKCz/40ORNt8q/v4/r4ft5rzfwyceruSjiW6/LHjX7p91HszZ9m/xocyL75XvZmz8X/Wh4fOeau79mcerfLyXu/j+vh+3mnEt1+TjX7p91DszZ9n/xocibb5Xy1NSGpL7Pa3DzX2+/hk49XeVfJ7KOEVL87NAJtyVjDe27adejs/SP9W8xedY/9cdluz+NHldedluy+NHnZWu04N7ekzJ+lO04LN3pMiPpC2b2Xo7/AA9m6/Kr2W69umj0unZbs/jR5XVvtOC/FJpbhHNO04N7ekzJ+kKd7L0OHs3X5VOy3ZfGjzsoKqlNNhu9rsXLuWj2nBZu9JkR3KntGqjqejwFxwg3uFlS2Sbfk1ZqYIpM0nxb9Js6hfSwufTxEmIEkvOql7MoLO/jQ5E23yp6DF1Ontj+Q3kp9+z+PV3lVG2S2s+Py8ffLfen8p91HszZ9m/xocyL75XvZmz8X/Wh4fOeau79mcereXJe7+P6+H7eax4luvyw41+6fdQ7M2fZ38aHIm2+U7M2fZv8aHMi++VoBry1/Fq7kvHB7cAOPVvJTxLdfk41+6fdl1mzqFlLM5lPEHNiJBDzqsHoIst0d3evqq/F1Oo4/kO5L5nuH9K7s1pmJ1dXYLTas6zqj6CK/CNOadDFnujv71J36dyXydl3lWnQR9BFluju706CK/CNOakvw5clNLA+JrHPtZ7bhRrEeDKKzMTMfwq9BFY7o7+9Q0tKanFZ7W4eatnQ/kpsa95rX0GixyWmtdYbMFIvkiso+y3Z/GjyunZbsvjR52WtvWfxau5JvWZxat5KpxrulymLoyey3Xt00el172W7P40eV1rb2L6uH05pvWdxau5Jx7nKYujI7Ldl8aPOygqqU02G72uxcu5bu9ZnFq3ks3bN7w3vodVsx5bWtpLTn2fHTHNohs0mzqF9LC59PEXGIEkvOql7M2fZ/wDGhyJtvlT0GLqdPbH8hvJT79n8eruSp2vbWfF4u+W+9P5T7qPZlBZv8aHMi++V72Zs/F/1oeHznmr2/ZnHq3y8k38f18P281HEt1+WHGv3T7qHZmz7O/jQ5E23ynZmz7N/jQ5kX3yr2/hk49Xck37M49W8uScS3X5ONfun3ZlXs6hZTTOZTxAiIkEPOql7M2fZn8aHMi++VYrsXVJ74vku5Kbfwx8ereXJTxLaefyznLfdj8p9/wCil2Zs/F/1oeG/Gea87M2fZ/8AGhyJtvlX9/H9fD9vNeb9pOPV3JRxLdfljxr90+6j2Zs+zf40OZF98r3sygxf9aHhvxnmru/ZnHq3lyXu/j+vh+3mnEt1+Uca/dPuodmUGF/8aHIm2+V8tTUpqS+z2tw819vv4ZOPV3l5L5TY98U9r92n5VjDe27adXZ+kf6uSa3nX/JRdluz+NHldOy3ZfGjzstbes/i1dyTeszi1HJTxrvScpi6Mnst17dNHpdOy3Z/GjyutbexfVw+nNN7f4tXck41zlMXRk9luy+NHnZQVVKabDd7XYuS3d6zOLVvJZu2b3hvfQ6rZjy2tbSWnPs+OmObRD6OhDepwZD5LfDJU9m2fkNXeGVBQlvVIM2/Jb4hU+Jtn5t1d4hVG3nLxF/VJZtmZN1b4ZSzcWjeH/EeaYm2Zm3VviFMTcWreH/IeahgWbZ+Q1d4ZSzbMyGrfDKYm2fm3V3iFMTbMzbq3xCoPFBWgdUnyHyXeGVFJtEvpuiEbRiw5hpUtcW9Unzb8l3iFZlxZuY1H1rbWNY8Xofo2z4ssWtkjWYmNHeN1+I8PlPNMbrO3jqfpK4uMWo4fP6pcWdmNT9ay0el3K9Id43WbvHUfSV50jsXEeHynmubizcxqPrS4xajh8/qmhuV6Q6xus7fOp+kr5wOc24DiAdbHVfQ3FnZjU/WvnDqrOz/AMqG2xEbugiIrSgIiICIiD62jr6JlLC11RCCIgCC3QqbtGgs7+VBmTbdXxiKrOy1mddVCdgpM66y+z7RoLN/lQZEX3U7RoMX/agth8vqvjETla9Ufb6d0vs+0qEB38qHU23Vjn/kcht/FiyssRFlXZqR5+LOmxYq+fi2JtvvmifGaaIY2Ftx3Kv08dhvjuWei3VpWvpWceOmPwrC/wBNHfjGi96ePPfGd1nosmxf6aPLfHcu62saWQiNwfZtj6LNRTE6RMGqx1t/lCUtUabHZjXYuarosZiJjSWVbTWdYX+1HZ/BjzunajsvgsysqCLDhU6NvM5erQ7Ude/Qx6WUVRXvmjwhjWb17t1VRFMY6x46InPkmNJl30snnd7rlznO4nE/krxFno1ay+2oQ3qcGQ+S3wyp7Ns/Iau8MqChLepwZj5LfEIU+Jtn5t1d4hXHt5y8xf1SWbZmQ1b4ZSzcejeH/GeaYm2Zm3VviHkmJuLVvD/kPNYsSzbPyGrvDKWbZmQ1b4ZTE2z826u8QpibZmbdW+IeSCGtA6pPkPku8MqazcLMhq3wzyUNa5vVJ8x8l3iFTYm2Zm3VviHksv4Zz6Y/zoWbi0bw/wCI80s2z8hq7wymJuLVvD/kPNAWkPzbq7xCoYFm2ZkNW+GUs3Ho3h/xnmmJtmZt1b4hTE3Fq3h/yHmgWbhfkNXeGV8rsnim07u66+qLm4X5t1d4hXyuySA6a9u7vsrGL0W/s7n0T9s/51aFhZ+mp+gpYWbpqPoKXFnZjU/WlxZumo+tYvWvbDF3cPkPNLCztNT9BXJe1rhdzRu97/VedJHZ2+zU+Ip0lGsOrCzNNR9BWdti14bW0P02V/pI7N32ajxFn7Wc13RYSDkdHXW3FE78K+0zHCl9PQY+p0/H8lvJT79n8ervKq1CB1ODIfJb4ZU9hZ+Q1d4ZVO3nLwN/VLrfszj1b5eS938f18P281xYWZkNW+GeS9sMeg4f8Z5qGD3fs/j1d5U37M49W+XkubCz8hq7wylhZmQ1b4ZQR12PqlRx/JdyWXvWZxat5LSrgOqT5D5LvDKy7CzchqPoW2nk9T9A9F/7Ot7F9XD6c03rP4tTyXlhi7uHyHmlhZ2mp+grJ6F7vWZxajkm9i+rh9Oa5sLN01H0FLC/dw+Q80HRxWdxau5L5o6r6Sws7TU/QV82dVZ2f+XP27/j/cREVpzhERAREQEREBERAREQEREBERAREQEREBERAREQEREH01Ltykip4mOklBbEGmzRqpf/AJBRWd8SbMn6AvlEVadmpKlOw4pnV9Z/8gorN+JNkR9AXT/+QUHQvLJJelwENuzvXyKJyuNHIYWr/wDIq/Pejzvfd5r3/wCRV+W9Hlb6VkotvCx9G/lsPbDUk29Wyscx7mYXtwmze5fSiqYWM+I/Kx0C+GX1TOBv4W7DseLNrvR5NeTZMVtPDT+jQ60zFfpH8NtBzXraqMYrvfmT3BUEW/7Xg/7auQxJ+tzZb2lrZJ1ua98XdbRQIrXKYOyPZv5fF2wmdWTBjt4Z3Oi+UgqpaYu6MgYtbhfSP4Hfgr5Q6qptOHHTSK101bsVYxTrjjRc7Uqc825+idqVOWbcvtVNFU4dOizx8vdKepq5arB0pBwaWCgRFlEREaQ12tNp1kREUsX2lFLEKSAGSO/Qt8Wym6aKz/iR6u8ZfKdGzyj2To2eUeypzssTOurmz9PiZ13n1nTRWZ8SPVvjJ00WL5kfD/m9V8n0bPKPZOjZ5R7JykdUfbo7n1fTRWf8SPV3jJ00VmfEj1b4y+U6NnlHsnRs8o9k5SOp9ujufS1ssRpJ7SR/JdpLdZvSMs3ebqPEWZ0bPKPZOjZ5R7LKNmiI01dPYf8AZxaI8dWl0jL8TeH/ACeqdIze3m6nxFm9Gzyj2To2eUeynl46uhz09rT6Rlm7zdR4i86Rl+JvD/k9Vm9Gzyj2To2eUeycvHU56e1pGRlnbzdT4i+eOqv9Gzyj2To2eUey248e4r58/F08NNFBFf6NnlHsnRs8o9lsV1BFf6NnlHsnRs8o9kFBFf6NnlHsnRs8o9kFBFf6NnlHsnRs8o9kFBFf6NnlHsnRs8o9kFBFf6NnlHsnRs8o9kFBFf6NnlHsnRs8o9kFBFf6NnlHsnRs8o9kFBFf6NnlHsnRs8o9kFBFf6NnlHsnRs8o9kFBFf6NnlHsnRs8o9kFBFf6NnlHsnRs8o9kFBFf6NnlHsnRs8o9kFBFf6NnlHsnRs8o9kFBFf6NnlHsnRs8o9kFBFf6NnlHsnRs8o9kFBfTMniwN+KzTzLI6NnlHsnRs8o9luxZuHr4GjZ6eH/Kz/7J08P+Vn/2WN0bPKPZOjZ5R7Lfzk9EaNnp4f8AKz/7J08P+Vn/ANljdGzyj2To2eUeyc5PQ0bD54sDvis08y+YOqv9Gzyj2To2eUey0Zc3E08E6KCK/wBGzyj2To2eUey0igiv9Gzyj2To2eUeyCgiv9Gzyj2To2eUeyDsC5A5rSn2FVwdK17d9j2saB9dzbL+1mA2NwtB22q5z3PMrbusfltsCDe4FsjfO6JQnZ9WKZ1SYSIW6uJHp/7CjdSzNpm1JjIhc7CH9xPJSVO0amqwCdzXNYSWtwAAXAGg/AXo2lUiJkQLOiY7E1nRggHnmPVB7DsutnDDFA52NuJuYzF7f7Up2FtMMxmkfhte9wuIdrVkAYI5GjAMIJY0mwNwL27jpyUjduVwpXwdIDidixFou3IjLlqh4OBsWvMBl6HIPEeEkXJJI/2CFCaCqFYKQwnpzoznlf8A0p5Nt18j8bpWl1wR8NuRBJBGWRuT7qHtCfrHTjo2y+YMA7rIJextoa9WdbncW/36Lp+xaqGPpagNij6PpC4nTO1rc7riHa9bBE2KOVuBotYsabjUXuM7HRO1aksLJcErHNLSHt1BN8yM8jmh4I6ihlhaJAWyRFgkDwdQTb/eSqq5LtGSWN0ZiiDDGI2tDeAA3y9b/wC1TQEREQIiICIiAiIg9AuQOa0p9hVcHSte3fY9rGgfXc2y/tZgNjcLQdtquc9zzK27rH5bbAg3uBbI3zuiUJ2fVimdUmEiFuriR6f+wo3UszaZtSYyIXOwh/cTyUlTtGpqsAnc1zWElrcAAFwBoPwF6NpVIiZECzomOxNZ0YIB55j1Qew7LrZwwxQOdjbibmMxe3+1KdhbTDMZpH4bXvcLiHa1ZAGCORowDCCWNJsDcC9u46clI3blcKV8HSA4nYsRaLtyIy5aoeDgbFrzAZehyDxHhJFySSP9ghQmgqhV9VMLumtfD6Wvf2U8m26+R+N0rS64I+G3IgkgjLI3J91G3aVQ2oNQ3oxKb3cIxpa1raWsgm7ErGtjfMGxMfiJc53CG5m/skuxqiKF0ryMAiErT5hu5eh3gue2asyB8pjlsXXD2DMOFiMu5cu2vWOYGGRpYGdHhwC1svT7Rn6IeDnsquEcshp3BkV8ZNsrXv8A6KhZSzvgM7YyY23u7uFrX/2PdT1W1aysjDKiQPZiDrYQM8+X5PuvWbSfHS9AyJjWukbI+2jsOgsgrU9PLVTCGBhkkdo0alTybOlbFHIxzXhxc0gGxa5ou4H+lHHWzQ1jqqAiKUkkFo4b8uSmdtWZxZeOGzcZIwZOc4WLiOf/AKQUUREQIiICIiD0C5A5rSn2FVwdK17d9j2saB9dzbL+1mA2NwtB22q5z3PMrbusfltsCDe4FsjfO6JQnZ9WKZ1SYSIW6uJHp/7CjdSzNpm1JjIhc7CH9xPJSVO0amqwCdzXNYSWtwAAXAGg/AXo2lUiJkQLOiY7E1nRggHnmPVB1T7JrKnD0cQs+MyNxOAuAvex6+1+ruADcWZGi5ZtSrZHFG17MEQcGtMbSN4WN8s8uaT7UrKiPo5JbtsBYNA0QcTUFVBUx00sLmzSWwt53Nv9qZ2xdot1pne4/PNcP2lUPqW1LujMzSHB+AAgg3v7r2Da9bBF0ccoDLk2LAb31Gmnogk7ErGtjfM1sTHYruceEN1v7KrPSuiginD2viluAR3EagqwNsVeMOkMcti7J7BYhwsRl3eiifXSOiZE1kbI2NLQ0Nvrqc75+qCRmx9oPcWtpnFwtcAi+a7dsLabGlzqR4AF9R/7XjNt18by9kzQ82u7o23JHfpr3XUp/wCQ7RMLGdK27CTiwC5BsbaaXCHgiOxK/BG7oL9JewuLi3NQQ0FVPUvp44XOljNnN5Z2/wBqY7ZrjJjdI0uuTcxtOosRppbuXNPtSppZHSQFjHu4iGDPO+n5CCR+xqmEF05ZGwRtkve98Wg/K9n2JVQPkjcAZWPDQwfWCSAQfyCohtSo6NzJBHI10Yj3mDQaG47wun7ZrnzumdK3G4tN8DcrG4tllnf3Q8EUuzquGm6xJC5sWW8SO8Aj/YUTqaZkAndGRE61ndxve3/+T7KWp2hUVcjX1DmvwZAYQAMgNB+ApTtR+CmjbFGI4MRDDmC52pQVqeknqsfQRl/RtxOt3DmpKiglhaJAWyRFgkD2nUE2/wB5LmkraiiLzTvwdIMLt0G45Z9ykl2lJLG6MxRBhjEbWhvAAb5et/8AaCmiIiBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQV+tfZ+0619n7X1NFBAaSAmKMkwt8K6n6vBZ/wo9XeCqk7VETpo50/UIidN18d1r7P2vetfZ+19h1eCzPhR6t8FOrQYvlR8P8Ah9VHNx0R9xr2vj+tfZ+0619n7X2HVoLP+FHq7wU6vBZnwo9W+CnNx0PuNe18d1r7P2nWvs/a+qraeAUk5EUYIhcRaKywGxR4RuN05Lfiy8SNdFrZ88ZomYjRU619n7TrX2ftXOij8jfZDEzCdxvstqwp9a+z9p1r7P2tjoorN3Gajw170UWLgZp/j9VX5iOi/wAjPcxutfZ+0619n7Wz0UVnfDZqfDWTRUjakyYnFuG2gusq5YmJlqybNasxWJ11cda+z9p1r7P2rvZUWfxX5E/SnZUVh8V+ZH0qeNQ5TL0UutfZ+0619n7XddSNpcGF7nYgdRZVFsraLRrCves0ndlY619n7TrX2ftV0UsVnrX2fteda+z9r6qip4TSQExRkmFvhXU/VoLP+FHq7wVUnaoidNHOn6hETpuvjutfZ+0619n7X2PVoLM+FHq3wU6tBi+VHw/4fVRzcdEfca9r4/rX2ftOtfZ+19h1aCz/AIUervBTq8FmfCj1b4Kc3HQ+417Xx3Wvs/ada+z9r6qtp4RSTkRRgiFx+VZYDYmYRuN05Lfiy8SNdFrZ88ZomYjRU619n7TrX2ftXOij8jfZOij8jfZbVhT619n7TrX2ftXOij8jfZOij8jfZBT619n7TrX2ftW3RR4TuN05KvQ0jarHicW4baC6i1orGss6Um87sOOtfZ+0619n7V3sqKzvivyJ+lOyorD4r8yPpWvjUb+Uy9FLrX2ftOtfZ+1d7KivbpX6X4U7Kiz+K/K/0pxqHKZeil1r7P2nWvs/au9lRWHxX5kfSqtdSNpcGFznYgdRZTXLW06Qwvs+Skb0uetfZ+151r7P2vqqKCE0kBMUZJhac4rqfq8Fn/Cj1d4K0TtUROmjjz9QiJ03Xx/Wvs/a8619n7X2PV4LM+FHq3wU6tBi+VHw/wCH1Uc3HRH3Gva+P619n7XnWvs/a+x6vBZ/wo9XeCnVoLM+FHq3wU5uOh9xr2vj+tfZ+151r7P2vqq2ngFJORFGCIXW+FZYDYo8I3G6clvxZeJGui1s+eM0TMRoqda+z9p1r7P2rnRR+RvsnRR+Rvstqwp9a+z9p1r7P2rnRR+RvsnRR+Rvsgp9a+z9p1r7P2rboo8J3G6clXoaVtTjxPLcNtBdRa0VjWWdKTed2HHWvs/ada+z9q72VFZ3xX5E/SnZUVh8V+ZH0rXxqN/KZeil1r7P2nWvs/au9lRXt0r9L8KdlRZ/Fflf6U41DlMvRS619n7TrX2ftXeyorN+K/O30qrXUjaXBhcXYgdRZTXLW06Qwvs+Skb0vq6FzepwZt+S3xCFPibZ+bdXeIVFQY+p09sfyG8lPv2fx6u8q5dvOXkr+qXOJtmZt1b4h5LqrqYmUcxbYPbESCHJv2Zx6t8vJQV+PqdRx/IdyURGswisRNofJ9q1+f8AKkz9U7Wr8v5UmXqqaLrblej0fCx9sey2/adbI0tfUvLXDCRfUK23hH4WStZvCPwsoiI8mda1r5Q9XjuE/herx3CfwpZNRgBY3IaBdWHILlnA38LpdfhY+2DiX6y5fkx1h3FY+yCLzXt3d9lsP4HfgrI2Pe81r6DRc/6hWtaeELGy2m2WNZaFxZ2Y1P1pcWbpqPrK63rP4tXcuSb1mcWo5clxXYZm17fBtyPfdZq09s3vDe+h1WYr2H0Q421ftkREW1XfbUJb1ODMfJb4hU+Jtn5t1d4hUVBj6nT2x/JbyU+/Z/Hq7yrj285eZv6pc4m2Zm3VviFMTcWreH/Iea637M49W+Xkvd/H9fD9vNYsHGJtn5t1d4hTE2zM26t8QrrfwycervKm/ZnHq3y8kFauc3qdRYt+S7xCvn28I/C+ir8fU6jj+Q7kvnW8I/Cv7J5S7H0702eoiK26IiIg8dwn8KPZBzmvbQd9lI7hP4XGxr/GtfQaLVm9ErOy/thoXFnZjU/WUu2zcxqPrXW9Z/Fq7km9ZnFqOSouw5uMWo4fP6pcWdmNT9a63sX1cPpzTes7i1dyQc4hZuY1H1rO2wR8G3I991p71mcWreSzds3vDe+h1stuH1wr7V+qX0dC5vU4M2/Jb4hU+Jtn5t1d4hUVBj6nT2x/IbyU+/Z/Hq7yqrbzl4C/qlzibZmbdW+IUxNxat4f8h5rrfszj1b5eS938f18P281iwcYm2fm3V3iFMTbMzbq3xCut+0nHq7ypv2Zx6t8vJBWrnN6nUZt+S7xCvn28I/C+ir8fU6i+P5DuS+dbwj8K/snpl2Pp3ps9REVt0RERB47hP4UeyLfGvyHfZSO4T+Fxse95rX0Gi1ZvRKzsv7YaFxZ2Y1P1pcWbmNR9a63rP4tXck3rM4tRy5Ki7Dm4xajh8/qlxZ2Y1P1rrexfVw+nNN6zuLV3JBzcWbmNR9aztsEEw2toe+6096zOLVvJZu2b3hvfQ62W3D64V9q/VL6OhA6nBkPkt8MqezbPyGrvDKgoSOpwZj5LfEKnuLP4dXeIVVt5y8Df1SWFmZDVvhlC1pdYhtsP+M80uLM4dW+IVFWPwU0zmOAcIiQRIcioiNZYxEzOjroYrP+FHq7wl70MVmfCj1b4S+Q7Vrs/wCVJn6p2rX5fypMvVWuWv1dDkMvc+oroYhRzkRsB6F1vhWWC0bo/CqP2nWvaWvqXlrhhIvqFbbwj8KzhxzSJ1Xdmw2xRMWnV7ZeOG6fwvV47hP4W5bXw42ZvHUfSV6C4uyceHynmpGAFjbjuC6GWiy+35u50edx9qI4yH5nU/SVlbI1mvyHddbLzZjvwVjbIOc1/TvstObZ74aTvTrqVzUy5a7saebSsLO01P0FLDC3TUfQlxZ2mp+sry4s3TUfWVQX2dtgD4NuR7rLNWltg/J/B77rNV7D6Icbav2yIiLarvtqEDqcGQ+S36Cp7Cz8hq7wyoKEjqcHD8lviEKe4s/NurvEK49vOXmb+qSwszIat8Mr2wx6Dh/xnmvLizOHVviFLjFq3h/ynmsWJYWfkNXeGUsLMyGrfDPJLiz+HV3iFLizM26t8QoIK4DqdRkPku8Mr59vCPwvoK4jqdRp8l3iFfPt4R+Ff2Xyl1/p/ps9REVt0RERB47hP4Uex9Zr8h3XUjuE/hR7H1m/A77LVm9ErOy/thpWFn5DU/QUsLNyGo+gry4s7TU/WlxZumo+sqi7D2wxd3D5DzSws7Ian6Clxi7uHznmlxvaan6yiCws3Iaj6Fm7YteG3I91lo3Fm6aj6ys7bBv0NuR77rbh9cK+1fql9JQgdTgyHyW+GVPYWfkNXeGVBQkdTg0+S3xCFPcWfm3V3iFVbecvA39UlhZmQ1b4Z5L2wxaDh/xnmvLizOHVviFLjFq3h/yHmsWJYWfkNXeGUsLMyGrfDKXFn8OrvEKXFmcOrfEKCCuA6nUZD5LvDK+fbwj8L6CuI6nUafJd4hXz7eEfhX9l8pdf6d6bPbJZEVt0SyWREHjhun8KPZFrzX5DuupHcJ/Cj2RrN+B32WrN6JWNl/bDRsLOyGp+gr2ws3Iaj6CvLiztNT9ZS4s3TUfWqLsFhi0HD5DzXthZ2Q1P0FeXGLu4fOeaXFnaan6ygWFm5DUfQs7bAsYbcj3WWjcWbpqPrKztsEEw25HvutuH1wr7V+qX01Bj6nT2x/Ib5VPv2fx6u8qzaOvomUsLX1MIIiAII0Km7RoLP/lQZk23VXtS2s+Dwl8V96fx+Fzfszj1b5eSgr8fU6i+P5DvKou0aCzf5UGRF91Q1lfRPpZmsqICTEQAG6lK0trHgUxX3o/H4fJIiLrPRi1m8I/CyVotqIsI3wgmXjuE/hR9Yi84XjqiIg74RLZZwN/C6VVtfShoBmbovev0v+Zq7PEr1YJ38DvwVkbHvea19BpZX319KWOAmbosWmqjTF9mNdi5rn7dpesRVv2e8UyRazdGKz+LU8k3rM4tRyWT2o/P4Med07Ufl8GPKy5HBu6fN4uqTbN/g3vodbLMViqqjU4bsa3DyVdW8dZrXSXNz3i+SbQIiLNpfcUGPqdPbH8hvlU+/Z/Hq7yrNo6+iZSwtdUQgiIAgjQqbtGgs7+TBqfpXJtS2s+Dzl8V96fx+Fzfszj1b5eSb+P6+H7eap9o0Fm/yYNRfdTtGgxf9qC2Hy+qx3LdPhjwr9vwub+GTj1d5UOOzOPVvl5Kn2jQWf8AyYMybbqdoUFm/wAmDUfSm5bocK/b8Ja/H1Oo4/kO8q+dbwj8LXrK+ifSzNbUQkmIgAN1KxG1EQaN8K9ssTETq6uwVmtZ1hKii6xF5wnWIvOFadBKii6xF5wnWIvOEEjuE/hcbGvea19BpZeGoiwnfCq0tU6mxYWNdi5rDJWbV0huwXimSLS3bOs/i1PJN6zOLVvJZPaj8/gx53TtR+XwY8rKpwbulzeLq1t7F9XD6c03rP4tXclk9qPv8mPSydqPz+DHndODc5vF1au9ZnFq3ks7bN/g3vodbKPtR+XwY8rKCqqnVOHExrcPJbMeK1bay059ox3xzWJfYUOPqdPbH8hvlU+/Z/Hq7y8lm0dfRMpYWvqIQREAQRoVN2hQWd/JgzJtuqnaltZ8Hir4r70/j8Lm/ZnHq3y8k38f18P281T7RoLN/kwai+6naNBi/wCzBw+X1WO5box4V+34XN+z+PV3lQ47M49W+XkqfaNBZ/8AJg1Nt1O0aCzf5MGovupuW6fBwr9vwlr8fU6i+P5DvKvnW8I/C16yvon0szWVEJJiIAA1Kw21EWEb4V7ZYmInV1dgrNazrCZFF1iLzhOsRecK06CVFF1iLzhOsRecIJHcJ/C42Pe81r6DSy8NRFhO+FVpao02KzGuxc1hkrNq6Q3YLxTJFpbu9Z/FqeSb1mcWo5LJ7Ufn8GPO6dqPy+DHlZVeDd0ubxdWtvYvq4fTmm9vcWp5LJ7Uff5Melk7Ufn8GPO6cG5zeLq1t6zeLUclm7ZveG99DrZR9pvy+DHlZQVVUanDdjW4eSzx4rVtrLTn2jHfHNYlZRetF3Acytqo/wCOSQukYJ2vd0jWRkDJ13Yc+Vj+Vac5iIr82yZ4KM1M0kLBewY51nO00FvUKN2z5m7PjrbsMUjsDQL3vnla3ogqItOm2FVVLYyx8LTI0PDXE3AJsL5d5KnP/Ga0RdJ0tORyDzcm17aa5IaMVFqv2DVMYMTmdI6VsbWX1Jc5t/dpUDtlVDdpGhxRmQC5cHboFr39kFFFru2G+CRgqZ2tbaQvLGkloaL5DK69m2EYqZ0/T4m9CJGANzvu3B5cQ5oaMdFpz7CqqeCWaeSCNsZLd5xBcRfIZZ8JVaKgmlon1TC0xsNnZ5g5WH93y/BQVUVihopa6rbSw4RK69g42uR3KSromwUzJmS47yOieMNrOba9uYzQU0REQIiICIiAi9aLuA5lbVR/xySF0jBO17ukayMgZOu7DnysfyiWIivzbJngozUzSQsF7BjnWc7TQW9Qo3bPmbs+OtuwxSOwNAve+eVreiCoi0YtjzSQRTdJGGyxve3UndF7HLVSP2BUxsxSTU7G4Q67nnXlpqgykV6p2TUU20YqGR0ZllLQMLrgEm2f9qy3/jtU5pcyemcwXu5ryQABmdO7RBkItY7EdE9nWJmtZaRzywYiA0Xy53VOppWR0kFTFI57JbtOJti1wtca+oQVUREQIiICIiAiIgIiICL1ou4DmVtVH/HJIXSME7Xu6RrIyBk67sOfKx/KJYiK/NsmeCjNTNJCwXsGOdZztNBb1Cjds+Zuz4627DFI7A0C9755Wt6IKiLRi2PNJBFN0kYbLG97dSd0XsctVI/YFTGzFJNTsbhDruedeWmqDKRXqnZNTTbQjonYHTSWw4XXGZsrJ2BKyWFkk8fxHlt23Nhha6/s5BkItpuw43NhlbUkwztODc3gd61xfSzT3qEbCqer9YfLBHFYHE95GoBHd9wQ0ZaK1TUE1VBLNGWYIc33OgsTf9W/Kip6d9RUsgZYPebDFlmgiRX6rZpp6V8hkvLE8MmZbhLgSLHv0N1QQEREQIiICIiArBrqtzi41U5cRYkyG5HJQNF3Acytqf8A47JC6Rgna93SNZGQMjd2HPlY/lEszrtV0zZXTve9puC84rH+106vrJZMRqJcV7jC62foApZtkzwUZqZpIWC9msc4hztNBb7guWUVTDSMr43hrL2Y5riHYrnIeuV0F+i2dteWKIwVZjjczEz4xAt3/sgH1XJg2jDSTwGttkx0kYcb55AX/v8ACihpqqNsL461zC6mfM0NcQQBe7R+cKnbQbQmxupqyQxxU7HXkeW5ObfCNf8A8ESS0O043OfPtINBLWYnTO3s3ZDLuLXKeOk2vM6SOGvG7O6K7yWvJF8725DuKqVtFVPbUydcfUNpCGymQm4dplzGZsfyq1LPtGqqBBHVz9ISXAOldm4D/aDRZQ7QjqoadlcX1bi57Y3EuZYXaTc6k2OVtPZY76qrDheolyBYLPIGHkPRT1VTXU73U8tQ4ujed8OzB7xi117uYVifZ9VJQxy1FQwsZAZIWi5JbiF+7m79IhPWUFVEWiWvklmbAKlrTvN3Sbi9+63JUanadaHCN07HYHtfiYwC5GncL2uoIKiow9WEuBstmEu7m3va+oFzdW9pUFUxklTWSA1HSNY+MDMXBtfu0HcgzzUTGd83SOEjyS5wNib6ruernqWMZK/E1l7C1szqTzJ5qXZdD2hXspi/owQSXWvawuoKqA09TLCST0b3Mva17GyCJFs1OwmwYGdZcZXwGVrejsMhcgm+WSxkBERECIiArBrqtzi41U5cRYkyG5HJV0RKx12q6Zsrp3ve03BecVj/AGunV9ZLJiNRLivcYXWz9AFVXUT5I5Gvic5rwci02IQbFLQ7UnbT9HWYTgDo2mRwLA44R3ZXvZSy7D2rLBikrI5GgYQ0zE917aei4otn7Vmgp3QV2Bj2ksHSvGHeAtkLXJIXLmVUVNUwx18nR9CycsNwX3sLfgYuahKY022IosJq7vfK2INJu7MuaN45jNhVRtNtGCv7Pjqi1zd+7JTgGVyfZWJqasZJKJtpymQvjiYcTiHE5i5vkAqtbU7UoZ3081S/pGEYnB1yDbTFqpQuCirqWeMz15DR0sj7Xfawz3Tre/7XNbsuZ1IJZKkOjjgD4mtjAtwkgju4wb5qq7tGnpKet6y58cjjgBcXXOYdcH0/3+VzRNr68y9FVOHQREkGQizO8C3d6IOp9hVVPBLNNJDG2Mlu84guIvkMs+EqtFQTS0T6thaYmGzs8wcsv7vl+Cpa99VDUNa6rknDWska4kkZgEHP8rRpXVNXsqSsnrCxkMzXEMiBLrWse69sSDGpKSWsn6CEAyEEhp77C9vypqygFNAyRsvSfEdE8WtZ7bXse8Z6rqeSp2btaZ4ka6cOJ6TCDe+dxyOarTVU00ccb3bkfC0Cw9T+fVBCiIiBERAREQFYNdVucXGqnLiLEmQ3I5KuiJWOu1XTNldO972m4Lzisf7XTq+slkxGolxXuMLrZ+gCqrqJ8kcjXxOc14ORabEINilodqTtp+jrMJwB0bTI4FgccI7sr3spZdh7VlgxSVkcjQMIaZie69tPRc0WzdsSxRGCsLI3MxR/GcBbv/ZAPquDBtGCkng67bJjpIg5188gCf77rhQlMabbEUWE1d5HytiAJu7MuaN45jNhVYx7UZtU0nXHOmB6Qu6UluTb3P8AXoppaHacbnPqNpBoJazG6V+9mchl3Frv/Cmjo9rzOkZDXi7Z3RXe5weSL53tyHcUEU1DWsnjNRXGNrekeeiv8PCLnC3LX+l5U7NqRRFz6+SSFsIkjZmR9N2kXy4hpddtodoMqoadtcX1bi97Y3OLmWzaTc6k2Pdp7LHfV1YcAaiXIFgs8gYeQ9PRSNisp9rCiqDV1FOyMHC69ml9ichYZ5tKoxtrJse1g6N5jdvFwBs4WAFrW78vwVcrKCqiIbNXySzNgFS1pOJu6TcXv3W5KjUbTrQ7o3VDX4HtfiY0C5GncL2uiHMRrNpOFI1wc67pMJFi92pvzOuqmbscCoggfUjHNia0sbiaHg2sTl7qnDWzw1bqpjh07iTjIBzOp/K8hraiH5crhYEC+eG+tuRPNBARYkHuRERAiIgIiICsGuq3OLjVTFxFiTIbkclXRErHXqrpmzOnkdI04gXuxWP9ruXaVXN0eKUhzHFwLBhOI9+Xfkqi6je+ORr4nOa8aFpsQg2qbY+1JYqcxVbWscAW/FcMGIC17Dvxd3Nex7M2nmxle2zYw0tbK42FrhvsCeS7p6LbVU2CSKswN6ICIiUtAbYcvWwPqqpoq6nZORW4ZMLDIxsjrm+gP/4qErNVT1E088UdXK+aOSOMAtDRIHWsXWOunsF66gqY9pwzSVLHVT5nB78Ic1tmh2IetjyGYUUmzdoROc6baLGC7WYjK7ezNgMu4tPtkpYtn7VmqHmHaBLm1BZaSRwJc3v0sch6oIRRS1MdO2WpPVJZCKd2G5L3HvH9Z8rhZhqamMlnTyjDdtg86aW/QW6afaHWKemZUwirs90bGRtDGi5BcCBxZE6aD+liVlIaYQHH0nSxdJcd2ZH/AIUoRyVEkkrJJCHOYGtFwLWAsBb+lI6vqnVDpjM7E54eQNCRplpkpqDZpq6WqqTJgjp2gkCxLsxoPwo62jFM92GVr2WY5txZzg5uIGyCIVU7al1Q2V7ZnEuL2mxudV5NUzTsYyV5eGEkX1uTcqJEFo7RqzG5jp3ODm4LuNyG8ge4KqiICIiIEREBF60XcBzK2p/+OyQmRona93SNZGQMjd2HPlY/lEsRdRvfHIHxOcx40c02IV2bZM0FGamaSFgvZrHOIc7TQW9QuW0dTDSM2hHIGsvZrmuIdiuch65XQaVDRbZkgb0NWGxSxbpdKbYcr25Z2BVWSirI4pA+sGKONrTGHuuGu+nl/WiloNl7SljgmpqxsYe04D0jhh3gLZDIkkLmSGoYJTHXzWhpWPaC4glrsNxroLqEpnUW0oWvE1fEGEtbie8m5BIGHK4ILXJ2XX1cb6ZlVG6NtS5mGRxxFzQRfTQNHcclHU0M7nyNqq973iRkcZJLgXEXFz3AC/uq1bLtKhqX081ZNjjcC7DK4gOt/tSLsOzq2Cogpoau9WcTmR3uwNGIXB58R00us2eKWgEPRTu+NEJCWEgakW/SsTnaLNmwVL6kuilecBDjjv3569wOvJc7ObXVYqBBUZRwuxiRxO4eK2qISbOhqaqhrZDUubDHE1rhfEXNDhkLnQapFFVU+0GUVHU3jlwyMx5MILQ4Fwz7v9KpUwvoKjBHKHgsY4uA3XAgHQ6j8rToBUTUs+1ZK6WMxyC4Yy9wABe17ZB39IOJdgV80kks1RAXG73vc46ZnFppkfZRdhvjY500zLGN7mdHvXLWlwueRA/YVevkrKWolpZKuV4aT9Zs6+d7et1AK2qAYBUzARizLPO6PTkgv9hySh76eRpijhZK50m6d5t7D11Xk2yI4pZ29YIDZWRRFzeIuF8+QsqLa2qa3C2pmAw4LB54eX4XQ2hVhrh1h7sVs3G5FtLE6W9EFqbZDomysfK1s8EjWStPC3Fob+neswixsrTNpVbGsDZjZjxILi5uNL87Ks5xc4uOZJuUHiIiIetF3Acytqf/AI6+EyNE7Xu6RrIyNHXdhz5EH8rEVg11Y55e6qnLiMJJkNyOSJTzbKmgozUzSRMbezWOdZztNBb7guW0dTBSM2hHIGsvZrmuIOK5yHrldRdeqjM2Z08j5GnEC92Kx55ruXaNZMWXlIcxxe0sAacR78u/JBfodmbSljgmpqxsYe04LSOGHeAtkMiSRkuZIahglMdfNaGlY9ocSCWuw3HKwv8A6U9DRbalgaIKzDFJFul0pthyuByzsCq0lDWxxPD6zejja0xh7smu+nl/WihKSooJ3PkbU17nvEjI4y4khziLi5Olv/KqVc20qOfoZqqcPaQ8jpTk4i+fqtA0e0qfEZ66AsuxvxXF4cQXAAAjuLXI/Zu0a8zxdaje01bgWvJBc8Xu7TkO4qRUmdtFmz4p31BMT3XY4O3iTe4vr3XI/Cip4KvazRG2RpFO3R5tZpOZ9zmtJtFXMlpKOKpjFUwPMTGtGANzBN+8mx/oLFjnlphLHGQOkGBzsOZHoUQsVbp6LFSsn6SN8TCXAatIxWBOdrlWIaLtCjqNp1dQ4FrwHBrASRkCf2FPUbPqKWGSB1a50j6QSFtrtLWk3be+VrLIp62opm4YpXBmLEWatJ9RoUCtpnUdXJA5wcWHUd41CgXcsr5pXSyuLnvN3OPeVwgIiIgREQEREBERAXUb3xPD43OY9uYc02IXKIJo6qoiLTHPKwtBAwvItfVeisqha1RLk3AN85N5fhQIiUstVUTG808khy4nk6aLvtCs/wD659cXzDrzVdEFk7Qq3RhhqHmxJBvvZ656rllbUsDgJn2cwxm5vunUZ6XUCIOo5Hxuuw2vqO4+hHeFNLXVMzZWySkiUhzxYZ209lXRARERAiIgIiICIiApzXVbnFxqpi4ixJkNyOSgREp+u1XTNlM8jpGnEC84rH+1JLtGsnLAZTdji5pYA04j35d+SqLqOR8Tw+NzmPbo5psQg3aGh2zJA0QVgbFLFul0pthyvblnYFVpKCtjieH1e9HG1pjD3GzXfTy/rRZ8dVURFpinlYWggYXkWvqvRW1QtapmybgG+cm8vwg130O0oA8TV8TWEsbie8m5BIGG4uCC13JOytoVUb6ZtUx0bakswyEglzQRfTQNHNY0tVUTG808shyO88nTT/ak7Qrf/wCufXF8w680GrFs2thqIKaGsvVnE5jL3YGi4uDz4u5Z08M1CIjFM/40XSExkgakW/S4O0at0YYah5sSQ6+9nrnquWVtSwOAmeQ5hjNzfdOoF9LoLuzoamqoa2Q1Lmwxxta5uTi5oIyF9ANV7FBVwV7KOjqT0cmGRpfkw3biBcMxp/pZkcj4nXYbX1HcRyI7wrDdp1jak1AncJCQctLgWGWmQQaEv/H6+aSSWaeAuJLnPc86ZnFppkfZRdhyRsLppmZxvczo96+FpcM8siB66hUHVtU5hY6pmLTe4LzY31QVlUAwColAjFmWed0enJBf7DllD308gMTIWSl0m6d4XsPXVeTbHZFLM3rBAbKyKIuZxFwvnyFvyqLa2qa3C2pmAw4LB54eX4XTa+qDXDp3uxW4jci2hBOlkFqbY7omysdK1s8EjGStI3W4tCD3+uXus0ixIVpm0qtjYwJjuPEguLkuGl+dlVcS5xccyTcoPEREQIiIC6ifIyRronOa8HItNiFyuo5HxPD43OY9ujmmxCJbdJsbarmwSQVTWB7QWkSuGAOtYGw77jRGbL2m42bXN3I7ENlcS0EXw+wvbRZh2lWkxkVUrejYGNwuLbD+vwFGyqqGY8E8rceTrPIxflBq1NBPNNPEKyeaaOWNrBJo7Fob3yIuVK6hr45YYn7SmxumcHlryWts0OxDPM2PosXrlSSSaiW5ABOM92nsuxtCsEmPrUxdiDs3k3PNBqxw1UzYcddIKeZ+GCW13lzicj3jMZ58llVb5iyFk0oeGtIaB9O8b3/tSna1XiYWuY3o8WDDG0BpOpGWvqq0tRLNHEyR12xAtYOQvdB6yrnZC+FshwPaGkel72v3C/coURARERAiIgIiICIiAiqdZk9PZOsyensgtoqnWZPT2TrMnp7ILaKp1mT09k6zJ6eyC2iqdZk9PZOsyensgtoqnWZPT2TrMnp7ILaKp1mT09k6zJ6eyC2iqdZk9PZOsyensgtoqnWZPT2TrMnp7ILaKp1mT09k6zJ6eyC2iqdZk9PZOsyensgtoqnWZPT2TrMnp7ILaKp1mT09k6zJ6eyC2iqdZk9PZOsyensgtoqnWZPT2TrMnp7ILaKp1mT09k6zJ6eyC2iqdZk9PZOsyensgtoqnWZPT2TrMnp7ILaKp1mT09k6zJ6eyC2iqdZk9PZOsyensgtoqnWZPT2TrMnp7ILaKp1mT09k6zJ6eyC2iqdZk9PZOsyensgtoqnWZPT2TrMnp7ILaKp1mT09k6zJ6eyC2iqdZk9PZOsyensgtoqnWZPT2TrMnp7ILaKp1mT09k6zJ6eyC2iqdZk9PZOsyensgtoqnWZPT2TrMnp7IN3YtXTUksrqoYmOZbDgDsXp6LzblXTVtcJaVmGPBa2HD3n/AMWH9LD6zJ6eydZk9PZE6oUREQIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiD/2Q==",
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"400\"\n",
       "            height=\"300\"\n",
       "            src=\"https://www.youtube.com/embed/3uiEyEKji0M\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.YouTubeVideo at 0x7f86d4fade48>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import YouTubeVideo\n",
    "YouTubeVideo(\"3uiEyEKji0M\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Locality is hard.  Potential trade-offs:\n",
    "- redundant computation to save memory bandwidth\n",
    "- sacrificing parallelism to get better reuse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Temporaries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "The issue of \"temporaries\" occurs when the result of a calculation is stored in a temporary variable in RAM, and then that variable is loaded to do another calculation on it. This is many orders of magnitude slower than simply keeping the data in cache or registers and doing all necessary computations before storing the final result in RAM. This is particularly an issue for us since numpy generally creates temporaries for every single operation or function it does. E.g. $a=b\\cdot c^2+ln(d)$ will create four temporaries (since there are four operations and functions)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### Scaling to multiple cores and nodes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "We have a separate section for scalability, but it’s worth noting that this is also important for speed - if we can't scale across all the computing resources we have, we'll be stuck with slower computation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Scalability / parallelization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Often we'll find that we have more data than we have memory to handle, or time to compute. In such a case we would like to be able to scale our algorithm across [multiple cores](http://www.makeuseof.com/tag/processor-core-makeuseof-explains-2/) (within one computer) or nodes (i.e. multiple computers on a network). We will not be tackling multi-node scaling in this course, although we will look at scaling across multiple cores (called parallelization). In general, scalable algorithms are those where the input can be broken up into smaller pieces, each of which are handled by a different core/computer, and then are put back together at the end."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.7.0 64-bit",
   "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.7.0"
  },
  "vscode": {
   "interpreter": {
    "hash": "c0f2ef85a013b3e049ebe5ad4c39b42bd34e85c4b0686b7bef6ce5ff62aa91e1"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}