diff --git a/docs/source/notebooks/simple_stochastic_optimization.ipynb b/docs/source/notebooks/simple_stochastic_optimization.ipynb new file mode 100644 index 0000000000..6580f4290c --- /dev/null +++ b/docs/source/notebooks/simple_stochastic_optimization.ipynb @@ -0,0 +1,794 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true + }, + "source": [ + "# Stochastic Optimization using inferred model" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true + }, + "source": [ + "Condider a problem where you have to perform optimization for decision making. You don't know about the world very much and thus you have to do bayesian inference to discover it.\n", + "\n", + "This formulation covers a wide range of practical tasks, here are some of them:\n", + " * Maximizing stochastic output of model $E_{p(\\theta|D)}[f(X)] \\rightarrow max X$\n", + " * Minimizing costs of inputs to archive target $E_{p(y|X, \\theta, D)}[L(y_{target}, X, y)] \\rightarrow min X$\n", + " * Minimizing Bayesian risk: $E_{p(y|X, \\theta, D)}[L(y_p, y)] \\rightarrow min y_p$\n" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": true, + "deletable": true, + "editable": true + }, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "import functools\n", + "import numpy as np\n", + "from theano import theano, tensor as tt\n", + "import matplotlib.pyplot as plt\n", + "import pymc3 as pm\n", + "np.random.seed(42)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true + }, + "source": [ + "## Quadratic problem\n", + "Consider a simple quadratic problem with unknown parameters. Here comes Bayesian inference to get model with that we need for our opimization" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": true, + "deletable": true, + "editable": true + }, + "outputs": [], + "source": [ + "def f(x, a, b, c):\n", + " return a*x**2 + b*x + c" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true + }, + "source": [ + "### True data\n", + "Generate true data" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false, + "deletable": true, + "editable": true + }, + "outputs": [], + "source": [ + "x_obs = np.random.uniform(-5, 5, size=(15,)).astype('float32')\n", + "a, b, c = 1, 2, 3\n", + "min_ = np.array([-b/2/a])\n", + "y_obs = f(x_obs, a, b, c) + np.random.normal(size=x_obs.shape).astype('float32')" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true + }, + "source": [ + "### Final task\n", + "Our task is to find th minimum of the unknown function without evaluating it. That differs a lot from traditional Bayesian optimization using Gaussian Process. Here we have to deal with uncertainty in loss.\n", + "\n", + "Now let's look at true optima" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false, + "deletable": true, + "editable": true + }, + "outputs": [ + { + "data": { + "text/plain": [ + "array([-1.])" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "min_" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true + }, + "source": [ + "And at the data" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false, + "deletable": true, + "editable": true + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEF1JREFUeJzt3X2IZFV+xvHnmXFYt6NBzRRmcOzuZSNZZJOdCZ3BYEiM\nL8uskdWFEDAVMUToDSgomBfX/mM3kIYNyWr+SNhQG40DqayIL7iIm+zECCIkJjXuOI7OJm7MdMdh\ndMoYUWlwGeeXP+5tZqa32rpVdW+9nPp+oLl1T93q+yttH0+dOvdcR4QAAJNvy6gLAACUg0AHgEQQ\n6ACQCAIdABJBoANAIgh0AEgEgQ4AiSDQASARBDoAJOKcYZ5s+/btMT8/P8xTAsDEO3DgwNsRUet2\n3FADfX5+Xq1Wa5inBICJZ3ulyHEMuQBAIgh0AEgEgQ4AiSDQASARBDoAJKJroNs+1/a/2X7J9iu2\n/zhvf8j2f9s+mP/sqr5cAJgszaY0Py9t2ZJtm83qzlVk2uKHkq6OiA9sb5P0vO3v5s/9QUQ8Wl15\nADC5mk1pcVFaW8v2V1ayfUmq18s/X9ceemQ+yHe35T/ctw4AulhaOh3m69bWsvYqFBpDt73V9kFJ\nJyTtj4gX8qeWbR+yfb/tT1RTIgBMptXV3toHVSjQI+KjiNglaaekPbY/K+krkj4j6RclXSTpjzq9\n1vai7ZbtVrvdLqlsABh/s7O9tQ+qp1kuEfGupGcl7Y2I4/lwzIeS/lbSnk1e04iIhYhYqNW6LkUA\nAMlYXpZmZs5um5nJ2qtQZJZLzfYF+eNPSrpO0g9s78jbLOkmSYerKREAJlO9LjUa0tycZGfbRqOa\nL0SlYrNcdkjaZ3ursv8BPBIRT9n+Z9s1SZZ0UNLvVVMiAEyuer26AN+oa6BHxCFJuzu0X11JRQCA\nvnClKAAkgkAHgEQQ6ACQCAIdABJBoANAIgh0AEgEgQ4AiSDQASARBDoAJIJAB4BEEOgAkAgCHQAS\nQaADQCIIdAAYULMpzc9LW7Zk22ZzNHUUWQ8dALCJZlNaXDx9M+iVlWxfGt466OvooQPAAJaWTof5\nurW1rH3YCHQAGMDqam/tVSLQAWAAs7O9tVeJQAeAASwvSzMzZ7fNzGTtw0agA8AA6nWp0ZDm5iQ7\n2zYaw/9CVCowy8X2uZKek/SJ/PhHI+Krtj8l6WFJPyXpgKRbIuJHVRYLAOOoXh9NgG9UpIf+oaSr\nI+JzknZJ2mv7Ckl/Kun+iPgZSf8n6bbqygQAdNM10CPzQb67Lf8JSVdLejRv3yfppkoqBAAUUmgM\n3fZW2wclnZC0X9J/SXo3Ik7mh7wh6ZJqSgQAFFEo0CPio4jYJWmnpD2SPlP0BLYXbbdst9rtdp9l\nAgC66WmWS0S8K+lZSb8k6QLb61+q7pR0bJPXNCJiISIWarXaQMUCADbXNdBt12xfkD/+pKTrJB1R\nFuy/kR92q6QnqyoSANBdkcW5dkjaZ3ursv8BPBIRT9l+VdLDtv9E0vclPVBhnQCALroGekQckrS7\nQ/vrysbTAQBjgCtFASARBDoAJIJAB4BEEOgAkAgCHQASQaADQCIIdABIBIEOAIkg0AEgEQQ6ACSC\nQAeARBDoAJAIAh0AEkGgA0AiCHQASASBDgCJINABIBEEOgAkgkAHgEQQ6ACQCAIdABLRNdBtX2r7\nWduv2n7F9p15+9dsH7N9MP+5vvpyAQCbOafAMScl3R0RL9o+X9IB2/vz5+6PiD+vrjwAQFFdAz0i\njks6nj9+3/YRSZdUXRgAoDc9jaHbnpe0W9ILedMdtg/ZftD2hSXXBgDoQeFAt32epMck3RUR70n6\npqRPS9qlrAf/jU1et2i7ZbvVbrdLKBkA0EmhQLe9TVmYNyPicUmKiLci4qOIOCXpW5L2dHptRDQi\nYiEiFmq1Wll1AwA2KDLLxZIekHQkIu47o33HGYd9SdLh8ssDABRVZJbLlZJukfSy7YN5272Sbra9\nS1JIOirpy5VUCAAopMgsl+clucNTT5dfDgCgX1wpCgCJINABIBEEOoCJ0WxK8/PSli3ZttkcdUXj\nhUAHMBGaTWlxUVpZkSKy7eLi5qE+jeFPoAOYCEtL0tra2W1ra1n7Rr2GfyoIdAATYXW1eHsv4Z8S\nAh3ARJidLd7eS/inhEAHMBGWl6WZmbPbZmay9o16Cf+UEOgAJkK9LjUa0tycZGfbRiNr36iX8E9J\nkUv/AWAs1OudA7zTcVI2Zr66mvXMl5eLvXaSEegAklQ0/FPCkAsAJIJAB4BEEOgAkAgCHQASQaAD\nQCIIdABIBIEOAIkg0AEgEQQ6ACSia6DbvtT2s7Zftf2K7Tvz9ots77f9Wr69sPpyAQCbKdJDPynp\n7oi4XNIVkm63fbmkeyQ9ExGXSXom3wcAjEjXQI+I4xHxYv74fUlHJF0i6UZJ+/LD9km6qaoiAQDd\n9TSGbnte0m5JL0i6OCKO50+9KeniUisDAPSkcKDbPk/SY5Luioj3znwuIkJSbPK6Rdst2612uz1Q\nsQCAzRUKdNvblIV5MyIez5vfsr0jf36HpBOdXhsRjYhYiIiFWq1WRs0AgA6KzHKxpAckHYmI+854\n6juSbs0f3yrpyfLLAwAUVeQGF1dKukXSy7YP5m33Svq6pEds3yZpRdJvVlMiAKCIroEeEc9L8iZP\nX1NuOQCAfnGlKAAkgkAHgEQQ6ACQCAIdABJBoANAIgh0AEgEgQ4AiSDQASARBDoAJIJAB4BEEOgA\nkAgCHQASQaADQCIIdABIBIEOAIkg0AEgEQQ6ACRi7AO92ZTm56UtW7JtsznqigBgPBW5p+jINJvS\n4qK0tpbtr6xk+5JUr4+uLgAYR2PdQ19aOh3m69bWsnYAwNnGOtBXV3trB4Bp1jXQbT9o+4Ttw2e0\nfc32MdsH85/rqyhudra3dgCYZkV66A9J2tuh/f6I2JX/PF1uWZnlZWlm5uy2mZmsHQBwtq6BHhHP\nSXpnCLX8mHpdajSkuTnJzraNBl+IAkAng4yh32H7UD4kc+FmB9letN2y3Wq32z2fpF6Xjh6VTp3K\ntoQ5AHTWb6B/U9KnJe2SdFzSNzY7MCIaEbEQEQu1Wq3P0wEAuukr0CPirYj4KCJOSfqWpD3llgUA\n6FVfgW57xxm7X5J0eLNjAQDD0fVKUdvflnSVpO2235D0VUlX2d4lKSQdlfTlCmsEABTQNdAj4uYO\nzQ9UUAsAYABjfaUoAKA4Ah0AEpFcoLPcLoBpNdbL5/aK5XYBTLOkeugstwtgmiUV6Cy3C2CaJRXo\nLLcLYJolFegstwtgmiUV6Cy3C2CaJTXLRcrCmwAHMI2S6qEDwDQj0AEgEQQ6ACSCQAeARBDoAJAI\nAh0AEkGgA0AiCHQASASBDqAw7jcw3pK7UhRANbjfwPjr2kO3/aDtE7YPn9F2ke39tl/LtxdWWyaA\nUeN+A+OvyJDLQ5L2bmi7R9IzEXGZpGfyfQAJ434D469roEfEc5Le2dB8o6R9+eN9km4quS4AY4b7\nDYy/fr8UvTgijueP35R0cUn1ABhT3G9g/A08yyUiQlJs9rztRdst2612uz3o6QCMCPcbGH/9Bvpb\ntndIUr49sdmBEdGIiIWIWKjVan2eDsCwfNzUxHpdOnpUOnUq2xLm46XfQP+OpFvzx7dKerKccgCM\n0vrUxJUVKeL01ETmm0+GItMWvy3pXyT9rO03bN8m6euSrrP9mqRr830AE46piZOtyCyXmyNiR0Rs\ni4idEfFARPxvRFwTEZdFxLURsXEWzNTiSjpMMqYmTjYu/S8RH1cx6ZiaONkI9BLxcRWTjqmJk41A\nLxEfVzHpmJo42Vicq0Szs9kwS6d2YFLU6wT4pKKHXiI+rgIYJQK9RHxcBTBKDLmUjI+rAEaFHjoA\nJIJAB4BEEOgAkAgCfchYGgCd8HeBMvCl6BBxk110wt8FyuLs/hTDsbCwEK1Wa2jnGzfz850vPJqb\ny9aWxnTi7wLd2D4QEQvdjmPIpWQf99GZpQHQCX8XKAuBXqJuqy2ykh064e8CZSHQS9RttUWWBkAn\n/F2gLAR6ibp9dO5laQBmPUwPloxAWfhStERlfbm1cdaDlPXY+I8cmE58KToCZX105kYZAPpBoJeo\nrI/OzHqYPgyxoQxcWFSyMlZb5EYZ04ULi1CWgXroto/aftn2QdvpDo4PGbMepgtDbChLGT30X4uI\nt0v4Pcit98qWlrJhltnZLMzpraWJITaUhSGXMcWNMqYHQ2woy6Bfioak79k+YHux0wG2F223bLfa\n7faApwPSwxAbyjJooP9yRPyCpC9Iut32r2w8ICIaEbEQEQu1Wm3A0wHp4cIilGWgIZeIOJZvT9h+\nQtIeSc+VURgwTRhiQxn67qHb/gnb568/lvR5SYfLKgwA0JtBeugXS3rC9vrv+fuI+IdSqgIA9Kzv\nQI+I1yV9rsRaAAAD4NJ/AEgEgQ4AiSDQASARBDoAJIJAB4BEEOgAkAgCHQASQaADQCIIdABIBIEO\nAIkg0AEgEQQ6ACSCQMdQNJvS/Ly0ZUu2bTZHXRGQHu4piso1m9Li4uk726+sZPsSN3UAykQPHZVb\nWjod5uvW1rJ2AOUh0FG51dXe2gH0h0BH5WZne2sH0B8CHZVbXpZmZs5um5nJ2gGUh0BH5ep1qdGQ\n5uYkO9s2GnwhCpSNWS4YinqdAAeqNlAP3fZe2/9h+4e27ymrKABA7/oOdNtbJf2VpC9IulzSzbYv\nL6swAEBvBumh75H0w4h4PSJ+JOlhSTeWUxYAoFeDBPolkv7njP038jYAwAhUPsvF9qLtlu1Wu92u\n+nQAMLUGmeVyTNKlZ+zvzNvOEhENSQ1Jst22vTLAOYdhu6S3R13ECPH+ef+8//EzV+QgR0Rfv932\nOZL+U9I1yoL83yX9VkS80tcvHBO2WxGxMOo6RoX3z/vn/U/u+++7hx4RJ23fIekfJW2V9OCkhzkA\nTLKBLiyKiKclPV1SLQCAAXDp/49rjLqAEeP9Tzfe/wTrewwdADBe6KEDQCII9E3Yvtt22N4+6lqG\nyfaf2f6B7UO2n7B9wahrGoZpX5fI9qW2n7X9qu1XbN856pqGzfZW29+3/dSoa+kXgd6B7UslfV7S\nNN5TZ7+kz0bEzyublvqVEddTOdYlkiSdlHR3RFwu6QpJt0/hP4M7JR0ZdRGDINA7u1/SH0qaui8Y\nIuJ7EXEy3/1XZReMpW7q1yWKiOMR8WL++H1lwTY1S3nY3inp1yX9zahrGQSBvoHtGyUdi4iXRl3L\nGPhdSd8ddRFDwLpEZ7A9L2m3pBdGW8lQ/YWyTtypURcyiKm8wYXtf5L00x2eWpJ0r7LhlmR93PuP\niCfzY5aUfQxvDrM2jJbt8yQ9JumuiHhv1PUMg+0bJJ2IiAO2rxp1PYOYykCPiGs7tdv+OUmfkvSS\nbSkbbnjR9p6IeHOIJVZqs/e/zvbvSLpB0jUxHfNaC61LlDrb25SFeTMiHh91PUN0paQv2r5e0rmS\nftL230XEb4+4rp4xD/1j2D4qaSEixnGxnkrY3ivpPkm/GhFTsTxmqusS9cJZD2afpHci4q5R1zMq\neQ/99yPihlHX0g/G0LHRX0o6X9J+2wdt//WoC6pa/iXw+rpERyQ9Mk1hnrtS0i2Srs7/vR/Me6yY\nIPTQASAR9NABIBEEOgAkgkAHgEQQ6ACQCAIdABJBoANAIgh0AEgEgQ4Aifh/OZkP5EpigTIAAAAA\nSUVORK5CYII=\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.plot(x_obs, y_obs, 'bo');" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true + }, + "source": [ + "For inference you can use every method you like. Here I use SVGD just to show that it's result can be used in inferece as well as any other tool in variational module." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false, + "deletable": true, + "editable": true, + "scrolled": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "WARNING (theano.tensor.blas): We did not found a dynamic library into the library_dir of the library we use for blas. If you use ATLAS, make sure to compile it with dynamics library.\n", + "100%|██████████| 4000/4000 [00:24<00:00, 166.11it/s]\n" + ] + } + ], + "source": [ + "sgd = functools.partial(pm.sgd, learning_rate=.001) # simple problem - simple optimizer\n", + "with pm.Model() as model:\n", + " abc = pm.Normal('abc', sd=1, shape=(3,))\n", + " x = theano.shared(x_obs, 'x_obs')\n", + " x2 = x**2\n", + " o = tt.ones_like(x)\n", + " X = tt.stack([x2, x, o]).T\n", + " y = X.dot(abc)\n", + " pm.Normal('y', mu=y, observed=y_obs)\n", + " histogram = pm.fit(4000, method='svgd', obj_optimizer=sgd)\n", + " # or like that if you want to customize SVGD\n", + " # svgd = pm.SVGD(n_particles=500, jitter=.1)\n", + " # histogram = svgd.fit(1000, obj_optimizer=sgd)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true + }, + "source": [ + "As a result we have a `Histogram` instance. Histogram is a linkage from raw samples to theano computational engine. If we use NUTS for example there is no direct way to use trace results in computational stuff except Histogram, we'll see it later." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "histogram" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": false, + "deletable": true, + "editable": true + }, + "outputs": [ + { + "data": { + "text/plain": [ + "array([ 1.0527544 , 1.99613583, 1.92554891], dtype=float32)" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "histogram.histogram.get_value().mean(0)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true + }, + "source": [ + "It contains particles inferred with SVGD, shared variable has exactly 100 samples for 3 variables." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false, + "deletable": true, + "editable": true + }, + "outputs": [ + { + "data": { + "text/plain": [ + "(theano.tensor.sharedvar.TensorSharedVariable, (100, 3))" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "type(histogram.histogram), histogram.histogram.get_value().shape" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true + }, + "source": [ + "Here is an option to use NUTS for inference and reuse it's trace directly in optimization. Let's do it!" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "collapsed": false, + "deletable": true, + "editable": true + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/root/macos/dev/pymc3/pymc3/sampling.py:188: UserWarning: Instantiated step methods cannot be automatically initialized. init argument ignored.\n", + " warnings.warn('Instantiated step methods cannot be automatically initialized. init argument ignored.')\n", + "100%|██████████| 10000/10000 [00:20<00:00, 487.37it/s]\n" + ] + } + ], + "source": [ + "with model:\n", + " step = pm.NUTS()\n", + " trace = pm.sample(10000, step=step)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA04AAACICAYAAADQ+bu1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvXmUXNd93/m5tXZ1VfW+r1gIgiBAgMTCnbQokrJER/Y4\nJ4rGntiO4liZSeLEJ3EcjzOJ7Th27CSaM2N7cmzJiRc59kxi2ZYtUwskkwRFUQQBEASIjdjR+1bV\nVdW1L3f+uPVqfbV1d/WCvp9z+nQtb7nvvfte/b7397u/n5BSotFoNBqNRqPRaDSaylg2uwEajUaj\n0Wg0Go1Gs9XRwkmj0Wg0Go1Go9FoaqCFk0aj0Wg0Go1Go9HUQAsnjUaj0Wg0Go1Go6mBFk4ajUaj\n0Wg0Go1GUwMtnDQajUaj0Wg0Go2mBlo4aTRbACHE3xVCfHuz26HRaDSanYn+HdJoaqOFk0aj0Wg0\nGo1Go9HUQAsnjUaj0Wg0Go1Go6mBFk4azQYihPhZIcRNIURICHFZCPGDxV+L3xRCBIQQV4UQLxZ8\n0SWE+F0hxLQQwi+E+PNNaL5Go9Fotjn6d0ijWT22zW6ARrPDuAk8B8wCnwL+UAjxQPa7J4A/AXqA\nvwn8qRBit5TSB3wRWAEOZv8/vdEN12g0Gs19gf4d0mhWiZBSbnYbNJodixDiPPDzQCfwK8CwzN6U\nQojTwG8A3wSmgG4ppX+z2qrRaDSa+w/9O6TR1I8O1dNoNhAhxI8KIc4LIZaFEMvAIdTIHsCULB7J\nuAsMAaOAT/9YaTQajWat6N8hjWb1aOGk0WwQQohx4AvAP0aN2nUAHwAiu8iwEEIUrDIGTAMTQJcQ\nomMj26vRaDSa+wv9O6TRrA0tnDSajcMNSGABQAjxGdRIn0Ef8E+EEHYhxKeAA8CrUsoZ4KvAfxZC\ndGa/f36D267RaDSa7Y/+HdJo1oAWThrNBiGlvAx8DngbmAMeAd4qWOQdYB+wCPwy8LeklEvZ734E\nSAJXgXngpzao2RqNRqO5T9C/QxrN2tDJITQajUaj0Wg0Go2mBtrjpNFoNBqNRqPRaDQ10MJJo9Fo\nNBqNRqPRaGqghZNGo9FoNBqNRqPR1EALJ41Go9FoNBqNRqOpgW2zG1APPT09cteuXZvdDI1Go9Gs\nM2fPnl2UUvZudjtqoX+HNBqN5v6kkd+hbSGcdu3axZkzZza7GRqNRqNZZ4QQdze7DfWgf4c0Go3m\n/qSR3yEdqqfRaDQajUaj0Wg0NdDCaYsjpSSeSm92MzQajUaj0Wg024RQIrTZTbgv0cJpi/Mrr17h\n8V/+FtPL0c1uikaj0Wg0Go1mC3Jr+RbfuPMNABaji7w5+SYToYlNbtX9hxZOW5hYMs0X3rxNIJrk\nm1fmNrs5Go1Go9FoNJotyFXfVVKZFAAriRUAgvHgZjbpvkQLpy3MVz+Yyb2+MqM7v0aj0Wg0Gs12\nJpKMcNV3tan7EEI0dfs7GS2ctiiXpgP83J9+wL4+D0fHOri1EN7sJmk0Go1Gs2ZuBW7x6q1XN7sZ\nGs2mcHbuLLeWb23IHCSJbPo+dhpaOG1R/uTsJGkp+d3PnGBvr4dbi1o4aTQajWZrkZEZYqlYQ+tc\nXVKj7VJqo06z89BiZnujhdMW5dJ0kMPD7Yx0tjLW1cpCKE4sqbPraTQajWbrcHHxIn99769JZ+r/\nfbIIZXpkZKZZzdqRhBIhfU63EYLmhdM1c9s7HS2ctijTy1FGu1oB6G9vAWA+GN/MJmk0Go1mG5DM\nJJkNz5Z9ns6kuea71pDIqcV8ZB6oLYIWo4t1eZiSmeS6tm+nEE/HeXPyTT5Y/GCzm6KpwUZ6Wo19\nbcZ9FU6GWYwubug+NwItnLYgmYxkNhBjMCuYjP8zAZ2SXKPRaDTVubR4iXNz58rmUNwN3eXm8k1u\nBW5taHt8MR+nZ07zof/DmsuevHOSt6bfqvz93ZOcmT2zns27LzCyqflj/lWtn5EZTk2eyglhTX0s\nRhe5E7iz2c2oyck7J3lz6s0N3ecbE29weub0hu5zI9DCaQsSiqdIZSTdHieQF06zwcbiyDUajUaz\n84im1CBbIp0o+twYfW5GOFe1eRvxlIqWiKQidW3LSKVsRjKd1MZ9E0ikE6wkVri4cHFT9j8bni3r\nr7WWv7l8s4ktUixGF7m8dLni96dnTlf9vhrNzHxnFqoXSdZ3/2mqo4XTFiQYTQLQ1mIDYKDdBcBs\nQAsnjUajKUUI4RJC7N/sdmwVrMIKlAskw1DbrHkwMyszZGRmS02ONzL8bZfwwPfm3zMNwzRY7bk1\n+kxKpla1/lpIpBOcmzvHmbn6PYnn5s5xzXetia1SnJ45ve4epXBy45J9baV77X6hacJJCPFfhRDz\nQogPCj7rEkKcFEJcz/7vbNb+tzMBQzi57AB4nDa8ThszWjhpNBpNEUKITwLnga9l3z8qhPiLzW3V\n5mIIpFKjyRiF3mhjqnB/adlcgZJIJzg1eapu4/R24Dag5oBsdd5feJ+ZlRnOzZ1r2j42I9OhIeQN\nT6lmHdC5IZpGMz1Ovwd8vOSznwW+JaXcB3wr+15TguFxas8KJ4CB9hbtcdJoNJpyfgF4HFgGkFKe\nB3ZvZoMaYSm61DRjdT22m8wk1+ShklJyJ3hnze2ol9nwLCuJFW4tb+w8ro1gKjS12U1YV6KpKOlM\nOies46k4786+u8mtuk/Ywo6mcDLMq7dezXlOLy9d5srSlYrL3wncwRfzbVTzatI04SSlPAWUHukP\nAL+fff37wP/UrP1vZ4IxI1SvWDjN6DlOGo1GU0pSShko+ayq2SCEaBFCnBZCvC+EuCSE+MUmtq8i\ni9FF3pl5Z01zNa4sXSmb81PJs7QaT9PJOyf52u2vrbp9M+EZlmPLpt9ttOdramWKidBEaSOawq3A\nraYmDbgf5nm9du81Ts8WJw9YiCxUXScQD6xL4djF6OKmz/nZqJTh37r3LdPP/TE/l5YubUgbClmO\nq+eBIZzuBO7kPL9mXF66zHenv7shbauHjZ7j1C+lnMm+ngX6Ky0ohPisEOKMEOLMwkL1G+l+IxhV\nMcbtrQXCqa2FWZ1VT6PRaEq5JIT4YcAqhNgnhPgN4Ds11okDH5VSHgEeBT4uhHiy2Q0FZfgZIUnx\ntEqasJKsnAyhkIzMcH7+fJHBdztwuyzLnFEnqVEWIgsNG5OFxt+t5Vu8MfFG0fdGtjeDufDcqtq2\nmvaU8v78+xuW/ODq0tW6kwZEkhEC8VLtX2P7vquradaWo94sgHPhOaZXpnlr6i3enFx7drjTM6d5\nfeL1upc3Qj8jyUjOA7tZiTQaxUjOUsrb029zN3B3g1uz/dm05BBSxRBUHOuRUn5eSnlcSnm8t7d3\nA1u2+QRKkkOAyqy3EIqTSuvidhqNRlPATwIHUWLoj4Eg8FPVVpAKQ63Ys39Nd3/MR+Z5a+otXrv3\nmtGOhtZfiCwwvTJtapCnM2nOzJ4hkozk5ziVbH8yNAlUFhfvzr676pTF1/3Xueq7Wja3qHRfFxYu\nVDzuRsXDWggnw1WTLKyGa75rq/JkvD7xOm9NmadgDyfDvDHxxpq9I1MrU3WJrUpewPWaA3Zr+RZf\nvf3Vhtc7O3eW8/Pn16UNqyGVSZHMJHl94nUuLl5ESlnuvcySSCeK5mvdDtzmbnDtAuXVW682JZRx\nM+a1bWc2WjjNCSEGAbL/t7+vuQkEY0ksAtyOvHAaaHeRkbCwoovgajQajYGUMiKl/FdSyhPZwbZ/\nJaWsGdcshLAKIc6jfodOSinfaXZb16v+kJlxOx+ZZz4yz1XfVdNQvaXoUlma7/Pz58uM2NVmlzMz\nDJdjy8TS9YWYL0YXK4qH1TARmqiade3U5KmiJAsZmeHd2XeLztFEcIKl6FLd+7y5fJPTM6dN52Mk\n0olVhWTeCd4hnAybhubNhmd5e/rtos8qieL3599f9byve8F7nLxzsmrCjXoLrF71XW26od6srJHG\n8dWaa/bNu9/MDY6ACqe9tNh4SJyUskyw1gplbGTbBuuZlOPt6be3RV2rtbDRwukvgB/Lvv4x4Msb\nvP9tQSCapM1lx2LJPwDzRXD1PCeNRqMxEEK8JoT469K/WutJKdNSykeBEeBxIcQhk203LWR8I2rQ\n1MpmN70yjZSSV2+9Wldx2kb5zvR3uO6/XteypaIumAhWXb4e70m1c1xqvF/zX2MhslA05+Pi4kXe\nmVkfPX1h4QLXfNdqTnIvnThveJpK5/WsJFZ4b/49/DF/0bHUmjd2ceGiqbAoXS+ZTrIYXQRgLjKX\n22clTt45ybenvl11380knUnzzsw7LEWX+Nrtr/Hm5Ju5UNiN2n+182MsE0vlbbha89TOL5zn5J2T\nZZ/fDd5tSNCbUXhuXp94fVWCLJKMlNXe8sf8dYeoxtPxbentqks4CSEeaXTDQog/Bt4G9gshJoUQ\nPw78KvCyEOI68FL2vaaEYDRZlBgCVHII0LWcNBqNpoSfBv5F9u9fo1KT1+3akVIuA69RngW2qSHj\njdSgqSeErbCY5vTKNEBVD0NpmFq1ydmghEahkWPM9WikaGk1MhQb88bcqEohamvJmmcmHGZWZkyW\nNEdKWdWrYSZeDM9BLUOxcOL8veC9nEFrFhZWbVv3gvfwxXyEEqGi5SZC1b1oUkqu+q5y8u5JTs+c\n5u3pt+tOY294pJKZZJHQk1Lii/kqhvsVemfMqOaJNObMLcWWWIoucW5eeRFDiRBn585W3S6o/jUR\nVOc2lUlxw38DKSWnJk+VnfNqBWvfm3+PU5OnqvaLr9/5On99Lz+mc3npctXwu0p98tLipZqCfimm\nrnHpgMkN/w1iqVjZsawmBPD1idd5baL6tavaxugS15frG1jZSthqLwLAfxZCOFEpxv+bSQajMqSU\nP1Thqxfr3OeOJRhL0eYqvjQDbdrjpNFoNKVIKUuto7eEEKdNF84ihOhFZeNbFkK4gJeBX2tWG9fC\nXHguZwB2u7pzn5+fP28qqAqN21AixGJ0kR5XT9Ey0yvTVYXSq7de5Xt3fS9WizX32Vdvf5UuVxdP\nDj7JbHiWc3PnaLW3rvq4Srm6VO5BCsQDvDX1Fge6D7C7fTfvzb9X9H0inUAgsFsLBhor2LaF4mGt\nab3PzJ1hIbLAK3teMf2+3qQBK4kV3HZ30WeFoXb+eH2JE0LJEN+eLPb2fLCYK6HJ3o69ZeukM2m+\ncfcbPNr7KIOewaLvCkWpP+anzdkGKHGdzCQZ9Y7mvk+kE2VetHdm3iEYD2Kz2Hhh9AXuhe5xzXeN\ndmd7XceTSCdwWB2599UGDoKJIP1u8zxjlUR9oYAzEkR0ubq4E7zD3cBdWu2trCRWiq5jrfA/Q6is\nJkzQzNvTiLfMF/OV7dcYPCn1VH7o/5D56DzH+o813M5CDE9yYXjmagr7LkYWK363klgpegZtFery\nOEkpnwP+F2AUOCuE+CMhxMtNbdkOJhBNFtVwAuhoteOwWlgI6TlOGo1GY5AtrG789QghvheoZaEN\nAq8JIS4A76LmOH2l6Y2twjXfNb4zVZ4MMJzKGyOFnoLplem6DBXDsCk0yM2MslLDy2xuki+qDGRj\nJH6tCQtqeTCM7RuZ10pH4L9595ucvFseymTGNX/ew3dxMW8Ql56LpehSTS9aqaErpcyFtVWidIR/\nJbHCqclT6xIiWSqaSjG8WIXMRmaRUvKh/0OS6WTVzI7BuAqb9Mf8ZaLwO9PfKZorNhGayC2fyqRY\nji/nQv0K59LUCmtrFLO5XZFkhHdm3inr22bZ9OYicyTTSlCZiZ8byzeKt50y7/uFfTqRTqw63bdx\nDuvhu9Pf5fRM1bGiIjIy03Aq9IzM5Pp4KpMyDcGtdQ80gi/m49TkqZreyEA8sCrBthbq9Tghpbwu\nhPg/UCEQvw48JtST4OeklH/arAbuRILRJH1eT9FnQgi6PQ4WdXIIjUajKeQsKiOeAFLAbeDHq60g\npbwAPNb8plWn0AvSrDlP9dZKWu+5Buu5vdnw7Jqzuhmirx7Oz5/n8cHHc++llNwO3Ga8bbxoBHwx\nukh3SzfXl69zw3/DbFNF2wBlzNssttzxlHqVDKG7FF0q8jDWSz1iNhAP5ARbRmZ4a/qtVYvg0vVK\nhVVGZnJ1vAoFaaWMdAb+mJ+3p9/muZHnqi5n9O/psPKwlPa7pegSK8kV2hxtuc8MgVRIIp3IeWnM\nhFM0FS0SC6VCyiwxxuWly7ltNkq1sEADX8xHV0tXw9sOxoMNCyejltszw8/UDOutRCQZ4U7wTtG1\nKORW4BZ72vcA5t6rYCJYtO5EcCI3CFLJ+9sM6hJOQojDwGeA7wNOAp+UUp4TQgyh5jFp4bSOBEzm\nOAF0exwsaeGk0Wg0OaSUuze7DavF+NGvJgrM6qw0MpF7tQImnAyXhZGBMtbqycJVj2CbCE3kDKVS\nAvFA0XHWm+1vPYqKBhNB7gXv5d5Ph6e56rtKPB3nQPeB3OenZ07T0dKBy+aquU3Da2Zc8ycGnzBd\nLpRUoVXxdHzdsp2VhtIVerkkkmiyeTUiV1Pk+Ozc2Vx4YL1ejGrhl4vRRdocbcRSMdMkKVB8TxWG\nORZS2CcqhQ+utU7Z7cBthj3DZf3YbE7fd6e/y5ND5eXnVuuBiSQjNcNvk5lkw9u/snSFRCaBP+Yn\nkoww3jZuutzVpauMekZZiJo/33xRX044RVPRIs/xRlKvx+k3gN9BeZdyd5iUcjrrhdKsE1JKlqNJ\nOtzlwqnH42RxZX0m4mo0Gs12RgjxN6t9v50iIQqNtpXECgvRBXa3Kz24VuN5NYYrqNTphvFayHen\nv1vX+u8vvF9zmatLVysKpytLV+raTy3SmTRWi7Uhj1UinSgyng3xGU/Hy8L4lmPLuDyVhVMyk8Ru\nKf89N8hkir0bhder3uKwtagmnlfbP9ZKNePbH/PXPR+qUthYIcuxZWinKDFDKaVzgUopLSpdKdTw\nwsKF3OtG5zv5Y36uLF1hKbpUllWyUhZJs+K29YTtmYUavj7xOg90PsC+jn0sRBfoa+0rWyYQDzRc\nb63UQ1UoXoPJ4uO8sXyD24Hb9LaWJ+MJJoIsRhdJppNl8x03knqF0/cBUSnV0QohLEBLtn7GF5vW\nuh1IJJEmkcrQ1eoo+67H4+TD2eo3t0aj0ewQPlnlO8k2jYQ4NXkKUJ6TXe27Gl6/1EiWSN6aemtV\nxWUbmWdRSr1Z6r499W2eHX521fupxdfvfJ0nh55c0zyIwrApMwO22rGevHOSl8crTwkvFS7tjvoE\nQzUaydhoFrZWL4UemEpUEg+1UnHXS2nNoGYJQatoPElBIwWWLyxcyAmV0nNTKwy0lHoEW6VMhTf8\nN3BYHFxeuozb7ubpoaeLkq+Y9S0pZV2hhWbtKx04MOYcGsk2CpkMTTIZmiwTsRtNvcLpm6j04YbE\nbgW+ATzdjEbtZHxhNZrV6S4XTmqOU6LhTqrRaDT3G1LKz2x2Gxqh0SQAl5cur044lRiOGZlZlWha\nC42kNl6LODO45rvGg50PIoQw/W1cay2mtYb/mSWwMGrdlF4vp9W5pn1BY/Pl1lIstlJIWyGVQuOa\nxWoLONdivbx/lZgMTVbcR7Vnh5lQXGv9KsMbFU6GOXn3JC+Nv1R1+Su+Kw0Vva2nz5UKqkbXbyb1\nCqcWKWXOLymlXBFCrF8eUk0Of0QJJzOPU6/HSSKdIRRPmc6B0mg0mp2IEOL7gINAi/GZlPLfbl6L\nyml01BhUoct6MYz71RSyXG8abUM9HhIjM5sZN5dvshBdIBgPmoZ4rTVRhRFqtJ7ejFqhYVuZS4uX\nONhzsKn7MAzxRrxn252tMiBeWuerVmbARkQTVPfErTaZxkZSr78rLIQ4arwRQhwDmjebcAdTy+ME\nsKhTkms0Gg0AQojfAj4N/CQqs96nAPPZx9uMRozr0gQABmsJxdoo6vGQXFqsbrwZnqtmeNcKt7le\nCRsMCj1uhWm9tzKNCPq1eoDWw7sgkes2Z26rUS2N/Gopfe5shcGYrUS9HqefAv6HEGIa9cM0gPqh\n0qwzOY+TiXDq8SgX/lI4wZ71LWKv0Wg025WnpZSHhRAXpJS/KIT4HPDVzW7UetCIUVQpRXClSeWa\n1VE6Gr+ezIZnGWsba9r215N605cbIYmbyVoz3W1lVuPJbpRUJtX0fWwn6hJOUsp3hRAPAfuzH12T\nUm79YaxtiC+sTmul5BCgPU4ajUZTgOECiGRLZCyhCtxuexqpPaRpPhsxV2y1NXI2mguLF2ovpGmI\n9S4K3Kxt7nTqLoALnAB2Zdc5KoRASvkHTWnVDsYfTmC1CLwt5ZcmF6qnazlpNBqNwVeEEB3AfwTO\noTLqfWFzm6S5H1ltkdhG2C5hUVrUa6rRSEbB7Ua9BXC/COwFzgNGwKoEtHBaZ3yRBJ2tdiyW8kmC\nXa0OhIAFXctJo9FoAJBS/lL25ZeEEF9BJTPa2DRyGo1Go8mxHeZWrpZ6PU7HgYflWlPTaGriDyfo\nNAnTA7BZLXS47PjC2uOk0Wg0AEKIC8D/C/x/UsqbgH5AajQazSYyEZrY7CY0jXqz6n2ASgihaTK+\ncMI0o55Bl9uRy7yn0Wg0Gj4JpID/LoR4Vwjx00KI7THDXqPRaO5DluPLm92EplGvcOoBLgshvi6E\n+Avjr5kN26n4IwnTxBAG3W4nSzpUT6PRaACQUt6VUv4HKeUx4IeBw8D2mGGv0Wg0mm1FvaF6v9DM\nRmjy+MJJjo1X9zjdXNBZUjQajcZACDGOKpHxadQ83J/Z3BZpNBqN5n6kLo+TlPIN4A5gz75+F5W9\nSLOOSCmVx8ltr7hMl0eH6mk0Go2BEOId4M9Qv2efklI+LqX83CY3a/2JLIKeZrw+JCMQnN7sVtRP\nJg2Bic27/skwZLbIZP9UDIJTzdm2zECl9N2xZdgS9YwkLN+FtLYDN4u6hJMQ4ieAPwF+O/vRMPDn\nzWrUTiUYS5HOyIrJIQC63Q78kQSZjP4B1Wg0GuBHpZRHpZS/KqW8tdmNaQpRPyzdhODkZrfk/mDu\nshIiqzWEZUZdj2rG68qcWs6MTFKJoUKSEWUQmxG4p4ReZHF17a1FJpUVkoV2hQTfbSWaZj+AuUu1\nt5NOQmgWlu9BMlp7+arbSqhtZVIw8Q6Es8e++CEEJiHVhBwwSzfVcWZS+euTjKh9LVxT+66FlOUC\nNzS3fu2Nh9R58W3Eo05CYGqLCMatQ71znP4R8AwQBJBSXgf6mtWonYo/60nqqpEcIiNhObpFRn80\nGo1mE5FSXtvsNqyJZATmryhDG5TxGZigyIg1RvvTSVi6rgzJakR9eUMzFlDvy5Bq36ttczyoXqfi\nysiMl4zUZ1IFhn6Ngb5UNGsISmUwh+fz34UX8/sySKwog67UQE2sKOMXqdpVaKymohDO1kiSWaN4\n6my5gDFDyuw5zO4v6lPHVknorMyB/44ytkvbDjB1Dqbfy79PRmD2YlYomPy2Z7ICLFyjxlMqqvpS\nPcdUiP+u6nPRgBJ76QSkstdhIXt7FZ3LePYY76o2GyzdUOckNAOzF6h53athbCuWTTJg7Mc4F5iI\n0nQSFq+ra1NPnanSc50Iqf/JGEydyR7HRZj7QH0eD9Xe5uwFmDxd3KblO7B4TfXN0Iy6TqVMn1f7\nKk3jnYxQdB7NvI6VRFnUZ97/lm7mz6tBPKieK4XnJBpQgzVTZ+sTwsFp8/2ZtStQknVPZvLPCykh\nWiW5RGBSeR6bIZ7roN45TnEpZUIIVVtICGFjTXeExgxfRAmnWln1AHzheFWBpdFoNJotSCYNmQTY\nXMpImb2oPo8HwdMPvptK7LR0KKPJYoPlrJGRihYLlNCcMsp696vlM2nlIVi8rr5398DCVfXa1Qlt\nQ+DwqPfBaWWA9D4Ejlaw2JXRlo6DsCjDNRmF0Sey+46BsKq2z2YNydaeYi/IwCGwu9V2jOOwu5Ux\n6elX7XN6lJG0fA/cfdC1G2YuqPULhUFGgt2ptgMw8Ei+HcbxhefVPi3Z8PbFD9W+C9tktH/2ojLI\nQiWFOVNRdU5W5tQ5GXqMnKgUNrBYleEI0DEO3oG8MIn4IPKOOq/to/nt+e/kr+l8MN+GQmQ6L5hK\nCU5DazfYnOoaG0ZrPKQM7+AUdD2gjPHBR8HqUG1enlD7jAfA1ZU9J9dVH+ranT2vafDfVte4a0++\nLUbbF6bVfnoeVJ+VGvKZFMxfKv7cO5D/rpDQLHgGlJDo3AWtXbB4QwnckRNq+VRMnf+lG+p970P5\ndgKkjbbFVPvSWWN54UPVJybfVe9Hn1DnLerLi6aBQ5BKQks7iOz5c7Zl2zaj+uDgYXUvFpKOqf+B\nyfLjyqTUfYBUfcY7oAx4i131lVQsv2xkSZ1nY73IIkRQ+x04DPbsfuMr6rjSwPQ56DuQ7dNS9Q+r\nHQaOqO0bxAJqsMHpVf/tLvUMcHVALAjtI/n7ZOgxtQ2y9UEji+qvsF8ay/puQ8coBKbBVmBjrsyD\nwwXO9qzHNHserDZ1/pLRvBgythv1q3tSCOg/pJYvfD5hAadb9c+Fq+qaL99T/T6+ovpCxAf2Vmjt\nVP08FVf93wjX7H1IXd8NpF7h9IYQ4ucAlxDiZeAfAn/ZvGbtTHIepxpZ9UAlkdBoNBrNNiEZhYWs\nN0BmlBFU6hmIBSARzr6ReQPcoFA0zX2QX3bhWrmIgbwXC5QRE/XnDbbcullhNfqEMtpKCc+rZVfm\ny78r3d/Ch+Xha8ZovdGWwnXC85XDz5bvFL83ExjphPLegDKSzYpuLl6Hnn35kfpSL9vcJejemz/X\nwWm13cJzl2vTXXMvU3AabC1KgMZMRtwn3lHXu/8gRYE+Zsfku6VG2wMT6lqVnofle9njynqCZs6X\nb2PxOnSMqf5liAhXuxInhR4Ld4/yUJVuG8zD0mY/UIav2fF5B8rP7fK9/Db9d4r7cyWv6eKH4O7N\nb8s4fpnOC2xQAsUQTaC8KJaSICpD4FtsSjgko9D3sDK+jXZF/BC4ULyeIXbMQi0NEW3cb8JSfp+C\nav/Sjfz70r7pv63EjdWphGghxjUxhFU6qTxgo4/nRS6ogQbDe5WMqr/QjHrfNpRfrtC7WYjZNTCe\nE6WsZAeCYseqAAAgAElEQVQchKg91y6TVs864xkjswJQWIrPqVnocTqZP1fByfwzz+gHTm/x8gtX\nzQcmmki9wulngR8HLgL/AHgV+J1mNWqn4qsjVK8zmzhCF8HVaDQaEEK0Av8cGJNS/oQQYh+wX0r5\nlSrrjAJ/APSjoic+L6X8v5va0NkS48wsnMoQMVBs0JqRKDFgzQSImUE3e0GN/JcauZUMWV8Dmd3N\n5vxUmudT7/f1YhjJpUR9tUMbl27mXwcmikf266XWnJNMieFficIQpdI+0wiFIggKRvkLqNXHSjET\nTQalnrzVUslwr0W1+V+ZVN5rNH+5+LvSkDFj+Xr3Z3aPgbkoLiQeqn3+S8PjJk6XL1MaImvQrPlw\n9SQomTpTYd0G73WzYzMLl6yU0KNJ1CWcpJQZ4AvZP02T8NcRqmd4nJZ0Zj2NRqMB+F3gLPBU9v0U\n8D+AisIJVTD3n0spzwkhvMBZIcRJKeXlKuvcP0TqmP+xk2l0jpDm/qKRwYKtyv1wDPUydwke2rjd\n1SWchBC3MZnTJKXcs+4t2sH4wkkcVgtuR+XRrpzHSRfB1Wg0GoC9UspPCyF+CEBKGRHGhNwKSCln\ngJns65AQ4goqW+zOEE4ajUajWRX1huodL3jdAnwK6Fr/5uxs/OEEnW471X7znTYrXqdNe5w0Go1G\nkRBCuMgO7gkh9gJ1xzILIXYBjwE14rk0Go1GsyVJJ7MJMJpPvQVwlwr+pqSU/xfwfU1u247DF0lU\nreFkoIvgajQaTY6fB74GjAoh/hvwLeBn6llRCOEBvgT8lJSybFa/EOKzQogzQogzCws1UkFrNBqN\n5r6n3lC9owVvLSgPVL3eKk2d+MOJulKMd7m1cNJoNBoAKeVJIcQ54ElUvt1/KqWsOTNaCGFHiab/\nJqX80wrb/jzweYDjx4/rEhwajUazw6lX/Hyu4HUKuAP87dXuVAhxBwihstanpJTHq6+xM/BFEhwY\nbKu5XLfbwdRyrOZyGo1Gc79SMqAH2TlLwJgQYkxKaZJbO7euAP4LcEVK+X82q42FtFldBNN1FJHU\naDQaTWPUk+1vnag3q94LTdj3C/WMCu4k/OFE1RpOBl1uBxenAhvQIo1Go9myfK7KdxL4aJXvnwF+\nBLgohDAK4fyclPLV9WpcKXuGn+D8vdebtfnGGX1CpeE2S1vs6TOv27Rd6dlnno7b3atSqMf072kZ\nTq956uetiq2luPjshuzTqWpCrRarwzyFv6Yxevdv2PwmqD9U759V+36jRuzuZ9IZyXI0WTUVuUGX\n24kvnEBKWTWRhEaj0dyvrGVAT0r5bVRY38bRuQfWKpysdjUJumsvOD2AyBdA7doDkSUlAkafUAbd\nzHlwdVaui9O9Ny+cBg8r41NKVaiyUDgNPFJcl6bv4fJ6OKCKrra0gf9u/UZ3+wgETAphmtG7XxX7\nLWT0idp1mlxdqtBousTIbR/JGq/JfPFfd686h+mEOk6nVxUVNeoiDR+H6FLj6Z6796paR9371HUR\nVhg+WlzE1aC0re0j6homwuaFjgtpG4H2YVUYVdjMixobdIzlj8voM4vXwNkOnt7yWkTDx/IFYLv3\n5utfWWxKbAeni5dvH1EFgeMmRYG7ducLK7v7VDFkVyd0P5A/Jx3jqqZScKq8YLSrU50nozDr4JFs\nodUL0DGqrnnhdSul74CqwWR15IVz4Xl3dRTX1CrEO6DaBsV9zzgOY5nQLLS0VxbmfQ+r2m2Fgq/v\ngKrvVGn/bUP58zz6uBoiWrhSfL+1DYHFAalI8X3s6VftMYob9zyo7nV7K8y8V+61qXcApfchcHhU\nv2xpV9vJpCCTNB+wqIW9VdWZc7ihpUNdf1B939VRXFgY1DIbaAs3klXvBPAX2fefBE4DqzgjgLrU\n3xBCSOC3s3HkO5pANImU0NVaWzV3ux0k05JQPEVby8apbE2D3PgmnP8jCM7A7ufgxN9XDyKNRrNu\nCCFagH8IPIv6bXkT+C0p5ZaKZxY2Bww+mhc6A4eUoepsh6XrythxthUbYp6BvGEIMFQanQj0H1SG\nkLtXGRYyW7zT5lTGsJSwfEcZsf67Slz1HyzY5mNqGZsz29CsATL4qDIindnw8dEnYO4DcHiVmOja\no4q+du5S39vdWTGHMv6iflXA0ulRgiy8qIzKTFoJP/9tte9MWgmn1q7i+lIOtzKsPf35fYAyFpMx\nVYzV7lafde5SRlvvw2qdeEAt0zmeX6/vIWXAu/uUmEgn899Z7UocRpagfVR9VpilyzuoDOFMRhnw\ntpbi9oRm1bajy+ZC1epU16a1J38uDQwRbGtRhmZoVl0rBPkqMALahtX3Fjt4+yGdUkakQdSnjNz2\n4ez1d6n/7SPqXEQWVZ9buqWM0q7d6lx4B/PbsDlh4LB6bVas1GJT/UNKcHUDWeE0fEz9d3WqPpGM\nKkHSNqyWNYRTxy7VF4VQ+3aj+q3Do9pjMHxUCUthgVhWPDg8qo9k0pCIqGtrsanrYfQ7IZSAMnBk\nP+8YU9dy6aYSCxa7WmfgsDqnC9eUoT94RBWZtdjUclKWC9vuvfnrWIq7W12HTEoJq/Yx1aaoTwmI\njnFw9+TFp82pBiyW7+ULCDvb1L1ntStBJayqPUZR2dYuJZzsLkCobtL3sLqmcx9A/yElPEAJNkP4\nONuK7yPjeuVed6s+0rVHneMWr7q/kpF8IdqWdnX/R/35YsjCoj4Hda+W0ndALSOs+YLOhvjrOwjz\nl9RnxrPR1anuiXgo378N4dS9V/1PxfLe4mw/30hHQr3CaQQ4KqUMAQghfgH4Kynl31nlfp+VUk4J\nIfqAk0KIq1LKU4ULCCE+C3wWYGxsbJW72T4YyR7q8zipZXwrCS2ctiKxIPzZ/wrX/kr9OHSMwhv/\nAd7+f+CZfwpP/SNlFGg0mvXgD1BzZn8j+/6HgS+iymZsLQxxAsooMQz/3oLqje5eCGcz+LV488LJ\n6TXfpsOTF1RCgCj5TRACOrNGafcD6q8Qa4XfHJuzuL2gjLLCdrp7zddFqBF/V0HVko4xZYTFAsog\ncp9Qn1tsyvBzuJURn4pDa7cynJau54VM4bbtrqzhmMXTnzfanJ68IV10PC3FgsfYloG9Fdpb8+/L\nQn+EMtIL13N1qdfewWIBgoTQvBookymqJjB2FxjhFnvJ8ZYYgpZsmxwmx1d6vg3askLKMDoHHqnc\nlkKEJS+8J0/njfH+gxBeUv2qa3exl8LhUf1r9mJ5W3seVEaxo7X498/sWCwF597hVf2wfcR82faR\nysfg9CqD3OjHZoLHYi/u1yMnjBNQ7MXo3a88G6W4upRgMs790FFyglcU9JOBQ/n7HQruLaEEVeH9\nbXxnL+iPoO4b49wUih5j2ZHHyz8DFapa2jecJfPpPX1KODnbiu/7voNKCPrvKuFlsebFeWKl8vPD\nbD+Djyqh6vAo28jmVH3B4cmLdptT9b3CQYHBI8V9wujTtfbdJOoVTv1AYSBmIvvZqpBSTmX/zwsh\n/gx4HDhVssyOymbkj6jTW1dWPY9aZimcYFePNsC3FMkYfPEHYfo9eOkX4cn/TT0IFq/Dt/4tvPbL\ncPb34Ef+TD2INRrNWjkkpXy44P1rQoitW8i2a48yHqt931VQW37wUWUg3A9R2d17IR5WYqkQw2h0\nePLGsdNr7mFbD7r2KC+XdZXJgZ1eJQSrCUdv1kQqFbLbDSFUeKLIij+7GzqydofbJILC7lKeUuP4\nXR3F4qXSAEAlLFa1/mopFf+1ECUit21Y3a9mogmUKClaX2B6sxaKpoFHyo1+M9FbSKGXcuix+ub0\nWO3F6+XWP1owCJDF6TVfFtT9agjvQsyEbDVsTsBZ8Jpi0Vo4GFK0XumARzkSidigh2RddZxQI3qn\nhRC/kPU2vQP8/mp2KIRwCyG8xmvgY8AHq9nW/UTO41RPcojsMjol+Rbkaz+rXOqf+j149qfyD4ee\nffDpL8Lf+7oKAfmjT+fd3xqNZi2cE0I8abwRQjwBnNnE9pjiNJ4F7t5iI6oWNmdlY2y7YbEXjyRv\nFi0dWYNtDefUO1guAO9XLNYG5pAIFUZXaOw2Kl62Eu0jtUVNo9hb19Z3rA7W1Het9nKBqKmbegvg\n/jLwGcCf/fuMlPJXVrnPfuDbQoj3UfOk/kpK+bVVbuu+YSGkJiT2ems/YHKheuE1ZHPRrD8zF5Q3\n6cl/BA9/v/kyY0/C3/59WL4LJ//NhjZPo7lPOQZ8RwhxJ1vq4m3ghBDiohDiwuY2LU+7o32zm7Bh\nvDDWjES8Go1GY47caunIs7QCQSnl7woheoUQu6WUDaaVASnlLeBIzQV3GPOhOEKoxA+16C4I1dNs\nIf7636nR1O/5merLjT8Nx38czvxXFcpX6urXaDSN8PHNboBGo9FoNg/JxgmnujxOQoifB/4l8L9n\nP7IDf9isRu1EFkIxut1ObNbal6TVYaPFbsG3ooXTlsF3C65/Ax7/B/WFonzPv1TxvG/8WvPbptHc\nx0gp7wJBoB3oNv6klHez320JdlLpiI2aa7BRPNT9UO2FNDuWj4x+ZLOboNlA6g1y/EHg+4EwgJRy\nGmhwhp+mGnPBOH11hOkZdGdrOWm2CO/+FxUHfuzv1re8pxce+xG49GcqXblGo1kVQohfAi4Av44q\nivs54D9taqM09xXDnuGG19lJQnkncXzgeNlnraXZ7zQbzpbzOAEJqQIIJeSSOmjWkflQjP62+oVT\nl9uhQ/W2CokIvPeH8NDfgLbB2ssbPP4Tql7C2d9tXts0mvufvw3slVJ+REr5Qvbvo5vdqO2O2+7m\n5fGXay63u313zWW2K12uLl7Z8wqWuk2lPC5bhQxh2windRsndWgSdss2z5KoWTP1Pg3+uxDit4EO\nIcRPAN8EvtC8Zu085oNx+ry1Uy4adLkd2uO0Vbj0p6pI3+OfbWy97r3w4PequU4pfS01mlXyAbAF\nUrVVZzuGr9nrSHlsM8kOZnhbOls6y77blmy/S7cujHir1EjaoWgx2Ry6Xd2rXrejpWNDBW29WfX+\nE/AnwJeA/cC/kVL+RvW1NPWSzkgWV+L0NeBx6tbCaWsgJZz+AvQeUEkfGuXE31fFLq/91fq3TaPZ\nGfx74D0hxNeFEH9h/G12o3YKHU5zzfrS+Es8MVhcF+bJoSdNl20WHx1bm+Nxq4td5yrSfDdi+JsZ\no3s7TOr5VOGJwSd4buS54u3WU4NoHWmxtWxpD+CL4y9udhNqcqD7QF3LHR84zpBnqKFtu2wuxtvG\nK37ftd7p4NdITeEkhLAKIV6TUp6UUv4LKeVPSylPbkTjdgpLK3EyEvraGvM4Lel05JvP1DmYOQ8n\nfryBOhcF7P0otI8pr5NGo1kNvw/8GvCr5Oc4fW5TW1QDyzaooVLvnAEzg9QqrDisjrLjNPNO1eLp\noVUMSAEWi4WWOgpnNrRNYWF/V75wuadKAdBkJrmqfQy4B+pedlfbroa3XypiqmE2T6vVVv98HiEE\n3a5uvI7iKfGHug9VWKM5tDna6gopfXbk2ZrLOEoL15rw1NBTdbXLwCbWvx5Yb2u+OPMjvY/UXL7W\noMZSdKnofaV5XQ6LY1Uitc3RZvr58YHjHO8vn1dWyEYPcNR8eksp00BGCLFzilBsMPPZGk6NJIfo\n8jiIJTNEEqlmNUtTD+/+jqqeffjTq1vfYoVjPwa3T8HijfVtm0azM4hIKX9dSvmalPIN42+zG7Xd\nMYyRQz2NG7mrEUigxE4hXoeXjpbVRWEe6DIfIa9mIB7tP8qejj2592bemcJje7j74YrbSqZXJ5we\n7Xu0aPR91DvKiYETDHoamD9bAbfdXZfhb2AmngsN8loevU/s/oTp5+txLI1Sj9fCY68shAFe2fNK\nXX271nZKWc9EIkf6yqv9lHoOX9nzStkyXS3Vz09pX6jkVRJCNJyoodryfa19q36eNIt6h71WgItC\niP8ihPh146+ZDdtJzAVjQGPCyaj3tKRTkm8e4SX44Etw5H+GFvPRkrp47EdUFXGdJEKjWQ1vCiH+\nvRDiKSHEUeNvsxtVSqFxVMtQqES94TJrwTCIjPaOtY2ZGloG65nN6uO78iW5nht5rszbZObhcdvN\nc1VV8upVGtkGZWA+1JVPPV74GgBRPLpdad9m1CsWLMLCwZ6DOS/NWNtYkVgppJ6in7XCpj6262NF\nHoIXx19kvD0r3Ew2X9iPq3n0aomBRr0ya8Vr99Lj6qm6zHp5LuxWO8+PPG96vx7sOWi6TrW+NNY2\nVvS+dIDBtA0V5vwc7V/do7GRArP1JozZrmnc6xVOfwr8a+AUcLbgT7MO5DxODYXqKZGl5zltIuf/\nENJxVcx2LXj7Yf8rcP6PIBlbn7ZpNDuHx4AngV9hm6QjrzTqv6t9V1lYUyFDnqEi8VBtXoAZr+x5\nhcf6Hsu9NwsLa3QOSyM0EpLkdXixWqxFnx3sLjc625xtZV6p/V37GfGUJzaoJEAqYSYM6k1NXjpK\n7rDU7+lZK4V95Ehv3gOxr7O82LrNYuOFsRdy751WJ1ahznuhKC68Fg91P8Qzw89UbcPzI89X/X4j\nE4e47W6EEEV9f70pFSQeh8dUQJjdswLBU0NPVTynpYJu3Fv5vje8XYUisVD0dLc0loThxMAJHut7\njHZn/UFn9cyj8zq8OXG9kanE14OqwkkIMQYgpfx9s7+NaeL9z3xQCadeT/0epx6PeggbokuzwWQy\nqnbT+DPQXzlco26O/z2I+uDKX659WxrNDqIgBfkL2yUd+VOD5qPtrbbWqkaE0+rkueHnit6vN622\nVjwOT9UwtELqnUf05NCTpvMinhp6Co/Dw6N9jwLK01RpronZ8QpEmWdqb8deU4/HiYETdbW1GoUC\nopYQfHLoyYrXqFQg7+/ab2o454xLkxH/eibhF56HeiftGwk/Cr1zhpgC2NO+J2dInxg4URQeZojY\n9fLelPZDm8XWcGibmWA0o3S7RhKLvR17i/qYWTjcWpNPOKwO2p3t9Lv717Sddmc7L42/VDEjYqNh\nb06bk0HPYJm3bnfbbvrd/RVLFtTy7j038lyujzRjjlczqeVx+nPjhRDiS01uy45lNhij2+3AYat/\nwvBwp7pJp5ejzWqWpho3vwXLd1VSiPVg9/dA526dJEKjWQVCiO8TQvyMEOLfGH+b3aZqVJpYLZGm\noWRPDj3Jx3d/3HT5StSTUU4IUWZgWS1Wnh95vqbhY1BvGuBKcyg6Wzp5fuT5nFHvdXgrhtM1e2S6\n1vZLxUCt7HBdLV1FHrzC5Q0j/Wj/UR4ffJy9HXurjuqbJmqo1I8aCKsyY8A9wAtjL9Tloett7S3y\nwh3rP8ahnkMNFYUd9Y5W/G5X+66i9y22Fh4feLzubUNB2GYdeqswDG6gdYCXx19mf9f+Iq9mu6N5\nU/6P9R8rm4/V5mxsKkC1eWzVRGc9AzEOq4MnBp/AbrVzrP9YxXugnkEKp9VJv7u/oXmUhclZDDa6\n2HQtS72wNXsqLqVZEzOBKIMdjWX/6XE7cVgtWjhtFqe/AJ5+eOiT67M9iwWO/V249x2Yv7o+29Ro\ndgBCiN8CPg38JOo361NAYzFsW4ghz1Bu9Nplc/HS+Et0tXSVzdkRQuQERmGYkDGiXOoJMhtpdtvd\nHO49XHUOU7PISLlmA9+MWhnEzMRROJ4inZG5874eAs0w5gxPjUVYeHGsPO20y+aqS6Q+3P0w4+3j\ndXsC18pqPShOq7NsTk4thr3Vwx+Nvn+k7whPDD5RV82fVnsrkUSKmUC5jWSxWMoEmSESC+fZgbk4\nrpYVs1qmRQOJJJZKs7gS58O5UJnhf6zvGM+OPMvTQ0/z1NBTjHpHa4ZGrgf1ZFx0291l5//F8Rdz\noss4lnrEjBCCY/3H6HZ102pvZXf77pohnqvJJLne1BJOssJrzToyG4gx0NbYQ8piEQx1tDCphdPG\n478D178BR38MbOsYt/7Y3wGLHc7+3vptU6O5/3laSvmjgF9K+YvAU8CDm9wmU4LRJB5rP4lUhmAs\nSTSZJpnOAJDMZHLLPTH4BKF4klRGlo0eJ9OSftc4J/qeot/dz9HeZxhwD5DKSIKxJE/2P8+B9qdY\njiQ4e9fPmbs+osk0Tw6oUCOJEiyRuGBP2wOm7Ywl06afpzMSXySBLxxHkm9bPJUmnZGkMpJ4qnhd\niSRTIJCWwnEiyTTn7vm5NBMs20c8lebiZIBbCyssRxK5/WakxOvw0unsZ8IfIZpMF5Xk8EUSLK7E\nGfGMsLgSJxRTWe0SKRuzhcazVG2Kp9JIJAIbV2aD3FoM5ybon73r5+5SpKhdZ+76mA2oOahttmHm\n/LZce4PRZO4446m02nb2mIe9w+zt2Mu+zn3I7PfGd1JKFkLF81pDsSTpTLG5JaXk/N0VBpwPlBn8\nC6E4Z+76sAoHiVSGlXiSCX+EM3d9XJ9fwR9O8NTQU2WC8vZiOHeM95YiLIXjRCtcd8iLhVKv24Qv\nwo35EMlMJtcHjfZLKYkm1DZ94TiXpgM5UXrmjo9wPIPL5qKrpYuXxz9OJJEikkgRjCUJRvNZCUda\nDzAfgCH3EHbhIJNRfTidkYTi6j4y+vn+zod5rO8YGZnh8kyQqQIbKZXOt+9A1wE+OvZR5kNxQvFk\nkYckmc4QSaRJpDMEouXZEVvtrRzpO0IynSEUV9/bLXZW4im6nMXeokd7j3O07wSJlDo/ANPLMT6Y\nCnBnKUwwVrz9TEZis9jw2Ly4bN7cXLB4wsniSly1K5Xh8nSQUDzJQiie224pqYwseq4Uks5IEqn8\nd7GkuoeT6QxL4Ti+cLzo/nJanVyeCZJOlXsSr8/GmF0uP0/RZDp3TdMZJRZX4qmKpXRccoy5gLqX\ng7Fk7tkIsBxNMBeMsRxJkcxkiKXUs3NxJc5yJEGmwjloBrUCC48IIYKoUTxX9jXZ91JKuYZUYhqD\n6eUoj+9uvMDXcKeLKb8WThvOmd8FkfUQrSfuHnj4++H9P4KXfh7sW7dgn0azhTAeghEhxBCwBGx8\nvuMaxJJpnImj+HyCr/pn+NAfKltmYXGO3hYHiXSMa8EQLlsaGZ4lnsqQkZJWhy1bgqKDybkw491w\ndykMrPC+3w+AK2HUW8kb/pemAzhiQV54qIXzE8ucnfXT4ejjry7OcGSkg1aHlTN3fXQ4+vhycAqA\n/QNeBtpauLmwQiyZYdB5jL+8UVDCcTGMx9ZJ2D/FRX8g97EFCz+Yla3X50KE/ftJyTjheAqHzcLt\nxXDROfny+Sn29XmJJlMMtrt4946v7LxEUnsQCN74cIHlyACh8CRzQbXPQHCJB9pT3FpYAeAv3p/O\nrbenx0Nq5QEmA99hcjlKMjiFlJI7YSvBpDpfY+6HgXlWYhmWIwnaWpSHYWElxoQvQqfbwUpcGXpT\n/hjfuDRLNNlDl7WHu0thLk4FMGO43cu3k4tZI9HD1XtzZGSGi8tqeZc1RTQdIuj147K10Wq3cn1+\nJech+TC4RDS9QtA/w76eFLPBGLPZDLzv+9U5emUPLAUdQJjz9wLEM2pdp1UZt4FoglPXF3JX5hzq\n2s5nojnj9cvn1We3/eq62KJTHB3r5Nw9P3t7PVycCvBgv4cx10Guhya4tRDj2qzaZrvLnhMW7/uX\nAZAr8wAcGm7ng4Jzcyu7/YBvF+dsfqaWo9jlowRi8GWfasNlvzIx+1rGmI/dU62OzBJP2egUR3jt\n2jyhmCrBksg8RCId4+bKe7l9XJoOYIvagRRXAvl+dGM+RIfLybdvzvPhcgiLsDJiV8LlXrZdH7QF\neKDPw9cvzeaOJRayMDk7n9tOr8fJgcG23Dk1ljvj8jG1HEWmH+L6pIubU1N0tjrwR/LJuy5klxXh\n6TIv2M2FFXq9Tt78cJGUidDp8SjRdCfb1gXfIpFkimuz6hky4YuSIUO/t4Xz9uXsM4HcM2HZN8vd\nsDofnXKJ+VAs993RsU4+mAqQDM6SyiS4FFgu2veriRlc1hD7+ryMOB9lccmTu14P9nvxhxMsrMSZ\nC6W56fOz4p/j4KDkxvwKl/wB3LZ2xtwP817g7aLtfvn8FH3eFuazAwdP7unm/Um1b7c8ysXlUwD8\nwD7jGqr7+9s3FnPnfbj1QaYiH7IYsPPRcYnDsjEhe1WFk5TSWu17zdoJx1MEY+oHo1GGO1y8fm2h\n9oKa9SMZg3N/APs/Ae31ZVdqiGOfUSnOL/0ZPPrD6799jeb+4ytCiA7gPwLnUNERX6i1khDivwJ/\nA5iXUja9IuepDxfqLnxrjMqnM5kiL0Bp3T7DQKqHjJR868oc/kSxh8MwVo50vlD0+bXZUM4wA1RR\nkgrs9hzm9soFAB7p/J6cMQ7gtLpw4uKbV+Yqrn99Xu1nssJAYKtNjdEaHqhWWzuRtFpHCMG3Kmz7\n1uIKdkvxvA0hBLs9j/C+/zUAPDY1om8IM4C9nsewWeycu6eMy1QmbwAXXo/zE8vsch/iTviDsn3P\nBWP0WIpH1islTDhjIhYLwwUrnZcvn59ipPVBrgbfMf2+Embha8Ot+2ixqDAz47hvLqzgsnRzcWoe\nS2cKC4NFfcLMG2PwQQVB2WptZ8KnRH2lcK5B114cFhcCUeTBNEQTgMPSgsNS3xSHy9MhhCjuwHdK\n7p2bCyvcXChexmsvHtBeWImzcL3c5jK8Wi3WfErxQtFUSMYkPLXSuTJYXKmeBMxl8xJOBZgLxbhr\nyx+X195FRhZ7EQ2hYhFWMjJNPPAwD1SJLjT67PX5EC5bcUKTD+fyfWHM/TDB5CJOa2tO5BxsfxaL\nsFZ87s0XeFu/eytfYNdSkIik8FlSSosl7/2ybZBogtoeJ02Tmcm6/gfbG69wPtThYj4UJ55K47Rp\njbshXP5zlf3u8Z9ozvZ3PQvd+5RXSwsnjaYmUspfyr78khDiK0CLlLK6JaL4PeA3gT9oVtsKqRQG\nZcFCBvNwmvXCzMBcr6xnAG32xlIcr5V2ew+L8UkAXNbK6dsbodCI99iL05tbslm/+lrKp87Va7yX\n0pTw7TgAABuqSURBVO0cYjJyDYdl/TMj1mLQtZeZ6M2iz3qc5lnYxlofZqT1IdPvCulw9LGcmK+5\nXL10O+vLAFhIX4v53KqNTh6wVdjjUdn/jOvSbs8n+9jf9jiJdAxbjTT5hWKwGnaLo+ya2QoSx3jt\nXYSS5QMEldjjOYy94N462P5MleekwLKBwqn+NG6apmCM/KxGOA13GJn1dO2fDePd31HCZvf3NGf7\nQsDxz8DkaZgtH8XUaDQKIcQJIcRAwfsfBf478EtCiJqxz1LKU0D9v+TrzEjrfvZ5jxV9ZhgpqxE1\n4+5DeGzFBr/htdjnPV62/GoSIDxosp3NwGPv5EjnCzzoPV6XgV1oMBbitrVTz/Rti7BwpPMFelsq\nZ38zaLFWziZXaMB3O4c40vlCTcO1GVQSGGYIIYpSkVdi3H2wzGu5EdhE3ji3kG+nWf82zr7TUn/G\nv2oMuR4ou+cq0W7vZdC1thpprVYvPU7zSJd93mPsb6uebbCw/zksLRUHCIx2Fp7btbLbfZjDHR+p\ne3mvvbtItNksjlUPUqw3WjhtMobHaaij8VC9vX3Kv3p9rjxWXtME7n4HJt+FE39fCZxmceSHwOqE\nd36refvQaLY/vw0kAIQQzwO/ivIeBYDPr8cOhBCfFUKcEUKcWVhYW1j0Jw8P0efN//B3O4dyIWig\nRoBLQ4MODbfjcdYXGPLo4C4e633c9LfEZs3/1Lfbe+lyDDLkMk8MUQ0jVMdhXbvp0O0cKmrDcEm7\nn9jdzaOjHXS0Onhs1LxYamHoUKsjm03QXm7k7/IcKjPqH25/Jjcir6jvmb63tySuSQhcdivWrNHp\ntXXR5XZU9XLYrRaOjNRncJuxv+1xRlpL0jKv4TfJmNdlxmqmEWwkD7adyIVbFtJuz2cqdNqstLvs\nWISVPZ7D7PEcBtQcODMebn+Gh9ufNv2ukN6WUfZ66yuqu8tzqKpgrXYNAMa6WtnXdpzh1nzem9HO\n/LVptbXlhEbhM2Csq36RaBEWXhz9OMeGKj8bzOqNHhsvP/9/43DxgMYnHhnM3ROt9vwz7ak9xd7q\navfNSKc6lv62FjpbHbizz8Ze78Z6bXWo3iYz6YtgEaojNMr+fvWjcW02xMcOlleA16wzb/wauHvh\n6I82dz+tXSrxxLu/A0//E+jdkgnCNJrNxiqlNDxGnwY+L6X8Eipk7/x67EBK+XmyIuz48eNrSttk\nsQie2tvNhC9CKJbiwKAXIQTJq220Oix8bHwvGSlw2CxcnV0kYnPjdth44YF+QrEkTpsVq0Vw+raP\nQ8NteFvsSCm5tRimx+2kvTVveF2bDXF1NshoZyv9bS18765hrBYr1+dCXJ4J8uKeE+zpdZNKS8KJ\nVM4Yimczf0XiKfqyv0kfTAXocjuwZAWCfVaJu9lAlH09Xexv7+bt7PyEfX1ejvR34LRbsAjB0kqC\n9lY78WSai1MB+rwtHHErwfCRkedwO6ykMjIndo4Db91YZLSzlYFsFMZ4tzIGx7qV0TThiyAltLls\ndLTmvTWv3lIi9Hv3DOAPJ0hmMvS4nQiRN8aml6PEkmlSGcloZytTy1FC8QiXAzDS4eKpPd1E4mkG\nO1q4MBlgNhDjow/10WK3cmFymQODbbTYrRwabmc5kuCNDxf45KHdvDs/w6GeQ7x+JUS/10nIGsEi\n7Lw4Noi9QGSmM5LWe/3s79rPWJubXT3q2FLpDGkpcyH386EYHcvdrCQd7Pd2kUo7Ge9yc30+xO3F\nMCPtXZzY9QB2q+Dmkg+7r50uVydv3LjDcIeLHreXlLTjslv52PggVovInYN0Rk3cf/ng32JiKcJ4\ndzs2i2AlnsrN8frIg320ufKFZid8EW4thnn2gZ7c1ABrQWhUIJLE02JDQO58J7IJTZw2CwsrcTIz\n7Vgsgld2D2G3WshkJHeWwvS3tRBLpgnGUqRcqi2f2D2Iw2Yhlkzn+kYgmiQUS9Lf1oLdaiEQSZJI\nZ+j1OrnmS3Jz+SbjbT0c7FEeGSmHmAuF6GptxWFTpm44nsLtHCaRypBIZ/A4bYTEHpKJVh4d6COT\nkQSiSXb1uLm5sEK320FHq4P5YIxer5N0RnJxKkCv18lIZyszgSinb/vo9Tp5eq8SaneXwsyHVDKU\nj+zvA+CdW0v0t7Xkrvet79gJxZM8NODFabPy8d0DZdMtwvEUU8tRdnW7czU+ezxOWh1Wri8vMRcL\nYbdY6G9z0dnSyVNDw6zEU/jDCUa7Wrm9GObC5DJOm5Xx7laWMw4eG+4kHbHz5J5uzt31s7vXTUZC\nn9dZ1E8B4uIwI95hBr2dJNMSl8NKIpXBYbPksg56nLZcPzBETTCWpKWkf/zAo+qafOLQILZ7nQgE\nzwz05Z5Znzw8xFI4kRNA8VSa5UgyZxcb85yOjXdybLyTV2+pfX1894OML0R5sHNjbSQtnDaZm4th\nRrtaGyp+a+B22hjrauXqrPY4NZ0rX4Fbr8PHfhkc6+Pmr8r3/Ayc/yP45s/DD/1x8/en0Ww/rEII\nm5QyBbwIfLbguy372zZaMgLscdrIyAwWi8CWnUQ91t3KrYgzF27kLRiNfmpvfoRWCFHuAUFlxNs/\n4GV65UnuBO5gtSij7IE+D0LA7h4PVovAaSM3agt5b02hl+vQsHmxz4F2Fw6bhb62Fn7g0WGeiX6c\naDLKaFs+vKZwQHBPtp1fu2NjzDtGu0sdU+n03GceqF7TqPT8mdHpNg9/K/XGPdDnIZqyMpPsosXW\nojyCWSfWiV3F3r/HxopH1TtaHTmD8COjHwHg+w93Ichw8h7s69xTZoxaLYKXd71c1i6b1VLUYfu8\nLVwNqJDN/vYWvA7VqMMjHRwu8VQ90NONzf4IA+4BMjJNNBXGbXcQTiZz2y5tw/4Btb39A/nr0+6y\ns6tbiTnj2hiMdrXmzrvh2SukULQbFNo0fd4WnAvqQhvnxGIRuT7hdtro9ji5ErQWrVvoPWx32Yva\nZbrPgmMVQjDQVpz42ejrDpslt4+nh4s9S0bfKbyvjEEEm1UU9YPBdleuDxiMd7tzYt/giRKvyice\nPMbdZR9pyyKA6Rx1t9PGg/3F8/eMazAft+NL5I/VqNHmcdpy9+54l6pjta/fw3hmhIi4yb6uUQZG\nlZh7usZ99tTo4fz2s80zzpndaqHdZW6zFnrOXjzQX5Ty3GGzcKjnEAuRhaLrZ7GIIq+R02alvy1/\nTo6OddLmKr/eFmHhsb76PH7ryZb9cdkp3FoIs6envsl3Zhweaeed2z4yGbmhk+N2FFE/vPovoO8g\nPPEPNmaf7h54/qeVcLr4J/DI39qY/Wo024c/Bt4QQiyiUpK/CSCEeAAVrrctMJuLYWShKq3h1ChD\nniGGPPmQGSEED/StLZnCoGeQSDJCIF58intcPVBHVFdpgdHNZj2L8CrD0rLhBYWN2k7PjjxNKpPi\n8tJlwsn6My6C6htHRlcfPriZbLfkDw/1PMBDPfDqrVdXtb5xvPu79pORGdOCwxaL4OCQGvSwW92b\nUuTa47RBSRTdeNs4422N1ScvHSwZ9Y7ij/vX2rxVo+c4bSKZjOT24kpu1GU1vHSgn4VQPJdSVrPO\nSAl/+U8hPA8/8JtgUkW8aTz1j2HkBPzVP4PgdO3lNZodhJTyl4F/jsqO96zMW8AW4CdrrS+E+GPg\nbWC/EGJSCPHjzWpro7TYWjjYc5Dj/VsjGUMhj/U9lhvl7W01T7qwnbBnM3/1t/ZvckvWjt1ix2Vz\ncaDrwGY3RbNB7Ovch9O68ZkZN5NHeh/h+ZHnN23/2uO0idxeChNLZniwf/XC6SP7e7EI+Our82Wh\nBJp14L0vwuUvw0u/CMNHN3bfVhv84G/Dbz0Lf/L34Mf+cmOFm+b/b+/eY+Sq7gOOf3/z3t1Z78sP\nsA3GWISkCY+AcZyWEmTMowHhKIaIlIZHg9I0omqTSA00VaT0n6RVkzZpqkYEKBAlgYQiYmhcRAgS\naSsw76eNMTgYuwabfb/msTO//nHPjGd3Z3ZmvbNzZ+/8PvLonnvn7vicc8+cueeec881TU5Vnyyz\nbW+Nf/vZ+seofuZ7VbaR2qPtbF23dcE9Ys0gGo6y5eQtgTr5bI82YDi5MS3KGk4+Kjz0buMpVWfO\nrai7PcbGdb08tvsIX73k9Op/YGr33quw82uw/gJvkoYZsvks6ak0qVyKyalJUlMp0rl0MZzKpaYv\np7z9MrkMefKoKnnNk9c8yvSw4E0BGw6FCZ9zOZG9jxJ64CrCGy4iEooU34uGoiTCCRKRBPFInLZw\nG4mIt54IJ2iLtNEebScZTRIPx5fckAZjgm5V+yreHX+3rs9VaoQgNJoKEpHmmOa41FkrzmLf0D6S\n0eO/sLrUbF69ua5DJ4Oq8N2L+TCVvbGGk68efukwq7sSC7rHCWDLh1by7Z17ODw82fRTh/olk8sw\nlB5iKD3EeHac8ew4E9kJbzk1weTU5PRtqQEm9z9BalUvqZ44qYe2z2oETelU9f94hpCEiIfjhCRE\niBAi4oUlhCDFpaLkNOe98jly3T1Mjb9J7uX95I/j+Svg3UDaGe0kGUuSjCbpjHWSjCZJxrxwV6yL\nnkQPvYleehI99CX66En00BXvqvjkb2PMwpy98mzSubRd1DDTdMW7OHfVudV3DJDexPFfRG4l65et\nJx6Os7pj/g8JNgtnDSefHOif4LdvvM+Xt35gwT+YW13D6bHdR/iTzc07vKOeJrIT9Kf66Z/sZyA1\nQH+qn4HJAQbTgwymBhlODzOYdsvUIBNTE1U/MyIR2qPtdITjtI8dpY0pEivPoLt9OW2RNuLh+LSe\nnMJ6W8Tr5YmH41640ANUWC/pAYqGosd3vHNZuGcbHHwOvfFX5E48k5zmyOay0xpzM3u5Cg3CsewY\no5lRxjJjjGa95Vh2jAOjBxjLjhXXywlJiO54N72J3mKjqid+rIFVbGzFvXB3vLs4i5cxZm4hCdEW\nsQteC7XpxE30T/b7HQ1TQUhC5DVffcfj0B33JrXoipefATJoRIQ1yfIPwjWLzxpOPrn36QOEBD5z\n3toFf9aGFUnW9bXz693vLbmGUy6fK57Ul75GMiOMZkYZSg8VG0X9Ka+RNJAaYHJqsuzndUQ7iifw\nvYleNnRtoDvRTXf82CsZTdIebfcaSdEOOiIdtEfbve7vPf8JO/7Ca6hc+ws4eXODc6SCcBSuvgt+\ntAX5ydVEbtxJZMUHiIfjJKnPUI5sPstweriYx4WGaP9kP4PpweL6noE9DKQGGM2UnwZfELriXcWG\nVaFRVTgOXfEuumJd3tK9lsWWFadUNcaY+Vrettyb2c80pfPXnM9AaqD6jsdhZftKtpy8pSmHXM4l\nHomTnkr7HQ0zT3am4oPJTI77nn6HLR9cWZehdSLC1g+t4sdPvs3R0XRDn6Kc1zzj2fFpDZ5Co6dc\nQ2jmtkq9HAVhCRd7Ovra+li3bF0x3JvopS/RR2+bWyZ6j3/cfXoUHv6KNxnECWfC9tthRZPdM5Zc\nCdf9Eu68FO64GK74Lnz40wt6YnypaCg6r5OPbD7LUGqIgdQAQ+khBlODDKSO9foNpAYYTA3y1tBb\nDKYHGUoPzXnFsTPaOa0xVdrAKja4SrZ3x7vpjHVa75Yxxsxw4UkXVrzA6IdkzBsWvliWWqMJ4II1\nF5DNZ/2OhpknXxpOInIZ8D0gDNyuqt/2Ix5++emuA/SPZ/izT2yo22de+7GTuet/f8e3du7mO1ef\nVXU4mKqSyWdm3etTCBfWZzZ4StdHMiOMZcbKPoekVOF+msJrdXI1nbFOlsWWTdteblsymlzc+2vy\nedjzEPzXrTB6GP7wq/CJWyDSpDdd9m2Azz8KP7/Om2nvie/AOZ+DM672nv3UQNFQlBXtK2qekjiv\necayYwynhhnODDOcHmYoPcRw2gsPZ46tj6RHODh6kOGMF56rjHXGOr2G1YxerELPYlukjfaI18NY\nadkWabP7uIwxgVEYVWGaVzQcJWoz5S450ugZTEQkDOwFLgYOAk8Dn1XV1yr9zcaNG/WZZ55Z0P+r\nqih6bIlSOBcrbinZp7B92t+Wvlf6OSXrhXDpjf1HRif47zeP0BkPMTCZ5q7/eZPTT+jgby4/ffoE\nAHOEp/JT5DVfDGfzWdK5NJlchkwuQzqX5vl33mf3uwOc0B1mTU+USCRHKpshR5apfIZM3tsvNZVi\nIjtR8+QGHdGOY42ZaOUGT7ntyWiyOXoE9v8WJvohPQKj78HwOzByCI7uheEDsOoMuOKf4KTz/I5p\nbfI5eOk+2HUb/N/zICHvmU9rz4Pe9dCxEhLLINoBU5OQmYDTLoZmOBbzVBjOWdrQGkoPMZIZmd7w\nSpc0xjLDjGfH5zWmPhaKEQ/HiUfixMNxYuEYiXCCWDhWdr3wKqwnItPfi4VjhCVMJBQhJCFvJsTC\nTIllwqFQiIhEEBEEObZ0YWD69nLbCuGZ627GtnKf7f2Tsp/dCCLyrKo23wOLZqjH75AxxpjmM5/f\nIT96nDYB+1T1LQARuRfYBlRsOC3U9h3b2TtY06M9GiK0Ft4Abnxk4Z9VONmLhWPEwjGW9wlHx5Uj\nh8KoRiAf5Yw1fZza3Tnt5K8j2lG84l4IF+71Kb33JxlLBuPekwe/5DWQCjpWwLI1sPos2PK38JHt\n3nOTlopQGM7+Y+/13mvw6gOw7zHY9SPIVRgz/bW3oW3pPRk+HAoXe5FOZvYT0isp9KpOZCeKvamT\nU5PF8LRt2QlSuRSZXKa4TOfSpKfS3jKXZjw7XrxYMW2fXHrRbnr2W6VGVkFXvIvHP/O4jzE0xhhj\nGsePHqergMtU9Sa3/jngY6p684z9vgB8AVgO9AGvNzSi/lkOvO93JBrE0hpcrZReS+vCrFPV2sZ7\n+khEjgJvL/BjWqms1MryZDbLk/IsX2azPJntePKk5t+hpr3Erqq3AbeJyDOqeorf8WkUl96mH7ZS\nD5bW4Gql9FpaW0M9GnetnH+VWJ7MZnlSnuXLbJYnsy12nvhxN/Qh4KSS9bVumzHGGGOMMcY0JT8a\nTk8Dp4nIehGJAdcAO3yIhzHGGGOMMcbUpOFD9VR1SkRuBh7Bm478TlV9dY4/ua0xMWsarZReS2tw\ntVJ6La2mVpZ/s1mezGZ5Up7ly2yWJ7Mtap40fHIIY4wxxhhjjFlq7ImPxhhjjDHGGFOFNZyMMcYY\nY4wxpoqmaTiJyGUi8rqI7BORW8q8f4OIHBWRF9zrJj/iWQ8icqeIHBGRVyq8LyLyfZcXL4nIOY2O\nY73UkNYLRWS45Lh+o9FxrBcROUlEHheR10TkVRH5yzL7BOLY1pjWIB3bhIjsEpEXXXq/WWafuIjc\n547tUyJySuNjunA1pjUw9XGjVPuNC4pKdYOI9IrIoyLyhlv2uO0V60QRud7t/4aIXO9XmupFRMIi\n8ryIPOzW17u6Yp+rO2Jue8W6RERuddtfF5FL/UlJ/YhIt4jcLyJ7RGS3iHy81cuKiHzZfXdeEZGf\nuTq5pcqKlDl3rGe5EJFzReRl9zffFxGhVqrq+wtvkog3gVOBGPAi8Hsz9rkB+IHfca1Tei8AzgFe\nqfD+J4GdgACbgaf8jvMipvVC4GG/41mntJ4InOPCncDeMuU4EMe2xrQG6dgKkHThKPAUsHnGPl8C\nfujC1wD3+R3vRUxrYOrjBuVp1d+4oLwq1Q3APwC3uO23AH/vwmXrRKAXeMste1y4x+/0LTBvvgL8\ntFAvAj8HrnHhHwJ/7sJl6xKXjy8CcWC9K1Nhv9O1wDy5G7jJhWNAdyuXFWANsB9oKykjN7RaWaHM\nuWM9ywWwy+0r7m//qNa4NUuP0yZgn6q+paoZ4F5gm89xWjSq+gQwMMcu24B71PMk0C0iJzYmdvVV\nQ1oDQ1UPq+pzLjwK7MarBEsF4tjWmNbAcMdrzK1G3WvmzDrb8E4CAO4HLprXVawmUWNazfy0zG/c\nHHVD6ffjbuBTLlypTrwUeFRVB1R1EHgUuKyBSakrEVkLXA7c7tYF2IJXV8DsPClXl2wD7lXVtKru\nB/bhla0lSUS68E6Q7wBQ1YyqDtHiZQVvxus2EYkA7cBhWqysVDh3rEu5cO8tU9Un1WtF3VPyWVU1\nS8NpDfBOyfpByp+EbXfdcPeLyEll3g+KWvMjKD4u3rCgnSLyYb8jUw+uu/yjeFfrSwXu2M6RVgjQ\nsXXDbF4AjuBVxhWPrapOAcNAX2NjWR81pBVapz6uh8B972sxo25YpaqH3VvvAqtcuFLeBC3P/hn4\nayDv1vuAIVdXwPT0VapLgpYn64GjwL+7IYy3i0gHLVxWVPUQ8I/AAbwG0zDwLFZWoH7lYo0Lz9xe\nk2ZpONXiIeAUVT0Tr9V4d5X9zdLwHLBOVc8C/gV40Of4LJiIJIH/AP5KVUf8js9iqpLWQB1bVc2p\n6tnAWmCTiHzE7zgtlhrSavWxmdNcdYO7ytsyvZgicgVwRFWf9TsuTSaCNxzr31T1o8A43hCsohYs\nKz14PSjrgdVAB0u792xR+FkumqXhdAgovWK51m0rUtV+VU271duBcxsUNz9UzY+gUNWRwrAgVf0V\nEBWR5T5H67iJSBTvZOEnqvpAmV0Cc2yrpTVox7bADSV5nNk/ZsVj64ZYdAH9jY1dfVVKa4vVx/UQ\nmO99LSrUDe8VhiW75RG3vVLeBCnP/gC4UkR+hzdMcwvwPbwhRRG3T2n6KtUlQcoT8K70Hyzp0b4f\nryHVymVlK7BfVY+qahZ4AK/8tHpZgfqVi0MuPHN7TZql4fQ0cJqbNSSGd4PbjtIdZtwHciXeuOmg\n2gFc52YK2QwMl3RPBoqInFC4D0RENuGVySV5sunScQewW1W/W2G3QBzbWtIasGO7QkS6XbgNuBjY\nM2O3HUBh1p6rgN+4q2JLSi1pbbH6uB6q/sYFxRx1Q+n343rglyXby9WJjwCXiEiPuwp/idu25Kjq\nraq6VlVPwTv2v1HVa/EuSlzldpuZJ+Xqkh3ANeLNpLYeOA3vJvclSVXfBd4RkdPdpouA12jhsoI3\nRG+ziLS771IhT1q6rDh1KRfuvRER2ezy+LqSz6pOm2D2DD02K8ZevJk/vu62/R1wpQt/C3gVb5aQ\nx4EP+h3nBaT1Z3hjV7N4V1w+D3wR+KJ7X4B/dXnxMrDR7zgvYlpvLjmuTwK/73ecF5DW8/G6jl8C\nXnCvTwbx2NaY1iAd2zOB5116XwG+4baX1lEJ4Bd4N+HuAk71O96LmNbA1McNzNdZv3FBfM1RN/QB\njwFvAL8Get3+FetE4E/d92kfcKPfaatT/lzIsVn1TnV1xT5Xd8Td9op1CfB1l1evM4+ZwJr1BZwN\nPOPKy4N4s5+1dFkBvol3seoV4Md4M+O1VFmh/Llj3coFsNHl75vADwCpNW7iPsAYY4wxxhhjTAXN\nMlTPGGOMMcYYY5qWNZyMMcYYY4wxpgprOBljjDHGGGNMFdZwMsYYY4wxxpgqrOFkjDHGGGOMMVVY\nw8kYY4wxxhhjqrCGkzHGGGOMMcZU8f/vWddeaQ3HCAAAAABJRU5ErkJggg==\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "pm.traceplot(trace[100:]);" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "trace" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "histogram = pm.Histogram(trace[100:], model=model)\n", + "histogram" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "(9900, 3)" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "histogram.histogram.get_value().shape" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true + }, + "source": [ + "# Inference results" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true + }, + "source": [ + "We don't have exact function as our observations were noisy" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": { + "collapsed": false, + "deletable": true, + "editable": true + }, + "outputs": [ + { + "data": { + "text/plain": [ + "array([ 1.04420257, 1.99506605, 2.05719256], dtype=float32)" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "histogram.histogram.get_value().mean(0)" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "isinstance(histogram, pm.variational.opvi.Approximation)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Basicly, any subclass of `pm.variational.opvi.Approximation` __must__ implement some abstract methods. One required is `.random(size)`. That allows to create unified interface for _any_ approximation we get from inference and reuse it just everywhere.\n", + "\n", + "The most frequent case is to create stochastic function and all we know about it is that it has PyMC3 variables with somehow inferred distribution. That is pretty easy with `Approximation` class. There are two great methods for that:\n", + "\n", + "* `.apply_replacements(node or list of nodes, ...)`\n", + "```\n", + " Replace variables in graph with variational approximation. By default, replaces all variables\n", + "\n", + " Parameters\n", + " ----------\n", + " node : Theano Variable(s)\n", + " node or nodes for replacements\n", + " deterministic : bool\n", + " whether to use point with highest density as initial point for distribution\n", + " if True - constant initial point will produce constant latent variables\n", + " include : list\n", + " latent variables to be replaced\n", + " exclude : list\n", + " latent variables to be excluded for replacements\n", + " more_replacements : dict\n", + " add custom replacements to graph, e.g. change input source\n", + "\n", + " Returns\n", + " -------\n", + " node(s) with replacements\n", + "```\n", + "\n", + "* `.sample_node(node or nodes, ...)`\n", + "\n", + "```\n", + " Samples given node or nodes over shared posterior\n", + "\n", + " Parameters\n", + " ----------\n", + " node : Theano Variable(s)\n", + " size : scalar\n", + " number of samples\n", + " more_replacements : dict\n", + " add custom replacements to graph, e.g. change input source\n", + "\n", + " Returns\n", + " -------\n", + " sampled node(s) with replacements\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": { + "collapsed": false, + "deletable": true, + "editable": true + }, + "outputs": [ + { + "data": { + "text/plain": [ + "array([ 18.1869278 , 16.50654984, 14.91315365, 13.40673256,\n", + " 11.98729515, 10.65483952, 9.4093647 , 8.25086784,\n", + " 7.17935371, 6.19481897, 5.29726696, 4.48669481,\n", + " 3.76310253, 3.12649107, 2.57685995, 2.11421013,\n", + " 1.73854065, 1.44985199, 1.24814403, 1.13341665,\n", + " 1.10566974, 1.16490352, 1.31111813, 1.54431319,\n", + " 1.86448908, 2.27164555, 2.76578259, 3.34690046,\n", + " 4.01499891, 4.77007771, 5.61213732, 6.54117775,\n", + " 7.557199 , 8.66020012, 9.85018349, 11.12714481,\n", + " 12.49108982, 13.94201374, 15.47992039, 17.10480499,\n", + " 18.81666946, 20.61551666, 22.50134277, 24.47415543,\n", + " 26.53393936, 28.68070984, 30.91445923, 33.23519516,\n", + " 35.64290237, 38.13759232], dtype=float32)" + ] + }, + "execution_count": 22, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "histogram.apply_replacements(y, deterministic=True).eval()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true + }, + "source": [ + "# Picture\n", + "Let's look at what we got from inference. There will be mean function." + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": { + "collapsed": false, + "deletable": true, + "editable": true + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztvXl03Od13/25AAaDGezERuzgAoqiJUtyWTWOHS+y5ai2\n3yhunLoW6zpNUjZt3GOnrhPbyknsNEqcZnHSk568ZWM76vvSdh3bOXYcL5GtXSIlkSJFiSIpUiQB\nEgRA7ACxzfb0jzuPfkMKJEESg1lwP+fMwWzAPDOY+c79fZ+7iHMOwzAMo/ApyfUCDMMwjJXBBN0w\nDKNIMEE3DMMoEkzQDcMwigQTdMMwjCLBBN0wDKNIMEE3DMMoEkzQDcMwigQTdMMwjCKhbDUfrLGx\n0fX09KzmQxqGYRQ8+/fvH3XONV3tfqsq6D09Pezbt281H9IwDKPgEZG+5dzPLBfDMIwiwQTdMAyj\nSFi2oItIqYgcEJHvpS9vEJFnROSEiPwfESnP3jINwzCMq3EtEfrHgSMZl/8I+KJzbjMwAfzKSi7M\nMAzDuDaWJegi0gG8D/jr9GUB7gK+mb7Lg8DPZ2OBhmEYxvJYboT+58BvAqn05QZg0jmXSF8+C7Qv\n9YsislNE9onIvpGRkRtarGEYhnF5riroIvJ+4Lxzbv/1PIBzbpdzbrtzbntT01XTKA3DMIzrZDl5\n6G8Bfk5E3gtUADXAXwB1IlKWjtI7gIHsLdMwDKOASSahpAREsvowV43QnXOfcc51OOd6gH8FPOyc\n2wE8AnwwfbePAt/J2ioNwzAKlVQK5udhFeY330ge+m8B/1lETqCe+pdWZkmGYRhFxPw8nDkDCwtZ\nf6hrKv13zj0KPJo+fxK4c+WXZBiGUSSkUjA0BCdPwvr1EI1m9eGsUtQwDCNbLCzA8eMwOAilpVl/\nuFVtzmUYhrFmSKVgeBhOnICpKVhczPpDWoRuGIaRDRYXVcwHBlTcs5zhAibohmEYK08qBefPq3c+\nNQXl5atiuZigG4ZhrDSxmEbn/f16uabGInTDMIyCw0fnr74K4+MQDquYp1JX/90bxATdMAxjJYnH\n1Wrp79cK0XBYrzPLxTAMo4BwDkZH1W4ZHYVIRP3zykqL0A3DMAqKzOjcOaio0FN3t/roWcYE3TAM\nYyVwDiYm4NgxzT8PhTQ6r6mB9na1X7KMCbphGMZKEI/rRmh/v3ZWjEbVctmwAdatg7Ls13GaoBuG\nYdwozmm++bFjQZl/WRlUV8PGjSrslrZoGIZRACQSuhF68qQKdyQCVVWwaRPU11vaomEYRkHgHExO\nahOuoSH1zcvK1DvfsEG99JISPWUZE3TDMIwbwUfnJ05oJF5RodH5xo3qnZeXq6ib5WIYhpHH+Oj8\nlVfUOy8r01NtrQp6KBT46avAcoZEV4jIsyLygogcFpHPp6//GxE5JSIH06fbs79cwzCM/GD3btjU\nnaC3+RSf+qUT7DnkNLOlslLFvLZWBX2VonNYXj/0ReAu59wFEQkBT4rID9K3fco5983sLc8wDCP/\n2L0bdv47h8xP8tMcIcUQf/tsBYsNId7xvnXqnYfDqxqdw/KGRDvn3IX0xVD6lP1pp4ZhGHnK/ffD\n4nyCVk7RxXFKcYwT5Xs/iGp0Xle36tE5LNNDF5FSETkInAcecs49k77pARE5JCJfFJHwZX53p4js\nE5F9IyMjK7RswzCM3NHf56hgkh6OUscwi4QppZRXWKeCHg5rZB4Kreq6liXozrmkc+52oAO4U0Ru\nAT4DbAX+KbAO+K3L/O4u59x259z2pqamFVq2YRhG7tjQqdF5D8coJckUlTgqiZf3XuydrzLXlOXi\nnJsEHgHucc4Npu2YReArwJ3ZWKBhGEZe4Rx/8JuT3JSOzueJUkYJM6zjV3+vJ8hDX0Xv3LOcLJcm\nEalLn48AdwNHRaQ1fZ0APw+8lM2FGoZh5AWJBB+68xSfes8R6klwgQg1VHPfJ3v5hY/WBE25csBy\nvkJagQdFpBT9AviGc+57IvKwiDQBAhwEfi2L6zQMw8g9qZTmnR89ytvbh3n7f4hCTSlsaIR3dgdF\nRKswzGIprirozrlDwB1LXH9XVlZkGIaRryQS2q/l8GHtrlhfr1Whvb1a6u8FPUdYpahhGMZy8NH5\nkSM6M7SyUsW7qUkHWFRU5DQ6BxN0wzCM5eH7nb/8spb8h8Mand98s0bnOdoIzcQE3TAM42qkUjA+\nDi+9pLNCw2Ftkbt+PXR2BjnnOYzOwQTdMAzj6sRi6p0fPariHono8IqbbtIoPYeZLZmYoBuGYVwJ\nH50fOqTReWWlinhHB/T0BGK+iiX+l8ME3TAM40rEYjpa7sgRtVS83dLbq+dzVBW6FCbohmEYlyOZ\nhLExePFFmJhQIY9GoatLI/RwOG+iczBBNwzDuDyxmGa1ZEbnNTWwZYsKex5F52CCbhiGsTTJpHrm\nhw7BzMzF0Xlrq0bmebARmklukyYNwzDyEedgYUErQo8d44lXyvjhcxW8yjoGqm7iY39UyYd+Kfd5\n55diEbphGMalpFJaDXrgAHv2zfC3z4UZoppTbODlC0187DdC7P7mkiMgcooJumEYRibOwfy8eucn\nTvD4gXIWqWGGdQywlQQRpmPl3P87uS0iWgoTdMMwjEwSCRgagoMHYXaWU1SQIsoJtjBLLXHCxCin\nvz/XC309JuiGYRieTO/81VchEmEdNYzQyBC9JIiQIASU0NWV68W+HhN0wzAMTyIBZ8/CgQOashiN\ncvf7owyylXkqSRIiTjnRKDzwQK4X+3pM0A3DMEA3Qmdn4fnn4fRpzTmvrOSd72/jU/9tI+1tEeKE\n6e4Wdu2CHTtyveDXc9WcGxGpAB4Hwun7f9M597sisgH4OtAA7Ac+4pyLZXOxhmEYWSMeh74+zTuP\nx7VnS3U1bNnCB2+u5IO/HoZo/hQRLcVyIvRF4C7n3G3A7cA9IvJTwB8BX3TObQYmgF/J3jINwzCy\nSCoFc3O6ETowoEVEVVVaQNTTo8Mr8qyIaCmuKuhOuZC+GEqfHHAX8M309Q+ig6INwzAKj1gMTpxQ\nuyWZVAGvq9PhFZWVar/kWRHRUizLQxeRUhE5CJwHHgJeBSadc4n0Xc4C7dlZomEYRhZJpWB6Gvbt\n02KiSET7tXR3Q3t7wUTnsExBd84lnXO3Ax3AncDW5T6AiOwUkX0ism9kZOQ6l2kYhpElFhZ0cMWL\nL+rlaBQaGnR4RWWlinmOJxEtl2vKcnHOTQKPAG8G6kTEH4N0AAOX+Z1dzrntzrntTU1NN7RYwzCM\nFSWR0MHPzzyjjbiqqtRq6e7W4c8FFJ3DMgRdRJpEpC59PgLcDRxBhf2D6bt9FPhOthZpGIax4jgH\ni4samR89qsJdWwuNjeqdV1erd15SONndy3H5W4EHRaQU/QL4hnPueyLyMvB1Efl94ADwpSyu0zAM\nY2VJJDQq37cPLlxQ37y6GjZtUmHPkzmh18JVBd05dwi4Y4nrT6J+umEYRmHhS/wPHYJXXtEMFh+d\nb92q3nlFRa5Xec0UzrGEYRjGShGPw+AgPPmkdlasqlJBv+UWPV9RURBpipdigm4YxtrCFxHt3auV\noeXlarW0tUFnZ5B3XoCYoBuGsbZYXNQiooMH1UevqtKMFp+mWGAboZkU5qoNwzCuh0RC54M++yyc\nOaMCvm6dpil2dOjlPBr6fK2YoBuGsTZwTkv8X35ZN0NLSjTn3G+EFmCa4qUU7soNwzCuhUQCJia0\niGh8PCgi6u3VytBIpKCjczBBNwxjLeCLiA4e1Ai9vBzq66GlBTZuVHEvsJzzpTBBNwyj+InHYXgY\n9uxRDz0SUe9861aN0gs0TfFSTNANwyhukkktInr2Wc1uCYU057y9XdMUvXdeBJigG4ZR3MTjOlJu\n714V9vp6WL8etmzRcv+KioLeCM2kOJ6FYRjGUiQS2uv86ad1+HNlpW6AdnaqqFdXF4V37jFBNwyj\nOPFpikeOwP79el1dHTQ3azfF2tqiEnMwQTcMo1iJxzVNce9eTVOsrlYx37RJN0QjkaLYCM3EBN0w\njOLDb4Tu3asRenm5RuetrZqmWFNTNBuhmZigG4ZRfMTjWtq/Zw9MTQVivnVrEJ0XyUZoJsX3jAzD\nWNskEtoS9+mn4dQpjc7XrYOuLu2oWFNT8BWhl2M5I+g6ReQREXlZRA6LyMfT139ORAZE5GD69N7s\nL9cwDOMK+I3Qo0c17zyV0jTF5mbYvFnPF+DgiuWynB2BBPBJ59zzIlIN7BeRh9K3fdE59yfZW55h\nGMY1EIupxfLUUzA2phuhvry/uVlL/IvQavEsZwTdIDCYPj8jIkeA9mwvzDAM45pIJtU7f+457aZY\nVha0xu3u1ui8yNIUL+WavqpEpAedL/pM+qqPicghEfmyiNSv8NoMwzCWTyymxUNPPaVRelWVlvdv\n2KADLKJREMn1KrPKsgVdRKqAbwGfcM5NA38FbAJuRyP4P73M7+0UkX0ism9kZGQFlmwYhnEJ8bgK\n+tNPw8mTmsXS3KwVoR0dWkRUZDnnS7EsQReRECrmu51z3wZwzg0755LOuRTwv4A7l/pd59wu59x2\n59z2pqamlVq3YRiG4jdCX3xR7ZbFRS3vb23V6HzduqLeCM1kOVkuAnwJOOKc+7OM61sz7vYB4KWV\nX55hGMZViMW0IvTJJ2FwUHPOm5o0q6WlRTdGi3gjNJPlHIO8BfgI8KKIHExf91ngwyJyO+CA08C/\nz8oKDcMwLkcyqYK+Zw+88EKwEdrTo3ZLU1PRb4RmspwslyeBpXYSvr/yyzEMw7gGYjHo61NBn5vT\n+aBdXZrV0tKiXvoaYm0chxiGUVDs3q1BdkmJ/ty9e4k7xWIwOwuPP66iXlGhvnlnp2a31NWtGavF\nU/zbvoZhFBS7d8POnRpwg2r1zp16fseO9J1SKc1sOXRIW+MmkyribW0q6M3NRdl862qYoBuGkVfc\nf38g5p65Ob3+NUGPxWB4mIf/5DFe+sEIZ6ljnnbu+o3NfOD/adeN0DWICbphGHlFf/9Vrk/nnH/v\n9/fw0g9eZpEwo7TSRwcPf7GNhVsa+fAvl67aevOJtWUwGYaR93R1XeF6b7UcOcL+rzyFY45p1jFC\nJ4N000crn/n82sg5XwoTdMMw8ooHHtAq/UyiUb2eWEynDz3+OOWcZYEow3QwQhfjtDNHLf1niru8\n/0qYoBuGkVfs2AG7dmnmoYj+3LULdnwoob3O9++HgweppIRx2jlPK2N0MsF6HGWXjfDXAuahG4aR\nd+zYkbEBCmq1LMS0T8vjj8PsLG96Uz17n+9imE0M00aMiiCSX6NYhG4YRv4Ti8HMDDzyCJw+DdEo\nb31vBzt+vYuqpnZmaaC7WzSS33HVv1a0WIRuGEZ+k0hbLc89B88/r5uiXV3Q2cn73t/J+77QBlUm\nZWARumEY+YzvpHjqlFot589fPLSiqwsqK3O9yrzBBN0wjPwlHtfy/scegxMndGhFR4eW+Hd1afOt\nIh9acS2YoBuGkZ8kEiroBw7Avn0q7M3NGpn7fi1rYGjFtWCCbhhG/uGtloEBjc6Hh4Nhz52dGp2v\n0fL+K2GCbhhG/hGLwcICPPEEHDmibXA7OoLT+vVmtSyBCbphGPmFz2p5/vmgz3lDQ2C1dHaa1XIZ\nljOCrlNEHhGRl0XksIh8PH39OhF5SESOp3/WZ22VqZS2xzQMo7jxVsvZs2q1DA7q0IotW1TIu7uh\npibXq8xblhOhJ4BPOue2AT8F/LqIbAM+DfzEOdcL/CR9OXukUnoyDKN48UMrHnkEjh3TnuZdXUFW\nS0uLWS1X4KqC7pwbdM49nz4/AxwB2oF7gQfTd3sQ+PlsLZKSEv0nWpRuGMVLZlbLnj0q7C0tQb55\nd7dZLVfhmjx0EekB7gCeAVqcc4Ppm4aAlhVd2aWUpvsbm6gbRvHhrZaTJ+HRR7WAqKEBenvVaunp\nsayWZbBsQReRKuBbwCecc9OZtznnHOAu83s7RWSfiOwbGRm5vlUmk3oqLTXrxTCKkVgMLlzQatBj\nx7Rfbne3ZrN0dOhoOeOqLEvQRSSEivlu59y301cPi0hr+vZW4PxSv+uc2+Wc2+6c297U1HR9qxTR\nb3DnzHoxjGIjPYGIZ5/Vfi2zs1o0lO7XwoYNZrUsk+VkuQjwJeCIc+7PMm76LvDR9PmPAt9Z+eW9\ntggV81QqmOJtom4YhU8yGVgtTzyhBUTNzbB5s4r5pk1mtVwDy4nQ3wJ8BLhLRA6mT+8FvgDcLSLH\ngXenL2eH+Xk9HINA1FMpFXnDMAoT75tPTcHDD8Px41pAtGGDZrX4n8ayuepxjHPuSeByeULvWtnl\nXIa5Of3mDoWgokIjdm+92KGYYRQmvhr0mWc0s2VxMYjM29qCkUXGsimMStHFRS00GBrS1CYfpTtn\n1othFCLJ5GvDnl+zWhob1WLp6FBhr1i7w56vl8IQdNCNkv5+mJi42E8368UwCgtvtYyOajVoX59W\nf/b2albLpk2af25cM4Uh6OPjcOYMjI3pz4WFIHXRsl4Mo7CIxdRGffRReOEFvdzTo355d7eeN66L\nwhD0+Xl9A4yNaW+Hs2dVxFOpIAPGRN0w8h+fonjgADz9tB5xr1+vG6CdnRqlh0K5XmXBUhiCHo2q\n5XL+vA6KPXtWo3ZvvYhYwZFh5Ds+RfHMGe3V0t+vE4e2btVN0Jtu0vFyxnVTGIJeVqZlwAsLunky\nNaWTvy9cCPxzs14MI3/xvvnkJDz0EBw9quPkNm1Sq6WnR4uJjBuiMAS9pETL/isqVMwHBzXj5cwZ\nzXrxFaSglw3DyC9iMbVO9+7VitD5eRVyX9a/ebOlIK8AhSHooRDU1upOuHNqt4yOwrlzGrFnZr34\n84Zh5AeJhAr6sWOa1TI6GlSDtrXBtm0arRs3TGF8JdbUQH29vil8utPwsEbskQhUVkJd3cW9Xnzx\nkWEYuSOVCqzShx+GU6f0s7pxo4r5li2af26sCIUh6M7pN/rcnAp6PK7++eCgbpj296u4V1Rc3OvF\nDuEMI3c4p0WB09Pqm7/4oh5t9/QELXG7u4PPrHHDFIbiVVXpt/r69Srmi4v6rT85qV56ebkK+4YN\nQWR+aTMvwzBWl3hcg7Bnn9UUxakp3QTt7NROijfdpJ9dY8UoDEEvL1cx9z2TfWHR+HhgvVRUqPC3\ntFw84cisF8NYfRIJ/Zx633x8XD+bGzeqoG/bZl0Us0BhCHpJiYq1F3V/qquDkREYGFAvPRpVP726\n2hp4GUauSKX0KHpwUH3zEyd0H8ynKG7dqhaqBVorTuEoXVmZbowuLMDiIt//0QLf3ZVijGp6medd\n9/bxrg+VqrDfdFPQldFXkfoRdoZhZA/vm09Owo9/DIcOXeyb+9mg9nnMCoUj6KBvjKYmvvn1GH+y\nK0YdCSLEOUuC73xnisXSId77kXTWS0+Pfgn4KlIR89MNI5t4MZ+b05a4e/eqRerngra1wS23WBfF\nLFJ4CldRwW//eTPDNDJCAwmqmaGBOUI88e1htV/6+zW10Y+tA43SrSujYWQPvwl6+LDOBh0fv3hQ\nxa23mm+eZZYzgu7LInJeRF7KuO5zIjJwyQSj7JHZfEuE44NVjNDGFK2cp5kE5UxTzwQprR4dGNDW\nADMzF1eRWmsAw8gOfhP03Dnt03LsmO5xbdigm6G33qq+uR0lZ5XlvLp/A9yzxPVfdM7dnj59f2WX\ndQm+8Va6ArSzu4R5apmkhUlamaaeRcKUUqtRwsmTKuqnTukhoHVlNIzs4TdBR0eDfPNoVKPy5uYg\nVdF886xzVUF3zj0OjK/CWi5PaWmQseIcDzwAFdEypmhkiibGaSVJDfd+MKpvpPl5bZp/+rSefL8X\nsK6MhrGSeN98YkInDz33nGagdXdrdL5pE9x8M4TDuV7pmuBGNkU/JiL/BtgHfNI5N7FCa1qasjIV\n5kSCHfeVAcL994cZ6mthY1OCnf96nvdsi8N0uqBoaCiYQVpZqZ3cMqtILT/dMG4M30FxZkaj8iee\n0POdndp0q6MDbrtNP3/GqnC9htZfAZuA24FB4E8vd0cR2Ski+0Rk38jIyHU+XBp/yJZMsmOHBt8L\nqQhPv9LEz/1Kp3p1kYi22q2o0DzYvj61YPzoOo9ZL4ZxYyQSmsVy9Cj88Ic6r8CPkGtvVzGvrzff\nfBW5rlfaOTfsnEs651LA/wLuvMJ9dznntjvntjc1NV3vOhURFfVLNkmprdVG+e3tKua+O2MqpYLe\n16d+emb/dPPTDeP6SSTU2jx3TouHTp/W4qH2dvXNt21TcTfffFW5LstFRFqdc4Ppix8AXrrS/VeU\nzBa5PrdcRCed+DadsZiKdX29jq07fVrvV1Gh3d3Ky1//NwzDWB7JpPrmIyO6CXrokH6murtV0Ht7\ntcTf8s1XnasKuoh8DXgH0CgiZ4HfBd4hIrcDDjgN/PssrvH1ZEbp3gsvK9Mo3efCzs+rUCcSuvs+\nMKBeXiSiRUdL/Q3DMK6Mb4c7OgpPPgn79+tnbONGjcg3btRN0Gg01ytdk1xV0J1zH17i6i9lYS3X\nRsYm6WsVoaGQHu75Frv+MDCZ1E3ScFh7wkQimlKVmZ9u/V4M48r4jJYLF7R46LHHtIPi5s0q5h0d\nmm9eVWUBUo4obBUrLVVBzxTkSCRo4hWPq39eXa2XBwb00DAc1sPB+nr9HeeCLwbDMF5PZkbLSy/B\n97+vdmZbmwZH7e1wxx36mTLfPGcUtoJ5qyUzUgeNENraVOhnZyGR4IdH6njq7weZoI9ZKrj3sxX8\n/H/eGpQiW/90w7g8PqPl5En4x3/UiuymJv2ctbRoj5b16/Uo2cgZhS3oEGS+JJNBV0Wf+ZJIQDzO\nPzwUZ/ffO5I00sIIpZzi//uDUlxpOR/41FaN1jOzXkzUDSMgHtfI/MwZ3QQ9cUKDptbWYCboxo1W\nPJQHFIdylZToKbMKtKREM1/a2vjiNzYxSRVxqhilgRrmqKWPv/2vL6vPHo9f3BrAmngZhpJMBjNB\nf/ITeOGFoB1uS4tmjW3dqlan+eY5p/AjdE9G0RGggl5aCo2NHJpJ0EyMHo4QIsV5kjQxSZyTcDis\nkUVXV/ClYJukhqGfg7k5TU988kkt6wfdBG1s1AKiW2/V7DE7qs0Liku1lkpFDIWo6Whm9GyScubp\n5lViQBjHZkbgZFiji3A4yHxJpWyT1Fjb+IZbY2OamvjUU/qZ6O7Wjc+uLq0EramxTdA8ovi+Vn0K\nY4Z18vkvhFkMNTHMJs7SRQlhFqniHe+u0Uq3w4e13ed4ugeZ3yC1SlJjLeJzzcfGtKz/xz8OerTU\n16tvftttQVW2kTcUn6DD69oD7NgBf/mVSqpa1zPCRlIl7fzb+yK8/Z01+gUwPKzNhY4dg+npwE+3\nzozGWsPnmo+P6+bnD3+olktLSxCZ33GHirptguYdxekpLJHOuGMH7LivCibbYMDBi4tabFRfr2/e\ngQF9g5aVaaVbVdXF/WJsw8codryYT01pssB3vwtnz6qYt7SoJfmGN6ioW1l/XlKcgg6XT2esqwsK\niVKpIJKfmNAipNJSLT7askXLlzP9dBN1o1hxTrO9pqc11/wf/kHFvK5O880bGjTXfMsWDXzss5CX\nFK+gw8X9z+HiHPXOThXreDwYgDE5qaLuI/zNm/XNa5kvRjHjA5zJSc01//GP4fhxPUrt7FQxv/lm\nuOkmDXJsEzRvKX6FyuzOCPpmLC3VyCMe1yZe8bjelkrpm/rUKf29UEgLJsrKgi8FE3Wj2PBifu6c\neuaHD6twd3Xp56S3F974Rq2qNjHPa9aGOvk3YWa73HSO+muWy0svBQOlR0dV1H3fl/b2wL7xVo5h\nFAPxuHrmg4Pw+OMq5mVlWjhUV6dHqbfequmJltGS96wNQYfXD8YoKQla7iaT+sZ+5ZXAM5+d1ctl\nZfDWt2qfCp8O6f+eYRQyiYSmI547p3nm+/frZ6SnR4uFenrUN29osIyWAmHtCDoEmS+X9lFvaQl8\n8hMngi6OU1OayhgKqag3NARtAWwwhlHIJBK6AXr2LOzZA888o0FNd7fuMW3cCLffru2oy8tzvVpj\nmaw9RfLZKn4jFFSwW1q0J0VPj/qH9fX6Rp6YgCNH9E3vC4+SyYt/3zAKCS/m587Bvn2wd6+mK/b0\nqE/e2Ql33MHXHm2lZ2sFJaVCTw/s3p3rhRtXY21F6J6lhmOUl6uo+w3Ukyf19pERFfVDh/R33/IW\nbfrlM2RCIUvhMgqHTJtl/37t0bK4qF55ZaXmmt96K197rI1f/ViYuQV9b/f1wc6d+id27Mjh+o0r\nctUIXUS+LCLnReSljOvWichDInI8/bM+u8vMApdG6pmi/oY3BIee69bp/aemtJr06adV4H0jL4vU\njUIhU8wPHNCJQwsLKuYVFbpPdNttsHEjn/18mLmFi+Vhbg7uvz9HazeWxXIsl78B7rnkuk8DP3HO\n9QI/SV8uPHxeeqao+yZdt96qaVv19SrszqmoHzqkfuPERPC7JupGvpNpsxw8qBkt8/Nqs5SXaybX\nLbdoB8WKCvrOLr3p39+/uss2ro3lzBR9XER6Lrn6XnRwNMCDwKPAb63gulaHzGrSTPvFi3pmb3Qv\n6BMT+oEAePObg0Ea/m+Z/WLkG36Df3BQ03OfeEKzuDZsUJulrU33j3p7NVIvLaWrS22WS+nqWv3l\nG8vnejdFW5xzg+nzQ0DL5e4oIjtFZJ+I7BsZGbnOh8simXnlmZF6RUUw9LarS3PWq6r09uFhPWTd\ns0c/KCLqp1t3RiPf8GI+MKCByKOPqu3S1QXV1fzo+Hre/Z9uJvrmm9l0a5Td39Bc8wce0NyATKJR\nvd7IX254U9Q550Tksn6Dc24XsAtg+/bt+elL+PTFpSL1zk69j4/UUymdrTg6qtNbkkn46Z8OKk/B\nqkmN/CAe1yPKoSENQJ54Qj3ztJX4o8Pr+PT/vpmjbCVOlJNnQq/b+Lz/frVZurpUzG1DNL+5XuUZ\nFpFW59ygiLQC51dyUTlhKfulpCSYZgQXTzSan4fz5/UQNpnUPHUTdSMf8L1ZxscvFvNYLCjnr6vj\nT7+1lWNxh5lqAAAVrElEQVTcTJxKkmhk7jc+d+wITkbhcL2Wy3eBj6bPfxT4zsosJ8dk+uDefikp\n0U2j7m61X7q71XOsrNTfGR7WculHH9WBABA0/EJzd3t69M9YLq+RdbyYj46qzbJ3byDm3d1qHdbV\nwc0389TsGy4Sc49tfBYuVw0jReRr6AZoo4icBX4X+ALwDRH5FaAP+JfZXOSqcrlIPRTS6MZXmIpo\nxsD8vOaq+w/SO9+pqY7xOF/9Kuz8D2XMzemftlxeI6v4fuYTE/re3LNHo/N4XKOJ2lrtydLbC9u2\n0dYZ5cSZ11eB2sZn4bKcLJcPX+amd63wWvKHJQZkXCTqHuf0kHZ6OhD1VAre/nZobubzvx1ncQ4y\nX+bMQ1rDWDGcU398ZETfk/v2BWK+YYNu6NfVaT/zbdugqorP/WGYnTt5LeAA2/gsdMzovRKZol5a\nGjT08pG6F/qysqCi9Phx/d2f+Rn6zrRQRhxwFx3W2iGtsaL4GaDnzwc2y8sv6/WbNqk92NioBUTb\ntunlcNg2PosQE/SrkdnQCwIB7+wMBN234x0e1gyYtKi/qfIO9s9uogz1072o2yGtsWKkUhpinzun\n77+9e7X3kIiKeTSqFmBvrw6pSIu5xzY+iwsT9OWQKerOqXiXlWnxUWmpWjFlZfrhGhnR+77yCp/7\nuXk+/rVF+thCCE17DEdDdkhrrAzJpAYQQ0Mq6I89prNAQyEV8FBIK51vukkLh6JRa4Nb5JigLxef\np545+ch3afRTkEIhvc/0NExP856u0/z3X4zxhb+PsXdhG5s74Xf+K9y3wwYFGDdIPB70Zenv1x5D\nZ86oV97VpcLd2KhtcLdu1cjc2uAWPSbo10Lm5CN/uaxMPzi33x50Xjx9WnOAp6e5uzvJ3Q8k4JZZ\neNObNEqKOftwGdeHH+Y8NqaReV+fivnwsForLS2BZ97bG3joNm1oTWCCfq34PHVvv5SV6am2Vpsb\nlZXph+fUKd2k8v7m4qJuXG3frvcFa71rXBuZmSyDgzp85cABLe1vatJTRYWK+saNugkajZqYryFM\n0K8HP6nI56p7y6WmRlvv+sg9FNLIaXY2sGsWFzVSb27Wy+Gwibpxdfzm59CQvqcOHtTNz7k5bXu7\nbp2+57q7NTLv6lJxNzFfU6y9iUUrhc928QVFoJcrKzVSv+UW9S7b2vRDdeGCHiYfPqytS0+ehMVF\nvvaVBXq6nVWSrjGuqYLYt77t79eRcc88oy0nFhe1gdy6dRqJb96sAUVXF0QiJuZrEIvQb4RLm3r5\nyBy0gCMcDgqSzpzRTSzQQ+VEgm/vnuA3/mQbkw4cYfr6SqySdA2wezcXFfRctoLYOS3Z9w22zp7V\nfvynT+t7zQ9zjkY1jXbLFo3W0y1wjbWHuFUczLB9+3a3b9++VXu8VSWz94tvHRCP87d/OcxXP3eY\n+dnT/FNe5WffPMlb/0kUqqr42Bda2MNWTvJG5qkjRhhHKd3d+pk1ipOenqV7jV/0f3dO20qMjqpf\nfuoUPP+8RurRqAp3JKJZLe3tmprY1KRBhIl50SEi+51z2692P4vQV4rMtMb0Zunurwr/8XfakPly\neojwLJWM7TlKqnSct908RowUG4lTzSxHuJUZWolRQX+/HSoXM5erFH7t+mRSj+aGhnQD9Phxrfyc\nnta88sbGYJB5R4dmstTXa2ReYi7qWsYEfSXJzIBJJLj/t0uZmS+llEZO8EZiVJAkzN89eZy33TZM\nD+OUEqOCBaJc4Bg3c55NrO+oAqwApFi54jQgn5I4OKibn8ePw4kT6pe3tEBDg0bhTU0amff2qu1S\nUWGb64Ztiq44frMUONevm6UJyligltNs4wQ3cwBtw/vWt0VoZIoqxmjkBLdygDexlwc+kU53TOe7\nWwve4mLJaUARxx/+zrymuPb16ab53r1w9Kj+4zs7g8i8rU2j8i1boLparRcTcwOL0LNDerO0sytB\nf3+SFCUkKSNGJefopbm+Em4J87byciR0lsd/MsZ5FmkhyT2/OM+7Ombh1Zuhs5Pdf1/Dzl8rsRa8\nRcSlTbF6OhL84Wem+dCbh6F/VKPyY8c0M8oXCUUiKt7NzZpj3tpqOebG67BN0Syyezf82r9LsjCf\nwiEkKaUqkuJ//mWc+352TKOwkyf1Azw2phurNTWahtbaClu2cNsvbOal8w2kLhlCYBunRYAvFBob\n08h8bAxefVW/tWdn9X3gxbymRm2WjRv1ZzRqm59rCNsUzQM0Eivl/s+WcK4/QU9Xgs//fin37SiH\nRLqqz6ednTql4dr4eFBVurBA1flxWtjKGG3EqAT00Npa8BY4ySRf3zXN//jtQabHx7iJM/za+17l\nrt5xtVja2zUir6nRyuKWltdmgRKJ2OansSQ3JOgichqYAZJAYjnfIGsNbU8q4NJZMC4JLp2bXlOj\nXmhVlX5ow2HNNZ6Y0HS1xUXeSJwok5xjA+fYzBRNOMqsBW+hko7Kv/FXY/zRJwcpY5QeTrOOPn74\nD/NwVwV33V2v74n6eh1K0d6u2Sy2+WlchZWI0N/pnBtdgb9T3GQWIWWkNlJZqR/WcDhIRTt9WrMc\nxsb48FsX+PKT64gwSx3jnGYLs6WdPPBAVa6fkXGt+IrP4WG+8rkRKjlLC2doZ4Q4wgANfP3heu76\nxZAKeWOjemtNTdZgy1gWZrmsNn7yka8uLSnRqMsXivjD7GgU+vt529Z5pHSU7zy2yKvE2FAyw7/4\n9AT/4t2bINVkh96FgO/DMjKiDdvGxymd6aOXM1SywAwhFqlhmlqGqIWW9EZoV5emKUYi5pcby+JG\nBd0B/ygiDvifzrldl95BRHYCOwG6zCdQLo3WRYJhBN5Xr6tTUR8e5mfCY/xM7yhEZqFxChqmYO+I\nbpB1d+uXgB2G5yd+aPPgoHZFHBiAwUFuY4xRUsxSyTSNTBMhQS310SrY0KZf8LW1ZrEY18SNCvpb\nnXMDItIMPCQiR51zj2feIS3yu0CzXG7w8YqLzEKkZFIvV1ZqhktlpfqoJ05oH5hz5zSNbWBAMyD8\ngIPhYR0C3NamkZyRH3h7ZWhIs1dGRvR/NTAAi4vc/bMhvvKjGsapZ5YaEkRIUsd993dAd7NNFzKu\nixsSdOfcQPrneRH5O+BO4PEr/5ZxEX4mqRd1EY3KfPvdmho97K6q0igvfcjO/LxGfNPTuoHqS8Cb\nm+3wPJekUvrFOzqqIu57sQwP6xdwSQlUVfGO91Qz11TFg1+vZyxRSaS6gY99rpUP/mqdfpnb/9C4\nDq5b0EWkEihxzs2kz78H+L0VW9law0frqZSe/NCMcFgj74aGIFofGtIofWBABX1qCiYnVTw2bAhs\nGBOF1cP75OPjKuATE/o/GRrS/0s8rmJeXa3/15oa3ntvBe/9aJ0eXTU3620WlRs3wI1E6C3A34n6\ne2XAV51zP1yRVa1VMqN1761HIiru0agKwbp1Gq2fP6+iMT8fpDrOzam4Dw2pSHR2mrBnG18cNDoa\n/E/GxvTnyIhG6yLBjE8v6BUVuk/S1qZ7J5WVQetlw7hOrvsd5Jw7Cdy2gmsxPJmZMM6p9VJdrXNI\nIxFNYztxQn310VEVjwsX9Lrx8eCQf3BQRb293Q7jVxrf3nZ8XIV8akptlakp/X/MzOjrHQ6/FpFT\nUxNkMjU0aLGQF3fDWAEsJMhXfCZMKhV4634KTXm5ikNbm1aY+jarfpP0wgUViqkpFZyzZ9WGWb9e\nI32LBK8f59TuGh3V13ZiQm2v8+f158yM3kdEj6R8GmpVlf7/fNVnba1eZ7nlxgpin+x851IbprQ0\n8Fq9t376tIr2yIiK+IULGp3PzAQNncbH9ZC/pUWzaKqq9IvBUuKWRyKhEfnIiJ6mp/U0Oamv88SE\nCrmfE1tfr699dbVG4NXVemTV0BBcZxgrjAl6oZBpw/hipLIy/VlVpVkuR44EmRXxeOCnj48HDb/G\nxlT8m5r0cl2dTbm5HM7p6zg9Hfji/qjnwoUgx9xvZJeU6Oua6ZP76t+mJr0uGrViMCNrmKAXEt6G\ncU6FvaxMhTgU0lNVlfrq/f0q7OGwbpTOzqqwT06qqKxfr0I0PKyC3tqqouPtmLUuOJkZK2NjQSQ+\nP6/n5+Y0Ks8cO+g3rL1PHonovkVTk762Zq8Yq4AJeiGS6a+nUmqdZAp7e7umNPb1BTZMLKZCNDam\nJx81NjbqZd93u6kp8Ht9KuVawDmNuCcn9cvOpx3Oz+vrtrCg52dmNGpPJvV191W9ftCEP/mI3L+W\na+V1NHKKCXoh4/11n+IYjarILC7qaKO2NhX1c+dU2OfmAmH32RkNDSo87e163cCARpmZsyvLy4sz\ncvci7gu0xsf1aObCBX2N5uf1tLCg90sm9fXzBV+1tfo6lZfr0VBlZdC33LdFLrbXzMhrTNCLgUxh\n9yPwystViDZtUktlcFCFfXw8yGmfnQ0i0qGhQKBaWtSHP3NGI8y6uoujzdLSwsyU8ZuWCwsq4Bcu\nBCI+M6OiHYvp7fG4irq/DvQ19V+APnPFv9b+tfM+ue1JGDmgAD+VxmXxwu5cILrhsEaUfhblyIhG\n4WNjeltNjUahXtjPn1fhr6sL/N/R0eBv1NYGghYOB3aPf+x8wu81xGKBaM/MqJgvLur5WExPqVQg\n5n7ASDKpfyccDr7UIpFAyEPpNrd1dYF3bj65kUNM0IsRkUBkvfDE4yq+vlXv2Jhuig4Pq3jV1Kio\nTU+rsA8MqLh7AWtpURthfFz/Riikt0WjeruP3r094x9fZHX840zxXlwMvG+/kbmwoJfn5oL7JhJ6\n33g8EPF4XP+ef34+6g6Hg+ftC70aG/U1MSE38gQT9GLGC7vfMC0vD6yEUEizMjo6VLhHR9WCCIVU\nwLwoTk7q9cPD+vtexHy0WlYW2A7eS/bDOqJRFcGSEr72rXI+93slnBoso6ujhM//fqmO6Fuu4Puq\nWVAh9qd4POg+6dc8OxuIczyuIp8p+D4aTySCSDweD45CfH6/98H9F1RpqV7vWxv7528YeYK9G9cK\nJSUquF7YvZh5Ae/o0M3BsTGNwhcWVPQaG/W+U1P6RTA4qOLo27v6yNxHquFwUABVWgrl5Xz3J6X8\nxR+HKCHEBkpJni3nD36pjIqBEL/wwZJgfZf6zr4Dpa+WzRRxb5UkEoFgJxLBNKhUKrivv21+Xv/O\n4qJ+icRi+rjeOmls1C8gv08QCunt0WiQyeK/pMwjN/IQE/S1hh+mEQoFkfTCQiDKDQ0qfL7J1Oys\nimAkEuRne3vC57j77pCRSCCI3t4Jhdj9x+WEKKFZF0AKoRTHV38vxC+8Ib0uH337aF0kEGYfXftN\n31RK7+MnPnlh9oVAJSWBeGdG6d4T9wVZNTUq1FVVQUZKKBS8Rj7bx6ciWmWtkeeYoK9lSkqCjU0v\nft62qKrSiHVxUS0X39mxsjKIfDM9aG/l+C8Db0WUlhKmgkbCOEpIECJJOQ5hZlFgKt0uNlMovZh7\nvBCLBMU8/rS4GETgS0X0vrmZb2xWVxd4/5FIMNsVguKszJJ9X5FrGAWAvVONoFDJZ8VUVgYZH/5n\nY2OQl+2zQzLF3FsgvujGX59IsJ55ypnF4RAc/m1XSzn0pT10nyXjPfXMCN3bG4mE/rw0Uvf+ui+y\n8n3H/dGCL8GvqAgsIR/Rl5eriNfU6PP2ewH+7xhGAWGCblyMF1ZvofhINx4PqiX9yVsvCwuBR11S\ncnHKXyLBOz88z1e+lmARBzjKiRPBcc8HYkFvGi/WmWvwEXbm9V7c0/78a31ovI3ihTsUCr6kfOGV\n74Loh3D7DB4fwYdC5o0bBY0JurE0khE5+8jdt3/1dktm3vb09Osj9XQ2yT//UAOx6iT/e1eKEZJ0\nSIIP/7Lj7ncQbGT6iNuLrrdCfITuM3Z8Gqbf5M28LXPtEGxslpUFIu4tpkzBNxE3ioQbEnQRuQf4\nC6AU+Gvn3BdWZFVGfpEp7qBi6C0O71s3NweZM96L9358PM699yW590PJwCbx4u03NDNP/j7eu/aX\nvSWTKfp+LV7cfdqkP3nx9huevoWB2SlGEXIjM0VLgf8B3A2cBZ4Tke86515eqcUZeYoXQx8BZ6YK\nVlVdnJni/XR/3ts3mbd5wc7MavGZM3BxBO4j80yB9idvyfijiswKVhNwYw1wIxH6ncCJ9Cg6ROTr\nwL2ACfpaI1Pg4eIsFC/OmaLtT95n95chsEu87+3JPJ+5gbrUTxNvY41yI4LeDpzJuHwW+GeX3klE\ndgI7Abq6um7g4YyC4UrVn164/U8v9JfefqW/5wXbcsIN4yKyvinqnNsF7ALYvn37Ep9WY02RWTgE\nFk0bxgpyI5+mAaAz43JH+jrDMAwjB9yIoD8H9IrIBhEpB/4V8N2VWZZhGIZxrVy35eKcS4jIx4Af\noWmLX3bOHV6xlRmGYRjXxA156M657wPfX6G1GIZhGDeA7UgZhmEUCSbohmEYRYIJumEYRpFggm4Y\nhlEkiFuqMi9bDyYyAvSt2gOuHI3AaK4XsYqstecL9pzXCoX6nLudc01Xu9OqCnqhIiL7nHPbc72O\n1WKtPV+w57xWKPbnbJaLYRhGkWCCbhiGUSSYoC+PXblewCqz1p4v2HNeKxT1czYP3TAMo0iwCN0w\nDKNIMEG/BkTkkyLiRKQx12vJNiLyxyJyVEQOicjfiUhdrteULUTkHhE5JiInROTTuV5PthGRThF5\nREReFpHDIvLxXK9pNRCRUhE5ICLfy/VasoUJ+jIRkU7gPUB/rteySjwE3OKceyPwCvCZHK8nK2TM\nxv3nwDbgwyKyLberyjoJ4JPOuW3ATwG/vgaeM8DHgSO5XkQ2MUFfPl8EfhNYE5sOzrl/dM4l0hf3\nogNMipHXZuM652KAn41btDjnBp1zz6fPz6Ai157bVWUXEekA3gf8da7Xkk1M0JeBiNwLDDjnXsj1\nWnLELwM/yPUissRSs3GLWtwyEZEe4A7gmdyuJOv8ORqQpa52x0Im6zNFCwUR+TGwfomb7gc+i9ot\nRcWVnrNz7jvp+9yPHqLvXs21GdlHRKqAbwGfcM5N53o92UJE3g+cd87tF5F35Ho92cQEPY1z7t1L\nXS8itwIbgBdEBxt3AM+LyJ3OuaFVXOKKc7nn7BGRXwLeD7zLFW9+65qcjSsiIVTMdzvnvp3r9WSZ\ntwA/JyLvBSqAGhH5/51z/zrH61pxLA/9GhGR08B251whNvhZNiJyD/BnwNudcyO5Xk+2EJEydNP3\nXaiQPwfcV8zjFEUjkweBcefcJ3K9ntUkHaH/F+fc+3O9lmxgHrpxOf4SqAYeEpGDIvL/5npB2SC9\n8etn4x4BvlHMYp7mLcBHgLvS/9uD6ejVKHAsQjcMwygSLEI3DMMoEkzQDcMwigQTdMMwjCLBBN0w\nDKNIMEE3DMMoEkzQDcMwigQTdMMwjCLBBN0wDKNI+L9p80VAaYWDVwAAAABJRU5ErkJggg==\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "x.set_value(np.linspace(-5, 5, dtype='float32'))\n", + "plt.plot(x_obs, y_obs, 'bo');\n", + "y_dist = histogram.sample_node(y).eval()\n", + "for line in y_dist:\n", + " plt.plot(np.linspace(-5, 5), line, 'r-', alpha=.02);" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true + }, + "source": [ + "We need an optima for real function but we know only we distribution over possible functions. To find the minima we can get the minimum of expected function. It can be done out of the box with PyMC3 tools." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true + }, + "source": [ + "# Optimization\n", + "At first we must create someting we wrt want to minimize and what to minimize. In our problem it is `x` coordinate. " + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": { + "collapsed": false, + "deletable": true, + "editable": true + }, + "outputs": [], + "source": [ + "# wrt what to minimize\n", + "# set x far from optima\n", + "x_ = theano.shared(pm.floatX(10.), 'x')\n", + "# and what to minimize\n", + "# Notice that we use random variable in out loss so it bocomes stochastic\n", + "y_ = f(x_ ,abc[0], abc[1], abc[2])" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": { + "collapsed": false, + "deletable": true, + "editable": true + }, + "outputs": [], + "source": [ + "N = 2000\n", + "hist_x = np.empty(N)\n", + "def cb(approx, current_loss, i):\n", + " hist_x[i] = x_.get_value()" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": { + "collapsed": false, + "deletable": true, + "editable": true + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "E_q[Loss] = 1.1633: 100%|██████████| 2000/2000 [00:00<00:00, 20250.94it/s]\n" + ] + } + ], + "source": [ + "sgd = functools.partial(pm.sgd, learning_rate=.001)\n", + "opt = pm.Optimizer(histogram, y_, [x_], optimizer=sgd)\n", + "opt.fit(N, callbacks=[cb])" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true + }, + "source": [ + "# Optimization results\n", + "It's all right!" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": { + "collapsed": false, + "deletable": true, + "editable": true + }, + "outputs": [ + { + "data": { + "text/plain": [ + "(array(-0.7868053913116455, dtype=float32), array([-1.]))" + ] + }, + "execution_count": 27, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "x_.get_value(), min_" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true + }, + "source": [ + "# Loss history" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": { + "collapsed": false, + "deletable": true, + "editable": true, + "scrolled": false + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6gAAAD8CAYAAAB6tolUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4XNWB9/HvmaZR75JlybLkgivuVNsYTAATamghIZRs\nAoFsks1mU0iFJdn3ZUOSzb6b3SWwAZINHUJJQjMQYwwGbAM27lW2ZVuyeh9pynn/mNF45CrbkmYk\n/T7Po2fuPffcO2d0nxnNT+fcc421FhEREREREZF4c8S7ASIiIiIiIiKggCoiIiIiIiIJQgFVRERE\nREREEoICqoiIiIiIiCQEBVQRERERERFJCAqoIiIiIiIikhAUUEVERERERCQhKKCKiIiIiIhIQlBA\nFRERERERkYTgincDAPLy8mxZWVm8myEiIiIiIiL9YNWqVbXW2vxj1TupgGqMeQi4FNhvrZ0aKcsB\nngTKgArgOmttw9GOU1ZWxsqVK0+mKSIiIiIiIpKgjDE7e1PvZIf4PgIsOqjsTuANa+144I3IuoiI\niIiIiMhRnVRAtdYuBeoPKr4C+H1k+ffAlSfzHInA5w/GuwkiIiIiIiJDXn9MklRord0XWa4CCvvh\nOQZMU7ufufe+yY+e/4SddW3xbo6IiIiIiMiQ1a+TJFlrrTHGHm6bMeY24DaA0tLS/mzGSekMBvnU\npEKeXLGbx97fxcWnFnHHgrFMLc6Md9NERERERCTC7/dTWVmJz+eLd1OGNa/XS0lJCW63+4T2N9Ye\nNj/2/gDGlAF/iZkkaRNwrrV2nzGmCFhirZ1wtGPMmTPHJvokSdXNPh5atoNH399Fa2eAeePyuH3B\nWOaOy8UYE+/miYiIiIgMazt27CA9PZ3cXH0/jxdrLXV1dbS0tFBeXt5jmzFmlbV2zrGO0R9DfF8E\nbo4s3wy80A/PMeAKM7x8/9OTePf7C/neoolsqm7hC797n8t/8w5/WbOXYOjkgr6IiIiIiJw4n8+n\ncBpnxhhyc3NPqhf7pAKqMeZxYDkwwRhTaYz5EnAvcIExZgvwqcj6kJHhdXPHuWN5+7vn8X+vOpXW\nzgBfe+wjzvvFEv53eQVtnYF4N1FEREREZFhSOI2/kz0HJ3UNqrX2c0fYdP7JHHcw8LqdfO70Uq6b\nM4rX1lVx/1vb+PEL67jv1U1cf3opN501mpLslHg3U0REREREZNDojyG+w4rTYbj41CKe//u5PHP7\nWcwfn8/vlu3gnJ//jTv+uIoVFfWc7HW+IiIiIiKS2Hbv3k15eTn19eG7cDY0NFBeXk5FRUVc2lNW\nVkZtbS0AZ5999gkf55FHHmHv3r191axjUkDtI8YY5pTl8J83zGLpd8/j1nPG8O62Oq69fzmX/WYZ\nz66qpDOg+6mKiIiIiAxFo0aN4o477uDOO+8E4M477+S2226jrKys3587EDj6ZYbvvvvuCR9bAXUI\nKM5K5vsXT2L59xfysyun0tEV5J+eXs3ce//Gvy3ezP4WTX0tIiIiIjLU/OM//iPvvfcev/71r1m2\nbBnf/va3D1vvD3/4A9OmTWP69OnceOONAFRUVLBw4UKmTZvG+eefz65du45afsstt3D77bdzxhln\n8N3vfpe6ujouvPBCpkyZwpe//OUeozjT0tIAWLJkCeeeey7XXHMNEydO5IYbbojWu+eeezjttNOY\nOnUqt912G9ZannnmGVauXMkNN9zAjBkz6OjoYNWqVSxYsIDZs2dz0UUXsW/fvj79HZ70bWb6wmC4\nzczJCIUsb2+t5ZF3dvC3TTW4nYZLTi3i5rPLmDEqSxdzi4iIiIicpA0bNjBp0iQA/vnP61i/t7lP\njz95ZAZ3XTblmPVeffVVFi1axGuvvcYFF1xwyPZ169bxmc98hnfffZe8vDzq6+vJycnhsssu45pr\nruHmm2/moYce4sUXX+T5558/Yvktt9xCbW0tL7zwAk6nk2984xvk5eXxk5/8hL/+9a9ceuml1NTU\nkJeXR1paGq2trSxZsoQrrriCdevWMXLkSObOnct9993HvHnzou0AuPHGG7nuuuu47LLLOPfcc/nF\nL37BnDlz8Pv9LFiwgBdeeIH8/HyefPJJXn31VR566KEerzH2XHTr7W1mTmqSJOkdh8Ow4JR8FpyS\nz47aNn7/bgXPrKrk+Y/3MrU4gxvPHM3l04tJ9jjj3VQRERERETkJL7/8MkVFRaxdu/awAfXNN9/k\n2muvJS8vDyAaCpcvX86f/vQnIBwQv/vd7x61HODaa6/F6QxniKVLl0brXXLJJWRnZx+2faeffjol\nJSUAzJgxg4qKCubNm8ff/vY3fv7zn9Pe3k59fT1Tpkzhsssu67Hvpk2beryuYDBIUVHRCfyWjkwB\ndYCV56Vy9+VT+PZFE3juoz38cflOvvfsJ/zsrxu4ZnYJXzhzNGPz0+LdTBERERGRQas3PZ394eOP\nP2bx4sW89957zJs3j+uvv77PA1ys1NTU494nKSkpuux0OgkEAvh8Pr761a+ycuVKRo0axd13333Y\ne5laa5kyZQrLly8/qXYfja5BjZO0JBc3njmaV745n6e+chbnTSjgj+/t5PxfvsXnH3yPlz/Zhz8Y\ninczRURERESkF6y13HHHHfz617+mtLSU73znO4e9BnXhwoU8/fTT1NXVAURn/T377LN54oknAHj0\n0UeZP3/+UcsPds455/DYY48B4V7choaGXre9O4zm5eXR2trKM888E92Wnp5OS0sLABMmTKCmpiYa\nUP1+P+vWrev18/SGelDjzBjD6eU5nF6eQ03LZJ5auZvH3t/FHY9+SEF6Ep87vZTPnV7KiExvvJsq\nIiIiIiJH8OCDD1JaWhod/vrVr36Vhx9+mLfeeosFCxZE602ZMoUf/vCHLFiwAKfTycyZM3nkkUf4\nj//4D774xS9y3333kZ+fz8MPPwxwxPKD3XXXXXzuc59jypQpnH322ZSWlva67VlZWdx6661MnTqV\nESNGcNppp0W3dU/GlJyczPLly3nmmWf4xje+QVNTE4FAgG9+85tMmdJ3PdaaJCkBBUOWv23czx/f\n38lbm2twGMN5Ewr47GmjOHdCPm6nOr5FRERERGIdbmIeiQ9NkjTEOB2GT00u5FOTC9lV186jH+zk\n2VV7eH1DNXlpSVw9q5hr54xiXIGuVRURERERkaFDATXBleam8P2LJ/HtCyewZFMNT63czf8s28Fv\nl25n9uhsrptTwiXTRpKWpFMpIiIiIiKDm1LNIOF2OrhgciEXTC5kf4uP5z/aw5MrdvO9Zz/hn/+8\nnktOLeK600YxZ3S27qsqIiIiIsOStVbfhePsZC8hVUAdhArSvdx2zlhunT+Gj3Y38tSK3fx59V6e\nXlVJeV4q184p4epZJRRmaGIlERERERkevF4vdXV15ObmKqTGibWWuro6vN4TzyGaJGmIaO8K8NIn\nVTy1cjcf7KjHYeC8CQVcO2cUCycW4HFpYiURERERGbr8fj+VlZWHvX+nDByv10tJSQlut7tHeW8n\nSVJAHYJ21Lbx9MrdPPthJdXNnWSluLnk1CKunFnM7NJsHA79R0lERERERAaOAqoQCIZ4e0stz320\nh9fWV+HzhyjOSubKmSO5ckYx4wvT491EEREREREZBhRQpYe2zgCvra/i+Y/28vaWGkIWJhdlcOXM\nkVw+vZgRmbpeVURERERE+ocCqhxRTUsnf1mzl+c/3svq3Y0YA2eNyeXKGcUsOnUEGV73sQ8iIiIi\nIiLSSwqo0is7att4/qM9vPDxHirq2vG4HJw/sYArZozk3AkFeN3OeDdRREREREQGOQVUOS7WWlZX\nNvH8R3v4y5q91LZ2kepxcsHkQi6ZNpJzTskjyaWwKiIiIiIix08BVU5YIBhi+fY6/rpmH6+sq6Kx\n3U96kosLJhdy6fQi5o3L121rRERERESk1xRQpU/4gyHe3VbHX1bv5dV1VTT7AmR4XVwweQQXTx3B\nvPF5GgYsIiIiIiJHpYAqfa4rEOKdrbX8Zc0+XltfRYsvQKrHybkTC1g0ZQTnTSwgLckV72aKiIiI\niEiC6W1AVZqQXvO4HJw3sYDzJhbQFTiV97bX8cq6Kl5bV8Vf1+zD43Qwf3weF00dwacmFZKT6ol3\nk0VEREREZBBRD6qctGDI8uGuBl5ZW8Ura6vY09iBw8AZ5bksmjqCC6cUUpSZHO9mioiIiIhInGiI\nr8SFtZZ1e5vDYXVdFVv3twIwZWQG508sYOGkQqYVZ+JwmDi3VEREREREBooCqiSErftbWby+mjc3\nVrNqZwMhC3lpSSycmM/CiYXMG5+n61ZFRERERIY4BVRJOA1tXby1uYY3Nu5nyab9tPgCeJwOzhiT\nw/kTCzh/UiGjclLi3UwREREREeljCqiS0PzBEKt2NvDGhmre2Lif7TVtAJxSmMbCiYWcP6mAmaOy\ncDl1v1URERERkcFOAVUGlR21bby5cT9vbqzm/e31BEKWrBQ3500oYOHEAs45JZ/MZHe8mykiIiIi\nIidAAVUGrWafn7c31/LGxmqWbKqhvq0Lh4Hpo7KYPy6PeePzmVmahVu9qyIiIiIig4ICqgwJwZDl\n492NLNm0n2Vba1m9u5GQhVSPkzPH5DJ3XB7zx+cxriANYzQzsIiIiIhIIlJAlSGpqcPP8m11LNta\nw7IttVTUtQMwIsMbDatzx+WRn54U55aKiIiIiEi3uAdUY0wF0AIEgcDRGqOAKidqd30772yt5e2t\ntby7tZaGdj8AE0ekM29cHvPG53FGeS7JHmecWyoiIiIiMnwlSkCdY62tPVZdBVTpC6GQZd3eZt7e\nWsM7W2tZUdFAVyCEx+lg9uhs5o3PY964PKYWZ+J0aDiwiIiIiMhAUUCVYa+jK8iKinqWba3l7S21\nbNjXDEBWipuzx+Yyb1w+88fn6d6rIiIiIiL9LBEC6g6gAbDAb621Dxy0/TbgNoDS0tLZO3fu7Jd2\niHSrbe3kna21LNtSy7Kttexr8gFQmpMSGQqcw+nlORRlJse5pSIiIiIiQ0siBNRia+0eY0wBsBj4\nurV26eHqqgdVBpq1lm01beHrV7fU8t72Olo7AwCUZCdzelkOp5XncFpZNmPzNUOwiIiIiMjJiHtA\nPagxdwOt1tpfHG67AqrEWyAYYsO+Fj6oqGdlRT0rKuqpbe0CICfVw5zR2ZxensOcshymjMzQPVhF\nRERERI5DbwOqq5+ePBVwWGtbIssXAvf0x3OJ9AWX08GpJZmcWpLJl+aVY61lR20bKyrqWVHRwIqK\nel5bXw1AisfJzNIsTivL4fSyHGaUZpHi6Ze3koiIiIjIsNJf36oLgeciwyJdwGPW2lf66blE+pwx\nhjH5aYzJT+Ozp5UCUN3sY0VFPSsrGvhgRz3//sYWrAWXwzClOJPTy7I5rSzcy5qT6onzKxARERER\nGXwGZIjvsWiIrwxGzT4/q3Y2sGJHOLR+XNlIVyAEwLiCtHAPa3k2c0bnUJKdrOtYRURERGTYSqhr\nUI9FAVWGAp8/yCd7mvhgR/ga1lUVDbREJl7KT09iekkmpxZnMW1UJtNLstTLKiIiIiLDRlyvQRUZ\njrxuJ6eV5XBaWQ4AwZBlU1ULKyrqWb27kTV7mnhj4366/ydUkp3MtJJMppVkMa0kk6nFmWR43XF8\nBSIiIiIi8aWAKtJPnA7D5JEZTB6ZES1r8flZu6eZNZXhwLqmspGXPqmKbh+Tn8r0SGCdVpLJlJGZ\neN3OeDRfRERERGTAKaCKDKB0r5uzxuZy1tjcaFl9WxdrKhv5pLKJ1ZVNvLO1luc+2gOEQ+4phelM\nK86MDg2eMCJdt7kRERERkSFJ16CKJKDqZl94WHBlE6srG/lkTxON7X4APC4Hk4oymB4ZHjy9JJMx\n+Wk4HZqESUREREQSkyZJEhlCrLXsru+IhtXVuxtZu6eJtq4gAKkeJ1OKM8MTMUVCa2lOimYOFhER\nEZGEoEmSRIYQYwyluSmU5qZw2fSRQHgSpu01raypDF/Lurqyid8v30lXYAcAWSluTi3OjE7ENGVk\nBsVZut2NiIiIiCQu9aCKDCFdgRCbq1uioXVNZRObqlsIhsLv8/QkFxOL0plUlMHEERlMLEpnQmE6\nqUn6X5WIiIiI9B8N8RURIHx/1nV7m9lY1czGfS1s2NfMxqoWWiP3aDUGRuekRAPrpKIMJo3IoCQ7\nGYeuaxURERGRPqAhviIChO/POnt0NrNHZ0fLrLVUNnREw2p3eH11fVX0Pq2pHicTRkR6W4symDQi\nnfEF6WSm6F6tIiIiItI/1IMqIlHtXQE2V7eycV8zG/Y1s6GqhY37mmn2BaJ1CtKTGF+YxviC9Ojj\nKYVpZKV44thyEREREUlk6kEVkeOW4nExY1QWM0ZlRcustexr8rGxqpkt1a1s2d/KluoWnlq5m/bI\nLMIAeWlJjC9I45TCNMYVpkeW08lJVXAVERERkd5RQBWRozLGMDIrmZFZySycWBgtD4Us+5p9bK5u\nYWt1K1v2t7C5upVnP9wTvb4VIDfVw7hIWB1fmMa4/DTG5KdRmJGkGYVFREREpAcFVBE5IQ6HoTgr\nmeKsZM6bUBAtt9ZS1exjc3W4p3Xr/lY2V7fw/Md7aIkZKpzicVKel0p5Xipj8tMYk5fKmPzwerpX\n17mKiIiIDEcKqCLSp4wxFGUmU5SZzIJT8qPl1lr2t3SydX8r22vb2F7Tyo7aNtZUNvHSJ/sIxVwO\nn5eWxJj81JjQmsaY/FRGZafgcTni8KpEREREZCAooIrIgDDGUJjhpTDDy9xxeT22dQaC7KprZ3tt\nGztiwuvi9dXUtXVF6zkdhlHZyYzJT4v0vKZSnptKaW4KRZnJOHVbHBEREZFBTQFVROIuyeVkfGE6\n4wvTD9nW1O5ne204sIbDaxvba9t4d1stPn8oWs/tNJRkp1Cak8Lo3O7HVEpzwsvJHudAviQRERER\nOQEKqCKS0DJT3MwszWZmaXaP8u5JmnbWtrGzvp2dde3srm9nZ30bH+5q6HG9K4RvjxMOrqmMzg2H\n2FE5KYzOSSEn1aMJm0REREQSgAKqiAxKsZM0nX3QNmstje1+dta3s6u+nV11beysa2dnfTvvbK3l\n2Q99PeqnJbkoyU6mJDsl8nhgeVR2ChnJLgVYERERkQGggCoiQ44xhuxUD9mpnh73dO3m8wepbAj3\nuu6sC4fYyoYOKhvaeW97XY/b5ACkJ7koVoAVERER6XcKqCIy7HjdTsYVpDOu4NBrXq21NHX4o4E1\n/Bhe3l3fzvJttbR1BXvsk+JxUpTpZWRWMiMyvBRlJTMyM/xYlOmlKNOrW+eIiIiI9IICqohIDGMM\nWSkeslI8TC3OPGR79/Dh2AC7t6mDfY0+9jV1sKmqhZrWTqztuV96kouiLC8jMiPhNTMSXrPCyyOz\nvKR49JEsIiIiw5u+DYmIHIfY4cOnlhwaYAG6AiGqm31UNfvY29jBviYfVU0HltfvbaK2teuQ/TK8\nLkZGel2jQTbSGzsiM3yLntQkfWyLiIjI0KVvOiIifczjcjAqJzxL8JF0BoJUN3Wyt6kjHF6jvbDh\nntjVlU3Utx0aYlM9TgozvOSnJ1GQ4aUwPYmCjCQK0r09HtOTdF2siIiIDD4KqCIicZDkclKam0Jp\n7pFDrM8f7BFe97d0sr8l8tjsY01lI/ubO+nwBw/Z1+t2UJjhpSD9oPCanhQuz0iiID2JzGS3gqyI\niIgkDAVUEZEE5XU7KctLpSwv9Yh1rLW0dgaobg6H15qWTqqbfexv7owG2g37mnlrc+chsxNDuLc3\nHGIPBNm8tCRy0zzkpYWX8yPrGl4sIiIi/U3fNkREBjFjDOleN+leN+MK0o5at60zEO19DYfX2GUf\nW2taWb69jqYO/2H3T3Y7yUv3kJsaCa7p4RCbm+ohLz0pGmjz0jzqmRUREZETooAqIjJMpCa5KE9y\nUX6UHlkIT/JU39ZFbWsnNa2d1LZ0UtfWRW1LJ7WtndS2dlHZ0M7Huxupb+skZA89hssRnkwqN9VD\ndoqHnDQPOSkeclI95KaFy3JTD5Rnp3pwOx399MpFRERksFBAFRGRHjwuByMiMwcfSzBkaWgPh9m6\n1kiobemkvq2L+rYu6tq6aGjrYsPeZurauo7YOwvhWYxzUj2RnySyU9xkp3rISnGTneIhK9lNVoqH\n7NTwemayG6/b2ZcvXUREROJMAVVERE6Y02GiQ3t7IxAM0dDujwbY+rYu6tu7qG/tor6tk/p2P/Vt\nnVQ2tLN2j5+G9i46A6EjHi/Z7SQ7xR25d20kyKa4Y5Y9ke1uMpPdZCS7yfAq2IqIiCQqBVQRERkw\nLqeD/PQk8tN7F2gBOrqCNLR30djup7G9i4Z2P40d4fWGtvB6U0f4cUNVM03tfho7/AQPN/Y4Isnl\niAbWzGQ3GV4Xmcnug8pilpMPbE/TLXxERET6jQKqiIgktGSPk2RPMiOzknu9TyhkaekM0BgJtg3t\n4eHFzb4AzR1+mjv8NEV+mn1+alo72VbTFl23R862OAzRntjMgwJsz1AbDr7hSaxcpCW5SPO6SPO4\ncDgUcEVERA6nXwKqMWYR8O+AE/gfa+29/fE8IiIih+NwmGh4HJ17fPuGQpbWrgBN7eGw2hQJtM0d\ngWiAjYbbyGNVsy9a1nWUIcnd0pLCgTXdGwmtkeX0JHfPda+LtEhZeLsrsuwmxe1U0BURkSGnzwOq\nMcYJ/CdwAVAJrDDGvGitXd/XzyUiItLXHA4T7gn1uk9of58/GA60Pj9NHQFaOwO0+Py0+rqXwz+t\nnf4e6/uafNF6bV3BYz6PMZDmcUXDa4rHRWqSk1SPi9QkFykeZ/QxLenA9p71IuseFylJTs2kLCIi\ncdcfPainA1uttdsBjDFPAFcACqgiIjLked1OvG4nBRnHngX5SIIhS2tnONC2+sIBtyW6HAm3vgDN\nvgMBuL0rSFtngLrWLtq6ArR3BmnrCuDzH7tHt5vH5SDVExNikyLhNRJ2u4Pt4cJuSkzoTe0OxB4n\nLoVeERE5Dv0RUIuB3THrlcAZ/fA8IiIiQ5IzZojyyQoEQ7T7g9HA2t4ZpLUzQHtXuKe2PRKE27uC\nPYJtW3dZZ4Da1s5IWXj9aDMrH8zjdOB1O8LXEkfCe/dystuJN2Y52RPZ7naSHNnHG7PtcPsne5wk\nuRyauEpEZIiI2yRJxpjbgNsASktL49UMERGRIc3ldJDhdJzwkOXDiQ290bDbGYyG3rbOQOQniC8Q\npKMriM8fpMMfXu7wh9ebOvzh5UhZuLz34TdWzxB7UCA+TMD1upwkuR14XQ6S3OGQ6408JkW3hR8P\n2eZy6PpfEZF+0h8BdQ8wKma9JFLWg7X2AeABgDlz5hxlvkQRERFJJP0ReruFQpbOQCgaWI8Ubjti\nQ23MckdXqEf9Fl+AmpbOHvt3dAUJHOU2RL3hcYaDa3e4DQfZcDjuDrJe94FA6z1svQPbPC4HHqcD\nj8uBO/KY5Oq57jno0amQLCJDUH8E1BXAeGNMOeFgej3w+X54HhERERliHA4TubWQs1+fJxAM0RUM\n0ekP4QsE6fSH6AyEw21nIERnINyb23msbbHlMfUb27sO2ie87PMHOclsHOV0GDxOB26nwRMJu+FA\naw4JvN3bwvUjyzFlnpiyHgH5oPruyPO5HOFHt9OBq/vRYXAdtF1Dr0XkePV5QLXWBowxXwNeJXyb\nmYestev6+nlERERETpTL6cDldJDiGfjnDgQPDrzhYNsVCIV/gqEey/6Y9c5ACH/QRraF9/EHw73O\n0foHHaPFH6AupqzH8SLL/cXpMOEg6wgHWZfTgTsSZF0x5bHBNzbwdgdgl8OBx3Xo9tjjhfc3uF2O\n43o+d6Re7PO5HQ7cLgVtkXjol2tQrbUvAS/1x7FFREREBrPucJyaFLepQHqw1hII2Who9UcC9OEC\nbVcwRCBoCYRCdAUtgci6PxR5DIYIhMLl/ph1f8x+/sh+/tCB/WOP2+EPP/Y8Xnedns93skO1e8vp\nMIcEZpfD4HQYHA5wmvCy02FwGIPLaXAagyOynyNmu9Nx0LbIeuyyM7J/7PHCxwCnwxHZFll2ED1+\nj2M4jv2csdsdkTYc+pw92+RwgCvShu7l2N+BwrycrMT4ZBQRERGRuDDGRIfrpibFuzXHpztcHxyS\nDw3EMaE22DNo+w+q173sjwTi7u3hYx/YHgqFnztkLcFQzI89dFv3PwCCB9eNqR+0lmAw8hiCYChE\nMGQJWQiEQoRCRPdPZA5Dj8DrcMSE56MF9cNsi4Znh8FpDgTycLkDpyG6v/PgYxz0TwOHw+Aw4TDv\niOwXXTbmMOUx22L2NTGvz2HC75/usB5djt3mOLAc+7sJ/67C+8XWNxxoj4m2J/LcHDhOtI7jwLYU\nj3NI/INAAVVEREREBqUD4RqS6d/rlhOFteHQGowE4EAk6PYmMMfuE4rZHozdFuw+RiQY9wjM9Azg\nB+1/SOg+Sptin7PntgPhPBiydAaCBC0H2nus19Qj6B/YnuC5vk+sv+ciUjyDP94N/lcgIiIiIjJM\nhHvp0CzOJ+DgcG8t4YBsLTZ0YLnHttCB5e59u49z8PagjdkW6l4/sG93yLYQfY7YY0XLsIRC9GhH\ndx3b/TpCMWWRR7fTEe9fcZ9QQBURERERkSFP4X5wGBoxW0RERERERAY9BVQRERERERFJCMba+F8x\nbIypAXbGux3HkAfUxrsR0oPOSWLSeUk8OieJSecl8eicJCadl8Sjc5KYEv28jLbW5h+rUkIE1MHA\nGLPSWjsn3u2QA3ROEpPOS+LROUlMOi+JR+ckMem8JB6dk8Q0VM6LhviKiIiIiIhIQlBAFRERERER\nkYSggNp7D8S7AXIInZPEpPOSeHROEpPOS+LROUlMOi+JR+ckMQ2J86JrUEVERERERCQhqAdVRERE\nREREEoICqoiIiIiIiCQEBdRjMMYsMsZsMsZsNcbcGe/2DBfGmFHGmL8ZY9YbY9YZY/4hUn63MWaP\nMebjyM+nY/b5fuQ8bTLGXBS/1g9txpgKY8wnkd//ykhZjjFmsTFmS+QxO1JujDH/L3Je1hhjZsW3\n9UOPMWZjDPvcAAAgAElEQVRCzPvhY2NMszHmm3qvDDxjzEPGmP3GmLUxZcf93jDG3Bypv8UYc3M8\nXstQcoTzcp8xZmPkd/+cMSYrUl5mjOmIed/cH7PP7Mhn39bIuTPxeD1DwRHOyXF/Zuk7Wt86wnl5\nMuacVBhjPo6U670yAI7yfXho/22x1urnCD+AE9gGjAE8wGpgcrzbNRx+gCJgVmQ5HdgMTAbuBr59\nmPqTI+cnCSiPnDdnvF/HUPwBKoC8g8p+DtwZWb4T+NfI8qeBlwEDnAm8H+/2D+WfyGdWFTBa75W4\n/P7PAWYBa2PKjuu9AeQA2yOP2ZHl7Hi/tsH8c4TzciHgiiz/a8x5KYutd9BxPoicKxM5dxfH+7UN\n1p8jnJPj+szSd7SBOS8Hbf8l8JPIst4rA3NOjvR9eEj/bVEP6tGdDmy11m631nYBTwBXxLlNw4K1\ndp+19sPIcguwASg+yi5XAE9YazuttTuArYTPnwyMK4DfR5Z/D1wZU/4HG/YekGWMKYpHA4eJ84Ft\n1tqdR6mj90o/sdYuBeoPKj7e98ZFwGJrbb21tgFYDCzq/9YPXYc7L9ba16y1gcjqe0DJ0Y4ROTcZ\n1tr3bPjb3h84cC7lOB3hvXIkR/rM0ne0Pna08xLpBb0OePxox9B7pW8d5fvwkP7booB6dMXA7pj1\nSo4ekqQfGGPKgJnA+5Gir0WGLTzUPaQBnauBZIHXjDGrjDG3RcoKrbX7IstVQGFkWedlYF1Pzy8P\neq/E3/G+N3R+Bt7fEe5x6FZujPnIGPOWMWZ+pKyY8LnopvPSP47nM0vvlYE1H6i21m6JKdN7ZQAd\n9H14SP9tUUCVhGaMSQOeBb5prW0G/hsYC8wA9hEebiIDa561dhZwMfD3xphzYjdG/mOq+1cNMGOM\nB7gceDpSpPdKgtF7I/EYY34IBIBHI0X7gFJr7UzgW8BjxpiMeLVvmNFnVmL7HD3/Aar3ygA6zPfh\nqKH4t0UB9ej2AKNi1ksiZTIAjDFuwm/GR621fwKw1lZba4PW2hDwIAeGJupcDRBr7Z7I437gOcLn\noLp76G7kcX+kus7LwLkY+NBaWw16rySQ431v6PwMEGPMLcClwA2RL3hEhpHWRZZXEb7G8RTC5yB2\nGLDOSx87gc8svVcGiDHGBVwFPNldpvfKwDnc92GG+N8WBdSjWwGMN8aUR3onrgdejHObhoXItQ6/\nAzZYa38VUx57/eJngO6Z5l4ErjfGJBljyoHxhC/Slz5kjEk1xqR3LxOeaGQt4d9/94xwNwMvRJZf\nBG6KzCp3JtAUMyRF+laP/27rvZIwjve98SpwoTEmOzLE8cJImfQhY8wi4LvA5dba9pjyfGOMM7I8\nhvD7Y3vk3DQbY86M/H26iQPnUvrACXxm6TvawPkUsNFaGx26q/fKwDjS92GG+N8WV7wbkMistQFj\nzNcIn0An8JC1dl2cmzVczAVuBD4xkSnNgR8AnzPGzCA8lKEC+AqAtXadMeYpYD3h4Vp/b60NDnir\nh75C4Lnw5yUu4DFr7SvGmBXAU8aYLwE7CU+kAPAS4RnltgLtwBcHvslDX+SfBRcQeT9E/FzvlYFl\njHkcOBfIM8ZUAncB93Ic7w1rbb0x5qeEv3wD3GOt7e1kMnIYRzgv3yc8K+ziyOfZe9ba2wnPYnqP\nMcYPhIDbY37/XwUeAZIJX7Mae92qHIcjnJNzj/czS9/R+tbhzou19nccOr8B6L0yUI70fXhI/20x\nkVEtIiIiIiIiInGlIb4iIiIiIiKSEBRQRUREREREJCEooIqIiIiIiEhCSIhJkvLy8mxZWVm8myEi\nIiIiIiL9YNWqVbXW2vxj1UuIgFpWVsbKlSvj3QwRERERERHpB8aYnb2ppyG+IiIiIiIikhAUUI/B\nHwyxoqKeJZv2x7spIiIiIiIiQ5oC6jFUNfm49v7l3PLwCt7aXBPv5oiIiIiIiAxZCXENaiIblZMS\nXX75k32keJxMK8kkyeWMY6tERERERCSW3++nsrISn88X76YMa16vl5KSEtxu9wntb6y1fdyk4zdn\nzhybyJMkPfHBLu780yfR9bH5qbz+rQUYY+LYKhERERER6bZjxw7S09PJzc3V9/Q4sdZSV1dHS0sL\n5eXlPbYZY1ZZa+cc6xga4tsL159eyuXTR0bXt9W08eyHe6hr7cQfDMWxZSIiIiIiAuDz+RRO48wY\nQ25u7kn1Yiug9tKiqSN6rH/76dXM/tnr/J+XNsSpRSIiIiIiEkvhNP5O9hwooPbSp08touLeS3j1\nm+f0KH99Q3WcWiQiIiIiIonCWsu8efN4+eWXo2VPP/00ixYtGvC2nHvuuXRfQvnpT3+axsbGEzrO\n888/z/r16/uyacekgHqcRmZ5e6zvru9g6/7WOLVGREREREQSgTGG+++/n29961v4fD5aW1v5wQ9+\nwH/+53/26/MGAoGjbn/ppZfIyso6oWMroA4C6V43L/z9XL4078BFv/+2eHMcWyQiIiIiIolg6tSp\nXHbZZfzrv/4r99xzDzfddBNjx47tUeeVV15h1qxZTJ8+nfPPPx+A+vp6rrzySqZNm8aZZ57JmjVr\njlp+9913c+ONNzJ37lxuvPFGOjo6uP7665k0aRKf+cxn6OjoiD5fWVkZtbW1VFRUMGnSJG699Vam\nTJnChRdeGK334IMPctpppzF9+nSuvvpq2tvbeffdd3nxxRf5zne+w4wZM9i2bRvbtm1j0aJFzJ49\nm/nz57Nx48Y+/x3qNjMnYPqorOjPNx7/iOpmHx/uamB6SRZOh8a9i4iIiIjE0z//eR3r9zb36TEn\nj8zgrsumHLPeXXfdxaxZs/B4PBx8p5KamhpuvfVWli5dSnl5OfX19dF9Zs6cyfPPP8+bb77JTTfd\nxMcff3zEcoD169ezbNkykpOT+dWvfkVKSgobNmxgzZo1zJo167Bt27JlC48//jgPPvgg1113Hc8+\n+yxf+MIXuOqqq7j11lsB+NGPfsTvfvc7vv71r3P55Zdz6aWXcs011wBw/vnnc//99zN+/Hjef/99\nvvrVr/Lmm2+e8O/0cBRQT8Ll00eyZON+/vTRHq76r3e5dX45P7xkcrybJSIiIiIicZKamspnP/tZ\n0tLSSEpK6rHtvffe45xzzonegiUnJweAZcuW8eyzzwKwcOFC6urqaG5uPmI5wOWXX05ycjIAS5cu\n5Rvf+AYA06ZNY9q0aYdtW3l5OTNmzABg9uzZVFRUALB27Vp+9KMf0djYSGtrKxdddNEh+7a2tvLu\nu+9y7bXXRss6OzuP/xd0DAqoJ2lE5oFrUp/7aC9fWTCWzkAIj9NBdoqbZl+AnFRPHFsoIiIiIjK8\n9Kansz85HA4cjv69mjI1NfW494kNzE6nMzrE95ZbbuH5559n+vTpPPLIIyxZsuSQfUOhEFlZWdEe\n3P6ia1BPUlFWcnS5trWTOT97nbn3vslp//I6P35hHbN+upiugO6VKiIiIiIy3J155pksXbqUHTt2\nAESH+M6fP59HH30UgCVLlpCXl0dGRsYRyw92zjnn8NhjjwHh3tDua1V7q6WlhaKiIvx+f/T5ANLT\n02lpaQEgIyOD8vJynn76aSA8a/Hq1auP63l6QwH1JJ1anHnEbY9/sAuAhvaugWqOiIiIiIgkqPz8\nfB544AGuuuoqpk+fzmc/+1kgPOnRqlWrmDZtGnfeeSe///3vj1p+sDvuuIPW1lYmTZrET37yE2bP\nnn1c7frpT3/KGWecwdy5c5k4cWK0/Prrr+e+++5j5syZbNu2jUcffZTf/e53TJ8+nSlTpvDCCy+c\n4G/iyIy1ts8PerzmzJljD76AeLDwB0OM/+HLR63z8j/MZ3RuCr96bTNfWziOrBQN+RURERER6Usb\nNmxg0qRJ8W6GcPhzYYxZZa2dc6x91YN6ktxOBxt/uog//N3pR6yzprKR19ZV8z/LdvDTv2wYwNaJ\niIiIiIgMHgqofcDrdnLOKflU3HsJFfdewhfnlvXY/r1nP2F1ZSMAm6r7drprERERERGRoUIBtR9c\nPasEgBEZB2b4ffidCgAqGzpIhGHVIiIiIiIiieaYAdUY85AxZr8xZm1MWY4xZrExZkvkMTtSbowx\n/88Ys9UYs8YYc/g7xA5xeWnh6ZuvmV3CY7eeES2/6azRNLb7+eN7O7n7xXV85r/eiVcTRURERESG\nHHUExd/JnoPe3Af1EeA3wB9iyu4E3rDW3muMuTOy/j3gYmB85OcM4L8jj8PKiEwvy753HkWZyTgd\nhitmjGRcfhq3nzuWdXub+fEL66J1G9q6yI7cJ3XVzgamlWTidqpjW0RERETkeHi9Xurq6sjNzcUY\nE+/mDEvWWurq6vB6vceufAS9msXXGFMG/MVaOzWyvgk411q7zxhTBCyx1k4wxvw2svz4wfWOdvzB\nPIvv8dpZ18aC+5ZE179yzhi+/+lJ7GnsYO69b3L1rBJ+ed30+DVQRERERGQQ8vv9VFZW4vP54t2U\nYc3r9VJSUoLb7e5R3ttZfHvTg3o4hTGhswoojCwXA7tj6lVGyo4aUIeT0bmpZHhdNPsCAPx26XZe\nW1/Njto2AJ79sJLzJxVQ09LJzWeXxbGlIiIiIiKDh9vtpry8PN7NkJN0ogE1ylprjTHHPdDYGHMb\ncBtAaWnpyTZjUHnyK2cRspat+1v5hyc+jobTbl999EMABVQRERERERlWTjSgVhtjimKG+O6PlO8B\nRsXUK4mUHcJa+wDwAISH+J5gOwalSUUZAEwuyiAz2c0tD6+Ic4tERERERETi70Rn43kRuDmyfDPw\nQkz5TZHZfM8Emo51/elwZozh3AkFPH37WYfdXt3sIxgaVtldRERERESGsWNOkmSMeRw4F8gDqoG7\ngOeBp4BSYCdwnbW23oSny/oNsAhoB75orT3m7EfDaZKkI6lr7eStzTU899Ee8tOT+NOH4Y7nry8c\nxz9dOCHOrRMRERERETlxvZ0kqVez+PY3BdSe/rx6L19//CMAHAa2/MunefyDXVw9q4RkjzPOrRMR\nERERETk+/T2Lr/SjMfmp0eWQhbE/eAmAD3c20Ozz8+NLJzM6N/VIu4uIiIiIiAxKJ3oNqvSjKSMz\n2fjTRVw1s7hH+Z8+2sPrG/bzL3/dEKeWiYiIiIiI9B8F1ATldTv55XXTufuyyWSnuEmNGdr72vpq\nfv9uRfwaJyIiIiIi0g8UUBOYMYZb5pbz0U8u5OEvnh4t/9SkAu56cR2Pf7ALnz8YxxaKiIiIiIj0\nHQXUQSIvzRNd/u8vzCY/PYnv/+kTJv74Ff7pqdX88LlP2FbTGscWioiIiIiInBwF1EEiLz0puux2\nOrjhjNLo+rMfVvLo+7u4+N/fZmNVczyaJyIiIiIictI0i+8gkZ7kIifVwzcWjgPg6wvHk5+exH/9\nbRt7GjsA6AqEuPq/3iUz2c1vbpjFrNLseDZZRERERETkuOg+qENAQ1sXM3+6+JDy2xeM5c6LJ8ah\nRSIiIiIiIgf09j6oGuI7BGSneg5bfv9b21i1s2GAWyMiIiIiInJiFFCHiP/90un8/OppzBiVxfcW\nHeg1/fXrm1m7pwlrLS99so/d9e2sqKiPY0tFREREREQOT9egDhHzx+cDcN1po7DWkpns5gfPfcLb\nW2p5e8uyQ+qvvutCMpPdA91MERERERGRI1IP6hBkjOHzZ5RSnpd6xDpPr9w9gC0SERERERE5NgXU\nIexbF5wCwPcvnkh+zG1qAH721w1c/ptlBEPxnyRLREREREQENIvvsLJubxPffOJjRuem8PqG/dHy\neePy+I/PzWR3QzvTSrLi2EIRERERERmKNIuvHGLKyEwWf2sBk4syepQv21rLzJ8u5vLfvMO2mlZW\n7aynqcMfp1aKiIiIiMhwpUmShqGvnjeO4uxkvvfsJ4dsO/+Xb0WXN9yziGSPcyCbJiIiIiIiw5h6\nUIchr9vJZ08rjV6X+tk5owBwmJ71ZtzzGu9uqx3o5omIiIiIyDClHtRh7I1/WoA/ECLd62b26Gwu\nnzGS6367nK5AiI1VLXQGQnz+wfcB+M5FE/jKOWNwOfU/DRERERER6R+aJEl68AdDAPz0L+v5w/Kd\nPbbdeOZoRuem8OX5Y+LRNBERERERGaR6O0mSelClB3ekh/SeK6Zy0ZQR3PA/70e3/e974cA6tTiT\n08tyCFobrS8iIiIiInKy1IMqR9XU7ueJFbuYWJTBr1/fzEe7GgEoyvTiMIYfXzqZs8flEgpZMrxu\nHAdfyCoiIiIiIsNeb3tQFVDluPyflzbwwNLtR9z+q+umc9WskgFskYiIiIiIJDoN8ZV+8b1FE5lW\nkklRZjLZKW4WxtyWBuBbT60mJ9XDy59U8c9XTMHr1m1qRERERESkdxRQ5bg4HYZLp42Mrp9RnsP7\nO+p71Lnl4RUAjC9MY2ZpFrNH5wxoG0VEREREZHDSEF85KYFgiCdW7KaxvYvN1a28uHrvIXV+duVU\n0r0uijKTmVaSqV5VEREREZFhRkN8ZUC4nA6+cOZoAKy1bKtpxet2smpnQ7TOj55fG13+wacncsvZ\n5Tz7YSVXzyqhwx8kM9k94O0WEREREZHEc1I9qMaYCqAFCAIBa+0cY0wO8CRQBlQA11lrG450DFAP\n6lDk8wdxOgyvrK3i649/dNg6pTkp7Kpv53c3z+H8SYUD3EIRERERERkoA9mDep61tjZm/U7gDWvt\nvcaYOyPr3+uD55FBpHsY72XTR7Jo6gje2FBNTUsnj76/i41VLQDsqm8H4NV1VYzKScHtdFCel4q1\nFmN0uxoRERERkeGmL3pQ58QGVGPMJuBca+0+Y0wRsMRaO+Fox1EP6vBireXVddU8uWIXf9tU02Pb\nlJEZjMpO4f4bZ8epdSIiIiIi0tcGqgfVAq8ZYyzwW2vtA0ChtXZfZHsVoLGb0oMxhkVTR7Bo6ghq\nWzv51lOrWbo5HFTX7W1m3d5m1u5pYu2eJp5cuZvirGR+8/lZcW61iIiIiIj0t5PtQS221u4xxhQA\ni4GvAy9aa7Ni6jRYa7MPs+9twG0ApaWls3fu3HnC7ZDB75evbeI/3tzK/V+Yze1/XHXI9iXfPpcO\nf5CqJh9et5OzxubGoZUiIiIiInIietuD2me3mTHG3A20AreiIb5ynLoCIVbtbOCssbl868mP+dNH\neyjMSGLGqCxeXVd9SP1371zI4vXVXDunhBSPJqMWEREREUlk/R5QjTGpgMNa2xJZXgzcA5wP1MVM\nkpRjrf3u0Y6lgCqxfP4gv3h1EzefXcaonBSe/2gP33zy48PWLUhP4lOTC7lmdgm+riBnjc3VBEsi\nIiIiIglmIALqGOC5yKoLeMxa+y/GmFzgKaAU2En4NjP1RzuWAqoci88fJMnloKHdz2/f2sbGqhbe\n2lxz2LpP334Wp5XlaDZgEREREZEEMeBDfE+GAqqciKZ2P4s3VNPc4efnr27E5w9Ft33+jFLe2lTD\nl+eXs6/Jx5fnlVOQ4Y1ja0VEREREhi8FVBlW/MEQ//jkx7y+obpHUI318j/MZ1JRBp2BIEku5wC3\nUERERERk+FJAlWFrf4uPV9ZW8ZMX1h2xzhfnlnHB5ELcTgftXUEWnJI/gC0UERERERleFFBl2Htn\nay276tu5fPpI/v2NLTywdPsR637tvHHcfHYZr6zdx0VTRmg4sIiIiIhIH1JAFTlIMGTpDAR54oPd\n/OmjStbuaT5i3YunjmD++HyKs5PVuyoiIiIicpIUUEWOYn+Lj1+8uokRGV5uPKuMz/52Odtr2yjK\n9LKvydejborHydljc/mnCycwqSgjTi0WERERERm8FFBFjoPPHyQYsqQmuXhw6XZC1vLg29upbe3q\nUS831UNOqoepxZl8aV452akeXllbxadPHUFRZnKcWi8iIiIiktgUUEX6wFMrdtMZCBKy8NA7O9hZ\n137EujecUcrnTi8lGLKML0wjxeMawJaKiIiIiCQuBVSRfrB4fTUG2NvUwZJNNby5cf9h643NT+W6\nOaMAmDsuj1E5KaR4nLidjgFsrYiIiIhIYlBAFRkANS2dpHicOB2Gv3tkBe9uqztq/f971anMLM0i\n3eumucNPTqqHQs0YLCIiIiJDnAKqSJz4/EF+tXgzLT4/72+vZ3tt21HrLzglny/NK2dqcSY5qZ4B\naqWIiIiIyMBRQBVJEHsbO1i7p4nqZh8+f4gWn5+H36mgpTNwSN3po7LAWlZXNgGwaMoIfn7tNDK8\n7oFutoiIiIhIn1FAFUlgNS2d+PxBCjO8NHZ0sW5PM29srOaDHfU0dwSoau55q5vsFDeZyW5mj85h\ndWUjAA/eNIfyvNR4NF9ERERE5LgooIoMUv5giN317aR73azf18z72+toaPfz6roq6tu6Dqk/Jj+V\na2aXkJvqYXxhOh/sqOec8flMHhm+Z6u1FmPMQL8MEREREZEoBVSRIcZaS01LJ3/9ZB/NHQH+7fXN\nR61fmpNCTUsnHf4gX5pXzs1nlWEMpHic5KYlDVCrRUREREQUUEWGvFU766ls6KAww8sHO+q5ZFoR\n22vaeOz9neysa6ckJ4Wlm2sOu29WipsRGV6SXA6unTOKFl+AcyfkMyY/lSSXc4BfiYiIiIgMdQqo\nIkJlQztdgRD+oOWP7+2kqcOPzx9kyeYaRmZ6cTsdbNnfGq2fl+bhtLIcOvxBijKTyUl1YzDUt3fx\n+dNLCVnLqcWZhCw4HRo2LCIiIiK9o4AqIkfU2N5FWpILf9Dyv+9VUJqTws66dl5dV0VThx+300FF\nXRs+f+iQfYsyvTS2+5kxKovpo7KYVJROQ1sXIzKTmTU6i7zUJBwKryIiIiISQwFVRE6KtZbF66up\nbvbRGQjRGQhP3rSiop7a1i6aOvyH3S8vLYkkl4P89CTSklwkuRyMyknBWstFU0cwOjeV4qzkAX41\nIiIiIhJPvQ2oroFojIgMPsYYLpwy4ojb9zV1sKehg0DI8uGuBnxdQdxOB+v2NtPS6ae5I8AHO+rp\nCh7ohf398p0ATByRTnOHn71NPq4/bRRj89NITXJxxpgcclI8JLkd/GX1PmaNzqY8L5XuDlnNRiwi\nIiIytCmgisgJKcpMpigz3BN65pjcw9bZ19RBXWsXBRlJVDZ08N72Opo6/Hy8q5Ha1k4A/rx6L21d\nwV49Z1luCjNLs5lWkklWipu/rqni9Q3VfGpSAYumFnFGeQ4l2ckYY3rcXicUshp2LCIiIjIIaIiv\niMSFtZZtNW2MzU+lrq2LqiYfG6taaGzvoq0zSFcwSH1bFz5/iOc+2tPr46Z4nMwszWLd3mYCQUuy\nx0lNSydzx+USCFqunl3CqcWZBIKWUTnJdAZC/Pr1zVwwuZDzJhSol1ZERESkH+gaVBEZUlp8fkIW\n/MEQy7bUsqexg/MnFbC9po3739pGIGhZMCGfqiYf6/Y2UZyVzPs76mnvZe8shGcmDobCn4nFWcnM\nKM0iyelgY1ULLqfhlrPLCAQtlY0dNLZ3cWpxJpdNH4nXHb41T1OHn7QkF06HUa+tiIiISAwFVBEZ\n9rpD4v4WH62+ALlpSfz69c2MK0hj2/42tta0kp7koq0rQHtXkOKs5B69tR6no8c1tEcycUQ6Lb4A\nexo7epRPLc5gfEE6GV4X00dlsbGqBZ8/yNTiTCYXZfCrxZu5YsZItu1v5ayxeSzdUsMlpxbhchom\nFKb36M2tbGinOCtZPbwiIiIyKCmgioicAGstnYEQSS4HABV17SzfVseiqSN44eM9NHX4OaUwncXr\nq/E4HVgsFXXteN1Olm6uASDZ7aTDH6QwI4nq5s4Tase0kkwWnJLP21tq2VLdQltXkPz0JOaNy2Pl\nzno6uoJ0+kMUZnq5elYJH+yoI2hhUlE618wqYXttG9kpHooyvSS5HHyyp4nJIzMoSPfiMGAtGAMh\nC80dfjKT3T16fDu6giR7nCf/CxURERFBAVVEZMB1h9vuIb8APn+Qv6zZx+SiDByO8KRQo3NTeWZl\nJd+7eCIbq5rZUdNGW1eAGaOyeG97PdXNPjZXt1Lb2smpxZmMzk3BYQwrK+rZ2+Q76XYmu53kpHoA\n2NPYgSMSVMcXpBG0lkDQsqu+HYDzJuSTneKhtTPAKYXp5KZ5eHtLLU0dfr523jg+3t3IS5/s4/rT\nS3E7DbNKs8nwugmEQmSneHA4DJnJbgBaO/9/e3ceHOdZH3D8+9t3L+0lrW7Zkm3Jlg8lcWInJCHE\nOYYkdiCFcrQNwwDlKEehLTDAQOmRaWdaoNN22hmm0EKAMBzhaEJoOJJSSOKQBCfxfcuyZN3XSqvV\n3sfTP95XQnYkxwq2tFF+n5kdvft439139+fn3ef3PM/7bIHneya4cUMtLpfQN5GiPuwnVywRdJLh\nXLHE2HSO2pAXYzjrs1RKKaXUy5cmqEop9TKWzhWZzhaoC/vOKi+WDEcHp3jq1Dgt1QH6JlKkckVe\nv7WJ3SfHiCVzFEuGVK7IVCZP58g0+3onAbi2tZreWIqqgBfBvmZ2PJllXU2QM7EU7fUhmior2NMd\nYzyZI+J3E/C6yRaKTKTm/93bc80ku+cSZ9R2xprqwGwSDFAd9JIvlkhkCmftt21NFde2VjOayHJq\nZBqALU0RUrkihVKJwXiGq1qqMMbuDCiWDJEKD2G/vUj9da013PdUN//6R1fRM57ih8/3MZnK0RwN\nUOGxWFsTYGAyzY3tdcTTebyWC6/bRaTCPbtKtTGGdL7Icz0TrK6qoK0uxOGBOLUhH892T/DaLfVn\nJdIH++JsqA/NOwLdN2GPtkf8dtLudUbqlVJKqZVuWRNUEdkF/BtgAV8xxnzufI/XBFUppS6dVK5A\nwHvhvyqWL5bIFkqEfPY+xhiODibIF0u85+t72Lamig/dsp5HjgxzzdpqDvZN0jk6zXtvbOXhA0OE\n/G7iqRx+r8WZ8RRHB6cI+z0c7I/PvkZTpZ9iyTCSuLAp0NGAx1nd+bfXBM9Nhs9NgM8V9rlJZAsL\nP9t8kjMAABK/SURBVGAe1UEv+UKJnPN5zJjvtTY3hskVSjRW+vn1qfF5X79QshNdgIjfzVSmwJu3\nr+ZAX5zOkWm2Nldy7bpqfnpoiPX1IQrFEhsbwoT9bpLZIiVj+Pqvu3n7dWuc0W7YczpGLJXj5o11\nNEb83Pd0D4lMnk/t3MxkKkcslSMa8HJscIqusSTN0QDb1lSxuTHMzw4Ncf+eXhor/Xzijk3UhLx8\n4WfHedO21UxnC+QKJZ45HeMPrmmmc2SamzfWURPycu/u00SDXjY1hEnmimRyRdbXh6gKeOgaTTI8\nleGNV9mLh+05HeM33TF2tNeRyhVojgZYVxOY91rqM+MpfrSvn5bqAG11QaIBL7FkjqqAB5/bYmw6\ny6bGMB7LRalkELE7WaoC3tnnMMZgDGdNV4+n8xRLhqDPIpMv8X/HhtnRXsfxoQRXtlQRcDoXHj85\nypamyOzsgowTq5DPPfs6xZLBOs/iZ6lcAWPAAG6XvOgMgMlU7qzjn/s+zv2MZtprMz/LNVM/iyXD\nr46PcOum+tn3nSuUfqfOj2S2QCJToLHS/5KfY6VYigXvRhIZfJZFZcDuuJov/i/mQvd5Kc/9SmeM\n3ekb9Omvc/6uli1BFRELOAHcDvQBe4C3GWOOLLSPJqhKKbWyzYxCBrzuFzSQMvkiz5yO4RLY0V5H\nMlvgxHCC1dEKwj4PlkvwWDLb8D8xnGAimWP72ihup+GYzBX55Pf3c2VLFXe/qoW+iTRPnhrj153j\ntDeE6BpN4ve42NgQZudljZSM4TenYzxyeJgtTWFOjSbJForsvKyR3Z1jFEuG5miAY0NTtNYGEYSm\nSj9TmTxul4tHjw6RzZdmE+zmaAV9E2mqAh4mFxht9loubthQw+6TYxTmG2ZeYaIBz7wj7xUei0iF\nG2NgbDpLNOClOujlpDNC/lL43C4qvBapXJGc05lw7iJnL9aJsRCPJeSL9o5Br0VrXZBktsjYdBav\n5aJjVYT+yTRdo8l59796bZSpdB4Ru9OjUDT0TqQIet10jdn7WC5hS1OY3liaxoif48OJs57jitWV\nHOyPE/K5mXY6Wq5qqaIh4qNzZJpTo0na6oI0RwO014f45lM9BH0WH7h5Pc90jdNY6WdTQ5gDfXH2\n9U7y+q1Nsx0pjRE/8XSeI4NTvHZzPXdduYqP3b+PWDLHrZvquGF9LQ/u6yeVKxLxu0lkClyzLsrx\n4Wn2906yvi6IMTCVyVMf9rO+PsTtHQ18dfdp9vdOzv4/qA/7GElk8VjC669oomTgzssbuXpdlAO9\ncb7269M82TnOlqYIW1dXki+VuGxVJaurKni6a5y9ZyYIeN081TXOjRtqiVS4ebLT/m3tHe21bGup\nAuDYUIJCybChPsTzPROcHJmmLuzjzssbyRVKiAghn0W+aEjlCrhE8Fgu9vdN0j2eJJ7Kc3tHIz6P\ni7baIPfuPk0yV+Ty1REaIxU0RHwMTKZJ5oq8+4Z17OmeYFWVnxPDCYansmxtruT4UILNTRGCXovx\nZI6HDwxy25Z6WuuCdI0m+d+jI3SOJGis9DOWyM12XP3Jjlbi6Tzfe7YPgNdd0ciuy5vIFUo8fmKU\nD9zcxvM9E5wYniYa8DCZzrPr8ka+9mQ3Z8ZTfGLnJn60r58zsRSXrYpweGCK27c0cGw4Qcjr5raO\nBu556DDXrIvSHK1gd+c40YCH69tq6GiKkMgUKBnDr46P8mxPjFs31bP3zAQul/DhWzYwOp3lkcND\nvHp9Da9uq2Vf7wSPHBlmKp3H63bZnUu1QX51fJSOVRHa60N43S6S2SJnYincLkHE7kB6uivG+3e0\nsbtzjPb6ELVhH8cGp/jy41189nVb7Prms3j85BgRv4f6sG/2EpOpdJ7W2iC7Lm/ky4910VZnzwKK\nJXN87LaN/OzwED94ro/btjSwtbmSXxwd5o7LGtnTHWNbSxS3JTRHK6gKeHm2O8ZPDg6yoT7Eof4p\n1tYE2NgQZn19iOlMgTOxFEcG4uzvi/PJnZt4/RVNPLivn/ue6uGjt7XTPZYi5HcT8bu5rrWGwXia\nurCP/sk0+WKJHzzXx5Od47xlezNnYklu2VTPXVub2N05xiOHh7nz8kZu2ljHF3/ZSf9kmtFElpKB\n0USGj99ux9Njubh6bZRiyeD32OecZ7piNET8HB9K8Obtq3n85ChttSHecnXz4k9yS2g5E9RXA/cY\nY3Y69z8DYIz5x4X20QRVKaXUy1H/ZJqw303E72FsOktlhQdLhKIxDMUz1Ed8HOqfYk11gLDfjd9j\nzSboE8kc0aCXTqcBHfG7KZYM/3NgkHTeTrR6YymCzk8XRfxuHjsxyl/d1UFlhYcH9/bj81hsXV1J\nLJnjB8/1cevmejbUh3j8xChVTsMzkSmQyOQJeN1EKtw8tG+AkUSWYsngtVxsaQqTKZSwRAj63Eyk\nclzZXEXQZ+FzW3z/uV48lovBeIbjQ1OcGJ5mQ32Im9rrODY0RSJToFAyjE9ned+OVvaemeSnh4b4\nvStXMZrI0BINUB30srW5itFEhsMDUwzGM1QGPAzHM1R4LYyBSIUbj+Xi2GACEUjlilzZUsWP9w8s\n+PkHvBYBr0XHqkpqQ14CXovByQxPdY2/4Cem3rxtNYcG4pwYnuZV66Jsbozw8MFBYsnceWN8bmK7\nudEe1d7TPUFLtT0NPJktks4VZxONC1Ed9J712pUVHuLphafS+z0uMvkXX1V8sbxu12xS/1LMTZqV\nUsvriU/dSkt1YLkPY0HLmaC+FdhljHmfc/8dwHXGmI+c87j3A+8HWLNmzdU9PT0X9TiUUkop9fLX\nP5mmJujFa7lI5gqEnet3z2dmCvDx4QRul9BaG5p3am4mX8QlgtftIpUr4LFcCPaIpojQG0vRWOnH\nY71wumy+aCf1M9M/Z9pTuWIJn9uaPYbpbIGRRJZD/XHu6GgkmSvQN5HmqpYqEpk8LqdjAOwpwpZL\n8LhczEwymEjl8VhC2G8nsPFUnr7JFO31YU6PJSmUSjRVVhDxuxERvvlUD7d11JPMFhlNZFlV5Sdf\nNOQKJa5rs6/n3tMdw3IJuy5rxG25iKfzjCayuF3Cb07HaKsLEqmwZy+cHE7Q3hCmNuTjYF+coakM\nO9prqQl6KRqDz23ROTJNPJ0j4vfwZOcYV62J4ve4qA/7efjAAOtq7dHdsN9NtlCiIeyjZOAffnIU\nr9vFmuoAuy5vBGA0kaW1Nsj+3knGkzkqPBYt1QECXouSsS8LcIkwPJVhXU2QJ06O0tFkjxauqQng\ntVykckWuWF3JeDJLlbOa+eGBOHvPTHLTxjpODk9zZDDO5sYIa6oD9E6kaIkGqAp48Fgu3JZwuH+K\nkUSGkM9DIpMn6HPz59/dy4duXo+I3aGwqqqCa1ur+dKvupjK5LnnDZcxlc5zdHCKgck0GxvCrK0J\ncmLY/i3tgckMa2sCTCRzTKbyjCQybGmyV1ff1BjmTCzFj/cP0BDxc31bNR7L5bxfeOzEKOl8kZvb\n66gL+xibzlEd9PLwgQHetL2Z6oCXB/b2kSuWuGvrKr7/bB+v2VBDe73doXLfU930xFJ84Kb1VAU8\n3L+nl4aIj74J++fREpkC0YAHgz1if2RwioN9cdrqQty8sY6xaXu2yEgiS/9EmoaIj0iFhxvba3mm\nK0aFx+JVrVGS2SI940kM0BtLkc4VGYhnZkd0333DOoI+N6fHknzliS5u62jAYwmD8QzRgJeOpghn\nYim8bhf7eye5bFUl17VV8+DefjpHp2mtCXLnFfaaCyeGE7xmQy3Zgt2h11Id4FB/nCtWV2KAWDLH\n+HSOK1sq+fJjXVyzLkpbbYgKr4tEpsAvjo4Q8rvZ0hTBGHvWzE8ODlId9OL3WNSEvJwamaY+4mdN\ndYBnu2P0xlK864Z1PHFyjMMDcXa01zE8leGJk2Nsbgpz66Z6Ehl7FlDE78ZtuXjr1c2kckUSmTxf\neeI0LhFu2VRHfcTHL46OkM3bl0mE/W6++MtT/PVdHaRyBTD2+W9VVQXPdI3TsaqSoM/i4QODeN32\n8/aMp/iz7+zlkzs3cX1bNVevrb7wk+syKPsEdS4dQVVKKaWUUuVMr99U6ndzoQnqpVg+sB9omXO/\n2SlTSimllFLqZUmTU6WWxqVIUPcA7SLSKiJe4G7goUvwOkoppZRSSimlVpCLvl6yMaYgIh8Bfo79\nMzP3GmMOX+zXUUoppZRSSim1slyS30Fd9EGIjALlvkpSLTC23AehzqIxKU8al/KjMSlPGpfyozEp\nTxqX8qMxKU/lHpe1xpi6F3tQWSSoLwci8uyFXNSrlo7GpDxpXMqPxqQ8aVzKj8akPGlcyo/GpDyt\nlLhcimtQlVJKKaWUUkqpRdMEVSmllFJKKaVUWdAE9cL953IfgHoBjUl50riUH41JedK4lB+NSXnS\nuJQfjUl5WhFx0WtQlVJKKaWUUkqVBR1BVUoppZRSSilVFjRBfREisktEjotIp4h8ermP55VCRFpE\n5JcickREDovIXzjl94hIv4jsc26vm7PPZ5w4HReRnct39CubiHSLyEHn83/WKasWkUdF5KTzN+qU\ni4j8uxOXAyKyfXmPfuURkU1z6sM+EZkSkY9qXVl6InKviIyIyKE5ZYuuGyLyLufxJ0XkXcvxXlaS\nBeLyTyJyzPnsHxCRKqd8nYik59SbL83Z52rn3NfpxE6W4/2sBAvEZNHnLG2jXVwLxOX+OTHpFpF9\nTrnWlSVwnvbwyv5uMcbobYEbYAGngDbAC+wHOpb7uF4JN6AJ2O5sh4ETQAdwD/CJeR7f4cTHB7Q6\ncbOW+32sxBvQDdSeU/YF4NPO9qeBzzvbrwN+CghwPfDMch//Sr4556whYK3WlWX5/G8CtgOH5pQt\nqm4A1UCX8zfqbEeX+729nG8LxOUOwO1sf35OXNbNfdw5z/MbJ1bixO7O5X5vL9fbAjFZ1DlL22hL\nE5dz/v2fgb9xtrWuLE1MFmoPr+jvFh1BPb9rgU5jTJcxJgd8F3jjMh/TK4IxZtAY87yznQCOAqvP\ns8sbge8aY7LGmNNAJ3b81NJ4I/ANZ/sbwO/PKb/P2J4GqkSkaTkO8BXitcApY0zPeR6jdeUSMcY8\nDsTOKV5s3dgJPGqMiRljJoBHgV2X/uhXrvniYox5xBhTcO4+DTSf7zmc2ESMMU8bu7V3H7+NpVqk\nBerKQhY6Z2kb7SI7X1ycUdA/BL5zvufQunJxnac9vKK/WzRBPb/VQO+c+32cP0lSl4CIrAO2Ac84\nRR9xpi3cOzOlAY3VUjLAIyLynIi83ylrMMYMOttDQIOzrXFZWndzduNB68ryW2zd0PgsvfdgjzjM\naBWRvSLymIjscMpWY8dihsbl0ljMOUvrytLaAQwbY07OKdO6soTOaQ+v6O8WTVBVWROREPBD4KPG\nmCngP4D1wFXAIPZ0E7W0bjTGbAfuBD4sIjfN/Uenx1SXB19iIuIF3gB83ynSulJmtG6UHxH5LFAA\nvuUUDQJrjDHbgI8D3xaRyHId3yuMnrPK29s4uwNU68oSmqc9PGslfrdognp+/UDLnPvNTplaAiLi\nwa6M3zLG/DeAMWbYGFM0xpSA/+K3UxM1VkvEGNPv/B0BHsCOwfDM1F3n74jzcI3L0rkTeN4YMwxa\nV8rIYuuGxmeJiMgfA3cBb3caeDjTSMed7eewr3HciB2DudOANS4X2Us4Z2ldWSIi4gbeDNw/U6Z1\nZenM1x5mhX+3aIJ6fnuAdhFpdUYn7gYeWuZjekVwrnX4KnDUGPMvc8rnXr/4JmBmpbmHgLtFxCci\nrUA79kX66iISkaCIhGe2sRcaOYT9+c+sCPcu4EfO9kPAO51V5a4H4nOmpKiL66zeba0rZWOxdePn\nwB0iEnWmON7hlKmLSER2AZ8C3mCMSc0prxMRy9luw64fXU5spkTkeuf76Z38NpbqIngJ5yxtoy2d\n24BjxpjZqbtaV5bGQu1hVvh3i3u5D6CcGWMKIvIR7ABawL3GmMPLfFivFK8B3gEcFGdJc+AvgbeJ\nyFXYUxm6gQ8AGGMOi8j3gCPY07U+bIwpLvlRr3wNwAP2+RI38G1jzM9EZA/wPRF5L9CDvZACwE+w\nV5TrBFLAu5f+kFc+p7Pgdpz64PiC1pWlJSLfAW4BakWkD/hb4HMsom4YY2Ii8vfYjW+AvzPGXOhi\nMmoeC8TlM9irwj7qnM+eNsZ8EHsV078TkTxQAj445/P/U+DrQAX2Natzr1tVi7BATG5Z7DlL22gX\n13xxMcZ8lReubwBaV5bKQu3hFf3dIs6sFqWUUkoppZRSalnpFF+llFJKKaWUUmVBE1SllFJKKaWU\nUmVBE1SllFJKKaWUUmVBE1SllFJKKaWUUmVBE1SllFJKKaWUUmVBE1SllFJKKaWUUmVBE1SllFJK\nKaWUUmVBE1SllFJKKaWUUmXh/wGCTi/FL42njgAAAABJRU5ErkJggg==\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig, ax = plt.subplots(2, 1, figsize=(16, 4))\n", + "ax[0].plot(hist_x, label='X coordinate');\n", + "ax[1].plot(opt.hist, label='Y coordinate');\n", + "ax[0].legend();\n", + "ax[1].legend();" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.4.3" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/pymc3/__init__.py b/pymc3/__init__.py index 9f9afda7e4..259cefc240 100644 --- a/pymc3/__init__.py +++ b/pymc3/__init__.py @@ -27,6 +27,7 @@ from .tests import test from .data import * +from .optimization import * import logging _log = logging.getLogger('pymc3') diff --git a/pymc3/optimization.py b/pymc3/optimization.py new file mode 100644 index 0000000000..90b7303083 --- /dev/null +++ b/pymc3/optimization.py @@ -0,0 +1,83 @@ +import tqdm +import numpy as np +import theano +from theano.configparser import change_flags +import pymc3 as pm + +__all__ = [ + 'Optimizer' +] + + +class Optimizer(object): + """ + Optimization with posterior replacements + + Parameters + ---------- + approx : Approximation + loss : scalar + wrt : shared params + more_replacements : other replacements in the graph + optimizer : callable that returns updates, pm.adam by default + """ + @change_flags(compute_test_value='off') + def __init__(self, approx, loss, wrt, more_replacements=None, optimizer=None): + if optimizer is None: + optimizer = pm.adam + self.optimizer = optimizer + self.wrt = wrt + self.approx = approx + self.loss = approx.apply_replacements(loss, more_replacements=more_replacements) + updates = self.optimizer(self.loss, self.wrt) + self.step_function = theano.function([], self.loss, updates=updates) + self.hist = np.asarray(()) + + @change_flags(compute_test_value='off') + def refresh(self, kwargs=None): + """ + Recompile step_function and reset updates + + This can be needed sometimes when + you use stateful optimizers like ADAM + and want to reuse the function for new problem + + Parameters + ---------- + kwargs : kwargs for theano.function + """ + updates = self.optimizer(self.loss, self.wrt) + self.step_function = theano.function([], self.loss, updates=updates, **kwargs) + self.hist = np.asarray(()) + + def fit(self, n=5000, callbacks=()): + """ + Perform optimization steps + + Parameters + ---------- + n : int + number of iterations + callbacks : list[callable] + list of callables with following signature + f(Approximation, loss_history, i) -> None + """ + progress = tqdm.trange(n) + scores = np.empty(n) + scores[:] = np.nan + i = 0 + try: + for i in progress: + scores[i] = self.step_function() + for callback in callbacks: + callback(self.approx, scores[:i + 1], i) + if i % ((n+1000)//1000) == 0: + progress.set_description( + 'E_q[Loss] = %.4f' % scores[max(0, i - n//50):i+1].mean() + ) + except (KeyboardInterrupt, StopIteration): + self.hist = np.concatenate((self.hist, scores[:i])) + else: + self.hist = np.concatenate((self.hist, scores)) + finally: + progress.close() diff --git a/pymc3/variational/opvi.py b/pymc3/variational/opvi.py index 71b2e6dda9..0a594b62c0 100644 --- a/pymc3/variational/opvi.py +++ b/pymc3/variational/opvi.py @@ -593,11 +593,11 @@ def apply_replacements(self, node, deterministic=False, Parameters ---------- - node : Theano Variables (or Theano expressions) + node : Theano Variable(s) node or nodes for replacements deterministic : bool - whether to use zeros as initial distribution - if True - zero initial point will produce constant latent variables + whether to use point with highest density as initial point for distribution + if True - constant initial point will produce constant latent variables include : list latent variables to be replaced exclude : list @@ -623,7 +623,7 @@ def sample_node(self, node, size=100, Parameters ---------- - node : Theano Variables (or Theano expressions) + node : Theano Variable(s) size : scalar number of samples more_replacements : dict