diff --git a/Lab.ipynb b/Lab.ipynb new file mode 100644 index 00000000..2d7849ef --- /dev/null +++ b/Lab.ipynb @@ -0,0 +1,683 @@ +{ + "cells": [ + { + "attachments": {}, + "cell_type": "markdown", + "id": "91129cb1", + "metadata": {}, + "source": [ + "# No-glue-code" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "896323ee", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "\u001b[32m\u001b[1m Activating\u001b[22m\u001b[39m project at `~/Cambdrige`\n" + ] + } + ], + "source": [ + "using Pkg\n", + "Pkg.activate(\"..\")" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "baed58e3", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "\u001b[36m\u001b[1m[ \u001b[22m\u001b[39m\u001b[36m\u001b[1mInfo: \u001b[22m\u001b[39mPrecompiling AdvancedHMC [0bf59076-c3b1-5ca4-86bd-e02cd72cde3d]\n", + "\u001b[36m\u001b[1m[ \u001b[22m\u001b[39m\u001b[36m\u001b[1mInfo: \u001b[22m\u001b[39mPrecompiling Turing [fce5fe82-541a-59a6-adf8-730c64b5f9a0]\n", + "WARNING: Method definition sample(Random.AbstractRNG, AbstractMCMC.AbstractModel, AbstractMCMC.AbstractSampler, AbstractMCMC.AbstractMCMCEnsemble, Integer, Integer) in module AbstractMCMC at /home/jaimerz/.julia/packages/AbstractMCMC/bE6VB/src/sample.jl:81 overwritten in module Inference at /home/jaimerz/Cambdrige/Turing.jl/src/inference/Inference.jl:210.\n", + " ** incremental compilation may be fatally broken for this module **\n", + "\n", + "WARNING: Method definition kwcall(Any, typeof(StatsBase.sample), Random.AbstractRNG, AbstractMCMC.AbstractModel, AbstractMCMC.AbstractSampler, AbstractMCMC.AbstractMCMCEnsemble, Integer, Integer) in module AbstractMCMC at /home/jaimerz/.julia/packages/AbstractMCMC/bE6VB/src/sample.jl:81 overwritten in module Inference at /home/jaimerz/Cambdrige/Turing.jl/src/inference/Inference.jl:210.\n", + " ** incremental compilation may be fatally broken for this module **\n", + "\n" + ] + } + ], + "source": [ + "using Random\n", + "using LinearAlgebra\n", + "using PyPlot\n", + "\n", + "#What we are tweaking\n", + "using Revise\n", + "using AdvancedHMC\n", + "using Turing\n", + "using DynamicPPL" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "3d76390f", + "metadata": {}, + "source": [ + "## Model" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "a7d6f81c", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "funnel (generic function with 2 methods)" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Just a simple Neal Funnel\n", + "d = 21\n", + "@model function funnel()\n", + " θ ~ Uniform(-1, 1) #Normal(0, 3)\n", + " z ~ MvNormal(zeros(d-1), exp(θ)*I)\n", + " x ~ MvNormal(z, I)\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "5f408f2b", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Model{typeof(funnel), (), (), (), Tuple{}, Tuple{}, ConditionContext{NamedTuple{(:x,), Tuple{Vector{Float64}}}, DefaultContext}}(funnel, NamedTuple(), NamedTuple(), ConditionContext((x = [1.2142074831535152, 1.23371919965455, -0.8480146960461767, 0.1600994648479841, 1.9180385508479283, -3.401523464506408, -0.0957684186471088, 0.6734622629464286, -3.2749467689509633, -1.6760091758453226, 1.9567202902549736, 0.1136169088905351, 0.11117896909388916, -0.5373922347882832, -0.12436857036298687, -1.2901071061088532, 1.702584517514787, -0.44460133117954226, 1.0818722439221686, 1.2208011493237483],), DefaultContext()))" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "Random.seed!(1)\n", + "(;x) = rand(funnel() | (θ=0,))\n", + "funnel_model = funnel() | (;x)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "d852c160", + "metadata": {}, + "source": [ + "## Sampling" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "486d475d", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "AdvancedHMC.HMCSampler{Nothing, Nothing, Nothing, Nothing}(AdvancedHMC.NUTS_alg(500, 0.95, 10, 1000.0, 0.1), nothing, nothing, nothing, nothing)" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "nadapts=500 \n", + "TAP=0.95\n", + "ϵ=0.1\n", + "nuts = AdvancedHMC.NUTS(nadapts, TAP; ϵ=ϵ)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "9e114ad8", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "AdvancedHMC.HMCSampler{Nothing, Nothing, Nothing, Nothing}(AdvancedHMC.HMC_alg(0.1, 20), nothing, nothing, nothing, nothing)" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ϵ=0.1\n", + "n_leapfrog=20\n", + "hmc = AdvancedHMC.HMC(ϵ, n_leapfrog)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "1f729dc6", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "AdvancedHMC.HMCSampler{Nothing, Nothing, Nothing, Nothing}(AdvancedHMC.HMCDA_alg(500, 0.95, 1.0, 0.1), nothing, nothing, nothing, nothing)" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "n_adapts = 500\n", + "TAP = 0.95\n", + "λ = 0.1 * 10\n", + "ϵ=0.1\n", + "hmcda = AdvancedHMC.HMCDA(n_adapts, TAP, λ; ϵ=ϵ)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "b0193663", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "\u001b[33m\u001b[1m┌ \u001b[22m\u001b[39m\u001b[33m\u001b[1mWarning: \u001b[22m\u001b[39m[DynamicPPL] attempt to link a linked vi\n", + "\u001b[33m\u001b[1m└ \u001b[22m\u001b[39m\u001b[90m@ DynamicPPL ~/.julia/packages/DynamicPPL/jjVG9/src/varinfo.jl:791\u001b[39m\n", + "\u001b[33m\u001b[1m┌ \u001b[22m\u001b[39m\u001b[33m\u001b[1mWarning: \u001b[22m\u001b[39m[DynamicPPL] attempt to link a linked vi\n", + "\u001b[33m\u001b[1m└ \u001b[22m\u001b[39m\u001b[90m@ DynamicPPL ~/.julia/packages/DynamicPPL/jjVG9/src/varinfo.jl:791\u001b[39m\n", + "\u001b[32mSampling: 100%|█████████████████████████████████████████| Time: 0:00:02\u001b[39m\n" + ] + }, + { + "data": { + "text/plain": [ + "Chains MCMC chain (5000×34×1 Array{Real, 3}):\n", + "\n", + "Iterations = 1:1:5000\n", + "Number of chains = 1\n", + "Samples per chain = 5000\n", + "parameters = param_1, param_2, param_3, param_4, param_5, param_6, param_7, param_8, param_9, param_10, param_11, param_12, param_13, param_14, param_15, param_16, param_17, param_18, param_19, param_20, param_21\n", + "internals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size, is_adapt\n", + "\n", + "Summary Statistics\n", + " \u001b[1m parameters \u001b[0m \u001b[1m mean \u001b[0m \u001b[1m std \u001b[0m \u001b[1m mcse \u001b[0m \u001b[1m ess_bulk \u001b[0m \u001b[1m rhat \u001b[0m \u001b[1m ess_per_se\u001b[0m ⋯\n", + " \u001b[90m Symbol \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Real \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Missin\u001b[0m ⋯\n", + "\n", + " param_1 0.1027 0.4682 0.0125 1316.8261 1.0006 missin ⋯\n", + " param_2 0.6380 0.7443 0.0088 7305.3358 1.0007 missin ⋯\n", + " param_3 0.6571 0.7388 0.0087 7222.3134 0.9999 missin ⋯\n", + " param_4 -0.4590 0.7424 0.0081 8600.6777 0.9998 missin ⋯\n", + " param_5 0.0827 0.7254 0.0078 8658.7613 1.0009 missin ⋯\n", + " param_6 1.0204 0.7597 0.0109 4919.8215 0.9999 missin ⋯\n", + " param_7 -1.7932 0.8261 0.0145 3273.3659 1.0001 missin ⋯\n", + " param_8 -0.0484 0.7195 0.0071 10192.8327 1.0002 missin ⋯\n", + " param_9 0.3575 0.7262 0.0076 9149.6800 1.0002 missin ⋯\n", + " param_10 -1.7292 0.8133 0.0135 3701.3245 0.9999 missin ⋯\n", + " param_11 -0.8752 0.7379 0.0093 6376.3368 1.0004 missin ⋯\n", + " param_12 1.0242 0.7599 0.0103 5479.1056 1.0000 missin ⋯\n", + " param_13 0.0675 0.7458 0.0079 8945.0993 1.0009 missin ⋯\n", + " param_14 0.0668 0.7140 0.0072 9814.6348 1.0006 missin ⋯\n", + " param_15 -0.2908 0.7255 0.0076 9112.4223 0.9998 missin ⋯\n", + " param_16 -0.0508 0.7068 0.0070 10008.6090 1.0001 missin ⋯\n", + " param_17 -0.6693 0.7322 0.0087 7073.7412 0.9999 missin ⋯\n", + " param_18 0.8904 0.7460 0.0093 6393.5556 1.0004 missin ⋯\n", + " param_19 -0.2438 0.7394 0.0079 8715.5189 1.0000 missin ⋯\n", + " param_20 0.5602 0.7217 0.0082 7751.5157 1.0000 missin ⋯\n", + " param_21 0.6376 0.7380 0.0084 7807.3097 1.0011 missin ⋯\n", + "\u001b[36m 1 column omitted\u001b[0m\n", + "\n", + "Quantiles\n", + " \u001b[1m parameters \u001b[0m \u001b[1m 2.5% \u001b[0m \u001b[1m 25.0% \u001b[0m \u001b[1m 50.0% \u001b[0m \u001b[1m 75.0% \u001b[0m \u001b[1m 97.5% \u001b[0m\n", + " \u001b[90m Symbol \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m\n", + "\n", + " param_1 -0.8352 -0.2336 0.1326 0.4584 0.9071\n", + " param_2 -0.7920 0.1390 0.6151 1.1211 2.1535\n", + " param_3 -0.7435 0.1493 0.6307 1.1539 2.1429\n", + " param_4 -1.9727 -0.9420 -0.4536 0.0433 0.9569\n", + " param_5 -1.3355 -0.4084 0.0832 0.5671 1.4991\n", + " param_6 -0.3763 0.4823 1.0017 1.5315 2.5668\n", + " param_7 -3.4720 -2.3401 -1.7762 -1.2272 -0.2403\n", + " param_8 -1.4292 -0.5439 -0.0520 0.4395 1.3675\n", + " param_9 -1.0777 -0.1229 0.3547 0.8448 1.7903\n", + " param_10 -3.4370 -2.2589 -1.6796 -1.1744 -0.2281\n", + " param_11 -2.4021 -1.3726 -0.8447 -0.3686 0.5171\n", + " param_12 -0.4100 0.5065 0.9956 1.5327 2.5705\n", + " param_13 -1.4140 -0.4160 0.0706 0.5537 1.5270\n", + " param_14 -1.3651 -0.4031 0.0653 0.5342 1.4844\n", + " param_15 -1.7440 -0.7812 -0.2779 0.1959 1.0957\n", + " param_16 -1.3863 -0.5423 -0.0520 0.4442 1.3074\n", + " param_17 -2.1487 -1.1499 -0.6642 -0.1710 0.6959\n", + " param_18 -0.5586 0.3798 0.8693 1.3917 2.4085\n", + " param_19 -1.7016 -0.7273 -0.2266 0.2350 1.2082\n", + " param_20 -0.8251 0.0794 0.5603 1.0190 1.9974\n", + " param_21 -0.7633 0.1377 0.6239 1.1340 2.1128\n" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "\u001b[33m\u001b[1m┌ \u001b[22m\u001b[39m\u001b[33m\u001b[1mWarning: \u001b[22m\u001b[39mTail ESS calculation failed: OverflowError(\"4750 * 4503599627370496 overflowed for type Int64\")\n", + "\u001b[33m\u001b[1m└ \u001b[22m\u001b[39m\u001b[90m@ MCMCChains ~/.julia/packages/MCMCChains/OVsxE/src/stats.jl:319\u001b[39m\n" + ] + } + ], + "source": [ + "nuts_samples = sample(funnel_model, nuts, 5000)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "f610b909", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "\u001b[33m\u001b[1m┌ \u001b[22m\u001b[39m\u001b[33m\u001b[1mWarning: \u001b[22m\u001b[39m[DynamicPPL] attempt to link a linked vi\n", + "\u001b[33m\u001b[1m└ \u001b[22m\u001b[39m\u001b[90m@ DynamicPPL ~/.julia/packages/DynamicPPL/jjVG9/src/varinfo.jl:791\u001b[39m\n", + "\u001b[33m\u001b[1m┌ \u001b[22m\u001b[39m\u001b[33m\u001b[1mWarning: \u001b[22m\u001b[39m[DynamicPPL] attempt to link a linked vi\n", + "\u001b[33m\u001b[1m└ \u001b[22m\u001b[39m\u001b[90m@ DynamicPPL ~/.julia/packages/DynamicPPL/jjVG9/src/varinfo.jl:791\u001b[39m\n", + "\u001b[32mSampling: 100%|█████████████████████████████████████████| Time: 0:00:02\u001b[39m\n" + ] + }, + { + "data": { + "text/plain": [ + "Chains MCMC chain (5000×32×1 Array{Real, 3}):\n", + "\n", + "Iterations = 1:1:5000\n", + "Number of chains = 1\n", + "Samples per chain = 5000\n", + "parameters = param_1, param_2, param_3, param_4, param_5, param_6, param_7, param_8, param_9, param_10, param_11, param_12, param_13, param_14, param_15, param_16, param_17, param_18, param_19, param_20, param_21\n", + "internals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, numerical_error, step_size, nom_step_size, is_adapt\n", + "\n", + "Summary Statistics\n", + " \u001b[1m parameters \u001b[0m \u001b[1m mean \u001b[0m \u001b[1m std \u001b[0m \u001b[1m mcse \u001b[0m \u001b[1m ess_bulk \u001b[0m \u001b[1m rhat \u001b[0m \u001b[1m ess_per_se\u001b[0m ⋯\n", + " \u001b[90m Symbol \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Real \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Missin\u001b[0m ⋯\n", + "\n", + " param_1 0.1116 0.4844 0.0126 1412.2510 1.0030 missin ⋯\n", + " param_2 0.6409 0.7630 0.0056 18494.8500 1.0003 missin ⋯\n", + " param_3 0.6563 0.7341 0.0054 18494.8500 1.0023 missin ⋯\n", + " param_4 -0.4489 0.7738 0.0057 18494.8500 1.0013 missin ⋯\n", + " param_5 0.0916 0.7387 0.0054 18494.8500 1.0008 missin ⋯\n", + " param_6 1.0122 0.7602 0.0068 13709.0981 1.0030 missin ⋯\n", + " param_7 -1.7991 0.8076 0.0124 4323.3788 1.0009 missin ⋯\n", + " param_8 -0.0475 0.7271 0.0053 18494.8500 1.0059 missin ⋯\n", + " param_9 0.3593 0.7176 0.0053 18494.8500 0.9999 missin ⋯\n", + " param_10 -1.7389 0.8314 0.0122 4786.2571 1.0019 missin ⋯\n", + " param_11 -0.8884 0.7405 0.0064 17067.3833 1.0013 missin ⋯\n", + " param_12 1.0324 0.7586 0.0068 12775.6485 1.0027 missin ⋯\n", + " param_13 0.0612 0.7115 0.0052 18494.8500 1.0026 missin ⋯\n", + " param_14 0.0576 0.7049 0.0052 18494.8500 1.0025 missin ⋯\n", + " param_15 -0.2848 0.7059 0.0052 18494.8500 0.9999 missin ⋯\n", + " param_16 -0.0663 0.7493 0.0055 18494.8500 1.0001 missin ⋯\n", + " param_17 -0.6799 0.7329 0.0054 18494.8500 1.0002 missin ⋯\n", + " param_18 0.9009 0.7595 0.0060 16083.8415 1.0022 missin ⋯\n", + " param_19 -0.2384 0.7235 0.0053 18494.8500 0.9999 missin ⋯\n", + " param_20 0.5663 0.7420 0.0055 18494.8500 1.0001 missin ⋯\n", + " param_21 0.6437 0.7433 0.0055 18494.8500 1.0003 missin ⋯\n", + "\u001b[36m 1 column omitted\u001b[0m\n", + "\n", + "Quantiles\n", + " \u001b[1m parameters \u001b[0m \u001b[1m 2.5% \u001b[0m \u001b[1m 25.0% \u001b[0m \u001b[1m 50.0% \u001b[0m \u001b[1m 75.0% \u001b[0m \u001b[1m 97.5% \u001b[0m\n", + " \u001b[90m Symbol \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m\n", + "\n", + " param_1 -0.8729 -0.2411 0.1414 0.4873 0.9276\n", + " param_2 -0.7746 0.1213 0.6172 1.1341 2.1519\n", + " param_3 -0.7742 0.1636 0.6370 1.1344 2.1403\n", + " param_4 -1.9930 -0.9673 -0.4454 0.0924 1.0236\n", + " param_5 -1.3644 -0.4021 0.0955 0.5800 1.5932\n", + " param_6 -0.4151 0.4951 0.9882 1.5015 2.5778\n", + " param_7 -3.4943 -2.3275 -1.7638 -1.2414 -0.2990\n", + " param_8 -1.4757 -0.5401 -0.0424 0.4405 1.3866\n", + " param_9 -1.0262 -0.1276 0.3563 0.8391 1.8048\n", + " param_10 -3.4816 -2.2922 -1.6942 -1.1709 -0.2376\n", + " param_11 -2.4214 -1.3706 -0.8625 -0.3788 0.5300\n", + " param_12 -0.4144 0.5254 1.0030 1.5337 2.5786\n", + " param_13 -1.3274 -0.4277 0.0578 0.5478 1.4726\n", + " param_14 -1.3147 -0.4071 0.0520 0.5357 1.4133\n", + " param_15 -1.7091 -0.7450 -0.2665 0.1876 1.0607\n", + " param_16 -1.5507 -0.5647 -0.0675 0.4274 1.4156\n", + " param_17 -2.1845 -1.1587 -0.6694 -0.1713 0.6950\n", + " param_18 -0.5178 0.3903 0.8748 1.4069 2.4258\n", + " param_19 -1.6924 -0.6976 -0.2310 0.2270 1.1589\n", + " param_20 -0.8190 0.0547 0.5392 1.0695 2.0687\n", + " param_21 -0.8290 0.1653 0.6314 1.1214 2.1541\n" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "\u001b[33m\u001b[1m┌ \u001b[22m\u001b[39m\u001b[33m\u001b[1mWarning: \u001b[22m\u001b[39mTail ESS calculation failed: OverflowError(\"4750 * 4503599627370496 overflowed for type Int64\")\n", + "\u001b[33m\u001b[1m└ \u001b[22m\u001b[39m\u001b[90m@ MCMCChains ~/.julia/packages/MCMCChains/OVsxE/src/stats.jl:319\u001b[39m\n" + ] + } + ], + "source": [ + "hmc_samples = sample(funnel_model, hmc, 5000)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "88df45a3", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "\u001b[33m\u001b[1m┌ \u001b[22m\u001b[39m\u001b[33m\u001b[1mWarning: \u001b[22m\u001b[39m[DynamicPPL] attempt to link a linked vi\n", + "\u001b[33m\u001b[1m└ \u001b[22m\u001b[39m\u001b[90m@ DynamicPPL ~/.julia/packages/DynamicPPL/jjVG9/src/varinfo.jl:791\u001b[39m\n", + "\u001b[33m\u001b[1m┌ \u001b[22m\u001b[39m\u001b[33m\u001b[1mWarning: \u001b[22m\u001b[39m[DynamicPPL] attempt to link a linked vi\n", + "\u001b[33m\u001b[1m└ \u001b[22m\u001b[39m\u001b[90m@ DynamicPPL ~/.julia/packages/DynamicPPL/jjVG9/src/varinfo.jl:791\u001b[39m\n", + "\u001b[32mSampling: 100%|█████████████████████████████████████████| Time: 0:00:01\u001b[39m\n" + ] + }, + { + "data": { + "text/plain": [ + "Chains MCMC chain (5000×32×1 Array{Real, 3}):\n", + "\n", + "Iterations = 1:1:5000\n", + "Number of chains = 1\n", + "Samples per chain = 5000\n", + "parameters = param_1, param_2, param_3, param_4, param_5, param_6, param_7, param_8, param_9, param_10, param_11, param_12, param_13, param_14, param_15, param_16, param_17, param_18, param_19, param_20, param_21\n", + "internals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, numerical_error, step_size, nom_step_size, is_adapt\n", + "\n", + "Summary Statistics\n", + " \u001b[1m parameters \u001b[0m \u001b[1m mean \u001b[0m \u001b[1m std \u001b[0m \u001b[1m mcse \u001b[0m \u001b[1m ess_bulk \u001b[0m \u001b[1m rhat \u001b[0m \u001b[1m ess_per_sec\u001b[0m ⋯\n", + " \u001b[90m Symbol \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Real \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Missing\u001b[0m ⋯\n", + "\n", + " param_1 0.0979 0.4865 0.0229 427.6675 1.0077 missing ⋯\n", + " param_2 0.6547 0.7415 0.0160 2189.7809 1.0004 missing ⋯\n", + " param_3 0.6347 0.7416 0.0140 2846.6874 1.0009 missing ⋯\n", + " param_4 -0.4482 0.7324 0.0148 2459.9117 1.0002 missing ⋯\n", + " param_5 0.0916 0.7201 0.0128 3150.8292 1.0022 missing ⋯\n", + " param_6 0.9939 0.7645 0.0163 2285.0805 1.0002 missing ⋯\n", + " param_7 -1.7991 0.8208 0.0261 1001.8156 1.0031 missing ⋯\n", + " param_8 -0.0504 0.7234 0.0136 2815.2275 1.0008 missing ⋯\n", + " param_9 0.3700 0.7229 0.0132 3028.1210 0.9998 missing ⋯\n", + " param_10 -1.7251 0.8101 0.0261 966.5697 1.0029 missing ⋯\n", + " param_11 -0.8600 0.7541 0.0168 2021.1769 1.0020 missing ⋯\n", + " param_12 1.0075 0.7484 0.0167 2050.6918 1.0005 missing ⋯\n", + " param_13 0.0569 0.7187 0.0117 3750.8085 1.0008 missing ⋯\n", + " param_14 0.0608 0.7254 0.0134 2916.2452 1.0003 missing ⋯\n", + " param_15 -0.2655 0.7254 0.0126 3303.5375 1.0016 missing ⋯\n", + " param_16 -0.0366 0.7243 0.0128 3216.3677 1.0016 missing ⋯\n", + " param_17 -0.6590 0.7431 0.0154 2371.9178 1.0009 missing ⋯\n", + " param_18 0.8751 0.7536 0.0160 2242.7235 1.0004 missing ⋯\n", + " param_19 -0.2233 0.7202 0.0123 3419.9118 1.0002 missing ⋯\n", + " param_20 0.6038 0.7478 0.0142 2803.1610 1.0011 missing ⋯\n", + " param_21 0.6409 0.7377 0.0137 2922.8470 1.0005 missing ⋯\n", + "\n", + "Quantiles\n", + " \u001b[1m parameters \u001b[0m \u001b[1m 2.5% \u001b[0m \u001b[1m 25.0% \u001b[0m \u001b[1m 50.0% \u001b[0m \u001b[1m 75.0% \u001b[0m \u001b[1m 97.5% \u001b[0m\n", + " \u001b[90m Symbol \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m \u001b[90m Float64 \u001b[0m\n", + "\n", + " param_1 -0.8759 -0.2497 0.1436 0.4725 0.9130\n", + " param_2 -0.7777 0.1435 0.6347 1.1346 2.1667\n", + " param_3 -0.7896 0.1384 0.6153 1.1279 2.1692\n", + " param_4 -1.9185 -0.9338 -0.4423 0.0496 0.9832\n", + " param_5 -1.3330 -0.3826 0.0886 0.5713 1.4915\n", + " param_6 -0.4397 0.4663 0.9664 1.4970 2.5635\n", + " param_7 -3.4716 -2.3299 -1.7589 -1.2145 -0.2936\n", + " param_8 -1.4562 -0.5463 -0.0707 0.4393 1.3843\n", + " param_9 -1.0222 -0.1147 0.3627 0.8514 1.8522\n", + " param_10 -3.3582 -2.2815 -1.6821 -1.1519 -0.2374\n", + " param_11 -2.3854 -1.3465 -0.8462 -0.3597 0.6050\n", + " param_12 -0.4173 0.4949 0.9801 1.4995 2.5221\n", + " param_13 -1.3876 -0.4168 0.0545 0.5379 1.4619\n", + " param_14 -1.3516 -0.4284 0.0526 0.5433 1.4733\n", + " param_15 -1.7321 -0.7393 -0.2599 0.2137 1.1228\n", + " param_16 -1.4597 -0.5141 -0.0427 0.4371 1.4198\n", + " param_17 -2.1839 -1.1502 -0.6285 -0.1511 0.7155\n", + " param_18 -0.6034 0.3688 0.8616 1.3863 2.3647\n", + " param_19 -1.6165 -0.7083 -0.2238 0.2560 1.1999\n", + " param_20 -0.7926 0.0770 0.5904 1.0971 2.1209\n", + " param_21 -0.7719 0.1303 0.6252 1.1271 2.1344\n" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "\u001b[33m\u001b[1m┌ \u001b[22m\u001b[39m\u001b[33m\u001b[1mWarning: \u001b[22m\u001b[39mTail ESS calculation failed: OverflowError(\"4750 * 4503599627370496 overflowed for type Int64\")\n", + "\u001b[33m\u001b[1m└ \u001b[22m\u001b[39m\u001b[90m@ MCMCChains ~/.julia/packages/MCMCChains/OVsxE/src/stats.jl:319\u001b[39m\n" + ] + } + ], + "source": [ + "hmcda_samples = sample(funnel_model, hmcda, 5000)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "bbf0131e", + "metadata": {}, + "source": [ + "### Plotting" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "9c61e0ab", + "metadata": {}, + "outputs": [], + "source": [ + "theta_nuts = Vector(nuts_samples[\"param_1\"][:, 1])\n", + "x10_nuts =Vector(nuts_samples[\"param_11\"][:, 1]);" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "0b0923f1", + "metadata": {}, + "outputs": [], + "source": [ + "theta_hmc = Vector(hmc_samples[\"param_1\"][:, 1])\n", + "x10_hmc =Vector(hmc_samples[\"param_11\"][:, 1]);" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "fec8ace5", + "metadata": {}, + "outputs": [], + "source": [ + "theta_hmcda = Vector(hmcda_samples[\"param_1\"][:, 1])\n", + "x10_hmcda =Vector(hmcda_samples[\"param_11\"][:, 1]);" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "8869229b", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAArEAAAL3CAYAAACd2x1cAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABY2klEQVR4nO3dd3wVVf7/8fcNKYQQQgu9R6z0IoIFFBGEVbGyq0gTdnXFhg0URSxfrNhWxYIglrWzCCjKIkRUZFXASkDQCIQOkgSEJCTz+4NfZs5N7k1uQpKbk7yej0ceTuZOOXcSru/M+cw5PsdxHAEAAAAWiQh3AwAAAICSIsQCAADAOoRYAAAAWIcQCwAAAOsQYgEAAGAdQiwAAACsQ4gFAACAdQixAAAAsA4hFgAAANYhxFZjnTt3ls/nU0xMjPbs2XNUx/L5fPL5fGXUsuonNTVVPp9Pbdq0KfRa/rVdtmxZkcdo06aNfD6fZs+e7bd+9uzZ7jGio6O1c+fOoMfIyspSgwYN3O3vv//+oNseOHBATz31lAYNGqRmzZopJiZGtWvX1nHHHafhw4dr3rx5ysvL89unX79+8vl8Sk1NLfK9HI386+Dz+fTuu+8G3e7ss88OeL1GjRoln8+nUaNGFXme/Oua/zNbtmyZe96SfN1zzz3uMfPy8jR79mwNGDBAjRo1UlRUlOrXr69jjz1W559/vh5++OFyvXamnJwcLVmyRLfeeqt69uypunXrKioqSk2aNNH555+vhQsXBt138+bNev755/X3v/9d3bt3V0xMjHw+n8aOHXtUbbrnnnsKXb+aNWuqUaNG6ty5s0aNGqXXX39dhw4dKtXxzZ9hkyZNdODAgYDbbdmypVJ85uVfD/N3SPLeR3G/w0fD/FwJ9lW3bt1yO39lk/+5U1H/PnFEZLgbgPD4+uuv9f3330uSsrOz9dprr+mGG24Ic6tQ3nJycvTqq6/q5ptvDvj63LlztXfv3mKP88knn2j48OHatWuXIiMj1b17d51++uk6fPiwNm7cqNdff12vv/66evbsqf/9739l/TZCduedd2ro0KGKjCz/j7omTZpo5MiRhdavWbNG3333nRo3bqxBgwYVer1Lly6SjvxRcN5552np0qWSpG7duumMM85QjRo19Ouvv2rRokWaP3++atWqpfHjx5fre5Gk5ORkDRgwQNKR93baaacpLi5OP//8s+bPn6/58+fr73//u2bMmFEozL333nu66aabyq1t5rXMzc1Venq6UlJS9Morr+iVV17RjTfeqKefflp//etfS32OHTt26LHHHtPdd99dVs2ukuLi4nTJJZcEfK1WrVoV3BpUN4TYamrmzJmSpObNmystLU0zZ84kxFZxnTp10tq1azVr1qygIfbll1+WJPXs2VNff/11wG0WLlyoCy64QLm5uRozZoymTZumRo0a+W2zadMm/d///Z/efvvtsn0TJVCrVi2tX79eL730kq6++upyP9/xxx9f6K6udORu2XfffRf0dXO7pUuXqlmzZvroo4/UqVMnv9fT09P13nvvqWnTpmXc8sAiIiJ08cUX64YbbtDpp5/u99pbb72lK664Qi+88IJOPfVUjRgxwu/1tm3b6rrrrlO3bt3UrVs3vf3223rggQfKrG3BruXGjRt1zz336LXXXtPf/vY37d27V//85z9LfPzY2FgdOnRIjz76qK655holJiaWQaurpoYNGxb5ew2UJ8oJqqE///xT//73vyVJr776qmrXrq0ffvghaGhB1ZCYmKjzzjtPP/30k1auXFno9U2bNmnJkiXq1auXTjzxxIDH2LNnj4YPH67c3Fxdf/31mjlzZqEAK0mtWrXSjBkz9J///Kes30bI8v8ou/fee/Xnn3+GrR2hevPNNyVJU6ZMKRRgJSkhIUFjxozRueeeWyHtOeuss/Tuu+8WCrCSNGzYMLeres6cOYVev+CCC/TUU09p1KhR6tSpU4XcCZekpKQkvfrqq7r11lslHfkd+PXXX0t8nGbNmumSSy5RZmZmkSU1AMKLEFsNvfPOO8rIyFCHDh105plnatiwYZK8u7PBrFixQueee67q1q2r2rVrq0ePHu6du4JSUlLk8/lUr169IuvTevToIZ/Pp3nz5rnrfv75Z02ZMkWnnnqqmjdvrujoaDVo0EBnn3120Dt7+TVg/fr1U05Ojh566CGddNJJio2NVYMGDXTRRRdp7dq1Qdvxxx9/6N5771WPHj2UkJCg2NhYtWvXTpdddpk++uijQtsfPnxYL730kvr166f69esrJiZGbdu21TXXXKPNmzcHPc+CBQvUt29fxcfHKyEhQaeffrrfey9vY8aMkaSAP7dZs2YpLy/P3SaQf/3rX9q3b58aNWqkhx9+uNjznXHGGSG1Kz09XZMnT1bHjh0VFxenmJgYNWvWTKeeeqruvvtu5eTkhHQc0+DBg9W3b19t27ZNjz/+eIn3r2g7duyQpIB/FFRGXbt2laQif9/D5YEHHlCzZs10+PDhUv/sH3jgAUVGRmrGjBn67bffSrRvaT4f3n//fY0dO1YdOnRQvXr1VLNmTbVt21ZjxozRunXrSvUeAvnvf/+r8847T40bN1ZUVJTq1aun9u3ba/jw4frss8/K7Dymomr+8wWrKTXXL126VOecc47q1aun2NhYdevWLeAfUZJXf79s2TKtWbNGF110kRo2bKiYmBideOKJeuyxx+Q4TtD2LFmyRBdddJGaNm2q6OhoNWrUSBdeeKFWrFhRmkuAckKIrYbyw2p+WMn/75tvvqmDBw8G3Oedd97R6aefrkWLFqlly5Y6//zzFRsbq7Fjxwbsmj7++OPVu3dv7du3L+jduB9++EHffvutGjdurCFDhrjrp0+frnvvvVd79+5Vx44dddFFF+m4447T0qVLNWzYME2YMCHoe8vJydHgwYN17733qlWrVhoyZIji4uI0d+5c9enTJ2DR/XfffaeOHTtqypQp2rBhg0477TRdcMEFatKkiRYsWKCHHnrIb/vMzEwNGDBA48aN07fffqtOnTrp/PPPV0xMjGbMmKGuXbtq9erVhc7z+OOP67zzztNnn32mE088UUOGDNGhQ4c0dOhQPf3000HfU1nKfwir4M/acRzNmjVLtWrVKrKOMD9wDxs2TDExMWXSpj///FOnnXaaHnjgAe3YsUP9+/d3f+a//vqr7rvvvqAP2BQn/2f38MMPH/XDi+WtVatWkqQZM2YoKysrzK0p3i+//CJJFVbeUBJRUVHuH+eLFy8u1THat2+vcePGKTs7W5MnTw55v9J+Plx22WX697//rdjYWJ111lkaOHCgIiIiNGvWLHXv3l1ffvllqd6H6ZVXXtE555yjhQsXqm3btrr44ot1xhlnqE6dOnrzzTf1/vvvH/U5ysvLL7+s/v37a+/evRo0aJC6dOmi1atXa+TIkXriiSeC7vfxxx+rV69eSklJ0YABA9S7d2+tX79et9xyS9C67VtuuUVnn3225s2bp1atWmno0KFq166d5s2bp9NPP12zZs0qp3eJEnNQraxbt86R5ERFRTk7d+501x9//PGOJGfOnDmF9tm2bZsTHx/vSHKmT5/u99p///tfp2bNmo4kp+Cv04svvuhIcgYOHBiwLTfddJMjybn55pv91i9btszZuHFjoe1TUlKcFi1aOJKclStX+r22dOlStw1du3Z1tm3b5r528OBBZ+DAgY4k5+9//7vffvv373datmzpSHJGjBjhZGZm+r2+b98+Z/HixX7rLr/8ckeS85e//MXZsWOH32uPP/64I8lp3769c/jwYXf9d99959SoUcOJiIhw3nnnHb99XnvtNcfn8zmSnNatWxd63/nva+nSpYVeM7Vu3dqR5MyaNctv/axZsxxJTv/+/R3HcZxJkyYV+lkvXrzYvQaO4zgjR450JDn33Xefu01OTo4TERER9PektF555RVHknPuuec62dnZfq/l5uY6y5Ytc7KyskI+Xv51WL58ueM4jnPRRRc5kpybbrrJb7v+/fsHvF75733kyJFFnif/ugb6mZmmTJniSHL69u1b5Hb5vzuSnMaNGzvjxo1zZs6c6axatcrvd6ky2LZtm5OQkOBIcp566qlit8+/BlddddVRnTfUa+k4R/5d5V/PnJyckI6f/zmSlJTkOM6R9xkXF+f4fD5n9erV7nabN28O+JnnOKX7fHAcx3nzzTed/fv3+63Ly8tznnnmGUeSc9JJJzl5eXl+r+dfjylTpoT0/tq2bev3b8O0Y8cOZ9WqVSEdx3FC//13HMf57bffit02/9/tb7/9FnB9VFSUM3/+/IBtSEhIcP7880+/1/r27ev+jGbMmOH32pIlSxyfz+fUqFHD2bx5s99rL7zwgiPJOeaYY5zvvvvO77Xk5GQnPj7eiY6OdtavXx9S+1G+CLHVzO233+5Ici6++GK/9Q8//HDQ/zncf//9jiTnlFNOCXjMG264IeAHekZGhlOrVi0nIiLC2bJli99r2dnZTmJioiPJ+fHHH0Nu//PPP+9Icm699Va/9fn/8/H5fM6aNWsK7ffVV185kpx27dr5rX/iiSccSU6XLl1CCgo///yz4/P5nGbNmjkZGRkBtxk8eLAjye8Dd+zYsY4kZ9iwYQH3ueCCC4oNsaF+FRdi169f70hy+vXr527z17/+1ZHkLFu2zHGcwCF2+/bt7jkWLVpU7LUKVf7vXsE/kEqrYIhNSUlxIiMjnZiYGCc1NdXdrrKFWMdxnAceeMCJi4sr9DONj493RowY4aSkpBR7jPKWk5PjXruOHTuG9AdGOELsokWL3OtXMEwGUzDEOo7jTJ48udAf48FCbGk/H4rTu3dvR5Lz008/+a0vaYitVauWk5CQEPJ5i5L/+1/UV/4f3mURYidMmBBwv/wbMJ999pnf+vwQe9FFFwXcb9CgQYX+IM/NzXWaNWvmSHK++eabgPvlf14VvPlCiA0PRieoRg4fPqxXXnlFkgrVPY4YMUJ33HGHPvvsM23cuFFJSUnua/njk15xxRUBjzty5Eg9+eSThdbHx8frkksu0Zw5czRnzhxNmjTJfW3hwoXatWuXTj75ZJ100kmF9t2/f78++ugjrV69Wrt371Z2drYkadu2bZIUtEasVatW6ty5c6H1J5xwgiQpLS3Nb/2iRYskSVdddZVq1KgR8JimDz/8UI7j6Nxzz1V8fHzAbfr166cPP/xQX375pf7yl79I8q7h8OHDA+4zcuTIYmtjBw4cqCZNmgR9/d133w2p2719+/Y6/fTTlZycrF9//VX16tXTf/7zHyUlJYVcw1qWevbsKelIl3+DBg30l7/8RfXr1y+z4x933HEaM2aMXnjhBd11111Ba+gqgzvuuEP//Oc/NW/ePCUnJ2vVqlX68ccflZmZqTlz5uidd97Ru+++q8GDB4etjVdffbWWLFmiBg0a6N1331V0dHTY2lIUc4zioxnP9dZbb9WMGTP08ccfa+nSpTrzzDODblvaz4d8GzZs0KJFi7RhwwZlZmYqNzdXklcvvW7duqAPXYbi5JNP1rJlyzRixAjdcMMN6tq1qyIijq6qsKghtor6vCqp8847L+D6E044QSkpKYU+20PZb9GiRX77rV69Wlu3blVSUpK6d+8ecL9+/fpJUpmUd+DoEWKrkYULF2r79u1q3ry5Bg4c6Pda48aNNXjwYH3wwQd6+eWX/YbD2bJli6Qjw+YEEmy9dCQsz5kzR7Nnz/YLsfk1RaNHjy60z/z58zV69OgiaxgzMjICrs+vKyyoTp06klSo1vD333+XdKSGNxT5TzrPnDmz2Afhdu3a5S4fzTXMN3HiRPcDNJBly5aFXDs6ZswYLV++XLNmzVKTJk106NAhjR49usj/2Tdo0EARERHKy8srcsKEkurXr59uv/12PfLIIxo5cqR8Pp/at2+vU089VRdccIHOO++8o/4fbf6wS6+//rpuueWWgE//58u/Bk4RD32Yr5f1gPd169bVyJEj3TFn//jjD82dO1eTJ0/Wtm3bNHLkSP3+++/FjsG5e/du3XLLLYXWH3/88Zo4cWKp2nbDDTdo5syZqlevnhYvXqxjjz22VMcJ5KWXXtLnn39eaP3EiRND/vdp2r17tyS5D5iWVp06dTR58mTdeOONuv322wOO7JGvtJ8Pubm5Gj9+vJ5//vkif++Cfe6F6tlnn9Vf/vIXvfrqq3r11VcVHx+vnj176qyzztKVV14Z9POzKBU1xFZxn+3BHiAuyX75P7+NGzcW++/a/PkhfAix1Uj+h+qhQ4fUt2/fQq/n/0U6e/Zs3XvvvSHdmSzOGWecoaSkJK1fv15ffvml+vTpo507d+rDDz9UzZo1Cz1ElJaWpmHDhungwYO67bbbdMUVV6hNmzaqXbu2IiIi9Mknn2jgwIFBP+iPNuwUJ//uTpcuXQLe8TX16tWrXNtyNC699FJdf/31euWVV9xwGmigflNkZKQ6deqkNWvW6Ouvv9aVV15ZZu158MEHdfXVV2v+/Pn6/PPP9cUXX2jWrFmaNWuWevbsqaVLlyouLq7Ux2/atKluuOEGTZs2TZMmTSpytqn88xT3B8H+/fslSbVr1y51u0JRr149jRkzRl27dlW3bt20e/duffHFF+5EBEW1L7/nxdS3b99Shdibb75ZTz31lOrWratPPvnEHZ2grHz++ecB2ztq1KhShdhVq1ZJOhLaj3aIr2uuuUZPPPGEvv76a7377rvq3bt3wO1K+/nw5JNPasaMGWrSpImmT5+uPn36qHHjxqpZs6Yk6fLLL9e///3vYv+wKs4JJ5ygdevW6ZNPPtGnn36qL7/8UsuXL9enn36qe++9VzNnzgzaW1TeCs7uV1BpP9tLsl9+G5o0aVLoRk9BDRs2LFV7ULYIsdXEtm3b9OGHH0o6MtbnF198EXTbrVu3atGiRe6IAc2bN1dKSkrQ6fSKmmYvf+rDu+66S7NmzVKfPn302muv6fDhw7rssssKTUs4f/58HTx4UBdeeGGhUQEk74nostKqVSutXbtWKSkpOvvss4vdvmXLlpKkU089Vf/6179CPk/z5s21ceNGpaamBiyfqOipCuPi4nTZZZdp5syZ2rx5swYNGqQWLVoUu98FF1ygNWvW6K233tIjjzxSZiMUSEeG0rnuuut03XXXSToyq9zw4cP19ddf6+GHH9bUqVOP6vi33367XnjhBX344YdFDiWUf+dmw4YNRR4v/3exNHevSqNr165q2LChdu/e7d5lLEqbNm2OOvTku+222zR9+nQlJCTok08+UY8ePcrkuKbZs2eX2R29nJwcdzi+c84556iPFx0drfvuu09XXnml7rzzTn3yyScBtyvt50N+W59//nmdf/75hV4vy8+9yMhIDR482C1JycjI0PTp0zV16lT94x//0IUXXnhUfzAGkl9ykpmZGfD1nJwct1QsnPJ/fg0aNGACB0swxFY1MXv2bOXm5qpXr15yjjzQF/Drtttuk+Q/Zmz+XdvXX3894LGLqzEcNWqUIiIi9Pbbb+vPP/8sspQgf8rT1q1bF3rNcRy98cYbIbzb0OVPXfnyyy+79WdFyR9o/oMPPijR/OxHew3Lw9ixY9WgQQM1aNBA48aNC2mf6667TgkJCdq5c6duv/32Yrdfvnx5qdvXs2dPd7alNWvWlPo4+RISEnTHHXdIkvt7HshZZ50lSfr++++DBtmcnBx98MEHftsfreIC5759+9zu5FD+4CgrEydO1COPPKKEhAQtXrzYrWGuzO68805t3bpVUVFRZTb97RVXXKHOnTvrl19+0Ysvvhhwm9J+PhT1uffTTz+Vye9/MHXq1NE999yjunXr6s8//9T69evL/ByJiYmKjo7W3r17A5Yiffzxxzp8+HCZn7ekevbsqYYNG+rnn3/WTz/9FO7mIASE2Goif3D74rqM86ePXLBggVvzc9VVV6l27dpasWKFnnrqKb/tly1bphkzZhR5zBYtWmjAgAHKyMjQHXfcoR9//FGtWrUK+D///Aew3n33Xb+/zHNzc3X33XeXeTH92LFj1aJFC61evVrjxo0r1IWckZGh//73v+73Xbt21cUXX6zNmzfroosuCngH9cCBA3r99dfdhzGkI+GvRo0aevvttzV37ly/7d98882wzGx1yimnuHf1LrroopD2adCggebMmaOIiAg9+eSTGjt2bMD/KaWlpWn8+PEaOnRoscecO3euPvvss0LdiTk5Oe6Dd4H+514a1157rVq1aqWVK1cGHbS8R48e6t+/vxzH0fDhwwvdITp48KCuueYa/f7770pMTNRVV11VJm07+eST9eyzz7qBxrR9+3aNHDlS2dnZat26ddDu7LI2efJkPfTQQ6pbt64VAfbXX3/ViBEj9Mgjj0g6MjlHWf3u+Hw+TZs2TZKCjkta2s+H/M+9Z555xu/fwbZt2zRixIgyCXh//vmnpk+fHrCWc/ny5dq3b59q1KhRLn8gRUVFuQ+NTp482e89fvfddxo/fnyZn7M0oqKiNGXKFDmOowsvvDBgjXZubq4+/fRTffXVV2FoIQqinKAaSE5O1oYNGxQTE1PkQPaSdNJJJ6lbt25atWqV5syZo5tvvlnNmjXTiy++qOHDh+uGG27QSy+9pA4dOigtLU3Lly/XjTfeWOysOKNHj9bHH3/sjmKQf3e2oPPOO0/du3fXt99+q2OPPVZ9+/ZVXFycVq5cqa1bt+r2228PWGZQWrVr19YHH3ygwYMHa9asWZo7d65OPfVU1a5dW5s3b9bq1at18skn+5UazJo1S/v27dNHH32k4447Tp07d1bbtm3lOI5SU1P13XffKTs7W2vXrlXjxo0lHamRmzZtmm677TZddNFF6tWrl5KSkvTLL7/o66+/1k033WTFrFKSdP7552vBggUaMWKEZs6cqVdeeUU9evRQ69atdfjwYW3cuFHfffedHMfRKaecUuzxkpOT9eSTT6phw4bq2rWrGjVqpMzMTH311VfauXOnmjdvXuSd05KIiYnRvffeq1GjRhU5Fe1rr72mc845RytXrlS7du3Up08fNW3aVOnp6VqxYoX27Nmj+vXr67333itUElNav/zyi6699lpdf/316tixo5KSkhQZGam0tDStXLlSOTk5ql+/vt58880Kmcb1gw8+cB/wPOaYY/TMM88E3K5hw4Z69NFH/dZt27ZNF154oft9/oONH3zwgd/vxLPPPqtu3bqVuG0pKSnutLd5eXlKT09XSkqKfvnlFzmOo8TERP3rX//SZZddVuJjF+Xcc89Vv3793NFGAinN58Mdd9yhRYsW6cUXX9TSpUvVrVs3ZWRkKDk5We3atdOFF15Y6I/fksrOztbNN9+sW2+9VR07dlT79u0VFRWl1NRUN5DdeeedSkxMPKrzBHP//ffrs88+04svvqjk5GR16tRJaWlp+uabb3T55Zdr2bJl7oO24TR+/Hht2rRJjzzyiE4//XSddNJJOuaYYxQbG6vt27drzZo12rdvn5577rmQPt9QzipwOC+EyZVXXulIci655JKQts8fO/WEE07wW798+XJn4MCBTp06dZxatWo5Xbt2dZ5//nnHcbyxTIM5dOiQU79+fXcs119//TXotpmZmc4dd9zhHHfccU7NmjWdRo0aOUOHDnW++eYbdxzHguNEBltvKqqNu3btciZPnux07NjRiYuLc2JjY5127do5w4YNCzgmam5urvPGG284gwcPdho3buxERUU5DRo0cDp06OCMHj3amTt3bqGB+x3HcebNm+ecdtppTlxcnFO7dm2nT58+zrvvvlvkOIoqMOZiMKFOdhCKQOPEFpSZmek8/vjjzoABA5wmTZo40dHRTq1atZxjjz3WGT58uLNgwYJCg7MHsnr1amfixInOaaed5jRv3tyJjo52EhMTne7duzv/93//5+zevTvkdjtO4XFiC8rNzXU6duwYdFzdfH/++afz9NNPO/369XMaNGjgREZGOnXq1HG6du3qTJo0yW9CjaKEOrbpDz/84Dz++OPOeeed5xx//PFO3bp1ncjISKd+/fpOnz59nKlTpzq7du0K6ZxlIZRxQIP9zub/Phf3VdzvdEH519L8io6Odho2bOh06tTJGTFihPP66687Bw8eLNV7DjRObEErV670O38gpfl8+P77753zzz/fadq0qVOzZk2nffv2zm233eZkZGS4/x4L/q6WZJzYnJwcZ8aMGc7f/vY35/jjj3cSEhKc2NhYJykpybn44oudJUuWFHsMU0kmO8i3YsUK55xzznHq1KnjxMbGOp07d3aeffZZJy8vr9hxYoONvxrs2uSPExvsd6y4a/fFF184V1xxhdO6dWsnJibGiY+Pd4499lhn6NChzksvveTs3bu3RO1E+fA5ThlV/gMAAAAVhJpYAAAAWIcQCwAAAOsQYgEAAGAdQiwAAACsQ4gFAACAdQixAAAAsA4hFgAAANYhxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwDiEWAAAA1iHEAgAAwDqEWAAAAFiHEAsAAADrEGIBAABgHUIsAAAArEOIBQAAgHUIsQAAALAOIRYAAADWIcQCAADAOoRYAAAAWIcQCwAAAOsQYgEAAGAdQiwAAACsQ4gFAACAdQixAAAAsA4hFgAAANYhxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwDiEWAAAA1iHEAgAAwDqEWAAAAFiHEAsAAADrEGIBAABgHUIsAAAArEOIBQAAgHUIsQAAALAOIRYAAADWIcQCAADAOoRYAAAAWIcQCwAAAOsQYgEAAGAdQiwAAACsQ4gFAACAdQixAAAAsA4hFgAAANYhxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwTmS4G1CR8vLytHXrVsXHx8vn84W7OQBw1BzHUWZmppo1a6aICO5LAKg+qlWI3bp1q1q2bBnuZgBAmdu8ebNatGgR7mYAQIWpViE2Pj5e0pEP+zp16oS5NQBw9DIyMtSyZUv38w0AqotqFWLzSwjq1KlDiAVQpVAiBaC6oYAKAAAA1iHEAgAAwDqEWAAAAFiHEAsAAADrEGIBAABgHUIsAAAArEOIBQAAgHUIsQAAALAOIRYAAADWIcQCAADAOoRYAAAAWCcy3A0AEF5tJi50l1MfHBLGlgAAEDruxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwDiEWAAAA1iHEAgAAwDqEWAAAAFiHEAsAAADrEGIBAABgHUIsAAAArEOIBQAAgHUIsQAAALAOIRYAAADWIcQCAADAOoRYAAAAWIcQCwAAAOsQYgEAAGAdQiwAAACsQ4gFAACAdQixAAAAsA4hFgAAANYhxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAVQrDYTF6rNxIXhbgYAAC5CLAAAAKxDiAUAAIB1IsPdAAAVj9IAAIDtuBMLAAAA6xBiAQAAYB3KCYAqyCwXSH1wSLnvBwBAReNOLAAAAKzDnVgApcJdWwBAOHEnFgAAANYhxAIAAMA6lBMACIixZAEAlRl3YgEAAGAdQiwAAACsQ4gFAACAdQixAAAAsA4hFgAAANYhxAJVXJuJCxlpAABQ5RBiAQAAYB3GiQWqCe7GAgCqEu7EAgAAwDqEWAAAAFiHEAsAAADrEGIBAABgHUIsAAAArEOIBQAAgHUYYgtAyBimCwBQWXAnFgAAANYhxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwDiEWAAAA1iHEAgAAwDqEWAAAAFiHEAsAAADrEGIBAABgHUIsAAAArEOIBQAAgHUIsQAAALBOZLgbAKBk2kxc6C6nPjgkjC0BACB8uBMLAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwDiEWAAAA1iHEAgAAwDqEWAAAAFiHEAsAAADrMGMXUEWYM3kBAFDVcScWAAAA1uFOLGAx7r4CAKor7sQCAADAOoRYAAAAWIcQCwAAAOsQYgEAAGAdQiwAAACsQ4gFAACAdQixAAAAsA7jxAKVmDkObOqDQ8LYEgAAKhfuxAIAAMA6hFgAAABYhxALAAAA61ATC+CoUbsLAKho3IkFAACAdQixAAAAsA4hFgAAANYhxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKzDZAeAJcwJBQAAqO64EwsAAADrEGIBAABgHUIsAAAArEOIBQAAgHUIsQAAALAOIRYAAADWIcQCAADAOowTC6BMmePZpj44JIwtAQBUZYRYoBIg+AEAUDKUEwAAAMA6hFgAAABYhxALAAAA6xBiAZSbNhMX+tX7AgBQVgixAAAAsA4hFgAAANYhxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwDiEWAAAA1okMdwMA+GOaVgAAisedWAAAAFiHEAsAAADrEGIBAABgHUIsAAAArEOIBQAAgHUIsQAAALAOIRYAAADWIcQCAADAOoRYAAAAWIcQCwAAAOsQYgEAAGAdQiwAAACsQ4gFAACAdQixAAAAsA4hFgAAANYhxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwTmS4GwBUV20mLgx3EyqM+V5THxwSxpYAAKoKQiyACkWgBQCUBcoJAAAAYB1CLAAAAKxDiAUAAIB1CLEAAACwDiEWAAAA1iHEAgAAwDqEWAAAAFiHEAsAAADrEGIBAABgHUIsAAAArEOIBQAAgHUIsQAAALAOIRYAAADWIcQCAADAOoRYAAAAWCcy3A0AqpM2ExeGuwkAAFQJ3IkFAACAdQixAAAAsA4hFgAAANYhxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1mLELqADM1AUAQNniTiwAAACsQ4gFAACAdQixAAAAsA4hFgAAANbhwS4AYWM+8Jb64JAwtgQAYBvuxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwDiEWAAAA1iHEAgAAwDqEWAAAAFiHaWcBVApMQQsAKAnuxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwDiEWAAAA1iHEAgAAwDqEWAAAAFiHEAsAAADrEGIBAABgHUIsAAAArEOIBQAAgHUIsQAAALAOIRYAAADWIcQCAADAOoRYAAAAWIcQCwAAAOsQYgEAAGAdQiwAAACsQ4gFAACAdSLD3QAAKKjNxIXucuqDQ8LYEgBAZcWdWAAAAFiHEAsAAADrUE4AlBOzSxwAAJQt7sQCAADAOoRYAAAAWIcQCwAAAOtQEwscJYaDAgCg4nEnFgAAANYhxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwDiEWAAAA1mHGLgBWYGY0AICJO7EAAACwDndigTJk3i0EAADlhzuxAAAAsA4hFgAAANYhxAIAAMA61MQCqNSoMwYABMKdWAAAAFiHEAsAAADrEGIBAABgHUIsAAAArEOIBQAAgHUIsQAAALAOIRYAAADWIcQCAADAOoRYAAAAWIcZu4BSYBYpAADCizuxAAAAsA53YoFicNcVAIDKhzuxAAAAsA4hFgAAANahnABAlWGWfqQ+OCSMLQEAlDfuxAIAAMA63IkFYB3uuAIAuBMLAAAA6xBiAQAAYB3KCYD/jy5qOzGOLwBUT9yJBQAAgHUIsQAAALAOIRYAAADWIcQCAADAOoRYAAAAWIcQCwAAAOswxBaqHYbSAgDAfoRYIADGHgUAoHKjnAAAAADWIcQCAADAOoRYAAAAWIcQCwAAAOvwYBeqDR7WAgCg6iDEoloj2FZd+T9bhlEDgKqJcgIAAABYhxALAAAA6xBiAQAAYB1qYgFUaUwzDABVE3diAQAAYB1CLAAAAKxDOQGqNIbQQjCUGQCA3QixAKoN/qgBgKqDcgIAAABYhzuxqHK42wYAQNXHnVgAAABYhxALAAAA61hTTjBt2jS9//77SklJUWxsrPr06aOHHnpIxx13XLibhkqAEgIcjWC/P4xaAACVlzUhNjk5Wddee6169uypw4cP64477tA555yjn3/+WXFxceFuHsoZIRUAAJisCbGLFi3y+3727Nlq1KiRvv32W51xxhlhahXKE8EVlQV3agGg8rEmxBaUnp4uSapfv36YW4KSyA8D5v/8GXQelVVxf0jxuwsA4WNliM3Ly9ONN96oU089VR06dAi6XVZWlrKystzv//jjD0nSli1bVKdOnXJvJwo7nLFb0pGfQcF1Ra0HKjvzd7ci7du3T5K0d+/esJwfAMqa4zjKzMxUs2bNFBERfAwCn+M4TgW2q0xcc801+uijj/T555+rRYsWQbe75557NHXq1ApsGQAAAMrC5s2bi8x51oXY8ePHa968efrss8/Utm3bIrcNdCe2TZs2Ok2DFamo8m4qAJS7gzqglfqvNm/eTA8TgCohIyNDLVu21L59+5SQkBB0O2vKCRzH0XXXXae5c+dq2bJlxQZYSYqJiVFMTEyh9ZGKUqSPEAvAflFOtCSpTp06hFgAVYrP5yvydWtC7LXXXqs33nhD8+bNU3x8vLZv3y5JSkhIUGxsbJhbBwAAgIpkzYxdzz33nNLT09WvXz81bdrU/XrrrbfC3TQAAABUMGvuxFpWugsAAIByZM2dWAAAACAfIRYAAADWIcQCAADAOoRYAAAAWIcQCwAAAOsQYgEAAGAdQiwAAACsQ4gFAACAdQixAAAAsA4hFgAAANYhxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwDiEWAAAA1iHEAgAAwDqEWAAAAFiHEAsAAADrEGIBAABgHUIsAAAArEOIBQAAgHUIsQAAALAOIRYAAADWIcQCAADAOoRYAAAAWIcQCwAAAOsQYgEAAGAdQiwAAACsQ4gFAACAdQixAAAAsA4hFgAAANYhxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwjlUh9rPPPtN5552nZs2ayefz6T//+U+4mwQAAIAwsCrEHjhwQJ07d9YzzzwT7qYAAAAgjCLD3YCSOPfcc3XuueeGuxmoQJGNG7nLh3fsDGmfGh2Oc5d9u/4o8f5lxWy7KVg7QtnefG9S8PdnHstJrBdw+2BK2r6C58hpEOcuRySvKvZYwdoX7P0U1cZg56jonz0AoPxZFWJLKisrS1lZWe73GRkZYWwNAAAAyopV5QQlNW3aNCUkJLhfLVu2DHeTAAAAUAaq9J3YSZMmacKECe73GRkZBFnLlKYbOPfHde6y2aUcrDs8+/gW7nJ0ypZij29uX1DUngPusmOsN7vJzZKA3SfXd5cb/m9vwO3z+nbzlgucL6t9XXc5Zq/XruxgDTS6+s33al5n83zm+zlsXFfJ/32Y1zx9TG93uV5ysIZ4zH39yiWMNhW85tHGclFlB+WBMgUAqByqdIiNiYlRTExMuJsBAACAMlalywkAAABQNVl1J3b//v3asGGD+/1vv/2mNWvWqH79+mrVqlUYW4aKVlS3sdnFG8qT+X5P0BvHDVY2YHaxS/7d4bnGerNbPlg3flZCg4DnyOzT1l2O2ZsTcJuCr5nvY9PUPu5yToJXhFD/e5+73HCPd23+GJLkrTfKGkwFR0YwRyGQ8V4TF2701hvX03xPtX/Z5y77fQgFKbtQgVEODhvLJR0F4miFMnoDZQYAUP6sCrHffPONzjzzTPf7/HrXkSNHavbs2WFqFQAAACqaVSG2X79+chyn+A0BAABQpVkVYmGPknatBhv8PrfAE/ElPneQkQqCtclvvVFOYJYQFGyT2e293xgtIP7L39xl8z2ZZQrNF+8JuG+wEoLfLoz2+94sD1CSNypA23e845ojICRs9MZNNs9X8w+v5MAsE0hP8h6MzEowziUpIdUrnjhUzyuvT1Dw0RtKIuioBfIvDQnl9yXYBBimoy0BoIQAACoWIRYAqoAOUz5WREytcDdDqQ8OCXcTAFQTjE4AAAAA63AnFuWipF2rftuXols2lPOZ3c4KYYB8xyghMLvYowtsv8vorjefzN9lPPFvduObXfQyls1tooKcOyrd/+/OrARv2ezeXz/We6/HvuSNNpA2oEHA7bec7y03XO61ySwzqFmgF94seYgJPKCBH3NEAr/raXTvmz9HswTAbySEAvsEmxjCbzKHIKUJoZSrFPz9oGygaG0mLgx3E8oFd5iByoc7sQAAALAOIRYAAADWoZwAlUKwkQOCrS+qSzfoRAhBup2dICMP1DDKD8zufb+yBPmXAQRjjipQ/3tvfb2XV3jnC9J9bp47el9Nv+Nmdsp2l2PSo4xXvKHozGOZoyGYkmZ722w5y5gQ4X/73GVzlANJSm/jtSXYSAXBJk4Idj3Nax7sGhRk/i5EB93KE2x0glCOj+qrosokKFsAQsedWAAAAFiHEAsAAADrUE6ASiGkyQcMRT0xHmyfgxf2cpfNiQiCPRFvdjtn9mkbcF9JUoEn5/OZIxXs7dQu4Dbbb+rjncMoDWjxQQ13eeOoeHf5hIkb/PbPXuVNLHDWvz51lz8df6q77F/KEHg0hJ3dvNKApn3S3OXfO3jnjvzcf7KDmHSvZCG9TY2A681raE7yEIw5okBUEduZx4oyShDM+fxCKRsoaRkLUN4oWwBCx51YAAAAWIcQCwAAAOtQTgAAQDVTVmULlCUgnAixqFCh1BqWpk4x2CxM5lBasXNXejsYxzJrZc1tzGOas00VHGLLFKzuM+nGr9zlP8b0dpfN+tED6d4/x11dvE6S2Dre8FKN5h3yO+6qt7xa1tfW93SX44yZwJole7Numcc96e9efW3WQye4y9u+bO6tb+XV6XrVsYVl1zW/82pnf/3nMe5y23e84b38Zuwy6mDNa1tUPW1E8ipvO+Nnae4frNb5cAizdAWb3Q0AUHlQTgAAAADrcCcWAACEhPIBVCaEWFSoUIYrKulwW5J/CYHJnOnJF2QmL7OEwCxZyDa6vM3u64LDe2W1r+sum2UH+431MUZZw4Ehme5y4zrecs1nmxhH9WbASt9Xx13+Kf0kv3PHGINKHdhc22tjgtelv7eTtxyV7u275Xavq3/XWd4QWfG/e8eM3ucNz+VfMiAdaO1t1+qkre7ypp+ausv1v/fObc74tfv0HHe54fIkd9mcwSzozGvFvBZwmxIOtxXK9gCA8CLEAgCAkFTUOLblhTvJVQs1sQAAALAOd2JRpQUrMzC7joM+uW48ER8Z5An4gvYHKS1YP9bbp9lr3lxUe9p4pQJq4y02X+w9yX+ontcNn9naf9asaO8UykvwuugbrfJGJPittVcSkLjGK1PYer038kD0594oB2ZXf2ydQ8Y2RlslZbU67C6bIxo0M87hP4eWp+Fy7xrU/MNrqzmDWYs3vNETzBEMpOCjGASbpSvYaAPMzAUA9iLEAgCAaqE8yiEoUQgfygkAAABgHe7EosoJ9uS62V3sN4lCCMc0t88rMPB+zF6v+90cDcGUl+B1vae3CfzEf0SHDHc5pZM3tUDMJq+EICfB63qXpOh93t+h8d97x23x0I/ucm6m15U+oH+KuzxzyZnu8vAxS9zluY+fZZzB6/bffXq2TB+f9aS7/I/1l7vLz496w12+5Ilb3WXzvZplEAkbveuXYBx/1xBv1IKG/9vrd+5gZQPBBNs+lBKCgr9PJd2HMgUAKB+EWAAAgBKijCD8KCcAAACAdbgTi2opz5h8YFdSjLtc72Wv69d8Ij7b2L6oJ+VNZhd20uw44xXvif8WD3lP4P9udPvvMUYCyOzkdeObT/VL0kl/98oGln91YsDlVidtc5dfW9/TXW6W7I0cMFNeaYE6eetjWu73ln/0H51g6Df/cJfjFnrlD0OHeOsbpBqjIfQ1R1bw/n5ON66/WTZQ29g6p4F5/aRo49oWNVpEvmDlBMFKT8xjHg4ywkVRgpUQUGYAAGWHEAsAAFAA5QKVH+UEAAAAsA53YmGNorpiS9pNa5YE1Ev2tjcnPsg2urAjkle5y46xjeQ/oYK5v9klbU4s0OwpbxSBLbcf4x1osjfBgVlCUJTktce6yzHpgf8m3ZHhdfUfzPAmNUi6fa3XDuM4vw6c6S5P2+Otf/Pz/n7H7dZ8k7v8k05yl83SgkNGT3+zZG9khdq/eO9198neZA4ms4Qgq75/GYUKjBCRz/w5BRPS70qIXf0l/b2jhAAAyg4hFgAAoIDymBihKqkM5RaUEwAAAMA63ImFNYrqivWbvCBIF6/Z1X84SAmAWRoQbRzHLCEo2P2deWkfd7nds95oA5l92rrLbf75m7u89sFW7vIJE70u+V+/9EoLIoxJDcyJEg4MyfQ7d1+jSz9ZxyqQiM3ec/59T/nZXTZHQ4it442Y8J8DXje+OfFBjLxRCyT/UoYIY0SDtnO9UoiNowJPxnConncNEzZmucvmtU1cuNFdzjKupeRfNhAZ5OdkjkhQXt34lAcAsEFluGtaHrgTCwAAAOsQYgEAAGAdyglgjaMdKN582t0sFZDR7Xzwwl7ucvyXXgmA2TWduNB/4PxEYznbeGr+UD3vb8Rd//RKBRou97reMwt0k+er/703MUBma68bPu53/6f0l3fyJjWINSYmMEchMCc7+OrjDu5y0z5p7nLjOv5lCvlO+vtPAfeV/EsQIr/3JkLYcpZ37hiv2sFvxIXsTd57ymztbe9nSJK7aE6CIEk5xuQTMkaaMMsREjZ6P+8oY6QIsygiz/h5hTKyAQBURlW1XKA43IkFAACAdQixAAAAsA7lBKjUSlNCEGykArPbOfix2gbcxu84BZglBFF7DrjLWd28bvLENbnucnqbGu5yZmtvue073gQAKTd7EwbEGF3vBZ1ujDZwQQOvO3zKsyPc5VSz4KGV16WfusVbH5HufRSsbdLcXTZHMCjILFmIN9bH/+512O81Ri0wzxG9z9u+0SqvLOG3C7332vB/XtlGwREhzPIC8+eUsDFweYB39f1HozAnvTiswAr+7BmRAEBlU9SYtlW51IA7sQAAALAOIRYAAADWsa6c4JlnntEjjzyi7du3q3Pnznr66ad18sknh7tZKCdl2XUb7Fhmd3HtX/a5yz5jvTmKQMzeHL/9s+p7IwakJ3nd3i3e8CY+MEsOGhn7b73e6973Pet1n8fW8Y5pjkiQ2dobtUCSNj50grs8pY03ekB2XQXUpsUud3n/v5u5y2a3//oDjd1ls+Sg78Af/Y61/CtvZATzfJmdvPdnnm/TT03d5Zh073zm9TNHZTBHkzAnPpD8fx6xRrPM8gCZkyAYoxOYI02Y64N9GJrbSJIoJwBgkXBOn1vepQxW3Yl96623NGHCBE2ZMkWrVq1S586dNXDgQO3cyf9UAAAAqhOrQuz06dM1btw4jR49WieeeKJmzJihWrVq6eWXXw530wAAAFCBrCknyM7O1rfffqtJkya56yIiInT22WdrxYoVAffJyspSVpY3L3tGRka5txOVS0nLEXJ/XBdwfezc4MeJN7qtzaf0zRICs8vc1HqKN5rB2gdbucstXvNGLThk9GbnJOT57Z90+1p3edVbXjlBljEKgdmlX+P+Bu7y7lH+ZRH5zEkNYjt4/2bM40tSq8HeZAlmacIpA9cH3McrDpCyErxlc4QGcxQHc6QH81oWFGwEC3O9+XM115ulBUHLTYKeGQDsVFVGLLDmTuzu3buVm5urxo0b+61v3Lixtm/fHnCfadOmKSEhwf1q2bJlRTQVAAAA5axK32SYNGmSJkyY4H6fkZFBkAUAAFVaVbnTWhxrQmzDhg1Vo0YN7dixw2/9jh071KRJk4D7xMTEKCYmpiKaBwuYg9yb3cuhjFpgbmMeR5IU5Gl380n5KGO9+dS9qeFyr+RgVxfvKf1TjFEBfnrhJL99khOOdZdbpJpD+nuTBpiTHTRM8s4RYYwQEJXudcqYJQttn/KOs/V6/3Icv9EGhmS6y8fGef9G/3eaVyIR+Xkdb3vj3Adae8uH6nnt2HWpV/pgTo4g+U+QYF7zGsby4SAlBMEE+/04WqWZsAMAUDxrygmio6PVvXt3LVmyxF2Xl5enJUuWqHfv3mFsGQAAACqaNXdiJWnChAkaOXKkevTooZNPPllPPPGEDhw4oNGjR4e7aQAAAH6qS7d+uFgVYocNG6Zdu3bp7rvv1vbt29WlSxctWrSo0MNeAAAAqNqsCrGSNH78eI0fPz7czYCFgtU55vXt5i5HJK9yl4MN13S4iHpJ8x+UuX+eMUSUOXSUOcRTzT/qussJG73hr1bt84aparTRqwU9sp23vHGU+YpXHzv9tDfd5XkneO/VHErLrINtluzVqOZO3uMuD2+S4nfuNz/v77Uj2avzfW14T3c5bqE56Jh33N2ne+/PrAU2Zw6L+92rC07wq/cNzryeoQy9FWxfU1F1rKHU0VIHCwDlw7oQCwAAYIPSTPlKCULorHmwCwAAAMjHnVhUaaEMb2SWEJR0X6nAkFtGl7S5voZRQmAOsZV+cn13ueH/9rrLaQO84aUOtPa6+rNS/Wf+itnrdcu3+MCb+Sq9jbc8adMIb39jJi9z8LlHh7zmLk9I+Ku73MbY5rX1PWWKM4bJMmcOO9cYYus1eftkba7tbT/be0/pSd4xzVIGsyRiV5caMiWu8ZYjkosfSsv8WQQrB/H7OYZYAlCWQ3EBAEqGEAsAAHCUKAOoeJQTAAAAwDrciUWVdjRPhpdmX799gjwRbxYENDTKDNaP9Wabivs98PHNGa0kaWvf6IDbBdv/+Me8mbV+M2bEmvSGV3LQpk+au7z/383c5YPGiAKSZM479r80b2auHVO8mbwiB3izdLVe7I10YJ47cY1XNmC+v4SNWe5yi0/9z20KVipgXnO/WdKCjEYRrDSgqNm+QvkdYcYuACgfhFgAAICjFMpIBJQclC3KCQAAAGAd7sSiyimP7luzy7rgoPjBurCdRK88wHzm3tw/25gEof73xuD+xqQGUYu9koPdxmgGBfep+Yf3xP+uLt420fu85d+nesUM8Qu9VmW29o6TuiXRXY4wJh+I2eRfupCVYCwbIw/8dqn3t3FTszQhtYm73O7ZDe6yeQ0O1TPHTPCYE0RI/tcwlIkMolO2BFx/2FgO9fempL9TlBAAFY87ntUDd2IBAABgHUIsAAAArEM5AaqE8n4CvDSD2hcsO8hndp+b3dwJ8tZn1fe6/bec5T3J3/Yd7wl/Sdrfvq67bD7Zb3bXZ/Zp6y5vbe11+2uIN1KBfvRGEeh7wnp3+fdMryTCLDOQpPjfvTZmp3vnHnHREnd5zvv93eVEc/KCId4MB+YkD1ndarrLO43lFm/4lwOY7ylmr3fdHKPswLz+ZmlHSX+WRf0+MfIAUD4oB0AouBMLAAAA6xBiAQAAYB2f4zhO8ZtVDRkZGUpISFA/XaBIX1TxO6DaOdruYXMUg5IOnu/fRe4N7m8+mZ82oIGCabTKG9EgPcl7yt8ctaDgZAnuuY3RCeJ/9z4S9pqjE7Tc77dP3ML4gPsHm7zAbIf5/szSiWDra/+yL2C7Jf/rXNTEBG5bjesc/+Vv7rKt5QCHnD/1uT5UyxvfVkRMrXA3BwgryhCqhvy8lp6erjp16gTdjjuxAAAAsA4hFgAAANZhdALAEEqXcsEua3OfYCUEZplBsMkRgnWlS3HuUkJqrkzmPsFGOtg4yvtbNf77Gu5ydl3vOGYJwNa+XmlA27nZxtn8JzvYOMo7d4sPvOOmt/GWzfaaJQHmyAHpxkgF5mQHZvlBTgPvGkj+JRZ/jOntnW9jltda43qYP6Pav9QLuD5UjEiA6ohuelRG3IkFAACAdQixAAAAsA7lBEAJFexCDjYigV+3s7He3D7b6CY3u8gjkgMfJ6LAuYN1bZtd9Mc/5k0mEGxyBLOEIMqYuMAsazDLBCSp4XJv5IL4Lze6y+bkA+Z7Ms+9dazXpV//e++Y5mgIbecGHqGhIHOyhGATHMi4NsEmoQi1TIASAlQ1lArAVtyJBQAAgHW4EwsAgAW4Ywr4I8Si2iivp8qDlhAEOYfZnR1h7Lv/wl7ucvwu7zhmt3je8V5XvSQFGzug3ssr3OW0m/oEbIc5OkHbud5T/eZECaaCIyOYtlx+TMDjtnvWGyHg0Mn1iz3fcdN+9Y5jvNeCJQDmhAXmqAfmtTL3CWUSBMoEAMAulBMAAADAOtyJBQCgkqF0ACgeIRbVRkV0F4dyDrPL2/wHGP/lbwG3MbvFs4wn/CX/yQ6ClR00WnXIXTa77psv9p7qLziZQKDjF/TbhWYBgzcxgTm6wbpJ7dzlZsl5CiRxoTeygd+IAkUINnGCyfxZMEEBAFQ9lBMAAADAOtyJBQAgDCgZAI5OqULsli1b9MEHH2jTpk3Kzs72e2369Oll0jCgqjK7v82u7by+3dzliORV7rLZFW52oxc8lvk0f3SKNyqAXxe9UU6w2xgtwJwwoGGQCQoKlhaYIwmYoxOYoxjEzl3pLpuTPJjlC2b7zJEeoosqAQhSKuB3DYzNK7q0gPIFACh/JQ6xS5Ys0fnnn6927dopJSVFHTp0UGpqqhzHUbdu3Yo/AAAAAHCUShxiJ02apFtuuUVTp05VfHy83nvvPTVq1EhXXHGFBg0aVB5tBACg1Oi2B6qmEofYtWvX6t///veRnSMjdfDgQdWuXVv33nuvLrjgAl1zzTVl3kjAdqF0L5slBCZzYH+ze17y76I3SwjMbnXzuAkNvN6Snd1qusuhjE5gHl8KXkLgN8qC0b60AQ0Cbi8Z5zZKKhyjrKHgB1VI1zCECQ7KCyUEAFD+Sjw6QVxcnFsH27RpU23c6A2Ps3v37rJrGQAAABBEie/EnnLKKfr88891wgknaPDgwbr55pv1ww8/6P3339cpp5xSHm0EAKDU2kxcWGgdJQaA/UocYqdPn679+/dLkqZOnar9+/frrbfeUvv27RmZAAgiWPdyZAhd3n4D+xfY/rDxNL95LL+uf3MfY33zPd6oAMFGKjBHPyg4EUFMuuMumyUEZvmD2fbmi/cEPJ8UFfA4wUYUCPR9oPVF7V+coso/GHkAACqHEofYdu28GXji4uI0Y8aMMm0QAAAlwV1VoHoqcU1su3bttGfPnkLr9+3b5xdwAQAAgPJS4hCbmpqq3NzcQuuzsrKUlpZWJo0CAAAAihJyOcEHH3zgLn/88cdKSEhwv8/NzdWSJUvUpk2bMm0cUJ0Eq+E061ILMofYClYfG2xGLPMff+LCIDN/GedeP9a/JvbYl/YG3MesazXXm2r+kecum8OGZV7Yy9gq+NBiwWpcg9WolrR2tajtqYMNH8oGAJhCDrFDhw6VJPl8Po0cOdLvtaioKLVp00aPPfZYmTbO9MADD2jhwoVas2aNoqOjtW/fvnI7FwAAACq3kENsXt6ROydt27bV119/rYYNG5ZbowLJzs7WpZdeqt69e2vmzJkVem4AAABULiUeneC337yuwkOHDqlmzZpFbF12pk6dKkmaPXt2hZwPqAhH2/1do8CwV4H2N/+Rm93wZld/Vn1vmCtzli5zm2bJ3pBakv8sX+ZMWdnGrFvmUF9mWUOU8WyoOWxYwbKBQO2WgpdeMPxV5UY5AICyVOIHu/Ly8nTfffepefPmql27tn799VdJ0l133cUdUgAAAFSIEofY+++/X7Nnz9bDDz+s6Ohod32HDh300ksvlWnjjlZWVpYyMjL8vgAAAGC/EpcTzJkzRy+88IL69++vq6++2l3fuXNnpaSklOhYEydO1EMPPVTkNmvXrtXxxx9f0mZKkqZNm+aWIQBVkTnaQDDButUjjPWxIZwrlG0k/9KCw+YLR9G9z2gBxftx6kDVqVMn3M0AgApT4hCblpamY445ptD6vLw85eTkBNgjuJtvvlmjRo0qcpujmUBh0qRJmjBhgvt9RkaGWrZsWerjAQAAoHIocYg98cQTtXz5crVu3dpv/bvvvquuXbuW6FiJiYlKTEwsaRNCFhMTo5iYmHI7PgAAAMKjxCH27rvv1siRI5WWlqa8vDy9//77WrdunebMmaMFCxaURxslSZs2bdLevXu1adMm5ebmas2aNZKkY445RrVr1y638wIAAKDyKfGDXRdccIHmz5+v//73v4qLi9Pdd9+ttWvXav78+RowYEB5tFHSkfDctWtXTZkyRfv371fXrl3VtWtXffPNN+V2TgAAAFROPsdxnOI3qxoyMjKUkJCgfrpAkb6o4ncAgErukPOnPteHSk9P58EuAFVCfl4r7nOtxOUE+bKzs7Vz5053Jq98rVq1Ku0hAQAAgJCUOMT+8ssvGjNmjL788ku/9Y7jyOfzKTc3t8waBwAAAARS4hA7atQoRUZGasGCBWratKl8Pl95tAsAAAAIqsQhds2aNfr2229LPQEBAAAAcLRKPDrBiSeeqN27d5dHWwAAAICQhBRiMzIy3K+HHnpIt912m5YtW6Y9e/b4vZaRkVHe7QUAAABCKyeoW7euX+2r4zjq37+/3zY82AUAAICKElKIXbp0qbucmpqqli1bqkaNGn7b5OXladOmTWXbOgAAACCAEk92UKNGDW3btk2NGjXyW79nzx41atSoUt+JZbIDAFUNkx0AqGpCneygxA925ZcNFLR//37VrFmzpIcDAAAASizkIbYmTJggSfL5fLrrrrtUq1Yt97Xc3FytXLlSXbp0KfMGAgAAAAWFHGJXr14t6cid2B9++EHR0dHua9HR0ercubNuueWWsm8hAAAAUEDIITb/4a7Ro0frySefpPYKAAAAYVPiGbtmzZpVHu0AAAAAQlbiB7sAAACAcCPEAgAAwDqEWAAAAFiHEAsAAADrEGIBAABgHUIsAAAArEOIBQAAgHUIsQAAALAOIRYAAADWIcQCAADAOoRYAAAAWIcQCwAAAOsQYgEAAGAdQiwAAACsQ4gFAACAdQixAAAAsA4hFgAAANYhxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwDiEWAAAA1iHEAgAAwDqEWAAAAFjHihCbmpqqq666Sm3btlVsbKySkpI0ZcoUZWdnh7tpAAAACIPIcDcgFCkpKcrLy9Pzzz+vY445Rj/++KPGjRunAwcO6NFHHw138wAAAFDBrAixgwYN0qBBg9zv27Vrp3Xr1um5554jxAIAAFRDVoTYQNLT01W/fv0it8nKylJWVpb7fUZGRnk3CwAAABXAiprYgjZs2KCnn35a//jHP4rcbtq0aUpISHC/WrZsWUEtBAAAQHkKa4idOHGifD5fkV8pKSl++6SlpWnQoEG69NJLNW7cuCKPP2nSJKWnp7tfmzdvLs+3AwAAgAoS1nKCm2++WaNGjSpym3bt2rnLW7du1Zlnnqk+ffrohRdeKPb4MTExiomJOdpmAgAAoJIJa4hNTExUYmJiSNumpaXpzDPPVPfu3TVr1ixFRFhZCQEAAIAyYMWDXWlpaerXr59at26tRx99VLt27XJfa9KkSRhbBgAAgHCwIsQuXrxYGzZs0IYNG9SiRQu/1xzHCVOrAAAAEC5W9MmPGjVKjuME/AIAAED1Y0WIBQAAAEyEWAAAAFiHEAsAAADrEGIBAABgHUIsAAAArEOIBQAAgHUIsQAAALAOIRYAAADWIcQCAADAOoRYAAAAWIcQCwAAAOsQYgEAAGAdQiwAAACsQ4gFAACAdQixAAAAsA4hFgAAANYhxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwDiEWAAAA1iHEAgAAwDqEWAAAAFiHEAsAAADrEGIBAABgHUIsAAAArEOIBQAAgHUIsQAAALAOIRYAAADWIcQCAADAOoRYAAAAWIcQCwAAAOsQYgEAAGAdQiwAAACsQ4gFAACAdQixAAAAsA4hFgAAANaxJsSef/75atWqlWrWrKmmTZvqyiuv1NatW8PdLAAAAISBNSH2zDPP1Ntvv61169bpvffe08aNG3XJJZeEu1kAAAAIg8hwNyBUN910k7vcunVrTZw4UUOHDlVOTo6ioqLC2DIAAABUNGtCrGnv3r16/fXX1adPnyIDbFZWlrKystzvMzIyKqJ5AAAAKGfWlBNI0u233664uDg1aNBAmzZt0rx584rcftq0aUpISHC/WrZsWUEtBQAAQHkKa4idOHGifD5fkV8pKSnu9rfeeqtWr16tTz75RDVq1NCIESPkOE7Q40+aNEnp6enu1+bNmyvibQEAAKCc+ZyiUmA527Vrl/bs2VPkNu3atVN0dHSh9Vu2bFHLli315Zdfqnfv3iGdLyMjQwkJCeqnCxTpo44WgP0OOX/qc32o9PR01alTJ9zNAYCjlp/XivtcC2tNbGJiohITE0u1b15eniT51bwCAACgerDiwa6VK1fq66+/1mmnnaZ69epp48aNuuuuu5SUlBTyXVgAAABUHVY82FWrVi29//776t+/v4477jhdddVV6tSpk5KTkxUTExPu5gEAAKCCWXEntmPHjvr000/D3QwAAABUElbciQUAAABMhFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwDiEWAAAA1iHEAgAAwDqEWAAAAFiHEAsAAADrEGIBAABgHUIsAAAArEOIBQAAgHUIsQAAALAOIRYAAADWIcQCAADAOoRYAAAAWIcQCwAAAOsQYgEAAGAdQiwAAACsQ4gFAACAdQixAAAAsA4hFgAAANYhxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwDiEWAAAA1iHEAgAAwDqEWAAAAFiHEAsAAADrEGIBAABgHUIsAAAArEOIBQAAgHUIsQAAALAOIRYAAADWIcQCAADAOtaF2KysLHXp0kU+n09r1qwJd3MAAAAQBtaF2Ntuu03NmjULdzMAAAAQRlaF2I8++kiffPKJHn300XA3BQAAAGEUGe4GhGrHjh0aN26c/vOf/6hWrVoh7ZOVlaWsrCz3+/T0dEnSYeVITrk0EwAqVI6yJUkZGRlhbgkAlI38zzPHKTqsWRFiHcfRqFGjdPXVV6tHjx5KTU0Nab9p06Zp6tSphdZ/rg/LuIUAEF4tW7YMdxMAoExlZmYqISEh6Os+p7iYW44mTpyohx56qMht1q5dq08++URvv/22kpOTVaNGDaWmpqpt27ZavXq1unTpEnTfgndi9+3bp9atW2vTpk1FXhR4MjIy1LJlS23evFl16tQJd3OswXUrOa5Z6fzxxx9q06aNUlNTVa9evXA3BwCOmuM4yszMVLNmzRQREbzyNawhdteuXdqzZ0+R27Rr106XXXaZ5s+fL5/P567Pzc1VjRo1dMUVV+iVV14J6XwZGRlKSEhQeno6/5MMEdesdLhuJcc1Kx2uG4DqKqzlBImJiUpMTCx2u6eeekr333+/+/3WrVs1cOBAvfXWW+rVq1d5NhEAAACVkBU1sa1atfL7vnbt2pKkpKQktWjRIhxNAgAAQBhZNcTW0YqJidGUKVMUExMT7qZYg2tWOly3kuOalQ7XDUB1FdaaWAAAAKA0qtWdWAAAAFQNhFgAAABYhxALAAAA6xBiAQAAYJ1qHWIXLlyoXr16KTY2VvXq1dPQoUPD3SRrZGVlqUuXLvL5fFqzZk24m1Nppaam6qqrrlLbtm0VGxurpKQkTZkyRdnZ2eFuWqXzzDPPqE2bNqpZs6Z69eql//3vf+FuUqU1bdo09ezZU/Hx8WrUqJGGDh2qdevWhbtZAFChqm2Ife+993TllVdq9OjR+u677/TFF1/o8ssvD3ezrHHbbbepWbNm4W5GpZeSkqK8vDw9//zz+umnn/T4449rxowZuuOOO8LdtErlrbfe0oQJEzRlyhStWrVKnTt31sCBA7Vz585wN61SSk5O1rXXXquvvvpKixcvVk5Ojs455xwdOHAg3E0DgApTLYfYOnz4sNq0aaOpU6fqqquuCndzrPPRRx9pwoQJeu+993TSSSdp9erV6tKlS7ibZY1HHnlEzz33nH799ddwN6XS6NWrl3r27Kl//etfkqS8vDy1bNlS1113nSZOnBjm1lV+u3btUqNGjZScnKwzzjgj3M0BgApRLe/Erlq1SmlpaYqIiFDXrl3VtGlTnXvuufrxxx/D3bRKb8eOHRo3bpxeffVV1apVK9zNsVJ6errq168f7mZUGtnZ2fr222919tlnu+siIiJ09tlna8WKFWFsmT3S09Mlid8rANVKtQyx+XfA7rnnHk2ePFkLFixQvXr11K9fP+3duzfMrau8HMfRqFGjdPXVV6tHjx7hbo6VNmzYoKefflr/+Mc/wt2USmP37t3Kzc1V48aN/dY3btxY27dvD1Or7JGXl6cbb7xRp556qjp06BDu5gBAhalSIXbixIny+XxFfuXXKErSnXfeqYsvvljdu3fXrFmz5PP59M4774T5XVS8UK/b008/rczMTE2aNCncTQ67UK+ZKS0tTYMGDdKll16qcePGhanlqGquvfZa/fjjj3rzzTfD3RQAqFCR4W5AWbr55ps1atSoIrdp166dtm3bJkk68cQT3fUxMTFq166dNm3aVJ5NrJRCvW6ffvqpVqxYUWiO9h49euiKK67QK6+8Uo6trFxCvWb5tm7dqjPPPFN9+vTRCy+8UM6ts0vDhg1Vo0YN7dixw2/9jh071KRJkzC1yg7jx4/XggUL9Nlnn6lFixbhbg4AVKgqFWITExOVmJhY7Hbdu3dXTEyM1q1bp9NOO02SlJOTo9TUVLVu3bq8m1nphHrdnnrqKd1///3u91u3btXAgQP11ltvqVevXuXZxEon1GsmHbkDe+aZZ7p3/CMiqlQHyFGLjo5W9+7dtWTJEneYu7y8PC1ZskTjx48Pb+MqKcdxdN1112nu3LlatmyZ2rZtG+4mAUCFq1IhNlR16tTR1VdfrSlTpqhly5Zq3bq1HnnkEUnSpZdeGubWVV6tWrXy+7527dqSpKSkJO4CBZGWlqZ+/fqpdevWevTRR7Vr1y73Ne4yeiZMmKCRI0eqR48eOvnkk/XEE0/owIEDGj16dLibVilde+21euONNzRv3jzFx8e7tcMJCQmKjY0Nc+sAoGJUyxArHRnmKDIyUldeeaUOHjyoXr166dNPP1W9evXC3TRUIYsXL9aGDRu0YcOGQkG/Go5uF9SwYcO0a9cu3X333dq+fbu6dOmiRYsWFXrYC0c899xzkqR+/fr5rZ81a1axZS4AUFVUy3FiAQAAYDeK8wAAAGAdQiwAAACsQ4gFAACAdQixAAAAsA4hFgAAANYhxAIAAMA6hFgAAABYhxALAAAA6xBigTKwbds2XX755Tr22GMVERGhG2+8MeB277zzjo4//njVrFlTHTt21IcfflixDQUAoIogxAJlICsrS4mJiZo8ebI6d+4ccJsvv/xSf/vb33TVVVdp9erVGjp0qIYOHaoff/yxglsLAID9mHYWCMGuXbvUsWNHXX/99brjjjskHQml/fr100cffaT+/fu72/br109dunTRE0884XeMYcOG6cCBA1qwYIG77pRTTlGXLl00Y8aMCnkfAABUFdyJBUKQmJiol19+Wffcc4+++eYbZWZm6sorr9T48eP9AmxRVqxYobPPPttv3cCBA7VixYryaDIAAFVaZLgbANhi8ODBGjdunK644gr16NFDcXFxmjZtWsj7b9++XY0bN/Zb17hxY23fvr2smwoAQJXHnVigBB599FEdPnxY77zzjl5//XXFxMSEu0kAAFRLhFigBDZu3KitW7cqLy9PqampJdq3SZMm2rFjh9+6HTt2qEmTJmXYQgAAqgdCLBCi7OxsDR8+XMOGDdN9992nsWPHaufOnSHv37t3by1ZssRv3eLFi9W7d++ybioAAFUeNbFAiO68806lp6frqaeeUu3atfXhhx9qzJgx7mgDa9askSTt379fu3bt0po1axQdHa0TTzxRknTDDTeob9++euyxxzRkyBC9+eab+uabb/TCCy+E6y0BAGAthtgCQrBs2TINGDBAS5cu1WmnnSZJSk1NVefOnfXggw/qmmuukc/nK7Rf69at/coO3nnnHU2ePFmpqalq3769Hn74YQ0ePLii3gYAAFUGIRYAAADWoSYWAAAA1iHEAgAAwDqEWAAAAFiHEAsAAADrEGIBAABgHUIsAAAArEOIBQAAgHUIsQAAALAOIRYAAADWIcQCAADAOoRYAAAAWIcQCwAAAOv8P2WdkDrOyZ0TAAAAAElFTkSuQmCC", + "text/plain": [ + "Figure(PyObject
)" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig, axis = plt.subplots(2, 2, figsize=(8,8))\n", + "fig.suptitle(\"AdvancedHMC's NUTS - 21-D Neal's Funnel\", fontsize=16)\n", + "\n", + "fig.delaxes(axis[1,2])\n", + "fig.subplots_adjust(hspace=0)\n", + "fig.subplots_adjust(wspace=0)\n", + "\n", + "axis[1,1].hist(x10_nuts, bins=100, range=[-6,2])\n", + "axis[1,1].set_yticks([])\n", + "\n", + "axis[2,2].hist(theta_nuts, bins=100, orientation=\"horizontal\", range=[-4, 2])\n", + "axis[2,2].set_xticks([])\n", + "axis[2,2].set_yticks([])\n", + "\n", + "axis[2,1].hist2d(x10_nuts, theta_nuts, bins=100, range=[[-6,2],[-4, 2]])\n", + "axis[2,1].set_xlabel(\"x10\")\n", + "axis[2,1].set_ylabel(\"theta\");" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "fe4c8b70", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAArEAAAL3CAYAAACd2x1cAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABLxklEQVR4nO3deZyVZf0//vewDfsICgqyo5VboIL7gqJhUkqp2SdbMOWRftQyNMNyLfuoqdlHW8xvqVn5cclwQzPFQE3NBUkxNUBHBARZdAZQB2Tm94c/zrkH5owzMMOZa3g+Hw8ej2vO3Pd93ucM4mvu632uq6SmpqYmAAAgIW2KXQAAADSWEAsAQHKEWAAAkiPEAgCQHCEWAIDkCLEAACRHiAUAIDlCLAAAyRFiAQBIjhALm9GgQYOipKQkbrrppnqPGzVqVJSUlMRFF11U6/Fp06ZFSUlJ7s/zzz9f73V22WWX3LEnn3xywePWrFkTN954Y4wbNy4GDBgQnTp1is6dO8eQIUPi2GOPjT/96U+xevXqhr7MTfL222/HzTffHF/5yldixx13jI4dO0bnzp3jU5/6VHz729+O8vLyguc++uij8T//8z9xzDHH5N7rkpKSePzxxzeppuy1SkpKok2bNtGtW7fo169fHHLIIXH22WfH008/vdHXHz9+fO7ap59+esHjLrnkkigpKYnx48dv9HM1hXXvx/o/i3WvY9q0ac323Ov+26jvz5lnntlsz9+SrPv3YNSoUcUuBYqiXbELADbeDTfcENdee22d33vqqafi3//+98deY8aMGXHsscfG66+/HiUlJTFs2LDYa6+9ok2bNlFeXh533XVX3HnnnfHDH/4w/v3vf0fnzp2b+mXUMnHixPjTn/4Ubdq0iV133TWOOuqoWLVqVTzzzDNx7bXXxg033BCTJ0+Oww8/fINzv/3tb8e//vWvZqtt//33jx122CEiIt5///1YunRpPP/88zFt2rS46qqr4uCDD44bbrghhgwZstHPcf3118d3v/vdGDp0aFOV3SoNGzYshg8fXuf39tprr81bDFAUQiwkaMCAAfHBBx/ELbfcEldeeWWUlpZucMwNN9wQEREjR46MZ555ps7rzJgxIw488MB477334nOf+1xcc801MXjw4FrHLFmyJK6++uq46qqrYvXq1c0eYnv27BkXX3xxnHTSSbH99tvnHl+5cmVMmDAhbr311vjyl78cc+bMiR49etQ69/DDD48vfOELsccee8Qee+wR+++/f7zxxhtNVtvJJ5+8wV3QmpqaeOCBB+LMM8+M6dOnx3777RdPPvnkBu9jQ3Tu3Dnee++9+OEPfxi33nprE1XdOo0bN26DmQpgy6KdABLUvn37+OpXvxrLly+Pu+66a4Pvv/fee3HrrbfG9ttvH2PGjKnzGmvWrInjjjsu3nvvvRg3blzcfffddQavXr16xf/8z//E448/XmdYbmrXXHNNXHDBBbUCbERE165d43e/+11069Ytli9fHlOmTNng3CuuuCIuvPDC+PznP7/B+c2lpKQkjjzyyHj66adjxx13jMWLF9fbulGfk046Kbp27Rq33357zJgxo4krBWhdhFhI1De/+c2IyN9xzbrjjjtixYoV8fWvfz3atm1b5/m33HJLvPbaa9GhQ4f49a9/HW3a1P/PwciRI6NTp06bXvgm6Ny5c3zyk5+MiIg333yzqLWsb6uttoqf//znERHxyCOPxHPPPdfoa/Tu3TvOOuusqKmpie9///uNPn/hwoUxceLE2GmnnaJz587RrVu3GDlyZPziF7+IDz/8cIPjlyxZEtdcc00ceeSRMXjw4OjUqVN07949RowYEZdffnl88MEHja6hLlVVVXHFFVfEnnvuGd26dYsOHTrEdtttFyNHjoxzzjknli9f3iTPs76LLrqozt7ydQr1lGYfX7NmTVx++eWxyy67RKdOnWLrrbeOL37xi/Hyyy9vcL3y8vIoKSmJQYMGRU1NTVx//fWx5557RpcuXaKsrCw+85nPxJNPPlmw3vfffz+uuuqq2GeffWKrrbaKjh07xic/+ck455xzYtmyZZvyVkCrJMRConbZZZfYa6+94uGHH94g0P3ud7+LiIgTTzyx4Pl33313RESMGTMmtttuu+YrtAmtWbMm92GiPn36FLeYOnz2s5+Nnj17RkTEQw89tFHXOPvss6N3797x8MMPN+oajz76aOy6665x9dVXxwcffBCHH3547L///jF37tw444wzYuzYsbFmzZpa5zz44IPxne98J1544YUYOHBgjBs3Lvbaa6949dVXY9KkSXHooYdGVVXVRr2Odaqrq2Ps2LFxzjnnxJw5c+LAAw+MY489NnbbbbdYsmRJXHHFFTFv3rxNeo7msmbNmjjyyCPjRz/6UQwYMCDGjh0bXbp0icmTJ8d+++1X74cMTzzxxDj99NNjq622is997nOx3XbbxUMPPRSHHHJI/POf/9zg+IULF8bee+8dZ599dsyePTtGjhwZRx55ZO4XgBEjRjRpawy0BkIsJOykk06K6urqWqsdzJ49Ox577LE46KCDYscddyx47rPPPhsRH91hTcXvfve7WLp0aXTq1Ck++9nPFrucDZSUlMQee+wREREvvfTSRl2ja9eucd5550VExKRJk6KmpuZjz1m0aFF88YtfjHfffTd+9atfxdy5c+Puu++Ov/71rzF79uw49NBD429/+1tceumltc7bc88948knn4z58+fHtGnT4v/+7//i4Ycfjnnz5uXuGl5zzTUNrv2mm26KmpqaWnc2H3/88Zg6dWrsvvvuMX/+/Lj//vvjlltuiYceeijKy8vjmWeeif79+zf4OTanJ554IpYtWxZz586NBx54IP785z/Hf/7znxgzZky8++67G7yf67zxxhsxbdq0mDVrVkydOjVuu+22+Pe//x3f/OY3o6qqKi644IJax9fU1MSXvvSlePHFF+Okk06K8vLyeOihh+Ivf/lLzJkzJ84666woLy+v95dS2BIJsVAEJ554Yr1LBE2fPr1B1/nyl78cnTt3zoWHiHx7wbp2g0KWLFkSER9NYafgxRdfjO9973sREXH++efHtttuW+SK6rbNNttERGzS9O8pp5wSQ4YMiRkzZsRtt932scf//Oc/j2XLlsVpp50Wp556aq3WkK233jpuvvnmaN++ffziF7+oFYp32mmn2GeffTa4Xo8ePXKrXtxxxx0b/ToiIhYvXhwREQceeGB069Ztg++PGDEitt5660Zf9+KLL67zv51BgwZtUr1ZJSUlceONN9aaqejYsWNcfPHFERHx8MMPFzz32muvjU984hO5r9u2bRs/+clPIiJi+vTpte6KP/jgg/GPf/wjhg8fHtddd12t96ldu3bx05/+NHbdddf4+9//HrNmzWqy1wepszoBFEF2qaa6/PWvf839z78+3bt3j2OOOSb+8Ic/xLRp0+Kggw6Km2++Obp16xbHHXdcU5ZcVPPnz4/Pf/7zsXLlyjjqqKNi0qRJxS6poOrq6oj4KABtrPbt28cll1wSX/nKV+K8886LY445Jtq3b1/w+HUfcjv++OPr/P72228fO+64Y/z73/+O2bNn1wpXa9eujWnTpsUTTzwRb731Vrz//vtRU1OTC7uvvvrqRr+OiIg99tgj2rZtGzfccEN84hOfiC9+8YtN0gpSaImtdb9ENIUBAwbEsGHDNnh8p512ioiIBQsW1Hleu3bt4ogjjtjg8e222y569OgR77zzTixbtiwXjtf9/I455pho127D/y23adMmDjrooJg1a1Y88cQTseuuu270a4LWRIiFIqhrqaasUaNGNSjERnx0x/UPf/hD3HDDDfHee+/FwoUL4+STT/7YpbB69eoVb775Zrz99tuNKf1j3XXXXXWumHDyySfHAQcc0OjrLVq0KEaPHh1vvPFGjBkzJm6//fZNCojrq+vnsM0228SVV165UddbunRpRESuN3ZjffnLX44rrrginn/++fjNb35T7yYIr732WkR8dLfz4yxZsiQXYmfPnh1f+MIX6m19qKysbGTltQ0dOjSuvvrq+N73vhenn356nH766TFw4MDYd99943Of+1wcd9xx0aFDh0Zfd3MssTVgwIA6H+/evXtERMF+4T59+hT8paN79+7xzjvv1PrQ3Lqf3/nnnx/nn39+vTWtm0EBhFhI3sEHHxxDhw6NO++8MxYtWhQRH99KEPFRP+Sbb75ZcA3ZjTVz5sz4/e9/v8Hjo0aNanSIffvtt+PQQw+N//znP3HYYYfFXXfd1eTLfNVV68CBAzcqxNbU1OR2Udttt902qa6SkpK47LLLYsyYMfHjH/+43l961t39PfbYY6NLly71Xjc7dX/sscfGSy+9FJ/73OfinHPOiZ133jm6d+8e7du3j9WrVzfZe33GGWfEl770pbjnnnvi8ccfj8cffzxuvfXWuPXWW+PCCy+Mxx57rCgf1Fv3vhXycSt2NNV56+o44IADPnaTi1122WWjaoLWSIiFxK3bhvT888+Phx9+OHbaaafYd999P/a8o48+Ou6666548MEHY/HixU3WY3rRRRc1yR2yJUuWxKGHHhovv/xyjB49Ou65557o2LHjphe4noZ8cKqh7r///njnnXciIuIzn/nMJl/vM5/5TIwePTqmTp0aV111VcHl0vr37x+zZ8+O73//+zFixIgGXfuVV16JF154IXr37h2TJ0/eYBp79uzZm1x/1rbbbhsTJkyICRMm5J7/m9/8Zjz55JMxadKkOn+Z2FTr7vCuWLGizu+3lE/7r/tg29FHHx1nn312kauBdPhgF7QC48ePj169esXWW28d3/rWtxp0zgknnBCDBg2K1atXx6mnnvqxd6Wee+65eP/995ui3I+1dOnSOPTQQ+Oll16K0aNHx7333lv0NWo/TkVFRXz3u9+NiI92Diu0JWpjXX755VFSUhJXXXVVwankdSs13H777Q2+7rq1Wfv27VtnH+Yf//jHjai24T71qU/l1sKdOXNmszzHug0v6lrTNSLq3DCjGNb9/O64444m/aUKWjshFlqBfv36xdtvvx1Lly6N73znOw06p3379nH77bdHx44dY/LkyTFu3Lh4/fXXNzhu+fLlcf7558f++++/yWuGNsTy5ctj9OjRMWvWrDjssMNafIBdt+3sXnvtFbNnz44+ffrE//t//6/Jrr/nnnvGcccdFytWrIjf/va3dR7zve99L7baaqv42c9+ltseeH2vv/56rWD6iU98Itq2bRsvvvhiTJs2rdax9957b1x99dVNUv8jjzwS999//wZr1NbU1MR9990XER+1bzSHQw89NNq0aRMPPvhgrRU/ampq4pprrok777yzWZ63sY4++ugYOXJkPP3003HiiSfW+cvKO++8E9ddd12dm1bAlko7AWzBRo4cGY8++mgcd9xxce+998Z9990Xu+++ewwZMiTatGkTb7zxRjz77LOxdu3aGDJkyGbZdvbkk0+OF154IUpKSqJnz55x6qmn1nncuHHjYty4cbUe++1vf1sr6L311lsREfGtb30rt2xRnz59YvLkyRtV229/+9tc4KuqqoqlS5fGjBkzcnc1R40aFTfccEOTh7Kf/OQnMXny5Hjvvffq/H6/fv3i7rvvjmOOOSbOPvvs3JJMffr0iYqKinj55Zdj7ty5sffee8dXv/rViPjow2unn356/O///m+MHj06DjzwwOjbt2+8+uqrMWPGjDjvvPPikksu2eTaX3jhhfjud78b3bt3jz322CP69u0b77//fsyYMSPeeOONKCsrix/96Eeb/Dx16d+/f5xxxhm1XmPPnj3jX//6V8ybNy8mTZoUl112WbM8d2O0adMm7rrrrhg7dmz8/ve/jz//+c8xbNiwGDBgQKxevTpee+21ePHFF2Pt2rUxfvz4Ou+cw5bIfwmwhRs5cmT85z//iT/84Q9x9913x4wZM+Kll16KkpKS6NOnT3zhC1+IY4455mOXeWoq6wJhTU1NvdPjgwYN2iDEzp8/v87dkP7973/nxpsSMP/xj3/EP/7xj4iI3Faiu+22W4wYMSKOP/74Zts4YocddogJEybEr371q4LHHHTQQfHSSy/FL37xi5gyZUo888wzUVVVFb17944BAwbEV7/61TjmmGNqnXP11VfHpz/96fjVr34Vzz33XMycOTN22223uPXWW+P4449vkhD7+c9/PioqKuKxxx6L2bNnx1NPPRWdOnWK/v37x6RJk+K0006Lfv36bfLzFHL11VfHgAED4re//W088cQT0bVr19h///3j9ttvj8rKyhYRYiM+aut46qmn4qabborbbrstXnjhhXj66aejZ8+e0bdv3zjllFPiqKOOapa+cEhVSY0GHAAAEqMnFgCA5AixAAAkR4gFACA5QiwAAMkRYgEASI4QCwBAcoRYAACSI8QCAJAcIRYAgOQIsQAAJEeIBQAgOUIsAADJEWIBAEiOEAsAQHKEWAAAkiPEAgCQHCEWAIDkCLEAACRHiAUAIDlCLAAAyRFiAQBIjhALAEByhFgAAJIjxAIAkBwhFgCA5AixAAAkR4gFACA5QiwAAMkRYgEASI4QCwBAcoRYAACSI8QCAJAcIRYAgOQIsQAAJEeIBQAgOUIsAADJEWIBAEiOEAsAQHKEWAAAkiPEAgCQHCEWAIDkCLEAACRHiAUAIDlCLAAAyRFiAQBIjhALAEByhFgAAJIjxAIAkBwhFgCA5AixAAAkR4gFACA5QiwAAMkRYgEASI4QCwBAcoRYAACSI8QCAJAcIRYAgOQIsQAAJEeIBQAgOUIsAADJEWIBAEiOEAsAQHLaFbuAzam6ujoWLlwY3bp1i5KSkmKXA7DJampqYsWKFdG3b99o08Z9CWDLsUWF2IULF0b//v2LXQZAk3vzzTejX79+xS4DYLPZokJst27dIuKjf+y7d+9e5GoANl1lZWX0798/9+8bwJZiiwqx61oIunfvLsQCrYoWKWBLo4EKAIDkCLEAACRHiAUAIDlCLAAAyRFiAQBIjhALAEByhFgAAJIjxAIAkBwhFgCA5AixAAAkR4gFACA57YpdALD5DZo0JTcuv2xsESsBgI3jTiwAAMkRYgEASI4QCwBAcoRYAACSI8QCAJAcIRYAgOQIsQAAJEeIBQAgOUIsAADJEWIBAEiOEAsAQHKEWAAAkiPEAgCQHCEWAIDkCLEAACRHiAUAIDntil0AUFyDJk3JjcsvG1vESgCg4dyJBQAgOUIsAADJEWIBAEiOEAsAQHKEWAAAkiPEAgCQHCEWAIDkCLEAACRHiAUAIDlCLAAAyRFiAQBIjhALAEByhFgAAJLTrtgFAC3ToElTcuPyy8YWsRIA2JA7sQAAJEeIBQAgOUIsAADJEWIBAEiOEAsAQHKEWAAAkiPEAgCQHCEWAIDkCLEAACRHiAUAIDm2nYUtRHYbWQBInRALrZzwCkBrpJ0AAIDkCLEAACRHiAUAIDlCLAAAyRFiAQBIjhALAEByLLEFrZBltQBo7dyJBQAgOe7EQsKyd1zLLxtbxEoAYPMSYoEGE5oBaCmEWGglmqIPVi8tAKnQEwsAQHKEWAAAkiPEAgCQHCEWAIDkCLEAACRHiAUAIDlCLAAAyRFiAQBIjhALAEBy7NgFfCw7eQHQ0rgTCwBAcoRYAACSI8QCAJAcIRYAgOQIsQAAJEeIBQAgOUIsAADJsU4ssMmy68iWXza2iJUAsKVwJxYAgOQIsQAAJEeIBQAgOUIsAADJEWIBAEiOEAsAQHKEWAAAkiPEQiIGTZpSaz1WANiSCbEAACRHiAUAIDlCLAAAyRFiAQBIjhALAEBy2hW7ACBNVkoAoJjciQUAIDlCLAAAyRFiAQBIjhALAEByhFgAAJIjxAIAkBwhFgCA5AixAAAkR4gFACA5QiwAAMkRYgEASI4QCwBAcoRYAACS067YBQCNM2jSlGKXAABF504sAADJEWIBAEiOEAsAQHKEWAAAkiPEAgCQHCEWAIDkCLEAACRHiAUAIDlCLAAAyRFiAQBIjhALAEByhFgAAJIjxAIAkBwhFgCA5AixAAAkR4gFACA5QiwAAMkRYgEASI4QCwBAcoRYAACS067YBQCty6BJU3Lj8svGFrESAFozd2IBAEiOEAsAQHKEWAAAkiPEAgCQHCEWAIDkCLEAACRHiAUAIDlCLAAAyRFiAQBIjhALAEByhFgAAJIjxAIAkBwhFgCA5LQrdgFAbYMmTcmNyy8bW8RKNt2615L66wCg5XEnFgCA5AixAAAkR4gFACA5QiwAAMkRYgEASI4QCwBAciyxBS1YdrktACDPnVgAAJIjxAIAkBwhFgCA5AixAAAkR4gFACA5VicAml12lYXyy8YWsRIAWgt3YgEASI47sdACWA8WABpHiIUiEVwBYONpJwAAIDlCLAAAyRFiAQBIjhALAEByhFgAAJIjxAIAkBwhFgCA5AixAAAkR4gFACA5QiwAAMkRYgEASI4QCwBAcoRYAACSI8QCAJAcIRYAgOQIsQAAJEeIBQAgOUIsAADJEWIBAEhOu2IXAFuSQZOmFLsEAGgV3IkFACA5QiwAAMnRTgBsVtmWivLLxhaxEgBS5k4sAADJEWIBAEiOEAsAQHKEWAAAkiPEAgCQHCEWAIDkCLEAACRHiAUAIDlCLAAAyRFiAQBIjhALAEByhFgAAJIjxAIAkBwhFgCA5AixAAAkR4gFACA5QiwAAMkRYgEASI4QCwBAcoRYAACSI8QCAJAcIRYAgOQIsQAAJEeIBQAgOe2KXQBsCQZNmlLsEgCgVXEnFgCA5AixAAAkR4gFACA5QiwAAMkRYgEASI4QCwBAciyxBc3EslobL/velV82toiVANBSCbFA0QirAGws7QQAACRHiAUAIDnaCaAJ6YPdeN47ABrDnVgAAJIjxAIAkBwhFgCA5AixAAAkR4gFACA5QiwAAMkRYgEASI4QCwBAcoRYAACSI8QCAJAcIRYAgOQIsQAAJEeIBQAgOUIs0KINmjQlBk2aUuwyAGhhhFgAAJIjxAIAkBwhFgCA5AixAAAkR4gFACA5QiwAAMkRYgEASI4QCwBAcoRYAACSI8QCAJAcIRYAgOQIsQAAJKddsQuAFA2aNCU3Lr9sbBErAYAtkzuxAAAkR4gFACA5QiwAAMkRYgEASI4QCwBAcoRYAACSY4ktaITs0loAQPEIsbCJBFsA2Py0EwAAkBwhFgCA5AixAAAkR4gFACA5QiwAAMmxOgGQhOwqEOWXjf3YxwFo3dyJBQAgOUIsAADJEWIBAEiOEAsAQHKEWAAAkiPEAgCQHCEWAIDkCLEAACTHZgfwMbKL6dMy+JkA4E4sAADJEWIBAEiOdgL4/5miBoB0uBMLAEBy3IkFWo3s3fTyy8YWsRIAmps7sQAAJEeIBQAgOUIsAADJ0RPLFs2KBACQJndiAQBIjhALAEByhFgAAJIjxAKt0qBJU/Q8A7RiQiwAAMkRYgEASI4QCwBAcqwTyxZHnyQApE+IZYshvAJA66GdAACA5LgTC2wxsnfjyy8bW8RKANhU7sQCAJAcd2JpFdxhA4AtizuxAAAkx51YWjUrEuDvAEDrJMTS6ggtAND6aScAACA5QiwAAMnRTkBy1rULWIWApmJ1C4D0JHMn9tJLL42RI0dGt27donfv3jFu3Lh49dVXi10WAABFkEyInT59epx22mnx1FNPxUMPPRRr1qyJz3zmM7Fq1apilwa0UoMmTan3g4Lrvu/DhACbXzLtBH/9619rfX3TTTdF796947nnnouDDjqoSFXRnAQDmlOhv191Pa7dAKDlSeZO7PoqKioiIqJnz55FrgQAgM0tmTuxWdXV1XHmmWfG/vvvH7vuumvB46qqqqKqqir39TvvvBMREfPnz4/u3bs3e51saJ//mRoREU/9YPQGj63/+IeVS+u9Vr///n0TVwcfb/78+blx9u9o9u9j9u9xc3v33XcjImL58uWb7TkBmlNNTU2sWLEi+vbtG23aFL7fWlJTU1OzGetqEqeeemo88MAD8fjjj0e/fv0KHnfRRRfFxRdfvBkrAwCgKbz55pv15rzkQuzpp58ed999dzz66KMxePDgeo+t607soEGD4oA4MtpF++YuFaDZvR+r4p/xcLz55ptmmIBWobKyMvr37x/vvvtulJWVFTwumXaCmpqaOOOMM2Ly5Mkxbdq0jw2wERGlpaVRWlq6wePton20KxFigfS1r+kQERHdu3cXYoFWpaSkpN7vJxNiTzvttLjlllvi7rvvjm7dusWiRYsiIqKsrCw6depU5OoAANicklmd4Ne//nVUVFTEqFGjok+fPrk/t912W7FLAwBgM0vmTmxirbsAADSjZO7EAgDAOkIsAADJEWIBAEiOEAsAQHKEWAAAkiPEAgCQHCEWAIDkCLEAACRHiAUAIDlCLAAAyRFiAQBIjhALAEByhFgAAJIjxAIAkBwhFgCA5AixAAAkR4gFACA5QiwAAMkRYgEASI4QCwBAcoRYAACSI8QCAJAcIRYAgOQIsQAAJEeIBQAgOUIsAADJEWIBAEiOEAsAQHKEWAAAkiPEAgCQHCEWAIDkCLEAACRHiAUAIDlCLAAAyRFiAQBIjhALAEByhFgAAJIjxAIAkBwhFgCA5AixAAAkR4gFACA5SYXYRx99ND7/+c9H3759o6SkJO66665ilwQAQBEkFWJXrVoVw4YNi1/+8pfFLgUAgCJqV+wCGuOzn/1sfPazny12GSSk7a6fzI3Xzno1N263be/c+MPFb2/Sc2SvtfpT/XLjNtNn5MbVB++RG7dftio3LlnyTp11FKrv/S/sXeu5uz3xep3P3eGV+Y26bqHXU+jc+jTk/dzU979QLZv6swQgHUmF2MaqqqqKqqqq3NeVlZVFrAYAgKaSVDtBY1166aVRVlaW+9O/f/9ilwQAQBNo1Xdizz333Jg4cWLu68rKSkF2C1BrqjkzXZ9V06tH/osCU9DrT92vk53Cj4hYsd/g3Ljr7Hfz38i0MrQt0EKwZOzQ3LjjO/nrlC5fkxvP/+8dcuM1ZdW1nrtnj/z5VWUluXHvyLcWRKbNYEXP9pmz66579dZdcuMOmaOz7QoRhdsiFn13v3wdMz6o8/jszyXbapFtg8iq9fOKiA8zrSEN0ZTtI5ty3eaqA2BL1KpDbGlpaZSWlha7DAAAmlirbicAAKB1SupO7MqVK2POnDm5r19//fWYOXNm9OzZMwYMGFDEymguDf2kfKGp2ew0dHVmOjw7vZ9dwSA7LZ5tG8heJ9sCEBHR44Yn819krlVIdlq+4zvVdR5TlZn27/dIgSn5iFi6V8/cuKx8bZ3nf9Aj/7tqtuUga8nwrXPjbm/U5K+ZaUtY/7kLtUVs/9CyOo8pZP3rNkTbAu9zY1eg2NTp/caeo4UAoOkkFWKfffbZOOSQQ3Jfr+t3/cY3vhE33XRTkaoCAGBzSyrEjho1Kmpqaj7+QAAAWrWkQiytV6Fp3ULTr+s/Xmjx++z0crYBfG1m3NjF/bd5enmtr0sKrIaQXbVgyfC2ufHgO/LT7e3zw1pePy4/vd/h3fy5pRW1P6iYbQ/omJm5z6420DVz/JrMygNVtVYqyMuujJC1csetan2dvW6htojse1CoPSMr+3j2Z7f+ChKFVjRobHuA6X2AdAmxAK3Arhc+GG1KOxe7jFrKLxtb7BKAVszqBAAAJMedWFqEhkzjN+SY9Y9rqufOWrveQvvZTRGy0/jZVQGyKwxkp/QrhubbA8rm5rdI7pC/TJRW5PvA119dILuZwOtfyG9N0De2yo0XHpw/p+/0/LWy9WWVZrol5h/aMTfOrloQERGZ9oLstbplWio+yKye0CnzPr9Ta5OHfCtCtuUg+7PYoP1g+ozc8MO6XkQDFVqZQptB0xg0aUqxSyDcEaf1cicWAIDkCLEAACRHOwEtzqZ+qrzQQvjZ6eKCn44vcM0Fh+dXC+i99R61vpdtIaj16f/M1H+htoHsSgfZFQl6zcyun5C3YmDbWl9nVxjo8kb+d9IlwyPzePaMzIYImfqyLQtv71F3C8GqsStqPfc2F+ZXMci2DWQ3c1j+6fz5Ky7eLzeubwOHuqy/aUJDWk4adPx6rSF1HV/fShjaDkjBx7V1aDcgVe7EAgCQHCEWAIDkaCegRcu2Bqy/KkAh6089r5NtISjUWlDo8X63zCn4fNnp82zbQHblgOyU+dLM1HtZ5jrZKfbsdbKf3q8asLrWc5c+kv/eHt9/OTee/vIncuMO7+ZXLZh/VL6dYJvH8r/DrhhYkjk+f/1sy8HWf8y3Snz0OvLnZ1sksi0Ogyfn681uSjD/Kzvkzy3PH1+RadvY/qH8ThDr/0yz73n2vc3+g5ad6i/0M94Ym9JCoBWB5qIlgC2RO7EAACRHiAUAIDnaCQAgcS1pYwmtDWwuQiwtWqE+2I3asSvzePXB+WWy2mR2f8ruvpXdPSrbU5ldRmv982PovnU+de0+zHxPbLZ/NNuv2iaz5FX7ivyyWqXzai+xNXd8vuf07dt2zY23yZzf8Z38dfc4/j+58dPdB+TGX//EM7nxI6fvn7l+frJm+4dqL7GVXVYrK/s6+t2Tf32vnzskNx48Od//m+2V7TQ5/zNamflZrP9c2R7cQn3MbQsso1Zoh7bs483Vr6oPFqDpaCcAACA57sQCAE2mJbU2rKPFoXUSYmkWzbGUUH0tBA2RXa6rXWYKOzLXLV2+JuqSbSHYYIepzHV7TZmbG6/Yb3BuXDF2aG6cXTKrYlC+PaDbC/lxWXl+Sj67+1abXStrPXXprO658aqB+euWvpBfGuuzP5qWG9/8l9H558vsxvW7Nw/JjQec91Zu3O9X2+XG2V3LIgrv8jX0pg8yR+Vr6vJG/phsG0VV5n3qOjvfAtBp8j9z4yWZ3b4iInpNyf/8ssttZVs7Ci3Pln08u2OX5a8A0iLEAgCtmq13Wyc9sQAAJMedWJpFoenYQlO2DZnK3dQp3uyn2LPXyk4vZ1sFsp90z05Tr8msbLD+OdmdqLI7dlUMyk+ld8wsVLB6q7prXTI831qQ3UHrw/WOy+7gVTovvzNXVWYrsFtvyLQQZFoAsrt0HbjPS7nxY0/tnBt3GVT499ylB9bderFqYL6OvtPzz5d9rR/0yF93m6eX58avnNUtNx56U/59HvKrwjumZVc3qMn8LAsptOJFQ/7O1ndcY2lfANg0QiwAsEXbHB9G07LQ9LQTAACQHHdi2ayaqlWg7XrTxoU+fV57k4G87FRuTZ1H1D43e83q9Y5bmlmIPztlnv3Efu3VBmpvWLBOn/0W5MZvPbF9ncdUvdm11tdDJ+fbCbIbH3R4Nz+l/+EB+RUNSqfkp+v3GTMrNz5663y7xFMV+U0TVnw6f81+99Sue8Wn8+NarQyZFoeus/MbJLT773x7xbJ386/vPyfn2zaG3pQ/N2v9vx/ZzSqy7QTrb0SxTva39ezPvtDKBvVN9TdVG4AWAoBN404sAEAz0krQPIRYAACSo52AJBVqE4go/OnzrOyU8Modt8qNuz3xem6cXZ2gPtkNDjpmNzjIbGTQdfa7mcfzmwZkNyiY91Kf/EXL8o+32XVlblxdmW9RiIhY+O3M9Hvme1/+5tTcOLvBQcd38m0NTy8YkBt/osvi3Djb1lA+v1duvH4bxIQR0/LPMS+zAsIL+daCNy7Ob2qw9a+2yo1XD89fp2dmY4bsJgjZjSfWbx9pW2AViawOBTa0yE7jt9mI9hZtAAAtgxALANCMmmv1gy29TUE7AQAAyXEnli1GrSnpTDtCdqq/1gYMmXNXFGgTiIjYbvrbme/VvdlB9lPzhT7xn90AoKosP8Xee3J+ev71L9T+T7bLY/nVBqo+nV9nIdtCkN0s4dFfXp8b7/LkCbnxH/8zMjfe+o/5WvOf3Y+IWFvrq1vfzT9Hm8wKCPF493xNmdUUSpdnV0/Iv+5sO0ah1QJivXaC7PuZ3WyiVgtBM7FJAUDLIMQCACRoc2zSsLE2R6uDdgIAAJLjTixJ2pgF6AutaJB9vNB1Ok3Oj0szC+2vf052U4Os7DT3p67Kf5q+1qf3M9P4VWX56faKoaW5cZc38m0GHx2XH580+pHc+D+rts2NsxsZ7HD7t+qs76TRf8+Nf3fwIblxaf/8yghdMhslRNReWaFLpoUg2y7RpiL/T8zrX8i3RfSdnn+fsq0a2RUJ5l28X248+I5ltZ47+342ZEq/1mYVhVoWMrQMAGyaQZOmNPvdWHdiAQBIjhALAEBytBOQpOx0b0ThKd+GTAtXZ9oDCk0v13dMTWaqOrvSQfYT9K/99w5Rl15/rLv9YPVW+XFpRX6cnaqPiOjUPb8Cwu+mHhJ1eSx2zo3bV2RWQBiQv9atN+RXGrjyv2/Ojc+95ev55x5Y+7rVZfmp/9Vb5VsFtnmsfdRl6YH547vOXpEbZzebyOr3yAcFj+lWoAUku/FBdtOLbMtIVeZa2dYQKxsAm8uWvr5rU3EnFgCA5AixAAAkRzsBLVqh6deGrk6QnV6OAtO32faA6gLTy9lJ/Or1Vidom1lsv9C0dbc3aqIu2U/jv71Hx9y418y62wzWl10xYNv/Wpgbl8/vlRv/7IBbc+NryvNtAwO75Wt9at6uufHLH2yfG7fZNb+JwfuV+foiIkrn5VsI+uy3IDde+Ubf3Hh5ZgOGoTflVzNYcPjWuXG2daLfI/n3I/tz6breZgeFVjTInlOozaDT5H/WeUxDp/cbclxD/94CW6bNub5ra25dcCcWAIDkCLEAACSnpKampu55zhbql7/8ZVxxxRWxaNGiGDZsWFx77bWx1157NejcysrKKCsri1FxdLQrqfsT1KRh/dUJsrJTts3xafD1n3t1ZvH8bAtCdgo7u1JBdvOCFQPzmxdkWw6qympvarBOdoOBhsquSJBtD6h6s2tuXF32YW6cXfFg/RaCrGw7wT5jZuXGb6zIv+55L/XJjbu8ka8juylExaD8xg7Z1oL63o/trn4iN87+PAr9LLIa2w6wMedvTh/UvBePx/3R/8zbo01p52KXAxRJa2obWJfXKioqonv37gWPS+pO7G233RYTJ06MCy+8MGbMmBHDhg2LMWPGxNtvt6z/qQAA0LySCrE/+9nPYsKECXHiiSfGzjvvHNddd1107tw5brjhhmKXBgDAZpTM6gSrV6+O5557Ls4999zcY23atInDDjssnnzyyTrPqaqqiqqqqtzXlZWVdR5Hy9KQFoBN/SR5oeco9In27KoDtVY8iPU2P8icn20hqOqZb1/p+E6+JaBsbv6T9dk2g9KK/FR6rylzc+PsdHlExMJv59dNyK5UkN1YICry0/7Vs/LTMgdmWgCyftDngdz4pDMn5sZDv/9yreOe7j4gN37qwfzqBtlNFLZ5Id8G0PGdfAtBtydez42XDM9vBJHd4KBQa0ZEREnmfZ7/lfz52TaFDvHx2mZWPchujpDV0toHAOqyqSsepNiOkMyd2KVLl8batWtj2223rfX4tttuG4sWLarznEsvvTTKyspyf/r37785SgUAoJklcyd2Y5x77rkxcWL+TlJlZaUgCwBskVK821qfZELsNttsE23bto3FixfXenzx4sWx3Xbb1XlOaWlplJaW1vk9Wq7NMX2bnZ7O/kdQq7WgwPHry05Jr860EGSnw9tkrvv+F/bOjbMtBNs8vTw3zrYiZGua/9/5qfOIiKrK/GTK+5kWggkjHsuN/9h9ZP74zIoEWTNuy7cDHLVVfrzP9/MtB0dvPaPWOYtPzq88sGbrfBtAtsVh1dj841WPZz9hmt+sINtCUOv9WJZ/z1dmNo6IiIjM1/1umZMb19roYr0NEtap9fMu0EIAQMuXTDtBhw4dYs8994ypU6fmHquuro6pU6fGvvvuW8TKAADY3JK5ExsRMXHixPjGN74RI0aMiL322it+/vOfx6pVq+LEE08sdmkAABultU3zby5Jhdjjjz8+lixZEhdccEEsWrQohg8fHn/96183+LAXAACtW3I7dm0KO3Zt2RqydFd9OzU1RKHe2ewSXVmFlodaeHB+aaq+0/P/iX7Qo3YHUHa5rmVfXZUbb9t9Rf7x+7fPjRuyI9aHB+SXosteZ3FlfgmviIh2mR7XFZ/O98F2eyG/uFVDni+7nFi2L7g+2fcz+x5md/KqVWuBn2trWD7Ljl2QPndia2uVO3YBAECEEAsAQIKS6omFTVFo6ji7RFYUmPYvtKvX+rJLQZUuzy95VZV5vOvsd3PjbAtB9vHBy/NLbLVflm8TWHZx7TaY7LJV2en9ZZEfZ9sD+mTaA96KfJvBmrJ8W8KVn74rNz57ylfzr6f/yiikU/f8MlmlFdka820DKwbmx9nWgrK5+V31sm0C2d3Jau2KFrV/Ztn3sPrgPXLj7PtW6OdaSENaT5riHAA2nhALANBA+ldbDu0EAAAkx51YtnhrG7BrU31TxdnvdXvi9dw4u1JBh1fy09lLxg7NjbNT7KUFWgiyU+xb/zG/01VExAc96p6Wf3uPjrlx9axMm8G7+XF+DYGIfcb8Oze+8Fdfz9eUaUVYX7ZNIfscZXPzrQX9Ls/vpvXS9bvUeZ1srdtndunKvgcrM7ucRdRu1ciOs7ukbYqNaQfQQgCweQmxAAANNGjSlAYfq/WgeWknAAAgOe7EQgPUtzpB9nvZT80XalPo+E6+JaDXlHz7QbbNoFdmWnzFfrVbCLKWf7omM843CPR8If/4ik/np9u7vVH3Jh/TX/5EbtxmYH6lgp5T8hscZDcoiKi9ScGqsfnWgvmZlRHmP7hrbjxkSr614LX/zm9QkF0ZYcHhW+fG/W7JH1+ReTwiIga1zQ17z/gg6pL9ubyfaUfItnw0JasTQBrcHW093IkFACA5QiwAAMnRTsAWqdDUb/bx7OoC2daA+qaKs8cVeo7spgbZx7d5Ov982YX+KzJT59kp/IiIwZPrXpGg4zv5DQCG3pSfrl973sLceOX/9c2ND97pP7nxY0/tnBsvPTCzCsC87HoGEb1n5J87u+lCdvOBrGxbxIALn8iNF313vzrPnf+VfMvB9g8tq3Wt7IoNtX4emfe80KoRWdljsu95dmWEhqxesUEdQIujjaD1cScWAIDkCLEAACRHOwFbpEJTv7Ueb+D0cEM+lV7o0/HVB++Rv05mRYIOmenyVV8Ykhv3u/S1WtfNrmiQld0AoKpnfkWC0kvyn/Lv9crc3HjuOzvlxj17ZFchyJ9bVVb7ObJT7qUVpbnxwoPz51eXfZiv/Z58W0R2FYfs6gLZzQq6zs63V2TbByJqr24w5FfRKIXaRNpkfnZ1N0TUz+oE0Dim99lU7sQCAJAcIRYAgORoJ4BGqm+zg8IrEvSIumSnz7PHZ9sMhp75VP6Y9c7vcUPdC/pnVyrIrmiw8OD8CgPVZQNy424v1L15QMXQfJvAik+vrvXc/xmYf019p+dXQPjkpfl2iVfPzbdCfJBpU1hyXGZTg0fq3qygPtlzsu0BBVctyGrAz2tjaCEgJabyaQ3ciQUAIDlCLAAAydFOAAUUmmqub9q40PcasglCrYX3M+cWWtkgovZUela2JSDbWtDzhfwxSw+s89SYf2j++G5v5FsR2lTU/uei5wv59oCKQfnfh0uX5zcN6Ds9f37X2dkVBnrmRtmWiuyGCNnXuv7rLLQZQbsCmx00aDWKjWBFAlKhfYDWyJ1YAACS404sACTMXVa2VEIsFNBc08PZqfHqT+Wn3j+cPiM3bp85pv2y/LmrM8d/9L1VUZfsBgdl5fml+z/okZ982WnSvILXXSc71f+Dc5+o9b0L3/h6bpxdASHbvrD9Q/niV+64Vb6+snwrQva5KwblV0n4oEd+I4ceNzxZ67mzmyVkFWqvaNvYFQwaSAsBQPFoJwAAIDlCLAAkbNCkKcUuAYpCOwFsZtnp7DYFPllfEw2T/WR+p1n5x7PT7Uv3yq8E0PGd/KYE2Wn8bPtBVrat4ewpX631vZ4VdVfZ75Y5ufHLl+U3VOh3T/6Y0gLnZtsP1mzdJTdef4OJWPJO1KXQSgXZNoNsC0D2fcqeC0DL504sAADJcScWAJqZFQSg6W1UiJ0/f37cc889MW/evFi9uvZ+6j/72c+apDBorQptcJBVaFp8g0n/zHR4dvr9/Ux7QK8pc3Pj7GYCbTKrIURmQ4WFB+dXDvjkE/lp+77Tt6r11AsPzo+7vJGf1Hntv3fIjXealG8tyLYvdJ1d92YF2deTXXmh0OoJEbVbHgq1BJQ0sv3AJggALV+jQ+zUqVPjqKOOiiFDhsQrr7wSu+66a5SXl0dNTU3ssccezVEjAADU0ugQe+6558bZZ58dF198cXTr1i3uvPPO6N27d5xwwglxxBFHNEeNANCiaReAza/RIfbll1+O//u///vo5Hbt4v3334+uXbvGj370ozj66KPj1FNPbfIioSXblKnjhhxfq/1gve9lP3WfbQ/oWmAzgE6T/5kbVx9c98zJ4Mn5FqElY4fWeUxERJc38m0HvWd8UPC4dbIbJxTalCC7IUK21g4F2i7qu1ahVQjaFTimKaf9tRAANL9Gr07QpUuXXB9snz59Yu7cfL/d0qVLm64yAAAooNF3YvfZZ594/PHHY6eddoojjzwyzjrrrHjxxRfjL3/5S+yzzz7NUSMAfCxT+rBlaXSI/dnPfhYrV66MiIiLL744Vq5cGbfddlvsuOOOViZgi9RUU8eFFt4vtIJBROFP3Wen5btljsm2EGSn96PAxgdlc6ty4+xqARG1N1EoJNuOUOha72dWRug6+93ceG09182uspA9p1BrR0PeTysKAKSl0SF2yJAhuXGXLl3iuuuua9KCAGBjNHb7VXduIW2N7okdMmRILFu2bIPH33333VoBFwAAmkujQ2x5eXmsXbvhRF9VVVUsWLCgSYoCAID6NLid4J577smNH3zwwSgrK8t9vXbt2pg6dWoMGjSoSYuD1qIh/ZaFdpvKqq9XM/scH/TI/37aKXNOdnerQstnZXtXszbsv833xK7fL5u7VoFas9fqVqCvN9sj/OF67023J+q+bqH3p76+4o87l5ZNSwBsuRocYseNGxcRESUlJfGNb3yj1vfat28fgwYNiquuuqpJi8v6yU9+ElOmTImZM2dGhw4d4t1332225wIAoGVrcIitrq6OiIjBgwfHM888E9tss02zFVWX1atXx3HHHRf77rtv/O53v9uszw0AQMvS6NUJXn/99dz4gw8+iI4dOzZpQYVcfPHFERFx0003bZbng6bUVNPc9R2fPa5X5kPaNZlp+ew0/jZP132dNVt3yY3f3iP/3/f2y2rvjLXN08tz4+xyWx3fqc6Ns20N2cerMktkZY8p1MpQX/N+dseuhvyDpm2gdWnsigR10ZIAaWr0B7uqq6vjxz/+cWy//fbRtWvXeO211yIi4vzzz3eHFACAzaLRIfaSSy6Jm266KX76059Ghw4dco/vuuuu8dvf/rZJi9tUVVVVUVlZWesPAADpa3Q7wc033xzXX399jB49Ok455ZTc48OGDYtXXnmlUdeaNGlSXH755fUe8/LLL8enPvWpxpYZERGXXnpprg0BWqLGTm039PhaxxU6p8Dj2d9st5ueH9e3g1aPWXU/3qmecxpzzPoa9Pq2MLMuHhPdu3cvdhkAm02jQ+yCBQtihx122ODx6urqWLNmTaOuddZZZ8X48ePrPWZTNlA499xzY+LEibmvKysro3///ht9PQAAWoZGh9idd945HnvssRg4cGCtx//85z/H7rvv3qhr9erVK3r16tXYEhqstLQ0SktLm+36AAAUR6ND7AUXXBDf+MY3YsGCBVFdXR1/+ctf4tVXX42bb7457rvvvuaoMSIi5s2bF8uXL4958+bF2rVrY+bMmRERscMOO0TXrl2b7XkBAGh5Gv3BrqOPPjruvffeePjhh6NLly5xwQUXxMsvvxz33ntvHH744c1RY0R8FJ533333uPDCC2PlypWx++67x+677x7PPvtssz0nAAAtU0lNTU1NsYvYXCorK6OsrCxGxdHRrqR9scsB2GQf1LwXj8f9UVFR4YNdQKuwLq993L9rjW4nWGf16tXx9ttv53byWmfAgAEbe0kAAGiQRofY2bNnxze/+c144oknaj1eU1MTJSUlsXZtfQvxAADApmt0iB0/fny0a9cu7rvvvujTp0+UlJQ0R10AAFBQo0PszJkz47nnntvoDQgAAGBTNXp1gp133jmWLl3aHLUAAECDNCjEVlZW5v5cfvnlcc4558S0adNi2bJltb5XWVnZ3PUCAEDD2gm22mqrWr2vNTU1MXr06FrH+GAXAACbS4NC7N///vfcuLy8PPr37x9t27atdUx1dXXMmzevaasDAIA6NHqzg7Zt28Zbb70VvXv3rvX4smXLonfv3i36TqzNDoDWxmYHQGvT0M0OGv3BrnVtA+tbuXJldOzYsbGXAwCARmvwElsTJ06MiIiSkpI4//zzo3PnzrnvrV27Nv75z3/G8OHDm7xAAABYX4ND7PPPPx8RH92JffHFF6NDhw6573Xo0CGGDRsWZ599dtNXCAAA62lwiF334a4TTzwx/vd//1fvFQAARdPoHbtuvPHG5qgDAAAarNEf7AIAgGITYgEASI4QCwBAcoRYAACSI8QCAJAcIRYAgOQIsQAAJEeIBQAgOUIsAADJEWIBAEiOEAsAQHKEWAAAkiPEAgCQHCEWAIDkCLEAACRHiAUAIDlCLAAAyRFiAQBIjhALAEByhFgAAJIjxAIAkBwhFgCA5AixAAAkR4gFACA5QiwAAMkRYgEASI4QCwBAcpIIseXl5XHSSSfF4MGDo1OnTjF06NC48MILY/Xq1cUuDQCAImhX7AIa4pVXXonq6ur4zW9+EzvssEPMmjUrJkyYEKtWrYorr7yy2OUBALCZJRFijzjiiDjiiCNyXw8ZMiReffXV+PWvfy3EAgBsgZIIsXWpqKiInj171ntMVVVVVFVV5b6urKxs7rIAANgMkuiJXd+cOXPi2muvjW9961v1HnfppZdGWVlZ7k///v03U4UAADSnoobYSZMmRUlJSb1/XnnllVrnLFiwII444og47rjjYsKECfVe/9xzz42KiorcnzfffLM5Xw4AAJtJUdsJzjrrrBg/fny9xwwZMiQ3XrhwYRxyyCGx3377xfXXX/+x1y8tLY3S0tJNLRMAgBamqCG2V69e0atXrwYdu2DBgjjkkENizz33jBtvvDHatEmyEwIAgCaQxAe7FixYEKNGjYqBAwfGlVdeGUuWLMl9b7vttitiZQAAFEMSIfahhx6KOXPmxJw5c6Jfv361vldTU1OkqgAAKJYk5uTHjx8fNTU1df4BAGDLk0SIBQCALCEWAIDkCLEAACRHiAUAIDlCLAAAyRFiAQBIjhALAEByhFgAAJIjxAIAkBwhFgCA5AixAAAkR4gFACA5QiwAAMkRYgEASI4QCwBAcoRYAACSI8QCAJAcIRYAgOQIsQAAJEeIBQAgOUIsAADJEWIBAEiOEAsAQHKEWAAAkiPEAgCQHCEWAIDkCLEAACRHiAUAIDlCLAAAyRFiAQBIjhALAEByhFgAAJIjxAIAkBwhFgCA5AixAAAkR4gFACA5QiwAAMkRYgEASI4QCwBAcoRYAACSk0yIPeqoo2LAgAHRsWPH6NOnT3zta1+LhQsXFrssAACKIJkQe8ghh8Ttt98er776atx5550xd+7cOPbYY4tdFgAARdCu2AU01He/+93ceODAgTFp0qQYN25crFmzJtq3b1/EygAA2NySCbFZy5cvjz/96U+x33771Rtgq6qqoqqqKvd1ZWXl5igPAIBmlkw7QUTE97///ejSpUtsvfXWMW/evLj77rvrPf7SSy+NsrKy3J/+/ftvpkoBAGhORQ2xkyZNipKSknr/vPLKK7njv/e978Xzzz8ff/vb36Jt27bx9a9/PWpqagpe/9xzz42KiorcnzfffHNzvCwAAJpZSU19KbCZLVmyJJYtW1bvMUOGDIkOHTps8Pj8+fOjf//+8cQTT8S+++7boOerrKyMsrKyGBVHR7sSfbRA+j6oeS8ej/ujoqIiunfvXuxyADbZurz2cf+uFbUntlevXtGrV6+NOre6ujoiolbPKwAAW4YkPtj1z3/+M5555pk44IADokePHjF37tw4//zzY+jQoQ2+CwsAQOuRxAe7OnfuHH/5y19i9OjR8clPfjJOOumk+PSnPx3Tp0+P0tLSYpcHAMBmlsSd2N122y0eeeSRYpcBAEALkcSdWAAAyBJiAQBIjhALAEByhFgAAJIjxAIAkBwhFgCA5AixAAAkR4gFACA5QiwAAMkRYgEASI4QCwBAcoRYAACSI8QCAJAcIRYAgOQIsQAAJEeIBQAgOUIsAADJEWIBAEiOEAsAQHKEWAAAkiPEAgCQHCEWAIDkCLEAACRHiAUAIDlCLAAAyRFiAQBIjhALAEByhFgAAJIjxAIAkBwhFgCA5AixAAAkR4gFACA5QiwAAMkRYgEASI4QCwBAcoRYAACSI8QCAJAcIRYAgOQIsQAAJEeIBQAgOUIsAADJSS7EVlVVxfDhw6OkpCRmzpxZ7HIAACiC5ELsOeecE3379i12GQAAFFFSIfaBBx6Iv/3tb3HllVcWuxQAAIqoXbELaKjFixfHhAkT4q677orOnTs36JyqqqqoqqrKfV1RURERER/GmoiaZikTYLNaE6sjIqKysrLIlQA0jXX/ntXU1B/WkgixNTU1MX78+DjllFNixIgRUV5e3qDzLr300rj44os3ePzxuL+JKwQorv79+xe7BIAmtWLFiigrKyv4/ZKaj4u5zWjSpElx+eWX13vMyy+/HH/729/i9ttvj+nTp0fbtm2jvLw8Bg8eHM8//3wMHz684Lnr34l99913Y+DAgTFv3rx63xTyKisro3///vHmm29G9+7di11OMrxvjec92zjvvPNODBo0KMrLy6NHjx7FLgdgk9XU1MSKFSuib9++0aZN4c7XoobYJUuWxLJly+o9ZsiQIfGlL30p7r333igpKck9vnbt2mjbtm2ccMIJ8fvf/75Bz1dZWRllZWVRUVHhf5IN5D3bON63xvOebRzvG7ClKmo7Qa9evaJXr14fe9w111wTl1xySe7rhQsXxpgxY+K2226LvffeuzlLBACgBUqiJ3bAgAG1vu7atWtERAwdOjT69etXjJIAACiipJbY2lSlpaVx4YUXRmlpabFLSYb3bON43xrPe7ZxvG/AlqqoPbEAALAxtqg7sQAAtA5CLAAAyRFiAQBIjhALAEBytugQO2XKlNh7772jU6dO0aNHjxg3blyxS0pGVVVVDB8+PEpKSmLmzJnFLqfFKi8vj5NOOikGDx4cnTp1iqFDh8aFF14Yq1evLnZpLc4vf/nLGDRoUHTs2DH23nvvePrpp4tdUot16aWXxsiRI6Nbt27Ru3fvGDduXLz66qvFLgtgs9piQ+ydd94ZX/va1+LEE0+Mf/3rX/GPf/wjvvKVrxS7rGScc8450bdv32KX0eK98sorUV1dHb/5zW/ipZdeiquvvjquu+66+MEPflDs0lqU2267LSZOnBgXXnhhzJgxI4YNGxZjxoyJt99+u9iltUjTp0+P0047LZ566ql46KGHYs2aNfGZz3wmVq1aVezSADabLXKJrQ8//DAGDRoUF198cZx00knFLic5DzzwQEycODHuvPPO2GWXXeL555+P4cOHF7usZFxxxRXx61//Ol577bVil9Ji7L333jFy5Mj4xS9+ERER1dXV0b9//zjjjDNi0qRJRa6u5VuyZEn07t07pk+fHgcddFCxywHYLLbIO7EzZsyIBQsWRJs2bWL33XePPn36xGc/+9mYNWtWsUtr8RYvXhwTJkyIP/zhD9G5c+dil5OkioqK6NmzZ7HLaDFWr14dzz33XBx22GG5x9q0aROHHXZYPPnkk0WsLB0VFRUREf5eAVuULTLErrsDdtFFF8V5550X9913X/To0SNGjRoVy5cvL3J1LVdNTU2MHz8+TjnllBgxYkSxy0nSnDlz4tprr41vfetbxS6lxVi6dGmsXbs2tt1221qPb7vttrFo0aIiVZWO6urqOPPMM2P//fePXXfdtdjlAGw2rSrETpo0KUpKSur9s65HMSLihz/8YRxzzDGx5557xo033hglJSVxxx13FPlVbH4Nfd+uvfbaWLFiRZx77rnFLrnoGvqeZS1YsCCOOOKIOO6442LChAlFqpzW5rTTTotZs2bFrbfeWuxSADardsUuoCmdddZZMX78+HqPGTJkSLz11lsREbHzzjvnHi8tLY0hQ4bEvHnzmrPEFqmh79sjjzwSTz755AZ7tI8YMSJOOOGE+P3vf9+MVbYsDX3P1lm4cGEccsghsd9++8X111/fzNWlZZtttom2bdvG4sWLaz2+ePHi2G677YpUVRpOP/30uO++++LRRx+Nfv36FbscgM2qVYXYXr16Ra9evT72uD333DNKS0vj1VdfjQMOOCAiItasWRPl5eUxcODA5i6zxWno+3bNNdfEJZdckvt64cKFMWbMmLjtttti7733bs4SW5yGvmcRH92BPeSQQ3J3/Nu0aVUTIJusQ4cOseeee8bUqVNzy9xVV1fH1KlT4/TTTy9ucS1UTU1NnHHGGTF58uSYNm1aDB48uNglAWx2rSrENlT37t3jlFNOiQsvvDD69+8fAwcOjCuuuCIiIo477rgiV9dyDRgwoNbXXbt2jYiIoUOHugtUwIIFC2LUqFExcODAuPLKK2PJkiW577nLmDdx4sT4xje+ESNGjIi99torfv7zn8eqVavixBNPLHZpLdJpp50Wt9xyS9x9993RrVu3XO9wWVlZdOrUqcjVAWweW2SIjfhomaN27drF1772tXj//fdj7733jkceeSR69OhR7NJoRR566KGYM2dOzJkzZ4OgvwWublfQ8ccfH0uWLIkLLrggFi1aFMOHD4+//vWvG3zYi4/8+te/joiIUaNG1Xr8xhtv/Ng2F4DWYotcJxYAgLRpzgMAIDlCLAAAyRFiAQBIjhALAEByhFgAAJIjxAIAkBwhFgCA5AixAAAkR4iFJvDWW2/FV77ylfjEJz4Rbdq0iTPPPLPO4+6444741Kc+FR07dozddtst7r///s1bKAC0EkIsNIGqqqro1atXnHfeeTFs2LA6j3niiSfiv/7rv+Kkk06K559/PsaNGxfjxo2LWbNmbeZqASB9tp2FBliyZEnstttu8e1vfzt+8IMfRMRHoXTUqFHxwAMPxOjRo3PHjho1KoYPHx4///nPa13j+OOPj1WrVsV9992Xe2yfffaJ4cOHx3XXXbdZXgcAtBbuxEID9OrVK2644Ya46KKL4tlnn40VK1bE1772tTj99NNrBdj6PPnkk3HYYYfVemzMmDHx5JNPNkfJANCqtSt2AZCKI488MiZMmBAnnHBCjBgxIrp06RKXXnppg89ftGhRbLvttrUe23bbbWPRokVNXSoAtHruxEIjXHnllfHhhx/GHXfcEX/605+itLS02CUBwBZJiIVGmDt3bixcuDCqq6ujvLy8Uedut912sXjx4lqPLV68OLbbbrsmrBAAtgxCLDTQ6tWr46tf/Wocf/zx8eMf/zhOPvnkePvttxt8/r777htTp06t9dhDDz0U++67b1OXCgCtnp5YaKAf/vCHUVFREddcc0107do17r///vjmN7+ZW21g5syZERGxcuXKWLJkScycOTM6dOgQO++8c0REfOc734mDDz44rrrqqhg7dmzceuut8eyzz8b1119frJcEAMmyxBY0wLRp0+Lwww+Pv//973HAAQdERER5eXkMGzYsLrvssjj11FOjpKRkg/MGDhxYq+3gjjvuiPPOOy/Ky8tjxx13jJ/+9Kdx5JFHbq6XAQCthhALAEBy9MQCAJAcIRYAgOQIsQAAJEeIBQAgOUIsAADJEWIBAEiOEAsAQHKEWAAAkiPEAgCQHCEWAIDkCLEAACRHiAUAIDn/H9uDeh4KGghBAAAAAElFTkSuQmCC", + "text/plain": [ + "Figure(PyObject
)" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig, axis = plt.subplots(2, 2, figsize=(8,8))\n", + "fig.suptitle(\"HMC - 21-D Neal's Funnel\", fontsize=16)\n", + "\n", + "fig.delaxes(axis[1,2])\n", + "fig.subplots_adjust(hspace=0)\n", + "fig.subplots_adjust(wspace=0)\n", + "\n", + "axis[1,1].hist(x10_hmc, bins=100, range=[-6,2])\n", + "axis[1,1].set_yticks([])\n", + "\n", + "axis[2,2].hist(theta_hmc, bins=100, orientation=\"horizontal\", range=[-4, 2])\n", + "axis[2,2].set_xticks([])\n", + "axis[2,2].set_yticks([])\n", + "\n", + "axis[2,1].hist2d(x10_hmc, theta_hmc, bins=100, range=[[-6,2],[-4, 2]])\n", + "axis[2,1].set_xlabel(\"x10\")\n", + "axis[2,1].set_ylabel(\"theta\");" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "2c9052ab", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAArEAAAL3CAYAAACd2x1cAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABOfElEQVR4nO3deXiU1f3//9eEQEggCaSGTQJJkOLagGyyVEFAEBGxSLVaWQRcPkhL0VpwYVEr7rVSflqVRREqIrKDiiggCq2yqMiiIBEEAoGQhDUBMt8//OW+z4SZMCHL5CTPx3Xl8uTOfd9zMiH44pz3fY7H6/V6BQAAAFgkLNQdAAAAAIqKEAsAAADrEGIBAABgHUIsAAAArEOIBQAAgHUIsQAAALAOIRYAAADWIcQCAADAOoRYAAAAWIcQC5xDYmKiPB6Ppk2bVuh5nTp1ksfj0bhx43yOr1ixQh6Px/nYsGFDofe57LLLnHOHDBkS8LxTp05p6tSp6tOnjxo1aqTIyEhFRUUpOTlZt9xyi2bMmKHc3FyfawYOHOjTF4/Ho8jISNWtW1etW7fW3XffrQULFuj06dOF9tGUkpIij8ejiIgIHTp0KOjrStKBAwf01ltv6fbbb1fTpk1VvXp1RUVF6eKLL9af/vQnpaamBrx21apVeuqpp9S3b1/nZ+3xeLR69epi9cm8l8fjUVhYmKKjo9WwYUN17txZDz74oP73v/8V6zUCGTdunPO6vXr1Cnje22+/LY/Ho06dOpVKP4KV/7uzYsUKn+P538e5fveKw9/vRMGPPn36lNrrlyepqanyeDxKTEwMdVeAoISHugNAZTNlyhRNnDjR79fWrl2rzZs3n/Me69ev1y233KKdO3fK4/EoJSVFbdq0UVhYmFJTUzVv3jzNmTNHjzzyiDZv3qyoqCif65s0aaKOHTtKkk6fPq3Dhw9r06ZNev311/X666+rcePGmjx5srp06VJoP7788kt98803kqTc3Fy9/fbb+vOf/xzM21CiRo4cqRkzZigsLEyXX365evfurWPHjunLL7/UxIkTNWXKFM2dO1fdunU769o//elP+vrrr0utbx06dNBFF10kSTpx4oQOHjyoDRs2aMWKFXrhhRd0zTXXaMqUKUpOTi6V11+8eLFWrVqlq6++ulTuX1GYvxMFXXnllWXcGwDBIMQCZaRRo0Y6efKkZs6cqeeff14RERFnnTNlyhRJUuvWrfXll1/6vc/69ev129/+VsePH1evXr308ssvKykpyeec9PR0/eMf/9ALL7yg3Nzcs0Jsx44d/Y5uff311xo1apQ++OADde/eXXPnztWNN94Y8HuaPHmyJOnCCy/Unj17NHny5JCE2Li4OI0fP16DBw/WhRde6Bw/evSohg4dqnfeeUe33Xabtm/frtq1a/tc261bN91888268sordeWVV6pDhw766aefSqxvQ4YM0cCBA32Oeb1eLV26VCNGjNDKlSvVvn17rVmz5qyfY3FFRUXp+PHj+tvf/qY1a9aU6L0rmkC/EwDKL8oJgDJStWpV/fGPf1RGRobmzZt31tePHz+ud955RxdeeKG6d+/u9x6nTp1Sv379dPz4cfXp00fz58/3G3zi4+P11FNPafXq1X7DciApKSlasmSJbr31Vp05c0YDBgxQdna233OPHz+u//znP5Kk6dOnq2bNmvr2228Dhu/S9PLLL2vMmDE+AVaSatasqcmTJys6OloZGRlavHjxWdc+99xzGjt2rG688cazri8tHo9HPXv21P/+9z81bdpU+/fvL7R05HzdfPPNSkhI0Nq1azV37twSvz8AhBIhFihDd911lyR3xNU0e/ZsHTlyRP3791eVKlX8Xj9z5kz9+OOPqlatml555RWFhRX+K9y6dWtFRkYWqY8ej0eTJk1SZGSkDh8+rNdff93vebNnz1Z2drYuv/xyde7cWbfeeqskd3S2vIiKilKzZs0kSbt37w5xb3zVqlVLL730kiTpk08+0bp160r0/tWrV9fjjz8uSXr44Yd15syZIl1/+PBhjR07Vs2bN1d0dLSioqJ0xRVX6Mknn9Tx48fPOv/IkSN6/fXX9bvf/U5NmzZVjRo1VKNGDV1xxRV65JFHlJmZWRLflvLy8vTaa6+pQ4cOqlWrlqpWrao6deooJSVFw4cPL7QGujimTZsmj8dz1sh6vkA1peZxr9er1157TS1btlSNGjUUGxur6667LuBIeX5driTNmTNHHTt2VExMjGrUqKEOHTpoyZIlAft7+vRpvfHGG+rUqZPi4uIUERGhpKQk3XfffeXudwE4H4RYoAxddtllatOmjT7++OOz/ieSH/4GDRoU8Pr58+dLkrp376569eqVWj9/9atfqUePHpKkZcuW+T0nv7/5wTz/v++8845OnDhRan0rqlOnTjmhpn79+qHtjB/XX3+94uLiJAV+r4ujf//+uvzyy7V161a//3gKZPPmzUpJSdHjjz+uAwcOqGPHjuratavS09P12GOPqUOHDsrKyvK55uuvv9bdd9+t1atXq169errxxhvVsWNH7du3T0899ZRat25dIg//DRkyRPfcc4/Wr1+v1q1bq1+/frryyit14sQJ/etf/9LGjRuL/RqlZdCgQbr//vtVq1Yt9erVS/Xq1dOyZcvUuXNn/fe//w143dixY9WvXz9JUs+ePdW0aVN98cUX6tWrl99R9iNHjqhbt24aOnSo1q1bp9/85jfq3bu3IiIi9Oqrr6pFixbnfMgUKO8IsUAZGzx4sPLy8nzq73744Qd99tlnuvrqq9W0adOA13711VeSfhlhLW0tW7aUJH333Xdnfe3777/XZ5995pRISFL79u118cUXKysrS++9916p9y9YkydP1sGDBxUZGanrr78+1N05i8fjcR4c8vdeF1dYWJieeuopSdL48eOD+gfGiRMn1Lt3b+3evVuPPvqoUlNTtXTpUi1YsEA//vij/vCHP2jjxo36y1/+4nNdYmKiPv74Y6Wlpemzzz7TO++8ow8//FC7du1S//79tX37do0ZMybovo8bN05er9dn5HPXrl2aOnWqGjZsqJ07d2rZsmWaOXOmPvjgA33//ffavHlzuX0Q66efftKKFSu0adMmLV++XLNmzdLmzZt11113KScnp9D35uWXX9aaNWu0evVqvfPOO9q4caPGjh0rr9erUaNGnXX+vffeqxUrVqhXr17asWOHVqxYodmzZ2vr1q36xz/+oUOHDjllQ4CtCLFAkAYNGlToMjwrV64M6j633XaboqKiNG3aNHm9XklueUH+aGYg6enpkqQ6deoU4zsJzgUXXCBJfkfO8vvbu3dvxcfHO8fz+19eSgq+/fZb/fWvf5UkPfbYY6pbt26Ie+RfYe91Sbjxxhv129/+Vnv27NE///nPc57/5ptvaseOHerVq5eeeOIJVatWzflaVFSUXnvtNdWpU0fTp0/X4cOHna81bNhQXbp0OavMJSoqSq+88orCw8M1e/bsYn0v+/fvl/TLigH+ZiMuueQSNWrUqMj3ffPNNwP+bpekiRMn6te//rXzeZUqVfT3v/9dkrRy5UqdOnXK73WPP/642rZt63Ns9OjRio2N1ffff+8zs7Nlyxb95z//UYMGDTRz5syz/r4YMWKEevbsqR9++EFLly4tqW8NKHOsTgAEyVwqyZ8PPvjA+R9sYWJiYtS3b19Nnz5dK1as0NVXX6233npL0dHRznRheZCXlydJZ/1P/PTp03rzzTclnR26+/fvr4cfflirVq3Sjh071KRJk7LprB8///yzbrzxRh09elS9e/f2O1pVXgR6r0vSM888o/bt2+uZZ57R3Xff7ZQw+JP/AFx+nXNBNWvWVKtWrbRkyRJ9+eWXuu6663y+/sUXX+izzz7Trl27dPz4cecfa9WqVVN6eroOHz581ioRwbr44osVHR2tJUuW6O9//7tuv/32ElnVobAltkpKeHi4U6ZjqlevnmrXrq3Dhw/r0KFDfsO5v1VCIiIilJycrA0bNmjPnj1KSEiQJC1ZskRer1fXX3+9oqOj/falU6dOWrJkiVOSANiIEAsEyd9SSaZOnToFFWKlX8Lf9OnTNWXKFB0/flx79+7VkCFDzloKq6D4+Hjt3r1bBw4cKErXz8vBgwcl6ayws3jxYqWlpfldRaFu3brq2bOnFixYoClTpjgjTOcyb948vys2DBky5LyCRVpamrp06aKffvpJ3bt317vvvluiAdHfn4MLLrhAzz///HndL9B7XZLatWunPn36aN68eXrqqacK7euPP/4oSbrzzjt15513Fnrf/NkB6ZdNJ/r27XvOjSKys7PPO8RGR0dr6tSpGjRokB599FE9+uijql+/vq666ir16NFDt99+u2rWrFnk+5bFElv169dX1apV/X4tJiZGhw8f1smTJ/1+PdDockxMjCT5XJf/85s8efI5Z0XMnx9gG0IsEALXXHONmjRpojlz5igtLU3SuUsJpF/qVHfv3l0my1itX79eknTFFVf4HM//n+LJkyd1zTXXnHXdnj17JP3yJPfjjz8ecKUF08aNG53RXVOnTp2KHGIPHDiga6+9Vt9//726du2qefPmFWmZsWD462vjxo3PK8R6vV7nAZuC73VJe+qpp7Rw4UJNmjSp0PV880eGe/Tocc4SjMaNGzvtIUOGaPXq1WrXrp3Gjx+vlJQU1a5d2wluDRo00L59+5yR2fPVt29fde3aVQsWLNBnn32mzz//XHPnztXcuXM1ZswYLVu2rNTfS3/y37dAzrWaSEldm9+P5s2bKyUlpdBzC5YoADYhxAIhkL9Mz2OPPaaPP/5Yl1xyidq1a3fO62666SbNmzdPH374ofbv319qNZ4HDx7Uhx9+KEk+U8X79u1zlvQ5dOiQPv/884D32Lt3rz744APdcMMN53y9cePGnbVd7/lIT0/Xtddeqy1btqhLly5asGCBqlevXuz7FlTcEGZasmSJU1dacFq+pF1yySUaOHCgJk+erDFjxgTckS0hIUFbt27V4MGDdcsttwR172PHjmnJkiUKCwvTkiVLVKtWrbO+nv8PtpIQGxvrM1K8e/duDR8+XPPnz9f9998fdI16UeTXBh85csTv10tyk4ziyC8r6NChg/71r3+FuDdA6eHBLiBEBg4cqPj4eP3qV7/SPffcE9Q1d9xxhxITE5Wbm6v77rvvnCM/69atK/JyV16vV/fff79OnDihuLg4DR482PnatGnTdObMGbVt21Zerzfgx0MPPSSpbB/wOnjwoK699lp999136tKlixYuXFjkNXLLWlZWlvOEf7du3dS8efNSf83x48crMjJSb731VsDVEPJXcXj33XeDvm9WVpbOnDmjmJiYswKsJL399tslGv4LSkhI0Pjx4yWp1JbYyt8MY+vWrX6/7m8zjVDI//ktWLAgYHkCUBEQYoEQadiwoQ4cOKCDBw8GvVVr1apV9e6776p69eqaO3eu+vTpo507d551XkZGhrOWZ05OTtB9+uabb9SzZ0/NmjVLVapU0dtvv+3zYEj+qgQDBgwo9D79+/eXJC1atKhMau4yMjLUpUsXbdq0SV27di33ATZ/29k2bdrohx9+UP369QNuKlHSLrzwQg0fPlx5eXl6+eWX/Z5z9913q3Hjxpo9e7b+9re/+R15TEtL8+lz3bp1Vbt2bWVmZmr69Ok+565du1ajR48ukf5v2LBBs2bN8vuPs4ULF0ryLXEoSW3atFFMTIw2b9581vc4e/bsgO9nWWvRooX69u2r3bt363e/+53fzR+OHTumGTNmBF3HD5RHlBMAlmndurVWrVqlfv36aeHChVq0aJFatGih5ORkhYWF6aefftJXX32lM2fOKDk52W896OrVq52Hk06fPq3MzExt2rTJmQ5NSkrS5MmT1blzZ+ealStXavv27YqIiNBtt91WaB8vu+wyXXnllVq/fr3eeustPfDAAyX3BvgxZMgQffPNN/J4PIqLi9N9993n97w+ffqoT58+PsfeeOMNvfHGG87n+/btkyTdc889ToCvX7/+eW/b+sYbb2jFihWSpJycHB08eFDr169XRkaGpF/qfqdMmVJqwcuf0aNH6/XXX/dZHstUo0YNLV68WL169dKzzz6r1157Tb/5zW/UsGFDHT9+XN9//722bNmiOnXqaOjQoZJ+WSpqzJgx+stf/qL+/ftr0qRJSk5O1q5du/TFF1/oj3/8o1atWlXsKfeffvpJt912myIjI3XllVcqISFBp0+f1rfffqtt27apWrVqevbZZ4v1GoFERkZq/Pjxzvf4yiuv6MILL9SWLVu0efNmPfroo3riiSdK5bWLaurUqcrMzNTSpUvVrFkzpaSkKCkpSV6vV6mpqfr666+Vm5urLVu2lNul54BzIcQCFmrdurW+//57TZ8+XfPnz9f69ev13XffyePxqH79+rr55pvVt29f9e3b1+/T0Dt27NCOHTsk/bJMT2xsrBISEnTdddepV69e6tmzp8LDff96yC8NuPHGG4N6srx///5av369Jk+eXOohNj8Qer3eQqfAExMTzwqxP//8s9+dkjZv3uy0ixMwP//8c6d2OH+b0SuuuEKtWrXSrbfeWiYbVxRUq1YtjR492in78Oeyyy7TN998o1dffVVz587VN998ozVr1uiCCy5Qw4YN9eCDD+rmm2/2uWbEiBFKSkrSs88+q82bN+u7777TxRdfrEmTJunee+8tkaWwrrrqKj399NNatWqVtmzZog0bNig8PFwNGzbUsGHDNHz4cGeb4dIwYsQIxcXF6Z///Kc2bNig7777Tq1atdJLL72kiy66qNyE2OjoaH300UeaNWuW3n77ba1bt04bN25UTEyM6tevrzvuuEO9e/cO6TJ4QHF5vKVZpAQAAACUAmpiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwDiEWAAAA1iHEAgAAwDqEWAAAAFiHEAsAAADrEGIBAABgHUIsAAAArEOIBQAAgHUIsQAAALAOIRYAAADWIcQCAADAOoRYAAAAWIcQCwAAAOsQYgEAAGAdQiwAAACsQ4gFAACAdQixAAAAsA4hFgAAANYhxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwDiEWAAAA1iHEAgAAwDqEWAAAAFiHEAsAAADrEGIBAABgHUIsAAAArEOIBQAAgHUIsQAAALAOIRYAAADWIcQCAADAOoRYAAAAWIcQCwAAAOsQYgEAAGAdQiwAAACsQ4gFAACAdQixAAAAsA4hFgAAANYhxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwDiEWAAAA1iHEAgAAwDqEWAAAAFiHEAsAAADrEGIBAABgHUIsAAAArEOIBQAAgHUIsQAAALAOIRYAAADWIcQCAADAOuGh7kBZysvL0969exUdHS2PxxPq7gBAsXm9Xh05ckQNGjRQWBjjEgAqj0oVYvfu3auEhIRQdwMAStzu3bvVsGHDUHcDAMpMpQqx0dHRkn75yz4mJibEvQGA4svOzlZCQoLz9xsAVBaVKsTmlxDExMQQYgFUKJRIAahsKKACAACAdQixAAAAsA4hFgAAANYhxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwDiEWAAAA1gkPdQcAlL3EUYuddurTN4SwJwAAnB9GYgEAAGAdQixQwSWOWuwz8goAQEVAiAUAAIB1CLEAAACwDiEWAAAA1iHEAgAAwDqEWAAAAFiHdWIBOFg/FgBgC0ZiAQAAYB1CLAAAAKxDiAUAAIB1qIkFKjl28wIA2IiRWAAAAFiHEAsAAADrEGIBAABgHUIsAAAArEOIBQAAgHUIsQAAALAOIRYAAADWIcQCAADAOoRYAAAAWIcQCwAAAOsQYgEAAGAdQiwAAACsQ4gFAACAdQixAAAAsA4hFgAAANYhxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwDiEWAAAA1iHEAgAAwDqEWAAAAFgnPNQdAFA2EkctDnUXAAAoMYzEAjinxFGLCcEAgHKFEAsAAADrEGIBAABgHUIsAAAArEOIBQAAgHUIsQAAALAOIRYAAADWYZ1YoIIwl8BKffqGEPYEAIDSx0gsAAAArEOIBQAAgHUIsQAAALAOIRYAAADWIcQCAADAOoRYAAAAWIcQCwAAAOsQYgEAAGAdQiwAAACsw45dQAVk7t4VynsAAFBaGIkFAACAdQixAAAAsA4hFgAAANYhxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwDiEWAAAA1iHEAgAAwDqEWAAAAFiHEAsAAADrEGIBAABgnfBQdwBA0SSOWuy0U5++IYQ9AQAgdBiJBQAAgHUYiQUsZo7KAgBQmTASCwAAAOsQYgEAAGAdQiwAAACsQ4gFAACAdQixAAAAsA4hFgAAANYhxAIAAMA6hFgAAABYhxALAAAA67BjF4CgmTuEpT59Qwh7AgCo7BiJBQAAgHUIsQAAALAOIRYAAADWIcQCAADAOoRYAAAAWIcQCwAAAOuwxBaA88JyWwCAUGIkFgAAANYhxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1WCcWQLGxZiwAoKwxEgsAAADrEGIBAABgHUIsAAAArEOIBQAAgHUIsQAAALAOIRYAAADWIcQCAADAOoRYAAAAWIcQCwAAAOsQYgEAAGAdQiwAAACsEx7qDgAILHHU4lB3AQCAcokQC5QzBFcAAM6NcgIAAABYhxALAAAA6xBiAQAAYB1CLIBSkzhqMTW+AIBSQYgFAACAdQixAAAAsA4hFgAAANYhxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1wkPdAQAVCzt0AQDKAiOxAAAAsA4hFgAAANYhxAIAAMA6hFgAAABYhwe7AJQ682Gv1KdvCGFPAAAVBSOxAAAAsA4hFgAAANYhxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwDiEWAAAA1iHEAgAAwDqEWAAAAFiHEAsAAADrEGIBAABgHUIsAAAArEOIBQAAgHUIsQAAALAOIRYAAADWIcQCAADAOoRYAAAAWIcQCwAAAOsQYgEAAGAdQiwAAACsQ4gFAACAdQixAAAAsA4hFgAAANYhxAIAAMA6hFgAAABYhxALAAAA64SHugMApMRRi0PdBQAArMJILAAAAKxDiAUAAIB1KCcAQoQSAgAAzh8jsQAAALAOIRYAAADWIcQCAADAOoRYAAAAWIcQCwAAAOsQYgEAAGAdQiwAAACsQ4gFAACAdQixAAAAsA4hFgAAANYhxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwDiEWAAAA1gkPdQcAVC6JoxY77dSnbwhhTwAANmMkFgAAANYhxAIAAMA6hFgAAABYhxALAAAA6/BgF4BygQe+AABFQYgFEDJmcAUAoCgoJwAAAIB1CLEAAACwDiEWAAAA1qEmFihD1IACAFAyGIkFAACAdQixAAAAsA4hFgAAANYhxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwDtvOAmWA7WYBAChZjMQCAADAOoRYAAAAWIcQCwAAAOtQEwug3DFriFOfviGEPQEAlFeMxAIAAMA6jMQCJYhVCAAAKBuMxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwDiEWAAAA1iHEAgAAwDqEWAAAAFiHEAsAAADrEGIBlGuJoxYrcdTiUHcDAFDOEGIBAABgnfBQdwCwHaOEAACUPUZiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwDiEWAAAA1iHEAgAAwDqEWAAAAFiHEAsAAADrEGIBAABgHUIsAAAArMO2s8B5YKtZAABCi5FYAAAAWIcQCwAAAOsQYgEAAGAdQiwAAACsQ4gFAACAdQixAAAAsA4hFgAAANYhxAIAAMA6hFgAAABYhxALAAAA67DtLAArBNrqN/XpG8q4JwCA8oCRWAAAAFiHEAsAAADrUE4AnIM5jc3UNQAA5QMjsQAAALAOIRYAAADWIcQCAADAOoRYAAAAWIcQCwAAAOuwOgEAq7F6BABUTozEAgAAwDqEWAAAAFiHEAsAAADrEGIBAABgHUIsAAAArEOIBQAAgHUIsQAAALAO68QCRWCuSQoAAEKHkVgAFVLiqMX8owMAKjBGYgE/CD8AAJRvhFgAFQb/+ACAyoNyAgAAAFiHEAsAAADrEGIBAABgHUIsAAAArMODXcD/j4eCAACwByOxAAAAsA4hFgAAANahnACVGiUEAADYiRALoNIw/9GS+vQNIewJAKC4CLEAKjRG2wGgYqImFgAAANZhJBaVDiNzAADYjxCLCo0aSAAAKiZCLCocRloBAKj4CLGoNAi3AABUHDzYBQAAAOswEgugUqJeGgDsxkgsAAAArMNILIBKL9CobP5xRmoBoPwhxKJC4KEtAAAqF8oJAAAAYB1CLAAAAKxjTTnBhAkT9P7772vr1q2KjIxU+/bt9cwzz6hZs2ah7hqACsRfaUqgchVqZQEgdKwZiV25cqWGDRumtWvXatmyZTp16pSuu+46HTt2LNRdQxlIHLXY+QAAALBmJPaDDz7w+XzatGmqU6eO1q1bp6uvvjpEvUIoEGQBAIA1IbagrKwsSVJcXFyIe4J8JbF4PAEVNjnXn1fKDQCg9FgZYvPy8jRixAh16NBBl19+ecDzcnJylJOT43x++PBhSdLPP/+smJiYUu9nZXM6+6DT/vnnn4t9D8B25/t7UBSZmZmSpIyMjFJ/LQAoC16vV0eOHFGDBg0UFha48tXj9Xq9ZdivEnHfffdp6dKlWr16tRo2bBjwvHHjxmn8+PFl2DMAAACUhN27dxea86wLsffff7/mz5+vVatWKSkpqdBz/Y3EJiYmqqN6KlxVS7urAFDqTuiY/quPtXv3bmaYAFQI2dnZSkhIUGZmpmJjYwOeZ005gdfr1fDhwzV37lytWLHinAFWkiIiIhQREXHW8XBVVbiHEAvAflW91SRJMTExhFgAFYrH4yn069aE2GHDhmnmzJmaP3++oqOjlZaWJkmKjY1VZGRkiHsHAACAsmTNOrGvvPKKsrKy1KlTJ9WvX9/5mDVrVqi7BgAAgDJmzUisZaW7AAAAKEXWjMQCAAAA+QixAAAAsA4hFgAAANYhxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwDiEWAAAA1iHEAgAAwDqEWAAAAFiHEAsAAADrEGIBAABgHUIsAAAArEOIBQAAgHUIsQAAALAOIRYAAADWIcQCAADAOoRYAAAAWIcQCwAAAOsQYgEAAGAdQiwAAACsQ4gFAACAdQixAAAAsA4hFgAAANYhxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwDiEWAAAA1iHEAgAAwDqEWAAAAFiHEAsAAADrEGIBAABgHUIsAAAArEOIBQAAgHWsCrGrVq3SjTfeqAYNGsjj8WjevHmh7hIAAABCwKoQe+zYMaWkpGjSpEmh7goAAABCKDzUHSiK66+/Xtdff32ou4EKILxBfad9eu++Il0b1uJSn8/zNmwu0n2Lc07B1w7bf7hI9zIV9fxgrzeZ/TXfp2AE+zMqzs+yLJT3/gGArawKsUWVk5OjnJwc5/Ps7OwQ9gYAAAAlxapygqKaMGGCYmNjnY+EhIRQdwkAAAAloEKPxI4ePVojR450Ps/OzibIVgLBTIfn1a3tfmJM8Z7P1O+p61o57bD04047u3Oi045MP+X3HAUoFTgRH+W0dz3U2Gk3e/Ynn9c+cfmFTjttSJLTrrc21+1HfFX3/N9nOm3Px8b5q93jezvW8ns889JYn9eutdl9D81SAfP7yDG+DxnvU+SmPU77dID3P1iBfpblBSUEAFA6KnSIjYiIUERERKi7AQAAgBJWocsJAAAAUDFZNRJ79OhRbd++3fl8586d2rhxo+Li4tSoUaMQ9gwl6Xym9ItaQmA+1R8W4Npg+3G4WTWn3cCYJpcx/V71o6+cdo4xrR5RYLUBf/esmeT2dadRMlBQhHuadvV0f7Xjv3SP1x/ncfsR75YcmKUCZgmBTzlAAT/dWMu9Jt79nszv1SwhMI9n3NHOacd96p4e7M/bLFkwf5Z5QV3tn/nzDvRnRQquj6xIAAClz6oQ+9VXX6lz587O5/n1rgMGDNC0adNC1CsAAACUNatCbKdOneT1ekPdDQAAAISYVSEWyFewfKCw6V/nHHOx/QDT0ebT/uYT9IUt2l9P7tf29vX/xH+mz/R56jnPrzPxC6e9v+slZ30vzjXGKgRmCUFeLXc1hJNx7sONZtmAuWJCVhOP0a7ltBsvdPsUke772nGfuu9bhrESw+Hh7f32tbbc0oJam7Pcvpo/u7r+VzwoLYHKUHxWWziPFRMoIQCA0keIBYAK4PKxHyosInAN8/lIffqGEr0fAJQkVicAAACAdRiJRZGUxVPXwSx+f9Zrm5sGBNFHc7rYfKI90jwnQImCubpAQdUz/Ndsm1P35tR7gzk7/b5eygb335efPV/LaV/4qXu+5LtaQdM3M532vnFuP3Jqu+UEtbcZ/Rh2zGlHfOyuhhCsbcYmDE1GrHHaccb7b5ZnRBibPARTwhHsSgPm+2b+hRbMn8/z+rNWDLatWpA4anGou1ChMdINFA8jsQAAALAOIRYAAADW8Xgr0ZpV2dnZio2NVSfdpHBP1XNfAGsUdbMDs5zglLn5gDHlbS70b05zmysKSL4rCZj9MMsGYmaskT8HAjzJf6RJcJPpjZacdtrm6gTRO9x/n9be5q5gkHaVWzYQu8P91TdXJzCPp7d22/FfuudIvisM/DAgxm+fTsS7v2dmSYV53HwN89pA5QdS4BIEU6CfdzBT+jZN+5/0HtdqLVHCiHdL/MEu2I1yBdgqP69lZWUpJiYm4HmMxAIAAMA6hFgAAABYh9UJEDLFeRq8sPKBQNebC+mb14cHmI6uai62b5xvlg8UvJc5ze3z2sZmCeaGAzlulxThf48Gn9KA82GWEJivceL3mU775OEo4xx3NYOaSe45+2v5TlXv7+au5VD9R7ePZhmAuUGC+X2bpQjprWOM891rA60gIUknjFKNqkX8eQejuCUENpUjoOIq7uoSlCOgvGMkFgAAANYhxAIAAMA6lBMAAFDJUCqAioAQi5ApTq1gwWsD1Tz61CYWUmPp73yzjtU8v+AvzZanGjjtusvcpaPMuk9zuS7T6DtmO+0JM/o57WoZ7iSJWTebG+fb8yNNzv0rbF5zMtld5uqWRtuc9rKPr3La3e5a67SviNrttF/8l9u/gv1KesPdSczcyctcMivQ+9FkhPt6Mt7zwupKzeW6TOY15jsVaPe1QIpb00odLACUPsoJAAAAYB1GYgEAqGSKu3JBSaO8AeeDEIty53ymcs3p4vAAxxVoKSZjCtssOThg7Mxl7np1uGMtn+vrLnN3nIr7NNXva0QY7dTb3aWpxi/t67Qvudadkv8u1S1RSJzp7pSVMeyYz32np0x12htzEpz2//f4LU47q6U79R67zu3JlsS6Ttvb1Z1i35LlHl82xS0zKMzOIe571WSEuwRZ9h3tnPaJeHeJrZNx7vdUW+6OaebuWz4lAObPUb5LcSnAzy/gz9tcXq0Yy7wVFOjPLcttAUDpIMQCAICQKsrIMKO2yEdNLAAAAKzDSCzKnfOZcg30xLl53Hxa3SwhCHROgznu9P5eo7TAPC757tKVYewkdcHQVKe985NaTjt2nXtt9Qy3FGFHhvsaYcaKAifi3an3gvq/MNLvvczX/r/6XzntdxPdqXuzZCGmtjs9f0nsfqe9pYm70oC5YkJBJ5NznPbuMe2dtrlDWFZL95xLHt7rtM33LKJA2UA+c7cvKXDZhllaYO745bNqQRFXKghWoD+3lBAAQOkgxAIAAGuU5kNplCrYhXICAAAAWIeRWFQIgaZsfcoGjGlkc7H9SGNKOdC0c52J7hP3GcYT9wX99sH/Ou1v77rYaTdWpt/XNmU1qeb3+Infu9dmH/a99rK+bmmDWQbw/so2TnusUTYQiLnKwa1T3BKFaGO2/a0HXvS55qYVw5x23WVu30/GueeYqx7UfbeW39c2Ny4INL1f8EqzhMNc0cDneuNnn2GUI5ibLpj3qfqRW3YRaNUCKbiVB1iRAABKHyEWAABApVOqQIlC6aGcAAAAANZhJBblWjDTtQW/FhbgCXXzeKBF9c3zT5lT1sbrmdPRku8U9lfpVzrtCLlP/P8wIMZpR+9w/+2Y4/9hfI29fo7Tfnefu6LAjnUFLkh0m+YmBZc0/8lpm6sQjG+3wO99H0692WnnGisjnEx2p/pvnjfC56VvueZ/TnvRj+6mCLE73FUSIifVcNppV7mrLKS3dlc9iP/SPb7rIfd40zeznfZPN9byee2kN9wyioDT9cbxuP3+/xxVNU4P9s+aqTgrElByAADFQ4gFAAAoJedTokAJQnAoJwAAAIB1GIlFuVbcjQ/yAhz3ua/RPnWdO8VulhyYCq4ucMJ48t180j6Qeqsz/d4r9XZ3Gn7CjH5+r40o8PC+WSpwWaK7gcDurFpO+5aUDU7bLCEwmSsbPNVnrtO+8+tBTvtopu+GA4H6ldXEY7T9r7jQaMlpp13jYfd9rvUXt1zCLCGotzY34GubAq0uccJ4nyOC2Owg2Kn+QKUGwfy5pYQAAIqHEAsAAFCOFKUEoTKXHlBOAAAAAOswEgtrBDvFG8wT40V9Mtycmo5IP+7ztRPxsX6/Zq50EN2xvdPONMoPzKn3xJnulPnhZu79a29zj2cMO+bbMWPzgy0b3Sf7aya5KygsmuuuHDDrLnfDAnOzgu/kliW8pxZO22c1g1jfUgRzQwU1cQs3Aq2MEFPbfW92xbnvQbVPkpx20n531YHcOHdFhxPx5joCviUB5ioSO3q6f6WZJQvmRgZmiYm5IUKg0oLCBLOKAasQAChNlXnUlpFYAAAAWIcQCwAAAOtQTgBrBPuUeEk9GW6WEJgKrk4QyI6X3Gn8RkvckoC0q9wn9s2NAcwp85H3z3ba3x5PcNrvfe1O9Uu+KxI8leiuKnDrlJF++2RuamCuWmDe1zx+Z/RBpz1hhjvtL0k1W/qfft9hlAdUN45H7jBWKvh9ptPMzTBWETBKA6pluP/GPhnn96XOYm6QYJZtRBhlA6bzKSEI9GetOOUqAFAWEkctrlAlBYzEAgAAwDqEWAAAAFiHcgJYqeCT4MWZsjXvZW5VYC62n/SG+9S8CpQTmBscmFPYjZa4x82NE2LjE532yTh3dYLqGW5pwfilfZ22OUWuP/v/HiTpoZvvctpN/uH29/f13Sfzx67p7bS3ZLqrGczv85LTNksOVNdtnkzO8Xm9za3fcdoXvXuP+4U4/ysVHFvrlgro3VpO02OUCpglFeZ7vnOIbymD+T7X2uyuxGCWepjHTeaqEWEBVhQojE+ZSRArYQAASgcjsQAAALAOIRYAAADWsa6cYNKkSXruueeUlpamlJQUTZw4UW3atDn3hahQirtSgXmOOT1sLoqftMn/fSIL6cuJO9o5bXPjg4zOiU7bnOY2SxaqZ7j3NJ/Mv2LKVqf9w0rfP+s/z3Gn2Y8McKfxa2a5pQnvyt2kINCKBKbdWW6fWky4z2lf1nenz3mXvup+Ld5YZWF/N3clhmNPuSUEgVZl6HbXWqdtbsxwMs793uqtde9ZGLNsw2T+jMOMlQpOGOUH5s/edFbpilGOEOg8ygkAlFdF2RyhuEp7JQSrRmJnzZqlkSNHauzYsVq/fr1SUlLUvXt3HThwINRdAwAAQBmyKsS++OKLGjp0qAYNGqRLL71Ur776qqKiojRlypRQdw0AAABlyOP1er3nPi30cnNzFRUVpffee099+vRxjg8YMECZmZmaP3/+Wdfk5OQoJ8d9ojo7O1sJCQnqpJsU7ql61vmA6dR17jS8WRpQcIF8c4H+QNPZJnNq23zK3pyGj10X4bTNVQsKymri8Xu8ybU7/R43XRK732mbpQXj2y1w2lOH3+S0zXKAwl7P3Owg11ipoNmz7koFBVcb8McsOSi40oD5Mwi0KYUpJ0DZQEUoATjpPa7VWqKEEe8qLCK4jTgAoDxvepCdna3Y2FhlZWUpJiYm4HnWjMQePHhQZ86cUd26dX2O161bV2lpaX6vmTBhgmJjY52PhIQEv+cBAADALtY92FUUo0eP1siR7vab+SOxAAAAtijPo6ahZE2IveCCC1SlShXt37/f5/j+/ftVr149v9dEREQoIiLC79dgt4JPjJtKalrYnHYubFF8s4Qg0NS2OZ1tnh9nTItHprtlCSfig6vyMTcguCxxr9M2VwXIGHbMaU9Pmeq07/x6kNM2Swhe2NrV7VN8cGU3Wza6GydEG9UW9daeds95qoHTbvbPTKedE+9/Ctws4UjrWMvna9Uz3DIMc7OJQOUckcb7fNrvGcGrCCUIAFARWFNOUK1aNbVs2VLLly93juXl5Wn58uVq165dIVcCAACgorFmJFaSRo4cqQEDBqhVq1Zq06aNXnrpJR07dkyDBg0698UAAADFwLR++WJViL311luVnp6uMWPGKC0tTc2bN9cHH3xw1sNeAAAAqNisWWKrJOQv2cASWxVbUWsWi3N+wWvM3aACMWtAzWWrzF2pzONJb7jLV+3t67s0VVZLtya27jL3GnO5ruo/unXhZg2tKdCSXhcMTXXa5k5ekpR92P0+zNcwmUtsBRL/pbtMmLlkmPl+HG7mu7xXnYlf+L1XoJ3Y8oxdtgL9vAPVWRd3d7jSxhJbQNlhJLZsBLvEllUjsQAAAMEgcFZ81jzYBQAAAORjJBYVTqBp3UDTxT7LYgWYXg52StmcPDeXhaq3OtPvteaUubmkV0Sz9k7b3BGs4O5d3truMlQnfu+2EyfVcNqHm7nnNzGW4TJ37Fq27iqn/X9j3nPa3x5311X+eY5vKUNMV3fZqqNxbnlO9A7338YRh912ToCNtWptznTaJ+Nq+T2nwRzfHchOGLup+TB34zKPG2Uep4MoLShMeSghAAAwEgsAACqgxFGLQ90FlDJCLAAAAKxDOQEqjYBTvwGOBzo/o3Oiz+cxM9Y4bfMXqvY290lxc0UCcycqk7krWO1tbpnBrp7uXaN3eHyuObrT3bnKfMo/9Xb3+rBMt8jB3Flri9x2TaM0wCwheO/rFk77sr6+U/rmvWomZblf2OHWDXS7a63TXjTXLVlocq17r4M7Ep22+X2bCu6EdsLYSazWZve1A+2sFmbs2BUWxOoEhZUJUEIAAOUDIRYAAFRIgUoKWLmgYqCcAAAAANZhswMggPNZnSDQeeYmCObUtjlNbq5mYD6Nb5YvRKaf8rmvOa1ubhRgMjc4MDdEMJmbI5gbH5grCkQclg9zo4XEme5rm+UPJnPVAlOglRvMEoyCzO/bFPdpqt/j5vsc6P0PZkOE8ojNDoCSxSht6AW72QEjsQAAALAOIRYAAADW4cEuIIBgp5HN88yyAXN62md1AuNac2q7wRy3bd4zZob/+0vSiXh3dYLYHW5lkFlaYJYQmE/y/zDAnaK5JWWD0/6otrs7Qq9G25z2+yvb+Ly2WULgu+KCe19zxYTIdP8bO2Te0c5v/yI37XHa5oYPv9zLLatIu8r9/uI+lV/mz8LckOKUcV+zQKG8lxAA5RHT8ChrjMQCAADAOoRYAAAAWIdyAqAEBXrCXcb0uYySAHPa+tR1rZx2pHFPc3WCk3G+KxCYqxiYspok+T0e/g+3ZCH+dbcUYUvzuk47d537xP77O90SArM0QJJSb3fLF2LX1XLazZ51+7TtIXdDhKZvuiUHZgmBuTJCVhP3Po0Xuq9VcDUCs5wg4fEvnPbe4e2dts+qBwE2NfD5uQCVBNP+qCgYiQUAAIB1CLEAAACwDuUEQAnymaoOcDwvwPEwnyf8XWYJQcHyAfOp/cPN3Kf0q7V0ywZOHnbLA75LbeC0E40pefN4rLGpQbWW7moBWRnGzgeFCJtpfLLRbQbavMBcPSFmhlsakBaoNKDg6xnvoXmeWdphKq2VB2zaIAGVW6CtWAui7ADlHSOxAAAAsA4hFgAAANahnAAoQYGmkc2pZp9pdaNdNcAKBuYUecH7nzBWLjBLDfJWu1P/eRvcKfqcMe4UveSuCmBuXJB2lXtGvUk1/J4vSWnGtg05ASoNmj37k9tXo/Th/8a857Rn3tbNvcD4vqtnuKsfFORTKhBgg4lApR2lNdVfnPtSioDyyCw7oLQA5REjsQAAALAOI7EAAIQQo5zA+SHEAsUUzFSwebyq0c42Fv2PMKbFzZKDyE17nHZYgSlyc9H/vLrunH5ax1pOu57c+0YYKw/s6un++tdMclch0Dp3tQBzQ4PqP7rlA5KUG+eus1Atw53UMVc6uER7nba5YcH4pX2ddvylvpso5DO/t4IrG0QEKCEwmeULVYu4qUFZT+9TQgAARUc5AQAAAKzDSCwAACEU7LqtJYHSBVQkhFigmMyp4KJOQ8fMWOO0zU0QzOnywkQE2CDBFLbfrSGonhHrtOutznba2/7sTtcnrnVXIUhNdu9Tb22B1QmucssOYne4ZQcRh92yA7PEIauJWzbQaIlbKmCWGZglBGYZxd6+ST6vXS/A922+/+HG9af9nh345xXoZ1oQZQAAEDqUEwAAAMA6jMQCAFDKmMYHSt55hdiff/5ZCxYs0K5du5Sb6zvF+OKLL5ZIxwAbBbPZgTnFHmhxfrO0wGQ+cS/5lhOYT/AfaZLn9xpzut7UdOA6p33qulZOO3adWxpwuJnvNUlvGJsrBFgZwWw3Xpjp93jtbe7fIWYJgdlvc8OHgsz+KohVCMLMUg2j1CKQYEsG2LAAAMpWkUPs8uXL1bt3byUnJ2vr1q26/PLLlZqaKq/XqyuvvLI0+ggAAAD4KHKIHT16tB588EGNHz9e0dHRmjNnjurUqaM77rhDPXr0KI0+AgBQLlEmAIROkUPsli1b9J///OeXi8PDdeLECdWsWVOPP/64brrpJt13330l3knAdj7Ty0E8+R4WYJq74GoEZjlCpHGvZpvcc8xpeXODg0ZLjA0EjCl58zWOGOebGxpIUkbnRKd9Ms5decAsDzDvlXlp7DnPMftqHjfLDySpwRy3lOGEcd/IAFP6Pu+z8d4GU/5R8JxgVjQI9l4AgPNX5NUJatSo4dTB1q9fXzt27HC+dvDgwZLrGQAAABBAkUdir7rqKq1evVqXXHKJevbsqQceeEDffvut3n//fV111VWl0UcAAMolc6MCSguAslXkEPviiy/q6NGjkqTx48fr6NGjmjVrlpo2bcrKBIAhUKlAMIvqB5yaLvB5jlEGEBZgpQJzWr5aRi2nXdV4kt98Yt8sZcirFem06y1xSwakwCsJBCohMEsOam12z/EpnTD6bV5bcHUCczWEWpuznHag9y3QihCBmPcJK7DxxOkgrg90LwBAySlyiE1OdrfwqVGjhl599dUS7RAAAGWNUVTAPkWuiU1OTtahQ4fOOp6ZmekTcAEAAIDSUuQQm5qaqjNnzpx1PCcnR3v27PFzBQAAAFCygi4nWLBggdP+8MMPFRvr1qudOXNGy5cvV2JiYol2DrBZUWshA51v1mQW3MnLrD8NtNyWuRRWIIGW9Eqc6fF7vCBz6a74L/3XwZrLYm17qLHf8+M+TXXabjXu2TuVBarnNfnUGAdRx+rzPhvnF3xvzJ9BMHXPKN8oIwDsFXSI7dOnjyTJ4/FowIABPl+rWrWqEhMT9cILL5Ro50x///vftXjxYm3cuFHVqlVTZmZmqb0WAAAAyregQ2xe3i/jD0lJSfryyy91wQUXlFqn/MnNzVW/fv3Url07TZ48uUxfGwAAAOVLkVcn2LnTnRI8efKkqlevXqIdCmT8+PGSpGnTppXJ6wFlIZjp6LBCdpgyrzenuc3z4j6V37YClBzs7+buptXsn747hJl2Dkly2tHunieqtTnT7ZMxLZ9nTNc3fTM74H3zmSUEZvmAFPj79iktCFAiEWgpM58SAuM+hZUiUDZQflAWAFQ+RX6wKy8vT0888YQuvPBC1axZUz/++KMk6bHHHmOEFAAAAGWiyCH2ySef1LRp0/Tss8+qWrVqzvHLL79cb7zxRol2rrhycnKUnZ3t8wEAAAD7ebxer7coF1x00UX697//rS5duig6Olpff/21kpOTtXXrVrVr106HD/ufwvNn1KhReuaZZwo9Z8uWLbr44oudz6dNm6YRI0YE9WDXuHHjnDIEUyfdpHBP1aD7CQDl1Unvca3WEmVlZSkmJibU3QGAYsvOzlZsbOw5/14rck3snj17dNFFF511PC8vT6dOnSrSvR544AENHDiw0HOKs4HC6NGjNXLkSOfz7OxsJSQknPf9AAAAUD4UOcReeuml+uyzz9S4cWOf4++9955atGhRpHvFx8crPj6+qF0IWkREhCIiIkrt/gAAAAiNIofYMWPGaMCAAdqzZ4/y8vL0/vvva9u2bXrrrbe0aNGi0uijJGnXrl3KyMjQrl27dObMGW3cuFHSL+UNNWvWLLXXBQAAQPlT5Ae7brrpJi1cuFAff/yxatSooTFjxmjLli1auHChunXrVhp9lPRLeG7RooXGjh2ro0ePqkWLFmrRooW++uqrc18MAACACqXID3bZLL9QmAe7AFQUPNgFoKIptQe78uXm5urAgQPOTl75GjVqdL63BAAAAIJS5BD7ww8/6K677tIXX3zhc9zr9crj8ejMmTMl1jkAAADAnyKH2IEDByo8PFyLFi1S/fr15fF4SqNfAAAAQEBFDrEbN27UunXrfDYgAAAAAMpSkVcnuPTSS3Xw4MHS6AsAAAAQlKBCbHZ2tvPxzDPP6KGHHtKKFSt06NAhn69lZ2eXdn8BAACA4MoJatWq5VP76vV61aVLF59zeLALAAAAZSWoEPvpp5867dTUVCUkJKhKlSo+5+Tl5WnXrl0l2zsAAADAjyJvdlClShXt27dPderU8Tl+6NAh1alTp1yPxLLZAYCKhs0OAFQ0wW52UOQHu/LLBgo6evSoqlevXtTbAQAAAEUW9BJbI0eOlCR5PB499thjioqKcr525swZ/fe//1Xz5s1LvIMAAABAQUGH2A0bNkj6ZST222+/VbVq1ZyvVatWTSkpKXrwwQdLvocAAABAAUGH2PyHuwYNGqR//vOf1F4BAAAgZIq8Y9fUqVNLox8AAABA0Ir8YBcAAAAQaoRYAAAAWIcQCwAAAOsQYgEAAGAdQiwAAACsQ4gFAACAdQixAAAAsA4hFgAAANYhxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwDiEWAAAA1iHEAgAAwDqEWAAAAFiHEAsAAADrEGIBAABgHUIsAAAArEOIBQAAgHUIsQAAALAOIRYAAADWIcQCAADAOoRYAAAAWIcQCwAAAOtYEWJTU1M1ePBgJSUlKTIyUk2aNNHYsWOVm5sb6q4BAAAgBMJD3YFgbN26VXl5efr3v/+tiy66SJs2bdLQoUN17NgxPf/886HuHgAAAMqYFSG2R48e6tGjh/N5cnKytm3bpldeeYUQCwAAUAlZEWL9ycrKUlxcXKHn5OTkKCcnx/k8Ozu7tLsFAACAMmBFTWxB27dv18SJE3XPPfcUet6ECRMUGxvrfCQkJJRRDwEAAFCaQhpiR40aJY/HU+jH1q1bfa7Zs2ePevTooX79+mno0KGF3n/06NHKyspyPnbv3l2a3w4AAADKSEjLCR544AENHDiw0HOSk5Od9t69e9W5c2e1b99er7322jnvHxERoYiIiOJ2EwAAAOVMSENsfHy84uPjgzp3z5496ty5s1q2bKmpU6cqLMzKSggAAACUACse7NqzZ486deqkxo0b6/nnn1d6errztXr16oWwZwAAAAgFK0LssmXLtH37dm3fvl0NGzb0+ZrX6w1RrwAAABAqVszJDxw4UF6v1+8HAAAAKh8rQiwAAABgIsQCAADAOoRYAAAAWIcQCwAAAOsQYgEAAGAdQiwAAACsQ4gFAACAdQixAAAAsA4hFgAAANYhxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwDiEWAAAA1iHEAgAAwDqEWAAAAFiHEAsAAADrEGIBAABgHUIsAAAArEOIBQAAgHUIsQAAALAOIRYAAADWIcQCAADAOoRYAAAAWIcQCwAAAOsQYgEAAGAdQiwAAACsQ4gFAACAdQixAAAAsA4hFgAAANYhxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwjjUhtnfv3mrUqJGqV6+u+vXr684779TevXtD3S0AAACEgDUhtnPnznr33Xe1bds2zZkzRzt27NAtt9wS6m4BAAAgBMJD3YFg/eUvf3HajRs31qhRo9SnTx+dOnVKVatWDWHPAAAAUNasCbGmjIwMzZgxQ+3bty80wObk5CgnJ8f5PDs7uyy6BwAAgFJmTTmBJP3tb39TjRo19Ktf/Uq7du3S/PnzCz1/woQJio2NdT4SEhLKqKcAAAAoTSENsaNGjZLH4yn0Y+vWrc75f/3rX7VhwwZ99NFHqlKlivr37y+v1xvw/qNHj1ZWVpbzsXv37rL4tgAAAFDKPN7CUmApS09P16FDhwo9Jzk5WdWqVTvr+M8//6yEhAR98cUXateuXVCvl52drdjYWHXSTQr3UEcLwH4nvce1WkuUlZWlmJiYUHcHAIotP6+d6++1kNbExsfHKz4+/ryuzcvLkySfmlcAAABUDlY82PXf//5XX375pTp27KjatWtrx44deuyxx9SkSZOgR2EBAABQcVjxYFdUVJTef/99denSRc2aNdPgwYP1m9/8RitXrlRERESouwcAAIAyZsVI7BVXXKFPPvkk1N0AAABAOWHFSCwAAABgIsQCAADAOoRYAAAAWIcQCwAAAOsQYgEAAGAdQiwAAACsQ4gFAACAdQixAAAAsA4hFgAAANYhxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwDiEWAAAA1iHEAgAAwDqEWAAAAFiHEAsAAADrEGIBAABgHUIsAAAArEOIBQAAgHUIsQAAALAOIRYAAADWIcQCAADAOoRYAAAAWIcQCwAAAOsQYgEAAGAdQiwAAACsQ4gFAACAdQixAAAAsA4hFgAAANYhxAIAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwDiEWAAAA1rEuxObk5Kh58+byeDzauHFjqLsDAACAELAuxD700ENq0KBBqLsBAACAELIqxC5dulQfffSRnn/++VB3BQAAACEUHuoOBGv//v0aOnSo5s2bp6ioqKCuycnJUU5OjvN5VlaWJOm0TkneUukmAJSpU8qVJGVnZ4e4JwBQMvL/PvN6Cw9rVoRYr9ergQMH6t5771WrVq2Umpoa1HUTJkzQ+PHjzzq+WktKuIcAEFoJCQmh7gIAlKgjR44oNjY24Nc93nPF3FI0atQoPfPMM4Wes2XLFn300Ud69913tXLlSlWpUkWpqalKSkrShg0b1Lx584DXFhyJzczMVOPGjbVr165C3xS4srOzlZCQoN27dysmJibU3bEG71vR8Z6dn8OHDysxMVGpqamqXbt2qLsDAMXm9Xp15MgRNWjQQGFhgStfQxpi09PTdejQoULPSU5O1u9//3stXLhQHo/HOX7mzBlVqVJFd9xxh958882gXi87O1uxsbHKysrif5JB4j07P7xvRcd7dn543wBUViEtJ4iPj1d8fPw5z3v55Zf15JNPOp/v3btX3bt316xZs9S2bdvS7CIAAADKIStqYhs1auTzec2aNSVJTZo0UcOGDUPRJQAAAISQVUtsFVdERITGjh2riIiIUHfFGrxn54f3reh4z84P7xuAyiqkNbEAAADA+ahUI7EAAACoGAixAAAAsA4hFgAAANYhxAIAAMA6lTrELl68WG3btlVkZKRq166tPn36hLpL1sjJyVHz5s3l8Xi0cePGUHen3EpNTdXgwYOVlJSkyMhINWnSRGPHjlVubm6ou1buTJo0SYmJiapevbratm2r//3vf6HuUrk1YcIEtW7dWtHR0apTp4769Omjbdu2hbpbAFCmKm2InTNnju68804NGjRIX3/9tT7//HPdfvvtoe6WNR566CE1aNAg1N0o97Zu3aq8vDz9+9//1nfffad//OMfevXVV/Xwww+HumvlyqxZszRy5EiNHTtW69evV0pKirp3764DBw6Eumvl0sqVKzVs2DCtXbtWy5Yt06lTp3Tdddfp2LFjoe4aAJSZSrnE1unTp5WYmKjx48dr8ODBoe6OdZYuXaqRI0dqzpw5uuyyy7RhwwY1b9481N2yxnPPPadXXnlFP/74Y6i7Um60bdtWrVu31r/+9S9JUl5enhISEjR8+HCNGjUqxL0r/9LT01WnTh2tXLlSV199dai7AwBlolKOxK5fv1579uxRWFiYWrRoofr16+v666/Xpk2bQt21cm///v0aOnSopk+frqioqFB3x0pZWVmKi4sLdTfKjdzcXK1bt05du3Z1joWFhalr165as2ZNCHtmj6ysLEnizxWASqVShtj8EbBx48bp0Ucf1aJFi1S7dm116tRJGRkZIe5d+eX1ejVw4EDde++9atWqVai7Y6Xt27dr4sSJuueee0LdlXLj4MGDOnPmjOrWretzvG7dukpLSwtRr+yRl5enESNGqEOHDrr88stD3R0AKDMVKsSOGjVKHo+n0I/8GkVJeuSRR9S3b1+1bNlSU6dOlcfj0ezZs0P8XZS9YN+3iRMn6siRIxo9enSouxxywb5npj179qhHjx7q16+fhg4dGqKeo6IZNmyYNm3apHfeeSfUXQGAMhUe6g6UpAceeEADBw4s9Jzk5GTt27dPknTppZc6xyMiIpScnKxdu3aVZhfLpWDft08++URr1qw5a4/2Vq1a6Y477tCbb75Zir0sX4J9z/Lt3btXnTt3Vvv27fXaa6+Vcu/scsEFF6hKlSrav3+/z/H9+/erXr16IeqVHe6//34tWrRIq1atUsOGDUPdHQAoUxUqxMbHxys+Pv6c57Vs2VIRERHatm2bOnbsKEk6deqUUlNT1bhx49LuZrkT7Pv28ssv68knn3Q+37t3r7p3765Zs2apbdu2pdnFcifY90z6ZQS2c+fOzoh/WFiFmgAptmrVqqlly5Zavny5s8xdXl6eli9frvvvvz+0nSunvF6vhg8frrlz52rFihVKSkoKdZcAoMxVqBAbrJiYGN17770aO3asEhIS1LhxYz333HOSpH79+oW4d+VXo0aNfD6vWbOmJKlJkyaMAgWwZ88ederUSY0bN9bzzz+v9PR052uMMrpGjhypAQMGqFWrVmrTpo1eeuklHTt2TIMGDQp118qlYcOGaebMmZo/f76io6Od2uHY2FhFRkaGuHcAUDYqZYiVflnmKDw8XHfeeadOnDihtm3b6pNPPlHt2rVD3TVUIMuWLdP27du1ffv2s4J+JVzdLqBbb71V6enpGjNmjNLS0tS8eXN98MEHZz3shV+88sorkqROnTr5HJ86deo5y1wAoKKolOvEAgAAwG4U5wEAAMA6hFgAAABYhxALAAAA6xBiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwDiEWAAAA1iHEAiVg3759uv322/XrX/9aYWFhGjFihN/zZs+erYsvvljVq1fXFVdcoSVLlpRtRwEAqCAIsUAJyMnJUXx8vB599FGlpKT4PeeLL77QH/7wBw0ePFgbNmxQnz591KdPH23atKmMewsAgP3YdhYIQnp6uq644gr96U9/0sMPPyzpl1DaqVMnLV26VF26dHHO7dSpk5o3b66XXnrJ5x633nqrjh07pkWLFjnHrrrqKjVv3lyvvvpqmXwfAABUFIzEAkGIj4/XlClTNG7cOH311Vc6cuSI7rzzTt1///0+AbYwa9asUdeuXX2Ode/eXWvWrCmNLgMAUKGFh7oDgC169uypoUOH6o477lCrVq1Uo0YNTZgwIejr09LSVLduXZ9jdevWVVpaWkl3FQCACo+RWKAInn/+eZ0+fVqzZ8/WjBkzFBEREeouAQBQKRFigSLYsWOH9u7dq7y8PKWmphbp2nr16mn//v0+x/bv36969eqVYA8BAKgcCLFAkHJzc/XHP/5Rt956q5544gkNGTJEBw4cCPr6du3aafny5T7Hli1bpnbt2pV0VwEAqPCoiQWC9MgjjygrK0svv/yyatasqSVLluiuu+5yVhvYuHGjJOno0aNKT0/Xxo0bVa1aNV166aWSpD//+c+65ppr9MILL+iGG27QO++8o6+++kqvvfZaqL4lAACsxRJbQBBWrFihbt266dNPP1XHjh0lSampqUpJSdHTTz+t++67Tx6P56zrGjdu7FN2MHv2bD366KNKTU1V06ZN9eyzz6pnz55l9W0AAFBhEGIBAABgHWpiAQAAYB1CLAAAAKxDiAUAAIB1CLEAAACwDiEWAAAA1iHEAgAAwDqEWAAAAFiHEAsAAADrEGIBAABgHUIsAAAArEOIBQAAgHUIsQAAALDO/wN6ia7y6q6icAAAAABJRU5ErkJggg==", + "text/plain": [ + "Figure(PyObject
)" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig, axis = plt.subplots(2, 2, figsize=(8,8))\n", + "fig.suptitle(\"HMCDA - 21-D Neal's Funnel\", fontsize=16)\n", + "\n", + "fig.delaxes(axis[1,2])\n", + "fig.subplots_adjust(hspace=0)\n", + "fig.subplots_adjust(wspace=0)\n", + "\n", + "axis[1,1].hist(x10_hmcda, bins=100, range=[-6,2])\n", + "axis[1,1].set_yticks([])\n", + "\n", + "axis[2,2].hist(theta_hmcda, bins=100, orientation=\"horizontal\", range=[-4, 2])\n", + "axis[2,2].set_xticks([])\n", + "axis[2,2].set_yticks([])\n", + "\n", + "axis[2,1].hist2d(x10_hmcda, theta_hmcda, bins=100, range=[[-6,2],[-4, 2]])\n", + "axis[2,1].set_xlabel(\"x10\")\n", + "axis[2,1].set_ylabel(\"theta\");" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "843becb3", + "metadata": {}, + "source": [] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "91baadc8", + "metadata": {}, + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 1.9.0", + "language": "julia", + "name": "julia-1.9" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.9.0" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/src/AdvancedHMC.jl b/src/AdvancedHMC.jl index 28bc440f..df123320 100644 --- a/src/AdvancedHMC.jl +++ b/src/AdvancedHMC.jl @@ -169,6 +169,7 @@ include("sampler.jl") export sample include("abstractmcmc.jl") +include("constructors.jl") ## Without explicit AD backend function Hamiltonian(metric::AbstractMetric, ℓ::LogDensityModel; kwargs...) @@ -263,4 +264,4 @@ function __init__() end end -end # module +end # module \ No newline at end of file diff --git a/src/abstractmcmc.jl b/src/abstractmcmc.jl index e491b53b..b4344ef0 100644 --- a/src/abstractmcmc.jl +++ b/src/abstractmcmc.jl @@ -1,30 +1,3 @@ -""" - HMCSampler - -A `AbstractMCMC.AbstractSampler` for kernels in AdvancedHMC.jl. - -# Fields - -$(FIELDS) - -# Notes - -Note that all the fields have the prefix `initial_` to indicate -that these will not necessarily correspond to the `kernel`, `metric`, -and `adaptor` after sampling. - -To access the updated fields use the resulting [`HMCState`](@ref). -""" -struct HMCSampler{K,M,A} <: AbstractMCMC.AbstractSampler - "Initial [`AbstractMCMCKernel`](@ref)." - initial_kernel::K - "Initial [`AbstractMetric`](@ref)." - initial_metric::M - "Initial [`AbstractAdaptor`](@ref)." - initial_adaptor::A -end -HMCSampler(kernel, metric) = HMCSampler(kernel, metric, Adaptation.NoAdaptation()) - """ HMCState @@ -53,140 +26,39 @@ struct HMCState{ adaptor::TAdapt end -""" - $(TYPEDSIGNATURES) - -A convenient wrapper around `AbstractMCMC.sample` avoiding explicit construction of [`HMCSampler`](@ref). -""" -function AbstractMCMC.sample( - model::LogDensityModel, - kernel::AbstractMCMCKernel, - metric::AbstractMetric, - adaptor::AbstractAdaptor, - N::Integer; - kwargs..., -) - return AbstractMCMC.sample( - Random.GLOBAL_RNG, - model, - kernel, - metric, - adaptor, - N; - kwargs..., - ) -end - -function AbstractMCMC.sample( - rng::Random.AbstractRNG, - model::LogDensityModel, - kernel::AbstractMCMCKernel, - metric::AbstractMetric, - adaptor::AbstractAdaptor, - N::Integer; - progress = true, - verbose = false, - callback = nothing, - kwargs..., -) - sampler = HMCSampler(kernel, metric, adaptor) - if callback === nothing - callback = HMCProgressCallback(N, progress = progress, verbose = verbose) - progress = false # don't use AMCMC's progress-funtionality - end - - return AbstractMCMC.mcmcsample( - rng, - model, - sampler, - N; - progress = progress, - verbose = verbose, - callback = callback, - kwargs..., - ) -end - -function AbstractMCMC.sample( - model::LogDensityModel, - kernel::AbstractMCMCKernel, - metric::AbstractMetric, - adaptor::AbstractAdaptor, - parallel::AbstractMCMC.AbstractMCMCEnsemble, - N::Integer, - nchains::Integer; - kwargs..., -) - return AbstractMCMC.sample( - Random.GLOBAL_RNG, - model, - kernel, - metric, - adaptor, - N, - nchains; - kwargs..., - ) -end - -function AbstractMCMC.sample( - rng::Random.AbstractRNG, - model::LogDensityModel, - kernel::AbstractMCMCKernel, - metric::AbstractMetric, - adaptor::AbstractAdaptor, - parallel::AbstractMCMC.AbstractMCMCEnsemble, - N::Integer, - nchains::Integer; - progress = true, - verbose = false, - callback = nothing, - kwargs..., -) - sampler = HMCSampler(kernel, metric, adaptor) - if callback === nothing - callback = HMCProgressCallback(N, progress = progress, verbose = verbose) - progress = false # don't use AMCMC's progress-funtionality - end - - return AbstractMCMC.mcmcsample( - rng, - model, - sampler, - parallel, - N, - nchains; - progress = progress, - verbose = verbose, - callback = callback, - kwargs..., - ) -end - function AbstractMCMC.step( rng::AbstractRNG, - model::LogDensityModel, - spl::HMCSampler; + model::AbstractMCMC.LogDensityModel, + spl::AbstractMCMC.AbstractSampler; init_params = nothing, kwargs..., -) - metric = spl.initial_metric - κ = spl.initial_kernel - adaptor = spl.initial_adaptor +) + # Unpack model + logdensity = model.logdensity + vi = logdensity.varinfo - if init_params === nothing - init_params = randn(rng, size(metric, 1)) - end + # Define metric + metric = make_metric(spl, logdensity) # Construct the hamiltonian using the initial metric hamiltonian = Hamiltonian(metric, model) + # Define integration algorithm + # Find good eps if not provided one + integrator = make_integrator(rng, spl, hamiltonian, init_params) + + # Make kernel + κ = make_kernel(spl, integrator) + + # Make adaptor + n_adapts, adaptor = make_adaptor(spl, metric, integrator) + # Get an initial sample. h, t = AdvancedHMC.sample_init(rng, hamiltonian, init_params) # Compute next transition and state. - state = HMCState(0, t, h.metric, κ, adaptor) - + state = HMCState(0, t, metric, κ, adaptor) + # Take actual first step. return AbstractMCMC.step(rng, model, spl, state; kwargs...) end @@ -194,7 +66,7 @@ end function AbstractMCMC.step( rng::AbstractRNG, model::LogDensityModel, - spl::HMCSampler, + spl::AbstractMCMC.AbstractSampler, state::HMCState; nadapts::Int = 0, kwargs..., @@ -302,4 +174,4 @@ function (cb::HMCProgressCallback)(rng, model, spl, t, state, i; nadapts = 0, kw elseif verbose && isadapted && i == nadapts @info "Finished $nadapts adapation steps" adaptor κ.τ.integrator metric end -end +end \ No newline at end of file diff --git a/src/constructors.jl b/src/constructors.jl new file mode 100644 index 00000000..9d33eff6 --- /dev/null +++ b/src/constructors.jl @@ -0,0 +1,205 @@ +abstract type AbstractHMCSampler <:AbstractMCMC.AbstractSampler end + +########## +# Custom # +########## +""" + HMCSampler + +A `AbstractMCMC.AbstractSampler` for kernels in AdvancedHMC.jl. + +# Fields + +$(FIELDS) + +# Notes + +Note that all the fields have the prefix `initial_` to indicate +that these will not necessarily correspond to the `kernel`, `metric`, +and `adaptor` after sampling. + +To access the updated fields use the resulting [`HMCState`](@ref). +""" +Base.@kwdef struct CustomHMC{I,K,M,A} <: AbstractMCMC.AbstractSampler + "[`integrator`](@ref)." + integrator::I=Leapfrog + "[`AbstractMCMCKernel`](@ref)." + kernel::K=nothing + "[`AbstractMetric`](@ref)." + metric::M=nothing + "[`AbstractAdaptor`](@ref)." + adaptor::A=nothing +end + +######## +# NUTS # +######## +""" + NUTS(n_adapts::Int, δ::Float64; max_depth::Int=10, Δ_max::Float64=1000.0, init_ϵ::Float64=0.0) + +No-U-Turn Sampler (NUTS) sampler. + +Usage: + +```julia +NUTS() # Use default NUTS configuration. +NUTS(1000, 0.65) # Use 1000 adaption steps, and target accept ratio 0.65. +``` + +Arguments: + +- `n_adapts::Int` : The number of samples to use with adaptation. +- `δ::Float64` : Target acceptance rate for dual averaging. +- `max_depth::Int` : Maximum doubling tree depth. +- `Δ_max::Float64` : Maximum divergence during doubling tree. +- `init_ϵ::Float64` : Initial step size; 0 means automatically searching using a heuristic procedure. + +""" +Base.@kwdef struct NUTS_alg <: AbstractMCMC.AbstractSampler + n_adapts::Int # number of samples with adaption for ϵ + δ::Float64 # target accept rate + max_depth::Int=10 # maximum tree depth + Δ_max::Float64=1000.0 # maximum error + init_ϵ::Float64=0.0 # (initial) step size + integrator_method=Leapfrog # integrator method + metric_type=DiagEuclideanMetric # metric type +end + +####### +# HMC # +####### +""" + HMC(ϵ::Float64, n_leapfrog::Int) + +Hamiltonian Monte Carlo sampler with static trajectory. + +Arguments: + +- `ϵ::Float64` : The leapfrog step size to use. +- `n_leapfrog::Int` : The number of leapfrog steps to use. + +Usage: + +```julia +HMC(0.05, 10) +``` + +Tips: + +- If you are receiving gradient errors when using `HMC`, try reducing the leapfrog step size `ϵ`, e.g. + +```julia +# Original step size +sample(gdemo([1.5, 2]), HMC(0.1, 10), 1000) + +# Reduced step size +sample(gdemo([1.5, 2]), HMC(0.01, 10), 1000) +``` +""" +Base.@kwdef struct HMC_alg <: AbstractMCMC.AbstractSampler + init_ϵ::Float64 # leapfrog step size + n_leapfrog::Int # leapfrog step number + integrator_method=Leapfrog # integrator method + metric_type=DiagEuclideanMetric # metric type +end + +######### +# HMCDA # +######### +""" + HMCDA(n_adapts::Int, δ::Float64, λ::Float64; ϵ::Float64=0.0) + +Hamiltonian Monte Carlo sampler with Dual Averaging algorithm. + +Usage: + +```julia +HMCDA(200, 0.65, 0.3) +``` + +Arguments: + +- `n_adapts::Int` : Numbers of samples to use for adaptation. +- `δ::Float64` : Target acceptance rate. 65% is often recommended. +- `λ::Float64` : Target leapfrog length. +- `ϵ::Float64=0.0` : Initial step size; 0 means automatically search by Turing. + +For more information, please view the following paper ([arXiv link](https://arxiv.org/abs/1111.4246)): + +- Hoffman, Matthew D., and Andrew Gelman. "The No-U-turn sampler: adaptively + setting path lengths in Hamiltonian Monte Carlo." Journal of Machine Learning + Research 15, no. 1 (2014): 1593-1623. +""" +Base.@kwdef struct HMCDA_alg <: AbstractMCMC.AbstractSampler + n_adapts::Int # number of samples with adaption for ϵ + δ::Float64 # target accept rate + λ::Float64 # target leapfrog length + init_ϵ::Float64=0.0 # (initial) step size + integrator_method=Leapfrog # integrator method + metric_type=DiagEuclideanMetric # metric type +end + +export CustomHMC, HMC_alg, NUTS_alg, HMCDA_alg +######### +# Utils # +######### + +function make_integrator(rng, spl::Union{HMC_alg, NUTS_alg, HMCDA_alg}, + hamiltonian, init_params) + init_ϵ = spl.init_ϵ + if iszero(init_ϵ) + init_ϵ = find_good_stepsize(rng, hamiltonian, init_params) + @info string("Found initial step size ", init_ϵ) + end + return spl.integrator_method(init_ϵ) +end + +function make_integrator(rng, spl::CustomHMC, hamiltonian, init_params) + return spl.integrator +end + +######### + +function make_metric(spl::Union{HMC_alg, NUTS_alg, HMCDA_alg}, logdensity) + d = LogDensityProblems.dimension(logdensity) + return spl.metric_type(d) +end + +function make_metric(spl::CustomHMC, logdensity) + return spl.metric +end + +######### + +function make_adaptor(spl::Union{NUTS_alg, HMCDA_alg}, metric, integrator) + adaptor = StanHMCAdaptor(MassMatrixAdaptor(metric), + StepSizeAdaptor(spl.δ, integrator)) + n_adapts = spl.n_adapts + return n_adapts, adaptor + end + +function make_adaptor(spl::HMC_alg, metric, integrator) + return 0, NoAdaptation() + end + + function make_adaptor(spl::CustomHMC, metric, integrator) + return spl.n_adapts, spl.adaptor + end + +######### + +function make_kernel(spl::NUTS_alg, integrator) + return HMCKernel(Trajectory{MultinomialTS}(integrator, GeneralisedNoUTurn())) +end + +function make_kernel(spl::HMC_alg, integrator) + return HMCKernel(Trajectory{EndPointTS}(integrator, FixedNSteps(spl.n_leapfrog))) +end + +function make_kernel(spl::HMCDA_alg, integrator) + return HMCKernel(Trajectory{EndPointTS}(integrator, FixedIntegrationTime(spl.λ))) +end + +function make_kernel(spl::CustomHMC, integrator) + return spl.kernel +end \ No newline at end of file diff --git a/src/sampler.jl b/src/sampler.jl index 7d1b7eb5..d8b63ce8 100644 --- a/src/sampler.jl +++ b/src/sampler.jl @@ -246,4 +246,4 @@ function sample( @info "Finished $n_samples sampling steps for $n_chains chains in $time (s)" h κ EBFMI_est average_acceptance_rate end return θs, stats -end +end \ No newline at end of file