{ "cells": [ { "cell_type": "markdown", "id": "af33f23a", "metadata": {}, "source": [ "# Adaptation of Hamiltonian Samplers\n", "\n", "In Hamiltonian samplers, we construct proposal by emulating Hamiltonian dynamics. This involves using an integrator (in this case, we must use a symplectic integrator to preserve detailed balance). \n", "\n", "As part of the symplectic integrator, you need to choose the `step_size` and number of steps `L`. It can be a little tricky to tune these by hand... Luckily we have some built-in scheme to adapt these parameters\n", "\n", "Ones which find optimal hyperparameters during warmup, and sets them statically for the rest of the run:\n", "- **Initial Step Size**: tunes the initial guess for the `step_size` to get a target acceptance rate and have no numerical error.\n", "- **Dual Averaging**: tunes the `step_size` to reach a target acceptance rate.\n", "- **ChEES (Change in the Estimator of the Expected Square)**: tunes the number of integration steps `L`.\n", "\n", "And those which adjusts the hyperparameters throughout sampling \n", "- **NUTS (No-U-Turn Sampler)**: Sets the number of steps `L` dyamically\n", "- **WALNUTS (Within-Orbit Adaptive Leapfrog No-U-Turn-Sampler)**: [Coming soon....] Sets the `step_size` and number of steps `L` dynamically " ] }, { "cell_type": "code", "execution_count": 1, "id": "49ef4f2a", "metadata": {}, "outputs": [], "source": [ "import hemcee\n", "import jax\n", "import jax.numpy as jnp\n", "\n", "import corner\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "id": "3f0aa647", "metadata": {}, "source": [ "For this tutorial, we'll sample the Rosenbrock potential. It is known for having extremely narrow regions of high probability, meaning you have to tune the step-size and integration length very carefully." ] }, { "cell_type": "code", "execution_count": 2, "id": "579a6081", "metadata": {}, "outputs": [], "source": [ "def log_prob(x: jnp.ndarray) -> float:\n", " \"\"\"Rosenbrock potential\"\"\"\n", " x0 = x[..., 0]\n", " x1 = x[..., 1]\n", " return -1 * ((1.0 - x0) ** 2 + 100.0 * (x1 - x0 ** 2) ** 2)\n", "\n", "total_chains = 100\n", "dim = 2\n", "\n", "keys = jax.random.split(jax.random.PRNGKey(0), 2)\n", "initial_states = jax.random.normal(keys[0], (total_chains, dim))" ] }, { "cell_type": "markdown", "id": "09a57baa", "metadata": {}, "source": [ "## Default Behavior of `hemcee`\n", "\n", "The `hemcee.HamiltonianEnsembleSampler` and `hemcee.HamiltonianSampler` by default having Initial Step Size, Dual Averaging, and ChEES all on!" ] }, { "cell_type": "code", "execution_count": null, "id": "46457f53", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/clarkmiyamoto/nyu/h-emcee/src/hemcee/backend/backend.py:47: UserWarning: Explicitly requested dtype float64 requested in zeros is not available, and will be truncated to dtype float32. To enable more dtypes, set the jax_enable_x64 configuration option or the JAX_ENABLE_X64 shell environment variable. See https://github.com/jax-ml/jax#current-gotchas for more.\n", " self.accepted: jnp.ndarray = jnp.zeros(self.nwalkers, dtype=self.dtype, device=self.device)\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Using 100 total chains: Group 1 (50), Group 2 (50)\n", "Using inital step size of 0.004999999888241291\n", "Starting warmup...\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "100%|██████████| 2000/2000 [00:40<00:00, 49.73it/s]\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Warmup complete.\n", "Found step size 0.010159287601709366 and integration length: 25\n", "Starting main sampling...\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "100%|██████████| 1000/1000 [00:19<00:00, 52.20it/s]\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Main sampling complete.\n", "Autocorrelation time: [41.01967 41.53285]\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAc4AAAHZCAYAAAAPNduTAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAOlBJREFUeJzt3Ql4VNX9xvHDDgFZAyhLQAHZBCSAiuKCe6GAWGTRUBCL1kpBqai4gLhA3UFa61oKjVWpiiKComURcWORTUIFBQICQpAdBIT7f97T/80z2ecmk8y9M9/P8+QBJmO8mWTmnXPO7/xOKcdxHAMAAMJSOry7AQAAghMAAI8YcQIA4AHBCQCABwQnAAAeEJwAAHhAcAIA4AHBCQCAB2VNAJw8edJs27bNnHLKKaZUqVLRvhzAV9TD5MCBA6ZevXqmdGneCwPFLRDBqdBs2LBhtC8D8LUtW7aYBg0aRPsygJgXiODUSNN9YahatWq0Lwfwlf3799s3lu7zBEDxCkRwutOzCk2CE8j/eQKgeLEgAgCABwQnAAAeEJwAAHhAcAIA4AHBCQCABwQnAAAeEJwAAHhAcAIA4AHBCQCABwQnAAAeEJwAAHhAcAIA4AHBCQCABwQnAACxdqwYSk56errJyMgo8H6JiYkmKSmpRK4JAPyE4ESW0GzZsqU5fPhwgY9KQkKCSUtLIzwBxB2CE5k00lRopqam2gDNiwIzJSXF3p9RJ4B4Q3AiB4VmcnIyjwwA5ILiIAAAPCA4AQDwgKlaFJrWOgtC9S2AWENwwjOFoapqVSBUEKpvAcQaghOeqZJWo82C9ntSfQsgFhGcKHR4shUFQDyiOAgAAA8ITgAAPCA4AQDwgOAEAMADghMAAA8ITgAAPCA4AQDwgOAEAMADghMAAA8ITgAAPCA4AQDwgOAEAMADghMAAA8ITgAAPCA4AQDwgOAEAMADDrKOE+np6SYjIyPf+6SlpZXY9QBAUBGccRKaLVu2NIcPHy7wvgkJCSYxMbFErgsAgojgjAMaaSo0U1NTbYDmR6GZlJRUYtcGAEFDcMYRhWZycnK0LwMAAo3iIAAAPCA4AQDwgOAEAMADghMAAA8ITgAAPCA4AQDwgOAEAMADghMAAA8ITgAAPCA4AQDwgJZ7KHbhnLpCj1wAQUFwotgoDHXaSkpKSoH31f0UsDSYB+B3BCeKjUJQYRjOOaAKV92P4ATgdwQnipWCkDAEEEsoDgIAwAOCEwAADwhOAAA8IDgBAPCA4AQAwAOCEwAADwhOAAA8IDgBAPCA4AQAwAOCEwAADwhOAAA8IDgBAPCA4AQAwAOCEwAADwhOAAA8IDgBAPCA4AQAwAOCEwAADwhOAAA8IDgBAPCA4AQAwAOCEwAADwhOAAA8IDgBAPCA4AQAwAOCEwAADwhOAAA8KOvlzvCf9PR0k5GRke990tLSSux6ACDWEZwBD82WLVuaw4cPF3jfhIQEk5iYWCLXBQCxjOAMMI00FZqpqak2QPOj0ExKSiqxawOAWEVwxgCFZnJycrQvAwDiAsVBAAB4QHACAOABwQkAgAcEJwAAHhCcAAB4QHACAOABwQkAgAcEJwAAHhCcAAB4QHACAOABwQkAgAf0qoVvhHP8Gc3qAUQbwYmoUxjq2LOUlJQC76v7KWA56QVAtBCciDqFoMIwnAO5Fa66H8EJIFoITviCgpAwBBAEFAcBAOABwQkAgAcEJwAAHhCcAAB4QHACAOABwQkAgAcEJwAABCcAAMWDEScAAB4QnAAAeEBwAgDgAcEJAIAHBCcAAB4QnAAAeEBwAgDgAcEJAIAHBCcAAB4QnAAAeEBwAgDgAcEJAIAHBCcAAB4QnAAAeEBwAgDgAcEJAIAHBCcAAB4QnAAAeFDWy51RctLT001GRka+90lLSyux6wEA/A/B6dPQbNmypTl8+HCB901ISDCJiYklcl0AAILTlzTSVGimpqbaAM2PQjMpKanErg0A4h0jTh9TaCYnJ0f7MgAAISgOAgDAA4ITAAAPCE4AADwgOAEA8IDgBADAA4ITAAAPCE4AADwgOAEA8IDgBADAA4ITAAAPaLmHwAnnVBh6+AIoLgQnAkNhqNNgUlJSCryv7qeApQE+gEgjOBEYCkGFYTjnlCpcdT+CE0CkEZwIFAUhYQggmigOAgDAA4ITAAAPCE4AADwgOAEA8IDgBADAA4ITAAAPCE4AADwgOAEA8IDgBADAA4ITAAAPCE4AADwgOAEA8IDgBADAA4ITAAAPCE4AADzgPM4Slp6eHtZBzAAAfyI4Szg0W7ZsaQ4fPlzgfRMSEkxiYmKJXBcAIHwEZwnSSFOhmZqaagM0PwrNpKSkErs2AEB4CM4oUGgmJydH438NACgiioMAAPCA4AQAwAOCEwAADwhOAAA8IDgBAPCAqlrErHAaSbDtB4BXBCdijsJQDSRSUlIKvK/up4BlzyyAcBGciDkKQYVhOK0NFa66H8EJIFwEJ2KSgpAwBFAcKA4CAMADghMAAA8ITgAAPCA4AQDwgOAEAMADghMAAA/YjhIh6enpYe0bBAAEG8EZodDU4dSHDx8Oq1ONOtsAAIKJ4IwAjTQVmqmpqTZA80NvVAAINoIzghSaycnJkfySAACfoTgIAAAPGHEi7nH8GAAvCE7ELY4fA1AYBCfiFsePASgMghNxjePHAHhFcBaAxgYAgFAEZz5obIBQFBEBIDgLQGMDCEVEAELF7YjTyxQsjQ3iG0VEAGI6OMMJxF27dplrr72W3rIoliKiSDbzp0Uj4D+BCE7HceyfixcvNpUrV87zfgrMlJQUc+TIkQK/ZqVKlcxbb71VYMP1WrVqmerVq5v9+/cX4soRTypUqGB/r/Q7GCn6euqBnN/v6aFDh7I8TwAUr1JOAJ5tW7duNQ0bNoz2ZQC+tmXLFtOgQYNoXwYQ8wIRnCdPnjTbtm0zp5xyiilVqlShvoZGjApfvbhUrVrVBBHfgz/47eegp/CBAwdMvXr1TOnStJ8Gilsgpmr1YhCpd9J6ofPDi11R8D34g59+DtWqVYv2JQBxg7enAAB4QHACAOBB6XiqeBw7dqz9M6j4HvwhFn4OAGK8OAgAAL+ImxEnAACRQHACABBr21EisY8TiFXh7uPkeQREZj90IIJToUnnIKBonYN4HgGR6cAViODUSFP80qkF8JMdO3aY5s2bZz5P8sLzCCi4I1hBz6PABKc7PeunTi2A3xS0jMHzCChYOMuBFAcBAOABwQkAgAcEJwAAHhCcAAB4QHACABDN4NywYYOZMWOGOXbsWKS/NID/d+jQIbNz5077J4AAB+eqVavM+eefb+bMmWMyMjIK/XWOHj1q99SEfgD4H4Xlpk2b7POC4AQCHJzp6emmR48eZvDgwebFF1+0bYuyC/cglgkTJtgT7d0PugYB/6OgXLt2rR1t7tmzx1SuXJmHBghqcGq0edZZZ5nHH3/cHD9+3Nx///2md+/eZujQoWbatGmZG0vDCc/Ro0ebffv2ZX6oYxAAYwNTz4ft27ebhIQEghOIgoh1Dlq+fLn56aef7N+7detmfvnlF9OuXTv77njp0qVm3bp1Zvz48WF1ZdABwRwSDOSkJZAff/zRzsQoOAEEODi1trlw4ULzyiuv2HBMTU019evXtyPGSZMm2XVPhWirVq0i9b8E4o7enOqUkzJlypg6depE+3KAuFToqdoTJ05k+be6yWtU+fTTT9vpWIWm6J3xjTfeaKdyV65cWfQrBuJ4fXPv3r3279WrV2eaFghScH777bdm4sSJdp3F1aJFC1sUpM8pJD///PPMz9WtW9ecd955pmbNmpG5aiAO1zbffPNN8/XXX9sagoKOPQLgo6la7dPs3LmzrejbvXu3GTlypElMTLSf6969u/nnP/9pbrjhBjNu3DhbYduxY0c7favRKNO0QOFs3rzZLF682B4hVrFiRdO4cWMeSiAIwampIm0V6dmzp+nUqZMZNmyYLQK66667MsOzf//+pnbt2uaBBx4wt99+u6lRo4Zdk5k1axbbSoBCUiGQTqVX0VyVKlWKNE27YsUK+zXyo+dzUlJSof8fQCzzFJx64nbo0MHUqlXL9OvXzz65FJQSGp6XXXaZOfvss20hg8JW00ru5wAULjibNm1qC4LatGlTpIfw4osvDuv/l5aWRngCRQ3OSpUqmUGDBmW+2+3bt68tBBowYID985577rGhqlHogQMHTLNmzbx8eQB50N5NPf9UJ3DmmWcW6XFSLYLeAOdFgZmSkmK3vjDqBCKwxumGpqpqNQLVyFOhef3119ttKJqeffLJJ+2ajBof6J1rOHs3AeRdGKSKdO3f1JvRom5Dad68uUlOTubhBkp6H6f2kSkwtX6p6VqF48CBA83MmTPNd999Z5YsWUK5PBABqqTViLN8+fLm9NNP53kFBLnlnsLSbaOnkeeFF15odu3aZbsIaY0TQNFHm9ripWlabevSti8AAe8cpODUtO2oUaPM/PnzbcVeUYsXAPyvil3buE499VRTtWpVc9FFF9EtCIilJu+tW7e2I822bdtG6ksCJt6DUwU669evt1O0+gAQI71qtd45ZMgQioCACNKyh3o863zaU045xY44AcTQiJPKWSCy1ClIo80ffviBhxaIxdNRAETOxo0bbb2AOm/poIRLL72UhxfwCYIT8CEdlKDALFeunG1GoBoCADE2VQsgcr755hvbwUdbUNQXGoB/EJyAD6dpNeLUNi/1ewbgLwQn4DNfffWVbSqiczfZEw34D2ucgM/s27fPNjxQT1m2oAD+w4gT8Nk07bZt2+y5mzqOr6gN3QFEHsEJ+GyaVqegqJpWfWmLcmA1gOJBcAI+auiuU4XU8EDduGjoDvgTwQn4hBq6q4q2YsWKplGjRkzTAj5FcRDgk4bu27dvtz1p1S2oXbt20b4kAHkgOAGfBKdCs2XLlvYUlA4dOkT7kgDkgeAEfGDz5s22MCghIcHUrFmToiDAxwhOwAfWrl1rp2rLli1rjhw5Eu3LAZAPioOAKFMl7ZdffmkPrdZos1atWtG+JAD5IDiBKK9tvv/++2bPnj32NJQrrrjCNG7cmJ8J4GMEJxDlvZtq5q5OQa1atbJFQTQ9APyNNU4gijQ9W69ePdOwYUM72iQ0Af9jxAlEeZp23rx5Zv/+/XYbCgD/Y8QJRLGS9tNPP7UBWr16dX4OQEAQnECULFy40E7VHjt2jNEmECBM1QJR3Iaiszc12uzSpQs/ByAgCE4gSqG5fv16U7p0adstiJNQgOAgOIEoSE1NtWub5cuXN7169eIkFCBACE6ghH344Ydm9uzZdpq2Tp06pl+/fvwMgAChOAgoYU888YTZunWrKVeunGndujWjTSBgGHECJdwpaPXq1eb48eOmVKlSdpoWQLAQnEAJ0Zrm448/bqdoT548aUeaVNMCwUNwAiUYnK+//ro5evSo/fcFF1xAiz0g3tc4N27caN555x27fnPOOecUuuhBLyzui4uoHRkQdJ9//rnZsWOH/bvO3RwyZEi0LwlANEecWre56KKLbO/NL774wlx//fW2CKIwJkyYYI9Ycj/UABsIujFjxtiTUKRJkyb2+QIgToNz8+bN5tprr7VhOXfuXLN48WLz0ksvmaeeespu8vZq9OjRdh3I/diyZUskLhOImk8++cSsWrUq898DBw7kpwHE61Stihy0btO0aVNz77332k4o0qlTJ1tur897pbMJ9QHEigcffDDz72qxN2DAgKheD4AojjgVlJ07dzZnn322nVZ1aX+a1nG2b99e1P8FEPjR5vz58zP/fdNNN9HUHYjHEafWasqUKWP/rrUad73GcRy7P030p/aruf7zn/+Ytm3bmtq1axf9yoGAGD9+fJZ/33XXXVG7FgBRGnF+++23ZuLEiVlGkwpMNyx/+eUXc+TIERusVatWtbdrGlcn3IcGKRAPzdzVYs9VpUoVOgUB8Tbi3LBhg52a3bNnj9m9e7cZOXKkSUxMzBxlutO3Ck2FqaZrH374YfPss8+aL7/80tSrVy/S3wPgW9lHlxQFAXEWnNrAra0iPXv2tMU/w4YNs6NLvTgoPEODs2LFina0eeutt5qVK1faStuOHTsWx/cA+JKeLwsWLMizSMjv0tLSCryPnvdJSUklcj1AIINTgdihQwdTq1Yt29xAT5r+/fvbz4WGp9Y/tY3k+++/NwcPHjRff/21adOmTfF8B4BPvfrqq1n+rZkatdnzOz2PdUZoSkpKgffV/RSwhCfiiafgrFSpkhk0aFBmm7C+ffva6ViV1uvPe+65x4aq/q5tKG+88YZp0KCBrbAF4s3f//73LP8eO3asCQKFoMIwIyMj3/vpPgpX3Y/gRDzxvMbphqZGlRqBauSpoFTzA61z3n777ebJJ580mzZtsof16h0pEI9FQVrTd6mS/KqrrjJBoSAkDIEIb0dxi380stR0rUJThQ8zZ860BURLly4lNBG3tLYfqlu3blG7FgA+aoCgsNSHAlQjzwsvvNDs2rXLrmmqIQIQr0VBy5Yty3IbnYKA2FHklnsKTk3bjho1ynZHWbFiBYVAiGt/+9vfsvxbU7RBmqYFUEKno6gAaPny5bYzEBCvdLTe/fffn+U2Nf8AEDsich6n1jt1tmBoEwQgXvvShp4lq4MOtIULQOyI2IiT0ES809pm9tFljx49MivRAcSGiAUnEO9mzZpltm3bluU2bc0CEFsITiBCo83hw4dnue3iiy/m+DAgBhGcQASsXbvW7Ny5M8ttTzzxBI8tEIMITqAYGh6oU5AOQgAQewhOIALt9bI3PJg8eTKPKxCjIrIdBYh17kHtucntFBF10gIQmxhxAkX07bffZvl39iIhALGF4ASKYNy4cTluGz9+PI8pEMMITqCQDh8+nCM41XqShgdAbCM4gUJ69dVXc9w2b948Hk8gxhGcQCHdcsstOdpOahsKgNhGcAKFMHfu3By3jRkzhscSiAMEJ1AIV199dY7bdCYtgNhHcAIeLVq0KMdtV155pUlISOCxBOIAwQl49Otf/zrHbW+//TaPIxAn6BwEeLB06VJz8ODBLOfPtmjRgtEmEEcYcQIedO3aNcdtzz//PI8hEEcITsDDaFNND0JVr17dXHjhhTyGQBwhOIEw/epXv8pxW/ZTUQDEPoITCMOmTZvMTz/9lOW2ihUrmsaNG/P4AXGG4ATC0KFDhxy3vfbaazx2QBwiOIECfPLJJ2bv3r05br/88st57IA4RHACBejWrVuuo00aHgDxieAE8vHNN9+YQ4cO5bi9b9++PG5AnCI4gTwoMHPbt/nII4/wmAFxjM5BiGuO4+T5ueXLl5s9e/aYsmXL2k5Bv/zyi7199OjRJXiFAPyG4ATyMGjQoCwBqw9V14a22wMQfyIenFu3bjWfffaZfZferFkz06ZNm0j/L4AS2be5ZcuWHLfPnj2bRx+IcxENztWrV5sePXqY2rVr2xedc845xzzzzDOmSZMmnr7O0aNH7Ydr//79kbxMoEC5tdHr0qWLqVOnDo8eEOciVhy0efNm25JswIABZsGCBWbKlClmyZIlZvfu3Z6/1oQJE0y1atUyPxo2bBipywQKtHbtWvPjjz/muP2tt97i0QMQueD88MMP7dTs+PHjTeXKlW2IJicnmxUrVphp06aZ+fPnh/21VHyxb9++zI/cpsyA4tK/f/8ctyUmJjLaBBDZqVoVTqSnp9ugbN++vXn00UfNnDlzzLFjx2z4aUT62GOPmcGDBxf4tSpUqGA/gJK2a9cus27duhy3b9iwgR8GgMiOOK+88kpz6qmn2o3hffr0MQ888ICZMWOGmTt3rpk1a5Z9Fz916lQ7dZvfFgAgmpKSknLcpjV7ugQBiPiI8/TTTzepqal2XVNrRCrZ79Wrl/2cCirq1atnFi5caKdxKeeHX0ebJ06cyHH7Rx99FJXrARAHnYMUnhpxNmjQwBw5csRO07pUbKEjmHJ7YQL8oH79+rn+Trdq1Soq1wMgjhognH/++ebOO+80kyZNstO3a9assVW2OmVCI06guIW7HHDy5Ek7A5KWlmbKly+fa9Gb+7WYKQFQbMGpd+ha3xw6dKgpXbq0fSevaVqaIcBvFIb6UAV4djVq1OCgagAl13JPzbG/+uorc/z4cVshW7169eL6XwGFpiWFv/zlL7l+TkVtAFCivWpr1qxZnF8eiMgJKGPGjMlxu97sqS8tAGTHsWKIazNnzsz19vfff7/ErwVAMHA6CuK6kfttt92W4/a6devavrQIjwqrCqLOS7ntkQWCiOBE3MqttZ6okA0mrDBUY4iUlJQC76v7KWAJT8QCghNxO9pU8Vr2LSj6t/Ybo2AKQYVhRkZGvvfTfRSuuh/BiVhAcCIunXvuubnezgko3igICUPEG4qDEHemT59u2+tlpwC44ooronJNAIKDESfirnPQzTffbKpWrZo5NeselK4DCQCgIAQn4kr2Zgfqp6yP1q1b2760AFAQpmoRNw4fPmyPu8vNvHnzSvx6AAQTwYm4MWzYsFxv155NHX0HAOEgOBEX1q1bZw8eyM20adNK/HoABBfBibgwYMCAXG/X1C1rmwC8IDgR8z777DPb8CA3I0eOLPHrARBsBCdiXvfu3XO9Pa9CIQDID8GJmPbiiy/m+TlGmwAKg+BETG8/ufvuu3P93MqVK0v8egDEBhogIGY6AmW/n/rRuh2CQqlbEP1VARQWwYmY9N5775mtW7fm+rn169ebUqVKlfg1AYgNTNUi5hw6dCjXA6qlVatWpnbt2iV+TQBiB8GJmJNXaLpbUwCgKAhOxBQdF5ZXh6C8CoUAwAuCEzGladOmeX5u9OjRJXotAGITwYmYMW7cuDw/x+knACKF4ETMePrpp3O9PTk52X4AQCQQnIgJ+TVqZ7QJIJIITgTek08+aX766adcP9enT58Svx4AsY0GCAiUgwcPZvn30aNHzRtvvGHOPPPMLLdXrlzZrFixwrz88sslfIUAYh3BiUDJ3vEnJSXF9qTNTrdpzyYdggBEGlO1CKxXX33VpKen5/q5hIQE07p16xK/JgCxj+BEIO3Zs8dMmjQpz8/n1QQBAIqK4EQgXXPNNXl+7t577zU1atQo0esBED8ITgTOxx9/bI4cOVKoUAWAoiI4ESjbtm2zI8q8vPPOOyV6PQDiT7FU1W7fvt2uQekIJyCS8htNqsK2Xr16POA+lZaWVuB9EhMTOWQc8RecP/zwg2nXrp256KKL7MigY8eOkf5fIE5dcskl+X5++PDhJXYtCJ/CUFXOemNTEN1PAZuUlMRDjPgJzvXr15t9+/bZj8mTJ5sRI0Zk9gl1HCesfXXa1K4P1/79+zOPjKpatWqkLxkB8OKLL5qFCxeali1b5vr5Tz75pMSvCeFRCCoMMzIy8r2f7qNw1f0ITsRVcLZt29Z069bNdO/e3bzwwgu28baOc9KeunCDc8KECbmedKEnVJMmTSJ9yfCBLVu2ZPn3zz//bPbu3WtOPfVUs3PnTjNlyhRz3nnnmTPOOMOsWbMmy3179+5tKlasWMJXDC8UhIQhYkVEg/PEiRP2Y926dea5554ztWvXtiGo/XbffPONOe2008ybb75Z4NdR0I4cOTLLiLNhw4Z2ygfxEaL9+/e30/4PPvigmT17dpbP6fcr1KhRo+wUHwAELjhLly5tw7JTp052VKCRQIUKFcygQYPs1OvQoUPD+jr6b/SRnb42YpuCsV+/fmbr1q3232PGjDENGjQw9evXz3W2YvHixYQmgOBuR3Ff2MqUKWMWLFhg//7222/bUahGjIsWLTJfffVVkf4fhw4dMhs3brQf+jtih8JSI039qWPCfve732XevnnzZjvVH2rIkCGmQ4cOUbpaAPEqoiNOdw3z0ksvtcH2hz/8wU6zLVu2zJ5UoSm18uXL23XQwq5JKSwPHDhgm3jro3HjxvYkDASbwlEjTY04FZo68UTrm3rDNXbsWLNjxw5z/PjxLMVBua2DA0CggtMdceqF78YbbzR169Y1s2bNsv/Whz6vrSpFKeRQSJ5yyil26lcj2U2bNhGeAffFF1+YXr162appvRFyQ1P0e1SzZk271WT37t12f7BonZOpewAx0zmoc+fO9hzEDz/80E6luVNs2ryuAC0KBae+hporKDw1+lTVJYLptddes/szFZotWrTIEpounbWp36GyZcvan39qaqoNWACImc5B5cqVM4MHD7bFQlIcZyJqylbt1/TnsWPHTJ06dZiyDZCTJ0/aitmHH37Y/vvyyy83zz77rKlSpUqO+06dOtX+qT28CtU+ffqU+PUCQLEfZO2GZnHRCEX/D62NqdlCpUqVbDUv/E9r06q0drcmae1b6+EqKstOezndI8IUmn/5y19K/HoBoESCs7hpfUsvwOqLq2m8r7/+2t7eqFEjO/qEP2lvptYzVTCmmQk1ydA6ZvYGCC6Fq05C0T5NhSsARFtgg1PhqA81RdAWFxUJabuLpm817UelrT+sXbs2S2hqCl8/I52XqcYYWgPXfVT0k/0Nj6Zz//GPf9i/d+3a1Vx55ZW2oxAARFNgg9PlFhtp9Lly5crMjfOEp7+EhqYKe55//nm71cSl9Ut1lgql4rL09HQ7Jf/Xv/7V1KpVq9jWzAEgrs7j1EjlnHPOMdWrV7drnzNnzrRTgfBnaKrvbGho5tfYXTRF64YmAERb4EecomlZNZHXaFMnaKjLjHrj3nHHHebXv/4107Y+C81w1qC/++47M3fuXDu6dDsIAYAfxMSIM3S/n5ojqBLzv//9r/n73/9u1z4RHQq/woSmvPTSS/bPK664ghNxAPhKTAWnXpR1gLY6E6kSU0VDOkybsxqj09hA57AWJjQ13a77y80331zMVwoAcRycmrL9zW9+Y7p06WJ74mp/p6b7dEzZkiVLon15cUGP+cCBA831119vj4M7++yzPYWmTJw40RZ7tW/f3lbSAoCfxFRwil6gVVGrCk3t79T2BYXmH//4R/Puu+9yokox0hFfCkq1xFMzAzVhV9cfL6Gp0aY7TavZAipoAfhNTBQHZR91qieupghVlakGCTpV48svv7QvxMJWlcjSG5THHnvM3HfffXbvpbYIvfrqq7Znceg+znA888wzdpq9Y8eO5qqrrorwlSII0tLSCryP9m8nJSWVyPUAMR+cohGOtjBoulab7N0m8HoR1wuzRqFU20YuNNVz9qGHHrL/Vis99ZzVvkyv9HPS4QCi6XVGm/FFYagOUSkpKQXeV/dTwBKeiIaYDE43PEeMGGEOHjxoR0MaCYm2q6hLjaZy1bVGT8CC8AKe01tvvWVDc/r06ebf//53Zmj26NHDfPTRR5n30zpnOHs2NbWrGQJ3tKlqWsQXhaDCMCMjI9/76T4KV92P4EQ0xGxwutO2Q4cOtR1oli9fnnn7qlWrzLBhw2zD8AsvvDCq1xhU2UNTxUAKPPcwc5f60WbvCJQbvQi+8sor9u86MUXbihB/FISEIfwu5oqDstN6m6ZnL7vssiy3r1692txwww22v60qOFH40FQl83/+8x9bgHXLLbfY9ngqFNJZqeHSnluNNs877zzWNgH4WswHp2hvpyo1L7300iy3q9NQ//79bSEL4eltTTM0NPXm48cff7T//umnn8z8+fPtm5UhQ4bY+2pkr9NrfvnllzwraXWAtej+TI0D8LOYnqrNPvJU4ckll1xiG4e7VHV76623mhUrVpgnnngirDXPeHXixAkzduxY8+ijj9p/9+7d24aminp0VqaqahWgaravx1NHhal7kD7U1F0FQ+eee6751a9+Zdcw3YDUFO3Ro0dtFS77NgH4XdwEpxueKkDp27evLVpxqXDob3/7m2nZsqUdJRGeOc2bN8+MHDnShqKcddZZZtasWXarj0JTI0VVRWo9U3s5VSi0e/du271JI/vPP//cPuYqHNKH1p3r169vPvjgA9tlSLTvk9EmAL+Li6naUNobePfdd+f6ueHDh9sA1dQh/kejRR08rTVihWalSpVs4c6aNWtsaCpAFXgKzex0oolG+E899ZQ9V1OtEEXrmApYhaZ+FhrJau+t9tcCgN/F1YjTpdM2tE1lwoQJOT43atQou+apkamqROOZCny031VN83Umpg6f1ihSFHy//e1v7WNU0ChRJ6T8/ve/t/tn27RpY/fWqhViaGgy2gQQFHE34nT3eGo9Tut1udH6XJ8+feyLe7wWDWkaViNAhaYCU9PZCk2dearCn6efftp06tSpwNDUf68tQVoHbdKkiR3Ra100e2hqHycABEFcBqdoHVPBqfXO3KiA6LrrrrOjz0OHDpl4Mm3aNBtoGiGWLVvWNoxQsGm7yYYNG8ztt99u92cWRI/b448/bh9LrWdqFK/Wh4QmgCCLy6naUGPGjLHbJN5+++0cn9M+RDVKWLdund2Un1/RUPaN//mJRgGM2zmpINlHfu4WEo0OJ0+ebD9kwIABdptPXlQlq9Z7qqzV+qe2A+mMVEITQNDFbXC64dW6dWs7wqpbt66dRsyN9iEqWO+//34bGLFYdavgV3/YcGnd8pxzzsn1cyoaUgXu999/bwuJ1BxBx43pNgWxugzpMWd6FkVBM3hES9wGZ/bWfNpOoe0od9xxhx1dZafTVtQ4Xlsq1B0nlsJTYabvSZWvRaXHTqP4zz77zFbRqhBIt6nISFO/3bt3t/8fQhOFRTN4RBvBGVIw1K9fP9snU6348lrXvOeee8ycOXNsRa6axAed2txpFP3ee+/ZUbhGnkU9XkxFVVobVUMJFRNpC5DecGhqVx2HwlkfBfJCM3hEG8GZLTy1Z1EFQfmFp05YOf/8822zBBXMtGrVygSRpk979uxpPv30UzsCzG2k7YWmtHVqigL4kUcesVW0eoy0L7Z9+/Zm5syZdh8oUFQ0g0c0xW1VbX7U+k2Nyi+44IICG5Nr1KmpSe15DNLWFYWk+vQqNHVuqf5d2OlTFVGpetad6tWB4e3atbPTvzt27DBnnnmmbXZQrVq1CH8XAFDyCM48Rp7ax/nPf/4z38pRt1pV05MaqeZVXOQHx44ds8UUKnLSNLNGmppSVWjqc6KCHS9UBPT666/brStuk3ZtVenSpYu5+eabbUVtvXr1bIs9PaYAEAuYqs2nYEi9bd99913zwAMP2GnI/Gha96677jITJ060fVpVcaqRa+3atU1J0qhX7fAUktpG4/6p1nm5nU6ijkBy55132hGopqjDoelXvblQH1pp1KiRDc0WLVpkCU3t3eR8RQCxhOAsgKpntYdT65ja01nQfkhV32pEp7W85557zralK076/2maWFWs+jO/47uqVKlig6158+b2TcHUqVNtwHXt2jXzxBOvWwHUj1ZTshpx60ix7KEZziHWABAkBGeY4XnTTTfZEaTO9FQnnXCqVW+88UYbtgpPbWUpahGR1iE1mnRDUh+bNm3KcT/tSdX/SyGpLTb6UwHmHvm1fPlyWwClgGvQoIH517/+ZatgvdBIVa30UlJS7OOjtUz1oyU04Tfs90SkEZweKIxWr15tpyTffPPNLJ9T4YvW/HKjYNJHzZo1zeDBg+32D7Wgy48C0d22oUYCOs9Sa4VqTh9KFawKvzPOOMNWsepD/50CTB9qcac/czvxRfdTk/a8evbmR1PRV199tV0z/eSTT2x/XwU7I034Bfs9UVxKOUXZuFdCtAdQwaTtEzoMOdrUsFx7OdUJR9OTEm7FqIJN4afmACo80qHOqkBVqGYPTgWepnu1x9KdItbaq867VDVvcnKyvY/WUdVcYMmSJTbEdFtu9N+qYbs+9P/TC0uFChWy3EfFQqoWlvnz59u9rRkZGfbfKijSuuaqVatsE3itnYZS03ft09R6J/z3/HDvp+1UBRW9xQr1SXZ/f/MbkWrmZNmyZfY5hfi030POMOIsBFWIqgBIoaeOQ4WpplXQqapVH6Iga9u2rTn33HPtCFLFSO+//35mYOo8TK0fKmRDt41oJPrxxx/bkaW+Zuh0bcOGDe2HinM0RavRYDj0XkpFTip20ihSQah1UW1d0V5Ml6Z3L774YtOjRw/7oVEv4Cfs90RxIDiLGKAaEWotU8UxRdnHqZGrQvS1116zFbqhgakmC+qpG9pAXffV9O/SpUszb9fIUyMJrcWqEMjrOo+osGjgwIH2a7uj1M2bN9sP0Wi1W7duNig1VcveTADxpliCUy/6GrXESz9STVFqje/Pf/5z5jSnFworrUGGBqYqXUeMGJElMPV5jTCnT5+eOUWsIh2NVBWYKgJyt5cUhoJf03gqftLX0c9P/0+tx2qrisJSHZPcQqJonPICFBeKiBC14Fy7dq0ZP368rbJs1qyZHb3oxTbWab3woYcesvsgNRrUVoxwqm9FoanuO6Im6Oq8o5Fd9p6u6syj4HSnYrX+qLl4Te1GggqfdM0auWodVpXBmo5OTU211wPEIoqIENXg1HmLCkm92GoUpgIaTSUqPIcPHx7219FUpD5CF22DQNsy1JNVHzp3Uq34CmqcIO4oU284VLEruW0zcStqNS2s0a1Gfu+8807Ert99zNXI/k9/+pP9u6pmsxcQAfHcNH7RokV2m1dBYUzjj9gVseDU1KxatukkDK3TiUZOOsx4ypQptnBFxSbhUAOBcePGmaDTCFRvGGbMmGEfg+xVqNmFu16oJ6XXfZdehE7BFmXqF4ilIiJGpnCVjeSLrbrYaIrWpUOMFRzaeqGeplorC6elmw5U1laP0BGnqkODSE82NQrQh0rj1f9WfV3z2vMJwJ8YmSKiwanRpoJTe6DWr19vp2y1fcENTx0tpdtUgdq7d+8CD4HW1GAsTg/qiacOQno8tm/fbvdc6hiuH374IdqXBiBKI1Mth5R0T+twMeVcjMHpTu1pm4KmJ3XE1KRJk+yWCIWqCkvUKF37ARUW2sbghdujwc9rnQX1sHXpJBJNs2oErdG3PtSSTwVFKsZxv0etZ2afjnVPMdGfbjGRpsDDmbZV5W5Bo1z3ewjdD6rryevgaapq/cH9nSmol4n7eVVK+/m5FHQ6vP2rr74yu3fvzvd+WlNVuHp9PSxJ6rmt4kAFaKw79P/nL4fVE8iJsHnz5jkVKlRwbrvtNmfXrl2Zt2/fvt1p166d89lnn3n+mlu2bNF3wgePAb8D+fwO6HnC84jfEV4nTLE+jyTiFSbaf6i2a9ddd52djuzbt6/dZ6jCIbWqK8xapTreqI2cpn0LO8px10n1dfzQtq8w+B78wW8/B71D1gxEQZ2hsj+P/PZ9hCuo1x3kaw/qdXu59nCfR1IspZnaKK8TPFTgo20ZmkrUZnq1kCvMnkNVdkZqr6IeuKD94LPje/AHP/0cwqnIzut55Kfvw4ugXneQrz2o1x3utYe7s6HY9jSoUEh9TdXhRimuY63iYZ4cABDbirVXbZDfnQAAkJu42d2u7S06dzLI21z4HvwhFn4OQf4+gnrdQb72oF53cV17IM7jBADAL+JmxAkAQCQQnAAAeEBwAgDgAcEJAIAHBCcAAB4QnCFNfYNaYMx1wws18z9x4kTgHjS18Fy7dm20LwMgONesWWPuvPNOe3pIEE772Lp1q/nwww9tP+DNmzfb23Td4Z7O4gdpaWn2IwiPd6xR8Og0Hh04f+utt9rWmEGgo/fatGlj7r//frN06VITFHq+Tp8+3R4dtnr16mhfTlzasGGDmTFjRubpUpEQ1yPOlStXmg4dOtj+hO7RXBq9+XUEpydex44d7RFtAwYMMH369LEHhbt9SIMQnqtWrTKtW7c2s2bNMkGls2X1Aq6fwZQpU8yyZctMUK77/PPPt6PNTp06mc8//9yMGDHCPPvss8bvdM7vvn377MfkyZPN8uXLMz/n5+drly5dzBNPPGHP4b3vvvvMd999Z4Jm48aN5plnnjF/+tOfzBtvvGGCRK83+p2fM2eOPcYtYpw4tXLlSqdy5crOnXfe6QTB3r177bFst99+u/371q1bnYcfftg566yznO7du2fe78SJE45frVixwqlUqZJz9913O0H1zTffODVq1HB69erlXH755U7r1q2ds88+25k2bZrjZydPnnTuvfdep2/fvpm37d+/33nkkUfs9T/22GOOn+3evdvp2bOn88ILLzjJycnODTfc4KxZs8a3v/ObNm1y6tev79xzzz3OwYMHndmzZzunnnqq8+WXXzpBsmrVKqdBgwbOZZdd5px//vlO6dKlnccff9wJgs2bNztJSUnOqFGj8n1eFEZcBqce0Fq1ajn9+vWz//7ll1+ccePGOQMHDnSuvvpq5/3333d++uknx2/XfOaZZ2Y5z/TAgQPO9OnTnebNmzvXXXed42fffvutU6pUKeehhx7KfMz//e9/23/re/j6668dv9M133jjjc6gQYMyn3BLlixxhg8f7tSsWdN5+eWXHT8bPHiwc9FFF2W5TeH55JNPOh07dnRSU1Mdvz7uO3futL//esP49ttvO506dXKGDh1qX8x/85vfOH6jgL/kkkuyvDB369bN3j516lR7brHfKfybNm3q3HXXXZlvTl555RWnbt269vnsd++99559zOXYsWPOfffd51xzzTXO7373O/szKEp4FmuTd7/S+prOZ9P0pk5q1xSK1jjr1q1rp2w1BTd69Gg7DZqQkGD8QGcoHj9+3K5Jde7c2d5WpUoV07NnT3PkyBHz1FNPmRdeeMHccsstxm/0Bu3TTz+1f2/WrJn98/LLLzd79+41Bw8etJ+vUaOGnf7s1auX8Stdp9ZL2rVrl7k+q6nz2rVrm/Lly5sHH3zQ/l0/E79dt65XJxZpylNTts2bN8/8vRoyZIi97bnnnjO9e/f2ze+8S89TPa6aXlZNgq5RfUcHDRpkjh49aoYOHWr8Ro95enq6WbFihWnfvr159NFH7XSh1tk03az6hMcee8wMHjzY+JGWfV5//XXTtGlTc++999qfgehnUK5cuUAsC2k6X6dzSbdu3exrvJ67WufXOvm6devM+PHjC1dr4cSpGTNm2Kk2TbtdddVV9h2t3tmKpkATEhLstJxf/Pzzz3akoxGxpk9CHTp0yE5j9e/f3/ErjY41stGoU1NYGiX897//tZ/T9NWAAQOcrl27Ojt27HD8TNM++n3Ztm1bltv1vWgGQyN//Tz8aMOGDU5iYqIzZMgQ+/MIfbednp5ufzZz5sxx/Oq3v/2tnfqUm266yT53W7VqZb8fv02Bfv/993Y0rBGbftf12L7zzjv28f7xxx/tLIVGpBkZGYWeLixuCxcuzHy8XRp5Nm7c2Jk/f77jdx999JFz6aWX2pmgK664ws5WiJa6NMN43nnnFfo1Pu6CM3Q95K233rJhtHjxYvvv0F9gvcBMnDjR8ZPVq1fbaRKtU+lFMNRTTz1l1378+qItR44csdd54YUXOkuXLs3xRqZixYp27dnPNL18+umnO5MnT84MH9e//vUvu26+ceNGx680RVihQgXntttuc3bt2pV5+/bt2+0aeuhSgF+4z8t//OMfztixY51bb73VOe2002w4adq2SZMmzu9//3v7++Unur433njDXnOfPn2yfO7Pf/6zfbz9ds3u4CGvn4FeP08//XRn7ty5mZ/7+OOP7cDDb9eelpbm1KtXz7650iAplN4oanCk52xhxF1wZg9IPbgazYXSC59+qUN/OfwS+F988YV9cdYTMXSdROs9GnUePXrU8bN9+/bZIiH3Ot3vSy/Y+gX3c+i4tOajJ91LL71ki1Zc69ats9+D/vSzmTNn2vC89tprnddff91Zu3atHVkojLZs2eL4lUZAGrmpyCb0jZfedCmk/Eq/JyrgC31u3nHHHbbATIVDfqFZE80Khc6mhL5WHj9+3F5v06ZN7euQjB492v5MfvjhB8dv1y6zZs1yypYt69SpUyfLm0L9LDQa/eCDDwr1/4vL4JT8pkfuv/9+W63qDu1LkoIk+zsnN1zc2/WioUpIjTAV8HoCVq1a1QZSUGkKtHPnzs6ePXucIMxW/PGPf7QFQapU/eqrr2yAqkJbo5/QkZxfLVu2zLn44oudRo0a2WtW4c3y5csdP1OBh4pT3FkJv05xZqfpwGrVqtlqVFVf641X9erVcyy5RNP69evt77NCUGGY2++wfv81Qm7SpIl9DVJhn97E6/ffz9f+2muv2WpgLbHo77q/3ihqNKqRZ2HEfHDmNfWQG037DBs2zP6SR6PKU08wldmr9FtTT3q3lP37cP9Ula17vdpKoJFz0B5v0ehM7761XuWXaVqFoB5PVQ5mH8GHfn+abrvgggvs6K19+/Z2xOb38Mk++tcIXy/gQQh7v249CYdmhxQ4zZo1s2ubfvldF40itU6squu//vWvNoD0Rjav34n27dvbquby5cvbqvJoCvfaNZ2sN+Za6mrRokWR3yjGbHCGFpkUNG/vevrpp+2DG413ggoQBbYKfPRuSCNJbRHQvk2X+yLu13faeU2XhAq9dq3Z3nLLLfaJ6JfRsq5J19OmTRsbiCoUy/77oykrl97ALFiwwE4jRmOGAsGhN2R6XfLbrMrhw4dt6GjaXrQum1sA6Xmg76FatWpOmTJlfDFiDvfaRYVYejOsQVFR3yjGZHBqtFCuXDmnR48eYY2EQtepovFLXdDmdK1fhlJ1nirz/CScqZ7cAl9ThipM8QON+LW/V1Ou+rtbBRw6nRPUEQ+Qn+xrrQoi/e7ruaDAcd8w6nn9wQcfZDafCNK1R7J+IuaCU6MdlYGr1PiMM86wBRD5hefIkSNtuXi0iwsK2pw+YcIEe5umb9XJQ5t5/fIi7nWqR2s9Dz74oOMnulY9/iNGjMgS9Nr+o6ICvUsNLZyZNGmSM2XKlChdLVA89BrpvsHVeqD7XFbxj5ZUevfu7dvK/YKuXVmg16pIzNjFXHCqvFgPkKbP3nzzTbvnKDQ8Q6fZRPfRSClaVWHuD/HZZ5+162XZKzLVwcjtkOJO1Y4ZM8b57rvvHL/wMl2i0b32O5577rmZ7wb9QNcyfvz4LB1RVPyg70Ojfr1ZUXHBokWL7PegNR6FqtYJgVii1yT3Tbme05q9U3cyTc/6vcPXyTyuXZW1kbz2mAtOvRt69913M0NS7dyyh6f7zsR9gLPvx/Pr5nS1kPKrcKZL9LhrKlzBk986aLRohO9y37HqTYCuV2uYCkvtyROt72h9E4hFeu1xX3+0bUODCz+safrl2mMqOHObulT5tDatZw/PF198MXN055dimyBuTvcyXaJCJ/WKzL5v1q99OrX+Gkp78fThl98XoLify5ri1HPYT1XAfrj2mOpV6/ZTDFWxYkXTvXt3249Q527qKK4GDRrYo5TcI378ci5k165d7Tmb1113nT20t2/fvqZt27Zm2rRpZufOnba/rt+VKVPG9ulUL8v+/fvbx3bgwIFm5syZ9vFWb2D1GfW7Ro0a2Q/R96Ieo+oNrJ+HX35fgOKmIwDV81W/90HTuhivvZTS08QQNfJ1z9YMbXD9888/m3fffdc2cFdD8Y8++sg2vfYj/bBHjhxpNm3aZL8XhZEaLqtZdFC4v1Z67C+77DLb7HrBggX2MOIgGjNmjJk6dar5+OOPMxvVA7HOff0MIqcYrz2mRpw6oFdBo8BZuHChPT3BfeA08tSLnk5+WLRokWnVqpXxKwW6Rmjq7H/gwAFz2mmnmcTERBMketz18xg1apSZP3++Dc4ghqZmAPS7pDcuerNFaCKeBDU0i/vac85tBnikqZGZQlNHJs2bNy/L52fPnm1fAPXh59B0Va1a1TRu3NiGTdBCM1amekS/K7t27bJvtoI04gdQfGJiqtadnlVoarSm8/p0NmXolK3OrNQ5eKeeempUrzXeBHmqx6VzUHUGIQDERHBmD00dIvzyyy9nCU0Vd+RWOAQAQFwFp9bQ3OnZvEITAIBICvQwTKG5efNmu452zTXXmFdeeYXQBAAUq8CPOG+++Wa7hvb8888TmgCAYhfo4JQ9e/aYatWqsYYJACgRgQ9OAABKUqDXOAEAKGkEJwAAHhCcAAB4QHACAOABwQkAgAcEJwAAHhCcAAB4QHACAOABwQkAgAcEJwAAHhCcAAB4QHACAOABwQkAgAnf/wG0O5BBTcyyDAAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sampler = hemcee.HamiltonianEnsembleSampler(total_chains, \n", " dim, \n", " log_prob, \n", " step_size = 0.01,\n", " adapt_inital_step_size=True, # Default to `True`, uses default parameters\n", " adapt_step_size=True, # Default to `True`\n", " adapt_length=True, # Default to `True`\n", " )\n", "samples = sampler.run_mcmc(keys[1], \n", " initial_states, \n", " num_samples=10**4, \n", " warmup=10**5,\n", " thin_by=5,\n", " show_progress=True)\n", "\n", "try:\n", " tau = hemcee.autocorr.integrated_time(samples)\n", " print(f'Autocorrelation time: {tau}')\n", "except Exception as e:\n", " print(e)\n", "\n", "corner.corner(samples.reshape(-1, dim).__array__())\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "02ce0500", "metadata": {}, "source": [ "## Modifying the Adapters\n", "\n", "Incases where the adapter behavior fails, you can modify it for your purposes.\n", "\n", "### Dual Averaging\n", "\n", "Dual averaging implementes an optimization scheme to adjust the step-size to reach a target acceptance ratio\n", "\n", "- `target_accept`: Target acceptance rate.\n", "- `stepsize_inter`:\n", "- `t0`: Controlls how quickly adaptation begins. Large $t_0$ means adaptation is slower early on, and smaller means adaptation is more reactive.\n", "- `gamma`: Controls how strongly the running dual average is pulled toward the instantaneous noisy update\n", "- `kappa`: Controls the weight given to the running average vs. most recent log‐stepsize. In otherwords, how quickly the optimization scheme \"trust\" new information. Large $\\kappa$ means a more timid algorithm, and vice-versa.\n", "- `agg`: Method to aggregate the acceptance probabilites. Do we do a harmonic average, or a regular mean?\n", "\n", "Since the the scheme is affine-invariant, the acceptance ratios are affine invariant. Thus dual averaging itself is also affine invariant." ] }, { "cell_type": "code", "execution_count": 7, "id": "ffa69753", "metadata": {}, "outputs": [], "source": [ "from hemcee.adaptation import DAParameters\n", "\n", "# These are the default / recommended parameters.\n", "da_parameters = DAParameters(\n", " target_accept = 0.651, # Recommended for ChEES, 0.8 is recommend for NUTS\n", " stepsize_inter = 0.9,\n", " t0 = 10.0, # How quickly adaptation begins\n", " gamma = 0.05,\n", " kappa = 0.75, #Controls the weight given to the running average vs. most recent log‐stepsize\n", " agg = 'harmonic' # How to average 'mean' or 'harmonic'\n", ")" ] }, { "cell_type": "markdown", "id": "7aac2f36", "metadata": {}, "source": [ "### ChEES\n", "\n", "The ChEES adapter maximizes the “Change in the Estimator of the Expected Square\"\n", "$$\n", "\\text{ChEES}(T) \\equiv \\frac{1}{4} \\mathbb E\\left[ \\big( ( \\tilde \\theta(T) - \\mathbb E[\\theta])^2 - (\\theta - \\mathbb E[\\theta] )^2 \\big)^2\\right]\n", "$$\n", "where $\\theta \\in \\mathbb R^d$ are the current positions of your walkers, and $\\tilde \\theta(T) \\in \\mathbb R^d$ is the position of your walkers after being time-evolved (via Leapfrog integrator) for $T (= \\epsilon L)$ total time (step size * number of steps). To maximize this quantity, we use the ADAM optimizer.\n", "\n", "The current metric is invariant under orthogonal transformations, but we are interested in being invariant to all affine invariant transformations. If you view the vector norm as an inner product, then you can ask, what's the affine invariant inner product. You get\n", "$$\n", "\\text{ChEES}_{\\Sigma^{-1}}(T) \\equiv \\frac{1}{4} \\mathbb E\\left[ \\big( ( \\tilde \\theta(T) - \\mathbb E[\\theta])^2_{\\Sigma^{-1}} - (\\theta - \\mathbb E[\\theta] )^2_{\\Sigma^{-1}} \\big)^2\\right]\n", "$$\n", "where $(\\vec v)^2_{\\Sigma^{-1}} = \\vec v^T \\Sigma^{-1} \\vec v$, and $\\Sigma^{-1}$ is the empirical covariance of the complement ensemble.\n", "\n", "Criterion parameters\n", "- `T_min`: Minimum allowed integration time.\n", "- `T_max`: Maximum allowed integration time.\n", "- `T_interpolation`: \n", "\n", "ADAM optimizer parameters\n", "- `lr_T`: Learning rate of ADAM optimizer.\n", "- `beta1`: 1st-moment smoothing (beta1 in ADAM), this is the momentum term.\n", "- `beta2`: 2nd-moment smoothing (beta2 in ADAM).\n", "- `regularization`: Avoids division by zero in the ADAM algorithm.\n" ] }, { "cell_type": "code", "execution_count": 8, "id": "b709a435", "metadata": {}, "outputs": [], "source": [ "from hemcee.adaptation import ChEESParameters\n", "\n", "# These are the default / recommended parameters.\n", "chees_parameters = ChEESParameters(\n", " T_min = 0.25, # Minimum allowed integration time\n", " T_max = 10.0, # Maximum allowed integration time\n", " T_interpolation = 0.9, \n", "\n", " # ADAM optimizer parameters\n", " lr_T = 0.025, # Learning rate of optimizer\n", " beta1 = 0.0, # 1st-moment smoothing (beta1 in ADAM)\n", " beta2 = 0.95, # 2nd-moment smoothing (beta2 in ADAM)\n", " regularization = 1e-7 # Avoids division by zero\n", ") " ] }, { "cell_type": "markdown", "id": "ce88bdfd", "metadata": {}, "source": [ "We can adjust the behavior of our Hamiltonian sampler by this into the `adapt` arguments" ] }, { "cell_type": "code", "execution_count": 10, "id": "c9998af6", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Using 100 total chains: Group 1 (50), Group 2 (50)\n", "Using inital step size of 0.004999999888241291\n", "Starting warmup...\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "100%|██████████| 2000/2000 [00:38<00:00, 52.01it/s]\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Warmup complete.\n", "Found step size 0.010159287601709366 and integration length: 25\n", "Starting main sampling...\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "100%|██████████| 1000/1000 [00:18<00:00, 54.98it/s]\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Main sampling complete.\n", "Autocorrelation time: [41.01967 41.53285]\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAc4AAAHZCAYAAAAPNduTAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAOlBJREFUeJzt3Ql4VNX9xvHDDgFZAyhLQAHZBCSAiuKCe6GAWGTRUBCL1kpBqai4gLhA3UFa61oKjVWpiiKComURcWORTUIFBQICQpAdBIT7f97T/80z2ecmk8y9M9/P8+QBJmO8mWTmnXPO7/xOKcdxHAMAAMJSOry7AQAAghMAAI8YcQIA4AHBCQCABwQnAAAeEJwAAHhAcAIA4AHBCQCAB2VNAJw8edJs27bNnHLKKaZUqVLRvhzAV9TD5MCBA6ZevXqmdGneCwPFLRDBqdBs2LBhtC8D8LUtW7aYBg0aRPsygJgXiODUSNN9YahatWq0Lwfwlf3799s3lu7zBEDxCkRwutOzCk2CE8j/eQKgeLEgAgCABwQnAAAeEJwAAHhAcAIA4AHBCQCABwQnAAAeEJwAAHhAcAIA4AHBCQCABwQnAAAeEJwAAHhAcAIA4AHBCQCABwQnAACxdqwYSk56errJyMgo8H6JiYkmKSmpRK4JAPyE4ESW0GzZsqU5fPhwgY9KQkKCSUtLIzwBxB2CE5k00lRopqam2gDNiwIzJSXF3p9RJ4B4Q3AiB4VmcnIyjwwA5ILiIAAAPCA4AQDwgKlaFJrWOgtC9S2AWENwwjOFoapqVSBUEKpvAcQaghOeqZJWo82C9ntSfQsgFhGcKHR4shUFQDyiOAgAAA8ITgAAPCA4AQDwgOAEAMADghMAAA8ITgAAPCA4AQDwgOAEAMADghMAAA8ITgAAPCA4AQDwgOAEAMADghMAAA8ITgAAPCA4AQDwgOAEAMADDrKOE+np6SYjIyPf+6SlpZXY9QBAUBGccRKaLVu2NIcPHy7wvgkJCSYxMbFErgsAgojgjAMaaSo0U1NTbYDmR6GZlJRUYtcGAEFDcMYRhWZycnK0LwMAAo3iIAAAPCA4AQDwgOAEAMADghMAAA8ITgAAPCA4AQDwgOAEAMADghMAAA8ITgAAPCA4AQDwgJZ7KHbhnLpCj1wAQUFwotgoDHXaSkpKSoH31f0UsDSYB+B3BCeKjUJQYRjOOaAKV92P4ATgdwQnipWCkDAEEEsoDgIAwAOCEwAADwhOAAA8IDgBAPCA4AQAwAOCEwAADwhOAAA8IDgBAPCA4AQAwAOCEwAADwhOAAA8IDgBAPCA4AQAwAOCEwAADwhOAAA8IDgBAPCA4AQAwAOCEwAADwhOAAA8IDgBAPCA4AQAwAOCEwAADwhOAAA8IDgBAPCA4AQAwAOCEwAADwhOAAA8KOvlzvCf9PR0k5GRke990tLSSux6ACDWEZwBD82WLVuaw4cPF3jfhIQEk5iYWCLXBQCxjOAMMI00FZqpqak2QPOj0ExKSiqxawOAWEVwxgCFZnJycrQvAwDiAsVBAAB4QHACAOABwQkAgAcEJwAAHhCcAAB4QHACAOABwQkAgAcEJwAAHhCcAAB4QHACAOABwQkAgAf0qoVvhHP8Gc3qAUQbwYmoUxjq2LOUlJQC76v7KWA56QVAtBCciDqFoMIwnAO5Fa66H8EJIFoITviCgpAwBBAEFAcBAOABwQkAgAcEJwAAHhCcAAB4QHACAOABwQkAgAcEJwAABCcAAMWDEScAAB4QnAAAeEBwAgDgAcEJAIAHBCcAAB4QnAAAeEBwAgDgAcEJAIAHBCcAAB4QnAAAeEBwAgDgAcEJAIAHBCcAAB4QnAAAeEBwAgDgAcEJAIAHBCcAAB4QnAAAeFDWy51RctLT001GRka+90lLSyux6wEA/A/B6dPQbNmypTl8+HCB901ISDCJiYklcl0AAILTlzTSVGimpqbaAM2PQjMpKanErg0A4h0jTh9TaCYnJ0f7MgAAISgOAgDAA4ITAAAPCE4AADwgOAEA8IDgBADAA4ITAAAPCE4AADwgOAEA8IDgBADAA4ITAAAPaLmHwAnnVBh6+AIoLgQnAkNhqNNgUlJSCryv7qeApQE+gEgjOBEYCkGFYTjnlCpcdT+CE0CkEZwIFAUhYQggmigOAgDAA4ITAAAPCE4AADwgOAEA8IDgBADAA4ITAAAPCE4AADwgOAEA8IDgBADAA4ITAAAPCE4AADwgOAEA8IDgBADAA4ITAAAPCE4AADzgPM4Slp6eHtZBzAAAfyI4Szg0W7ZsaQ4fPlzgfRMSEkxiYmKJXBcAIHwEZwnSSFOhmZqaagM0PwrNpKSkErs2AEB4CM4oUGgmJydH438NACgiioMAAPCA4AQAwAOCEwAADwhOAAA8IDgBAPCAqlrErHAaSbDtB4BXBCdijsJQDSRSUlIKvK/up4BlzyyAcBGciDkKQYVhOK0NFa66H8EJIFwEJ2KSgpAwBFAcKA4CAMADghMAAA8ITgAAPCA4AQDwgOAEAMADghMAAA/YjhIh6enpYe0bBAAEG8EZodDU4dSHDx8Oq1ONOtsAAIKJ4IwAjTQVmqmpqTZA80NvVAAINoIzghSaycnJkfySAACfoTgIAAAPGHEi7nH8GAAvCE7ELY4fA1AYBCfiFsePASgMghNxjePHAHhFcBaAxgYAgFAEZz5obIBQFBEBIDgLQGMDCEVEAELF7YjTyxQsjQ3iG0VEAGI6OMMJxF27dplrr72W3rIoliKiSDbzp0Uj4D+BCE7HceyfixcvNpUrV87zfgrMlJQUc+TIkQK/ZqVKlcxbb71VYMP1WrVqmerVq5v9+/cX4soRTypUqGB/r/Q7GCn6euqBnN/v6aFDh7I8TwAUr1JOAJ5tW7duNQ0bNoz2ZQC+tmXLFtOgQYNoXwYQ8wIRnCdPnjTbtm0zp5xyiilVqlShvoZGjApfvbhUrVrVBBHfgz/47eegp/CBAwdMvXr1TOnStJ8Gilsgpmr1YhCpd9J6ofPDi11R8D34g59+DtWqVYv2JQBxg7enAAB4QHACAOBB6XiqeBw7dqz9M6j4HvwhFn4OAGK8OAgAAL+ImxEnAACRQHACABBr21EisY8TiFXh7uPkeQREZj90IIJToUnnIKBonYN4HgGR6cAViODUSFP80qkF8JMdO3aY5s2bZz5P8sLzCCi4I1hBz6PABKc7PeunTi2A3xS0jMHzCChYOMuBFAcBAOABwQkAgAcEJwAAHhCcAAB4QHACABDN4NywYYOZMWOGOXbsWKS/NID/d+jQIbNz5077J4AAB+eqVavM+eefb+bMmWMyMjIK/XWOHj1q99SEfgD4H4Xlpk2b7POC4AQCHJzp6emmR48eZvDgwebFF1+0bYuyC/cglgkTJtgT7d0PugYB/6OgXLt2rR1t7tmzx1SuXJmHBghqcGq0edZZZ5nHH3/cHD9+3Nx///2md+/eZujQoWbatGmZG0vDCc/Ro0ebffv2ZX6oYxAAYwNTz4ft27ebhIQEghOIgoh1Dlq+fLn56aef7N+7detmfvnlF9OuXTv77njp0qVm3bp1Zvz48WF1ZdABwRwSDOSkJZAff/zRzsQoOAEEODi1trlw4ULzyiuv2HBMTU019evXtyPGSZMm2XVPhWirVq0i9b8E4o7enOqUkzJlypg6depE+3KAuFToqdoTJ05k+be6yWtU+fTTT9vpWIWm6J3xjTfeaKdyV65cWfQrBuJ4fXPv3r3279WrV2eaFghScH777bdm4sSJdp3F1aJFC1sUpM8pJD///PPMz9WtW9ecd955pmbNmpG5aiAO1zbffPNN8/XXX9sagoKOPQLgo6la7dPs3LmzrejbvXu3GTlypElMTLSf6969u/nnP/9pbrjhBjNu3DhbYduxY0c7favRKNO0QOFs3rzZLF682B4hVrFiRdO4cWMeSiAIwampIm0V6dmzp+nUqZMZNmyYLQK66667MsOzf//+pnbt2uaBBx4wt99+u6lRo4Zdk5k1axbbSoBCUiGQTqVX0VyVKlWKNE27YsUK+zXyo+dzUlJSof8fQCzzFJx64nbo0MHUqlXL9OvXzz65FJQSGp6XXXaZOfvss20hg8JW00ru5wAULjibNm1qC4LatGlTpIfw4osvDuv/l5aWRngCRQ3OSpUqmUGDBmW+2+3bt68tBBowYID985577rGhqlHogQMHTLNmzbx8eQB50N5NPf9UJ3DmmWcW6XFSLYLeAOdFgZmSkmK3vjDqBCKwxumGpqpqNQLVyFOhef3119ttKJqeffLJJ+2ajBof6J1rOHs3AeRdGKSKdO3f1JvRom5Dad68uUlOTubhBkp6H6f2kSkwtX6p6VqF48CBA83MmTPNd999Z5YsWUK5PBABqqTViLN8+fLm9NNP53kFBLnlnsLSbaOnkeeFF15odu3aZbsIaY0TQNFHm9ripWlabevSti8AAe8cpODUtO2oUaPM/PnzbcVeUYsXAPyvil3buE499VRTtWpVc9FFF9EtCIilJu+tW7e2I822bdtG6ksCJt6DUwU669evt1O0+gAQI71qtd45ZMgQioCACNKyh3o863zaU045xY44AcTQiJPKWSCy1ClIo80ffviBhxaIxdNRAETOxo0bbb2AOm/poIRLL72UhxfwCYIT8CEdlKDALFeunG1GoBoCADE2VQsgcr755hvbwUdbUNQXGoB/EJyAD6dpNeLUNi/1ewbgLwQn4DNfffWVbSqiczfZEw34D2ucgM/s27fPNjxQT1m2oAD+w4gT8Nk07bZt2+y5mzqOr6gN3QFEHsEJ+GyaVqegqJpWfWmLcmA1gOJBcAI+auiuU4XU8EDduGjoDvgTwQn4hBq6q4q2YsWKplGjRkzTAj5FcRDgk4bu27dvtz1p1S2oXbt20b4kAHkgOAGfBKdCs2XLlvYUlA4dOkT7kgDkgeAEfGDz5s22MCghIcHUrFmToiDAxwhOwAfWrl1rp2rLli1rjhw5Eu3LAZAPioOAKFMl7ZdffmkPrdZos1atWtG+JAD5IDiBKK9tvv/++2bPnj32NJQrrrjCNG7cmJ8J4GMEJxDlvZtq5q5OQa1atbJFQTQ9APyNNU4gijQ9W69ePdOwYUM72iQ0Af9jxAlEeZp23rx5Zv/+/XYbCgD/Y8QJRLGS9tNPP7UBWr16dX4OQEAQnECULFy40E7VHjt2jNEmECBM1QJR3Iaiszc12uzSpQs/ByAgCE4gSqG5fv16U7p0adstiJNQgOAgOIEoSE1NtWub5cuXN7169eIkFCBACE6ghH344Ydm9uzZdpq2Tp06pl+/fvwMgAChOAgoYU888YTZunWrKVeunGndujWjTSBgGHECJdwpaPXq1eb48eOmVKlSdpoWQLAQnEAJ0Zrm448/bqdoT548aUeaVNMCwUNwAiUYnK+//ro5evSo/fcFF1xAiz0g3tc4N27caN555x27fnPOOecUuuhBLyzui4uoHRkQdJ9//rnZsWOH/bvO3RwyZEi0LwlANEecWre56KKLbO/NL774wlx//fW2CKIwJkyYYI9Ycj/UABsIujFjxtiTUKRJkyb2+QIgToNz8+bN5tprr7VhOXfuXLN48WLz0ksvmaeeespu8vZq9OjRdh3I/diyZUskLhOImk8++cSsWrUq898DBw7kpwHE61Stihy0btO0aVNz77332k4o0qlTJ1tur897pbMJ9QHEigcffDDz72qxN2DAgKheD4AojjgVlJ07dzZnn322nVZ1aX+a1nG2b99e1P8FEPjR5vz58zP/fdNNN9HUHYjHEafWasqUKWP/rrUad73GcRy7P030p/aruf7zn/+Ytm3bmtq1axf9yoGAGD9+fJZ/33XXXVG7FgBRGnF+++23ZuLEiVlGkwpMNyx/+eUXc+TIERusVatWtbdrGlcn3IcGKRAPzdzVYs9VpUoVOgUB8Tbi3LBhg52a3bNnj9m9e7cZOXKkSUxMzBxlutO3Ck2FqaZrH374YfPss8+aL7/80tSrVy/S3wPgW9lHlxQFAXEWnNrAra0iPXv2tMU/w4YNs6NLvTgoPEODs2LFina0eeutt5qVK1faStuOHTsWx/cA+JKeLwsWLMizSMjv0tLSCryPnvdJSUklcj1AIINTgdihQwdTq1Yt29xAT5r+/fvbz4WGp9Y/tY3k+++/NwcPHjRff/21adOmTfF8B4BPvfrqq1n+rZkatdnzOz2PdUZoSkpKgffV/RSwhCfiiafgrFSpkhk0aFBmm7C+ffva6ViV1uvPe+65x4aq/q5tKG+88YZp0KCBrbAF4s3f//73LP8eO3asCQKFoMIwIyMj3/vpPgpX3Y/gRDzxvMbphqZGlRqBauSpoFTzA61z3n777ebJJ580mzZtsof16h0pEI9FQVrTd6mS/KqrrjJBoSAkDIEIb0dxi380stR0rUJThQ8zZ860BURLly4lNBG3tLYfqlu3blG7FgA+aoCgsNSHAlQjzwsvvNDs2rXLrmmqIQIQr0VBy5Yty3IbnYKA2FHklnsKTk3bjho1ynZHWbFiBYVAiGt/+9vfsvxbU7RBmqYFUEKno6gAaPny5bYzEBCvdLTe/fffn+U2Nf8AEDsich6n1jt1tmBoEwQgXvvShp4lq4MOtIULQOyI2IiT0ES809pm9tFljx49MivRAcSGiAUnEO9mzZpltm3bluU2bc0CEFsITiBCo83hw4dnue3iiy/m+DAgBhGcQASsXbvW7Ny5M8ttTzzxBI8tEIMITqAYGh6oU5AOQgAQewhOIALt9bI3PJg8eTKPKxCjIrIdBYh17kHtucntFBF10gIQmxhxAkX07bffZvl39iIhALGF4ASKYNy4cTluGz9+PI8pEMMITqCQDh8+nCM41XqShgdAbCM4gUJ69dVXc9w2b948Hk8gxhGcQCHdcsstOdpOahsKgNhGcAKFMHfu3By3jRkzhscSiAMEJ1AIV199dY7bdCYtgNhHcAIeLVq0KMdtV155pUlISOCxBOIAwQl49Otf/zrHbW+//TaPIxAn6BwEeLB06VJz8ODBLOfPtmjRgtEmEEcYcQIedO3aNcdtzz//PI8hEEcITsDDaFNND0JVr17dXHjhhTyGQBwhOIEw/epXv8pxW/ZTUQDEPoITCMOmTZvMTz/9lOW2ihUrmsaNG/P4AXGG4ATC0KFDhxy3vfbaazx2QBwiOIECfPLJJ2bv3r05br/88st57IA4RHACBejWrVuuo00aHgDxieAE8vHNN9+YQ4cO5bi9b9++PG5AnCI4gTwoMHPbt/nII4/wmAFxjM5BiGuO4+T5ueXLl5s9e/aYsmXL2k5Bv/zyi7199OjRJXiFAPyG4ATyMGjQoCwBqw9V14a22wMQfyIenFu3bjWfffaZfZferFkz06ZNm0j/L4AS2be5ZcuWHLfPnj2bRx+IcxENztWrV5sePXqY2rVr2xedc845xzzzzDOmSZMmnr7O0aNH7Ydr//79kbxMoEC5tdHr0qWLqVOnDo8eEOciVhy0efNm25JswIABZsGCBWbKlClmyZIlZvfu3Z6/1oQJE0y1atUyPxo2bBipywQKtHbtWvPjjz/muP2tt97i0QMQueD88MMP7dTs+PHjTeXKlW2IJicnmxUrVphp06aZ+fPnh/21VHyxb9++zI/cpsyA4tK/f/8ctyUmJjLaBBDZqVoVTqSnp9ugbN++vXn00UfNnDlzzLFjx2z4aUT62GOPmcGDBxf4tSpUqGA/gJK2a9cus27duhy3b9iwgR8GgMiOOK+88kpz6qmn2o3hffr0MQ888ICZMWOGmTt3rpk1a5Z9Fz916lQ7dZvfFgAgmpKSknLcpjV7ugQBiPiI8/TTTzepqal2XVNrRCrZ79Wrl/2cCirq1atnFi5caKdxKeeHX0ebJ06cyHH7Rx99FJXrARAHnYMUnhpxNmjQwBw5csRO07pUbKEjmHJ7YQL8oH79+rn+Trdq1Soq1wMgjhognH/++ebOO+80kyZNstO3a9assVW2OmVCI06guIW7HHDy5Ek7A5KWlmbKly+fa9Gb+7WYKQFQbMGpd+ha3xw6dKgpXbq0fSevaVqaIcBvFIb6UAV4djVq1OCgagAl13JPzbG/+uorc/z4cVshW7169eL6XwGFpiWFv/zlL7l+TkVtAFCivWpr1qxZnF8eiMgJKGPGjMlxu97sqS8tAGTHsWKIazNnzsz19vfff7/ErwVAMHA6CuK6kfttt92W4/a6devavrQIjwqrCqLOS7ntkQWCiOBE3MqttZ6okA0mrDBUY4iUlJQC76v7KWAJT8QCghNxO9pU8Vr2LSj6t/Ybo2AKQYVhRkZGvvfTfRSuuh/BiVhAcCIunXvuubnezgko3igICUPEG4qDEHemT59u2+tlpwC44ooronJNAIKDESfirnPQzTffbKpWrZo5NeselK4DCQCgIAQn4kr2Zgfqp6yP1q1b2760AFAQpmoRNw4fPmyPu8vNvHnzSvx6AAQTwYm4MWzYsFxv155NHX0HAOEgOBEX1q1bZw8eyM20adNK/HoABBfBibgwYMCAXG/X1C1rmwC8IDgR8z777DPb8CA3I0eOLPHrARBsBCdiXvfu3XO9Pa9CIQDID8GJmPbiiy/m+TlGmwAKg+BETG8/ufvuu3P93MqVK0v8egDEBhogIGY6AmW/n/rRuh2CQqlbEP1VARQWwYmY9N5775mtW7fm+rn169ebUqVKlfg1AYgNTNUi5hw6dCjXA6qlVatWpnbt2iV+TQBiB8GJmJNXaLpbUwCgKAhOxBQdF5ZXh6C8CoUAwAuCEzGladOmeX5u9OjRJXotAGITwYmYMW7cuDw/x+knACKF4ETMePrpp3O9PTk52X4AQCQQnIgJ+TVqZ7QJIJIITgTek08+aX766adcP9enT58Svx4AsY0GCAiUgwcPZvn30aNHzRtvvGHOPPPMLLdXrlzZrFixwrz88sslfIUAYh3BiUDJ3vEnJSXF9qTNTrdpzyYdggBEGlO1CKxXX33VpKen5/q5hIQE07p16xK/JgCxj+BEIO3Zs8dMmjQpz8/n1QQBAIqK4EQgXXPNNXl+7t577zU1atQo0esBED8ITgTOxx9/bI4cOVKoUAWAoiI4ESjbtm2zI8q8vPPOOyV6PQDiT7FU1W7fvt2uQekIJyCS8htNqsK2Xr16POA+lZaWVuB9EhMTOWQc8RecP/zwg2nXrp256KKL7MigY8eOkf5fIE5dcskl+X5++PDhJXYtCJ/CUFXOemNTEN1PAZuUlMRDjPgJzvXr15t9+/bZj8mTJ5sRI0Zk9gl1HCesfXXa1K4P1/79+zOPjKpatWqkLxkB8OKLL5qFCxeali1b5vr5Tz75pMSvCeFRCCoMMzIy8r2f7qNw1f0ITsRVcLZt29Z069bNdO/e3bzwwgu28baOc9KeunCDc8KECbmedKEnVJMmTSJ9yfCBLVu2ZPn3zz//bPbu3WtOPfVUs3PnTjNlyhRz3nnnmTPOOMOsWbMmy3179+5tKlasWMJXDC8UhIQhYkVEg/PEiRP2Y926dea5554ztWvXtiGo/XbffPONOe2008ybb75Z4NdR0I4cOTLLiLNhw4Z2ygfxEaL9+/e30/4PPvigmT17dpbP6fcr1KhRo+wUHwAELjhLly5tw7JTp052VKCRQIUKFcygQYPs1OvQoUPD+jr6b/SRnb42YpuCsV+/fmbr1q3232PGjDENGjQw9evXz3W2YvHixYQmgOBuR3Ff2MqUKWMWLFhg//7222/bUahGjIsWLTJfffVVkf4fhw4dMhs3brQf+jtih8JSI039qWPCfve732XevnnzZjvVH2rIkCGmQ4cOUbpaAPEqoiNOdw3z0ksvtcH2hz/8wU6zLVu2zJ5UoSm18uXL23XQwq5JKSwPHDhgm3jro3HjxvYkDASbwlEjTY04FZo68UTrm3rDNXbsWLNjxw5z/PjxLMVBua2DA0CggtMdceqF78YbbzR169Y1s2bNsv/Whz6vrSpFKeRQSJ5yyil26lcj2U2bNhGeAffFF1+YXr162appvRFyQ1P0e1SzZk271WT37t12f7BonZOpewAx0zmoc+fO9hzEDz/80E6luVNs2ryuAC0KBae+hporKDw1+lTVJYLptddes/szFZotWrTIEpounbWp36GyZcvan39qaqoNWACImc5B5cqVM4MHD7bFQlIcZyJqylbt1/TnsWPHTJ06dZiyDZCTJ0/aitmHH37Y/vvyyy83zz77rKlSpUqO+06dOtX+qT28CtU+ffqU+PUCQLEfZO2GZnHRCEX/D62NqdlCpUqVbDUv/E9r06q0drcmae1b6+EqKstOezndI8IUmn/5y19K/HoBoESCs7hpfUsvwOqLq2m8r7/+2t7eqFEjO/qEP2lvptYzVTCmmQk1ydA6ZvYGCC6Fq05C0T5NhSsARFtgg1PhqA81RdAWFxUJabuLpm817UelrT+sXbs2S2hqCl8/I52XqcYYWgPXfVT0k/0Nj6Zz//GPf9i/d+3a1Vx55ZW2oxAARFNgg9PlFhtp9Lly5crMjfOEp7+EhqYKe55//nm71cSl9Ut1lgql4rL09HQ7Jf/Xv/7V1KpVq9jWzAEgrs7j1EjlnHPOMdWrV7drnzNnzrRTgfBnaKrvbGho5tfYXTRF64YmAERb4EecomlZNZHXaFMnaKjLjHrj3nHHHebXv/4107Y+C81w1qC/++47M3fuXDu6dDsIAYAfxMSIM3S/n5ojqBLzv//9r/n73/9u1z4RHQq/woSmvPTSS/bPK664ghNxAPhKTAWnXpR1gLY6E6kSU0VDOkybsxqj09hA57AWJjQ13a77y80331zMVwoAcRycmrL9zW9+Y7p06WJ74mp/p6b7dEzZkiVLon15cUGP+cCBA831119vj4M7++yzPYWmTJw40RZ7tW/f3lbSAoCfxFRwil6gVVGrCk3t79T2BYXmH//4R/Puu+9yokox0hFfCkq1xFMzAzVhV9cfL6Gp0aY7TavZAipoAfhNTBQHZR91qieupghVlakGCTpV48svv7QvxMJWlcjSG5THHnvM3HfffXbvpbYIvfrqq7Znceg+znA888wzdpq9Y8eO5qqrrorwlSII0tLSCryP9m8nJSWVyPUAMR+cohGOtjBoulab7N0m8HoR1wuzRqFU20YuNNVz9qGHHrL/Vis99ZzVvkyv9HPS4QCi6XVGm/FFYagOUSkpKQXeV/dTwBKeiIaYDE43PEeMGGEOHjxoR0MaCYm2q6hLjaZy1bVGT8CC8AKe01tvvWVDc/r06ebf//53Zmj26NHDfPTRR5n30zpnOHs2NbWrGQJ3tKlqWsQXhaDCMCMjI9/76T4KV92P4EQ0xGxwutO2Q4cOtR1oli9fnnn7qlWrzLBhw2zD8AsvvDCq1xhU2UNTxUAKPPcwc5f60WbvCJQbvQi+8sor9u86MUXbihB/FISEIfwu5oqDstN6m6ZnL7vssiy3r1692txwww22v60qOFH40FQl83/+8x9bgHXLLbfY9ngqFNJZqeHSnluNNs877zzWNgH4WswHp2hvpyo1L7300iy3q9NQ//79bSEL4eltTTM0NPXm48cff7T//umnn8z8+fPtm5UhQ4bY+2pkr9NrfvnllzwraXWAtej+TI0D8LOYnqrNPvJU4ckll1xiG4e7VHV76623mhUrVpgnnngirDXPeHXixAkzduxY8+ijj9p/9+7d24aminp0VqaqahWgaravx1NHhal7kD7U1F0FQ+eee6751a9+Zdcw3YDUFO3Ro0dtFS77NgH4XdwEpxueKkDp27evLVpxqXDob3/7m2nZsqUdJRGeOc2bN8+MHDnShqKcddZZZtasWXarj0JTI0VVRWo9U3s5VSi0e/du271JI/vPP//cPuYqHNKH1p3r169vPvjgA9tlSLTvk9EmAL+Li6naUNobePfdd+f6ueHDh9sA1dQh/kejRR08rTVihWalSpVs4c6aNWtsaCpAFXgKzex0oolG+E899ZQ9V1OtEEXrmApYhaZ+FhrJau+t9tcCgN/F1YjTpdM2tE1lwoQJOT43atQou+apkamqROOZCny031VN83Umpg6f1ihSFHy//e1v7WNU0ChRJ6T8/ve/t/tn27RpY/fWqhViaGgy2gQQFHE34nT3eGo9Tut1udH6XJ8+feyLe7wWDWkaViNAhaYCU9PZCk2dearCn6efftp06tSpwNDUf68tQVoHbdKkiR3Ra100e2hqHycABEFcBqdoHVPBqfXO3KiA6LrrrrOjz0OHDpl4Mm3aNBtoGiGWLVvWNoxQsGm7yYYNG8ztt99u92cWRI/b448/bh9LrWdqFK/Wh4QmgCCLy6naUGPGjLHbJN5+++0cn9M+RDVKWLdund2Un1/RUPaN//mJRgGM2zmpINlHfu4WEo0OJ0+ebD9kwIABdptPXlQlq9Z7qqzV+qe2A+mMVEITQNDFbXC64dW6dWs7wqpbt66dRsyN9iEqWO+//34bGLFYdavgV3/YcGnd8pxzzsn1cyoaUgXu999/bwuJ1BxBx43pNgWxugzpMWd6FkVBM3hES9wGZ/bWfNpOoe0od9xxhx1dZafTVtQ4Xlsq1B0nlsJTYabvSZWvRaXHTqP4zz77zFbRqhBIt6nISFO/3bt3t/8fQhOFRTN4RBvBGVIw1K9fP9snU6348lrXvOeee8ycOXNsRa6axAed2txpFP3ee+/ZUbhGnkU9XkxFVVobVUMJFRNpC5DecGhqVx2HwlkfBfJCM3hEG8GZLTy1Z1EFQfmFp05YOf/8822zBBXMtGrVygSRpk979uxpPv30UzsCzG2k7YWmtHVqigL4kUcesVW0eoy0L7Z9+/Zm5syZdh8oUFQ0g0c0xW1VbX7U+k2Nyi+44IICG5Nr1KmpSe15DNLWFYWk+vQqNHVuqf5d2OlTFVGpetad6tWB4e3atbPTvzt27DBnnnmmbXZQrVq1CH8XAFDyCM48Rp7ax/nPf/4z38pRt1pV05MaqeZVXOQHx44ds8UUKnLSNLNGmppSVWjqc6KCHS9UBPT666/brStuk3ZtVenSpYu5+eabbUVtvXr1bIs9PaYAEAuYqs2nYEi9bd99913zwAMP2GnI/Gha96677jITJ060fVpVcaqRa+3atU1J0qhX7fAUktpG4/6p1nm5nU6ijkBy55132hGopqjDoelXvblQH1pp1KiRDc0WLVpkCU3t3eR8RQCxhOAsgKpntYdT65ja01nQfkhV32pEp7W85557zralK076/2maWFWs+jO/47uqVKlig6158+b2TcHUqVNtwHXt2jXzxBOvWwHUj1ZTshpx60ix7KEZziHWABAkBGeY4XnTTTfZEaTO9FQnnXCqVW+88UYbtgpPbWUpahGR1iE1mnRDUh+bNm3KcT/tSdX/SyGpLTb6UwHmHvm1fPlyWwClgGvQoIH517/+ZatgvdBIVa30UlJS7OOjtUz1oyU04Tfs90SkEZweKIxWr15tpyTffPPNLJ9T4YvW/HKjYNJHzZo1zeDBg+32D7Wgy48C0d22oUYCOs9Sa4VqTh9KFawKvzPOOMNWsepD/50CTB9qcac/czvxRfdTk/a8evbmR1PRV199tV0z/eSTT2x/XwU7I034Bfs9UVxKOUXZuFdCtAdQwaTtEzoMOdrUsFx7OdUJR9OTEm7FqIJN4afmACo80qHOqkBVqGYPTgWepnu1x9KdItbaq867VDVvcnKyvY/WUdVcYMmSJTbEdFtu9N+qYbs+9P/TC0uFChWy3EfFQqoWlvnz59u9rRkZGfbfKijSuuaqVatsE3itnYZS03ft09R6J/z3/HDvp+1UBRW9xQr1SXZ/f/MbkWrmZNmyZfY5hfi030POMOIsBFWIqgBIoaeOQ4WpplXQqapVH6Iga9u2rTn33HPtCFLFSO+//35mYOo8TK0fKmRDt41oJPrxxx/bkaW+Zuh0bcOGDe2HinM0RavRYDj0XkpFTip20ihSQah1UW1d0V5Ml6Z3L774YtOjRw/7oVEv4Cfs90RxIDiLGKAaEWotU8UxRdnHqZGrQvS1116zFbqhgakmC+qpG9pAXffV9O/SpUszb9fIUyMJrcWqEMjrOo+osGjgwIH2a7uj1M2bN9sP0Wi1W7duNig1VcveTADxpliCUy/6GrXESz9STVFqje/Pf/5z5jSnFworrUGGBqYqXUeMGJElMPV5jTCnT5+eOUWsIh2NVBWYKgJyt5cUhoJf03gqftLX0c9P/0+tx2qrisJSHZPcQqJonPICFBeKiBC14Fy7dq0ZP368rbJs1qyZHb3oxTbWab3woYcesvsgNRrUVoxwqm9FoanuO6Im6Oq8o5Fd9p6u6syj4HSnYrX+qLl4Te1GggqfdM0auWodVpXBmo5OTU211wPEIoqIENXg1HmLCkm92GoUpgIaTSUqPIcPHx7219FUpD5CF22DQNsy1JNVHzp3Uq34CmqcIO4oU284VLEruW0zcStqNS2s0a1Gfu+8807Ert99zNXI/k9/+pP9u6pmsxcQAfHcNH7RokV2m1dBYUzjj9gVseDU1KxatukkDK3TiUZOOsx4ypQptnBFxSbhUAOBcePGmaDTCFRvGGbMmGEfg+xVqNmFu16oJ6XXfZdehE7BFmXqF4ilIiJGpnCVjeSLrbrYaIrWpUOMFRzaeqGeplorC6elmw5U1laP0BGnqkODSE82NQrQh0rj1f9WfV3z2vMJwJ8YmSKiwanRpoJTe6DWr19vp2y1fcENTx0tpdtUgdq7d+8CD4HW1GAsTg/qiacOQno8tm/fbvdc6hiuH374IdqXBiBKI1Mth5R0T+twMeVcjMHpTu1pm4KmJ3XE1KRJk+yWCIWqCkvUKF37ARUW2sbghdujwc9rnQX1sHXpJBJNs2oErdG3PtSSTwVFKsZxv0etZ2afjnVPMdGfbjGRpsDDmbZV5W5Bo1z3ewjdD6rryevgaapq/cH9nSmol4n7eVVK+/m5FHQ6vP2rr74yu3fvzvd+WlNVuHp9PSxJ6rmt4kAFaKw79P/nL4fVE8iJsHnz5jkVKlRwbrvtNmfXrl2Zt2/fvt1p166d89lnn3n+mlu2bNF3wgePAb8D+fwO6HnC84jfEV4nTLE+jyTiFSbaf6i2a9ddd52djuzbt6/dZ6jCIbWqK8xapTreqI2cpn0LO8px10n1dfzQtq8w+B78wW8/B71D1gxEQZ2hsj+P/PZ9hCuo1x3kaw/qdXu59nCfR1IspZnaKK8TPFTgo20ZmkrUZnq1kCvMnkNVdkZqr6IeuKD94LPje/AHP/0cwqnIzut55Kfvw4ugXneQrz2o1x3utYe7s6HY9jSoUEh9TdXhRimuY63iYZ4cABDbirVXbZDfnQAAkJu42d2u7S06dzLI21z4HvwhFn4OQf4+gnrdQb72oF53cV17IM7jBADAL+JmxAkAQCQQnAAAeEBwAgDgAcEJAIAHBCcAAB4QnCFNfYNaYMx1wws18z9x4kTgHjS18Fy7dm20LwMgONesWWPuvPNOe3pIEE772Lp1q/nwww9tP+DNmzfb23Td4Z7O4gdpaWn2IwiPd6xR8Og0Hh04f+utt9rWmEGgo/fatGlj7r//frN06VITFHq+Tp8+3R4dtnr16mhfTlzasGGDmTFjRubpUpEQ1yPOlStXmg4dOtj+hO7RXBq9+XUEpydex44d7RFtAwYMMH369LEHhbt9SIMQnqtWrTKtW7c2s2bNMkGls2X1Aq6fwZQpU8yyZctMUK77/PPPt6PNTp06mc8//9yMGDHCPPvss8bvdM7vvn377MfkyZPN8uXLMz/n5+drly5dzBNPPGHP4b3vvvvMd999Z4Jm48aN5plnnjF/+tOfzBtvvGGCRK83+p2fM2eOPcYtYpw4tXLlSqdy5crOnXfe6QTB3r177bFst99+u/371q1bnYcfftg566yznO7du2fe78SJE45frVixwqlUqZJz9913O0H1zTffODVq1HB69erlXH755U7r1q2ds88+25k2bZrjZydPnnTuvfdep2/fvpm37d+/33nkkUfs9T/22GOOn+3evdvp2bOn88ILLzjJycnODTfc4KxZs8a3v/ObNm1y6tev79xzzz3OwYMHndmzZzunnnqq8+WXXzpBsmrVKqdBgwbOZZdd5px//vlO6dKlnccff9wJgs2bNztJSUnOqFGj8n1eFEZcBqce0Fq1ajn9+vWz//7ll1+ccePGOQMHDnSuvvpq5/3333d++uknx2/XfOaZZ2Y5z/TAgQPO9OnTnebNmzvXXXed42fffvutU6pUKeehhx7KfMz//e9/23/re/j6668dv9M133jjjc6gQYMyn3BLlixxhg8f7tSsWdN5+eWXHT8bPHiwc9FFF2W5TeH55JNPOh07dnRSU1Mdvz7uO3futL//esP49ttvO506dXKGDh1qX8x/85vfOH6jgL/kkkuyvDB369bN3j516lR7brHfKfybNm3q3HXXXZlvTl555RWnbt269vnsd++99559zOXYsWPOfffd51xzzTXO7373O/szKEp4FmuTd7/S+prOZ9P0pk5q1xSK1jjr1q1rp2w1BTd69Gg7DZqQkGD8QGcoHj9+3K5Jde7c2d5WpUoV07NnT3PkyBHz1FNPmRdeeMHccsstxm/0Bu3TTz+1f2/WrJn98/LLLzd79+41Bw8etJ+vUaOGnf7s1auX8Stdp9ZL2rVrl7k+q6nz2rVrm/Lly5sHH3zQ/l0/E79dt65XJxZpylNTts2bN8/8vRoyZIi97bnnnjO9e/f2ze+8S89TPa6aXlZNgq5RfUcHDRpkjh49aoYOHWr8Ro95enq6WbFihWnfvr159NFH7XSh1tk03az6hMcee8wMHjzY+JGWfV5//XXTtGlTc++999qfgehnUK5cuUAsC2k6X6dzSbdu3exrvJ67WufXOvm6devM+PHjC1dr4cSpGTNm2Kk2TbtdddVV9h2t3tmKpkATEhLstJxf/Pzzz3akoxGxpk9CHTp0yE5j9e/f3/ErjY41stGoU1NYGiX897//tZ/T9NWAAQOcrl27Ojt27HD8TNM++n3Ztm1bltv1vWgGQyN//Tz8aMOGDU5iYqIzZMgQ+/MIfbednp5ufzZz5sxx/Oq3v/2tnfqUm266yT53W7VqZb8fv02Bfv/993Y0rBGbftf12L7zzjv28f7xxx/tLIVGpBkZGYWeLixuCxcuzHy8XRp5Nm7c2Jk/f77jdx999JFz6aWX2pmgK664ws5WiJa6NMN43nnnFfo1Pu6CM3Q95K233rJhtHjxYvvv0F9gvcBMnDjR8ZPVq1fbaRKtU+lFMNRTTz1l1378+qItR44csdd54YUXOkuXLs3xRqZixYp27dnPNL18+umnO5MnT84MH9e//vUvu26+ceNGx680RVihQgXntttuc3bt2pV5+/bt2+0aeuhSgF+4z8t//OMfztixY51bb73VOe2002w4adq2SZMmzu9//3v7++Unur433njDXnOfPn2yfO7Pf/6zfbz9ds3u4CGvn4FeP08//XRn7ty5mZ/7+OOP7cDDb9eelpbm1KtXz7650iAplN4oanCk52xhxF1wZg9IPbgazYXSC59+qUN/OfwS+F988YV9cdYTMXSdROs9GnUePXrU8bN9+/bZIiH3Ot3vSy/Y+gX3c+i4tOajJ91LL71ki1Zc69ats9+D/vSzmTNn2vC89tprnddff91Zu3atHVkojLZs2eL4lUZAGrmpyCb0jZfedCmk/Eq/JyrgC31u3nHHHbbATIVDfqFZE80Khc6mhL5WHj9+3F5v06ZN7euQjB492v5MfvjhB8dv1y6zZs1yypYt69SpUyfLm0L9LDQa/eCDDwr1/4vL4JT8pkfuv/9+W63qDu1LkoIk+zsnN1zc2/WioUpIjTAV8HoCVq1a1QZSUGkKtHPnzs6ePXucIMxW/PGPf7QFQapU/eqrr2yAqkJbo5/QkZxfLVu2zLn44oudRo0a2WtW4c3y5csdP1OBh4pT3FkJv05xZqfpwGrVqtlqVFVf641X9erVcyy5RNP69evt77NCUGGY2++wfv81Qm7SpIl9DVJhn97E6/ffz9f+2muv2WpgLbHo77q/3ihqNKqRZ2HEfHDmNfWQG037DBs2zP6SR6PKU08wldmr9FtTT3q3lP37cP9Ula17vdpKoJFz0B5v0ehM7761XuWXaVqFoB5PVQ5mH8GHfn+abrvgggvs6K19+/Z2xOb38Mk++tcIXy/gQQh7v249CYdmhxQ4zZo1s2ubfvldF40itU6squu//vWvNoD0Rjav34n27dvbquby5cvbqvJoCvfaNZ2sN+Za6mrRokWR3yjGbHCGFpkUNG/vevrpp+2DG413ggoQBbYKfPRuSCNJbRHQvk2X+yLu13faeU2XhAq9dq3Z3nLLLfaJ6JfRsq5J19OmTRsbiCoUy/77oykrl97ALFiwwE4jRmOGAsGhN2R6XfLbrMrhw4dt6GjaXrQum1sA6Xmg76FatWpOmTJlfDFiDvfaRYVYejOsQVFR3yjGZHBqtFCuXDmnR48eYY2EQtepovFLXdDmdK1fhlJ1nirz/CScqZ7cAl9ThipM8QON+LW/V1Ou+rtbBRw6nRPUEQ+Qn+xrrQoi/e7ruaDAcd8w6nn9wQcfZDafCNK1R7J+IuaCU6MdlYGr1PiMM86wBRD5hefIkSNtuXi0iwsK2pw+YcIEe5umb9XJQ5t5/fIi7nWqR2s9Dz74oOMnulY9/iNGjMgS9Nr+o6ICvUsNLZyZNGmSM2XKlChdLVA89BrpvsHVeqD7XFbxj5ZUevfu7dvK/YKuXVmg16pIzNjFXHCqvFgPkKbP3nzzTbvnKDQ8Q6fZRPfRSClaVWHuD/HZZ5+162XZKzLVwcjtkOJO1Y4ZM8b57rvvHL/wMl2i0b32O5577rmZ7wb9QNcyfvz4LB1RVPyg70Ojfr1ZUXHBokWL7PegNR6FqtYJgVii1yT3Tbme05q9U3cyTc/6vcPXyTyuXZW1kbz2mAtOvRt69913M0NS7dyyh6f7zsR9gLPvx/Pr5nS1kPKrcKZL9LhrKlzBk986aLRohO9y37HqTYCuV2uYCkvtyROt72h9E4hFeu1xX3+0bUODCz+safrl2mMqOHObulT5tDatZw/PF198MXN055dimyBuTvcyXaJCJ/WKzL5v1q99OrX+Gkp78fThl98XoLify5ri1HPYT1XAfrj2mOpV6/ZTDFWxYkXTvXt3249Q527qKK4GDRrYo5TcI378ci5k165d7Tmb1113nT20t2/fvqZt27Zm2rRpZufOnba/rt+VKVPG9ulUL8v+/fvbx3bgwIFm5syZ9vFWb2D1GfW7Ro0a2Q/R96Ieo+oNrJ+HX35fgOKmIwDV81W/90HTuhivvZTS08QQNfJ1z9YMbXD9888/m3fffdc2cFdD8Y8++sg2vfYj/bBHjhxpNm3aZL8XhZEaLqtZdFC4v1Z67C+77DLb7HrBggX2MOIgGjNmjJk6dar5+OOPMxvVA7HOff0MIqcYrz2mRpw6oFdBo8BZuHChPT3BfeA08tSLnk5+WLRokWnVqpXxKwW6Rmjq7H/gwAFz2mmnmcTERBMketz18xg1apSZP3++Dc4ghqZmAPS7pDcuerNFaCKeBDU0i/vac85tBnikqZGZQlNHJs2bNy/L52fPnm1fAPXh59B0Va1a1TRu3NiGTdBCM1amekS/K7t27bJvtoI04gdQfGJiqtadnlVoarSm8/p0NmXolK3OrNQ5eKeeempUrzXeBHmqx6VzUHUGIQDERHBmD00dIvzyyy9nCU0Vd+RWOAQAQFwFp9bQ3OnZvEITAIBICvQwTKG5efNmu452zTXXmFdeeYXQBAAUq8CPOG+++Wa7hvb8888TmgCAYhfo4JQ9e/aYatWqsYYJACgRgQ9OAABKUqDXOAEAKGkEJwAAHhCcAAB4QHACAOABwQkAgAcEJwAAHhCcAAB4QHACAOABwQkAgAcEJwAAHhCcAAB4QHACAOABwQkAgAnf/wG0O5BBTcyyDAAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sampler = hemcee.HamiltonianEnsembleSampler(total_chains, \n", " dim, \n", " log_prob, \n", " step_size = 0.01,\n", " adapt_step_size=da_parameters,\n", " adapt_length=chees_parameters,\n", " )\n", "\n", "samples = sampler.run_mcmc(keys[1], \n", " initial_states, \n", " num_samples=10**4, \n", " warmup=10**5,\n", " thin_by=5,\n", " show_progress=True)\n", "\n", "try:\n", " tau = hemcee.autocorr.integrated_time(samples)\n", " print(f'Autocorrelation time: {tau}')\n", "except Exception as e:\n", " print(e)\n", "\n", "corner.corner(samples.reshape(-1, dim).__array__())\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "771b85d6", "metadata": {}, "source": [ "## Turning off adaptation\n", "\n", "**I want to preface, this is something you should NOT generally do**, and this should only be done in extreme scenarios. But, if for some reason adaptation is giving you a headache, you can disable it. Just plug in `False` into the `adapt` parameters." ] }, { "cell_type": "code", "execution_count": null, "id": "2661dc24", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/clarkmiyamoto/nyu/h-emcee/src/hemcee/backend/backend.py:47: UserWarning: Explicitly requested dtype float64 requested in zeros is not available, and will be truncated to dtype float32. To enable more dtypes, set the jax_enable_x64 configuration option or the JAX_ENABLE_X64 shell environment variable. See https://github.com/jax-ml/jax#current-gotchas for more.\n", " self.accepted: jnp.ndarray = jnp.zeros(self.nwalkers, dtype=self.dtype, device=self.device)\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Using 100 total chains: Group 1 (50), Group 2 (50)\n", "Starting warmup...\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "100%|██████████| 2000/2000 [00:23<00:00, 84.03it/s]\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Warmup complete.\n", "Found step size 0.01 and integration length: 10\n", "Starting main sampling...\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "100%|██████████| 1000/1000 [00:11<00:00, 87.17it/s]\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Main sampling complete.\n", "The chain is shorter than 50 times the integrated autocorrelation time for 2 parameter(s). Use this estimate with caution and run a longer chain!\n", "N/50 = 200;\n", "tau: [210.34383 205.32188]\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdMAAAHZCAYAAADdZdHLAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAPMVJREFUeJzt3QmcjXX///GvLYwiTJQ0UlIj3DWlonBX2kvllqihKIp2RdbSqlWlX5u4laYFrVJu7futkiUxKYUhKZPsO9f/8f7ej2v+Z2bOzJyZ6yzXdc7r+Xich5lzjjPXWd/nu32+lRzHcQwAAKiwyhX/rwAAgDAFACAKaJkCAOARYQoAgEeEKQAAHhGmAAB4RJgCAOARYQoAgEdVTQDs2bPHrFq1yuyzzz6mUqVKiT4cwFdUd2Xjxo2mUaNGpnJlvh8DiRCIMFWQHnTQQYk+DMDXVqxYYRo3bpzowwBSUiDCVC1S98Oidu3aiT4cwFc2bNhgv2y67xMA8ReIMHW7dhWkhClQ+vsEQPwxwAIAgEeEKQAAHhGmAAB4RJgCAOARYQoAgEeEKQAAHhGmAAB4RJgCAOARYQoAgEeEKQAAHhGmAAB4RJgCAOARYQoAgEeEKQAAqbAFG+InLy/P5OfnR+320tPTTUZGRtRuDwD8iDBFoSDNzMw0W7ZsidqjkpaWZnJzcwlUAEmNMEUBtUgVpDk5OTZUvVKIZmdn29uldQogmRGmKEZBmpWVxSMDABFiAhIAAB4RpgAAeESYAgDgEWEKAABhCgBAYtEyBQDAI8IUAACPCFMAADwiTAEA8IgwBQDAI8IUAACPCFMAADwiTAEA8IgwBQDAI8IUAACPCFMAADxic/AUkZeXZ/Lz80u9Tm5ubtyOBwCSCWGaIkGamZlptmzZUuZ109LSTHp6elT/fiQhrb+ZkZER1b8LAPFCmKYAtUgVpDk5OTZU4xVqui2Fc3Z2dpnX1fUUugQqgCAiTFOIgjQrKytuf0/BqICMpHtZgavrEaYAgogwRUwpHAlIAMmO2bwAAHhEmAIA4BFhCgCAR4QpAAAeEaYAAHhEmAIA4BFhCgCAR4QpAAAeEaYAAHhEmAIA4BFhCgCAR4QpAAAeEaYAAHhEmAIA4BFhCgCAR4QpAAAeEaYAAHhEmAIA4BFhCgCAR4QpAAAeEaYAAHhEmAIA4BFhCgCAR4QpAAAeEaYAAHhEmAIA4BFhCgCAR4QpAAAeEaYAAHhU1esNILHy8vJMfn5+qdfJzc2N2/EAQCoiTAMepJmZmWbLli1lXjctLc2kp6cbP4sk9HUfMjIy4nI8ABApwjTA1CJVkObk5NhQDWoI6dgU9tnZ2WVeV9dT6Pr1vgBITYRpElCQZmVlmaBSMCogI+muVuDqeoQpAD8hTOELCkcCEkBQMZsXAACPCFMAADwiTAEA8IgwBQDAI8IUAACPCFMAADwiTAEA8IgwBQDAI8IUAACPCFMAADwiTAEA8IgwBQDAI8IUAACPCFMAADwiTAEA8IgwBQDAI8IUAACPCFMAADwiTAEA8IgwBQDAI8IUAACPCFMAADwiTAEA8IgwBQDAI8IUAACPCFMAADwiTAEA8IgwBQDAI8IUAACPCFMAADwiTAEA8IgwBQDAI8IUAACPCFMAADwiTAEA8IgwBQDAo6pebwCIt9zc3DKvk56ebjIyMuJyPABAmPpUXl6eyc/P9xwqyUQBmZaWZrKzs8u8rq6nx4dABRAPhKlPgzQzM9Ns2bIlotBQyKQCBaMCMpIvGQpcXY8wBRAPhKkPKQQUpDk5OTZUS5Nq3Zm6r6l0fwEEA2HqYwrSrKysRB8GAKAMzOYFAMAjwhQAAI8IUwAAPCJMAQDwiDAFAMAjwhQAAI8IUwAAPCJMAQDwiDAFAMAjwhQAAI8IUwAAPCJMAQDwiDAFAMAjwhQAAI8IUwAAPCJMAQDwiDAFAMAjwhQAAI8IUwAAPCJMAQDwiDAFAMAjwhQAAI8IUwAAPCJMAQDwiDAFAMAjwhQAAI8IUwAAPCJMAQDwiDAFAMCjql5vAOWTl5dn8vPzS71Obm4uDysABAhhGucgzczMNFu2bCnzumlpaSY9PT0ux5WsIvlSosc4IyMjLscDIHkRpnGkFqmCNCcnx4ZqafiQrzg9dvoykp2dXeZ1dT2FLoEKwAvCNAEUpFlZWYn40ylBwaiAjKQ7XYGr6xGmALwgTJGUFI4EJIB4YTYvAAAeEaYAAHhEmAIA4BFhCgCAR4QpAAAeEaYAAHhEmAIA4BFhCgCAR4QpAAAeEaYAAHhEmAIA4BFhCgCAR4QpAAAeEaYAAHhEmAIA4BFhCgCAR4QpAAAeEaYAAHhEmAIA4BFhCgCAR1W93gD+Jy8vz+Tn55f6cOTm5vJw+VAkz0t6errJyMiIy/EACB7CNEpBmpmZabZs2VLmddPS0uwHMxJPz4Oej+zs7DKvq+spdAlUAOEQplGgFqmCNCcnx4ZqaWjh+IeCUQEZSY+CAlfXI0wBhEOYRpGCNCsrK5o3iRhTOBKQALxiAhIAAB4RpgAAeESYAgDgEWEKAIBHhCkAAB4RpgAAeESYAgDgEWEKAIBHhCkAAB4RpgAAeESYAgDgEbV5gQixVRuAkhCmZWCfUrBVG4CyEKalYJ9SCFu1ASgLYVoK9imFi63aAJSGMI0A+5QCAErDbF4AADwiTAEA8IhuXiDKWEIDpB7CFIgSltAAqStlw5T1o4g2ltAAqSslw5T1o/DDEhq6g4HkkZJhyvpRJBLdwUDyCVSYzps3z+y9996eb8dtEbB+FEHoDv7888/ta7UkmzZtisFRAki6MHUcx/7bsWPHqN1mzZo1TfXq1c2GDRuidptApPbdd197Ko1en3qdKlDL8z4BEH+BCNONGzdG/Ta3bt1qWrZsGfXbBRL5PqlTpw5PAJAAlZwAfJ3ds2ePWbx4sWnRooVZsWKFqV27tklmai0fdNBBKXFfU+3+xuK+6i2sIG3UqJGpXJk6LEAiBKJlqg+IAw880P6sD6Bk/8B1pdJ9TbX7G+37SosUSCy+xgIA4BFhCgBAqoSpZjbefvvt9t9kl0r3NdXubyrdVyCVBGICEgAAfhaYlikAAH5FmAIAkApLY7TOdNWqVWafffYxlSpVSvThAL4S6TpT3kdA7NZrByJMFaRa6A6gZCoE0bhxY95HQAzfR4EOU7VIJRUq5ADlsWzZMns6//zzC94nJeF9BJRdnays91Ggw9Tt2k2lCjlAWTZv3my2b99esJNSWUMgvI+AslV0KDEQYQogfJjWrVvXjoUCSCzCFAioWrVqFfoXQJKE6dKlS82bb75pVq5caY477jhz8cUXV+h21HWlk4s9R4HCvv32WzN37lzTpEkT07p1ax4eIFnWmS5YsMB06NDBvPPOO2bWrFnmkksuMQ8++GCFbmv06NF2Fwz3xExeoLCcnBwzadIk8/zzz5s1a9bw8ADJEKbLly83Xbp0sQH63nvvmS+//NI8++yz5uGHHzY///xzuW9v6NChZv369QUnzeIF8D9//vmn3d937dq15q+//uJhAZKhm1eTH1555RXTrFkzM2zYsILFrm3atDHVqlWr0OQIFQGnEDgQ3owZM8yOHTvszPaOHTua/fbbj4cKCHrLVOHZtm1bc9RRRxXaoPjII480VatWNb///rvXPwEgxMcff2xbp1WqVDGnnXYaE5CAIIfp7t27C37WWKnGOSV0Exqt19m5c2fB7x9++CHjO4DH5TCrV68269ats70+W7du5fEEghqmP/30k3n00UcLtTrdEFWA7tq1y77J9c3ZLbKgLmB9iw4NVwDlM336dPu+0xBKgwYNmJwHBHXMdMmSJbZb9++//7aTHwYOHGjS09MLVY1Q16+CVAGrrt677rrLjB071nz99de2iDCA8vvss8/MiBEjbJjqPacJf02bNmXpGBC0MFUXk7pzO3fubCcYXXvttbYVOnjwYPvmDg3TGjVq2FZp//79zfz58+0M32OPPTYW9wFIenrvaaJfXl6enXyk+Qknn3xyog8LQEXCVCF5zDHHmPr169uCDArQ7t2728tCA1XjqVrS8uuvv5pNmzbZxeWtWrUqz58CEELF7GfOnGmDVNTDo25eAAEM05o1a5rLLrusYPZgt27dbFdujx497L9DhgyxQaufNTli8uTJdisbzewFUHEaXtGXU1dFq4sB8MmYqRukan2qpao3tcJT4zcaN73xxhvNQw89ZL9Jq0pLWlpaLI4bSBlaBqOxUpfmI/Tq1SuhxwQgSkUb3AlGaoGqq1dB2rNnTzNt2jT7LXr27NkEKRAFqnf9ww8/FPx+wAEHeOrinTdvXsG2bSXRkE1GRkaF/waQajxVQHJn8CpU1UIdN26cfaMyRgpEz+OPP17od31p9UJVk8qiHqXc3FwCFYhXOUEFqrp8Bw0aZCuzKEyZbAREb3eY0Fap9O3b19Nt6kuvJhKWRCGanZ1t8vPzCVMg3luwaZLRnDlz2A4KiKIpU6YU+v3ggw/2PIv38MMPN1lZWR6PDEDUw1Tjp3369ClUuAGAd6+//nqxHZXYDBxI4v1MCVIg+oUaQpfDKEQvvfRSHmYgmcMUQHQVXf5y9dVX0yoFfIowBXxo6dKlxbp4KdQA+BdhCvjQ+++/X+h31eJt0qRJwo4HQOkIU8CHnn766UK/n3DCCXTxAj5GmAI+7OJV4ZNQvXv3JkwBHyNMAZ+54447Cv3eokULc+655ybseADEsWgDgMipBGc4a9asMS+//HKh8+69915apYDP0TIFfGTSpEkFe5ZKw4YNzfnnn5/QYwJQNsIU8IktW7aYp556qtB5w4cPT9jxAIgcYQr4xIsvvlio4pGoTCcA/yNMAZ94/vnnC/1et25d9gQGAoIwBXzSxbto0aJi5QMBBANhCvjAF198YdatW1foPIraA8FBmAI+aJUWbYU2a9bM7l0KIBgIU8AHE4+WLVtW6LyHHnqI8VIgQAhTIMHGjBlT6Pdq1aqZzp07J+x4AJQfYQokkCYdLV68uNB57dq1S9jxAKgYwhRIoNtuu63YedThBYKH2rxAgqgOb9ENwKtXr27OOussXzwnubm5ZV4nPT3dZGRkxOV4AD8jTIEEmTFjRrHzxo4da3eJSSQFZFpamsnOzi7zurqeQpdARaojTIEEGTZsWLHz+vbtaxJNwaiAzM/PL/V6uo4CV9cjTJHqCFMgQV28q1atKnRex44dffNcKBwJSCByTEACEiBcC3TKlCk8F0BAEaZAnC1dutRMmzat0HlNmjQx++23H88FEFCEKRBnTzzxRLHz+vXrx/MABBhhCsTZ448/Xuy8K6+8kucBCDAmIAFR5DhOmROPdu7caSpVqlRw3VatWtHFCwRcTMJ0z5499t/KlWn4AmV15w4fPtyGK4DgqhqLWqP33XefWblypd1G6pRTTjHdu3eP9p8BAkfvjaITj2rXrm1OPvnkhB0TgOiIatPxxx9/NCeddJLZa6+9bH3RvLw8M3LkSHPdddeV63a2b99uNmzYUOgEBN0tt9xS7LwHH3zQNGjQICHHA8CHYaoAvPvuu03Pnj3N+PHjzcCBA82bb75p9tlnHzt78ZJLLon4tkaPHm3q1KlTcDrooIOidZhAQmisdObMmcXOv+CCCxJyPAB8GqYq0L169WpTr149+/u2bdtMjRo1zGmnnWa6dOlit5nShseRGDp0qFm/fn3BacWKFdE6TCAhRowYUew8TTyiVQokh6iEqWYlbtmyxezYscP88ssvZteuXTZIf/vtNzN58mRzzjnn2OLd7777bsTBrLGk0BMQVHpvqLemqAkTJiTkeAD4dAKSZiJq9wh1z3bo0MEsX77cVnTR9lI9evQwvXv3Nm3atLGbHquF2rx5c2YvImX8+9//LnZelSpV7HsCQHKI6gSkE0880cyaNcsWyFbr8oEHHjDPPvusvezXX381jRs3Nvvvvz9BipSiYYuiwrVUAQRX1JfG6Nv2pEmTigXm559/bho2bEiQIqXoda9u3qK6du2akOMBEKCiDaFBumDBAvP000+bnJwc89lnnzH+iZSqgHTttdcWK17SsmVLU7NmzTgfGYDAlhPUcpklS5aYtWvX2m/orVu3juWfAxIu9IuklsMsXLiw2HVefPFFemiAJBPTMNW46dlnn21OP/10U6tWrVj+KcB3nnrqqWLn1a1b185sB5BcYl7oXoGqE5BK1Cq98847i51/8803J+R4AMQWleiBGAhXoETLYa644goebyAJEaZAlGn27sMPP1zs/EGDBrHVGpCkCFMgyqZPnx72/GHDhvFYA0mKMAWi3CoNN1bavn17WyUMQHIiTIEo71mqrQiL0s5JAJIXYQpE0fXXX1/sPG1DyHIYILnFfGkMkIyVjUpqlc6ePbvY+drnN9z/L1pyE0BwEaZAlIwcOdLs2bPHnlzVqlUz11xzDY8xkOTo5gWiVKThrbfeKnb+1KlTeXyBFECYAlFQUjGGzp078/gCKYAwBaKwHObdd98tdn7z5s15bIEUQZgCHr300kthz3/ggQd4bIEUQZgCHlulV199ddhWKV28QOogTAEPJk6cGPb8nJwcHlcghRCmgIdW6Y033ljs/Pr161OkAUgxhClQQd99913YYgzaMYY6vEBqoWgDEIGiobl161bTrVs3U6NGjWLXPeuss3hMgRRDmAIRKFr6T0th1q9fH3bP0gYNGvCYAimGbl6gAu66666w59966608nkAKIkyBcsrNzTU///xzsfMPOeQQxkqBFEWYAuWswdupU6ewl40aNYrHEkhRhClQDm+//bb5+++/i51/6KGHmo4dO/JYAimKMAXKsa70hhtuCHvZlClTzH777cdjCaQowhSI0CuvvGJ2795d7Py99trLZGZm8jgCKYwwBSJ0yy23hD1/zJgxPIZAiiNMgQirHW3bti3sZd27d+cxBFIcRRuQ0sKVAwwnOzvbVKtWLew2azVr1ix0O0ULPKTCUqGypKenm4yMjLgcD5AIhClQhi+++MKsXLnS7Nixo9hlPXv2TLnwDA1I1SDWF42y6HoKXQIVyYowBcqYwXvBBReEvWzEiBEpXaRBwaiAzM/PL/V6uo4CV9cjTJGsCFOgFHPmzDGbNm0K28UbblPwVKNwJCABJiABpTr77LPDnn/TTTexrhRAAVqmQAndu19++WXYcVKhVQogFEtjgDA2b95sunbtGvaxueqqq8zBBx/M4wagAGEKhGmVfvrpp2b79u1hHxtapQDiEqb6MCqpewwIQqtUS17CycrKMi1atIj7MQFIsTD94YcfTLdu3cysWbNK/GYP+NnUqVNLLObw8ssvx/14AKTYBKSFCxea9u3bm4svvtg0bdrUVK9evUK3oxAODeINGzZE8SiB/2/Pnj3FHo577rnH1K9fv9B5WhrTpEkTxkoBxDZM1TU2cOBA06NHD/Pkk0/a83788Udbz7RevXrlWos2evRoc8cdd0Tr0IASFa1epPKA4YYodN7dd9+dstWOAMSpm7dq1ap2rLRv3752m6ozzzzT9OrVy3To0MG2VCdMmBDxbQ0dOtSsX7++4LRixYpoHSZQIr1+H3roobCXHX/88eaII47g0QMQ25bpunXrzOLFi23JsEGDBtnzxo8fb1atWmU++ugjW3qtTp06JS43CKXu4Yp2EQMVNW7cuBIv+7//+z/ToEEDHlwAsQ1TfdCceuqpZtq0aWbZsmW2Qkzr1q3tqWXLlub33383H374obnwwgtN5cqV6S6Dryxfvtzce++9YS9r06aNHS8FgJh382os6eabbzYTJ04077zzTqFxp8aNG5uGDRuaRYsWEaTwpd69e5d4mcZKASBuS2OOPfZYM2PGjIIuM83ude3cudM0b97c7Nq1K5p/EvBMQxNa0hVOp06dzNFHH82jDCC+tXm1NOaTTz6xs3r79OljWrVqZVup6v7VvpDhdt8A/Noqve666+J6LACCKSYVkDSDV5OOTj/9dDsWtffee9sg1dgp4Cdz5841X3/9dYmt0rZt28b9mAAET8x2jTn88MPNXXfdVbAoXpOOAD/R2ujOnTuXePltt90W1+MBEFwx34KNEIVfKhsV9eKLL9qCIlrWVdSAAQPsF8KSygoCQCj2M0XKTjoaNmyYqVmzplm7dm2xy4cMGZKQ4wIQTPS9IiWpbGBJbr/99rgeC4DgI0yRcjQp7rnnnivx8muuuSauxwMg+AhTpBxtpFCSmTNnmrS0tLgeD4DgI0yRcmOlr7/+eomXn3jiiXE9HgDJgTBFSu0KU1o1o2eeeSauxwMgeRCmSBmqGx266XyoKlWq2E0YAKAiCFOkTKu0tA3n3333XcZKAVQYYYqUkJ2dXeJlBx10EMXsAXhC0QYkpa1btxb8rApH2k/30EMPLXa9GjVq2I0ZAMALwhRJSfvryrZt2+zuReFKBso///lPNqqPk9zc3DKvk56ebjIyMuJyPEA0EaZIatqQftmyZSVeTjH72FNAau1uaV3tLl1PoUugImgIUyQttUqvuOKKEi/v0qWLqVu3blyPKRUpGBWQWuNbGl1HgavrEaYIGsIUSWv8+PGlXl5a0CK6FI4EJJIZs3mRlFatWmUmTJhQ4uW9evUyjRo1iusxAUhehCmSkrpwS3P55ZfH7VgAJD/CFEln1KhRJVY6crt3tY8pAEQLYYqksnDhwlIrHWnc7tJLL7XrSwEgWghTJI0///zTrhstjfYqZQYvgGhjNi8CZe3ataXu+qLSgDo1bdrULF68uNDlZ555pmnTpk0cjhJAqiFMEWh///23LWJfq1Yt88YbbxScv3LlSvPzzz8Xuu7XX39NMXsAMUGYIrAWLFhgt03buHGjadWqVanX/eCDDwhSADHDmCkC6YcffrBBqpbprl27zNy5c80ff/wR9ronnXSSycrKivsxAkgdhCkCWW/XDVKFZN++fQsKNah713GcQtefNm0arVIAMUU3LwK3tdpll11mJyIdddRR5tVXXzV16tSxk45UtH7NmjV2tu5ee+1lrz9ixAiCFEDM0TJFoIwZM8b8+uuvZv/99y8IUnfJS7169ezPlSv/72VdrVo1M3DgwIQeL4DUQJgiUOOkY8eOtT/ff//9hdaLaocYd9mMQlS06be29AKAWCNMEQh79uwx/fr1s5ONzjrrLHPuuecWunz16tUFm4JXqVLFtlyPOeaYBB0tgFRDmCIQVJDhv//9r11PqlZpUZp8JBorVaAOHTo0AUcJIFUxAQm+mVhUEgXlkCFD7M8333yz2Xfffc3mzZsLXWfZsmX233322cf8+9//NvXr14/xEQPA/0eYwvduueUWs2HDBnPssceaq6++Ouw46Lp16womH6mUoKiFCgDxQDcvfG369OnmzTffNFWrVjVPPPGEHQ8Nx+3mPe+88+J8hABAyxQ+pjKBN910k/35hhtusCUDdV5pYXrYYYfF9RgRfbm5uWVeJz093W6nB/gF3bzwpd9//93O3v3tt99st21ZE4p0fTnggAPidISINgWkuvCzs7PLvK6up9AlUOEXhCl8Jy8vz5x44onmr7/+sr9rbWlp60U/++wzM3/+fPszH67BpedOAZmfn1/q9XQdBa6ux/ONpA1TLZzXJs0a22rSpElBWTcgEjt27LAflApStVTuvvtuc+qpp5YapBdffLHZvn27OeOMM2yJQQSXwpGAhEn1CUiqUNOpUyfTrVs3O771wAMPmN27d0fzTyDJDRs2zMyePdtWN/r8889Nr169Srzupk2bbJ1eLas5/fTTzfPPP88MXgDBDlPt5PHPf/7TtiJeeeUVc88999jC4+7EkPJQK0NLIUJPSH7a3PvJJ5+0Pz/77LO2Z6M0L730kt055pBDDjGTJk0yNWrUiNORAkAMunk1dtG/f3/bPffggw/a8zIzM+2GzNoSS112WkSvnT0iMXr0aHPHHXdE49AQEEuWLLFrSN3CDGeffXaZ5QWffvpp+/OAAQMIUgDBD1Mtjj/zzDNN165dC87TWNfMmTNtzVSF7ZFHHmm3w9JGzWXRzM3Q3T7UMo00iOEvCr2yqJv2/PPPt8teWrZsaU455RQ7FhqOijIcd9xx5j//+Y/dPUbVkHr06BGDIweAOIepWp3XXnutLeUm6ua9/fbb7b8aQ9VYqqrYfPjhhxGFafXq1e0JqUFrSJcuXWq3U9PrRAFcWgjreuPGjbM/X3XVVaZRo0ZxPFoAiOFsXjdIpW3btnYSSVZWlv29Q4cOpkGDBua7776L1p9DktBY54QJE2zvhoI0kpq6c+fOtdurqSqSvsQBQFKuM9XEEXfyiFoYWu6w9957m9atW8fizyGg1GOh8U5RV22ky1oeeeQR++9FF11kGjduHNNjBABf1ObVGNe9995rt8/Shx8g69evt+tDNV562mmn2eVUkdBkNg0fSOi4OgAkbZhOnTrVdsNpuYOKlVM3FTJt2jS7DvnHH380Bx54oHnhhRfsl66ybNu2zTz88MNm586dduxdu8gAQNKHaYsWLcyaNWvs4vujjz46ln8KAbBr1y7Tu3dvc+GFF9qau4ceeqj9krXffvtFtPZYM8QXLFhgx+cfffTRuBwzACS8Nq+Ww+Tk5Jhq1arF8s8gADR23qdPH/Piiy/aUpNaS6qiHjVr1izz/2rMXUVAVH9X19eymGOOOSYuxw0Avih0T5DCcRxzzTXX2CDVDFx1/3fu3DmiB0Zduhpz1wxeLZfSkqt27drxoALwlUBtDr5s2TKzefPmRB8GyhmkWvKidaFa/qKlMOUJ0vvuu88uqdKGCQpS9XYAgN8Eags2jb1q38rmzZvbdau1atVK9CGldEhGQpOE5s2bZ39Wi/K9996zJwVrKG2xpl1fQsdXtVHCN998Y4N05MiRtjpSpH8XyY9NxOEnVYPWMtWsT5WQ27Jlizn44IMJVB9TGLpBqglommj09ttv28lEKuQRWuhDS2TOPfdc+7MCU3WeZ82aZbt233rrLbsrDCBsIg4/ClSYNmzY0KxYscLWZFUhfbVmmjZtmujDQhhaDjVkyBD7s4p16Hn6+OOP7S4v8umnn9odhsJNQNKktZdfftmOt7/22msEKQphE3H4UaDCVONl+hBWV6/WHB5//PGJPiSEoX1F3TJ/Wh6l9cV63hSkamkqJLUXqc47+eSTC9Vh1pZ9qtUro0aNKnP3GKQmNhGH3wRqApLKzWnvSrVmNF6qcnSTJ0+2u9P8+eefiT48GGOmTJlirrjiCvtYKBTVg/DFF1/YnYMUoh07drT73mrvUVVB0ji4xkfd7l1tw7Zu3To71jpo0CAeUwCBEKgwVYBecMEFtvqNZnq++uqrdoNodQVq6QQSR0GoJSyqsas1pQrUMWPGmD/++MOeRAFZt25d+zyecMIJBeUBVcBBtEXf9OnT7YSjiRMn2mU0ABAEgfu00gex9jdVKTqd1LpREQD9rBaNJrEwyze+1OV++eWX21apqHX5+OOP2xm7mnSkiUbaq/T777+3k0cUku7EJI17ayx80aJFBTV3x44dyxIYAIESuDDVOlN1F6qbUK0hzQLVzF4F7EMPPWRn+qrlSqDGh77AqDygxj/1vChE+/XrV3C5Wpnq1tXkI42T6l+dp/+n51CXaVKZglbuv//+Qv8fAIIgUN28opBUKTkt/NdyC7Vs3NaR9lBVgQB1FTKGGnuaLKQxUAVp7dq1bZm/cEGo50ihqedOgbp27VobqPq/GkudM2eOvd7w4cMZJwUQSIEM0zZt2pjLLrvM9OrVq1gBfU1KGjp0qJ2YRLWk2Fm8eLE58cQTbeH5/fff304CczeDL+l508xd/asg1TpTBeu3335rL1chjjvvvDOGRwwAsRO4bt5QzZo1s0snFKBaLuNaunSpncyiyS2aEeq2XosqWoUH/yvhVxZVJVJXuvsYqnXZtm1b+7u6ejWG7dLsa223JgpSPV+7d++2/0d73KqrXmtQ9aWI5wNAUAU6TFUB6fzzz7ddus8991yhQNUY6mOPPWZbQVdeeWVE23yhbDNmzLAzdl0KQ3dpixvGoYGcl5dnv+yE0lIZlQ7UrF9tGK8Z2QpgwhRAUAWumzeUWjpq1Vx66aV2NqkmtITSJJennnrKFlfXvqrwRo9jly5d7ISvilLRes241sSxs846y24MHtqSBYAgCnTLNLQy0o033mhbQdrmK9TKlSvNXXfdZWeRXn/99ZSmqwC1Ph988EHbdS768lL0cY6EZu2ec845ttdAk4+0Tlg9B0AsURAf8ZAUYepSCTu1RjW7VJNbXPrwfvfdd+1Y6tNPP22OOOII26pl+UzZ1BWrcWcteZGBAwcW6wEojdsFrNnW//rXv2yX/D/+8Q9bvD6SjcGBiqIgPuIpacJUwajJLtrzUmNwajkV7drVN1TN9NWOJKrGo1nBKE5jngsXLrRdsm+88YbdMs3dBUZjnOXZXEA7xKhbWGtRtXSpXr169ja1lAaIJQriI56SKkzdkyYcHXjggbZrsuh606+++squj1SoEqbGzqzVMhctUVHYaabu/PnzbQiGGj9+vF2KpC5frfGdNm1aRM+LZvdqBu9nn31mn5vXX3/dThwD4oGC+IiXpAlTl5bBaKcSzd5VcQCFQNEWqvZFveqqq2wgDB482K55TCUKxGeffda23lUwIbRL3FWnTh1bHEMnhaHWlIpm3KrsX0nLjcJ18ypIVVJQXe3u7QBAMkm6MHUpTPv27WsnuGic1C22Huqdd96xu5aoak+3bt1SoqWqINVEIhWlL1pVSl3f2plH/x566KElLlUpTwF6jbmquL2qI6XC4wsgNSVtmIq6EzV5pn79+nZvTLVUi3Jr+mpccNiwYUldKF9Bqsfj4Ycftr+PHDnSdO/e3Rx++OEFy1NC14jq+j///LPtqq3oJuwK5I8++shOOgKAZJXUYSrqjlRAqqCDxlDd2aVFW1fatWTAgAE2XFUIwi3yoECJtJhAIooO6PgioQLyb7/9th0zFt1HfWnQeTq51IWrmbfayFszozV2qvulrnNtIlBemrFLkAJIdkkXpuECTa0qrUNVoYAnn3zS7NixI+z/VZAqUFVN6ZJLLrFjhU2aNDFBp67WN99803z99df2d82sPf744wsu0xeNn376yZ60fChcYGsmdKNGjUz16tULaulq2zuXbkMTmdS1/ssvv5jGjRubDz/80Bx22GFUNgKQ9JIuTEuiVpgCUrN8tcRDLa+SaEarThpTVbiqKERQZ6Bqtq7GhBWk+qKhtZ4aE1XgafauunGLTkBSt65alGrVa03p6tWrbWtVgamQDDdmquUvqkKlwhn68qIgDepjhtREcQd4kVJhqgkwOqnwusYOy3rzvP/++3afTXWJav1q0Or7qku7T58+JicnxwapJlmpoPzcuXPtRt5uF7EmaWnCkVqbKj6/995729a7dt1xv3QoQHV7ClR9IQktAajau6q1q8t0Gx988IENXSAIKO6AaEiZMA2lHU806UjFG0procrGjRttt6/GGlXgXeOv6vqNdGlIomgikbap0xiogu/iiy82rVu3toUYVMZPQaovFSeccIK9PwrL9evX2yDURK3Q4vVu4FauXNmGrFqqbliqxavucD2OasHrC4i2ZAOCguIOiIaUDFO1Uk877TTz2muv2QIPv/76a5n/R+OJd9xxh3nmmWfsGlV1aep2/Biq6pLVLF2to1WXrQJV3bkaH9V5ClItgVFLVQHp0mVF94BViCpAdVLVIo0ra+xZJy0rUq1ddRMrlLUpuyocAUFDcQek9K4x0Wihqgu0ZcuWEf8ftcoUqtroWoUP1FXsZReVaFNLWgGnYNR4p74waMKRKDi1DEa0n2hR7du3L1Z43p2spZarblt0HQWzdn1RkJ566qm2RUqQAkhVKdkyDaXJOOPGjTOjR48utESkLKqiNHz4cHvSbdxzzz12vFBiXURfk4q0G45a1EuWLLH/uie1QNVdq4pDClTtzhJK3dRqZev/qzJRu3btCgJUE4Y0fqQZu2rdqvXpflFQt7FatLquijBofFS/a6xUJQIpWg8glaV8mIrCUC04zW7VhJ3QJR+R0P9TqGgSj1prXbt2jVq1H3WrajNtBZ9q5iowFeQlLe8RjVkqSHW/ilJXrcY4FYCqSqSTAlKTjvQF4K+//rK/a5xVJ5UV1HkKToWsfle9Y/2uFvDUqVPLtYsMEHTM+kU4hGkIhY/GF7Um9ZNPPjHltWDBAnvS0htN9lFFpeOOO67g8gYNGth/NS5ZUutVhRI0Fqlt5HRS7Vx1qRalsVC1JBXgWoqif3XSzjlqIbvrQUu6n2qZ6rY10UjBrElHboWookXuJTRIRceuLyDsR4pUwaxflKaSE2kJnQRS60wf5Oq+jNbWXSXdbXVrKjCeeuopu4enQqPozNaSJjCEKxjvUutNRQ/UYlXLUeO0mhH722+/2epLanXOmjXLbn1W9Nj0/1RFKDMz0/6sk5bpqOU4efLkiMYqdZ+KFqDQ33GXwOjY9a+K0et21ZWsk+67AlNh6m4YoMdEx6AiDX6tDJVKIn1/uNfTl7RU29whWrSOOtx8g6ItV60U0Mz5rKysuB0bEpszKdsyLekDXkGhlp7WlWpizSOPPGJry5ZFY4xlXa4uWo0tqqCBxiA1mUmTeorOoHVbmQpetXBLW9+qVmu4lmtRaq1q5nIkioazunK1LEg0S/iFF14oV7F7IFkw6xcl4ROxBApVjYNqxq+6Q6+//nrbavRKk3r07VbfftwAVtesZtK2bdvWLjFp2LCh3V9U45iRUAtSt6mJRZqQpC5gzTZWKT8vFPYvv/yyLXChv6Fv2xMnTiRIgQgwtppaYhKmaimpCzG0Sk5QaR2pAnXmzJm221ct1dIm/5RELVEFqLp13f+vtZ433XSTXfNa3q5QBbEC9+OPPzYrVqwoNs6pCUbqFu7UqVPEa2F1XCrCoCpImpikrmdX79697VKgZHhOgVhibDU1RT1MNf6nvTLVqlHLqGfPnnb5RdCpq/XOO++0pQU1nqqJSqHblZVEY46qDqQWqVp3otanJjlp3Kq8IaqwU81gBWlogCow9XjrpMlFup6+GS9fvty2sPU3S/rio/vi7u1atMtZk5XUtavQDy3wACA6FZX0vtMX37ICWreLFAlTTUhRcGp5iMb7ZsyYYZeNKFDVTRophURoUGhg2C+0ifbYsWNN//797axd7cZSGgWpOzFJrURVT1JLtyLUstUMYXeCkrqDFfIaV1W9XDfs3Dq7mkykN7SOUccdjiY9XXfddYVmHCt8zzzzTNtiDlo9YiAoY6u0YJNL1MJUH/CTJk2yH8QaZxN98Ct4NM6mbsnBgwdHdFsqoKAqQ36lVqDCSXuEapKS1nSqNehWCArlTg7S5uRDhgwpWFpSEWrdut3n2kquWbNmtlUZbjbvAQccYHr16mXGjBlj/09JrWj3i4qW2agGserr6k0OILZowSaXqIWpuitXrVplu3ddqsKjFqmWhehDX62nSy+9tMzbGjp0qBk4cGChD/yDDjrI+I1acSrVp64aFVJQkXidwk1UimapPYVpJJOLyjO+Wb9+fVv4HkD80IJNHlEJU7V8FKZaU6VydurudWvAKlBVVUjnqTWl8ClrQozWNJZWdMAvdD/c+6IvCiocr5OWwKi1qiUwf/zxR6IPE0CAxaIFGynGauMcpu4kGpWp0yQdjSU+9thjdmmHgla1XEeOHGmLBqgsnsbjko1a3259WrVCNXFHk4z04o5knSoAxKsFGyndnlYGBH3uRHocJnBFdQKSJr5o02lNQFKwaJzQHX/T2kdNlFGFifJyJ9wkYiJSpAWiNGGq6MxcFTZwW64a79Txa1xVj0VZVImp6O25Red1TO6sWy1nCVf+T0KLOeg64R4/9zY109i9nIIMweI+b2W9Vt3L9drx06Q+RMe+++5rvvnmG1tLOxrUElYwJ0Pjp2bNmnaHsNLmg7ifqRUuCujEwLRp05zq1as7Xbp0cV555RVn0aJFzpAhQ5wDDjjAWbFiRblvT/9Hh8qJx4DXQMmvgbLeW7yPeP/w/jFlPgYVySiJWW1eVQ3SJCJNzFFLR5NhNAnp6KOPLvdtqYWlMdcWLVrYAgXRqs/rV+6Eq1S4r6l2f2NxX/UWVo+H6iWXthZY7yNNEtQ8BrfXI8iPPcfOYx/N102k76O4lxPUZCRNwtFOJDpALdWo6JIL3TFN8BHd+aC96Ssqle5rqt3faN/XSIZP9D7S5grxOJ544th57KP1uqnIMGRcavMG+UUOAECkqA8HAECqhKnWnWpbtCCsP/Uqle5rqt1fv91Xvx1PeXDsPPZ+et0EYnNwAAD8LDAtUwAA/IowBQDAI8IUAACPCFMAADwiTAEASPUwdScjp8KkZO4jokFlBbWxAeJPG0toc4qgvm5CN89AEoXpDz/8YG655Raza9euYjusBN3KlSvNzJkzzdSpU83y5cvtebqPyfpi1l6MOiXb8+g3ixYtMr169TJnnHGG6d+/v/nqq69MkKg86Y8//mj3TQ5aKOnzSvsdz5o1q8Sdnvz8urn88stNp06dTL9+/WyddSRJmM6fP98cc8wxtpaiu2WYWm7J0HpbsGCB3Q9Ve8D26NHDdO3a1Vx//fUF9VWTLVC///57c+SRR5rp06ebZKcNG0aMGGGf14kTJ5rvvvsurn+7Xbt2tlXapk0b89///tfccMMNZuzYsSYoYaQPcwVSq1at7L7JQWlhL1y40LRv397WRm7atGmgimToy8tJJ51k9tprL3PuueeavLw8+9l03XXXmaBZunSpeeSRR8zNN99sJk+eHN0bdwJo/vz5Tq1atZxbbrnFSTbr1q1z/vGPfzg33nij/XnlypXOXXfd5bRs2dI555xzCq63e/duJxnMmzfPqVmzpnPrrbc6yW7hwoVO3bp1nfPPP9/p1KmTc+SRRzpHHXWUM2nSpJj/7T179jjDhg1zunXrVnDehg0bnLvvvtsew/333+/4/bGrX7++fc/r54ceesipVKmSk5eX5/jdpk2bnNNPP93p379/wXm5ubnO3LlzneXLlzt+tm3bNufSSy91rr/++oLztm7d6hx99NH28e/Ro4cTFN9//73TuHFj59RTT3XatWvnVK5c2XnggQeidvuBC1O9+PSmuvjii+3vu3btcu644w6nZ8+ezplnnum88847ztq1a52g0v1r3ry589VXXxWct3HjRmfKlCnO4Ycf7lx00UVOsvjpp5/sG/LOO+8seC6nTp1qf9f91YdNstB96927t3PZZZfZYJNvv/3WfkjVq1fPGT9+fMyP4fLLL3c6dOhQ6DwFqoLp2GOPdXJychw/WrNmjT3uG264oeA8PYZ6v+t9oteJn0NVgXTSSSc5c+bMsa+DM844w2nTpo2zzz77OCeccEJcnnsvFD6jRo0qCFIZPHiw869//cvJyspyHnzwQcfvli1b5jRr1swet9sQmTBhgtOwYUP7ORQNgevm1bia9qJTd6d2ldcu8B9//LEdP1F3r7rPnnnmGTvQH0Taa3Lnzp2FxrL23ntv07lzZzNs2DDbVaf7F3T6IvfFF1/Ynw877DD7r7rw7rnnHjNp0iQzdOhQ07dvX/PWW2+ZZKD7u2TJkkJ7iaorX3v+9unTx4waNcpuWRirv+1ui6huUb2GXDoe/X3tM/zkk0/68n2jx0vv82uuuabgvLvvvtvOKRgwYIA577zz7GvFfT35zbp16+xjnp+fbwYNGmTPGz9+vJkyZYrt+lW3/6uvvmr8Rq8bd8LUL7/8Yuem1KhRw/z222+2i/Scc86xe0y/++67xs80LKYx3mbNmtnPUHevUg11VKtWLXrDZk4AvfHGG7abTF1m+pb3559/2m98oi7RtLQ02xUURPoWq9aLvnWrWyLU5s2bnc6dOzvdu3d3koFa3G533YEHHmi/6S5evNhe9vXXX9supJNPPtlZvXq1kwwGDRpkX6+rVq0qdL7us3pa1Oug5zhWlixZ4qSnpzt9+vSxj724rWS17PQ8zJgxw/EjtaBdL7/8sj3WyZMnO3/99Zfz6aef2pae23ryGz3Ges9ee+21zrnnnuv85z//KbhsxYoVTnZ2tnP11VfbzzD3+fCTL774wnaJqndAPYAaYrvyyivtZQsWLLAt7B9//NGXx+7Sa2TIkCGFzlML9eCDD3Y+/vjjqPyNQIVp6Djha6+9ZkPnyy+/tL+HPpH6wHj00UedoNILVN0PGt/SB2Cohx9+2HatxPJDN57UbaT71L59e2f27NnFvjTVqFHDjpEnA3VhN23a1Hn88ccLwsz10ksv2Q+ppUuXxvQYPvroI6d69erONddcY7tPXb///rsdqw8dXvBzl913331X6DzNJzjvvPMcv1KXvp5ffQmYNm1aoctuvvlmG1R+DqNvvvnGhr5C9Iknnig4/6233nIyMzPt/A6/cRtYRbmPs/JE78f33nuv4LIPPvjANs4qIqabg0ebmuf6AqBuny5dutguBs2ME7frbNmyZebAAw+0lwWRuhxatmxpuzdPPfVU+7u6sk4++eSCmXWaEejOYA46dRtdeeWV9r5mZmba83Sf9Vw3bNjQHHLIIUmzwbxmZX/77bfm1ltvtfdbr+F69eoVdME2adIk5ksm9DrScquLLrrI/P7773ZmbOvWrW3X+p9//mmHUPxOj5NO7mtF3ZAaCtH98Ct16c+YMcN07NjRjBs3zr6uNYNdNKzTvHlz242qbkc/UpeoXiNFl659/vnn9n1ayWdL2n766Sfz9ttvm0suucQccMAB9jw3O3TSY633WpUqVQo+X9QFfN9999lliRXiBFBp3+BGjBhhZ75qFqyf6VtR0W9ObsvbPV8tNc20VEtUrQbNAq1du7adAZsK1C3atm1b5++//3aCLrRX5brrrrOTjjS7Vt/41VWpWaqHHnpoodZiLKll17FjR6dJkyb272rSmybIBNHIkSOdjIyMqE0kiXV3Y6NGjZzjjjvOueKKK2y3aZ06dWxvVJBoCGrAgAG+/Dz6+eef7ftLvQBDhw4N+57S+1G9Ynrt63NWkx7Vc6D3Y0X5ej9TTZbQN4dIvPHGG+ajjz4yL7zwgvnkk0/MUUcdZfy8APree+81q1evtpNvtHZLg/mh99n9V2u6tBZR902tBk1EOuKII0yQlOd5FHeS1XPPPWefSz+3OMIVFVALT/dXrSetzQv3ONx///32m/Ps2bNtL4peC++8846dCBQvGzZssMe7ceNG++09PT3dBIla2J9++qmdXPL+++/H9bHzQq/vnJwcW7xB73/1PKk3KijUotOkI02g0kRBP70/N2/ebNfkq8dCrelrr73WFvYZPHhw2Ne3eoTUy6e6BV9++aXtQagwx4dCJ5yU1e/tGjNmjG3FFJ204zcaqNc3UU1I0IC4WpxalqB1pa7t27fbf/08hhIpTa7RJKOik25Chd5PfUO/6qqr7Do2v33jLYuOXcfdqlUrOy6pyXBFX787d+4stAzqk08+sa0Vv/ek+NEPP/xg5xUsWrTICSK1joK6XlwTJbV+1m+2bNlix3RfeeUV+7smqamFql6u0Baq3pfqEdJncZUqVaKSG74LUy1mrlatWqHJBCUFqugBcfm9O7CshfN9+/YtdP0333zT+eOPP5ygiqS7JdwXBnVBakJMkERSVCCoH5x+tmPHjkQfAnxmU5GQV7Dqvaj3Zn5+fsGXWn0eaWa1vpRFg6/CVK0XVabQQuZDDjnE6dKlS6mBOnDgQLuc4tdff3WCoqyF86NHj7bnTZ8+3VbrGD58eCA/hPWC1hIM3V99Uwz37TCUKpH4dWlDNIoKaAmE67HHHnMmTpyYoKMFUsOukKVG7nIqfQb99ttvzk033eRceOGFUV0V4asw1fIABai6vl599VW7Big0UEO7yETXUctHD47fuU/q2LFjnRNPPNF294ZS1Sa1TPVlwu3mve2225xffvnFCaJIu1vc3gWtszz++OMLvjkGiY753nvvLTQBRhMadH/V46AvRVpf+vnnn9v7qjWRCtr169cn9LiBZLdnz56Cxog+i9TrqUpy6tqNdoU1X4WpviVo3ZIbnCopVzRQ3W8b7gNUdL2e30WycP7tt992UqW7Rc+nuucVMqWNq/pdJEUFbr/9dnu5xmf8XpMVSBZ79uwp+Iw95ZRTbAMsFnNrfBOm4boyNXVZC92LBuq4ceMKWnZBnKSTDAvno9XdoolXF1xwgZ3QkCxKKiqgUxBfr0AyfAbddNNN9rMnVkVgfLPy362XGEoL27VkRItsNb1Zi95VsEBbRqlWpPhtsXCqLJwvDy0H0Rc3TVfv3r27fc569uxpa9HqeVSN5SBtSeWlqEAQX69AMjjyyCPNnDlzYraUx1frTFWVIrSyj1uxYtu2bbYikIrY161b164p0/qgoNMTq0Lnqtqk+63Q0Zq5oKyXKy/3pabnVBWP5s2bZ9eRam/KZHbbbbeZ559/3nzwwQcFRf0BxJebJ7Him5apFrQrUBQsWoh92WWXFdxxtVD1QZSWlmbLVwW1VGBR+kKg1lmQF86Xh55PPc/aOUM7/ShMkzlIixYVIEiBxIl1r1Blv7RI1SpTkB5++OG22k8oVdvQh5JOyRKkLtWFPPjgg22oJHOQxrO7xS/0Wl2zZo39ApisvQ0AfNLN63btKkjVUrvwwgttKbnQ7t6tW7ea9evXm/333z+Rh4qAdLf4iYqY+7V4OYAkCdOiQaq6s9o0NzRI3R1EAADwq4SFqVv0u7QgBQAgCBLW5FOQLl++3I6fXXDBBWbChAkEKQAgkBLaMu3Xr58dO3v66acJUgBAYCV0zPTvv/82derUYUwUABBoCZ/NCwBA0DFNFgAAjwhTAAA8IkwBAPCIMAUAwCPCFAAAjwhTAAA8IkwBAPCIMAUAwCPCFAAAjwhTAAA8IkwBAPCIMAUAwCPCFAAA483/A60XwXDnPdVDAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sampler = hemcee.HamiltonianEnsembleSampler(total_chains, \n", " dim, \n", " log_prob, \n", " step_size = 0.01\n", " adapt_inital_step_size=False, # <- Not recommended, but can be done.\n", " adapt_step_size=False,\n", " adapt_length=False,\n", " )\n", "\n", "samples = sampler.run_mcmc(keys[1], \n", " initial_states, \n", " num_samples=10**4, \n", " warmup=10**5,\n", " thin_by=5,\n", " show_progress=True)\n", "\n", "try:\n", " tau = hemcee.autocorr.integrated_time(samples)\n", " print(f'Autocorrelation time: {tau}')\n", "except Exception as e:\n", " print(e)\n", "\n", "corner.corner(samples.reshape(-1, dim).__array__())\n", "plt.show()" ] } ], "metadata": { "kernelspec": { "display_name": "hemcee", "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.13.7" } }, "nbformat": 4, "nbformat_minor": 5 }