Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

MultiNest #1319

Open
wants to merge 61 commits into
base: main
Choose a base branch
from
Open

MultiNest #1319

wants to merge 61 commits into from

Conversation

ben18785
Copy link
Collaborator

@ben18785 ben18785 commented Mar 13, 2021

Created MultiNest sampler and simplified ellipsoidal nested sampling by creating an Ellipsoid class shared by both samplers. Things done:

  • Created Ellipsoid class used by both methods
  • Added MultiNest which uses a EllipsoidTree class to form a binary tree of ellipsoid leaves which cover the sampling domain
  • Corrected some parts of ellipsoidal sampling which meant that ellipsoidal updating wasn't being triggered
  • Changed name of sampler -> controller in rejection sampling notebook to avoid ambiguity
  • Made a notebook that plots lots of pretty pictures of the ellipsoids produced within MultiNest and shows how the approach works really quite well for the tricky Goodwin Oscillator model

- created new getter
- used getter in notebooks and test files
- also added points to ellipsoid
- looks like way too many ellipsoids are being created though
- also often get a "singular matrix" error (likely because ellipsoids are too small)
- also added f_s function
- changed n_points > 50
- added getter for ellipsoidtree
- corrected pseudocode to be more representative of actual approach
- added enlargement factor capacity
- made ellipsoid gap used
- corrected hyperparameters
- also added max number of tries for while loop for kmeans

addresses #282
- Error handling for singular matrix issue when forming bounding ellipse
- Added back in recursion as it was used in the Matlab implementation we were following
- Cleaned up the notebook and removed the logistic example as Goodwin is better
- Improved docstrings for algo pseudocode

Addresses #282
@ben18785 ben18785 marked this pull request as ready for review March 13, 2021 13:54
@codecov
Copy link

codecov bot commented Mar 13, 2021

Codecov Report

Merging #1319 (9d949dc) into master (9052e5d) will not change coverage.
The diff coverage is 100.00%.

Impacted file tree graph

@@            Coverage Diff             @@
##            master     #1319    +/-   ##
==========================================
  Coverage   100.00%   100.00%            
==========================================
  Files           87        88     +1     
  Lines         8933      9190   +257     
==========================================
+ Hits          8933      9190   +257     
Impacted Files Coverage Δ
pints/_log_priors.py 100.00% <ø> (ø)
pints/_nested/_rejection.py 100.00% <ø> (ø)
pints/__init__.py 100.00% <100.00%> (ø)
pints/_nested/__init__.py 100.00% <100.00%> (ø)
pints/_nested/_ellipsoid.py 100.00% <100.00%> (ø)
pints/_nested/_multinest.py 100.00% <100.00%> (ø)
pints/_mcmc/_mala.py 100.00% <0.00%> (ø)
pints/_mcmc/_nuts.py 100.00% <0.00%> (ø)
pints/_mcmc/_dream.py 100.00% <0.00%> (ø)
pints/_mcmc/__init__.py 100.00% <0.00%> (ø)
... and 12 more

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 9052e5d...9d949dc. Read the comment docs.

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 Ben!
Some minor comments on the first few files. Will look at the rest at some later point in time :-)

docs/source/nested_samplers/index.rst Show resolved Hide resolved
docs/source/nested_samplers/nested_multinest_sampler.rst Outdated Show resolved Hide resolved
],
"source": [
"n = 400\n",
"log_pdf = pints.toy.AnnulusLogPDF()\n",
Copy link
Member

Choose a reason for hiding this comment

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

This is very confusing. You sample from an annulus distribution, and then want to convert it to be within a certain range.
So why not set the parameters for the Annulus so that they end up in that range? Or if they are fixed just scale everything e.g. draws = (draws - lower) / (upper - lower) ?
I don't understand where the gaussian distribution comes in

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Thanks. Yep, I think I just need to explain this better to be honest. The Gaussian distribution comes in as a prior but I hadn't explained this before.

"metadata": {},
"outputs": [],
"source": [
"from pints._nested.__init__ import Ellipsoid\n",
Copy link
Member

Choose a reason for hiding this comment

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

If these are useful for the user, we should make them public.
Users should never import anything from a hidden _module.

(Also, if you were going to do it, the first line should just be from pints._nested import .... The __init__ file is what determines what lives inside the pints._nested namespace)

Copy link
Member

Choose a reason for hiding this comment

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

Having said all that, I'm not sure where these should live :D

  • pints.EllipsoidTree seems to general?
  • pints.MultiNestSampler.EllipsoidTree is possible, but maybe not nice from a code organisation point of view?
  • Some new module pints.multinest where these two can live?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Yeah, I wasn't sure this. The notebook using them for pedagogical purposes (although as your comment above suggests, these bits need to be improved), but, for running MultiNest, they're not necessary. That's why I kept them non-public... I'm happy for them to be public but I just wasn't sure about their utility for a user?

Copy link
Member

Choose a reason for hiding this comment

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

Wait, the multinest class doesn't use them?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

MultiNest uses them but a user wouldn't need ever to call EllipsoidTree or Ellipsoid (or any of their methods) themselves.

Copy link
Member

Choose a reason for hiding this comment

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

I dunno. If I was using multinest I think Gary would make me visualise the ellipses to check what it was doing :D

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Hmm, re: location, I'm not sure either. Ellipsoid is required by ellipsoidal nested sampling and multinest; ellipsoidtree is only required by multinest...

Copy link
Member

Choose a reason for hiding this comment

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

We could change from pints._nested to pints.nested, but only import the samplers?
So then you'd do

import pints
pints.MultiNest

but

import pints.nested
e = pints.nested.Ellipsoid

?

Copy link
Member

Choose a reason for hiding this comment

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

You'd still be able to do pints.nested.MultiNest, so that would be slightly messy

},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAACVvklEQVR4nOydd3hUVdrAf2f6pPeekNCLdKSIomABxYJYse9a114Whf3su6vYy+rq2lbXtWBFFBVWUVEUaQm9B0glCenJ9Jnz/TGZMJlMTSbU+T1PniQz55577sy973nPe94ipJREiBAhQoQjH8WhHkCECBEiRAgPEYEeIUKECEcJEYEeIUKECEcJEYEeIUKECEcJEYEeIUKECEcJqkN14pSUFJmfn3+oTh8hQoQIRyRr1qzZL6VM9fbeIRPo+fn5rF69+lCdPkKECBGOSIQQe329FzG5RIgQIcJRQkSgR4gQIcJRQkSgR4gQIcJRQkSgR4gQIcJRQkSgR4gQIcJRwiHzcokQwRcLCst5avE2KhqMZCXomT11ADNGZh/qYUWIcNgTEegRDisWFJYz97MNGK12AMobjMz9bANARKhHiBCAoAS6EGIa8AKgBN6QUs7zeL8X8BaQCtQBV0gpy8I81ghHGlJC8z6wmcBYB621YKgFqwHsFufrNrOznVINSg17ftjNeQ4FDYoYGoihTsZSZ43jqW+3RgR6hAgBCCjQhRBK4GXgdKAMWCWEWCil3OzW7GngP1LKd4QQU4DHgSt7YsARDkNsZqjaBNVboG4X1O5y/t63IeSu7gRQd37daNLAy70hoRck9YaM4yBjGKQOBJUmpHOEYtKxOqzsN+ynylBFnamOJksTTeYm529LE82WZsx2Mxa7BYvdgtluZm31WgCGpQxDCIFCKBA4f+tUOqJUUUSro4lSR7X/nahLJFmXTJIuiWS987dOpQvtw4twzBOMhj4W2CmlLAYQQnwInAe4C/TBwN1tf/8ALAjjGCMcTkgJ9Xtgz89QuhIq1zkFucPqfF8oIbEXJPWBXnFQsw1OvAuS+0JUMkQlgSYaVFpQap2/Ec7j7RamPbuUxsZGEkUzCaKFRFpIEY0M1NUzK9kB9Xud57YanOdTqCFzOOSfCAUnQd4EZ/8+6GzSMTD3i1/Y3RJLTnorexv3UtJcwr7WfVQbqqkz1SHpXARGIIjRxBCniUOr1KJValEr1WiV2vY2sZpYHNKBAwdSSmwOG/Wmespt5RisBuePzYBd2r2ONV4bT05MDtkx2WTHZpMTk0N+XD79EvuRqEvs0tcX4egmGIGeDZS6/V8GjPNosw6YidMscz4QK4RIllLWhmWUEQ4tpibYsQR2LYXdy6Cx7XbQJ0LmCDjhVufv9OOcwlzpRcUOhMIp3G+aNoa5n22g0pqMS47q1UoeP3souLRohx3qip2Tyb71ULICfnsJlj8PSg30mQKDz4MBZzrH2IbVbmXe999jjSpGq6tAqatAoa1EKC28VQwUg0ahIS8uj4zoDAYnDyY9Kp20qDTSotJI0icRp4kjThNHrCYWhei+k5iUEpPdRL2pnlpjLXWmOmpNtdQaa6lsraS8pZxt9dtYWroUm8PWflyKPoX+if3pl9CPISlDGJE6gozoDIQQ3R5ThCOXcG2K/hl4SQhxDbAMKAc6qR1CiBuAGwDy8vLCdOoIPYKxATZ/AVu+hOIfnRq0PsmpCU+8AwomQUp/CLMAcZk+/JpEFEpI6ceC0iieWpNERcM4escLHh3VwkRZBFsWwvZvMSrVrO93MqvT+7LGXMP6/esxp5rRA9KuwW7Owto4GoclDWlO4ee7LyE9Oj0sgjpYhBDoVXr0MXqyYrJ8trM77NQYayhuLGZH/Q62129nR/0OPtz2IebNZgDSotIYkTqCUemjOCn7JPLiIs/YsYYIVFNUCDEBeFhKObXt/7kAUsrHfbSPAbZKKXP89TtmzBgZSc51mOFwOM0Zhe86BbnNBIn5MPBsGHQO5BzvFKaHAZ6mE3Bq8rOnJ6ON28ZPu75iVd1mbEgUUjLQoWR01jg+2dKPmro0pDUJ9zCM7AQ9y+dMOQRX0j1sDhs76ndQVFNEYXUh66rXUdFaAUBebB4nZp/IybknMzZjLCpFxKntaEAIsUZKOcbre0EIdBWwHTgVp+a9CrhMSrnJrU0KUCeldAgh/g7YpZQP+us3ItAPI2xmWD8ffv0H7N8OungYehGMuByyRoZdCw8HE+ctpbzBCIBQNaCOL0QVV4RSVwVA7/jeTMqZxLiUEYzYv4eYFf+C/dtojOvP3fUX8b11SHtferWSx2cOPWq8aEqbSvm5/Gd+Kf+FVftWYbKbSNYlM61gGtMLpnNcynER08wRTLcEelsHZwHP43RbfEtK+XchxKPAainlQiHEhTg9WyROk8stUkqzvz4jAr3nsTskNc1mKhqN1DSbMVntWGwONCoF8Xo16VGCviXzUf/2D2jZBxlDYcJtMPhcUOuBwzfIp2DOQpSxm1AnrkAVXQyAzdALe9NQvv/TzeTG5XY8wGGHjZ/Bj49BXTFLFCcx13A5uoT0Hrmmw+VzM9lM/FL+C1/v/pqfSn/C4rAwOHkwVwy6gqn5U9EoQ/MQinDo6bZA7wkiAj28SCnZW2vg1121bChvYFNFE1v3NWOxOby15hzFb8xWzSdPUUOhahhb+17PkInnMCz3wCaiL7PGodRmDVYDH237iGdXvolUNuCwJGFtHI21cQTSmhzYdGI1wS/PwS/PQlQKXPQ25Hnu8XePnv7cujpZNFma+Hb3t7y35T2KG4tJ06dx4/AbOb/f+agVXdjIjnBIiAj0oxSLzcHynftZvGkfP+/Y326CiNerGZIVx5CsOHolR5OVoCMtVodOrSSqtYyE7/9MVNnPNMYP5Lvsm/m0cQArd9dhc0hG5SVw77SBjO+d3MGs4c6hsDdbHVY+3f4pr657lVpTLQXRw9m5YySGxv64bOEhCc3KdfDRVdBYBuf/C4ZeGLax9uTnFo7JQkrJbxW/8er6VymsLqRXXC/mjp3LxOyJ3RpbhINDRKAfRUgp+X13HZ+vLefbTftoNFqJ1aqY0CeZk/qlMLFvCgUp0Z1tpFLCmrdh8V+cvuKnPwKj/wAKpzBsMln5bE0Zr/+8m/IGI5eMyWX+6tLOAwAEsHve9J69UDc21Gzgod8eYkf9Dkanj+bOUXcyIm1Et80ai1ZuIeubPzLcsYW/au5g+Fk3hEWDLpizyIvneng+t3BOFtJhZ1nJUp4tfJHipj3M7DeT+46/jyh1VLfGGKFn8SfQI9veRwiNBiufri3jvd/3squmlRititMHpzN9aCYn9U9Bq/LjfWIxwKJ7YN370HsynPcSxHd0QorTqblmYgGXjs3j+e928OpPu9AoFVjsnU02WQn6cF9eBw4IagPJ2b9gifuG1KgUXpj8ApNzJ7dPVjNGZndZAC8oLGful3uR1tm8pX6KOZaXuOazJODCgH0GmkiyEvRehW44PrcKL/36fN1mgfLVztVI1San737zPmipBmsrQjo4GZiEwKLSUF/yDyp/fpPsnAnoUvo7PZwyhjmjcv0Ea0U4fIgI9MOc0joDry0r5uM1pZisDkbmJfD0RcOZPjQTvSYIF0JjPbx3MZStglPmwqR727Vyb+jUSuacOZCh2fHc9sFaFAIcbuqmWiloNdsomLPIqzDrjta8oLCcR77cRL3BCtjQZc/HHLcBR/MIrh/yAFPy+gfVTzA8tXhbm9lCy5+sd/KF5gGeUrzI1d/28zveQMnDFhSW02q2dTpOr1Yye+oAv/0G87n5mizi9W02cJsZtn0N6z+C4p/A2gqAWZPIFmsmpbZUjJqBDO+Xy4CsRFAoETYLWnMzon4H1WW/oS5ZRu7unxA2l1+DgJR+zijc3idDwcks2G4+LDZ9I3QkItAPU3ZUNfPKj7v4Yl0FCgHnj8zmqgn5HJcdH3wnhjp451zYvw0ufscZPRkk04dl0mQaytzPNhCnU9FsspEQpabFZKPB6Azz9ybMupopseOxDnTZH6CO24Sp6kysdZP4x3dlXDomfALdXSg2EsOfrTfyifZRzm/9kInzdD4F1YGJ4ABGq52nFm8D6GTfBlCIjm08P4tgJgmX8IzXqztNsgAWs4FnHr6Nq+UCUkQTBl0GUcMvhT6T+aYhl7u/rsRobVttWUG/Rcnjgzra3dOB2trNXLz4j+REZ/P+SU+jqd7ijMYtXwubPoe17wDQR/bmbNs4FolxlDWkRTJiHiZEBPphRnmDkWeXbOezwjJ0KiVXT8jn+kkFZMaHuFy3muDDy5x+5bM+hL6nhjyWS4/P5fst1SzfuZ+1D5zO2f/4pU17PoBLUM0YmR1Q2PnT6NyP1aQsdQrzfWdjrT8R8G1q6AoLCssR0MHOvVoO5Cv7OK5S/o9/NpyHROd1QvJn8vB2/XBA+Pqa4B75clPQk4RrMnVnlNjOs6pXyKeKZY6hvGk/i9X2Efw9ezgzBmXzt3lLDwhzj/49BfDg5MHMO2kety29jTdLv+VPw/8EA89yvmm3QWUR/3r7TcZZVjBX/QFz+YAiRx8+sE/hpW/tEYF+iIlULDpMaDBY+PuizUx++ke+XF/BDSf1ZvmcKTx4zuDQhTnA4rlQ8huc/0qXhDk4w9LvnTYAo9XOuyv2BrTf+nrfJcjKG4xIt/8XFJZ3aAMg1LVoUpZibRyBtf6A10U47fZPLd7mddPyXdsZxAkDpyiK2l9zF6z+xuHLFOKJZ38LCss7TZIu/E0SLmYqlvGR5lGUOLjM8heuss7lJ8dwWq2y/Twh2d2BU3JPYXLuZP67+b8YbW5tlCrIGcO8lrOZYfkbJ5qf5zHrLPSYeUL9Ot+ZL4OH46Gx3Gu/EXqeiEA/xDgckvd/L+GUp3/kjV92c97wLH788ynMPWsQSdFdDPrY+R2sfgsm3ArHXdDlsS0oLOcP/14FwIvf7yAhyruvskvI+RJ2SiG8aqD3fLSOBYXl7RozgCZpOUgF5qrp0PZqIPtzqPgSZKtlf1qkjnGKLT7bz546ALWioweRWiGYPDCVYGMv3ftzF+6eZCXo/a5MpilW8qzmVVY4BnGm5XF+dRzn9Tz+JiFfXNj/QposTayvWe/zuDKZxmv2c5hqeYJLLfcfaPDCcFh4mzMzZoSDSkSgH0I2ljcy85Vf+cvnG+ifHss3d5zEUxcN75426rDD4vud6WqnPNDlblx2XZfWaXNImo1W1MqOYstd2M6eOgC9uuNGrQDsPlxj7VIy97MNPPLlpnaNWRW9HXtrP6Q9FoAEvTrkgJwFheVMnLeUgjmLmDhvaYeVAPgWZFIoKZHpZIv9HV5XCNGxD0/JLWDR+kqvWr833M/vT2DPnjrAtzCmlmfUr7DW0Zc/Wu+lhc6uhq5jfX0vkwem+jx3v4R+AJQ2d3Zd7dyfYJ1yKAvO2wx3rIfRV8O6+fDS8fDdw85snREOChEb+iHAYLHx5Lfb+M9ve0iK1vLcJcOZMSI7PPk1tn0DNVvgwrdA3bFAgqcnxeSBqfywtcarXdvbUt8mIUGjIlqr8nrMjJHZrN5bx3srStqFWyAhZ7TaO5xHqBtwtAxq///hc4eELMwDbczOnjqA2Z+sw2r3GJ0EExp0WDq87Jp4wPm5eB5ntUufZhNPPFcbvkw1CXp1+3i9BRK9mLoYVZ2D2y23YvFWEQTntU+ct5TZUwdwwejsTt/L/JWlLFpfSYPB2um7dJlavKUGCJgRc/ozcNI9lHw8h7xfnqPm57d4XnMjx591TcTG3sNEBPpBZvWeOu75eB0ldQauHN+Le84YcMDlLBys/xCi02BQR48Wb4LuvytK2t/3FHy+NMcGo5Wih87wefofttYEral6Qzp0oDxw7jvnF/HU4m1Bu8X525h1n3geXrip0wajA0hRNLPB3qtTv64+Qt2cTYxSE6XxPgGCc3LxJrAfPndI+1hd1+XqY86peYxZ8gOMmoXc3Av8jMn1vWpVik7fi9VxYCLy/P5X7lsJQGV1MhM/Xtpp/IFiABbskszdcxn9bGP4m/ot/m59kq8+X84iy1NMHzfE53ERukdEoB8kTFY7zyzZxhu/7CYnUc+H149nXO/k8J5EStjzizPdrbLjVxtocw06Cj5fmqPAOTn4epi7643iMOagit6BGQcui2Ao7o/BbgA2evEWiaeFPPYx3zHJZx/+NGqzzdFJMD90jv8VRlD53z3H2bTV6V/e/0wqfg1uIzbQd+9q99TibZw5LIV3N79Lpq4Pz3/d1O4hE8r34Lrf1tOHmZZH+JNyIberPmf/N+dB3meQOSzgeCKETsSGfhDYUdXMuS/9wus/72ZC72Rsdsmlr63wat/tFi1VzkCijM4PS7CC1tXOl31VEngjrztYG8egUDeiTljV4XVP7xBf+FrteI7L2zinKZ3nXKce4bMPb/Zol0b9+MyhZCfoEThD8YO1/c8Ymc3yOVPYPW86y+dM6XDM/Qs2cNf8og4eQp//1PbZJOaHPWq3osHIM6ufoaS5hMbyqV7dHe/6qMjn/oR7Py5sqPiHfSYXWh4CaYc3z2Dlorf87nNE6BoRDb2H+WRNGQ8s2EiURsmNk3rzn9/2hhx4E3T0pbnZ+Vuf0OmteL3aqw+zJy4B8cPWGp9t3G2zrnG4xljeYOzk4x0KtuYh2FoL0KYvwm7shcOc0eG8BXMWEa9XIwSdbL8LCstptXSO0lQrBLOnDugUoKNWinZ7uAIH1yi/ZYcjm9/t+agVAqtb9I7L9h1Io+7K9+jvdXe7twuzzQEaQDq8mmy6jiQ5ezkfbP2KqwZfxT8/9V7xyLXH7e/+9baSWSf7cr3uGd6JfpHRK+9mvPVGPmVSSJp/uDDaHbTaHUgkiSoVKkUY9q8OAyICvYcwWGw8+MUmPllTxvjeSbxw6Uhm/vPXgPZdT7zZvmd/so6HF26i0eixmeVKquQS7G59eBN0nrhv2AXS6MsbjNw5v4g75xd1ek9Cu1DPTtCTn6zn1111AYV8u725YhZRBf9An/cGxpLrOgh1ScfgGndh4G3DEiBG57zNPQN01ArRHnV5gXIZgxSl3Ga5FasDEqNUPm3fXckh42uzdvXeOj5dU+51kvflL18q05x/1O5kxsgZ7W0rGozt0byek5FWpfA/oQsb2rRFmON+w9Y4goayM4jXVwVUAnzdv772Bq6bNo6Lv72XhxyP8YzmVewWBQscJwZ8DsBpJqtvtdBksmJzSLQqBdEaFRnxzkyi7nhOkpef1gd7qo6f65vZ2mqi2uN5yNaqGREXxaTEWM5KjSdVc2SmE44I9B6gtM7A9f9ZzbaqZm4/tR93nNoPpUKEHOAB3m3fVrv0Hn4/PANUeqjZ2qkPb4IuSq0gMVrrVWgFGyjjC5cwd2UA9NSOWy22DmNytzc7E2dZUGS9RlSvVzBVzsTWPNznuQJtWDYYrN4/xzahlyNqeED1LisdA/jKMb79mMIHfW/+hoqvzdoPfi/t5NYZ6Hq2yxwM6IjatRSGzOg0wXjT+KGzt4xaIYjRqWi0l6LL+gilrgJL7UmYq8/kvYoyNMrgtFbPcbrOb7TaUQqBXUqy3e6vu+YXcT338G+e5An1a1RYklkpB3Xop9lkZeXuOlYU17KurJFd1S3Utlo8T91ORpyO47LjGZOfiM3u4KWlOzHaHDgy9RTnRvNocz00w+BoHZOT4uit1xKjclqc66w2dhrMFDYZWFTTyNztZZyVGs/NeWmMijuykpIFJdCFENOAF3BWLHpDSjnP4/084B0goa3NHCnl1+Ed6pHB78W1/Om9tdjsDt7+w1hO7n/AFt2VLHzB2L47aDcFJ8H2xTDtifYkXL76MFgdJHp53VeCqVBxP28wQqej2eI0nvhfNI2x/0af8wHWxq2Yq89E2uJ8nsvf5+vrM4jFwOvqZwC423ozsm1bKdy2aV+Toy8ffX/XY0FDbc6pqDZ8zvRNZ7CzUQTtgeL+md96ahYl8gve3fw+0q7DUHoV9pbBB87jRQnwhvtn5bkSsUvZyVzlvC64yXonn2se4iXNPzjd/CTR8Sl8tb6CL4oq+GlbDRa7A41SwXHZcZw2KJ3eqdHsrTXwyZqyDllAVQpBTqKe4poWvtviLEEolQJ7r2hsfeMQJjuqrQ1kGyRL75ns8zqklGxtNfFJVT3vVuznq5pGZqYn8kCfTDK1R0Zlp2Bqiipx1hQ9HSjDWVN0lpRys1ub14BCKeUrQojBwNdSynx//R6N+dA/WFnCAws2kpcUxRtXj6F3akz7ewsKy726ygUqTuAr/7Un7bm2N34Kn/wRLnkPBp3ttw9PW7dereSC0dkdTADdIRwFHawOK+P/ORdz9HcglVhqJ2OpnwCOjj72Lg3QV/EHl33fnWiMvKl5mtFiO3+w3ssvjqEdjgmnPbfP3K99Cm9v+LoeAVw+Po/T4io4ZdnFPGO9kH/YZ4Y07npTPe9ufpcPtn5Aq7UVc/1YzDVngD10bdTznL7utcQodfuKx13oDxZ7+ELzAIscE3hIdQeNRivpcVrOHpbFqQPTGNUrsYM5JVA++D31Bk76zwqEyY6i3oJsu8kFoeWjb7HZebmkmn+WVqNXKHh2YC5npSaE+vH0CN3Nhz4W2CmlLG7r7EPgPGCzWxsJuFSneKCi68M98nA4JPO+3cpry4qZ1D+Vf8wa2cHbwluVGXDe5O5mBm8aq88gGA/ataRB50FSH1hyvzOHi9q3YPDs0ZcJAECjFCgViqAFvWsj0huhpNhVK9Q8MPEu5n45HJK+RJu2GE3yT1jqJ2Ctn4C0xQW9Yen+GSTRxL81TzJE7OFu683twjxBrw45mCkY/AlzvVrZaRIKdD0T59Vwv/14blF9wVeOCeyWmR1Wat4+4365DXy0/SO+2f0NJpuJ03qdRn/N+cxb6D2S02Vw8TVypRCdJhBfK6F6g7Xd3dXV/slvt7K5MZ/X7Gdzi+oLVqZezumnnsakfqkofWxS+jNbbm81cdW2PViHJqLc3YKotyDcBp8So/VxJZ2JUSm5r3cmF2UkcdPmPfxx4x5uyk3lwT5ZKA7jAtvBCPRswD3+twzwLML4MLBECHEbEA2c5q0jIcQNwA0AeXned9CPNKx2B/d+sp7PC8u5akIvHjx7MCplR29QXz7gURpVUGlnvWn27nSIPlSq4Jzn4Z1z4Os/w7kveRUMoZoALHbJ8xd613S9EaNTeRWKoabYbbfHGpJQGq/BrC0lOm0Z2uQf0ST/iNo8mIsGzGT6sLT2Prz14/4ZpDWu52XNiyTRxI3Wu/jeMbq9XbTW+7i7S7aPz9ylifszP3kbT0WDkQe5hv9p7+UF9UtcbHkQE879EPfPWKjrqFFs4P5VTyDWV6BX6Tmr4CyuHHwlfRL6MHHeUp9jFl7S9LrjkDIo7xYX7pueAzNjSYjSUNFoYnX2ldjrf+CxtO9hwGXt7b1NSj7jAHrFMX3NdnRKBX+OTeTNPVW4txJAncHCf1fs5YrxnQPHfNE7SstXo/rx0M4KXi2todJs5cVBeWj91BQ4lIRrU3QW8LaU8hkhxATgXSHEcVLKDk6sUsrXgNfAaXIJ07kPGQaLjT/9dy0/ba/hz2f055bJfb2G7wfaDA0U3egtCMaFAC4Y7fHQF0xyFrJY9iTE58HJ93YSDL6Wrq5NLE8y43Qd+lhQWM7/fb6BVot3jb3BRyh8MJGcLrzaY+29+OuEJxnVx87nOz/ni51f8MHev/JVxfNMyp3ElNwpnJh9otcyajOGpTOj9SNY+ndKbQlcZH2IDbJ3hzbhTNPrji9zkLvdOxRcdug/W2/iNfWzPK1+ldust5GZoGfe0iXY4zcQFbMZpd7p32035hDdeCHf33Q3sZrY9n78Xa8/YQ7eff5nTx3g1fPJdS4pJW8t38MT32wlTq/ilctHMe24DMSiC2Hdh2BpBU20z4nfm0lQmRlF9cBYBuo1vDO0Nzk6DX1V6g6Twc2n9OG7LVXcv2AjjUYrt0zu6//i3NAoFDzWL5scnYa/7qrA7HDw+pCCTknaDgeCEejlQK7b/zltr7lzLTANQEr5mxBCB6QA1eEY5OFIo8HK1f9eyfqyBubNHMqlY32vOHxpFa6kT4EEvj+tR+LDZ/yUudBYCj8+BuYmOP1RUBywRfoSMN4eGIWA+84c2KF7lxAa+egSr3lMfG0qhuLp40/4L58zhTtG3cEtI27h14pfWbJnCT+V/cSi4kWoFWqGpQ5jbMZYjs84nmGpw9BWbYaFtzuLNQw6h+uLL2arpXPFp54qr9eViFB/uL6/76wjmCfO4i/KRdijynk8O5ZWRytanELcVHUmtuahSGsSRuggzKF73kzeLA/+VpSZ8Tr+/PF6Pl1bxmmD0nnigqEku8wg/c90ZgitKIL8iT6/+x+21rTviVQ0GIkriKOmfywj4qL4cHgf4tpKMXqbJC85PpfZn6znqcXbiNGquPqE/BCuVXBLXho6heD/dpRz19YSXhyUd9iZX4IR6KuAfkKIApyC/FLgMo82JcCpwNtCiEGADvAdmXKE02iwcsWbv7NtXzOvXDGaqUMy/Lb3Ffxhl5K75hf5tFG6Z8vzFzziVUgqFHDeP0EbC7+9BFUbYebrEHPANAHeBcyYXkkdTCuXj+vlU/A8dM4Qn5qnr2sK1tMnGOGvUqiYlDOJSTmTsDlsFFYX8lPpT6yqWsWr617llXWvoEEwwGRiMCoGT76DQcfN4ro9eh5YsDXocYeDrmji7ljsFoobi9lRv4Ni+w76Dl/NnuZtfKCwEFcfz60NZfRp7sufzLMo298HaY/pcLy3z9jXxB7Qbx3fq7CHz+18TwBUNJr4dG0Z04Zk8M/LR6Fw13DT2hSG/dshf6Lf7971Oa5oaOHiol2MitXzwfA+xPqrqwuolAqevmg4zSYbj361mQEZsYwPMf3GtTmpNNvszNu9j95RWu7O9//sH2wCCnQppU0IcSuwGKdL4ltSyk1CiEeB1VLKhcA9wOtCiLtwKo3XyEDuM0cojQYrl7+5gu37WvjXlaOZPDAt4DGuh/iej9Z1Mmf4+pDchYu/48H7g3rAq+ZkLlLCX4v/jXhxLNppf4URl4NC4dfePGNkNtf/ZzWFJfXcf/agTm08ry1YzdOf6cHbdYXi5qlSqDg+43iOzzge9u+gafnzrN2+gNU6HZtS8vjKYWD+ns9hz+eohIqMIRk0NCZgaE0kTpnFrFGjGNnHhtFmRK/q2ULYvjDZTFS2VlLWXEZZS5nzd3MZu5t2U9JUgl06Pze1Qs2AxAFcnncRw1OHMzx1OHL9pxy35AG+jF/ITMWtFLvJU1+fsa/vD7yX03PH1/fgyrrpa4P9p+01LFxX0fEecZnI7Nb2vr199xKnufCq0/rwvLGJPL2G/w7rHVCYu1AqBM9dMpxzX1rO7E/WseTOk4OrzevGHb3S2WU08+TufQyJ0TM1JYSykD1MQLfFnuJIdFvcs7+VU57+EY1SEbQwd6dgzqKgQuITo9RISadIUG/eMt5c1RYUljP743UdogX7iTIeV7/JGMU26hKG8tfW81nQPICshCivAnhvbSuTn/6Rm07uw73TOppbukuwXi7BXm87DgfsWQa/v+YslKzUwIjL4OR7IS4Lh3RQ2lzK5trNbKvbxt6mvexp2kNJUwkWR8eglVhNLOlR6aRFpRGvjSdOE9f+E6uJJVodjVqpRq1Qo1FqnL8VGhRCgU3asDvs2KXzx+FwYLQbMVgNtFhbaLW20mptpdnSTJ2pjlpjLbWmWmqNtbRYWzqMQ6vUkhOTQ25cLv0S+tE/sT/9EvuRF5eHWuElmnH7Yvj0OiwO+BvX8W7zaJ/fcbDfk7d0Dv6+B19eXe50cmmtXA//OsmZ9vm4C/z2IZUC64Q0omPUfDduIL30wXuvuFhRXMulr63g9il9ufuM0FdlRruDGYU72GO08OPYAQfVT92f22JEoAdJi9nGcQ8tBuCpC4dx0ZjcAEccwP3BCAZvbmyuhycYYejb79zBlfpfuckxnyxRy1pHX960ncXPynE8OnNkh37u/WQdXxRV8PN9k0mL1XXq62ARlPCvK4aiD2DdB859A30SHH8djL2+3cTkD4d0UNVaRUlzCVWGKqoN1VS1On/XGGtoMDfQbGmm2dLcriF3F6VQEqOJIUmXRLIumRR9Csn6ZJJ1yWREZ5ATm0NOTA4p+pTQ8+TX7oLPb4SyVTB4BkybB3GZ3RpvKK6mwcROePqEb/z4rxy36WlOMj9HqUwnMUrN9GGZ/LC1pkNfErAOTcSRqSdrSzNrb/GeGTMYbn5vDT9v38/yuVOI04Ue6r/bYGbKqm2Mi4/mg+G9w1PPIAgiAr2bWGwO/vD2SpbvrOWmk/sw58zgNdZgtBV3fHmZhBKkE2gloMHKBcpl/Em5kDxFDftlHEtUk7nsj3dA1ihK641MfvpHrhjfqz0v92GFlFC5DrYucmriVRsBAX0mO81JA6eDOvwmEyklBpuBJnMTrdZWrA4rFocFq9352+ZwauYqhQqlQolStP0olOiUOqLV0e0/WqW2ZwWA3Qa/vgA/zgOFGibdA+Nv6VT0JFSCEezBrETd7+eFq3cz6svTqJKJXGB5pL2NWil46sLhHfaZbFlR2IYmotrZhHpXc6dAoVAmng1ljZzz0i88fM5grplYENLn4OKd8v3ct72MJ/rncHV2Spf6CJXuBhYd00gpue/T9SzfWcszFw3ngtE5QR0XqlYO/su1heJOF8hzwYKaD+ynMt8+mZMUG5ilXMrFti/h9S8gLoe9qnGcpijgpnEjgz5njyKlc7Nsz8/OfO97foHWGhAKyB0PZ/wNhsyE+J7N1CeEaBfIhz1KFZx0Dww5H5Y8AN8/CqvehIl3wKirujTh+Usw5l75KlBmT097fstXc8kR+7nXekOHdla75J6P1pEQpabeYEVqFdgGxSPqzCh3NXey4Yca4zA0J56BGbEs2lDZZYF+VVYyX9U08FhxJdNTE0jRHFqRGhHoAXh6yTY+Lyznz2f0D0mYdyWlqcS3hh6KO93sqQM62dDBqfFEa1TtD5sDBT85hvOTYzgD4218e6aBpsLPGLNnISeqrPDKM5B+HOQe7/ydMRTSBoM2xttpw4OpCep3Q/UWp121ch3s2wDmRuf7cdnQZwoUnAz9p0L0wdGKjliSesOl70HxT/Dj4/DNvbDsaRh3A4y8CmLTg+7Klyuhe4rf8gYjauWBLJYu3LNvtmvNUrJ5/oNcxje8aTuzU5FrcCo4jUar09QyOAEEqDfWE+Vlk9efm6vrfU/NfeqQDF5cuoNGo7VLlcOEEPy9Xw5TVm1lXnElTw8M3hTbE0QEuh8+WFnCyz/sYtbYvJACEYKpDuQNXxp6qO50Lm3E3R/YlWYAvNeovGnaCOTwLK77PZ8S5TV8f2k00RW/w95fYMMnTh9h1yhjMyEhFxLyID7XaafWJTjzsOsSQBMFQgkKlfMHCVYj2ExgNYDFAIb9Ti27dT+0VENDidMWbnAr0KzSQ/oQGHohZI2A/BMhscC7A3QE//Q+2fmzZzksewqW/s1pjhlwFoy62vmesqNA8zRf+IuFcMeVpkKnUmC2ObybPlprYdFdDN76BZ/aT+Ixm6cn9AEcEhzpehxpelRbG1Ea7VwwPq+T1u1rFevS1L1p7iPzEpASNlU0ckKfrikH/aN1XJuTymulNVyTncxxsZ2D2g4WEYHug9V76njwi42c3D+Vv543JCR7Z1eiDX0VhfCWLyMYAvk8e9NWvlxXwcrddfxtxnFE9+8F/U8B7nOaPBpKnLbqqk1Qt9u5+Vi6EjZ9Do5uZGbUxDq17Pgcp+07qTckFUBKf0ju16mUXoRukj/R+bN/B6x5G4rehy0LnRPxgLNg0DmQfyILtjR3EoKhFi6J06upaTZ3fLG11qkc/PoiWI3Ms17Kv+xnt2e59IZUgHVAHKLJgrKkxWcwna9JRymET839s5tPAGBndUuXBTrAPfkZfFhZx1N79vHO0N6BD+ghIk+LF6qaTPzpvbVkJ+h5cdbITrlZAuFPm/H2UCT4sTl6y5fRXbwJ+/pWCw8v3MSwnHguPd5j2SgEJPZy/gz0yFbnsIOpEUwNzvJ3xganNi7tTkHvcADSabNV6Z2/1XqISoboVKc2H+Hgk9IPpv4dTn0Qdi2FzV84N5nXvQ9CSV/6cKscSJGiL5sdvSgnBYkISahXtwlzS0MlKz7/H8ev2kl21Y9gt8CA6XDqg3z5VgUygAJk7xUDehXqDTXtyba8KU2+Yhx8rZbLG4wkRTvdDev85FoPhjiVkptyU3li9z6KmgyMiDs093VEoHtgsTn403/X0Gq28d9rx3XJruYvstPbwxCtVRGtVflMD1AwZ1G3Q8UD8devNtNotPLe9eNCm8AUSohKcv5EOPJQaWHAmc4fmwVKfoPdyzD9tJAblItQq5z3cKOMYqfMplymUCFTqZSJNEs9regwoUWJHTU29FhIFQ2ki3oKxD4GKfaSJeoAqK1IgHHXwuirIc0ZrDZ7aqzf/SapUWDrHYuiyoii/oDQ9ban5CtIypdzggAWra9Ep1aEJf//dW1ml6f37OO/ww6Nlh4R6B48+e1W1pY08PJloxiQERv4AC+4bqxHvtzkNc+JJxUNRi4fn8d/V5R0es9lU+/Juos/bK3ms8Jybj+1HwMzvBeQiHAMoNK029rvWDWRuoZ6BopSBiv2MljspUBUMlpVzDTHSjTC/x6RQWoplan87hjEZkcvfnUcxxaZR/GZ53Ro522/x30VIPrGgUKg2n4gxW+oe0qzpw7wmmJD4nzezTZHp8LfXSFWpeS6nFSe2rOPnQYTfaMOfvxGRKC7sWx7DW/8spurJvRi+rDQAzE8N5GCdfHPStD7LcrswltWwlD8br21PaFvMrM/WcfAjFhumdwnuAFHOOpxrTILrf0otPcD2gLczh3K099uwdhYTbQwEY0JPWasqEiMjSY5IZYfywX1dh0HMqo7yfbjqWW2HUjM6npspEaBOSuK/jYFdo2aCoONrAQ9kwem8siXm9qzOrpy2ANeNz8fnznUp5mootEEQGwXAou8cVV2Mi/sreKNsv3M6x+cV1w4iQj0Nva3mLn7o3X0T4/hL2f5zl3ijr86maFEhbo0iGBwtx2G4nfrre2cT9fTKzmaFrOND64fjzbIfBgRjn4CFwuxUethq35s2nE88MUmRhQksGZvfdCJz3x5hdnyY0ABpg31/MXNdOK5km0wWpn98TpidCqfm5++8tGnxmipaTGTnxKe2IJUjZrz0xP5aF8dcwsyiFcfXBF7eGZpP8hIKZn72QaaTFZeuHRkpwri3nAJyPIGY3sl+kBVhTzJTtC3e7AkRAWnIbjbDgP53brjra3J5mBbVTMPnTOEfuldMy9FOHqZMTKb5XOmsHvedJbPmdKh6MbjM4eSnaBHcOA+njIonRazjZP7p3p933NlOXHeUgrmLPKehEutwJ4bjaLSSHVVa/uz5gurQ/o0b1Y0GJk9dUAns4pereSUAc6av/3SwhdbcX1OCga7g4/21Yetz2A5pjV0z2jOc4dnMSizsw3Zm6miq77m4CzP9tRFw9tzs0yctzQoW7unlhNKbnF/rpSdvFoiRAiAN0+pXTXOxGKpsVq/brPBBN7Zc6JApUC1u9mr22EoZCXofa44lmzeR2a8jl7J4fNKOS42imGxej7eV8f1uamBDwgjx6xA93ZTLdm0r73uoa92nkEKgUhs07zdBXa0VuVzDN4Q4NU+Hkp6WV9tM+N1By2pUISjG9eekSJAJZ9AypAUYMuNRlFrItosMYaQb8pbVkj3NNTuz4/RYuf+BRuZPjQz7M/AxRlJ3L+jnC0tRgbFHLxUzMesycWXCcLTXOHLrKEM4gbQq5U8dM4QHjpnSIflXoPRytzPNvDIl5sCCvPsBH2nJa8LX8tIb7ZKb221KgX3hTk1boRjF1WbILfYHH7bBQq8c6TqQK8irdbabroJFsmBrVhvph53vt1USYvZxvmjwu8KPCMtEZWAjw+y2SUoDV0IMQ14AWeBizeklPM83n8OmNz2bxSQJqVMCOM4w06w5gpf7exSdgpaUCsEMToVDYaOecwnzlvqdVIIJMwDuWeFUlxixshsHA7JXxZswGR1kByt4YGzB/eYX3uEY4/0OKeb3r7Gjs+Mp8nSlWzLF/a8aFRmO6uvP7F9kghlVezKGeMvO6mUkjd/2U1+chRj88MfQ5GiUTElKY4F1fU80Cf8KwBfBBToQggl8DJwOlAGrBJCLJRSbna1kVLe5db+NuAwSdPnm8x4XbvLkjue5gpfpopA1drd6UoqgOwALogugi1rJqWksLQBk9XBQ+cM5g9dzC4XIYIv9Bol6XFadtW0tr/mzWSpVgjUSuHVicChV+JI1nG2LqpdmLsrLuUNxvYEdr48V8D7M+c+sSRGa6hrtfDkhcMCmoi6ylmp8SypbWJ9i5HhBym/SzAa+lhgp5SyGEAI8SFwHrDZR/tZwEPhGV7PcfbwLF5bVtzhNW8acVertbvfPAofGRQT9GrMNkfwFXm6wT9/3MW7K/Zy46TeEWEeoccY3SuRlbvrkFIihPBqsrQ6JAl6NdFaVbvLrxDOGqXRBXFYgAdH5nc4xtez5quYRqDUunWtFgQEZTrtKqcnx6OglG9rGg8rgZ4NlLr9XwaM89ZQCNELKACW+nj/BuAGgLy8vJAGGk6klKworiU5WoNWpaCy0eRTw+5KtXbPm8dXBkVXMES4KsH74qPVpc6ApBFZEZt5hB5lQp8Uvt6wjy2VzQzOivO5Om00Wil66IwOr0kpOeH3LQzSasjRHSjp5i94Ltgatd4mFgk8+7/tQafFDpVkjYpxCdF8u7+R+3p3r2JUsITby+VS4BMpvdfpklK+BrwGzopFYT530CzfWcv6skYenzmUWWMDTyyhVmv3tYuvFAKHlJ1uyp60Y39eWMZ9n67npH4pPHnh8B5bXkaIADB9aCaPfrmJj1aX8vC5Q3wWu/CWI2ltk4HdRgu39zqQoz1Q8FywClcoLr7hZFpKPA/trKDEaCavC7VPQyUYgV4OuDsq57S95o1LgVu6O6ie5uUfdpIep2VmD+xug++bxCFlp5JZPcnCdRXc89E6xhck89qVY9CojlmnpggHiaRoDWcNzeSj1aXcMrmvz9T13l7/oroBrUJwdmpC+2v+gufcFaJASlEoLr7h5OQkZ8DeLw0tXHYQBHowT/gqoJ8QokAIocEptBd6NhJCDAQSgd/CO8TwsqWyid+Ka7n2xIIeC3X3dZP09M3jzqL1ldw1v4gx+Um8ec0Y9JpIWH+Eg8Mdp/bDbHPw3HfbafDhzeL5upSSb/Y3Mikxlli357IrmrV7FOrEeUtZUFjODZN64zmHhJrkqysMiNKRqlHxS31Lj57HRUCBLqW0AbcCi4EtwEdSyk1CiEeFEOe6Nb0U+FAeqqrTQfLhyhI0SgUXje656MhQ/MN7gk/XlHH7h4WMzE3g39ccT9QhrnMY4diid2oM15yQz/u/l7TnG/fEU7nZ1GKk1GThzJR4v+0Cve6ZkqO8wch9n67n5R92olEpSI3R+kxH0BMIITgxIYZf6ps5GKIxqCddSvk18LXHaw96/P9w+IbVM5isdj4vLGfacRkk+rjRwkFXNlLDxZu/7OavX21mYt9k/nXlmPao1AgRDiazpw5g2fYaKhqNaNtK0blwV25cG557k1XQJxZbpQGykjv0E8ympwtvJhqzzUFNi5mPbpzA8T3gcx6IExNj+by6gR0GM/2jezal7jH1tP+yYz9NJhsXhrCrHUp6WndC3UjtLlJKnv3fdv6xdCfThmTwwqwRkeyJEQ4ZOrWS164awwWv/IpWBQl6QXWzucMz5L7haR+Uimiw8PgPm4hTKDo5DAT7DPoyxUjJIRHmAOMSnJkcVze1RgR6OFmyeR+xOhXjeycHboz3HfbZH6/jkYUbMba00itacu1JBZw5Io/oxEQUCu8C1H1SSIhSY7LaMVqdGoureHN3hL/V7uDBLzbywcpSLhmTy2Mzh6KMeLNEOMQUpETz9h+O56q3VmKX8OnNJzAqL7H9fZc2LTUKZJwG1fZGrzn/fSlH3pQtX5ufoq39oYiM7q3XEqdSUNRk4LLM4GRPVzlmBLqUku+3VDN5QFrQ3h7tyzcpyTTvo8CwhwxTFamWGjTSWbKqdIvTD1OhVJGYmUXukKH0GTWWvGEjUCiUnSYFz5DneoOV2Z+sA7rmvthgsHDze2v5dVctt07uyz1n9I8k24pw2DAsJ4FPbjqBP7y9kote/Y3bpvTl5lP6olEp2rVpR7LT+0NR66xBWtFgDLgy9pXfv296jPd0vNBponCnqyvxYFAIwYjYKAqbDGHpzx/HjEDfvb+V2lYLJ/QJfoasq63n+KaNHNe8mRh7K3YUVGtT2RI7iCZVLEalHgmk6QQ3jk6kevcuNv34PUWLFxGbksrx517A02u1AXNQWO3S783mi+KaFq59ZzVl9QaeuWh4jwVIRIjQHfqmxfDVbSfx0Bcbef67HXy2tpx7zujfnn7DkaQFqwPR5FR24vXqgIVbfCXX21jehC98mWNCKRTTVUbGRfNSSRVGuwN9iEXnQ+GYEehrSxoAGNUr0X9DwGIysnLBx1xT9hkqh429+lyWx4xnd1Q+VkXnzdQdwIUDR/Ds3lz2ZY5hNBWcbNvK0rde5RR1PP9LmUKVLsPvOYOtcOTip+013P5BIUqF4P3rxx8y+2CECMEQr1fz/KUjOX9UDo9/vYU7PiwiXq9GqRCYtQoUtSYEzg1PIQjoe+7veVEIcHhxKPFVRCYYX/fuMjI2CruEjS1Gjo8PT3UkbxwzAn1LZRN6tZK+qf4rk+xY9Rvfv/kKrfV1JAw+nrdb+1KpSPB7TEKUm0YhlKwkl5WqHIbkDWZs5Q9cWLmA5UkTKIof7rOPYHNK2B2SF77fwT+W7qB/WiyvXzWGvDAm548QoSc5uX8qJ/ZNYenWat5dsZdl22vQrq1DKgRalYIT+6Xwv81VXo8tbzBy7dur2FThWwsXeBfmgM8av6H4unfVNDM4xrkZurU1ItDDQkmdgdwkvc/Qd6vJxA/vvMaGpUtIy+/DuXfPJav/INL91A0Fp0YhZWeNAiHYpMxmT87FnFz9AyfV/Uq0rZXlSRO8hsl5y/fiSW2LmTs+LOKXnfu5YFQOf5txXCRgKMIRh1IhOH1wOqcPTueNXft46LddTBM6yqpb+X6Ld2HuorTewIQ+ySgEfLW+soM7pGdxC08avaQggOCjSLtjmsnRadArFGxv7ZzhNZwcMwJ9X6OJzHjvwQjNtfv5/MlHqdm7m7EzLuKEiy5DqXIuzzx32L3N0P4KPLei4bde06H6Z0btL8Ki1LAqYUyndoGS+P+6cz93f7SOOoOFJy4YysVjciObnxGOeLbbrETlxfL6SUNRCIHV7uC9FXt5/JutHYS1TqXgsfOHMtNtn+ikfqkdnsVAZktfwUiTB6by3ooSn5WOXHTHNKMQgn7RWra3mv226y7HjEA3We1Eaztrs9V7ivl83sNYTEbOv+9Beo883m8/LgHvEux3zS/ymR7XRYPRxssv/ZVv//kcLFtKqy6ZzboDKWz9BUqYrHae/HYbby3fTe+UaN685gSGZMV7bRshwpHG2iYDI+OiULQpJ2qlgmsmFpAQpQlo2vBUtvrM/drnc+jrGVtQWM6na8o7CHMBXDC6s6tkdxN89Y/SsbyhZ1MAHDMC3eaQKBUdd5er9xTz8V//D5VWy6WPPkVqXn5QNrJg0uO6kxCl5ouiCp5uHMoE7UZOrPqBxr6ZVFh1fu1wG8sbuXN+ETurW7h6Qi/mnDkoYmKJcNRgsDvY0mrktrz0Tu91JTDP33PoK8zfV1rdH7bWdGrb3QRfA6J1fFJVT5PNTlwPBf0dM+n34nSqDja02rISPv7b/ai0Wi55aF67MPfMAzH3sw0sKOyYXNJXelxfBpBGo5XZH6+jrMnCktTTENLOqIqfeO6SEV5rhRotduZ9s5UZLy+nxWTj3WvH8sh5EXt5hKOL9c0G7BJGxYVnU9+X2TI7QR9yZKm317uboym/LdtiibHnzC7HjEBPjdVR2fYlGRob+GzeIyiVSi5+8DES0p0uhf5sZO74stVJnFWIPHFIZ5UWgEZ1PGviR5HfUsy/P/muU9sftlZz+nM/8epPu5g5Kptv7zyJk/qlhny9ESIc7mxodj5HI8JUzacrAjeU5F8zRma3F63uSoIvV9GOcrPveqrd5ZgxuQzOjGXp1iqaW4189cxjGBrqueTheSRmZLW3CXa2VvqwmSuF8LmT7k5h/HCGNm+id9ly4GoASusMPP7NFr7esI++aTHMv2E844JMURAhwpHIllYjSWolqWHKBhpM3hdPk+rkgal8uqY86ORf3cnRlKNzKnulJkuXjg+GY0agj8lPwiHhs9deZ/+2zZx9531k9O3foU2wNjJftjq7lD6FfYd2ChXr4oYysX4Fxdt3MH+XjX8v34MA/nxGf26Y1CdSjCLCUc+WFhODovVh9dby5rTw1OJt7QLa0+3w0zXlXDA6mx+21vR4ZtQUtQqdQlAWEejdZ0KfZAY5Ktm/YgnDTz+LARNO6tQmmFSdCwrLffq7CrwLe7VSgJvZBWBX4hAmNK1h3nP/5n8JE5k5Moc/T+3v07UyQoSjCYeUbDOYuCwz/BHOCwrLmf3JuvZ4kfIGI7M/WUe0RuXVpPrD1hqWz5nitZ9w5ncRQpCt1UQEejiQFhOn1SylTp1I2mkXem0TzJLtqcXbfAYveHtdKQRPXTi8/djyBiNxOhUOVBRrc+jdUszCv8xmaG7glAQRIhwtlJosGOwOBkaHX4F55MtNHYL/wJkvyVttU/AdEdoT+V0ytWqqzLYuHx+IoNb1QohpQohtQoidQog5PtpcLITYLITYJIR4P7zD7D6/fvw+GJtZlXsGc77YhtnmPWHWjJHZLJ8zhd3zpnv1QAm1qKxDSmaMzOaEvslccnwuSdEamkw2huXEc9Y5U1FbWkg27uvydUWIcCSyrS1ickAP5Af3zGgaCIUQHcrVQfAOEqGSrFFRZ+05gR5QQxdCKIGXgdOBMmCVEGKhlHKzW5t+wFxgopSyXgiR1lMD7gr7S/aw9puFDJ1yBkMnns4N767hzg+L+MeskagCZD7zXHYlRKm93jC+bOcpMVru/LCQRRsqsdolpwxI5bYpfRndKwlDUyPb579K2ZZNZPUfFLbrjRDhcKekzeyQr++5ymHe0KuVnQS167l118K7G0TkiyR1zwr0YDT0scBOKWWxlNICfAic59HmeuBlKWU9gJSyOrzD7B6/zH8XjU7PSbOu5owhGdw/fRDfbNzH1f9eyf4W3z6h3vzSW0w2p03cDb1ayaxxuZ1cpgRQ02Lm+y3VXDG+F0vvOZm3/zCW0b2cdsOouHgSs3Io37op3JccIcJhzV6jmSilghR1+K2+3lyHXa+7ux16S4jn0sJ7qtB7klpJvc2OzVcGsW4SjEDPBkrd/i9re82d/kB/IcRyIcQKIcQ0bx0JIW4QQqwWQqyuqekcidUT7Nu1g12rf2f02TPQx8YBcN1JvXnywmGs2lPPlKd/5PVlxTSbOmvd3pZdVockWqPq4Iv6txnHcf7IbE7ql9JB2PdOjeax84fy219O5aFzhtDbI9PjgsJy1hljWbduS4flXoQIRzt7jRZ66TQ9ko/o4XOHoPZIwqdWCB4+d0gHk6rDhzdaRYPRq0+7wKnUdedZTW6bwOptPaOlh2t6VAH9gFOAHGCZEGKolLLBvZGU8jWcBX4YM2ZMz5fABn779AN0MbGMOrPjouLiMbmMykvgkS838/evt/Dcd9s5bVA6k/qnMiwnnvzkaJ/LqwajlfNHZqPTKNlQ1sjDX26i2WRDIWBkXiJnDE7nnOFZfmdzl/Z/nIwhz97CvrrmsCfVjxDhcGWvyUJBD5lbPJ0b4vVqhKCDC+OMkdkB3ZS1KkUHhc4lsLqzQZrkEuhWO6ka7yuJ7hCMQC8Hct3+z2l7zZ0y4HcppRXYLYTYjlPArwrLKLtIw75KiteuYvzMS9BGdY5G65sWy7vXjmNdaQPv/b6XpVtrWLiuAnAmyRfCdw7lf/+6B4WAwVlxnDM8i4l9UjixbwrxPpLoe+LS/hvVzlVDnK2JekVSWJPqR4hwOCKlpMRo5uTE2B47h7s/ui9vFV9uypMHpnZ63ZOuFsCIatuzM9gdAVp2jWAE+iqgnxCiAKcgvxS4zKPNAmAW8G8hRApOE0xxGMfZJYr+9zUKhYLhp53pt93w3ASG5ybw2doy5n2zlepmM9EaFZkJOnZVt/oMFMqI0/HVbZ392YPBpf0blU5tQOcwd3g9QoSjlUabHaNDkqUNv4bqiS9vlXs+WodDSuL1anRqBQ0Ga7ubsq9cTZ505Vl1CXSjo2cEekAbupTSBtwKLAa2AB9JKTcJIR4VQpzb1mwxUCuE2Az8AMyWUtb2yIiDxG6zsemn7+l7/ARikgKH0C8oLOf/Pt9IdbNTsDabbeyuaUWr8m3jq2zserJ617LOrHAm7NG2CfTubrpEiHA4s6CwnNNf/gWAfy7e3uP7Rr6Erl1KJE7zqcnq6JAoL1hB3ZVnVd+W8dXYQxp6UH7oUsqvpZT9pZR9pJR/b3vtQSnlwra/pZTybinlYCnlUCnlhz0y2hAo2VCEqbmJwZMmB9Xe1waower7g++O8HVtuliFU0tRO6whZW6LEOFIw2X+2NeWnKq+3ns203ASzDPq6V8ezDFdfVb1h1pDP1LZ9tvPaKOi6TVsVFDtQ10+dVf4ujK3ZcQ6BXpCtC6kzG0RIhxptCtNbYVmhNkRlmAdf3jzVvGGy3ulYM4iDBZbJy8Zd3wVwAgGnaJnbehHpUCXDge71qykz+ixqNROgbmgsLz9C/PmdhSKth1q2kxfzBiZzdtXjwbgsQuHR4R5hKMal9IktU6xIyz2Dq/3BC7FKRhc8Sb1BisI3/7svgpgBIPLgmsLooZwVzgqBXr1nmJMLc30Gu7UzoMpXBHsTJ6doPeaEqCrWM1OO7xaow1LfxEiHK64lCapVrYVCZAdXu8pZozMDliz1xNXLhhfenpXJ6GergN8VAr0/73+EgB5Q4YBwedl0KkPfBx6taJTRKhaKWg123xq+V2htaEBgOjE8GedixDhcKJdaVIJsEsE3Tdd+sN9VR7IjOKNBqOVeB9a+uHqvHBUCvSq4p0A7d4tgfIyuDT4jjlaBJccn0uim1+5K2ObS8u/c34RIx5Z0i3B3tpQB0B0QiTbYoSjG5f5Qx+lRtgcYTNdesNzVe5uRglFrDeZrJ0mgsPZeeGoFOixKan0Hzex/f9AeRl8afCL1ldi8uPlAs5ZvDs79Q1Vlai1OvQxPRdkESHC4cKMkdmcNCiNgSkxYTVdeuLVa80uidaq2D1vetD9OCQdJoKenITCwVEn0K0mE837a0jtVdD+WqBag740+HqDNagAg+7s1NeW7iU5Nw+hOOq+iggRvNJssxPbQ1XvXfh6pl3eLFHq4J8394mgu5OQ7KHNUBdHXYGLukqnppyUndP+WqDCFb5yOoSCryT5/oplSCmpKdlLn9HjunXuCBGOJH5paOnxc/h7pssbjKgVAoVo08CDIFyeOOa2E+p7SIE76gR6U00VAPFpGR1e91fc1VdOB61K4bPKiSeeZp1gKp40Vu3D2NRIeu++QZ0jQoQIweHtmXbH6pDo1QqSorUdCkZ/8Hup11Qf4doEdQUU6QPUYegqR51ANzY3AaCPiw/6GF8aPBAwSQ943yTx51njOl/pFqeAzx18XNBjjRDhSKdvlJbBMT3vqggHyj56w2h1dFo1j+mVFLCucHdwhfxHNPQgMTa5BHpcSMf50+BdN4WrKlFilBopodFo9Vk8NpiKJ6WbNqCPiycpO9dr2wgRjkasDok6zP7YvsybM0ZmM3HeUp9C3TNjYjB1hbtDREMPEavZjBCKsAXqzBiZzeq9dby3oqR9KVZvsKJWCuL1aioajO0bou5feqBcy3ableK1K+kzamyPBxtEiHA4YZMSVRjv+UDmzdlTB3Dn/CKvx5Y3GFlQWN5JqPeUF8sBDb1nnvmjTqADvsO7gsBzpp88MJX3VpTgaVVzryLuzT7uyy7vWrqVblyPubWVfuNP7PpgI0Q4ArHLAyHw4SCQeXPGyGwe+XKTz+LRrmfX1VdPaOYumtsEekwPefkcnb5yXfQM8pYiwJsw94an66IriMK9VJ27/+rW335GrdOTP2xk1wYbIcIRitIZKBo2gjFvPnTOEJ+pPYxWOw8v3BQwPUg4qLM4S88l9UAtVTgKNXS1ToeUDqxmE2qtLqRjvc30odx3njeWr6WbqbWFbb/+zKATT0alObhVzyNEONQohQhrcqpA5k04sHL2ZXrx5s3W1apE/qiz2lAJiO0hG3pQvQohpgkhtgkhdgoh5nh5/xohRI0Qoqjt57rwDzU4otq8W1zeLqHQXV/TYF2bNi/7AZvFzPDTz+rW+SJEOBJRhVmgBwocdCdUS0+4M0HWWe0kqVU9tm8WUEMXQiiBl4HTcdYOXSWEWCil3OzRdL6U8tYeGGNIuNwVWxvqiUtJC+nYUAKMPIMSgnVtctjtFH67kIy+/SP+5xGOScIt0IP1THlq8baQrbHhTsJVb7OR2EPmFgjO5DIW2CmlLAYQQnwInAd4CvTDgsTMLADqK8rJ7OtbwHpzcwoUjODi+UtGAF3bQNm6/Cca9lVy7p//GPxFRYhwFKEKsw0dgvNMOdhFbLxRa7GRFESa7q4SjEDPBkrd/i8DvMWqXyCEmARsB+6SUpZ6NhBC3ADcAJCXlxf6aIMgIT0ThVLF/rISn218uTk9PnMoj88c6jcYITtB337jhGpbc9jtrPhsPqm9Cug7ZnxIx0aIcLSgUyp6rKYm+PZJD3YFLqDHvFzKzVbGxkeHtU93wmWZ/xLIl1IOA/4HvOOtkZTyNSnlGCnlmNTU1DCduiNKlYqU3F7s27ndZ5tAbk7L50zh+UtGBG2XC5aiJYuoryznhIsuj/ieRzhmiVUqabYHTnrXFfwVswnm2c1O0IclCZc37FJSabaQrfWeYz0cBCPQywH3UMacttfakVLWSinNbf++AYwOz/C6Ru6Q46jYvgWbxeL1fW9Lr0S7IKXKwk8fbOPb1zai/nU/s+OTOcuhZ6BFSZZWjVal4K75RV0qbmFobODXj96j17CR9BkTScYV4dglRqWg2dYzGnogZc29voEnApg8sGcUTYAaiw2bhCxdz3m2BWNyWQX0E0IU4BTklwKXuTcQQmRKKSvb/j0X2BLWUYZI7pDhrFn0BWWbN5A/ovPc0r70kjDAqmSsWUWG3Tm3bV9ZRXS8BqVagaXBzMBmGIIGDFChdLBeo2RLfedAokD8+J83sJpNTL7mhoh2HuGYJlappKWbGrovs4q/tLkFcxYRr1ejVor2EnPuSODTNeWM6ZXUI5Gi5SangnlINXQppQ24FViMU1B/JKXcJIR4VAhxbluz24UQm4QQ64DbgWt6asDB0GvoCLRR0Wz9dZnX92dPHUCGQsnlLVrONWhQSViqs/BRmo2UK3sTdV4uj9samKds4vl4I/+JMfGzzopawjSjhhuadAxqgae/DS4H+rbffmbLLz8y7vxLSI7kbYlwjBOnUtJk67pA92dW8eeVImnzN5f41NS7U9sgEGXmNoF+iDV0pJRfA197vPag299zgbnhHVrXUWk09Bt3Att++4XJ19yINiqqw/uj9XquatZicjj4OsrCJrXdud6y2Ln3k/UgDhSJdQioUkmqVDZWaG3k2hRMMKs41aihttRB1Z4m0vN9JwJr2l/Nd6+/TGbfAYyfeUlPXnaECEcEsSoFLXYHNodE1YWcJv7MKsF4qlkdkiiNigaD1asbY7h9z13sMTgFel4PCvSjM/QfGH76WVhNRtZ//22H18u21vHVS+tITI3im2zBJo29Q7SB1SG9LscAEFCqdvBRtIXPo83ohILPnlzDmm/3eK1EYjWZWPDU35BScuatd6NQ9myVlggRjgRSNU7teL/V1qXj/YX6e6bc8NdHoNKU4WaHwUSWVt1jeVzgKBboGX36kXfcMNYsWoDF5LwBGmsMfPvaRuJTozj/nlHsbOniTCygPErQ98q+9B6ZyooFxXz3783Y3eqPOhx2vv3nc9Ts3c3022eTmHl41iCMEOFgk65xGgaqLcEVj/EkkCB2eartnjedbD9tQ4kwDQc7DWb6RoUnC6wvjlqBDnDCxVfSWl/H759/hHRIlv5nKwDTbx6GLlrd5ZnYlWhr5vg8zrhuCOPO6832lVV8+dI6rBY70uHgu9dfZvvvyzn58j9QMHJMOC8rQoQjmvQ2Db3K3DWBHoog9tfWWwK9C0Zn89TibRTMWdQlbzZfSCnZaTDRNyq0/FKhctQl53Ine8AgBk+awuovP0cTNYCKHQYmXzmQ+FSnIJ89dQCzP16H1S2GXwH4cqhKjFJT+OAZHV4TQjDmzHxiErV8/84Wvn5lHVFRv7Fh6RLGz7yEMefM7KGrixDhyCS1zcuj2tI1k0soRSgCtXWPMA2mbGRXqbLYaLE7elxDP6oFOsApV11Hycb1/Pbxy6T1u55BJ2R2bOBhaFMqBQ4fNnRf+ZQBBo7PxGo2893rL+Cw7mDMOTM54eIrujv8CBGOOtLaTC5VXTS5QGhFKIJtG0zZyK6yrdUEQP/ontXQj2qTC4A+No7R59yI3dqIufFzbBZz+3tPLd7WaQPUapcoffiJC/C5BGuqqWbD//6Bw7oTlf5kopMnR/zNI0TwglahIE2jotTkPfDvUOHPh7275pf1zQYAjuvhWqpHvUAHqN8XR1TS2dRV7OSzxx/G2NIM+P4C7VJ63SGX0MlHVUrJ1l+X8Z/7bqO2rIRz7p7DoJPO4veFuynbVh/mK4kQ4eigQK9lt8EcuGGYWVBYzsR5S73ayP3tqXW34MW6ZgN5Og0JPZhpEY4Bge6wO9i9fj8DTjiJs277M5U7tvLB/fdQuXObzy8wO0HvM82m+yRQV1HOZ48/xKIXniQxM5srn/gH/cdNZPIVTjv90ne2YDF2zU4YIcLRTC+9hr0HWUP3F5AE3jdQ3elO0NGGZiPDYntWO4djQKDXlrdiNdnJGZjIoIknc+EDf8dmsfLB/bO5RqwlSZg6tHftgPtzd6reU8yiF5/i7bv/RMX2rSSfejGvxJzFyOfWMHHeUhZt3sdp1wympd7Er5/vOhiXGSHCEUW+Xkul2YqhB7MueuLPRg4dy0b6oitBRw1WG3tNFobHRgVu3E2O+k3Rqt2NAGT0dha+yBk4hKuffomf33+bDUuXcAW/UBrbm/WaAhzpvblj+oG6n64db4W0k2reT761khNbS3n3vjLUWh2jz55BTf4EHli8F6PVuXx0T8U7dHIOG34o47hJWaTkxB6aDyBChMOQAr3T22Ov0cygMNqVfeV4geBqj7o2UCfOWxqwrF2wrGt29jP0IGjoR71Ab6gxolIriE06sLusjYrmtOtuYcw5F7Bm0QK0vy4jp3o7VMPGrTEsV0ShVGu4Wq/C3NSA3tyEss2ZMbnvAAaedSODTjoFfUwsE+ct9Tnrf3/bSWz7fR+/fLyT8+4cEdkkjRChjT5t7nvbDaawCfRAbocJUWqvnmrehLS3FAJdzca4oqEFBTA6rufyoLs46gV6834TsSl6r8I0IT2DU/94E03DzuSVj5cS21hGorWBaFsrSqudKpuDtPQ8JozsT1pBX3IGDSE6IbFDH/5mfV20muPPKuCXj3dQubORrH4JPXGJESIccfSP0qEUsLnFxHmhVYr0iS+Typ3zi5j72XrMXlL2qpXCa0DSjJHZrN5bx3srStr307qajXFFYwtDY/XE9mDIv4ujXqCbjTZ0Ub4vc0FhOX/5YgtGkQYJne8sAQwaNIKTfHyBgSqODz4pizXf7mHt4r0RgR4hQhs6pYK+UTo2dTX9hhf82beNVu+2+miNyqdw/mFrTSfniFD90k12B2ubDFyTnRJU++5y1G+K2q12VBrfl+ltVnfHm6uiO4HCkNUaJcMm57B3Yy11la2hDT5ChKOYITF6NodRoHfFvt1o9B3cFIzNPRBFzQbMDskJCTEhj60rHPUC3eEA/Niug/lyXIEF3nxXveWDeHzm0A4z+OATsxEKwbYVlV56jxDh2GRwtI4Ks5W6LmZd9CSQ26E3/E0C4cjGuLy+BQE9WkfUnaAEuhBimhBimxBipxBijp92FwghpBDisMlGpdEpsZp83zDBfjm+fFehY3Y3b7UIo+I09Douma0r9uFwhLnceYQIRyguN77CJkNY+nMpV74ivT0JlFkxHNkYl9Y1MSI2isQeDihyEVCgCyGUwMvAmcBgYJYQYrCXdrHAHcDv4R5kd9DoVZgNvgV6V2b1rgQY9B+bjqHRQtXuppCOixDhaGVUXBRKAasau26KlFJiMBioqKhg69at9FbWMneshn7KGgoUteQoGkgQxnYvNYXA50rak2BW3/6otdhY22Tg1GTfBXDCTTDTxlhgp5SyGEAI8SFwHrDZo91fgSeA2WEdYTeJSdRSsqkWKaVXTxf3bGzeNjd9EchU4+kPe88pfRECSjbVktknPrSLiBDhKCRapWRItN6nQPfmUy4dDt78dhUq435yNEbSVQbsls4pBCZ6qTDXJLXk9+rFyaOHMGDAAHS6wImyQkkC5skPdU1I4LTDTKBnA6Vu/5cBHcrWCyFGAblSykVCiMNKoCekRWGzODA0WohO6Ji60vOGSdSrMDe1kGBuRm+zoLFb0TqsqO027EKBXaHAJpRYlSqikxOx1dejjItDeFQi8uYP+3+LNnN3ejwlm+sYd27vg3b9ESIczhwfH837lXVYHRK1Wzk692dIIKGpkgWfbyRLNDJG2HEoocGmZ4ctjhOH9ubEIb2Ij49Hq9WiVjul+TfryvjPLzswtTSRFWVnVAoYqkv4/PPtKJVK+vfvz/jx48nLy+uRGJHva5tIUasOSsi/i24bdoQQCuBZgigMLYS4AbgBIC8vr7unDoqEDKedrra8pV2gOywWlnz+Iz8tWMa5DfvIbakmo7WORHMzOnvw+SV2fPYIAIrYWFRpaagzM1FnZbF5h5HxilhKY9MojUnDrNJgtNopMhs5br8Ju82BUnXU70dHiBCQ4+OjebN8PxuaDYxy2zh8avE2pNXESFUVfZX7iRZWTFLFHnsiZY4EKh1xWHEqUrt2abnhok5WYC6dFM+lk4Z0eM3hcFBeXs7GjRtZv349W7ZsITc3l6lTp5KTkxO267I4HCyta2ZqShyKgxhQGIxALwfcS9XntL3mIhY4DvixbZbLABYKIc6VUq5270hK+RrwGsCYMWMOyu5gen4cQkDpymL0P35I6/JfMW3cSC+LhRsAg0pLaUwaW5PyqNPGYoxNpFEXS61DhVmpxqxUY1MoUUiJStpROuxoHDairUaS7CbGpqioKa9G21hH5ua9pKwq4nxzS/v5HQiqohIpiU2nMn0MjqTxLPh8LTMvHB2JHI1wzHNiYiwC+KGuuV2gNzU1kdOylcnaGpRIyhzx/G5PpcwRj8PLtl8oboQKhYLc3Fxyc3M59dRTKSoqYtmyZbzxxhuMHj2aqVOnotF0v4jzj3XNNNrsnJuWGLhxGAlGoK8C+gkhCnAK8kuBy1xvSikbgXaveSHEj8CfPYX5ocBWU0PTp58iZX+Kfm8m6aeX0Q0dSuLll3PXJgfbE3LZr4/v4NYo8B0s5I3PAOmx6a2xW8lorSW3uZq85ip6NVeR11zFoJ3/Y/XY8SS+9A82PrOLhFHD0Y8YQfS4seiGDEGojvo4rwgROpCicZokfqxr5vacZH799Vd++eUXBqhs7LIls8GeQZP0b7LoailJjUbD2LFjGT58OD/++CO//fYbJSUlXHzxxaSmhh7i784X1Q0kqpRMSjw4/ucuAkoQKaVNCHErsBhQAm9JKTcJIR4FVkspF/b0IEPFUlbG/pdepvGrr8Bmg1NeBqDX/34gKicdgL3zlrLfR4RnKDO+t2WGRammJC6DkrgMlru9rrdZubUFlvU6gdxaOKW4mJalS6kBFNHR6MeMJnrsOKInjEc7aFBEg49wTDAlKY4X91bxzCv/wlS3n0GDBkH2MD5evBejPBD0p1YKkHQoGRmOos5arZapU6fSr18/Pv30U9566y2uuOIKsrO7thlqsDv4dn8j56clolEcXNNqUCqhlPJr4GuP1x700faU7g+raziMRmr+8RJ1776LUChIvGwWibNmkWBP5LOn11JaLhnQZibzlnzHdXOE6vESLEaVGqOQlCUU8E5WDtfOm46tthbDypW0/v47hhW/U/3TMgBUaWnEnHwyMZMnEz1hPAq93m8muQgRjkQsFgva7Zuwa1MoiU/mz2efRe/eTqcBbXRcp/sdgqsl2hV69+7Ntddey3/+8x/eeecd/vjHP5KRkRFyP9/VNtFqdzAjPSEs4woFIeWhCXQZM2aMXL06fFYZ48ZNlN9zN9a9JcRfeAGpt92GOt2pjUuH5D//9yvJOTGcfcvw9mM8BeTkgan8sLWG8gYjAu/ad3e5vklLudLBumwVy+dM6fS+taqK1uW/0vLjj7T+8gsOgwGh1dI8eARvqXrzU9pgDGrnElOvVobkFxshwuFEXV0dH374Ifuqa3h/0tlMSU3ktaGH3gOssbGRN954AyEE119/PbGxoaW+nrVuF1tbTayeMDjoIKdQEEKskVJ6Dd48Koy2TV9/TcXcv6BMTiLv7beJHt/BqxKhEPQfm07h/0pprjO1p9L1V/FbQrtQj1IrMHhJ7tMvLZqyepPfXDCe2AG1wnuGNwB1ejoJM88nYeb5OCwWDKtW0fLjT1R//jW3t/zOnxRKVqUP4qfsEfyeMSgsBWwjRDjYFBcX8/HHHyOl5KorLsdk1/JJVT1GuwO9MvxmilBWt/Hx8Vx22WW8+eabLFy4kMsuuyxo82eJ0cyPdc3clZ/eI8I8EEe8QG/88ksq7r0P/ahR5Lz4AqrkZK/thkzKpnBJCRt+LOOEmX07ve8tSZdLQ0+M1jIyWc+K4nrs0llEeta4XP42Y2iHG0WnVvjM6uZCpVAwPDs2KCGs0GiImTiRmIkTmdA0jH4NpZxcVsSk8iJOqNyISanm18yhtE6OIWrs8RGbe4Qjgg0bNvDZZ5+RkpLCrFmzSEpK4py6Zv5TUcvSuiampyaE9XyB8qR7IzMzk9NOO41vv/2WdevWMWLEiKDO9X5lHQK4LNO7HOppjmhn6NYVK6iYM5eo448n7803fApzgLhkPb1HpLL5lwrMXpLc+9sILW8w8uuuOsb3TiQ7QY9DSn7YWsOCwvIOeVySorU++wCnxp+hV9M7O/TqRVmJUWxPzOP1oedy1dT7uffEP/F97mjGV22h5Oqr2TV1Gvv/9RrWquqQ+44Q4WCxZs0aPv30U/Ly8rj22mtJSkoCYEJCDElqJV9UN4T9nIFKz/li7NixZGdns3TpUmy2wAnErA7JB5W1TEmOI0fXfdfHrnDECnRbbS3lf56NJj+fnH/+E0UQYbyjz8rHbLBRuKSk03uBXJ8ksHxXXYckXXfNL+L+BRva2wTyjsmK12NutaGN9hKXHAD3nDNSKNiQ0oc3j7+Esjc+IevJJ1Cnp1Pz3HPsnDKF0ltupXXFCg7V/kiECN4oKiriyy+/pG/fvlx++eUdQu9VCsH5aYl8W9MYUvbFBYXlPjOhuvD1XPrLogpOn/UpU6bQ1NREYWFhwLF8UV1PlcXG1VmHRjuHI1ig1zz/AvaGBrKffRZlTHCpKVNzY+l3fDrrlpbSUt8x/0NXknRJ4L0VJe03gr9JQa9WcteEAhwOSXxq6H6zvhIFnTeuN/Hnnkuvd/9Dn2+/IfmPf8BYWEjJNX9g93kzaPjkExwmU8D+I0ToSbZv384XX3xBQUEBl156qdfgncuzkrFIyaf76oPq02VK8ZcJFXw/lwL/WVTB6fmSmZnJmjVr/I5FSsnLJdX0j9Id1GRcnhyRAt1SUkLDJ5+QdPll6Ab0D+nYcef2Rkr4ef72Dq/PGJnNBaOzCdUK7V4Aw9ekkBil5vGZQxnVluQ+KatruZEDpenV5OeTds899P1hKZl//zsIQeX9D7Bz8hSqX3gBW11dl84bIUJ3qK6u5uOPPyYjI4NLLrkElY8AusExekbGRvFeZW1Qq8tgTSnenktvXmzejhVCMHLkSPbt20dVVZXPsfxQ18yWVhO35KUd1FB/T45IgV7/4XxQKkn647UhHxufqmfs2QUUF9Wwa21He7O3klPB4FrSedOin79kBIUPnsGMkdlUbG9AqVKQnN2z0WMKrZaEC2ZSsOBz8t5+G/2oUdS++i92TjmVqsfnRezsEQ4aZrOZjz76CI1Gw6xZswJmOLw8K5mtrSZ+98jA6M20EmxFIW/Ppa/n3FufAwcOBGDXrl1ej5FS8o+SKjK1as4/BL7n7hxxXi5SSpqXLCHmxBNRp3etuuzw03LZuaaaH/67ldRescQlO5dkoUSIuuO+pPOVblNKScnmWrL6xaPW9HyxWHBqF9HjxxE9fhzm4mJq//Uadf/9L/Xvv0/8hReQfO11aHIiLo8Reo6vv/6a2tparrrqKuLiApsizk9P4LHiCl4prWZ824rWl5dKQpSaei8ODt5MLJ7P5cR5S/3WAnYnLi6O5ORk9uzZwwknnNDp/WX1LfzW0Mrf+mUf9MhQT444Dd1aXoG1rIzok07sch9KpYIzrhuCdEg+fKGQkx53zvz+lkrZCXom9knqZJIJNvS4pqSZ+n0GCoZ3L0dEV9H27k3WE/Po8+03xM+cScMnn7Jr2jT2PfootpqaQzKmCEc327dvZ926dZx00kkUFBQEdUy0Usk12Sks3t/E9lbn3o8v04qUBFVRyJt2H2o1oszMTKqrO69sHVLy9+IKcnRqrjyEm6EujjiBbtldDICubRnUVRLSooiZlI612sSIChtIsHux2+nVSp6/ZATL50zhvesn8NwlI7pUwWTjT+WoNAr6j8sIame+p9Dk5pL5yMP0/d8SEi66kPqPPmbnGVOpfuEF7C0tgTuIECEIzGYzX375JWlpaUyaNCmkY/+YnYpOIXipxGmz9rVybjRaO5hSEvRqdGoFd80van+ufG2cAiFVI0pOTqahoaGT++JXNY2sbzZyb0Em2kOsncMRaHK
Copy link
Member

Choose a reason for hiding this comment

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

From the annulus example, I was expecting the ellipses to encompass the data, but below that doesn't seem to be the case. I guess they're just some characteristic shape of a distribution, e.g. its std dev?

Warrants some comment in the notebook I think

pints/_nested/__init__.py Outdated Show resolved Hide resolved
pints/_nested/__init__.py Outdated Show resolved Hide resolved
pints/_nested/__init__.py Show resolved Hide resolved
pints/_nested/__init__.py Show resolved Hide resolved
@@ -414,6 +419,8 @@ def _initialise_logger(self):
self._logger.add_time('Time m:s')
self._logger.add_float('Delta_log(z)')
self._logger.add_float('Acceptance rate')
if self._sampler._multiple_ellipsoids:
Copy link
Member

Choose a reason for hiding this comment

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

In the other controllers we handle this in a more general way, by letting each sampler/optimiser define a method _log_init (that does nothing by default) which can update the logging

https://github.com/pints-team/pints/blob/master/pints/_mcmc/__init__.py#L623-L624

and _log_write:

https://github.com/pints-team/pints/blob/master/pints/_mcmc/__init__.py#L820-L821

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

2 participants