{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Hamiltonian Monte Carlo\n", "\n", "Having introduced first-generation MCMC methods previously, we will now turn our attention to a more sophisticated class of algorithm, which improve upon random-walk jump algorithms by using information about the posterior distribution to inform candidate transitions: namely, gradient information.\n", "\n", "In order to implement these methods in code, we require more powerful mathematical software tools that allow for automated gradient calculation. There are several open source toolboxes that can support gradient-based Monte Carlo methods, and we will look at one of them in detail." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction to Mathematical Expressions with Theano\n", "\n", "Theano is a Python library that allows you to define, optimize, and evaluate mathematical expressions involving multi-dimensional arrays efficiently. Theano features:\n", "\n", "* __tight integration with numpy__ – Use numpy.ndarray in Theano-compiled functions.\n", "* __transparent use of a GPU__ – Perform data-intensive calculations up to 140x faster than with CPU.(float32 only)\n", "* __efficient symbolic differentiation__ – Theano does your derivatives for function with one or many inputs.\n", "* __speed and stability optimizations__ – Get the right answer for log(1+x) even when x is really tiny.\n", "* __dynamic C code generation__ – Evaluate expressions faster.\n", "* __extensive unit-testing and self-verification__ – Detect and diagnose errors.\n", "\n", "Theano is part programming language, part compiler. It is often used to build machine learning, though it is not in itself a machine learning toolkit; think of it as a **mathematical toolkit**.\n", "\n", "After a brief introduction to the Theano package, we will use it to implement a modern MCMC algorithm, *Hamiltonian Monte Carlo (HMC)*." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Adding Two Scalars\n", "\n", "To get us started with Theano and get a feel of what we're working with, \n", "let's make a simple function: add two numbers together. Here is how you do\n", "it:\n", "\n", "### Step 1 - Declaring Variables" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "from theano import function, shared\n", "import theano.tensor as tt\n", "import theano\n", "\n", "x = tt.dscalar('x')\n", "y = tt.dscalar('y')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In Theano, all symbols must be typed. In particular, tt.dscalar\n", "is the type we assign to \"0-dimensional arrays (scalar) of doubles\n", "(d)\". It is a Theano type." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "theano.tensor.var.TensorVariable" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "type(x)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "TensorType(float64, scalar)" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x.type" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "TensorType(float64, scalar)" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "tt.dscalar" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice that we have created objects of the type TensorVariable. A **tensor** is a generalization of an array to (potentially) multiple dimensions. Thus, everything from a scalar to a 5-dimensional hyper-matrix can be accomodated with the same abstraction. All expressions defined in Theano are performed by associating tensors with operations and with one another." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Step 2 - Symbolic Expressions\n", "\n", "The second step is to combine *x* and *y* into their sum *z*:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "z = x + y" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*z* is yet another *Variable* which represents the addition of\n", "*x* and *y*. You can use the pp function to *pretty-print* out the computation associated to *z*.\n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(x + y)\n" ] } ], "source": [ "from theano.printing import pp\n", "print(pp(z))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Step 3 - Compiling a Function\n", "\n", "The last step is to create a function taking *x* and *y* as inputs\n", "and giving *z* as output:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "f = function([x, y], z)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The first argument to function() is a list of Variables\n", "that will be provided as inputs to the function. The second argument\n", "is a single Variable *or* a list of Variables. For either case, the second\n", "argument is what we want to see as output when we apply the function. *f* may\n", "then be used like a normal Python function.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can call the function:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "5.0\n" ] } ], "source": [ "print(f(2, 3))" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "28.4\n" ] } ], "source": [ "print(f(16.3, 12.1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you are following along and typing into an interpreter, you may have\n", "noticed that there was a slight delay in executing the function\n", "instruction. Behind the scenes, *f* was being **compiled into C code**." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Internally, Theano builds a graph structure composed of interconnected Variable nodes, op nodes and apply nodes. \n", "\n", "An op node encapsulates a particular mathematical operation, such as an arithmetic operation or a transformation.\n", "\n", "An apply node represents the application of an op to some variables. It is important to draw the difference between the definition of a computation represented by an op and its application to some actual data which is represented by the apply node. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A Variable is the main data structure you work with when\n", "using Theano. By calling tt.dscalar with a string argument, you create a\n", "Variable representing a floating-point scalar quantity with the\n", "given name. If you provide no argument, the symbol will be unnamed. Names\n", "are not required, but they can help debugging." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Adding Two Matrices\n", "\n", "If we want to work with matrices instead of scalars, the only change\n", "from the previous example is that you need to instantiate *x* and\n", "*y* using the matrix Types:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "x = tt.dmatrix('x')\n", "y = tt.dmatrix('y')\n", "z = x + y\n", "f = function([x, y], z)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here is the expression graph corresponding to the addition of x and y:\n", "\n", "![expression graph](images/expression_graph.png)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "dmatrix is the Type for matrices of doubles. Then we can use\n", "our new function on 2D arrays:\n" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[11., 22.],\n", " [33., 44.]])" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "f([[1, 2], [3, 4]], [[10, 20], [30, 40]])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following types are available:\n", "\n", "* **byte**: bscalar, bvector, bmatrix, brow, bcol, btensor3, btensor4\n", "* **16-bit integers**: wscalar, wvector, wmatrix, wrow, wcol, wtensor3, wtensor4\n", "* **32-bit integers**: iscalar, ivector, imatrix, irow, icol, itensor3, itensor4\n", "* **64-bit integers**: lscalar, lvector, lmatrix, lrow, lcol, ltensor3, ltensor4\n", "* **float**: fscalar, fvector, fmatrix, frow, fcol, ftensor3, ftensor4\n", "* **double**: dscalar, dvector, dmatrix, drow, dcol, dtensor3, dtensor4\n", "* **complex**: cscalar, cvector, cmatrix, crow, ccol, ctensor3, ctensor4" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "An example of a slightly more interesting function is the logistic curve:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "x = tt.dmatrix('x')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The logistic transformation:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "s = 1 / (1 + tt.exp(-x))" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[0.5 0.73105858]\n", " [0.26894142 0.11920292]]\n" ] } ], "source": [ "logistic = function([x], s)\n", "print(logistic([[0, 1], [-1, -2]]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Theano supports functions with multiple outputs. For example, we can\n", "compute the elementwise difference, absolute difference, and\n", "squared difference between two matrices *a* and *b* at the same time." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "a, b = tt.dmatrices('a', 'b')\n", "diff = a - b\n", "abs_diff = abs(diff)\n", "diff_squared = diff ** 2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "When we use the function f, it returns the three computed results as a list." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[array([[ 1., 0.],\n", " [-1., -2.]]), array([[1., 0.],\n", " [1., 2.]]), array([[1., 0.],\n", " [1., 4.]])]" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "f = function([a, b], [diff, abs_diff, diff_squared])\n", "\n", "f([[1, 1], [1, 1]], [[0, 1], [2, 3]])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Setting a Default Value for an Argument\n", " \n", "Let's say you want to define a function that adds two numbers, except that if you only provide one number, the other input is assumed to be one. In Python, the default value for parameters achieves this effect.\n", "\n", "In Theano we make use of the [In](http://deeplearning.net/software/theano/library/compile/io.html#function-inputs) class, which allows you to specify properties of your function's parameters with greater detail. Here we give a default value of 1 for y by creating an In instance with its value field set to 1. Inputs with default values must follow inputs without default values (like Python's functions). There can be multiple inputs with default values. These parameters can be set positionally or by name, as in standard Python." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "from theano import In\n", "\n", "x, y, w = tt.dscalars('x', 'y', 'w')\n", "z = (x + y) * w\n", "g = function([x, In(y, value=1), In(w, value=2, name='w_by_name')], z)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "g(33) = 68.0\n" ] } ], "source": [ "print('g(33) = {}'.format(g(33)))" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "g(33, 0, 1) = 33.0\n" ] } ], "source": [ "print('g(33, 0, 1) = {}'.format(g(33, 0, 1)))" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "g(33, w_by_name=1) = 34.0\n" ] } ], "source": [ "print('g(33, w_by_name=1) = {}'.format(g(33, w_by_name=1)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Maintaining State with Shared Variables\n", "\n", "It is also possible to make a function with an internal state. For example, let’s say we want to make an accumulator: at the beginning, the state is initialized to zero. Then, on each function call, the state is incremented by the function’s argument.\n", "\n", "First let’s define the accumulator function. It adds its argument to the internal state, and returns the old state value." ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "state = shared(0)\n", "inc = tt.iscalar('inc')\n", "accumulator = function([inc], state, updates=[(state, state+inc)])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This code introduces a few new concepts. The shared function constructs so-called shared variables. \n", "\n", " state = shared(0)\n", "\n", "These are hybrid symbolic and non-symbolic variables whose value may be shared between multiple functions. Shared variables can be used in symbolic expressions but they also have an internal value that defines the value taken by this symbolic variable in all the functions that use it. It is called a shared variable because its value is shared between many functions. The value can be accessed and modified by the get_value and set_value methods.\n", "\n", "The other new thing in this code is the updates parameter of function. \n", "\n", " updates=[(state, state+inc)\n", "\n", "updates must be supplied with a list of pairs of the form (shared-variable, new expression). It can also be a dictionary whose keys are shared-variables and values are the new expressions. Here, the accumulator replaces the state‘s value with the sum of state and the increment amount inc." ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0\n" ] } ], "source": [ "print(state.get_value())" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0\n" ] } ], "source": [ "print(accumulator(1))" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1\n" ] } ], "source": [ "print(state.get_value())" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1\n" ] } ], "source": [ "print(accumulator(300))" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "301\n" ] } ], "source": [ "print(state.get_value())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is possible to reset the state. Just use the set_value method:" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [], "source": [ "state.set_value(-1)" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-1\n" ] } ], "source": [ "print(accumulator(3))" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2\n" ] } ], "source": [ "print(state.get_value())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we mentioned above, you can define more than one function to use the same shared variable. These functions can all update the value." ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [], "source": [ "decrementor = function([inc], state, updates=[(state, state-inc)])" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2\n" ] } ], "source": [ "print(decrementor(2))" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0\n" ] } ], "source": [ "print(state.get_value())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You might be wondering why the updates mechanism exists. You can always achieve a similar result by returning the new expressions, and working with them in NumPy as usual. While the updates mechanism can be a syntactic convenience, it is mainly there for *efficiency*. Updates to shared variables can sometimes be done more quickly using in-place algorithms (e.g. low-rank matrix updates). Also, Theano has more control over where and how shared variables are allocated, which is one of the important elements of getting good performance on the GPU." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise: Create and manipulate Theano objects\n", "\n", "To give you some practice with basic Theano data structures and functions, try making the operations below work by implementing the functions that are needed." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def make_vector():\n", " \"\"\"\n", " Create and return a new Theano vector.\n", " \"\"\"\n", "\n", " pass\n", "\n", "def make_matrix():\n", " \"\"\"\n", " Create and return a new Theano matrix.\n", " \"\"\"\n", "\n", " pass\n", "\n", "def elemwise_mul(a, b):\n", " \"\"\"\n", " a: A theano matrix\n", " b: A theano matrix\n", " \n", " Calcuate the elementwise product of a and b and return it\n", " \"\"\"\n", "\n", " pass\n", "\n", "def matrix_vector_mul(a, b):\n", " \"\"\"\n", " a: A theano matrix\n", " b: A theano vector\n", " \n", " Calculate the matrix-vector product of a and b and return it\n", " \"\"\"\n", "\n", " pass\n", "\n", "a = make_vector()\n", "b = make_vector()\n", "c = elemwise_mul(a, b)\n", "d = make_matrix()\n", "e = matrix_vector_mul(d, c)\n", "\n", "f = function([a, b, d], e)\n", "\n", "rng = np.random.RandomState([1, 2, 3])\n", "a_value = rng.randn(5).astype(a.dtype)\n", "b_value = rng.rand(5).astype(b.dtype)\n", "c_value = a_value * b_value\n", "d_value = rng.randn(5, 5).astype(d.dtype)\n", "expected = np.dot(d_value, c_value)\n", "\n", "actual = f(a_value, b_value, d_value)\n", "\n", "assert np.allclose(actual, expected)\n", "print(\"SUCCESS!\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: Logistic regression\n", "\n", "Here is a non-trivial example, which uses Theano to estimate the parameters of a logistic regression model using gradient information. We will use the bioassay example as a test case:" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [], "source": [ "rng = np.random\n", "\n", "dose = np.array([-0.86, -0.3 , -0.05, 0.73])\n", "deaths = np.array([0, 1, 3, 5])\n", "training_steps = 1000" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We first declare Theano symbolic variables:" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Initial model: 1.0 0.0\n" ] } ], "source": [ "x = tt.vector(\"x\")\n", "y = tt.vector(\"y\")\n", "w = theano.shared(1., name=\"w\")\n", "b = theano.shared(0., name=\"b\")\n", "\n", "print(\"Initial model:\", w.get_value(), b.get_value())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "... then construct the expression graph:" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [], "source": [ "# Probability that target = 1\n", "p_1 = 1 / (1 + tt.exp(-(x*w + b))) \n", "\n", "# The prediction threshold\n", "prediction = p_1 > 0.5 \n", "\n", "# Cross-entropy loss function\n", "xent = -y * tt.log(p_1) - (5-y) * tt.log(1-p_1) \n", "\n", "# The cost to minimize\n", "cost = xent.mean() \n", "\n", "# Compute the gradient of the cost\n", "gw, gb = tt.grad(cost, [w, b]) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compile Theano functions:" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [], "source": [ "step = theano.shared(10., name='step')\n", "train = theano.function(\n", " inputs=[x, y],\n", " outputs=[prediction, xent],\n", " updates=((w, w - step * gw), (b, b - step * gb), (step, step * 0.99)))\n", "predict = theano.function(inputs=[x], outputs=prediction)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Train model:" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Final model: 7.795152457520207 0.8536420249656105\n" ] } ], "source": [ "for i in range(training_steps):\n", " pred, err = train(dose, deaths)\n", " \n", "w, b = w.get_value(), b.get_value()\n", "\n", "print(\"Final model:\", w, b)" ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8VPW9//HXZyYbYYcEgYQQNhXEBUwR17oL9hZs6wLFalsttmpv77W3ai+9ba8trdr763avvS1a64bg0lZpS7Eu9apVhCiIgCAxISGsYU3IPjPf3x8z4BASMsBMzszk/Xw85jFnzjkz8+Zk8ubknDnnmHMOERFJLz6vA4iISPyp3EVE0pDKXUQkDancRUTSkMpdRCQNqdxFRNKQyl1EJA2p3EVE0pDKXUQkDWV49cZ5eXmuuLjYq7cXEUlJ77zzzk7nXH5n83lW7sXFxZSWlnr19iIiKcnMKmOZT5tlRETSkMpdRCQNqdxFRNKQyl1EJA2p3EVE0lCn5W5mD5vZDjNb3cF0M7NfmlmZma0ys4nxjyki0on586G4GHy+8P38+V4n8lQsa+6PAFOOMH0qMCZymw387/HHEhE5CvPnw+zZUFkJzoXvZ8/u1gXf6ffcnXOvmVnxEWaZDjzmwtfrW2pm/cxsiHNua5wyikg3Fgo5WoIhWoIhWgMhgiFHa8gRCIZoDToCoRAj7/o2WQ0Nhz6xoYGmO+/mrTMvI+QcwZAj5CDkXOQG7sBwKDzeRaY7Pn7sIo8PGyb8mMhjAIeLGubgNMeh833yxHzGF/RN4FKLz0FMBcCmqMfVkXGHlbuZzSa8dk9RUVEc3lpEkoVzjqbWEHsbW9jX2Mrehlb2NYZv9c0BGlqC7G8OUN8cYH9zgIbmII2tQZoO3kI0BYI0tgTDZR4I3wKhzq/zXL65ut3xWVs286VHlsf7n3rc+vbITIlyt3bGtfvTcM7NA+YBlJSU6MrcIiliX2Mrm3Y3UL2nke21TdTUNVNT18zO/c3U7A8P76pvoSUQOuLrZPqNntkZ9MzKIDfLT26Wn+xMP/1ys8jJ9JGT6Scnw092po8sv4+sDB+ZkfvsyHCG38j0he8z/D4yfUbzYwX02Hp4wbcWFPDHW8/B7zN8Frn5+HjYPh42A2vnsREZBswsch8ez8F5wtOAg9PDw5FxUS1pBhm+xH+XJR7lXg0Mi3pcCGyJw+uKSBdqbAmyYUcd67fV8eH2Oip3hct8054G6poCh8zr9xkDe2aR3zub/N7ZnHRCbwb0yqJfjyz65WbSt0cm/Xpk0qdHeLhndgY9s/1kZ/gTE/4n94a3sUdvmsnNJfu+e5lQ1D8x75nk4lHui4DbzWwhcBawT9vbRZJbbVMr71bu4d2qvazbWhsu890NB7cJZ2f4KBqQS2H/HpQU92dY/1yGDehBYf9chvTNoX9uFj5fe3+0e2TWrPD9nDlQVQVFRTB37sfju6FOy93MFgAXAnlmVg18D8gEcM79GlgMXAmUAQ3AlxIVVkSOzbZ9TSzfuDty28O6bbU4Bz6DEXk9GTe0D5+ZUMhJg3tx4gm9GT6wJ/5kKu9YzJrVrcu8rVi+LTOzk+kOuC1uiUTkuDnnWLu1lhdWb+OFNdtZv70OgNwsPxOL+vMvl5zIJ4r7c0ZRP3KzPDs5rCSQfqoiaSIYcrxbtYcXVm9jyZptVO9pxGdQUjyAOVeO5ayRAxg3pA8Zfh2Y3h2o3EVS3J76Fp4q3cTjb1WyeW8jWX4f544eyNcvHs2lY09gYK9sryOKB1TuIilq7ZZaHn1zI8+t3ExzIMRZIwZw55STuPjkQfTOyfQ6nnhM5S6SQpxzvLh2Ow+9UcGyit3kZPr47MQCbji7mLFD+ngdT5KIyl0kRbxfvY8f/GUtyyp2U9i/B/9+5clcWzKMfrlZXkeTJKRyF0ly22ubuH/Jev6wopoBuVnM/cx4risZph2jckQqd5Ek1dgSZN5r5fz6/z4iGHLMvmAkt100mj7ani4xULmLJKFlFbv516dWsnlvI1PHD+bbU8dSNDDX61iSQlTuIkkkFHL86tUyfvrihxQNyOWp2ZM5a+RAr2NJClK5iySJmrpm7nh6Ja9v2Mm004fyo8+eSq9s/YrKsdEnRyQJvFm2k288tZLaxlZ+/NlTmfGJYQdPIStyLFTuIh4Khhy/fHkDv3xlAyPyevLYlyfp++oSFyp3EY8EgiG++cx7PL9yC5+dUMAPrhpPT22GkTjRJ0nEAy2BEP+8YAVL1mzjW1ecxK0XjtJmGIkrlbtIF2tqDXLr/Hd5Zd0O/uOfxnHTeSO8jiRpSOUu0oUaW4LMfryU1zfs5IdXjef6ycO9jiRpSuUu0kX2Nwf48iPLKd24m59cfRrXlAzr/Ekix0jlLtIF9jW28sXfLWNV9T5+PmMC004f6nUkSXMqd5EEaw2GuOXxUlZv3scDn5/IlPGDvY4k3YDKXSTBfvDntSwt383/u+Z0Fbt0GZ0zVCSBFiyr4rG3Krn5vBF87sxCr+NIN6JyF0mQ5Rt3893nV3P+mDzunnqy13Gkm1G5iyTA5r2NfO2Jdyjsn8v/zJyoC2tIl9M2d5E4a2wJcsvjpTS1hlg4+0z65uriGtL1VO4iceSc487fr2LNlloeuqGE0YN6ex1Juin9rSgSR/NeK+dP723h3y4/iUvGnuB1HOnGVO4icbJuWy3/9bf1TB0/mFsvHOV1HOnmVO4icRAIhrjz2VX0yclk7mdO1RkexXPa5i4SBw++XsGq6vARqAN6ZnkdR0Rr7iLHq2zHfn720odMOWUwV56qI1AlOcRU7mY2xczWm1mZmd3dzvQiM/u7ma0ws1VmdmX8o4okn2DIceez75Gb5eeeq07R5hhJGp2Wu5n5gQeAqcA4YKaZjWsz23eAp51zE4AZwK/iHVQkGT3y5kberdrL9z49jkG9c7yOI3JQLGvuk4Ay51y5c64FWAhMbzOPAw5c1bcvsCV+EUWS08ad9fzkhXVcfPIgrjqjwOs4IoeIZYdqAbAp6nE1cFabeb4P/M3Mvg70BC6NSzqRJBUKOe76/SoyfT5+pG/HSBKKZc29vU+ta/N4JvCIc64QuBJ43MwOe20zm21mpWZWWlNTc/RpRZLE/LcrebtiN9/5p7EM7qvNMZJ8Yin3aiD6emCFHL7Z5SbgaQDn3FtADpDX9oWcc/OccyXOuZL8/PxjSyzisZ37m7lvyXrOG53HtbpUniSpWMp9OTDGzEaYWRbhHaaL2sxTBVwCYGZjCZe7Vs0lLf3ipQ00tgb5z+n6dowkr07L3TkXAG4HXgA+IPytmDVmdo+ZTYvM9k3gK2b2HrAA+KJzru2mG5GU91HNfp5cVsXMScMYld/L6zgiHYrpCFXn3GJgcZtx340aXgucG99oIsnnJ0vWk53h4xuXnOh1FJEj0hGqIjF6p3I3S9Zs45YLRpHfO9vrOCJHpHIXiYFzjh8tXkd+72xuPn+E13FEOqVyF4nBC2u2807lHu647ER6Zut8e5L8VO4inWgNhrh/yTpGD+rFNWcWeh1HJCYqd5FOLFy+ifKd9dw95WRd6FpShj6pIkewvznAL176kEkjBnDJ2EFexxGJmTYeihzBvNfK2bm/hYduHKsDliSlaM1dpAM1dc08+Fo5nzptCGcM6+d1HJGjonIX6cBv36igORDkm5fpgCVJPSp3kXbsa2zliaWVTD11CCN1mgFJQSp3kXY8sbSS/c0BvvbJUV5HETkmKneRNppag/zuHxVccGI+4wv6eh1H5Jio3EXaeLp0Ezv3t3DrhVprl9SlcheJ0hoM8Zv/K2diUT/OGjHA6zgix0zlLhLlT+9tYfPeRm69cLS+1y4pTeUuEhEKOf731Y846YTeXHyyjkaV1KZyF4l4ed0ONuzYz9cuHIXPp7V2SW0qdxHC52v/1atlFPbvwT+dNsTrOCLHTeUuAiwt382Kqr3ccsFInflR0oI+xSLAr14tI69XFteUDPM6ikhcqNyl21u9eR+vb9jJl88bQU6m3+s4InGhcpdu7+E3KuiZ5ef6ycO9jiISNyp36dZ27m/mz6u28rkzC+mTk+l1HJG4UblLt7ZwWRUtwRA3nF3sdRSRuFK5S7fVGgzxxNIqzh+Tx+hBOq2vpBeVu3Rbf1uznW21TVprl7Skcpdu69G3NlLYv4dONSBpSeUu3dIHW2tZVrGbG84ejl+nGpA0pHKXbunRNzeSk+njWh20JGlK5S7dzt6GFp5buZmrziigX26W13FEEkLlLt3O06WbaGoNceM5xV5HEUmYmMrdzKaY2XozKzOzuzuY51ozW2tma8zsyfjGFImP0BPzmTb9HMrv+zRjzz4N5s/3OpJIQmR0NoOZ+YEHgMuAamC5mS1yzq2NmmcM8G3gXOfcHjPT1w8k+cyfj/vKVxjc1Bh+XFkJs2eHh2fN8i6XSALEsuY+CShzzpU751qAhcD0NvN8BXjAObcHwDm3I74xReJgzhz8B4r9gIYGmDPHmzwiCRRLuRcAm6IeV0fGRTsRONHM/mFmS81sSnsvZGazzazUzEpramqOLbHIMXJVVe1P6Gi8SAqLpdzb+xKwa/M4AxgDXAjMBB4ys36HPcm5ec65EudcSX5+/tFmFTku+/IGtz+hqKhrg4h0gVjKvRqI/jJwIbClnXmed861OucqgPWEy14kKTS2BPnx+V+gJSvn0Am5uTB3rjehRBIolnJfDowxsxFmlgXMABa1mec54CIAM8sjvJmmPJ5BRY7Hn1dt4akxF1B1789h+HAwC9/Pm6edqZKWOv22jHMuYGa3Ay8AfuBh59waM7sHKHXOLYpMu9zM1gJB4FvOuV2JDC5yNJ5cVsWo/J6M+pfZ8K+3eB1HJOE6LXcA59xiYHGbcd+NGnbAHZGbSFL5YGstK6r28p1PjcVM55GR7kFHqEraW7CsiqwMH5+bWOh1FJEuo3KXtNbYEuSP727myvGD6d9T55GR7kPlLmntT6u2UNcc4PNn6eLX0r2o3CWtLVhWxehBvfhEcX+vo4h0KZW7pK0DO1JnTirSjlTpdlTukrY+3pHa9mwZIulP5S5pqaElwB/f3cynTh2iC3JIt6Ryl7T051VbqWsOMHOSzhsj3ZPKXdLSk29rR6p0byp3STtrt9SycpN2pEr3pnKXtKMdqSIqd0kzDS0BnlsRPiJVO1KlO1O5S1o5sCNVR6RKd6dyl7SyMHJqX+1Ile5O5S5pY/22Ot7VEakigMpd0siCZVVk+X18Vqf2FVG5S3poag3yh3eruWL8YAbo1L4iKndJD4vf30ptU4CZk4Z1PrNIN6Byl7SwcNkmigfmcvbIgV5HEUkKKndJeWU76li2cTcztCNV5CCVu6S8Bcs2kek3rj5TO1JFDlC5S0o7sCP18nGDyeuV7XUckaShcpeU9sKabexpaGWGdqSKHELlLilt4bJNDBvQg3NH5XkdRSSpqNwlZZXX7Oet8l3M+EQRPp92pIpEU7lLynpq+Sb8PuMa7UgVOYzKXVJScyDIs+9Uc+nYQQzqk+N1HJGko3KXlLRk9TZ21bfoGqkiHVC5S0qav7SK4QNzuWBMvtdRRJKSyl1SzrpttSzbuJvrzxquHakiHYip3M1sipmtN7MyM7v7CPNdbWbOzEriF1HkUE8srSQrw6cjUkWOoNNyNzM/8AAwFRgHzDSzce3M1xv4Z+DteIcUOaCuqZU/vruZT582lP46ta9Ih2JZc58ElDnnyp1zLcBCYHo78/0AuB9oimM+kUM8t2Iz9S1BvnC2rpEqciSxlHsBsCnqcXVk3EFmNgEY5pz7cxyziRzCOccTS6s4taAvpxf29TqOSFKLpdzb22PlDk408wE/A77Z6QuZzTazUjMrrampiT2lCLB84x7Wb6/jC5OH69S+Ip2IpdyrgeizMhUCW6Ie9wbGA6+a2UZgMrCovZ2qzrl5zrkS51xJfr6+wiZH5/GllfTJyeDTpw/1OopI0oul3JcDY8xshJllATOARQcmOuf2OefynHPFzrliYCkwzTlXmpDE0i3tqGtiyeqtXFMyjB5Zfq/jiCS9TsvdORcAbgdeAD4AnnbOrTGze8xsWqIDigA8vXwTrUHHrLN0RKpILDJimck5txhY3GbcdzuY98LjjyXysWDI8eTbVZw3Oo+R+b28jiOSEnSEqiS9V9btYMu+Jq6frK8/isRK5S5J7/GllQzuk8OlYwd5HUUkZajcJalV7KzntQ9r+PxZRWT49XEViZV+WySpPfxGBVl+n66RKnKUVO6StPbUt/DMO5u4asJQBvXWBTlEjobKXZLWE0sraWoNcfP5I72OIpJyVO6SlJpagzz61kYuOimfE0/o7XUckZSjcpek9NyKzezc38JXtNYuckxU7pJ0QiHHQ29UcMrQPpw9aqDXcURSkspdks6rH+6gbMd+Zl8wUmd/FDlGKndJOvNeK2dI3xyuPHWI11FEUpbKXZLK+9X7WFq+my+fO4JMHbQkcsz02yNJ5cHXy+mVncF1OmhJ5Lio3CVpbN7byF/e38rMScPok5PpdRyRlKZyl6TxuzcqMOBL547wOopIylO5S1LY19jKgmVVfOq0IQzt18PrOCIpT+UuSeGJpZXUtwR10JJInKjcxXO1Ta3Me62ci08exPiCvl7HEUkLKnfx3G9fr2BfYyt3XHai11FE0obKXTy1p76F375RwdTxg7XWLhJHKnfx1G9eK6e+JcC/aq1dJK5U7uKZHXVNPPJmBdNPH6rT+orEmcpdPPOrv39Ea9DxjUu11i4Sbyp38cSWvY08+XYVV08sZEReT6/jiKQdlbt44r9fKcPh+Polo72OIpKWVO7S5Sp31fNM6SZmTiqisH+u13FE0pLKXbrcL17egN9n3H6R1tpFEkXlLl2qbEcdz63YzI3nFDOoT47XcUTSlspdutR9S9bTI9PPLRfoHDIiiaRyly7z8gfbeXHtdm67eDQDe2V7HUckrancpUs0tgT53qI1jBnUi5vP01q7SKLFVO5mNsXM1ptZmZnd3c70O8xsrZmtMrOXzWx4/KNKSpk/H4qLweeD4mJe+o+fUr2nkR9cNZ6sDK1TiCRap79lZuYHHgCmAuOAmWY2rs1sK4AS59xpwLPA/fEOKilk/nyYPRsqK8E5qKzkkp/+Bz+sf4/JIwd6nU6kW4hlFWoSUOacK3fOtQALgenRMzjn/u6ca4g8XAoUxjempJQ5c6Ch4ZBRuYFmZj7/G48CiXQ/sZR7AbAp6nF1ZFxHbgL+2t4EM5ttZqVmVlpTUxN7SkktVVXtjvZXb2p3vIjEXyzlbu2Mc+3OaHY9UAL8pL3pzrl5zrkS51xJfn5+7CkltRQVHd14EYm7WMq9GhgW9bgQ2NJ2JjO7FJgDTHPONccnnqSkuXMht81pBXJzw+NFpEvEUu7LgTFmNsLMsoAZwKLoGcxsAvAbwsW+I/4xJaXMmkXVvT+nuk8+zgyGD4d582DWLK+TiXQbGZ3N4JwLmNntwAuAH3jYObfGzO4BSp1ziwhvhukFPGNmAFXOuWkJzC1JrDkQ5GvuZHbc+SQvf/OT9MnJ9DqSSLfTabkDOOcWA4vbjPtu1PClcc4lKezHi9exZkst875wpopdxCM6mkTi6q/vb+WRNzdy03kjuPyUwV7HEem2VO4SN5W76rnz2VWcPqwfd0052es4It2ayl3iojkQ5LYn38UMHvj8BJ1iQMRjMW1zF+nMj/7yAas31/LgDSW6upJIEtDqlRy3v6zayqNvVXLzeSO4bNwJXscREVTucpw27qznrt+vYkJRP+6aqu3sIslC5S7HrKElwG1PvovfZ/z3zAlk+vVxEkkW+m2UY9LUGuQrj5XywdZafnbd6drOLpJktENVjlprMMTtT67gH2W7+K9rTufik7WdXSTZaM1djkow5Pi3Z97jpQ+285/TTuHqM3XqfpFkpHKXmDnn+M5zq3l+5Ra+dcVJ3HhOsdeRRKQDKneJiXOOHy3+gAXLqrj1wlHcdtForyOJyBGo3CUmv3y5jAdfr+CGs4fzrStO8jqOiHRCO1TliIIhx31L1jHvtXI+N7GQ73/6FCKndRaRJKZylw7VNrXyzwtW8Or6Gr4weTjf+/Q4fD4Vu0gqULlLuyp21nPTo8up2tXA3M+MZ9ZZw72OJCJHQeUuh3l9Qw23zQ8fefrEzWcxeeRAryOJyFFSuctBzjl+94+N/PAvaznxhN48eEMJwwboyFORVKRyFwC21zbxvefXsGTNNq445QR+eu0Z9MzWx0MkVem3t5sLhhzz367k/iXraQ2GuGvKydxywUjtOBVJcSr3buyDrbV8+w/vs3LTXs4fk8cPrxrP8IE9vY4lInGgcu+GGluC/PzlD3no9Qr69cjk59edwfQzhur76yJpROXejdQ1tTL/7Sp++0YFNXXNXFtSyL9fOZZ+uVleRxOROFO5dwO79jfzu39s5LG3NlLbFOD8MXn8atZEPlE8wOtoIpIgKvc0tnlvIw++Vs7C5VU0B0JcMW4wt140itMK+3kdTUQSTOWeZvY1tPLX1Vt5buVm3q7Yjd+MqyYU8NVPjmT0oN5exxORLqJyTwNNrUFeWbeD51Zs5tX1NbQEQ4zI68k3LhnDNSXDKOjXw+uIItLFVO4pyDnH+u11/KNsF2+W7WRp+S7qW4Lk987m+snDmX7GUE4r7Ktvv4h0Yyr3FBAIhviopp6Vm/aEC/2jXezc3wxA8cBcpk8oYOr4wZw9ciAZfp2iX0RU7kmnoSVA2Y79rN5cy5ot+1i9pZZ1W2tpDoQAyO+dzXmjB3LO6DzOHZ2nTS4i0q6Yyt3MpgC/APzAQ865e9tMzwYeA84EdgHXOec2xjdqegiFHDvrm9lR28ym3Q1U7Kqncmfkflc922ubD87bOyeD8UP78oXJwxlf0JfxBX0Zld9Tm1tEpFOdlruZ+YEHgMuAamC5mS1yzq2Nmu0mYI9zbrSZzQDuA65LROCEmD8f5syBqiooKoK5c2HWrJieGgw56lsC1DUF2FPfwp6GFvY0tH48XN/CjrpmttU2sX1fEzvqmgmE3CGvkdcrm+KBuZw/Jp/igbmMzO/F+KF9GTagh4pcRI5JLGvuk4Ay51w5gJktBKYD0eU+Hfh+ZPhZ4H/MzJxzh7ZYHARDjtZgiJBzBEOOUAiCB4adIxByBIOOQChEMBR5HHlOa9ARCIZoCYYIBMPj8v/0e874wZ1kNDWG36Cyktabbmbxe1tYcd6VNLQEaGwN0dgSoLE1SENLkLqmAPubAuxvDt+OpE9OBoP65DC4Tw6TRw1kcJ8cBvfN4YQ+ORT068Hwgbn0zsmM92ISkW4ulnIvADZFPa4GzupoHudcwMz2AQOBnfEIGe3B18u596/r4vZ6b/zvPR8Xe0RmcxMlv/kvvtvzdHKz/PTI9NMjct8rO4MhfXPonZ1Jr5wMemVn0DsnfOuXm0X/3CwG9MykX24W/XpkageniHgilnJvb7tA2zXyWObBzGYDswGKiopieOvDTRoxgG9dcRJ+n+E3w+cz/AZ+X3g4w2f4fT4yfEaGP/zYZ0am30em30eG/8CwkeHzUXB/+///FNTV8N73Lj+mjCIiXoul3KuBYVGPC4EtHcxTbWYZQF9gd9sXcs7NA+YBlJSUHNMmm4lF/ZlY1P9Yntq+oiKorGx/vIhIioplm8FyYIyZjTCzLGAGsKjNPIuAGyPDVwOvJGJ7e0LMnQu5bS4ll5sbHi8ikqI6LXfnXAC4HXgB+AB42jm3xszuMbNpkdl+Cww0szLgDuDuRAWOu1mzYN48GD4czML38+bF/G0ZEZFkZF6tYJeUlLjS0lJP3ltEJFWZ2TvOuZLO5tNXOURE0pDKXUQkDancRUTSkMpdRCQNqdxFRNKQZ9+WMbMaoJ2jh2KSRwJObRAHynV0lOvoJWs25To6x5NruHMuv7OZPCv342FmpbF8FairKdfRUa6jl6zZlOvodEUubZYREUlDKncRkTSUquU+z+sAHVCuo6NcRy9ZsynX0Ul4rpTc5i4iIkeWqmvuIiJyBElb7mZ2jZmtMbOQmXW4V9nMppjZejMrM7O7o8aPMLO3zWyDmT0VOV1xPHINMLMXI6/7opkddnJ5M7vIzFZG3ZrM7KrItEfMrCJq2hldlSsyXzDqvRdFjfdyeZ1hZm9Fft6rzOy6qGlxXV4dfV6ipmdH/v1lkeVRHDXt25Hx683siuPJcQy57jCztZHl87KZDY+a1u7PtItyfdHMaqLe/+aoaTdGfu4bzOzGts9NcK6fRWX60Mz2Rk1L5PJ62Mx2mNnqDqabmf0yknuVmU2Mmhbf5eWcS8obMBY4CXgVKOlgHj/wETASyALeA8ZFpj0NzIgM/xr4Wpxy3Q/cHRm+G7ivk/kHEL5wSW7k8SPA1QlYXjHlAvZ3MN6z5QWcCIyJDA8FtgL94r28jvR5iZrnVuDXkeEZwFOR4XGR+bOBEZHX8XdhrouiPkNfO5DrSD/TLsr1ReB/2nnuAKA8ct8/Mty/q3K1mf/rwMOJXl6R174AmAis7mD6lcBfCV+9bjLwdqKWV9KuuTvnPnDOre9ktoMX73bOtQALgelmZsDFhC/WDfAocFWcok2PvF6sr3s18FfnXEOc3r8jR5vrIK+Xl3PuQ+fchsjwFmAH0OlBGseg3c/LEfI+C1wSWT7TgYXOuWbnXAVQFnm9LsnlnPt71GdoKeEroiVaLMurI1cALzrndjvn9gAvAlM8yjUTWBCn9z4i59xrtHMVuijTgcdc2FKgn5kNIQHLK2nLPUbtXby7gPDFufe68IVGosfHwwnOua0AkftBncw/g8M/WHMjf5L9zMyyuzhXjpmVmtnSA5uKSKLlZWaTCK+NfRQ1Ol7Lq6PPS7vzRJbHgYu9x/LcROaKdhPhtb8D2vuZdmWuz0V+Ps+a2YFLcibF8opsvhoBvBI1OlHLKxYdZY/78orlGqoJY2YvAYPbmTTHOfd8LC/Rzjh3hPHHnSvW14i8zhDgVMJXsTrg28A2wgU2D7gLuKcLcxU557aY2UjgFTN7H6htZz6vltfjwI3OuVBk9DEvr/beop1xsV7s/bg+U52I+bXN7HqgBPhk1OjDfqbOuY/ae34Ccv0JWOAVOZDKAAACdUlEQVScazazrxL+q+fiGJ+byFwHzACedc4Fo8YlannFoss+X56Wu3Pu0uN8iY4u3r2T8J87GZG1r/Yu6n1Mucxsu5kNcc5tjZTRjiO81LXAH51zrVGvvTUy2GxmvwP+rStzRTZ74JwrN7NXgQnA7/F4eZlZH+AvwHcif64eeO1jXl7tOJ6Lvcfy3ETmwswuJfwf5iedc80HxnfwM41HWXWayzm3K+rhg8B9Uc+9sM1zX41DpphyRZkB3BY9IoHLKxYdZY/78kr1zTLtXrzbhfdQ/J3w9m4IX7w7lr8EYhF9MfDOXvewbX2RgjuwnfsqoN296onIZWb9D2zWMLM84FxgrdfLK/Kz+yPhbZHPtJkWz+V1PBd7XwTMsPC3aUYAY4Blx5HlqHKZ2QTgN8A059yOqPHt/ky7MNeQqIfTCF9nGcJ/rV4eydcfuJxD/4JNaK5ItpMI75x8K2pcIpdXLBYBN0S+NTMZ2BdZgYn/8krUXuPjvQGfIfy/WTOwHXghMn4osDhqviuBDwn/zzsnavxIwr98ZcAzQHaccg0EXgY2RO4HRMaXAA9FzVcMbAZ8bZ7/CvA+4ZJ6AujVVbmAcyLv/V7k/qZkWF7A9UArsDLqdkYilld7nxfCm3mmRYZzIv/+ssjyGBn13DmR560Hpsb5895ZrpcivwcHls+izn6mXZTrx8CayPv/HTg56rlfjizHMuBLXZkr8vj7wL1tnpfo5bWA8Le9Wgn3103AV4GvRqYb8EAk9/tEfRMw3stLR6iKiKShVN8sIyIi7VC5i4ikIZW7iEgaUrmLiKQhlbuISBpSuYuIpCGVu4hIGlK5i4ikof8PFuDT/aSd4BsAAAAASUVORK5CYII=\n", "text/plain": [ "