From 8e1662f5ce9f4849d48cb98581a149aa573721f6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C4=B0smail=20=C5=9Een=C3=B6z?= Date: Mon, 13 Jan 2025 15:01:31 +0100 Subject: [PATCH 1/8] add binomial model tests --- test/models/regression/binomialreg_tests.jl | 97 +++++++++++++++++++++ 1 file changed, 97 insertions(+) create mode 100644 test/models/regression/binomialreg_tests.jl diff --git a/test/models/regression/binomialreg_tests.jl b/test/models/regression/binomialreg_tests.jl new file mode 100644 index 000000000..cfede8878 --- /dev/null +++ b/test/models/regression/binomialreg_tests.jl @@ -0,0 +1,97 @@ +@testitem "Linear regression" begin + using BenchmarkTools, Plots, Dates, LinearAlgebra, StableRNGs + + include(joinpath(@__DIR__, "..", "..", "utiltests.jl")) + + function generate_synthetic_binomial_data( + n_samples::Int, + n_features::Int, + true_beta::Vector{Float64}; + seed::Int=42 + ) + rng = StableRNG(seed) + + # Generate design matrix X + X = randn(rng, n_samples, n_features) + + # Generate number of trials for each observation + n_trials = rand(rng, 5:20, n_samples) + + # Compute logits and probabilities + logits = X * true_beta + probs = 1 ./ (1 .+ exp.(-logits)) + + # Generate binomial outcomes + y = [rand(rng, Binomial(n_trials[i], probs[i])) for i in 1:n_samples] + + return X, y, n_trials + end + n_samples = 1000 + n_features = 2 + true_beta = [-1.0 , 0.6] + n_iterations = 100 + n_sims = 20 + + + @model function binomial_model(prior_xi, prior_precision, n_trials, X, y) + β ~ MvNormalWeightedMeanPrecision(prior_xi, prior_precision) + for i in eachindex(y) + y[i] ~ BinomialPolya(X[i], n_trials[i], β) where { + dependencies = RequireMessageFunctionalDependencies(β = MvNormalWeightedMeanPrecision(prior_xi, prior_precision)) + } + end + end + + function binomial_inference(binomial_model, iterations, X, y, n_trials) + return infer( + model = binomial_model(prior_xi = zeros(n_features), prior_precision =diageye(n_features),), + data = (X=X, y=y,n_trials=n_trials), + iterations = iterations, + free_energy = true, + options = (limit_stack_depth = 100,), + ) + end + + function run_simulation(n_sims::Int, n_samples::Int, n_features::Int, true_beta::Vector{Float64};iterations = n_iterations) + # Storage for results + coverage = Vector{Vector{Float64}}(undef, n_sims) + fes = Vector{Vector{Float64}}(undef, n_sims) + for sim in 1:n_sims + # Generate new dataset + X, y, n_trials = generate_synthetic_binomial_data(n_samples, n_features, true_beta, seed = sim) + X = [collect(row) for row in eachrow(X)] + + # Run inference + results = binomial_inference(binomial_model, iterations, X, y, n_trials) + # Extract posterior parameters + post = results.posteriors[:β][end] + m = mean(post) + v = var(post) + estimates = map((x,y) -> Normal(x,sqrt(y)), m, v) + coverage[sim] = map((d,b) -> cdf(d, b), estimates,true_beta) + fes[sim] = results.free_energy + + end + + return coverage, fes + end + + function in_credible_interval(x, lwr=0.025, upr=0.975) + return x >= lwr && x <= upr + end; + + + coverage, fes = run_simulation(n_sims, n_samples, n_features, true_beta) + for i in 1:n_sims + @test fes[i][end] < fes[i][1] + end + coverages = Vector{Float64}(undef, n_features) + for i in 1:n_features + coverages[i] = sum(in_credible_interval.(getindex.(coverage,i))) / n_sims + @test coverages[i] >= 0.8 + end + + @test_benchmark "models" "binomialreg" binomial_inference(binomial_model, $n_iterations, $X, $y, $n_trials) + + +end \ No newline at end of file From 60252b8434d38329ddd681710020bcc3f86977bc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C4=B0smail=20=C5=9Een=C3=B6z?= Date: Mon, 13 Jan 2025 16:47:28 +0100 Subject: [PATCH 2/8] add binomial regression example --- examples/.meta.jl | 6 + .../basic_examples/Binomial Regression.ipynb | 473 ++++++++++++++++++ 2 files changed, 479 insertions(+) create mode 100644 examples/basic_examples/Binomial Regression.ipynb diff --git a/examples/.meta.jl b/examples/.meta.jl index 7d8c78ffc..6e733ca05 100644 --- a/examples/.meta.jl +++ b/examples/.meta.jl @@ -21,6 +21,12 @@ return ( description = "An extensive tutorial on Bayesian linear regression with RxInfer with a lot of examples, including multivariate and hierarchical linear regression.", category = :basic_examples ), + ( + filename = "Binomial Regression.ipynb", + title = "Binomial Regression", + description = "An example of Bayesian inference in Binomial regression with Expectation Propagation.", + category = :basic_examples + ), ( filename = "Kalman filtering and smoothing.ipynb", title = "Kalman filtering and smoothing", diff --git a/examples/basic_examples/Binomial Regression.ipynb b/examples/basic_examples/Binomial Regression.ipynb new file mode 100644 index 000000000..93fd5a144 --- /dev/null +++ b/examples/basic_examples/Binomial Regression.ipynb @@ -0,0 +1,473 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Bayesian Binomial Regression " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This notebook is an introductory tutorial to Bayesian binomial regression with `RxInfer`." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "\u001b[32m\u001b[1m Activating\u001b[22m\u001b[39m project at `~/.julia/dev/RxInfer/examples`\n" + ] + } + ], + "source": [ + "# Activate local environment, see `Project.toml`\n", + "import Pkg; Pkg.activate(\"..\"); Pkg.instantiate(); " + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "using RxInfer, ReactiveMP,Random, Plots, StableRNGs, LinearAlgebra, StatsPlots, LaTeXStrings" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Likelihood Specification\n", + "For observations $y_i$ with predictors $\\mathbf{x}_i$, Binomial regression models the number of successes $y_i$ as a function of the predictors $\\mathbf{x}_i$ and the regression coefficients $\\boldsymbol{\\beta}$\n", + "\n", + "$$\\begin{equation}\n", + "y_i \\sim \\text{Binomial}(n_i, p_i)\\,,\n", + "\\end{equation}$$\n", + "where:\n", + "$y_i$ is the number of successes, $n_i$ is the number of trials, $p_i$ is the probability of success. The probability $p_i$ is linked to the predictors through the logistic function:\n", + "$$\\begin{equation}\n", + "p_i = \\frac{1}{1 + e^{-\\mathbf{x}_i^T\\boldsymbol{\\beta}}}\n", + "\\end{equation}$$\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Prior Distributions\n", + "We specify priors for the regression coefficients:\n", + "\n", + "$$\\begin{equation}\n", + "\\boldsymbol{\\beta} \\sim \\mathcal{N}_{\\xi}(\\boldsymbol{\\xi}, \\boldsymbol{\\Lambda})\n", + "\\end{equation}$$\n", + "as a Normal distribution in precision-weighted mean form.\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Model Specification\n", + "\n", + "The likelihood and the prior distributions form the probabilistic model\n", + "$$p(y, x, \\beta, n) = p(\\beta) \\prod_{i=1}^N p(y_i \\mid x_i, \\beta, n_i),$$\n", + "where the goal is to infer the posterior distributions $p(\\beta \\mid y, x, n)$. Due to logistic link function, the posterior distribution is not conjugate to the prior distribution. This means that we need to use a more complex inference algorithm to infer the posterior distribution. Before dwelling into the details of the inference algorithm, let's first generate some synthetic data to work with." + ] + }, + { + "cell_type": "code", + "execution_count": 61, + "metadata": {}, + "outputs": [], + "source": [ + "function generate_synthetic_binomial_data(\n", + " n_samples::Int,\n", + " true_beta::Vector{Float64};\n", + " seed::Int=42\n", + ")\n", + " n_features = length(true_beta)\n", + " rng = StableRNG(seed)\n", + " \n", + " X = randn(rng, n_samples, n_features)\n", + " \n", + " n_trials = rand(rng, 5:20, n_samples)\n", + " \n", + " logits = X * true_beta\n", + " probs = 1 ./ (1 .+ exp.(-logits))\n", + " \n", + " y = [rand(rng, Binomial(n_trials[i], probs[i])) for i in 1:n_samples]\n", + " \n", + " return X, y, n_trials\n", + "end\n", + "\n", + "\n", + "n_samples = 10000\n", + "true_beta = [-3.0 , 2.6]\n", + "\n", + "X, y, n_trials = generate_synthetic_binomial_data(n_samples, true_beta);\n", + "X = [collect(row) for row in eachrow(X)];" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We generate `X` as the design matrix and `y` as the number of successes and `n_trials` as the number of trials. Next task is to define the graphical model. RxInfer provides a `BinomialPolya` factor node that is a combination of a Binomial distribution and a PolyaGamma distribution introduced in [1]. The `BinomialPolya` factor node is used to model the likelihood of the binomial distribution. \n", + "\n", + "Due to non-conjugacy of the likelihood and the prior distribution, we need to use a more complex inference algorithm. RxInfer provides an Expectation Propagation (EP) [2] algorithm to infer the posterior distribution. Due to EP's approximation, we need to specify an inbound message for the regression coefficients while using the `BinomialPolya` factor node. " + ] + }, + { + "cell_type": "code", + "execution_count": 62, + "metadata": {}, + "outputs": [], + "source": [ + "@model function binomial_model(prior_xi, prior_precision, n_trials, X, y) \n", + " β ~ MvNormalWeightedMeanPrecision(prior_xi, prior_precision)\n", + " for i in eachindex(y)\n", + " y[i] ~ BinomialPolya(X[i], n_trials[i], β) where {\n", + " dependencies = RequireMessageFunctionalDependencies(β = MvNormalWeightedMeanPrecision(prior_xi, prior_precision))\n", + " }\n", + " end\n", + "end" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Havin specified the model, we can now utilize the `infer` function to infer the posterior distribution." + ] + }, + { + "cell_type": "code", + "execution_count": 80, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "\u001b[32mProgress: 100%|█████████████████████████████████████████| Time: 0:00:02\u001b[39m\u001b[K\n" + ] + }, + { + "data": { + "text/plain": [ + "Inference results:\n", + " Posteriors | available for (β)\n", + " Free Energy: | Real[21992.9, 16235.8, 13785.0, 12519.7, 11800.1, 11366.2, 11094.1, 10918.7, 10803.3, 10726.3 … 10558.2, 10558.2, 10558.2, 10558.1, 10558.1, 10558.1, 10558.1, 10558.1, 10558.1, 10558.1]\n" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "n_features = length(true_beta)\n", + "n_max_iterations = 50\n", + "strategy = StopEarlyIterationStrategy(1e-5, [])\n", + "results = infer(\n", + " model = binomial_model(prior_xi = zeros(n_features), prior_precision = diageye(n_features),),\n", + " data = (X=X, y=y,n_trials=n_trials),\n", + " iterations = n_max_iterations,\n", + " free_energy = true,\n", + " options = (limit_stack_depth = 100,),\n", + " showprogress = true,\n", + " \n", + ")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can now plot the free energy to see if the inference algorithm is converging." + ] + }, + { + "cell_type": "code", + "execution_count": 81, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3deXwUVb738VPV3Vk7ndUkkLAlQcANhaCyCAwQ9QpBQRgQwXFDHF/3+uCjw+iAOj6My1zEG/W6juCGuIDKCIqKgGIQIqCIsgkGQSCQhJCt03vV80dhG0ISOqZJdXd93n/4qjpUqn9ddve365yq05KqqgIAAKOS9S4AAAA9EYQAAEMjCAEAhkYQAgAMjSAEABgaQQgAMDSCEABgaAQhAMDQCEIAgKERhAAAQwunIHz22Wf37NnjX/V4PDoWE444Ym3l9XqZg7BNfD6foih6VxFOFEXx+Xx6VxFmvF5vcHcYTkG4cuXK3bt3+1edTqeOxYQjjlhbuVwugrBNPB4PH+tt4vP5+IbaJoqiuFyu4O4znIIQAICgIwgBAIZGEAIADI0gBAAYGkEIADA0ghAAYGgEIQDA0AhCAIChRWwQfnJQPeLQuwgAQMiL2CBcsFtZV8ZUTwCA0zDrXcCZkhglatx6FwEAequtrb399tuDPj9nx7NarQsXLjwTe47YIEyKEtUEIQDDq6ioWLNmzdNPP613Ie01adKkl1566UzsOWKDMDFKqnEzXTIACKvVOnHiRL2raK/JkyefoT1H7BghXaMAgEAQhAAAQ4voIORHvgAApxOxQZgUJVW7GCMEAJxGxAYhXaMAgEBEdBDSNQoAOJ2IvX2CrlEACCOHDx9et27d9u3bc3Nzb7zxxo586IgNQluUqPMIVQhJ70oAAKe1aNGitWvX2u32b7/9toODMGK7Rk2SiDOLOnpHASCU7NixY+XKlf7VXbt2ffjhh0KIWbNmrVy5cuzYsR1fUsQGoWByGQAIPfHx8dddd11dXZ22ev/99+/cuVPfkiK2a1T8euFol3i96wCAUPKLXb1lna9jzhKiZLFkpDmuUdR069Zt0KBBb7311vTp048ePfrxxx8/88wzHVJLiyI/CAEAjXWKlf7a16R2SBJGm0TcKTnz5z//ee7cudOnT1+4cGFhYWF6enpHlNKySA5CfoACAE5llsXIznpeRzh69Og777xzy5YtCxcuPEM/KNEmjBECADqULMu33nrrzTffbDabhw4dqnc5ER6EdI0CQCi67bbbdu/efdttt0nSiXPTf//737m5uY899tjatWtzc3PvvffeDiuGrlEAQEdzOBxRUVE33HCDv+U//uM/Gp8dRkdHd1gxkRyEiVFSNV2jABBili1b9vTTT998882pqan+xqioqKioKF3qiewgFPvr9S4CAHCy7du3FxYWzpgxQ+9CTojwIGSMEABCzezZs/Uu4SSRfLFMEl2jAIDTieQg5IwQAHBaBCEAwNAIQgCAoUXyxTKMEQKAyWQ6ePBgfn6+3oW0lyRJkiSpZ2CO1EgOQqtFOH3CqwhzJJ/3AkBrunfvXlJS4vV69S6kveLi4gjCNpOESLCIWo9I6bgJCgAg5PTt21fvEkJahJ8rMe82AKB1kR6EFq6XAQC0JsKDMCmaebcBAK2J8CCkaxQA0LpID0K6RgEArYr0IOSeegBAqyI8CBkjBAC0LsKDkDFCAEDrIj0IGSMEALQq0oOQMUIAQKsiPAiTopl3GwDQmggPQrpGAQCti/QgpGsUANCqCA9CukYBAK2L8CCkaxQA0LoID8JYs1CFcPn0rgMAEKoiPAiFEDZOCgEALYv8IGRyGQBAK4wQhKLGo3cRAIBQFflBmBQlql16FwEACFWRH4R0jQIAWmGEIKRrFADQosgPQrpGAQCtiPwgpGsUANCKQIPQ6XRu2LBh5cqVBw8ebNyuKMrGjRs/+uijmpqaxu2VlZUrVqz45ptvmuxny5YtK1asOHbsWOPG6urqjz76qKSkRFGUtj+F06BrFADQioCCcPfu3Z06dZo5c+azzz573nnnzZs3T2v3+XyFhYXTp09//vnne/Xq9cMPP2jtxcXFffr0WbBgweTJk//0pz/59zN16tTrrrtuwYIFffr0+eqrr7TGbdu29erV6/nnn7/llluuvvpqny/I08Aw7zYAoDVqACorK3/66Sdtef369bIsV1VVqaq6bNmynj17NjQ0qKo6e/bs8ePHa9sMHjz4ySefVFX1+PHjGRkZGzZs0P4wIyPj+PHjqqoWFRUNHTpU2/iaa66ZM2eOqqoNDQ15eXkffPBBS2WMGTNm+fLl/tXa2tpAin9vn2/cKm8gW0a8AI8Y/Orr630+n95VhBOHw+F2u/WuIpy43W6Hw6F3FeHE5/PV19cHd58BnRGmpqbm5ORoy7169VIUxW63CyHef//98ePHx8bGCiGmTJnywQcfeL3eo0ePrl+/fsqUKUKIpKSkq6666v3339c2Hj16dFJSkrbxunXrKisrPR7PihUrrr/+eiFEbGzsuHHjtI2DKDFKqnYxRggAaJ65rX/wxBNPDB06NDs7Wwhx8ODBAQMGaO3dunXzer1HjhwpLy+Pi4tLS0vT2rt27bpnzx4hxC+//NK7d2+t8ayzzoqNjT148KDD4fB6vV27dvVvvHXr1pYeura2tri4uKGhQVuVZXn8+PGnLTjBrNZ4xJkYfQw7iqJwHNqEI9ZWiqJIksRBC5zyK70LCRttPWKSJEmS1Po2bQvCN99885VXXvnyyy+1VbfbbbFYtGVtweVyuVwus/m33VosFqfT2WRjIURUVJS2sRDCv71/42YdP3583bp1WqwKIWw221VXXXXammNUUe1qbbfG4XK5Gv8vwGk5nU5JkmQ58i+uDhan02kymYI+0h/BPB4Ph6tNFEXRXmYBbh8VFdU4kprVhiB877337r777k8//dTfTZqZmVlZWaktV1RUCCE6deoky3JdXZ3L5YqOjtbaO3XqpP2Tf2OXy1VbW9upU6eUlBQhRGVlZefOnbWNtYVmdevWbcaMGWPGjNFW6+rq4uLiTlt2J1nUuD2BbBnxfD4fx6FNVFWNjY0lCAMny7LJZOL7VuC0IIyJidG7kLCh9ToE96Ms0Hf4ypUr77jjjuXLl5933nn+xkGDBq1du1ZbXrt27UUXXRQXF9e1a9fs7OwvvvhCa//8888HDx7cZOPPP/9c28xqtfbt27fxTgYNGhSUJ+aXGCVquX0CANCCgM4Id+7cOW7cuOHDhy9ZsmTJkiVCiNtvv7179+5/+tOf/vnPf/7lL3/p06fP7Nmzi4qKhBAmk+mee+654447HnjggeLiYrvdro3kTZgw4cEHH5wxY8bgwYMfeuihe+65R/uiPWvWrLvvvtvtdm/fvn3nzp3vvvtucJ+hRRYWWdi9Ir7N46EAgMgXUDjEx8c/9NBDJ/2Z2SyESE5O3rhx43PPPff1118vWLDAP2J35513ZmZmfvbZZxkZGcXFxVofaXR0dHFx8TPPPFNcXPzYY49NnDhR23jKlCmJiYnLly9PSkrauHGjdllpcCVGiRq3Gm8+zXgpAMCAJFUNm1sLCgsLm4wRJiQkBPKHfZZ63xtl6pNk9CAM/IhBY7fbGSNsE+0qBsYIA8cYYVspiuJwOOLj44O4T0O8w5OYXAYA0AJDBGFilKgmCAEAzTFIEPIDFACA5hkkCOkaBQA0zxBBmETXKACgBYYIQrpGAQAtMUgQ0jUKAGgeQQgAMDRDBGFSlKimaxQA0BxDBGFilMQZIQCgWQYJQrpGAQDNM0QQcvsEAKAlhghCbp8AALTEEEFos4h6j1CIQgDAKQwRhLIk4s2ijt+pBwCcwhBBKOgdBQC0wDhByIWjAIBmGCUIk6K5cBQA0AyjBGGihTNCAEAzDBOEjBECAJpjlCCkaxQA0CyjBCFdowCAZhkmCOkaBQA0xzhByBkhAKAZRgpCZpYBAJzCKEGYFCVVu+gaBQA0ZZQgpGsUANAsghAAYGhGCUJ+mxcA0CyjBCG3TwAAmmWUIIy3CJciPIredQAAQoxRglASwmYRtdxBAQA4mVGCUNA7CgBojqGCkAtHAQBNGSgIuXAUAHAqAwUhXaMAgFMZKgjpGgUANGWgIEyKEtUuvYsAAIQYAwUhP0ABADiVoYKQMUIAQFOGCkLGCAEATRkoCLl9AgBwKgMFIV2jAIBTGSoI6RoFADRloCBMIggBAKcwUBAmRknVdI0CAE5mqCDkjBAA0JSBgjDGJIQQTp/edQAAQomBglBwUggAOIXRgpA7KAAAJzFWEHJPPQCgCWMFIV2jAIAmjBaEdI0CAE5irCDknnoAQBPGCsJExggBACczWhDSNQoAOInBgtBC1ygA4CQGC0LGCAEAJzNWECZFM0YIADiJsYKQMUIAQBMGC0LGCAEAJzNWENI1CgBowlhBSNcoAKCJjg5Cp9NZUVFxantlZaXD4TjTj55oEbUeQRICAPwCDcIHH3ywX79+KSkpCxYs8DfOnDkzpZE+ffpo7aNHj/Y3Dh482L99UVFRZmbmhRde2K9fv/3792uN5eXlgwcPPv/88zt16vTQQw8F6Xk1zyyLaFnYPWf0QQAA4STQIOzZs+e8efN69+7tdDr9jUVFRVW/uuqqq0aNGqW119XVvfzyy1r7+vXrtcY9e/Y88MADJSUlhw4duuyyy+655x6t/f777+/Ro0dZWdkPP/zwv//7vyUlJcF7ds2gdxQA0Jg5wO2mTp0qhHjkkUea/deampr333+/uLi4lT288cYbV155Za9evYQQM2fO7NWrV11dXXx8/OLFi9esWSOEyM7OnjBhwqJFiy655JK2PYm2SIwSNR6RdeYeAAAQVoIzRrh48eK8vLyLLrrI3zJlypTY2NgLLrjg448/1lpKS0u1FBRC9OjRQ5KkX375pby8vL6+3t/eq1evffv2tfQoXq/3yJEjpb86cuTI7yg1KUpUu37H3wEAIlOgZ4StW7BgwfTp0/2rTz75ZJ8+fUwm06uvvnrttdd+//33OTk5NTU1cXFx/m3i4+Orq6vNZrMQwt+uNbb0KHv37p0zZ05sbKy2mpiY2Po5aLOsJsuRGnd9vNLWP4wAdrtdkiS9qwgnDQ0NPp9Plo11cXV7OJ1Ok8lksVj0LiRseDwen8/n9Xr1LiRsKIridDpVNdARrpiYGC1oWhGEIPz+++9/+OGHyZMn+1v8p4a33nrrggUL1qxZk5OTk56e7g85RVFqa2szMjJSUlKEEDU1NampqUKI6urq9PT0lh6od+/eM2bMGDNmjLZaV1dntVrbWm1KrM9tslitRvxoU1X1dxwxI5MkKTY2liAMnNlsJgjbRAvCmJgYvQsJG4qimEym+Pj4IO4zCO/wl156afz48Wlpac3+a11dnfb/+Pzzz9+0aZPWuGXLloSEhOzs7OTk5OzsbH/7pk2bzj///PaX1IokfpIQANBIoGeEJSUlBw4cKC8v//bbb5csWTJw4MDs7GwhhNvtXrx48ZtvvunfsqKi4rXXXhs6dKjZbF64cGFFRcUVV1whhJg6deqDDz740ksvDRky5K9//estt9wSHR0thPjzn/88Z86cLl26bN++/ZNPPpk3b94ZeJq/4QcoAACNtSEIi4uL+/TpU1tbu2TJki5dumhBuHPnzjFjxowYMcK/pcVi2bZt21tvvaUoSt++fb/88suzzjpLCJGSkrJy5coHH3zw6aefvvzyy+fOnattP2vWLJfLdf3116empr733nvdunUL9nM8CbdPAAAakwIfctRdYWFhkzHChISEtu7kmR3Kjmr1mUGmYFcXBn7fETMyu93OGGGbcLFMWzFG2FaKojgcjpAbIwwv3D4BAGjMcEFI1ygAoDEDBqGoYa5RAMCvDBmEXDUKAPiV4YKQMUIAQGOGC0LGCAEAjRkuCBMsosEnfEQhAEAIYcAglCURbxZ1XC8DABBCGDAIBb2jAIBGjBiEzLsNAPAzYhByBwUAwM+wQUjXKABACKMGocQZIQBAY8QgZIwQAOBnxCBkjBAA4GfMIOT2CQDACcYMQs4IAQAnGDEIGSMEAPgZMQjpGgUA+BkzCOkaBQCcYMQgpGsUAOBnxCBkZhkAgJ8xg5CZZQAAJxgxCOPNwqMIt6J3HQCAEGDEIBRC2CyilpNCAIBhg5A7KAAAGoMGYVI0F44CAIQwbBAmWriVEAAghHGDkK5RAIAQwrBBSNcoAEBj0CBkljUAgMaoQcgYIQBACGHYIOwcLx20M0YIADBqEObZpL21BCEAwLhBKPbW6l0EACAEGDQIu8RLx1yqw6t3HQAAvRk0CGVJdLdKP9XROwoARmfQIBRC5NoEw4QAAOMGYZ5N+olhQgAwPOMGYa5N+okzQgAwPOMGIXdQAACEsYOQOygAAAYOwu5W6XCD6vLpXQcAQFfGDUKzLLrESz/X0zsKAIZm3CAU9I4CAAwfhNLeGs4IAcDQDB2EuTYmlwEAozN0EHIHBQDA4EHIGCEAGJ2hg7BHgnSgXvUoetcBANCPoYMw2iQ6xUkHuIMCAAzM0EEo6B0FAMMzehDmJnC9DAAYmuGDkN+gAABjM3oQ5tnET3V6FwEA0A9ByOQyAGBoBKG0r15ViEIAMCqjB2GsWaRESwftJCEAGJTRg1BwBwUAGBtByIyjAGBoBCF3UACAoRGEdI0CgKERhHSNAoChEYSip036qVYlCQHAmAINwiVLltx+++0FBQUffPCBv3Hjxo0FjXz99ddae319/W233Zabmztw4MC1a9f6t1+9evWll16am5s7Y8YMu92uNbrd7rvuuisvLy8/P7/xzjuM1SKsFlHWQBQCgBGZA9xuy5Yt3bt3Ly4u/uWXX/yNFRUVZWVlRUVF2mpOTo62cO+99x46dGjdunVffvnluHHj9u7dm5aWVlFRMX78+BdffHHIkCHTp0+/7777nnrqKSHEI4888vXXX69evXrHjh2TJk369ttvc3Nzg/ocTy/PJu2tFZ3jOvhhAQD6CzQIH3vsMSHEqlWrmrQnJSWNGjWqcYvD4Xj11Ve//PLLrKysyZMnL1y4cNGiRTNnznz99dcvueSSSZMmCSEefvjh4cOH//d//3d0dPQLL7zw+uuvd+vWrVu3bmPGjFm4cOHDDz8cjKfWBrk2aW+tOjRT6uDHBQDorr1jhLt37x44cODo0aNfe+01VVWFEPv373c6nX379tU26Nev344dO4QQ27dvz8/P1xr79u1rt9t/+eWXqqqqI0eO9O/fv8nGHYw7KADAsAI9I2xWbm7uM888k5ubu2vXrnvuuaeqqmrmzJnHjh2zWq2SdOLsKikpadeuXUKIysrKvLw8rVGWZZvNVlFRoWVnQkKC1p6YmFheXt7Sw+3YsWPy5MkWi0VbTU1N3bp1a3vq98uOkj86ZKqvdwRlbyHLbrf7/78gEA0NDT6fT5a5pixQTqfTZDL536Q4LY/H4/P5vF6v3oWEDUVRnE6nGvAFjjExMWbzaZKuXUF4zjnnnHPOOUKI/v37u93uoqKimTNnJiUl1dfXq6qqfebW1tampqYKIZKTk+vr67U/VFW1vr4+JSUlKSlJCFFfX68t1NXVaRs3q3fv3g8//PAVV1yhrbpcLqvV2p76/c5LV5/d47Nao4Oyt5ClqmqwjphBSJIUGxtLEAbObDYThG2iBWFMTIzehYQNRVFMJlN8fHwQ9xm0d3hycrLD4RBCdOnSRZKkn376SWvftWuXdhFNbm7uzp07tca9e/dKkpSdnZ2Wlmaz2bRTxsYbN1+rLFut1uRfBfFA9LRJe/gxJgAwpECD0G63Hz9+3OPxOByO48ePu91uIURxcfGxY8eEEAcOHHjkkUe0czWbzTZu3Lh58+YpirJ169ZVq1Zdf/31QoipU6euWrVq69atiqLMmzfv2muvtVqtsixPmzbt8ccf93q9e/fufffdd2+44YYz9mRblBwtzLKocHb8IwMAdBZoEM6ZMyc/P//QoUPPPfdcfn7+6tWrhRCffPJJ9+7drVZr3759L7rookcffVTb+H/+53927dqVkpJSUFDw1FNPde/eXQjRo0ePoqKiUaNGpaSk/Pjjj/Pnz9c2njt3rt1uT0tLu/jii+fMmdOvX7/gP8sAML8MABiTFPiQY0scDkdsbGw7251OZ3R0dOuXchQWFs6YMWPMmDHaal1dnf8qm/a7fq3vyi7StLxIHg0K7hEzArvdzhhhm3CxTFsxRthWiqI4HI7gjhG262IZTbNp19Z23V8HeTbBHRQAYEB81T0h1ybxGxQAYEAE4QmMEQKAMRGEJ+TZpL3cQQEAxkMQnpAeK7yqOO7Suw4AQMciCH+TS+8oABgPQfgbhgkBwIAIwt/kJgguHAUAoyEIf8OPMQGAARGEv8mzST/VEYQAYCwE4W/ybII7KADAaAjC33SOl+o8os6jdx0AgA5EEP5GEiKHYUIAMBiC8CTcQQEARkMQniTPxh0UAGAsBOFJuIMCAIyGIDwJXaMAYDQE4UnoGgUAoyEIT9IlXjrmUh1evesAAHQUgvAksiS6W5lfBgAMhCBsimFCADAUgrAphgkBwFAIwqa4gwIADIUgbIogBABDIQibomsUAAyFIGyqu1U63KC6fHrXAQDoEARhU2ZZdImXfq6ndxQADIEgbAa9owBgHARhM/JsEj9VDwAGQRA2I9fG5DIAYBQEYTPOTpR2VROEAGAIBGEzBqZLJeVcOAoAhkAQNiM5WvROkr4q56QQACIfQdi8y7OlVQcVvasAAJxxBGHzCrLkTw9xRggAkY8gbN6gdOmnWrXCqXcdAIAzjCBsnlkWl2XKaw7TOwoAEY4gbFFBlrSK3lEAiHQEYYsKsqRPDhKEABDhCMIW9U6SzLLgznoAiGwEYWtGdpa4dhQAIhtB2JqCLGnVIa6XAYBIRhC2ZlSW/EUZc60BQCQjCFuTGi16JUobmWsNACIXQXgal2fTOwoAkYwgPI2CLJm7CQEgghGEpzEoQ9pVrVa59K4DAHBmEISnESWLIZnSauZaA4AIRRCeHr2jABDBCMLTuzxb+pS51gAgQhGEp3dOkqSo4scashAAIhBBGJBR/BIFAEQogjAg/CQTAEQqgjAgBVny52WKh0tHASDiEIQBSYsRuTaphLnWACDiEISB4pcoACAiEYSBKsiS+W1CAIg8BGGgLsuUdhxnrjUAiDQEYaCiZDEoQ/q8jN5RAIgoBGEbMNcaAEQegrANLs+WPmauNQCILARhG5ybLLl94qdashAAIgdB2AaSEKOyJK4dBYBIQhC2DXOtAUCEMQe43Y4dOzZu3Lhnz56xY8cOHDhQa9yzZ88777yzc+fOhISESZMmDR8+XGt//vnnf/75Z205IyPjrrvu0pb3798/f/78ioqKUaNG3XzzzZIkae1vvvnm8uXLk5KSZs6cefbZZwfpqZ0RBVnyf37l8SgmC18hACAiBPpxPmvWrBUrVrzxxhvffPONv3H+/Pnl5eVXXXVVXl5eYWHhsmXLtPbFixeXlZUlJycnJyfbbDat0eFwXHbZZRaLZeLEiY8//vj8+fO19oULF957773jxo1LT08fMmTIsWPHgvfsgi8jVnS3SpsqOCkEgAgR6BnhihUrhBAjR45s3Pjcc8/5z+oqKirefPPNa665RlsdP3781Vdf3Xjjt99+Oz09Xcu/pKSkG2644a677jKZTPPnz583b97EiRMnTpy4cePGV1555e67727nszqjLs+WVh1SB2VIehcCAAiCdnXw+VNQCHH48OGMjAz/6htvvHHnnXf+61//crlOzMXy9ddfDxs2TFseMmTI4cOHDx48WFdXt2PHjqFDh2rtQ4cOLSkpaU9JHaAgS2bSUQCIGIGeEbZu9erVy5cv/+6777TVYcOGJSUlmc3mBQsWvPjii+vXr4+KiiorK8vPz9c2iIqKstlsZWVlHo9HCJGamqq1p6WlHTlypKVHKS0tnT179hNPPKGtJiYmvv7660Gpv00usoofa6K2ltnzEsKsg9Rutzf+7oLTamho8Pl8ssyAcKCcTqfJZLJYLHoXEjY8Ho/P5/N6vXoXEjYURXE6naoa6MdvTEyM2XyapAtCEG7atGnKlCnvvPNO165dtZa5c+dqC7fddlvPnj0//PDDcePGxcXFud1u/1+5XK64uLjY2FhtWXvnOJ3OuLi4lh4oIyNjxIgRl1xyib+llY3PnDgh/vMc9ek98r+GhNnno8/n0+WIhS9VVWNjYwnCwMmyTBC2iRaEMTExehcSNhRFkSQp8I+yQN6/7Q3CrVu3jh079qWXXiooKDj1X2NjY/Py8srKyoQQ2dnZ+/fv19qPHj3qcrmys7NtNpvFYtm/f/+5554rhNi/f392dnZLjxUfH3/hhRf6H6iurk6vT6j/c57Ie8fzQIPUzRpOJ1iyLPOZ3ibyr/QuJGxwxNpKlmVVVTlibRL011i79vX9999feeWVRUVFhYWF/kan01ldXe3fYNOmTVqP6IQJEz788MPKykohxKuvvvqHP/whJSXFbDZfffXVr776qhCirq7uvffemzBhQntK6hiJUeKWXvIT3zNSCADhTw3MzJkzc3JyYmNj09LScnJyPvroI1VVx4wZExsbm/OrCRMmqKq6f//++Pj4Sy+9dNCgQVarde7cuf6d3H777dnZ2SNHjszMzNy8ebPWuHPnzuzs7OHDh/fo0WPy5Mk+n6+lGsaMGbN8+XL/am1tbYDFnwlHGtTk19xlDTqW0Gb6HrFwVF9f38oLEqdyOBxut1vvKsKJ2+12OBx6VxFOfD5ffX19cPcpqYENOZaXl9fX1/tXMzIy4uPjjxw50tDQ4G+MiYnp3LmzEKKmpmbHjh1CiN69eycnJzfez969e8vKyvr16xcfH+9vdDqdmzdvTk1N7dOnTys1FBYWzpgxY8yYMdpqXV1dQkJCIMWfIf/1lS/BIh4ZYNKxhjbR/YiFHbvdzhhhm3CxTFsxRthWiqI4HI7GCdJ+gY4Rpqenp6enN2nMzMxsduPExET/7DNN5OXl5eXlNWmMiYkZMmRIgJWEjll95X7ve2f1NSVF6V0KAOD34qvu79clXhrTVX5mByOFABDGCMJ2mfG9ij4AAA97SURBVH2h/NR2X71H7zoAAL8XQdgueTZpWKb80m5OCgEgXBGE7TX7Inn+94qbKASA8EQQtlffFOmCFPHaHpIQAMISQRgE919kenSr4iUKASAMEYRBcGm61MUqluwjCQEg/BCEwXFfX9M/vlWUMPs5CgAAQRgkV2RL8Rax4gAnhQAQZgjCoLm3r/zwVoIQAMIMQRg047rLDV6x+jDdowAQTgjCoJGEmNVXfnSrT+9CAABtQBAG05RceX+9+OooJ4UAEDYIwmAySeL/ni//cxsjhQAQNgjCILvpbPmbSnUNI4UAECYIwiCLMYnXhpuuW+vdV0cWAkAYIAiD7w+dpFkXmMZ/5nN49S4FAHA6BOEZcff58gUp0oxiriAFgFBHEJ4pzw82ba9Wn9vJhTMAENIIwjMl1izeHWn6f9/4vjzCYCEAhC6C8AzqniC9Ntx83VrfITtZCAAhiiA8swqypNt7yxNX+/gJewAITQThGTf7IjkrXvq/G7lwBgBCEUF4xklCLBxqWntYXfgjZ4UAEHIIwo6QYBHvFZju2+TbVMFgIQCEFoKwg/RKlF4cYpq42lfh1LsUAEAjBGHHubqbPCVXmrzG66WLFABCBkHYoebmm2JM4o9rfLUevUsBAAghCMIOZpLEvwvMZyeKAcu8248zXggA+iMIO5pZFo8NMD1wkTziI+/SfXSSAoDOzHoXYFDX58nnJkvXfubbXKk+nG8ySXoXBABGxRmhbi5Mlb6+xvxtpTrqI2+5Q+9qAMCoCEI9pUaLj640D+skXfxvL7cYAoAuCEKdmSTx936mpwbKhZ96mXoGADoeY4QhYWw3uVeSNH6Vr6RcfXqQKYrvJwDQUfjEDRW9EqWvxpqPOkT+Mu+y/Qr9pADQMQjCEJIYJZYVmJ64xPTIVqXve94l+4hDADjjCMKQMypL+vpq8xOXmB77Thn4gXf5AQYOAeAMIghD1KgsafM15tkXyg9sUQYRhwBwxhCEoUsSorCrvOUa813ny38pUYYs9645TF8pAAQZV42GOlkSE3vI47rJb/ykTP/S19Uqrs+Tx3WXU6P1rgwAIgJnhOHBLIs/9ZR3TTT/17nyZ4fU3Lc9V37sXbBbqXLpXRkAhDnOCMOJRRbju8vjuwunz7TqkLKkVL2nxHNusjSxh3xdrpweq3d9ABCGCMKwFGMShV3lwq7C7jV9eEB5Z5/64DeewRnStT3kYZlSro05vAEgUARheIs3iz/myH/MEXavafl+Zdl+9cEtikdRL02XB6ZLgzKk/mlSHP+TAaBlfEZGiHizmJwrT84VQoiDdvWro+qGcnXW18q2KvXcZGlgunRpunRBvHRugt6FAkCIIQgjUHa89Mcc6Y85Qgjh9IktleqGcnXJPvXuo9H1Xs/ZidLZiVLvJKlXotCW43kVADAwPgIjXIxJDM6QBmdI4nxRV9fgi07YXa3urlF316hL94kfa5Qfa9SzYqSzE8XZiVKOTeoUK7Ljpcw40SWePlUAhsBHnbEkRYlL0qVL0n+7mkZRxYF69cda8WON+nOduvWY+KVeOeIQB+2qRRZZcVLneNE5TsqKExmxUnK0SImWUqJFyq8LZm7AARDmCEKjkyXRPUHqniAuz2p6rWmNWxxqUA/ZRVmDetAuSuvUqkpR5VKqXKLKJapcapVLWM0iJUZKjRbJ0cJmkeLNIt4iEqNEgrZsFknRwmqW4szCahFRsojX/muWok2CM04AoYCPIrQoMUokRknnJAkhWrwfo9otqlzqMaeocYsat2r3CrtX1LpFrUetcIoGr6h2iXqv0uAV9R7hVoTdI1yKaPCqLp9o8IoYk4g1ixiTiDVJQghblDBJQhIiKVoIISyysJolIUSUSfgHMi2ysFpOLJslkWD5rTaTLGyWk8rTdn4qq1mytHwu6w9ph0OKjlZl+TQz2yV31Cw/JknYLCF9b4zLLcmysDAzbsA8HqEoUnRI/18NLYoiZEXEB3WfBCHaJSlKJEVJOScuRm3zu9npEw6vcPhUp08IIWrcQlGFoooatxBCuBVh96pCCLdP2L0n/sSjiHrPiWWvKo67f0spnyL2eBrvXji8QttzE/VexdPyh7UW0kIIRTHLsipEc7to5HhHze/jU0WtJ6Tnm1VVSZKEEJ7Tbgk/VZUkiSPWBjflmuYPCuYOCULoKcYkYkwiubUE1fOrst1uj42NlWUGQgPldDpNJpPFYjn9phBCCOHxeHw+JSYmRu9CwoaiKA6HO7j75B0OADA0ghAAYGgEIQDA0MI4CIuKitzuIPcUR7YFCxZUVlbqXUU4Wbp06d69e/WuIpx89tlnmzZt0ruKcLJly5ZVq1bpXUU4KS0tXbJkSXD3GcZB+Mwzz1RVVeldRThZtGhRaWmp3lWEk2XLln333Xd6VxFOVq1atX79er2rCCfr16//9NNP9a4inHz33Xfvv/9+cPcZxkEIAED7EYQAAEMjCAEAhiapakhPVNFY7969a2trY2NjtdUDBw5kZ2dzs3PgDh8+nJaWFhUVpXchYaO8vNxqtcbFxeldSNioqqoymUyJiYl6FxI2amtrPR5Pamqq3oWEjYaGhvr6+vT09AC3nzJlyty5c1vfJpyC8NixY9XV1doMTkIIl8sVHd1RkzxGBI5YW7ndbovF4n/J4bS8Xq8kSSZTcxO8ojk+n09VVbOZSb4Cpaqqx+MJ/At9p06d/KdPLQmnIAQAIOjoVwQAGBpBCAAwNIIQAGBoBCEAwNDC8lIlVVXffPPNzZs35+Tk3HrrrfyUV7Oqqqq2bNmyb9++Sy+99IILLvC3l5eXL1y4sKKiYvTo0SNGjNCxwlCzf//+FStWlJaWZmZmTps2LTMzU2tXFGXRokVbt27t2bPnLbfcwv0nfqWlpStWrPj5559tNtvo0aMHDBjg/6d33313/fr1Xbt2nT59enx8cH9OPBLs3bt3zZo1o0ePzsrK0lq+/fbbt956Kzo6+sYbb8zJydG3vNBRU1Pz9ttv+1cHDhx4/vnna8slJSVLly61Wq033XRT165d2/MoYXlG+Le//e3RRx/t2bPnBx98MH78eL3LCVHXXHPN3/72t7///e+Np/RtaGgYOHDg7t27u3XrNmXKlMavMEyaNGnLli1dunTZvn37Oeecs2/fPq39rrvuKioq6tmz5zvvvHP99dfrW2RIWbduXWlpaY8ePdxu98iRI5cuXaq1P/bYY/fdd19eXt7atWuvuOIKLk1vwuv1Tps2bebMmTt37tRaSkpKhg0blpyc7Ha7BwwYcODAAX0rDB1Hjx6dOXNm6a9qamq0du2llZmZWV1dPWDAgPLy8nY9jBpuampqrFbrDz/8oKqq3W5PTEz85ptv9C4qFGn3J40ePfrxxx/3Ny5YsGDAgAGKoqiqunjx4gsuuEC3+kKPw+HwL48YMeLhhx9WVbWioiImJqa0tFRV1Zqamri4uF27dulWYgh74IEHxo8fr6qq0+lMS0tbv369qqput7tz585r167VubgQ88gjj8yePTsrK2vVqlVay4QJEx544AFteerUqffee69+1YWW3bt3p6Wlndp+xRVXzJs3T1u++uqrtXfr7xZ+Z4Rbtmyx2WznnnuuECIuLm7w4MFffPGF3kWFombn3Fm3bl1BQYF2h3hBQcG2bduOHz/e4aWFqMZ97C6XKyEhQQhRUlLSpUuXHj16CCFsNtvFF1+8bt063UoMVQ6HY/PmzVqf1fbt210u18CBA4UQFotl+PDhvEMb271791tvvTVnzpzGjevWrbv88su15YKCAo5YYy6Xq6io6Nlnn921a5e/8YsvvgjiEQu/IDxy5MhZZ53lX83IyDh8+LCO9YSXsrIy/9FLTU01m81lZWX6lhSC3nrrrX379k2bNk3wejudDRs29OjRIyUlRVEU7cP9yJEjaWlp/ul4OGKNKYoyffr0J554ovG3Lo/HU1lZ6X+Zpaen8670M5vNw4YNq6ys3LRpU35+/uLFi4UQx48fdzqdQTxi4ReEZrPZ5/P5Vz0eD9OGBc5sNnu9Xm3Z5/P5fD4u/WiiuLj4zjvvfPvtt5OSkkRzrzeOWGMDBgzYsmVLcXFxTU3NvffeK05+jQneoSd78skn+/TpM3LkyMaNJpNJlmX/QfN6vbzG/HJycpYvX/6Pf/zj5ZdffvbZZ++55x4hhMViEUIE8YiFXxB27ty5rKxMURRt9dChQ506ddK3pDCSlZXl/3quLXD0GisuLr722muXLFkyZMgQraVz586HDh3yb3Do0KHOnTvrVF0oMpvNKSkp/fv3v//++7Vrrzp37lxRUeF2u7UNeIc2tnjx4nXr1uXn5+fn55eXl99xxx0vvviiLMuZmZn+lxmvsZYMHjy4rKzM4XBYrVabzdb4iLXzNRZ+QXjxxRdbLJa1a9cKIQ4fPrxx48bRo0frXVTYKCwsXL58udPpFEIsXbp0xIgRXNrut3HjxnHjxr3yyivDhg3zNw4ZMsRut2/YsEEIsW/fvm3btl155ZX61RhaGhoa/MubN2/WLmE/55xzsrKyVqxYIYSoqqpas2ZNYWGhbiWGmNdee23RokUvvPDCCy+8kJycfNddd40ZM0YIUVhYuGTJEiGEqqpLly7liPlpH1aa5cuX5+bmajNojx07VjtiiqK89957Y8eObdfDtOdKG728/PLL6enpN910U05Ozt133613OSHqoYce6t+/f2JiYnZ2dv/+/T/66CNVVb1eb0FBQf/+/adNm5aamvrVV1/pXWYIyc3N1U5uNP/4xz+09meffTYjI+Pmm2/u1q2b/9I+qKp6+eWXjxgxYtq0aUOGDMnIyNiwYYPWvnTp0rS0tBtvvLFXr1633nqrvkWGrMZXjZaWlnbu3Pnaa68dOXLkeeedV11drW9toeP+++/Pz8+fOnXq8OHDU1JSPv30U619586d6enpkyZNGjp0aH5+vt1ub8+jhOuvT+zevfubb77Jzc29+OKL9a4lRO3fv7+ystK/ql3RIITw+Xyff/55ZWXlsGHD/PeMQwjx/fff+zv0hBBpaWndunXTlnfs2PHdd9+dffbZ/fv316m6UGS320tKSo4ePZqenj5w4MDGP9xYWlpaUlLStWvXwYMH61hhKNu2bVv37t1tNpu2WlNT89lnn0VHR48aNYpJQvwcDsemTZsOHTqUmpp68cUXayP3mqqqqtWrVyckJIwYMaKdY4ThGoQAAARF+I0RAgAQRAQhAMDQCEIAgKERhAAAQyMIAQCGRhACAAyNIAQAGBpBCAAwNIIQAGBoBCEAwNAIQgCAof1/ye1XNulvviUAAAAASUVORK5CYII=", + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n" + ], + "text/html": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plot(results.free_energy[1:end])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Free energy is converging to a stable value, indicating that the inference algorithm is converging. Let's visualize the posterior distribution and how it compares to the true parameters." + ] + }, + { + "cell_type": "code", + "execution_count": 83, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdeXxTZaI38HOy72nTNE3SPd2hpS2llL1QoCwDyqDgMsogOio4jHrVV1yur3ccF0bH8c5cN1DBBdHr4DggvkrBQQoUWcrSFrrQfUnXNG32/f0jTqdSKC1Ne3KS3/cPPu1Jmvya0vx6znnO85Aej4cAAAAIVgyqAwAAAFAJRQgAAEENRQgAAEENRQgAAEENRQgAAEENRQgAAEENRQgAAEENRQgAAEENRQgAAEENRQgAAEHNX4rQ4/E89thjVKegmMPhoDoCXIfT6aQ6AgwHPyD/54dvdKSfzDXqcDiEQqHdbqc6CJUMBoNYLKY6BVyTx+Mxm81CoZDqIHBNJpOJz+czGP7yJz4M5YdvdPjvAgAAQQ1FCAAAQQ1FCAAAQQ1FCAAAQQ1FCAAAQQ1FCAAAQQ1FCAAAQQ1FCAAAQY1FdQAA6rlcLpvN5p3w4lr/ulwuBoNhtVolEglBEGw2myAIFotFkiSTyWQwGCwWi8vlslj4nQKgGfzSQlBwu90Wi8Vms3n/NZutJpPNZPrpX4fDw2RymUy2202wWGyXi2Cxfvaxy0UwmUyXy26329lsPZNJOJ0OBoNwuZwk6XG7nQThcbudLpeVxSL4fK5AwBUKuUIhTyDgcP+Fx+Px+XyqXwkAuBKKEAKQx+MxmUwGg6G/36DTGbq7Df39FpLkMZlcJpNPkhwWi8fhSLlcXkgINyKCx2SO6BfB4/FYrdbhy8ztdtntNrvdZjLZdDqr02l3u/vdbpvbbXM4zEymMzRUFBIiDA0ViURCkUgkFAqxEzkB/vjHP54+fZrqFEAQBOFyuZhM5g184YoVK9atW+fzPATmGvUrfjgFH104nU6dTqfX9+l0Bp3OoNebSZLHYolZLDGfLxYKxTyecOzzT46kCIfncjktFpPZbLRajU6nyeEwOp0moZAdEiIKCRGGhYlDQkIkEgmmyrxh15prdM6cOUuXLk1JSaEkFYzdkSNH+vv7P/zww/F4cPwpCnRlt9t1Ol1XV09rq66nx8Rmh7BYIXx+hEiUGB4u8s8uYTJZIpFUJJIO3mizWcxmY1ubqa6u3+lscrtNcrlYpQqRyUJCQkJEIhFVaQPMggULZs+eTXUKuEEWi+XQoUPj9OAoQqATq9U6UH56vZXNDuVyw8Ti9JSUEJIkqU53g7hcPpfLDw0N937qdruMxr6GBv2lS512ezWb7VAqQyIifupFLpdLbVqAwIMiBBowGAxtbdq6uvaeHiuXK+PxwiSSGIVCQt/yGwaDwZRIZBKJzPupw2E3GPSVlXqHo8Fu18tkvOjocKUyXCaT3diJFgC4AooQ/JfBYGht1dbWtun1Lh5PGRIyKTU1LCDLbxhsNkcmU8hkCuKnQUD99fXdFRW1TudpuVwUHx+hVEZIJIH5NwHAxEARgn/xeDw6na61VVtX12G1MrlclUyWrVRKr/+VQYAkyX+dYkxwuZx9fT1lZd2nTp1ls+3R0fLo6HC5XI4rNABGC0UI/sJisTQ0NF282GS38wQCVVhYnkCAcSLXxGSyZLIImSyCIAibzaLXd7e2dtlsl+RyfmKiWq1WCQQCqjMC0AOKECjm8Xh6enpqahrr67vZbJVSOUMgwDUko8Pl8iMiogki2uPxGAy9Fy5oT5w4JpNxEhJUkZFqjDu9qrvv/vXXX++bmOdisViffrpr8eLFE/N0MFooQqCMzWZrbm6pqGg0m9kSSWxSUhaDgdEfY0KSpHegjcczyWDoLS/Xnj59IjSUjUYcqrGx6bnnPszKmjMBz7V164Pt7e1XbNy5c+eRI0f0en1MTMzGjRtHdY1jaWnphx9+2NzcrFAoFi5ceMstt4zqeiGHw/HDDz8sWrRo5F8S2FCEQIGenp7q6ob6+m4ORx0enhMVhVOAPnZFI5aVtXkbMTFRHRUViaOmXkKhWCIJnYAn4nCuctFLQ0NDfn5+WFjY4cOHZ8yYUVFRoVarR/JoO3bs2Lx586OPPnrPPfd0dHS88847nZ2dDz300MjzGI3GxYsXu91ujLHyQhHChOru7r5woVqrtUul8YmJmSOc2wxu2KBGnGww9JaVaU+ePBodLU1KilYqlf457UCQeP75570frFixYu/evSdOnFi9ejVBEDqdrri4uKOjw3trXFxcYWHhwFe1t7dv2rTps88+u/nmm71b7r333q6uLu/H3377bUlJiUqluuOOO6RSKUEQxcXFYrG4ra3t2LFjWVlZa9asIQjiq6++Ighi27ZtJEmuXr1aLpc3Njbu2bPHarUWFhZOmzaNIIi2traBSARB7Nq1a/ny5aGhoUVFRVFRUeXl5WfOnHnqqadMJtOnn37a3d0dHR19yy23KJXKcX/hxgF+DWCC6HS6w4dLiorKbbaYlJR8lSoWLTiRvI0YFzc5JWWR1RpTXNz6xRcHz5690N/fT3W0oGaz2Q4cOKDX63NzcwmCuHz5clpa2ldffdXQ0LBly5bXXnvtwoULg+//9ddfq1SqgRYkCIIkSYVCQRDEM8888/jjj6tUqpMnT+bk5Hh/srt27brjjjv27t2rUqmefvrpv/71r0MzVFRUTJs2TafTCQSClStXfv755wRBVFZWvvjiiwP3efLJJ9va2giCeOedd1atWvX9998rlUqz2Txjxoy+vr7Jkyd3dXWVlZWNy2s0/vBOBONOp9NduFCl1dpCQhKTkiJxNIZaDAZDLlfJ5Sq73drc3Hrx4umwMPakSbFqtRrTf0+whx56aMeOHQ6H4+23346OjiYIYuvWrb/85S/feecdgiCmT5/++OOPP/7444O/pLGxUaPRDH2onp6e119/vaKiQqPRPPjgg4sXL96+fftjjz1GEMTkyZO9DyiXy7dt27Z58+ZVq1Zt2LDh/vvv9/4yPvLII/fee+8f/vAHgiBiYmKefPLJ2267bZjYM2fOfPvttwmCqK6uttlsTz/9NN0v2sEeIYyjrq6uf/7z+HfflVmtsUlJ+QpFFFrQf3A4vMjIhJSUBUxmSklJ5//+76EzZy709vZSnSuIvPnmm2az+fjx408//bR3Is2KiooFCxZ4by0oKKitrdXpdIO/hM/nGwyGoQ9VXV0tl8sHOnLOnDnl5eXej7OysrwfREVFDRxEHayiomJgFtbZs2fX1taazeZhYnt3XgmCSExMnDt3rkqluvPOO//3f//XT5ZwuAEoQhgXBoPh8OGSoqKLNltcUtK88HA1KtA/kSQpkykSE6fFx89vbRV+8825774rbmtro++bGu3k5uYuXrzYW4QCgcBisXi3WywWkiSv2NnKzc29ePHiFe14xRd6v3ZgSNR1p+Lj8/mDn5TFYnE4HBaL5XQ6B+4zuBoHJrxlMBh/+9vfzp8/P2vWrGeffXbwoVR6QRGCj7lcrsrK6v37TxgMyuRkVCBtcDjcyMiE1NQFLFbasWOtX355sKqq2uFwUJ0rMFmt1oGds56enpKSkkmTJhEEMX/+/N27d7tcLoIgPvrooxkzZlxRhAsXLkxLS9u4caPJZPJuOXXq1LfffpuSksJms7/55huCIIxG4549e+bPn3+tZxeLxUwms6enx/tpfn7+rl273G6390nnzJnDYrFiYmLq6+v1ej1BEIcOHbrqoQKDwWC1WmNjY3/729/+5je/uXjx4lhfF4rglAD4UkdHx8mTFVarLDZ23lWHjIP/CwmRh4TIjca+Cxfqz5z5PiMjKjFRQ/eTQP6mt7d38uTJcXFxAoGgrKxs7dq1d955J0EQjz76aHFxcVZWVkRERG1trXd452AMBmPv3r333HNPVFRUenp6Z2cnSZI7duzg8XgffPDBvffem5mZWVlZWVBQcOutt17r2Vks1qZNm9LT06Oioj7++OMnn3xy1apVOTk5ISEhLS0t//jHPwiCiIuLW716dXZ2dlJSklgsjoiIGPo4lZWVK1asyMnJYbFY58+f946yoSMszOtHaL0wr8lkKi0tb2y0qFTpISFyquOMi7EvzEs7druts7PRaGyIjw+dNCkpJCSE6kTXMczCvFu3bh28HuG8eQuEQpVSGTMBqY4e/fq55568++67B2+02WxVVVUWiyUhIUEu//evjMfjqa+v7+/vnzRpEofDudZjtra2NjY2KpXK2NjYgYOfJpOpsrJSqVRGRkZ6t5jNZgaDwePxCIJwOp1ms1kikXhvstvtJpNJIpEwmUyPx1NXV2cymdLS0ths9kCSyspKh8ORkZHR398vEomYTKbRaGSz2QNHRw0GQ1VVFYPBSEtLG9dfjY8++ujQoUNYmBf8lMvlqqqqOX++SSRKSk2Nw4HQQMLhcKOikl0ujVbbuHfv6eho0eTJGu9gfbp75pkt586dm5jnSk+/Z+HChVds5HK5U6ZMGXpnkiSvOi70CpGRkQNtN0AoFObk5AzeMnjyBBaLNdCCBEFwOJyBoiVJMiEhYWiStLQ078feqxIJgrhifiKxWOy97pDWUIQwJnq9vri41GoNjYvLx7HQQMVksqKiEiIjNV1dbUVFleHhVVlZKXSvwyVLlixZsoTqFOAXMFgGblxtbf3+/SfZ7FSNJhstGPBIklQoIlNT55FkclFR1cGDxwZmPwl4PT09zz33f5OTs4XC0JCQiAULln/55ZdUhwKfQRHCjbBarT/8UPLjj+0azbzw8BFNkAgBQyaLSEmZ43Jpioouff/98YHBh4GqtLQ0OTnj1VeLa2qeMJsP9vXt+eGH2XfeuWnlyltveFiDTqd74IEH3njjjYEtJpPpwQcffOWVV3yU+kaUlZW1tLQMfCoSicZ4XemRI0cGRrf6MxQhjFp7e/v+/cW9vWHJyTM4HB7VcYACJEnK5arU1HyHI/7//b8Lhw+f8I6zDzx9fX2LF6/Q6x+0Wr8niDsJIocg5ng8z9hsZQcPXn7ssSdv7GGNRuO2bdtefvnlgddtz549u3fvHjpMdCK99NJLe/fuHfj0L3/5i1AoHMsDrl27tqmpacy5xh3OEcIouFyuiorK8+c7oqJyJBIZ1XGAYt46DAtTdnW17t9fGh0tyMxMGxhVERjefPMtiyXB7X5uyC3hVusnb7899dlnt1z10oKRuOmmmz777LMHH3yQIIgdO3bcdtttg2cW/eqrr0pKSsLCwu65557w8HCCIOrr67/66qumpia1Wr1+/XrvxjNnzrS1tUml0q+++io6Onrjxo3eMaKDOZ3OXbt2lZWVxcfHb9iwwTu8s6am5tNPP+3p6YmOjr777rtbWlouXbpktVoNBkNmZubSpUsNBoP3soLPPvssIyOjuLi4srJy5cqVCxcu/Pvf/15cXDx9+vTbb7+dIAiXy7V3795Tp065XK758+cvW7aMIIgvv/zSZDK99957CoVi5cqVkyZNMplMH374YW1t7ZQpU+66667rXuk/YbBHCCNlMBi++664qsqRkjIPLQgDSJJUKKKSk+f39UV8/fWpH38sHX6CLnr58stvLZZfXePGdB4v4+DBgzf84L/+9a+91wM0NDQ0NjbOmzdv4KZNmzb99a9/zcjIMJlM06dP9+44fvPNN06nc+bMmZ2dnbm5ud6jjsePH9+8efO7776bmZm5b98+b60O5na7Fy9efPDgwalTp5aVlRUUFLjdbp1ON2/ePIlEMm/ePIfDUV9fz+VyORyOQCAIDQ31jjV99NFHvU+xffv2tWvX9vT0aDSa1atXP/DAA4cPH05PT9+yZcvu3bsJgujr6/v6669TU1NTUlKeeOKJt956iyAIoVDIYDAkEkloaCiHwzEajbm5uTU1NTk5OXv37r3iYhJqYY8QRqSrq+v7789JJGkaTRTVWcAfMRgMtTpeqYzRautra49mZcUkJycGwCzeWm0bQcRd61aXK37wSbXR8i7dcOnSpc8//3zdunUDlx5dunRpz549DQ0N3l235ubmDz/88OGHH/YuOuhwOBYvXnz69OlDhw7ddNNNBEGIxeJPPvmEJMnc3NzBber197//3WKxfPTRRyRJ3nnnnXl5eYcOHQoLC+Nyuffff//gg58JCQmzZ8++//77h0a95ZZbnnnmGYIgTp8+3dra+vXXXxME0d/f//XXX99xxx0ymez9998nCMJgMCiVyueee27Tpk1Llizh8/lr1671XoPxxz/+MTs7+89//jNBEGvWrImNja2vr4+Pj7/hV8+HaP/fFCZAY2PT0aNVajUOh8J1MBjMyMhEuz3q4sXqS5cO5+YmR0dH0/rSUoFASBDXHO7BZBquuK5utNatW7dz5849e/Z89913J06c8G48f/682WyeO3eu91OtVuttxG+//da7oASfz6+vr29tbfXeIT093fsiq9Xqnp4ej8cz+DU/d+5cbW3twEzZ9fX1NTU1BQUF06ZNU6lUhYWFq1evvv3224dfnDI9Pd37gUKhGFhAWKFQeIdKWa3W+++///Dhw+Hh4Tab7arjp86dO3fkyJGBiw77+/trampQhEADHo+noqLy3Ln2uLhZfP6YTptD8OBwePHxU4zGvuLi8vDwxtzc9NDQiVgIfjzMnJlTX3/Y5brlajdaHI4fc3L+ayyPv27duuTk5Nzc3ISEhIEiFAqFGo2mqKho4G7emVzWrVu3d+/eGTNmEASxfPly75SkxPWm1RaJRAUFBd6VmLwEAgGTyfzb3/6m1Wr37t37/PPPt7W1XbHe0xUGapIkyaFP99577+l0uvr6eiaTWVJSsmrVqqvGuOuuu5588snBW4Z5xomEc4RwTS6Xq6TkzIULfUlJc9CCMFoikTQ1dbbHk7h/f+mJE2cGr41AIxs33keSHxBE5dCbGIw/xMfHT58+fSyPr1ard+7c+ac//WnwxlmzZrW1tZ09ezY0NNR7xs5sNjudzr6+Pu/AnJqamsOHD4/wKZYuXXrw4MGenh7vozEYDIfD0dvbazabVSrVAw88sHbt2suXLxMEERoa2tnZeQPfRU9PT3h4uHeqtvfee29g++AHXLFixRdffOFyubwxbDbb8PugEwl7hHB1Vqu1uPiUXi9NTs6h9aEtoJZcrpLJFK2ttfX1xdnZccnJif7z9jcSM2fOfOSRzX/9a4HN9hZB3EwQ3t8FHYv1Bzb7/c8++2Hs387Q2bHDwsI+++yzDRs2qNVqDodTX1+/ffv2wsLCJ554Ys6cOZMnT9br9d79wpHIzMz805/+NHfu3NTUVIfD0dzc/P3337e2tq5Zs2by5MkEQbS0tHgv21i/fv1tt9328ccf33bbbS+99NLIv4V169bNmzdv/vz5er0+Ly9vYPujjz56++23CwSC119//eabby4rK0tPT588ebLRaOzu7i4vL/eTmXsx6bYf8Z9Jt/v6+v75z9NMZpxafeX0g8EsCCfd9iGr1dzScpHP78/Lm3zD1xtc18gn3R45j8fz1ltvP/3083Y7g8lMIQiz1XohOztv5863vUVyA9xud19f3xVHjO12u81mG3gTcLvdDQ0NNptNo9EMTHLd2tqq1+tTU1OtViuLxeJyuTabzel0ese8eDwevV5/1QPRDoejtraWxWLFx8d7j21aLJb6+noWi5WQkDD4aKd30K93N9Q7fNRgMHjHlHq/iiAI72+B3W632+3eI5xWq7WmpkalUoWGhhoMhsHzs/f39/N4PO+X22y2mpoasVgcExMzqr+wx3XSbRShH/GTIuzr6ztw4KRUmhEWpqQ6i39BEY6dXt+t1ZYlJ4dkZk4aeHP3ofEoQi+bzXb8+PH6+no+n5+Tk5OcnDy2pDA6WH0CJk5/f39REVoQxktIiFwimdfaWltff2TmzJSYmIlYBcknuFzuggULFixYQHUQ8D06HayH8dbf33/gwI8SCVoQxhGDwYyOTlap8oqLG4uLf6TFIJr29vbHHn8sNiGOxWYLhIKZc2d+8sknfnI4DcYOe4TwE28LikST0YIwAYRCSUrKnPb2hq++Ks7NTUhI0PjtmKwTJ04s/cUyfrwo7GZVtkrjtrm6Kns2/cdDH37y0b6v9g6dz2wkenp6HnnkkaHb//CHP8TGxo4x8B133OFyudhsdkxMzJ133pmRkTHGB7wBvb29JSUly5cvn/invgEoQiCIQS2IpSRgwpAkqVLFh4ZGnDp1obGxffr0Kf5wjvwKOp1u2Yrl8iXqyF/8e+CYKF4aMS/67Gulv3vk4W3vvHsDD8vn870lYbVaN2zY8Je//MW7SP3ghXNv2J49e1588cXY2NjS0tJp06adOHEiOzt77A87KvX19Q899FB9ff0EP++NQRECYTQai4pOCoWT0IIw8Xg8QXLyjI6Opr17S7KyotPSUvzq+oo333qTo+INbkEvlpCtuX/yB0+9//xz/3dgppWREwgEd9xxB0EQBoNhw4YNN910U2xsbHFx8dmzZ3t6eg4dOrR58+ZLly6lpaV5B6ZeunSprKxs7dq1BEHY7fYdO3ZcuHAhNjb2wQcfvGp3Ll68OCsra+3atefOnfvb3/4WHR29e/fuyspKkUi0Zs0a7/QuHR0dX3755Zw5cz744IOcnJxf/vKXu3fvvnDhAovFWr58+aJFiwiCsNlsb7zxxtq1a99++22Hw/HQQw8pFIq//OUvnZ2dd95558AlHJcvX961a1dPT09BQYH3avqPP/64r69v69atBEE8/PDDPB7v4sWLu3fvNhgMS5Ys8c7KXV9fX1xcnJKSsmvXrkWLFhUWFr733nsVFRUcDmfx4sUrVqwY7at6w/zoPxxQwmQyHThwQiicpFBEUp0FgldERIxGM6+szFhUdNRgMFAd59/27t8rzQu/6k18lUgWLz906JCvnquoqGjdunWHDh3Kz8/n8/k7d+48e/as96bz589/8MEHBEG43e5FixadOHGisLCwu7t7zpw5DodjmMe0WCwsFqu8vLyvr2/p0qVxcXG/+MUvzpw5QxBEW1vbli1b/uM//iMrKys6Orq5ubm2tragoCArK+uBBx7wrjxss9m2bNny0EMPZWZmEgRRWFh4zz33KBQKjUZTWFjoXZm5tLS0oKAgNDR0/vz5W7duffXVVwmCiIyMZLFYGo1Go9EwGIwjR44sW7ZMrVbPmTPnqaee2r59O0EQ1dXVjz/++H/9139Nnz49IiLiiSee+P7772+66abZs2dXV1f76lUdCewRBjWHw3H48EkeLxn7gkA5DoeXlJTb0dG8d2/J7NkpcXFjPVXmE+3tHaHzrvk3IkfOa2tr8+HTxcfHD54Lbai9e/fabLYdO3YQBHHzzTefPHnyu+++G7rztGfPnmPHjp04caK0tHT79u3Jycnz58+32Wxms7mtre2zzz7LyckhCMJkMn388cdK5U/DAl5++WWHw2E0Gu12+yeffLJ69eqB7ZmZmXfccceHH364cuXK9evXe5McOXJkzZo1zz///NNPP+1d9SIrK2vGjBlPPPFEQUHBm2++uWbNGu8jPPvss1u3bvWu2aTRaG6//fbf/OY3BEFYrdbPP//ce0j8qaeeuuOOO7w7ixMMRRi8PB5PSckZu11JoyHsEPAiIqKl0rDjx8+2tLTn5maNx7WGoyIUCZ2Wa+5yuS0u357X9O54DePChQt1dXUDU1c3NDR4Z0e7Qmtrq8vlysvL27p1q1qtvnTp0rp160wmU0RERHt7+8ApQ4VCMdCCHR0dv/rVrxobG9VqtcFgYLPZA4/mXT6CwWCEhYWlpqZ6N8rlcp1O541UXV09MLNaV1dXf3//0NgvvPDCa6+9RhCEy+Wqq6tzu90EQWg0moEXcMuWLb/+9a9fe+215cuXb968WaPRjOw18wEUYfA6f768tZWRlJRKdRCAn+HxBMnJs1pba/btOzJvXqZCoaAwzKy8WUWXDofnXeWQidvu0lV3Dazq4BPe6Ve8WCyW0+n0fmw0Gr0fSCSSgoKCzz//fPjH+d3vfpeVlTXw6ZYtW9atW7d582aCIF588cXy8vKhT/fKK6+kp6d7l1fcvXu3d70kr8FnbQfPvu29gEQsFr/22mtLliwZJo9EIvnggw8Gz742NEBhYWFra+upU6d27tw5e/bs5ubmCVvGC+cIg1RdXX15uS4hYarfjlmHYEaSZFRUskw29cCBsrKyCu/eAyU2Pbixo7jF1HTlLg5BEE17alJSUwd2znwuISHh2LFjBEHYbLZPP/3Uu3HZsmVFRUXnz5/3ftrd3e1ds3d4/f393onQ9Hr9xx9/PPx9bDab9xzeCK1ater1118fuB7Uu4caFham0+msVuvAfV599dWBucOuuhd7+fJlBoORl5f3+9//vqura+BrJwCKMBh1dXUdP345Pj6XycQhAfBfUmlYUtK8S5dsBw4UUzWCZtq0aVue3HJp6+nO460e109X0Nv7bHU7K3qPtn/60a7x+1Ny8+bNBw8ezMjImDp16sCMbikpKdu2bVu5cuX06dOzs7Pz8vJGsl7E//k//+fxxx+fO3fu7NmzrzVb929/+9v3339/9uzZmZmZKSkpI8/5zDPPxMbGajSauXPnJiQkPP300wRBxMbGrlmzRqPRJCQk6PX6F198kc/ne+8THx/vHU16hY0bNyYmJnoXSnzxxRcncpEmzDXqRyZmrlGj0bh/f4lCMU0ioesScVTBXKNU6exs6e29OGNGkkZznXVcx2mu0Z07dz759JN9/f0StdTtcPc29eQX5G97e1tiYuKNPeBgV6yjO5jT6WxpaYmIiBj6v8575FCpVI6wiU0mk1arjYmJGXw08go2m62lpUWpVA5etn6EvF8bHh4+zKWQFoulra1tmMfX6XQ6nU6tVnsn+x4Mc42Cz9jt9u+/PymRTEILAo0oFFEiUcjx46Wdnbpp0zIn7NTRgPXr1991112lpaV1dXXeSbejoqJ89eDDNBmLxYqLi7vqTdHR0aN6FqFQeN3a5nK5CQk3uODMSL6Wz+cPfx+ZTCaTyW4swFigCIPL6dMXnE5lZCQuGQSaEQhEKSlzGhsrdLpj8+blTPzi5pA5slsAACAASURBVCwWa/r06WNchhf8E84RBpGmpqa6OlN0NIaJAi0xGIz4+AySTNi/v6S9vZ3qOBA4UITBwmw2Hz9eGRs71a/mrwIYLYUiSqmcfvDgxfLyS34yxAHoDu+JQcHj8ZSUnBWJkgUCv5vUGGC0RCJpYuKcCxf6jh07NfwEYwAjgSIMCpWV1Z2dLJUqjuogAL7BZnOSk/O6usTffutfc5MCHaEIA19vb++ZM01xcVnXvysAfZAkGROTxman7NtX0tLSQnUcoLGRjhq1WCw//PBDY2NjVFRUYWHh4GnoBqupqTl8+DCHw5kzZ87AMFmtVltZWalWq0d1kSb4hNPpPHbsXHh4OodD8ZyNAOMhPFwtEIj++c/TGRm6zMxrrkDLZrPXrl17Y4vogm8Nc93kMIxG48qVK8cjDzHyC+pTUlLUanVaWtqpU6fsdvvRo0eHXvr95ptvPv/888uXL2cymQwGwzsH64YNG7744gsOh3PXXXf993//97UeHxfUE+NzQX1p6YXaWkKjmeLbhw1OuKDebzmdjrq60uhoMiMjRSwWDx0R1tvb29vbS0k2uILJZLqBC/YJglAoFON02cxIi7C5udl7/abT6czIyHjsscfuu+++wXeoqqqaOnXq6dOnvfOUD+js7AwLC/vd737HYrFQhMPzeRH29fXt23cyOXkBplLzCRShP/N4PI2NFUxm6+LFc27sfRYmxsRMoTUqIz1HODCLAYvFkkqlQ+fA/fLLLwsLC8Vi8b59+yorKwe2KxQKJpPpk6wwWqdPl4eFpaEFIRiQJBkXl+50xnz77fGhywABDGPUb5H79u2rr68fWLBxQF1dXVNT04oVK7Kzs7/77ruNGzf+53/+56ge2e12v/TSSwOfTp06deHChaONR2sOh8OHY8FbWlra2lzJycqBlVxgjDwej9PpxOvpz+TyKKNRuG/f8QULpoSHX31leaCWb9/orst7qm74+4yuCM+cOXPfffd9/PHHcrn8ipu8M67W1taKRKLKysqMjIz77rtPpVKN6vEHH8Tv7++ncO0VSrjdbl99y06n8+TJKqVyarC9huPN4/HgJfVnHo9HJovgcLhFRWfnzEmOisJsgn7Hh290IzGSKURGUYQXLlxYsWLFu+++W1hYOPRWtVo9ZcoU75nM1NTUkJCQ6urqURUhg8F49dVXR37/wGO32321Hnd1dS2TqZbJqFzRNPB4PB6XyzXM5P1AOZfLxWazw8IiBIK5J06czMlxpqYmUx0KfsaHb3S+MtJzhFVVVcuXL3/99ddXrVo1eHtzc7N3hEthYeHly5e9R43a29t7e3tjYmJ8HhdGwmg0lpW1RkVhTlEIXny+UKOZdeZMR2npeczEBsMb6ajRuLg4Lpc7f/5876dLlixZvXq1x+NhMBgnTpzIy8vzeDyLFi3i8/kFBQW7du3KycnZtm0bQRBff/31vn37jh49ymAwZs2adfPNNy9fvnzo42PUKOG7wVTFxT/q9Qq1+jort8FoYdSo/7NYLFwud+BomNvtqq0tVatds2ZNm/jFm+Cq/HDU6EiL8Kuvvhp8ejMlJWXKlCkEQXzyySfLli0LCwsjCMJut3/22WeNjY05OTkDbXfx4sWKioqBL0xPT7/i+govFCHho/8fHR0dBw9WpqTMG7+Fs4MWitD/XVGEBEF4PJ6GhnKpVJ+fn4fD2v6AxkU43lCEhI/+f3z3XTGDkSyTRfgkEgyGIvR/Q4vQq6mpksfrKCiY4W9np4KQHxYh5hoNKJ2dnd3dHrQgwBViYlKdzqhDh0qsVivVWcDvoAgDSnn5Zbk8ieoUAP5IrU5wu+MOHDhuNpupzgL+BUUYOLq7u7Vae1iYkuogAH5KpYpjsRK/++44Vm6CwVCEgaOi4rJMlogxMgDDiIiI4fPTiop+RBfCABRhgNDr9S0t5vBwzKMBcB3h4ZECwaRvvz2h1+upzgJ+AUUYIMrLq6XSBOwOAoxEeLg6JGRKUdEprM0EBIowMPT19TU29kdERFMdBIA2ZLKIkJDMAwdO63Q6qrMAxVCEgaCqqk4i0YxkblkAGCCTKeTy7KKiMzhGGuTw1kl7TqeztrYzPDyK6iAA9BMSIg8LyyoqOoUlDIMZipD2WlpaWCwFm425owBuRGhouFSaceDAj0ajkeosQA0UIe1VVjbLZDg7CHDjwsKUQmHqwYM/WiwWqrMABVCE9Nbf39/T45BKw6gOAkBvCkU0k6n5/vsTNpuN6iww0VCE9FZf3yQSReOqCYCxU6vj7Xb14cM/Dl5pB4IBipDG3G53VVUbhskA+EpMTIrJJD9y5KR3jXEIEihCGtNqtR5PCJeLVYEAfCY2dpJOJz5+/LTb7aY6C0wQFCGNVVc3h4RgmAyAj8XHZ2i17JMnz/rJcq0w3lCEdGW1Wtva+rH0IIDPkSSp0WTX1TnKyy9SnQUmAoqQrjo6OjgcBWaTARgPDAYjMXHa+fPd9fUNVGeBcYe3UbpqbOwQixVUpwAIWEwmKz4+9/jxyx0dHVRngfGFIqQll8vV0qILDQ2nOghAIOPxBFFRuYcPX8AEbIENRUhLXV1dTGYIi8WmOghAgBOJpCEh6YcOnbRarVRngfGCIqQlrbZTIMAwGYCJIJermMy44uJTLpeL6iwwLlCEtFRX1xkaihOEABMkMjKxry+kpOQMLqgISChC+tHr9TYbi88XUh0EIIjExaU3NbkrKiqpDgK+hyKkn/b2Dj4fx0UBJhRJkomJ086d62xoaKQ6C/gYipB+Gho6Q0JQhAATzXtBxbFj1Tqdjuos4EsoQppxuVw9PSaRKITqIADBiMcTKJXZhw+XYhBpIEER0oxer2exJJhQBoAqISFyBiOmpAQzkQYOvJ/STG9vL5uN3UEAKkVFJbW3MzFwJmCgCGmmo0MvEKAIAajknZX73DmtVqulOgv4AIqQZtrb9RJJKNUpAIIdi8WOiZlWXFxuNBqpzgJjhSKkE7PZbLF4sBIvgD8QCiVCYfLRo2cw4wzdoQjpRK/Xc7nYHQTwF0plrMEQUlp6geogMCYoQjrp6dFzOChCAD8SF5dRXW1sbGyiOgjcOBQhnWi1vWIxRsoA+BEGgxEbO/X48aq+vj6qs8ANQhHShsfj6e7uF4mkVAcBgJ/h84UyWfqxY2dxspCmUIS0YbVaPR42k8miOggAXEkuVxmNIeXll6gOAjcCRUgbJpOJxcKKEwB+Ki4uvayss7Ozk+ogMGooQtowmUxMJooQwE8xmazIyKzi4vM2m43qLDA6KELaMBrNLJaA6hQAcE0SiYzBiCotLaM6CIwOipA29HoTFuMF8HNRUSm1tebm5haqg8AooAhpA0UI4P8YDEZsbHZJySWLxUJ1FhgpFCFt9PWZeTwcGgXwdwKBmMfTnDiBdZpoA0VID1ar1e1m4doJAFpQqzVaLXn5ch3VQWBEUIT0YDKZ2GwcFwWgB5IkY2MzT5+uMxgMVGeB60MR0oPFYmEycVwUgDZ4PIFUmvLjj+dxgNT/oQjpwW63MxgcqlMAwChERMR0djIxH7f/QxHSg93uIEk21SkAYHSio9N//LHKarVSHQSGgyKkB5vNwWKhCAFoRiAQczgx589fpDoIDAdFSA9WK4oQgJaiopKqq/UdHR1UB4FrQhHSg8ViZ7NxjhCAfhgMpkqVceJEORZp8lsoQnrAoVEA+goNDbdaQysra6gOAleHIqQHFCEArcXETD5/vqm/v5/qIHAVKEJ6wDlCAFrjcLhSacrJkxdwWaEfQhHSA/YIAeguIiKmo4NsamqmOghcCUVIA06nkySZJElSHQQAbhxJktHRGadOVTmdTqqzwM+gCAEAJohQKPF4IqqqLlMdBH4GRQgAMHGiopIvXGjCaoV+BUUIADBxOBwenx9XVlZJdRD4txGtb+dwOIqKig4fPmw0GjMzM9evX8/lcofezWAwbNu2raqqKjw8/Ne//nVycjJBEDab7c0336ysrJw0adLGjRuv+oUAAD7R3a0NDQ0TiaRUBxmOWq2pqTmcktInlfp1zuAxoj3CsrKyZ599NjQ0NDs7e8eOHStWrBh6H71en5eXV1JSkpubKxQKq6urvdvXrVu3f//+/Pz8f/zjHxs2bPBldgAAgiAIwuPx7P7fN1asiX3kxWV3/Tbn7t9Mq6o6S3Woa2IyWSEhyaWlmIDUX5AjuajF6XSyWD/tO3Z0dCiVysbGxpiYmMH3efLJJysrK//xj38M3lhXVzdp0qS2tjaZTNbT0xMZGVlVVRUbGzv0KRwOh1AotNvtY/heaM9gMIjF4qHbnU7nZ58dTEtbOvGRYDCPx2O1Wvl8PtVB4ErvvP/cobo9kbfFMThMgiAsWmPztrq3thbFxCRRHe3qPB5PZeWRxYtTIyIiqM4y0a71RkehEe0RDrQgQRC9vb1MJnPoHv2BAwduvfXW999//4UXXigpKfFuLCkpyczMlMlkBEGEhYWlp6cP3AQA4BN2u23fdzuj7tR4W5AgCL5KpLhVvf3j31MbbBgkSUZEpJ05cwnX1/uDEZ0jHOByuR566KHNmzcPLcKGhobf//73t956a0hIyM033/z666/fddddWq1WLpcP3EehUGi12mEe/JZbbhn4dP78+ffee++o4tGdxWJhMplDtzudTpvNjiXNKOfdI8QFnf6mqalaqBaTzJ/9XKQpssp/lPrzb41AILl8mVFZWRUbG3P9eweQa73RjRMOhzN4X+6qRlGEbrd7w4YNLBbrlVdeGXorm81etmzZyy+/TBBERETEyy+/fNddd3G5XIfDMXAfm83G4/Gu9fgMBuO2224b+DQ5OXmYOwckh8Nx1W/Z6XSy2WwOB6tPUM/j8eAH4W9EIonLfuXCDi6ri8vl+fkPKzZ2Snn5yYQEzXXfqQPJtd7oxgmDcf0DnyN99T0ez4MPPtjU1LR///6rjvyMiopKSvrpcHxycnJbW5t3Y3Pzv+cTamlpiYyMvNZTkCS5du3aEeYJSAwG46o/MwaDwWCQI/lxwrjyeDwkiR+E31Gr4zx6t73PxpH++61JV9Ixf/YqP/9hicUhXV2KurqG1NRkqrNMnGu90VFoRGk8Hs/mzZsrKir27t0rEAgGtre3t3/33Xfej2+99dYffvjB+/Hhw4enTJlCEMTChQtbW1tLS0sJgjh9+nRHR0dBQYGPv4MgwGQyPR4XziUAXMszj7/b8Ealvrzb7XA7DPaOb5o9Zz133f441bmuT61OLitrxKRr1BrRHuHRo0fffPPNlJSUBQsWeLe89957WVlZJ06ceOCBB7wrL2/atGnBggX5+fkymezkyZN79+4lCEIikbz44ovLli3Lz88/cuTIyy+/LBKJxu+bCVQkSbJYTLfbxWQG0fETgJHLmbpg++tH3tnx3KW9pwUiydI5v7rjdw/TYp56Hk9AEPL6+oakpESqswSvEV0+YbFYrhjkolareTye0WjUarUDR0QdDkdJSYnD4Zg+ffrg0bG1tbWXLl2aNGmSRqO51lPg8gli2FHFf//7IYViJo8nuOqtMDFw+YT/s1gsXC7X3468Dc9sNnZ0lKxatSBIzhT64eUTI3rd+Xz+VTtMJBINtCBBEGw2e968eUPvlpCQkJCQcMMRgSAIHo/jdDqufz8AoBuBQOR2y5qamjWaeKqzBCk6/d0UzHg8NooQIFAplUlnz152ua4c+woTA0VID1wuihAgYAmFEocjpLm5heogQQpFSA/YIwQIbEpl0rlzl91uN9VBghGKkB6wRwgQ2MTiEItF1NraSnWQYIQipAcul+1yBfWQWoCAp1Qmnz1bgyuGJx6KkB7YbLbbjT1CgEAmkYSaTHzvtFwwkVCE9MDhcDwe7BECBDi5PKG8vI7qFEEHRUgPAoHA5bJQnQIAxpdMpujudvb29lIdJLigCOlBIBA4HCaqUwDAuJNIYmtqGqhOEVxQhPTAZrM5HNLhwNFRgACnUERfvtxps9moDhJEUIS0IZUKrFYz1SkAYHyxWGw2W9XQ0ER1kCCCIqQNFCFAkFAq4ysqGnEdxYRBEdKGRIIiBAgKAoHYZhO0t7dTHSRYoAhpQyQSOJ0oQoCgEBoad+lSA9UpggWKkDYEAoHbjSIECApyuUqrNRkMBqqDBAUUIW0IhULsEQIECZIkBYKY2toGqoMEBRQhbfD5fLfbisnpAYKEUhlbWdnmdDqpDhL4UIS0QZKkVCqwWnFZPUBQ4HC4DEYYhsxMABQhnSgUUqOxj+oUADBBJBJ1bS0WZhp3KEI6kcslZjOKECBYhIUpW1r0VquV6iABDkVIJyEhIU4nihAgWDAYDA5H2dqKhZnGF4qQTqRSqcPRj/kmAIKHTBZZU4Ojo+MLRUgnLBZLLOZYLBgvAxAspNKwnh4bLigcVyhCmomICMF4GYDgQZIkn69uacHR0XGEIqSZ8HCpxYIiBAgiYWGR1dU4OjqOUIQ0I5VK7XY91SkAYOKIRFKzmYll68cPipBmpFKp04nxMgDBhc9XNzVhp3C8oAhphs1mY7wMQLAJD4+sqdHiL+BxgiKkH5UqtL9fR3UKAJg4PJ7Abufo9TgtMi5QhPSjUsnN5m6qUwDAhOLxFO3tnVSnCEwoQvqRy+U2Ww/VKQBgQoWGRjQ2ogjHBYqQfvh8vkTCMptxgW2Qam9vqqwsNZuNVAeBCSUWh/b0WDDv6HhgUR0AbkRMjLy5uVsgEFMdBCbUxYun//M/HzEYwj2ecI/n3Ny5s5588iUeT0B1LpgIJElyOPLOzs6YmBiqswQa7BHSkkIRZrHgNGFw6ehofvjhBxobP9Lp/t7bu02v//Hbb5OfeupBqnPBxBGJIlpacHTU91CEtCSXy+32HoylDiofffSOXv80QWj+tYF0ODaVl7d2dWHyrWAhkykaG7vdbjfVQQINipCWOBxOWJjAYMBY6iBSWVnp8Uy7YqPTOa2+/hIleWDisVhskhTrdLh6ysdQhHQVFSXv68PR0SAiEokJ4spJtphMnUgkoSQPUILHU2i1HVSnCDQoQrqKiJBbrSjCILJixQqh8P2fb+tms0+npGRTEwioIJNF1NWhCH0MRUhXMpnM5epzu11UB4EJsmjRLdnZerH4PoIoJQgtSX4ZGrr0uef+yGRi7HcQEQolBoPLbDZTHSSg4FeIrlgsVlRUqE7XKZerqM4CE4EkyTfe+PiHH77et+/trq6O9PSM9ev3hoerqc4FE43LlfX29goEuGzGZ1CENBYfr/rxRy2KMKjk56/Iz19BdQqgEpcb2t3dGxkZSXWQwIFDozSmVCrt9i4cHQUIKhKJrK0NA0d9CUVIYxwOR6WS6PUYMgMQREQiqU5ndjqdVAcJHChCetNoVHq9luoUADBxSJJkMiVYsN6HUIT0plKpbLYOzDQBEFS4XJlOhyL0GRQhvXG53IgIUV8fVmUCCCIiUShOE/oQipD2EhJwdBQguEgksvZ2PWYb9hUUIe0plUqrtR2/EgDBg8ViEwTPYMCipL6BIqQ9gUAgl/P7+3GcBCCIsNkyzL7tKyjCQJCYqNbpsBYPQBDh8aQ6XT/VKQIEijAQREdH2e1aXFkPEDyEQnF3Nw6N+gaKMBBwudzY2NDubgyZAQgWAoFYp0MR+gaKMEAkJETr9U1UpwCACcJisV0uptVqpTpIIEARBoiIiAgOx2w2G6kOAgAThM0WY+CoT6AIAwRJkqmpkd3dzVQHAYAJwmCIjEb87esDKMLAERcXYzK14IJCgCDB44l7e7FH6AMowsAhFArVaqFO10F1EACYCAKBSKfDHqEPoAgDSnJyjF6Po6MAQUEgEPf0YI/QB1CEAUWlUnk8vXY7BpIBBD42m+NwkDabjeogtMcayZ1MJtNf//rX77//vr+/f/Lkyc8++2x8fPwV99mzZ8/nn38+8Om7774bGhpKEERZWdkLL7xQV1eXnZ390ksvhYeH+zA9XIHJZKakqBobW6KiEqnOAgDjzjtwlMvlUh2E3ka0R9jR0XHx4sVHHnnkrbfeYjKZCxcuHPo3yMWLF/v7+9f8C4/H835hfn7+ggULdu/ezefzb7/9dt9/B/Bz8fExRmMjhswABAMGg2+xWKhOQXsj2iPUaDQfffSR9+M333xTIBBUVVVNmTLlirslJiauWbNm8JZvvvkmKSlp48aNBEG89tprUqm0rKwsIyPDF8nh6qRSaVSUsKurTaGIpDoLAIwvJpOHQ6NjN+pzhOXl5SwWKyYmZuhNhw4dWrp06X333Xf69GnvFrPZLBAIvB9zOBwOh3P+/PmxxIWRSEvT6HS1VKcAgHHHZvNMJowJGKsR7REOMBqN69evf+6550JCQq64adq0afHx8QqFoqSkZO7cuf/85z9nzJiRn5//xBNPnDlzJicnZ8eOHf39/R0d1xzc73Q6B596vP3225955plRxaM7k8lEkuTYH0cgEPD5Vq22WSoNG/ujwWBWqxWHnf2Z1Wp1u90++T2iBZfL3dmpo9dl9b56oxshHo/HYl2n6UZRhBaL5aabbpo+ffqWLVuG3rps2TLvB4WFhVqtdtu2bTNmzEhPT3/rrbdWr17tdDpnzZqVnZ09zGAZJpN56NChgU/DwsJEItHI4wUAj8fjq285Nzft2LE2lSraJ48GXh6PhyRJPp9PdRC4JpIkuVwugxEs4+GdzlCHQ0uvt0ofvtH5ykiL0Gq1rlq1Kioq6t13371umavV6oFDoOvXr1+/fj1BECaTKTIyMj09/VpfRZKkRqMZYR4YXmRkJIdTZTL1C4USqrMAwHjhcnl6PQ6NjtWI/m6y2+1r166VSCQffPDB4D+16urq/vCHP3g/PnnypNvt9m7csWPHggULvNtbWloIgnA6nU888cSUKVOmTp3q4+8ArobBYGRkxHV21lMdBADGEZvNNZkwWGasRlSEFRUVR48ePXTokEKhkMlkMpnMewyzsbHxjTfe8N7nmWeeEYvFMTExU6ZM+eUvf+kdKUoQxJIlSyIjI6VSaW1t7RdffDFO3wYMFRcX63S24+J6gADGYDAIgmW326kOQm+kD8/8G41GvV6vVquvOEDf3t7O5XK919dfi8PhEAqFQf7jNBgMYrHYhw9YVlZx+TIzJibVh48ZzDwej9VqxTlCf2axWILqHCFBEDU1R5Yvz5JIaHMSxOdvdGPny/8uIpEoKipq6H9BpVI5fAvCOElIiDeZGl0uJ9VBAGC8MBhcXEo4RkH0d1MQEggESUnhHR1YuR4gYKEIxw5FGOAmTUrq66vFTiFAoCJJpsvlojoFvaEIA5xYLE5ODmtvb6A6CACME4Z3xD7cMBRh4Js8OcVgqHM6HVQHAQDfwx7h2KEIA59QKExLi9Bq66gOAgC+R5JMpxNFOCYowqCQlpZsNjfY7TijDhBoGAyGy4VDo2OCIgwKfD4/PT1Sq8WSFACBhsHAHuFYoQiDRUpKks3WjIlmAAIMg8FwOFCEY4IiDBZcLnfKlJjW1hqqgwCALzEYTBwaHSMUYRBJTk50u9utVjPVQQDAZ5hMHBodKxRhEGGz2VOmxLa1VVMdBAB8hsFg4tDoGKEIg0tioobB6DYa+6gOAgC+4V0vmuoU9IYiDC4sFmv69JTW1nKqgwCAb7jdLjabSXUKekMRBp3o6GiFwtPZ2UJ1EADwAbfbzWTinXxM8PIFo9zcjJ6eS5iJGyAAYI9w7FCEwUgqlaalKXApBUAAcLtdLBaKcExQhEEqPT3Vbm82m41UBwGAMcGh0bHDyxekuFzutGmJra0VVAcBgDHBodGxQxEGL40mXiKx6nQdVAcBgLFwMxh4Jx8TvHzBiyTJ6dMnd3RUYFVPAPpyu11MJvYIxwRFGNTkcnlCggRLFQLQl8fjwh7hGOHlC3aZmZPM5jpMQApAUx6Pk81mU52C3lCEwU4gEOTmJjY1nac6CADcCLfbyuPxqE5BbyhCIDSa+IgIT3t7I9VBAGDUXC4U4VihCME7amZKX18VDpAC0A6KcOxQhEAQBCESiXJzNU1NF6gOAgCj4HDYeTwWBsuMEV4++EliYkJ4uLOzs5nqIAAwUna7VSjE7uBYoQjhJyRJzpiRpddX2mwWqrMAwIjYbChCH0ARwr+JRKLs7NjGxjKqgwDAiNjtVrEYRThWKEL4mZSUpLAwW2dnK9VBAOD6cGjUJ1CE8DMkSc6YkanTVdjtVqqzAMB1uFxWPh9FOFYoQriSRCKZOhUHSAFowO028/l8qlPQHooQriI1NTk83N7WVk91EAAYjsNhkEgkVKegPRQhXAVJkrNmTbVYLhuNfVRnAYCrs9ttHI6Hy+VSHYT2UIRwdXw+f+7c9ObmUpfLSXUWALgKk6lfLsfuoA+gCOGaVCrV5Mlh9fWYbgbAH5lM/XK5mOoUgQBFCMPJzEyXSo2dnS1UBwGAK9nthpAQ7BH6AIoQhsNgMGbPztbrL5rNRqqzAMDPuFz9GCnjEyhCuA6xWDxzZkpj4xm32011FgD4icfjcTpNIpGI6iCBAEUI1xcXF5ucLGluvkR1EAD4icVilEr5TCaT6iCBAEUII5Kdnc5md+h0HVQHAQCCIAijESNlfAZFCCPCZrPnzMnu6rpgsZiozgIAhNncq1CEUJ0iQKAIYaRCQ0PnzEltaDiFKwsBKGe398jlcqpTBAgUIYxCTEx0Rob88uUzHo+H6iwAwctutzEYVgwZ9RUUIYzOlCmTIyPdLS3VVAcBCF56fXd0dBhJklQHCRAoQhgdkiRnzsxhMFq7utqozgIQpEymHrUax0V9BkUIo8bhcObPn9bbW44puQEoYbd3h4WFUZ0icKAI4UZIJJIFC6a0tJxxOOxUZwEILjabhcNxicW4dsJnUIRwg5RKZXa2uq7uNGacAZhIen13dDSOi/oSihBuXFpaSlwcu6npItVBAIKIydStUuG4qC+hCOHGkSQ5fXo2n9/d0dFEdRaAoODxeGy2blxB6FsoQhgTFos1FnjAigAAHf1JREFUf/50m60as68BTID+fp1czhMIBFQHCSgoQhgrgUBQUJDb03O+v19HdRaAAKfTtSUmqqlOEWhQhOADUql00aKpWu0ZLFsIMH48Ho/N1q5Wq6gOEmhQhOAbcrk8P39SY+OPNpuF6iwAgUmv746IEOC4qM+hCMFnIiMj8/Li6upOOp0OqrMABKDe3jaNBsdFfQ9FCL6UlJSQmRleW4uLCwF8zO12OxwdOC46HlCE4GMZGZOSk/n19eewQgWAD/X2diqVYh6PR3WQAIQiBN+bOnWKSuVoaCinOghA4Ojr0yYk4LjouBhREXZ1dd16660ymYzH4+Xk5Bw5cmToff785z8nDNLb2+vdvmvXLo1GI5fL4+LiduzY4cvs4K8YDMbMmTlicW9r62WqswAEApfL6XB0qlQ4LjouRlSEFotl/vz5lZWVJpPp7rvvvvnmm81m8xX36e3tnT9/ftG/eFeM1Gq199xzz86dO7u7uz///PMHHnigoaHB598D+CEWi5WfP53JbNZq66jOAkB7nZ0tGk04h8OhOkhgGlERxsTE/Pa3v1UoFEwmc+PGjX19fY2NjUPvJpVKNf/CZDIJgujo6GCxWHPnziUIIi8vTywWa7Va334D4Ld4PN7ChTNJsrGtrZbqLAD01t/fkJIST3WKgDXqc4TffPNNREREQkLC0Js++ugjhUKRk5Pz0UcfebdkZmYuXLhw06ZN+/fvf/jhh3NycqZPnz7Mg/cOgqEWAWCgC7XaBoqjANCWTtchlzNDQ0OpDhKwWKO6d21t7caNG999992he+irVq2688475XL50aNH161bFxoaunLlSpIkCwsL/+d//qelpaW6uvo3v/mNd0/xqpxOp0ajGfj0V7/61SuvvDKqeHRnMplIkqQ6he/l5WUcPnyqvt4SERFLdZaxslqt+BPNn1mtVrfbHWC/R62tl2bNUhuNATJt0wS/0fF4PBbrOk1Hjvy3uqmpKT8/f8uWLQ888MDw93zsscd6enp27tx5+PDhtWvXVldXh4SEGI3G1NTU7du3L1u2bOiXOBwOoVBotwf1Kq8GgyFQF9u0WCxFRSUMhkaliqM6y43zeDxWq5XP51MdBK7JYrFwuVwGI3DGw5vNhq6uH1etKgiYb8oP3+hG+sq2tLQUFBRs2rTpui1IEITH4/EWfnV1dWJiYkhICEEQIpEoNTW1urp6LHGBpvh8/uLFM93uuvb2q5xdBoBraW+vnzw5NmBa0D+N6MXt7u4uKCiYOXNmQUHBmTNnzpw5YzAYCII4ffr07bff7r3PBx98cOnSpc7Ozi+++GL79u233HILQRAzZ848d+7c/v37rVZrUVFRSUnJnDlzxu+bAX/m7UKnsxZdCDBCTqfD4dDGxcVQHSTAjegcYUdHR2xsbHt7+5YtW7xbXn311aysLKfT2d/f791y5syZrVu39vX1JSQkbN++fcWKFQRBZGRkfPzxxy+88ML9998fGRn5/vvv5+TkjNN3Av6Pz+cvWpRXVFTS0UFGROB3G+A62tsbU1KUXC6X6iABbhTnCMcVzhESfnnofDyYTKYDB0q43GTadSHOEfq/QDpH6Ha7qqv/uXLldO9l2QHDD9/oAuG/C9CLUCgsLJzpcl3G9YUAw9BqGzSa0ABrQf+EIgQKCIXCwsLZHE5rQ0OFnxyTAPArTqfDYKjNyEihOkhQQBECNbhcbkHBTJmsr64O61QAXKm19fKkSSqRSER1kKCAIgTKsNnsefPy1GrH5cun3W4X1XEA/IXdbrXbm9PSkqgOEixQhEAlJpM5e3ZuUhKnuvqEwxHUQ6UABrS0VGdkxGDpwQmDIgSKkSSZk5OZmSm7fLnEbrdSHQeAYhaLyeNpT06+ynzOME5QhOAX0tPTZsyIqq09brGYqM4CQKXW1sqpUxPYbDbVQYIIihD8RWJiwty5iY2NJSZTP9VZAKjR39/L4fTGx8dRnCPIoAjBj8TGxixcmN7W9mNPTzvVWQAmmtvtbm29kJc3aZhVemA8oAjBvyiVyuXL88zmiuZmzM8OwaW19XJcHF+tVlMdJOigCMHvSCSSZcvmiESdtbWluKwCgoTZbLTbG3NyMqgOEoxQhOCPuFxuQcGs+HiypuaE3W6jOg7A+PJ4PM3NF6ZPT8ZMtpRAEYKfYjAYubnZubnKurqjRmMf1XEAxpFW26BUEnFxsVQHCVIoQvBrSUkJBQWT2ttPYvgMBCqbzWIw1EyfPoXqIMELRQj+TqVSLVuG4TMQsJqby6dN02BaUQqhCIEGJBLJ0qWzRaLO2tqzGD4DgaSjozkkxJKUhHlkqIQiBHrg8XgLFsxMSCCrq4+azUaq4wD4gNlsMBgqZ83KJkmS6ixBDUUItMFkMnNysubNS2hpOd7e3kR1HIAxcbmcjY1nZs1K87fl2oMQihBoJjo6auXKWQxGfV3dORwmBfqqrz+fkSGPjo6iOgigCIGGRCLR4sWz4+M91dVHzWYD1XEARq2trS4szJKRMYnqIEAQKEKgKRaLNW1a9rx5Ca2tJzo7W6mOAzAK/f29NlvtrFlTGQy8A/sF/BiAxqKjo37xixkEcRmHSYEu7HZbW1tpfn6WQCCgOgv8BEUI9CYWixcvnh0T466uPm61mqmOAzAcj8dTX1+akxMdHh5OdRb4NxQh0B6LxcrLmzprVnRj49GODowmBf9VX38hNpaVkpJEdRD4GRQhBIj4+LiVK2ey2Y3V1T/a7Vaq4wBcqanpUmioYcaMqbhq0N+gCCFwiMXiRYvmTJsmr68/ggsNwa9otQ1cbmd+fh4W3fVDKEIIKCRJJiUlrFw5k8VqqK4+iSWcwB90drY6nbULFkxns9lUZ4GrQBFCABKLxYWFc7OyQmprj2DZCqCWTtdpNF5ctCgPaw36LRQhBCaSJNPSkn/xi2lOZ2VNzRmn00F1IghGBoO+p+fcokXTsLiEP0MRQiALDQ1dsmTupEnc2tojvb1dVMeB4GI2G1pbTy1cmB0aGkp1FhgOihACHJPJnDIlvbAw02Qqq609i7OGMDEMBn1T04kFCybjkkH/hyKEoCCXy1esyJ8yRdjYeKStrd7j8VCdCAJZX19PW9vJhQsz1Go11Vng+lCEECyYTGZqavLKlbOk0s6qqqP9/b1UJ4LApNN1dHeXLl06TalUUp0FRgRFCMFFKBTOm5dXUJCo052urT3ncNipTgQBpbOzpb+/bOnSPJlMRnUWGCkUIQQjlUp1003z09LYtbWHcek9+EpbW73dXr106UyJREJ1FhgFFCEEKTabnZExecWKPC63qarqONY1hDFqa6tlMhuXLJklFAqpzgKjgyKEoCaVShcunD1zprqtraSx8SIuN4Qb4Ha7amvP8vnaRYtm8Xg8quPAqLGoDgBAMZIk4+Pj1GpVZWXNxYuHRSKNUhmPFVNhhCwWU0PD6dRUaXb2TMwjSlMoQgCCIAgul5uZmZ6YGF9WVllT80+pNCkiIhqrBMDwdLrOjo5zM2YkJCYmUJ0FbhyKEODfhELhjBk5qan6s2cvVVbWh4cny+UqqkOBn2prq3W5Gn7xi1xMHEN3KEKAK4WEhCxYMLO7u/vUqYs1NY0qVdr/b+/Og5so3weAb5JNNleTNk2TpmnTNj1S7IHQYmkFSynUsxwV/0Hw1kEUT5QZGcf7BEcdFUQUGVQc7QAFFXEsCC1UKy1TBIHeoWmSpm3aNPdmN7vfP9Zffp0WsAg5aJ7POMzbd9/gE95mn7zvvruvWCwNd1AggpAk0dV1QqXyl5bOwTAs3OGAKwWJEIALk8vlt9wy12AwnDhxfGAgPilJh2GwewBAnM7R3t6W669PzMubBpPnUwMkQgAuisViaTQatVrd2dnd2nqUzVbGxSXDZjpRi6Zpo7HT6+2pqMiDZ6dNJZAIAfgXHA5Hp8vKyEjv6dE3NTVarUlJSdlCYUy44wIh5XY7entbNRpeUdFc+DI0xUAiBGBSUBTNzMxQKBKs1uHW1iaCkCqVWTExseGOCwQdTdNmc7fb3TVnTk5qqibc4YCrDxIhAJcBRVGtNj0tLbWvr+/EiRajUZCYmB0bKw93XCBYXC67wXAyNRVbuPAmuFl+qoJECMBlY7PZGo0mOTnZZDK1tp4aGMDk8gyZTBnuuMDV9H8Dwe65c3UaDQwEpzJIhAD8R2w2Ozk5Wa1Wm0ymkyfPtbV1xMWly+UqeCrNFGC19g8MnM3MlMyYcRPcIDHlQSIE4IqwWCy1Wp2UlDQwMNDeru/oOCMQaJRKDdxrcY2y20fM5rMyGXnzzXmwuXyUgEQIwFXAYrGUSqVSqXS5XD09vWfPNhCEJCEhHeZLryE47jEa29jsoTlzslNS4AF7UYRF03S4Y0AQBCEIQiQS+XxRvUuqw+GIiYFF+ZGLpmm32z2ZTXYoijKZTKdOdY+M+MVijVKpQVFuCCIEHo8Hw7DLnZ0mScJk6vT5DHl5KTpdForCCCGIIvBEB/0NwNXHXD5MTk62Wq3t7fqurk4MUysUqXD3YaShKL/ZrLfbu/Lzk6ZNm8fj8cIdEQgDGBFGkAj8ogTGmvyIcByv16vX954504vjfLFYLZcn8Xiw/iIoJj8i9Pnw/v4et7s3M1Oel6eD3XRDJgJPdDAiBCDo+Hx+Tk62Tpc1MjKi1/d1dBxmsaQxMcnx8YkcDnwGQ83jcQ0M6L3ePp1OqdOVisXicEcEwgw+hACECIvFkslkMpns+uvzBgcHu7v7OjpOc7kJcXHJcXEKWJoRAnb7sMXSxeHY8vNTtdpymAgFDJgajSAROGMAxvrPU6MXg+O42Wzu6OgbGvLy+UkymRr2e7pCF5wapShqaMhktXbLZHR+vlatVsO9nmEUgSc6GBECEDYYhqWlpaWlpblcLoPB2N5+oq+PwjBFbGyiVBoPJ+srRNO03T48PGzEcXNKSmxhYY5CoQh3UCASwYgwgkTgFyUw1lUfEU7kdrv7+y29vRaj0cbhxIpESrlcxePBIy4nixkRer1uq9XkcvXFxXEyMlQaTQrsFxE5IvBEByNCACKIUCjUatO12nSCIIaGhgwGy/nz7QTBxzCFTKaUSGThDjCikSRhsRhwfADDvNnZKo2mUCqFqWbw7y5jRDgwMGA0GjUaTXx8/MSjBEE4nc7Aj2KxmMvl0jRts9nGNhMIBBd8gjuMCJGI/KIExgrBiPCC/1Or1Wo2W7q7LU4njWFyoVAmkcgEAlju/w+Xyz4yMuD1DtK0Xa2W6HQZCQkJsPgoYkXgiW6yI8LZs2d3dHRotdr29vYHHnjg/fffH9dg//79y5YtC7y9mpqaioqKkZGRzMzMQBubzfbxxx+vXr36qoQOQDRgsVhyuVwul+fn5zqdzuHhYYtlyGhsd7loHk/G58tiY+Oj8D59gvDZbEMOx6DHMyCVolqtQqnMiI+P93q9AoEAsiC4LJNNhOvXr7/99tvZbPb58+cLCgoWL148b968cW1KSkrq6+vH1shksuHhYabc3t6el5dXXV19xTEDEKXEYrFYLNZoNLNmIW63e3h4eGBg2GzWGww4lyvDsHiJRCYSSabqKhuKopzOUZtt0OsdRBBHcnL8tGmKhIQsoVAY7tDAtW2yibCqqooppKamJicnm83mCzYzGAzx8fEX/L384osvqqqqEhMT/1ugAICxhEKhUChMTk5GEATH8eHh4cFBq9HY19fn4nCEKCrh8SQikUQkkly7T7GhadrtdjgcNo9nlCRtJOmUy8XZ2XKlUieTyaZqvgehd9mLZRoaGkwmU3l5+cRDf/75Z1lZmdlsrqys/PLLL2Wy/7+wT5LkV199tXXr1kv8zTRN19XVBX7U6XQpKSmXGx4AUQjDMJVKpVKpCgoQiqIcDofdbrfZHFZrl9Fo9/lYKCpBUYlAECMSSQQCccSmEIqivF6X0znqctkIYpQkHbGxAoVCqlTGSqXJUqk0YiMH17TLu32iu7u7rKzsnXfeWb58+bhDIyMjGIYJhUKbzVZdXZ2RkTE27e3du3fVqlUGg+Fij3UnCILP55eVlQVqFixYsGbNmst5L9c8p9MJT3uKZDRNezyea24izuv1OhwOp9NptdqHhpyjoy4Wi8dmC1BUwGYLeDwBhgkwTIhh/FA+740gcK/X4/W6cdxNkm6/3+33u2kal0iECQkxcrlUKpXGxMRc7kYQbrebz+dDvoxkIT7R8fl8Lvdf9n65jERoMBjKysrWrl37r6tdamtrn3vuuY6OjkDN4sWLc3Nz33zzzYu9BFaNIhG5mAqMFZZVo8GA47jb7fZ4PB6Px+Fwj456HA6Pw+GmKA6HI2CzMRaLy2KhbDaXw+GiKJfDQQN/oij30mmGpmm/nyRJgiB8JEkw/1EUQVEETRM0Tfj9OEm6+XyORCKUSoVSqVAkEopEIqFQyOfzr3Cdi8vlEggEkAgjWQSe6Cb7bctisVRWVq5atWoyaz6NRmNcXNzY1x44cGDjxo3/MUYAwFWFYRiGYWM/pAyfz+fxeHAcJwiCSWU47sNxl89Her0EjhNOJ4njBE1Tl/77eTwuhnExjCsS/VPAMC6XK+RyuVwul5k6gj3/QOSY1O+i3+8vKyuLiYmJjY397LPPEAQpLi6ePn261WrNy8tramrSaDQbN26USqUajebMmTOvvvrq2Psrtm/fXlJSkpWVFaw3AQC4Gng8HjyHGkShSSVCmqaZq3ctLS1MTWpq6vTp0zEMu+uuu5iZIp1Ot2vXrj179iQlJX3//fcLFy4MvJzFYq1fvz4IwQMAAABXCp41GkEicOocjDVlrhFOYXCNMPJF4IkOfl0AAABENUiEEeTDDz/EcTzcUYCLOn/+/LfffhvuKMClfPfddz09PeGOAlwUQRATn9AZdpAII8imTZusVmu4owAXdfr06V27doU7CnAptbW1f/31V7ijABdls9k++uijcEcxHiRCAAAAUQ0SIQAAgKgGiRAAAEBUi5TbJ/x+P5/P12g04Q4knHp7e9VqNYfDCXcg4MI8Ho/dblcqleEOBFyUxWKJiYm55p4HGz0oiurr6wvlqX758uWvvfbapdtEylOOOBxOV1cXSZLhDiSccBzHsGt1x5woAX0U4Xw+H5fLhY15I1mIP0Qqlepf20TKiBAAAAAIC7hGCAAAIKpBIgQAABDVIBECAACIapAIAQAARLVIWTUabRobGw8fPjw0NJSYmLhixYqkpKSJbYxG47Zt20ZHR6urq0tLS8cewnF8x44dubm54+rBVfTLL7/88ccfo6OjaWlpK1eunLiNLYIg586d+/rrrymKuvvuu3NzcxEEIUny8OHDx44dc7vds2fPXrJkCaxgDJLu7u4ff/xRr9dLJJLbb7991qxZE9vgOP755593dXUVFhYuX7480BenTp3auXMniqIrV67Mzs4ObeDRgqbp5ubmuro6q9Wam5u7fPnyCy4W/f3333/88UemLzIzM5lKl8u1Y8eOzs5OlUp1zz33KBSKoIYKI8Lw2L17t8fj0Wq1Z8+ezcvL0+v14xoMDw/PmjXLYrGo1eo77rjjwIEDY4++/vrra9eu/e6770IXcfTZuXMni8XSarUHDx4sLCy02+3jGrS1tRUXF3M4HIFAUFpaevr0aQRB6urqnn32WZIkFQrF888//8gjj4Qj9qjQ0NDQ3d2dnp7u8/kqKipqamomtlm2bNmePXuysrLefffddevWMZWtra033nijWCymabq4uLirqyu0gUcLk8m0bNmy4eHhlJSUrVu3VlRUTLxBbufOndXV1YmJiWw2u7S0tLOzE0EQiqLKy8v3799fUFDQ1tY2c+bMkZGR4MZKg3ArLS3dtGnTuMqNGzdWVlYy5Y8//risrCxwqLW1taio6MEHH3ziiSdCFmQ08/v9KpXqwIED4+pXr169atUqpvzMM8888MADNE17vd5AgxMnTqAo6na7QxZq1HrppZeWLl06rvLkyZNisdjhcNA03dbWJhQKR0ZGaJpesWLF888/z7R56KGHnnzyyRBHGyUIgiAIgik7HA4+n9/S0jKuzXXXXbdjxw6m/OSTTz722GM0TXd3dyMIwnQcTdOpqak//fRTUEOFEWGYGY1GvV6fl5c3rr6+vn7hwoVMubKy8ujRo36/H0EQkiQfeeSRTz/9lMvlhjrWaHXmzBm3263T6cbV19fXV1ZWMuWFCxceOXIEQZCxkz9er5fH40FPBZvX621ubs7Pzx9XX19fX1paKhaLEQTJzs5WKBTNzc3IRToOXHUoiqLoP1ffSJIkSZLpi7FMJpNWq2XKGRkZv/32G4IgSqVSLpe3trYiCNLb22uz2YI9fQ2JMGw2bNiQkpKSnp6+Zs2auXPnjjtqNpsTEhKYskKh8Pv9AwMDCIK8/fbbZWVlhYWFoQ43Kj3xxBNqtbqoqOiTTz5JS0sbd3RcH5nN5rFHCYJ45pln1q5dGzgXgKvujz/+0Gq1cXFxJEm++OKL446O7SAEQRQKhclkomm6v7//Eh0HguGpp55atGjRxHyWn59/6NAhBEFomj548KDJZEIQRCgU7t27t7q6Oicnp6CgYPPmzYFrh0ECiTBYVq9ejU4wY8aMQIPHH3+8paVl7969H3zwwf79+8e9HEXRwHw6U+DxeOfOnfvmm29eeumlkL2Lqe2WW26Z2Ed33nlnoMFbb73V3Ny8bdu2NWvWTNzlblwf8Xi8wCG/33/vvfdKpdL169eH4I1MVXv27JnYQSiKWiwWpkFRUVFzc/PRo0cdDkfgEmAAl8tl5lEYBEHweDwWi3WJjgOXa9GiRRM7qKqqamybl19++cSJE1u3bp348vfee2/z5s233nprSUmJy+USCAQIgtjtdmb6et++fZs2bXr66afb2tqC+zaCOvEKJuPZZ5+97777xlUuW7bs1VdfZcqtra18Pp+iqFdeeUWlUhUWFhYWFiYkJCgUihUrVoQ83mi0ZMmSQHcEzJw58+uvv2bKu3fvvu6665iy3+9fuXJlZWWlx+MJaZRR7Oeff1ar1eMqt2zZMnfuXKZMUZRSqTxy5AhN05mZmfv27WPqt2/fXlxcHMpQo81bb701bdq0/v7+izVwOp2NjY09PT1btmy56aabaJquqanJzs4ONFiyZMkrr7wS1CBhRBgGNE17PB6mTFFUS0sL8yx2HMcPHTqE4ziCIFVVVXv27GG+t9bU1Nxxxx0sFuvhhx/+4YcftmzZwvzGzJ8//4UXXgjjG5nCfD5fYNDg8Xj+/vvv1NRUBEFsNtvhw4eZ+qqqqsBKxZqaGuZbME3Tq1ev1uv1u3fv5vP5YQg9arjd7kC5ubk5sKFBU1MTM9t52223HT9+3GAwIAjS0NBA0/Ts2bORi3QcCIYPPvhg27ZtdXV1Y/dssdvtzLVAhkgkKikpkclkH3300d13340gSHx8/ODgoMPhQBCEoqienh65XB7UOOGh22FAEIRSqZwzZ05sbGxTU1NMTMyvv/4aFxdnMBg0Go3BYEhOTvZ6vfPnz2exWKmpqXV1dQcPHhy3FuDRRx/l8XgffvhhuN7F1Nbe3l5eXl5cXCwQCOrr6wsKCmpra7lc7pEjRyorK5kvK4ODg6WlpVlZWSiKnjp1qrGxUaVS1dbWLl26dNq0aYGdgHbv3h3l+4sFyc0330ySpFqt1uv1bW1ttbW1JSUlCILk5OSsW7fu/vvvRxBk3bp1NTU18+bN279//xtvvPHggw8iCNLX11daWjpz5kyv16vX648dOxYfHx/mNzMVdXZ2ZmVlpaWlBf55N2zYUF5e3tjYeOONNzKpZ8uWLTU1NQqF4ujRo/Pnz9+2bRubzaYoavHixR0dHRUVFS0tLR6Pp6GhQSKRBC9USITh0d/f39zc7HQ609PTb7jhBuY+X5/P19raOmPGDGadIUEQhw4dGh0draiomPhB7e3tZbFYKSkpYYg+Ouj1+pMnT/p8vuzs7OnTpzOVDoejra2tqKiI+dHlctXV1dE0vWDBAmZF3MjICLP4OyA3NxeGhsHgcrmamposFotCoZg9e7ZIJGLqT506pVKpAmOI48ePd3Z2zpgxIycnJ/Bah8NRV1eHouiCBQuY61LgqvN6vX///ffYmoyMjNjYWKfTefbsWeYBCG63+9ixY0NDQzk5OWOXUCAIcvz48a6uLpVKNWfOnGDv0gqJEAAAQFSDa4QAAACiGiRCAAAAUQ0SIQAAgKgGiRAAAEBUg0QIAAAgqkEiBAAAENUgEQIAAIhqkAgBAABENUiEAAAAohokQgAAAFENEiEAAICo9j/C85EPyf3j0QAAAABJRU5ErkJggg==", + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n" + ], + "text/html": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "m = mean(convert(MvNormalMeanCovariance, results.posteriors[:β][end]))\n", + "Σ = cov(convert(MvNormalMeanCovariance, results.posteriors[:β][end]))\n", + "\n", + "p = plot()\n", + "# i = 2\n", + "covellipse(m, Σ, aspect_ratio=1,n_std=3, label=\"3σ Contours\", color=:blue, fillalpha=0.2)\n", + "scatter!([m[1]], [m[2]], label=\"Mean estimate\", color=:blue)\n", + "scatter!([true_beta[1]], [true_beta[2]], label=\"True Parameters\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# References\n", + "[1] Polson, N. G., Scott, J. G., & Windle, J. (2013). Bayesian inference for logistic models using Polya-Gamma latent variables. *Journal of the American Statistical Association*, 108(1), 136-146.\n", + "\n", + "[2] Minka, T. (2001). Expectation Propagation for approximate Bayesian inference. *Uncertainty in Artificial Intelligence*, 2, 362-369." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 1.11.0", + "language": "julia", + "name": "julia-1.11" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.11.0" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} From 797b04cce7ed5c1a88f10c381ea63f095ae95da6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C4=B0smail=20=C5=9Een=C3=B6z?= Date: Mon, 13 Jan 2025 16:50:17 +0100 Subject: [PATCH 3/8] remove unnecesary argument --- test/models/regression/binomialreg_tests.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/test/models/regression/binomialreg_tests.jl b/test/models/regression/binomialreg_tests.jl index cfede8878..6786a4945 100644 --- a/test/models/regression/binomialreg_tests.jl +++ b/test/models/regression/binomialreg_tests.jl @@ -5,12 +5,11 @@ function generate_synthetic_binomial_data( n_samples::Int, - n_features::Int, true_beta::Vector{Float64}; seed::Int=42 ) rng = StableRNG(seed) - + n_features = length(true_beta) # Generate design matrix X X = randn(rng, n_samples, n_features) @@ -42,7 +41,7 @@ end end - function binomial_inference(binomial_model, iterations, X, y, n_trials) + function binomial_inference(binomial_model, iterations, X, y, n_trials, n_features) return infer( model = binomial_model(prior_xi = zeros(n_features), prior_precision =diageye(n_features),), data = (X=X, y=y,n_trials=n_trials), @@ -52,17 +51,18 @@ ) end - function run_simulation(n_sims::Int, n_samples::Int, n_features::Int, true_beta::Vector{Float64};iterations = n_iterations) + function run_simulation(n_sims::Int, n_samples::Int, true_beta::Vector{Float64};iterations = n_iterations) # Storage for results + n_features = length(true_beta) coverage = Vector{Vector{Float64}}(undef, n_sims) fes = Vector{Vector{Float64}}(undef, n_sims) for sim in 1:n_sims # Generate new dataset - X, y, n_trials = generate_synthetic_binomial_data(n_samples, n_features, true_beta, seed = sim) + X, y, n_trials = generate_synthetic_binomial_data(n_samples, true_beta, seed = sim) X = [collect(row) for row in eachrow(X)] # Run inference - results = binomial_inference(binomial_model, iterations, X, y, n_trials) + results = binomial_inference(binomial_model, iterations, X, y, n_trials, n_features) # Extract posterior parameters post = results.posteriors[:β][end] m = mean(post) @@ -81,7 +81,7 @@ end; - coverage, fes = run_simulation(n_sims, n_samples, n_features, true_beta) + coverage, fes = run_simulation(n_sims, n_samples, true_beta) for i in 1:n_sims @test fes[i][end] < fes[i][1] end From bc74270fb8d4aa3350511509eef930b20a1ca06b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C4=B0smail=20=C5=9Een=C3=B6z?= Date: Mon, 13 Jan 2025 17:25:52 +0100 Subject: [PATCH 4/8] formatter --- test/models/regression/binomialreg_tests.jl | 55 +++++++++------------ 1 file changed, 22 insertions(+), 33 deletions(-) diff --git a/test/models/regression/binomialreg_tests.jl b/test/models/regression/binomialreg_tests.jl index 6786a4945..0e9bcca08 100644 --- a/test/models/regression/binomialreg_tests.jl +++ b/test/models/regression/binomialreg_tests.jl @@ -3,55 +3,48 @@ include(joinpath(@__DIR__, "..", "..", "utiltests.jl")) - function generate_synthetic_binomial_data( - n_samples::Int, - true_beta::Vector{Float64}; - seed::Int=42 - ) + function generate_synthetic_binomial_data(n_samples::Int, true_beta::Vector{Float64}; seed::Int = 42) rng = StableRNG(seed) n_features = length(true_beta) # Generate design matrix X X = randn(rng, n_samples, n_features) - + # Generate number of trials for each observation n_trials = rand(rng, 5:20, n_samples) - + # Compute logits and probabilities logits = X * true_beta probs = 1 ./ (1 .+ exp.(-logits)) - + # Generate binomial outcomes y = [rand(rng, Binomial(n_trials[i], probs[i])) for i in 1:n_samples] - + return X, y, n_trials end n_samples = 1000 n_features = 2 - true_beta = [-1.0 , 0.6] + true_beta = [-1.0, 0.6] n_iterations = 100 n_sims = 20 - - @model function binomial_model(prior_xi, prior_precision, n_trials, X, y) + @model function binomial_model(prior_xi, prior_precision, n_trials, X, y) β ~ MvNormalWeightedMeanPrecision(prior_xi, prior_precision) for i in eachindex(y) - y[i] ~ BinomialPolya(X[i], n_trials[i], β) where { - dependencies = RequireMessageFunctionalDependencies(β = MvNormalWeightedMeanPrecision(prior_xi, prior_precision)) - } + y[i] ~ BinomialPolya(X[i], n_trials[i], β) where {dependencies = RequireMessageFunctionalDependencies(β = MvNormalWeightedMeanPrecision(prior_xi, prior_precision))} end end function binomial_inference(binomial_model, iterations, X, y, n_trials, n_features) return infer( - model = binomial_model(prior_xi = zeros(n_features), prior_precision =diageye(n_features),), - data = (X=X, y=y,n_trials=n_trials), + model = binomial_model(prior_xi = zeros(n_features), prior_precision = diageye(n_features)), + data = (X = X, y = y, n_trials = n_trials), iterations = iterations, free_energy = true, - options = (limit_stack_depth = 100,), - ) + options = (limit_stack_depth = 100,) + ) end - function run_simulation(n_sims::Int, n_samples::Int, true_beta::Vector{Float64};iterations = n_iterations) + function run_simulation(n_sims::Int, n_samples::Int, true_beta::Vector{Float64}; iterations = n_iterations) # Storage for results n_features = length(true_beta) coverage = Vector{Vector{Float64}}(undef, n_sims) @@ -60,26 +53,24 @@ # Generate new dataset X, y, n_trials = generate_synthetic_binomial_data(n_samples, true_beta, seed = sim) X = [collect(row) for row in eachrow(X)] - + # Run inference - results = binomial_inference(binomial_model, iterations, X, y, n_trials, n_features) + results = binomial_inference(binomial_model, iterations, X, y, n_trials, n_features) # Extract posterior parameters post = results.posteriors[:β][end] m = mean(post) v = var(post) - estimates = map((x,y) -> Normal(x,sqrt(y)), m, v) - coverage[sim] = map((d,b) -> cdf(d, b), estimates,true_beta) + estimates = map((x, y) -> Normal(x, sqrt(y)), m, v) + coverage[sim] = map((d, b) -> cdf(d, b), estimates, true_beta) fes[sim] = results.free_energy - end - + return coverage, fes end - function in_credible_interval(x, lwr=0.025, upr=0.975) + function in_credible_interval(x, lwr = 0.025, upr = 0.975) return x >= lwr && x <= upr - end; - + end coverage, fes = run_simulation(n_sims, n_samples, true_beta) for i in 1:n_sims @@ -87,11 +78,9 @@ end coverages = Vector{Float64}(undef, n_features) for i in 1:n_features - coverages[i] = sum(in_credible_interval.(getindex.(coverage,i))) / n_sims + coverages[i] = sum(in_credible_interval.(getindex.(coverage, i))) / n_sims @test coverages[i] >= 0.8 end @test_benchmark "models" "binomialreg" binomial_inference(binomial_model, $n_iterations, $X, $y, $n_trials) - - -end \ No newline at end of file +end From 62734c58444790a259bfa8b44f67b4c92852ed80 Mon Sep 17 00:00:00 2001 From: Bagaev Dmitry Date: Tue, 14 Jan 2025 13:25:44 +0100 Subject: [PATCH 5/8] move example out --- .../basic_examples/Binomial Regression.ipynb | 473 ------------------ 1 file changed, 473 deletions(-) delete mode 100644 examples/basic_examples/Binomial Regression.ipynb diff --git a/examples/basic_examples/Binomial Regression.ipynb b/examples/basic_examples/Binomial Regression.ipynb deleted file mode 100644 index 93fd5a144..000000000 --- a/examples/basic_examples/Binomial Regression.ipynb +++ /dev/null @@ -1,473 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Bayesian Binomial Regression " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "This notebook is an introductory tutorial to Bayesian binomial regression with `RxInfer`." - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "\u001b[32m\u001b[1m Activating\u001b[22m\u001b[39m project at `~/.julia/dev/RxInfer/examples`\n" - ] - } - ], - "source": [ - "# Activate local environment, see `Project.toml`\n", - "import Pkg; Pkg.activate(\"..\"); Pkg.instantiate(); " - ] - }, - { - "cell_type": "code", - "execution_count": 10, - "metadata": {}, - "outputs": [], - "source": [ - "using RxInfer, ReactiveMP,Random, Plots, StableRNGs, LinearAlgebra, StatsPlots, LaTeXStrings" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Likelihood Specification\n", - "For observations $y_i$ with predictors $\\mathbf{x}_i$, Binomial regression models the number of successes $y_i$ as a function of the predictors $\\mathbf{x}_i$ and the regression coefficients $\\boldsymbol{\\beta}$\n", - "\n", - "$$\\begin{equation}\n", - "y_i \\sim \\text{Binomial}(n_i, p_i)\\,,\n", - "\\end{equation}$$\n", - "where:\n", - "$y_i$ is the number of successes, $n_i$ is the number of trials, $p_i$ is the probability of success. The probability $p_i$ is linked to the predictors through the logistic function:\n", - "$$\\begin{equation}\n", - "p_i = \\frac{1}{1 + e^{-\\mathbf{x}_i^T\\boldsymbol{\\beta}}}\n", - "\\end{equation}$$\n", - "\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Prior Distributions\n", - "We specify priors for the regression coefficients:\n", - "\n", - "$$\\begin{equation}\n", - "\\boldsymbol{\\beta} \\sim \\mathcal{N}_{\\xi}(\\boldsymbol{\\xi}, \\boldsymbol{\\Lambda})\n", - "\\end{equation}$$\n", - "as a Normal distribution in precision-weighted mean form.\n", - "\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Model Specification\n", - "\n", - "The likelihood and the prior distributions form the probabilistic model\n", - "$$p(y, x, \\beta, n) = p(\\beta) \\prod_{i=1}^N p(y_i \\mid x_i, \\beta, n_i),$$\n", - "where the goal is to infer the posterior distributions $p(\\beta \\mid y, x, n)$. Due to logistic link function, the posterior distribution is not conjugate to the prior distribution. This means that we need to use a more complex inference algorithm to infer the posterior distribution. Before dwelling into the details of the inference algorithm, let's first generate some synthetic data to work with." - ] - }, - { - "cell_type": "code", - "execution_count": 61, - "metadata": {}, - "outputs": [], - "source": [ - "function generate_synthetic_binomial_data(\n", - " n_samples::Int,\n", - " true_beta::Vector{Float64};\n", - " seed::Int=42\n", - ")\n", - " n_features = length(true_beta)\n", - " rng = StableRNG(seed)\n", - " \n", - " X = randn(rng, n_samples, n_features)\n", - " \n", - " n_trials = rand(rng, 5:20, n_samples)\n", - " \n", - " logits = X * true_beta\n", - " probs = 1 ./ (1 .+ exp.(-logits))\n", - " \n", - " y = [rand(rng, Binomial(n_trials[i], probs[i])) for i in 1:n_samples]\n", - " \n", - " return X, y, n_trials\n", - "end\n", - "\n", - "\n", - "n_samples = 10000\n", - "true_beta = [-3.0 , 2.6]\n", - "\n", - "X, y, n_trials = generate_synthetic_binomial_data(n_samples, true_beta);\n", - "X = [collect(row) for row in eachrow(X)];" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We generate `X` as the design matrix and `y` as the number of successes and `n_trials` as the number of trials. Next task is to define the graphical model. RxInfer provides a `BinomialPolya` factor node that is a combination of a Binomial distribution and a PolyaGamma distribution introduced in [1]. The `BinomialPolya` factor node is used to model the likelihood of the binomial distribution. \n", - "\n", - "Due to non-conjugacy of the likelihood and the prior distribution, we need to use a more complex inference algorithm. RxInfer provides an Expectation Propagation (EP) [2] algorithm to infer the posterior distribution. Due to EP's approximation, we need to specify an inbound message for the regression coefficients while using the `BinomialPolya` factor node. " - ] - }, - { - "cell_type": "code", - "execution_count": 62, - "metadata": {}, - "outputs": [], - "source": [ - "@model function binomial_model(prior_xi, prior_precision, n_trials, X, y) \n", - " β ~ MvNormalWeightedMeanPrecision(prior_xi, prior_precision)\n", - " for i in eachindex(y)\n", - " y[i] ~ BinomialPolya(X[i], n_trials[i], β) where {\n", - " dependencies = RequireMessageFunctionalDependencies(β = MvNormalWeightedMeanPrecision(prior_xi, prior_precision))\n", - " }\n", - " end\n", - "end" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Havin specified the model, we can now utilize the `infer` function to infer the posterior distribution." - ] - }, - { - "cell_type": "code", - "execution_count": 80, - "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "\u001b[32mProgress: 100%|█████████████████████████████████████████| Time: 0:00:02\u001b[39m\u001b[K\n" - ] - }, - { - "data": { - "text/plain": [ - "Inference results:\n", - " Posteriors | available for (β)\n", - " Free Energy: | Real[21992.9, 16235.8, 13785.0, 12519.7, 11800.1, 11366.2, 11094.1, 10918.7, 10803.3, 10726.3 … 10558.2, 10558.2, 10558.2, 10558.1, 10558.1, 10558.1, 10558.1, 10558.1, 10558.1, 10558.1]\n" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "n_features = length(true_beta)\n", - "n_max_iterations = 50\n", - "strategy = StopEarlyIterationStrategy(1e-5, [])\n", - "results = infer(\n", - " model = binomial_model(prior_xi = zeros(n_features), prior_precision = diageye(n_features),),\n", - " data = (X=X, y=y,n_trials=n_trials),\n", - " iterations = n_max_iterations,\n", - " free_energy = true,\n", - " options = (limit_stack_depth = 100,),\n", - " showprogress = true,\n", - " \n", - ")\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We can now plot the free energy to see if the inference algorithm is converging." - ] - }, - { - "cell_type": "code", - "execution_count": 81, - "metadata": {}, - "outputs": [ - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3deXwUVb738VPV3Vk7ndUkkLAlQcANhaCyCAwQ9QpBQRgQwXFDHF/3+uCjw+iAOj6My1zEG/W6juCGuIDKCIqKgGIQIqCIsgkGQSCQhJCt03vV80dhG0ISOqZJdXd93n/4qjpUqn9ddve365yq05KqqgIAAKOS9S4AAAA9EYQAAEMjCAEAhkYQAgAMjSAEABgaQQgAMDSCEABgaAQhAMDQCEIAgKERhAAAQwunIHz22Wf37NnjX/V4PDoWE444Ym3l9XqZg7BNfD6foih6VxFOFEXx+Xx6VxFmvF5vcHcYTkG4cuXK3bt3+1edTqeOxYQjjlhbuVwugrBNPB4PH+tt4vP5+IbaJoqiuFyu4O4znIIQAICgIwgBAIZGEAIADI0gBAAYGkEIADA0ghAAYGgEIQDA0AhCAIChRWwQfnJQPeLQuwgAQMiL2CBcsFtZV8ZUTwCA0zDrXcCZkhglatx6FwEAequtrb399tuDPj9nx7NarQsXLjwTe47YIEyKEtUEIQDDq6ioWLNmzdNPP613Ie01adKkl1566UzsOWKDMDFKqnEzXTIACKvVOnHiRL2raK/JkyefoT1H7BghXaMAgEAQhAAAQ4voIORHvgAApxOxQZgUJVW7GCMEAJxGxAYhXaMAgEBEdBDSNQoAOJ2IvX2CrlEACCOHDx9et27d9u3bc3Nzb7zxxo586IgNQluUqPMIVQhJ70oAAKe1aNGitWvX2u32b7/9toODMGK7Rk2SiDOLOnpHASCU7NixY+XKlf7VXbt2ffjhh0KIWbNmrVy5cuzYsR1fUsQGoWByGQAIPfHx8dddd11dXZ22ev/99+/cuVPfkiK2a1T8euFol3i96wCAUPKLXb1lna9jzhKiZLFkpDmuUdR069Zt0KBBb7311vTp048ePfrxxx8/88wzHVJLiyI/CAEAjXWKlf7a16R2SBJGm0TcKTnz5z//ee7cudOnT1+4cGFhYWF6enpHlNKySA5CfoACAE5llsXIznpeRzh69Og777xzy5YtCxcuPEM/KNEmjBECADqULMu33nrrzTffbDabhw4dqnc5ER6EdI0CQCi67bbbdu/efdttt0nSiXPTf//737m5uY899tjatWtzc3PvvffeDiuGrlEAQEdzOBxRUVE33HCDv+U//uM/Gp8dRkdHd1gxkRyEiVFSNV2jABBili1b9vTTT998882pqan+xqioqKioKF3qiewgFPvr9S4CAHCy7du3FxYWzpgxQ+9CTojwIGSMEABCzezZs/Uu4SSRfLFMEl2jAIDTieQg5IwQAHBaBCEAwNAIQgCAoUXyxTKMEQKAyWQ6ePBgfn6+3oW0lyRJkiSpZ2CO1EgOQqtFOH3CqwhzJJ/3AkBrunfvXlJS4vV69S6kveLi4gjCNpOESLCIWo9I6bgJCgAg5PTt21fvEkJahJ8rMe82AKB1kR6EFq6XAQC0JsKDMCmaebcBAK2J8CCkaxQA0LpID0K6RgEArYr0IOSeegBAqyI8CBkjBAC0LsKDkDFCAEDrIj0IGSMEALQq0oOQMUIAQKsiPAiTopl3GwDQmggPQrpGAQCti/QgpGsUANCqCA9CukYBAK2L8CCkaxQA0LoID8JYs1CFcPn0rgMAEKoiPAiFEDZOCgEALYv8IGRyGQBAK4wQhKLGo3cRAIBQFflBmBQlql16FwEACFWRH4R0jQIAWmGEIKRrFADQosgPQrpGAQCtiPwgpGsUANCKQIPQ6XRu2LBh5cqVBw8ebNyuKMrGjRs/+uijmpqaxu2VlZUrVqz45ptvmuxny5YtK1asOHbsWOPG6urqjz76qKSkRFGUtj+F06BrFADQioCCcPfu3Z06dZo5c+azzz573nnnzZs3T2v3+XyFhYXTp09//vnne/Xq9cMPP2jtxcXFffr0WbBgweTJk//0pz/59zN16tTrrrtuwYIFffr0+eqrr7TGbdu29erV6/nnn7/llluuvvpqny/I08Aw7zYAoDVqACorK3/66Sdtef369bIsV1VVqaq6bNmynj17NjQ0qKo6e/bs8ePHa9sMHjz4ySefVFX1+PHjGRkZGzZs0P4wIyPj+PHjqqoWFRUNHTpU2/iaa66ZM2eOqqoNDQ15eXkffPBBS2WMGTNm+fLl/tXa2tpAin9vn2/cKm8gW0a8AI8Y/Orr630+n95VhBOHw+F2u/WuIpy43W6Hw6F3FeHE5/PV19cHd58BnRGmpqbm5ORoy7169VIUxW63CyHef//98ePHx8bGCiGmTJnywQcfeL3eo0ePrl+/fsqUKUKIpKSkq6666v3339c2Hj16dFJSkrbxunXrKisrPR7PihUrrr/+eiFEbGzsuHHjtI2DKDFKqnYxRggAaJ65rX/wxBNPDB06NDs7Wwhx8ODBAQMGaO3dunXzer1HjhwpLy+Pi4tLS0vT2rt27bpnzx4hxC+//NK7d2+t8ayzzoqNjT148KDD4fB6vV27dvVvvHXr1pYeura2tri4uKGhQVuVZXn8+PGnLTjBrNZ4xJkYfQw7iqJwHNqEI9ZWiqJIksRBC5zyK70LCRttPWKSJEmS1Po2bQvCN99885VXXvnyyy+1VbfbbbFYtGVtweVyuVwus/m33VosFqfT2WRjIURUVJS2sRDCv71/42YdP3583bp1WqwKIWw221VXXXXammNUUe1qbbfG4XK5Gv8vwGk5nU5JkmQ58i+uDhan02kymYI+0h/BPB4Ph6tNFEXRXmYBbh8VFdU4kprVhiB877337r777k8//dTfTZqZmVlZWaktV1RUCCE6deoky3JdXZ3L5YqOjtbaO3XqpP2Tf2OXy1VbW9upU6eUlBQhRGVlZefOnbWNtYVmdevWbcaMGWPGjNFW6+rq4uLiTlt2J1nUuD2BbBnxfD4fx6FNVFWNjY0lCAMny7LJZOL7VuC0IIyJidG7kLCh9ToE96Ms0Hf4ypUr77jjjuXLl5933nn+xkGDBq1du1ZbXrt27UUXXRQXF9e1a9fs7OwvvvhCa//8888HDx7cZOPPP/9c28xqtfbt27fxTgYNGhSUJ+aXGCVquX0CANCCgM4Id+7cOW7cuOHDhy9ZsmTJkiVCiNtvv7179+5/+tOf/vnPf/7lL3/p06fP7Nmzi4qKhBAmk+mee+654447HnjggeLiYrvdro3kTZgw4cEHH5wxY8bgwYMfeuihe+65R/uiPWvWrLvvvtvtdm/fvn3nzp3vvvtucJ+hRRYWWdi9Ir7N46EAgMgXUDjEx8c/9NBDJ/2Z2SyESE5O3rhx43PPPff1118vWLDAP2J35513ZmZmfvbZZxkZGcXFxVofaXR0dHFx8TPPPFNcXPzYY49NnDhR23jKlCmJiYnLly9PSkrauHGjdllpcCVGiRq3Gm8+zXgpAMCAJFUNm1sLCgsLm4wRJiQkBPKHfZZ63xtl6pNk9CAM/IhBY7fbGSNsE+0qBsYIA8cYYVspiuJwOOLj44O4T0O8w5OYXAYA0AJDBGFilKgmCAEAzTFIEPIDFACA5hkkCOkaBQA0zxBBmETXKACgBYYIQrpGAQAtMUgQ0jUKAGgeQQgAMDRDBGFSlKimaxQA0BxDBGFilMQZIQCgWQYJQrpGAQDNM0QQcvsEAKAlhghCbp8AALTEEEFos4h6j1CIQgDAKQwRhLIk4s2ijt+pBwCcwhBBKOgdBQC0wDhByIWjAIBmGCUIk6K5cBQA0AyjBGGihTNCAEAzDBOEjBECAJpjlCCkaxQA0CyjBCFdowCAZhkmCOkaBQA0xzhByBkhAKAZRgpCZpYBAJzCKEGYFCVVu+gaBQA0ZZQgpGsUANAsghAAYGhGCUJ+mxcA0CyjBCG3TwAAmmWUIIy3CJciPIredQAAQoxRglASwmYRtdxBAQA4mVGCUNA7CgBojqGCkAtHAQBNGSgIuXAUAHAqAwUhXaMAgFMZKgjpGgUANGWgIEyKEtUuvYsAAIQYAwUhP0ABADiVoYKQMUIAQFOGCkLGCAEATRkoCLl9AgBwKgMFIV2jAIBTGSoI6RoFADRloCBMIggBAKcwUBAmRknVdI0CAE5mqCDkjBAA0JSBgjDGJIQQTp/edQAAQomBglBwUggAOIXRgpA7KAAAJzFWEHJPPQCgCWMFIV2jAIAmjBaEdI0CAE5irCDknnoAQBPGCsJExggBACczWhDSNQoAOInBgtBC1ygA4CQGC0LGCAEAJzNWECZFM0YIADiJsYKQMUIAQBMGC0LGCAEAJzNWENI1CgBowlhBSNcoAKCJjg5Cp9NZUVFxantlZaXD4TjTj55oEbUeQRICAPwCDcIHH3ywX79+KSkpCxYs8DfOnDkzpZE+ffpo7aNHj/Y3Dh482L99UVFRZmbmhRde2K9fv/3792uN5eXlgwcPPv/88zt16vTQQw8F6Xk1zyyLaFnYPWf0QQAA4STQIOzZs+e8efN69+7tdDr9jUVFRVW/uuqqq0aNGqW119XVvfzyy1r7+vXrtcY9e/Y88MADJSUlhw4duuyyy+655x6t/f777+/Ro0dZWdkPP/zwv//7vyUlJcF7ds2gdxQA0Jg5wO2mTp0qhHjkkUea/deampr333+/uLi4lT288cYbV155Za9evYQQM2fO7NWrV11dXXx8/OLFi9esWSOEyM7OnjBhwqJFiy655JK2PYm2SIwSNR6RdeYeAAAQVoIzRrh48eK8vLyLLrrI3zJlypTY2NgLLrjg448/1lpKS0u1FBRC9OjRQ5KkX375pby8vL6+3t/eq1evffv2tfQoXq/3yJEjpb86cuTI7yg1KUpUu37H3wEAIlOgZ4StW7BgwfTp0/2rTz75ZJ8+fUwm06uvvnrttdd+//33OTk5NTU1cXFx/m3i4+Orq6vNZrMQwt+uNbb0KHv37p0zZ05sbKy2mpiY2Po5aLOsJsuRGnd9vNLWP4wAdrtdkiS9qwgnDQ0NPp9Plo11cXV7OJ1Ok8lksVj0LiRseDwen8/n9Xr1LiRsKIridDpVNdARrpiYGC1oWhGEIPz+++9/+OGHyZMn+1v8p4a33nrrggUL1qxZk5OTk56e7g85RVFqa2szMjJSUlKEEDU1NampqUKI6urq9PT0lh6od+/eM2bMGDNmjLZaV1dntVrbWm1KrM9tslitRvxoU1X1dxwxI5MkKTY2liAMnNlsJgjbRAvCmJgYvQsJG4qimEym+Pj4IO4zCO/wl156afz48Wlpac3+a11dnfb/+Pzzz9+0aZPWuGXLloSEhOzs7OTk5OzsbH/7pk2bzj///PaX1IokfpIQANBIoGeEJSUlBw4cKC8v//bbb5csWTJw4MDs7GwhhNvtXrx48ZtvvunfsqKi4rXXXhs6dKjZbF64cGFFRcUVV1whhJg6deqDDz740ksvDRky5K9//estt9wSHR0thPjzn/88Z86cLl26bN++/ZNPPpk3b94ZeJq/4QcoAACNtSEIi4uL+/TpU1tbu2TJki5dumhBuHPnzjFjxowYMcK/pcVi2bZt21tvvaUoSt++fb/88suzzjpLCJGSkrJy5coHH3zw6aefvvzyy+fOnattP2vWLJfLdf3116empr733nvdunUL9nM8CbdPAAAakwIfctRdYWFhkzHChISEtu7kmR3Kjmr1mUGmYFcXBn7fETMyu93OGGGbcLFMWzFG2FaKojgcjpAbIwwv3D4BAGjMcEFI1ygAoDEDBqGoYa5RAMCvDBmEXDUKAPiV4YKQMUIAQGOGC0LGCAEAjRkuCBMsosEnfEQhAEAIYcAglCURbxZ1XC8DABBCGDAIBb2jAIBGjBiEzLsNAPAzYhByBwUAwM+wQUjXKABACKMGocQZIQBAY8QgZIwQAOBnxCBkjBAA4GfMIOT2CQDACcYMQs4IAQAnGDEIGSMEAPgZMQjpGgUA+BkzCOkaBQCcYMQgpGsUAOBnxCBkZhkAgJ8xg5CZZQAAJxgxCOPNwqMIt6J3HQCAEGDEIBRC2CyilpNCAIBhg5A7KAAAGoMGYVI0F44CAIQwbBAmWriVEAAghHGDkK5RAIAQwrBBSNcoAEBj0CBkljUAgMaoQcgYIQBACGHYIOwcLx20M0YIADBqEObZpL21BCEAwLhBKPbW6l0EACAEGDQIu8RLx1yqw6t3HQAAvRk0CGVJdLdKP9XROwoARmfQIBRC5NoEw4QAAOMGYZ5N+olhQgAwPOMGYa5N+okzQgAwPOMGIXdQAACEsYOQOygAAAYOwu5W6XCD6vLpXQcAQFfGDUKzLLrESz/X0zsKAIZm3CAU9I4CAAwfhNLeGs4IAcDQDB2EuTYmlwEAozN0EHIHBQDA4EHIGCEAGJ2hg7BHgnSgXvUoetcBANCPoYMw2iQ6xUkHuIMCAAzM0EEo6B0FAMMzehDmJnC9DAAYmuGDkN+gAABjM3oQ5tnET3V6FwEA0A9ByOQyAGBoBKG0r15ViEIAMCqjB2GsWaRESwftJCEAGJTRg1BwBwUAGBtByIyjAGBoBCF3UACAoRGEdI0CgKERhHSNAoChEYSip036qVYlCQHAmAINwiVLltx+++0FBQUffPCBv3Hjxo0FjXz99ddae319/W233Zabmztw4MC1a9f6t1+9evWll16am5s7Y8YMu92uNbrd7rvuuisvLy8/P7/xzjuM1SKsFlHWQBQCgBGZA9xuy5Yt3bt3Ly4u/uWXX/yNFRUVZWVlRUVF2mpOTo62cO+99x46dGjdunVffvnluHHj9u7dm5aWVlFRMX78+BdffHHIkCHTp0+/7777nnrqKSHEI4888vXXX69evXrHjh2TJk369ttvc3Nzg/ocTy/PJu2tFZ3jOvhhAQD6CzQIH3vsMSHEqlWrmrQnJSWNGjWqcYvD4Xj11Ve//PLLrKysyZMnL1y4cNGiRTNnznz99dcvueSSSZMmCSEefvjh4cOH//d//3d0dPQLL7zw+uuvd+vWrVu3bmPGjFm4cOHDDz8cjKfWBrk2aW+tOjRT6uDHBQDorr1jhLt37x44cODo0aNfe+01VVWFEPv373c6nX379tU26Nev344dO4QQ27dvz8/P1xr79u1rt9t/+eWXqqqqI0eO9O/fv8nGHYw7KADAsAI9I2xWbm7uM888k5ubu2vXrnvuuaeqqmrmzJnHjh2zWq2SdOLsKikpadeuXUKIysrKvLw8rVGWZZvNVlFRoWVnQkKC1p6YmFheXt7Sw+3YsWPy5MkWi0VbTU1N3bp1a3vq98uOkj86ZKqvdwRlbyHLbrf7/78gEA0NDT6fT5a5pixQTqfTZDL536Q4LY/H4/P5vF6v3oWEDUVRnE6nGvAFjjExMWbzaZKuXUF4zjnnnHPOOUKI/v37u93uoqKimTNnJiUl1dfXq6qqfebW1tampqYKIZKTk+vr67U/VFW1vr4+JSUlKSlJCFFfX68t1NXVaRs3q3fv3g8//PAVV1yhrbpcLqvV2p76/c5LV5/d47Nao4Oyt5ClqmqwjphBSJIUGxtLEAbObDYThG2iBWFMTIzehYQNRVFMJlN8fHwQ9xm0d3hycrLD4RBCdOnSRZKkn376SWvftWuXdhFNbm7uzp07tca9e/dKkpSdnZ2Wlmaz2bRTxsYbN1+rLFut1uRfBfFA9LRJe/gxJgAwpECD0G63Hz9+3OPxOByO48ePu91uIURxcfGxY8eEEAcOHHjkkUe0czWbzTZu3Lh58+YpirJ169ZVq1Zdf/31QoipU6euWrVq69atiqLMmzfv2muvtVqtsixPmzbt8ccf93q9e/fufffdd2+44YYz9mRblBwtzLKocHb8IwMAdBZoEM6ZMyc/P//QoUPPPfdcfn7+6tWrhRCffPJJ9+7drVZr3759L7rookcffVTb+H/+53927dqVkpJSUFDw1FNPde/eXQjRo0ePoqKiUaNGpaSk/Pjjj/Pnz9c2njt3rt1uT0tLu/jii+fMmdOvX7/gP8sAML8MABiTFPiQY0scDkdsbGw7251OZ3R0dOuXchQWFs6YMWPMmDHaal1dnf8qm/a7fq3vyi7StLxIHg0K7hEzArvdzhhhm3CxTFsxRthWiqI4HI7gjhG262IZTbNp19Z23V8HeTbBHRQAYEB81T0h1ybxGxQAYEAE4QmMEQKAMRGEJ+TZpL3cQQEAxkMQnpAeK7yqOO7Suw4AQMciCH+TS+8oABgPQfgbhgkBwIAIwt/kJgguHAUAoyEIf8OPMQGAARGEv8mzST/VEYQAYCwE4W/ybII7KADAaAjC33SOl+o8os6jdx0AgA5EEP5GEiKHYUIAMBiC8CTcQQEARkMQniTPxh0UAGAsBOFJuIMCAIyGIDwJXaMAYDQE4UnoGgUAoyEIT9IlXjrmUh1evesAAHQUgvAksiS6W5lfBgAMhCBsimFCADAUgrAphgkBwFAIwqa4gwIADIUgbIogBABDIQibomsUAAyFIGyqu1U63KC6fHrXAQDoEARhU2ZZdImXfq6ndxQADIEgbAa9owBgHARhM/JsEj9VDwAGQRA2I9fG5DIAYBQEYTPOTpR2VROEAGAIBGEzBqZLJeVcOAoAhkAQNiM5WvROkr4q56QQACIfQdi8y7OlVQcVvasAAJxxBGHzCrLkTw9xRggAkY8gbN6gdOmnWrXCqXcdAIAzjCBsnlkWl2XKaw7TOwoAEY4gbFFBlrSK3lEAiHQEYYsKsqRPDhKEABDhCMIW9U6SzLLgznoAiGwEYWtGdpa4dhQAIhtB2JqCLGnVIa6XAYBIRhC2ZlSW/EUZc60BQCQjCFuTGi16JUobmWsNACIXQXgal2fTOwoAkYwgPI2CLJm7CQEgghGEpzEoQ9pVrVa59K4DAHBmEISnESWLIZnSauZaA4AIRRCeHr2jABDBCMLTuzxb+pS51gAgQhGEp3dOkqSo4scashAAIhBBGJBR/BIFAEQogjAg/CQTAEQqgjAgBVny52WKh0tHASDiEIQBSYsRuTaphLnWACDiEISB4pcoACAiEYSBKsiS+W1CAIg8BGGgLsuUdhxnrjUAiDQEYaCiZDEoQ/q8jN5RAIgoBGEbMNcaAEQegrANLs+WPmauNQCILARhG5ybLLl94qdashAAIgdB2AaSEKOyJK4dBYBIQhC2DXOtAUCEMQe43Y4dOzZu3Lhnz56xY8cOHDhQa9yzZ88777yzc+fOhISESZMmDR8+XGt//vnnf/75Z205IyPjrrvu0pb3798/f/78ioqKUaNG3XzzzZIkae1vvvnm8uXLk5KSZs6cefbZZwfpqZ0RBVnyf37l8SgmC18hACAiBPpxPmvWrBUrVrzxxhvffPONv3H+/Pnl5eVXXXVVXl5eYWHhsmXLtPbFixeXlZUlJycnJyfbbDat0eFwXHbZZRaLZeLEiY8//vj8+fO19oULF957773jxo1LT08fMmTIsWPHgvfsgi8jVnS3SpsqOCkEgAgR6BnhihUrhBAjR45s3Pjcc8/5z+oqKirefPPNa665RlsdP3781Vdf3Xjjt99+Oz09Xcu/pKSkG2644a677jKZTPPnz583b97EiRMnTpy4cePGV1555e67727nszqjLs+WVh1SB2VIehcCAAiCdnXw+VNQCHH48OGMjAz/6htvvHHnnXf+61//crlOzMXy9ddfDxs2TFseMmTI4cOHDx48WFdXt2PHjqFDh2rtQ4cOLSkpaU9JHaAgS2bSUQCIGIGeEbZu9erVy5cv/+6777TVYcOGJSUlmc3mBQsWvPjii+vXr4+KiiorK8vPz9c2iIqKstlsZWVlHo9HCJGamqq1p6WlHTlypKVHKS0tnT179hNPPKGtJiYmvv7660Gpv00usoofa6K2ltnzEsKsg9Rutzf+7oLTamho8Pl8ssyAcKCcTqfJZLJYLHoXEjY8Ho/P5/N6vXoXEjYURXE6naoa6MdvTEyM2XyapAtCEG7atGnKlCnvvPNO165dtZa5c+dqC7fddlvPnj0//PDDcePGxcXFud1u/1+5XK64uLjY2FhtWXvnOJ3OuLi4lh4oIyNjxIgRl1xyib+llY3PnDgh/vMc9ek98r+GhNnno8/n0+WIhS9VVWNjYwnCwMmyTBC2iRaEMTExehcSNhRFkSQp8I+yQN6/7Q3CrVu3jh079qWXXiooKDj1X2NjY/Py8srKyoQQ2dnZ+/fv19qPHj3qcrmys7NtNpvFYtm/f/+5554rhNi/f392dnZLjxUfH3/hhRf6H6iurk6vT6j/c57Ie8fzQIPUzRpOJ1iyLPOZ3ibyr/QuJGxwxNpKlmVVVTlibRL011i79vX9999feeWVRUVFhYWF/kan01ldXe3fYNOmTVqP6IQJEz788MPKykohxKuvvvqHP/whJSXFbDZfffXVr776qhCirq7uvffemzBhQntK6hiJUeKWXvIT3zNSCADhTw3MzJkzc3JyYmNj09LScnJyPvroI1VVx4wZExsbm/OrCRMmqKq6f//++Pj4Sy+9dNCgQVarde7cuf6d3H777dnZ2SNHjszMzNy8ebPWuHPnzuzs7OHDh/fo0WPy5Mk+n6+lGsaMGbN8+XL/am1tbYDFnwlHGtTk19xlDTqW0Gb6HrFwVF9f38oLEqdyOBxut1vvKsKJ2+12OBx6VxFOfD5ffX19cPcpqYENOZaXl9fX1/tXMzIy4uPjjxw50tDQ4G+MiYnp3LmzEKKmpmbHjh1CiN69eycnJzfez969e8vKyvr16xcfH+9vdDqdmzdvTk1N7dOnTys1FBYWzpgxY8yYMdpqXV1dQkJCIMWfIf/1lS/BIh4ZYNKxhjbR/YiFHbvdzhhhm3CxTFsxRthWiqI4HI7GCdJ+gY4Rpqenp6enN2nMzMxsduPExET/7DNN5OXl5eXlNWmMiYkZMmRIgJWEjll95X7ve2f1NSVF6V0KAOD34qvu79clXhrTVX5mByOFABDGCMJ2mfG9ij4AAA97SURBVH2h/NR2X71H7zoAAL8XQdgueTZpWKb80m5OCgEgXBGE7TX7Inn+94qbKASA8EQQtlffFOmCFPHaHpIQAMISQRgE919kenSr4iUKASAMEYRBcGm61MUqluwjCQEg/BCEwXFfX9M/vlWUMPs5CgAAQRgkV2RL8Rax4gAnhQAQZgjCoLm3r/zwVoIQAMIMQRg047rLDV6x+jDdowAQTgjCoJGEmNVXfnSrT+9CAABtQBAG05RceX+9+OooJ4UAEDYIwmAySeL/ni//cxsjhQAQNgjCILvpbPmbSnUNI4UAECYIwiCLMYnXhpuuW+vdV0cWAkAYIAiD7w+dpFkXmMZ/5nN49S4FAHA6BOEZcff58gUp0oxiriAFgFBHEJ4pzw82ba9Wn9vJhTMAENIIwjMl1izeHWn6f9/4vjzCYCEAhC6C8AzqniC9Ntx83VrfITtZCAAhiiA8swqypNt7yxNX+/gJewAITQThGTf7IjkrXvq/G7lwBgBCEUF4xklCLBxqWntYXfgjZ4UAEHIIwo6QYBHvFZju2+TbVMFgIQCEFoKwg/RKlF4cYpq42lfh1LsUAEAjBGHHubqbPCVXmrzG66WLFABCBkHYoebmm2JM4o9rfLUevUsBAAghCMIOZpLEvwvMZyeKAcu8248zXggA+iMIO5pZFo8NMD1wkTziI+/SfXSSAoDOzHoXYFDX58nnJkvXfubbXKk+nG8ySXoXBABGxRmhbi5Mlb6+xvxtpTrqI2+5Q+9qAMCoCEI9pUaLj640D+skXfxvL7cYAoAuCEKdmSTx936mpwbKhZ96mXoGADoeY4QhYWw3uVeSNH6Vr6RcfXqQKYrvJwDQUfjEDRW9EqWvxpqPOkT+Mu+y/Qr9pADQMQjCEJIYJZYVmJ64xPTIVqXve94l+4hDADjjCMKQMypL+vpq8xOXmB77Thn4gXf5AQYOAeAMIghD1KgsafM15tkXyg9sUQYRhwBwxhCEoUsSorCrvOUa813ny38pUYYs9645TF8pAAQZV42GOlkSE3vI47rJb/ykTP/S19Uqrs+Tx3WXU6P1rgwAIgJnhOHBLIs/9ZR3TTT/17nyZ4fU3Lc9V37sXbBbqXLpXRkAhDnOCMOJRRbju8vjuwunz7TqkLKkVL2nxHNusjSxh3xdrpweq3d9ABCGCMKwFGMShV3lwq7C7jV9eEB5Z5/64DeewRnStT3kYZlSro05vAEgUARheIs3iz/myH/MEXavafl+Zdl+9cEtikdRL02XB6ZLgzKk/mlSHP+TAaBlfEZGiHizmJwrT84VQoiDdvWro+qGcnXW18q2KvXcZGlgunRpunRBvHRugt6FAkCIIQgjUHa89Mcc6Y85Qgjh9IktleqGcnXJPvXuo9H1Xs/ZidLZiVLvJKlXotCW43kVADAwPgIjXIxJDM6QBmdI4nxRV9fgi07YXa3urlF316hL94kfa5Qfa9SzYqSzE8XZiVKOTeoUK7Ljpcw40SWePlUAhsBHnbEkRYlL0qVL0n+7mkZRxYF69cda8WON+nOduvWY+KVeOeIQB+2qRRZZcVLneNE5TsqKExmxUnK0SImWUqJFyq8LZm7AARDmCEKjkyXRPUHqniAuz2p6rWmNWxxqUA/ZRVmDetAuSuvUqkpR5VKqXKLKJapcapVLWM0iJUZKjRbJ0cJmkeLNIt4iEqNEgrZsFknRwmqW4szCahFRsojX/muWok2CM04AoYCPIrQoMUokRknnJAkhWrwfo9otqlzqMaeocYsat2r3CrtX1LpFrUetcIoGr6h2iXqv0uAV9R7hVoTdI1yKaPCqLp9o8IoYk4g1ixiTiDVJQghblDBJQhIiKVoIISyysJolIUSUSfgHMi2ysFpOLJslkWD5rTaTLGyWk8rTdn4qq1mytHwu6w9ph0OKjlZl+TQz2yV31Cw/JknYLCF9b4zLLcmysDAzbsA8HqEoUnRI/18NLYoiZEXEB3WfBCHaJSlKJEVJOScuRm3zu9npEw6vcPhUp08IIWrcQlGFoooatxBCuBVh96pCCLdP2L0n/sSjiHrPiWWvKo67f0spnyL2eBrvXji8QttzE/VexdPyh7UW0kIIRTHLsipEc7to5HhHze/jU0WtJ6Tnm1VVSZKEEJ7Tbgk/VZUkiSPWBjflmuYPCuYOCULoKcYkYkwiubUE1fOrst1uj42NlWUGQgPldDpNJpPFYjn9phBCCOHxeHw+JSYmRu9CwoaiKA6HO7j75B0OADA0ghAAYGgEIQDA0MI4CIuKitzuIPcUR7YFCxZUVlbqXUU4Wbp06d69e/WuIpx89tlnmzZt0ruKcLJly5ZVq1bpXUU4KS0tXbJkSXD3GcZB+Mwzz1RVVeldRThZtGhRaWmp3lWEk2XLln333Xd6VxFOVq1atX79er2rCCfr16//9NNP9a4inHz33Xfvv/9+cPcZxkEIAED7EYQAAEMjCAEAhiapakhPVNFY7969a2trY2NjtdUDBw5kZ2dzs3PgDh8+nJaWFhUVpXchYaO8vNxqtcbFxeldSNioqqoymUyJiYl6FxI2amtrPR5Pamqq3oWEjYaGhvr6+vT09AC3nzJlyty5c1vfJpyC8NixY9XV1doMTkIIl8sVHd1RkzxGBI5YW7ndbovF4n/J4bS8Xq8kSSZTcxO8ojk+n09VVbOZSb4Cpaqqx+MJ/At9p06d/KdPLQmnIAQAIOjoVwQAGBpBCAAwNIIQAGBoBCEAwNDC8lIlVVXffPPNzZs35+Tk3HrrrfyUV7Oqqqq2bNmyb9++Sy+99IILLvC3l5eXL1y4sKKiYvTo0SNGjNCxwlCzf//+FStWlJaWZmZmTps2LTMzU2tXFGXRokVbt27t2bPnLbfcwv0nfqWlpStWrPj5559tNtvo0aMHDBjg/6d33313/fr1Xbt2nT59enx8cH9OPBLs3bt3zZo1o0ePzsrK0lq+/fbbt956Kzo6+sYbb8zJydG3vNBRU1Pz9ttv+1cHDhx4/vnna8slJSVLly61Wq033XRT165d2/MoYXlG+Le//e3RRx/t2bPnBx98MH78eL3LCVHXXHPN3/72t7///e+Np/RtaGgYOHDg7t27u3XrNmXKlMavMEyaNGnLli1dunTZvn37Oeecs2/fPq39rrvuKioq6tmz5zvvvHP99dfrW2RIWbduXWlpaY8ePdxu98iRI5cuXaq1P/bYY/fdd19eXt7atWuvuOIKLk1vwuv1Tps2bebMmTt37tRaSkpKhg0blpyc7Ha7BwwYcODAAX0rDB1Hjx6dOXNm6a9qamq0du2llZmZWV1dPWDAgPLy8nY9jBpuampqrFbrDz/8oKqq3W5PTEz85ptv9C4qFGn3J40ePfrxxx/3Ny5YsGDAgAGKoqiqunjx4gsuuEC3+kKPw+HwL48YMeLhhx9WVbWioiImJqa0tFRV1Zqamri4uF27dulWYgh74IEHxo8fr6qq0+lMS0tbv369qqput7tz585r167VubgQ88gjj8yePTsrK2vVqlVay4QJEx544AFteerUqffee69+1YWW3bt3p6Wlndp+xRVXzJs3T1u++uqrtXfr7xZ+Z4Rbtmyx2WznnnuuECIuLm7w4MFffPGF3kWFombn3Fm3bl1BQYF2h3hBQcG2bduOHz/e4aWFqMZ97C6XKyEhQQhRUlLSpUuXHj16CCFsNtvFF1+8bt063UoMVQ6HY/PmzVqf1fbt210u18CBA4UQFotl+PDhvEMb271791tvvTVnzpzGjevWrbv88su15YKCAo5YYy6Xq6io6Nlnn921a5e/8YsvvgjiEQu/IDxy5MhZZ53lX83IyDh8+LCO9YSXsrIy/9FLTU01m81lZWX6lhSC3nrrrX379k2bNk3wejudDRs29OjRIyUlRVEU7cP9yJEjaWlp/ul4OGKNKYoyffr0J554ovG3Lo/HU1lZ6X+Zpaen8670M5vNw4YNq6ys3LRpU35+/uLFi4UQx48fdzqdQTxi4ReEZrPZ5/P5Vz0eD9OGBc5sNnu9Xm3Z5/P5fD4u/WiiuLj4zjvvfPvtt5OSkkRzrzeOWGMDBgzYsmVLcXFxTU3NvffeK05+jQneoSd78skn+/TpM3LkyMaNJpNJlmX/QfN6vbzG/HJycpYvX/6Pf/zj5ZdffvbZZ++55x4hhMViEUIE8YiFXxB27ty5rKxMURRt9dChQ506ddK3pDCSlZXl/3quLXD0GisuLr722muXLFkyZMgQraVz586HDh3yb3Do0KHOnTvrVF0oMpvNKSkp/fv3v//++7Vrrzp37lxRUeF2u7UNeIc2tnjx4nXr1uXn5+fn55eXl99xxx0vvviiLMuZmZn+lxmvsZYMHjy4rKzM4XBYrVabzdb4iLXzNRZ+QXjxxRdbLJa1a9cKIQ4fPrxx48bRo0frXVTYKCwsXL58udPpFEIsXbp0xIgRXNrut3HjxnHjxr3yyivDhg3zNw4ZMsRut2/YsEEIsW/fvm3btl155ZX61RhaGhoa/MubN2/WLmE/55xzsrKyVqxYIYSoqqpas2ZNYWGhbiWGmNdee23RokUvvPDCCy+8kJycfNddd40ZM0YIUVhYuGTJEiGEqqpLly7liPlpH1aa5cuX5+bmajNojx07VjtiiqK89957Y8eObdfDtOdKG728/PLL6enpN910U05Ozt133613OSHqoYce6t+/f2JiYnZ2dv/+/T/66CNVVb1eb0FBQf/+/adNm5aamvrVV1/pXWYIyc3N1U5uNP/4xz+09meffTYjI+Pmm2/u1q2b/9I+qKp6+eWXjxgxYtq0aUOGDMnIyNiwYYPWvnTp0rS0tBtvvLFXr1633nqrvkWGrMZXjZaWlnbu3Pnaa68dOXLkeeedV11drW9toeP+++/Pz8+fOnXq8OHDU1JSPv30U619586d6enpkyZNGjp0aH5+vt1ub8+jhOuvT+zevfubb77Jzc29+OKL9a4lRO3fv7+ystK/ql3RIITw+Xyff/55ZWXlsGHD/PeMQwjx/fff+zv0hBBpaWndunXTlnfs2PHdd9+dffbZ/fv316m6UGS320tKSo4ePZqenj5w4MDGP9xYWlpaUlLStWvXwYMH61hhKNu2bVv37t1tNpu2WlNT89lnn0VHR48aNYpJQvwcDsemTZsOHTqUmpp68cUXayP3mqqqqtWrVyckJIwYMaKdY4ThGoQAAARF+I0RAgAQRAQhAMDQCEIAgKERhAAAQyMIAQCGRhACAAyNIAQAGBpBCAAwNIIQAGBoBCEAwNAIQgCAof1/ye1XNulvviUAAAAASUVORK5CYII=", - "image/svg+xml": [ - "\n", - "\n", - "\n", - " \n", - " \n", - " \n", - "\n", - "\n", - "\n", - " \n", - " \n", - " \n", - "\n", - "\n", - "\n", - " \n", - " \n", - " \n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n" - ], - "text/html": [ - "\n", - "\n", - "\n", - " \n", - " \n", - " \n", - "\n", - "\n", - "\n", - " \n", - " \n", - " \n", - "\n", - "\n", - "\n", - " \n", - " \n", - " \n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "plot(results.free_energy[1:end])" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Free energy is converging to a stable value, indicating that the inference algorithm is converging. Let's visualize the posterior distribution and how it compares to the true parameters." - ] - }, - { - "cell_type": "code", - "execution_count": 83, - "metadata": {}, - "outputs": [ - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdeXxTZaI38HOy72nTNE3SPd2hpS2llL1QoCwDyqDgMsogOio4jHrVV1yur3ccF0bH8c5cN1DBBdHr4DggvkrBQQoUWcrSFrrQfUnXNG32/f0jTqdSKC1Ne3KS3/cPPu1Jmvya0vx6znnO85Aej4cAAAAIVgyqAwAAAFAJRQgAAEENRQgAAEENRQgAAEENRQgAAEENRQgAAEENRQgAAEENRQgAAEENRQgAAEENRQgAAEHNX4rQ4/E89thjVKegmMPhoDoCXIfT6aQ6AgwHPyD/54dvdKSfzDXqcDiEQqHdbqc6CJUMBoNYLKY6BVyTx+Mxm81CoZDqIHBNJpOJz+czGP7yJz4M5YdvdPjvAgAAQQ1FCAAAQQ1FCAAAQQ1FCAAAQQ1FCAAAQQ1FCAAAQQ1FCAAAQQ1FCAAAQY1FdQAA6rlcLpvN5p3w4lr/ulwuBoNhtVolEglBEGw2myAIFotFkiSTyWQwGCwWi8vlslj4nQKgGfzSQlBwu90Wi8Vms3n/NZutJpPNZPrpX4fDw2RymUy2202wWGyXi2Cxfvaxy0UwmUyXy26329lsPZNJOJ0OBoNwuZwk6XG7nQThcbudLpeVxSL4fK5AwBUKuUIhTyDgcP+Fx+Px+XyqXwkAuBKKEAKQx+MxmUwGg6G/36DTGbq7Df39FpLkMZlcJpNPkhwWi8fhSLlcXkgINyKCx2SO6BfB4/FYrdbhy8ztdtntNrvdZjLZdDqr02l3u/vdbpvbbXM4zEymMzRUFBIiDA0ViURCkUgkFAqxEzkB/vjHP54+fZrqFEAQBOFyuZhM5g184YoVK9atW+fzPATmGvUrfjgFH104nU6dTqfX9+l0Bp3OoNebSZLHYolZLDGfLxYKxTyecOzzT46kCIfncjktFpPZbLRajU6nyeEwOp0moZAdEiIKCRGGhYlDQkIkEgmmyrxh15prdM6cOUuXLk1JSaEkFYzdkSNH+vv7P/zww/F4cPwpCnRlt9t1Ol1XV09rq66nx8Rmh7BYIXx+hEiUGB4u8s8uYTJZIpFUJJIO3mizWcxmY1ubqa6u3+lscrtNcrlYpQqRyUJCQkJEIhFVaQPMggULZs+eTXUKuEEWi+XQoUPj9OAoQqATq9U6UH56vZXNDuVyw8Ti9JSUEJIkqU53g7hcPpfLDw0N937qdruMxr6GBv2lS512ezWb7VAqQyIifupFLpdLbVqAwIMiBBowGAxtbdq6uvaeHiuXK+PxwiSSGIVCQt/yGwaDwZRIZBKJzPupw2E3GPSVlXqHo8Fu18tkvOjocKUyXCaT3diJFgC4AooQ/JfBYGht1dbWtun1Lh5PGRIyKTU1LCDLbxhsNkcmU8hkCuKnQUD99fXdFRW1TudpuVwUHx+hVEZIJIH5NwHAxEARgn/xeDw6na61VVtX12G1MrlclUyWrVRKr/+VQYAkyX+dYkxwuZx9fT1lZd2nTp1ls+3R0fLo6HC5XI4rNABGC0UI/sJisTQ0NF282GS38wQCVVhYnkCAcSLXxGSyZLIImSyCIAibzaLXd7e2dtlsl+RyfmKiWq1WCQQCqjMC0AOKECjm8Xh6enpqahrr67vZbJVSOUMgwDUko8Pl8iMiogki2uPxGAy9Fy5oT5w4JpNxEhJUkZFqjDu9qrvv/vXXX++bmOdisViffrpr8eLFE/N0MFooQqCMzWZrbm6pqGg0m9kSSWxSUhaDgdEfY0KSpHegjcczyWDoLS/Xnj59IjSUjUYcqrGx6bnnPszKmjMBz7V164Pt7e1XbNy5c+eRI0f0en1MTMzGjRtHdY1jaWnphx9+2NzcrFAoFi5ceMstt4zqeiGHw/HDDz8sWrRo5F8S2FCEQIGenp7q6ob6+m4ORx0enhMVhVOAPnZFI5aVtXkbMTFRHRUViaOmXkKhWCIJnYAn4nCuctFLQ0NDfn5+WFjY4cOHZ8yYUVFRoVarR/JoO3bs2Lx586OPPnrPPfd0dHS88847nZ2dDz300MjzGI3GxYsXu91ujLHyQhHChOru7r5woVqrtUul8YmJmSOc2wxu2KBGnGww9JaVaU+ePBodLU1KilYqlf457UCQeP75570frFixYu/evSdOnFi9ejVBEDqdrri4uKOjw3trXFxcYWHhwFe1t7dv2rTps88+u/nmm71b7r333q6uLu/H3377bUlJiUqluuOOO6RSKUEQxcXFYrG4ra3t2LFjWVlZa9asIQjiq6++Ighi27ZtJEmuXr1aLpc3Njbu2bPHarUWFhZOmzaNIIi2traBSARB7Nq1a/ny5aGhoUVFRVFRUeXl5WfOnHnqqadMJtOnn37a3d0dHR19yy23KJXKcX/hxgF+DWCC6HS6w4dLiorKbbaYlJR8lSoWLTiRvI0YFzc5JWWR1RpTXNz6xRcHz5690N/fT3W0oGaz2Q4cOKDX63NzcwmCuHz5clpa2ldffdXQ0LBly5bXXnvtwoULg+//9ddfq1SqgRYkCIIkSYVCQRDEM8888/jjj6tUqpMnT+bk5Hh/srt27brjjjv27t2rUqmefvrpv/71r0MzVFRUTJs2TafTCQSClStXfv755wRBVFZWvvjiiwP3efLJJ9va2giCeOedd1atWvX9998rlUqz2Txjxoy+vr7Jkyd3dXWVlZWNy2s0/vBOBONOp9NduFCl1dpCQhKTkiJxNIZaDAZDLlfJ5Sq73drc3Hrx4umwMPakSbFqtRrTf0+whx56aMeOHQ6H4+23346OjiYIYuvWrb/85S/feecdgiCmT5/++OOPP/7444O/pLGxUaPRDH2onp6e119/vaKiQqPRPPjgg4sXL96+fftjjz1GEMTkyZO9DyiXy7dt27Z58+ZVq1Zt2LDh/vvv9/4yPvLII/fee+8f/vAHgiBiYmKefPLJ2267bZjYM2fOfPvttwmCqK6uttlsTz/9NN0v2sEeIYyjrq6uf/7z+HfflVmtsUlJ+QpFFFrQf3A4vMjIhJSUBUxmSklJ5//+76EzZy709vZSnSuIvPnmm2az+fjx408//bR3Is2KiooFCxZ4by0oKKitrdXpdIO/hM/nGwyGoQ9VXV0tl8sHOnLOnDnl5eXej7OysrwfREVFDRxEHayiomJgFtbZs2fX1taazeZhYnt3XgmCSExMnDt3rkqluvPOO//3f//XT5ZwuAEoQhgXBoPh8OGSoqKLNltcUtK88HA1KtA/kSQpkykSE6fFx89vbRV+8825774rbmtro++bGu3k5uYuXrzYW4QCgcBisXi3WywWkiSv2NnKzc29ePHiFe14xRd6v3ZgSNR1p+Lj8/mDn5TFYnE4HBaL5XQ6B+4zuBoHJrxlMBh/+9vfzp8/P2vWrGeffXbwoVR6QRGCj7lcrsrK6v37TxgMyuRkVCBtcDjcyMiE1NQFLFbasWOtX355sKqq2uFwUJ0rMFmt1oGds56enpKSkkmTJhEEMX/+/N27d7tcLoIgPvrooxkzZlxRhAsXLkxLS9u4caPJZPJuOXXq1LfffpuSksJms7/55huCIIxG4549e+bPn3+tZxeLxUwms6enx/tpfn7+rl273G6390nnzJnDYrFiYmLq6+v1ej1BEIcOHbrqoQKDwWC1WmNjY3/729/+5je/uXjx4lhfF4rglAD4UkdHx8mTFVarLDZ23lWHjIP/CwmRh4TIjca+Cxfqz5z5PiMjKjFRQ/eTQP6mt7d38uTJcXFxAoGgrKxs7dq1d955J0EQjz76aHFxcVZWVkRERG1trXd452AMBmPv3r333HNPVFRUenp6Z2cnSZI7duzg8XgffPDBvffem5mZWVlZWVBQcOutt17r2Vks1qZNm9LT06Oioj7++OMnn3xy1apVOTk5ISEhLS0t//jHPwiCiIuLW716dXZ2dlJSklgsjoiIGPo4lZWVK1asyMnJYbFY58+f946yoSMszOtHaL0wr8lkKi0tb2y0qFTpISFyquOMi7EvzEs7druts7PRaGyIjw+dNCkpJCSE6kTXMczCvFu3bh28HuG8eQuEQpVSGTMBqY4e/fq55568++67B2+02WxVVVUWiyUhIUEu//evjMfjqa+v7+/vnzRpEofDudZjtra2NjY2KpXK2NjYgYOfJpOpsrJSqVRGRkZ6t5jNZgaDwePxCIJwOp1ms1kikXhvstvtJpNJIpEwmUyPx1NXV2cymdLS0ths9kCSyspKh8ORkZHR398vEomYTKbRaGSz2QNHRw0GQ1VVFYPBSEtLG9dfjY8++ujQoUNYmBf8lMvlqqqqOX++SSRKSk2Nw4HQQMLhcKOikl0ujVbbuHfv6eho0eTJGu9gfbp75pkt586dm5jnSk+/Z+HChVds5HK5U6ZMGXpnkiSvOi70CpGRkQNtN0AoFObk5AzeMnjyBBaLNdCCBEFwOJyBoiVJMiEhYWiStLQ078feqxIJgrhifiKxWOy97pDWUIQwJnq9vri41GoNjYvLx7HQQMVksqKiEiIjNV1dbUVFleHhVVlZKXSvwyVLlixZsoTqFOAXMFgGblxtbf3+/SfZ7FSNJhstGPBIklQoIlNT55FkclFR1cGDxwZmPwl4PT09zz33f5OTs4XC0JCQiAULln/55ZdUhwKfQRHCjbBarT/8UPLjj+0azbzw8BFNkAgBQyaLSEmZ43Jpioouff/98YHBh4GqtLQ0OTnj1VeLa2qeMJsP9vXt+eGH2XfeuWnlyltveFiDTqd74IEH3njjjYEtJpPpwQcffOWVV3yU+kaUlZW1tLQMfCoSicZ4XemRI0cGRrf6MxQhjFp7e/v+/cW9vWHJyTM4HB7VcYACJEnK5arU1HyHI/7//b8Lhw+f8I6zDzx9fX2LF6/Q6x+0Wr8niDsJIocg5ng8z9hsZQcPXn7ssSdv7GGNRuO2bdtefvnlgddtz549u3fvHjpMdCK99NJLe/fuHfj0L3/5i1AoHMsDrl27tqmpacy5xh3OEcIouFyuiorK8+c7oqJyJBIZ1XGAYt46DAtTdnW17t9fGh0tyMxMGxhVERjefPMtiyXB7X5uyC3hVusnb7899dlnt1z10oKRuOmmmz777LMHH3yQIIgdO3bcdtttg2cW/eqrr0pKSsLCwu65557w8HCCIOrr67/66qumpia1Wr1+/XrvxjNnzrS1tUml0q+++io6Onrjxo3eMaKDOZ3OXbt2lZWVxcfHb9iwwTu8s6am5tNPP+3p6YmOjr777rtbWlouXbpktVoNBkNmZubSpUsNBoP3soLPPvssIyOjuLi4srJy5cqVCxcu/Pvf/15cXDx9+vTbb7+dIAiXy7V3795Tp065XK758+cvW7aMIIgvv/zSZDK99957CoVi5cqVkyZNMplMH374YW1t7ZQpU+66667rXuk/YbBHCCNlMBi++664qsqRkjIPLQgDSJJUKKKSk+f39UV8/fWpH38sHX6CLnr58stvLZZfXePGdB4v4+DBgzf84L/+9a+91wM0NDQ0NjbOmzdv4KZNmzb99a9/zcjIMJlM06dP9+44fvPNN06nc+bMmZ2dnbm5ud6jjsePH9+8efO7776bmZm5b98+b60O5na7Fy9efPDgwalTp5aVlRUUFLjdbp1ON2/ePIlEMm/ePIfDUV9fz+VyORyOQCAIDQ31jjV99NFHvU+xffv2tWvX9vT0aDSa1atXP/DAA4cPH05PT9+yZcvu3bsJgujr6/v6669TU1NTUlKeeOKJt956iyAIoVDIYDAkEkloaCiHwzEajbm5uTU1NTk5OXv37r3iYhJqYY8QRqSrq+v7789JJGkaTRTVWcAfMRgMtTpeqYzRautra49mZcUkJycGwCzeWm0bQcRd61aXK37wSbXR8i7dcOnSpc8//3zdunUDlx5dunRpz549DQ0N3l235ubmDz/88OGHH/YuOuhwOBYvXnz69OlDhw7ddNNNBEGIxeJPPvmEJMnc3NzBber197//3WKxfPTRRyRJ3nnnnXl5eYcOHQoLC+Nyuffff//gg58JCQmzZ8++//77h0a95ZZbnnnmGYIgTp8+3dra+vXXXxME0d/f//XXX99xxx0ymez9998nCMJgMCiVyueee27Tpk1Llizh8/lr1671XoPxxz/+MTs7+89//jNBEGvWrImNja2vr4+Pj7/hV8+HaP/fFCZAY2PT0aNVajUOh8J1MBjMyMhEuz3q4sXqS5cO5+YmR0dH0/rSUoFASBDXHO7BZBquuK5utNatW7dz5849e/Z89913J06c8G48f/682WyeO3eu91OtVuttxG+//da7oASfz6+vr29tbfXeIT093fsiq9Xqnp4ej8cz+DU/d+5cbW3twEzZ9fX1NTU1BQUF06ZNU6lUhYWFq1evvv3224dfnDI9Pd37gUKhGFhAWKFQeIdKWa3W+++///Dhw+Hh4Tab7arjp86dO3fkyJGBiw77+/trampQhEADHo+noqLy3Ln2uLhZfP6YTptD8OBwePHxU4zGvuLi8vDwxtzc9NDQiVgIfjzMnJlTX3/Y5brlajdaHI4fc3L+ayyPv27duuTk5Nzc3ISEhIEiFAqFGo2mqKho4G7emVzWrVu3d+/eGTNmEASxfPly75SkxPWm1RaJRAUFBd6VmLwEAgGTyfzb3/6m1Wr37t37/PPPt7W1XbHe0xUGapIkyaFP99577+l0uvr6eiaTWVJSsmrVqqvGuOuuu5588snBW4Z5xomEc4RwTS6Xq6TkzIULfUlJc9CCMFoikTQ1dbbHk7h/f+mJE2cGr41AIxs33keSHxBE5dCbGIw/xMfHT58+fSyPr1ard+7c+ac//WnwxlmzZrW1tZ09ezY0NNR7xs5sNjudzr6+Pu/AnJqamsOHD4/wKZYuXXrw4MGenh7vozEYDIfD0dvbazabVSrVAw88sHbt2suXLxMEERoa2tnZeQPfRU9PT3h4uHeqtvfee29g++AHXLFixRdffOFyubwxbDbb8PugEwl7hHB1Vqu1uPiUXi9NTs6h9aEtoJZcrpLJFK2ttfX1xdnZccnJif7z9jcSM2fOfOSRzX/9a4HN9hZB3EwQ3t8FHYv1Bzb7/c8++2Hs387Q2bHDwsI+++yzDRs2qNVqDodTX1+/ffv2wsLCJ554Ys6cOZMnT9br9d79wpHIzMz805/+NHfu3NTUVIfD0dzc/P3337e2tq5Zs2by5MkEQbS0tHgv21i/fv1tt9328ccf33bbbS+99NLIv4V169bNmzdv/vz5er0+Ly9vYPujjz56++23CwSC119//eabby4rK0tPT588ebLRaOzu7i4vL/eTmXsx6bYf8Z9Jt/v6+v75z9NMZpxafeX0g8EsCCfd9iGr1dzScpHP78/Lm3zD1xtc18gn3R45j8fz1ltvP/3083Y7g8lMIQiz1XohOztv5863vUVyA9xud19f3xVHjO12u81mG3gTcLvdDQ0NNptNo9EMTHLd2tqq1+tTU1OtViuLxeJyuTabzel0ese8eDwevV5/1QPRDoejtraWxWLFx8d7j21aLJb6+noWi5WQkDD4aKd30K93N9Q7fNRgMHjHlHq/iiAI72+B3W632+3eI5xWq7WmpkalUoWGhhoMhsHzs/f39/N4PO+X22y2mpoasVgcExMzqr+wx3XSbRShH/GTIuzr6ztw4KRUmhEWpqQ6i39BEY6dXt+t1ZYlJ4dkZk4aeHP3ofEoQi+bzXb8+PH6+no+n5+Tk5OcnDy2pDA6WH0CJk5/f39REVoQxktIiFwimdfaWltff2TmzJSYmIlYBcknuFzuggULFixYQHUQ8D06HayH8dbf33/gwI8SCVoQxhGDwYyOTlap8oqLG4uLf6TFIJr29vbHHn8sNiGOxWYLhIKZc2d+8sknfnI4DcYOe4TwE28LikST0YIwAYRCSUrKnPb2hq++Ks7NTUhI0PjtmKwTJ04s/cUyfrwo7GZVtkrjtrm6Kns2/cdDH37y0b6v9g6dz2wkenp6HnnkkaHb//CHP8TGxo4x8B133OFyudhsdkxMzJ133pmRkTHGB7wBvb29JSUly5cvn/invgEoQiCIQS2IpSRgwpAkqVLFh4ZGnDp1obGxffr0Kf5wjvwKOp1u2Yrl8iXqyF/8e+CYKF4aMS/67Gulv3vk4W3vvHsDD8vn870lYbVaN2zY8Je//MW7SP3ghXNv2J49e1588cXY2NjS0tJp06adOHEiOzt77A87KvX19Q899FB9ff0EP++NQRECYTQai4pOCoWT0IIw8Xg8QXLyjI6Opr17S7KyotPSUvzq+oo333qTo+INbkEvlpCtuX/yB0+9//xz/3dgppWREwgEd9xxB0EQBoNhw4YNN910U2xsbHFx8dmzZ3t6eg4dOrR58+ZLly6lpaV5B6ZeunSprKxs7dq1BEHY7fYdO3ZcuHAhNjb2wQcfvGp3Ll68OCsra+3atefOnfvb3/4WHR29e/fuyspKkUi0Zs0a7/QuHR0dX3755Zw5cz744IOcnJxf/vKXu3fvvnDhAovFWr58+aJFiwiCsNlsb7zxxtq1a99++22Hw/HQQw8pFIq//OUvnZ2dd95558AlHJcvX961a1dPT09BQYH3avqPP/64r69v69atBEE8/PDDPB7v4sWLu3fvNhgMS5Ys8c7KXV9fX1xcnJKSsmvXrkWLFhUWFr733nsVFRUcDmfx4sUrVqwY7at6w/zoPxxQwmQyHThwQiicpFBEUp0FgldERIxGM6+szFhUdNRgMFAd59/27t8rzQu/6k18lUgWLz906JCvnquoqGjdunWHDh3Kz8/n8/k7d+48e/as96bz589/8MEHBEG43e5FixadOHGisLCwu7t7zpw5DodjmMe0WCwsFqu8vLyvr2/p0qVxcXG/+MUvzpw5QxBEW1vbli1b/uM//iMrKys6Orq5ubm2tragoCArK+uBBx7wrjxss9m2bNny0EMPZWZmEgRRWFh4zz33KBQKjUZTWFjoXZm5tLS0oKAgNDR0/vz5W7duffXVVwmCiIyMZLFYGo1Go9EwGIwjR44sW7ZMrVbPmTPnqaee2r59O0EQ1dXVjz/++H/9139Nnz49IiLiiSee+P7772+66abZs2dXV1f76lUdCewRBjWHw3H48EkeLxn7gkA5DoeXlJTb0dG8d2/J7NkpcXFjPVXmE+3tHaHzrvk3IkfOa2tr8+HTxcfHD54Lbai9e/fabLYdO3YQBHHzzTefPHnyu+++G7rztGfPnmPHjp04caK0tHT79u3Jycnz58+32Wxms7mtre2zzz7LyckhCMJkMn388cdK5U/DAl5++WWHw2E0Gu12+yeffLJ69eqB7ZmZmXfccceHH364cuXK9evXe5McOXJkzZo1zz///NNPP+1d9SIrK2vGjBlPPPFEQUHBm2++uWbNGu8jPPvss1u3bvWu2aTRaG6//fbf/OY3BEFYrdbPP//ce0j8qaeeuuOOO7w7ixMMRRi8PB5PSckZu11JoyHsEPAiIqKl0rDjx8+2tLTn5maNx7WGoyIUCZ2Wa+5yuS0u357X9O54DePChQt1dXUDU1c3NDR4Z0e7Qmtrq8vlysvL27p1q1qtvnTp0rp160wmU0RERHt7+8ApQ4VCMdCCHR0dv/rVrxobG9VqtcFgYLPZA4/mXT6CwWCEhYWlpqZ6N8rlcp1O541UXV09MLNaV1dXf3//0NgvvPDCa6+9RhCEy+Wqq6tzu90EQWg0moEXcMuWLb/+9a9fe+215cuXb968WaPRjOw18wEUYfA6f768tZWRlJRKdRCAn+HxBMnJs1pba/btOzJvXqZCoaAwzKy8WUWXDofnXeWQidvu0lV3Dazq4BPe6Ve8WCyW0+n0fmw0Gr0fSCSSgoKCzz//fPjH+d3vfpeVlTXw6ZYtW9atW7d582aCIF588cXy8vKhT/fKK6+kp6d7l1fcvXu3d70kr8FnbQfPvu29gEQsFr/22mtLliwZJo9EIvnggw8Gz742NEBhYWFra+upU6d27tw5e/bs5ubmCVvGC+cIg1RdXX15uS4hYarfjlmHYEaSZFRUskw29cCBsrKyCu/eAyU2Pbixo7jF1HTlLg5BEE17alJSUwd2znwuISHh2LFjBEHYbLZPP/3Uu3HZsmVFRUXnz5/3ftrd3e1ds3d4/f393onQ9Hr9xx9/PPx9bDab9xzeCK1ater1118fuB7Uu4caFham0+msVuvAfV599dWBucOuuhd7+fJlBoORl5f3+9//vqura+BrJwCKMBh1dXUdP345Pj6XycQhAfBfUmlYUtK8S5dsBw4UUzWCZtq0aVue3HJp6+nO460e109X0Nv7bHU7K3qPtn/60a7x+1Ny8+bNBw8ezMjImDp16sCMbikpKdu2bVu5cuX06dOzs7Pz8vJGsl7E//k//+fxxx+fO3fu7NmzrzVb929/+9v3339/9uzZmZmZKSkpI8/5zDPPxMbGajSauXPnJiQkPP300wRBxMbGrlmzRqPRJCQk6PX6F198kc/ne+8THx/vHU16hY0bNyYmJnoXSnzxxRcncpEmzDXqRyZmrlGj0bh/f4lCMU0ioesScVTBXKNU6exs6e29OGNGkkZznXVcx2mu0Z07dz759JN9/f0StdTtcPc29eQX5G97e1tiYuKNPeBgV6yjO5jT6WxpaYmIiBj6v8575FCpVI6wiU0mk1arjYmJGXw08go2m62lpUWpVA5etn6EvF8bHh4+zKWQFoulra1tmMfX6XQ6nU6tVnsn+x4Mc42Cz9jt9u+/PymRTEILAo0oFFEiUcjx46Wdnbpp0zIn7NTRgPXr1991112lpaV1dXXeSbejoqJ89eDDNBmLxYqLi7vqTdHR0aN6FqFQeN3a5nK5CQk3uODMSL6Wz+cPfx+ZTCaTyW4swFigCIPL6dMXnE5lZCQuGQSaEQhEKSlzGhsrdLpj8+blTPzi5pA5slsAACAASURBVCwWa/r06WNchhf8E84RBpGmpqa6OlN0NIaJAi0xGIz4+AySTNi/v6S9vZ3qOBA4UITBwmw2Hz9eGRs71a/mrwIYLYUiSqmcfvDgxfLyS34yxAHoDu+JQcHj8ZSUnBWJkgUCv5vUGGC0RCJpYuKcCxf6jh07NfwEYwAjgSIMCpWV1Z2dLJUqjuogAL7BZnOSk/O6usTffutfc5MCHaEIA19vb++ZM01xcVnXvysAfZAkGROTxman7NtX0tLSQnUcoLGRjhq1WCw//PBDY2NjVFRUYWHh4GnoBqupqTl8+DCHw5kzZ87AMFmtVltZWalWq0d1kSb4hNPpPHbsXHh4OodD8ZyNAOMhPFwtEIj++c/TGRm6zMxrrkDLZrPXrl17Y4vogm8Nc93kMIxG48qVK8cjDzHyC+pTUlLUanVaWtqpU6fsdvvRo0eHXvr95ptvPv/888uXL2cymQwGwzsH64YNG7744gsOh3PXXXf993//97UeHxfUE+NzQX1p6YXaWkKjmeLbhw1OuKDebzmdjrq60uhoMiMjRSwWDx0R1tvb29vbS0k2uILJZLqBC/YJglAoFON02cxIi7C5udl7/abT6czIyHjsscfuu+++wXeoqqqaOnXq6dOnvfOUD+js7AwLC/vd737HYrFQhMPzeRH29fXt23cyOXkBplLzCRShP/N4PI2NFUxm6+LFc27sfRYmxsRMoTUqIz1HODCLAYvFkkqlQ+fA/fLLLwsLC8Vi8b59+yorKwe2KxQKJpPpk6wwWqdPl4eFpaEFIRiQJBkXl+50xnz77fGhywABDGPUb5H79u2rr68fWLBxQF1dXVNT04oVK7Kzs7/77ruNGzf+53/+56ge2e12v/TSSwOfTp06deHChaONR2sOh8OHY8FbWlra2lzJycqBlVxgjDwej9PpxOvpz+TyKKNRuG/f8QULpoSHX31leaCWb9/orst7qm74+4yuCM+cOXPfffd9/PHHcrn8ipu8M67W1taKRKLKysqMjIz77rtPpVKN6vEHH8Tv7++ncO0VSrjdbl99y06n8+TJKqVyarC9huPN4/HgJfVnHo9HJovgcLhFRWfnzEmOisJsgn7Hh290IzGSKURGUYQXLlxYsWLFu+++W1hYOPRWtVo9ZcoU75nM1NTUkJCQ6urqURUhg8F49dVXR37/wGO32321Hnd1dS2TqZbJqFzRNPB4PB6XyzXM5P1AOZfLxWazw8IiBIK5J06czMlxpqYmUx0KfsaHb3S+MtJzhFVVVcuXL3/99ddXrVo1eHtzc7N3hEthYeHly5e9R43a29t7e3tjYmJ8HhdGwmg0lpW1RkVhTlEIXny+UKOZdeZMR2npeczEBsMb6ajRuLg4Lpc7f/5876dLlixZvXq1x+NhMBgnTpzIy8vzeDyLFi3i8/kFBQW7du3KycnZtm0bQRBff/31vn37jh49ymAwZs2adfPNNy9fvnzo42PUKOG7wVTFxT/q9Qq1+jort8FoYdSo/7NYLFwud+BomNvtqq0tVatds2ZNm/jFm+Cq/HDU6EiL8Kuvvhp8ejMlJWXKlCkEQXzyySfLli0LCwsjCMJut3/22WeNjY05OTkDbXfx4sWKioqBL0xPT7/i+govFCHho/8fHR0dBw9WpqTMG7+Fs4MWitD/XVGEBEF4PJ6GhnKpVJ+fn4fD2v6AxkU43lCEhI/+f3z3XTGDkSyTRfgkEgyGIvR/Q4vQq6mpksfrKCiY4W9np4KQHxYh5hoNKJ2dnd3dHrQgwBViYlKdzqhDh0qsVivVWcDvoAgDSnn5Zbk8ieoUAP5IrU5wu+MOHDhuNpupzgL+BUUYOLq7u7Vae1iYkuogAH5KpYpjsRK/++44Vm6CwVCEgaOi4rJMlogxMgDDiIiI4fPTiop+RBfCABRhgNDr9S0t5vBwzKMBcB3h4ZECwaRvvz2h1+upzgJ+AUUYIMrLq6XSBOwOAoxEeLg6JGRKUdEprM0EBIowMPT19TU29kdERFMdBIA2ZLKIkJDMAwdO63Q6qrMAxVCEgaCqqk4i0YxkblkAGCCTKeTy7KKiMzhGGuTw1kl7TqeztrYzPDyK6iAA9BMSIg8LyyoqOoUlDIMZipD2WlpaWCwFm425owBuRGhouFSaceDAj0ajkeosQA0UIe1VVjbLZDg7CHDjwsKUQmHqwYM/WiwWqrMABVCE9Nbf39/T45BKw6gOAkBvCkU0k6n5/vsTNpuN6iww0VCE9FZf3yQSReOqCYCxU6vj7Xb14cM/Dl5pB4IBipDG3G53VVUbhskA+EpMTIrJJD9y5KR3jXEIEihCGtNqtR5PCJeLVYEAfCY2dpJOJz5+/LTb7aY6C0wQFCGNVVc3h4RgmAyAj8XHZ2i17JMnz/rJcq0w3lCEdGW1Wtva+rH0IIDPkSSp0WTX1TnKyy9SnQUmAoqQrjo6OjgcBWaTARgPDAYjMXHa+fPd9fUNVGeBcYe3UbpqbOwQixVUpwAIWEwmKz4+9/jxyx0dHVRngfGFIqQll8vV0qILDQ2nOghAIOPxBFFRuYcPX8AEbIENRUhLXV1dTGYIi8WmOghAgBOJpCEh6YcOnbRarVRngfGCIqQlrbZTIMAwGYCJIJermMy44uJTLpeL6iwwLlCEtFRX1xkaihOEABMkMjKxry+kpOQMLqgISChC+tHr9TYbi88XUh0EIIjExaU3NbkrKiqpDgK+hyKkn/b2Dj4fx0UBJhRJkomJ086d62xoaKQ6C/gYipB+Gho6Q0JQhAATzXtBxbFj1Tqdjuos4EsoQppxuVw9PSaRKITqIADBiMcTKJXZhw+XYhBpIEER0oxer2exJJhQBoAqISFyBiOmpAQzkQYOvJ/STG9vL5uN3UEAKkVFJbW3MzFwJmCgCGmmo0MvEKAIAajknZX73DmtVqulOgv4AIqQZtrb9RJJKNUpAIIdi8WOiZlWXFxuNBqpzgJjhSKkE7PZbLF4sBIvgD8QCiVCYfLRo2cw4wzdoQjpRK/Xc7nYHQTwF0plrMEQUlp6geogMCYoQjrp6dFzOChCAD8SF5dRXW1sbGyiOgjcOBQhnWi1vWIxRsoA+BEGgxEbO/X48aq+vj6qs8ANQhHShsfj6e7uF4mkVAcBgJ/h84UyWfqxY2dxspCmUIS0YbVaPR42k8miOggAXEkuVxmNIeXll6gOAjcCRUgbJpOJxcKKEwB+Ki4uvayss7Ozk+ogMGooQtowmUxMJooQwE8xmazIyKzi4vM2m43qLDA6KELaMBrNLJaA6hQAcE0SiYzBiCotLaM6CIwOipA29HoTFuMF8HNRUSm1tebm5haqg8AooAhpA0UI4P8YDEZsbHZJySWLxUJ1FhgpFCFt9PWZeTwcGgXwdwKBmMfTnDiBdZpoA0VID1ar1e1m4doJAFpQqzVaLXn5ch3VQWBEUIT0YDKZ2GwcFwWgB5IkY2MzT5+uMxgMVGeB60MR0oPFYmEycVwUgDZ4PIFUmvLjj+dxgNT/oQjpwW63MxgcqlMAwChERMR0djIxH7f/QxHSg93uIEk21SkAYHSio9N//LHKarVSHQSGgyKkB5vNwWKhCAFoRiAQczgx589fpDoIDAdFSA9WK4oQgJaiopKqq/UdHR1UB4FrQhHSg8ViZ7NxjhCAfhgMpkqVceJEORZp8lsoQnrAoVEA+goNDbdaQysra6gOAleHIqQHFCEArcXETD5/vqm/v5/qIHAVKEJ6wDlCAFrjcLhSacrJkxdwWaEfQhHSA/YIAeguIiKmo4NsamqmOghcCUVIA06nkySZJElSHQQAbhxJktHRGadOVTmdTqqzwM+gCAEAJohQKPF4IqqqLlMdBH4GRQgAMHGiopIvXGjCaoV+BUUIADBxOBwenx9XVlZJdRD4txGtb+dwOIqKig4fPmw0GjMzM9evX8/lcofezWAwbNu2raqqKjw8/Ne//nVycjJBEDab7c0336ysrJw0adLGjRuv+oUAAD7R3a0NDQ0TiaRUBxmOWq2pqTmcktInlfp1zuAxoj3CsrKyZ599NjQ0NDs7e8eOHStWrBh6H71en5eXV1JSkpubKxQKq6urvdvXrVu3f//+/Pz8f/zjHxs2bPBldgAAgiAIwuPx7P7fN1asiX3kxWV3/Tbn7t9Mq6o6S3Woa2IyWSEhyaWlmIDUX5AjuajF6XSyWD/tO3Z0dCiVysbGxpiYmMH3efLJJysrK//xj38M3lhXVzdp0qS2tjaZTNbT0xMZGVlVVRUbGzv0KRwOh1AotNvtY/heaM9gMIjF4qHbnU7nZ58dTEtbOvGRYDCPx2O1Wvl8PtVB4ErvvP/cobo9kbfFMThMgiAsWmPztrq3thbFxCRRHe3qPB5PZeWRxYtTIyIiqM4y0a71RkehEe0RDrQgQRC9vb1MJnPoHv2BAwduvfXW999//4UXXigpKfFuLCkpyczMlMlkBEGEhYWlp6cP3AQA4BN2u23fdzuj7tR4W5AgCL5KpLhVvf3j31MbbBgkSUZEpJ05cwnX1/uDEZ0jHOByuR566KHNmzcPLcKGhobf//73t956a0hIyM033/z666/fddddWq1WLpcP3EehUGi12mEe/JZbbhn4dP78+ffee++o4tGdxWJhMplDtzudTpvNjiXNKOfdI8QFnf6mqalaqBaTzJ/9XKQpssp/lPrzb41AILl8mVFZWRUbG3P9eweQa73RjRMOhzN4X+6qRlGEbrd7w4YNLBbrlVdeGXorm81etmzZyy+/TBBERETEyy+/fNddd3G5XIfDMXAfm83G4/Gu9fgMBuO2224b+DQ5OXmYOwckh8Nx1W/Z6XSy2WwOB6tPUM/j8eAH4W9EIonLfuXCDi6ri8vl+fkPKzZ2Snn5yYQEzXXfqQPJtd7oxgmDcf0DnyN99T0ez4MPPtjU1LR///6rjvyMiopKSvrpcHxycnJbW5t3Y3Pzv+cTamlpiYyMvNZTkCS5du3aEeYJSAwG46o/MwaDwWCQI/lxwrjyeDwkiR+E31Gr4zx6t73PxpH++61JV9Ixf/YqP/9hicUhXV2KurqG1NRkqrNMnGu90VFoRGk8Hs/mzZsrKir27t0rEAgGtre3t3/33Xfej2+99dYffvjB+/Hhw4enTJlCEMTChQtbW1tLS0sJgjh9+nRHR0dBQYGPv4MgwGQyPR4XziUAXMszj7/b8Ealvrzb7XA7DPaOb5o9Zz133f441bmuT61OLitrxKRr1BrRHuHRo0fffPPNlJSUBQsWeLe89957WVlZJ06ceOCBB7wrL2/atGnBggX5+fkymezkyZN79+4lCEIikbz44ovLli3Lz88/cuTIyy+/LBKJxu+bCVQkSbJYTLfbxWQG0fETgJHLmbpg++tH3tnx3KW9pwUiydI5v7rjdw/TYp56Hk9AEPL6+oakpESqswSvEV0+YbFYrhjkolareTye0WjUarUDR0QdDkdJSYnD4Zg+ffrg0bG1tbWXLl2aNGmSRqO51lPg8gli2FHFf//7IYViJo8nuOqtMDFw+YT/s1gsXC7X3468Dc9sNnZ0lKxatSBIzhT64eUTI3rd+Xz+VTtMJBINtCBBEGw2e968eUPvlpCQkJCQcMMRgSAIHo/jdDqufz8AoBuBQOR2y5qamjWaeKqzBCk6/d0UzHg8NooQIFAplUlnz152ua4c+woTA0VID1wuihAgYAmFEocjpLm5heogQQpFSA/YIwQIbEpl0rlzl91uN9VBghGKkB6wRwgQ2MTiEItF1NraSnWQYIQipAcul+1yBfWQWoCAp1Qmnz1bgyuGJx6KkB7YbLbbjT1CgEAmkYSaTHzvtFwwkVCE9MDhcDwe7BECBDi5PKG8vI7qFEEHRUgPAoHA5bJQnQIAxpdMpujudvb29lIdJLigCOlBIBA4HCaqUwDAuJNIYmtqGqhOEVxQhPTAZrM5HNLhwNFRgACnUERfvtxps9moDhJEUIS0IZUKrFYz1SkAYHyxWGw2W9XQ0ER1kCCCIqQNFCFAkFAq4ysqGnEdxYRBEdKGRIIiBAgKAoHYZhO0t7dTHSRYoAhpQyQSOJ0oQoCgEBoad+lSA9UpggWKkDYEAoHbjSIECApyuUqrNRkMBqqDBAUUIW0IhULsEQIECZIkBYKY2toGqoMEBRQhbfD5fLfbisnpAYKEUhlbWdnmdDqpDhL4UIS0QZKkVCqwWnFZPUBQ4HC4DEYYhsxMABQhnSgUUqOxj+oUADBBJBJ1bS0WZhp3KEI6kcslZjOKECBYhIUpW1r0VquV6iABDkVIJyEhIU4nihAgWDAYDA5H2dqKhZnGF4qQTqRSqcPRj/kmAIKHTBZZU4Ojo+MLRUgnLBZLLOZYLBgvAxAspNKwnh4bLigcVyhCmomICMF4GYDgQZIkn69uacHR0XGEIqSZ8HCpxYIiBAgiYWGR1dU4OjqOUIQ0I5VK7XY91SkAYOKIRFKzmYll68cPipBmpFKp04nxMgDBhc9XNzVhp3C8oAhphs1mY7wMQLAJD4+sqdHiL+BxgiKkH5UqtL9fR3UKAJg4PJ7Abufo9TgtMi5QhPSjUsnN5m6qUwDAhOLxFO3tnVSnCEwoQvqRy+U2Ww/VKQBgQoWGRjQ2ogjHBYqQfvh8vkTCMptxgW2Qam9vqqwsNZuNVAeBCSUWh/b0WDDv6HhgUR0AbkRMjLy5uVsgEFMdBCbUxYun//M/HzEYwj2ecI/n3Ny5s5588iUeT0B1LpgIJElyOPLOzs6YmBiqswQa7BHSkkIRZrHgNGFw6ehofvjhBxobP9Lp/t7bu02v//Hbb5OfeupBqnPBxBGJIlpacHTU91CEtCSXy+32HoylDiofffSOXv80QWj+tYF0ODaVl7d2dWHyrWAhkykaG7vdbjfVQQINipCWOBxOWJjAYMBY6iBSWVnp8Uy7YqPTOa2+/hIleWDisVhskhTrdLh6ysdQhHQVFSXv68PR0SAiEokJ4spJtphMnUgkoSQPUILHU2i1HVSnCDQoQrqKiJBbrSjCILJixQqh8P2fb+tms0+npGRTEwioIJNF1NWhCH0MRUhXMpnM5epzu11UB4EJsmjRLdnZerH4PoIoJQgtSX4ZGrr0uef+yGRi7HcQEQolBoPLbDZTHSSg4FeIrlgsVlRUqE7XKZerqM4CE4EkyTfe+PiHH77et+/trq6O9PSM9ev3hoerqc4FE43LlfX29goEuGzGZ1CENBYfr/rxRy2KMKjk56/Iz19BdQqgEpcb2t3dGxkZSXWQwIFDozSmVCrt9i4cHQUIKhKJrK0NA0d9CUVIYxwOR6WS6PUYMgMQREQiqU5ndjqdVAcJHChCetNoVHq9luoUADBxSJJkMiVYsN6HUIT0plKpbLYOzDQBEFS4XJlOhyL0GRQhvXG53IgIUV8fVmUCCCIiUShOE/oQipD2EhJwdBQguEgksvZ2PWYb9hUUIe0plUqrtR2/EgDBg8ViEwTPYMCipL6BIqQ9gUAgl/P7+3GcBCCIsNkyzL7tKyjCQJCYqNbpsBYPQBDh8aQ6XT/VKQIEijAQREdH2e1aXFkPEDyEQnF3Nw6N+gaKMBBwudzY2NDubgyZAQgWAoFYp0MR+gaKMEAkJETr9U1UpwCACcJisV0uptVqpTpIIEARBoiIiAgOx2w2G6kOAgAThM0WY+CoT6AIAwRJkqmpkd3dzVQHAYAJwmCIjEb87esDKMLAERcXYzK14IJCgCDB44l7e7FH6AMowsAhFArVaqFO10F1EACYCAKBSKfDHqEPoAgDSnJyjF6Po6MAQUEgEPf0YI/QB1CEAUWlUnk8vXY7BpIBBD42m+NwkDabjeogtMcayZ1MJtNf//rX77//vr+/f/Lkyc8++2x8fPwV99mzZ8/nn38+8Om7774bGhpKEERZWdkLL7xQV1eXnZ390ksvhYeH+zA9XIHJZKakqBobW6KiEqnOAgDjzjtwlMvlUh2E3ka0R9jR0XHx4sVHHnnkrbfeYjKZCxcuHPo3yMWLF/v7+9f8C4/H835hfn7+ggULdu/ezefzb7/9dt9/B/Bz8fExRmMjhswABAMGg2+xWKhOQXsj2iPUaDQfffSR9+M333xTIBBUVVVNmTLlirslJiauWbNm8JZvvvkmKSlp48aNBEG89tprUqm0rKwsIyPDF8nh6qRSaVSUsKurTaGIpDoLAIwvJpOHQ6NjN+pzhOXl5SwWKyYmZuhNhw4dWrp06X333Xf69GnvFrPZLBAIvB9zOBwOh3P+/PmxxIWRSEvT6HS1VKcAgHHHZvNMJowJGKsR7REOMBqN69evf+6550JCQq64adq0afHx8QqFoqSkZO7cuf/85z9nzJiRn5//xBNPnDlzJicnZ8eOHf39/R0d1xzc73Q6B596vP3225955plRxaM7k8lEkuTYH0cgEPD5Vq22WSoNG/ujwWBWqxWHnf2Z1Wp1u90++T2iBZfL3dmpo9dl9b56oxshHo/HYl2n6UZRhBaL5aabbpo+ffqWLVuG3rps2TLvB4WFhVqtdtu2bTNmzEhPT3/rrbdWr17tdDpnzZqVnZ09zGAZJpN56NChgU/DwsJEItHI4wUAj8fjq285Nzft2LE2lSraJ48GXh6PhyRJPp9PdRC4JpIkuVwugxEs4+GdzlCHQ0uvt0ofvtH5ykiL0Gq1rlq1Kioq6t13371umavV6oFDoOvXr1+/fj1BECaTKTIyMj09/VpfRZKkRqMZYR4YXmRkJIdTZTL1C4USqrMAwHjhcnl6PQ6NjtWI/m6y2+1r166VSCQffPDB4D+16urq/vCHP3g/PnnypNvt9m7csWPHggULvNtbWloIgnA6nU888cSUKVOmTp3q4+8ArobBYGRkxHV21lMdBADGEZvNNZkwWGasRlSEFRUVR48ePXTokEKhkMlkMpnMewyzsbHxjTfe8N7nmWeeEYvFMTExU6ZM+eUvf+kdKUoQxJIlSyIjI6VSaW1t7RdffDFO3wYMFRcX63S24+J6gADGYDAIgmW326kOQm+kD8/8G41GvV6vVquvOEDf3t7O5XK919dfi8PhEAqFQf7jNBgMYrHYhw9YVlZx+TIzJibVh48ZzDwej9VqxTlCf2axWILqHCFBEDU1R5Yvz5JIaHMSxOdvdGPny/8uIpEoKipq6H9BpVI5fAvCOElIiDeZGl0uJ9VBAGC8MBhcXEo4RkH0d1MQEggESUnhHR1YuR4gYKEIxw5FGOAmTUrq66vFTiFAoCJJpsvlojoFvaEIA5xYLE5ODmtvb6A6CACME4Z3xD7cMBRh4Js8OcVgqHM6HVQHAQDfwx7h2KEIA59QKExLi9Bq66gOAgC+R5JMpxNFOCYowqCQlpZsNjfY7TijDhBoGAyGy4VDo2OCIgwKfD4/PT1Sq8WSFACBhsHAHuFYoQiDRUpKks3WjIlmAAIMg8FwOFCEY4IiDBZcLnfKlJjW1hqqgwCALzEYTBwaHSMUYRBJTk50u9utVjPVQQDAZ5hMHBodKxRhEGGz2VOmxLa1VVMdBAB8hsFg4tDoGKEIg0tioobB6DYa+6gOAgC+4V0vmuoU9IYiDC4sFmv69JTW1nKqgwCAb7jdLjabSXUKekMRBp3o6GiFwtPZ2UJ1EADwAbfbzWTinXxM8PIFo9zcjJ6eS5iJGyAAYI9w7FCEwUgqlaalKXApBUAAcLtdLBaKcExQhEEqPT3Vbm82m41UBwGAMcGh0bHDyxekuFzutGmJra0VVAcBgDHBodGxQxEGL40mXiKx6nQdVAcBgLFwMxh4Jx8TvHzBiyTJ6dMnd3RUYFVPAPpyu11MJvYIxwRFGNTkcnlCggRLFQLQl8fjwh7hGOHlC3aZmZPM5jpMQApAUx6Pk81mU52C3lCEwU4gEOTmJjY1nac6CADcCLfbyuPxqE5BbyhCIDSa+IgIT3t7I9VBAGDUXC4U4VihCME7amZKX18VDpAC0A6KcOxQhEAQBCESiXJzNU1NF6gOAgCj4HDYeTwWBsuMEV4++EliYkJ4uLOzs5nqIAAwUna7VSjE7uBYoQjhJyRJzpiRpddX2mwWqrMAwIjYbChCH0ARwr+JRKLs7NjGxjKqgwDAiNjtVrEYRThWKEL4mZSUpLAwW2dnK9VBAOD6cGjUJ1CE8DMkSc6YkanTVdjtVqqzAMB1uFxWPh9FOFYoQriSRCKZOhUHSAFowO028/l8qlPQHooQriI1NTk83N7WVk91EAAYjsNhkEgkVKegPRQhXAVJkrNmTbVYLhuNfVRnAYCrs9ttHI6Hy+VSHYT2UIRwdXw+f+7c9ObmUpfLSXUWALgKk6lfLsfuoA+gCOGaVCrV5Mlh9fWYbgbAH5lM/XK5mOoUgQBFCMPJzEyXSo2dnS1UBwGAK9nthpAQ7BH6AIoQhsNgMGbPztbrL5rNRqqzAMDPuFz9GCnjEyhCuA6xWDxzZkpj4xm32011FgD4icfjcTpNIpGI6iCBAEUI1xcXF5ucLGluvkR1EAD4icVilEr5TCaT6iCBAEUII5Kdnc5md+h0HVQHAQCCIAijESNlfAZFCCPCZrPnzMnu6rpgsZiozgIAhNncq1CEUJ0iQKAIYaRCQ0PnzEltaDiFKwsBKGe398jlcqpTBAgUIYxCTEx0Rob88uUzHo+H6iwAwctutzEYVgwZ9RUUIYzOlCmTIyPdLS3VVAcBCF56fXd0dBhJklQHCRAoQhgdkiRnzsxhMFq7utqozgIQpEymHrUax0V9BkUIo8bhcObPn9bbW44puQEoYbd3h4WFUZ0icKAI4UZIJJIFC6a0tJxxOOxUZwEILjabhcNxicW4dsJnUIRwg5RKZXa2uq7uNGacAZhIen13dDSOi/oSihBuXFpaSlwcu6npItVBAIKIydStUuG4qC+hCOHGkSQ5fXo2n9/d0dFEdRaAoODxeGy2blxB6FsoQhgTFos1FnjAigAAHf1JREFUf/50m60as68BTID+fp1czhMIBFQHCSgoQhgrgUBQUJDb03O+v19HdRaAAKfTtSUmqqlOEWhQhOADUql00aKpWu0ZLFsIMH48Ho/N1q5Wq6gOEmhQhOAbcrk8P39SY+OPNpuF6iwAgUmv746IEOC4qM+hCMFnIiMj8/Li6upOOp0OqrMABKDe3jaNBsdFfQ9FCL6UlJSQmRleW4uLCwF8zO12OxwdOC46HlCE4GMZGZOSk/n19eewQgWAD/X2diqVYh6PR3WQAIQiBN+bOnWKSuVoaCinOghA4Ojr0yYk4LjouBhREXZ1dd16660ymYzH4+Xk5Bw5cmToff785z8nDNLb2+vdvmvXLo1GI5fL4+LiduzY4cvs4K8YDMbMmTlicW9r62WqswAEApfL6XB0qlQ4LjouRlSEFotl/vz5lZWVJpPp7rvvvvnmm81m8xX36e3tnT9/ftG/eFeM1Gq199xzz86dO7u7uz///PMHHnigoaHB598D+CEWi5WfP53JbNZq66jOAkB7nZ0tGk04h8OhOkhgGlERxsTE/Pa3v1UoFEwmc+PGjX19fY2NjUPvJpVKNf/CZDIJgujo6GCxWHPnziUIIi8vTywWa7Va334D4Ld4PN7ChTNJsrGtrZbqLAD01t/fkJIST3WKgDXqc4TffPNNREREQkLC0Js++ugjhUKRk5Pz0UcfebdkZmYuXLhw06ZN+/fvf/jhh3NycqZPnz7Mg/cOgqEWAWCgC7XaBoqjANCWTtchlzNDQ0OpDhKwWKO6d21t7caNG999992he+irVq2688475XL50aNH161bFxoaunLlSpIkCwsL/+d//qelpaW6uvo3v/mNd0/xqpxOp0ajGfj0V7/61SuvvDKqeHRnMplIkqQ6he/l5WUcPnyqvt4SERFLdZaxslqt+BPNn1mtVrfbHWC/R62tl2bNUhuNATJt0wS/0fF4PBbrOk1Hjvy3uqmpKT8/f8uWLQ888MDw93zsscd6enp27tx5+PDhtWvXVldXh4SEGI3G1NTU7du3L1u2bOiXOBwOoVBotwf1Kq8GgyFQF9u0WCxFRSUMhkaliqM6y43zeDxWq5XP51MdBK7JYrFwuVwGI3DGw5vNhq6uH1etKgiYb8oP3+hG+sq2tLQUFBRs2rTpui1IEITH4/EWfnV1dWJiYkhICEEQIpEoNTW1urp6LHGBpvh8/uLFM93uuvb2q5xdBoBraW+vnzw5NmBa0D+N6MXt7u4uKCiYOXNmQUHBmTNnzpw5YzAYCII4ffr07bff7r3PBx98cOnSpc7Ozi+++GL79u233HILQRAzZ848d+7c/v37rVZrUVFRSUnJnDlzxu+bAX/m7UKnsxZdCDBCTqfD4dDGxcVQHSTAjegcYUdHR2xsbHt7+5YtW7xbXn311aysLKfT2d/f791y5syZrVu39vX1JSQkbN++fcWKFQRBZGRkfPzxxy+88ML9998fGRn5/vvv5+TkjNN3Av6Pz+cvWpRXVFTS0UFGROB3G+A62tsbU1KUXC6X6iABbhTnCMcVzhESfnnofDyYTKYDB0q43GTadSHOEfq/QDpH6Ha7qqv/uXLldO9l2QHDD9/oAuG/C9CLUCgsLJzpcl3G9YUAw9BqGzSa0ABrQf+EIgQKCIXCwsLZHE5rQ0OFnxyTAPArTqfDYKjNyEihOkhQQBECNbhcbkHBTJmsr64O61QAXKm19fKkSSqRSER1kKCAIgTKsNnsefPy1GrH5cun3W4X1XEA/IXdbrXbm9PSkqgOEixQhEAlJpM5e3ZuUhKnuvqEwxHUQ6UABrS0VGdkxGDpwQmDIgSKkSSZk5OZmSm7fLnEbrdSHQeAYhaLyeNpT06+ynzOME5QhOAX0tPTZsyIqq09brGYqM4CQKXW1sqpUxPYbDbVQYIIihD8RWJiwty5iY2NJSZTP9VZAKjR39/L4fTGx8dRnCPIoAjBj8TGxixcmN7W9mNPTzvVWQAmmtvtbm29kJc3aZhVemA8oAjBvyiVyuXL88zmiuZmzM8OwaW19XJcHF+tVlMdJOigCMHvSCSSZcvmiESdtbWluKwCgoTZbLTbG3NyMqgOEoxQhOCPuFxuQcGs+HiypuaE3W6jOg7A+PJ4PM3NF6ZPT8ZMtpRAEYKfYjAYubnZubnKurqjRmMf1XEAxpFW26BUEnFxsVQHCVIoQvBrSUkJBQWT2ttPYvgMBCqbzWIw1EyfPoXqIMELRQj+TqVSLVuG4TMQsJqby6dN02BaUQqhCIEGJBLJ0qWzRaLO2tqzGD4DgaSjozkkxJKUhHlkqIQiBHrg8XgLFsxMSCCrq4+azUaq4wD4gNlsMBgqZ83KJkmS6ixBDUUItMFkMnNysubNS2hpOd7e3kR1HIAxcbmcjY1nZs1K87fl2oMQihBoJjo6auXKWQxGfV3dORwmBfqqrz+fkSGPjo6iOgigCIGGRCLR4sWz4+M91dVHzWYD1XEARq2trS4szJKRMYnqIEAQKEKgKRaLNW1a9rx5Ca2tJzo7W6mOAzAK/f29NlvtrFlTGQy8A/sF/BiAxqKjo37xixkEcRmHSYEu7HZbW1tpfn6WQCCgOgv8BEUI9CYWixcvnh0T466uPm61mqmOAzAcj8dTX1+akxMdHh5OdRb4NxQh0B6LxcrLmzprVnRj49GODowmBf9VX38hNpaVkpJEdRD4GRQhBIj4+LiVK2ey2Y3V1T/a7Vaq4wBcqanpUmioYcaMqbhq0N+gCCFwiMXiRYvmTJsmr68/ggsNwa9otQ1cbmd+fh4W3fVDKEIIKCRJJiUlrFw5k8VqqK4+iSWcwB90drY6nbULFkxns9lUZ4GrQBFCABKLxYWFc7OyQmprj2DZCqCWTtdpNF5ctCgPaw36LRQhBCaSJNPSkn/xi2lOZ2VNzRmn00F1IghGBoO+p+fcokXTsLiEP0MRQiALDQ1dsmTupEnc2tojvb1dVMeB4GI2G1pbTy1cmB0aGkp1FhgOihACHJPJnDIlvbAw02Qqq609i7OGMDEMBn1T04kFCybjkkH/hyKEoCCXy1esyJ8yRdjYeKStrd7j8VCdCAJZX19PW9vJhQsz1Go11Vng+lCEECyYTGZqavLKlbOk0s6qqqP9/b1UJ4LApNN1dHeXLl06TalUUp0FRgRFCMFFKBTOm5dXUJCo052urT3ncNipTgQBpbOzpb+/bOnSPJlMRnUWGCkUIQQjlUp1003z09LYtbWHcek9+EpbW73dXr106UyJREJ1FhgFFCEEKTabnZExecWKPC63qarqONY1hDFqa6tlMhuXLJklFAqpzgKjgyKEoCaVShcunD1zprqtraSx8SIuN4Qb4Ha7amvP8vnaRYtm8Xg8quPAqLGoDgBAMZIk4+Pj1GpVZWXNxYuHRSKNUhmPFVNhhCwWU0PD6dRUaXb2TMwjSlMoQgCCIAgul5uZmZ6YGF9WVllT80+pNCkiIhqrBMDwdLrOjo5zM2YkJCYmUJ0FbhyKEODfhELhjBk5qan6s2cvVVbWh4cny+UqqkOBn2prq3W5Gn7xi1xMHEN3KEKAK4WEhCxYMLO7u/vUqYs1NY0qVdr/b+/Og5so3weAb5JNNleTNk2TpmnTNj1S7IHQYmkFSynUsxwV/0Hw1kEUT5QZGcf7BEcdFUQUGVQc7QAFFXEsCC1UKy1TBIHeoWmSpm3aNPdmN7vfP9Zffp0WsAg5aJ7POMzbd9/gE95mn7zvvruvWCwNd1AggpAk0dV1QqXyl5bOwTAs3OGAKwWJEIALk8vlt9wy12AwnDhxfGAgPilJh2GwewBAnM7R3t6W669PzMubBpPnUwMkQgAuisViaTQatVrd2dnd2nqUzVbGxSXDZjpRi6Zpo7HT6+2pqMiDZ6dNJZAIAfgXHA5Hp8vKyEjv6dE3NTVarUlJSdlCYUy44wIh5XY7entbNRpeUdFc+DI0xUAiBGBSUBTNzMxQKBKs1uHW1iaCkCqVWTExseGOCwQdTdNmc7fb3TVnTk5qqibc4YCrDxIhAJcBRVGtNj0tLbWvr+/EiRajUZCYmB0bKw93XCBYXC67wXAyNRVbuPAmuFl+qoJECMBlY7PZGo0mOTnZZDK1tp4aGMDk8gyZTBnuuMDV9H8Dwe65c3UaDQwEpzJIhAD8R2w2Ozk5Wa1Wm0ymkyfPtbV1xMWly+UqeCrNFGC19g8MnM3MlMyYcRPcIDHlQSIE4IqwWCy1Wp2UlDQwMNDeru/oOCMQaJRKDdxrcY2y20fM5rMyGXnzzXmwuXyUgEQIwFXAYrGUSqVSqXS5XD09vWfPNhCEJCEhHeZLryE47jEa29jsoTlzslNS4AF7UYRF03S4Y0AQBCEIQiQS+XxRvUuqw+GIiYFF+ZGLpmm32z2ZTXYoijKZTKdOdY+M+MVijVKpQVFuCCIEHo8Hw7DLnZ0mScJk6vT5DHl5KTpdForCCCGIIvBEB/0NwNXHXD5MTk62Wq3t7fqurk4MUysUqXD3YaShKL/ZrLfbu/Lzk6ZNm8fj8cIdEQgDGBFGkAj8ogTGmvyIcByv16vX954504vjfLFYLZcn8Xiw/iIoJj8i9Pnw/v4et7s3M1Oel6eD3XRDJgJPdDAiBCDo+Hx+Tk62Tpc1MjKi1/d1dBxmsaQxMcnx8YkcDnwGQ83jcQ0M6L3ePp1OqdOVisXicEcEwgw+hACECIvFkslkMpns+uvzBgcHu7v7OjpOc7kJcXHJcXEKWJoRAnb7sMXSxeHY8vNTtdpymAgFDJgajSAROGMAxvrPU6MXg+O42Wzu6OgbGvLy+UkymRr2e7pCF5wapShqaMhktXbLZHR+vlatVsO9nmEUgSc6GBECEDYYhqWlpaWlpblcLoPB2N5+oq+PwjBFbGyiVBoPJ+srRNO03T48PGzEcXNKSmxhYY5CoQh3UCASwYgwgkTgFyUw1lUfEU7kdrv7+y29vRaj0cbhxIpESrlcxePBIy4nixkRer1uq9XkcvXFxXEyMlQaTQrsFxE5IvBEByNCACKIUCjUatO12nSCIIaGhgwGy/nz7QTBxzCFTKaUSGThDjCikSRhsRhwfADDvNnZKo2mUCqFqWbw7y5jRDgwMGA0GjUaTXx8/MSjBEE4nc7Aj2KxmMvl0jRts9nGNhMIBBd8gjuMCJGI/KIExgrBiPCC/1Or1Wo2W7q7LU4njWFyoVAmkcgEAlju/w+Xyz4yMuD1DtK0Xa2W6HQZCQkJsPgoYkXgiW6yI8LZs2d3dHRotdr29vYHHnjg/fffH9dg//79y5YtC7y9mpqaioqKkZGRzMzMQBubzfbxxx+vXr36qoQOQDRgsVhyuVwul+fn5zqdzuHhYYtlyGhsd7loHk/G58tiY+Oj8D59gvDZbEMOx6DHMyCVolqtQqnMiI+P93q9AoEAsiC4LJNNhOvXr7/99tvZbPb58+cLCgoWL148b968cW1KSkrq6+vH1shksuHhYabc3t6el5dXXV19xTEDEKXEYrFYLNZoNLNmIW63e3h4eGBg2GzWGww4lyvDsHiJRCYSSabqKhuKopzOUZtt0OsdRBBHcnL8tGmKhIQsoVAY7tDAtW2yibCqqooppKamJicnm83mCzYzGAzx8fEX/L384osvqqqqEhMT/1ugAICxhEKhUChMTk5GEATH8eHh4cFBq9HY19fn4nCEKCrh8SQikUQkkly7T7GhadrtdjgcNo9nlCRtJOmUy8XZ2XKlUieTyaZqvgehd9mLZRoaGkwmU3l5+cRDf/75Z1lZmdlsrqys/PLLL2Wy/7+wT5LkV199tXXr1kv8zTRN19XVBX7U6XQpKSmXGx4AUQjDMJVKpVKpCgoQiqIcDofdbrfZHFZrl9Fo9/lYKCpBUYlAECMSSQQCccSmEIqivF6X0znqctkIYpQkHbGxAoVCqlTGSqXJUqk0YiMH17TLu32iu7u7rKzsnXfeWb58+bhDIyMjGIYJhUKbzVZdXZ2RkTE27e3du3fVqlUGg+Fij3UnCILP55eVlQVqFixYsGbNmst5L9c8p9MJT3uKZDRNezyea24izuv1OhwOp9NptdqHhpyjoy4Wi8dmC1BUwGYLeDwBhgkwTIhh/FA+740gcK/X4/W6cdxNkm6/3+33u2kal0iECQkxcrlUKpXGxMRc7kYQbrebz+dDvoxkIT7R8fl8Lvdf9n65jERoMBjKysrWrl37r6tdamtrn3vuuY6OjkDN4sWLc3Nz33zzzYu9BFaNIhG5mAqMFZZVo8GA47jb7fZ4PB6Px+Fwj456HA6Pw+GmKA6HI2CzMRaLy2KhbDaXw+GiKJfDQQN/oij30mmGpmm/nyRJgiB8JEkw/1EUQVEETRM0Tfj9OEm6+XyORCKUSoVSqVAkEopEIqFQyOfzr3Cdi8vlEggEkAgjWQSe6Cb7bctisVRWVq5atWoyaz6NRmNcXNzY1x44cGDjxo3/MUYAwFWFYRiGYWM/pAyfz+fxeHAcJwiCSWU47sNxl89Her0EjhNOJ4njBE1Tl/77eTwuhnExjCsS/VPAMC6XK+RyuVwul5k6gj3/QOSY1O+i3+8vKyuLiYmJjY397LPPEAQpLi6ePn261WrNy8tramrSaDQbN26USqUajebMmTOvvvrq2Psrtm/fXlJSkpWVFaw3AQC4Gng8HjyHGkShSSVCmqaZq3ctLS1MTWpq6vTp0zEMu+uuu5iZIp1Ot2vXrj179iQlJX3//fcLFy4MvJzFYq1fvz4IwQMAAABXCp41GkEicOocjDVlrhFOYXCNMPJF4IkOfl0AAABENUiEEeTDDz/EcTzcUYCLOn/+/LfffhvuKMClfPfddz09PeGOAlwUQRATn9AZdpAII8imTZusVmu4owAXdfr06V27doU7CnAptbW1f/31V7ijABdls9k++uijcEcxHiRCAAAAUQ0SIQAAgKgGiRAAAEBUi5TbJ/x+P5/P12g04Q4knHp7e9VqNYfDCXcg4MI8Ho/dblcqleEOBFyUxWKJiYm55p4HGz0oiurr6wvlqX758uWvvfbapdtEylOOOBxOV1cXSZLhDiSccBzHsGt1x5woAX0U4Xw+H5fLhY15I1mIP0Qqlepf20TKiBAAAAAIC7hGCAAAIKpBIgQAABDVIBECAACIapAIAQAARLVIWTUabRobGw8fPjw0NJSYmLhixYqkpKSJbYxG47Zt20ZHR6urq0tLS8cewnF8x44dubm54+rBVfTLL7/88ccfo6OjaWlpK1eunLiNLYIg586d+/rrrymKuvvuu3NzcxEEIUny8OHDx44dc7vds2fPXrJkCaxgDJLu7u4ff/xRr9dLJJLbb7991qxZE9vgOP755593dXUVFhYuX7480BenTp3auXMniqIrV67Mzs4ObeDRgqbp5ubmuro6q9Wam5u7fPnyCy4W/f3333/88UemLzIzM5lKl8u1Y8eOzs5OlUp1zz33KBSKoIYKI8Lw2L17t8fj0Wq1Z8+ezcvL0+v14xoMDw/PmjXLYrGo1eo77rjjwIEDY4++/vrra9eu/e6770IXcfTZuXMni8XSarUHDx4sLCy02+3jGrS1tRUXF3M4HIFAUFpaevr0aQRB6urqnn32WZIkFQrF888//8gjj4Qj9qjQ0NDQ3d2dnp7u8/kqKipqamomtlm2bNmePXuysrLefffddevWMZWtra033nijWCymabq4uLirqyu0gUcLk8m0bNmy4eHhlJSUrVu3VlRUTLxBbufOndXV1YmJiWw2u7S0tLOzE0EQiqLKy8v3799fUFDQ1tY2c+bMkZGR4MZKg3ArLS3dtGnTuMqNGzdWVlYy5Y8//risrCxwqLW1taio6MEHH3ziiSdCFmQ08/v9KpXqwIED4+pXr169atUqpvzMM8888MADNE17vd5AgxMnTqAo6na7QxZq1HrppZeWLl06rvLkyZNisdjhcNA03dbWJhQKR0ZGaJpesWLF888/z7R56KGHnnzyyRBHGyUIgiAIgik7HA4+n9/S0jKuzXXXXbdjxw6m/OSTTz722GM0TXd3dyMIwnQcTdOpqak//fRTUEOFEWGYGY1GvV6fl5c3rr6+vn7hwoVMubKy8ujRo36/H0EQkiQfeeSRTz/9lMvlhjrWaHXmzBm3263T6cbV19fXV1ZWMuWFCxceOXIEQZCxkz9er5fH40FPBZvX621ubs7Pzx9XX19fX1paKhaLEQTJzs5WKBTNzc3IRToOXHUoiqLoP1ffSJIkSZLpi7FMJpNWq2XKGRkZv/32G4IgSqVSLpe3trYiCNLb22uz2YI9fQ2JMGw2bNiQkpKSnp6+Zs2auXPnjjtqNpsTEhKYskKh8Pv9AwMDCIK8/fbbZWVlhYWFoQ43Kj3xxBNqtbqoqOiTTz5JS0sbd3RcH5nN5rFHCYJ45pln1q5dGzgXgKvujz/+0Gq1cXFxJEm++OKL446O7SAEQRQKhclkomm6v7//Eh0HguGpp55atGjRxHyWn59/6NAhBEFomj548KDJZEIQRCgU7t27t7q6Oicnp6CgYPPmzYFrh0ECiTBYVq9ejU4wY8aMQIPHH3+8paVl7969H3zwwf79+8e9HEXRwHw6U+DxeOfOnfvmm29eeumlkL2Lqe2WW26Z2Ed33nlnoMFbb73V3Ny8bdu2NWvWTNzlblwf8Xi8wCG/33/vvfdKpdL169eH4I1MVXv27JnYQSiKWiwWpkFRUVFzc/PRo0cdDkfgEmAAl8tl5lEYBEHweDwWi3WJjgOXa9GiRRM7qKqqamybl19++cSJE1u3bp348vfee2/z5s233nprSUmJy+USCAQIgtjtdmb6et++fZs2bXr66afb2tqC+zaCOvEKJuPZZ5+97777xlUuW7bs1VdfZcqtra18Pp+iqFdeeUWlUhUWFhYWFiYkJCgUihUrVoQ83mi0ZMmSQHcEzJw58+uvv2bKu3fvvu6665iy3+9fuXJlZWWlx+MJaZRR7Oeff1ar1eMqt2zZMnfuXKZMUZRSqTxy5AhN05mZmfv27WPqt2/fXlxcHMpQo81bb701bdq0/v7+izVwOp2NjY09PT1btmy56aabaJquqanJzs4ONFiyZMkrr7wS1CBhRBgGNE17PB6mTFFUS0sL8yx2HMcPHTqE4ziCIFVVVXv27GG+t9bU1Nxxxx0sFuvhhx/+4YcftmzZwvzGzJ8//4UXXgjjG5nCfD5fYNDg8Xj+/vvv1NRUBEFsNtvhw4eZ+qqqqsBKxZqaGuZbME3Tq1ev1uv1u3fv5vP5YQg9arjd7kC5ubk5sKFBU1MTM9t52223HT9+3GAwIAjS0NBA0/Ts2bORi3QcCIYPPvhg27ZtdXV1Y/dssdvtzLVAhkgkKikpkclkH3300d13340gSHx8/ODgoMPhQBCEoqienh65XB7UOOGh22FAEIRSqZwzZ05sbGxTU1NMTMyvv/4aFxdnMBg0Go3BYEhOTvZ6vfPnz2exWKmpqXV1dQcPHhy3FuDRRx/l8XgffvhhuN7F1Nbe3l5eXl5cXCwQCOrr6wsKCmpra7lc7pEjRyorK5kvK4ODg6WlpVlZWSiKnjp1qrGxUaVS1dbWLl26dNq0aYGdgHbv3h3l+4sFyc0330ySpFqt1uv1bW1ttbW1JSUlCILk5OSsW7fu/vvvRxBk3bp1NTU18+bN279//xtvvPHggw8iCNLX11daWjpz5kyv16vX648dOxYfHx/mNzMVdXZ2ZmVlpaWlBf55N2zYUF5e3tjYeOONNzKpZ8uWLTU1NQqF4ujRo/Pnz9+2bRubzaYoavHixR0dHRUVFS0tLR6Pp6GhQSKRBC9USITh0d/f39zc7HQ609PTb7jhBuY+X5/P19raOmPGDGadIUEQhw4dGh0draiomPhB7e3tZbFYKSkpYYg+Ouj1+pMnT/p8vuzs7OnTpzOVDoejra2tqKiI+dHlctXV1dE0vWDBAmZF3MjICLP4OyA3NxeGhsHgcrmamposFotCoZg9e7ZIJGLqT506pVKpAmOI48ePd3Z2zpgxIycnJ/Bah8NRV1eHouiCBQuY61LgqvN6vX///ffYmoyMjNjYWKfTefbsWeYBCG63+9ixY0NDQzk5OWOXUCAIcvz48a6uLpVKNWfOnGDv0gqJEAAAQFSDa4QAAACiGiRCAAAAUQ0SIQAAgKgGiRAAAEBUg0QIAAAgqkEiBAAAENUgEQIAAIhqkAgBAABENUiEAAAAohokQgAAAFENEiEAAICo9j/C85EPyf3j0QAAAABJRU5ErkJggg==", - "image/svg+xml": [ - "\n", - "\n", - "\n", - " \n", - " \n", - " \n", - "\n", - "\n", - "\n", - " \n", - " \n", - " \n", - "\n", - "\n", - "\n", - " \n", - " \n", - " \n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n" - ], - "text/html": [ - "\n", - "\n", - "\n", - " \n", - " \n", - " \n", - "\n", - "\n", - "\n", - " \n", - " \n", - " \n", - "\n", - "\n", - "\n", - " \n", - " \n", - " \n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "m = mean(convert(MvNormalMeanCovariance, results.posteriors[:β][end]))\n", - "Σ = cov(convert(MvNormalMeanCovariance, results.posteriors[:β][end]))\n", - "\n", - "p = plot()\n", - "# i = 2\n", - "covellipse(m, Σ, aspect_ratio=1,n_std=3, label=\"3σ Contours\", color=:blue, fillalpha=0.2)\n", - "scatter!([m[1]], [m[2]], label=\"Mean estimate\", color=:blue)\n", - "scatter!([true_beta[1]], [true_beta[2]], label=\"True Parameters\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# References\n", - "[1] Polson, N. G., Scott, J. G., & Windle, J. (2013). Bayesian inference for logistic models using Polya-Gamma latent variables. *Journal of the American Statistical Association*, 108(1), 136-146.\n", - "\n", - "[2] Minka, T. (2001). Expectation Propagation for approximate Bayesian inference. *Uncertainty in Artificial Intelligence*, 2, 362-369." - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Julia 1.11.0", - "language": "julia", - "name": "julia-1.11" - }, - "language_info": { - "file_extension": ".jl", - "mimetype": "application/julia", - "name": "julia", - "version": "1.11.0" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} From a5daa65d4410702b89fcda079223e3141a43ed5a Mon Sep 17 00:00:00 2001 From: Bagaev Dmitry Date: Tue, 14 Jan 2025 13:26:14 +0100 Subject: [PATCH 6/8] format ext folder --- scripts/format.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/format.jl b/scripts/format.jl index 5e55a137c..5f5facd43 100644 --- a/scripts/format.jl +++ b/scripts/format.jl @@ -10,7 +10,7 @@ s = ArgParseSettings() end commandline_args = parse_args(s) -folders_to_format = ["scripts", "src", "test"] +folders_to_format = ["scripts", "src", "test", "ext"] overwrite = commandline_args["overwrite"] formatted = all(map(folder -> JuliaFormatter.format(folder, overwrite = overwrite, verbose = true), folders_to_format)) From 5637b711ac5e1bcdc6958df9029fa00a939aad58 Mon Sep 17 00:00:00 2001 From: Bagaev Dmitry Date: Tue, 14 Jan 2025 13:32:54 +0100 Subject: [PATCH 7/8] rename test --- test/models/regression/binomialreg_tests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/models/regression/binomialreg_tests.jl b/test/models/regression/binomialreg_tests.jl index 0e9bcca08..6421e4bf7 100644 --- a/test/models/regression/binomialreg_tests.jl +++ b/test/models/regression/binomialreg_tests.jl @@ -1,4 +1,4 @@ -@testitem "Linear regression" begin +@testitem "Linear regression with BinomialPolya node" begin using BenchmarkTools, Plots, Dates, LinearAlgebra, StableRNGs include(joinpath(@__DIR__, "..", "..", "utiltests.jl")) From 2891509f4d92763b6f60f771a57709e8a1eafffe Mon Sep 17 00:00:00 2001 From: Bagaev Dmitry Date: Tue, 14 Jan 2025 13:50:45 +0100 Subject: [PATCH 8/8] bump version and ReactiveMP dep --- Project.toml | 4 ++-- codemeta.json | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index fa668524a..add7b452f 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "RxInfer" uuid = "86711068-29c9-4ff7-b620-ae75d7495b3d" authors = ["Bagaev Dmitry and contributors"] -version = "3.8.4" +version = "3.9.0" [deps] BayesBase = "b4ee3484-f114-42fe-b91c-797d54a0c67e" @@ -42,7 +42,7 @@ MacroTools = "0.5.6" Optim = "1.0.0" ProgressMeter = "1.0.0" Random = "1.9" -ReactiveMP = "~4.4.4" +ReactiveMP = "~4.5.0" Reexport = "1.2.0" Rocket = "1.8.0" Static = "0.8.10, 1" diff --git a/codemeta.json b/codemeta.json index a3e3598e2..8c0dba9ac 100644 --- a/codemeta.json +++ b/codemeta.json @@ -9,12 +9,12 @@ "downloadUrl": "https://github.com/reactivebayes/RxInfer.jl/releases", "issueTracker": "https://github.com/reactivebayes/RxInfer.jl/issues", "name": "RxInfer.jl", - "version": "3.8.4", + "version": "3.9.0", "description": "Julia package for automated, scalable and efficient Bayesian inference on factor graphs with reactive message passing. ", "applicationCategory": "Statistics", "developmentStatus": "active", "readme": "https://reactivebayes.github.io/RxInfer.jl/stable/", - "softwareVersion": "3.8.4", + "softwareVersion": "3.9.0", "keywords": [ "Bayesian inference", "message passing",