Skip to content

Conversation

@lorcandelaney
Copy link
Contributor

@lorcandelaney lorcandelaney commented Aug 7, 2019

See #772

Includes:

  1. Hyperrectangles-based methods for Slice Sampling: Adaptive and Non-Adaptive. Both methods are in pints/pints/_mcmc/_slice_hyperrectangles.py.
  2. Tests for Hyperrectangles-based methods.
  3. Notebooks for Hyperrectangles-based methods.

@codecov
Copy link

codecov bot commented Aug 7, 2019

Codecov Report

Merging #895 into master will not change coverage.
The diff coverage is 100%.

Impacted file tree graph

@@          Coverage Diff          @@
##           master   #895   +/-   ##
=====================================
  Coverage     100%   100%           
=====================================
  Files          55     56    +1     
  Lines        5673   5772   +99     
=====================================
+ Hits         5673   5772   +99
Impacted Files Coverage Δ
pints/_mcmc/_slice_hyperrectangles.py 100% <100%> (ø)

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 06e64cc...e43c4f9. Read the comment docs.

@MichaelClerx
Copy link
Member

Hey @lorcandelaney ! For these and other PRs, please check how your docstrings render, and if you're using the correct syntax (:math:`1 + x` for maths, ``1 + x`` for code) everywhere!

See CONTRIBUTING.md for instructions on locally building and inspecting the docs

Copy link
Collaborator

@ben18785 ben18785 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Great stuff Lorcan -- looking really good! Apart from the inline comments, just one ask from me:

  • Can you go through the text in the notebook examples once more? I noticed a few spelling mistakes in there (these are hard to comment on in Github's interface).

- [Slice Sampling: Stepout MCMC](./sampling-slice-stepout-mcmc.ipynb)
- [Slice Sampling: Doubling MCMC](./sampling-slice-doubling-mcmc.ipynb)
- [Slice Sampling: Overrelaxation MCMC](./sampling-slice-overrelaxation-mcmc.ipynb)
- [Slice Sampling: Hyperrectangles MCMC](./sampling-slice-hyperrectangles-mcmc.ipynb)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can we sort this so it's alphabetical please (sorry)?

"metadata": {},
"source": [
"# Inference: Slice Sampling with Adaptive Hyperrectangles\n",
"This example shows you how to perform Bayesian inference on multiple toy problems, using\n",
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Change to, "This example shows how to perform Bayesian inference on multiple toy problems..."

"source": [
"# Inference: Slice Sampling with Adaptive Hyperrectangles\n",
"This example shows you how to perform Bayesian inference on multiple toy problems, using\n",
"Slice Sampling with Adaptive Hyperrectangles, as described in [1].\n",
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

no capitalisation on slice sampling and adaptive hyperrectangles

" horizontal “slice”: $S = {x: y < f (x)}$. Note that $x_0$ is\n",
" always within $S$.\n",
"3. Find a hyperrectangle ($H = (L_1, R_1) ×···× (L_n, R_n)$) around\n",
" $x_0$, which preferably contains at least a big part of the slice.\n",
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This bit doesn't seem to be parseing for me. Can you check the formulae?

"\n",
"The implementation uses estimates (``width``) of the relative scales of the variables to randomly position a hyperrectangle with such dimensions uniformly over positions containing $x_0$ that lead to $H$.\n",
"\n",
"When a sample is rejected, the algorithm shrinks adaptively the hyperrectangle. Specifically, only the axis corresponding to the variable $x_i$ is shrunk, where $i$ maximises: $(R_i - L_i) |G_i|$, with $G$ being the gradient of $f(x)$ evaluated at the last rejectedvsample. By multiplying the magnitude of the component $i$ of the gradient by the width of the hyperrectangle in this direction, we get an estimate of the amount by which log $f(x)$ changes along axis $i$. The axis for which this change is thought to be largest is likely to be the best one to shrink in order to eliminate points outside the slice.\n",
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

rejectedvsample -- think a spelling issue

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Where is R_i set?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is rejection sampling used for step 4? For example, if a point is drawn that is in the rectangle not in the slice, then it is rejected, right? Can we say this explicitly?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is it the gradient of f(x) or of log f(x)?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Again, is it the gradient of log f(x) of f(x)?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What happens if adaptation is turned off? Is None returned in grad?

}
],
"source": [
"import os\n",
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please install PINTS, using e.g. pip install -e .[dev,test] (see the README and CONTRIBUTING).
After that these lines import os and os.chdir won't be necessary anymore and can be removed

"mcmc = pints.MCMCController(log_pdf, 3, xs, method=pints.SliceHyperrectanglesMCMC)\n",
"\n",
"# Add stopping criterion\n",
"mcmc.set_max_iterations(2000)\n",
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do we really need 2000 iterations for a unimodal Gaussian? Keep in mind these notebooks should run quickly for users, and are tested regularly, so the faster the better!

"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAEYCAYAAABBWFftAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzsvXmYHNV56P073T2LRqN9hlVCAoEAE2EWYWHkYAUbLiYx3NhJsJ/gx4kDOCb4ftf3+sln53qJifj4nI8Y7NiykY2RwCy2BQgj4wACLSwSkmY0aKSRZt/3nul9r6rz/VFV3VXd1bNpZjTC9XseabqrTp1zaunu9z3vJqSUuLi4uLi4uLi4uLi4uJw6ntM9ARcXFxcXFxcXFxcXlw8KroLl4uLi4uLi4uLi4uIyTbgKlouLi4uLi4uLi4uLyzThKlguLi4uLi4uLi4uLi7ThKtgubi4uLi4uLi4uLi4TBOuguXi4uLi4uLi4uLi4jJNuAqWi4uLi4uLi4uLi4vLNOEqWC4uLi4uLi4uLi4uLtOEq2C5uLi4uLi4uLi4uLhME76Z6lgI8UvgL4AhKeWfGNv+FbgHGDaa/YuU8hVj3zeBfwBU4H9IKV8db4yqqiq5atWq6Z+8i4uLi8tppaamxi+lrD7d85gq7u+TM5rUUKVKiafkdE/FxWVuEBnQ/y445/TOw2VCTPS3acYULGAr8GPgybztj0gpH7ZuEEJ8CPgccAVwHrBLCLFGSqmONcCqVas4fPjw9M3YxcXFxWVOIIToPN1zOBXc3ydnftP4G/wJP/deeS8+z0yKIC4uZwi7H9L//tk3HXdn1AzRTJQl5UtmcVIuxZjob9OMuQhKKfcBoxNsfgfwnJQyJaVsB1qAj8zU3FxcXFxcXFxASsnhgcNE09FZGc+f8AOQ0TKzMp6Ly4ww2g6aNr19Sum4eWfbTp49+ez0jjVHkFKias7nfaZzOmKw7hdCHBVC/FIIYarj5wPdljY9xrYChBD3CiEOCyEODw8POzVxcXFxcXE5Y8moGRpGGhhJjMz4WKPJUQ4OHOTVjnG98qcVWUSYdHGZ8wQ64P3noPOd6e23yGeiP9Y/vePMIXadGOJHbzSf7mnMCLNtn/8p8G+ANP7+B/AlQDi0dXzSpJRbgC0A69atc7+hXf5oeOT1pqL7vnbzmlmciYuLy0yiobGnew/Xn3s9y+Ytm9GxpPFT61qUXFwmSMqw9iYm6qQ1QaTGH1vuuWO9IUBfcBHCSRU4c5lVBUtKOWi+FkL8HNhpvO0BVliaLgf6ZnFqLi4uLkXJZDL09PSQTCZP91Q+cJSXl7N8+XJKStykByZl3jLm+eYRTodP91RcXFxmCzm2y+EHUQkxUTWJzzu1c0ukVfY2DXHTZWdT6ps7CuqsKlhCiHOllKat8y+BY8br3wHPCCF+gJ7k4hLg4GzOzcVlLrC5bnPRfbVh3V3omoV3ztZ0XAx6enpYsGABq1at+sD+wJ0OpJSMjIzQ09PDhRdeeLqnM6fwCA/aOALXdCAMBxLp7DRySiiagkS6GQNdpoyUktZgKxctvgiPmCXhWUr9n2eWhfVZ+LzPNTxCoEmJokl83qn1caB9hBP9Ec5aWM41F8ydRCAz9vQIIZ4F9gOXCiF6hBD/APy7EKJeCHEU+DPgawBSyuPAb4AG4L+Afxovg6CLi4vLbJFMJlm2bJmrXE0zQgiWLVvmWgYd8AjPGR+ntO34Nn5+9OeO+2ZCoXP54NEUaOK1zteo99eP3zgd1zPyjbSe2qAnd8Le759aH1NhPAvWB/Az4zW0kOlIdDHXfp1nzIIlpfy8w+bHx2j/IPDgTM3HxcXF5VRwlauZwb2uzniEB40ze0U7pabGbaNpkrSqUV4yxeVrl5mlfjtUng0X/ulpGT6uxAEmluUyakShdB+EZaunPujAsfHbTBJVUznQf4B156yjzFvm3GgCLoJzTos4RTweAapEPcMXk5xwi1C4uHzAGMvN0OS+q+6bhZm4uLhMFYGYERfBUCrEorJFtAZbKfeVM883D9CFt4a+MEd7gnzuIxdM+7j5mNa5d1r9HO4IcP9NF1PidXaqGYgNUD2vGq9nZpQwRVPcmlzF8Dfr/06TgjU1F9a5J6w3BZp4f/h9VKly4/IbnRvNcRdBTZNseauNGy+p5kPnLZyWPr3GApuqzr17dqq43yguLi4uk2SsjI5TYSJZIB988EGeeeYZvF4vHo+Hxx57jPXr10/rPKxs3LiRhx9+mHXr1k25j5MnT/L3f//31NbW8uCDD/L1r399Gmf4wWYmXAR7o7281PISn7jgE7zR9QYAn7vsc9n9rx4fmNbxJkLrkG6ZCCUyVFUWruz7E35eaH6Bq6qv4obzb5j28Qdjgzzf/Dx/cdFfcMHCmVcsp4v+UIISr8fxmk0XUspZNZh0hDp4pf0V/u6Kv6OipGLyHcxha7ipICqaMkajueMiGE5mWFDms3kYpFWNRFplT9PQpBWs+uF6IunIjHyG5yquguXi8sdA+1v294GQ/X2RCvIuc4P9+/ezc+dOamtrKSsrw+/3k06nT/e0xmXp0qX86Ec/YseOHad7KmccHkBLTS2LYDSloGqSRfPsySVGjbTSg/FsQt8ZTXIxFjI6CIvKdBchQCuiTMYzuovYSHJmaoINxYcA6Ah3cMHCC+gJxIkkFS4/dyHhZAYpKbiOc4HnDuqlQ0+5RIemQmQAFhWWHv3p+z/lolgnt85feWpjTJBjI7pr3nBimJUl+pjmc+mZTMqA6VqYkHLalLYJfc5Og4IlpeRoT4hLz1mQddMdiiR5+kAXGy+t5mpL0gjFiJPyeRyuiZKCwWNw3jWO1+ytXl0GKaZgndKZzVHjl6tgubicYdSGf12wbXPd5GrlbA4etW/Icyt0XQjnFv39/VRVVVFWpq9WV1VVZfc98MADvPzyyyQSCW644QYee+wxhBBs3LiRq6++mpqaGoaHh3nyySd56KGHqK+v584772TTpk10dHRw6623sn79eo4cOcKaNWt48sknqaiwrx6/9tprfPe73yWVSrF69WqeeOIJKisr+cY3vsHvfvc7fD4ft9xyCw8//LDtuLPOOouzzjqL3//+9zN/kT5giJEWZKAbVt4M8yaXGevn+9qAqQvfM50OWqoZAu88QeWF1yK42thY2C6WUhiOJoAJxuql45AMwcJz8waUEOqGRSsKhL8Sr648ZbQMmtR45O2XObv0Mi4/98M8/lY7cHrqDHaPxjl/8bysAmqj9U0+NHiEhrM/feoDte+Frvdg3Zdgwdm6oNy6G+2ijQC0ZeyLcf6En9rBWj658pMzltXParmd3AKAea1OTeJWpIYHcUpZ4DpCHZR6Szmv8jx9ZuZzN9bUxlMMJ3laUkpah2OsWlaBr5j7bTjJmyeH6A7E+Ysr9bmGE3pNvJ5AwqZgRZO69c3j9Fls+i8YbID51bB4eizBb/e+TTwT55ZVtzjubxmKEE3l8uHNtXjeuZMw3sXFZcrsbx2h+8hrjv8KrFcuZxy33HIL3d3drFmzhvvuu4+9e/dm991///0cOnSIY8eOkUgk2LlzZ3ZfaWkp+/bt4x//8R+54447+MlPfsKxY8fYunUrIyO6RaCxsZF7772Xo0ePsnDhQjZvtivbfr+fTZs2sWvXLmpra1m3bh0/+MEPGB0d5cUXX+T48eMcPXqUb33rW7NzMf5I8MRH0ZCQSRBOZmjoc7ZmNQ9GONwxsYKnpgDiJMBamYaEXmMyEknSPBShr6s5q/DkDymlZMu+Nl5+v89wVZuA8HTkKajZWrh9+CRq7VPQX1ewyyf0deaMlqEr3MVQuonuVI2tzdHho/zy2C8ncmoFPH3iaV7teHVSx3SPxtle08PBYve16z0Wpvqd9wEcfxE63p7YYFHdgoeZRKL7Peg7gtLzXlbgD6hJyOiK7hudb9ASbGE0mZtbR6iDp088nbU2ZqnfDi27io+taXDiZYgMwu6HEMFe4BQsNRMRsDNJPdNg35GiTbaEjvFSrO2UYqJeaX+FHS05y/10WLAmS+twjJff76O2K1i0TTKjj5lRrWM7W5VrOgMARJIObo5p496P5QLpwFh3+ujwUVqCLUX3v/x+P7tPDk1qvNnEVbBcXFxc5jiVlZXU1NSwZcsWqqurufPOO9m6dSsAu3fvZv369axdu5Y333yT48ePZ4+7/fbbAVi7di1XXHEF5557LmVlZVx00UV0d+suRitWrGDDhg0A3HXXXbz9tl0wO3DgAA0NDWzYsIGrrrqKbdu20dnZycKFCykvL+fuu+/mhRdeKLB6fRARQniFEEeEEDvHb32KY0lpZA3zsv1wD68eH0BRCwWwnUf7eavZP7E+J2gNmOn08GkjoD2dzljG1P8eHznOQGwgq+SZc+kMd7K5bjNpdQzX2Lgh9Gv26zQQaOOx0DG6R5uLHtoV7mJXl64MaJYqMUKqvH3sVyQnksXOgVAqRGtwcmnDQ4YFwfw7aYZOnvLCmtRUTPH32UgTpPTzd1LSa4dqCaVChFJ5ruf+Zug+VHyQSJ+esa/uV3rfI/r9cXo+zW31PSH2Ng0DkFEzdIe7HSY/xvNrut32HLZt3ly3mX09+7Lv+5XY9LkaWqc2WQVL07LXfrKK53BEz+KZcfjeMFGNz4rXUvPLNJrmn/7KZfp3/EXV8yc1j4lwppekcMJVsFxcXFzOALxeLxs3buR73/seP/7xj3n++edJJpPcd999bN++nfr6eu655x5bTSnTpdDj8WRfm+8VRV9pzHeryH8vpeTmm2+mrq6Ouro6GhoaePzxx/H5fBw8eJDPfvaz7Nixg1tvvXWmTn0u8X8BJ2ZjIAF6knbhIZ7W79UpiyBZwWnsngosWDE/DOTVIcokslYNJ6xj5Ne4EcI4F6lZnLr0Nnu79/JC8wvZ1fN8oTKWiY05d31Ae3r4vqSugHanCi1C1uQDpvJmHbMq1gLBbgh0zpoQaJ67Y6xLzKJMT8t88sYQHmJphUQ6bXvgxjr3jKorgpO3PBljK2ljaCNqxexm4BgiEbKNv+vEILWGJWVfzz5ebnuZQDJQeC5qpkDR1jvSiKcVounCfcf8xwraThdOimkkHbFnCnUar32vbm1Lx4pf33CfHkeXhzmW1+k5MnCKqzJdAPMtWOZPg8+xALPZdm656Z1OXAXLxcVlfHY/NPY/lxmlsbGR5ubc6ntdXR0rV67MKlNVVVVEo1G2b98+6b67urrYv38/AM8++ywf+9jHbPuvv/563nnnHVpadFeNeDxOU1MT0WiUUCjEbbfdxqOPPkpdXaH71QcJIcRy4M+BX8zGeB6kLlAJT1Y4K5YIQpUZth7bSluobcrjWQW/gnGO/hpO7NTjc0zeflT/h56+uT9kV7bMpBShRIYfvdHMYNhSTNrs3iJQlnS+pdcvKpiDXayckBCfZ+WSxjjCIWbI2VqSm9fq0X3Zuc5E2nwnTKE3mnJwtzJc+lIoSOlszatLDbM5eBRVUx33j4nwcLQnxO6Tg1g1LJlXk8163TKarmBN+voIwfH0KK/Gu4yhffa+T7yMPPGS46FSSkJpXflKKvmFyiXsexjqf+NwoMb7PSEOtAcK9xU2HreFJrWiSn9a0bKLI2aSDvPcoukoTzU8xcGBgzSmA7wULeKSGO7T/xpKLIFO/Tc3nHMRzRx8gtSBwq+liag8ilqohHmyymBef6ZV2XJdGkcb8ScKLeia1D6QVqnJ4Ca5cHFxcZkksx30Ho1G+epXv0owGMTn83HxxRezZcsWFi9ezD333MPatWtZtWoV11133aT7vvzyy9m2bRtf/vKXueSSS/jKV75i219dXc3WrVv5/Oc/TyqlC9ibNm1iwYIF3HHHHSSTSaSUPPLIIwV9DwwMsG7dOsLhMB6Ph0cffZSGhgYWLpyeGiqzzKPAPwMLijUQQtwL3AtwwQWnFugtpCnGWITcIvKKIpPElTh7uvdw0aKLivc51YyBpgCbSYCvjKFIkmpLIox3W0c41DHK366/gLMWlgO5wrCjsTTnAT2BOGcb+4yzQUgNIcCrpSjtPQAjJeCNQWmlozCXnY6msLd7L9efdz3zSxzclfIF1bEULKeLajteklQ0OgfCHOkaZd2q6sL2Bj2BOIFYhrXLFxVtMx6KppBSdGG6bdhBcDeu+aueDmKJd4G1eXOX1KR0F7r07v+HeSs+om9f898mNgHjGiXTGezPnnFNkiGI6MJ9SlHRtFzq8fGeq3gmzonRE1xz1jUIIYinVXZFuynxCKi4AESROmfJEFK1u0u+3jCIr1wXYRVpKKIiz7dttL2wr6xldHxLy2hilKULzh2zzdu9b3PMf4y7195NqbfUtq+mM8BjgTa+dvOaAgtWQtEXJLoj3dTGDTdHRwXVOBch9GP9RomQUE82mUttVwBNwkfzjjQXKRwTpVjalCpRlox2AOeYQ9mOH4s3ut4AJc19kUR2ngA/e/9nXLz44qIJKvI5FVVstjOgThRXwXJxmSYmUuDXxWUqXHvttbz77ruO+zZt2sSmTZsKtu/Zsyf7euPGjWzcuLFgX0dHBx6Ph5/97GdjHn/TTTdx6FBhLMXBgwcLtlk555xz6OnpGbPN6UYI8Qcp5afGafMXwJCUskYIsbFYOynlFmALwLp1607pV19gCA5SUq4EKU0l0KRFefI3w1ADcDnmGvV4K8ZOCpb1tRC6/JkvWGmaRlRLsxCJP5ri6QNdfDac4IKlekyGP5oiqgzxTs8wt6y+jvllpVmrhgDi6iigKybRpGIZUxcor+vZhjzXULpHWo1zyc0QmT1FPTNasJXGQCMaGjevvHnMcybQiUyb8UPjWLCkpERNIL3280+kVSSwp2mI8xZXcp43DPOrwFL4OKNq/Paw/qyvXarq8WBVF485tdrBWgKpAJ+44BPZbU+feJrGoWGqud35oETO8hJSerOvg8kgg/FBLl2sj6kBdT1BLk8dYPG80kkrWMsSbQTk2dnN5rMlug6AEkdKyeNvt5PKaHRpo1xYXYqUMpug48sfv4gKICVVfJqK1+NlX88+2kJtnDP/HM6vPJ9nD/UwEE6yYrFe7Nq8PwXPcf9RUCSsvCm76XhfmDWr9etfWFtq/Dgn6ZQBUVPgZC7jaWI8d9T4CO0jR2DBWaTVNIMhheoFZdl051ammuSioTeIP5amyjzWtEpanr20kfUwn/5Q0hhbp2Uowu/f7+crG1dTWuJF0yRdo3Eu9b9GVSQCyfVQvgiBxKcmkHKe8zTzT8HfBGUrCtq1BFu4hYkpWCavd75OS6CFr1xlLPRpqm4591nqvWUSjtdqrjknugqWi4vL+HWyXFzOYIQQ1xTbBVw1gS42ALcLIW4DyoGFQohfSSnvmq45FpAVICRre7ejqgqadn1uf/12wskM+C7HJlBqGh5NQfMU/rw7pTFuGowQTyssKNUz9UlkQQxWTXKQQ/Fe/jYVIqHpgk4kabcoNMV303QCdrXU8x+3/R1NAX2lfSiSol97nY+mqwjE5nOyq5/FyR4kS2CM7IDWGCxVSiIJJVuPyimexYZ1e90z9Ax0MKqmHcey9RHu4+xYA4omcgkzjFmYc9mx/zhf9PyB14LL2XDb5zlrgW6V+9WBzlzzg1tIaQrajf/beX4GB/oPANgUrPFizNTWPdmzEJa5/7rx16hS5dJFqwFQVI2UotE1Emfx8lKHnopgKB4+NZ61VIFhwUrZE32kjAx08eAAYU1DW63xlpGA4pn3uviSlDwWOMaazte59cJbsy6EKSNGTsl/0LKZzK0LAAaJwt8kr/AipSSRMV0lTdNLoWtkINJLIhXmPN98o1/7szAUTjLc04x6jvUcx1kjiQxCuA3mLyWSTLO9xs/KZRV85prlhW0nEv+YpzT4E34OjvaQUFXLuRnKpEXBesXTRqUsYUNed70Bs8SB/v6dlhGu7HuOzNuVlP7Z/+JgxyjNg1Gu1pLGio4+ftnwUdb1Ps/A/L8FcopTzqqch03BHVvNsZaAkFJmE/eYfTcHmrP7pJQEmt4lOhSk8tZv5zoxXJPhs2OOdbpxY7BcXFxc/khZtWoVx44dG7/hmc8h4GHgP/L+PQwsHu9gKeU3pZTLpZSrgM8Bb86ocqUPaiSCkHjQBUbP8e025eG4U+r2hhf5SM84KcUtEtIbJwZ5v1sXXkURIbBXiQAQycTwegQ9ySO8quYUCvO4RcleSoPv8stjv6QzrO83E1yk1DiJjMqyRBtlaqTgWID6VC6WI2vjkipNAxEa+sJZYcxUlFqHIzzyehPt/hiPvN5ELO2cIjqQSBPLqIhQD/TaU7Br1tgiw9K1NN4M7z2Wm4vMvfCpSUZjaRakB2kZygnjwbhF4ZSSx0PHeeLYE9lNb/W85WBpmTwH20c51mfcL8uNVM3Mh5bEISCJpScZh2W4uUlkNvkEoFtL3/3PgqQkAEvjTZSFO9HQyCj69YwkFVqGovSHk7SMnATIutBlVOfsiFnB2+aa6JzsBGA4mqY/lOT/27cjq/APRZLEEvkxWfDsG//MjrcesCgx+lhmlseeYAIQ2fkDNAaaearhqQnEEgn2GZk8O0dyqerLlAjXd22ByAAikwQlVXAetr7zFKzfNP6Getllb5tVsOwLKFGRgdY3ebdjF/v79ttnJwSb6zbTHDnIa1ot7XE9rmsk6hzDVxLqpEH4icftLpZTccVLG9czmlR0pcnSR113kIzq3GdG04t8xxIJjvXmlOvu0XhB0py5iqtgubi4uLh80DkBfFlK+Wf5/4CJ5TifbbLLxTmB2TPSCkqKtJqmKZ2rbWMTN4absi9jeUkSxnJTahoMZwWXQsNCbsXZ5xFEYwcIyUIhdkF6gHLVuV5XznpWpO4VkvctCpZpwdJQsoqTNNqZ82kajOLRMhze/RL+1El6o3l1mPIYiHTTd+glRqKWZB0OspqaLwRb5jhV6v31Y9b0yY4lJT2B4tkZFTTeSOsxOwLnmB0ppU3hS2TylKxkyJZQxHa0UXjZ7MtEC3YWtK2P7qAt8bYxb33uVvfSbYkGfV9MT3jiNWKssjFTec+CCPdD3xFk7VMFM5DpmJ7NEVgab8OjKRzrCeGPplBkmrihSLYOxzjUNuxwYsa1UvWxNaP+2W+bfsurHa/iEQKJsM3/ZKCRSDrieN/b/FG6R/Xnra47SIe/0PI4L2O4c4Z6ofYp6D5o62swnKQvmLCsbIzh9pZtY7oIOjigdb1H3cnnOTJkr/Fl9jGQ1K1Dh9N6ohTVVg9PR5MavxzeR4sIcijm7JYupSSSjjAcL7zOndFh/v2NvYxEU4TiGX6yu4UDbSN0tDUy2n2C2sHa3HS7Cp8pE1Wqeh1Acs9AMJbm14fbafMXlkw4J1KPL1aYSfF04roIurh8wOkOFv+xNn3fp8Lm4NHcmyLxZ/dddd+U+3dxmUb+leILil+dTEdSyj3AnlObzviIrChmBiBJhtQEqVArrdFu2uJdrKbKnFTB8fNTQ4SOtzL/mk/pLlOWmBPTfc9KXzDJ0kqzO3t/5pHdkR7WpMtZmmhHjhE473w+znVgs5tkTpBKZFQShsAcyvTYTu83jb/h4ys+bhyicV7kKIvDBxhVJG95QqyhMIGKeXhXJkLV6EKO1vTw2Y8swiM8jsJzTGQKjhfA/PQgy5IDUCQXw3hMpFiyP5pGShjJdLCsZJVt33AkxSAxhoQu2AtDIE+MdJMJ9FCyZDlS00grGilL7aOCZAX12/VshNWXFoyvafmpVXRC8RTxkRjCGgqjpQhqvRiPjeHW5XSORvxWfgmIvOshELobojdlSUZh6eHIr6hMfYw1/l0Mzb+UtorSrBDu8wgwrCECDYlESocPveGeaLrQmqn5zcf5zYQ1ZjTvKtQ9qys4517FYFjvJ1qpkPKqyHIKyCrAHm/u0yxz16JtOEaFt5QPZ4dzUuCz6qf+x3R/LJIQJJVRUfIWVnLWYDNhhv5HtaWxF8Y2NftBTaRVMqpGideMjcu1fqrhKayohkvf082vcCIoKC/xUOWNs0iE6BqZx6JUL2heDg4cZN0568DfzMWdzxGafyMj8y8u+BzqWQjt53ZkqI73I6/y8XR+llBYFdjPssajcNn39I2RAag82/al83rn66TVNH9+0Z87XrvpZkIWLCHEn8z0RFxcXFxcXGYCKeV2KWVjkX07Zns+EyJrwZJZIeHFWAuvd+4iktZd7FQjCYZ0EGDXDu6grO89/c3ef4cG53TXYw1tYioGb/cdQbOkxO5Vxii+Gx3UY1TMPh361eds2Q9kNIk/lub1hkEiyhBRdSRrUUuFBoklksTSFmuBlGhIPFIhLVVCaorfd73OOy2DPPJ6E5rMiW6mcKtKyW8af8NzJ5/LxgNN5IJcNrST88LvAzAgQ0TSzqm+i1m6fIZQ3xfty7pQ5mMKwYHgDt29zBCC2/0xfnWg06ak+TR97nU7fshARwOMtvPz4+O4hwLDgRAH2kbY1TBIJq9WlGaJ/QM9eYcEajtH6A8ls5YiJxLDraSzMTX25CEFdO6nOt5sa5tzP8SiYNmP9RrJU8rUGCBIZKx1pMx+JP2hJO+1j2bTpGfPLx3nqBgmjF2J1lOTC5ozOcswaV2RDcbTdAYHaR16P2tFy55zJleQ2UoyrVKmGp+Pxj9kt/dEiif9GcsVUTM/5yOtRFMKcdMq2Vtra9czGuNoj34O80q9Rr9mH7m+ANRC/crWHmBPo26lOu4/TnPoOM5IajsDHOoIZK9EeYmX+cGTXN3/a8qz87B0HPMTlEkOxl4irAygSY10uB8S+tzTSobOEcvn/PgO2o0yFFHGKDYOEOym+8B/kujQY8sToV6ah4+RVJIOKf1njolasH4mhCgFtgLPSCmD47R3cflAMS0ZAtvfOkVr0odOfQ4uLi5nBMJw9aL2yWysTb6i4Oxsl8Mmrw2dQJzn9B0i8/5aVroN5U4A8YzKyf4IDb6cG99L0TbuyyQZTtqFTsDiqvihbJ+qZheXzSxnVkJxXXhKZlTSMidg+bQkgc7jJL0LERf+qd4204fmuQTIWXPeTvbTGa6kLnSSCs41zsV+/Uo8nqx4fXJUjw8iGcrWmHK6QkLaFYs62UlJ/0v8t8u/btsezPSwZ3gAHL5Q6CQcAAAgAElEQVTONakRSoXY0VJcpzcv/YLUAHARaBnwlBGIOwmV0i60hnpQFp1vV1QcaB2OIoHGwSgXeGOsKYFIIkMinCQS6+cN0cmVspp4WjISSbF4XglmxkdF1c1C5rjlSk4cPP7Wdg4uuZCV5R9hsS+X6CGazLC/dQRK9ePe7wpy6egezgsH9L6MdomMqiduqSicvQRDWdK3h2Wc0UzEtt/836MpdIZTDJEkmlSoKM2JuqONL9MhQni0Nmh5Q0/wdOGfIgRoQthHHTwOF1zPE++0UxfZzkc9ndy3+MoiV9VOc2sTi7Wc8l7czVRaXmkONk59f38olb3m9b0hfP2/4LrP/R9oetXWulSLUZkeYjCczFqBs3FsVqW16TVU7TL7UJkklFbaNoUS+idlb89e+kMJlvJp/RnNm2hBwhILXuH8TTWMLg+FlD4ODe6n+5hRy+vCP6VpMMJwJE32KRo6gbrwfMqVEFJA+0gM8qohmFNSEgFeirQSOvQ891at5803vkOvR+XsS2+m3OtgapwhJqRgSSk/JoS4BPgScFgIcRB4Qkr5+ozOzsXFJcvycM34jaYJm/vfRMjPQmgSCMGfffPUJzTXmO7iyhO4Rg8++CDPPPMMXq8Xj8fDY489xvr166d3HhY2btzIww8/zLp166bcx9NPP833v/99ACorK/npT3/Khz/84XGOcgFs2pFZBFSXYaRFwRIIBzG6nSAJobJazrf1k7V+SAnBLlh4nuPQmgQ63oH2fXDj1xFm8L9PMhCzx10kMnGOBvdQquasWaETb7KwvER3qyvzEkupSAkv1Oe+V7JWt2x8l74tYQTFe4RAkznrg6lAeWUGj+HuqEkNFa+hAGk2645XCAwDn8VFKnc9TQXLZ8Ti0F/8O09KmKcGkQ7uh2AvCNyWeId0qosV8+bZrI+gpxN/+sTTDv1LWkeGeLMvZ2XMphFXM9kU1YpM0SpyCs3SRAfd+/OKi6vOiTT8CT9d4S6uOfsa2zMR0ZLsinfTWtNCvMLDRUtaiYkMIZkirejCaEbVsu5u+QpAVbw1+zohFDSp0pc6yiLf+aZnK92jMQ60jXDNZYKRaIq+UJCuuD1eTko9ZXhIGGn886w5EkkslVNy/YTAcE4sURN6Qd4K3VLilRlepYOUR+MGs/aAgWL57ERb3iWUzLBIc4plMw+wLwLUJoe40nDTGyRGHOfrXaLmFlNlviIMxRM1SOPcPR5jnkZ7KdHQ2D7QxhIqqNBKoNM5RmpeJsyuE3brsW0IgN4a1IrL8KlJypQognKUQ7+k8fJ/tLXVHOYZTA9BWcHmgrEEMESc0rzyZPpOfaNXKnikQnvfuzaFpCTlL4gxVLQMi5PdvCMUjoSHYaFKUOmmL7iGi4F4WuUXL+/lmvR+UqpGIB3lVwc6SQQjRESaZZqC8M1eMvcJx2BJKZuFEN8CDgM/Aq4Wuj/Cv0gpX5ipCbq4uLj8sbN//3527txJbW0tZWVl+P1+0ulx3CTmABdeeCF79+5lyZIl/OEPf+Dee+/lvffeOy1zEXqRneullM5SyRym1OshjqoLO1aXQCDnXJej3qNbmaR0LnacivjxjHRQoiRZFltGRWaETFnO4iClpLdzH0u0DBWZOCpaLgQkTxY96q9HAGfFch6Y4aRCZZkvu3Kt9wmH/a9zdr6IbgnCSio5AVoIPYOgEx41AyndeqEJUxA1XNPQlcdy5RJkyTm2K2MKuSVegSk2ezzOkRLdSpTBof9EJZf0wUkM31y3mdTwJ0lpMcJKv31npB8ycSMzn0Bd2KXH0Hjs8TMnByI8UfMOC5fkLHbSjODQcq5sXcnDVAq7wO8bqre972lroDSeAgExMpSTS+YglQxXL7o4ex0WJXvpVFvolCPEvQNUsAbzWdLyxfJssWZzfoUWJhWpK8Ca3e3Suggwr28/56cG0Mpy11VKfTyB1dKR+1+3tkVY4q2EJfoxHkuEy9mxBsqP/oqRq26nUYyyRi7JWmFLho/B4deybSNp/XouTg1xpNdH2KMQC0QM4d8pSNBrqzd2IDmAGtSf9fc8/czDi9Nn0IqiWq5WOg5dB9jZayivMmFJ4a5B3TMQ7KLzyv9hJNHQjxwKp2jyd1OTHKZazOOj8nzof99xPFX4iFuU0fy4OvOvqmlc6s9dm3Z/jGeOHKDSWORYmmhnec/LxD/8Rf1SKCmu79rCb6vLWLX8XBJplQXl1qQoZmIdc5ugQfi5Or+sQnwU2vYCUJEZYXngbTyZMNK8pWqGBU2/ZklCtyinUOjIhJFaNR6pogiJXyTQ4gfxBF7Er9zFxUBvOM4V4ZcJo8ICEBIWJnuzn3XNsZDzzDHRGKwrhRCPoGdiugn4tJTycuP1IzM4PxcXF5c/evr7+6mqqqKsTF82rKqq4rzzdOvDAw88wHXXXcef/MmfcO+992Z/xDZu3MjXvvY1brzxRi6//HIOHTrEZz7zGS655BK+9a1vAXqh4csuu4wvfvGLXHnllfzVX/0V8XhhJrbXXnuNj370o1xzzTX89V//NdGobq34xje+wYc+9CGuvPJKvv71rxccd8MNN7BkiS4RXX/99ae16LCUUkNPzX5GYBX1FDQ6CFksMbqg4CdBbfi3DKZPZttarSmK1Hj0wK9pDASNPgWBWJoTfSEGIimQGiuD77E00WEbWwIvRZrZHm0FTUWxWIdknoZVM3zEVo+p8DzM4Pniws0AMWqTwwxYXAY1CRKLwmURYD3Nr0NfHZrUa0FJwGMVnmLDXDKkC46D4aQlwYC+22tRqmTH29BTo7ugWbLtvZXoRaS7CywUQZL8KnVCn5NxkzJagobY7+lO1hhzNVCSEO6HQCcEOlCaXtWLsuahZ/zTLTTZuWYtWEq2T1UWpjf35CWOsGYMfMfTSxqLm1hfLeE9j2TT0a8e3ctoOMhQNEVaJmlPvIvWp5+DBEoSVpdJI912SqE7mCClFCq/aVTOjRxleagGsNRTMzNCtr+FlFCRCeDkAmh7LXMKs2bE2aUsKdS9DuLrs62v0ihGUSy9+UYabRk3Gwb0lN8eBO959PMbHexlnv8YEl0ZSasa4aR53yXtCT3tuaJJgokMmbg98Wi+WqZYgptUNEa0RG5GvTXQupvBoP6sZwxl1B9NcaQroFuWgRdqe/OySUoC8WTBtbJSlzSsy0LYEsqoeYqF+R2hapJSSxxlJJ2hO1lryyJaGe1ka/3PASiN6HNbFDpBT0cLx3rDKJpmM00NRSzKdV7ZB7NV9PBzjEaT2e+GEjVhTwAjVXYm26nI6LXoDoh+Xol1IDNRPBarti+jZw1cEDnBCAl+72nDTyK7OFCRCfKhod/nvoOkOqFEM9PFRC1YPwZ+jm6tyt5xKWWfYdVycXE5RcaKzzpT2Rw8WjTDILhZBifKLbfcwgMPPMCaNWv45Cc/yZ133snHP/5xAO6//36+853vAPCFL3yBnTt38ulPfxqA0tJS9u3bxw9/+EPuuOMOampqWLp0KatXr+ZrX/saAI2NjTz++ONs2LCBL33pS2zevNmmLPn9fjZt2sSuXbuYP38+3//+9/nBD37A/fffz4svvsjJkycRQhAMjh2a+/jjj/OpT31qJi7PZHhNCPFZ4AU5fnGb00fe1I5JP0c9IyylBMPhSN/u8VMVTzIsPMAykkqSHb2d2aXTfi1ObV8jMtTNpUsWI4RgIJy0KSvmKyFVru/aQufi9QyEzgIhiGpp0FQymoIQsDA1wFn+A3lztfZSiBm7lK3T5CDgHPT0szhTUmBtsr5flmgzjtYQ8aCeTEHYLSkBmaRK6osQpoBpFXbNlj6vvrNlMIwyMoxvnsQfTZFUCtv2CGucj6RXFCb2qI/+zn7OBS10MlKDmN88QV2J8Hgp8ZoCoMSXd/S7x5pJLNJYOt+5WHB+Msf8zHwaeja9eFqhQknx2+F2rjCzT5K7fn1KB0kW0Y1+flLY7+mADLGAChJp/bqHE4WucQrW62d1TTXiiLrbCaQzltK1ZlsYVhN0CSPFv+kmR/EnK1/BklJyrDdgGLhyo4vRFnbFc7WkzBT8XgRDxBDAopS+8BNnKSOWWLcF5bm7UaaEGUokUSUkh9rsIXZmnJ+UDI/8nO/87rLszI+KYY5EM3ys4rzs+QxFklBiP7OEopE4+gLywirbdutntSAbZB6/HtA/IxUZPS2+P93GYLqBj+SV7BtKJNkX60XxSeJCUit6+bRcRdzBvdRqKy9JjiDRCzWXhDtgwTKIByEdw6mcoOEhmlXYJJBSVA63j1CmRJEWN0N7vOhJrERJs4BySvoOO86uXmun1JjjkIizUi5ECF2Jhpwlaba/8idaB+s29OQWCdBdLYQQehiilE+NeaSLi4uLyylRWVlJTU0NW7Zsobq6mjvvvJOtW7cCsHv3btavX8/atWt58803OX48l+np9ttvB2Dt2rVcccUVnHvuuZSVlXHRRRfR3a0nJlixYgUbNmwA4K677uLtt9+2jX3gwAEaGhrYsGEDV111Fdu2baOzs5OFCxdSXl7O3XffzQsvvEBFRUXR+e/evZvHH388G491GvlfwG+BtBAiLISICCGcCzedTjQlK5wAWStEIq0SjNldQ8tUu8B/3JNbXX8jbS9Smi9gyFiurRnjtDL4Hm+eHCKpGG0bdjBPKyGtaMxTAgRUp8yBxQWX6kAtIAuK7GYFYMt7ay96jEpui1dLZ+c5EkszaKyUe6SexU0iSKRVhhTdAnuSUbxaCmFJXGCev8/QSoLhsJ5UAWzKlZVmkcuMZrbIn1sBeUpPOKlbfepTuuBLIgD+5mwMjTnH/mCSnqEASxId+DTdWjFw/En+q9WMcS0cMz/1uSfPrVICPYE473eHyKhaNgNbnAwjJPJaFnsHJ8Qgg+Ss257u/ZQr9o9Ov8i5OF469EpBbxGLdTUhFQ4J3aVyJJriSDKvppJxHvlubQtSZmyR/byfCzU7zj2T59iZNCySPV6lSB2xvFmn9Oe9Ot5sZoG3jSzQY+GQEq9MI5MtaKE/MN9QckYNB7WM1O0qyYxK63CMsJpzJ1UtLnzvtY/yfk8ITerzM597AJEKZcccDyEEXclDpLSYrphJe191KT+VkRaaxSjDIkF9YtSx6G9B7FjekyGCHfkj248nVxtNSknTYNTxozMcTpExF0NSEds+gSCZUbOFkXNzyt0QUxVW0AwLVs5WZbVgzSYTVbB2Yc+JU2FsK4oQ4pdCiCEhxDHLtqVCiNeFEM3G3yXGdiGE+JEQokUIcVQIcc1kT8TFxcXlg4zX62Xjxo1873vf48c//jHPP/88yWSS++67j+3bt1NfX88999xDMplzszJdCj0eT/a1+V5RDNejPAGtoFaNlNx8883U1dVRV1dHQ0MDjz/+OD6fj4MHD/LZz36WHTt2cOuttzrO++jRo9x999289NJLLFu2bFquxVSRUi6QUnqklCVSyoXGe+fMBacTtdAVDPRV7uahCFqyUCcsujprbD7QPsqv3m0jGM/kWTmchaqsvhHzEx5RiWdUFsky6o26zJboKsdEG9b+F6QGCaSGHUerSPRn21n3KkWUGIEkZbjBlahxpOkQJPQU7yYBklzd84T9TCUcE8NUhuphVF/tzxcYTUIO6dtNu4iiSkq0BL74ILWNHQXtNGm3NphKXDBjCLjJsJ7GHiARNO6d3n5JojMrnAMc8gwwkDqBEIKEWmglzrdglal24VRDZt3drIkVdosu3vH0Fr1zEsmopVgxUpJGzV7Pno63qYo3W3fbWJjqy77Ofz4kkqOaP6uQZTRJwDKWhgRDIY9rmWwCFIC4OkQKlQWpwVwhX6AtEcmL4dJR8tzj0lmlqlBNsc4ZIBBLM9ReT8FzmPe2TI3gU+OODczny4NgNJYmEM3QR5S+ZC52ri+kK7rHxDAJFBIZlaF0k01ZLlOjzG8xk6CMr2JVpIazGR5VqSL9rZwTyY0ZTSnUD25m1IgvS2U0ToiRsbuX8IqnzWkzIVLkx+WZ3WharlC4LzHCPu1YVtEFaBchEkJhKJoqLIptmW8+ImvlFPgMdSaFkrWkmhYsa6H02WSiCla5lDK7bGW8Lr5cqbMVyP/F/QbwhpTyEuAN4z3Ap4BLjH/3Aj+d4LxcXFzmIN3BRPbf/tYR2z+XydPY2Ehzc06YqaurY+XKlVllqqqqimg0yvbt24t1UZSuri7279djDJ599lk+9rGP2fZff/31vPPOO7S0tAAQj8dpamoiGo0SCoW47bbbePTRR6mrq3Ps+zOf+QxPPfUUa9asmfTcphtjMe8uIcS3jfcrhBAfOd3zKsBXhjj/Wqujk2237DpQeAyFAkQyoyKQDIk4/USN2Bg7YcOiIYRVEJVZH7uYliEpckJPgZCSihRK17Y56dan9qjd7addhNCQ+EYPmCPaxLODgWdJaTE8UsnLoCoJxXVh6+zYCSpS/dkp5/OKpw3V4lIpgTYRYqRvO4R0wdKM1/COIbOmVWfXt8xQM9XxpqwVzUpvKJldkTf3BPNSrcfSCj0nD+If3E9EGSSfTnKK9EC8G0VOPrGNhnS8P6rhAtguQsYW053KOaV2vpK0L2ZXRkLJnII0jn3PNo5JJGU/Hk0loSk0GrXGzP5qlHr2i15SKCxLtBWxQuWstUOxFKqUhBKZPDf8whnmW//iGdVWsNlEABnyLIWWQt7CotSZz5dZ6gBgWOSSV1jpFhFe83QQJElfqt62ELI42V1wSCKj5qw+eVT1PJXN8Hho+E1SEbuFMJLS1ZC05XPfLezKOUAok0FNJ/QU7gYZVbPFaWVUjb2ebsunzE5FVPeWkBKWxFpJCIXXPB20SnsdOU2C37DQOz0/+fFTo9pQdk+f4bobI8Nbnp4CK6PepyxYQJxJJhqDFRNCXCOlrAUQQlwLjBkwIqXcJ4RYlbf5DmCj8XobsAf4v43tTxo+8QeEEIuFEOdKKfNS8ri4uLjMAWY59Xw0GuWrX/0qwWAQn8/HxRdfzJYtW1i8eDH33HMPa9euZdWqVVx33XWT7vvyyy9n27ZtfPnLX+aSSy7hK1/5im1/dXU1W7du5fOf/zyplL6qv2nTJhYsWMAdd9xBMplESskjjxTmO3rggQcYGRnhvvv0WDufz8fhw/l+9LPKZnRPr5uAfwOiwE+AyV+4mcTjhUUroHt/gUJzTvQ4cmmh0U1i1CiykFQ0PEZihBYRYHlmAZTZj20QukXKmxnBXHOVFgFxW/gE3SJhjGHp35RThk4gygrra9lq5Rgp0xekBijR9L66RJhKWUJNeFR3j8k3Emgqo5lOlkfsmdKElAQtmfTeVd7j4yyjQhm1tzP+vpnpoiMrOBoughZh2BSAPUKgjrPCnUGlTYQozVubPi9SR6ykisC8lTaBv3UoymVOS9EJfa6hRAYJ1LQ8QazcLF+Qm8P7nlySiZNDA1Mq1dErIpQb8TFjnZ1XS7M8XINnoT5hazr4tJJTFNNifDerpFDolKHs+8r0EKVK1Cb0RpKKbd0grqi5TJVSQriXSDpNIJ6msrwEa+bMsEhz3LC2nBPJOkllu0tZlJ+TwRD+khRp47NhxpaNlZgln3zlUiA5IfR7mHNWc+5P5qkdqpAoUuP8SB0j81aTKCmMXdrn6QHOLthuVdYA6rqDugXTQWd4Cz2ubFGyh4YBhcuLpO/PpyIzQkk65zpcIwZQOnuoLPPlziH/s+pw6lJCLK1Sjo/l3Tv1bZgeEvoBUZwt9YF42lbewMS0UmWvuUUhHhS6BTEscnX0YkLR14rGON+ZZKIK1v8EfiuEMJcszgXunMJ4Z5tKk5SyXwhxlrH9fMBaqbDH2FagYAkh7kW3cnHBBc4paF1cZoxi9Z5MjAKYLpNgvJpSH8Q6WpPk2muv5d13nbOLb9q0iU2bNhVs37NnT/b1xo0b2bhxY8G+jo4OPB4PP/vZz8Y8/qabbuLQoUMFbQ4ePDjmvH/xi1/wi1/8Ysw2s8x6KeU1QogjAFLKgBDCOXvAaeatFj/VkSQHhkd1QcGaZAuHlOEWIWee9JEQukB1TlSPyQuIFMud3AGNv0vjbTSylNUsplQNE5e6hWVeSS6leKYweTdAESuCGeAuqUwPEU+uYFGq17ZfEVp2lT7f6lGiJZAOTjYCjQ5b2JwsmFUkpWTdBQelpR6R0SyWUJHzcnMYiqRs7oXFGDUVOweJskIJEGAlixO5uLf56SGoyEVXZDQ9HsWj6kKgqdSpUjI/7SftrSiIqTMZCjtbr5yER+u2JhFgVCaBpfo5UFjcGaDUcHHLtyyZ8wY4IfyoQhvfPAWc9NgV3qWJDqxqrsewR5hdWe+hRMKJnRzrixD1qJT6PKhIhii0/HgtmRXNme/xdGefnAyaYxFcYRT1nYrTmEBk52u6tMWVkYJsfdbzsipgZqKUhal+EiWLCo7BOCIf8z4IYIQEoyRZIsvH1CAWpAfxaWkCSqHLqxXzei1NdCAtqdeHRYLFlDgflIfV3TaUyLDA+JyMVYTYiWi6UIlPKhqVDm0B5mf8BduKuRrOuSyCUspDQojLgEvRb+VJKR3yhU4dpzN2vCNSyi3AFoB169bNrkOli4tB0Yx/R15jxeJ5zvtcXFxONxkhhFm4BiFENc7ljU474aTC/JSzkKC7fdm3DYViLPLrFg5TuSo4TgkAOYEu/we00TOKKjVWBA/zlpKiVM2AZZU5IRQ82YPGFlSC8XS2aDDA/GhnoduZFDblTLG4KwkpOSt+YswxAHxaypbIoCeYsI1jjd0w3akE+gq3eQpObmBO5ITlQoRU8UjVFj/l1DYQz6B5oNxjd7ucnxlBUNw65FxSGlvqcpNknpXJLxIsD9cQL/GSGUdBSsrilo600PAW3Ts2FYFdtIlcnJJATyVu6qoqVrc6Xbk3lZLReAYEHBB9CLmk6BgZrVCJUPM+3h1Zt8tiNqdCRL7iJAoV0e7kfqpjTQWfClVIW7FrKyVanMr0kMMeva6cv2K1bVtXMA4eva93PL2OxzkxTwmM2yY5xrMH+j1xLLWgpnKWLdu1zimWfanYxAOSHDC/F3KWK/3V/HShYmUf3SR3V+aiiyDoLhSrjGOuFkIgpXxykuMNmq5/QohzAfPJ6gFb5s7lQF/B0S4uLi4u08aqVas4duzY+A0/OPwIeBE4SwjxIPBXwLdP75Sc0QVqmX1tpWs0RmmZVWKRBHubyWTGXqXWgGXxVuIlS83DdIHbQrGkD9Y+QI9LSqsapV4PJelQQbtEnuCfTKUpy2szIhLMUwz3wzyr0OJkNxPlsGcg+3qs2ccsK+P+WBoj4/SESVBc+Uh5F7AgNVCwPZywr0UrmkZSUWyKTFLRwDu19eL6kaAtPfpYxIus6usYAvE403BINDchDnjsIl2LCNhk7qDIPbsSSbsSolHYrWDWeeZT7vPYlHXzPJS89qaL42TE7PPy3FQ9FOoL89P+MZO9tAwXxjcBVKSdzhFK1Rhezf7sREXx+KRTJd/imM9QOFlg5a2ONcPCnO1ZQaNETTIvEyJaWg1AXCjsERP/LI+F+X1oPoOmu7Fj2zlgfpmQgiWEeApYDdRBVs2VwGQVrN8BXwT+X+PvS5bt9wshngPWAyE3/srlTOVMqmc11lynyxJXEDfQbuk3UCiYubjMFFLKp4UQNcAn0GWs/y6lHN9MMsfQpD3l9fz0CPMzw2McoSORzFOCOQXLQVRrFyEWkpqQEDcYSbFi8TxK06Ei0RQ5nKwzflH8+2cs4WmqmOfULPQV/QXpQoVoLEzXLqdro3jKHPeMF9dVMDkHlodrSKk4Omrt9uguiQtladFU82ca7yUHUTTJgCh0Z5ynFK+3Z1VIcjFP9oQMFLSYGp48FSs/QYYVVUr93kzSeHJ2zP7VlFv8mF3tIZVRHV1oy9QwaOU2ZXZJUn8eJ/L5nawtaeI2R3tbe8KLuWfBWgd8aDKFGYUQz6IntKgSQvQA30VXrH4jhPgHoAv4a6P5K+i1tlqAOPD3Ex3HxcXFxcVlIgghnpJSfgE46bDtjCH/h3hJsnNKxznVvQEIy8kqN+OLBqVqbMz9sykymhn0KjLFV+0nOx9rFjkrsbx4krGlqOLC38JE65iZxXx4GI46x1dNBNP6Mhv3YbwxpjIHKQtT1ENxi2yB298kEGqacum13S7PGP0pqjaGcjBxgd90wZvI9ZnO+5hvjS6G1UVwMklEJkpAjG2ht3GmWLCAY8A5OCSdKIaU8vNFdn3Coa0E/mmifbu4uMw8s2HdcnGZZa6wvjHisa49TXMZB0GMDL/ztBTsmcxKbl6XAJQYRWz9sZRtl8y+nnj/oUSGqOIrcP+bNHNAILLSM8b3n9Nas5AaFFGyrGQ0CcLpZMe+AKVqgsQYsviog7Vnasz8jRj3+ZVMycerXCn0iLBmQ5wuPELgLVCMiisho/EMxSLHSrWxFx6smMpiflyZE5NNLDFVIkkla02fyLyszOQcrTMx71QwnqacSShpp8hEFawqoEEIcRBys5NS3j4js3JxcZnTnElukC4uQohvAv8CzBNChMn95qYxkibNJUwBPulb4Lwq7yigj49Zu8cpO59VIPFpiaLZ7PIJpxTKcI4vmQxzTL8aE6e5OmUyK4aTO5lnnCQDUSZfA2symNa82YhdySUrcN4/FEmybP7kkntONFHJdJE/9fE+L1NeFLFgxtpNypIzJSY+V6ursjXerUwtLIae3/dEMneeKtbslP2hJLHYxL7XpoOJKlj/OpOTcHFxcTmT2Fy3eVr7u++q+8Zt8+CDD/LMM8/g9XrxeDw89thjrF+/flrnYWXjxo08/PDDrFu3bsp9vPTSS3z729/G4/Hg8/l49NFHCwoZzwZSyoeAh4QQD0kpJ533XwixAj3m+Bx0XWSLlPKH0zzNLJo0UjFXXFQQYH8qtGWLyubSPTvRk2ma9eIxYydgmGPMgFzoc8iAZ2WqSvVkyXdpnAnGUzYymmQgMnuWhslyMDXIgGdyFjJApNwAACAASURBVMPUOAr0RDBrPM1VmhyTkpxedolOFqEr65WhJjzLzpm1sSeapn2vEGIlcImUcpcQogKmnK3TxcXFBYDNwaPjthlf9fjgs3//fnbu3EltbS1lZWX4/X7S6bn9YwvwiU98gttvvx0hBEePHuVv/uZvOHny5PgHzhz/RwhxF3ChlPLfDM
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Again, based on the trace plots it looks like 2000 iterations is massive overkill

"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAEYCAYAAABBWFftAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzsvXmYHNV56P073T2LRqN9hlVCAoEAE2EWYWHkYAUbLiYx3NhJsJ/gx4kDOCb4ftf3+sln53qJifj4nI8Y7NiykY2RwCy2BQgj4wACLSwSkmY0aKSRZt/3nul9r6rz/VFV3VXd1bNpZjTC9XseabqrTp1zaunu9z3vJqSUuLi4uLi4uLi4uLi4uJw6ntM9ARcXFxcXFxcXFxcXlw8KroLl4uLi4uLi4uLi4uIyTbgKlouLi4uLi4uLi4uLyzThKlguLi4uLi4uLi4uLi7ThKtgubi4uLi4uLi4uLi4TBOuguXi4uLi4uLi4uLi4jJNuAqWi4uLi4uLi4uLi4vLNOEqWC4uLi4uLi4uLi4uLtOEq2C5uLi4uLi4uLi4uLhME76Z6lgI8UvgL4AhKeWfGNv+FbgHGDaa/YuU8hVj3zeBfwBU4H9IKV8db4yqqiq5atWq6Z+8i4uLi8tppaamxi+lrD7d85gq7u+TM5rUUKVKiafkdE/FxWVuEBnQ/y445/TOw2VCTPS3acYULGAr8GPgybztj0gpH7ZuEEJ8CPgccAVwHrBLCLFGSqmONcCqVas4fPjw9M3YxcXFxWVOIIToPN1zOBXc3ydnftP4G/wJP/deeS8+z0yKIC4uZwi7H9L//tk3HXdn1AzRTJQl5UtmcVIuxZjob9OMuQhKKfcBoxNsfgfwnJQyJaVsB1qAj8zU3FxcXFxcXFxASsnhgcNE09FZGc+f8AOQ0TKzMp6Ly4ww2g6aNr19Sum4eWfbTp49+ez0jjVHkFKias7nfaZzOmKw7hdCHBVC/FIIYarj5wPdljY9xrYChBD3CiEOCyEODw8POzVxcXFxcXE5Y8moGRpGGhhJjMz4WKPJUQ4OHOTVjnG98qcVWUSYdHGZ8wQ64P3noPOd6e23yGeiP9Y/vePMIXadGOJHbzSf7mnMCLNtn/8p8G+ANP7+B/AlQDi0dXzSpJRbgC0A69atc7+hXf5oeOT1pqL7vnbzmlmciYuLy0yiobGnew/Xn3s9y+Ytm9GxpPFT61qUXFwmSMqw9iYm6qQ1QaTGH1vuuWO9IUBfcBHCSRU4c5lVBUtKOWi+FkL8HNhpvO0BVliaLgf6ZnFqLi4uLkXJZDL09PSQTCZP91Q+cJSXl7N8+XJKStykByZl3jLm+eYRTodP91RcXFxmCzm2y+EHUQkxUTWJzzu1c0ukVfY2DXHTZWdT6ps7CuqsKlhCiHOllKat8y+BY8br3wHPCCF+gJ7k4hLg4GzOzcVlLrC5bnPRfbVh3V3omoV3ztZ0XAx6enpYsGABq1at+sD+wJ0OpJSMjIzQ09PDhRdeeLqnM6fwCA/aOALXdCAMBxLp7DRySiiagkS6GQNdpoyUktZgKxctvgiPmCXhWUr9n2eWhfVZ+LzPNTxCoEmJokl83qn1caB9hBP9Ec5aWM41F8ydRCAz9vQIIZ4F9gOXCiF6hBD/APy7EKJeCHEU+DPgawBSyuPAb4AG4L+Afxovg6CLi4vLbJFMJlm2bJmrXE0zQgiWLVvmWgYd8AjPGR+ntO34Nn5+9OeO+2ZCoXP54NEUaOK1zteo99eP3zgd1zPyjbSe2qAnd8Le759aH1NhPAvWB/Az4zW0kOlIdDHXfp1nzIIlpfy8w+bHx2j/IPDgTM3HxcXF5VRwlauZwb2uzniEB40ze0U7pabGbaNpkrSqUV4yxeVrl5mlfjtUng0X/ulpGT6uxAEmluUyakShdB+EZaunPujAsfHbTBJVUznQf4B156yjzFvm3GgCLoJzTos4RTweAapEPcMXk5xwi1C4uHzAGMvN0OS+q+6bhZm4uLhMFYGYERfBUCrEorJFtAZbKfeVM883D9CFt4a+MEd7gnzuIxdM+7j5mNa5d1r9HO4IcP9NF1PidXaqGYgNUD2vGq9nZpQwRVPcmlzF8Dfr/06TgjU1F9a5J6w3BZp4f/h9VKly4/IbnRvNcRdBTZNseauNGy+p5kPnLZyWPr3GApuqzr17dqq43yguLi4uk2SsjI5TYSJZIB988EGeeeYZvF4vHo+Hxx57jPXr10/rPKxs3LiRhx9+mHXr1k25j5MnT/L3f//31NbW8uCDD/L1r399Gmf4wWYmXAR7o7281PISn7jgE7zR9QYAn7vsc9n9rx4fmNbxJkLrkG6ZCCUyVFUWruz7E35eaH6Bq6qv4obzb5j28Qdjgzzf/Dx/cdFfcMHCmVcsp4v+UIISr8fxmk0XUspZNZh0hDp4pf0V/u6Kv6OipGLyHcxha7ipICqaMkajueMiGE5mWFDms3kYpFWNRFplT9PQpBWs+uF6IunIjHyG5yquguXi8sdA+1v294GQ/X2RCvIuc4P9+/ezc+dOamtrKSsrw+/3k06nT/e0xmXp0qX86Ec/YseOHad7KmccHkBLTS2LYDSloGqSRfPsySVGjbTSg/FsQt8ZTXIxFjI6CIvKdBchQCuiTMYzuovYSHJmaoINxYcA6Ah3cMHCC+gJxIkkFS4/dyHhZAYpKbiOc4HnDuqlQ0+5RIemQmQAFhWWHv3p+z/lolgnt85feWpjTJBjI7pr3nBimJUl+pjmc+mZTMqA6VqYkHLalLYJfc5Og4IlpeRoT4hLz1mQddMdiiR5+kAXGy+t5mpL0gjFiJPyeRyuiZKCwWNw3jWO1+ytXl0GKaZgndKZzVHjl6tgubicYdSGf12wbXPd5GrlbA4etW/Icyt0XQjnFv39/VRVVVFWpq9WV1VVZfc98MADvPzyyyQSCW644QYee+wxhBBs3LiRq6++mpqaGoaHh3nyySd56KGHqK+v584772TTpk10dHRw6623sn79eo4cOcKaNWt48sknqaiwrx6/9tprfPe73yWVSrF69WqeeOIJKisr+cY3vsHvfvc7fD4ft9xyCw8//LDtuLPOOouzzjqL3//+9zN/kT5giJEWZKAbVt4M8yaXGevn+9qAqQvfM50OWqoZAu88QeWF1yK42thY2C6WUhiOJoAJxuql45AMwcJz8waUEOqGRSsKhL8Sr648ZbQMmtR45O2XObv0Mi4/98M8/lY7cHrqDHaPxjl/8bysAmqj9U0+NHiEhrM/feoDte+Frvdg3Zdgwdm6oNy6G+2ijQC0ZeyLcf6En9rBWj658pMzltXParmd3AKAea1OTeJWpIYHcUpZ4DpCHZR6Szmv8jx9ZuZzN9bUxlMMJ3laUkpah2OsWlaBr5j7bTjJmyeH6A7E+Ysr9bmGE3pNvJ5AwqZgRZO69c3j9Fls+i8YbID51bB4eizBb/e+TTwT55ZVtzjubxmKEE3l8uHNtXjeuZMw3sXFZcrsbx2h+8hrjv8KrFcuZxy33HIL3d3drFmzhvvuu4+9e/dm991///0cOnSIY8eOkUgk2LlzZ3ZfaWkp+/bt4x//8R+54447+MlPfsKxY8fYunUrIyO6RaCxsZF7772Xo0ePsnDhQjZvtivbfr+fTZs2sWvXLmpra1m3bh0/+MEPGB0d5cUXX+T48eMcPXqUb33rW7NzMf5I8MRH0ZCQSRBOZmjoc7ZmNQ9GONwxsYKnpgDiJMBamYaEXmMyEknSPBShr6s5q/DkDymlZMu+Nl5+v89wVZuA8HTkKajZWrh9+CRq7VPQX1ewyyf0deaMlqEr3MVQuonuVI2tzdHho/zy2C8ncmoFPH3iaV7teHVSx3SPxtle08PBYve16z0Wpvqd9wEcfxE63p7YYFHdgoeZRKL7Peg7gtLzXlbgD6hJyOiK7hudb9ASbGE0mZtbR6iDp088nbU2ZqnfDi27io+taXDiZYgMwu6HEMFe4BQsNRMRsDNJPdNg35GiTbaEjvFSrO2UYqJeaX+FHS05y/10WLAmS+twjJff76O2K1i0TTKjj5lRrWM7W5VrOgMARJIObo5p496P5QLpwFh3+ujwUVqCLUX3v/x+P7tPDk1qvNnEVbBcXFxc5jiVlZXU1NSwZcsWqqurufPOO9m6dSsAu3fvZv369axdu5Y333yT48ePZ4+7/fbbAVi7di1XXHEF5557LmVlZVx00UV0d+suRitWrGDDhg0A3HXXXbz9tl0wO3DgAA0NDWzYsIGrrrqKbdu20dnZycKFCykvL+fuu+/mhRdeKLB6fRARQniFEEeEEDvHb32KY0lpZA3zsv1wD68eH0BRCwWwnUf7eavZP7E+J2gNmOn08GkjoD2dzljG1P8eHznOQGwgq+SZc+kMd7K5bjNpdQzX2Lgh9Gv26zQQaOOx0DG6R5uLHtoV7mJXl64MaJYqMUKqvH3sVyQnksXOgVAqRGtwcmnDQ4YFwfw7aYZOnvLCmtRUTPH32UgTpPTzd1LSa4dqCaVChFJ5ruf+Zug+VHyQSJ+esa/uV3rfI/r9cXo+zW31PSH2Ng0DkFEzdIe7HSY/xvNrut32HLZt3ly3mX09+7Lv+5XY9LkaWqc2WQVL07LXfrKK53BEz+KZcfjeMFGNz4rXUvPLNJrmn/7KZfp3/EXV8yc1j4lwppekcMJVsFxcXFzOALxeLxs3buR73/seP/7xj3n++edJJpPcd999bN++nfr6eu655x5bTSnTpdDj8WRfm+8VRV9pzHeryH8vpeTmm2+mrq6Ouro6GhoaePzxx/H5fBw8eJDPfvaz7Nixg1tvvXWmTn0u8X8BJ2ZjIAF6knbhIZ7W79UpiyBZwWnsngosWDE/DOTVIcokslYNJ6xj5Ne4EcI4F6lZnLr0Nnu79/JC8wvZ1fN8oTKWiY05d31Ae3r4vqSugHanCi1C1uQDpvJmHbMq1gLBbgh0zpoQaJ67Y6xLzKJMT8t88sYQHmJphUQ6bXvgxjr3jKorgpO3PBljK2ljaCNqxexm4BgiEbKNv+vEILWGJWVfzz5ebnuZQDJQeC5qpkDR1jvSiKcVounCfcf8xwraThdOimkkHbFnCnUar32vbm1Lx4pf33CfHkeXhzmW1+k5MnCKqzJdAPMtWOZPg8+xALPZdm656Z1OXAXLxcVlfHY/NPY/lxmlsbGR5ubc6ntdXR0rV67MKlNVVVVEo1G2b98+6b67urrYv38/AM8++ywf+9jHbPuvv/563nnnHVpadFeNeDxOU1MT0WiUUCjEbbfdxqOPPkpdXaH71QcJIcRy4M+BX8zGeB6kLlAJT1Y4K5YIQpUZth7bSluobcrjWQW/gnGO/hpO7NTjc0zeflT/h56+uT9kV7bMpBShRIYfvdHMYNhSTNrs3iJQlnS+pdcvKpiDXayckBCfZ+WSxjjCIWbI2VqSm9fq0X3Zuc5E2nwnTKE3mnJwtzJc+lIoSOlszatLDbM5eBRVUx33j4nwcLQnxO6Tg1g1LJlXk8163TKarmBN+voIwfH0KK/Gu4yhffa+T7yMPPGS46FSSkJpXflKKvmFyiXsexjqf+NwoMb7PSEOtAcK9xU2HreFJrWiSn9a0bKLI2aSDvPcoukoTzU8xcGBgzSmA7wULeKSGO7T/xpKLIFO/Tc3nHMRzRx8gtSBwq+liag8ilqohHmyymBef6ZV2XJdGkcb8ScKLeia1D6QVqnJ4Ca5cHFxcZkksx30Ho1G+epXv0owGMTn83HxxRezZcsWFi9ezD333MPatWtZtWoV11133aT7vvzyy9m2bRtf/vKXueSSS/jKV75i219dXc3WrVv5/Oc/TyqlC9ibNm1iwYIF3HHHHSSTSaSUPPLIIwV9DwwMsG7dOsLhMB6Ph0cffZSGhgYWLpyeGiqzzKPAPwMLijUQQtwL3AtwwQWnFugtpCnGWITcIvKKIpPElTh7uvdw0aKLivc51YyBpgCbSYCvjKFIkmpLIox3W0c41DHK366/gLMWlgO5wrCjsTTnAT2BOGcb+4yzQUgNIcCrpSjtPQAjJeCNQWmlozCXnY6msLd7L9efdz3zSxzclfIF1bEULKeLajteklQ0OgfCHOkaZd2q6sL2Bj2BOIFYhrXLFxVtMx6KppBSdGG6bdhBcDeu+aueDmKJd4G1eXOX1KR0F7r07v+HeSs+om9f898mNgHjGiXTGezPnnFNkiGI6MJ9SlHRtFzq8fGeq3gmzonRE1xz1jUIIYinVXZFuynxCKi4AESROmfJEFK1u0u+3jCIr1wXYRVpKKIiz7dttL2wr6xldHxLy2hilKULzh2zzdu9b3PMf4y7195NqbfUtq+mM8BjgTa+dvOaAgtWQtEXJLoj3dTGDTdHRwXVOBch9GP9RomQUE82mUttVwBNwkfzjjQXKRwTpVjalCpRlox2AOeYQ9mOH4s3ut4AJc19kUR2ngA/e/9nXLz44qIJKvI5FVVstjOgThRXwXJxmSYmUuDXxWUqXHvttbz77ruO+zZt2sSmTZsKtu/Zsyf7euPGjWzcuLFgX0dHBx6Ph5/97GdjHn/TTTdx6FBhLMXBgwcLtlk555xz6OnpGbPN6UYI8Qcp5afGafMXwJCUskYIsbFYOynlFmALwLp1607pV19gCA5SUq4EKU0l0KRFefI3w1ADcDnmGvV4K8ZOCpb1tRC6/JkvWGmaRlRLsxCJP5ri6QNdfDac4IKlekyGP5oiqgzxTs8wt6y+jvllpVmrhgDi6iigKybRpGIZUxcor+vZhjzXULpHWo1zyc0QmT1FPTNasJXGQCMaGjevvHnMcybQiUyb8UPjWLCkpERNIL3280+kVSSwp2mI8xZXcp43DPOrwFL4OKNq/Paw/qyvXarq8WBVF485tdrBWgKpAJ+44BPZbU+feJrGoWGqud35oETO8hJSerOvg8kgg/FBLl2sj6kBdT1BLk8dYPG80kkrWMsSbQTk2dnN5rMlug6AEkdKyeNvt5PKaHRpo1xYXYqUMpug48sfv4gKICVVfJqK1+NlX88+2kJtnDP/HM6vPJ9nD/UwEE6yYrFe7Nq8PwXPcf9RUCSsvCm76XhfmDWr9etfWFtq/Dgn6ZQBUVPgZC7jaWI8d9T4CO0jR2DBWaTVNIMhheoFZdl051ammuSioTeIP5amyjzWtEpanr20kfUwn/5Q0hhbp2Uowu/f7+crG1dTWuJF0yRdo3Eu9b9GVSQCyfVQvgiBxKcmkHKe8zTzT8HfBGUrCtq1BFu4hYkpWCavd75OS6CFr1xlLPRpqm4591nqvWUSjtdqrjknugqWi4vL+HWyXFzOYIQQ1xTbBVw1gS42ALcLIW4DyoGFQohfSSnvmq45FpAVICRre7ejqgqadn1uf/12wskM+C7HJlBqGh5NQfMU/rw7pTFuGowQTyssKNUz9UlkQQxWTXKQQ/Fe/jYVIqHpgk4kabcoNMV303QCdrXU8x+3/R1NAX2lfSiSol97nY+mqwjE5nOyq5/FyR4kS2CM7IDWGCxVSiIJJVuPyimexYZ1e90z9Ax0MKqmHcey9RHu4+xYA4omcgkzjFmYc9mx/zhf9PyB14LL2XDb5zlrgW6V+9WBzlzzg1tIaQrajf/beX4GB/oPANgUrPFizNTWPdmzEJa5/7rx16hS5dJFqwFQVI2UotE1Emfx8lKHnopgKB4+NZ61VIFhwUrZE32kjAx08eAAYU1DW63xlpGA4pn3uviSlDwWOMaazte59cJbsy6EKSNGTsl/0LKZzK0LAAaJwt8kr/AipSSRMV0lTdNLoWtkINJLIhXmPN98o1/7szAUTjLc04x6jvUcx1kjiQxCuA3mLyWSTLO9xs/KZRV85prlhW0nEv+YpzT4E34OjvaQUFXLuRnKpEXBesXTRqUsYUNed70Bs8SB/v6dlhGu7HuOzNuVlP7Z/+JgxyjNg1Gu1pLGio4+ftnwUdb1Ps/A/L8FcopTzqqch03BHVvNsZaAkFJmE/eYfTcHmrP7pJQEmt4lOhSk8tZv5zoxXJPhs2OOdbpxY7BcXFxc/khZtWoVx44dG7/hmc8h4GHgP/L+PQwsHu9gKeU3pZTLpZSrgM8Bb86ocqUPaiSCkHjQBUbP8e025eG4U+r2hhf5SM84KcUtEtIbJwZ5v1sXXkURIbBXiQAQycTwegQ9ySO8quYUCvO4RcleSoPv8stjv6QzrO83E1yk1DiJjMqyRBtlaqTgWID6VC6WI2vjkipNAxEa+sJZYcxUlFqHIzzyehPt/hiPvN5ELO2cIjqQSBPLqIhQD/TaU7Br1tgiw9K1NN4M7z2Wm4vMvfCpSUZjaRakB2kZygnjwbhF4ZSSx0PHeeLYE9lNb/W85WBpmTwH20c51mfcL8uNVM3Mh5bEISCJpScZh2W4uUlkNvkEoFtL3/3PgqQkAEvjTZSFO9HQyCj69YwkFVqGovSHk7SMnATIutBlVOfsiFnB2+aa6JzsBGA4mqY/lOT/27cjq/APRZLEEvkxWfDsG//MjrcesCgx+lhmlseeYAIQ2fkDNAaaearhqQnEEgn2GZk8O0dyqerLlAjXd22ByAAikwQlVXAetr7zFKzfNP6Getllb5tVsOwLKFGRgdY3ebdjF/v79ttnJwSb6zbTHDnIa1ot7XE9rmsk6hzDVxLqpEH4icftLpZTccVLG9czmlR0pcnSR113kIzq3GdG04t8xxIJjvXmlOvu0XhB0py5iqtgubi4uLh80DkBfFlK+Wf5/4CJ5TifbbLLxTmB2TPSCkqKtJqmKZ2rbWMTN4absi9jeUkSxnJTahoMZwWXQsNCbsXZ5xFEYwcIyUIhdkF6gHLVuV5XznpWpO4VkvctCpZpwdJQsoqTNNqZ82kajOLRMhze/RL+1El6o3l1mPIYiHTTd+glRqKWZB0OspqaLwRb5jhV6v31Y9b0yY4lJT2B4tkZFTTeSOsxOwLnmB0ppU3hS2TylKxkyJZQxHa0UXjZ7MtEC3YWtK2P7qAt8bYxb33uVvfSbYkGfV9MT3jiNWKssjFTec+CCPdD3xFk7VMFM5DpmJ7NEVgab8OjKRzrCeGPplBkmrihSLYOxzjUNuxwYsa1UvWxNaP+2W+bfsurHa/iEQKJsM3/ZKCRSDrieN/b/FG6R/Xnra47SIe/0PI4L2O4c4Z6ofYp6D5o62swnKQvmLCsbIzh9pZtY7oIOjigdb1H3cnnOTJkr/Fl9jGQ1K1Dh9N6ohTVVg9PR5MavxzeR4sIcijm7JYupSSSjjAcL7zOndFh/v2NvYxEU4TiGX6yu4UDbSN0tDUy2n2C2sHa3HS7Cp8pE1Wqeh1Acs9AMJbm14fbafMXlkw4J1KPL1aYSfF04roIurh8wOkOFv+xNn3fp8Lm4NHcmyLxZ/dddd+U+3dxmUb+leILil+dTEdSyj3AnlObzviIrChmBiBJhtQEqVArrdFu2uJdrKbKnFTB8fNTQ4SOtzL/mk/pLlOWmBPTfc9KXzDJ0kqzO3t/5pHdkR7WpMtZmmhHjhE473w+znVgs5tkTpBKZFQShsAcyvTYTu83jb/h4ys+bhyicV7kKIvDBxhVJG95QqyhMIGKeXhXJkLV6EKO1vTw2Y8swiM8jsJzTGQKjhfA/PQgy5IDUCQXw3hMpFiyP5pGShjJdLCsZJVt33AkxSAxhoQu2AtDIE+MdJMJ9FCyZDlS00grGilL7aOCZAX12/VshNWXFoyvafmpVXRC8RTxkRjCGgqjpQhqvRiPjeHW5XSORvxWfgmIvOshELobojdlSUZh6eHIr6hMfYw1/l0Mzb+UtorSrBDu8wgwrCECDYlESocPveGeaLrQmqn5zcf5zYQ1ZjTvKtQ9qys4517FYFjvJ1qpkPKqyHIKyCrAHm/u0yxz16JtOEaFt5QPZ4dzUuCz6qf+x3R/LJIQJJVRUfIWVnLWYDNhhv5HtaWxF8Y2NftBTaRVMqpGideMjcu1fqrhKayohkvf082vcCIoKC/xUOWNs0iE6BqZx6JUL2heDg4cZN0568DfzMWdzxGafyMj8y8u+BzqWQjt53ZkqI73I6/y8XR+llBYFdjPssajcNn39I2RAag82/al83rn66TVNH9+0Z87XrvpZkIWLCHEn8z0RFxcXFxcXGYCKeV2KWVjkX07Zns+EyJrwZJZIeHFWAuvd+4iktZd7FQjCYZ0EGDXDu6grO89/c3ef4cG53TXYw1tYioGb/cdQbOkxO5Vxii+Gx3UY1TMPh361eds2Q9kNIk/lub1hkEiyhBRdSRrUUuFBoklksTSFmuBlGhIPFIhLVVCaorfd73OOy2DPPJ6E5rMiW6mcKtKyW8af8NzJ5/LxgNN5IJcNrST88LvAzAgQ0TSzqm+i1m6fIZQ3xfty7pQ5mMKwYHgDt29zBCC2/0xfnWg06ak+TR97nU7fshARwOMtvPz4+O4hwLDgRAH2kbY1TBIJq9WlGaJ/QM9eYcEajtH6A8ls5YiJxLDraSzMTX25CEFdO6nOt5sa5tzP8SiYNmP9RrJU8rUGCBIZKx1pMx+JP2hJO+1j2bTpGfPLx3nqBgmjF2J1lOTC5ozOcswaV2RDcbTdAYHaR16P2tFy55zJleQ2UoyrVKmGp+Pxj9kt/dEiif9GcsVUTM/5yOtRFMKcdMq2Vtra9czGuNoj34O80q9Rr9mH7m+ANRC/crWHmBPo26lOu4/TnPoOM5IajsDHOoIZK9EeYmX+cGTXN3/a8qz87B0HPMTlEkOxl4irAygSY10uB8S+tzTSobOEcvn/PgO2o0yFFHGKDYOEOym+8B/kujQY8sToV6ah4+RVJIOKf1njolasH4mhCgFtgLPSCmD47R3cflAMS0ZAtvfOkVr0odOfQ4uLi5nBMJw9aL2yWysTb6i4Oxsl8Mmrw2dQJzn9B0i8/5aVroN5U4A8YzKyf4IDb6cG99L0TbuyyQZTtqFTsDiqvihbJ+qZheXzSxnVkJxXXhKZlTSMidg+bQkgc7jJL0LERf+qd4204fmuQTIWXPeTvbTGa6kLnSSCs41zsV+/Uo8nqx4fXJUjw8iGcrWmHK6QkLaFYs62UlJ/0v8t8u/btsezPSwZ3gAHL5Q6CQcAAAgAElEQVTONakRSoXY0VJcpzcv/YLUAHARaBnwlBGIOwmV0i60hnpQFp1vV1QcaB2OIoHGwSgXeGOsKYFIIkMinCQS6+cN0cmVspp4WjISSbF4XglmxkdF1c1C5rjlSk4cPP7Wdg4uuZCV5R9hsS+X6CGazLC/dQRK9ePe7wpy6egezgsH9L6MdomMqiduqSicvQRDWdK3h2Wc0UzEtt/836MpdIZTDJEkmlSoKM2JuqONL9MhQni0Nmh5Q0/wdOGfIgRoQthHHTwOF1zPE++0UxfZzkc9ndy3+MoiV9VOc2sTi7Wc8l7czVRaXmkONk59f38olb3m9b0hfP2/4LrP/R9oetXWulSLUZkeYjCczFqBs3FsVqW16TVU7TL7UJkklFbaNoUS+idlb89e+kMJlvJp/RnNm2hBwhILXuH8TTWMLg+FlD4ODe6n+5hRy+vCP6VpMMJwJE32KRo6gbrwfMqVEFJA+0gM8qohmFNSEgFeirQSOvQ891at5803vkOvR+XsS2+m3OtgapwhJqRgSSk/JoS4BPgScFgIcRB4Qkr5+ozOzsXFJcvycM34jaYJm/vfRMjPQmgSCMGfffPUJzTXmO7iyhO4Rg8++CDPPPMMXq8Xj8fDY489xvr166d3HhY2btzIww8/zLp166bcx9NPP833v/99ACorK/npT3/Khz/84XGOcgFs2pFZBFSXYaRFwRIIBzG6nSAJobJazrf1k7V+SAnBLlh4nuPQmgQ63oH2fXDj1xFm8L9PMhCzx10kMnGOBvdQquasWaETb7KwvER3qyvzEkupSAkv1Oe+V7JWt2x8l74tYQTFe4RAkznrg6lAeWUGj+HuqEkNFa+hAGk2645XCAwDn8VFKnc9TQXLZ8Ti0F/8O09KmKcGkQ7uh2AvCNyWeId0qosV8+bZrI+gpxN/+sTTDv1LWkeGeLMvZ2XMphFXM9kU1YpM0SpyCs3SRAfd+/OKi6vOiTT8CT9d4S6uOfsa2zMR0ZLsinfTWtNCvMLDRUtaiYkMIZkirejCaEbVsu5u+QpAVbw1+zohFDSp0pc6yiLf+aZnK92jMQ60jXDNZYKRaIq+UJCuuD1eTko9ZXhIGGn886w5EkkslVNy/YTAcE4sURN6Qd4K3VLilRlepYOUR+MGs/aAgWL57ERb3iWUzLBIc4plMw+wLwLUJoe40nDTGyRGHOfrXaLmFlNlviIMxRM1SOPcPR5jnkZ7KdHQ2D7QxhIqqNBKoNM5RmpeJsyuE3brsW0IgN4a1IrL8KlJypQognKUQ7+k8fJ/tLXVHOYZTA9BWcHmgrEEMESc0rzyZPpOfaNXKnikQnvfuzaFpCTlL4gxVLQMi5PdvCMUjoSHYaFKUOmmL7iGi4F4WuUXL+/lmvR+UqpGIB3lVwc6SQQjRESaZZqC8M1eMvcJx2BJKZuFEN8CDgM/Aq4Wuj/Cv0gpX5ipCbq4uLj8sbN//3527txJbW0tZWVl+P1+0ulx3CTmABdeeCF79+5lyZIl/OEPf+Dee+/lvffeOy1zEXqRneullM5SyRym1OshjqoLO1aXQCDnXJej3qNbmaR0LnacivjxjHRQoiRZFltGRWaETFnO4iClpLdzH0u0DBWZOCpaLgQkTxY96q9HAGfFch6Y4aRCZZkvu3Kt9wmH/a9zdr6IbgnCSio5AVoIPYOgEx41AyndeqEJUxA1XNPQlcdy5RJkyTm2K2MKuSVegSk2ezzOkRLdSpTBof9EJZf0wUkM31y3mdTwJ0lpMcJKv31npB8ycSMzn0Bd2KXH0Hjs8TMnByI8UfMOC5fkLHbSjODQcq5sXcnDVAq7wO8bqre972lroDSeAgExMpSTS+YglQxXL7o4ex0WJXvpVFvolCPEvQNUsAbzWdLyxfJssWZzfoUWJhWpK8Ca3e3Suggwr28/56cG0Mpy11VKfTyB1dKR+1+3tkVY4q2EJfoxHkuEy9mxBsqP/oqRq26nUYyyRi7JWmFLho/B4deybSNp/XouTg1xpNdH2KMQC0QM4d8pSNBrqzd2IDmAGtSf9fc8/czDi9Nn0IqiWq5WOg5dB9jZayivMmFJ4a5B3TMQ7KLzyv9hJNHQjxwKp2jyd1OTHKZazOOj8nzof99xPFX4iFuU0fy4OvOvqmlc6s9dm3Z/jGeOHKDSWORYmmhnec/LxD/8Rf1SKCmu79rCb6vLWLX8XBJplQXl1qQoZmIdc5ugQfi5Or+sQnwU2vYCUJEZYXngbTyZMNK8pWqGBU2/ZklCtyinUOjIhJFaNR6pogiJXyTQ4gfxBF7Er9zFxUBvOM4V4ZcJo8ICEBIWJnuzn3XNsZDzzDHRGKwrhRCPoGdiugn4tJTycuP1IzM4PxcXF5c/evr7+6mqqqKsTF82rKqq4rzzdOvDAw88wHXXXcef/MmfcO+992Z/xDZu3MjXvvY1brzxRi6//HIOHTrEZz7zGS655BK+9a1vAXqh4csuu4wvfvGLXHnllfzVX/0V8XhhJrbXXnuNj370o1xzzTX89V//NdGobq34xje+wYc+9CGuvPJKvv71rxccd8MNN7BkiS4RXX/99ae16LCUUkNPzX5GYBX1FDQ6CFksMbqg4CdBbfi3DKZPZttarSmK1Hj0wK9pDASNPgWBWJoTfSEGIimQGiuD77E00WEbWwIvRZrZHm0FTUWxWIdknoZVM3zEVo+p8DzM4Pniws0AMWqTwwxYXAY1CRKLwmURYD3Nr0NfHZrUa0FJwGMVnmLDXDKkC46D4aQlwYC+22tRqmTH29BTo7ugWbLtvZXoRaS7CywUQZL8KnVCn5NxkzJagobY7+lO1hhzNVCSEO6HQCcEOlCaXtWLsuahZ/zTLTTZuWYtWEq2T1UWpjf35CWOsGYMfMfTSxqLm1hfLeE9j2TT0a8e3ctoOMhQNEVaJmlPvIvWp5+DBEoSVpdJI912SqE7mCClFCq/aVTOjRxleagGsNRTMzNCtr+FlFCRCeDkAmh7LXMKs2bE2aUsKdS9DuLrs62v0ihGUSy9+UYabRk3Gwb0lN8eBO959PMbHexlnv8YEl0ZSasa4aR53yXtCT3tuaJJgokMmbg98Wi+WqZYgptUNEa0RG5GvTXQupvBoP6sZwxl1B9NcaQroFuWgRdqe/OySUoC8WTBtbJSlzSsy0LYEsqoeYqF+R2hapJSSxxlJJ2hO1lryyJaGe1ka/3PASiN6HNbFDpBT0cLx3rDKJpmM00NRSzKdV7ZB7NV9PBzjEaT2e+GEjVhTwAjVXYm26nI6LXoDoh+Xol1IDNRPBarti+jZw1cEDnBCAl+72nDTyK7OFCRCfKhod/nvoOkOqFEM9PFRC1YPwZ+jm6tyt5xKWWfYdVycXE5RcaKzzpT2Rw8WjTDILhZBifKLbfcwgMPPMCaNWv45Cc/yZ133snHP/5xAO6//36+853vAPCFL3yBnTt38ulPfxqA0tJS9u3bxw9/+EPuuOMOampqWLp0KatXr+ZrX/saAI2NjTz++ONs2LCBL33pS2zevNmmLPn9fjZt2sSuXbuYP38+3//+9/nBD37A/fffz4svvsjJkycRQhAMjh2a+/jjj/OpT31qJi7PZHhNCPFZ4AU5fnGb00fe1I5JP0c9IyylBMPhSN/u8VMVTzIsPMAykkqSHb2d2aXTfi1ObV8jMtTNpUsWI4RgIJy0KSvmKyFVru/aQufi9QyEzgIhiGpp0FQymoIQsDA1wFn+A3lztfZSiBm7lK3T5CDgHPT0szhTUmBtsr5flmgzjtYQ8aCeTEHYLSkBmaRK6osQpoBpFXbNlj6vvrNlMIwyMoxvnsQfTZFUCtv2CGucj6RXFCb2qI/+zn7OBS10MlKDmN88QV2J8Hgp8ZoCoMSXd/S7x5pJLNJYOt+5WHB+Msf8zHwaeja9eFqhQknx2+F2rjCzT5K7fn1KB0kW0Y1+flLY7+mADLGAChJp/bqHE4WucQrW62d1TTXiiLrbCaQzltK1ZlsYVhN0CSPFv+kmR/EnK1/BklJyrDdgGLhyo4vRFnbFc7WkzBT8XgRDxBDAopS+8BNnKSOWWLcF5bm7UaaEGUokUSUkh9rsIXZmnJ+UDI/8nO/87rLszI+KYY5EM3ys4rzs+QxFklBiP7OEopE4+gLywirbdutntSAbZB6/HtA/IxUZPS2+P93GYLqBj+SV7BtKJNkX60XxSeJCUit6+bRcRdzBvdRqKy9JjiDRCzWXhDtgwTKIByEdw6mcoOEhmlXYJJBSVA63j1CmRJEWN0N7vOhJrERJs4BySvoOO86uXmun1JjjkIizUi5ECF2Jhpwlaba/8idaB+s29OQWCdBdLYQQehiilE+NeaSLi4uLyylRWVlJTU0NW7Zsobq6mjvvvJOtW7cCsHv3btavX8/atWt58803OX48l+np9ttvB2Dt2rVcccUVnHvuuZSVlXHRRRfR3a0nJlixYgUbNmwA4K677uLtt9+2jX3gwAEaGhrYsGEDV111Fdu2baOzs5OFCxdSXl7O3XffzQsvvEBFRUXR+e/evZvHH388G491GvlfwG+BtBAiLISICCGcCzedTjQlK5wAWStEIq0SjNldQ8tUu8B/3JNbXX8jbS9Smi9gyFiurRnjtDL4Hm+eHCKpGG0bdjBPKyGtaMxTAgRUp8yBxQWX6kAtIAuK7GYFYMt7ay96jEpui1dLZ+c5EkszaKyUe6SexU0iSKRVhhTdAnuSUbxaCmFJXGCev8/QSoLhsJ5UAWzKlZVmkcuMZrbIn1sBeUpPOKlbfepTuuBLIgD+5mwMjTnH/mCSnqEASxId+DTdWjFw/En+q9WMcS0cMz/1uSfPrVICPYE473eHyKhaNgNbnAwjJPJaFnsHJ8Qgg+Ss257u/ZQr9o9Ov8i5OF469EpBbxGLdTUhFQ4J3aVyJJriSDKvppJxHvlubQtSZmyR/byfCzU7zj2T59iZNCySPV6lSB2xvFmn9Oe9Ot5sZoG3jSzQY+GQEq9MI5MtaKE/MN9QckYNB7WM1O0qyYxK63CMsJpzJ1UtLnzvtY/yfk8ITerzM597AJEKZcccDyEEXclDpLSYrphJe191KT+VkRaaxSjDIkF9YtSx6G9B7FjekyGCHfkj248nVxtNSknTYNTxozMcTpExF0NSEds+gSCZUbOFkXNzyt0QUxVW0AwLVs5WZbVgzSYTVbB2Yc+JU2FsK4oQ4pdCiCEhxDHLtqVCiNeFEM3G3yXGdiGE+JEQokUIcVQIcc1kT8TFxcXlg4zX62Xjxo1873vf48c//jHPP/88yWSS++67j+3bt1NfX88999xDMplzszJdCj0eT/a1+V5RDNejPAGtoFaNlNx8883U1dVRV1dHQ0MDjz/+OD6fj4MHD/LZz36WHTt2cOuttzrO++jRo9x999289NJLLFu2bFquxVSRUi6QUnqklCVSyoXGe+fMBacTtdAVDPRV7uahCFqyUCcsujprbD7QPsqv3m0jGM/kWTmchaqsvhHzEx5RiWdUFsky6o26zJboKsdEG9b+F6QGCaSGHUerSPRn21n3KkWUGIEkZbjBlahxpOkQJPQU7yYBklzd84T9TCUcE8NUhuphVF/tzxcYTUIO6dtNu4iiSkq0BL74ILWNHQXtNGm3NphKXDBjCLjJsJ7GHiARNO6d3n5JojMrnAMc8gwwkDqBEIKEWmglzrdglal24VRDZt3drIkVdosu3vH0Fr1zEsmopVgxUpJGzV7Pno63qYo3W3fbWJjqy77Ofz4kkqOaP6uQZTRJwDKWhgRDIY9rmWwCFIC4OkQKlQWpwVwhX6AtEcmL4dJR8tzj0lmlqlBNsc4ZIBBLM9ReT8FzmPe2TI3gU+OODczny4NgNJYmEM3QR5S+ZC52ri+kK7rHxDAJFBIZlaF0k01ZLlOjzG8xk6CMr2JVpIazGR5VqSL9rZwTyY0ZTSnUD25m1IgvS2U0ToiRsbuX8IqnzWkzIVLkx+WZ3WharlC4LzHCPu1YVtEFaBchEkJhKJoqLIptmW8+ImvlFPgMdSaFkrWkmhYsa6H02WSiCla5lDK7bGW8Lr5cqbMVyP/F/QbwhpTyEuAN4z3Ap4BLjH/3Aj+d4LxcXFzmIN3BRPbf/tYR2z+XydPY2Ehzc06YqaurY+XKlVllqqqqimg0yvbt24t1UZSuri7279djDJ599lk+9rGP2fZff/31vPPOO7S0tAAQj8dpamoiGo0SCoW47bbbePTRR6mrq3Ps+zOf+QxPPfUUa9asmfTcphtjMe8uIcS3jfcrhBAfOd3zKsBXhjj/Wqujk2237DpQeAyFAkQyoyKQDIk4/USN2Bg7YcOiIYRVEJVZH7uYliEpckJPgZCSihRK17Y56dan9qjd7addhNCQ+EYPmCPaxLODgWdJaTE8UsnLoCoJxXVh6+zYCSpS/dkp5/OKpw3V4lIpgTYRYqRvO4R0wdKM1/COIbOmVWfXt8xQM9XxpqwVzUpvKJldkTf3BPNSrcfSCj0nD+If3E9EGSSfTnKK9EC8G0VOPrGNhnS8P6rhAtguQsYW053KOaV2vpK0L2ZXRkLJnII0jn3PNo5JJGU/Hk0loSk0GrXGzP5qlHr2i15SKCxLtBWxQuWstUOxFKqUhBKZPDf8whnmW//iGdVWsNlEABnyLIWWQt7CotSZz5dZ6gBgWOSSV1jpFhFe83QQJElfqt62ELI42V1wSCKj5qw+eVT1PJXN8Hho+E1SEbuFMJLS1ZC05XPfLezKOUAok0FNJ/QU7gYZVbPFaWVUjb2ebsunzE5FVPeWkBKWxFpJCIXXPB20SnsdOU2C37DQOz0/+fFTo9pQdk+f4bobI8Nbnp4CK6PepyxYQJxJJhqDFRNCXCOlrAUQQlwLjBkwIqXcJ4RYlbf5DmCj8XobsAf4v43tTxo+8QeEEIuFEOdKKfNS8ri4uLjMAWY59Xw0GuWrX/0qwWAQn8/HxRdfzJYtW1i8eDH33HMPa9euZdWqVVx33XWT7vvyyy9n27ZtfPnLX+aSSy7hK1/5im1/dXU1W7du5fOf/zyplL6qv2nTJhYsWMAdd9xBMplESskjjxTmO3rggQcYGRnhvvv0WDufz8fhw/l+9LPKZnRPr5uAfwOiwE+AyV+4mcTjhUUroHt/gUJzTvQ4cmmh0U1i1CiykFQ0PEZihBYRYHlmAZTZj20QukXKmxnBXHOVFgFxW/gE3SJhjGHp35RThk4gygrra9lq5Rgp0xekBijR9L66RJhKWUJNeFR3j8k3Emgqo5lOlkfsmdKElAQtmfTeVd7j4yyjQhm1tzP+vpnpoiMrOBoughZh2BSAPUKgjrPCnUGlTYQozVubPi9SR6ykisC8lTaBv3UoymVOS9EJfa6hRAYJ1LQ8QazcLF+Qm8P7nlySiZNDA1Mq1dErIpQb8TFjnZ1XS7M8XINnoT5hazr4tJJTFNNifDerpFDolKHs+8r0EKVK1Cb0RpKKbd0grqi5TJVSQriXSDpNIJ6msrwEa+bMsEhz3LC2nBPJOkllu0tZlJ+TwRD+khRp47NhxpaNlZgln3zlUiA5IfR7mHNWc+5P5qkdqpAoUuP8SB0j81aTKCmMXdrn6QHOLthuVdYA6rqDugXTQWd4Cz2ubFGyh4YBhcuLpO/PpyIzQkk65zpcIwZQOnuoLPPlziH/s+pw6lJCLK1Sjo/l3Tv1bZgeEvoBUZwt9YF42lbewMS0UmWvuUUhHhS6BTEscnX0YkLR14rGON+ZZKIK1v8EfiuEMJcszgXunMJ4Z5tKk5SyXwhxlrH9fMBaqbDH2FagYAkh7kW3cnHBBc4paF1cZoxi9Z5MjAKYLpNgvJpSH8Q6WpPk2muv5d13nbOLb9q0iU2bNhVs37NnT/b1xo0b2bhxY8G+jo4OPB4PP/vZz8Y8/qabbuLQoUMFbQ4ePDjmvH/xi1/wi1/8Ysw2s8x6KeU1QogjAFLKgBDCOXvAaeatFj/VkSQHhkd1QcGaZAuHlOEWIWee9JEQukB1TlSPyQuIFMud3AGNv0vjbTSylNUsplQNE5e6hWVeSS6leKYweTdAESuCGeAuqUwPEU+uYFGq17ZfEVp2lT7f6lGiJZAOTjYCjQ5b2JwsmFUkpWTdBQelpR6R0SyWUJHzcnMYiqRs7oXFGDUVOweJskIJEGAlixO5uLf56SGoyEVXZDQ9HsWj6kKgqdSpUjI/7SftrSiIqTMZCjtbr5yER+u2JhFgVCaBpfo5UFjcGaDUcHHLtyyZ8wY4IfyoQhvfPAWc9NgV3qWJDqxqrsewR5hdWe+hRMKJnRzrixD1qJT6PKhIhii0/HgtmRXNme/xdGefnAyaYxFcYRT1nYrTmEBk52u6tMWVkYJsfdbzsipgZqKUhal+EiWLCo7BOCIf8z4IYIQEoyRZIsvH1CAWpAfxaWkCSqHLqxXzei1NdCAtqdeHRYLFlDgflIfV3TaUyLDA+JyMVYTYiWi6UIlPKhqVDm0B5mf8BduKuRrOuSyCUspDQojLgEvRb+VJKR3yhU4dpzN2vCNSyi3AFoB169bNrkOli4tB0Yx/R15jxeJ5zvtcXFxONxkhhFm4BiFENc7ljU474aTC/JSzkKC7fdm3DYViLPLrFg5TuSo4TgkAOYEu/we00TOKKjVWBA/zlpKiVM2AZZU5IRQ82YPGFlSC8XS2aDDA/GhnoduZFDblTLG4KwkpOSt+YswxAHxaypbIoCeYsI1jjd0w3akE+gq3eQpObmBO5ITlQoRU8UjVFj/l1DYQz6B5oNxjd7ucnxlBUNw65FxSGlvqcpNknpXJLxIsD9cQL/GSGUdBSsrilo600PAW3Ts2FYFdtIlcnJJATyVu6qoqVrc6Xbk3lZLReAYEHBB9CLmk6BgZrVCJUPM+3h1Zt8tiNqdCRL7iJAoV0e7kfqpjTQWfClVIW7FrKyVanMr0kMMeva6cv2K1bVtXMA4eva93PL2OxzkxTwmM2yY5xrMH+j1xLLWgpnKWLdu1zimWfanYxAOSHDC/F3KWK/3V/HShYmUf3SR3V+aiiyDoLhSrjGOuFkIgpXxykuMNmq5/QohzAfPJ6gFb5s7lQF/B0S4uLi4u08aqVas4duzY+A0/OPwIeBE4SwjxIPBXwLdP75Sc0QVqmX1tpWs0RmmZVWKRBHubyWTGXqXWgGXxVuIlS83DdIHbQrGkD9Y+QI9LSqsapV4PJelQQbtEnuCfTKUpy2szIhLMUwz3wzyr0OJkNxPlsGcg+3qs2ccsK+P+WBoj4/SESVBc+Uh5F7AgNVCwPZywr0UrmkZSUWyKTFLRwDu19eL6kaAtPfpYxIus6usYAvE403BINDchDnjsIl2LCNhk7qDIPbsSSbsSolHYrWDWeeZT7vPYlHXzPJS89qaL42TE7PPy3FQ9FOoL89P+MZO9tAwXxjcBVKSdzhFK1Rhezf7sREXx+KRTJd/imM9QOFlg5a2ONcPCnO1ZQaNETTIvEyJaWg1AXCjsERP/LI+F+X1oPoOmu7Fj2zlgfpmQgiWEeApYDdRBVs2VwGQVrN8BXwT+X+PvS5bt9wshngPWAyE3/srlTOVMqmc11lynyxJXEDfQbuk3UCiYubjMFFLKp4UQNcAn0GWs/y6lHN9MMsfQpD3l9fz0CPMzw2McoSORzFOCOQXLQVRrFyEWkpqQEDcYSbFi8TxK06Ei0RQ5nKwzflH8+2cs4WmqmOfULPQV/QXpQoVoLEzXLqdro3jKHPeMF9dVMDkHlodrSKk4Omrt9uguiQtladFU82ca7yUHUTTJgCh0Z5ynFK+3Z1VIcjFP9oQMFLSYGp48FSs/QYYVVUr93kzSeHJ2zP7VlFv8mF3tIZVRHV1oy9QwaOU2ZXZJUn8eJ/L5nawtaeI2R3tbe8KLuWfBWgd8aDKFGYUQz6IntKgSQvQA30VXrH4jhPgHoAv4a6P5K+i1tlqAOPD3Ex3HxcXFxcVlIgghnpJSfgE46bDtjCH/h3hJsnNKxznVvQEIy8kqN+OLBqVqbMz9sykymhn0KjLFV+0nOx9rFjkrsbx4krGlqOLC38JE65iZxXx4GI46x1dNBNP6Mhv3YbwxpjIHKQtT1ENxi2yB298kEGqacum13S7PGP0pqjaGcjBxgd90wZvI9ZnO+5hvjS6G1UVwMklEJkpAjG2ht3GmWLCAY8A5OCSdKIaU8vNFdn3Coa0E/mmifbu4uMw8s2HdcnGZZa6wvjHisa49TXMZB0GMDL/ztBTsmcxKbl6XAJQYRWz9sZRtl8y+nnj/oUSGqOIrcP+bNHNAILLSM8b3n9Nas5AaFFGyrGQ0CcLpZMe+AKVqgsQYsviog7Vnasz8jRj3+ZVMycerXCn0iLBmQ5wuPELgLVCMiisho/EMxSLHSrWxFx6smMpiflyZE5NNLDFVIkkla02fyLyszOQcrTMx71QwnqacSShpp8hEFawqoEEIcRBys5NS3j4js3JxcZnTnElukC4uQohvAv8CzBNChMn95qYxkibNJUwBPulb4Lwq7yigj49Zu8cpO59VIPFpiaLZ7PIJpxTKcI4vmQxzTL8aE6e5OmUyK4aTO5lnnCQDUSZfA2symNa82YhdySUrcN4/FEmybP7kkntONFHJdJE/9fE+L1NeFLFgxtpNypIzJSY+V6ursjXerUwtLIae3/dEMneeKtbslP2hJLHYxL7XpoOJKlj/OpOTcHFxcTmT2Fy3eVr7u++q+8Zt8+CDD/LMM8/g9XrxeDw89thjrF+/flrnYWXjxo08/PDDrFu3bsp9vPTSS3z729/G4/Hg8/l49NFHCwoZzwZSyoeAh4QQD0kpJ533XwixAj3m+Bx0XWSLlPKH0zzNLJo0UjFXXFQQYH8qtGWLyubSPTvRk2ma9eIxYydgmGPMgFzoc8iAZ2WqSvVkyXdpnAnGUzYymmQgMnuWhslyMDXIgGdyFjJApNwAACAASURBVMPUOAr0RDBrPM1VmhyTkpxedolOFqEr65WhJjzLzpm1sSeapn2vEGIlcImUcpcQogKmnK3TxcXFBYDNwaPjthlf9fjgs3//fnbu3EltbS1lZWX4/X7S6bn9YwvwiU98gttvvx0hBEePHuVv/uZvOHny5PgHzhz/RwhxF3ChlPLfDM
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Comment doesn't match code (and again, do we need 2000?)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Never say something's obvious! (or trivial etc.)
If readers get it you're wasting space, if they don't you're making them feel bad

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Again please have a look at the number of iterations in all these examples

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

See Ben's comments on previous notebook

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should be :math:`stuff` , same goes for rest of code!

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please remove these lines, as the debug variable isn't used anywhere.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please remove

Copy link
Member

@MichaelClerx MichaelClerx left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @lorcandelaney , have made some initial suggestions to go with @ben18785 's comments.
Overall, it's looking very good!

One question: it seems the two notebooks are very similar, and mostly showcase the same code (just with the adaptive option switched on and off). Would it be better to replace this with a single notebook that shows how the adaptive method is better/worse in selected cases?

@MichaelClerx MichaelClerx mentioned this pull request Sep 3, 2019
49 tasks
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This means if (type(w) == int) or float)
and if float == True...

use np.isscalar(w) instead

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

if np.any(w < 0)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants