{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Inference plots - Pairwise scatterplots\n", "\n", "This example builds on the [adaptive covariance MCMC example](https://github.com/pints-team/pints/blob/master/examples/sampling-adaptive-covariance-mcmc.ipynb), and shows you a different way to plot the results.\n", "\n", "Inference plots:\n", "* [Predicted time series](https://github.com/pints-team/pints/blob/master/examples/plot-mcmc-predicted-time-series.ipynb)\n", "* [Trace plots](https://github.com/pints-team/pints/blob/master/examples/plot-mcmc-trace-plots.ipynb)\n", "* [Autocorrelation](https://github.com/pints-team/pints/blob/master/examples/plot-mcmc-autocorrelation.ipynb)\n", "* __Pairwise scatterplots__\n", "* [Pairwise scatterplots with KDE](https://github.com/pints-team/pints/blob/master/examples/plot-mcmc-pairwise-kde-plots.ipynb)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Setting up an MCMC routine\n", "\n", "See the adaptive covariance MCMC example for details." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from __future__ import print_function\n", "import pints\n", "import pints.toy as toy\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "# Load a forward model\n", "model = toy.LogisticModel()\n", "\n", "# Create some toy data\n", "real_parameters = [0.015, 500]\n", "times = np.linspace(0, 1000, 100)\n", "org_values = model.simulate(real_parameters, times)\n", "\n", "# Add noise\n", "noise = 50\n", "values = org_values + np.random.normal(0, noise, org_values.shape)\n", "real_parameters = np.array(real_parameters + [noise])\n", "\n", "# Get properties of the noise sample\n", "noise_sample_mean = np.mean(values - org_values)\n", "noise_sample_std = np.std(values - org_values)\n", "\n", "# Create an object with links to the model and time series\n", "problem = pints.SingleOutputProblem(model, times, values)\n", "\n", "# Create a log-likelihood function (adds an extra parameter!)\n", "log_likelihood = pints.GaussianLogLikelihood(problem)\n", "\n", "# Create a uniform prior over both the parameters and the new noise variable\n", "log_prior = pints.UniformLogPrior(\n", " [0.01, 400, noise*0.1],\n", " [0.02, 600, noise*100]\n", " )\n", "\n", "# Create a posterior log-likelihood (log(likelihood * prior))\n", "log_posterior = pints.LogPosterior(log_likelihood, log_prior)\n", "\n", "# Perform sampling using MCMC, with a single chain\n", "x0 = real_parameters * 1.1\n", "mcmc = pints.MCMCController(log_posterior, 1, [x0])\n", "mcmc.set_max_iterations(6000)\n", "mcmc.set_log_to_screen(False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plotting 1d and 2d histograms\n", "\n", "We can now run the MCMC routine and plot the histograms of the inferred parameters." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Running...\n", "Done!\n" ] } ], "source": [ "print('Running...')\n", "chains = mcmc.run()\n", "print('Done!')" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Select chain 0 and discard warm-up\n", "chain = chains[0]\n", "chain = chain[3000:]" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEKCAYAAADXdbjqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAGVxJREFUeJzt3XuwZWV55/HvT1AuGgWkIT0INhi8oKONHokO4qB4QTQixguMUVS0JcFxLDM1ImaUMrGGJCLRcYK2JREcQFRESYkXZFRi1XBpEBsQkIsdbelqWnQABSHAM3+sdWDTnHN6n7POvnV/P1W79lrvXpenV+9znvOu913vm6pCkqQuHjHqACRJk89kIknqzGQiSerMZCJJ6sxkIknqzGQiSerMZCJJ6sxkIknqzGQiSeps61EH0MXOO+9cy5YtG3UYkjRRLrvssl9V1ZLFPOZEJ5Nly5axatWqUYchSRMlyb8u9jG9zSVJ6sxkIknqzGQiSerMZCJJ6sxkIknqzGQiSerMZCJJ6sxkIknqzGQiSepsop+AlybFsmO/MWP5mhNeOeRIpMGwZiJJ6mxgySTJ7km+l+SaJFcn+S9t+U5Jzk9yffu+Y1ueJJ9MckOS1UmePajYJEmLa5A1k3uBv6yqpwHPA45Jsg9wLHBBVe0NXNCuA7wC2Lt9rQBOHmBskqRFNLBkUlXrqurydvkO4BpgN+BQ4NR2s1OB17TLhwKnVeMiYIckSwcVnyRp8QylAT7JMmBf4GJg16paB03CSbJLu9luwC96dlvblq0bRozSKNgwr83FwBvgkzwGOBt4b1XdPtemM5TVDMdbkWRVklUbNmxYrDAlSR0MNJkkeSRNIjm9qr7aFq+fvn3Vvt/Slq8Fdu/Z/QnAzRsfs6pWVtVUVU0tWbKoE4VJkhZokL25AnwOuKaqPt7z0bnAke3ykcDXe8rf0vbqeh5w2/TtMEnSeBtkm8n+wJuBK5Nc0ZYdB5wAfCnJUcDPgde3n50HHALcANwJvG2AsUmSFtHAkklV/ZCZ20EADpph+wKOGVQ8kqTB8Ql4SVJnJhNJUmcO9CjpIWZ79gV8/kWzs2YiSerMZCJJ6sxkIknqzGQiSerMZCJJ6sxkIknqzGQiSerMZCJJ6sxkIknqzGQiSerMZCJJ6sxkIknqzGQiSepskNP2npLkliRX9ZSdleSK9rVmegbGJMuS3NXz2acHFZckafENcgj6zwOfAk6bLqiqN04vJzkRuK1n+xuravkA45EkDcggp+29MMmymT5LEuANwIsHdX5J0vCMqs3kAGB9VV3fU7Znkh8l+UGSA0YUlyRpAUY10+IRwJk96+uAParq1iTPAb6W5OlVdfvGOyZZAawA2GOPPYYSrCRpbkOvmSTZGngtcNZ0WVXdXVW3tsuXATcCT55p/6paWVVTVTW1ZMmSYYQsSdqEUdzmeglwbVWtnS5IsiTJVu3yXsDewE0jiE2StACD7Bp8JvB/gackWZvkqPajw3noLS6AFwKrk/wY+ApwdFX9elCxSZIW1yB7cx0xS/lbZyg7Gzh7ULFIkgbLJ+AlSZ2ZTCRJnZlMJEmdmUwkSZ2ZTCRJnZlMJEmdmUwkSZ2NamwuaeiWHfuNGcvXnPDKIUcibX6smUiSOrNmIk2Q2WpXYA1Lo2XNRJLUmclEktSZyUSS1JltJtIWaq72F2m+rJlIkjozmUiSOhvkTIunJLklyVU9Zccn+WWSK9rXIT2ffSDJDUmuS/LyQcUlSVp8g6yZfB44eIbyk6pqefs6DyDJPjTT+T693ecfp+eElySNv4Elk6q6EOh3HvdDgS9W1d1V9TPgBmC/QcUmSVpco2gzeXeS1e1tsB3bst2AX/Rss7YtkyRNgGEnk5OBJwHLgXXAiW15Zti2ZjpAkhVJViVZtWHDhsFEKUmal6Emk6paX1X3VdX9wGd58FbWWmD3nk2fANw8yzFWVtVUVU0tWbJksAFLkvoy1GSSZGnP6mHAdE+vc4HDk2yTZE9gb+CSYcYmSVq4gT0Bn+RM4EBg5yRrgQ8DByZZTnMLaw3wLoCqujrJl4CfAPcCx1TVfYOKTZK0uAaWTKrqiBmKPzfH9h8FPjqoeCRJg+MT8JKkzkwmkqTOTCaSpM5MJpKkzkwmkqTOnBxL2sw5CZaGwWSiiTTXL8g1J7xyiJFIAm9zSZIWgclEktSZyUSS1JnJRJLUmclEktSZyUSS1FlfySTJMwYdiCRpcvVbM/l0kkuS/EWSHQYakSRp4vSVTKrqBcCbaKbWXZXkjCQvHWhkkqSJ0XebSVVdD/wV8H7gPwKfTHJtktfOtH2SU5LckuSqnrK/b/dZneSc6VpOkmVJ7kpyRfv6dLd/liRpmPptM3lmkpOAa4AXA39SVU9rl0+aZbfPAwdvVHY+8IyqeibwU+ADPZ/dWFXL29fR8/g3SJJGrN+ayaeAy4FnVdUxVXU5QFXdTFNbeZiquhD49UZl36mqe9vVi4AnLChqSdJY6TeZHAKcUVV3ASR5RJLtAarqCws899uBb/as75nkR0l+kOSABR5TkjQC/SaT7wLb9axv35YtSJIPAvcCp7dF64A9qmpf4H3AGUkeO8u+K5KsSrJqw4YNCw1BkrSI+k0m21bVb6dX2uXtF3LCJEcCrwLeVFXVHu/uqrq1Xb4MuBF48kz7V9XKqpqqqqklS5YsJARJ0iLrN5n8Lsmzp1eSPAe4a74nS3IwTW+wV1fVnT3lS5Js1S7vBewN3DTf40uSRqPfybHeC3w5yc3t+lLgjXPtkORM4EBg5yRrgQ/T9N7aBjg/CcBFbc+tFwIfSXIvcB9wdFX9esYDS5LGTl/JpKouTfJU4ClAgGur6t82sc8RMxR/bpZtzwbO7icWSdL4mc+0vc8FlrX77JuEqjptIFFJkiZKX8kkyReAJwFX0NyGAijAZCINwFxz3EvjqN+ayRSwz3TvK2lLMNsv9DUnvHLIkUjjr9/eXFcBfzjIQCRJk6vfmsnOwE+SXALcPV1YVa8eSFSS5s1bYxqlfpPJ8YMMQpI02frtGvyDJE8E9q6q77bjcm012NCk4fAveqm7foegfyfwFeAzbdFuwNcGFZQkabL02wB/DLA/cDs8MFHWLoMKSpI0WfpNJndX1T3TK0m2pnnORJKkvpPJD5IcB2zXzv3+ZeCfBxeWJGmS9JtMjgU2AFcC7wLOY5YZFiVJW55+e3PdD3y2fUmS9BD9js31M2ZoI6mqvRY9IknSxJnP2FzTtgVeD+y0+OFIkiZRX20mVXVrz+uXVfUPwIsHHJskaUL0e5vr2T2rj6CpqfzBQCKSNHEcYVn93uY6sWf5XmAN8IZN7ZTkFOBVwC1V9Yy2bCfgLJqJttYAb6iq36SZx/cTwCHAncBbq+ryPuOTJI1Qv7e5XtTzemlVvbOqrutj188DB29UdixwQVXtDVzQrgO8Ati7fa0ATu4nNknS6PV7m+t9c31eVR+fpfzCJMs2Kj4UOLBdPhX4PvD+tvy0dgKui5LskGRpVa3rJ0ZJ0uj0+9DiFPDnNAM87gYcDexD024y37aTXacTRPs+PcbXbsAverZb25Y9RJIVSVYlWbVhw4Z5nlqSNAjzmRzr2VV1B0CS44EvV9U7FjGWzFA207MtK4GVAFNTU44PJkljoN9ksgdwT8/6PTQN6Auxfvr2VZKlwC1t+Vpg957tngDcvMBzaDOxkLlGnJ9k/Nn7a/PT722uLwCXJDk+yYeBi4HTFnjOc4Ej2+Ujga/3lL8ljecBt9leIkmTod+xuT6a5JvAAW3R26rqR5vaL8mZNI3tOydZC3wYOAH4UpKjgJ/TPE0PzeCRhwA30HQNfts8/h2SpBHq9zYXwPbA7VX1T0mWJNmzqn421w5VdcQsHx00w7ZFMwmXJGnC9Dtt74dpuu9+oC16JPC/BxWUJGmy9FszOQzYF7gcoKpuTuJwKpLmZGeILUe/DfD3tLehCiDJowcXkiRp0vSbTL6U5DPADkneCXwXJ8qSJLX67c31sXbu99uBpwAfqqrzBxqZJGlibDKZJNkK+HZVvQQwgUiSHmaTyaSq7ktyZ5LHVdVtwwhKGmdzNSr7BLe2VP325vo9cGWS84HfTRdW1XsGEpUkaaL0m0y+0b4kSXqYOZNJkj2q6udVdeqwApIkTZ5NdQ3+2vRCkrMHHIskaUJtKpn0zjGy1yADkSRNrk0lk5plWZKkB2yqAf5ZSW6nqaFs1y7TrldVPXag0UmSJsKcyaSqthpWINLmwIENtaXqd2wuSZJmNZ/JsRZFkqcAZ/UU7QV8CNgBeCewoS0/rqrOG3J4kqQFGHoyqarrgOXwwLhfvwTOoZmm96Sq+tiwY5IkdTPq21wHATdW1b+OOA5JUgejTiaHA2f2rL87yeokpyTZcaYdkqxIsirJqg0bNsy0iSRpyNJMoDiCEyePAm4Gnl5V65PsCvyK5nmWvwaWVtXb5zrG1NRUrVq1avDBamTsHbVlcdTl4UhyWVVNLeYxR1kzeQVweVWtB6iq9VV1X1XdTzOL434jjE2SNA+jTCZH0HOLK8nSns8OA64aekSSpAUZem8ugCTbAy8F3tVT/HdJltPc5lqz0WeSpDE2kmRSVXcCj9+o7M2jiEWS1N2oe3NJkjYDJhNJUmcmE0lSZyYTSVJnJhNJUmcmE0lSZyYTSVJnJhNJUmcjeWhRm7/ZBmh0ID9p82TNRJLUmclEktSZyUSS1JltJhoq21KkzZM1E0lSZyYTSVJnJhNJUmcjazNJsga4A7gPuLeqppLsBJwFLKOZbfENVfWbUcWo4ZmtLUXSZBh1zeRFVbW8qqba9WOBC6pqb+CCdl2SNOZGnUw2dihwart8KvCaEcYiSerTKJNJAd9JclmSFW3ZrlW1DqB932Vk0UmS+jbK50z2r6qbk+wCnJ/k2n52ahPPCoA99thjkPFJkvo0sppJVd3cvt8CnAPsB6xPshSgfb9lhv1WVtVUVU0tWbJkmCFLkmYxkppJkkcDj6iqO9rllwEfAc4FjgROaN+/Por4JI3GXL36HCVhvI3qNteuwDlJpmM4o6q+leRS4EtJjgJ+Drx+RPFJkuZhJMmkqm4CnjVD+a3AQcOPSJLUxbh1DZYkTSCTiSSpM4eglzQRnL5gvFkzkSR1ZjKRJHVmMpEkdWYykSR1ZjKRJHVmMpEkdWYykSR1ZjKRJHVmMpEkdWYykSR1ZjKRJHVmMpEkdWYykSR1NvRRg5PsDpwG/CFwP7Cyqj6R5HjgncCGdtPjquq8YccnabI4mvB4GMUQ9PcCf1lVlyf5A+CyJOe3n51UVR8bQUySpA6Gnkyqah2wrl2+I8k1wG7DjkPS5s0ay3CNtM0kyTJgX+DitujdSVYnOSXJjiMLTJI0LyNLJkkeA5wNvLeqbgdOBp4ELKepuZw4y34rkqxKsmrDhg0zbSJJGrKRJJMkj6RJJKdX1VcBqmp9Vd1XVfcDnwX2m2nfqlpZVVNVNbVkyZLhBS1JmtXQk0mSAJ8Drqmqj/eUL+3Z7DDgqmHHJklamFH05tofeDNwZZIr2rLjgCOSLAcKWAO8awSxaQY2ZEralFH05vohkBk+8pkSSZpQPgEvSerMZCJJ6sxkIknqzGQiSepsFL25tJmYrZeXpC2PNRNJUmcmE0lSZyYTSVJnJhNJUmcmE0lSZ/bmmnCLOW6WvbO0JXCsucGwZiJJ6syaiSSNsUmpSZlMtkDezpIebiE/F7P9Qp+UBLCYTCYTwgQgaZzZZiJJ6mzskkmSg5Ncl+SGJMeOOh5J0qaN1W2uJFsB/wt4KbAWuDTJuVX1k0Gcb3O+r+ltMWnw/Dl70LjVTPYDbqiqm6rqHuCLwKEjjkmStAljVTMBdgN+0bO+FvjjEcXyMMN4QHBzqBVJ2vKkqkYdwwOSvB54eVW9o11/M7BfVf3nnm1WACva1acA183zNDsDv1qEcIfJmIfDmIfDmIdjrpifWFVLFvNk41YzWQvs3rP+BODm3g2qaiWwcqEnSLKqqqYWuv8oGPNwGPNwGPNwDDvmcWszuRTYO8meSR4FHA6cO+KYJEmbMFY1k6q6N8m7gW8DWwGnVNXVIw5LkrQJY5VMAKrqPOC8AZ5iwbfIRsiYh8OYh8OYh2OoMY9VA7wkaTKNW5uJJGkSVdXYv4CDaboA3wAcO8Pn2wBntZ9fDCxryx8PfA/4LfCpjfb5FvBj4Grg08BWbflZwBXtaw1wRVu+DLir57NPDzvmnn3PBa7qWd8JOB+4vn3fsS0P8Mn2HKuBZ49RzH8PXNvGdQ6wwwRc5+OBX/bEdkjPZx9oz3EdTff2cYl5bL/PwPfbY07HsMtcxxqH6zxHzO8DfkLzfb6Apuvt9D739Wx/7hjF/FZgQ0/5O3r2OZLm98n1wJFzxfzAPv1sNMoXTUP8jcBewKNoEsA+G23zF9M/DDQ9wM5qlx8NvAA4eoYL/Nj2PcDZwOEznPtE4EM9P3xXjTLm9vPXAmfw0F8Yfzf9xQOOBf62XT4E+Gb7b3wecPEYxfwyYOt2+W97Yh7n63w88F9n2Haf9tzbAHu2MW01DjGP8/eZ5pfc1Aznm+1YI7/Oc8T8ImD7dvnPp4/Vrv92TK/zW2f5Hu0E3NS+79gu77ip+CfhNlc/Q6wcCpzaLn8FOChJqup3VfVD4PcbH7Sqbm8Xt6b5j3tI41GSAG8AzhyXmJM8huYvoL+Z41inAq/pKT+tGhcBOyRZOg4xV9V3quredvUimmeK5mvY13k2hwJfrKq7q+pnNH817jdOMY/j93kOMx6LMbjOs6mq71XVne3qWH2fF+DlwPlV9euq+g3N3Y6DN7XTJCSTmYZY2W22bdpfULfRVPvmlOTbwC3AHTT/Mb0OANZX1fU9ZXsm+VGSHyQ5YAQx/zXNX5d3blS+a1Wta4+1DthlHnGMKuZeb6epQU0b1+sM8O4kq5OckmTHecQxyphhPL/PAP+U5Iok/71NGHMdaxyu82wx9zqKh36ft02yKslFSV4zw/ajjPlP2+/zV5JMPzA+n+v8gElIJjP9Z23cBa2fbR6+QdXLgaU01eYXb/TxETz0r7h1wB5VtS/NX4BnJHnssGJOshz4o6o6Z7ZtFhjHfLZd9JiTfBC4Fzi9LRrn63wy8CRgeRvniQs4x6i+G2P1fW69qar+PU2iOwB48yaONdLr3Jot5uaAyZ8BUzRtgtP2qOZJ9P8E/EOSJ41JzP9M0+byTOC7PFjjWdDv00lIJpscYqV3myRbA48Dft3Pwavq9zSNlg9UJ9tjvJamoWt6u7ur6tZ2+TKae5tPHmLMzweek2QN8EPgyUm+3362fvr2Vft+yzziGFXMJDkSeBXNl71gvK9zVa2vqvuq6n7gszx4i2Xcr/M4fp+pql+273fQtPU87HpudKxRX+e5YibJS4APAq+uqrt79rm5fb+Jpv1i33GIuapu7Ynzs8Bz5hHHw0xCMulniJVzaXofALwO+D/Tv5xmkuQxPb98t6ZpqL62Z5OXANdW1dqefZa0862QZC9gb5qGqaHEXFUnV9W/q6plNA1tP62qA2c41pHA13vK35LG84Dbpm+HjTrmJAcD76f5wXvg1sw4X+eN2psOA67qOcfhSbZJsmcb8yXjEHNr7L7PSbZOsnO7/EiaPyp6r+dMxxrpdZ4r5iT7Ap+h+T7f0rPPjkm2aZd3Bvan6fU1DjH3fp9fDVzTLn8beFkb+440nWW+Pds5HlB99DIY9Yvml/1Paf56+mBb9hGa/ziAbYEv0zTIXQLs1bPvGprM/VuajLsPsGv7H7eapmvw/6TtWdTu83ng6I1i+NN22x8DlwN/MsyYNzr2Mh7ay+jxNN0Rr2/fd2rLQzPZ2I3AlczQo2OEMd9Ac1/2IV1Tx/w6f6G9jqtpfqiX9nz2wTaG64BXjEvM4/p9pul9dBkP/gx+gge75891rJFd503E/F1gPRt1AQb+Q/ud+XH7ftQYxfw/er4D3wOe2nOst7fnuAF421wxT798Al6S1Nkk3OaSJI05k4kkqTOTiSSpM5OJJKkzk4kkqTOTiTZ7Se5rh5K4KsmXk2w/6pgAkhy3CMd4fZKrk9yfZKLmKNfmxWSiLcFdVbW8qp4B3EMzsmpfph/sG5B5J5MZ4rmK5un2CxclImmBTCba0vwL8EcASb6W5LL2L/sV0xsk+W2SjyS5GHh+kg8lubSt2aycHigvyfeTnJTkwiTXJHlukq8muT7J3/Qc78+SXNLWjj6TZKskJwDbtWWnz7bdTPH0/mOq6pqqum7QF03aFJOJthjt0DmvoHkSGeDtVfUcmoH53pNkevTVR9M8Rf7H1Qzr/amqem5bs9mOZkiKafdU1QtpJlj7OnAM8AzgrUken+RpwBuB/atqOc1ESW+qqmN5sMb0ptm2myUeaexsPeoApCHYLskV7fK/AJ9rl9+T5LB2eXeasZ5upflFfnbP/i9K8t+A7WkmDLqaZsRVeHDspCuBq6sd+yzJTe0xX0AzgN6lbYVmOx4ciLPXQXNst3E80tgxmWhLcFf71/4DkhxIMwDi86vqzjSj7G7bfvz7qrqv3W5b4B9pxjX7RZLje7YDmB519f6e5en1rWnGRzu1qj6wiRjn2u6BeKRx5W0ubakeB/ymTSRPpZnWeCbTieNXaWYzfN08z3MB8LokuwAk2SnJE9vP/q0dyXVT20ljz2SiLdW3gK2TrKaZpfCimTaqqv9HM9fDlcDXaEab7ltV/QT4K+A77bnOp5mQDWAlsDrJ6ZvYblZJDkuylqZh/htpZg+Vhs5RgyVJnVkzkSR1ZjKRJHVmMpEkdWYykSR1ZjKRJHVmMpEkdWYykSR1ZjKRJHX2/wGvXQbDq73/HgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.figure()\n", "plt.xlabel('Parameter 1')\n", "plt.ylabel('Frequency')\n", "bins = np.linspace(np.min(chain[:,0]), np.max(chain[:,0]), 50)\n", "plt.hist(chain[:,0], bins=bins)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plotting a histogram of two variables (showing their correlation) is a bit more involved. We can use Matplotlib's [hist2d method](http://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.hist2d), but we need to do some extra work to get a fixed aspect ratio." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAARMAAAEKCAYAAAAxcLHrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHKxJREFUeJzt3Xu0lNWZ5/Hv4w1ERRQQUa4i2uMlUSHGS7Al6fEe6bRmgmMiRlds02a6p9OJ0THJ2Kbt1Y6J3U5cjSHdSbSjLXaMShx1YhLNxYgJICqoICoKwiAQRW6C4DN/7F1SHM45ex/Yb50qzu+z1lmn6q2n3r2rTtVz3st+92PujojIjtqluzsgIjsHJRMRKULJRESKUDIRkSKUTESkCCUTESlCyUREilAyEZEilExEpIjdursDO2LAgAE+fPiIhrT1XkZMs2Xm9zIGN+9i1fdDWtusWTNXuPvAVFylycTMFgKrgc3AJncfa2bfACYQvp9vABe7+xIzM+Bm4CxgXVw+q7P1Dx8+gsefnFHlS3jfxk3pdLLHbs2VTtZv3JyM2XOPXRvQE2lle+5ur+bENeLTP97dj3H3sfH+je7+AXc/BngA+HpcfiYwOv5cBkxuQN9EpJCG/yt197fr7u4F1DbGJwC3ezAd6GdmgxvdPxHZPlUnEwd+amYzzeyy2kIzu97MFgEXsmXL5GBgUd1zF8dlWzGzy8xshpnNWL5ieYVdF5GuqDqZnOzuxxF2Ya4ws1MA3P0adx8K3AF8Ica2dyhwm0OI7j7F3ce6+9iBA5LHhESkQSpNJu6+JP5+A7gXOL5NyJ3AefH2YmBo3WNDgCVV9k9EyqksmZjZXma2T+02cBowx8xG14WdC7wQb08DLrLgBGCVuy+tqn8iUlaVp4YHAfeGM77sBtzp7g+b2T1mdjjh1PCrwOUx/kHCaeEFhFPDn62wbyJSWGXJxN1fBj7YzvLz2gnHw/yRV1TVHxGpVkuPgG2kZhuQtjljeGupAWka/CY5musbIiItS8lERIpQMhGRIpRMRKQIJRMRKULJRESKUDIRkSKUTESkCA1aa1G7NnC+RQ1IkxzaMhGRIpRMRKQIJRMRKULJRESKUDIRkSKUTESkiO4ownUj8HFgI/AS8Fl3f8vMRgDPA/Pi06e7++XbrFREmlIjxpmMd/cVdfcfAa52901mdgNwNfCV+NhLsTiXiLSYhg9ac/ef1t2dDpzf6D7IFqXKnubM/JarkQPycrRiadju0C1FuOpcAjxUd3+kmT1lZr80s3HtrVBFuESaU9VbJifHouQHAI+Y2Qvu/isAM7sG2EQoxAWwFBjm7ivNbAxwn5kd2aacKO4+BZgCMGbM2HL/DkVkh3RLES4zmwScA1wYZ6XH3Te4+8p4eybh4OxhVfZPRMrpjiJcZxAOuJ7r7uvq4gea2a7x9iHAaODlqvonImV1RxGuBUAvwm4PbDkFfApwnZltIpxKvtzd/1Bh/0SkoO4ownVoB/H3APdU1R8RqZbOZ4lIEUomIlKEZlprUUvefCcZ06dXeoa0fn12T8aUHLT10rK1yZgRA/skY0oNbNOAtHL0LolIEUomIlKEkomIFKFkIiJFKJmISBFKJiJShJKJiBShZCIiRWjQ2k6sT0ZZz5xBW+s2bk7G5A7sGjVor6y4EnJmf2u2Wd1ambZMRKQIJRMRKULJRESKUDIRkSIqTSZmttDMnjWz2WY2Iy670cxeMLNnzOxeM+tXF3+1mS0ws3lmdnqVfRORshqxZTLe3Y9x97Hx/iPAUe7+AWA+oQgXZnYEMBE4EjgD+OfanLAi0vwavpvj7j91903x7nRgSLw9AbgrzlL/CrCAOJu9iDS/ZirCdTCwqO6xxXHZVlSES6Q5NVMRrvZGD20z6khFuIIB++xRZD3zl65JxkyZsSgZc9Wpo7LaW/LW+mRM397p2d8OG7x3VnspK9dsTMYc0LdXkbZ2dk1ThIuwJTK07ulDgCVV9k9EymmaIlzANGCimfUys5GEIly/q6p/IlJW0xThcve5ZnY38Bxh9+cKd09fFCIiTaFpinDFx64Hrq+qTyJSHY2AFZEilExEpAglExEpQslERIrQTGuZSs3aVWpmsxw5ff7+7NeTMSeP7JuMGbRvuYFdOWVNF61MD37rt1d68Fv/vdOD/9Zn/D32zJjVbmenLRMRKULJRESKUDIRkSKUTESkCCUTESlCyUREilAyEZEilExEpAgNWsuUMyAtZ3BTThnNXhkxl/777GTMQ/enp4P50d+fl4wZd+iAZMwds15LxgDMWJSe2e2sw/snYz40bP9kzJzXVyVjThyVbqtkedSUVi5Xqi0TESmiO+rmfNLM5prZe2Y2ti52hJmtj7GzzezWKvsmImU1YjdnvLuvqLs/B/gz4DvtxL7k7sc0oE8iUljDj5m4+/MAccpGEdlJdHfdnLZGmtlTZvZLMxtXcd9EpKBuq5vTjqXAMHdfaWZjgPvM7Eh3f7s+KCalywCGDhtWaedFJF+HWyZmdrSZTTezRWY2xcz2q3ssqwRFR3VzOojd4O4r4+2ZwEvAYe3ETXH3se4+duCAgTndEJEG6Gw3ZzJwLXA0ocD4b8ysVrYtOetMR3VzOokfWCtUbmaHEOrmvJzxGkSkCXS2m7O3uz8cb3/TzGYCD5vZZ2inbGc7Oqqb8wng28BA4P+Y2Wx3Px04BbjOzDYBm4HL3f0P2/eyukfObFulZmy7c9KYZMzK844u0tZNv3opGTPh8AOTMQA/n5f+kz63fG0y5oaH5idjxh+d7tNRB++bjMkZkNbKg81K6SyZmJnt6+6rANz9UTM7D7gHSA4/7KRuzr2EXZ62y++J6xaRFtRZyr0B+E/1C9z9GeBjwI+r7JSItJ4Ot0zc/c4Olr8GfK6yHolIS9K1OSJShJKJiBTRaTIxs13N7K8b1RkRaV2dJhN33wxMaFBfRKSF5Qynf9zMbgGmAu8PAHD3WZX1SkRaTk4yOSn+vq5umQMfLd+dnd/KNRuTMU++tjIZ8/GjDkrG5JS+/JufPJ+Mee61N5MxFx07NBkD8LHD0zOk7dc7Xdbzq2f/UTJm317p9eQMIsyRU640x9D+exZZT3dIJhN3H9+IjohIa0uezTGzQWb2r2b2ULx/hJldWn3XRKSV5Jwa/gHwf4HadvV84L9X1SERaU05yWSAu98NvAfg7rUL8URE3peTTNaaWX/ilcJmdgKQriEgIj1KztmcLwLTgFFm9jhh6oBPVtorEWk5OclkLvDHwOGAAfPQMHwRaSMnKTzh7pvcfa67z3H3d4Enqu6YiLSWzuaAPTBO7LynmR1rZsfFn1OBPjkr70oRrvjY1Wa2wMzmmdnpO/C6RKTBOtvNOR24GBgC3FS3/G3gf3ShjawiXGZ2BDAROJJwGvpnZnZYvD6o2+XUEc6ZtrFfn/SozKMH9UvGPP1q+hj4/fOWJWMWLH07GbNoUbqt55el1wPwocHpEbAL3lydjLl95uvJmKtOPTQZM31herTxotXp0a0TP5geAbxP7527tHdnkyPdBtxmZufFKRWL6KQI1wTgLnffALxiZgsIs9lrl0qkBeQcM3l8B0bAdqUI18HAorr7i+MyEWkBOcnk+2z/CNiT3f044EzgCjM7pZPY9qb33uYqLDO7zMxmmNmM5SuWZ3ZDRKpW6QjYrhThImyJ1O94DgGWtLNOFeESaUKVjYDtahEuwsC4iWbWy8xGEopwZVUOFJHut70jYM/PeF6XinC5+1wzuxt4DtgEXNEsZ3JEJC1nPpNZZrbVCNg4cC31vC4V4YqPXQ9cn1q3iDSfZDKJ9X/PAkbE+NPMDHe/qdMnikiPkrOb8xPgHeBZ4kHYniin3mypOsL990lPt5jTn+OH9E3G7N0rPdDu5k+kaxbnTn94/wv/Lxnz3NI1yZj5C9NTSZ50/v9Mxjz54+uSMUcOStcjXr1+UzKmT8agxhzNWtc4J5kMcfcPVN4TEWlpOWdzHjKz0yrviYi0tJwtk+mEszK7AO8SDsK6u6e3oUWkx8hJJt8CTgSedfcydQFEZKeTs5vzIjBHiUREOpOzZbIUeCxe6LehtlCnhkWkXk4yeSX+7BF/RES2kTMC9m8b0RERaW05I2AHAlcSZkDrXVvu7i1Ra7hULdmNm9Lj9XJmWivVVo7dd0kfEvvch4cnY15dvi4Z8/e/eDGrTw8/+HQy5ht/eWoy5ojBeydjfnvgBcmYnAFgdzy9OBnzVx85JBmz+p30wLacz2tOzAF9eyVjSss5AHsH8AIwEvhbYCHw+wr7JCItKCeZ9Hf3fwXedfdfuvslwAkV90tEWkzOAdjaFcJLzexswoRFQ6rrkoi0opxk8ndmti/wN4R5SPoCf11pr0Sk5XSaTOL0A6Pd/QHC7GrjG9IrEWk5nR4ziTOdnbu9K++gCNf+ZvaImb0Yf+8Xl59qZqti7Gwz+/r2tisijZezm/NbM7sFmAqsrS1091mZbbQtwnUV8HN3/wczuyre/0p87Nfufk7mekWkieQkk5Pi7/pZZBzY3nEmE4BT4+3bgMfYkkxEpEXljIDdkeMktSJcDnzH3acAg9x9aVz3UjM7oC7+RDN7mnDG6EvuPncH2gbKzUqVMyAtp4Rojv57p69aWPLmO8mYcYcOSMbcOyddZvPxV9KlP7/zqW2m+23XPw7eJxmzYVN6UNbKdekBYHOeS5dHff2UEcmYG75yczLmc4+mL1XL+bu2sqzip/GUcNsRsOn57kIRriUxYTxiZi90EjsLGO7ua8zsLOA+QrmLtn25DLgMYOiwYTndF5EGSA5aM7NbgU8B/40wMdIngfQYbDoswrXMzAbHdQ8G3ogxb7v7mnj7QWB3M9vmX6uKcIk0p5wRsCe5+0XAm/GivxPZuvJeuzopwjUNmBTDJgH3x5gDLRbZMbPjY9/SJepFpCnk7Oasj7/XmdlBhC/4yIzndVSE6/fA3bH4+WuELR0Ihb0+b2abYpsTNSGTSOvISSYPmFk/4EbCcQ0Hvpt6UidFuFYCH2tn+S3ALRn9EZEmlHM25xvx5j1m9gDQ292TtYZFpGfJmc+kN/AXwEcIWyW/MbPJ7p4+NykiPUbObs7twGrCRX4AFwD/xpZjHSIiWcnkcHevP/bxaBxYJm00cqa1PhllPZ9/fXUyZsXaZA16Lj0uPePEgy8sTcYA/Ms9s5MxU6/c5pDaNobvv1cy5p6fpWd/e+yVdJnRW7+bHqD947npwX+TxqRHVOSUfW1WOT1/yszenwzJzD4MPF5dl0SkFeVsmXwYuMjMXov3hwHPm9mzhMp+qkMsIlnJ5IzKeyEiLS/n1PCrjeiIiLS21j3aIyJNRclERIpQMhGRIpRMRKQIa+ULc8eMGeuPPzmju7vRJTmzseUMfssZ2LZs1YZkzJqMkpVrN6Zjjjy4bzIG8gZlrVyzMRnzvRmvJWNWb0i/R/f/8uVkzNcmHpWMOWloela7QfumS3aWGrSW8/nIKVcKMGS/XjPdfWwqTlsmIlKEkomIFKFkIiJFVJpMuliEy8zsf5vZAjN7xsyOq7JvIlJWI7ZMxrv7MXUHcGpFuEYDP4/3Ac4kzEY/mjD7/OQG9E1ECumO3ZwJhOJbxN9/Wrf8dg+mA/1qs9iLSPOrOpnUinDNjPVuoE0RLqBWhOtgYFHdcxfHZSLSArKKcO2ArhThaq/03jaDYFSES6Q5VZpM6otwmdlWRbhiadD3i3ARtkTq6/EMIZQJbbvOKcAUCIPWqux/FXIGJW1+r8zLGtp/z2RMzsClPoUG0QGcM/mJZMxtnxmTjPniKaOSMdNfSZdduiRjFrmcErOr1qVnrMuZHS/n89ErIyanz6XLlVa2m9PVIlxx+UXxrM4JwKra7pCINL8qt0y6WoTrQeAsYAGwDvhshX0TkcIqSybbUYTLgSuq6o+IVEsjYEWkCCUTESlCyUREilAyEZEilExEpIiqR8BKGzmDiRq5npzBZjmDpOYvXZPV3o3npmcteyujZGnO4K5jh+6XjJnz+qpkzAF7907G5MyOlzMYMee93lDob1bqM1SjLRMRKULJRESKUDIRkSKUTESkCCUTESlCyUREilAyEZEilExEpIiWHrT2nqfLbeYMJtpZ5ZQizZltK2c9+/bZPatPOQO3+u2VXlepgXTD998rGZPT55xZ7V5atjYZc0DfxpUQLa05eyUiLafyZGJmu5rZU2b2QLz/UTObZWZzzOw2M9stLj/VzFbFgl2zzezrVfdNRMppxJbJXwHPA5jZLoRaORPd/SjgVbbMBwvw61iw6xh3v64BfRORQqouDzoEOBv4l7ioP7DB3efH+48A51XZBxFpjKq3TP4JuBKoXea4AtjdzGqlQs9n6/IWJ5rZ02b2kJkdWXHfRKSgKktdnAO84e4za8vipNETgX80s98Bq4Fa4ZZZwHB3/yDwbeC+DtZ7mZnNMLMZK1Ysr6r7ItJFVW6ZnAyca2YLgbuAj5rZD939CXcf5+7HA78CXgRw97fdfU28/SBhC2ZA25W6+xR3H+vuYwcMGFhh90WkKypLJu5+tbsPcfcRhK2RX7j7p2OpUMysF/AV4NZ4/0CLRXbM7PjYt3RJNhFpCt0xaO3LcRdoF2Cyu/8iLj8f+LyZbQLWE874dDpaaBfruYPScgZS5bw3OeVB9+md/pgctF96NjLI6/eK1RuLrGdwRp9KldFctHJ9kbZaWUOSibs/BjwWb38Z+HI7MbcAtzSiPyJSnkbAikgRSiYiUoSSiYgUoWQiIkUomYhIEUomIlKEkomIFNHSM601m5xSm6VmySo1ACpnQFqOt9alS3pC3kxiOfplzOy25M13kjEjBvZJxpSaaS1HqUGEOZ/F0rRlIiJFKJmISBFKJiJShJKJiBShZCIiRSiZiEgRSiYiUoSSiYgUUfmgNTPbFZgBvO7u55jZR4FvAnsAM4FL3X1TnLLxZuAsYB1wsbvPqrp/JTVr2cZGyB2MNmfR28mYUYPSJTtz3uucAWk5GjlD2sqMWeZyBq11x2exmYpwnQmMjj+XAZMb0DcRKaSZinBNAG73YDrQz8wGV9k/ESmnmYpwHQwsqnvu4rhMRFpAMxXham/HdJsrrOqLcC1XES6RplHlAdhaEa6zgN5A31iE69PAOAAzOw04LMYvZutSoUOAJW1X6u5TgCkAY8aMTV/OKSIN0TRFuIBpwEUWnACscvelVfVPRMpqpiJcDxJOCy8gnBr+bDf0TUS2UzMV4XLgikb0R0TK00xrUkTODGEARw3tW6S99Rs3J2NKlY7NmWktJ6aRA+1y+lN6MF7PHbIpIkUpmYhIEUomIlKEkomIFKFkIiJFKJmISBFKJiJShJKJiBRhYeBpazKz5YQJlrrDAMKUCt1F7ffs9hvZh+HuPjAV1NLJpDuZ2Qx3H5uOVPtqf+ftQz3t5ohIEUomIlKEksn2m6L21X43a4Y+vE/HTESkCG2ZiEgZ7t5jfoAzgHmE2dyuaufxXsDU+PiTwIi4vD/wKLAGuKWDdU8D5tTd359QyuPF+Hu/uvbfIpzSewY4rqL2bwReiG3cC/SL7b9EqBawBJgN3FpR+9cCr8c2ZhNm0au9/pXx9c8DTq+o/al1bS8EZsflk+Lr3xDfg1vbrKfLfSBM/DWvrr0DOlpXFe9BJ+1/EXgufgZ+TjjFW3vO5rr4aUW+X939BW/UD7Br/CIdQqgm+DRwRJuYv6h9uAjz1k6Nt/cCPgJc3t6HGfgz4M42H+b/RUxYwFXx/kvAJcDDsf2JwJMVtX8asFu8fUNd++OAuQ14/dcCX2rn/T8tfrifAf5zXLZr6fbbPP4t4OuxD6/GL16xzwDhyzy2nXbbruvuKt6DTtofD/SJtz9fW1e8v6b0d6wn7eYcDyxw95fdfSNwF6HwV70JhIqDAD8CPmZm5u5r3f03wDttV2pmexP+A/xdJ+u6DfgU4b/Nh+P9u4CRbF1srFj77v5Td69NfzYd+GBsfxGhhEjVr7+t42P7Y4B/jz9j47Ljq2o/lp39L7G94wlbKe+W/Ax0ou26TqOC96Aj7v6ou6+Ld6cTKj5Upiclk5wiX+/HxC/iKsLmZWe+QfjPt67N8kEeZ9ePv/vHddfaqLVf34+S7de7hLDLU3v9I4E/B75gZuPq4kq3/wUze8bMvkcoadIdr38csMzdX4zrXwqMNLOngIsIX+Z629MHgO+b2Wwz+1pMYO2t6x3gDcq/Bx21X+9S4KG6+71j/anpZvanGetP6knJJKfIV1YhsPeDzY4BDnX3e7ezH97md/H2zewaQqGz38RFS4FhwNeAXwN3mlltYtaS7U8GRgHHxDZrNaUb+vqBCwhbALX1rwOGufuxhF2j0+pef5f7EF3o7kcTEtc44DOdrKvt8h16DxLthxWafZqQNG+sWzzMw+jZ/wr8k5mNSrSR1JOSSU6Rr/djzGw3YF/gD52s80RgjJktJHxZDzOzx+Jjy2q7L/H3yrjuWhu19uv7UbJ9zGwScA5wYW3d7r7B3VfGdp8m7K9vUwhtR9t392Xuvtnd3wO+Cwzvhte/G+F4ytS69R8cXz+EJLu87vVvTx9w99fj79WEBFXbZWm7rt7AAYXfg87ax8z+BLgGONfdN9Q9Z0n8/TLhmMuxnbWRoyclk98Do81spJntQTi4Na1NzDS2/Ac9n1A4rMP/Cu4+2d0P8lBo7CPAfHc/tZ11TSIcfBtNOEI/Kba/kK2LjRVr38zOIBQ5OzfuN9de/xgz6x3bnxX79HIF7dcXnf8EMCO2NZOwtXBBvD0a+F3p9qM/AV5w98Xx/u+Bw81sVPwMXEQ4sPly3XO61Acz283MBsTbuxOS95wO1vVI6fegs/bN7FjgO4TPwBt1z9kvFsEjPvdkwlmfHVP6iG4z/xBOT84n/De+Ji67Lr7ZEP5z/AfhgNjvgEPqnruQ8B9iDeG/R9uzACPY+mxGf8LpuBfj7/3r2l9F2FJ5lvBfu4r2awdba6f/bo3tLwE2EnY9ZhF2Aapo/9/i63uG8AUZXPf6V8afecAdVbQfl/0AuLzNsusJp4Vrp4Y/viOfAUIymhlf51zgZracmdlmXaXfg0T7PwOW0eYUMHBS/Ns8HX9fWuL7pRGwIlJET9rNEZEKKZmISBFKJiJShJKJiBShZCIiRSiZSFMxs4vN7KAdXMcfmdkTZrbBzL5Uqm/SOSUT6bI4MrMqFwNdSibt9OcPwF8C3yzUJ8lQ5YdCmpSZjSBMg/AkYRj1fOAid19nZl8nDOTaE/gt8Ofu7nGY+m8JoyWnmdl84KuES/lXEq4PWWZm1xIuJBxMGKb+ReAE4EzC/CYfd/d3zWwMcBOwN2Fej4vjuscCd5jZesJw+SPaxrn70rb9IVzsB4CH0Z5vmNnZJd836Zy2THquw4Ep7v4B4G3CPBoQ5sr4kLsfRUgo59Q9p5+7/7G7f4twLcwJHi6Yuwu4si5uFHA24XL6HwKPergQbT1wdhz2/W3gfHcfA3wPuN7df0QYdn+hux9DuHZmm7gO+iPdTFsmPdcid3883v4hW3YLxpvZlUAfwiUAc4GfxLipdc8fAkyN1+DsAbxS99hDcevjWcKERA/H5c8Shr0fDhwFPBKvlt+VMLy/rVTc1HaeI91EyaTnansdhccLAP+ZMGvXorjL0rsuZm3d7W8DN7n7NDM7lTCzWs0GAHd/z8ze9S3XbLxH+MwZMNfdT0z0MRW3toPl0g20m9NzDTOz2pf0AsJuSy1xrIgzmJ3fyfP3JRwDgS1XueaaBwystW9mu5vZkfGx1cA+GXHSZJRMeq7ngUlm9gxhd2ayu79FuIr5WeA+wiX7HbkW+A8z+zVdrHfrYcrE84EbzOxpwhWtJ8WHfwDcamazCbs1HcV1yMwONLPFhIO/XzWzxW0mQJIK6KrhHiiezXkgHmQVKUJbJiJShLZMRKQIbZmISBFKJiJShJKJiBShZCIiRSiZiEgRSiYiUsT/ByIl7o+m2iujAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(1, 1)\n", "ax.set_xlabel('parameter 1')\n", "ax.set_ylabel('parameter 2')\n", "xbins = np.linspace(np.min(chain[:,0]), np.max(chain[:,0]), 25)\n", "ybins = np.linspace(np.min(chain[:,1]), np.max(chain[:,1]), 25)\n", "histplot = ax.hist2d(chain[:,0], chain[:,1], bins=(xbins, ybins), cmap=plt.cm.Blues)\n", "\n", "# Fix the aspect ratio\n", "ax.set_aspect((xbins[-1] - xbins[0]) / (ybins[-1] - ybins[0]))\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## A matrix of scatterplots\n", "\n", "We now have all the ingredients to create a matrix of parameter distribution plots." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/chonlei/anaconda3/envs/matplotlib_build/lib/python3.7/site-packages/matplotlib/axes/_axes.py:6521: MatplotlibDeprecationWarning: \n", "The 'normed' kwarg was deprecated in Matplotlib 2.1 and will be removed in 3.1. Use 'density' instead.\n", " alternative=\"'density'\", removal=\"3.1\")\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAt0AAALSCAYAAAABYW0nAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3XmcXXWZ4P/Pk4QQICCYFItZKMSIIkuQElHQDsJAFBRbZXFmBLfO4DK2dqvACK2tQ4s/e7EVBieIDfSAkBk6SsvSELtpWzRCAlHCEggYSJkYkpJVspDk+f1xT8FNUnVuJXVP3Vo+79frvu493+ecc59bVUme+uZ7nhOZiSRJkqTqjGp1ApIkSdJwZ9EtSZIkVcyiW5IkSaqYRbckSZJUMYtuSZIkqWIW3ZIkSVLFLLolSZKkill0S5IkSRWz6JYkSZIqNqbVCVRl4sSJ2d7e3uo0NAQs6VoCwEETDmpxJtVauHDhmsxsa3UekiSNRMO26G5vb2fBggWtTkNDwIwrZwBwx4fvaGkeVYuIx1udgyRJI5XLSyRJkqSKWXRLkiRJFbPoliRJkio2bNd09+TFF1+ks7OTdevWtTqVYWncuHFMnjyZnXbaqdWpSJIkDSojquju7Oxk9913p729nYhodTrDSmbS1dVFZ2cnBxxwQKvTkSRJGlRG1PKSdevWMWHCBAvuCkQEEyZM8H8RJEmSejCiim7AgrtCfm0lSZJ6NqKWl2jgtZ93U2l82cUnD1AmkiRJrTPiZrpbpauri+nTpzN9+nT23XdfJk2a9NL2hg0bWpbXvHnzeO9739uy95ckSRoJnOkeIBMmTGDRokUAfOUrX2H8+PF8/vOf32KfzCQzGTXK34UkSZKGk8qK7oj4PnAK8GRmHlKMXQ8cVOyyJ/B0Zk6PiHbgQWBJEZufmecUxxwJXAnsAtwM/GlmZn/z++ytn2XR7xb19zRbmL7vdL4181vbdczSpUt573vfy7HHHssvf/lLfvjDH3L44Yfz9NNPA3Ddddcxb948vve977Fq1So+8YlP8MQTTzBq1Ci+/e1vc/TRR29xvo6ODq655hoOOqj2ZT722GO59NJLWbt2LZ/73OdYt24du+66K1deeSXTpk3b4tgLLriAiRMn8tnPfhaA173udcybN4/Jkydz1VVXcemll7Jhwwbe+ta3cskll/jLgSRJUh9VWTVdCcysH8jMMzJzemZOB24A/qku/Gh3rLvgLlwGzAKmFY8tzjkcPPDAA3zsYx/j3nvvZdKkSb3u95nPfIYvfvGLLFiwgDlz5vDxj398m33OOOMM5syZA9RaJHZ1dXH44Yfz+te/np/97Gfce++9XHjhhVxwwQV9zm/x4sXMnTuXn//85yxatIiNGzdy3XXXbf8HlSRJGqEqm+nOzJ8WM9jbiFqbi9OBd5SdIyL2A/bIzF8U21cD7wVu6W9+2zsjXaUDDzyQN73pTQ33mzdvHkuWLHlp+6mnnmLt2rXssssuL42dfvrpvPvd7+bCCy/k+uuv5/TTTwfg6aef5qyzzuLRRx/d7vzmzZvH3XffTUdHBwBr165lypQp230eSZKkkapVa7rfBqzKzEfqxg6IiHuBZ4ELMvM/gElAZ90+ncVYjyJiFrVZcaZOndr0pKuy2267vfR61KhR1K+eqe97nZncddddjB07ttdz7b///owfP54HHniA66+/niuvvBKAL33pS5x00kl88pOfZOnSpcycue1/GIwZM4bNmzdv896ZyUc/+lG+9rWv7fBnlCRJGslatSj3g8AP6rZXAlMz8wjgz4BrI2IPoKfGz72u587M2ZnZkZkdbW1tTU14oIwaNYq99tqLRx55hM2bNzN37tyXYieccAKXXnrpS9vdF2Zu7YwzzuDrX/8669ev5+CDDwbgmWeeeWnpSnchvrX29nYWLlwIwF133cXy5ctfet85c+awZs0aoNaJ5YknnujfB5UkSRpBBrzojogxwPuA67vHMnN9ZnYVrxcCjwKvpTazPbnu8MnAioHLtjW+8Y1vMHPmTI4//ngmT37541966aXceeedHHbYYRx88MFcfvnlPR5/2mmnce211760tATg3HPP5Qtf+ALHHHNMr+972mmnsWrVKo444giuuOIKXv3qVwNw6KGH8uUvf5kTTjiBww47jBNPPJFVq1Y16dNKkiQNf9GERiC9n7y2pvvH3d1LirGZwPmZ+Ud1Y23A7zNzU0S8GvgP4NDM/H1E3A38d+CX1LqXfCczb2703h0dHblgwYItxh588EFe//rX9/+DqVdbf42Hws1xZlw5A4A7PnxHS/OoWkQszMyOVuchSdJIVNlMd0T8APgFcFBEdEbEx4rQmWy5tATg7cCvI+JXwP8DzsnM3xexTwDfA5ZSmwHv90WUkiRJ0kCqsnvJB3sZ/3APYzdQayHY0/4LgEN6ikmSJElDwYi7I2VmUutYqGarcqnS1obCshVJkqRuI+qWguPGjaOrq2tAi8ORIjPp6upi3LhxrU5FkiRp0BlRM92TJ0+ms7OT1atXtzqVYWncuHFbdFuRJElSzYgqunfaaScOOOCAVqchSZKkEWZELS+RJEmSWsGiW5IkSaqYRbckSZJUMYtuSZIkqWIW3ZIkSVLFLLolSZKkill0S5IkSRWz6JYkSZIqZtEtSZIkVcyiW5IkSaqYRbckSZJUMYtuSZIkqWKVFd0R8f2IeDIiFteNfSUifhsRi4rHu+pi50fE0ohYEhEn1Y3PLMaWRsR5VeUrSZIkVaXKme4rgZk9jP9dZk4vHjcDRMTBwJnAG4pj/ldEjI6I0cClwDuBg4EPFvtKkiRJQ8aYqk6cmT+NiPY+7n4qcF1mrgd+ExFLgaOK2NLMfAwgIq4r9n2gyelKkiRJlWnFmu5PR8Svi+UnexVjk4Dldft0FmO9jfcoImZFxIKIWLB69epm5y1JkiTtkIEuui8DDgSmAyuBvynGo4d9s2S8R5k5OzM7MrOjra2tv7lKkiRJTVHZ8pKeZOaq7tcRcTnw42KzE5hSt+tkYEXxurdxSZIkaUgY0JnuiNivbvOPge7OJjcCZ0bEzhFxADANuAu4G5gWEQdExFhqF1veOJA5S5IkSf1V2Ux3RPwAmAFMjIhO4MvAjIiYTm2JyDLgvwFk5v0RMYfaBZIbgU9l5qbiPJ8G/gUYDXw/M++vKmdJkiSpClV2L/lgD8NXlOx/EXBRD+M3Azc3MTVJkiRpQHlHSkmSJKliFt2SJElSxSy6JUmSpIpZdEuSJEkVs+iWJEmSKmbRLUmSJFXMoluSJEmqmEW3JEmSVDGLbkmSJKliFt2SJElSxSy6JUmSpIpZdEuSJEkVs+iWJEmSKmbRLUmSJFXMoluSJEmq2JhWJ6DBr/28m0rjyy4+eYAykSRJGpoqm+mOiO9HxJMRsbhu7JsR8VBE/Doi5kbEnsV4e0SsjYhFxeO7dcccGRH3RcTSiPh2RERVOUuSJElVqHJ5yZXAzK3GbgcOyczDgIeB8+tij2bm9OJxTt34ZcAsYFrx2PqckiRJ0qBW2fKSzPxpRLRvNXZb3eZ84ANl54iI/YA9MvMXxfbVwHuBW5qarAalRstaJEmShopWrun+KHB93fYBEXEv8CxwQWb+BzAJ6Kzbp7MY61FEzKI2K87UqVObnrCaz8JakiSNBC3pXhIRXwI2AtcUQyuBqZl5BPBnwLURsQfQ0/rt7O28mTk7Mzsys6Otra3ZaUuSJEk7ZIdmuiNifGY+v4PHng2cAhyfmQmQmeuB9cXrhRHxKPBaajPbk+sOnwys2JH3lSRJklplR2e6H9iRgyJiJnAu8J7MfKFuvC0iRhevX03tgsnHMnMl8FxEHF10LTkL+NEO5ixJkiS1RK8z3RHxZ72FgPGNThwRPwBmABMjohP4MrVuJTsDtxed/+YXnUreDnw1IjYCm4BzMvP3xak+Qa0Tyi7ULqD0IkpJkiQNKWXLS/4K+Ca1tddbazhDnpkf7GH4il72vQG4oZfYAuCQRu8nSZIkDVZlRfc9wA8zc+HWgYj4eHUpSZIkScNLWdH9EaCrl1hHBblIkiRJw1KvRXdmLimJraomHbWK/bIlSZKq05I+3ZIkSdJIYtEtSZIkVay06I6I0RHxuYFKRpIkSRqOSovuzNwEnDpAuUiSJEnDUl9uA39nRFwCXA/8oXswM++pLCtJkiRpGOlL0f3W4vmrdWMJvKP56UiSJEnDT8OiOzOPG4hEJEmSpOGqYfeSiNgnIq6IiFuK7YMj4mPVpyZJkiQND31pGXgl8C/Aq4rth4HPVpWQJEmSNNz0peiemJlzgM0AmbkR2FRpVpIkSdIw0pei+w8RMYHaxZNExNHAM5VmJUmSJA0jfele8mfAjcCBEXEn0AacVmlWkiRJ0jDSl6L7fuCPgIOAAJbg7eM1hLWfd9MW278b27XF+LKLTx7wnCRJ0vDWl+L5F5m5MTPvz8zFmfki8IuqE5MkSZKGi16L7ojYNyKOBHaJiCMi4o3FYwawa19OHhHfj4gnI2Jx3dgrI+L2iHikeN6rGI+I+HZELI2IX0fEG+uOObvY/5GIOHuHP60kSZLUAmXLS04CPgxMBv62bvxZ4H/08fxXApcAV9eNnQf8JDMvjojziu1zgXcC04rHm4HLgDdHxCuBLwMd1C7mXBgRN2bmU33MQZIkSWqpXovuzLwKuCoi3p+ZN+zIyTPzpxHRvtXwqcCM4vVVwB3Uiu5TgaszM4H5EbFnROxX7Ht7Zv4eICJuB2YCP9iRnCRJkqSB1pc13Xc2+Y6U+2TmSoDiee9ifBKwvG6/zmKst/FtRMSsiFgQEQtWr17djxQlSZKk5ulL0f0PDMwdKaOHsSwZ33Ywc3ZmdmRmR1tbW1OTkyRJknZUK+5IuapYNkLx/GQx3glMqdtvMrCiZFySJEkaElpxR8obge4OJGcDP6obP6voYnI08Eyx/ORfgBMjYq+i08mJxZgkSZI0JOzoHSk/0JeTR8QPqF0IOTEiOql1IbkYmFOsC3+Cl+9ueTPwLmAp8ALwEYDM/H1EfA24u9jvq90XVUqSJElDQcOiOzPviYgt7khZ3CCnocz8YC+h43vYN4FP9XKe7wPf78t7SpIkSYNNw6I7IkZTm4FuL/Y/MSLIzL8tPVBqoa1v9S5JktRKfVle8s/AOuA+iospJUmSJPVdX4ruyZl5WOWZSJIkScNUX4ruWyLixMy8rfJsNCS5lEOSJKlcX4ru+cDciBgFvEjtYsrMzD0qzUySJEkaJvpSdP8N8BbgvqLDiCRJkqTt0Jeb4zwCLLbgliRJknZMX2a6VwJ3RMQtwPruQVsGSpIkSX3Tl6L7N8VjbPGQJEmStB36ckfKvxyIRCRJkqThqi93pGwDvgi8ARjXPZ6Z76gwL0mSJGnY6MvykmuA64FTgHOAs4HVVSal5rOXtiRJUuv0pXvJhMy8AngxM/89Mz8KHF1xXpIkSdKw0ZeZ7heL55URcTKwAphcXUqSJEnS8NKXovt/RsQrgD8HvgPsAXyu0qwkSZKkYaS06I6I0cC0zPwx8Axw3IBkJUmSJA0jpWu6M3MT8J5mvmFEHBQRi+oez0bEZyPiKxHx27rxd9Udc35ELI2IJRFxUjPzkSRJkqrWl+UlP4+IS6h1MPlD92Bm3rMjb5iZS4Dp8NJM+m+BucBHgL/LzL+u3z8iDgbOpNay8FXAvIh4bfELgSRJkjTo9aXofmvx/NW6sQSa0af7eODRzHw8Inrb51TgusxcD/wmIpYCRwG/aML7S5IkSZXryx0pq1zHfSbwg7rtT0fEWcAC4M8z8ylgEjC/bp/OYmwbETELmAUwderUShKWJEmStldf+nQTESdHxBcj4i+6H/1944gYS229+P8thi4DDqS29GQl8Dfdu/ZwePZ0zsycnZkdmdnR1tbW3xQlSZKkpmhYdEfEd4EzgP9OrQA+Ddi/Ce/9TuCezFwFkJmrMnNTZm4GLqe2hARqM9tT6o6bTK1XuCRJkjQk9GWm+62ZeRbwVGb+JfAWtiyCd9QHqVtaEhH71cX+GFhcvL4RODMido6IA4BpwF1NeH9JkiRpQPTlQsq1xfMLEfEqoAs4oD9vGhG7Av8J+G91w/9fREyntnRkWXcsM++PiDnAA8BG4FN2LpEkSdJQ0pei+8cRsSfwTeAeakXx5f1508x8AZiw1diHSva/CLioP+8pSZIktUpfupd8rXh5Q0T8GBiXmc9Um5YkSZI0fDQsuiNiHPBJ4Fhqs9w/i4jLMnNd1clJkiRJw0FflpdcDTwHfKfY/iDwj9S6mEiSJElqoC9F90GZeXjd9r9FxK+qSkiSJEkabvrSMvDeiDi6eyMi3gzcWV1KkiRJ0vDSl5nuNwNnRcQTxfZU4MGIuA/IzDyssuwkSZKkYaAvRffMyrOQJEmShrG+tAx8fCASkSRJkoarvqzpliRJktQPFt2SJElSxSy6JUmSpIpZdEuSJEkVs+iWJEmSKmbRLUmSJFXMoluSJEmqmEW3JEmSVLGWFd0RsSwi7ouIRRGxoBh7ZUTcHhGPFM97FeMREd+OiKUR8euIeGOr8pYkSZK2V6tnuo/LzOmZ2VFsnwf8JDOnAT8ptgHeCUwrHrOAywY8U0mSJGkHNbwN/AA7FZhRvL4KuAM4txi/OjMTmB8Re0bEfpm5siVZDlLt593U6hQkSZLUg1bOdCdwW0QsjIhZxdg+3YV08bx3MT4JWF53bGcxtoWImBURCyJiwerVqytMXZIkSeq7Vs50H5OZKyJib+D2iHioZN/oYSy3GcicDcwG6Ojo2CYuSZIktULLiu7MXFE8PxkRc4GjgFXdy0YiYj/gyWL3TmBK3eGTgRUDmrBGjLJlOssuPnkAM5EkScNFS5aXRMRuEbF792vgRGAxcCNwdrHb2cCPitc3AmcVXUyOBp5xPbckSZKGilbNdO8DzI2I7hyuzcxbI+JuYE5EfAx4Ajit2P9m4F3AUuAF4CMDn7IkSZK0Y1pSdGfmY8DhPYx3Acf3MJ7ApwYgNUmSJKnpWt2nW5IkSRr2LLolSZKkill0S5IkSRWz6JYkSZIqZtEtSZIkVcyiW5IkSaqYRbckSZJUMYtuSZIkqWIW3ZIkSVLFLLolSZKkill0S5IkSRWz6JYkSZIqZtEtSZIkVcyiW5IkSaqYRbckSZJUMYtuSZIkqWIDXnRHxJSI+LeIeDAi7o+IPy3GvxIRv42IRcXjXXXHnB8RSyNiSUScNNA5S5IkSf0xpgXvuRH488y8JyJ2BxZGxO1F7O8y86/rd46Ig4EzgTcArwLmRcRrM3PTgGYtSZIk7aABn+nOzJWZeU/x+jngQWBSySGnAtdl5vrM/A2wFDiq+kwlSZKk5mjFTPdLIqIdOAL4JXAM8OmIOAtYQG02/ClqBfn8usM66aVIj4hZwCyAqVOnVpZ3q7Sfd1OrUxjxGn0Pll188gBlIkmShpKWXUgZEeOBG4DPZuazwGXAgcB0YCXwN9279nB49nTOzJydmR2Z2dHW1lZB1pIkSdL2a0nRHRE7USu4r8nMfwLIzFWZuSkzNwOX8/ISkk5gSt3hk4EVA5mvJEmS1B+t6F4SwBXAg5n5t3Xj+9Xt9sfA4uL1jcCZEbFzRBwATAPuGqh8JUmSpP5qxZruY4APAfdFxKJi7H8AH4yI6dSWjiwD/htAZt4fEXOAB6h1PvmUnUskSZI0lAx40Z2ZP6Pnddo3lxxzEXBRZUlJkiRJFfKOlJIkSVLFLLolSZKkill0S5IkSRWz6JYkSZIqZtEtSZIkVcyiW5IkSaqYRbckSZJUsVbcHEcattrPu6nX2LKLTx7ATCRJ0mDiTLckSZJUMYtuSZIkqWIuLxlkypYnSJIkaWhypluSJEmqmEW3JEmSVDGXl0gDpNHSIbubSJI0fDnTLUmSJFXMme4KeDGkJEmS6ll0S4NEf35Zc2mKJEmD25BZXhIRMyNiSUQsjYjzWp2PJEmS1FdDYqY7IkYDlwL/CegE7o6IGzPzgf6c11t2S5IkaSAMiaIbOApYmpmPAUTEdcCpQL+K7jKuy5YkSVKzDJWiexKwvG67E3jz1jtFxCxgVrH5fEQsGYDcGpkIrGl1Ek02LD/T45wyZD9TfKPH4a2/T/sPSDKSJGkbQ6Xojh7GcpuBzNnA7OrT6buIWJCZHa3Oo5n8TEPDcPxMkiQNVUPlQspOYErd9mRgRYtykSRJkrbLUCm67wamRcQBETEWOBO4scU5SZIkSX0yJJaXZObGiPg08C/AaOD7mXl/i9Pqq0G13KVJ/ExDw3D8TJIkDUmRuc3SaEmSJElNNFSWl0iSJElDlkW3JEmSVDGLbkmSJKliFt2SJElSxSy6JUmSpIpZdEuSJEkVs+iWJEmSKmbRLUmSJFXMoluSJEmqmEW3JEmSVDGLbkmSJKliFt2SJElSxSy6JUmSpIpZdEuSJEkVs+iWJEmSKmbRLUmSJFXMoluSJEmqmEW3JEmSVDGLbkmSJKliFt2SJElSxSy6JUmSpIpZdEuSJEkVs+iWJEmSKmbRLUmSJFXMoluSJEmqmEW3JEmSVDGLbkmSJKliFt2SJElSxSy6JUmSpIpZdEuSJEkVs+iWJEmSKmbRLUmSJFVsTKsTiIhlwHPAJmBjZnZExNeAU4HNwJPAhzNzRUQE8PfAu4AXivF7ys4/ceLE3H//9sry39wg3urfajZneXxUDEweQ8XDXUsAeO2Eg1qcSXPcc8/CNZnZ1uo8Gpk4cWK2t7e3Oo0BtaT4WTtomPysqbkWLhwaf3Yl9V3Li+7CcZm5pm77m5l5IUBEfAb4C+Ac4J3AtOLxZuCy4rlX++/fzp2/XFBJ0gAbNpaX3WPHtLbsXrthU2l8l7GjByiToeHEf5wBwG0fuqOleTTLLjvF463OoS/a29tZsKC6P6eD0YwrZwBwx4fvaGkeGpwihsafXUl91+qJ2B5l5rN1m7sB3fO1pwJXZ818YM+I2G/AE5QkSZK2w2AouhO4LSIWRsSs7sGIuCgilgP/hdpMN8AkYHndsZ3FmCRJkjRoDYai+5jMfCO1pSOfioi3A2TmlzJzCnAN8Oli355WIG+zajkiZkXEgohYsHrN6qryliRJkvqk5Wu6M3NF8fxkRMwFjgJ+WrfLtcBNwJepzWxPqYtNBlb0cM7ZwGyAI4/saHApoSSpNy+++CKdnZ2sW7eu1akMS+PGjWPy5MnstNNOrU5FUsVaWnRHxG7AqMx8rnh9IvDViJiWmY8Uu70HeKh4fSPw6Yi4jtoFlM9k5soBT1ySRojOzk5233132tvbqTWQUrNkJl1dXXR2dnLAAQe0Oh1JFWv1TPc+wNziL/IxwLWZeWtE3BARB1HryPc4tc4lADdTaxe4lFrLwI8MfMqSNHKsW7fOgrsiEcGECRNYvdplkNJI0NKiOzMfAw7vYfz9veyfwKeqzkuS9DIL7ur4tZVGjlbPdA95re7DvanB3W/624fbPt/S0NJ+3k2l8WUXnzxAmUiS6ll0S9IQ06iwHk66uro4/vjjAfjd737H6NGjaWur3ajxrrvuYuzYsS3Ja968eVxyySX88Ic/bMn7Sxp6LLolSUDPxfzl79mPFzufBuCwyXsOdEpMmDCBRYsWAfCVr3yF8ePH8/nPf36LfTKTzGTUqMHQBVeSeubfUJKkIWfp0qUccsghnHPOObzxjW9k+fLl7Lnny78UXHfddXz84x8HYNWqVbzvfe+jo6ODo446ivnz529zvo6ODpYsWfLS9rHHHsuvfvUr5s+fz1ve8haOOOIIjjnmGB555JFtjr3gggv41re+9dL26173Ojo7OwG46qqrOOqoo5g+fTqf/OQn2bx5c9O+BpKGFme6JUl98tlbP8ui3y1q6jmn7zudb838VuMde/DAAw/wD//wD3z3u99l48aNve73mc98hi9+8YscffTRLFu2jFNOOYXFixdvsc8ZZ5zBnDlzuPDCC+ns7KSrq4vDDz+cZ555hp/97GeMHj2aW2+9lQsuuIDrr7++T/ktXryYuXPn8vOf/5wxY8Ywa9YsrrvuOv7zf/7PO/R5JQ1tFt2SpCHpwAMP5E1velPD/ebNm7fFLPZTTz3F2rVr2WWXXV4aO/3003n3u9/NhRdeyPXXX8/pp58OwNNPP81ZZ53Fo48+ut35zZs3j7vvvpuOjg4A1q5dy5QpUxocJWm4suiWpEFoMF4suaMz0lXZbbfdXno9atQoal1la+rvoJmZDS+63H///Rk/fjwPPPAA119/PVdeeSUAX/rSlzjppJP45Cc/ydKlS5k5c+Y2x44ZM2aLZSPd752ZfPSjH+VrX/vaDn9GScOHa7olSUPeqFGj2GuvvXjkkUfYvHkzc+fOfSl2wgkncOmll7603X1h5tbOOOMMvv71r7N+/XoOPvhgAJ555hkmTZoE8FIhvrX29nYWLlwI1DqqLF++/KX3nTNnDmvWrAFqnVieeOKJ/n1QSUOWRfcQN3pUlD76a5exo0sfkjRYfOMb32DmzJkcf/zxTJ48+aXxSy+9lDvvvJPDDjuMgw8+mMsvv7zH40877TSuvfbal5aWAJx77rl84Qtf4Jhjjun1fU877TRWrVrFEUccwRVXXMGrX/1qAA499FC+/OUvc8IJJ3DYYYdx4oknsmrVqiZ9WklDTdT/d9xwdOSRHXnnLxe0Og0NESf+4wwAbvvQHS3No1l22SkWZmZHq/NopKOjIxcsGFl/TmdcOQOAOz58R4/xqpaXlN0cp7eWgftMrRWRrWgZOBI8+OCDvP71r99iLGJo/NmV1HfOdEuSJEkV80JKSRpBBuMFmpI0EjjTLUnqVZIM92WIreTXVho5Wl50R8SyiLgvIhZFxIJi7JsR8VBE/Doi5kbEnsV4e0SsLfZdFBHfbW32kjS8Pf70i2x84VmLwwpkJl1dXYwbN67VqUgaAINleclxmbmmbvt24PzM3BgR3wDOB84tYo9m5vQBz1CSRqDv/PIp/juw/55rePD5XVudzrAzbty4LTqtSBq+BkvRvYXMvK1ucz7wgVblIqkaETELmAUwderUFmej3jy7fjMX/bQLKO98Ikkq1/LlJUACt0XEwuIf4a19FLilbvuAiLg3Iv7yedgVAAAgAElEQVQ9It42MCkOXxs2bi59NLJpc/br0Wpbf95MyKTPn187LjNnZ2ZHZna0tbW1Oh1Jkio1GGa6j8nMFRGxN3B7RDyUmT8FiIgvARuBa4p9VwJTM7MrIo4EfhgRb8jMZ+tPWD+DNsUZNEmSJLVYy2e6M3NF8fwkMBc4CiAizgZOAf5LFlfwZOb6zOwqXi8EHgVe28M5X55Bm+gMmiRJklqrpUV3ROwWEbt3vwZOBBZHxExqF06+JzNfqNu/LSJGF69fDUwDHhv4zCVJkqS+a/Xykn2AuRHRncu1mXlrRCwFdqa23ARgfmaeA7wd+GpEbAQ2Aedk5u9bk7okSZLUNy0tujPzMeDwHsZf08v+NwA3VJ2XJEmS1EwtX9MtSZIkDXcW3ZIkSVLFWr2mW/204ql1pfFddx5dGt9z151K4416VY8dU/5726Or/lAab28rv8Pd6FFRGm9ke/OP6HlckiSpP6wsJEmSpIpZdEuSJEkVs+iWJEmSKmbRLUmSJFXMoluSJEmqmEW3JEmSVDGLbkmSJKli9uke5nYdW96nu1Ef6xc2bCqNN+pnfeA+u5XG+2vT5iyN97fPtyRJUjM40y1JkiRVzKJbkiRJqljLi+6IWBYR90XEoohYUIx9MyIeiohfR8TciNizbv/zI2JpRCyJiJNal7kkSZLUNy0vugvHZeb0zOwotm8HDsnMw4CHgfMBIuJg4EzgDcBM4H9FRPmiZUmSJKnFBkvRvYXMvC0zNxab84HJxetTgesyc31m/gZYChzVihwlSZKkvhoMRXcCt0XEwoiY1UP8o8AtxetJwPK6WGcxJkmSJA1ag6Fl4DGZuSIi9gZuj4iHMvOnABHxJWAjcE2xb0/937bpGVcU77MApkydWk3WkiRJUh+1vOjOzBXF85MRMZfacpGfRsTZwCnA8ZnZXVh3AlPqDp8MrOjhnLOB2QBHHtlR3sh5iJu4+9h+Hf/wyudL47MXLC+NnzfjwNL4iqfXlsb3GLdTafy1+40vjTfS9fyG0vjee+zcr/NLkiT1RUuXl0TEbhGxe/dr4ERgcUTMBM4F3pOZL9QdciNwZkTsHBEHANOAuwY6b0mSJGl7tHqmex9gbkR053JtZt4aEUuBnaktNwGYn5nnZOb9ETEHeIDaspNPZWb5LRMlSZKkFmtp0Z2ZjwGH9zD+mpJjLgIuqjIvSVJztZ93U6+xZRefPICZSFJrDIbuJZIkSdKwZtEtSZIkVcyiW5IkSaqYRbckSZJUsVZ3LxnyNm0ubwM+elRP9/N52YaNm0vjL2zoX3OWRvldtei3pfFjDtijNL7PK/rX53rXnUeXxpd3lff53nO38j7fE8aX9zFfu9XXd3PREr57fJex5flpx9XfxGqqN7GSJA1zznRLaonMnJ2ZHZnZ0dbW1up0JEmqlEW3JEmSVDGXl0iS+sRe25K045zpliRJkipm0S1JkiRVzKJbkiRJqphFtyRJklQxL6Tsp0Z9uLfuA721sWPKf+/ZuUH8Yz9YVBq/5Ud3lcb/31+9vzT+ttdMLI1fc88TpfEFy58vjb/roAml8TdNfWVpfPFvnymNv+XA8vNv3Qd90+Ytxxt9fxpp9PMhSZJGhpbPdEfEsoi4LyIWRcSCYuy0iLg/IjZHREfdvu0RsbbYd1FEfLd1mUuSJEl9M1hmuo/LzDV124uB9wH/u4d9H83M6QOTliRJktR//Sq6I+JQ4HJgEnALcG5mPlXE7srMo3bkvJn5YHGO/qQnSYPW/Me6gPLe15Kk4aO/y0suA74CHAo8DPwsIg4sYjv18RwJ3BYRCyNiVh/2PyAi7o2If4+It213xpIkSdIA6+/ykvGZeWvx+q8jYiFwa0R8iFox3RfHZOaKiNgbuD0iHsrMn/ay70pgamZ2RcSRwA8j4g2Z+Wz9TkXxPgtgytSp2/2hJEmSpGbq70x3RMQrujcy89+A9wP/COzflxNk5ori+UlgLtDrkpTMXJ+ZXcXrhcCjwGt72G92ZnZkZkfbxLbt+DiSJElS8/V3pvsbwOuB+d0DmfnriDgeuLDRwRGxGzAqM58rXp8IfLVk/zbg95m5KSJeDUwDHuvnZ5Ak9ZNr0yWpXL+K7sy8tpfxJ4A/6cMp9gHmFhdMjgGuzcxbI+KPge8AbcBNEbEoM08C3g58NSI2ApuAczLz9/35DFXbZezo0vimzeWrcBr1eb727CNL413vP7Rf5//bnz5aGj/1oH1L4z9ZUv7teWD1H0rj37jl4dL4cYeWv/8hk15RGt+6D3f3l6N73D7bkiSpGVraMjAzHwMO72F8LrWlJluP3wDcMACpSZIkSU3T8pvjSJIkScNdv4vuiBgdEZ9rRjKSJEnScNTvojszNwGnNiEXSZIkaVhq1pruOyPiEuB64KUr4zLzniadX5IkSRqymlV0v7V4rm/3l8A7mnR+SZIkachqStGdmcc14zySJEnScNSUojsi9gH+CnhVZr4zIg4G3pKZVzTj/MNZ1/MbSuO/fKKrNP7uQ15VGp8wfmxp/M//+cHS+ANPPFUaP+uIKaXx4w96ZWl8r3E7lcYvOPl1pfFX7Fx+fKM+6Fvr3rv7uOVda7fr+K1NmbBLv46XJEnDQ7NaBl4J/AvQXQE+DHy2SeeWJEmShrRmFd0TM3MOsBkgM7vvGClJkiSNeM0quv8QERMo/nc+Io4GnmnSuSVJkqQhrVndS/4MuBE4MCLuBNqA05p0bkmSJGlIa1bRfT/wR8BBQABL8BbzkiRJEtC8wvgXmbkxM+/PzMWZ+SLwiyadW5IkSRrS+lV0R8S+EXEksEtEHBERbyweM4Bd+3iOZRFxX0QsiogFxdhpEXF/RGyOiI6t9j8/IpZGxJKIOKk/+UuSJEkDob/LS04CPgxMBv62bvxZ4H9sx3mOy8w1dduLgfcB/7t+p6L/95nAG6i1J5wXEa/NzJZ1Slm7ofytdxk7ujS+567lfaYP3WfP0vivHi+/XvVHS1aVxpeufLY0vnx5+fkfXFV+/Jv2K+/TvfSp50rjVy/8bWn8vBmvKY3PX1be53z5c1v24V79/HoA/u99nQCceXh5H/LdxzVrhZYkSRrO+lUxZOZVwFUR8f7MvKFJOZGZDwJExNahU4HrMnM98JuIWAochUtZJEmSNIg1a033nRFxRUTcArUZ6Yj4WB+PTeC2iFgYEbMa7DsJWF633VmMSRpiImJWRCyIiAWrV69udTqSJFWqWUX3P7Djd6Q8JjPfCLwT+FREvL1k322mvnn5zt0v71T/j/ka/zGXBqPMnJ2ZHZnZ0dbW1up0JEmqVMvvSJmZK4rnJ4G51JaL9KYTqF9kOxlY0cM5X/7HfKL/mEuSJKm1WnpHyojYLSJ2734NnEjtIsre3AicGRE7R8QBwDTgrv4mL0mSJFWpyjtSfqAPx+0DzC0umBwDXJuZt0bEHwPfKc5zU0QsysyTMvP+iJgDPABsBD7Vys4lkiRJUl80pejOzHsiYos7UhY3yGl03GPA4T2Mz6W21KSnYy4CLupfxs0zdkz5fxZs2rzNkvMtjB7V0zL1l03YfWy/3v+oyXuUxsfvXN7S8O//+NDSeKPP96OHflcaf2Dl86Xxh5c9VRp/6we+XBr/5T99tTT+hn1escX2/7m/9vWe+Zp9AXhu7cbS43dt0BKyTKPvvSRJGj6aUnRHxGjgXUB7cc4TI4LM/NvSAyVJkqQRoFnLS/4ZWAfcR3ExpSRJkqSaZhXdkzPzsCadS5IkSRpWmtW95JaIOLFJ55IkSZKGlWbNdM+n1oVkFPAitYspMzPLr+KTJEmSRoBmFd1/A7wFuC8zy9tZSJIkSSNMs5aXPAIstuCWJEmSttWsme6VwB0RcQuwvntwKLQMbNRnupENG8ubtezSjz7OfTl/IzuNKv+96k/evH9p/PHVL5TG/+pfHymN33rzr0rjX/vMjNL4wfuNL43/fN8PlsYb9cK+5ledW2yv+sP6Lcb/9NhXlx7/3LryPt5lP1+Nfvb23mPn0rgkSRo6mlV0/6Z4jC0ekiRJkgrNuiPlXzbjPJIkSdJw1Kw7UrYBXwTeAIzrHs/MdzTj/JIkSdJQ1qwLKa8BHgIOAP4SWAbc3aRzS5IkSUNas4ruCZl5BfBiZv57Zn4UOLpJ55YkSZKGtGYV3S8Wzysj4uSIOAKY3JcDI2JZRNwXEYsiYkEx9sqIuD0iHime9yrGZ0TEM8W+iyLiL5qUvyRJklSZZnUv+Z8R8Qrgz4HvAHsAn9uO44/LzDV12+cBP8nMiyPivGL73CL2H5l5SjOSliRJkgZCv4vuiBgNTMvMHwPPAMf1Oys4FZhRvL4KuIOXi+6matTHuZFGfbjXbtjUr/NPGF/egXHFU+tK4297zcTS+NzFvy2N3/mbZ0vj//uMw0vjf7ff7qXx9RvLe1V3vVDeB3vxA6tK4799e3tp/Bvn/v2W+Rxb68/9jatr43/yb+Wt5ht9fyRJkqAJy0sycxPwnv6cArgtIhZGxKxibJ/MXFmcfyWwd93+b4mIX0XELRHxhn68ryRJkjQgmrW85OcRcQlwPfCH7sHMvKcPxx6TmSsiYm/g9oh4qGTfe4D9M/P5iHgX8ENg2tY7FcX7LIApU6dux8eQpOZpP++m3oP+J4kkjSjNKrrfWjx/tW4sgYZ9ujNzRfH8ZETMBY4CVkXEfpm5MiL2A54s9nm27ribI+J/RcTErdaDk5mzgdkARx7Z0b/7vEuSJEn91Kw7Uu7QOu6I2A0YlZnPFa9PpFa43wicDVxcPP+o2H9fYFVmZkQcRW15TFcTPoIkSZJUmWbNdBMRJ7PtHSm/2vsRAOwDzI2I7lyuzcxbI+JuYE5EfAx4Ajit2P8DwCciYiOwFjgzM53JliRJ0qDWrNvAfxfYlVrnku9RK47vanRcZj4GbNP+IjO7gON7GL8EuKS/+UqSJEkDqVk3x3lrZp4FPJWZfwm8BZjSpHNLkiRJQ1qzlpesLZ5fiIhXUVtnfUCTzj2kNerj3ciGjZtL47vuXH7+B3/7XGl8zR9eLI1/7I3lNxa9+aGVpfHv3bCoNH79F7f5D40t7P/K3UrjN8x7pDR+x2+eKo1/9/It279//RcLATj/7Nr4P91f3sf87CP3L42PHdOs32slSdJQ1qyi+8cRsSfwTWpt/RK4vEnnliRJkoa0ZnUv+Vrx8oaI+DEwLjOfaca5JUmSpKGuWRdSjgM+CRxLbZb7ZxFxWWaW36Nckoaw0pvfSJJUp1nLS64GngO+U2x/EPhHXm71J0mSJI1YzSq6D8rM+tZ//xYRv2rSuSUNQxExC5gFMHXq1BZnI0lStZrVWuHeiDi6eyMi3gzc2aRzSxqGMnN2ZnZkZkdbW1ur05EkqVLNmul+M3BWRDxRbE8FHoyI+4DMzMOa9D6SJEnSkNOsontmk84z4qzdsKk03t8+368cP7Y0/o728hnGP2zYWBp/z8GvKo1/4LLyZf1dz28ojX9/wROl8fefMK00PudfHy2NX3jmIaXxma/ZtzRepUY92p9bV/69kYaKRhekLrv45AHKRJKq06yWgY834zySJEnScOTt8iRJkqSKWXRLkiRJFWt50R0RyyLivohYFBELirFXRsTtEfFI8bxXMR4R8e2IWBoRv46IN7Y2e0mSJKmxlhfdheMyc3pmdhTb5wE/ycxpwE+KbYB3AtOKxyzgsgHPVJIkSdpOzepe0mynAjOK11cBdwDnFuNXZ2YC8yNiz4jYLzNX9naih7uWcOI/zugt3HKbM0vjoyJK4w0Ob9gBY9Pm8hNsavAGu+1c3l0lGuS/cVN5fiufW18ab5T/ml3Wlca//ovdtth+4tn7i/HTAXjFPTuVHj92TPnvrQ0+fqlG39uNDT67NFyUdTexs4mkoWIwzHQncFtELCzuUAewT3chXTzvXYxPApbXHdtZjG0hImZFxIKIWPDiiy9WmLokSZLU2GCY6T4mM1dExN7A7RHxUMm+Pc0bbjPdl5mzgdkARx7Zkbd96I6mJFqFRjO1/T2+0Uxso17PuzboE95oJv19l/+yNH7Nh44sje+5a/lM8/zfdJXGJ+2xa2l89Kgtf6TO+uday/mr330rAM83+Prst9e40njZ13/nBt+brXPb3jjALn/Sj6l2SZLUNC2f6c7MFcXzk8Bc4ChgVUTsB1A8P1ns3glMqTt8MrBi4LKVJEmStl9Li+6I2C0idu9+DZwILAZuBM4udjsb+FHx+kZqt5uPiDgaeKZsPbckSZI0GLR6eck+wNziYrsxwLWZeWtE3A3MiYiPAU8A3fcSvxl4F7AUeAH4yMCnLEmSJG2flhbdmfkYcHgP413A8T2MJ/CpAUhNkiRJappWz3RL0qBW1q5OkqS+avmFlJIkSdJwZ9EtSZIkVczlJS3Wl17LVR7fqM92o17SD698vjT+zfccUhp/+g/lNy9q1Gf8iCl7lcYX//aZ0vje47fss53FbSC7+5/v0qBPeaM+6WVfv/X9/Nr393svSZIGjjPdkiRJUsUsuiVJkqSKWXRLkiRJFbPoliRJkipm0S1JkiRVzKJbkiRJqphFtyRJklSxYd+ne3PC2g2beo036sM81JV9doAJ48f26/hX7LpTabxRH+s9dys/vr99wvd/5W6l8a3zGxW13tfd/cGnTNil9PhHV/2hNL73Hjv3GmvUg1ySJA0f/qsvSZIkVWxQFN0RMToi7o2IHxfb74iIeyJicURcFRFjivEZEfFMRCwqHn/R2swlSZKkxgZF0Q38KfAgQESMAq4CzszMQ4DHgbPr9v2PzJxePL468KlKkiRJ26fla7ojYjJwMnAR8GfABGB9Zj5c7HI7cD5wRWsylDSctZ93U6tTkCSNAINhpvtbwBeBzcX2GmCniOgotj8ATKnb/y0R8auIuCUi3tDTCSNiVkQsiIgFa9asrixxSZIkqS9aWnRHxCnAk5m5sHssMxM4E/i7iLgLeA7YWITvAfbPzMOB7wA/7Om8mTk7Mzsys2PixLZKP4OkHVP/y/Hq1f5yLEka3lo9030M8J6IWAZcB7wjIv5PZv4iM9+WmUcBPwUeAcjMZzPz+eL1zdRmxCe2KHdJ/VD/y3Fbm78cS5KGt5au6c7M86mt1yYiZgCfz8z/GhF7Z+aTEbEzcC619d5ExL7AqszMiDiK2i8NXWXvMSqGdy/uRn2wG33259ZtLI3vPq78R+RVe40rjTfKb81zG/p1/H4N3n/0qCiNb92nfOedtuzPvbxrbb/OL0mSBIPgQspefKFYejIKuCwz/7UY/wDwiYjYCKyl1uGkvCqTJEmSWmzQFN2ZeQdwR/H6C8AXetjnEuCSAU1MkiRJ6qdWr+mWJEmShj2LbkmSJKliFt2SJElSxSy6JUmSpIpZdEuSJEkVGzTdS4arDRs3l8bHjunf7z397RPdqA93I0+/8GJpfO89du7X+ffcdafS+Iqn1pXG29t2LY1v3Qc8txrv7te9o8r6oDf62jf62VHftZ93U6tTkCSNcM50S5IkSRWz6JYkSZIqZtEtSZIkVcyiW5IkSaqYRbckSZJUMbuXSJKGrbLONcsuPnkAM5E00ll0SxrybAkoSRrsBkXRHRGjgQXAbzPzlIh4B/DXwFhgIfCxzNwYEQH8PfAu4AXgw5l5T6vy7ov+9uEe7Br14V68/NnS+IH77FYab/T1a9SHu5Gt+5xHL+M7quu5Db3GGvXpHu4/O5IkjSSD5V/1PwUeBIiIUcBVwJmZeQjwOHB2sd87gWnFYxZw2cCnKkmSJG2flhfdETEZOBn4XjE0AVifmQ8X27cD7y9enwpcnTXzgT0jYr8BTViSJEnaTi0vuoFvAV8Euu95vQbYKSI6iu0PAFOK15OA5XXHdhZjW4iIWRGxICIWrF6zupqsJUmSpD5q6ZruiDgFeDIzF0bEDIDMzIg4E/i7iNgZuA3Y2H1ID6fJbQYyZwOzAY48smObuCRpePAiWklDRasvpDwGeE9EvAsYB+wREf8nM/8r8DaAiDgReG2xfycvz3oDTAZWDGC+kiRJ0nZr6fKSzDw/MydnZjtwJvCvmflfI2JvgGKm+1zgu8UhNwJnRc3RwDOZubIVuUuSJEl91eqZ7t58oVh6Mgq4LDP/tRi/mVq7wKXUWgZ+pEX5SZIkSX02aIruzLwDuKN4/QXgCz3sk8CnBjQxlXpu3cbS+CFT9ujX+ddu2FQa32Xs6H6df9PmLZf851bjW8e3VmUf8Ubv3axe4pIkqXqDpuiWpDJeMCdJGsoGQ8tASZIkaViz6JYkSZIq5vISSdKI1GjJ0rKLTx6gTCSNBBbdkiRtp/5cY2AxL41MLi+RJEmSKmbRLUmSJFUsaq2vh6+IWA08PkBvNxFYM0Dv1Sp+xqFl/8xsa3USPYmIWcCsYvMgYEkTTz9UvodDIU9zbI7tzXHQ/tmVtGP+//buPE6Ostr/+Pdk3xeSQAghCWvCIgQZwxJ2EMGwiUAEBJElelEUuawKiIIa7kUFAeEXdpEliIIIl0VUVMIaIEBYkgAJSYCsJCFkX87vj6cmzgwzTyXTU1PVPZ/365VXpvtUd5+q7uo+/fRTpyq+6G5OZjbe3avyziNLrCPKQbk8h+WQJzk2jXLIEUC2mF4CAAAAZIyiGwAAAMgYRXfTGpN3As2AdUQ5KJfnsBzyJMemUQ45AsgQc7oBAACAjDHSDQAAAGSMohsAAADIGEU3AAAAkDGKbgAAACBjFN0AAABAxii6AQAAgIxRdAMAAAAZo+gGAAAAMkbRDQAAAGSMohsAAADIGEU3AAAAkDGKbgAAACBjFN0AAABAxii6AQAAgIxRdAMAAAAZo+gGAAAAMkbRDQAAAGSMohsAAADIGEU3AAAAkDGKbgAAACBjFN0AAABAxii6AQAAgIxRdAMAAAAZo+gGAAAAMkbRDQAAAGSMohsAAADIGEU3AAAAkDGKbgAAACBjFN0AAABAxii6AQAAgIxRdAMAAAAZo+gGAAAAMtYm7wSy1rt3bx84cFDeabRYnhK3Zsli/U2eP0mStG2vwTln0jRefvmlee7eJ+881lel76+lvr7KbX/ChimX/bV3794+aNCgvNNAHZOS95fBFfL5VU5eemn99t2KL7oHDhykcc+PzzuNFmvN2niZ0LpVscqEg+/cT5L0xElP5ZpHU+nY1t7PO4cNUen7a6mvr3Lbn7BhymV/HTRokMaPr9z9tFztd/t+kqSnTnkq1zxaIrP123eZXgIAAABkjKIbAAAAyBhFNwAAAJCxws7pNrMekm6WtKPC8UOnSpokaaykQZKmSTrO3RfklCIAAEBhdGndRVOnTtXy5cvzTqUidejQQf3791fbtm0bdfvCFt2SrpH0mLsfY2btJHWS9ENJf3P30WZ2oaQLJV2QZ5IAAABFcPzA49W1a1cNGjRIZhxY3ZTcXfPnz9fMmTO1xRZbNOo+Cjm9xMy6SdpH0i2S5O4r3X2hpCMl3ZEsdoeko/LJEAAAoFg27bipevXqRcGdATNTr169SvoVoZBFt6QtJc2VdJuZvWJmN5tZZ0mbuPtHkpT8v3F9NzazUWY23szGz503t/myBrDB2F8BoGmYjII7Q6Vu26JOL2kj6fOSznL3583sGoWpJOvF3cdIGiNJu+5alXY+CWQorW/wvMUrovHeXds3ZTooIPbX9Zd1H272R6D5DbrwkQZj00aPaMZMkLWijnTPlDTT3Z9PLt+vUITPNrNNJSn5f05O+QEAAKCG+fPna+jQoRo6dKj69u2rzTbbbN3llStX5pbXk08+qaOOyn9GciFHut19lpnNMLPB7j5J0oGS3kz+fUPS6OT/P+eYJgAAABK9evXShAkTJEmXXXaZunTponPPPbfWMu4ud1erVkUd981Okdf4LEl3mdlrkoZK+rlCsf1FM5si6YvJZQAAABTUO++8ox133FHf/va39fnPf14zZsxQjx491sXvvfdenX766ZKk2bNn6+ijj1ZVVZWGDRum55577jP3V1VVpUmTJq27vNdee+nVV1/Vc889pz322EO77LKLhg8frilTpnzmthdffLGuvvrqdZeHDBmimTNnSpLuuOMODRs2TEOHDtWZZ56ptWvXNtk2kAo60i1J7j5BUlU9oQObOxcAAIBycvZjZ2vCrAlNep9D+w7V1Ydcnb5gPd58803ddtttuvHGG7V69eoGl/ve976n888/X7vvvrumTZumww47TBMnTqy1zMiRI3Xffffpkksu0cyZMzV//nztvPPOWrRokZ5++mm1bt1ajz32mC6++GKNHTt2vfKbOHGiHnjgAT3zzDNq06aNRo0apXvvvVcnnHBCo9a3PoUtugEAAFAZttpqK33hC19IXe7JJ5+sNYq9YMECLVu2TB07dlx33XHHHafDDz9cl1xyicaOHavjjjtOkrRw4UKdfPLJevfddzc4vyeffFIvvviiqqrCeO+yZcu0+eabb/D9xFB0AwAAVJjGjkhnpXPnzuv+btWqldz/06yqZu9rd9cLL7ygdu3aNXhfAwcOVJcuXfTmm29q7Nixuv322yVJP/rRj/SlL31JZ555pt555x0dcsghn7ltmzZtak0bqX5sd9epp56qyy+/vNHrmKbIc7oBAABQYVq1aqWePXtqypQpWrt2rR544IF1sYMOOkjXX3/9usvVB2bWNXLkSP3iF7/QihUrtP3220uSFi1apM0220yS1hXidQ0aNEgvvfSSJOmFF17QjBkz1j3ufffdp3nz5kkKnVimT59e2orWQdGNXPXu2j76D0DzYX8E0FyuvPJKHXLIITrwwAPVv3//dddff/31GjdunHbaaSdtv/32uummm+q9/bHHHqu777573dQSSbrgggt03nnnafjw4Q0+7rHHHqvZs2drl1120S233KItt9xSkvS5z31OP/7xj3XQQQdpp5120sEHH6zZs2c30doGTC8BAABRZjZK0ihJGjBgQM7ZoPD5xt0AACAASURBVBxcdtll6/7eeuutPzNiPXLkSI0cOfIzt+vTp4/uv//+1Pvv16+f1qxZU+u6vfbaS5MnT153+YorrpAURrEPOuggSWGay5NPPlnvfZ5wwglNeuBkXYx0AwCAKHcf4+5V7l7Vp0+fvNMByhJFNwAAAJAxim4AAIAK4PJaXUHQtErdthTdAAAAFeCjZR9p/vz5FN4ZcHfNnz9fHTp0aPR9cCAlAABABbjn/Xu05+Z7au7cuXmnUpE6dOhQq9PKhqLoBgAAqACfrvlUW2yxRd5poAEU3WVuzdr4T0hp8datrKR4mrTHT1Pq42et1O2L5pX181Xq6z3NytVro/F2beIzBvN+vZaaPwCUM97hAAAAgIwx0g0AAFCGBl34yLq/Z7WbX+u6aaNH5JITGsZINwAAAJAxim4AAAAgYxTdAAAAQMYougEAAICMUXQDAAAAGaN7ScGV2te26H2i0/JbvGx1NN61Y74v4aJvX9SW9/NV3Se7odd12uu56H240x6fPtwAWjLeAQEAAICMUXQDAAAAGaPoBgAAADJG0Q0AAABkjKIbAAAAyBhFNwAAAJAxim4AAAAgY/TpLlHWfXGz7mtbah/wUtd/3uIV0XindvGXaKl9vOvm73WuL/X5W7hkVYOxHp3blnTfKD/Vr6eGXpel7k9p+0Paa67Ux1+xKv5+kvZ+k7a/LliyMhrv3bV9NA4AeWKkGwAAAMgYRTcAAACQMYpuAAAAIGPM6QYAACjBoAsfaTA2bfSIZswERUbRDQAAUECxYh7lp7DTS8xsmpm9bmYTzGx8ct1lZvZBct0EM/ty3nkCAFDpzGyUmY03s/Fz587NOx2gLBV9pHt/d59X57pfu/tVuWQDAEAL5O5jJI2RpKqqqnhvSQD1KnrRXXhpfWvT+t6mWbpiTUm3X7Iy3rd3ZUpf3QG9O0Xjaev/2vRF0fjWm3SJxtu3jf8Yk7Z90uKfLK/dR3tV0kd4btI/vG/3DtHbp6EXd8uS1ie7uk/1rEXL641361Da6yWtz/X0eUuj8bQ+153at47GP1ywrKT7n5vSt3+jzu2i8TSlnpcAAEpR5HcYl/SEmb1kZqNqXP9dM3vNzG41s5713bDWz2Dz+BkMKDL2VwBAS1Dkonu4u39e0qGSvmNm+0i6QdJWkoZK+kjSL+u7obuPcfcqd6/q07tPsyUMYMOxvwIAWoLCFt3u/mHy/xxJD0ga5u6z3X2Nu6+VdJOkYXnmCAAAAKyPQhbdZtbZzLpW/y3pYEkTzWzTGot9RdLEPPIDAAAANkRRD6TcRNIDZiaFHO9298fM7E4zG6ow33uapG/llyIAAACwfgpZdLv7e5J2ruf6k3JIBwAAAChJIaeXAAAAAJWkkCPdTWmNx3s1p/WBTutDnabU26f13U2Ttn5pfXEvfXxSNH7m7gOj8SH9ukbjnyxbFY2vWRvvC1y3z3ZdC5fE49v0rd0nvG3rsL2q+wGn9Vkv9fktxbyU564lKvX5+tvbc6Lx/baNd1cptS972nOa1ie7fdt4H+u07ZO2P7VpHd9+pfb5/iClz3eatL76aX246eMNIEu8gwAAAAAZo+gGAAAAMkbRDQAAAGSMohsAAADIGEU3AAAAkDGKbgAAACBjFN0AAABAxiq+T3drS+8NW4pYD3Cp9MdevGx1NJ7Wd3fBkpXR+CYpfW27dygt/9vHvx+N795vo2i8a4f4S3RA707R+GY9O0bj0+ctrXW5uk/vrIXLJUl9e8S3T9rzE+vbPGvR8uht03oOp/VELkdrFe+VvHpN/PVe6v524JCNo/G0Ps5p++OqJP/5i+vfL3t1jffZnr0o3sf71Q8WReNDNon3zf/g43if7Ibyrla3731dac9f+5Q+2NPq7K91pe0z9OEGkCfeYQAAAICMNdtIt5lt5O4fN9fjAQAANIVBFz6SdwqoAJmMdJvZcDN7y8zeMLPdzOyvksab2Qwz2yOLxwQAAACKKquR7l9LOk5SF0mPSDrK3Z82s89LulbS8IweFwAAACicrIrutu7+uiSZ2Vx3f1qS3P1lM4sf2QYAAABUmKwOpKx5vxfVicUPzwcAAAAqTFZF9yVm1kmS3P3B6ivNbCtJv8voMQEAAIBCymR6ibs/1MD170r6nyweMytpfXfbtLZM779rx/hT9MGCeF/djil9i9P61p5WNTAaT8tvpz7do/Fv3flSNP7Fqs2i8dvunxCNn/31XePxvbeqdbld2/A9dLONwiyotz9cHL39Dv27ReNT5yxpMLbFxp2jt22JWineK/mDjxvenpI0aV78+dp+4/jzVf38NyStD3Ta/vz+gtBn+idPTKo3Pmbk0Ojt0/pYD+wR71u/cMmq+O1T+t7/7O/vRuN7DYpv3216xvuEb9Envk+k7TNp23/KrE+j8bT9uVRp77eoTHQ+QTX6dAMAAAAZo+gGAAAAMpZZ0W1mrc3sB1ndPwAAaB5mNsrMxpvZ+Llz5+adDlCWMiu63X2NpCOzun8AANA83H2Mu1e5e1WfPn3yTgcoS1mfBn6cmV0naaykdUdAufvLGT8uAAAAUBhZF917Jv//tMZ1LumAjB8XAABgvVVal5G09Zk2ekQzZYJqmRbd7r5/lvcPAAAAlINMi24z20TSzyX1c/dDzWx7SXu4+y1ZPm5Tat0q3oc7Lb50xZp4fOXqaLxbx7bReOd2pfXxfmnWgmj8hffjfY/P33fLaPzSR96Kxn/8le2j8arNe0bjW/eK903es3+vaPzB1z+odXnepytqXX/okE2jt4/14ZZK68Wd1nM47bVXjlzx9U7bngNS+kxPeH9hNJ7Wp/q16Yui8bQ+0107hP11t63ir+uGzFuyMhr/bkrf+4d/sHc0/s/3SjtA7uBt+0bj46bOi8YH9Ipv/wUp65/2fpl1H+40sR70ACpf1u8At0t6XFK/5PJkSWdn/JgAAABAoWRddPd29/skrZUkd18tKT70CwAAAFSYrA+kXGJmvRR+NZaZ7S4p/vssAAAAMhU70JKDLLORddF9jqSHJG1lZuMk9ZF0bMaPCQAAABRK1kX3G5L2lTRYkkmaJE49DwAAgBYm6wL4WXdf7e5vuPtEd18l6dmMHxMAAAAolExGus2sr6TNJHU0s10URrklqZukeE8oAAAAoMJkNb3kS5JOkdRf0q9qXP+JpB+uzx2Y2TRJixW6nax29yoz20jhlPKDJE2TdJy7xxtNF1zvru2j8cXL4n28bxn/fjR+9t5bReMzPlkajXfpEO/zfdH/xftwv/rie9H4wr03j8bP/MNr0fje28b7cL+/ML5+Vz08udblGVpW6/p2reM/Bh2+Y79ovBSV2Ic7jSm+3vMWr4jevmfndtF4Wh/uZ9+fH41Pnh9/PS16a1Y0nubOV2ZE49tsFM9/1CFbR+PH3/xCNH7NsTtH47sPiO9vT707Jxo/85p/R+OTrjsmGu/UunU0nnpehJR4147xj8S0fbIl9tYHsP4yKbrd/Q5Jd5jZV939jyXc1f7uXvNsChdK+pu7jzazC5PLF5SSKwAAAJC1rA+kHGdmt6jpzkh5pKT9kr/vkPSUUoruyfMn6eA794stkqm1KSMfrUocOfnwk+XR+GPTO0bjC5bGz/C2MGWkfcXq+MjRkgHxkcHRz/02Gp/3aTy/NybEz0DXqW18ZGyqap9RcpnekSRN0TmSpMuejo8sXvtKfGQVGyZtf121Zm309m1SfplYnXL7T5bHX+9LV8Vf72n765xlb0uSbn7txHrjrSz+ftCxbXz9VqyOr9/sFfFfCv7rifgZNdO276cr4ttvYd/4GUGPuDf+fmAp2ydt+6dp3Tp+/2nj1GmPzjg30LJlfSDlbWr8GSld0hNm9pKZjUqu28TdP5Kk5P+N67uhmY0ys/FmNn7VqlWNzx5A5thfAQAtQdYj3b3d/T4zu0gKZ6Q0s/U9I+Vwd//QzDaW9Fcze3t9H9Tdx0gaI0m77lrlT5z01Ibm3WTS5hB2ah8fiU2b0z3mhWnReNqc7r++PTsaf2pqfGRqxvwl0fiTr74ZjV94zL7R+J9eieeXNqd7+95dovGL/zix1uXqEe5tkkMRfrjXkOjts5zT3RQ6nlr8sbUN2V9LndO9YEn8l5OS53Qvj+/vf5h8siTp9J3uqjeedgxB2pzu6Z/Ej8G4Z1x8zvg1B8fndPfqGt++L38QP8QmbU73QxfH53S3SRmJTnu/TZP3nO5y2F8BNF7WI92NPiOlu3+Y/D9H0gOShkmabWabJve1qaT4UTsAAABAAWRddNc9I+XvJJ2VdiMz62xmXav/lnSwpInJfX0jWewbkv6cRdIAAABAU8p0eom7v2xmtc5ImZwgJ80mkh5IDpppI+lud3/MzF6UdJ+ZnSZpujilPAAAAMpApkW3mbWW9GWFvtptJB1sZnL3X8Vu5+7vSfrM5EJ3ny/pwAxSzUzaHMRZi+LdR7p1iHfn2LZXfI7n6L+/E43/8834DJ2bT9glGt9so3h3lKsHdo/GD9p6k2j8g0/ic3D/8vJH0fhJpw2Lxt9+5uVal1fssDhc/0a4vl/K+n+wID6HdrOe8e0TszKlE0VLlNbXPm3Od1p3j7Q5+i9Njc9Z7pDSLeeh98KPi5t1q39u9Lvz4/n/6sl3o/E+PTpE449/f69o/Kgxz0fjt54Y3x9+/uf4oTePXHJIND475f0wzYCUPuxpc75L7aNNH24AMVkfSPkXScslvS6JCgIAAAAtUtZFd3933ynjxwAAAAAKLesDKR81s4MzfgwAAACg0LIe6X5O4YDIVpJWKRxM6e7eLePHBQAAAAoj66L7l5L2kPS6u5d2fl4AAACgTGU9vWSKpIkU3AAAAGjJsh7p/kjSU2b2qKR1vbDSWgYCAAAAlSTrontq8q9d8q/ipPV97dQ+3re3b/d4X900/TrH+0B/1CXe9/fjj5dG42m9og+/8blofKeBPaPx3zwzLRq/9oa/RuMde24UjX/jrvhLfOLvv1Pr8gkPPiRJuvuicH3PTvGX7WszF0Xjnds1/Phpr412bbL+Iar8pPVFT5PWN/3fU+ZF4306x/uET/l4cTTevX3ou7/n5r3rjT/93nvR26ftr2l99Tc5+KfReJettovGVx73mdMn1HLYbptH49c99340/t3dB0bjS1fF32/T+nS3b8s+1VhmNkrSKEkaMGBAztkA5SnrM1L+JMv7BwAA2XP3MZLGSFJVVRVTRoFGyPqMlH0knS9pB0nrhnTd/YAsHxcAAAAokqx/a7tL0tuStpD0E0nTJL2Y8WMCAAAAhZJ10d3L3W+RtMrd/+nup0raPePHBAAAAAol6wMpVyX/f2RmIyR9KKl/xo8JAAAAFErWRfcVZtZd0n9LulZSN0k/yPgxAQAAgELJrOg2s9aStnH3hyUtkrR/Vo8FAAAAFFlmRbe7rzGzIyT9OqvHKIK0Xstp0vpgt25l0XjfHvE+370+aRuNn3vk4Gj80XfmROP/tU+8r+7k+fG+wsP7x/tsTz92eDTes0u8b/Lpu8ZnM9Xtk17dG7v6+rmL433ON+kaf/wenePbH7WtcdfiZasbjKf12U7rmz8v5fncZuMu0XiblP1xm76bRuOjnwv7+8Klq+qNj/pCvM/1xQdtHY2PmxbvM95rx6HR+AF7bRmN/+GND6Pxr27fNxrv1C6lN31KH+0pcz6NxtPeL9PipSr1vA0AKlvW00ueMbPrJI2VtKT6Snd/OePHBQAAAAoj66J7z+T/mqdBc0n06QYAAECLkfUZKZnHDQAAgBYv65FuJa0C656R8qcN3wIAAACoLJmeHMfMbpQ0UtJZkkzSsZLiR94BAAAAFSbrM1Lu6e4nS1rg7j+RtIek+OH5AAAAQIXJuuhelvy/1Mz6KZyhcouMHxMAAAAolKzndD9sZj0k/a+klxU6l9yU8WOWleq+0A157M1Z0fjg3l2j8asenhyNP3HOPtH41+98KRr/3j7x71BXnHt1NH7KxWdG4xcfuE00/vdpc6PxbfrG+y4/Nbn27av7J1df36dTvA/3Dv27ReOxvr1pPXvXrPVoPOuew3lobaauHRv/tpS2TdP6dA/o3SkaT+vDfNo9E6Lxzu3Duu00oHu98UsfnxS9fZrz9tkqGp9/ws7R+IMvx99vqvNvyL+mx/uE/+J38W6xL155WDSe1kd9+rz4eQHSnt+08yakvV/ThxtATNbdSy5P/vyjmT0sqYO7L8ryMQEAAOoadOEjeaeAFi7TotvMOkg6U9JeCqPcT5vZDe6+PMvHBQAAAIok6+klv5O0WNK1yeXjJd2p0MUEAAAAaBGyLroHu3vNSYT/MLNXM35MAAAAoFCyLrpfMbPd3f05STKz3SSNy/gxAQAA0Ehp89+njR7RTJlUlqyL7t0knWxm05PLAyS9ZWavS3J33ynjxwcAAAByl3XRfUjG9w8AAAAUXtYtA9/P8v6LoNReyml9YbffON4Hum+PDtH4FV/dMX7/3/tjNH7JacOi8XcXLonGr7nhvGh8z/69ovFNusfX7/q/xPuQb9WjczR+4JCNa13u+WLbWte/Nj3e4TLt+U17fZRy20rs0521tD7NC5esKun+Rx+2XTT+1T+EPt9vzPyk3vir7y+M3v7SL24bjQ/Y5+xoPG1/vPRLg6Pxv02dH42P2j3et3/XTXpE4+NnLIjG99mqTzTeqXtpfbKz3qfS+rwDqGxZn5ESAAAAaPEKXXSbWWszeyU5sY7M7HYzm2pmE5J/Q/POEQAAAEiT9ZzuUn1f0luSas6xOM/d788pHwAAAGCDFXak28z6Sxoh6ea8cwEAAABKUdiiW9LVks6XVPdIw5+Z2Wtm9msza1/fDc1slJmNN7Pxc+fNzTxRAI3H/goUX639dC77KdAYhSy6zewwSXPc/aU6oYskDZH0BUkbSbqgvtu7+xh3r3L3qj6940e7A8gX+ytQfLX20z7sp0BjFLLoljRc0hFmNk3SvZIOMLPfu/tHHqyQdJukeD87AAAAoAAKeSClu1+kMKotM9tP0rnu/nUz29TdPzIzk3SUpIk5pikp+76u3Tq2jcbT+r7ut218ROIfl8dP5frnSbOi8aO32zQa/879r0XjR2zXLxq//aXp0fgBwzaPxheuWBmNn/i72j+mvDVnca3rr/lKvM952vbv2rGQu1hhrXHX4mWrG4yXuj1j9y1JndrH+zzPWrg8Gm/dOqUv/5owW27Wp/Xfz5y58b73X7704Xj8rG9G41/ZYbNo/J5XZ0TjJ+wU31/T9oeBKX3St23TNRpv37a0caK8e9+nvb4AVLZyqwjuMrM+kkzSBEnfzjkfAAAAIFXhi253f0rSU8nfB+SaDAAAANAIRZ3TDQAAAFQMim4AAAAgYxTdAAAAQMYougEAAICMUXQDAAAAGSt895JK165N/HtPWjzNytVrS7r9kYP7RuOzF6+Ixn9z9Oei8VvGvx+N79CnSzR+/M7xvsP/fC9+uuJLDtq21uXTH+1U6/pO7eK7SJZ9d0t97stRa7NMe5uXet9zU17vQwf2iMa7dQiPP2zgRo16/P/76WHR+JSFi6Px6fOXRuNf3HLjaDxt+y1dGe+D3qZV/DWddV/7rPtwA0AMRTcAAADW26ALH2kwNm10/KR7LVnLG0oDAAAAmhlFNwAAAJAxim4AAAAgYxTdAAAAQMYougEAAICMUXQDAAAAGaNlYMGtWevReFrf2bRez1ts3Dkan7VoeTQ+qGOnaHzCBwvjt+/RIRo/ZPt4n/Cpc5ZE48MH9Y7G+3av/fid2oW+20P6dY3eDi3Trlv0LOn2rSzsr+3b1r9fPnHOPtHbp/XZ/touA6Lxe1+ZHo1v0yP+uv94ycpofLt+3aLxLPvaA0DRMdINAAAAZIyRbgAAUBFiJ20B8sZINwAAAJAxim4AAAAgYxTdAAAAQMYougEAAICMUXQDAAAAGaN7ScGl9eHO+v67dWgbjXftGH8JDW8X75Nd6vr17NyupPuv2+d7+ao1ta5P62NeisXLVkfjadsWza/U5yxp091g//wVq9bG779D/P6XrlgTjaf18U47L0CaKbM+jcbpfw+gJWOkGwAAAMgYQ2kAACDKzEZJGiVJAwbEfzFByxbrlT5t9IhmzKR4GOkGAABR7j7G3avcvapPnz55pwOUJYpuAAAAIGMU3QAAAEDGKLoBAACAjFF0AwAAABmje0kLt3J1Sl/gEntFl3r7eYtXROO9u7Yv6f47tW9d63L7tuHygN6dSrrf9UEf7vKT9XOW1le+7ut1Q+OlPn4a+nADQMMY6QYAAAAyRtENAAAAZIyiGwAAAMhYoYtuM2ttZq+Y2cPJ5S3M7Hkzm2JmY82sXd45AgAAAGkKXXRL+r6kt2pcvlLSr919G0kLJJ2WS1YAAADABihs+wQz6y9phKSfSTrHzEzSAZJOSBa5Q9Jlkm7IJUEAANDsBl34SN4pAI1S5JHuqyWdL6m6p10vSQvdfXVyeaakzeq7oZmNMrPxZjZ+7ry52WcKoNHYXwEALUEhR7rN7DBJc9z9JTPbr/rqehb1+m7v7mMkjZGkXXetqncZBO3aZPu9a83a+OZP6wuc1od71qLl0Xjf7h2i8c/k57WvL7VvMdKxv66/rPcHAEB2Cll0Sxou6Qgz+7KkDpK6KYx89zCzNslod39JH+aYIwAAALBeCjm9xN0vcvf+7j5I0tck/d3dT5T0D0nHJIt9Q9Kfc0oRAAAAWG+FLLojLlA4qPIdhTnet+ScDwAAAJCqqNNL1nH3pyQ9lfz9nqRheeYDAACADZfWeWba6BHNlEk+ym2kGwAAACg7FN0AAABAxii6AQAAgIwVfk43ylvWfa436tyupNvX7VNuVv/1QHPo1L51SbdP68O9cvXaaJzXPQBkh3dYAAAAIGOMdAMAgGYV62JR6R0s0HJRdAMAgMJIaysHlCuKbgAAAOSulC9csV9IitIfnDndAAAAQMYY6QYAAE2KKSLAZ5m7551DpsxsrqT3m+nhekua10yPlRfWsbwMdPc+eSexvhq5v5bb80W+2SrnfAu7v5rZKEmjkouDJU3K+CGL/DySW+NUcm7rte9WfNHdnMxsvLtX5Z1HllhHFE25PV/kmy3yrQxF3i7k1jjkxpxuAAAAIHMU3QAAAEDGKLqb1pi8E2gGrCOKptyeL/LNFvlWhiJvF3JrnBafG3O6AQAAgIwx0g0AAABkjKIbAAAAyBhFNwAAAJAxiu5mZmaWdw5ZMLNNzKxD3nk0FzNj3wGAJlb9GVmpn5UtFc9nwGngM2ZmX5S0v6Qpksa7++tm1srd1+acWpMxsxGSvinp25KW55xOJpJ1HCFplqTH3f35nFNCwsz2k7SxpDbufnfO6aQqw3y/KGlbSa3c/VozMy/wEfhlmO9+KqPXQzPYWNJshfpkVZE+L81sV0kL3P29vHOpy8w2V9hubd19ScG2224Kz+e4vHOpq7m3G6N1GTKzfSVdL2mBpAGSHjSzA9x9baWMlJrZoZJ+Julqd59XJ1YR32zNbHdJv5L0gqRPJD1sZofnmxUkycz2l3SPwv51jpn91sz65ZxWg8ow370k3a3wZXqkmV0rabiZFXLApgzzLavXQ9aSwY0HzGyMpJ+Y2aCifF6a2ZckjZXUpcZ1hfiMS7bbo5KulXSbmQ0u2Ha7QwUckMtju+X+hFS4bSQ96O7/6+4/kXSRpPvNbP/kiS3EDttYyYfDhZL+5u5Pm1lPMzvZzEaa2Y7u7uW+jonNJD3j7re7+9WSTpN0pZkdJhXnjbelSbb7oZL+x92vkrSXpO6SLjCzTWosUwjllm9imKTr3P0WSQdJWiTpGElfyDWrhpVNvmX6esiMmW2lUPxcJOlOSUskjTWzbfIuIJMvR9dJOsPdXzOzjkmodRLPJTcLNpc0WtJ3JV0q6XlJ/zCzHQqw3faSdKuk/3L3l8ysS3J9x+T/PHPbTNKVaubtRtGdrbmSelZfcPf7FKZg3GBm2xX5J8/1NE/SHyR9amb/LemvCh8c+0i6y8x2qYB1lKSpktaY2aaS5O4PKXww3GpmVRWyjmWjuhBJtvvLkgab2SbuvlzSGZI2kfTjGssUQrnlm3hV0u5mtm2S7+WSlko6Md+0GvS6pD3KId/kuX5F5fV6yNI8Sf9w939KelrSzyX9UdLvzGxgzlMlDlbYF54zswGSfmNmv5F0hZltnlduHsyQ9KykyZLmuPsvFYrwJ5L9IM/ttoPClJJ5ZjZQ0hgzu1HhOa3+MtXsXyyTon+epH+rmbcbRXcTM7O9kjl6kvQ3STub2a+r40nh/UdJ2+eQXpNI1vEAd18p6WaF+VBHSrrN3Ue5+3ckPSRpzzzzLIWZHWpmRyUX35DUVdJ5ZtYqmSP6Z0n/K2l4bkm2XH0kycxaSxqv8NzsZGYd3X2pwvEFu5nZETnmuI6ZbW5m7ZM3+mdV/Hxrfgi+JelNSXuZ2abuvkLSTyUNM7NT8sivLjPrl2zfLpL+JeltSXsXON9BZtbDzLopvB66qMCvh6wlnydfV5i6t62ZXZgUky7pKkn/J+kkM2vd3AWame1tZkdKukzSe5KukfSEwj4xTtIKST9MXn/NndvhZvYDM2srqZukU6q/pLn7b5Jcf2hmHXLI7QgzO0PhF4txkr6X/P+cwsj3y5KuM7Ouzf3FMnk+r5LUT9JGkr7ZnNutkPPcylHy5PSRdJOkrmZ2hrs/amaHSBqfFGpnJ4t3lLRdXrk2Vj3reLq7P2Zmt0h6xd2frbF4e4XiouyY2cEKc7j/S5LcfYWZjZL0oKRfJ/+mSeogqW9OabZIFo4hONfMxktapvAG+UdJ3w9he93dPzKzv0lak2OqktbNGbxS0jMK+8M5CnN4z1Zx8x1sZmPc/VN3/9DM/i3psBC2ce7+tpn9RVLuB2kl768/Vii0V0g6X6FIOy6JFy3fL0n6hULx0VXSKQqjumepgK+HLCU/33eSS4FMNwAAD85JREFU9P8UapFPFZ63J8xsmbtfk4yEviDpSHdvtu1RI7cbJbVV+My+WGEk9AZ3vyZZbh9JJyVf7ppN8hl1uaQL3H2VmV0o6V9mtsbdr0wWu0/SRckvKM2d208lXejuS5P64HRJz7v7zckyHypMv23u7bavwvvx9919qpmdJ+nfyevtV8limW43iu4mknxTmmNmv5fUWdJoM+vs7veb2S4Kc7lvTRbfXdJX88q1sepZxyvNrJO7/0lhxEaSZGYjJR0o6fh8Mm08M9tT0m8V5qA9ZWZdJXVJPgi/rDDn8LLk+iGSRuaYbotiZjsobP9vSlos6QRJ90s6OlnkxLCYfZDEbs4jT2ndF9T++s9cy7ckfUNh3uAekm6Q9PVk2dzzTfL4gsIHzkcKU8budfdP3P3BZMBnD0knm9kESV+TtF9uyWrdPNvfSDpVyYGTkka4+z3JqPduKl6+v1QosD+SdJ6kDu5+o5l9rIK9HrKW/Hz/qZndofAF4ziF6ZgHShpnZqvd/XpJmyp8Eewq6dPmGBmtJ7fDJXVy93PMrH2NRQdK2ix5vS1pjtySz6g7JR3u7i+YWW9JMyUdJekRM1sl6WGFX5p3NbOe7r4g67wayG0jhRkVd6r2gZT7StpS4YvNyubILbGrpJvd/XEL04S6KHyZ+q2ZLVeYnbCHMtxuFN1NJPlm7Aqjn5MVfkq53MyGSJqhsNMOVxgZHe3uk/PKtbEi67iNpDXufpWFrh5nSTq5HNdR4eemTyR9mKzX1ZJkZnMUpsycofANfWtJb7n71LwSbYHaSHra3f8tScmb5t4KxxUcozDncgdJQyUdmOfrL/nwnWFmNeda/k/ygfiMwhfvlxUO8Ns573wTHSV9RdIchV9z2prZnTUK73EKXzQHS7rW3d/JMVdJ2kXS5e7+tCSZ2QEKBds97v4nM/unwjS+ouS7naSz3P0fZjZI4SDKXyRfaH4m6SmFD/yivB6ay2qF7i23KLy/9leYmz/SzIYpHBx7nLsvzjG3WyWdkXzxXynpIjM7W9LJCiPdnzZjTvMlrZK0qZn1Unj/W60wDfJmhcJyG0lVClMnmqXgbiC3+xVGsxdJeszM7lLYZt+VdKK7L2zG3KSwndolf98r6UNJ7yq83g5WeK/YU1luN3fnXxP8k2TJ/7so/KwihZ9Ylkv6ed75NcM6/jS53FHSpnnnWuJ6Hq/whWKiwofA5gqjqLdL6pN3fi3tn8LBuccr/LrynqQfJNf/XGGu4OWSjsk7zxr5Hi7pBwo/S98r6Yd14hcpfIi3zzvXGvmeKckkbZRcN0zS3xU+HLsl13XOO9ckjyMUfq7uoFAQVb8v7Sbp3hrLdcw71xr5nlHjcpfkveQ8hS+JP1IYYeuad645bZ+tanye/LdC0XZpcrmdpN4Fym2ppOuTy7dJ2jGnvHZO3gtnJp9RrSSNUmhRvHmyTM8C5XaqQivPLRSmV22XU247SpqUvC9/M7lu2ySnI5tju3EgZQnMbCcz21n6zFHmWyeT9UcqzDk9Ppl7WHY2YB1PMrMR7r7M3T/KI9fGqrmOkuTu9yj8DHyju9/k4ejwhxS6CnRp4G7QxCwctNpFYc7nTxRGBg+V9G0zu1dhlPi3CqMVu+aWaA015lq+6e6rFFpqftvMLqix2D0KIy7N+bNqvWrk+44HH0uSu7+g8LPr0ZK+bGbflXSjmbVJps7kme9PJc1w9+XuPr3G+5IrfKjLzE6SdLHlcPBdTTXyfb/6Og+jold4aCX7hkLx9r5CsdkSLVOYPnKGQnevKxQOJP22u6/0Oud/yDm30ZIGmdlxkk5z94l5JOXuryocZ/GL5DNqrbuPUfgFtk+yWHOPIsdyu1VSL4UvyBe5+1s55TZR0rkKX9C3SK6brHBCpu7JYpluN6aXNJKFHs13SHrUzG5w93GS5O6vmNkShZGs0939ATN7RuEn5rLSiHXMZUcqRWQd/2ShO0a1AxV2yiU5pNki+WfnVR6lcNawwWbW3d0XSZKZrZXUxsxaezMebFXXBs61/LykHgonzipKvt2TnD6WtMLdnzGzbyoc9LdaYb706gLmO19hBHKmpKlmdqzCgaonFuz10F1h+toshTak1fZXGFHtqAKeQCRrHg7WnSHpEknfcfe/JPPf854O1FBuB0ia4jmf7dHd31TooiJJMrOvSuot6YMknlu7yQZy66Own+btUYWDry8zs+ovwzsr/HKa+Xaj6G4EM2sn6csKRza/p3CwjqoLNkl/lvQ7D83gq9vLlRXWUar+wDaz7ymcEOdEd5+TV74tWPW8ytskjUoO+FthZj9S+PnyPElH5VlgJYo817I+9c2/XKbQReIxhS/V/RU6ER2UjMrmqaF8l0j6i8KH6UEKo30nu/vbeSWaiG3fR83sTwrT1s6Q9PUCvB7ydJOkP7v7S8nlf+Zd1NZQN7enCpRb9UHb31QYwT3W3WfnnNI69eQ2K+eUlAwc/M7MJiocC9Re4f343eZ4fMvxy1BZS47KXaFwYORXFD5M7/bQ2L/mcpbnN85StOB1/L0nB+sly4yUNLEARUeLZOFMdce6+2gLJ2G6XNKt7v7d5AvR4+4+Kd8sg2Sa0gMKc1F/onBw2OkKIymj3X1Gc3YTSNNAvqdI+qLCvPQtFQ4CzX3UUUrN9xKF6UZnJyNtuYvke6DCWfDOkjSmKPnmrcifJUXNLSls95U0qwBfNGspcm55oejeAGY2VElfyZpzkpIuF0cqFGxXKnxQzShKIbAhWMfPrON0bzldBArJzPopdHd4RqEP850KB/rd6+6/zzO3+pjZ9pL299DurPq6xxV6v75ctA/vBvJ9TNK5ec1ZjWkg3ycUfo2aWaRtK0Xz/ZbT/QhoUZhesp4snJRjjMK0iv3M7JfufpskufsUM3tQ4fTn9ym01dolt2QbiXX8zDpup9B+Djkq8pzP+hR5rmV9IvMv8zyArUGR7buqaNtWiua7LLekAOSCojtF8vNIZ4WfAb/j7g+Z2e6Sfm9m7d39Rkly93fM7FSFN9Pd3H1KfllvGNaxwXUcVk7rWOGKPOezXkWea1mfIs6/jCFfAOWGojtFMnLyqYXTTnczs7bu/pyZfU3SH8xsubvfbqHTxRBJR5fb3F/WsTLWsZJ5aNs4o3pqRtEL7hreU3gtlct8RvLNVrnlC6AJMad7PZnZmQp9gb/v7p8k1+2lcMbCrxXlQKNSsI6VsY4AAKB4ODlOiuQnQbn7byV1UjhBRPdkpPRpSa8ptAUrW6xjZawjAAA1mdkpycHopdzHkWb2mplNMLPxyUAVGoGR7nqY2WCFkxiMl7S2Zv9fC2fCW6Zwwog2ks6RtK+7z8wj18ZiHStjHQEA5c3M2nhGJ54ys6cUOhGNb2w+Fs4MvMTd3cx2knSfuw9p+mwrH3O66zCzoxXOTPRB8m+8md1ePRXB3b+WHGjXT6H37hHlVqixjpWxjgCA/JnZIIUTSj2v0NVrssIJmpaa2aWSDlc44+gzCq0iPSmGn5E0XNJDZjZZ0sUKPd3nK5yMbbaZXaZwyvJNJW2rMEC0u6RDFT7bDnf3VWa2q6RfSeqi0HnolOS+qyTdZWbLJO0hafu6y7n7R3XzkfTL6vVz909rrG5nSYzWNhIj3TWYWVtJv5f0G3cfl7R22l2hp/P/enLa6RrLt3f3FTmk2misY2WsIwCgGJKie6qkvZLPnFslvenuV5nZRu7+cbLcnQqjxH9Jitw33f3MJNZT0sKkID9d0nbu/t9J0X2QpP0VCuZnJX3V3R81swck3SHpEUn/lHSku8+1cEK3L7n7qTVHupPPxthy6/KpZx2/IukXkjaWNMLdn23SjdhCMKf7s7opnBxFCmcSe1jhm+fxkmRmw8zs80l8ZfOn1yRYx8pYRwBAMcxw93HJ37+XVD3veX8ze97MXpd0gKQdatxmbI2/+0t6PFnuvDrLPeruqyS9Lqm1wqi6ksuDJA2WtKOkv5rZBIUR8/715Ji23Nh6biNJcvcHkiklRymcFRiNQNFdQ/Ki/pWko81s76Qt2dOSJkjax8w6Kvz08mGyfNn9TMA6VsY6AnloooOyTkwOynrNzJ6xcKp0oNzV/RxxM+sg6beSjnH3zymcb6BDjWWW1Pj7WknXJct9q85y1WdQXqvaJ4FaqzBN2CS94e5Dk3+fc/eD68kxbbkl9dym9kq5/0vSVmbWO21ZfBZF92f9W9ITkk4ys33cfY27360w97efu//ay/+kBqxjZawj8BlmluWxOqco7EPrrZ58pioctLyTwojZmKZJDcjVADPbI/n7eIWBnurCeV5yMOIxkdt3V3LWWknf2MDHniSpT/Xjm1lbM6seKV8sqet6LNcgM9u6ugNY8gtx9bxzbCAOpKzD3Zeb2V0K31ovMrMhCt8y+0j6NHrjMsE6VsY6ojK1gIOynqmxus+p/p/BgXLzlqRvmNn/kzRF0g3JPnuTwjSQaZJejNz+MoUTtX2gsF9ssb4P7O4rzewYSb8xs+4Ktd3Vkt6QdLtCi9zqfbah5WK+KulkM1ul0PVrJL8QNw4HUjbAzNopfGB8S9JySde4+yv5ZtW0WEegeFrCQVk11vVcSUPc/fQm2nxAs0v22YfdfcecU0HBMdLdAHdfKekfZvavcLFsTju93lhHoLDqHpT1PUlXKRyUdb7CCZ42Uhih+kuyXN2Dssaa2aYKo91Ta8QeTUaz1+egLCXLfFRPjmnLNXhQliSZ2f6STtN/DjgDgIpG0Z3Ca5xQpVKxjkDhxA7KqnL3GcmodeygrF+5+0Nmtp/CT9fV1h2UZWaxg7L2UFzacg0elGXhBBs3SzrU3ZkbirLm7tMUvoACURxICQDFU8kHZQ2Q9CdJJ7n75A3MDQDKFkU3ABRP9UFZrylMI7nB3RcqtBx7XdKDWr+Dsv6tcIDjekumZB0j6Uoze1Wh1eaeSfh2hYOyJihMJ2louZhLJfWS9Fszm2Bm6316agAoZxxICQAFwkFZAFCZGOkGAAAAMsZINwAAAJAxRrqRqyY6rfQQM3vWzFYkfX8BAAAKhaIbqcrgtNIf6z99jAEAAAqHPt0tQAs4rfQcSXPMbERTbjcAAICmwkh3yzFY0hh330nSJ5KqT898nbt/IemU0FHSYTVu08Pd93X3Xyr0Cd7d3XeRdK+k82sst5WkEZKOVDh73j/c/XOSlkkakZwu+lpJx7j7rpJulfQzd79f0niFAn6opNX1LddAPgAAAGWDke6Wo+JPKw0AAFBUFN0tR0WfVhoAAKDImF7SclTsaaUBAACKjqK75ajY00qbWV8zm6lwEOfFZjbTzLptSI4AAABZ4uQ4LQCnlQYAAMgXI90AAABAxhjpBgAAADLGSDcAAACQMYpuAAAAIGMU3QAAAEDGKLoBAACAjFF0AwAAABn7/23TAtT6pCfLAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Create figure with several subplots\n", "n_bins = 25\n", "n_param = log_likelihood.n_parameters()\n", "fig, axes = plt.subplots(n_param, n_param, figsize=(12, 12))\n", "for i in range(n_param):\n", " for j in range(n_param):\n", " ax = axes[i, j]\n", " \n", " # Create subplot\n", " if i == j:\n", " # Diagonal: Plot a 1d histogram\n", " bins = np.linspace(np.min(chain[:,i]), np.max(chain[:,i]), n_bins)\n", " ax.hist(chain[:,i], bins=bins, normed=True)\n", " ax.axvline(real_parameters[i], c='g', label='True value')\n", " ax.legend()\n", " elif i < j:\n", " # Upper right: No plot\n", " ax.axis('off')\n", " else:\n", " # Lower left: Plot a 2d histogram\n", " xbins = np.linspace(np.min(chain[:,j]), np.max(chain[:,j]), n_bins)\n", " ybins = np.linspace(np.min(chain[:,i]), np.max(chain[:,i]), n_bins)\n", " histplot = ax.hist2d(chain[:,j], chain[:,i], bins=(xbins, ybins), normed=True, cmap=plt.cm.Blues)\n", " # Fix the aspect ratio\n", " ax.set_aspect((xbins[-1] - xbins[0]) / (ybins[-1] - ybins[0]))\n", " ax.axhline(real_parameters[i], c='g')\n", " ax.axvline(real_parameters[j], c='g')\n", " \n", " # Customise tick labels\n", " if j > 0:\n", " # Only show y tick labels for the first column\n", " axes[i,j].set_yticklabels([])\n", " if i < n_param-1:\n", " # Only show x tick labels for the last row\n", " ax.set_xticklabels([])\n", " else:\n", " # Rotate the x tick labels to fit in the plot\n", " for tl in axes[i,j].get_xticklabels():\n", " tl.set_rotation(45)\n", " \n", " # Add labels for subplots at the edges\n", " axes[i,0].set_ylabel('parameter %d'%(i+1))\n", " axes[-1,i].set_xlabel('parameter %d'%(i+1))\n", "\n", "plt.show()" ] } ], "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.7.2" } }, "nbformat": 4, "nbformat_minor": 2 }